Add files via upload

1 files changed, 1652 insertions, 0 deletions

diff --git a/sourcecodes/k-best/src/UpdateHR.h b/sourcecodes/k-best/src/UpdateHR.h
new file mode 100644
index 00000000..19b93c1d
--- /dev/null
+++ b/sourcecodes/k-best/src/UpdateHR.h
@@ -0,0 +1,1652 @@
+#include<stdio.h>
+#include<stdlib.h>
+#include"math.h"
+
+#include"Arguments.h"
+#include"Model.h"
+#include"Engine.h"
+
+
+#define LOGADD_NEW(logA, logB) if(logA == MARK) logA = logB; else if (logA == LOG_ZERO) logA = logB; else if (logB == LOG_ZERO) logA = logA; else logA = logA + log(1+exp(logB-logA));
+    	
+//#define LOGMINUS(la,lb) if(la>8e99||lb-la>100) la=lb;else la+=log(1-exp(lb-la));
+
+#define LOGMINUS_NEW(logA, logB) if(logA == MARK) logA = logB; else if (logA == LOG_ZERO) logA = logB; else if (logB == LOG_ZERO) logA = logA; else logA = logA + log(1-exp(logB-logA));
+
+#define epsilon 1.0e5
+
+bool logSpecValEquals(double logX, double logSpecVal){
+	return (logX >= logSpecVal - epsilon && logX <= logSpecVal + epsilon);
+}
+
+//call by reference to the original value
+//the log sum result is stored in logA
+//replace LOGADD_NEW(logA, logB)
+//Pre: if neither logA nor logB is MARK/LOG_ZERO, then logA < 0 && logB < 0
+void logAdd_New(double & logA, double & logB){
+	    //if(logA == MARK){
+	    if(logSpecValEquals(logA, MARK)){
+    		logA = logB;
+    	}
+    	//else if (logA == LOG_ZERO){
+    	else if (logSpecValEquals(logA, LOG_ZERO)){
+    		logA = logB;
+    	}
+    	//else if(logB == MARK){
+    	else if(logSpecValEquals(logB, MARK)){
+    		logA = logA;
+    	}
+    	//else if(logB == LOG_ZERO){
+    	else if(logSpecValEquals(logB, LOG_ZERO)){
+    		logA = logA;
+    	}
+    	//Nither logA nor logB is MARK/LOG_ZERO
+    	//logA < 0 && logB < 0
+    	else{
+    		//1.1
+			//Gaurd logB-logA <= 100 to avoid that exp(logB-logA) go to inf
+    		if(logB - logA <= 100 ){
+    			logA = logA + log( 1 + exp(logB-logA) );
+
+    		} 
+    		//1.2		
+    		//if logB - logA > 100, then logB > 100 + logA so that B is much bigger than A (B > e^100 * A)
+    		//    So (A + B) almost = B 
+    		else{
+    			//logA = logB; 
+	    		//The following is the real formula,; but it is almost same as //logA = logB; 
+	    		logA = logB + log( exp(logA - logB) + 1);
+    		}
+    		
+    	}
+}// end of void logAdd_New(double & logA, double & logB)
+
+//Pre:
+//  (logA != MARK) and (if logB != LOG_ZERO then logA !=LOG_ZERO) and (logA < 0); (i.e. logA has been assigned some value)
+//  if logB != MARK/LOG_ZERO, then logA >= logB
+
+void logMinus_New(double & logA, double & logB){
+	    //if(logA == MARK){
+	    if(logSpecValEquals(logA, MARK)){
+    		cerr << "Error: logMinus_New(double & logA, double & logB): logA == MARK\n"; 
+    		exit(1);
+    	}
+    	//else if (logA == LOG_ZERO && logB != LOG_ZERO){
+    	else if (logSpecValEquals(logA, LOG_ZERO) && ! logSpecValEquals(logB, LOG_ZERO)){
+    		cerr << "Error: logMinus_New(double & logA, double & logB): logA == LOG_ZERO && logB != LOG_ZERO\n"; 
+    		exit(1);
+    	}
+    	//else if(logB == MARK){
+    	else if(logSpecValEquals(logB, MARK)){
+    		logA = logA;
+    	}
+    	//else if(logB == LOG_ZERO){
+    	else if(logSpecValEquals(logB, LOG_ZERO)){
+    		logA = logA;
+    	}
+    	//Nither logA nor logB is MARK/LOG_ZERO
+    	//logA >= logB
+    	else{
+    		if(logA == logB){
+    			logA = LOG_ZERO;
+    		}
+    		//logA > logB
+    		else{
+    			if(logA > logB){
+    				logA = logA + log(1 - exp(logB-logA));
+    			}
+    			//if logA < logB
+    			//Note: this is the case which only occurs in 
+    			//	sub_h_Fast_Details() and sub_RR_Fast_Details() for the testing; it should not happen generally
+    			else{
+    				cerr << "Warning: logMinus_New(double & logA, double & logB): logA < logB" << endl;
+    				logA = -( logB + log(1 - exp(logA-logB)) );
+    			}
+
+    			
+    		}
+    		
+    	}
+} // end of void logMinus_New(double & logA, double & logB)
+
+
+//Add: modified from void logMinus_New(double & logA, double & logB)
+
+//Only for void sub_RR_Fast_TSize() and void sub_h_Fast_TSize()
+//Pre:
+//  (logA != MARK) and 
+//  (logB != MARK) and 
+//   both logA >= logB and logA < logB are allowed
+// Post:
+//	if logA < logB, then return LOG_ZERO bound
+//  ow, return the correctly exact answer
+void logMinus_New_NonnegRes(double & logA, double & logB){
+	    //if(logA == MARK){
+	    if(logSpecValEquals(logA, MARK)){
+    		cerr << "Error: logMinus_New_NonnegRes(double & logA, double & logB): logA == MARK\n"; 
+    		exit(1);
+    	}
+    	//else if(logB == MARK){
+    	else if(logSpecValEquals(logB, MARK)){
+    		cerr << "Error: logMinus_New_NonnegRes(double & logA, double & logB): logB == MARK\n"; 
+    		exit(1);
+    	}
+    	//else if(logB == LOG_ZERO){
+    	else if(logSpecValEquals(logB, LOG_ZERO)){
+    		logA = logA;
+    	}
+    	//Nither logA nor logB is MARK/LOG_ZERO
+    	//logA >= logB
+    	else{
+    		if(logA == logB){
+    			logA = LOG_ZERO;
+    		}
+    		//logA > logB
+    		else{
+    			if(logA > logB){
+    				logA = logA + log(1 - exp(logB-logA));
+    			}
+    			//if logA < logB, return LOG_ZERO as the bound
+    			//Note: this is the case which only occurs in 
+    			//	sub_h_Fast_TSize() and sub_RR_Fast_TSize();
+    			else{
+    				logA = LOG_ZERO;
+    			}
+	
+    		}
+    		
+    	}
+} // end of void logMinus_New_NonnegRes(double & logA, double & logB)
+
+
+//call by reference to the original value
+//the log sum result is stored in logA
+void logAddComp(double & logA, double & logB){
+
+    	//if(logA == MARK){
+	    if(logSpecValEquals(logA, MARK)){
+    		logA = logB;
+    	}
+    	//else if (logA == LOG_ZERO){
+    	else if (logSpecValEquals(logA, LOG_ZERO)){
+    		logA = logB;
+    	}
+    	//else if(logB == MARK){
+    	else if(logSpecValEquals(logB, MARK)){
+    		logA = logA;
+    	}
+    	//else if(logB == LOG_ZERO){
+    	else if(logSpecValEquals(logB, LOG_ZERO)){
+    		logA = logA;
+    	}
+    	
+    	
+    	//Nither logA nor logB is MARK/LOG_ZERO
+    	else{
+    		//1
+    		if(logA < 0 && logB < 0 ){
+    			//1.1
+				//Gaurd logB-logA <= 100 to avoid that exp(logB-logA) go to inf
+	    		if(logB - logA <= 100 ){
+	    			logA = logA + log( 1 + exp(logB-logA) );
+	
+	    		} 
+	    		//1.2		
+	    		//if logB - logA > 100, then logB > 100 + logA so that B is much bigger than A (B > e^100 * A)
+	    		//    So (A + B) almost = B 
+	    		else{
+	    			//logA = logB; 
+		    		//The following is the real formula,; but it is almost same as //logA = logB; 
+		    		logA = logB + log( exp(logA - logB) + 1);
+	    		}
+    		}
+    		
+    		//2
+    		else if(logA > 0 && logB > 0){
+    			//always make exp(.): . > 0
+    			if(-logB + logA <= 100){
+    				//2.1
+    				logA = -( -logA + log( 1 + exp(-logB+logA) ) );
+    			} 			
+    			else{
+    				//2.2
+    				logA = -( -logB + log( 1 + exp(-logA+logB) ) );
+    			}
+    			
+    		}
+    		//3
+    		else if(logA < 0 && logB > 0){
+    			//3.1
+    			if(logA + logB > 0){   	
+    				if(logA + logB <= 100){
+    					//3.1.2
+    					logA = -logB + log( exp(logA + logB) - 1 );
+    				}	
+    				else{
+    					//3.1.1
+    					logA = logA + log( 1 - exp(-logB-logA) );
+    				}	
+	
+    			}
+    			//3.2
+    			else if(logA + logB < 0){
+    				if(-logA - logB <= 100){
+    					//3.2.2
+    					logA = -( logA + log( exp(-logB-logA)  - 1) );
+    				}
+    				else{
+    					//3.2.1
+    					logA = - ( -logB + log(1 - exp(logA+logB) ) );
+    				}
+
+    			}
+    			//3.3
+    			else{
+    				logA = LOG_ZERO;
+    			}
+    		}
+    		//4
+    		else if(logA > 0 && logB < 0){
+    			//4.1
+    			if(logA + logB > 0){
+    				if(logA + logB <= 100){
+    					//4.1.2
+    					logA = -logA + log( exp(logB+logA) - 1);
+    				}
+    				else{
+    					//4.1.1
+    					logA = logB + log( 1 - exp(-logA-logB) );
+    				}
+    				
+    			}
+    			//4.2
+    			else if(logA + logB < 0){
+    				if(-logA-logB <= 100){
+    					//4.2.2
+    					logA = -( logB + log( exp(-logA-logB) - 1) );
+    				}
+    				else{
+    					//4.2.1
+    					logA = -( -logA + log(1 - exp(logB+logA)) );
+    				}
+	
+    			}
+    			//4.3
+    			else{
+    				logA = LOG_ZERO;
+    			}			
+    		}
+		
+    	}
+	
+} // end of void logAddComp(double & logA, double & logB)
+
+    	
+    	
+//#include <vector>
+
+//forward
+double * hFunc;
+double * hFuncNume;
+
+
+//backward
+//AA_S = new double [S+1];
+double * AA_S;
+
+double * RRFunc;
+double * RRFuncNume;
+
+
+//mixIndegree
+double * Eta_U;
+double ** KFunc;
+
+//In_Out
+double * RRFuncInOut;
+
+//ofstream totalOfMixEta("debugEngineMixForEta_U.txt");
+
+void printIntArr(int * TSet, int leng){
+		for(int ind = 0; ind < leng; ind++){
+			cout << "Set[" << ind << "]=" << TSet[ind] << endl;
+		}	
+}
+
+
+int binIntLen(int S){
+	int SSetLen = 0;
+
+	while(S > 0){ 
+		if(S & 1){
+			SSetLen++;
+			
+		}
+		S >>= 1;
+	}
+	return SSetLen;
+	
+}
+
+
+void binIntToDecIntArr(int S, int ones, int *SSet){
+	int posFromLeft = 0;
+	int ind = 0;
+	while(S > 0){ 
+		if(S & 1){
+			SSet[ind++] = posFromLeft;
+			
+		}
+		S >>= 1;
+		posFromLeft++;
+	}
+	//printIntArr(SSet, ones);
+	
+}
+
+
+
+//HR: Compute F(S, T)
+double compute_f(int nh, int S, int T, int sizeofT, int *TSet, double * h, double ** alpha){
+	
+
+	double result_f = MARK;
+	double result_h = h[S-T];
+	//LOGADD(result_f, result_h);
+	result_f = result_h;
+	if (result_f == LOG_ZERO){
+		return result_f;
+	}
+
+	
+	for(int index = 0; index < sizeofT; index++){
+		int j = TSet[index];
+		double result_A = alpha[j][S-T];
+		//result_f += result_A;
+		if (result_A == LOG_ZERO){
+			return LOG_ZERO;
+		}
+		else{
+			result_f += result_A;
+		}
+	}
+	return result_f;
+}
+
+
+//HR: Compute: \sum_{k=1}^{|S|} (-1)^{k+1} \sum{T \subset S, |T|=k} F(S,T)
+void rec_compute_sum_F_S_T(
+	int d,
+	int nh, 
+	int S, 
+	int ones, 
+	int * SSet, 
+	int T,
+	int sizeofT,
+	int * TSet,
+	double * h, 
+	double ** alpha, 
+	double & sumOdd, 
+	double & sumEven)
+{
+
+	if(d < nh){	
+		rec_compute_sum_F_S_T(d+1, nh, S, ones, SSet, T, sizeofT, TSet, h, alpha, sumOdd, sumEven);
+		
+		//the binary repre of S has 1 in d'th position (starting from 0) from right
+		if(((S>>d) & 1) == 1){
+			TSet[sizeofT] = d;
+			rec_compute_sum_F_S_T(d+1, nh, S, ones, SSet, (T | (1<<d)), (sizeofT + 1), TSet, h, alpha, sumOdd, sumEven);
+		}
+	}
+	else{
+		if(sizeofT == 0){
+			//|T| == 0 does not count
+			return;
+			
+		}
+		else{	
+			double result_F = compute_f(nh, S, T, sizeofT, TSet, h, alpha);
+			//sizeofT is odd
+			if(sizeofT % 2 == 1){				
+				logAdd_New(sumOdd, result_F);
+				//LOGADD_NEW(sumOdd, result_F);
+			}
+			else{
+				logAdd_New(sumEven, result_F);
+				//LOGADD_NEW(sumEven, result_F);
+			}
+			
+		}
+		
+	}
+} // end rec_compute_sum_F_S_T()
+
+	
+//HR: Compute: \sum_{k=1}^{|S|} (-1)^{k+1} \sum{T \subset S, |T|=k} F(S,T)
+void compute_sum_F_S_T(int nh, int S, int ones, int * SSet, double * h, double ** alpha, double & sumOdd, double & sumEven){
+	//TSet's length is set to be SSet's length: it is fine since each TSetLen is set accordingly to each subset T;
+	int * TSet = new int[ones];
+	rec_compute_sum_F_S_T(0, nh, S, ones, SSet, 0, 0, TSet, h, alpha, sumOdd, sumEven);
+	delete[] TSet;
+	
+} // end compute_sum_F_S_T()
+
+
+void sub_h_Fast(int d, int nh, int S, int ones, int * SSet, double * h, double ** alpha){
+	if(d < nh){
+		
+		//this kind of ordering in the recursion is crucial; 
+		//it guarantee when you want to compute one S, all its subsets have already been computed
+		sub_h_Fast(d+1, nh, S, ones, SSet, h, alpha);
+		
+		SSet[ones] = d;
+		sub_h_Fast(d+1, nh, S|(1<<d), (ones+1), SSet, h, alpha);
+	}
+	//d == nh
+	else{
+		//S = 0 means that S is the empty set; 
+		if(S == 0){
+			h[0] = 0; //log(1) = 0;
+			//cout << "h[0] = 0" << endl;
+			//cout << endl;
+			return;
+
+		}
+		
+		//if S > 0
+		double sumOdd = MARK;
+		double sumEven = MARK;
+		
+		
+		compute_sum_F_S_T(nh, S, ones, SSet, h, alpha, sumOdd, sumEven);
+		
+		//if sumEven has been assigned some value instead of MARK
+		if(sumEven != MARK){
+			logMinus_New(sumOdd, sumEven);
+			//LOGMINUS_NEW(sumOdd, sumEven);
+		}
+						
+		h[S] = sumOdd;
+		//cout << "h[" << S << "] = "  << h[S] << endl;
+		//cout << endl;
+
+	}
+} // end sub_h_Fast()
+
+
+//HR:
+//	Compute h
+void compute_h_Fast(int nh, double * h, double ** alpha){
+	int * SSet = new int[nh];
+	
+	sub_h_Fast(0, nh, 0, 0, SSet, h, alpha);
+	
+	delete [] SSet;
+	
+} // end compute_h_Fast()
+
+
+
+
+
+//HR: Compute: \sum_{k=1}^{|S|} (-1)^{k+1} \sum{T \subset S, |T|=k} F(S,T)
+void compute_sum_F_S_T_Iter(int nh, int S, double * h, double ** alpha, double & sumOdd, double & sumEven){
+	
+	//all the subsets T of S, expect that T = 0 (empty set)
+	for(int T = 1; T <= S; T++){
+		//if T is the subset of S
+		if((S|T) == S){
+			
+			int TSetLen = binIntLen(T);
+			int * TSet = new int[TSetLen];
+			binIntToDecIntArr(T, TSetLen, TSet);
+			
+			double result_F = compute_f(nh, S, T, TSetLen, TSet, h, alpha);
+			//sizeofT is odd
+			if(TSetLen % 2 == 1){				
+				logAdd_New(sumOdd, result_F);
+				//LOGADD_NEW(sumOdd, result_F);
+			}
+			else{
+				logAdd_New(sumEven, result_F);
+				//LOGADD_NEW(sumEven, result_F);
+			}
+			
+			delete[] TSet;
+	
+		}
+	
+	}
+	
+}
+
+
+//HR:
+//	Compute h by iteration instead of recursion
+void compute_h_Iter(int nh, double * h, double ** alpha){
+	
+	//
+	h[0] = 0;
+
+	//h[S]: S: from 1 to (1<<nh)-1
+	for(int S = 1; S < (1<<nh); S++){
+		
+		double sumOdd = MARK;
+		double sumEven = MARK;
+		
+		//int ones = binIntLen(S);
+		//cout << "\nones = " << ones << endl;
+		//int * SSet = new int[ones];
+		//binIntToDecIntArr(S, ones, SSet);
+		//compute_sum_F_S_T(nh, S, ones, SSet, h, alpha, sumOdd, sumEven);
+		compute_sum_F_S_T_Iter(nh, S, h, alpha, sumOdd, sumEven);
+		
+		//if sumEven has been assigned some value instead of MARK
+		if(sumEven != MARK){
+			logMinus_New(sumOdd, sumEven);
+			//LOGMINUS_NEW(sumOdd, sumEven);
+		}
+						
+		h[S] = sumOdd;
+	}
+	
+}
+
+
+void rec_compute_AA_S(
+	int d,
+	int nh,
+	int S,
+	int ones, //leng of SSet
+    int * SSet,
+    int T,
+    int sizeofT,
+    int * TSet,
+    double * AA_S,
+    double ** alpha
+){
+
+	if(d < nh){	
+		rec_compute_AA_S(d+1, nh, S, ones, SSet, T, sizeofT, TSet, AA_S, alpha);
+		
+		//the binary repre of S has 1 in d'th position (starting from 0) from right
+		if(((S>>d) & 1) == 1){
+			TSet[sizeofT] = d;
+			rec_compute_AA_S(d+1, nh, S, ones, SSet, (T | (1<<d)), (sizeofT + 1), TSet, AA_S, alpha);
+		}
+	}
+	else{
+		if(T == 0){
+			//cout << "\nT == 0" << endl;
+			
+			AA_S[T] = 0;
+			
+			//cout << "AA_S[" << T << "]=" << AA_S[T] << endl;
+			//cout << endl;
+		}
+		else{
+			//TSet[0]: the right most element
+			//T - (1<<TSet[0]): the binary repre of the the left most biggest 
+			//                 subset of T (including all the elements of T except the right most element)
+			//(1<<nh) - 1 - S: V - S
+			//AA_S[T - (1<<TSet[0])] must have been computed
+			
+			//cout << "\nT = " << T << endl;
+			//cout << "sizeofT = " << sizeofT << endl;
+			//printIntArr(TSet, sizeofT);
+			//cout << "alpha[" << TSet[0] << "][" << (1<<nh) - 1 - S << "] = " 
+			//     << alpha[TSet[0]][(1<<nh) - 1 - S] << endl;
+			//cout << "AA_S[" << T - (1<<TSet[0]) << "] = " << AA_S[T - (1<<TSet[0])] << endl;
+			
+			//AA_S[T] = alpha[TSet[0]][(1<<nh) - 1 - S] + AA_S[T - (1<<TSet[0])];
+			
+			if(alpha[TSet[0]][(1<<nh) - 1 - S] == LOG_ZERO){
+				AA_S[T] = LOG_ZERO;
+			}
+			else{
+				AA_S[T] = alpha[TSet[0]][(1<<nh) - 1 - S] + AA_S[T - (1<<TSet[0])];
+			}
+			
+			
+			//cout << "AA_S[" << T << "]=" << AA_S[T] << endl;
+			//cout << endl;
+		}
+		
+	}
+
+
+}
+
+
+
+//	Given S \subset V,
+//  Compute all T \subset S: AA_S(T) at one time
+void compute_AA_S(
+	int nh,
+	int S,
+	int ones, //leng of SSet
+    int * SSet,
+    double * AA_S,
+    double ** alpha
+){
+	int * TSet = new int[ones];
+	rec_compute_AA_S(0, nh, S, ones, SSet, 0, 0, TSet, AA_S, alpha);
+	delete[] TSet;
+	
+}
+
+
+//	Given S \subset V,
+//  Compute all T \subset S: AA_S(T) at one time
+void compute_AA_S_Iter(
+	int nh,
+	int S,
+	//int ones, //leng of SSet
+    //int * SSet,
+    double * AA_S,
+    double ** alpha
+){
+
+    AA_S[0] = 0;
+    
+	//all the subsets T of S, expect that T = 0 (empty set)
+	for(int T = 1; T <= S; T++){
+		//if T is the subset of S
+		if((S|T) == S){
+			
+			int TSetLen = binIntLen(T);
+			int * TSet = new int[TSetLen];
+			binIntToDecIntArr(T, TSetLen, TSet);
+			
+			//TSet[TSetLen-1]: the left most element
+			//T - (1<<TSet[TSetLen-1]): the binary repre of the the right most biggest 
+			//                 subset of T (including all the elements of T except the left most element)
+			//(1<<nh) - 1 - S: V - S
+			//AA_S[T - (1<<TSet[TSetLen-1])] must have been computed
+			
+			AA_S[T] = alpha[TSet[TSetLen-1]][(1<<nh) - 1 - S] + AA_S[T - (1<<TSet[TSetLen-1])];
+			
+			delete[] TSet;
+	
+		}
+	
+	}
+	
+}
+
+
+
+
+//HR: Compute RF(S, T)
+double compute_RF(int nh, int S, int T, int sizeofT, int *TSet, double * RR, double ** alpha){
+	
+	double result_RF = MARK;
+	double result_RR = RR[S-T];
+	//LOGADD(result_f, result_h);
+	result_RF = result_RR;
+	if (result_RF == LOG_ZERO){
+		return result_RF;
+	}
+
+	//AA_S has been computed
+	if(AA_S[T] == LOG_ZERO){
+		return LOG_ZERO;
+	}
+	else{
+		result_RF += AA_S[T];
+	}
+	
+	return result_RF;
+}
+
+
+
+//HR: Compute: \sum_{k=1}^{|S|} (-1)^{k+1} \sum{T \subset S, |T|=k} RF(S,T)
+void rec_compute_sum_RF_S_T(
+	int d,
+	int nh, 
+	int S, 
+	int ones, 
+	int * SSet, 
+	int T,
+	int sizeofT,
+	int * TSet,
+	double * RR, 
+	double ** alpha, 
+	double & sumOdd, 
+	double & sumEven)
+{
+
+	if(d < nh){	
+		rec_compute_sum_RF_S_T(d+1, nh, S, ones, SSet, T, sizeofT, TSet, RR, alpha, sumOdd, sumEven);
+		
+		//the binary repre of S has 1 in d'th position (starting from 0) from right
+		if(((S>>d) & 1) == 1){
+			TSet[sizeofT] = d;
+			rec_compute_sum_RF_S_T(d+1, nh, S, ones, SSet, (T | (1<<d)), (sizeofT + 1), TSet, RR, alpha, sumOdd, sumEven);
+		}
+	}
+	else{
+		if(sizeofT == 0){
+			//|T| == 0 does not count
+			return;
+			
+		}
+		else{	
+			double result_RF = compute_RF(nh, S, T, sizeofT, TSet, RR, alpha);
+			//sizeofT is odd
+			if(sizeofT % 2 == 1){				
+				logAdd_New(sumOdd, result_RF);
+				//LOGADD_NEW(sumOdd, result_RF);
+			}
+			else{
+				logAdd_New(sumEven, result_RF);
+				//LOGADD_NEW(sumEven, result_RF);
+			}
+			
+		}
+		
+	}
+} // end rec_compute_sum_RF_S_T()
+
+	
+
+
+//HR: Compute: \sum_{k=1}^{|S|} (-1)^{k+1} \sum{T \subset S, |T|=k} RF(S,T)
+void compute_sum_RF_S_T(int nh, int S, int ones, int * SSet, double * RR, double ** alpha, double & sumOdd, double & sumEven){
+	
+	int * TSet = new int[ones];
+	rec_compute_sum_RF_S_T(0, nh, S, ones, SSet, 0, 0, TSet, RR, alpha, sumOdd, sumEven);
+	delete[] TSet;
+	
+} // end compute_sum_RF_S_T()
+
+
+void sub_RR_Fast(int d, int nh, int S, int ones, int * SSet, double * RR, double ** alpha){
+	if(d < nh){
+		sub_RR_Fast(d+1, nh, S, ones, SSet, RR, alpha);
+		
+		SSet[ones] = d;
+		sub_RR_Fast(d+1, nh, S|(1<<d), (ones+1), SSet, RR, alpha);
+	}
+	//d == nh
+	else{
+		//S = 0 means that S is the empty set; 
+		if(S == 0){
+			RR[0] = 0; //log(1) = 0;
+			//cout << "R[0] = 0" << endl;
+			//cout << endl;
+			return;
+
+		}
+		
+		//For given S, for all T \subset S, Compute AA_S[T]
+		AA_S = new double [S+1];
+		compute_AA_S(nh, S, ones, SSet, AA_S, alpha);	
+		//compute_AA_S_Iter(nh, S, AA_S, alpha);		
+		
+		double sumOdd = MARK;
+		double sumEven = MARK;
+		
+		
+		compute_sum_RF_S_T(nh, S, ones, SSet, RR, alpha, sumOdd, sumEven);
+			
+		//if sumEven has been assigned some value instead of MARK
+		if(sumEven != MARK){
+			logMinus_New(sumOdd, sumEven);
+			//LOGMINUS_NEW(sumOdd, sumEven);
+		}
+	
+		RR[S] = sumOdd;
+		delete[] AA_S;
+
+	}
+} // end sub_RR_Fast()
+
+
+
+
+
+
+
+//HR:
+//	Compute RR
+void compute_RR_Fast(int nh, double * RR, double ** alpha){
+	int * SSet = new int[nh];
+	
+	sub_RR_Fast(0, nh, 0, 0, SSet, RR, alpha);
+	
+	delete[] SSet;
+	
+} // end compute_RR_Fast()
+
+
+//HR: Compute: \sum_{k=1}^{|S|} (-1)^{k+1} \sum{T \subset S, |T|=k} RF(S,T)
+void rec_compute_sum_RF_S_T_TSize(
+	int d,
+	int nh, 
+	int S, 
+	int ones, 
+	int * SSet, 
+	int T,
+	int sizeofT,
+	int * TSet,
+	double * RR, 
+	double ** alpha, 
+	double & sumOdd, 
+	double & sumEven,
+	int TSize)
+{
+
+	if(d < nh){	
+		rec_compute_sum_RF_S_T_TSize(d+1, nh, S, ones, SSet, T, sizeofT, TSet, RR, alpha, sumOdd, sumEven, TSize);
+		
+		//the binary repre of S has 1 in d'th position (starting from 0) from right
+		if(((S>>d) & 1) == 1){
+			TSet[sizeofT] = d;
+			rec_compute_sum_RF_S_T_TSize(d+1, nh, S, ones, SSet, (T | (1<<d)), (sizeofT + 1), TSet, RR, alpha, sumOdd, sumEven, TSize);
+		}
+	}
+	else{
+		if(sizeofT == 0){
+			//|T| == 0 does not count
+			return;
+			
+		}
+		//HR: new for *_TSize
+		else if(sizeofT <= TSize){
+		//else{	
+			double result_RF = compute_RF(nh, S, T, sizeofT, TSet, RR, alpha);
+			//sizeofT is odd
+			if(sizeofT % 2 == 1){				
+				logAdd_New(sumOdd, result_RF);
+				//LOGADD_NEW(sumOdd, result_RF);
+			}
+			else{
+				logAdd_New(sumEven, result_RF);
+				//LOGADD_NEW(sumEven, result_RF);
+			}
+			
+		}
+		
+	}
+} // end rec_compute_sum_RF_S_T_TSize()
+
+
+//HR: Compute: \sum_{k=1}^{|S|} (-1)^{k+1} \sum{T \subset S, |T|=k} RF(S,T)
+void compute_sum_RF_S_T_TSize(int nh, int S, int ones, int * SSet, double * RR, double ** alpha, double & sumOdd, double & sumEven, int TSize){
+	
+	int * TSet = new int[ones];
+	rec_compute_sum_RF_S_T_TSize(0, nh, S, ones, SSet, 0, 0, TSet, RR, alpha, sumOdd, sumEven, TSize);
+	delete[] TSet;
+	
+	
+}// end compute_sum_RF_S_T_TSize()
+
+
+
+void sub_RR_Fast_TSize(int d, int nh, int S, int ones, int * SSet, double * RR, double ** alpha, int TSize){
+	if(d < nh){
+		sub_RR_Fast_TSize(d+1, nh, S, ones, SSet, RR, alpha, TSize);
+		
+		SSet[ones] = d;
+		sub_RR_Fast_TSize(d+1, nh, S|(1<<d), (ones+1), SSet, RR, alpha, TSize);
+	}
+	//d == nh
+	else{
+		//S = 0 means that S is the empty set; 
+		if(S == 0){
+			RR[0] = 0; //log(1) = 0;
+			//cout << "R[0] = 0" << endl;
+			//cout << endl;
+			return;
+
+		}
+		
+		//For given S, for all T \subset S, Compute AA_S[T]
+		AA_S = new double [S+1];
+		compute_AA_S(nh, S, ones, SSet, AA_S, alpha);	
+		//compute_AA_S_Iter(nh, S, AA_S, alpha);		
+		
+		double sumOdd = MARK;
+		double sumEven = MARK;
+		
+		
+		compute_sum_RF_S_T_TSize(nh, S, ones, SSet, RR, alpha, sumOdd, sumEven, TSize);
+		
+		
+		//if sumEven has been assigned some value instead of MARK
+		if(sumEven != MARK){
+			logMinus_New_NonnegRes(sumOdd, sumEven);
+			//logMinus_New(sumOdd, sumEven);
+			//LOGMINUS_NEW(sumOdd, sumEven);
+		}
+	
+		RR[S] = sumOdd;
+		delete[] AA_S;
+
+	}
+}// end sub_RR_Fast_TSize()
+
+
+//HR:
+//	Compute RR
+void compute_RR_Fast_TSize(int nh, double * RR, double ** alpha, int TSize){
+	int * SSet = new int[nh];
+	
+	sub_RR_Fast_TSize(0, nh, 0, 0, SSet, RR, alpha, TSize);
+	 
+	delete[] SSet;
+	
+} // end compute_RR_Fast_TSize()
+
+
+
+void rec_compute_sum_F_S_T_Details(
+	int d,
+	int nh, 
+	int S, 
+	int ones, 
+	int * SSet, 
+	int T,
+	int sizeofT,
+	int * TSet,
+	double * h, 
+	double ** alpha, 
+	double & sumOdd, 
+	double & sumEven,
+	double * FSum)
+{
+
+	if(d < nh){	
+		rec_compute_sum_F_S_T_Details(d+1, nh, S, ones, SSet, T, sizeofT, TSet, h, alpha, sumOdd, sumEven, FSum);
+		
+		//the binary repre of S has 1 in d'th position (starting from 0) from right
+		if(((S>>d) & 1) == 1){
+			TSet[sizeofT] = d;
+			rec_compute_sum_F_S_T_Details(d+1, nh, S, ones, SSet, (T | (1<<d)), (sizeofT + 1), TSet, h, alpha, sumOdd, sumEven, FSum);
+		}
+	}
+	else{
+		if(sizeofT == 0){
+			//|T| == 0 does not count
+			return;
+			
+		}
+		else{	
+			double result_F = compute_f(nh, S, T, sizeofT, TSet, h, alpha);
+			//sizeofT is odd
+			if(sizeofT % 2 == 1){				
+				logAdd_New(sumOdd, result_F);
+				//LOGADD_NEW(sumOdd, result_F);
+				
+			}
+			else{
+				logAdd_New(sumEven, result_F);
+				//LOGADD_NEW(sumEven, result_F);
+			}
+			
+			//FSum[sizeofT - 1] stores the sum_{T \subset S, |T| = sizeofT} F(S,T)
+			logAdd_New(FSum[sizeofT - 1], result_F);
+			//LOGADD_NEW(sumEven, result_F);
+			
+			
+			
+		}
+		
+	}
+}
+
+	
+
+void compute_sum_F_S_T_Details(int nh, int S, int ones, int * SSet, double * h, double ** alpha, double & sumOdd, double & sumEven, double * FSum){
+	//TSet's length is set to be SSet's length: it is fine since each TSetLen is set accordingly to each subset T;
+	int * TSet = new int[ones];
+	rec_compute_sum_F_S_T_Details(0, nh, S, ones, SSet, 0, 0, TSet, h, alpha, sumOdd, sumEven, FSum);
+	delete[] TSet;
+	
+	
+}
+
+
+void sub_h_Fast_Details(int d, int nh, int S, int ones, int * SSet, double * h, double ** alpha, double * FSum){
+	if(d < nh){
+		
+		//this kind of ordering in the recursion is crucial; 
+		//it guarantee when you want to compute one S, all its subsets have already been computed
+		sub_h_Fast_Details(d+1, nh, S, ones, SSet, h, alpha, FSum);
+		
+		SSet[ones] = d;
+		sub_h_Fast_Details(d+1, nh, S|(1<<d), (ones+1), SSet, h, alpha, FSum);
+	}
+	//d == nh
+	else{
+		//S = 0 means that S is the empty set; 
+		if(S == 0){
+			h[0] = 0; //log(1) = 0;
+			//cout << "h[0] = 0" << endl;
+			//cout << endl;
+			return;
+
+		}
+		
+		//if S > 0
+		double sumOdd = MARK;
+		double sumEven = MARK;
+		
+		if(S < (1 << nh) - 1 ){
+			compute_sum_F_S_T(nh, S, ones, SSet, h, alpha, sumOdd, sumEven);
+		}
+		//S == (1 << nh) - 1
+		else{
+			compute_sum_F_S_T_Details(nh, S, ones, SSet, h, alpha, sumOdd, sumEven, FSum);
+		}
+		
+		
+		//if sumEven has been assigned some value instead of MARK
+		if(sumEven != MARK){
+			logMinus_New(sumOdd, sumEven);
+			//LOGMINUS_NEW(sumOdd, sumEven);
+		}
+						
+		h[S] = sumOdd;
+		//cout << "h[" << S << "] = "  << h[S] << endl;
+		//cout << endl;
+
+	}
+}
+
+
+void compute_h_Fast_Details(int nh, double * h, double ** alpha, double * FSum){
+	int * SSet = new int[nh];
+	
+	sub_h_Fast_Details(0, nh, 0, 0, SSet, h, alpha, FSum);
+	
+	delete [] SSet;
+	
+}
+
+
+
+//HR: Compute: \sum_{k=1}^{|S|} (-1)^{k+1} \sum{T \subset S, |T|=k} F(S,T)
+void rec_compute_sum_F_S_T_TSize(
+	int d,
+	int nh, 
+	int S, 
+	int ones, 
+	int * SSet, 
+	int T,
+	int sizeofT,
+	int * TSet,
+	double * h, 
+	double ** alpha, 
+	double & sumOdd, 
+	double & sumEven,
+	int TSize)
+{
+
+	if(d < nh){	
+		rec_compute_sum_F_S_T_TSize(d+1, nh, S, ones, SSet, T, sizeofT, TSet, h, alpha, sumOdd, sumEven, TSize);
+		
+		//the binary repre of S has 1 in d'th position (starting from 0) from right
+		if(((S>>d) & 1) == 1){
+			TSet[sizeofT] = d;
+			rec_compute_sum_F_S_T_TSize(d+1, nh, S, ones, SSet, (T | (1<<d)), (sizeofT + 1), TSet, h, alpha, sumOdd, sumEven, TSize);
+		}
+	}
+	else{
+		if(sizeofT == 0){
+			//|T| == 0 does not count
+			return;
+			
+		}
+		//else {	
+		//HR: new for *_TSize
+		else if(sizeofT <= TSize){
+			double result_F = compute_f(nh, S, T, sizeofT, TSet, h, alpha);
+			//sizeofT is odd
+			if(sizeofT % 2 == 1){		
+			//if(sizeofT % 2 == 1 && sizeofT <= TSize){			
+				logAdd_New(sumOdd, result_F);
+				//LOGADD_NEW(sumOdd, result_F);
+			}
+			else{
+			//else if(sizeofT % 2 == 0 && sizeofT <= TSize){
+				logAdd_New(sumEven, result_F);
+				//LOGADD_NEW(sumEven, result_F);
+			}
+			
+		}
+		
+	}
+} // end rec_compute_sum_F_S_T_TSize()
+
+
+
+//HR: Compute: \sum_{k=1}^{|S|} (-1)^{k+1} \sum{T \subset S, |T|=k} F(S,T)
+void compute_sum_F_S_T_TSize(int nh, int S, int ones, int * SSet, double * h, double ** alpha, double & sumOdd, double & sumEven, int TSize){
+	//TSet's length is set to be SSet's length: it is fine since each TSetLen is set accordingly to each subset T;
+	int * TSet = new int[ones];
+	rec_compute_sum_F_S_T_TSize(0, nh, S, ones, SSet, 0, 0, TSet, h, alpha, sumOdd, sumEven, TSize);
+	delete[] TSet;
+	
+	
+}
+
+
+void sub_h_Fast_TSize(int d, int nh, int S, int ones, int * SSet, double * h, double ** alpha, int TSize){
+	if(d < nh){
+		
+		//this kind of ordering in the recursion is crucial; 
+		//it guarantee when you want to compute one S, all its subsets have already been computed
+		sub_h_Fast_TSize(d+1, nh, S, ones, SSet, h, alpha, TSize);
+		
+		SSet[ones] = d;
+		sub_h_Fast_TSize(d+1, nh, S|(1<<d), (ones+1), SSet, h, alpha, TSize);
+	}
+	//d == nh
+	else{
+		//S = 0 means that S is the empty set; 
+		if(S == 0){
+			h[0] = 0; //log(1) = 0;
+			//cout << "h[0] = 0" << endl;
+			//cout << endl;
+			return;
+
+		}
+		
+		//if S > 0
+		double sumOdd = MARK;
+		double sumEven = MARK;
+		
+		
+		compute_sum_F_S_T_TSize(nh, S, ones, SSet, h, alpha, sumOdd, sumEven, TSize);
+		
+		//if sumEven has been assigned some value instead of MARK
+		if(sumEven != MARK){
+			logMinus_New_NonnegRes(sumOdd, sumEven);
+			//logMinus_New(sumOdd, sumEven);
+			//LOGMINUS_NEW(sumOdd, sumEven);
+		}
+						
+		h[S] = sumOdd;
+		//cout << "h[" << S << "] = "  << h[S] << endl;
+		//cout << endl;
+
+	}
+} // end sub_h_Fast_TSize()
+
+
+//HR:
+//	Compute h
+
+void compute_h_Fast_TSize(int nh, double * h, double ** alpha, int TSize){
+	int * SSet = new int[nh];
+	
+	sub_h_Fast_TSize(0, nh, 0, 0, SSet, h, alpha, TSize);
+	
+	delete [] SSet;
+	
+} // end compute_h_Fast_TSize()
+
+
+
+
+void rec_compute_sum_RF_S_T_Details(
+	int d,
+	int nh, 
+	int S, 
+	int ones, 
+	int * SSet, 
+	int T,
+	int sizeofT,
+	int * TSet,
+	double * RR, 
+	double ** alpha, 
+	double & sumOdd, 
+	double & sumEven,
+	double * RFSum)
+{
+
+	if(d < nh){	
+		rec_compute_sum_RF_S_T_Details(d+1, nh, S, ones, SSet, T, sizeofT, TSet, RR, alpha, sumOdd, sumEven, RFSum);
+		
+		//the binary repre of S has 1 in d'th position (starting from 0) from right
+		if(((S>>d) & 1) == 1){
+			TSet[sizeofT] = d;
+			rec_compute_sum_RF_S_T_Details(d+1, nh, S, ones, SSet, (T | (1<<d)), (sizeofT + 1), TSet, RR, alpha, sumOdd, sumEven, RFSum);
+		}
+	}
+	else{
+		if(sizeofT == 0){
+			//|T| == 0 does not count
+			return;
+			
+		}
+		else{	
+			double result_RF = compute_RF(nh, S, T, sizeofT, TSet, RR, alpha);
+			//sizeofT is odd
+			if(sizeofT % 2 == 1){				
+				logAdd_New(sumOdd, result_RF);
+				//LOGADD_NEW(sumOdd, result_RF);
+			}
+			else{
+				logAdd_New(sumEven, result_RF);
+				//LOGADD_NEW(sumEven, result_RF);
+			}
+			//RFSum[sizeofT - 1] stores the sum_{T \subset S, |T| = sizeofT} RF(S,T)
+			logAdd_New(RFSum[sizeofT - 1], result_RF);
+			//LOGADD_NEW(RFSum[sizeofT - 1], result_RF);
+			
+		}
+		
+	}
+}
+
+
+
+void compute_sum_RF_S_T_Details(int nh, int S, int ones, int * SSet, double * RR, double ** alpha, double & sumOdd, double & sumEven, double * RFSum){
+	
+	int * TSet = new int[ones];
+	rec_compute_sum_RF_S_T_Details(0, nh, S, ones, SSet, 0, 0, TSet, RR, alpha, sumOdd, sumEven, RFSum);
+	delete[] TSet;
+	
+	
+}
+
+
+
+void sub_RR_Fast_Details(int d, int nh, int S, int ones, int * SSet, double * RR, double ** alpha, double * RFSum){
+	if(d < nh){
+		sub_RR_Fast_Details(d+1, nh, S, ones, SSet, RR, alpha, RFSum);
+		
+		SSet[ones] = d;
+		sub_RR_Fast_Details(d+1, nh, S|(1<<d), (ones+1), SSet, RR, alpha, RFSum);
+	}
+	//d == nh
+	else{
+		//S = 0 means that S is the empty set; 
+		if(S == 0){
+			RR[0] = 0; //log(1) = 0;
+			//cout << "R[0] = 0" << endl;
+			//cout << endl;
+			return;
+
+		}
+		
+		//For given S, for all T \subset S, Compute AA_S[T]
+		AA_S = new double [S+1];
+		compute_AA_S(nh, S, ones, SSet, AA_S, alpha);	
+		//compute_AA_S_Iter(nh, S, AA_S, alpha);		
+		
+		double sumOdd = MARK;
+		double sumEven = MARK;
+		
+		if(S < (1 << nh) - 1 ){
+			compute_sum_RF_S_T(nh, S, ones, SSet, RR, alpha, sumOdd, sumEven);
+		}
+		//if S == (1 << nh) - 1
+		else{
+			compute_sum_RF_S_T_Details(nh, S, ones, SSet, RR, alpha, sumOdd, sumEven, RFSum);
+		}
+		
+
+		
+		//if sumEven has been assigned some value instead of MARK
+		if(sumEven != MARK){
+			logMinus_New(sumOdd, sumEven);
+			//LOGMINUS_NEW(sumOdd, sumEven);
+		}
+	
+		RR[S] = sumOdd;
+		delete[] AA_S;
+
+	}
+}
+
+
+
+
+void compute_RR_Fast_Details(int nh, double * RR, double ** alpha, double * RFSum){
+	int * SSet = new int[nh];
+	
+	sub_RR_Fast_Details(0, nh, 0, 0, SSet, RR, alpha, RFSum);
+	
+	delete[] SSet;
+	
+}
+
+
+
+
+
+
+
+
+//HR: Compute: \sum_{k=1}^{|S|} (-1)^{k+1} \sum{T \subset S, |T|=k} RF(S,T)
+void compute_sum_RF_S_T_Iter(int nh, int S, double * RR, double ** alpha, double & sumOdd, double & sumEven){
+	
+	//all the subsets T of S, expect that T = 0 (empty set)
+	for(int T = 1; T <= S; T++){
+		//if T is the subset of S
+		if((S|T) == S){
+			
+			int TSetLen = binIntLen(T);
+			int * TSet = new int[TSetLen];
+			binIntToDecIntArr(T, TSetLen, TSet);
+
+			double result_RF = compute_RF(nh, S, T, TSetLen, TSet, RR, alpha);
+			//sizeofT is odd
+			if(TSetLen % 2 == 1){				
+				logAdd_New(sumOdd, result_RF);
+				//LOGADD_NEW(sumOdd, result_RF);
+			}
+			else{
+				logAdd_New(sumEven, result_RF);
+				//LOGADD_NEW(sumEven, result_RF);
+			}
+			
+			delete[] TSet;
+		}
+
+	}
+
+}
+
+//HR:
+//	Compute RR
+void compute_RR_Iter(int nh, double * RR, double ** alpha){
+	
+	
+	RR[0] = 0; //log(1) = 0;
+	
+	//RR[S]: S: from 1 to (1<<nh)-1
+	for(int S = 1; S < (1<<nh); S++){
+		//For given S, for all T \subset S, Compute AA_S[T]
+		AA_S = new double [S+1];
+		//compute_AA_S(nh, S, ones, SSet, AA_S, alpha);	
+		compute_AA_S_Iter(nh, S, AA_S, alpha);		
+		
+		double sumOdd = MARK;
+		double sumEven = MARK;
+				
+		compute_sum_RF_S_T_Iter(nh, S, RR, alpha, sumOdd, sumEven);
+		
+		//if sumEven has been assigned some value instead of MARK
+		if(sumEven != MARK){
+			logMinus_New(sumOdd, sumEven);
+			//LOGMINUS_NEW(sumOdd, sumEven);
+		}
+	
+		RR[S] = sumOdd;
+		//cout << "RR[" << S << "] = "  << h[S] << endl;
+		//cout << endl;
+
+		delete[] AA_S;
+		
+	}
+	
+}
+
+
+	
+	
+//	Given v \in V, given U \subset V-{v},
+//  Compute all T \subset (V-{v}-U) = TMax: Eta_U(T) at one time
+void rec_compute_Eta_U(
+	int d,
+	int nh,
+	int TMax,
+	int TMaxSetLen, //leng of SSet
+    int * TMaxSet,
+    int U,
+    int T,
+    int sizeofT,
+    int * TSet,
+    double * Eta_U,
+    double ** alpha
+){
+
+	if(d < nh){	
+		rec_compute_Eta_U(d+1, nh, TMax, TMaxSetLen, TMaxSet, U, T, sizeofT, TSet, Eta_U, alpha);
+		
+		//the binary repre of S has 1 in d'th position (starting from 0) from right
+		if(((TMax>>d) & 1) == 1){
+			TSet[sizeofT] = d;
+			rec_compute_Eta_U(d+1, nh, TMax, TMaxSetLen, TMaxSet, U, (T | (1<<d)), (sizeofT + 1), TSet, Eta_U, alpha);
+		}
+	}
+	else{
+		if(T == 0){
+			//cout << "\nT == 0" << endl;
+			
+			Eta_U[T] = 0;
+			
+			//cout << "AA_S[" << T << "]=" << AA_S[T] << endl;
+			//cout << endl;
+		}
+		else{
+	
+			Eta_U[T] = alpha[TSet[0]][U] + Eta_U[T - (1<<TSet[0])];
+			
+		}
+		
+	}
+
+
+}
+
+
+//	Given v \in V, given U \subset V-{v},
+//  Compute all T \subset (V-{v}-U) = TMax: Eta_U(T) at one time
+void compute_Eta_U(
+	int nh,
+	int TMax,
+	int TMaxSetLen, 
+    int * TMaxSet,
+    int U,
+    double * Eta_U,
+    double ** alpha
+){
+	int * TSet = new int[TMaxSetLen];
+	rec_compute_Eta_U(0, nh, TMax, TMaxSetLen, TMaxSet, U, 0, 0, TSet, Eta_U, alpha);
+	delete[] TSet;
+	
+}
+
+
+//	Given v \in V, given U \subset V-{v},
+//  Compute K_v(U)
+void rec_compute_K_v(
+	int d,
+	int nh,
+	int TMax,
+	int TMaxSetLen, //leng of SSet
+    int * TMaxSet,
+    int U,
+    int T,
+    int sizeofT,
+    int * TSet,
+    double * Eta_U,
+    double * K_v,
+    double ** alpha,
+    double & sumOdd,
+    double & sumEven
+){
+	if(d < nh){	
+		rec_compute_K_v(d+1, nh, TMax, TMaxSetLen, TMaxSet, U, T, sizeofT, TSet, Eta_U, K_v, alpha, sumOdd, sumEven);
+		
+		//the binary repre of S has 1 in d'th position (starting from 0) from right
+		if(((TMax>>d) & 1) == 1){
+			TSet[sizeofT] = d;
+			rec_compute_K_v(d+1, nh, TMax, TMaxSetLen, TMaxSet, U, (T | (1<<d)), (sizeofT + 1), TSet, Eta_U, K_v, alpha, sumOdd, sumEven);
+		}
+	}
+	else{
+			double abs_result = RRFunc[TMax - T] + Eta_U[T];
+			//sizeofT is odd
+			if(sizeofT % 2 == 1){				
+				logAdd_New(sumOdd, abs_result);
+				//LOGADD_NEW(sumOdd, abs_result);
+			}
+			else{
+				logAdd_New(sumEven, abs_result);
+				//LOGADD_NEW(sumEven, abs_result);
+			}
+		
+	}	
+	
+}
+
+
+
+//	Given v \in V, given U \subset V-{v},
+//  Compute K_v(U)
+void compute_K_v(
+	int nh,
+	double * K_v,
+	int v, 
+	int U,
+	int USetLen, 
+    int * USet,
+    double * Eta_U,
+    double ** alpha
+)
+
+
+{
+	//TMax = V - {v} - U
+	int TMax = (1 << nh) - 1 - (1 << v) - U;
+	int TMaxSetLen = nh - 1 - USetLen;
+	int * TMaxSet = new int[TMaxSetLen];
+	binIntToDecIntArr(TMax, TMaxSetLen, TMaxSet);
+	
+	//compute Eta_U
+	Eta_U = new double [TMax+1];
+	compute_Eta_U(nh, TMax, TMaxSetLen, TMaxSet, U, Eta_U, alpha);
+	
+//	totalOfMixEta << "\nEta_U U == " << U << ":" << endl;
+//	for(int index = 0; index < TMax+1; index++){
+//		totalOfMixEta << "Eta_U[" << index << "] = " << Eta_U[index] << endl;
+//	}
+	 
+	
+	
+	double result_K_v;
+	
+	double sumOdd = MARK;
+	double sumEven = MARK;
+	
+	int * TSet = new int[TMaxSetLen];
+	rec_compute_K_v(0, nh, TMax, TMaxSetLen, TMaxSet, U, 0, 0, TSet, Eta_U, K_v, alpha, sumOdd, sumEven);
+	
+	
+
+	// sum of itmes of Even sizeofT must have been assigned some value, since |T| = 0 must exist
+	// sum of itmes of odd sizeofT must have been assigned some value could be not assigned so that it is still MARK
+	if(sumOdd == MARK){
+		result_K_v = sumEven;	
+	}
+	else{
+		//if sum of itmes of Even sizeofT > sum of itmes of odd sizeofT, so that log(sumEven) > log(sumOdd) 
+		if(sumEven > sumOdd){
+			logMinus_New(sumEven, sumOdd);
+			//LOGMINUS_NEW(sumEven, sumOdd);
+			//result_K_v is log(prob) < 0
+			result_K_v = sumEven;
+		}
+		else if(sumEven == sumOdd){
+			//cerr<<"\n******Original K_v is zero!" << endl;
+			result_K_v = LOG_ZERO;
+		}
+		else{
+			//cerr<<"\n******Original K_v is negetive!" << endl;
+			logMinus_New(sumOdd, sumEven);
+			//LOGMINUS_NEW(sumOdd, sumEven);
+			//use result_K_v > 0 to represent negative prob
+			result_K_v = - sumOdd;
+			//cerr<<"K_v[" << U << "]" << result_K_v << endl;
+		}
+	}
+	K_v[U] = result_K_v;	
+	
+	delete[] TSet;
+	
+	delete[] Eta_U;			
+	
+
+}
+
+
+void init_all_K_v(int d, int S, double * K_v, double value, int ones, int nh){
+		//cerr<<" d: "<<d<<endl; 
+	
+	if (d < nh){
+		init_all_K_v(d+1, S, K_v, value, ones, nh);
+			
+		//HR: ones means the num of 1 bits, represeting the num of parents
+		//HR: S | (1 << d) reprents the S could have 1 bit in the d'th if current S does not have
+		//HR: if (ones < k) 
+		init_all_K_v(d+1, S | (1 << d), K_v, value, ones+1, nh);
+	}
+	else{
+		K_v[S] = value;
+		
+	}
+}
+	
+
+void rec_compute_all_K_v(
+	int d,
+	int nh,
+	int UMax,
+	int UMaxSetLen, //leng of SSet
+    int * UMaxSet,
+    int U,
+    int USetLen,
+    int * USet,
+    double * K_v,
+    int v,
+    double ** alpha
+){
+	if(d < nh){	
+		rec_compute_all_K_v(d+1, nh, UMax, UMaxSetLen, UMaxSet, U, USetLen, USet, K_v, v, alpha);
+		
+		//the binary repre of S has 1 in d'th position (starting from 0) from right
+		if(((UMax>>d) & 1) == 1){
+			USet[USetLen] = d;
+			rec_compute_all_K_v(d+1, nh, UMax, UMaxSetLen, UMaxSet, (U | (1<<d)), (USetLen + 1), USet, K_v, v, alpha);
+		}
+	}
+	else{			
+		compute_K_v(nh, K_v, v, U, USetLen, USet, Eta_U, alpha);
+	
+	}	
+	
+}
+
+
+void compute_all_K_v(
+	int nh,
+	double * K_v,
+	int v, 
+    double ** alpha
+){
+	
+	init_all_K_v(0, 0, K_v, LOG_ZERO, 0, nh);
+	
+	//UMax = V - {v} 
+	int UMax = (1 << nh) - 1 - (1 << v);
+	int UMaxSetLen = nh - 1;
+	int * UMaxSet = new int[UMaxSetLen];
+	binIntToDecIntArr(UMax, UMaxSetLen, UMaxSet);
+	
+	int * USet = new int[UMaxSetLen];
+	rec_compute_all_K_v(0, nh, UMax, UMaxSetLen, UMaxSet, 0, 0, USet, K_v, v, alpha);
+	
+	delete[] USet;
+	
+}
+
+