C++ codes for inferring the networks
发布时间：2018-01-21 点击次数：次
Source codes of inferring large-scale GRNs based on randomized singular value decomposition (IGRSVD)
Anjing Fan, Haitao Wang, Hua Xiang Xiufen Zou*
School of Mathematics and Statistics, Wuhan University, Wuhan, 430072, China
* Corresponding author: xfzou@whu.edu.cn.
We list C++ codes files for inferring gene regulatory networks using randomized algorithm based on singular value decomposition (IGRSVD). The detailed source codes and illustration of the codes are described as follows. Detailed function description has also been described in each major function. The users can copy and paste to obtain all the executable files
******************************************
If the readers are interested in all the source data and codes in this paper, they can contact with Prof.Xiufen Zou: xfzou@whu.edu.cn. We can provide all the executable files.
******************************************
------------------------The main function----------------------------
Function qr % QR factorization
Function SRFTsvd % SVD of small-scale matrix
Function solution % Get the solution
Function eval_fun % Evaluation the inferred network
Function main % Main code for inferring network
-------------------------------data-------------------------------------
We use the time series data Am*n with m time points and n genes.
------------------------------Input-------------------------------------
n_gene—the number of genes
n_time—the number of time points
FILE_NAME1—the time series data
FILE_NAME2—the gold standards data
For example, time series data with 2 time points and 3 genes is as follows:
4.5 6.7 6.4
5.4 3.5 6.7
------------------------------- Output------------------------------------
i_S1—the interaction network
i_SS—the interaction network with normalization
i_G—the interaction network by threshold q.
i—x13 and y13 for ROC.
All the experimental results are saved in the folder result.
******************************************
#include "stdafx.h"
#include "RSVD.h"
#include<math.h>
#include<stdlib.h>
#include<time.h>
#include<fstream>
//#include<dir.h>
#include<iostream>
using namespace std;
// global variable
int n_gene, n_time, n_k, n_c;
// generate a gauss random matrix.
double gaussrand()
{
static double V1,V2,S;
static int phase=0;
double X;
if(phase==0){
do{
double U1 = (double)rand()/RAND_MAX;
double U2 = (double)rand()/RAND_MAX;
V1 = 2 * U1 - 1;
V2 = 2 * U2 - 1;
S = V1 * V1 + V2 * V2;
}
while(S >= 1 || S == 0);
X= V1 * sqrt(-2 * log(S) / S);
}
else
X=V2*sqrt(-2*log(S)/S);
phase=1-phase;
return X;
}
bool mynull(double **A,double r,double *x)
{
int i,j,k,d;
double max,a;
for(i=0;i<r;i++)
{
d=i;
max=fabs(A[i][i]);
for(j=i+1;j<n_k;j++)
{
if(fabs(A[j][i])>max)
{
d=j;
max=fabs(A[j][i]);
}
}
for(j=i;j<n_k;j++)
{
a=A[i][j];
A[i][j]=A[d][j];
A[d][j]=a;
}
for(j=i+1;j<n_k;j++){
a=A[j][i]/A[i][i];
A[j][i]=0;
for(k=i+1;k<n_k;k++)
{
A[j][k]-=A[i][k]*a;
}
}
}
for(i=r;i<n_k;i++){
for(j=0;j<n_k;j++)
{
A[i][j]=0;
}
}
x[n_k-1]=1;
a=x[n_k-1]*x[n_k-1];
for(i=r-1;i>=0;i--)
{
x[i]=0;
for(j=i+1;j<n_k;j++){
x[i]-=x[j]*A[i][j];
}
x[i]=x[i]/A[i][i];
a+=x[i]*x[i];
}
a=sqrt(a);
for(i=0;i<n_k;i++)
x[i]/=a;
return true;
}
bool qr(double **Y,double **Q,double **R0)
{
int i,j,k,h,n_min=n_gene>n_k?n_k:n_gene;
double s,b,*u=new double[n_gene];
double **H=new double *[n_gene],**H0=new double *[n_gene],**H1=new double *[n_gene];
double **R=new double *[n_gene],**R1=new double *[n_gene];
cout<<"in qr"<<endl;
for(i=0;i<n_gene;i++)
{
H[i]=new double[n_gene];
H0[i]=new double[n_gene];
H1[i]=new double[n_gene];
R[i]=new double[n_k];
R1[i]=new double[n_k];
}
for(i=0;i<n_gene;i++)
{
for(j=0;j<n_gene;j++)
{
if(i==j)
H0[i][j]=1;
else
H0[i][j]=0;
}
}
for(i=0;i<n_gene;i++){
for(j=0;j<n_k;j++)
{
R[i][j]=Y[i][j];
}
}
for(i=0;i<n_min;i++)
{
s=0;
for(j=i;j<n_gene;j++)
{
u[j]=R[j][i];
s+=u[j]*u[j];
}
s=-(R[i][i]>0?1:-1)*sqrt(s);
u[i]-=s;
b=s*(s-R[i][i]);
for(j=0;j<n_gene;j++)
{
for(k=0;k<n_gene;k++)
{
if(j<i||k<i)
{
if(k==j)
H[j][k]=1;
else
H[j][k]=0;
}
else
{
if(k==j)
H[j][k]=1-u[j]*u[k]/b;
else
H[j][k]=-u[j]*u[k]/b;
}
}
}
for(j=0;j<n_gene;j++)
{
for(k=0;k<n_gene;k++)
{
H1[j][k]=H[j][0]*H0[0][k];
for(h=1;h<n_gene;h++)
{
H1[j][k]+=H[j][h]*H0[h][k];
}
}
}
for(j=0;j<n_gene;j++)
{
for(k=0;k<n_gene;k++)
{
H0[j][k]=H1[j][k];
}
}
R[i][i]=s;
for(j=i+1;j<n_gene;j++){
R[j][i]=0;
}
for(j=i;j<n_gene;j++){
for(k=i+1;k<n_k;k++)
{
R1[j][k]=H[j][i]*R[i][k];
for(h=i+1;h<n_gene;h++)
{
R1[j][k]+=H[j][h]*R[h][k];
}
}
}
for(j=i;j<n_gene;j++){
for(k=i+1;k<n_k;k++)
{
R[j][k]=R1[j][k];
}
}
}
for(i=0;i<n_gene;i++){
for(j=0;j<n_k;j++)
{
Q[i][j]=H0[j][i];
}
}
for(i=0;i<n_k;i++){
for(j=0;j<n_k;j++)
{
R0[i][j]=R[i][j];
}
}
for(i=0;i<n_gene;i++)
{
delete [] H[i];
delete [] H0[i];
delete [] H1[i];
delete [] R[i];
delete [] R1[i];
}
delete [] H;
delete [] H0;
delete [] H1;
delete [] R;
delete [] R1;
delete [] u;
return true;
}
bool qr1(double **Y,double **Q,double **R0)
{
int i,j,k,l=n_k;
double s,b,*u=new double[n_k];
double **H=new double *[n_k],**H0=new double *[n_k],**H1=new double *[n_k];
double **R=new double *[n_k],**R1=new double *[n_k];
for(i=0;i<n_k;i++)
{
H[i]=new double[n_k];
H0[i]=new double[n_k];
H1[i]=new double[n_k];
R[i]=new double[n_k];
R1[i]=new double[n_k];
}
for(i=0;i<n_k;i++)
{
for(j=0;j<n_k;j++)
{
if(i==j)
H0[i][j]=1;
else
H0[i][j]=0;
}
}
for(i=0;i<n_k;i++){
for(j=0;j<n_k;j++)
{
R[i][j]=Y[i][j];
}
}
for(i=0;i<l;i++)
{
s=0;
for(j=i;j<n_k;j++)
{
u[j]=R[j][i];
s+=u[j]*u[j];
}
s=-(R[i][i]>0?1:-1)*sqrt(s);
u[i]-=s;
b=s*(s-R[i][i]);
for(j=0;j<n_k;j++)
{
for(k=0;k<n_k;k++)
{
if(j<i||k<i)
{
if(k==j)
H[j][k]=1;
else
H[j][k]=0;
}
else
{
if(k==j)
H[j][k]=1-u[j]*u[k]/b;
else
H[j][k]=-u[j]*u[k]/b;
}
}
}
for(j=0;j<n_k;j++)
{
for(k=0;k<n_k;k++)
{
H1[j][k]=H[j][0]*H0[0][k];
for(l=1;l<n_k;l++)
{
H1[j][k]+=H[j][l]*H0[l][k];
}
}
}
for(j=0;j<n_k;j++)
{
for(k=0;k<n_k;k++)
{
H0[j][k]=H1[j][k];
}
}
R[i][i]=s;
for(j=i+1;j<n_k;j++){
R[j][i]=0;
}
for(j=i;j<n_k;j++){
for(k=i+1;k<n_k;k++)
{
R1[j][k]=H[j][i]*R[i][k];
for(l=i+1;l<n_k;l++)
{
R1[j][k]+=H[j][l]*R[l][k];
}
}
}
for(j=i;j<n_k;j++){
for(k=i+1;k<n_k;k++)
{
R[j][k]=R1[j][k];
}
}
}
for(i=0;i<n_k;i++){
for(j=0;j<n_k;j++)
{
Q[i][j]=H0[j][i];
}
}
for(i=0;i<n_k;i++){
for(j=0;j<n_k;j++)
{
R0[i][j]=R[i][j];
}
}
for(i=0;i<n_k;i++)
{
delete [] H[i];
delete [] H0[i];
delete [] H1[i];
delete [] R[i];
delete [] R1[i];
}
delete [] H;
delete [] H0;
delete [] H1;
delete [] R;
delete [] R1;
delete [] u;
return true;
}
bool svd(double **X,double **U,double *S,double **V)
{
int i,j,k,l;
double s,b,ss,*x;
double **A0,**A1,**H,**H0,**H1,**H2,*u;
double **A,**Q,**R;
x=(double *)calloc(n_k,sizeof(double));
A=(double **)calloc(n_k,sizeof(double *));
A0=(double **)calloc(n_k,sizeof(double *));
A1=(double **)calloc(n_k,sizeof(double *));
H=(double **)calloc(n_k,sizeof(double *));
H0=(double **)calloc(n_k,sizeof(double *));
H1=(double **)calloc(n_k,sizeof(double *));
H2=(double **)calloc(n_k,sizeof(double *));
Q=(double **)calloc(n_k,sizeof(double *));
R=(double **)calloc(n_k,sizeof(double *));
u=(double *)calloc(n_k,sizeof(double));
for(i=0;i<n_k;i++)
{
A[i]=(double *)calloc(n_k,sizeof(double));
A0[i]=(double *)calloc(n_k,sizeof(double));
A1[i]=(double *)calloc(n_k,sizeof(double));
H[i]=(double *)calloc(n_k,sizeof(double));
H0[i]=(double *)calloc(n_k,sizeof(double));
H1[i]=(double *)calloc(n_k,sizeof(double));
H2[i]=(double *)calloc(n_k,sizeof(double));
Q[i]=(double *)calloc(n_k,sizeof(double));
R[i]=(double *)calloc(n_k,sizeof(double));
}
cout<<"in svd"<<endl;
for(i=0;i<n_k;i++){
for(j=0;j<n_k;j++)
{
if(i==j)
H0[i][j]=1;
else
H0[i][j]=0;
A[i][j]=X[0][i]*X[0][j];
for(k=1;k<n_time;k++)
{
A[i][j]+=X[k][i]*X[k][j];
}
A0[i][j]=A[i][j];
}
}
cout<<"computing A is done."<<endl;
for(i=0;i<n_k-2;i++)
{
s=0;
for(j=i+1;j<n_k;j++)
{
u[j]=A[j][i];
s+=u[j]*u[j];
}
s=-(u[i+1]>0?1:-1)*sqrt(s);
b=s*(s-u[i+1]);
u[i+1]-=s;
for(j=0;j<n_k;j++){
for(k=0;k<n_k;k++)
{
if(j<=i||k<=i)
{
if(k==j)
H[j][k]=1;
else
H[j][k]=0;
}
else
{
if(k==j)
H[j][k]=1-u[j]*u[k]/b;
else
H[j][k]=-u[j]*u[k]/b;
}
}
}
for(j=0;j<n_k;j++)
{
for(k=0;k<n_k;k++)
{
H1[j][k]=H[j][0]*H0[0][k];
for(l=1;l<n_k;l++)
{
H1[j][k]+=H[j][l]*H0[l][k];
}
}
}
for(j=0;j<n_k;j++)
{
for(k=0;k<n_k;k++)
{
H0[j][k]=H1[j][k];
}
}
A[i+1][i]=s;
A[i][i+1]=s;
for(j=i+2;j<n_k;j++)
{
A[j][i]=0;
A[i][j]=0;
}
for(j=i+1;j<n_k;j++){
for(k=i+1;k<n_k;k++)
{
H1[j][k]=H[j][i+1]*A[i+1][k];
for(l=i+2;l<n_k;l++)
{
H1[j][k]+=H[j][l]*A[l][k];
}
}
}
for(j=i+1;j<n_k;j++){
for(k=i+1;k<n_k;k++)
{
H2[j][k]=H1[j][i+1]*H[i+1][k];
for(l=i+2;l<n_k;l++)
{
H2[j][k]+=H1[j][l]*H[l][k];
}
}
}
for(j=i+1;j<n_k;j++){
for(k=i+1;k<n_k;k++)
{
A[j][k]=H2[j][k];
}
}
}
A[n_k-1][n_k-2]=A[n_k-2][n_k-1];
cout<<" Three diagonalization is done. "<<endl;
ss=1;
do{
ss=0;
qr1(A,Q,R);
for(i=0;i<n_k;i++){
for(j=0;j<n_k;j++)
{
A[i][j]=R[i][0]*Q[0][j];
for(k=1;k<n_k;k++)
A[i][j]+=R[i][k]*Q[k][j];
if(i!=j)
ss+=fabs(A[i][j]);
}
}
}while(ss>1e-9);
cout<<" The calculated eigenvalues are completed. "<<endl;
for(i=0;i<n_k;i++){
for(j=0;j<n_k;j++)
{
if(i!=j)
A[i][j]=0;
else
A[i][j]=A[i][j]>0?A[i][j]:-A[i][j];
}
S[i]=sqrt(A[i][i]);
}
for(i=0;i<n_k;i++)
{
for(j=0;j<n_k;j++){
for(k=0;k<n_k;k++)
{
if(j==k)
A1[j][k]=A0[j][k]-A[i][i];
else
A1[j][k]=A0[j][k];
}
}
mynull(A1,n_k-1,x);
for(j=0;j<n_k;j++)
{
V[j][i]=x[j];
}
}
for(i=0;i<n_time;i++){
for(j=0;j<n_k;j++)
{
U[i][j]=X[i][0]*V[0][j];
for(k=1;k<n_k;k++)
U[i][j]+=X[i][k]*V[k][j];
//U[i][j]/=S[j];
}
}
for(i=0;i<n_k;i++)
{
if(S[i]>0){
for(j=0;j<n_time;j++){
U[j][i]/=S[i];
}
}
}
for(i=0;i<n_k;i++)
{
free(A[i]);
free(A0[i]);
free(A1[i]);
free(H[i]);
free(H0[i]);
free(H1[i]);
free(H2[i]);
free(Q[i]);
free(R[i]);
}
free(A);
free(A0);
free(A1);
free(H);
free(H0);
free(H1);
free(H2);
free(Q);
free(R);
return true;
}
bool solution(double **U,double *s,double **V,double *b,double *S1)
{
int i,j;
double *beta=new double[n_k],*zeta=new double[n_k];
for(i=0;i<n_k;i++)
{
beta[i]=U[0][i]*b[0];
for(j=0;j<n_time;j++)
{
beta[i]+=U[j][i]*b[j];
}
}
//compute zeta
for(i=0;i<n_k;i++)
{
zeta[i]=s[i]*beta[i];
}
for(i=0;i<n_gene;i++)
{
S1[i]=V[i][0]*zeta[0]/(s[0]*s[0]);
for(j=1;j<n_k;j++)
{
S1[i]+=V[i][j]*zeta[j]/(s[j]*s[j]);
}
}
delete [] zeta;
delete [] beta;
return true;
}
bool SRFTsvd(double **A,double k,double **U,double *s,double **V)
{
int i,j,g,h;
double **Omega=new double *[n_k],**Y=new double *[n_gene],**Q=new double *[n_gene],**R=new double *[n_k];
double **B=new double *[n_time],**W=new double *[n_k];
for(i=0;i<n_k;i++)
{
Omega[i]=new double[n_time];
R[i]=new double[n_k];
W[i]=new double[n_k];
}
for(i=0;i<n_gene;i++)
{
Y[i]=new double[n_k];
Q[i]=new double[n_k];
}
for(i=0;i<n_time;i++)
{
B[i]=new double[n_k];
}
//randomly produce Omega
for(i=0;i<n_k;i++){
for(h=0;h<n_time;h++)
{
Omega[i][h]=gaussrand();
}
}
//Y = Omega * A
for(i=0;i<n_k;i++){
for(h=0;h<n_gene;h++){
Y[h][i]=Omega[i][0]*A[0][h];
for(g=1;g<n_time;g++)
{
Y[h][i]+=Omega[i][g]*A[g][h];
}
}
}
//QR decompose for Y
cout<<"qr"<<endl;
qr(Y,Q,R);
cout<<"hqr"<<endl;
//compute B
for(i=0;i<n_time;i++){
for(j=0;j<n_k;j++){
B[i][j]=A[i][0]*Q[0][j];
for(h=1;h<n_gene;h++)
{
B[i][j]+=A[i][h]*Q[h][j];
}
}
}
cout<<"svd"<<endl;
//svd decompose for B
svd(B,U,s,W);
cout<<"hsvd"<<endl;
//compute V
for(i=0;i<n_gene;i++){
for(j=0;j<n_k;j++){
V[i][j]=Q[i][0]*W[0][j];
for(h=1;h<n_k;h++)
{
V[i][j]+=Q[i][h]*W[h][j];
}
}
}
for(i=0;i<n_k;i++)
{
delete [] Omega[i];
delete [] R[i];
delete [] W[i];
}
for(i=0;i<n_gene;i++)
{
delete [] Y[i];
delete [] Q[i];
}
for(i=0;i<n_time;i++)
{
delete [] B[i];
}
delete [] Omega;
delete [] R;
delete [] W;
delete [] Y;
delete [] Q;
delete [] B;
return true;
}
bool saureus_SVD01(double **A,double **B,double k,double **S1)
{
int i,j;
double *b=new double[n_time],*x1=new double[n_gene],*s=new double[n_k];
double **U,**V;
U=new double *[n_time];
V=new double *[n_gene];
for(i=0;i<n_time;i++)
{
U[i]=new double[n_k];
}
for(i=0;i<n_gene;i++)
{
x1[i]=0;
V[i]=new double[n_k];
}
cout<<"before SRFTsvd"<<endl;
SRFTsvd(A,k,U,s,V);
for(i=0;i<n_k;i++){
cout<<s[i]<<endl;
}
/* for(i=0;i<n_time;i++){
for(j=0;j<n_k;j++){
cout<<U[i][j]<<"t";
}
cout<<endl;
}
for(i=0;i<n_k;i++){
cout<<s[i]<<endl;
}
for(i=0;i<n_gene;i++){
for(j=0;j<n_k;j++){
cout<<V[i][j]<<"t";
}
cout<<endl;
}*/
cout<<"after SRFTsvd"<<endl;
for(j=0;j<n_gene;j++)
{
for(i=0;i<n_time;i++){
b[i]=B[i][j];
}
// cout<<"before solution"<<endl;
solution(U,s,V,b,x1);
for(i=0;i<n_gene;i++){
S1[i][j]=x1[i];
}
}
for(i=0;i<n_time;i++)
{
delete [] U[i];
}
for(i=0;i<n_gene;i++)
{
delete [] V[i];
}
delete [] U;
delete [] V;
delete [] b;
delete [] x1;
delete [] s;
return true;
}
bool rsvd(double **data,int k,int c,double l,bool **G,double **SS,int itertime)
{
int i,j,h;
char outfilename[100]="",str[10];
double min,max;
double **S_cen_diff,**S1;
S_cen_diff=(double **)calloc(n_time,sizeof(double *));
S1=(double **)calloc(n_gene,sizeof(double *));
ofstream fp;
//Calculate difference
//fp.open("S_cen_diff.txt");
cout<<"in rsvd"<<endl;
for(i=0;i<n_time;i++)
{
S_cen_diff[i]=(double *)calloc(n_gene,sizeof(double));
}
for(i=0;i<n_gene;i++)
{
S1[i]=(double *)calloc(n_gene,sizeof(double));
}
cout<<"before saureus"<<endl;
for(i=0;i<n_time;i++){
for(j=0;j<n_gene;j++)
{
if(i==0)
S_cen_diff[0][j]=0;
else
S_cen_diff[i][j]=data[i-1][j]-data[i][j];
}
}
for(i=0;i<c;i++)
{
cout<<"before saureus_SVD01"<<endl;
saureus_SVD01(data,S_cen_diff,k,S1);
for(j=0;j<n_gene;j++)
{
for(h=0;h<n_gene;h++)
{
SS[j][h]+=S1[j][h];
//SS[j][h]+=SS[j][h];
}
}
}
strcpy(outfilename,"result\");
itoa(itertime,str,10);
strcat(outfilename,str);
strcat(outfilename,"_S1.txt");
cout<<outfilename<<endl;
fp.open(outfilename);
for(j=0;j<n_gene;j++)
{
for(h=0;h<n_gene;h++)
{
SS[j][h]/=c;
fp<<SS[j][h]<<"t";
}
fp<<endl;
}
fp.close();
for(h=0;h<n_gene;h++){
min=SS[0][h];
max=min;
for(j=1;j<n_gene;j++)
{
if(SS[j][h]>max)
max=SS[j][h];
if(SS[j][h]<min)
min=SS[j][h];
}
for(j=0;j<n_gene;j++)
{
SS[j][h]=(SS[j][h]-min)/(max-min)*2-1;
if(SS[j][h]>=l||SS[j][h]<=-l)
G[j][h]=true;
else
G[j][h]=false;
}
}
cout<<"l="<<l<<endl;
strcpy(outfilename,"result\");
itoa(itertime,str,10);
strcat(outfilename,str);
strcat(outfilename,"_SS.txt");
fp.open(outfilename);
for(j=0;j<n_gene;j++)
{
for(h=0;h<n_gene;h++)
{
fp<<SS[j][h]<<"t";
}
fp<<endl;
}
fp.close();
strcpy(outfilename,"result\");
itoa(itertime,str,10);
strcat(outfilename,str);
strcat(outfilename,"_G.txt");
fp.open(outfilename);
for(j=0;j<n_gene;j++)
{
for(h=0;h<n_gene;h++)
{
fp<<G[j][h]<<"t";;
}
fp<<endl;
}
fp.close();
for(i=0;i<n_time;i++)
{
free(S_cen_diff[i]);
}
for(i=0;i<n_gene;i++)
{
free(S1[i]);
}
free(S_cen_diff);
free(S1);
return true;
}
bool eval_fun(bool **G,bool **Sdd,double &tpr,double &fpr)
{
int i,j;
int sd1=0,sd0=0,tp=0,fn=0,tn=0,fp=0;
for(i=0;i<n_gene;i++)
{
for(j=0;j<n_gene;j++)
{
sd1+=Sdd[i][j];
tp+=G[i][j]&&Sdd[i][j];
fn+=Sdd[i][j]&&(!G[i][j]);
tn+=!(G[i][j]||Sdd[i][j]);
fp+=G[i][j]&&(!Sdd[i][j]);
}
}
sd0=n_gene*n_gene-sd1;
tpr=1.0*tp/(tp+fn);
fpr=1.0*fp/(fp+tn);
return true;
}
int main(int argc, char* argv[])
{
int i,j,k,itertime,n_itertimes;
double n_l;
bool **G,**Sdd,**S13;
double **data,**SS;
ifstream fin;
ofstream fout;
clock_t t1,t2;
//tp,fp,tn,fn
int sd1=0,sd0=0,tp=0,fn=0,tn=0,fp=0;
double d,tpr,fpr,ppv,acc,mcc,spe,E_acc,E_error;
double x13[101],y13[101],AUC13=0;
char outfilename[100]="",str[10];
// {
char *FILE_NAME1="data92.txt",*FILE_NAME2="Sdd92.txt",*DIR_NAME="result";
n_itertimes=10;//
n_time=53;
n_gene=92;
n_k=10;
n_c=10;
n_l=0.1;
// }
G=(bool **)calloc(n_gene,sizeof(bool *));
SS=(double **)calloc(n_gene,sizeof(double *));
Sdd=(bool **)calloc(n_gene,sizeof(bool *));
S13=(bool **)calloc(n_gene,sizeof(bool *));
data=(double **)calloc(n_time,sizeof(double *));
for(i=0;i<n_gene;i++)
{
Sdd[i]=(bool *)calloc(n_gene,sizeof(bool));
G[i]=(bool *)calloc(n_gene,sizeof(bool));
S13[i]=(bool *)calloc(n_gene,sizeof(bool));
SS[i]=(double *)calloc(n_gene,sizeof(double));
}
for(i=0;i<n_time;i++)
{
data[i]=(double *)calloc(n_gene,sizeof(double));
}
fin.open(FILE_NAME1);
for(i=0;i<n_time;i++){
for(j=0;j<n_gene;j++)
{
fin>>data[i][j];
}
}
fin.close();
fin.open(FILE_NAME2);
for(i=0;i<n_gene;i++){
for(j=0;j<n_gene;j++)
{
fin>>Sdd[i][j];
}
}
fin.close();
system("md result");
/**********************begin n_itertimes**********************/
for(itertime=1;itertime<=n_itertimes;itertime++)
{
/*****************************************************************/
srand((unsigned)time(NULL));
cout<<"before rsvd"<<endl;
t1=clock();
rsvd(data,n_k,n_c,n_l,G,SS,itertime);
t2=clock();
for(i=0;i<n_gene;i++)
{
for(j=0;j<n_gene;j++)
{
sd1+=Sdd[i][j];
tp+=G[i][j]&&Sdd[i][j];
fn+=Sdd[i][j]&&(!G[i][j]);
tn+=!(G[i][j]||Sdd[i][j]);
fp+=G[i][j]&&(!Sdd[i][j]);
}
}
cout<<"tp="<<tp<<"fn="<<fn<<"tn="<<tn<<"fp="<<fp<<endl;
sd0=n_gene*n_gene-sd1;
tpr=1.0*tp/(tp+fn);
fpr=1.0*fp/(fp+tn);
ppv=1.0*tp/(tp+fp);
acc=1.0*(tp+tn)/(tp+fp+tn+fn);
//cout<<"mcc"<<endl;
//cout<<(tp*tn-fp*fn)<<endl;
//cout<<(tp+fp)<<"t"<<(tp+fn)<<"t"<<(tn+fp)<<"t"<<(tn+fn)<<endl;
//cout<<(tp+fp)*(tp+fn)*(tn+fp)*(tn+fn)<<endl;
mcc=(tp*tn-fp*fn)/sqrt(1.0*(tp+fp)*(tp+fn)*(tn+fp)*(tn+fn));
//cout<<"mcc="<<mcc<<endl;
spe=1.0*tn/(tn+fp);
E_acc=tp*sd0*1.0/(n_gene*n_gene)+tn*sd1*1.0/(n_gene*n_gene);
E_error=fp*sd0*1.0/(n_gene*n_gene)+fn*sd1*1.0/(n_gene*n_gene);
for(i=0;i<=100;i++)
{
d=0.01*i;
for(j=0;j<n_gene;j++)
{
for(k=0;k<n_gene;k++)
{
if(fabs(SS[j][k])>=d)
S13[j][k]=true;
else
S13[j][k]=false;
}
}
eval_fun(S13,Sdd,y13[100-i],x13[100-i]);
}
for(i=0;i<100;i++)
{
AUC13+=(x13[i+1]-x13[i])*(y13[i]+y13[i+1])/2;
}
strcpy(outfilename,"result\");
itoa(itertime,str,10);
strcat(outfilename,str);
strcat(outfilename,".txt");
fout.open(outfilename);
fout<<"Result:"<<endl;
fout<<"tpt"<<tp<<endl;
fout<<"fpt"<<fp<<endl;
fout<<"tnt"<<tn<<endl;
fout<<"fnt"<<fn<<endl;
fout<<"tprt"<<tpr<<endl;
fout<<"fprt"<<fpr<<endl;
fout<<"ppvt"<<ppv<<endl;
fout<<"acct"<<acc<<endl;
fout<<"mcct"<<mcc<<endl;
fout<<"spet"<<spe<<endl;
fout<<"E_acct"<<E_acc<<endl;
fout<<"E_errort"<<E_error<<endl;
fout<<"x13ty13"<<endl;
for(i=0;i<100;i++)
{
fout<<x13[i]<<"t"<<y13[i]<<endl;
}
fout<<"AUC13t"<<AUC13<<endl;
fout<<"TimeCost"<<(t2-t1)*1.0/1000<<endl;
fout.close();
/*****************************************************************/
}
/*********************end n_itertimes***********************/
for(i=0;i<n_time;i++)
{
free(data[i]);
}
free(data);
for(i=0;i<n_gene;i++)
{
free(SS[i]);
free(S13[i]);
free(Sdd[i]);
free(G[i]);
}
free(SS);
free(S13);
free(G);
free(Sdd);
return 0;
}