diff options
Diffstat (limited to 'sourcecodes/k-best/src/UpdateHR.h')
| -rw-r--r-- | sourcecodes/k-best/src/UpdateHR.h | 1652 |
1 files changed, 0 insertions, 1652 deletions
diff --git a/sourcecodes/k-best/src/UpdateHR.h b/sourcecodes/k-best/src/UpdateHR.h deleted file mode 100644 index 19b93c1d..00000000 --- a/sourcecodes/k-best/src/UpdateHR.h +++ /dev/null @@ -1,1652 +0,0 @@ -#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; - -} - - |