Add files via upload

1 files changed, 2278 insertions, 0 deletions

diff --git a/sourcecodes/k-best/src/Engine.h b/sourcecodes/k-best/src/Engine.h
new file mode 100644
index 00000000..2e16a701
--- /dev/null
+++ b/sourcecodes/k-best/src/Engine.h
@@ -0,0 +1,2278 @@
+#ifndef ENGINE_H
+#define ENGINE_H
+
+#include<stdio.h>
+#include<stdlib.h>
+
+#include"Arguments.h"
+#include"Model.h"
+
+#include <sys/time.h>
+#include <time.h> 
+
+#define MARK 9.0e99
+
+// Substitutes la += log(1+exp(lb-la)) <=> a = a + b.
+#define LOGADD(la,lb) if(la>8e99||lb-la>100) la=lb;else la+=log(1+exp(lb-la));
+
+#define LOG_ZERO -1.0e101
+
+using namespace std;
+
+#include"UpdateHR.h"
+
+
+
+//----------------
+// Fast Möbius transform routines.
+void sub_fumt(int j, int d, int S, int overlap, double *t, double *s, int n, int k){
+	if (d < n && overlap < k){
+		sub_fumt(j, d+1, S, overlap, t, s, n, k);
+		
+		
+		S |= (1 << d);
+		overlap += (d >= j+1);
+		sub_fumt(j, d+1, S, overlap, t, s, n, k); 
+	}
+	
+	else{
+		//base case: 
+		//since s[S] = MARK initially,  and the following 3 cases, s[S] must be assigned some score (neither MARK nor undecided)
+		s[S] = MARK;
+		
+		//jinS: j is in S
+		int jinS = ((S >> j) & 1);
+		if (overlap + jinS <= k) s[S] = t[S];
+		//if (jinS) LOGADD(s[S], t[S - (1 << j)]);
+		if (jinS) logAdd_New(s[S], t[S - (1 << j)]);
+	}
+}
+
+// Fast upward Möbius transform:
+// s(S) := \sum_{T \subseteq S : |T| \leq k} t(T).
+void fumt(double *t, double *s, int n, int k){
+	double *tmp, *sprev = new double[1 << n];
+	for (int T = 0; T < (1<<n); T ++) 
+		sprev[T] = t[T];
+	for (int j = 0; j < n; j ++){
+		//for each j
+		sub_fumt(j, 0, 0, 0, sprev, s, n, k);
+		//recyclely use sprev, s
+		tmp = sprev; 
+		sprev = s; 
+		s = tmp;
+	}
+	if (n % 2 == 0){ 
+		for (int S = 0; S < (1<<n); S ++){ 
+			s[S] = sprev[S];
+		}
+	}
+	else{ 
+		tmp = sprev; 
+		sprev = s; 
+		s = tmp;
+	}
+	delete [] sprev;
+}
+void test_fumt(){
+	int n = 4;
+	double t[16], s[16];
+	for (int S = 0; S < 16; S ++){
+		t[S] = log(S+1);
+	}
+	fumt(t, s, n, 2);
+	for (int S = 0; S < 16; S ++){
+		cerr<<" s["<<S<<"] = "<<exp(s[S])<<endl;
+	}
+}
+
+
+void sub_fdmt(int j, int d, int T, int overlap, double *s, double *t, int n, int k){
+	//cerr<<" T = "<<T<<":"; print_set(cerr, T); 
+	//cerr<<" overlap = "<<overlap<<endl; 
+	if (d < n && overlap <= k){
+		sub_fdmt(j, d+1, T, overlap, s, t, n, k);
+		T |= (1 << d);
+		overlap += (d <= j);
+		sub_fdmt(j, d+1, T, overlap, s, t, n, k); 
+	}
+	else if (d == n){
+		t[T] = s[T];
+		int jinT = ((T >> j) & 1);
+		//if (jinT == 0) LOGADD(t[T], s[T + (1 << j)]);
+		if (jinT == 0) logAdd_New(t[T], s[T + (1 << j)]);
+		if (overlap > k) t[T] = MARK;
+	}
+}
+// Fast downward Möbius transform:
+// t(T) := \sum_{T \subseteq S} s(S).
+void fdmt(double *s, double *t, int n, int k){
+	double *tmp, *tprev = new double[1 << n];
+	for (int S = 0; S < (1<<n); S ++) tprev[S] = s[S];
+	for (int j = 0; j < n; j ++){
+		
+		sub_fdmt(j, 0, 0, 0, tprev, t, n, k);	
+		tmp = t; t = tprev; tprev = tmp;
+	}
+	if (n % 2 == 0){ for (int T = 0; T < (1<<n); T ++) t[T] = tprev[T];}
+	else{ tmp = t; t = tprev; tprev = tmp;}
+	delete [] tprev;
+}
+
+
+void sub_fdmtLogAddComp(int j, int d, int T, int overlap, double *s, double *t, int n, int k){
+	//cerr<<" T = "<<T<<":"; print_set(cerr, T); 
+	//cerr<<" overlap = "<<overlap<<endl; 
+	if (d < n && overlap <= k){
+		sub_fdmtLogAddComp(j, d+1, T, overlap, s, t, n, k);
+		T |= (1 << d);
+		overlap += (d <= j);
+		sub_fdmtLogAddComp(j, d+1, T, overlap, s, t, n, k); 
+	}
+	else if (d == n){
+		t[T] = s[T];
+		int jinT = ((T >> j) & 1);
+		if (jinT == 0) logAddComp(t[T], s[T + (1 << j)]);
+		if (overlap > k) t[T] = MARK;
+	}
+}
+// Fast downward Möbius transform:
+// t(T) := \sum_{T \subseteq S} s(S).
+void fdmtLogAddComp(double *s, double *t, int n, int k){
+	double *tmp, *tprev = new double[1 << n];
+	for (int S = 0; S < (1<<n); S ++) tprev[S] = s[S];
+	for (int j = 0; j < n; j ++){
+		
+		sub_fdmtLogAddComp(j, 0, 0, 0, tprev, t, n, k);	
+		tmp = t; t = tprev; tprev = tmp;
+	}
+	if (n % 2 == 0){ for (int T = 0; T < (1<<n); T ++) t[T] = tprev[T];}
+	else{ tmp = t; t = tprev; tprev = tmp;}
+	delete [] tprev;
+}
+
+
+void sub_fdmtNoLog(int j, int d, int T, int overlap, double *s, double *t, int n, int k){
+	//cerr<<" T = "<<T<<":"; print_set(cerr, T); 
+	//cerr<<" overlap = "<<overlap<<endl; 
+	if (d < n && overlap <= k){
+		sub_fdmtNoLog(j, d+1, T, overlap, s, t, n, k);
+		T |= (1 << d);
+		overlap += (d <= j);
+		sub_fdmtNoLog(j, d+1, T, overlap, s, t, n, k); 
+	}
+	else if (d == n){
+		t[T] = s[T];
+		int jinT = ((T >> j) & 1);
+		//if (jinT == 0) LOGADD(t[T], s[T + (1 << j)]);
+		//if (jinT == 0) logAdd_New(t[T], s[T + (1 << j)]);
+		if (jinT == 0) t[T] += s[T + (1 << j)];
+		if (overlap > k) t[T] = MARK;
+	}
+}
+
+// Fast downward Möbius transform:
+// t(T) := \sum_{T \subseteq S} s(S).
+void fdmtNoLog(double *s, double *t, int n, int k){
+	double *tmp, *tprev = new double[1 << n];
+	for (int S = 0; S < (1<<n); S ++) tprev[S] = s[S];
+	for (int j = 0; j < n; j ++){
+		
+		sub_fdmtNoLog(j, 0, 0, 0, tprev, t, n, k);	
+		tmp = t; t = tprev; tprev = tmp;
+	}
+	if (n % 2 == 0){ for (int T = 0; T < (1<<n); T ++) t[T] = tprev[T];}
+	else{ tmp = t; t = tprev; tprev = tmp;}
+	delete [] tprev;
+}
+
+
+void test_fdmt(){
+	int n = 4;
+	double t[16], s[16];
+	for (int S = 0; S < 16; S ++){
+		s[S] = log(S+1);
+	}
+	fdmt(s, t, n, 3);
+	for (int T = 0; T < 16; T ++){
+		cerr<<" t["<<T<<"] = "<<exp(t[T])<<endl;
+	}
+}
+
+//----------------
+
+
+
+
+
+
+// NOTE: Engine only calls a veru limited set of interface functions
+// implemented in Model.
+class Engine {
+public:
+	Engine(){}
+	~Engine(){}	
+	
+	void init(Model *m){
+		model = m;
+		//test();
+	
+	}
+	void test(){
+		vector<int> T;
+		T.push_back(1); T.push_back(2);
+		
+		cerr<<" log p(x_0 | x_1,2): "<<model->log_lcp(0, T)<<endl;
+	
+		for (int h = 0; h < model->num_layers(); h ++){
+			int *Vh, *Vu, *Vl;
+			int nh, nu, nl;
+			model->layer(h, &Vh, &nh);	
+			model->upper_layers(h, &Vu, &nu);	
+			model->lower_layers(h, &Vl, &nl);
+			cerr<<"Layer "<<h<<":";
+			cerr<<endl<<" layer: ";
+			print_nodes(cerr, Vh, nh);
+			cerr<<endl<<" upper: ";
+			print_nodes(cerr, Vu, nu);
+			cerr<<endl<<" lower: ";
+			print_nodes(cerr, Vl, nl);
+			cerr<<endl;	
+		}
+		
+		//test_fumt();
+		test_fdmt(); 
+		exit(1);
+
+	}
+	// Computes all edge probabilities.
+	void compute_edge_probabilities(){
+//		int l = model->num_layers();
+	
+		//Arguments::option == 0 means the rebel method with new Dirichlet hyper-param (see sub_beta)
+		if(Arguments::option == 0){
+			compute_edge_probabilities(0);
+		}
+		//Arguments::option == 1 means the forward method
+		else if(Arguments::option == 1){
+			compute_edge_probabilities_forward(0);
+		}
+		//Arguments::option == 2 means the backward method
+		else if(Arguments::option == 2){
+			compute_edge_probabilities_backward(0);
+		}
+		//Arguments::option == 3 means the mix method in the UAI paper
+		else if(Arguments::option == 3){
+			compute_edge_probabilities_mixIndegree(0);
+		 	//compute_edge_probabilities_mixIndegreeNoLog(0);
+		}
+		//Arguments::option == 4 means the computation used in top-k
+		else if(Arguments::option == 4){
+			compute_edge_probabilities_top_k(0);
+		}
+		//Arguments::option == 5 means the computation required for in and out features
+		else if(Arguments::option == 5){
+			//For In_Out features
+			compute_edge_probabilities_in_out_features(0);
+		}
+		
+		
+	}
+	// Computes probabilities for edges pointing to the h-th layer.
+	void compute_edge_probabilities(int h){
+		cerr<<"\nREBEL Method:" << endl;
+		cerr<<" Compute edge probabilities for Layer "<<h<<":"<<endl;
+		
+		// Step 1: Compute beta[][];
+		// Step 2: Compute alpha[][];
+		int *Vh, *Vu; int nh, nu;
+		model->layer(h, &Vh, &nh);
+		model->upper_layers(h-1, &Vu, &nu);
+
+		
+		beta = new double*[nh];
+		alpha = new double*[nh];
+		
+		k = model->max_indegree();
+		
+		cerr<<" . "<<nh<<" nodes"<<endl;
+		
+		struct timeval startTime;
+		struct timeval endTime;
+		struct timezone myTimeZone;
+
+		
+		double betaTime;
+		double alphaTime;
+		double forwardTime;
+		double backwardTime;
+		double gammaTime;
+		
+		double step2BTime;
+	
+		gettimeofday(&startTime, &myTimeZone);		
+
+		//If Arguments::ADtree == 1, user wants to use ADtree, then make the whole tree first
+		if(Arguments::ADtree == 1){
+			model->makeADTree();
+		}
+		
+		for (int j = 0; j < nh; j ++){
+
+			beta[j] = new double[1 << nh];
+			//int i = Vh[j];
+			compute_beta(j, Vh, nh, Vu, nu);
+		
+		}
+		
+		//If Arguments::ADtree == 1, then delete the whole tree
+		if(Arguments::ADtree == 1){
+			model->freeADTree();
+		}
+		
+		gettimeofday(&endTime, &myTimeZone);		
+		betaTime = (endTime.tv_sec - startTime.tv_sec) + 
+		           (endTime.tv_usec - startTime.tv_usec) / 1000000.0 ;
+		           
+		cerr<<" . Tables beta are now ready."<<endl;
+		
+		gettimeofday(&startTime, &myTimeZone);
+		for (int j = 0; j < nh; j ++){
+
+			alpha[j] = new double[1 << nh];	
+			compute_alpha(j, Vh, nh);
+		
+			//cerr<<" . Tables beta and alpha computed for node "<<j<<"."<<endl; 
+		
+		}
+		gettimeofday(&endTime, &myTimeZone);		
+		alphaTime = (endTime.tv_sec - startTime.tv_sec) + 
+		           (endTime.tv_usec - startTime.tv_usec) / 1000000.0 ;
+		
+		cerr<<" . Tables alpha are now ready."<<endl;
+
+			
+		// Step 3: Compute g_forward[];
+		gf = new double[1 << nh];
+		
+		gettimeofday(&startTime, &myTimeZone);
+		
+		compute_g_forward(nh);
+		
+		gettimeofday(&endTime, &myTimeZone);		
+		forwardTime = (endTime.tv_sec - startTime.tv_sec) + 
+		           (endTime.tv_usec - startTime.tv_usec) / 1000000.0 ;
+		
+		//cerr<<" . gf[all] = "<<gf[(1<<nh)-1]<<endl;
+		
+		// Step 4: Compute g_backward[];
+		
+		gb = new double[1 << nh];
+		
+		gettimeofday(&startTime, &myTimeZone);
+		
+		compute_g_backward(nh);
+		
+		gettimeofday(&endTime, &myTimeZone);		
+		backwardTime = (endTime.tv_sec - startTime.tv_sec) + 
+		           (endTime.tv_usec - startTime.tv_usec) / 1000000.0 ;
+		
+		//cerr<<" . gb[all] = "<<gb[(1<<nh)-1]<<endl;
+	
+		// Step 5: Loop over end-point nodes j.
+		
+		gammaTime = 0.0;
+		step2BTime = 0.0;
+		
+		double *gamma = new double[1 << nh];
+		double *gfb = new double[1 << nh];
+		for (int j = 0; j < nh; j ++){
+			// Compute gfb[]. 
+			
+			//start counting gammaTime for each j
+			gettimeofday(&startTime, &myTimeZone);
+			
+			for (int S = 0; S < (1<<nh); S ++){
+				
+				//if j \in the binary repre of S
+				if ((S>>j) & 1){
+					gfb[S] = -MARK;
+				}
+				else{
+					int cS = (1 << nh) - 1 - S - (1 << j);
+					gfb[S] = gf[S] + gb[cS];
+					//cerr<<" j = "<<j<<": ";
+					//print_set(cerr, S);
+					//cerr<<"; ";
+					//print_set(cerr, cS);
+					//cerr<<": "<<gfb[S];
+					//cerr<<endl;
+				}
+			}
+			// Compute gamma_j() := gamma[].
+			fdmt(gfb, gamma, nh, k);
+			
+			gettimeofday(&endTime, &myTimeZone);		
+			gammaTime += (endTime.tv_sec - startTime.tv_sec) + 
+		           (endTime.tv_usec - startTime.tv_usec) / 1000000.0 ;
+			//end counting gammaTime for each j
+			
+			//cerr<<" gamma[0]: "<<gamma[0]<<endl;
+			//cerr<<" gfb[0]: "<<gfb[0]<<endl;
+			
+			cerr<<" . Incoming edges for node "<<j<<":";
+			// Loop over start-point nodes i.
+			
+			//start counting step2BTime for each j
+			gettimeofday(&startTime, &myTimeZone);
+			
+			for (int i = 0; i < nh; i ++){
+				if (i == j) continue;
+				cerr<<" "<<i;		
+				double log_prob = eval_edge(i,  j, gamma, beta[j], nh, k);
+		
+				model->print_edge_prob(
+					cout, Vh[i], Vh[j], exp(log_prob - gb[(1<<nh)-1])); //HR: - gb[(1<<nh)-1] mean /L(v)
+			
+			}
+			gettimeofday(&endTime, &myTimeZone);		
+			step2BTime += (endTime.tv_sec - startTime.tv_sec) + 
+		           (endTime.tv_usec - startTime.tv_usec) / 1000000.0 ;
+			//end counting step2BTime for each j
+			
+			
+			
+			cerr<<endl;
+			
+		} 
+		
+		
+		cout << "\n\n\nTime Statistics" << endl;
+		cout << "betaTime = " << betaTime <<  endl;				
+		cout << "alphaTime = " << alphaTime <<  endl;
+		cout << "forwardTime = " << forwardTime <<  endl;
+		cout << "backwardTime = " << backwardTime <<  endl;				
+		cout << "gammaTime = " << gammaTime <<  endl;
+		
+		cout << "step2BTime = " << step2BTime <<  endl;
+		
+		double totalTime = betaTime + alphaTime + forwardTime 
+						+ backwardTime + gammaTime + step2BTime;
+		cout << "totalTime = " << totalTime <<  endl;
+		
+		/*
+		ofResult << "\n\n\nTime Statistics" << endl;
+		ofResult << "betaTime = " << betaTime <<  endl;				
+		ofResult << "alphaTime = " << alphaTime <<  endl;
+		ofResult << "RRFuncTime = " << RRFuncTime <<  endl;
+		ofResult << "HFuncTime = " << HFuncTime <<  endl;
+		ofResult << "KFuncTime = " << KFuncTime <<  endl;		
+				
+		ofResult << "GammaTime = " << GammaTime <<  endl;
+		ofResult << "step2Time = " << step2Time <<  endl;
+		double totalTime = betaTime + alphaTime + RRFuncTime 
+						+ HFuncTime + KFuncTime + GammaTime + step2Time;
+		ofResult << "totalTime = " << totalTime <<  endl;
+		*/
+		
+		
+		delete [] gamma;
+		delete [] gfb;
+	
+		delete [] gf; delete [] gb;
+		for (int j = 0; j < nh; j ++){
+			delete [] beta[j]; delete [] alpha[j];
+		}
+		delete [] beta; delete [] alpha;
+		cerr<<" Edge probabilities now computed."<<endl; 
+	}
+	
+	
+	void compute_edge_probabilities_ind(int h){
+		cerr<<"\nREBEL Method with each individual edge:" << endl;
+		cerr<<" Compute edge probabilities for Layer "<<h<<":"<<endl;
+				
+		// Step 1: Compute beta[][];
+		// Step 2: Compute alpha[][];
+		int *Vh, *Vu; int nh, nu;
+		model->layer(h, &Vh, &nh);
+		model->upper_layers(h-1, &Vu, &nu);
+				
+		beta = new double*[nh];
+		alpha = new double*[nh];
+		
+		betaNume = new double*[nh];
+		alphaNume = new double*[nh];
+		
+		k = model->max_indegree();
+		
+		cerr<<" . "<<nh<<" nodes"<<endl;
+		
+		for (int j = 0; j < nh; j ++){
+
+			beta[j] = new double[1 << nh];
+
+			//int i = Vh[j];
+			compute_beta(j, Vh, nh, Vu, nu);
+			
+			alpha[j] = new double[1 << nh];	
+			compute_alpha(j, Vh, nh);
+			
+			//alphaNume[j] = new double[1 << nh];	
+			//compute_alphaNume(j, Vh, nh);
+		
+			cerr<<" . Tables beta and alpha computed for node "<<j<<"."<<endl; 
+		
+		}
+	
+	  for(int uu = 0; uu < nh; uu++){
+	   for(int vv = 0; vv < nh; vv++){
+	  	if(vv != uu){
+
+		cerr<<" Compute edge ("<< uu <<"<-" << vv <<")" <<endl;
+
+		for (int j = 0; j < nh; j ++){
+			betaNume[j] = new double[1 << nh];
+			compute_betaNume(j, Vh, nh, Vu, nu, uu, vv);
+			
+			alphaNume[j] = new double[1 << nh];	
+			compute_alphaNume(j, Vh, nh);
+		
+			cerr<<" . Tables betaNume and alphaNume computed for node "<<j<<"."<<endl; 
+		
+		}
+	
+		cerr<<" . Tables beta/betaNume and alpha/alphaNume are now ready."<<endl;
+
+			
+		// Step 3: Compute g_forward[];
+		gf = new double[1 << nh];
+
+		compute_g_forward(nh);
+
+		gfNume = new double[1 << nh];
+		compute_gNume_forward(nh);
+		
+		cerr << "\nalpha" << endl;
+		for (int j = 0; j < nh; j ++){		
+			for(int index=0; index < (1 << nh); index++){
+				cerr<< "alpha[" << j << "][" << index << "] =" <<  alpha[j][index] << endl;
+			}
+		}
+		
+		cerr << "\nalphaNume" << endl;
+		for (int j = 0; j < nh; j ++){		
+			for(int index=0; index < (1 << nh); index++){
+				cerr<< "alphaNume[" << j << "][" << index << "] =" <<  alphaNume[j][index] << endl;
+			}
+		}
+		
+		cerr << "\ngf:" << endl;
+		for(int index = 0; index < (1 << nh); index++){
+			cerr << "gf[" << index << "] = " << gf[index] << endl;
+		}
+		
+		cerr << "\ngfNume:" << endl;
+		for(int index = 0; index < (1 << nh); index++){
+			cerr << "gfNume[" << index << "] = " << gfNume[index] << endl;
+		}
+
+		cerr<<"gf[all] = "<<gf[(1<<nh)-1]<<endl;
+		cerr<<"gfNume[all] = "<<gfNume[(1<<nh)-1]<<endl;
+		
+		//For checking: Compute g_backward[];
+		
+		gb = new double[1 << nh];
+		compute_g_backward(nh);	
+		cerr<<" gb[all] = "<<gb[(1<<nh)-1]<<endl;
+				
+		cerr<<" . Incoming edges from node " << vv << " to node " << uu <<" :";
+		model->print_edge_prob(cout, Vh[vv], Vh[uu], exp(gfNume[(1<<nh)-1] - gf[(1<<nh)-1])); //HR: - gb[(1<<nh)-1] mean /L(v)		
+		cerr<<endl;
+		
+	  	}
+	   }
+	  }
+		
+		delete [] gf; 
+		delete [] gfNume; 
+		for (int j = 0; j < nh; j ++){
+			delete [] beta[j]; 
+			delete [] betaNume[j];
+			delete [] alpha[j];
+			delete [] alphaNume[j];
+		}
+		delete [] beta; 
+		delete [] betaNume; 
+		delete [] alpha;
+		delete [] alphaNume;
+		
+		cerr<<" Edge probabilities now computed."<<endl; 
+		
+		return;
+		
+		
+
+	}
+	
+		
+	//HR: new forward method, compute for each individual edge
+	void compute_edge_probabilities_forward(int h){
+		cerr<<"Forward method: " << endl;
+		cerr<<" Compute edge probabilities for Layer "<<h<<":"<<endl;
+				
+		// Step 1: Compute beta[][];
+		// Step 2: Compute alpha[][];
+		int *Vh, *Vu; int nh, nu;
+		model->layer(h, &Vh, &nh);
+		model->upper_layers(h-1, &Vu, &nu);
+				
+		//If Arguments::ADtree == 1, user wants to use ADtree, then make the whole tree first
+		if(Arguments::ADtree == 1){
+			model->makeADTree();
+		}
+
+	
+		beta = new double*[nh];
+		alpha = new double*[nh];
+		
+		betaNume = new double*[nh];
+		alphaNume = new double*[nh];
+		
+		k = model->max_indegree();
+		
+		cerr<<" . "<<nh<<" nodes"<<endl;
+		
+		for (int j = 0; j < nh; j ++){
+
+			beta[j] = new double[1 << nh];
+					
+			//int i = Vh[j];
+			compute_beta(j, Vh, nh, Vu, nu);
+
+			alpha[j] = new double[1 << nh];	
+			compute_alpha(j, Vh, nh);
+				
+			cerr<<" . Tables beta and alpha computed for node "<<j<<"."<<endl; 
+		
+		}
+		
+		hFunc = new double[1 << nh];
+		compute_h_Fast(nh, hFunc, alpha);
+		cerr<<"hFunc[all] = "<<hFunc[(1<<nh)-1]<<endl;
+
+		
+	    for(int uu = 0; uu < nh; uu++){
+	    	for(int vv = 0; vv < nh; vv++){
+	  			if(vv != uu){
+	  		
+					cerr<<"Compute edge ("<< uu <<"<-" << vv <<")" <<endl;
+	
+					for (int j = 0; j < nh; j ++){
+						betaNume[j] = new double[1 << nh];
+
+						compute_betaNume(j, Vh, nh, Vu, nu, uu, vv);
+			
+						alphaNume[j] = new double[1 << nh];	
+						compute_alphaNume(j, Vh, nh);
+		
+						//cerr<<" . Tables betaNume and alphaNume computed for node "<<j<<"."<<endl; 
+		
+					}
+		
+		
+					cerr<<" . Tables beta/betaNume and alpha/alphaNume are now ready."<<endl;
+
+		
+					hFuncNume = new double[1 << nh];
+					compute_h_Fast(nh, hFuncNume, alphaNume);
+
+					cerr<<"hFunc[all] = "<<hFunc[(1<<nh)-1]<<endl;
+					cerr<<"hFuncNume[all] = "<<hFuncNume[(1<<nh)-1]<<endl;
+				
+					cerr<<" . Incoming edges from node " << vv << " to node " << uu <<" :";
+					model->print_edge_prob(cout, Vh[vv], Vh[uu], exp(hFuncNume[(1<<nh)-1] - hFunc[(1<<nh)-1])); //HR: - gb[(1<<nh)-1] mean /L(v)		
+					cerr<<endl;
+		
+	  			}
+	   		}
+	  	}
+		
+		
+		//If Arguments::ADtree == 1, then delete the whole tree
+		if(Arguments::ADtree == 1){
+			model->freeADTree();
+		}
+		
+		delete [] hFunc; 
+		
+		//it better to move the following line inside the double for loop
+		delete [] hFuncNume; 
+		
+		for (int j = 0; j < nh; j ++){
+			delete [] beta[j]; 
+			delete [] betaNume[j];
+			delete [] alpha[j];
+			delete [] alphaNume[j];
+		}
+		delete [] beta; 
+		delete [] betaNume; 
+		delete [] alpha;
+		delete [] alphaNume;
+		cerr<<" Edge probabilities now computed."<<endl; 
+		
+		return;
+		
+		
+
+	}
+	
+	
+	//HR: new backward method, compute for each individual edge
+	void compute_edge_probabilities_backward(int h){
+		cerr<<"Backward method: " << endl;
+		cerr<<" Compute edge probabilities for Layer "<<h<<":"<<endl;
+				
+		// Step 1: Compute beta[][];
+		// Step 2: Compute alpha[][];
+		int *Vh, *Vu; int nh, nu;
+		model->layer(h, &Vh, &nh);
+		model->upper_layers(h-1, &Vu, &nu);
+		
+		//If Arguments::ADtree == 1, user wants to use ADtree, then make the whole tree first
+		if(Arguments::ADtree == 1){
+			model->makeADTree();
+		}
+	
+		beta = new double*[nh];
+		alpha = new double*[nh];
+		
+		betaNume = new double*[nh];
+		alphaNume = new double*[nh];
+		
+		k = model->max_indegree();
+		
+		cerr<<" . "<<nh<<" nodes"<<endl;
+		
+		for (int j = 0; j < nh; j ++){
+
+			beta[j] = new double[1 << nh];
+
+			compute_beta(j, Vh, nh, Vu, nu);
+
+			alpha[j] = new double[1 << nh];	
+			compute_alpha(j, Vh, nh);
+			
+			cerr<<" . Tables beta and alpha computed for node "<<j<<"."<<endl; 
+		
+		}
+		
+		RRFunc = new double[1 << nh];
+		compute_RR_Fast(nh, RRFunc, alpha);
+		cerr<<"RRFunc[all] = "<<RRFunc[(1<<nh)-1]<<endl;
+		
+					
+		//int uu;
+		//int vv;
+	    //uu = 2;
+		//vv = 1;	
+	    for(int uu = 0; uu < nh; uu++){
+	    	for(int vv = 0; vv < nh; vv++){
+	  			if(vv != uu){
+	  		
+					cerr<<"Compute edge ("<< uu <<"<-" << vv <<")" <<endl;
+	
+					for (int j = 0; j < nh; j ++){
+
+						betaNume[j] = new double[1 << nh];
+						compute_betaNume(j, Vh, nh, Vu, nu, uu, vv);
+			
+						alphaNume[j] = new double[1 << nh];	
+						compute_alphaNume(j, Vh, nh);
+		
+						cerr<<" . Tables betaNume and alphaNume computed for node "<<j<<"."<<endl; 
+		
+					}
+				
+					cerr<<" . Tables beta/betaNume and alpha/alphaNume are now ready."<<endl;
+		
+					RRFuncNume = new double[1 << nh];
+					compute_RR_Fast(nh, RRFuncNume, alphaNume);
+
+					cerr<<"RRFunc[all] = "<<RRFunc[(1<<nh)-1]<<endl;
+					cerr<<"RRFuncNume[all] = "<<RRFuncNume[(1<<nh)-1]<<endl;
+
+					cerr<<" . Incoming edges from node " << vv << " to node " << uu <<" :";
+					model->print_edge_prob(cout, Vh[vv], Vh[uu], exp(RRFuncNume[(1<<nh)-1] - RRFunc[(1<<nh)-1])); //HR: - gb[(1<<nh)-1] mean /L(v)		
+					cerr<<endl;
+		
+	  			}
+	   		}
+	  	}
+		
+		//If Arguments::ADtree == 1, then delete the whole tree
+		if(Arguments::ADtree == 1){
+			model->freeADTree();
+		}
+		
+		delete [] RRFunc;
+		
+		//it better to move the following line inside the double for loop 
+		delete [] RRFuncNume; 
+		
+		for (int j = 0; j < nh; j ++){
+			delete [] beta[j]; 
+			delete [] betaNume[j];
+			delete [] alpha[j];
+			delete [] alphaNume[j];
+		}
+		delete [] beta; 
+		delete [] betaNume; 
+		delete [] alpha;
+		delete [] alphaNume;
+		cerr<<" Edge probabilities now computed."<<endl; 
+		
+		return;
+		
+		
+
+	}
+	
+
+	void compute_edge_probabilities_in_out_features(int h){
+		cerr<<"Backward method for in and out features: " << endl;
+		//cerr<<" Compute poster probabilities for Layer "<<h<<":"<<endl;
+				
+		// Step 1: Compute beta[][];
+		// Step 2: Compute alpha[][];
+		int *Vh, *Vu; int nh, nu;
+		model->layer(h, &Vh, &nh);
+		model->upper_layers(h-1, &Vu, &nu);
+		
+		beta = new double*[nh];
+		alpha = new double*[nh];
+		
+		betaInOut = new double*[nh];
+		alphaInOut = new double*[nh];
+		
+		k = model->max_indegree();
+		
+		cerr<<" . "<<nh<<" nodes"<<endl;
+
+		//If Arguments::ADtree == 1, user wants to use ADtree, then make the whole tree first
+		if(Arguments::ADtree == 1){
+			model->makeADTree();
+		}
+
+		for (int j = 0; j < nh; j ++){
+
+			beta[j] = new double[1 << nh];
+			compute_beta(j, Vh, nh, Vu, nu);			
+
+			alpha[j] = new double[1 << nh];	
+			compute_alpha(j, Vh, nh);
+				
+			cerr<<" . Tables beta and alpha computed for node "<<j<<"."<<endl; 
+		
+		}
+		
+	
+		inFeatureEdges = new int[nh];
+		outFeatureEdges = new int[nh];	
+		//cerr << "call: getInOutFeatures(" << nh << ")" << endl;	
+		getInOutFeatures(nh);		
+		//cerr << "return: getInOutFeatures()"  << endl;
+		
+//		for(int i = 0; i < nh; i++){
+//			cerr << "inFeatureEdges[" << i << "] = " << inFeatureEdges[i] << endl;
+//		}
+//		for(int i = 0; i < nh; i++){
+//			cerr << "outFeatureEdges[" << i << "] = " << outFeatureEdges[i] << endl;
+//		}
+		
+		for (int j = 0; j < nh; j ++){
+			//cerr<<"\nj = " << j << endl;
+			betaInOut[j] = new double[1 << nh];
+			compute_betaInOut(j, Vh, nh, Vu, nu);
+
+			alphaInOut[j] = new double[1 << nh];			
+			compute_alphaInOut(j, Vh, nh);
+
+			cerr<<" . Tables betaInOut and alphaInOut computed for node "<<j<<"."<<endl; 
+
+		}
+				
+		cerr<<" Tables beta/betaInOut and alpha/alphaInOut are now ready."<<endl;
+
+		cerr << " Compute the postier probability for the specified In-Out features:" << endl;
+		RRFunc = new double[1 << nh];
+		compute_RR_Fast(nh, RRFunc, alpha);
+		//compute_h_Fast(nh, RRFunc, alpha); //To confirm that the result is the same for forward and backward
+
+		RRFuncInOut = new double[1 << nh];
+		compute_RR_Fast(nh, RRFuncInOut, alphaInOut);
+		//compute_h_Fast(nh, RRFuncInOut, alphaInOut); //To confirm that the result is the same for forward and backward
+
+		cerr<<"RRFunc[all] = "<<RRFunc[(1<<nh)-1]<<endl;
+		cerr<<"RRFuncInOut[all] = "<<RRFuncInOut[(1<<nh)-1]<<endl;
+
+		//cerr<<" . Incoming edges from node " << vv << " to node " << uu <<" :";
+		//model->print_edge_prob(cout, Vh[vv], Vh[uu], exp(RRFuncNume[(1<<nh)-1] - RRFunc[(1<<nh)-1])); //HR: - gb[(1<<nh)-1] mean /L(v)
+		cerr << "\nThe postier probability for the specified In-Out features: " <<  exp(RRFuncInOut[(1<<nh)-1] - RRFunc[(1<<nh)-1]);		
+		cerr<<endl;
+		
+		cout << "The postier probability for the specified In-Out features: " <<  exp(RRFuncInOut[(1<<nh)-1] - RRFunc[(1<<nh)-1]);		
+		cout<<endl;
+		
+		//If Arguments::ADtree == 1, then delete the whole tree
+		if(Arguments::ADtree == 1){
+			model->freeADTree();
+		}
+		
+		delete [] RRFunc; 
+		delete [] RRFuncInOut; 
+	
+		
+		delete [] inFeatureEdges;
+		delete [] outFeatureEdges;	
+		
+		for (int j = 0; j < nh; j ++){
+			delete [] beta[j]; 
+			delete [] betaInOut[j];
+			delete [] alpha[j];
+			delete [] alphaInOut[j];
+		}
+		delete [] beta; 
+		delete [] betaInOut; 
+		delete [] alpha;
+		delete [] alphaInOut;
+		cerr<<" Poster probabilities now computed."<<endl; 
+		
+		return;
+		
+		
+	}
+
+
+	void compute_edge_probabilities_top_k(int h){
+		cerr<<"For top-k: compute p(D): " << endl;
+		cerr<<" Compute edge probabilities for Layer "<<h<<":"<<endl;		
+		
+		// Step 1: Compute beta[][];
+		// Step 2: Compute alpha[][];
+		int *Vh, *Vu; int nh, nu;
+		model->layer(h, &Vh, &nh);
+		model->upper_layers(h-1, &Vu, &nu);
+		
+		
+		beta = new double*[nh];
+		alpha = new double*[nh];
+		
+		//betaNume = new double*[nh];
+		//alphaNume = new double*[nh];
+		
+		k = model->max_indegree();
+		
+		cerr<<" . "<<nh<<" nodes"<<endl;
+		
+		//If Arguments::ADtree == 1, user wants to use ADtree, then make the whole tree first
+		if(Arguments::ADtree == 1){
+			model->makeADTree();
+		}
+
+		for (int j = 0; j < nh; j ++){
+
+			beta[j] = new double[1 << nh];
+		
+			//cerr<<" Start compute_beta for node "<<j<<"."<<endl;
+			compute_beta(j, Vh, nh, Vu, nu);
+			//cerr<<" Tables beta computed for node "<<j<<"."<<endl;
+
+			alpha[j] = new double[1 << nh];	
+			compute_alpha(j, Vh, nh);
+				
+			cerr<<" . Tables beta and alpha computed for node "<<j<<"."<<endl; 
+		
+		}
+
+		//If Arguments::ADtree == 1, then delete the whole tree
+		if(Arguments::ADtree == 1){
+			model->freeADTree();
+		}
+			
+		RRFunc = new double[1 << nh];
+		
+		cerr<<"\n Compute_RR_Fast(): "<<endl;		
+		compute_RR_Fast(nh, RRFunc, alpha);
+		cerr<<"\n RRFunc[all] = "<<RRFunc[(1<<nh)-1]<<endl <<endl;
+		cout<<"\n RRFunc[all] = "<<RRFunc[(1<<nh)-1]<<endl <<endl;
+
+		fprintf(stderr, "RRFunc[all] = %18.8f\n", RRFunc[(1<<nh)-1]);
+		
+		//write to 
+		FILE * fp_PD = fopen("exactPD.txt", "a");
+		if(fp_PD){
+			fprintf(fp_PD, "%18.8f\n", RRFunc[(1<<nh)-1]);
+			fclose(fp_PD);
+			
+		}
+		else{
+			cerr <<"fopen exactPD.txt fails" << endl;
+			
+		}
+		
+	
+		delete [] RRFunc; 
+		//delete [] RRFuncNume; 
+		
+		for (int j = 0; j < nh; j ++){
+			delete [] beta[j]; 
+			//delete [] betaNume[j];
+			delete [] alpha[j];
+			//delete [] alphaNume[j];
+		}
+		delete [] beta; 
+		//delete [] betaNume; 
+		delete [] alpha;
+		//delete [] alphaNume;
+		//cerr<<" Edge probabilities now computed."<<endl; 
+		
+		return;
+
+	}
+
+		
+	//HR: new mix with max in-degree method, 
+	void compute_edge_probabilities_mixIndegree(int h){
+		cerr<<" Mix with max in-degree method: " << endl;
+		//cerr<<" Compute edge probabilities for Layer "<<h<<":"<<endl;
+				
+		// Step 1: Compute beta[][];
+		// Step 2: Compute alpha[][];
+		int *Vh, *Vu; int nh, nu;
+		model->layer(h, &Vh, &nh);
+		model->upper_layers(h-1, &Vu, &nu);
+		
+		//Step 1(a) (b)
+		
+		beta = new double*[nh];
+		alpha = new double*[nh];
+		
+		k = model->max_indegree();
+		
+		cerr<<" . "<<nh<<" nodes"<<endl;
+		
+		struct timeval startTime;
+		struct timeval endTime;
+		struct timezone myTimeZone;
+		
+		double betaTime;
+		double alphaTime;
+		double RRFuncTime;
+		double HFuncTime;
+		double KFuncTime;
+		double GammaTime;
+		
+		double step2Time;
+		
+		
+		//Step 1(a)
+
+		gettimeofday(&startTime, &myTimeZone);
+		
+		//If Arguments::ADtree == 1, user wants to use ADtree, then make the whole tree first
+		if(Arguments::ADtree == 1){
+			//cout << "Arguments::ADtree == 1"  << endl;
+			model->makeADTree();
+		}
+		else{
+			//cout << "Arguments::ADtree == 0"  << endl;
+		}
+
+		for (int j = 0; j < nh; j ++){
+			beta[j] = new double[1 << nh];
+					
+			//int i = Vh[j];
+			compute_beta(j, Vh, nh, Vu, nu);
+
+			//cerr<<" . Tables beta and alpha computed for node "<<j<<"."<<endl; 
+		
+		}
+		
+		//If Arguments::ADtree == 1, then delete the whole tree
+		if(Arguments::ADtree == 1){
+			model->freeADTree();
+		}
+		
+		gettimeofday(&endTime, &myTimeZone);		
+		betaTime = (endTime.tv_sec - startTime.tv_sec) + 
+		           (endTime.tv_usec - startTime.tv_usec) / 1000000.0 ;
+				
+		cerr<<" . Tables beta have been computed" << endl;
+
+
+		//HR: Add for Top-k inside void compute_edge_probabilities_mixIndegree(int h)
+		//HR: Use the Poster tool to compute the local scores for all the families of each variable (without max-indegree
+		//	restriction). So for each of n variables, there are 2^(n-1) family scores.
+		write_family_scores_to_files(nh);
+		
+
+		//step 1(b)
+		gettimeofday(&startTime, &myTimeZone);
+		for (int j = 0; j < nh; j ++){
+
+			alpha[j] = new double[1 << nh];	
+			compute_alpha(j, Vh, nh);
+			
+			//cerr<<" . Tables beta and alpha computed for node "<<j<<"."<<endl; 
+		
+		}
+		gettimeofday(&endTime, &myTimeZone);
+		alphaTime = (endTime.tv_sec - startTime.tv_sec) + 
+		           (endTime.tv_usec - startTime.tv_usec) / 1000000.0 ;
+		
+		
+		cerr<<" . Tables alpha has been computed" << endl;
+		
+
+		//Step 1(c)
+		//HR: backward method to compute: for all S \subset V RR(S)
+		RRFunc = new double[1 << nh];
+		
+		cerr<<" . To compute Table RR" << endl;
+		
+		gettimeofday(&startTime, &myTimeZone);
+		
+		compute_RR_Fast(nh, RRFunc, alpha);
+		//compute_RR_Iter(nh, RRFunc, alpha);
+		
+		gettimeofday(&endTime, &myTimeZone);		
+		RRFuncTime = (endTime.tv_sec - startTime.tv_sec) + 
+		           (endTime.tv_usec - startTime.tv_usec) / 1000000.0 ;
+		
+		//cerr<<"RRFunc[all] = "<<RRFunc[(1<<nh)-1]<<endl;
+		//fprintf(stderr, "RRFunc[all] = %18.8f\n", RRFunc[(1<<nh)-1]);
+		cerr<<" . Table RR computed." <<endl;
+
+
+		//HR: merge the function of opt 4 into here.
+		//write to exactPD.txt
+		//string directory_name = "./family_scores";
+		string str;
+		str.assign(Arguments::directoryName);
+		str.append("/");
+		str.append("exactPD.txt");
+
+		//FILE * fp_PD = fopen(str.c_str(), "a");
+		FILE * fp_PD = fopen(str.c_str(), "w");
+		if(fp_PD){
+			fprintf(fp_PD, "%18.8f\n", RRFunc[(1<<nh)-1]);
+			fclose(fp_PD);
+
+		}
+		else{
+			cerr <<"fopen exactPD.txt fails" << endl;
+
+		}
+		
+
+		//Step 1(d)
+		//HR: forward method to compute: for all S \subset V H(S)
+		
+		hFunc = new double[1 << nh];
+		
+		cerr<<" . To compute Table H" << endl;
+		
+		gettimeofday(&startTime, &myTimeZone);
+
+		compute_h_Fast(nh, hFunc, alpha);	
+		//compute_h_Iter(nh, hFunc, alpha);	
+		
+		gettimeofday(&endTime, &myTimeZone);
+		HFuncTime = (endTime.tv_sec - startTime.tv_sec) + 
+		           (endTime.tv_usec - startTime.tv_usec) / 1000000.0 ;
+		
+		
+		//cerr<<"hFunc[all] = "<<hFunc[(1<<nh)-1]<<endl;
+		//fprintf(stderr, "hFunc[all] = %18.8f\n", hFunc[(1<<nh)-1]);
+		cerr<<" . Table H computed" <<endl;
+
+			
+		//Step 1 (e)
+		KFunc = new double*[nh]; 
+		cerr<<" . To compute Table K" << endl; 
+		
+		gettimeofday(&startTime, &myTimeZone);
+		for(int i = 0; i < nh; i++){
+			//cerr<<" For v = " << i << endl;
+			KFunc[i] = new double [1<<nh];
+			compute_all_K_v(nh, KFunc[i], i, alpha);
+			
+		}
+		
+		gettimeofday(&endTime, &myTimeZone);	
+		KFuncTime = (endTime.tv_sec - startTime.tv_sec) + 
+		           (endTime.tv_usec - startTime.tv_usec) / 1000000.0 ;
+		
+		           
+		cerr<<" . Table K computed" <<endl;
+		
+	
+		// Step 1 (f): For all i \in V, for all Pa_i \subset V - {i} with |Pa_i| <= k, compute Gamma_i(Pa_i)
+		
+		double * Gamma = new double[1 << nh];
+		double * Gfb = new double[1 << nh];
+
+		GammaTime = 0.0;
+		step2Time = 0.0;
+
+		for (int v = 0; v < nh; v++){
+			//cerr<<" For v = " << v <<": " << endl;
+			
+			//start counting the GammaTime for each v
+			gettimeofday(&startTime, &myTimeZone);
+			
+			// Compute Gfb[]. 
+			for (int U = 0; U < (1<<nh); U++){				
+				//if v \in the binary repre of U
+				if ((U>>v) & 1){
+					Gfb[U] = LOG_ZERO;
+				}
+				else{
+
+					//Gfb[U] = hFunc[U] + KFunc[v][U];
+					if(KFunc[v][U] == LOG_ZERO){
+						Gfb[U] = LOG_ZERO;
+					}
+					//note that KFunc[v][U] == 0 is possible
+					else if(KFunc[v][U] <= 0){
+						Gfb[U] = hFunc[U] + KFunc[v][U];
+					}
+					else{
+						//Now if Gfb[U] > 0, then Gfb[U]
+				
+						Gfb[U] = -(hFunc[U] - KFunc[v][U]);
+					}
+				
+					//cerr<<" j = "<<j<<": ";
+					//print_set(cerr, S);
+					//cerr<<"; ";
+					//print_set(cerr, cS);
+					//cerr<<": "<<gfb[S];
+					//cerr<<endl;
+				}
+			}
+						
+			//cerr<< "Table Gfb has been computed" << endl;
+					
+			// Compute gamma_j() := gamma[].
+						
+			fdmtLogAddComp(Gfb, Gamma, nh, k);
+			
+			gettimeofday(&endTime, &myTimeZone);
+			 //end counting GammaTime for each v
+			GammaTime += (endTime.tv_sec - startTime.tv_sec) + 
+		           (endTime.tv_usec - startTime.tv_usec) / 1000000.0 ;
+		           
+			
+			//cerr<< "Table Gamma has been computed" << endl;
+						
+			//cerr<< " Gfb[0]: " <<Gfb[0]<<endl;
+			//cerr<< " Gamma[0]: " << Gamma[0]<<endl;
+			
+			
+			//Step 2
+			//cerr<<" . Incoming edges for node "<< v <<":";
+			// Loop over start-point nodes i.
+			
+			//start counting step2Time for each v
+			gettimeofday(&startTime, &myTimeZone);
+
+			for (int i = 0; i < nh; i++){
+				if (i == v) continue;
+				//cerr<<" "<< i;		
+				double log_prob_nume = eval_edge_mixIndegree(i,  v, Gamma, beta[v], nh, k);
+		
+				model->print_edge_prob(
+					cout, Vh[i], Vh[v], exp(log_prob_nume - RRFunc[(1<<nh)-1])); //HR: - gb[(1<<nh)-1] mean /L(v)
+							
+			}
+			
+
+			gettimeofday(&endTime, &myTimeZone);
+			//end counting step2Time for each v	
+			step2Time += (endTime.tv_sec - startTime.tv_sec) + 
+		           (endTime.tv_usec - startTime.tv_usec) / 1000000.0 ;
+		    	
+			//cerr<<endl;
+			
+		} 
+		
+		
+//		cout << "\n\n\nTime Statistics" << endl;
+//		cout << "betaTime = " << betaTime <<  endl;				
+//		cout << "alphaTime = " << alphaTime <<  endl;
+//		cout << "RRFuncTime = " << RRFuncTime <<  endl;
+//		cout << "HFuncTime = " << HFuncTime <<  endl;
+//		cout << "KFuncTime = " << KFuncTime <<  endl;		
+//				
+//		cout << "GammaTime = " << GammaTime <<  endl;
+//		cout << "step2Time = " << step2Time <<  endl;
+//		double totalTime = betaTime + alphaTime + RRFuncTime 
+//						+ HFuncTime + KFuncTime + GammaTime + step2Time;
+//		cout << "totalTime = " << totalTime <<  endl;
+				
+		
+		delete [] Gamma;
+		delete [] Gfb;
+
+		delete [] hFunc;
+		
+		delete [] RRFunc; 
+		
+		delete [] Eta_U;
+		
+		for (int j = 0; j < nh; j ++){
+			delete [] beta[j]; 
+			delete [] alpha[j];
+
+			delete [] KFunc[j]; 
+		
+		}
+		delete [] beta; 
+		delete [] alpha;
+
+		delete [] KFunc;
+		
+		
+		cerr<<" Edge probabilities now computed."<<endl; 
+		
+		return;
+
+	}
+	//HR: end mixIndegree
+	
+	
+	void compute_edge_probabilities_mixIndegreeNoLog(int h){
+		cerr<<"Mix with max in-degree method: " << endl;
+		cerr<<" Compute edge probabilities for Layer "<<h<<":"<<endl;
+		
+		
+		// Step 1: Compute beta[][];
+		// Step 2: Compute alpha[][];
+		int *Vh, *Vu; int nh, nu;
+		model->layer(h, &Vh, &nh);
+		model->upper_layers(h-1, &Vu, &nu);
+		
+		//Step 1(a) (b)
+		
+		beta = new double*[nh];
+		alpha = new double*[nh];
+		
+		k = model->max_indegree();
+		
+		cerr<<" . "<<nh<<" nodes"<<endl;
+
+		for (int j = 0; j < nh; j ++){
+
+			beta[j] = new double[1 << nh];
+					
+			//int i = Vh[j];
+			compute_beta(j, Vh, nh, Vu, nu);
+
+			alpha[j] = new double[1 << nh];	
+			compute_alpha(j, Vh, nh);
+				
+			cerr<<" . Tables beta and alpha computed for node "<<j<<"."<<endl; 
+		
+		}
+		
+		
+		//Step 1(c)
+		//HR: backward method to compute: for all S \subset V RR(S)
+		RRFunc = new double[1 << nh];
+		
+		cerr<<" . To compute Table RR" << endl;
+		
+		compute_RR_Fast(nh, RRFunc, alpha);
+		
+		cerr<<" . Table RR computed." <<endl;
+		
+		
+		//Step 1(d)
+		//HR: forward method to compute: for all S \subset V H(S)
+		
+		hFunc = new double[1 << nh];
+		
+		cerr<<" . To compute Table H" << endl;
+		
+		compute_h_Fast(nh, hFunc, alpha);	
+		
+		cerr<<" . Table H computed" <<endl;
+
+			
+		//Step 1 (e)
+		KFunc = new double*[nh]; 
+		cerr<<" . To compute Table K" << endl; 
+		for(int i = 0; i < nh; i++){
+			cerr<<" For v = " << i << endl;
+			KFunc[i] = new double [1<<nh];
+			compute_all_K_v(nh, KFunc[i], i, alpha);
+			
+		}
+		cerr<<" . Table K computed" <<endl;		
+		
+		// Step 2: Loop over end-point nodes v \in V
+		
+		//No Log transformation
+		double * Gamma = new double[1 << nh];
+		double * Gfb = new double[1 << nh];
+		ofstream ofNoLog("debugEngineNoLog.txt");
+		
+		for (int v = 0; v < nh; v++){
+			cerr<<" For v = " << v <<": " << endl;
+			// Compute Gfb[]. 
+			for (int U = 0; U < (1<<nh); U++){				
+				//if v \in the binary repre of U
+				if ((U>>v) & 1){
+					//Gfb[U] = LOG_ZERO;
+					//Gfb[U] = 0;
+					Gfb[U] = LOG_ZERO;//
+				}
+				else{
+					//Gfb[U] = hFunc[U] + KFunc[v][U];
+					
+					if(KFunc[v][U] == LOG_ZERO){
+						Gfb[U] = exp(hFunc[U]);
+					}
+					else if (KFunc[v][U] > 0){
+						Gfb[U] = - exp(hFunc[U] - KFunc[v][U]);
+						
+					}
+					else{
+						Gfb[U] = exp(hFunc[U] + KFunc[v][U]);
+					}
+					//
+					if(Gfb[U] == 0){
+						Gfb[U] = LOG_ZERO;
+					}
+					else{
+						Gfb[U] = log(Gfb[U]);
+					}
+					//
+										
+					//cerr<<" j = "<<j<<": ";
+					//print_set(cerr, S);
+					//cerr<<"; ";
+					//print_set(cerr, cS);
+					//cerr<<": "<<gfb[S];
+					//cerr<<endl;
+				}
+			}
+			cerr<< "Table Gfb has been computed" << endl;
+					
+			ofNoLog <<"For v = " << v <<": " << endl;
+			ofNoLog << "\nGfb: " << endl;
+			for(int index = 0; index < (1 << nh); index++){
+				ofNoLog << "Gfb[" << index << " ] =" << Gfb[index] << endl;
+			}
+			ofNoLog << endl;
+			
+			
+			// Compute gamma_j() := gamma[].
+			//fdmtNoLog(Gfb, Gamma, nh, k);
+			fdmt(Gfb, Gamma, nh, k);
+			
+			cerr<< "Table Gamma has been computed" << endl;
+			
+			cerr<< " Gamma[0]: " << Gamma[0]<<endl;
+			cerr<< " Gfb[0]: " <<Gfb[0]<<endl;
+			
+			cerr<<" . Incoming edges for node "<< v <<":";
+			// Loop over start-point nodes i.
+			for (int i = 0; i < nh; i++){
+				if (i == v) continue;
+				cerr<<" "<< i;		
+				double log_prob_nume = eval_edge_mixIndegree(i,  v, Gamma, beta[v], nh, k);
+		
+				model->print_edge_prob(
+					cout, Vh[i], Vh[v], exp(log_prob_nume - RRFunc[(1<<nh)-1])); //HR: - gb[(1<<nh)-1] mean /L(v)
+			
+			}
+			cerr<<endl;
+			
+		} 
+	
+		delete [] Gamma;
+		delete [] Gfb;
+
+		delete [] hFunc;
+		
+		delete [] RRFunc; 
+		
+		delete [] Eta_U;
+		
+		for (int j = 0; j < nh; j ++){
+			delete [] beta[j]; 
+			delete [] alpha[j];
+
+			delete [] KFunc[j]; 
+		
+		}
+		delete [] beta; 
+		delete [] alpha;
+
+		delete [] KFunc;
+				
+		cerr<<" Edge probabilities now computed."<<endl; 
+		
+		return;
+
+	}
+	
+	//HR: end mixIndegreeNoLog
+	
+	
+	// beta[j][T] := sum_W gamma[j][T \cup W] .
+	void compute_beta(int j, int *Vh, int nh, int *Vu, int nu){
+		//cerr<<"Compute beta, j = "<<j<<", i = "<< Vh[j];
+		//cerr<<", k = "<<k<<", nh = "<<nh<<":"<<endl;
+		
+		sub_init(0, 0, beta[j], MARK, 0, nh);
+
+		//cerr<<" Initialization done."<<endl;
+		
+		int *S = new int[k];		
+		sub_beta(-1, nh, beta[j], 0, S, Vh[j], Vu, nu, 0);
+		delete [] S;
+				
+	}// end void compute_beta()
+	
+	
+	// beta[j][T] := sum_W gamma[j][T \cup W] .
+	void compute_betaNume(int j, int *Vh, int nh, int *Vu, int nu, int uu, int vv){
+		//cerr<<"Compute betaNume, j = "<<j<<", i = "<< Vh[j];
+		
+		sub_init(0, 0, betaNume[j], MARK, 0, nh);
+		
+		//cerr<<" Initialization done."<<endl;
+
+		int *S = new int[k];
+		
+		//check whether uu == j
+		if(uu != j){
+			sub_beta(-1, nh, betaNume[j], 0, S, Vh[j], Vu, nu, 0);
+		}
+		else{
+			//cerr << "uu == j == " << uu << ": so call sub_betaNume()" << endl;
+			sub_betaNume(-1, nh, betaNume[j], 0, S, Vh[j], Vu, nu, 0, uu, vv);
+		}
+		
+		
+		
+		delete [] S;
+	}
+	
+	
+	void compute_betaInOut(int j, int *Vh, int nh, int *Vu, int nu){
+
+		//cerr << "call: sub_init(0, 0, betaInOut[j], MARK, 0, nh)" << endl;
+		sub_init(0, 0, betaInOut[j], MARK, 0, nh);
+		//cerr << "return: sub_init(0, 0, betaInOut[j], MARK, 0, nh)" << endl;
+
+		//cerr<<" Initialization done."<<endl;
+
+		int *S = new int[k];
+		
+		//cerr << "call: sub_betaInOut()" << endl;
+		
+		sub_betaInOut(-1, nh, betaInOut[j], 0, S, Vh[j], Vu, nu, 0);
+		//sub_beta(-1, nh, betaInOut[j], 0, S, Vh[j], Vu, nu, 0);
+		
+		//cerr << "return: sub_betaInOut()" << endl;
+
+		delete [] S;
+	
+	}//end of compute_betaInOut()
+	
+	
+
+	void sub_beta(int jprev, int nh, double *b, int d, int *S, int i, int *Vu, int nu, int T){
+		//cerr<<" S:"; print_nodes(cerr, S, d); 
+		//cerr<<"; T:"; print_nodes(cerr, T, Vu, nu); cerr<<endl;
+		
+		
+		double lb;	
+		if(Arguments::option < 1){
+		
+			//HR: for old prior for q(G_i)and new Dirichlet hyper-param
+			//used to compare rebel and out method		
+
+			//lb = model->log_prior(i, S, d) + model->log_lcpHR(i, S, d);
+			if(Arguments::ADtree == 0){
+				lb = model->log_prior(i, S, d) + model->log_lcpHR(i, S, d);
+			}
+			else{
+				lb = model->log_prior(i, S, d) + model->log_lcpHR_ADtree(i, S, d);
+			}
+		}
+		else if(Arguments::option >= 1) {					
+
+			//lb = 0 + model->log_lcpHR(i, S, d);
+			if(Arguments::ADtree == 0){
+				lb = 0 + model->log_lcpHR(i, S, d);
+			}
+			else{
+				lb = 0 + model->log_lcpHR_ADtree(i, S, d);
+			}
+		}
+
+		//LOGADD(b[T], lb);
+		logAdd_New(b[T], lb);
+		//
+		
+		if (d < k){
+			for (int j = jprev + 1; j < nu; j ++){	
+				if (Vu[j] != i){
+					S[d] = Vu[j];
+					
+					if (j < nh) {
+						sub_beta(j, nh, b, d+1, S, i, Vu, nu, T | (1 << j));
+					}
+					
+					else{
+						sub_beta(j, nh, b, d+1, S, i, Vu, nu, T);
+					}
+				}
+			}
+		}	
+	}
+	
+	                
+	void sub_betaNume(int jprev, int nh, double *b, int d, int *S, int i, int *Vu, int nu, int T, int uu, int vv){
+		//cerr<<" S:"; print_nodes(cerr, S, d); 
+		//cerr<<"; T:"; print_nodes(cerr, T, Vu, nu); cerr<<endl;
+		int vvInS = 0;
+	
+		for(int ind = 0; ind < d; ind++){
+			if(vv == S[ind]){
+				vvInS = 1;
+				break;
+			}
+		}
+		
+		if(vvInS == 1){
+			double lb;
+			if(Arguments::option < 1){
+				
+				//HR: for old prior for q(G_i)and new Dirichlet hyper-param
+				//used to compare rebel and out method
+
+				//lb = model->log_prior(i, S, d) + model->log_lcpHR(i, S, d);
+				if(Arguments::ADtree == 0){
+					lb = model->log_prior(i, S, d) + model->log_lcpHR(i, S, d);
+				}
+				else{
+					lb = model->log_prior(i, S, d) + model->log_lcpHR_ADtree(i, S, d);
+				}
+			}
+			else if(Arguments::option >= 1) {					
+				//HR: for 1 prior for q(G_i) and new Dirichlet hyper-param
+				//lb = 0 + model->log_lcpHR(i, S, d);
+				if(Arguments::ADtree == 0){
+					lb = 0 + model->log_lcpHR(i, S, d);
+				}
+				else{
+					lb = 0 + model->log_lcpHR_ADtree(i, S, d);
+				}
+				
+			}
+		
+			//double lb = model->log_prior(i, S, d) + model->log_lcp(i, S, d);
+			//double lb = 0 + model->log_lcpHR(i, S, d);
+
+			//LOGADD(b[T], lb);
+			logAdd_New(b[T], lb);
+		}
+		else{
+			//how to set? 0-no, MARK-no; log(0) = -inf = LOG_ZERO
+			//LOGADD(b[T], LOG_ZERO);
+			double tempLogB = LOG_ZERO;
+			logAdd_New(b[T], tempLogB);
+		}
+		
+
+		
+		if (d < k){
+			for (int j = jprev + 1; j < nu; j ++){	
+				if (Vu[j] != i){
+					S[d] = Vu[j];
+					
+					if (j < nh) {
+						sub_betaNume(j, nh, b, d+1, S, i, Vu, nu, T | (1 << j), uu, vv);
+					}
+					
+					else{
+						sub_betaNume(j, nh, b, d+1, S, i, Vu, nu, T, uu, vv);
+					}
+				}
+			}
+		}		
+	}
+	
+	
+
+	void getInOutFeatures(int no_vars){
+	
+		for(int i = 0; i < no_vars; i++){
+			inFeatureEdges[i] = 0;
+		}
+		for(int i = 0; i < no_vars; i++){
+			outFeatureEdges[i] = 0;
+		}
+			
+		FILE * fpIn = fopen(Arguments::inFeasFileName, "r");
+	
+		int fromNode;
+		int toNode;
+		//char c;
+		//int i=0;
+		while(!feof(fpIn)){
+			
+			fromNode = -1;
+			toNode = -1;
+			//c = -1;
+			
+			fscanf(fpIn,"%d->%d\n",&fromNode, &toNode);
+			
+			//cout << "fromNode = "<< fromNode << endl;
+			//cout << "toNode = "<< toNode << endl;
+	
+			//To guard against toNode == -1
+			if(toNode == -1){
+				break;
+			}
+			else{
+				inFeatureEdges[toNode] += (1 << fromNode);
+			}
+
+		}
+		
+		fclose(fpIn);
+		
+//		for(int i = 0; i < no_vars; i++){
+//			cout << "inFeatureEdges[" << i << "] " << inFeatureEdges[i] << endl; 
+//		}
+
+		
+		FILE *fp2 = fopen(Arguments::outFeasFileName, "r");
+
+		while(!feof(fp2)){
+
+			fromNode = -1;
+			toNode = -1;
+			
+			fscanf(fp2,"%d->%d\n",&fromNode, &toNode);
+			
+			//cout << "fromNode = "<< fromNode << endl;
+			//cout << "toNode = "<< toNode << endl;
+	
+			if(toNode == -1){
+				break;
+			}
+			else{
+				outFeatureEdges[toNode] += (1 << fromNode);
+			}
+			
+		
+		}
+		fclose(fp2);
+	
+//		for(int i = 0; i < no_vars; i++){
+//			cout << "outFeatureEdges[" << i << "] " << outFeatureEdges[i] << endl; 
+//		}
+
+	} //end of getInOutFeatures
+	
+	
+	void sub_betaInOut(int jprev, int nh, double *b, int d, int *S, int i, int *Vu, int nu, int T){
+//		cerr<< "Start sub_betaInOut()" << endl;
+//		cerr<<" S:"; print_nodes(cerr, S, d); 
+//		cerr<<" T: " << T << endl;
+//		cerr<<" d: " << d << endl;
+		//cerr<<"; T:"; print_nodes(cerr, T, Vu, nu); cerr<<endl;
+				
+		int compBetaNeeded = 0;
+		
+		//cerr << "inFeatureEdges[" << i << "] = " << inFeatureEdges[i] << endl;
+		//cerr << "outFeatureEdges[" << i << "] = " << outFeatureEdges[i] << endl;
+		
+		//if (each of inFeatureEdges[i] is in T ) && (none of outFeatureEdges[i] is in T )
+		if(((inFeatureEdges[i] & T) == inFeatureEdges[i]) && ((outFeatureEdges[i] & T) == 0) ){
+			compBetaNeeded = 1;		
+			//cerr << "inFeatureEdges[" << i << "] = " << inFeatureEdges[i] << endl;
+			//cerr << "outFeatureEdges[" << i << "] = " << outFeatureEdges[i] << endl;
+			//cerr << "T = " << T << endl;
+			//cerr<<" compBetaNeeded = 1;" << endl;
+			
+		}
+		
+		if(compBetaNeeded == 1){
+			double lb;
+			//Arguments::option must be 5
+			if(Arguments::option < 1){
+				
+				//HR: for old prior for q(G_i)and new Dirichlet hyper-param
+				//used to compare rebel and out method
+				//lb = model->log_prior(i, S, d) + model->log_lcpHR(i, S, d);
+				if(Arguments::ADtree == 0){
+					lb = model->log_prior(i, S, d) + model->log_lcpHR(i, S, d);
+				}
+				else{
+					lb = model->log_prior(i, S, d) + model->log_lcpHR_ADtree(i, S, d);
+				}
+			}
+			else if(Arguments::option >= 1) {					
+				//HR: for 1 prior for q(G_i) and new Dirichlet hyper-param
+				//cerr << "start: double lb = 0 + model->log_lcpHR(i, S, d)" << endl;
+				//lb = 0 + model->log_lcpHR(i, S, d);
+				if(Arguments::ADtree == 0){
+					lb = 0 + model->log_lcpHR(i, S, d);
+				}
+				else{
+					lb = 0 + model->log_lcpHR_ADtree(i, S, d);
+				}
+				//cerr << "end: double lb = 0 + model->log_lcpHR(i, S, d)" << endl;
+			}
+			
+		
+			//LOGADD(b[T], lb);
+			logAdd_New(b[T], lb);
+		}
+		else{
+			//LOGADD(b[T], LOG_ZERO);
+			double tempLogB = LOG_ZERO;
+			logAdd_New(b[T], tempLogB);
+		}
+
+		if (d < k){
+			for (int j = jprev + 1; j < nu; j ++){	
+				if (Vu[j] != i){
+					S[d] = Vu[j];
+					
+					if (j < nh) {						
+						sub_betaInOut(j, nh, b, d+1, S, i, Vu, nu, T | (1 << j));
+					}
+					
+					else{
+						sub_betaInOut(j, nh, b, d+1, S, i, Vu, nu, T);
+					}
+				}
+			}
+		}		
+	} //end of sub_betaInOut
+	
+	
+	void sub_init(int d, int S, double *a, double value, int ones, int nh){
+		//cerr<<" d: "<<d<<endl; 
+		a[S] = value;
+		
+		if (d < nh){
+			sub_init(d+1, S, a, value, ones, nh);
+
+			if (ones < k) sub_init(d+1, S | (1 << d), a, value, ones+1, nh);
+		}
+	}
+
+	//HR: From get_kbest_nets.cc in Top-k
+	//HR: add v itself into the binary repres, v posi should be 0
+	//Convert the parent set to the appropriate variable set
+	unsigned int parset2varset(int v, unsigned int set) {
+		unsigned int sinkleton = 1U << v;
+		unsigned int mask = sinkleton - 1;
+		unsigned int low = set & mask;
+		unsigned int high = set & ~mask;
+		return (high << 1) | low;
+	}
+
+
+	//HR: Add for Top-k inside void compute_edge_probabilities_mixIndegree(int h)
+	//HR: Use the Poster tool to compute the local scores for all the families of each variable (without max-indegree
+	//	restriction). So for each of n variables, there are 2^(n-1) family scores.
+	//	Then store them as file 0, 1, 2, ..., n-1 for variable 0, 1, ..., n-1.
+	void write_family_scores_to_files(int nh){
+		FILE* fout;
+		//nh is the number of variables
+		int num_of_vars = nh;
+		unsigned int num_of_parsets = 1U << (num_of_vars - 1);
+		//HR: for each variable
+		for (int v = 0; v < num_of_vars; v++) {
+			char str1[40];
+			sprintf(str1, "%d", v);
+			string str;
+			//string dirname = "./family_scores";
+			string dirname = Arguments::directoryName;
+			str.append(dirname);
+			str.append("/");
+			str.append(str1);
+
+			fout = fopen(str.c_str(), "wb");
+			unsigned int new_par_set;
+			for (unsigned int par_set = 0; par_set < num_of_parsets; par_set++) {
+				//HR: add v itself into the binary repres, v posi should be 0
+				new_par_set = parset2varset(v, par_set);
+				//This is fine.
+				//fwrite(beta[v] + new_par_set, sizeof(double), 1, fout);
+				//This is also fine.
+				fwrite(&(beta[v][new_par_set]), sizeof(double), 1, fout);
+
+			}
+			fclose(fout);
+
+		}
+//		//test
+//		//HR: for each variable
+//		FILE* fin;
+//		for (int v = 0; v < num_of_vars; v++) {
+//			char str1[40];
+//			sprintf(str1, "%d", v);
+//			string str;
+//			//string dirname = "./family_scores";
+//			string dirname = Arguments::directoryName;
+//			str.append(dirname);
+//			str.append("/");
+//			str.append(str1);
+//
+//			fin = fopen(str.c_str(), "rb");
+//			unsigned int new_par_set;
+//			for (unsigned int par_set = 0; par_set < num_of_parsets; par_set++) {
+//				//HR: add v itself into the binary repres, v posi should be 0
+//				new_par_set = parset2varset(v, par_set);
+//				double temp_family_score;
+//				fread(&temp_family_score, sizeof(double), 1, fin);
+//				if(temp_family_score != beta[v][new_par_set]){
+//					cerr << "temp_family_score = " << temp_family_score << endl;
+//					cerr << "beta[" << v << "][" << new_par_set << "] = " << beta[v][new_par_set] << endl;
+//					cerr << "Error: temp_family_score != beta["
+//							<< v << "][" << new_par_set << "]" << endl;
+//					exit(1);
+//				}
+//			}
+//			fclose(fin);
+//
+//		}
+//		//
+
+	}//write_family_scores_to_files(int nh)
+
+	
+	// alpha[j][S] = sum_T beta[j][T]  (...times the prior).
+	void compute_alpha(int j, int *Vh, int nh){
+		// Use fast truncated upward Möbius transform.
+		fumt(beta[j], alpha[j], nh, k); 
+	}
+	
+
+	void compute_alphaNume(int j, int *Vh, int nh){
+		// Use fast truncated upward Möbius transform.
+		fumt(betaNume[j], alphaNume[j], nh, k); 
+	}
+	
+
+	void compute_alphaInOut(int j, int *Vh, int nh){
+		// Use fast truncated upward Möbius transform.
+		fumt(betaInOut[j], alphaInOut[j], nh, k); 
+	}
+	
+	
+	void compute_g_forward(int nh){	
+		sub_gf(0, nh, 0);
+	
+	}
+	
+	
+	void compute_gNume_forward(int nh){	
+		sub_gfNume(0, nh, 0);
+	
+	}
+	
+
+	void sub_gf(int d, int nh, int S){
+		
+		if (d < nh){
+			sub_gf(d+1, nh, S);
+			sub_gf(d+1, nh, S | (1 << d)); 
+		}
+		else {
+			double sum = MARK; 
+			int T = S, J = 1;
+			
+			for (int j = 0; j < nh; j ++){
+				if (T & 1){
+					// Now S - J is a subset of S.
+					
+					double w = alpha[j][S - J] + gf[S - J];
+					//LOGADD(sum, w);
+					logAdd_New(sum, w);
+				} 
+					
+				T >>= 1; J <<= 1;
+			}
+			if (S == 0) 
+				gf[0] = 0; 
+			else 
+				gf[S] = sum;
+		}
+	}
+	
+	
+		void sub_gfNume(int d, int nh, int S){
+		
+		if (d < nh){
+			sub_gfNume(d+1, nh, S);
+			sub_gfNume(d+1, nh, S | (1 << d)); 
+		}
+		else {
+			double sum = MARK; 
+			int T = S, J = 1;
+			
+			for (int j = 0; j < nh; j ++){				
+				if (T & 1){
+					// Now S - J is a subset of S.
+					
+					double w = alphaNume[j][S - J] + gfNume[S - J];
+					//LOGADD(sum, w);
+					logAdd_New(sum, w);
+				} 
+	
+				T >>= 1; J <<= 1;
+			}
+			if (S == 0) 
+				gfNume[0] = 0; 
+			else 
+				gfNume[S] = sum;
+		}
+	}
+	
+	
+	
+	void compute_g_backward(int nh){
+		sub_gb(0, nh, 0);
+	
+	}
+	
+	
+	void sub_gb(int d, int nh, int S){
+		if (d < nh){
+			sub_gb(d+1, nh, S);
+			sub_gb(d+1, nh, S | (1 << d)); 
+		}
+		else {
+			double sum = MARK; 
+			int T = S, J = 1, complS = (1 << nh) - 1 - S;
+			for (int j = 0; j < nh; j ++){
+				if (T & 1){
+					// Now S - J is a subset of S.
+					double w = alpha[j][complS] + gb[S - J];
+					//LOGADD(sum, w);
+					logAdd_New(sum, w);
+				} 		
+				T >>= 1; J <<= 1;
+			}
+			if (S == 0) gb[0] = 0; else gb[S] = sum;
+		}
+	}
+	
+
+	double eval_edge(int i, int j, double *a, double *b, int nh, int k){
+		int T = 1 << i;  
+		return sub_eval_edge(-1, 1, T, i, j, a, b, nh, k);
+	}
+	
+
+	double sub_eval_edge(
+		int tprev, int d, int T, int i, int j, double *a, double *b, int nh, int k){
+		
+		//cerr<<" T:"; print_set(cerr, T); 
+		//cerr<<": "<<a[T]<<"; "<<b[T]<<endl;
+		
+		//if T, d = |G_v| = k	
+		if (d == k){
+			return a[T] + b[T];
+		}
+		
+		//if  T < k, can add more
+		double sum = a[T] + b[T];
+		for (int t = tprev + 1; t < nh; t ++){
+			if (t == i || t == j) continue;
+			
+			int Tnext = T | (1 << t);
+			double w = sub_eval_edge(t, d+1, Tnext, i, j, a, b, nh, k);
+			//LOGADD(sum, w);
+			logAdd_New(sum, w);
+		}
+		return sum;
+	}
+	
+	
+	double eval_edge_mixIndegree(int i, int j, double * Gamma, double * Beta_j, int nh, int k){
+		int T = 1 << i;  
+		return sub_eval_edge_mixIndegree(-1, 1, T, i, j, Gamma, Beta_j, nh, k);
+	}
+	
+	
+	double sub_eval_edge_mixIndegree(
+		int tprev, int d, int T, int i, int j, double * Gamma, double * Beta_j, int nh, int k){
+		
+		//cerr<<" T:"; print_set(cerr, T); 
+		//cerr<<": "<<a[T]<<"; "<<b[T]<<endl;
+		
+		//if T, d = |G_v| = k	
+		if (d == k){
+			double sumRes;
+			if(Beta_j[T] == MARK){
+				cerr<<"\n**** Beta_j[T] == MARK ";
+		   	}
+			
+			//return (Gamma[T] + Beta_j[T]);
+			if(Gamma[T] == LOG_ZERO){
+				sumRes = LOG_ZERO;
+			}
+					
+			else if(Gamma[T] <= 0){
+				sumRes = Beta_j[T] + Gamma[T];
+			}
+			else{
+				cerr<<"\n**** Gamma[T] > 0 ";
+				sumRes = - ( Beta_j[T] - Gamma[T] );
+
+			}
+			return sumRes;
+	
+		}
+		
+		//if  T < k, can add more
+		//double sum = Gamma[T] + Beta_j[T];
+		double sum;
+		if(Beta_j[T] == MARK){
+			cerr<<"\n**** Beta_j[T] == MARK ";
+		}
+		if(Gamma[T] == LOG_ZERO){
+				sum = LOG_ZERO;
+		}
+					
+		else if(Gamma[T] <= 0){
+			sum = Beta_j[T] + Gamma[T];
+		}
+		else{
+			cerr<<"\n**** Gamma[T] > 0 ";
+			sum = - ( Beta_j[T] - Gamma[T] );
+
+		}
+		
+	
+		for (int t = tprev + 1; t < nh; t ++){
+			if (t == i || t == j) continue;
+
+			int Tnext = T | (1 << t);
+			double w = sub_eval_edge_mixIndegree(t, d+1, Tnext, i, j, Gamma, Beta_j, nh, k);
+			logAddComp(sum, w);
+		}
+		return sum;
+	}
+	
+	
+	double eval_edge_mixIndegreeNoLog(int i, int j, double * Gamma, double * Beta_j, int nh, int k){
+		int T = 1 << i;  
+		return sub_eval_edge_mixIndegreeNoLog(-1, 1, T, i, j, Gamma, Beta_j, nh, k);
+	}
+	
+
+	double sub_eval_edge_mixIndegreeNoLog(
+		int tprev, int d, int T, int i, int j, double * Gamma, double * Beta_j, int nh, int k){		
+		//cerr<<" T:"; print_set(cerr, T); 
+		//cerr<<": "<<a[T]<<"; "<<b[T]<<endl;
+		
+		//if T, d = |G_v| = k	
+		if (d == k){
+			//Lose precision
+			//return (Gamma[T] * exp(Beta_j[T]));
+			double product_res;
+			//if Beta_j[T]  has not been initialized
+			if(Beta_j[T] == MARK){
+				cerr<<"***Beta_j[T]  has not been initialized" << endl;
+				product_res = 0;		
+			}
+			else{
+				if(Gamma[T] == 0){
+				 	product_res = exp(Beta_j[T]);
+				}
+				else if(Gamma[T] > 0){
+					product_res = exp(Beta_j[T] + log(Gamma[T]));
+				}
+				else{
+					product_res = - exp(Beta_j[T] + log(-Gamma[T]));
+				}
+				
+			}
+
+			return product_res;
+			
+		}
+		
+		//if  T < k, can add more
+		//double sum = Gamma[T] * exp(Beta_j[T]);
+		double sum;
+		double product_res;
+		
+			//if Beta_j[T]  has not been initialized
+			if(Beta_j[T] == MARK){
+				cerr<<"***Beta_j[T]  has not been initialized" << endl;
+				product_res = 0;		
+			}
+			else{
+				if(Gamma[T] == 0){
+				 	product_res = exp(Beta_j[T]);
+				}
+				else if(Gamma[T] > 0){
+					product_res = exp(Beta_j[T] + log(Gamma[T]));
+				}
+				else{
+					product_res = - exp(Beta_j[T] + log(-Gamma[T]));
+				}
+				
+			}
+			
+		sum = product_res;
+		for (int t = tprev + 1; t < nh; t ++){
+			if (t == i || t == j) continue;
+
+			int Tnext = T | (1 << t);
+			double w = sub_eval_edge_mixIndegreeNoLog(t, d+1, Tnext, i, j, Gamma, Beta_j, nh, k);
+			//LOGADD(sum, w);
+			sum += w;
+		}
+		return sum;
+	}
+	
+	
+	Model *model;
+	
+	double **alpha, **beta;
+	double **alphaNume, **betaNume;
+	
+	double *gf, *gb;
+	
+	double *gfNume;
+	int k;
+	
+	double **alphaInOut, **betaInOut;
+	int * inFeatureEdges; 
+	int * outFeatureEdges; 
+	
+	
+};
+
+
+
+#endif