/* EMBLEM and SLIPCASE Copyright (c) 2011, Fabrizio Riguzzi and Elena Bellodi This package uses the library cudd, see http://vlsi.colorado.edu/~fabio/CUDD/ for the relative license. */ #include <math.h> #include <stdlib.h> #include <stdio.h> #include <string.h> #include "cuddInt.h" #include "YapInterface.h" #define LOGZERO log(0.000001) #define CACHE_SLOTS 1 #define UNIQUE_SLOTS 1 typedef struct { int nVal,nRule; int firstBoolVar; } variable; typedef struct { DdNode *key; double value; } rowel; typedef struct { int cnt; rowel *row; } tablerow; tablerow * table; static variable * vars; static variable ** vars_ex; static int * bVar2mVar; static int ** bVar2mVar_ex; static double * sigma; static double ***eta; static double ***eta_temp; static double **arrayprob; static int *rules; static DdManager *mgr; static DdManager **mgr_ex; static int *nVars; static int *nVars_ex; static int nRules; double * probs; double * nodes_probs_ex; double ** probs_ex; static int * boolVars; static int * boolVars_ex; tablerow * nodesB; tablerow * nodesF; int ex,cycle; DdNode *** nodesToVisit; int * NnodesToVisit; double * example_prob; static int ret_prob(void); double Prob(DdNode *node,int comp_par); static int end_bdd(void); static int init_test(void); static int add_var(void); static int init(void); static int end(void); static int EM(void); static int Q(void); double ProbPath(DdNode *node, int comp_par); static int rec_deref(void); int indexMvar(DdNode *node); void Forward(DdNode *node); void GetForward(DdNode *node, double ForwProbPath); void UpdateForward(DdNode * node); double GetOutsideExpe(DdNode *root,double ex_prob); void Maximization(void); static double Expectation(DdNode **nodes_ex, int lenNodes); void init_my_predicates(void); FILE *open_file(char *filename, const char *mode); tablerow* init_table(int varcnt); double * get_value(tablerow *tab, DdNode *node); void add_or_replace_node(tablerow *tab, DdNode *node, double value); void add_node(tablerow *tab, DdNode *node, double value); void destroy_table(tablerow *tab,int varcnt); static int init(void) { int j,i; YAP_Term arg1,arg2,list; ex=0; cycle=0; arg1=YAP_ARG1; arg2=YAP_ARG2; nRules=YAP_IntOfTerm(arg1); vars_ex=NULL; nVars_ex=NULL; eta= (double ***) malloc(nRules * sizeof(double **)); eta_temp= (double ***) malloc(nRules * sizeof(double **)); rules= (int *) malloc(nRules * sizeof(int)); arrayprob=(double **) malloc(nRules * sizeof(double *)); probs_ex=NULL; bVar2mVar_ex=NULL; boolVars_ex=NULL; mgr_ex=NULL; nodes_probs_ex=NULL; list=arg2; for (j=0;j<nRules;j++) { rules[j]=YAP_IntOfTerm(YAP_HeadOfTerm(list)); list=YAP_TailOfTerm(list); eta[j]= (double **) malloc((rules[j]-1)*sizeof(double *)); eta_temp[j]= (double **) malloc((rules[j]-1)*sizeof(double *)); arrayprob[j]= (double *) malloc((rules[j]-1)*sizeof(double)); for (i=0;i<rules[j]-1;i++) { eta[j][i]=(double *) malloc(2*sizeof(double)); eta_temp[j][i]=(double *) malloc(2*sizeof(double)); } } return 1; } static int init_bdd(void) { mgr=Cudd_Init(0,0,UNIQUE_SLOTS,CACHE_SLOTS,5120); Cudd_AutodynEnable(mgr, CUDD_REORDER_GROUP_SIFT); Cudd_SetMaxCacheHard(mgr, 0); Cudd_SetLooseUpTo(mgr, 0); Cudd_SetMinHit(mgr, 15); mgr_ex=(DdManager **) realloc(mgr_ex, (ex+1)* sizeof(DdManager *)); mgr_ex[ex]=mgr; bVar2mVar_ex=(int **) realloc(bVar2mVar_ex, (ex+1)* sizeof(int *)); bVar2mVar_ex[ex]=NULL; bVar2mVar=bVar2mVar_ex[ex]; vars_ex=(variable **) realloc(vars_ex, (ex+1)* sizeof(variable *)); vars_ex[ex]=NULL; vars=vars_ex[ex]; nVars_ex=(int *) realloc(nVars_ex, (ex+1)* sizeof(int )); nVars=nVars_ex+ex; *nVars=0; probs_ex=(double **) realloc(probs_ex, (ex+1)* sizeof(double *)); probs_ex[ex]=NULL; probs=probs_ex[ex]; boolVars_ex=(int *) realloc(boolVars_ex, (ex+1)* sizeof(int )); boolVars=boolVars_ex+ex; *boolVars=0; return 1; } static int end_bdd(void) { bVar2mVar_ex[ex]=bVar2mVar; probs_ex[ex]=probs; vars_ex[ex]=vars; ex=ex+1; return 1; } static int init_test(void) { YAP_Term arg1; arg1=YAP_ARG1; nRules=YAP_IntOfTerm(arg1); mgr=Cudd_Init(0,0,UNIQUE_SLOTS,CACHE_SLOTS,0); Cudd_AutodynEnable(mgr, CUDD_REORDER_GROUP_SIFT); Cudd_SetMaxCacheHard(mgr, 1024*1024*1024); Cudd_SetLooseUpTo(mgr, 1024*1024*512); rules= (int *) malloc(nRules * sizeof(int)); bVar2mVar=NULL; probs=NULL; vars=NULL; nVars=(int *) malloc(sizeof(int )); *nVars=0; boolVars=(int *) malloc(sizeof(int )); *boolVars=0; return 1; } static int end_test(void) { free(bVar2mVar); free(vars); free(nVars); free(boolVars); Cudd_Quit(mgr); free(probs); free(rules); return 1; } static double Expectation(DdNode **nodes_ex,int lenNodes) { int i; double rootProb,CLL=0; for(i=0;i<lenNodes;i++) { if (!Cudd_IsConstant(nodes_ex[i])) { mgr=mgr_ex[i]; probs=probs_ex[i]; boolVars=boolVars_ex+i; nodesB=init_table(*boolVars); nodesF=init_table(*boolVars); bVar2mVar=bVar2mVar_ex[i]; vars=vars_ex[i]; Forward(nodes_ex[i]); rootProb=GetOutsideExpe(nodes_ex[i],example_prob[i]); if (rootProb<=0.0) CLL = CLL + LOGZERO*example_prob[i]; else CLL = CLL + log(rootProb)*example_prob[i]; nodes_probs_ex[i]=rootProb; destroy_table(nodesB,*boolVars); destroy_table(nodesF,*boolVars); } else if (nodes_ex[i]==Cudd_ReadLogicZero(mgr_ex[i])) CLL=CLL+LOGZERO*example_prob[i]; } return CLL; } static int end(void) { int r,i; for (i=0;i<ex;i++) { Cudd_Quit(mgr_ex[i]); free(bVar2mVar_ex[i]); free(vars_ex[i]); free(probs_ex[i]); fflush(stdout); } free(mgr_ex); free(bVar2mVar_ex); free(vars_ex); free(probs_ex); free(nVars_ex); free(boolVars_ex); free(nodes_probs_ex); for (r=0;r<nRules;r++) { for (i=0;i<rules[r]-1;i++) { free(eta[r][i]); free(eta_temp[r][i]); } free(eta[r]); free(eta_temp[r]); } free(eta); free(eta_temp); for (r=0;r<nRules;r++) { free(arrayprob[r]); } free(arrayprob); free(rules); return 1; } static int ret_prob(void) { YAP_Term arg1,arg2,out; DdNode * node; arg1=YAP_ARG1; arg2=YAP_ARG2; node=(DdNode *)YAP_IntOfTerm(arg1); if (!Cudd_IsConstant(node)) { table=init_table(*boolVars); out=YAP_MkFloatTerm(Prob(node,0)); destroy_table(table,*boolVars); } else { if (node==Cudd_ReadOne(mgr)) out=YAP_MkFloatTerm(1.0); else out=YAP_MkFloatTerm(0.0); } return(YAP_Unify(out,arg2)); } double Prob(DdNode *node,int comp_par) /* compute the probability of the expression rooted at node. table is used to store nodeB for which the probability has alread been computed so that it is not recomputed */ { int index,mVarIndex,comp,pos; variable v; double res; double p,pt,pf,BChild0,BChild1; double * value_p; DdNode *nodekey,*T,*F; comp=Cudd_IsComplement(node); comp=(comp && !comp_par) ||(!comp && comp_par); if (Cudd_IsConstant(node)) { if (comp) return 0.0; else return 1.0; } else { nodekey=Cudd_Regular(node); /* if (comp) nodekey=Cudd_Complement(nodefw); else nodekey=nodefw;*/ value_p=get_value(table,nodekey); if (value_p!=NULL) return *value_p; else { index=Cudd_NodeReadIndex(node); //Returns the index of the node. The node pointer can be either regular or complemented. //The index field holds the name of the variable that labels the node. The index of a variable is a permanent attribute that reflects the order of creation. p=probs[index]; T = Cudd_T(node); F = Cudd_E(node); pf=Prob(F,comp); pt=Prob(T,comp); BChild0=pf*(1-p); BChild1=pt*p; mVarIndex=bVar2mVar[index]; v=vars[mVarIndex]; pos=index-v.firstBoolVar; res=BChild0+BChild1; add_node(table,nodekey,res); return res; } } } static int add_var(void) { YAP_Term arg1,arg2,arg3,arg4,out,probTerm,probTerm_temp; variable * v; int i; DdNode * node; double p,p0; arg1=YAP_ARG1; arg2=YAP_ARG2; arg3=YAP_ARG3; arg4=YAP_ARG4; *nVars=*nVars+1; vars=(variable *) realloc(vars,*nVars * sizeof(variable)); v=&vars[*nVars-1]; v->nVal=YAP_IntOfTerm(arg1); v->nRule=YAP_IntOfTerm(arg3); v->firstBoolVar=*boolVars; probs=(double *) realloc(probs,(((*boolVars+v->nVal-1)* sizeof(double)))); bVar2mVar=(int *) realloc(bVar2mVar,((*boolVars+v->nVal-1)* sizeof(int))); probTerm=arg2; p0=1; for (i=0;i<v->nVal-1;i++) { node=Cudd_bddIthVar(mgr,*boolVars+i); p=YAP_FloatOfTerm(YAP_HeadOfTerm(probTerm)); bVar2mVar[*boolVars+i]=*nVars-1; probs[*boolVars+i]=p/p0; probTerm_temp=YAP_TailOfTerm(probTerm); probTerm=probTerm_temp; p0=p0*(1-p/p0); } *boolVars=*boolVars+v->nVal-1; rules[v->nRule]= v->nVal; out=YAP_MkIntTerm((YAP_Int)* nVars-1); return YAP_Unify(out,arg4); } static int equality(void) { YAP_Term arg1,arg2,arg3,out; int varIndex; int value; int i; variable v; DdNode * node, * tmp,*var; arg1=YAP_ARG1; //var arg2=YAP_ARG2; //value arg3=YAP_ARG3; //node varIndex=YAP_IntOfTerm(arg1); value=YAP_IntOfTerm(arg2); v=vars[varIndex]; i=v.firstBoolVar; tmp=Cudd_ReadOne(mgr); Cudd_Ref(tmp); node=NULL; for (i=v.firstBoolVar;i<v.firstBoolVar+value;i++) { var=Cudd_bddIthVar(mgr,i); node=Cudd_bddAnd(mgr,tmp,Cudd_Not(var)); Cudd_Ref(node); Cudd_RecursiveDeref(mgr,tmp); tmp=node; } if (!(value==v.nVal-1)) { var=Cudd_bddIthVar(mgr,v.firstBoolVar+value); node=Cudd_bddAnd(mgr,tmp,var); Cudd_Ref(node); Cudd_RecursiveDeref(mgr,tmp); } out=YAP_MkIntTerm((YAP_Int) node); return(YAP_Unify(out,arg3)); } static int one(void) { YAP_Term arg,out; DdNode * node; arg=YAP_ARG1; node = Cudd_ReadOne(mgr); Cudd_Ref(node); out=YAP_MkIntTerm((YAP_Int) node); return(YAP_Unify(out,arg)); } static int zero(void) { YAP_Term arg,out; DdNode * node; arg=YAP_ARG1; node = Cudd_ReadLogicZero(mgr); Cudd_Ref(node); out=YAP_MkIntTerm((YAP_Int) node); return(YAP_Unify(out,arg)); } static int bdd_not(void) { YAP_Term arg1,arg2,out; DdNode * node; arg1=YAP_ARG1; arg2=YAP_ARG2; node = (DdNode *)YAP_IntOfTerm(arg1); node=Cudd_Not(node); out=YAP_MkIntTerm((YAP_Int) node); return(YAP_Unify(out,arg2)); } static int and(void) { YAP_Term arg1,arg2,arg3,out; DdNode * node1, *node2,*nodeout; arg1=YAP_ARG1; arg2=YAP_ARG2; arg3=YAP_ARG3; node1=(DdNode *)YAP_IntOfTerm(arg1); node2=(DdNode *)YAP_IntOfTerm(arg2); nodeout=Cudd_bddAnd(mgr,node1,node2); Cudd_Ref(nodeout); out=YAP_MkIntTerm((YAP_Int) nodeout); return(YAP_Unify(out,arg3)); } static int or(void) { YAP_Term arg1,arg2,arg3,out; DdNode * node1,*node2,*nodeout; arg1=YAP_ARG1; arg2=YAP_ARG2; arg3=YAP_ARG3; node1=(DdNode *)YAP_IntOfTerm(arg1); node2=(DdNode *)YAP_IntOfTerm(arg2); nodeout=Cudd_bddOr(mgr,node1,node2); Cudd_Ref(nodeout); out=YAP_MkIntTerm((YAP_Int) nodeout); return(YAP_Unify(out,arg3)); } static int garbage_collect(void) { YAP_Term arg1,arg2,out; YAP_Int nodes,clearCache; arg1=YAP_ARG1; arg2=YAP_ARG2; clearCache=YAP_IntOfTerm(arg1); nodes=(YAP_Int)cuddGarbageCollect(mgr,clearCache); out=YAP_MkIntTerm(nodes); return(YAP_Unify(out,arg2)); } static int bdd_to_add(void) { YAP_Term arg1,arg2,out; DdNode * node1,*node2; arg1=YAP_ARG1; arg2=YAP_ARG2; node1=(DdNode *)YAP_IntOfTerm(arg1); node2= Cudd_BddToAdd(mgr,node1); out=YAP_MkIntTerm((YAP_Int) node2); return(YAP_Unify(out,arg2)); } static int create_dot(void) { char * onames[]={"Out"}; char ** inames; DdNode * array[1]; YAP_Term arg1,arg2; int i,b,index; variable v; char numberVar[10],numberBit[10],filename[1000]; FILE * file; arg1=YAP_ARG1; arg2=YAP_ARG2; YAP_StringToBuffer(arg2,filename,1000); inames= (char **) malloc(sizeof(char *)*(*boolVars)); index=0; for (i=0;i<*nVars;i++) { v=vars[i]; for (b=0;b<v.nVal-1;b++) { inames[b+index]=(char *) malloc(sizeof(char)*20); strcpy(inames[b+index],"X"); sprintf(numberVar,"%d",i); strcat(inames[b+index],numberVar); strcat(inames[b+index],"_"); sprintf(numberBit,"%d",b); strcat(inames[b+index],numberBit); } index=index+v.nVal-1; } array[0]=(DdNode *)YAP_IntOfTerm(arg1); file = open_file(filename, "w"); Cudd_DumpDot(mgr,1,array,inames,onames,file); fclose(file); index=0; for (i=0;i<*nVars;i++) { v=vars[i]; for (b=0;b<v.nVal-1;b++) { free(inames[b+index]); } index=index+v.nVal-1; } free(inames); return 1; } static int rec_deref(void) { YAP_Term arg1; DdNode * node; arg1=YAP_ARG1; node=(DdNode *) YAP_IntOfTerm(arg1); Cudd_RecursiveDeref(mgr, node); return 1; } double ProbPath(DdNode *node,int comp_par) { int index,mVarIndex,comp,pos,position,boolVarIndex; variable v; double res; double value,p,pt,pf,BChild0,BChild1,e0,e1; double * value_p,** eta_rule; DdNode *nodekey,*T,*F; comp=Cudd_IsComplement(node); comp=(comp && !comp_par) ||(!comp && comp_par); if (Cudd_IsConstant(node)) { value=Cudd_V(node); if (comp) { return 0.0; } else { return 1.0; } } else { nodekey=Cudd_Regular(node); value_p=get_value(nodesB,nodekey); if (value_p!=NULL) { return *value_p; } else { index=Cudd_NodeReadIndex(node); p=probs[index]; T = Cudd_T(node); F = Cudd_E(node); pf=ProbPath(F,comp); pt=ProbPath(T,comp); BChild0=pf*(1-p); BChild1=pt*p; value_p=get_value(nodesF,nodekey); e0 = (*value_p)*BChild0; e1 = (*value_p)*BChild1; mVarIndex=bVar2mVar[index]; v=vars[mVarIndex]; pos=index-v.firstBoolVar; eta_rule=eta_temp[v.nRule]; eta_rule[pos][0]=eta_rule[pos][0]+e0; eta_rule[pos][1]=eta_rule[pos][1]+e1; res=BChild0+BChild1; add_node(nodesB,nodekey,res); position=Cudd_ReadPerm(mgr,index); position=position+1; boolVarIndex=Cudd_ReadInvPerm(mgr,position); //Returns the index of the variable currently in the i-th position of the order. if (position<*boolVars) { sigma[position]=sigma[position]+e0+e1; } if(!Cudd_IsConstant(T)) { index=Cudd_NodeReadIndex(T); position=Cudd_ReadPerm(mgr,index); sigma[position]=sigma[position]-e1; } if(!Cudd_IsConstant(F)) { index=Cudd_NodeReadIndex(F); position=Cudd_ReadPerm(mgr,index); sigma[position]=sigma[position]-e0; } return res; } } } void Forward(DdNode *root) { int i,j; if (*boolVars) { nodesToVisit= (DdNode ***)malloc(sizeof(DdNode **)* *boolVars); NnodesToVisit= (int *)malloc(sizeof(int)* *boolVars); nodesToVisit[0]=(DdNode **)malloc(sizeof(DdNode *)); nodesToVisit[0][0]=root; NnodesToVisit[0]=1; for(i=1;i<*boolVars;i++) { nodesToVisit[i]=NULL; NnodesToVisit[i]=0; } add_node(nodesF,Cudd_Regular(root),1); for(i=0;i<*boolVars;i++) { for(j=0;j<NnodesToVisit[i];j++) UpdateForward(nodesToVisit[i][j]); } for(i=0;i<*boolVars;i++) { free(nodesToVisit[i]); } free(nodesToVisit); free(NnodesToVisit); } else { add_node(nodesF,Cudd_Regular(root),1); } } void UpdateForward(DdNode *node) { int index,position,mVarIndex; DdNode *T,*E,*nodereg; variable v; double *value_p,*value_p_T,*value_p_F,p; if (Cudd_IsConstant(node)) { return; } else { index=Cudd_NodeReadIndex(node); mVarIndex=bVar2mVar[index]; v=vars[mVarIndex]; p=probs[index]; nodereg=Cudd_Regular(node); value_p=get_value(nodesF,nodereg); if (value_p== NULL) { printf("Error\n"); return; } else { T = Cudd_T(node); E = Cudd_E(node); if (!Cudd_IsConstant(T)) { value_p_T=get_value(nodesF,T); if (value_p_T!= NULL) { *value_p_T= *value_p_T+*value_p*p; } else { add_or_replace_node(nodesF,Cudd_Regular(T),*value_p*p); index=Cudd_NodeReadIndex(T); position=Cudd_ReadPerm(mgr,index); nodesToVisit[position]=(DdNode **)realloc(nodesToVisit[position], (NnodesToVisit[position]+1)* sizeof(DdNode *)); nodesToVisit[position][NnodesToVisit[position]]=T; NnodesToVisit[position]=NnodesToVisit[position]+1; } } if (!Cudd_IsConstant(E)) { value_p_F=get_value(nodesF,Cudd_Regular(E)); if (value_p_F!= NULL) { *value_p_F= *value_p_F+*value_p*(1-p); } else { add_or_replace_node(nodesF,Cudd_Regular(E),*value_p*(1-p)); index=Cudd_NodeReadIndex(E); position=Cudd_ReadPerm(mgr,index); nodesToVisit[position]=(DdNode **)realloc(nodesToVisit[position], (NnodesToVisit[position]+1)* sizeof(DdNode *)); nodesToVisit[position][NnodesToVisit[position]]=E; NnodesToVisit[position]=NnodesToVisit[position]+1; } } return; } } } int indexMvar(DdNode * node) { int index,mVarIndex; index=Cudd_NodeReadIndex(node); mVarIndex=bVar2mVar[index]; return mVarIndex; } double GetOutsideExpe(DdNode *root,double ex_prob) { int i,j,mVarIndex,bVarIndex; double **eta_rule; double theta,rootProb, T=0; sigma=(double *)malloc(*boolVars * sizeof(double)); for (j=0; j<*boolVars; j++) { sigma[j]=0; } for (j=0; j<nRules; j++) { for (i=0; i<rules[j]-1; i++) { eta_temp[j][i][0]=0; eta_temp[j][i][1]=0; } } rootProb=ProbPath(root,0); if (rootProb>0.0) { for (j=0; j<*boolVars; j++) { T += sigma[j]; bVarIndex=Cudd_ReadInvPerm(mgr,j); if (bVarIndex==-1) { bVarIndex=j; } mVarIndex=bVar2mVar[bVarIndex]; eta_rule=eta_temp[vars[mVarIndex].nRule]; for (i=0; i<vars[mVarIndex].nVal-1;i++) { theta=probs[bVarIndex]; eta_rule[i][0]=eta_rule[i][0]+T*(1-theta); eta_rule[i][1]=eta_rule[i][1]+T*theta; } } for (j=0; j<nRules; j++) { for (i=0; i<rules[j]-1; i++) { eta[j][i][0]=eta[j][i][0]+eta_temp[j][i][0]*ex_prob/rootProb; eta[j][i][1]=eta[j][i][1]+eta_temp[j][i][1]*ex_prob/rootProb; } } } free(sigma); return rootProb; } void Maximization(void) { int r,i,j,e; double sum=0; double *probs_rule,**eta_rule; for (r=0;r<nRules;r++) { eta_rule=eta[r]; for (i=0;i<rules[r]-1;i++) { sum=(eta_rule[i][0]+eta_rule[i][1]); if (sum==0.0) { arrayprob[r][i]=0; } else arrayprob[r][i]=eta_rule[i][1]/sum; } } for(e=0;e<ex;e++) { nVars=nVars_ex+e; probs=probs_ex[e]; vars=vars_ex[e]; for (j=0;j<*nVars;j++) { r=vars[j].nRule; probs_rule=arrayprob[r]; for(i=0;i<rules[r]-1;i++) { probs[vars[j].firstBoolVar+i]=probs_rule[i]; } } } } static int randomize(void) { int i,j,e,rule; double * theta,p0; double pmass,par; double **Theta_rules; Theta_rules=(double **)malloc(nRules *sizeof(double *)); for (j=0;j<nRules;j++) { Theta_rules[j]=(double *)malloc(rules[j]*sizeof(double)); } for (j=0;j<nRules;j++) { theta=Theta_rules[j]; pmass=0; for (i=0;i<rules[j]-1;i++) { par=((double)rand())/RAND_MAX*(1-pmass); pmass=pmass+par; theta[i]=par; } theta[rules[j]-1]=1-pmass; } for(e=0;e<ex;e++) { nVars=nVars_ex+e; probs=probs_ex[e]; vars=vars_ex[e]; for (j=0; j<*nVars; j++) { rule=vars[j].nRule; theta=Theta_rules[rule]; p0=1; for (i=0; i<vars[j].nVal-1;i++) { probs[vars[j].firstBoolVar+i]=theta[i]/p0; p0=p0*(1-theta[i]/p0); } } } for (j=0;j<nRules;j++) { free(Theta_rules[j]); } free(Theta_rules); return 1; } static int EM(void) { YAP_Term arg1,arg2,arg3,arg4,arg5,arg6,arg7, out1,out2,nodesTerm,ruleTerm,tail,pair,compoundTerm; DdNode * node1,**nodes_ex; int r,lenNodes,i,iter; long iter1; double CLL0= -2.2*pow(10,10); //-inf double CLL1= -1.7*pow(10,8); //+inf double p,p0,**eta_rule,ea,er; double ratio,diff; arg1=YAP_ARG1; arg2=YAP_ARG2; arg3=YAP_ARG3; arg4=YAP_ARG4; arg5=YAP_ARG5; arg6=YAP_ARG6; arg7=YAP_ARG7; nodesTerm=arg1; ea=YAP_FloatOfTerm(arg2); er=YAP_FloatOfTerm(arg3); lenNodes=YAP_IntOfTerm(arg4); iter=YAP_IntOfTerm(arg5); nodes_ex=(DdNode **)malloc(lenNodes*sizeof(DdNode*)); nodes_probs_ex=(double *)malloc(lenNodes*sizeof(double)); example_prob=(double *)malloc(lenNodes*sizeof(double)); for (i=0;i<lenNodes;i++) { pair=YAP_HeadOfTerm(nodesTerm); node1=(DdNode *)YAP_IntOfTerm(YAP_HeadOfTerm(pair)); nodes_ex[i]=node1; pair=YAP_TailOfTerm(pair); example_prob[i]=YAP_FloatOfTerm(YAP_HeadOfTerm(pair)); nodesTerm=YAP_TailOfTerm(nodesTerm); } diff=CLL1-CLL0; ratio=diff/fabs(CLL0); if (iter==-1) iter1= 2147000000; else iter1=iter; while ( (diff>ea) && (ratio>er) && (cycle<iter1) ) { cycle++; for (r=0;r<nRules;r++) { for (i=0;i<rules[r]-1;i++) { eta_rule=eta[r]; eta_rule[i][0]=0; eta_rule[i][1]=0; } } CLL0 = CLL1; CLL1 = Expectation(nodes_ex,lenNodes); Maximization(); diff=CLL1-CLL0; ratio=diff/fabs(CLL0); } out2= YAP_TermNil(); for (r=0; r<nRules; r++) { tail=YAP_TermNil(); p0=1; for (i=0;i<rules[r]-1;i++) { p=arrayprob[r][i]*p0; tail=YAP_MkPairTerm(YAP_MkFloatTerm(p),tail); p0=p0*(1-arrayprob[r][i]); } tail=YAP_MkPairTerm(YAP_MkFloatTerm(p0),tail); ruleTerm=YAP_MkIntTerm(r); compoundTerm=YAP_MkPairTerm(ruleTerm,YAP_MkPairTerm(tail,YAP_TermNil())); out2=YAP_MkPairTerm(compoundTerm,out2); } out1=YAP_MkFloatTerm(CLL1); YAP_Unify(out1,arg6); free(nodes_ex); free(example_prob); return (YAP_Unify(out2,arg7)); } static int Q(void) { YAP_Term arg1,arg2,arg3,arg4,out,out1, term,nodesTerm,ruleTerm,tail,pair,compoundTerm; DdNode * node1,**nodes_ex; int r,lenNodes,i; double p1,p0,**eta_rule,CLL; arg1=YAP_ARG1; arg2=YAP_ARG2; arg3=YAP_ARG3; arg4=YAP_ARG4; nodesTerm=arg1; lenNodes=YAP_IntOfTerm(arg2); nodes_ex=(DdNode **)malloc(lenNodes*sizeof(DdNode*)); example_prob=(double *)malloc(lenNodes*sizeof(double)); for (i=0;i<lenNodes;i++) { pair=YAP_HeadOfTerm(nodesTerm); node1=(DdNode *)YAP_IntOfTerm(YAP_HeadOfTerm(pair)); nodes_ex[i]=node1; pair=YAP_TailOfTerm(pair); example_prob[i]=YAP_FloatOfTerm(YAP_HeadOfTerm(pair)); nodesTerm=YAP_TailOfTerm(nodesTerm); } for (r=0;r<nRules;r++) { for (i=0;i<rules[r]-1;i++) { eta_rule=eta[r]; eta_rule[i][0]=0; eta_rule[i][1]=0; } } CLL=Expectation(nodes_ex,lenNodes); out= YAP_TermNil(); for (r=0; r<nRules; r++) { tail=YAP_TermNil(); eta_rule=eta[r]; for (i=0;i<rules[r]-1;i++) { p0=eta_rule[i][0]; p1=eta_rule[i][1]; term=YAP_MkPairTerm(YAP_MkFloatTerm(p0), YAP_MkPairTerm(YAP_MkFloatTerm(p1),YAP_TermNil())); tail=YAP_MkPairTerm(term,tail); } ruleTerm=YAP_MkIntTerm(r); compoundTerm=YAP_MkPairTerm(ruleTerm,YAP_MkPairTerm(tail,YAP_TermNil())); out=YAP_MkPairTerm(compoundTerm,out); } free(nodes_ex); free(example_prob); out1=YAP_MkFloatTerm(CLL); YAP_Unify(out1,arg4); return (YAP_Unify(out,arg3)); } static int paths_to_non_zero(void) { double paths; YAP_Term arg1,arg2,out; DdNode * node; arg1=YAP_ARG1; arg2=YAP_ARG2; node=(DdNode *)YAP_IntOfTerm(arg1); paths=Cudd_CountPathsToNonZero(node); out=YAP_MkFloatTerm(paths); return(YAP_Unify(out,arg2)); } static int paths(void) { double paths; YAP_Term arg1,arg2,out; DdNode * node; arg1=YAP_ARG1; arg2=YAP_ARG2; node=(DdNode *)YAP_IntOfTerm(arg1); paths=Cudd_CountPath(node); out=YAP_MkFloatTerm(paths); return(YAP_Unify(out,arg2)); } static int dag_size(void) { int size; YAP_Term arg1,arg2,out; DdNode * node; arg1=YAP_ARG1; arg2=YAP_ARG2; node=(DdNode *)YAP_IntOfTerm(arg1); size=Cudd_DagSize(node); out=YAP_MkIntTerm(size); return(YAP_Unify(out,arg2)); } void init_my_predicates() /* function required by YAP for intitializing the predicates defined by a C function*/ { YAP_UserCPredicate("init",init,2); YAP_UserCPredicate("init_bdd",init_bdd,0); YAP_UserCPredicate("end",end,0); YAP_UserCPredicate("end_bdd",end_bdd,0); YAP_UserCPredicate("add_var",add_var,4); YAP_UserCPredicate("equality",equality,3); YAP_UserCPredicate("and",and,3); YAP_UserCPredicate("one",one,1); YAP_UserCPredicate("zero",zero,1); YAP_UserCPredicate("or",or,3); YAP_UserCPredicate("bdd_not",bdd_not,2); YAP_UserCPredicate("create_dot",create_dot,2); YAP_UserCPredicate("init_test",init_test,1); YAP_UserCPredicate("end_test",end_test,0); YAP_UserCPredicate("ret_prob",ret_prob,2); YAP_UserCPredicate("em",EM,7); YAP_UserCPredicate("q",Q,4); YAP_UserCPredicate("randomize",randomize,0); YAP_UserCPredicate("deref",rec_deref,1); YAP_UserCPredicate("garbage_collect",garbage_collect,2); YAP_UserCPredicate("bdd_to_add",bdd_to_add,2); YAP_UserCPredicate("paths_to_non_zero",paths_to_non_zero,2); YAP_UserCPredicate("paths",paths,2); YAP_UserCPredicate("dag_size",dag_size,2); } FILE * open_file(char *filename, const char *mode) /* opens a file */ { FILE *fp; if ((fp = fopen(filename, mode)) == NULL) { perror(filename); exit(1); } return fp; } tablerow* init_table(int varcnt) { int i; tablerow *tab; tab = (tablerow *) malloc(sizeof(rowel) * varcnt); for (i = 0; i < varcnt; i++) { tab[i].row = NULL; tab[i].cnt = 0; } return tab; } void add_node(tablerow *tab, DdNode *node, double value) { int index = Cudd_NodeReadIndex(node); tab[index].row = (rowel *) realloc(tab[index].row, (tab[index].cnt + 1) * sizeof(rowel)); tab[index].row[tab[index].cnt].key = node; tab[index].row[tab[index].cnt].value = value; tab[index].cnt += 1; } void add_or_replace_node(tablerow *tab, DdNode *node, double value) { int i; int index = Cudd_NodeReadIndex(node); for(i = 0; i < tab[index].cnt; i++) { if (tab[index].row[i].key == node) { tab[index].row[i].value=value; return; } } tab[index].row = (rowel *) realloc(tab[index].row, (tab[index].cnt + 1) * sizeof(rowel)); tab[index].row[tab[index].cnt].key = node; tab[index].row[tab[index].cnt].value = value; tab[index].cnt += 1; } double * get_value(tablerow *tab, DdNode *node) { int i; int index = Cudd_NodeReadIndex(node); for(i = 0; i < tab[index].cnt; i++) { if (tab[index].row[i].key == node) { return &tab[index].row[i].value; } } return NULL; } void destroy_table(tablerow *tab,int varcnt) { int i; for (i = 0; i < varcnt; i++) { free(tab[i].row); } free(tab); }