diff options
| author | ziejd2 | 2017-09-14 15:39:41 -0500 |
|---|---|---|
| committer | ziejd2 | 2017-09-14 15:57:50 -0500 |
| commit | c4f926438dcb8abe805e910399940f79ff643c4b (patch) | |
| tree | 3146d05ada5cf4b48d9bdd16c1baa498e8df7192 /sourcecodes/localscore/network_score.c | |
| parent | 57ebf49403b75dcf8482d174b59f7fd2a961d98e (diff) | |
| download | BNW-c4f926438dcb8abe805e910399940f79ff643c4b.tar.gz | |
Add files via upload
Diffstat (limited to 'sourcecodes/localscore/network_score.c')
| -rw-r--r-- | sourcecodes/localscore/network_score.c | 1441 |
1 files changed, 1441 insertions, 0 deletions
diff --git a/sourcecodes/localscore/network_score.c b/sourcecodes/localscore/network_score.c new file mode 100644 index 00000000..c50a520e --- /dev/null +++ b/sourcecodes/localscore/network_score.c @@ -0,0 +1,1441 @@ +#include<stdio.h> +#include <stdlib.h> +#include<string.h> +#define MATHLIB_STANDALONE 1 +#include "matrix.h" +#include "postc0.c" +#include "modified_postc.c" // only line that has been changed in this code +//#define OUT_DIR "./Result/" //WHERE BENE READ DATA FOR STRUCTURE LEARNING //change it to #define OUT_DIR "/var/www/html/compbio/BNW/bene-0.9-4/example/resdir/" for server. Also change from sprintf(line+9,"%d",i); to sprintf(line+52,"%d",i); + +struct white{ +int n; +int *list; +}; + +struct ban{ +int n; +int *list; +}; + +struct ndata{ +int idx,type,score_i; +char **ddata; +char **label; //for descrete node +int lindex; +double *prob,*cdata,*post_alpha; +double score_zeroparent; +double *score; +struct white wlist; +struct ban blist; +}; + + +struct nd{ +int num,nd,*discrete,nc,*continuous,row; +double score; +struct ndata *node; +}; + +void get_next_seq(int *,int *,int); +void learn_node(struct nd *); +double learn_parent(struct nd *,int,int); + +double run_strtof (const char * input) //Convert string to double +{ + double output; + char * end; + output = strtod (input, & end); + if (end == input) { + return 0.0; + } + else { + return output; + } +} + + +char* itoa(int val, int base){ //convert number to string + + static char buf[32] = {0}; + + int i = 30; + + for(; val && i ; --i, val /= base) + + buf[i] = "0123456789abcdef"[val % base]; + + return &buf[i+1]; + +} + +double standarizedata(double val,double mean,double std) +{ + return (val-mean)/std; + +} + +double mean(double *a,int l) +{ +int i; +double m; +m=0; +for(i=0;i<l;i++) +{ + m+=a[i]; +} +m=m/l; +return m; +} + + +double stdev(double *a,int l,double m) +{ +int i; +double sqsum=0.0; + +for(i=0;i<l;i++) +{ + sqsum+=(a[i]-m)*(a[i]-m); + +} + +return sqrt(sqsum/(l-1)); +} + + +main(int argc, char *argv[]) +{ + FILE * pFile,*banf,*whitef; + int i,s_i,j,k,l,m,id,ic,c,row,imatch,flag,bit,bit_count,ccount,maxccount; + int n = 0,ban_flag,white_flag,white_flag_all; + char line[2500],nf[250],nt[250],*file_name; + char **name; //name of variables + double tc,tcm,score_val; + struct nd network; + int MAX_PARENT; + FILE **inputFiles; + + + double min_score; + + + if(argc<=5) + { + printf("Use 5 arguments: 1. input data file name. 2. ban list file name. 3. white list file name. 4. integer value of the maximum number of parents are allowed and 5. output file directory\n"); + return; + } + + + + pFile=fopen (argv[1],"r"); + + + if (pFile==NULL) perror ("Error opening file"); + else + {//reading data from input file + + + row=0; + do{ row++; + } while (fgets (line,2500,pFile)!=NULL); + + row-=3; + network.row=row; + fclose(pFile); + pFile=fopen (argv[1],"r"); + ccount=0; + maxccount=0; + do { + c = fgetc (pFile); + if (c == '\t') + { + n++; + if(maxccount<ccount) + maxccount=ccount; + ccount=0; + } + else + ccount++; + + + } while (c != '\n'); + + + network.num=n+1; + network.discrete=(int *)malloc(sizeof(int)*network.num); + network.continuous=(int *)malloc(sizeof(int)*network.num); + network.node=(struct ndata *)malloc(sizeof(struct ndata)*network.num); + + MAX_PARENT=atoi(argv[4]); //maximum number of parents for each node + if(MAX_PARENT==0) //no restriction on parents + MAX_PARENT=network.num+1; + + name=(char **)malloc(sizeof(char *)*network.num); + + for(i=0;i<network.num;i++) + name[i]=(char *)malloc(sizeof(char)*(maxccount+1)); + + + + id=0; + ic=0; + for(i=0;i<=n;i++) + { + fscanf(pFile,"%d",&c); + network.node[i].idx=i; + network.node[i].type=c; + if(c==1) + { + network.continuous[ic]=i; + ic++; + network.node[i].prob=(double *)malloc(sizeof(double)*2); //probability array of size 2 for continuous node + network.node[i].prob[0]=0.0; //initialization + network.node[i].prob[1]=0.0; + network.node[i].cdata=(double *)malloc(sizeof(double)*row); //data matrix for continuous node + } + else + { + network.discrete[id]=i; + network.node[i].post_alpha=(double *)calloc(network.node[i].type,sizeof(double)); //prepare total count of each discrete data lebel initialize with zero by calloc which is used in learnnode + id++; + network.node[i].prob=(double *)malloc(sizeof(double)*c); //probability array of size equals to number of levels for discrete node + network.node[i].ddata=(char **)malloc(sizeof(char *)*row); //data matrix for continuous node + network.node[i].label=(char **)malloc(sizeof(char *)*c); //label array + for(j=0;j<c;j++) + network.node[i].prob[j]=1.0/c; + } + } + network.nd=id; network.nc=ic; + j=0; + + do{ //calculate average and prob for each node + for(i=0;i<=n;i++) //n is nimber of column in the data matrix + { + if(network.node[i].type==1) //if continuous + { + fscanf(pFile,"%s",line); + tc=run_strtof(line); + //network.node[i].prob[1]+=tc; //summation for average + network.node[i].cdata[j]=tc; // data values in data column vector + } + else //discrete + { + fscanf(pFile,"%s",line); + network.node[i].ddata[j]=(char *)malloc(sizeof(char)*strlen(line)); + if(j==0) //first iteration of do while + { + network.node[i].lindex=0; + network.node[i].label[0]=(char *)malloc(sizeof(char)*strlen(line)); + strcpy(network.node[i].label[0],line); //new label assigned to label list + network.node[i].post_alpha[0]=1; + + } + else //store different label for discrete node also count number of ocarance of each label for calculation of post alpha score + { + imatch=0; + for(k=0;k<=network.node[i].lindex;k++) + { + if(strcmp(network.node[i].label[k],line)==0) + { + imatch=1; + network.node[i].post_alpha[k]+=1; + break; + } + } + if(imatch==0) + { + network.node[i].lindex++; + k=network.node[i].lindex; + network.node[i].label[k]=(char *)malloc(sizeof(char)*strlen(line)); + strcpy(network.node[i].label[k],line); //new label assigned to label list + network.node[i].post_alpha[k]+=1; + } + } + strcpy(network.node[i].ddata[j],line); + } + + } + j++; + } while (fgets (line,250,pFile)!=NULL); //end of do while to read the input file + + //////////////normalize continuous data///////////// + for(i=0;i<=n;i++) + { + if(network.node[i].type==1) + { + + tcm=mean(network.node[i].cdata,row); + tc=stdev(network.node[i].cdata,row,tcm); + for(j=0;j<row;j++) + { + network.node[i].cdata[j]=standarizedata(network.node[i].cdata[j],tcm,tc); //dataval.data mean, standard deviation + } + } + } + + //////////////////adjust average and prob for normalized data + for(i=0;i<=n;i++) + { + if(network.node[i].type==1) //claculate average and prob for continuous nodes + { + network.node[i].prob[1]=mean(network.node[i].cdata,row); + + for(j=0;j<row;j++) + { + network.node[i].prob[0]+=(network.node[i].cdata[j]-network.node[i].prob[1])*(network.node[i].cdata[j]-network.node[i].prob[1]); + } + network.node[i].prob[0]*=1.0/(double)(row-1); + network.node[i].prob[0]*=(double)(row-1.0)/(double)row; + + } + } + + fclose (pFile); + }//end of else part for open input file + + pFile=fopen (argv[1],"r"); + + for(i=0;i<network.num;i++) + fscanf(pFile,"%s",name[i]); + + fclose(pFile); + + learn_node(&network); //call function to learn network and get local score whell 0 parent for all nodes + + + + //printf("%.2lf\n",learn_parent(&network,6,176)); + // printf("%.2lf\n",learn_parent(&network,6,47)); + // printf("%.2lf\n",learn_parent(&network,6,18)); + //printf("%.2lf\n",learn_parent(&network,7,80)); + ////////printing trylist to a file//////////////////////////////// + + //fprintf(outfile,"Number of nodes\n%d\n",network.num); + //fprintf(outfile,"Possible sets of parents per node\n"); + j=(int)pow(2.0,network.nd-1); //possible number of parents for discrete node + k=(int)pow(2.0,network.num-1); //possible number of parents for continuous nodes + + + for(i=0;i<network.num;i++) + { + network.node[i].score=(double *)malloc(sizeof(double)*(k+1)); + network.node[i].score_i=0; + + + network.node[i].wlist.n=0; + network.node[i].blist.n=0; + network.node[i].blist.list=(int *)malloc(sizeof(int)*network.num); + network.node[i].wlist.list=(int *)malloc(sizeof(int)*network.num); + // if(network.node[i].type>1) //discrete + // fprintf(outfile,"%d\t",j); + //else //continuous + // fprintf(outfile,"%d\t",k); + } + + banf=fopen (argv[2],"r"); //ban list file name + whitef=fopen (argv[3],"r"); //white list file name + + //read banlist +// fgets (line,2500,banf); //header + + +strcpy(nf,""); +strcpy(nt,""); + +while (fgets (line,2500,banf)!=NULL){ + + fscanf(banf,"%s",nf); //node from + fscanf(banf,"%s",nt); //node to + + for(i=0;i<network.num;i++) + { + if(strcmp(name[i],nt)==0) + { + j=network.node[i].blist.n; + for(l=0;l<network.num;l++) + { + if(strcmp(name[l],nf)==0) + { + ban_flag=0; + for(k=0;k<network.node[i].blist.n;k++) + { + if(network.node[i].blist.list[k]==l) + { + ban_flag=1; + break; + } + } + if(ban_flag==0) + { + network.node[i].blist.list[j]=l; + j++; + network.node[i].blist.n=j; + } + + } + } + } + } + +} + + + +//read whitelist +//fgets (line,2500,whitef); //header + +strcpy(nf,""); +strcpy(nt,""); + +while (fgets (line,2500,whitef)!=NULL){ + +fscanf(whitef,"%s",nf); //node from +fscanf(whitef,"%s",nt); //node to + +for(i=0;i<network.num;i++) + { + if(strcmp(name[i],nt)==0) + { + j=network.node[i].wlist.n; + for(l=0;l<network.num;l++) + { + if(strcmp(name[l],nf)==0) + { + white_flag=0; + for(k=0;k<network.node[i].wlist.n;k++) + { + if(network.node[i].wlist.list[k]==l) + { + white_flag=1; + break; + } + } + if(white_flag==0) + { + + network.node[i].wlist.list[j]=l; + j++; + network.node[i].wlist.n=j; + } + + } + } + } + } + +} + + //fprintf(outfile,"\nChild Parents Score\n"); + + m=(int)pow(2.0,network.num)-1; //number of combinations of parents + min_score=1.0; + for(i=0;i<network.num;i++) + { + + if(network.node[i].wlist.n==0) + { + network.node[i].score[0]=network.node[i].score_zeroparent; + if(min_score>network.node[i].score_zeroparent) + min_score=network.node[i].score_zeroparent; + + //fprintf(outfile,"%d\t0\t%lf\n",i+1,network.node[i].score_zeroparent); //printing score for zero parent + } + else + { + network.node[i].score[0]=1.0; + //fprintf(outfile,"%d\t0\t0.0\n",i+1); + + } + network.node[i].score_i=1; + + + for(j=0;j<m;j++) + { + + ban_flag=0; + + bit= (j+1 >> i) & 1; + + if(bit!=1) //check if the child itself is present in parent list + { + + bit_count=0; + //fprintf(outfile,"%d\t",i+1); + for(l=0;l<network.num;l++) + { + bit= (j+1 >> l) & 1; + if(bit==1) + { + //fprintf(outfile,"%d",l+1); + bit_count++; + + } + + } + + + if(network.node[i].type>1) //if discrete node checck if any continuous node is included in the parent list + { + flag=0; + ban_flag=0; + for(l=0;l<network.num;l++) + { + bit= (j+1 >> l) & 1; + + if(bit==1 && network.node[l].type==1) //check if continuous is parent of discrete + { + flag=1; + break; + } + if(bit==1) + { + for(k=0;k<network.node[i].blist.n;k++) + { + if(network.node[i].blist.list[k]==l) + { + ban_flag=1; + break; + } + } + + + if(ban_flag==1) + break; + } + + + + } + if(flag==0) //if no continuous parent for discrete node then learn the network with the list of parents + { + + if(bit_count>MAX_PARENT) + { + //fprintf(outfile,"\t0.0\n"); + + s_i=network.node[i].score_i; + network.node[i].score[s_i]=1.0; + network.node[i].score_i=s_i+1; + + + + } + else if(ban_flag==1) + { + //fprintf(outfile,"\t0.0\n"); + s_i=network.node[i].score_i; + network.node[i].score[s_i]=1.0; + network.node[i].score_i=s_i+1; + } + else // calculate score if white list check get satisfied + { + if(network.node[i].wlist.n==0) + { + score_val=learn_parent(&network,i,j); + + if(min_score>score_val) + min_score=score_val; + + + + //fprintf(outfile,"\t%lf\n",score_val); + s_i=network.node[i].score_i; + network.node[i].score[s_i]=score_val; + network.node[i].score_i=s_i+1; + } + else if(bit_count<network.node[i].wlist.n) //white list is greater than list of parents + { + //fprintf(outfile,"\tfor white\t0.0\n"); + s_i=network.node[i].score_i; + network.node[i].score[s_i]=1.0; + network.node[i].score_i=s_i+1; + } + else + { + white_flag=0; + white_flag_all=0; + + for(l=0;l<network.num;l++) + { + bit= (j+1 >> l) & 1; + + if(bit==1) + { + for(k=0;k<network.node[i].wlist.n;k++) + { + if(network.node[i].wlist.list[k]==l) + { + white_flag++; + //changed from white_flag=1; + break; + } + } + if(white_flag>=network.node[i].wlist.n) + { //!=1 + white_flag_all=1; + break; + + } + + } + } + if(white_flag_all==1) + { + score_val=learn_parent(&network,i,j); + if(min_score>score_val) + min_score=score_val; + + //fprintf(outfile,"\t%lf\n",score_val); + s_i=network.node[i].score_i; + network.node[i].score[s_i]=score_val; + network.node[i].score_i=s_i+1; + } + else + { + //fprintf(outfile,"\t0.0\n"); + s_i=network.node[i].score_i; + network.node[i].score[s_i]=1.0; + network.node[i].score_i=s_i+1; + } + + } + } + + } + else //fill up data when continuous is parent of discrete + { + //fprintf(outfile,"\t0.0\n"); + s_i=network.node[i].score_i; + network.node[i].score[s_i]=1.0; + network.node[i].score_i=s_i+1; + + } + } + else //if continuous node then learn the network with the parent set + { + + ban_flag=0; + + + for(l=0;l<network.num;l++) + { + bit= (j+1 >> l) & 1; + if(bit==1) + { + + for(k=0;k<network.node[i].blist.n;k++) + { + if(network.node[i].blist.list[k]==l) + { + ban_flag=1; + break; + } + } + if(ban_flag==1) + break; + } + + + } + if(bit_count>MAX_PARENT) + { + //fprintf(outfile,"\t0.0\n"); + s_i=network.node[i].score_i; + network.node[i].score[s_i]=1.0; + network.node[i].score_i=s_i+1; + } + else if(ban_flag==1) + { + //fprintf(outfile,"\t0.0\n"); + s_i=network.node[i].score_i; + network.node[i].score[s_i]=1.0; + network.node[i].score_i=s_i+1; + } + else + { + if(network.node[i].wlist.n==0) + { + score_val=learn_parent(&network,i,j); + if(min_score>score_val) + min_score=score_val; + + //fprintf(outfile,"\t%lf\n",score_val); + s_i=network.node[i].score_i; + network.node[i].score[s_i]=score_val; + network.node[i].score_i=s_i+1; + } + else if(bit_count<network.node[i].wlist.n) //white list is greater than list of parents + { + //fprintf(outfile,"\t0.0\n"); + s_i=network.node[i].score_i; + network.node[i].score[s_i]=1.0; + network.node[i].score_i=s_i+1; + } + else + { + white_flag=0; + white_flag_all=0; + + for(l=0;l<network.num;l++) + { + bit= (j+1 >> l) & 1; + + if(bit==1) + { + for(k=0;k<network.node[i].wlist.n;k++) + { + if(network.node[i].wlist.list[k]==l) + { + white_flag++; + break; + } + } + if(white_flag>=network.node[i].wlist.n) + { //!=1 + white_flag_all=1; + break; + + } + + } + } + if(white_flag_all==1) + { + score_val=learn_parent(&network,i,j); + if(min_score>score_val) + min_score=score_val; + + //fprintf(outfile,"\t%lf\n",score_val); + s_i=network.node[i].score_i; + network.node[i].score[s_i]=score_val; + network.node[i].score_i=s_i+1; + } + else + { + // fprintf(outfile,"\t0.0\n"); + s_i=network.node[i].score_i; + network.node[i].score[s_i]=1.0; + network.node[i].score_i=s_i+1; + } + + } + } + + + } + + + } + } + } + + +//print scores to separate files + +min_score*=10; +inputFiles= (FILE **) malloc(network.num * sizeof(FILE*)); + +printf("%s\n",argv[5]); + +for(i=0;i<network.num;i++) +{ + strcpy(line,argv[5]); + sprintf(line+10,"/%d",i); //sprintf(line+52,"%d",i); + + // strcat(line,".txt"); + + inputFiles[i] = fopen(line, "wb"); + + //print content + for(j=0;j<network.node[i].score_i;j++) + { + if(network.node[i].score[j]<1.0) + fwrite(&network.node[i].score[j],sizeof(double),1,inputFiles[i]); +// fprintf(inputFiles[i],"%lf\n",network.node[i].score[j]); + else + fwrite(&min_score,sizeof(double),1,inputFiles[i]); +//fprintf(inputFiles[i],"%lf\n",min_score); + } + + fclose(inputFiles[i]); +} + + + + + + +fclose(banf); +fclose(whitef); + + +// printf("Time elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC); + +} + + +/* learn the network for zero parent */ + +void learn_node(struct nd *net) +{ + int i,j,k,nk,l,m,i1,j1; + long num,temp; + double n,N,s2; + struct nd network=*net; //create a local copy + double *alpha,nu,rho,mu,phi,tau; + double post_nu,post_rho,post_mu,post_phi,post_tau; + double loglik,post_loglik; + + num=network.num; + + nk=1; + + for(i=0;i<num;i++) + { + if(network.node[i].type>1) + nk*=network.node[i].type; + } + + + + + for(i=0;i<num;i++) + { + if(network.node[i].type>1) + { + + alpha=(double *)malloc(sizeof(double)*network.node[i].type); + k=1; + + for(j=0;j<num;j++) + { + if(i!=j && network.node[j].type>1) + k*=network.node[j].type; + + } + N=0.0; n=0.0; + for(j=0;j<network.node[i].type;j++) + { + alpha[j]=2.0*k; + network.node[i].post_alpha[j]+=alpha[j]; + ////////update loglik from udisclik + n+=alpha[j]; + N+=network.node[i].post_alpha[j]; + } + + + loglik=0.0; + post_loglik=0.0; + + for(j=0;j<network.node[i].type;j++) + post_loglik+=lgammafn(network.node[i].post_alpha[j])-lgammafn(alpha[j]); + + post_loglik+=lgammafn(n); + post_loglik-=lgammafn(N); + + + } + else + { + nu=nk*2.0; + tau=nu; + rho=nk*2.0; + mu=network.node[i].prob[1]; + phi=network.node[i].prob[0]*nk; + loglik=0.0; + post_mu=mu; + post_tau=tau; + post_rho=rho; + post_phi=phi; + post_loglik=loglik; + + postc0(&post_mu,&post_tau,&post_rho,&post_phi,&post_loglik,network.node[i].cdata,&network.row); + + + } + network.node[i].score_zeroparent=post_loglik; // store the zero parent score + + } + + if(nk!=1) + free(alpha); + + *net=network; //coppied back to original +} + + +/* +learn node for a list of parents +*/ + +double learn_parent(struct nd *net,int child,int p) +{ + int i,j,jj,k,nk,l,m,i1,j1,k1,k2,cn_send; + long num,temp; + double n,N,s2; + struct nd network=*net; //create a local copy + int cn,ds,bit,*discrete,*continuous; + double alpha,*postalpha,alpha_parent,*postalpha_parent,nu,rho,mu,phi,tau; + double post_nu,post_rho,*post_mu,post_phi,post_tau; + double loglik,post_loglik=0,post_loglik_parent=0; + double *mu_mat,**tau_mat,tau11,tau12,tau22,div,phiinv; + char **label; + int *a,*l_list,flag,y_count,z_count; + double *y,mu_p,phi_p,*tau_send,*z_send,loglik_sum; + + + + num=network.num; + + continuous=(int *)malloc(sizeof(int)*network.nc); + discrete=(int *)malloc(sizeof(int)*network.nd); + + nk=1; + for(i=0;i<num;i++) + { + if(network.node[i].type>1) + nk*=network.node[i].type; + } + +//list of discrete snd continuous parents + cn=0; + ds=0; + alpha_parent=0; + alpha=2; + + k1=network.node[child].type; + k2=1; + + //printf("%d\n",child); + for(j=0;j<network.num;j++) //list the continuous parents and discrete parents + { + bit= (p+1 >> j) & 1; + + + if(bit==1 && j!=child) + { + + if(network.node[j].type>1) + { + discrete[ds]=j; + //printf("discrete parent=%d\n",j); + k1*=network.node[j].type; + k2*=network.node[j].type; + ds++; + } + else + { + continuous[cn]=j; + //printf("continuous parent=%d\n",j); + cn++; + } + } + else if(j!=child) + { + alpha*=network.node[j].type; + } + } + + if(network.nd<=2) + alpha=2; + + i=child; //i is the node with changes in its parrent nodes + + if(network.node[i].type>1 && cn==0) // no continuous node //prepare alpha and post alpha and then calculate local score + { + + k=k1; + postalpha=(double *)calloc(network.node[i].type,sizeof(double)); + + label=(char **)malloc(sizeof(char *)*(ds+1)); //label array + a=(int *)calloc((ds+1),sizeof(int)); + l_list=(int *)malloc(sizeof(int)*(ds+1)); + + for(j=0;j<=ds;j++) + label[j]=(char *)malloc(sizeof(char)*100); + + for(j=0;j<ds;j++) + { + l=discrete[j]; + l_list[j]=network.node[l].type; + + } + l_list[ds]=network.node[i].type; + + for(j=0;j<k;j++) + { + + get_next_seq(a,l_list,ds+1); // prepare a new combination of labels for parent + for(i1=0;i1<ds;i1++) + { + j1=a[i1]; + l=discrete[i1]; + strcpy(label[i1],network.node[l].label[j1]); + } + + j1=a[i1]; + strcpy(label[i1],network.node[child].label[j1]); + for(i1=0;i1<network.node[i].type;i1++) + postalpha[i1]=0; + for(i1=0;i1<network.row;i1++) // search through parents data and count how many number of matches for current label list. prepare post alpha from count+alpha + { + flag=0; + for(m=0;m<ds;m++) + { + j1=a[m]; + l=discrete[m]; + if(strcmp(label[m],network.node[l].ddata[i1])!=0) + { + flag=1; + } + } + j1=a[m]; + if(strcmp(label[m],network.node[child].ddata[i1])!=0) + { + flag=1; + } + if(flag==0) + { + j1=a[ds]; + postalpha[j1]+=1; + + } + } + for(m=0;m<network.node[i].type;m++) + { + if(postalpha[m]!=0) + { + postalpha[m]+=alpha; + post_loglik+=lgammafn(postalpha[m])-lgammafn(alpha); + } + } + } + + + k=k2; // repeat the last step of preparing alpha and post alpha for only parent set. this time exclude the child node label from computation. prepare alpha_parent and post_alphaparent + l=discrete[ds-1]; + postalpha_parent=(double *)calloc(network.node[l].type,sizeof(double)); + + for(j=0;j<ds;j++) + { + l=discrete[j]; + l_list[j]=network.node[l].type; + a[j]=0; + + } + for(j=0;j<k1/k2;j++) + alpha_parent+=alpha; + for(j=0;j<k;j++) + { + get_next_seq(a,l_list,ds); + + for(i1=0;i1<ds;i1++) + { + j1=a[i1]; + l=discrete[i1]; + strcpy(label[i1],network.node[l].label[j1]); + + } + + l=discrete[ds-1]; + for(i1=0;i1<network.node[l].type;i1++) + postalpha_parent[i1]=0; + + for(i1=0;i1<network.row;i1++) + { + flag=0; + for(m=0;m<ds;m++) + { + j1=a[m]; + l=discrete[m]; + if(strcmp(label[m],network.node[l].ddata[i1])!=0) + { + flag=1; + } + } + + if(flag==0) + { + j1=a[ds-1]; + postalpha_parent[j1]+=1; + + } + } + l=discrete[ds-1]; + for(m=0;m<network.node[l].type;m++) + { + if(postalpha_parent[m]!=0) + { + postalpha_parent[m]+=alpha_parent; + post_loglik_parent+=lgammafn(postalpha_parent[m])-lgammafn(alpha_parent); + } + } + } + + + loglik=post_loglik-post_loglik_parent; + free(postalpha); + free(postalpha_parent); + for(j=0;j<=ds;j++) + free(label[j]); + free(label); + free(a); + free(l_list); + + } + else if(network.node[i].type==1 && ds==0) // no discrete node //prepare continuous parameter set + { + z_count=0; + tau_mat=(double **)malloc(sizeof(double *)*(cn+1)); + + z_send=(double *)malloc(sizeof(double)*(cn+1)*network.row); + tau_send=(double *)malloc(sizeof(double)*(cn+1)*(cn+1)); + + mu_mat=(double *)malloc(sizeof(double)*(cn+1)); + nu=nk*2.0; rho=nk*2.0+cn; + mu_mat[0]=network.node[i].prob[1]; + for(j=0;j<cn;j++) + { + mu_mat[j+1]=0; + tau_mat[j]=(double *)malloc(sizeof(double)*(cn+1)); + + } + + tau_mat[j]=(double *)malloc(sizeof(double)*(cn+1)); + + for(j=0;j<network.row;j++) + { + z_send[z_count]=1; + z_count++; + for(j1=0;j1<cn;j1++) + { + l=continuous[j1]; + z_send[z_count]=network.node[l].cdata[j]; + z_count++; + } + + } + + phi=network.node[i].prob[0]*nk; + + //preparing tau + + //**tau_mat + for(j=0;j<=cn;j++) + for(m=0;m<=cn;m++) + { + tau_mat[j][m]=1.0; + } + + for(j=0;j<cn;j++) + { + l=continuous[j]; + mu_p=network.node[l].prob[1]; + phi_p=network.node[l].prob[0]*nk; + phiinv=1/phi_p; + tau11 = 1/nu + mu_p*phiinv*mu_p; + tau22 = phiinv; + tau12 = -mu_p*phiinv; + div=tau11*tau22-tau12*tau12; + if(j==0) + tau_mat[0][0]=tau22/div; + tau_mat[0][j+1]=-tau12/div; + tau_mat[j+1][0]=tau_mat[0][j+1]; + tau_mat[j+1][j+1]=tau11/div; + } + + for(j=0;j<=cn;j++) + for(m=0;m<=cn;m++) + { + if(tau_mat[j][m]==1.0) + { + tau_mat[j][m]=tau_mat[0][j]*(tau_mat[0][m]/tau_mat[0][0]); + + } + j1=j*(cn+1)+m; + tau_send[j1]=tau_mat[j][m]; + + } + + loglik=0.0; + + cn_send=cn+1; + + postc(mu_mat,tau_send,&rho,&phi,&loglik,network.node[i].cdata,z_send,&network.row,&cn_send); + + + free(z_send); + free(tau_send); + for(j=0;j<=cn;j++) + free(tau_mat[j]); + free(tau_mat); + free(mu_mat); + + } + else //hybrid cases + { + if(cn==0) //if parents are all discrete and child is continuous + { + + ///////////////////prepare y array and run c code//////////////////////////////// + loglik_sum=0; + k=k2; + y=(double *)calloc(network.row,sizeof(double)); + label=(char **)malloc(sizeof(char *)*(ds)); //label array + a=(int *)calloc((ds+1),sizeof(int)); + l_list=(int *)malloc(sizeof(int)*(ds)); + + for(j=0;j<ds;j++) + label[j]=(char *)malloc(sizeof(char)*100); + + for(j=0;j<ds;j++) + { + l=discrete[j]; + l_list[j]=network.node[l].type; + + } + + + for(j=0;j<k;j++) + { + + nu=(nk*2.0)/k2; + tau=nu; + rho=nu; + mu=network.node[i].prob[1]; + phi=network.node[i].prob[0]*nu/2; + loglik=0.0; + + get_next_seq(a,l_list,ds); // prepare a new combination of labels for parent + for(i1=0;i1<ds;i1++) + { + j1=a[i1]; + l=discrete[i1]; + strcpy(label[i1],network.node[l].label[j1]); + } + y_count=0; + for(i1=0;i1<network.row;i1++) // search through parents data and count how many number of matches for current label list. prepare post alpha from count+alpha + { + flag=0; + for(m=0;m<ds;m++) + { + + l=discrete[m]; + if(strcmp(label[m],network.node[l].ddata[i1])!=0) + { + flag=1; + } + } + + if(flag==0) + { + y[y_count]=network.node[i].cdata[i1]; + y_count++; + + } + } + + postc0(&mu,&tau,&rho,&phi,&loglik,y,&y_count); + loglik_sum+=loglik; + + + } + + loglik=loglik_sum; + + free(y); + for(j=0;j<ds;j++) + free(label[j]); + free(label); + free(a); + free(l_list); + + + }//end of continuous node discrete parent check + else //all the other hybrid cases when parents are both continuous and discrete + { + + + tau_mat=(double **)malloc(sizeof(double *)*(cn+1)); + z_send=(double *)malloc(sizeof(double)*(cn+1)*network.row); + tau_send=(double *)malloc(sizeof(double)*(cn+1)*(cn+1)); + mu_mat=(double *)malloc(sizeof(double)*(cn+1)); + + for(j=0;j<cn;j++) + { + tau_mat[j]=(double *)malloc(sizeof(double)*(cn+1)); + } + + + tau_mat[j]=(double *)malloc(sizeof(double)*(cn+1)); + + ///////////////////prepare parameters and run c code//////////////////////////////// + loglik_sum=0.0; + k=k2; + y=(double *)calloc(network.row,sizeof(double)); + label=(char **)malloc(sizeof(char *)*(ds)); //label array + a=(int *)calloc((ds+1),sizeof(int)); + l_list=(int *)malloc(sizeof(int)*(ds)); + + for(j=0;j<ds;j++) + label[j]=(char *)malloc(sizeof(char)*100); + + for(j=0;j<ds;j++) + { + l=discrete[j]; + l_list[j]=network.node[l].type; + + } + + + for(jj=0;jj<k;jj++) + { + + nu=(nk*2.0)/k2; + rho=nu+cn; + mu_mat[0]=network.node[i].prob[1]; + + for(j1=0;j1<cn;j1++) + mu_mat[j1+1]=0; + + phi=network.node[i].prob[0]*nu/2; + //preparing tau + //**tau_mat + for(j=0;j<=cn;j++) + for(m=0;m<=cn;m++) + { + tau_mat[j][m]=1.0; + } + + for(j=0;j<cn;j++) + { + l=continuous[j]; + mu_p=network.node[l].prob[1]; + phi_p=network.node[l].prob[0]*(nk/k2); + + phiinv=1/phi_p; + tau11 = 1/nu + mu_p*phiinv*mu_p; + tau22 = phiinv; + tau12 = -mu_p*phiinv; + div=tau11*tau22-tau12*tau12; + if(j==0) + tau_mat[0][0]=tau22/div; + + tau_mat[0][j+1]=-tau12/div; + tau_mat[j+1][0]=tau_mat[0][j+1]; + tau_mat[j+1][j+1]=tau11/div; + + + } + for(j=0;j<=cn;j++) + for(m=0;m<=cn;m++) + { + if(tau_mat[j][m]==1.0) + { + tau_mat[j][m]=tau_mat[0][j]*(tau_mat[0][m]/tau_mat[0][0]); + } + + j1=j*(cn+1)+m; + tau_send[j1]=tau_mat[j][m]; + + //printf("%.2lf\t",tau_mat[j][m]); + } + //printf("\n"); + loglik=0.0; + + + cn_send=cn+1; + + + get_next_seq(a,l_list,ds); // prepare a new combination of labels for parent + for(i1=0;i1<ds;i1++) + { + j1=a[i1]; + l=discrete[i1]; + strcpy(label[i1],network.node[l].label[j1]); + } + y_count=0; + z_count=0; + for(i1=0;i1<network.row;i1++) // search through parents data and count how many number of matches for current label list. prepare post alpha from count+alpha + { + flag=0; + for(m=0;m<ds;m++) + { + + l=discrete[m]; + if(strcmp(label[m],network.node[l].ddata[i1])!=0) + { + flag=1; + } + } + + if(flag==0) + { + y[y_count]=network.node[i].cdata[i1]; + + //prepare z send matrix + z_send[z_count]=1; + z_count++; + for(j1=0;j1<cn;j1++) + { + l=continuous[j1]; + z_send[z_count]=network.node[l].cdata[i1]; + z_count++; + } + y_count++; + } + } + + postc(mu_mat,tau_send,&rho,&phi,&loglik,y,z_send,&y_count,&cn_send); + + loglik_sum+=loglik; + + + } + + loglik=loglik_sum; + + free(y); + //free(z); + free(z_send); + free(tau_send); + for(j=0;j<=cn;j++) + free(tau_mat[j]); + free(tau_mat); + free(mu_mat); + for(j=0;j<ds;j++) + free(label[j]); + free(label); + free(a); + free(l_list); + + }//end of else part for hybrid when parents are mixed + + }//end of else part for hybrid + + free(continuous); + free(discrete); + return loglik; + + +} + +void get_next_seq(int *a,int *l,int size) +{ + int j,idx; + idx=size-1; + a[idx]+=1; + if(a[idx]==l[idx]) + { + a[idx]=0; + for(j=idx-1;j>=0;j--){ + a[j]+=1; + if(a[j]==l[j]) + a[j]=0; + else + break; + } + } +} + + |
