From c4f926438dcb8abe805e910399940f79ff643c4b Mon Sep 17 00:00:00 2001
From: ziejd2
Date: Thu, 14 Sep 2017 15:39:41 -0500
Subject: Add files via upload

---
 sourcecodes/localscore/network_score.c | 1441 ++++++++++++++++++++++++++++++++
 1 file changed, 1441 insertions(+)
 create mode 100644 sourcecodes/localscore/network_score.c

(limited to 'sourcecodes/localscore/network_score.c')

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;
+      }
+  }
+}
+
+
-- 
cgit 1.4.1