www.pudn.com > Genecircus20070919.rar > calculate_mi.cpp
// calculate_mi.cpp: implementation of the calculate_mi class.
//
//////////////////////////////////////////////////////////////////////
#include "calculate_mi.h"
#include "mainframe.h"
using namespace boost;
#define TWO_SQRTPI 1.77245385090551602729816
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
calculate_mi::calculate_mi()
{
prob_table = NULL;
norm_1D_table = NULL;
norm_2D_table = NULL;
name_for_now_gene = "";
law = 1;
step = 1;
pvalue = 0.0;//0.066027;//0.42875;
ok_readandwrite_index = true;
sigma = 0.141885;
other_name_index = -100;
}
calculate_mi::~calculate_mi()
{
if(prob_table != NULL){
for ( int i = 0; i < number; i++ )
{
delete[] prob_table[i];
delete[] norm_2D_table[i];
}
delete[] prob_table;
delete[] norm_2D_table;
delete[] norm_1D_table;
}
delete Integr;
}
double calculate_mi::get_GenePairs_Mi_and_save(int i ,int j)
{
number = this->number;
if(gene_matrix_for_calcu_mi->gene_pairs.size() != 0){
gene_matrix_for_calcu_mi->gene_pairs.clear();
}
gene_matrix_for_calcu_mi->resturct_genepairs(i,j);
double mi = 0.0;
rank_transformation(gene_matrix_for_calcu_mi->gene_pairs);
for (unsigned int a = 0; a < this->number; a++ )
{
double v = calculateprob( gene_matrix_for_calcu_mi->gene_pairs, a );
mi += log( v );
}
mi = mi/(this->number);
return mi;
}
void calculate_mi::get_GenePairs_Mi(int i ,int j,bool p)
{
number = this->number;
if(gene_matrix_for_calcu_mi->gene_pairs.size() != 0){
gene_matrix_for_calcu_mi->gene_pairs.clear();
}
gene_matrix_for_calcu_mi->resturct_genepairs(i,j);
double mi = 0.0;
rank_transformation(gene_matrix_for_calcu_mi->gene_pairs);
for (unsigned int a = 0; a number; a++ )
{
double v = calculateprob( gene_matrix_for_calcu_mi->gene_pairs, a );
mi += log( v );
}
mi = mi/(this->number);
if(p)
setmimatrix(i,j,gene_matrix_for_calcu_mi->micro_set[i].name,gene_matrix_for_calcu_mi->micro_set[j].name,mi);
else{
setmimatrix(i,j,gene_matrix_for_calcu_mi->micro_set[i].name,gene_matrix_for_calcu_mi->micro_set[j].name,mi,false);
}
}
void calculate_mi::setmimatrix(int i,int j,wxString name1,wxString name2,double mi,bool p)
{
if(p){
if(i != j){
gene_matrix_for_calcu_mi->mi_matrix->utype = 1;
gene_matrix_for_calcu_mi->mi_matrix->mi = mi;
gene_matrix_for_calcu_mi->mi_matrix->name = name1;
connect_matrix * mi_matrix_temp = new connect_matrix();
gene_matrix_for_calcu_mi->mi_matrix->lnext = mi_matrix_temp;
gene_matrix_for_calcu_mi->mi_matrix = mi_matrix_temp;
if(mi > pvalue && j != this->other_name_index){
gene_matrix_for_calcu_mi->index_temp.push_back(j);//in order to rebuild networks,
//save nodes index,
//j is center gene's round gene
}
}
else{
gene_matrix_for_calcu_mi->mi_matrix->utype = 1;
gene_matrix_for_calcu_mi->mi_matrix->mi = 0.0;
gene_matrix_for_calcu_mi->mi_matrix->name = name1;
connect_matrix * mi_matrix_temp = new connect_matrix();
gene_matrix_for_calcu_mi->mi_matrix->lnext = mi_matrix_temp;
gene_matrix_for_calcu_mi->mi_matrix = mi_matrix_temp;
}
}
else{
if(i != j){//gene_matrix_for_calcu_mi->pvalue[j]
gene_matrix_for_calcu_mi->pvalue[j].push_back(mi);
}
else{
gene_matrix_for_calcu_mi->pvalue[j].push_back(0.0);
}
}
}
void calculate_mi::Getresult(int obj_index,double max )
{
if(law == false){
for(int i = 0;i< gene_matrix_for_calcu_mi->micro_set.size();i++){
vector::const_iterator it = max_element(gene_matrix_for_calcu_mi->pvalue[i].begin(),gene_matrix_for_calcu_mi->pvalue[i].end());
double temp = *it;
if(max < temp){
max = temp;
}
}
}
else{
for(int i = 1;i< gene_matrix_for_calcu_mi->micro_set.size();i++){
for(int j = 0; jpvalue[i].size();j++){
gene_matrix_for_calcu_mi->pvalue[0].push_back(gene_matrix_for_calcu_mi->pvalue[i][j]);
}
}
sort(gene_matrix_for_calcu_mi->pvalue[0].begin(),gene_matrix_for_calcu_mi->pvalue[0].end());
int size = gene_matrix_for_calcu_mi->pvalue[0].size();
int size1 = size-50;
p p1;
for(int t = size1;t < size;t++){
p1.push_back(gene_matrix_for_calcu_mi->pvalue[0][t]);
}
if(dirname != "")
save_to_pvalue_temp(size-size1,p1,name_for_now_gene,dirname);
else
save_to_pvalue_temp(size-size1,p1,name_for_now_gene);
max = gene_matrix_for_calcu_mi->pvalue[0][size-10];
pvalue = max;
}
string temp1,temp2;
int count;
if(law){
using boost::lexical_cast;
temp1 = "Max = " + boost::lexical_cast(max) + '\r'+'\n';
gene_matrix_for_calcu_mi->mi_matrix = gene_matrix_for_calcu_mi->head_mi_matrix;
for(int i = 0;i< gene_matrix_for_calcu_mi->micro_set.size();i++){
if(gene_matrix_for_calcu_mi->mi_matrix->mi > max){
double tr = gene_matrix_for_calcu_mi->mi_matrix->mi;
using boost::lexical_cast;
temp2 = boost::lexical_cast(tr);
temp1 += temp2 +'\t' + (string)(gene_matrix_for_calcu_mi->micro_set1[i].name) + '\r'+'\n';
}
gene_matrix_for_calcu_mi->mi_matrix = gene_matrix_for_calcu_mi->mi_matrix->lnext;
}
}
else{
}
count = temp1.length();
const char * t = temp1.c_str();
wxString temp_dir = dirname;
dirname = dirname + gene_matrix_for_calcu_mi->micro_set[obj_index].name+".tmp";
wxFile file;
file.Create(dirname,true);
file.Open(dirname,wxFile::read_write);
file.Write(t,count);
file.Close();
dirname = temp_dir;
}
void calculate_mi::save_to_pvalue_temp(int size,p p1,wxString name,wxString dirname)
{
string temp1,temp2;
for(int i =0;i(p1[i]);
temp1 = temp1 + temp2 + '\r'+'\n';
}
if(dirname != ""){
int count = temp1.length();
const char * t = temp1.c_str();
wxString temp_dir = dirname;
dirname = dirname + "_pvalue" + name +".tmp";
wxFile file;
file.Create(dirname,true);
file.Open(dirname,wxFile::read_write);
file.Write(t,count);
file.Close();
dirname = temp_dir;
}
}
double calculate_mi::calculateprob(Pair_Vector pairs,int i)
{
double fxy = 0.0;
double fx = 0.0;
double fy = 0.0;
int ix;
int iy;
unsigned int size = pairs.size();
double qd = 2*sigma*sigma;
for ( unsigned int j = 0; j < size; j++ )
{
int mu_x = pairs[j].xi;
int mu_y = pairs[j].yi;
//int ix = abs( pairs[i].xi - mu_x );
//int iy = abs( pairs[i].yi - mu_y );
if(pairs[i].xi <= mu_x){
ix = -(pairs[i].xi - mu_x);
}
else{
ix = pairs[i].xi - mu_x;
}
if(pairs[i].yi <= mu_y ){
iy = -( pairs[i].yi - mu_y );
}
else{
iy = pairs[i].yi - mu_y;
}
double dx = pairs[i].x - pairs[j].x;
double dy = pairs[i].y - pairs[j].y;
if ( prob_table[ix][iy] == -1.0 )
{
if ( prob_table[ix][0] == -1.0 )
{
prob_table[ix][0] = exp(-(dx*dx)/qd);
prob_table[0][ix] = prob_table[ix][0];
}
if ( prob_table[0][iy] == -1.0 )
{
prob_table[0][iy] = exp(-(dy*dy)/qd);
prob_table[iy][0] = prob_table[0][iy];
}
prob_table[ix][iy] = prob_table[ix][0]*prob_table[0][iy];
prob_table[iy][ix] = prob_table[ix] [iy];
}
fx += prob_table[ix][0]/norm_1D_table[mu_x];//mu_x是高斯核的中心
fy += prob_table[0][iy]/norm_1D_table[mu_y];
fxy += prob_table[ix][iy]/norm_2D_table[mu_x][mu_y];
/*
//for(unsigned int j = 0;jqgaus(up,low);
ss = 2*ss*(1/TWO_SQRTPI);
return ss;
}