From c4f926438dcb8abe805e910399940f79ff643c4b Mon Sep 17 00:00:00 2001 From: ziejd2 Date: Thu, 14 Sep 2017 15:39:41 -0500 Subject: Add files via upload --- sourcecodes/k-best/index.html | 6 + sourcecodes/k-best/src/Arguments.h | 195 +++ sourcecodes/k-best/src/Engine.h | 2278 +++++++++++++++++++++++++ sourcecodes/k-best/src/GenerateGraphs.sh | 7 + sourcecodes/k-best/src/Model.h | 1002 +++++++++++ sourcecodes/k-best/src/UpdateHR.h | 1652 ++++++++++++++++++ sourcecodes/k-best/src/UpdateHR2.h | 511 ++++++ sourcecodes/k-best/src/arcToNet.cc | 233 +++ sourcecodes/k-best/src/arcs2dot.c | 60 + sourcecodes/k-best/src/buildk.sh | 23 + sourcecodes/k-best/src/buildk_poster.sh | 9 + sourcecodes/k-best/src/cfg.h | 9 + sourcecodes/k-best/src/clean.sh | 6 + sourcecodes/k-best/src/comp_pHat_D.cc | 213 +++ sourcecodes/k-best/src/data.sh | 27 + sourcecodes/k-best/src/data2net.sh | 35 + sourcecodes/k-best/src/data2netk.sh | 55 + sourcecodes/k-best/src/data2netk_poster.sh | 56 + sourcecodes/k-best/src/data2netk_poster.sh.bk | 56 + sourcecodes/k-best/src/datak.sh | 23 + sourcecodes/k-best/src/drawgraph.sh | 8 + sourcecodes/k-best/src/files.c | 59 + sourcecodes/k-best/src/files.h | 12 + sourcecodes/k-best/src/files.o | Bin 0 -> 3728 bytes sourcecodes/k-best/src/get_best_net.c | 63 + sourcecodes/k-best/src/get_best_order.c | 54 + sourcecodes/k-best/src/get_best_parents.c | 114 ++ sourcecodes/k-best/src/get_best_sinks.c | 78 + sourcecodes/k-best/src/get_kbest_nets | Bin 0 -> 122641 bytes sourcecodes/k-best/src/get_kbest_nets.cc | 1552 +++++++++++++++++ sourcecodes/k-best/src/get_kbest_parents | Bin 0 -> 25198 bytes sourcecodes/k-best/src/get_kbest_parents.cc | 249 +++ sourcecodes/k-best/src/get_local_scores.h | 15 + sourcecodes/k-best/src/graph.sh | 13 + sourcecodes/k-best/src/ilogi.c | 37 + sourcecodes/k-best/src/ilogi.h | 9 + sourcecodes/k-best/src/ls_BDe.c | 163 ++ sourcecodes/k-best/src/ls_BDe.h | 12 + sourcecodes/k-best/src/ls_LOO.c | 152 ++ sourcecodes/k-best/src/ls_LOO.h | 10 + sourcecodes/k-best/src/ls_NML.c | 77 + sourcecodes/k-best/src/ls_NML.h | 12 + sourcecodes/k-best/src/ls_XIC.c | 82 + sourcecodes/k-best/src/ls_XIC.h | 13 + sourcecodes/k-best/src/main.cc | 51 + sourcecodes/k-best/src/net2parents.c | 50 + sourcecodes/k-best/src/parents2arcs.c | 37 + sourcecodes/k-best/src/reg.c | 320 ++++ sourcecodes/k-best/src/reg.h | 7 + sourcecodes/k-best/src/score_net.c | 39 + sourcecodes/k-best/src/score_nets.c | 54 + sourcecodes/k-best/src/varpar.c | 19 + sourcecodes/k-best/src/varpar.h | 9 + sourcecodes/k-best/src/varpar.o | Bin 0 -> 1360 bytes sourcecodes/k-best/src/xtab.h | 80 + 55 files changed, 9906 insertions(+) create mode 100644 sourcecodes/k-best/index.html create mode 100644 sourcecodes/k-best/src/Arguments.h create mode 100644 sourcecodes/k-best/src/Engine.h create mode 100644 sourcecodes/k-best/src/GenerateGraphs.sh create mode 100644 sourcecodes/k-best/src/Model.h create mode 100644 sourcecodes/k-best/src/UpdateHR.h create mode 100644 sourcecodes/k-best/src/UpdateHR2.h create mode 100644 sourcecodes/k-best/src/arcToNet.cc create mode 100644 sourcecodes/k-best/src/arcs2dot.c create mode 100644 sourcecodes/k-best/src/buildk.sh create mode 100644 sourcecodes/k-best/src/buildk_poster.sh create mode 100644 sourcecodes/k-best/src/cfg.h create mode 100644 sourcecodes/k-best/src/clean.sh create mode 100644 sourcecodes/k-best/src/comp_pHat_D.cc create mode 100644 sourcecodes/k-best/src/data.sh create mode 100644 sourcecodes/k-best/src/data2net.sh create mode 100644 sourcecodes/k-best/src/data2netk.sh create mode 100644 sourcecodes/k-best/src/data2netk_poster.sh create mode 100644 sourcecodes/k-best/src/data2netk_poster.sh.bk create mode 100644 sourcecodes/k-best/src/datak.sh create mode 100644 sourcecodes/k-best/src/drawgraph.sh create mode 100644 sourcecodes/k-best/src/files.c create mode 100644 sourcecodes/k-best/src/files.h create mode 100644 sourcecodes/k-best/src/files.o create mode 100644 sourcecodes/k-best/src/get_best_net.c create mode 100644 sourcecodes/k-best/src/get_best_order.c create mode 100644 sourcecodes/k-best/src/get_best_parents.c create mode 100644 sourcecodes/k-best/src/get_best_sinks.c create mode 100644 sourcecodes/k-best/src/get_kbest_nets create mode 100644 sourcecodes/k-best/src/get_kbest_nets.cc create mode 100644 sourcecodes/k-best/src/get_kbest_parents create mode 100644 sourcecodes/k-best/src/get_kbest_parents.cc create mode 100644 sourcecodes/k-best/src/get_local_scores.h create mode 100644 sourcecodes/k-best/src/graph.sh create mode 100644 sourcecodes/k-best/src/ilogi.c create mode 100644 sourcecodes/k-best/src/ilogi.h create mode 100644 sourcecodes/k-best/src/ls_BDe.c create mode 100644 sourcecodes/k-best/src/ls_BDe.h create mode 100644 sourcecodes/k-best/src/ls_LOO.c create mode 100644 sourcecodes/k-best/src/ls_LOO.h create mode 100644 sourcecodes/k-best/src/ls_NML.c create mode 100644 sourcecodes/k-best/src/ls_NML.h create mode 100644 sourcecodes/k-best/src/ls_XIC.c create mode 100644 sourcecodes/k-best/src/ls_XIC.h create mode 100644 sourcecodes/k-best/src/main.cc create mode 100644 sourcecodes/k-best/src/net2parents.c create mode 100644 sourcecodes/k-best/src/parents2arcs.c create mode 100644 sourcecodes/k-best/src/reg.c create mode 100644 sourcecodes/k-best/src/reg.h create mode 100644 sourcecodes/k-best/src/score_net.c create mode 100644 sourcecodes/k-best/src/score_nets.c create mode 100644 sourcecodes/k-best/src/varpar.c create mode 100644 sourcecodes/k-best/src/varpar.h create mode 100644 sourcecodes/k-best/src/varpar.o create mode 100644 sourcecodes/k-best/src/xtab.h (limited to 'sourcecodes/k-best') diff --git a/sourcecodes/k-best/index.html b/sourcecodes/k-best/index.html new file mode 100644 index 00000000..429668c7 --- /dev/null +++ b/sourcecodes/k-best/index.html @@ -0,0 +1,6 @@ + + +Bayesian Network Web Server + + + diff --git a/sourcecodes/k-best/src/Arguments.h b/sourcecodes/k-best/src/Arguments.h new file mode 100644 index 00000000..574b2bf7 --- /dev/null +++ b/sourcecodes/k-best/src/Arguments.h @@ -0,0 +1,195 @@ +#ifndef ARGUMENTS_H +#define ARGUMENTS_H + +#include +#include + +using namespace std; + +class Arguments { +public: + static char* datafile; + static char* layeringfile; + static char* maxindegree; + static char* model; + static char* task; + static char* maxnumrecords; + + //HR: Add for in_out features + static char* inFeasFileName; + static char* outFeasFileName; + // + + //HR: Add for in_out features + static int option; + // + + //HR: Add for ADtree option + static int ADtree; + // + + //HR: Add for Topk dir option + static char* directoryName; + // + + + + static void init(int argc, char **args){ + for(int i = 1; i < argc; i ++){ + + if(args[i][0]=='-' && args[i][1]=='d' && args[i][2]=='\0'){ + int j = i + 1; + while (j < argc && args[j][0]!='-') { + Arguments::datafile = args[j]; + j ++; + } + } + else if(args[i][0]=='-' && args[i][1]=='l' && args[i][2]=='\0'){ + int j = i + 1; + while (j < argc && args[j][0]!='-') { + Arguments::layeringfile = args[j]; + j ++; + } + } + else if(args[i][0]=='-' && args[i][1]=='m' && args[i][2]=='\0'){ + int j = i + 1; + while (j < argc && args[j][0]!='-') { + Arguments::maxindegree = args[j]; + j ++; + } + } + else if(args[i][0]=='-' && args[i][1]=='u' && args[i][2]=='\0'){ + int j = i + 1; + while (j < argc && args[j][0]!='-') { + Arguments::maxnumrecords = args[j]; + j ++; + } + } + else if(args[i][0]=='-' && args[i][1]=='M' && args[i][2]=='\0'){ + int j = i + 1; + while (j < argc && args[j][0]!='-') { + Arguments::model = args[j]; + j ++; + } + } + else if(args[i][0]=='-' && args[i][1]=='T' && args[i][2]=='\0'){ + int j = i + 1; + while (j < argc && args[j][0]!='-') { + Arguments::task = args[j]; + j ++; + } + } + //HR: Add for in_out feature + else if(args[i][0]=='-' && args[i][1]=='i' && args[i][2]=='n' && args[i][3]=='\0'){ + int j = i + 1; + while (j < argc && args[j][0]!='-') { + Arguments::inFeasFileName = args[j]; + j ++; + } + } + else if(args[i][0]=='-' && args[i][1]=='o' && args[i][2]=='u' && args[i][3]=='t' && args[i][4]=='\0'){ + int j = i + 1; + while (j < argc && args[j][0]!='-') { + Arguments::outFeasFileName = args[j]; + j ++; + } + } + // + //HR: Add for in_out feature + else if(args[i][0]=='-' && args[i][1]=='o' && args[i][2]=='p' && args[i][3]=='t' && args[i][4]=='\0'){ + int j = i + 1; + while (j < argc && args[j][0]!='-') { + Arguments::option = atoi(args[j]); + j ++; + } + } + // + //HR: Add for ADtree feature + else if(args[i][0]=='-' && args[i][1]=='a' && args[i][2]=='d' && args[i][3]=='\0'){ + int j = i + 1; + while (j < argc && args[j][0]!='-') { + Arguments::ADtree = atoi(args[j]); + j ++; + } + } + // + //HR: Add for Topk dir + else if(args[i][0]=='-' && args[i][1]=='d' && args[i][2]=='i' && args[i][3]=='r' && args[i][4]=='\0'){ + int j = i + 1; + while (j < argc && args[j][0]!='-') { + Arguments::directoryName = args[j]; + j++; + } + } + // + } + + print_arguments(stderr); + + } + static void print_arguments(FILE *f){ + fprintf(f, " -d Data file:\n"); + fprintf(f, " %62s\n", Arguments::datafile); + //fprintf(f, " -l Layering file:\n"); + //fprintf(f, " %62s\n", Arguments::layeringfile); + fprintf(f, " -m Maximum indegree:\n"); + fprintf(f, " %62s\n", Arguments::maxindegree); + fprintf(f, " -u Maximum number of data records read:\n"); + fprintf(f, " %62s\n", Arguments::maxnumrecords); + //fprintf(f, " -M Model:\n"); + //fprintf(f, " %62s\n", Arguments::model); + //fprintf(f, " -T Task (Infer=I, Generate=G):\n"); + //fprintf(f, " %62s\n", Arguments::task); + + //HR: Add for in_out + fprintf(f, " -opt Option:\n"); + fprintf(f, " %62d\n", Arguments::option); + if(Arguments::option == 5){ + fprintf(f, " -in In_Feature File:\n"); + fprintf(f, " %62s\n", Arguments::inFeasFileName); + fprintf(f, " -out Out_Feature File:\n"); + fprintf(f, " %62s\n", Arguments::outFeasFileName); + } + // + //HR: Add for ADtree + fprintf(f, " -ad ADtree:\n"); + fprintf(f, " %62d\n", Arguments::ADtree); + // + //HR: Add for Topk dir + fprintf(f, " -dir directory name:\n"); + fprintf(f, " %62s\n", Arguments::directoryName); + // + } +}; + +char* Arguments::datafile = "testdata.dat"; +char* Arguments::layeringfile = "%"; +char* Arguments::maxindegree = "3"; +char* Arguments::model = "M"; +char* Arguments::task = "I"; +char* Arguments::maxnumrecords = "999999"; + + +//HR: For test +//char* Arguments::datafile = "cases/iris.idt"; +//char* Arguments::layeringfile = "%"; +//char* Arguments::maxindegree = "4"; +//char* Arguments::model = "M"; +//char* Arguments::task = "I"; +//char* Arguments::maxnumrecords = "150"; + + +//HR: Add for in_out +char* Arguments::inFeasFileName = "inFeature.txt"; +char* Arguments::outFeasFileName = "outFeature.txt"; +int Arguments::option = 5; + +//HR: Add for ADtree +int Arguments::ADtree = 0; + + +//HR: Add for Topk dir option +char* Arguments::directoryName = "./"; +// + +#endif 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 +#include + +#include"Arguments.h" +#include"Model.h" + +#include +#include + +#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<log_lcp(0, T)<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 "<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 "<layer(h, &Vh, &nh); + model->upper_layers(h-1, &Vu, &nu); + + + beta = new double*[nh]; + alpha = new double*[nh]; + + k = model->max_indegree(); + + cerr<<" . "<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."<print_edge_prob( + cout, Vh[i], Vh[j], exp(log_prob - gb[(1<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<<" . "<print_edge_prob(cout, Vh[vv], Vh[uu], exp(gfNume[(1<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<<" . "<print_edge_prob(cout, Vh[vv], Vh[uu], exp(hFuncNume[(1<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."<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<<" . "<print_edge_prob(cout, Vh[vv], Vh[uu], exp(RRFuncNume[(1<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."<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<<" . "<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 "<print_edge_prob(cout, Vh[vv], Vh[uu], exp(RRFuncNume[(1<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."<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<<" . "<makeADTree(); + } + + for (int j = 0; j < nh; j ++){ + + beta[j] = new double[1 << nh]; + + //cerr<<" Start compute_beta for node "<freeADTree(); + } + + RRFunc = new double[1 << nh]; + + cerr<<"\n Compute_RR_Fast(): "<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<<" . "<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 "<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 "<>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 = "<print_edge_prob( + cout, Vh[i], Vh[v], exp(log_prob_nume - RRFunc[(1<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<<" . "<>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 = "<print_edge_prob( + cout, Vh[i], Vh[v], exp(log_prob_nume - RRFunc[(1<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<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<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: "<>= 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<<": "< 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<<": "< 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 diff --git a/sourcecodes/k-best/src/GenerateGraphs.sh b/sourcecodes/k-best/src/GenerateGraphs.sh new file mode 100644 index 00000000..3427eb4d --- /dev/null +++ b/sourcecodes/k-best/src/GenerateGraphs.sh @@ -0,0 +1,7 @@ +#!/bin/bash +k=100 +rdir=./../iris3/resdir +for (( i = 0 ; i < k; i++ )) +do + dot ${rdir}/dot_$i -Tps -o ${rdir}/dotRes_$i.ps +done \ No newline at end of file diff --git a/sourcecodes/k-best/src/Model.h b/sourcecodes/k-best/src/Model.h new file mode 100644 index 00000000..1ba7a0dd --- /dev/null +++ b/sourcecodes/k-best/src/Model.h @@ -0,0 +1,1002 @@ +#ifndef MODEL_H +#define MODEL_H + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include"Arguments.h" + +using namespace std; + +#include"UpdateHR2.h" + +#define MAX_COUNT 1000 +#define MAX_RQ 8*10000 +static double lr[MAX_RQ][MAX_COUNT]; +static double lr2[MAX_RQ][MAX_COUNT]; + +//ofstream totalOf1("debugModelTotal1.txt"); + +double log_gammaratio(int a, int z){ + if(z >= MAX_COUNT || a >= MAX_RQ){ +// if(z >= MAX_COUNT){ +// cerr << " Counts larger than " << MAX_COUNT << " occured. " << endl; +// cerr << " Increase MAX_COUNT. Exit now.\n"; +// +// } +// else if(a >= MAX_RQ ){ +// cerr << " a larger than " << MAX_RQ << " occured. " << endl; +// cerr << " Increase MAX_RQ. Exit now.\n"; +// } +// +// exit(1); + + double log_gam_res = 0; + //seperate 2 cases: + //if only z >= MAX_COUNT but a < MAX_RQ + if(a < MAX_RQ){ + //can compute it faster based on lr + log_gam_res += lr[a][MAX_COUNT - 1]; + for(int k = MAX_COUNT - 1; k <= z-1; k++){ + log_gam_res += log(1.0/a + k); + } + + } + //o.w.: start from the beginning + else{ + for(int k = 0; k <= z-1; k++){ + log_gam_res += log(1.0/a + k); + } + } + +// totalOf1 << "log_gam_res = " << log_gam_res << endl << endl; + return log_gam_res; + + } + else if (z < MAX_COUNT && a > 0){ + return lr[a][z]; + } + else if (z < MAX_COUNT && a < 0){ + return lr2[-a][z]; + } + else{ + cerr << "Error: a = 0" << endl; + exit(1); + } +} + + +double log_gammaratio(double z){ + return 0.50 * log(M_PI) + (z - 0.50) * log(z) - z + 1.0 / (12 * z); +} + +template void printvec(ostream& f, const vector& v){ + int size= v.size(); + int i; + for(i=0; i void printvecs(ostream& f, const vector >& v){ + int size= v.size(); + int i; + for(i=0; i 0){ + if (T & 1){ + f<<" "<>= 1; + } +} +void print_set(ostream& f, int T){ + int j = 0; + while (T > 0){ + if (T & 1){ + f<<" "<>= 1; + } +} + + +// Reads and stores data. +class Data{ +public: + Data(){} + ~Data(){} + void init(){ + read_data(); + downcode(); + } + void read_data(){ + ifstream ifs(Arguments::datafile, ios::in); + if (!ifs){ + fprintf(stderr, " Cannot read file %s.\n", Arguments::datafile); + exit(1); + } + fprintf(stderr, " Reading file %s...\n", Arguments::datafile); + + char* buffer= new char[10000]; + ifs.getline(buffer, 10000); + char* pch= strtok(buffer,", \t"); + while (pch != NULL){ + string tempstring(pch); + //cerr << "Node = "<< tempstring << endl; + heads.push_back(tempstring); + pch = strtok(NULL, ", \t"); + } + numattributes = heads.size(); + + fprintf(stderr, " Heading read: %d attributes.\n", numattributes); + + dm.clear(); + vector temp; + while (true){ + temp.clear(); + ifs.getline(buffer, 10000); + pch = strtok(buffer,", \t"); + while (pch != NULL){ + temp.push_back(atoi(pch)); + pch = strtok(NULL, ", \t"); + } + if ((int)temp.size() != numattributes) break; + + dm.push_back(temp); + + if ((int)dm.size() >= atoi(Arguments::maxnumrecords)) break; + } + numrecords = dm.size(); + fprintf(stderr, " Data read: %d lines.\n", numrecords); + delete [] buffer; + } + void downcode(){ + //fprintf(stderr, "downcode() starts. \n"); + + int inuse[4096][256]; + for (int v = 0; v < 255; v ++){ + for (int i = 0; i < numattributes; i ++){ + inuse[v][i] = -1; + } + } + //Purpose: For each i attribute, + for (int t = 0; t < numrecords; t ++){ + for (int i = 0; i < numattributes; i ++){ + inuse[dm[t][i]][i] = 1; + } + } + + //ofstream of1("debugDowncode.txt"); + //of1 << "\n inuse check" << endl; + for (int v = 0; v < 30; v ++){ + + +// for (int i = 0; i < numattributes; i++){ +// of1 << " inuse[" << v << "][" << i << "] =" << inuse[v][i] << endl; +// } +// of1 << "\n"; + + + /* + FILE * fp1 = fopen("debug1.txt", "w"); + fprintf(fp1, "\n inuse check\n"); + for (int i = 0; i < numattributes; i++){ + fprintf(fp1, " inuse[%d][%d] = %d\n", v, i, inuse[v][i]); + } + fclose(fp1); + */ + + + /* + fprintf(stderr, "\n inuse check\n"); + for (int i = 0; i < numattributes; i++){ + fprintf(stderr, " inuse[%d][%d]=%d\n", v, i, inuse[v][i]); + } + fprintf(stderr, "\n"); + */ + } + + arities.clear(); + + maxarity = 0; + + for (int i = 0; i < numattributes; i ++){ + int count = 0; + for (int v = 0; v < 255; v ++){ + if (inuse[v][i] == 1){ + inuse[v][i] = count ++; + } + } + arities.push_back(count); + if (count > maxarity) maxarity = count; + } + +// of1 << "maxarity = " << maxarity << endl; +// printvec(of1, arities); + + for (int t = 0; t < numrecords; t ++){ + for (int i = 0; i < numattributes; i ++){ + dm[t][i] = inuse[dm[t][i]][i]; + + } + } + }//end downcode() + + + int get_index(string s){ + int i = 0; + while (i < (int)heads.size() && heads[i] != s){ i ++; } + return i; + } + void print_data(ostream & f){ + for (int i = 0; i < numattributes; i ++){ + f << " " << heads[i]; + } + f << endl; + for (int t = 0; t < numrecords; t ++){ + for (int i = 0; i < numattributes; i ++){ + f << " " << dm[t][i]; + } + f << endl; + } + } + + vector heads; + + vector< vector > dm; + + vector arities; + + int numattributes; + + int numrecords; + + int maxarity; +}; + + +// NOTE: We assume that children come first, that is, +// V[0] is the grandest child and V[n-1] is the grandest parent. +// Layers respectively from 0 to numlayers. +// Note that this is reverse to the input ordering. +class Layering{ +public: + Layering(){} + ~Layering(){ delete [] V;} + + void init(Data & data){ + set_layers(data); + + //print_layers(); + + } + void set_layers(Data & data){ + +// ofstream of1("debugSet_layers.txt"); + + V = new int[data.numattributes]; + int j = 0; + cnh.clear(); nh.clear(); cnh.clear(); + + if (Arguments::layeringfile[0] == '%'){ + numlayers = 1; + int layersize = data.numattributes; + nh.push_back(layersize); + cnh.push_back(j); + + for (int l = 0; l < layersize; l ++){ + int i = j; + V[j++] = i; + } + + edgeswithin.clear(); + edgeswithin.push_back(true); + + return; + } + + parselayeringfile(Arguments::layeringfile, true); + numlayers = nodesinlayers.size(); + + for (int h = numlayers-1; h >= 0; h --){ + int layersize = nodesinlayers[h].size(); + nh.push_back(layersize); + cnh.push_back(j); + for (int l = 0; l < layersize; l ++){ + int i = data.get_index(nodesinlayers[h][l]); + V[j++] = i; + } + } + + } + void print_layers(){ + for (int h = 0; h < numlayers; h ++){ + cerr<<"Layer "< tempstringvec; //will hold the tokenized line + while(pch != NULL){ + string tempstring(pch); + tempstringvec.push_back(tempstring); + pch = strtok(NULL, ", \t"); + } + + if(tempstringvec.size() < 3){ + cerr<<" ERROR: in reading file '"< tempnodesinlayer; + for(i = 2; i < (int)tempstringvec.size();i ++){ + tempnodesinlayer.push_back(tempstringvec.at(i)); + } + nodesinlayers.push_back(tempnodesinlayer); + + ifs>>tempchar; + ifs.putback(tempchar); + if(ifs.eof()){ + break; + } + } + delete [] buffer; + buffer = 0; + pch = 0; + cerr.width(13); + cerr<<"Layer Name"; + cerr.width(24); + cerr<<"Edges Within Allowed"; + cerr.width(14); + cerr<<"Node Names"< layernames; + vector edgeswithin; + vector > nodesinlayers; + int *V; + int numlayers; + vector nh; + vector cnh; +}; + + + +// Blocks for computing sufficient statistics. +class Tnode{ +public: + double pseudocount; + short int arity; + short int count; + short int depth; + Tnode** children; + + Tnode(const int a, const int d){ + pseudocount = 1; + arity = a; + count = 0; + depth = d; + children = new Tnode*[arity]; + for(int v = 0; v < arity; v ++){ + children[v] = NULL; + } + } + ~Tnode(){ + for(int v = 0; v < arity; v ++){ + if (children[v] != NULL){ + delete children[v]; + } + } + delete [] children; + //free_memory(); + } + void free_memory(){ + //cerr<<" Deleting d="<print(f); + } + } + } + + double evaluate(int d, int a, int r){ + double sum = 0; + + //if no parent, directly log[gamma(1+count) / gamma(1) ] (1/a = 1/1 = 1) + if (depth == d){ + sum = log_gammaratio(a, count); + //cerr << "count1: " << count << endl; + + } + else if (depth <= d - 1){ + for (int v = 0; v < arity; v ++){ + if (children[v] != NULL){ + double value = children[v]->evaluate(d, a, r); + sum += value; + } + } + if (depth == d - 1){ + sum -= log_gammaratio(r, count); + //cerr << "count2: " << count << endl; + + } + } + return sum; + } + + + double evaluateHR(int d, int alpha_v_pa, int r){ + //cerr << "start: root.evaluateHR(d =" + //<< d << ", alpha_v_pa=" << alpha_v_pa << ", r=" << r << ")" << endl; + double sum = 0; + + //if no parent, directly log[gamma(1+count) / gamma(1) ] (1/a = 1/1 = 1) + if (depth == d){ + //sum = log_gammaratio(a, count); + sum = log_gammaratio(alpha_v_pa, count); + //cerr << "count1: " << count << endl; + + } + else if (depth <= d - 1){ + for (int v = 0; v < arity; v ++){ + if (children[v] != NULL){ + double value = children[v]->evaluateHR(d, alpha_v_pa, r); + sum += value; + } + } + if (depth == d - 1){ + int alpha_pa = alpha_v_pa / r; + + sum -= log_gammaratio(alpha_pa, count); +// totalOf1 << "\nalpha_pa: " << alpha_pa << endl; +// totalOf1 << "count2: " << count << endl; +// totalOf1 << "-log_gammaratio(alpha_pa, count) = " << -log_gammaratio(alpha_pa, count) << endl; + //cerr << "count2: " << count << endl; + + } + } + return sum; + } + +}; + + +class Model{ +public: + Model(){} + ~Model(){} + + void makeADTree(vector< vector > & dm, vector & arities){ + //Create dmIndex from dm, dmIndex stroes {0, 1, R-1} + vector dmIndex; +// cout << "dm.size() = " << dm.size() << endl; + for(int i = 0; i < (int) dm.size(); i++){ + dmIndex.push_back(i); + } +// cout << "dmIndex:" << endl; +// for(int i = 0; i < (int) dm.size(); i++){ +// cout << "dmIndex[" << i << "] = "<< dmIndex[i] << endl; +// } + + // M: the num of attr +// int M = dm[0].size(); +// cout << "M = " << M << endl; + + //cout << "Start MakeADTree()" << endl; + //ADTreeNode * ADTreeNodeP1 = MakeADTree(0, dmIndex, dm, arities); + ADTreeNodeP1 = MakeADTree(0, dmIndex, dm, arities); + //cout << "End MakeADTree()" << endl; + + //cout << "Start PrintADTreeNode()" << endl; + //PrintADTreeNode(ADTreeNodeP1, 0, M, arities); + //cout << "End PrintADTreeNode()" << endl; + + } + + void makeADTree(){ + makeADTree(data.dm, data.arities); + } + + void freeADTree(vector< vector > & dm, vector & arities){ + int M = dm[0].size(); + //cout << "M = " << M << endl; + //cout << "Start FreeADTreeNode()" << endl; + FreeADTreeNode(ADTreeNodeP1, M, arities); + //cout << "End FreeADTreeNode()" << endl; + + } + + void freeADTree(){ + freeADTree(data.dm, data.arities); + } + + + void MakeADContiTable(vector< vector > & dm, vector & arities){ + cout << "\nStart MakeADContiTable(vector< vector > & dm, vector & arities)" << endl; + // M: the num of attr + int M = dm[0].size(); + cout << "M = " << M << endl; + + deque inquiryAttrs; + + for(int i = 0; i < (int) inquiryAttrs.size(); i++){ + cout << "inquiryAttrs[" << i << "] = " << inquiryAttrs[i] << endl; + } + + vector givenAttrs; + + vector givenVals; + + cout << "Call MakeContab() " << endl; + CondiContiTable resultTable = MakeContab(inquiryAttrs, ADTreeNodeP1, givenAttrs, givenVals, arities); + cout << "Return MakeContab() " << endl; + for(unsigned int i = 0; i < resultTable.contiTableCounts.size(); i++){ + cout << "resultTable.contiTableCounts[" << i << "] = " << resultTable.contiTableCounts[i] << endl; + } + + cout << "End MakeADContiTable(vector< vector > & dm, vector & arities)" << endl; + + + } //end void MakeADContiTable() + + void MakeADContiTable(){ + MakeADContiTable(data.dm, data.arities); + } + + vector MakeADContiTable(deque inquiryAttrs, vector< vector > & dm, vector & arities){ + //cout << "\nStart MakeADContiTable(vector< vector > & dm, vector & arities)" << endl; + // M: the num of attr + //int M = dm[0].size(); + //cout << "M = " << M << endl; + + //deque inquiryAttrs; + + //Full inquiry +// for(int i = 0; i < M; i++){ +// inquiryAttrs.push_back(i); +// } + + +// for(int i = 0; i < (int) inquiryAttrs.size(); i++){ +// cout << "inquiryAttrs[" << i << "] = " << inquiryAttrs[i] << endl; +// } + + vector givenAttrs; + + vector givenVals; + + //cout << "Call MakeContab() " << endl; + CondiContiTable resultTable = MakeContab(inquiryAttrs, ADTreeNodeP1, givenAttrs, givenVals, arities); + //cout << "Return MakeContab() " << endl; +// +// for(unsigned int i = 0; i < resultTable.contiTableCounts.size(); i++){ +// cout << "resultTable.contiTableCounts[" << i << "] = " << resultTable.contiTableCounts[i] << endl; +// } +// +// cout << "End MakeADContiTable(vector< vector > & dm, vector & arities)" << endl; + return resultTable.contiTableCounts; + + } //end vector MakeADContiTable() + + + + void init(){ + data.init(); + m = data.numrecords; + n = data.numattributes; + + layering.init(data); + + for (int a = 1; a < MAX_RQ; a ++){ + double aa = 1.0 / a; + lr[a][0] = 0; lr2[a][0] = 0; + + for (int k = 1; k < MAX_COUNT; k ++){ + lr[a][k] = lr[a][k-1] + log(aa + k - 1); + lr2[a][k] = lr2[a][k-1] + log(a + k - 1); + } + } + + } + + + + + + + + void test(){ + for (int i = 0; i < 1; i ++){ + vector T; + T.clear(); + T.push_back(1); + T.push_back(2); + + } + exit(1); + } + + + void testHR(){ + for (int i = 0; i < 1; i ++){ + vector T; + T.clear(); + T.push_back(1); + T.push_back(2); + T.push_back(3); + T.push_back(4); + + } + } //end void testHR() + + + void testHR_ADtree(){ + for (int i = 0; i < 1; i ++){ + deque T; + T.clear(); + T.push_back(1); + T.push_back(2); + T.push_back(3); + T.push_back(4); + + } + } //end void testHR_ADtree() + + + + //--- Interface functions --- below ---------- + // + // Log of local conditional probability, log p(xi | xS). + double log_lcp(int i, vector & T){ + + return log_lcp_mult(i, T); + } + double log_lcp(int i, int *T, int d){ + vector TT; + for (int j = 0; j < d; j ++) TT.push_back(T[j]); + return log_lcp_mult(i, TT); + } + + + double log_lcpHR(int i, vector & T){ + + return log_lcp_multHR(i, T); + } + + double log_lcpHR(int i, int *T, int d){ + vector TT; + for (int j = 0; j < d; j ++) TT.push_back(T[j]); + //cerr << "start log_lcp_multHR(i, TT)" << endl; + double result = log_lcp_multHR(i, TT); + return result; + } + + + double log_lcpHR_ADtree(int i, int *T, int d){ + deque TT; + for (int j = 0; j < d; j ++) TT.push_back(T[j]); + double result = log_lcp_multHR_ADtree(i, TT); + return result; + } + + + double log_prior(int i, int *T, int d){ + return lr[1][d] + lr[1][n-1-d] - lr[1][n-1]; + } + int num_layers(){ + return layering.numlayers; + } + void layer(int h, int** Vh, int* nh){ + *Vh = &(layering.V[layering.cnh[h]]); + *nh = layering.nh[h]; + } + void upper_layers(int h, int** Vu, int* nu){//similar to above but upper means parent + if (h < layering.numlayers){ + *Vu = &(layering.V[layering.cnh[h+1]]); + *nu = n - layering.cnh[h+1]; + } + else{ *Vu = NULL; *nu = 0;} + } + void lower_layers(int h, int **Vl, int *nl){//similar to above but lower means children + if (h > 0){ + *Vl = layering.V; + *nl = layering.cnh[h]; + } + else{ *Vl = NULL; *nl = 0;} + } + int max_indegree(){ + return atoi(Arguments::maxindegree); + } + bool edges_within(int h){ + return layering.edges_within(h); + } + int num_nodes(){ + return n; + } + + void print_edge_prob(ostream& f, int i, int j, double p){ + //f< "< "< & T){ + vector S; + + S.push_back(i); + + for (int k = 0; k < (int)T.size(); k ++){ + S.push_back(T[k]); + } + + //cerr << " Evaluation yields: " << evaluate(S) << " , log_gammaratio(1, 2) = " + // << log_gammaratio(1, 2) << " , a = " << data.arities[S[0]] << endl; + + return evaluate(S); + } + + double evaluate(vector & S){ + + Tnode root(data.maxarity+1, 0); + vector u; + + //for each record t + for (int t = 0; t < m; t ++){ + u.clear(); + for (int k = 0; k < (int)S.size(); k ++){ + int j = S[k]; + u.push_back(data.dm[t][j]); + } + //cerr << t+1 << "th update: "; + + root.update(S, u); + //cerr << endl; + + } + //root.print(cerr); + + + return root.evaluate(S.size(), 1, -data.arities[S[0]]); + + } + + + double log_lcp_multHR(int i, vector & T){ + vector S; + + S.push_back(i); + + for (int k = 0; k < (int)T.size(); k ++){ + S.push_back(T[k]); + } + + //cerr << " Evaluation yields: " << evaluate(S) << " , log_gammaratio(1, 2) = " + // << log_gammaratio(1, 2) << " , a = " << data.arities[S[0]] << endl; + double result = evaluateHR(S); + return result; + + + }//end log_lcp_multHR() + + + double log_lcp_multHR_ADtree(int i, deque & T){ + //totalOf1 << "Inside log_lcp_multHR()" << endl; + + double result2 = evaluateHR_ADtree(T); + + int insertIndex = -1; + for(int k = 0; k < (int)T.size(); k++){ + if(i < T[k]){ + insertIndex = k; + break; + } + } + if(insertIndex == -1){ + insertIndex = (int)T.size(); + } + + deque S; + + for(int k = 0; k <= insertIndex - 1; k++){ + S.push_back(T[k]); + } + S.push_back(i); + for(int k = insertIndex; k < (int)T.size(); k++){ + S.push_back(T[k]); + } + + double result1 = evaluateHR_ADtree(S); + + //cerr << " Evaluation yields: " << evaluate(S) << " , log_gammaratio(1, 2) = " + // << log_gammaratio(1, 2) << " , a = " << data.arities[S[0]] << endl; + double result = result1 - result2; + return result; + + } + + + // Dirichlet-multinomial model + double evaluateHR(vector & S){ + Tnode * root = new Tnode (data.maxarity+1, 0); + vector u; + + //for each record t + for (int t = 0; t < m; t ++){ + u.clear(); + for (int k = 0; k < (int)S.size(); k ++){ + int j = S[k]; + u.push_back(data.dm[t][j]); + } + //cerr << t+1 << "th update: "; + + (*root).update(S, u); + // + //cerr << endl; + + } + //root.print(cerr); + + int alpha_v_pa = 1; + for(int index = 0; index < (int)S.size(); index++){ + alpha_v_pa *= data.arities[S[index]]; + } + int r = data.arities[S[0]]; + //return root.evaluateHR(S.size(), alpha_v_pa, r); + + double result = (*root).evaluateHR(S.size(), alpha_v_pa, r); + //totalOf1 << " result of root.evaluateHR(" << S.size() <<", " + // << alpha_v_pa << ", " << r << ") = " << result << endl; + delete root; + return result; + + } //end double evaluateHR(vector & S) + + + double evaluateHR_ADtree(deque & S){ + + int alpha = 1; + for(int index = 0; index < (int)S.size(); index++){ + alpha *= data.arities[S[index]]; + } + + vector countVec = MakeADContiTable(S, data.dm, data.arities); + + double sum = 0; + for(int i = 0; i < (int) countVec.size(); i++){ + //totalOf1 << "countVec[" << i <<" ]" << countVec[i] << endl; + sum += log_gammaratio(alpha, countVec[i]); + } + +// totalOf1 << "evaluateHR_ADtree(deque & S) ends" << endl; + + return sum; + + } //end double evaluateHR_ADtree() + + + Data data; + int n; + int m; + Layering layering; + + ADTreeNode * ADTreeNodeP1; + +}; //end of class Model + + + + + +#endif 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 +#include +#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 + +//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< 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<>d) & 1) == 1){ + TSet[sizeofT] = d; + rec_compute_sum_RF_S_T_TSize(d+1, nh, S, ones, SSet, (T | (1<>d) & 1) == 1){ + TSet[sizeofT] = d; + rec_compute_sum_F_S_T_Details(d+1, nh, S, ones, SSet, (T | (1< 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< 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< 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) & 1) == 1){ + USet[USetLen] = d; + rec_compute_all_K_v(d+1, nh, UMax, UMaxSetLen, UMaxSet, (U | (1< +#include + +#include +#include +#include +#include +#include +#include + +#include"Arguments.h" +#include"Model.h" + +using namespace std; + + +struct ADTreenode{ + int count; + + //i.e. start attr which = *; + int startVaryNodeInd; + + //struct ADVarynode * ADVaryNodePs[]; + struct ADVarynode ** ADVaryNodePs; +}; + +typedef struct ADTreenode ADTreeNode; + +struct ADVarynode{ + //Represent Vary Index + int index; + + int mcv; + + //struct ADTreenode * ADTreeNodePs[]; + struct ADTreenode ** ADTreeNodePs; + +}; + +typedef struct ADVarynode ADVaryNode; + + +typedef struct { + //the count for each arrangment of contiTableAttrs + //contiTableCounts.size() == the product of the arity of each element in contiTableAttrs + //do not explicitly record the whole contiTable + vector contiTableCounts; + + //invariant: in the decreasing order + vector contiTableAttrs; + + //condiAttrs.size() == condiVals.size() + //invariant: in the increasing order + vector condiAttrs; + vector condiVals; + + //invariant: contiTableAttrs.size() + condiAttrs.size() == the whole original inquiry length + // finally: condiAttrs.size() == 0 +} CondiContiTable; + + + + + +ADTreeNode * MakeADTree(int, vector &, vector< vector > &, vector &); +ADVaryNode * MakeVaryNode(int, vector &, vector< vector > &, vector &); +void PrintADTreeNode(ADTreeNode *, int, const int, vector &); +void PrintADVaryNode(ADVaryNode *, int, const int, vector &); +void FreeADTreeNode(ADTreeNode *, const int, vector &); +void FreeADVaryNode(ADVaryNode *, const int, vector &); + + +CondiContiTable MakeContab(deque, ADTreeNode *, + vector &, vector &, const vector &); + +void MinusContab(CondiContiTable &, CondiContiTable *, int, int); + +void ConcatContab(CondiContiTable *, int, int, CondiContiTable &); + + + + + +//HR: +// a_i is in {0, ..., M-1} +// dm: the data matrix: R * M +// dmIndex: initially stores {0, 1, ..., R - 1}, it is just the light-weight index for dm +// arities: stores the arity of each attr +// could use vector & dmIndex instead of vector dmIndex +ADTreeNode * MakeADTree(int a_i, vector & dmIndex, vector< vector > & dm, vector & arities){ + //M = num of attr + int M = dm[0].size(); + + ADTreeNode * ADTreeNodeP1 = (ADTreeNode *) malloc(sizeof(ADTreeNode)); + ADTreeNodeP1->count = dmIndex.size(); + ADTreeNodeP1->startVaryNodeInd = a_i; + //base + if(a_i >= M){ + ADTreeNodeP1->ADVaryNodePs = NULL; + } + else{ + ADTreeNodeP1->ADVaryNodePs = (ADVaryNode **) malloc(sizeof(ADVaryNode*) * + (M - 1 - a_i + 1)); + } + for(int a_j = a_i; a_j <= M-1; a_j++){ + ADTreeNodeP1->ADVaryNodePs[a_j - a_i] = MakeVaryNode(a_j, dmIndex, dm, arities); + } + return ADTreeNodeP1; + +}//end ADTreeNode * MakeADTree() + + + +//HR: +// a_i is in {0, ..., M-1} +// dm: the data matrix: R * M +// dmIndex: initially stores {0, 1, ..., R - 1} +// arities: stores the arity of each attr +// could use vector & dmIndex instead of vector dmIndex +ADVaryNode * MakeVaryNode(int a_i, vector & dmIndex, vector< vector > & dm, vector & arities){ + ADVaryNode * ADVaryNodeP1 = (ADVaryNode *) malloc(sizeof(ADVaryNode)); + ADVaryNodeP1->index = a_i; + int n_i = arities[a_i];//may be used as global var + ADVaryNodeP1->ADTreeNodePs = (ADTreeNode **) malloc(sizeof(ADTreeNode *) * n_i); + vector dmIndexSub[n_i]; + for(unsigned int j = 0; j <= dmIndex.size() - 1; j++){ + int v = dm[dmIndex[j]][a_i]; + dmIndexSub[v].push_back(dmIndex[j]); + } + ADVaryNodeP1->mcv = 0; + int maxCount = dmIndexSub[0].size(); + for(int j = 1; j <= n_i - 1; j++){ + if((int) dmIndexSub[j].size() > maxCount){ + maxCount = dmIndexSub[j].size(); + ADVaryNodeP1->mcv = j; + } + } + for(int j = 0; j <= n_i - 1; j++){ + if(dmIndexSub[j].size() == 0 || j == ADVaryNodeP1->mcv){ + ADVaryNodeP1->ADTreeNodePs[j] = NULL; + } + else{ + ADVaryNodeP1->ADTreeNodePs[j] = MakeADTree(a_i+1, dmIndexSub[j], dm, arities); + } + } + //not necessary in fact + //dmIndex.clear(); + //HR: Neither necessary nor correct + //delete dmIndex; + return ADVaryNodeP1; + + +}//end ADVaryNode * MakeVaryNode() + + +//HR: +// ADTreeNodeP1: pointers to the ADTreeNode +// level: starting from 0 +// M: the num of attr +// arities: stores the arity of each attr +// initial call: PrintADTreeNode(ADTreeNodeP1, 0, M, arities) in Model.h +// pre-order to print the ADtree +void PrintADTreeNode(ADTreeNode * ADTreeNodeP1, int level, const int M, vector & arities){ + string strSpace = ""; + for(int i = 0; i < level; i++){ + strSpace += " "; + } + cout << strSpace << "ADTreeNode" << endl; + cout << strSpace << "Count = " << ADTreeNodeP1->count << endl; + cout << strSpace << "startVaryNodeInd = " << ADTreeNodeP1->startVaryNodeInd << endl; + if(ADTreeNodeP1->ADVaryNodePs == NULL){ + return; + } + else{ + for(int i = 0; i < M - 1 - ADTreeNodeP1->startVaryNodeInd + 1; i++){ + PrintADVaryNode(ADTreeNodeP1->ADVaryNodePs[i], level+1, M, arities); + } + } +} //end void PrintADTreeNode() + + +//HR: +// ADVaryNodeP1: pointers to the ADVaryNode +// level: starting from 0 +// M: the num of attr +// arities: stores the arity of each attr +// pre-order to print the ADtree +void PrintADVaryNode(ADVaryNode * ADVaryNodeP1, int level, const int M, vector & arities){ + string strSpace = ""; + for(int i = 0; i < level; i++){ + strSpace += " "; + } + cout << strSpace << "ADVaryNode" << endl; + cout << strSpace << "index = " << ADVaryNodeP1->index << endl; + cout << strSpace << "mcv = " << ADVaryNodeP1->mcv << endl; + int n_i = arities[ADVaryNodeP1->index]; //may be changed to global var + for(int j = 0; j <= n_i - 1; j++){ + if(j == ADVaryNodeP1->mcv){ + cout << strSpace << " mcv" << endl; + } + else if(ADVaryNodeP1->ADTreeNodePs[j] == NULL){ + cout << strSpace << " NULL" << endl; + } + else{ + PrintADTreeNode(ADVaryNodeP1->ADTreeNodePs[j], level+1, M, arities); + } + } + +}//end void PrintADVaryNode() + + + +//HR: +// ADTreeNodeP1: pointers to the ADTreeNode +// M: the num of attr +// arities: stores the arity of each attr +// initial call: FreeADTreeNode(ADTreeNodeP1, M, arities) in Model.h +// pre-order to free the ADtree +void FreeADTreeNode(ADTreeNode * ADTreeNodeP1, const int M, vector & arities){ + + if(ADTreeNodeP1->ADVaryNodePs != NULL){ + for(int i = 0; i < M - 1 - ADTreeNodeP1->startVaryNodeInd + 1; i++){ + FreeADVaryNode(ADTreeNodeP1->ADVaryNodePs[i], M, arities); + } + free(ADTreeNodeP1->ADVaryNodePs); + } + + free(ADTreeNodeP1); +} //end void FreeADTreeNode() + + +//HR: +// ADVaryNodeP1: pointers to the ADVaryNode +// M: the num of attr +// arities: stores the arity of each attr +// pre-order to free the ADtree +void FreeADVaryNode(ADVaryNode * ADVaryNodeP1, const int M, vector & arities){ + + int n_i = arities[ADVaryNodeP1->index]; //may be changed to global var + for(int j = 0; j <= n_i - 1; j++){ + if(j == ADVaryNodeP1->mcv){ + ; + } + else if(ADVaryNodeP1->ADTreeNodePs[j] == NULL){ + ; + } + else{ + FreeADTreeNode(ADVaryNodeP1->ADTreeNodePs[j], M, arities); + } + } + free(ADVaryNodeP1->ADTreeNodePs); + free(ADVaryNodeP1); + +}//end void FreeADVaryNode() + + + +//Error: if use ofstream of1 instead of ofstream & of1 +//In file included from Model.h:19, +// from main.cc:6: +//UpdateHR2.h: In copy constructor `std::basic_ios >::basic_ios(const std::basic_ios >&)': +///usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../include/c++/3.4.6/bits/ios_base.h:781: error: `std::ios_base::ios_base(const std::ios_base&)' is private +//UpdateHR2.h:220: error: within this context +//UpdateHR2.h: In copy constructor `std::basic_filebuf >::basic_filebuf(const std::basic_filebuf >&)': +///usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../include/c++/3.4.6/streambuf:769: error: `std::basic_streambuf<_CharT, _Traits>::basic_streambuf(const std::basic_streambuf<_CharT, _Traits>&) [with _CharT = char, _Traits = std::char_traits]' is private +//UpdateHR2.h:220: error: within this context +//UpdateHR2.h: In function `void print_CondiContiTable(CondiContiTable&, std::ofstream)': +//UpdateHR2.h:220: error: initializing argument 2 of `void print_vec(std::vector >, std::ofstream)' + + +void print_vec(vector & vec1, ofstream & of1){ + for(int i = 0; i < (int) vec1.size(); i++){ + of1 << vec1[i] << " "; + } + of1 << endl; +} + +void print_deque(deque & v1, ofstream & of1){ + for(int i = 0; i < (int) v1.size(); i++){ + of1 << v1[i] << " "; + } + of1 << endl; +} + +void print_CondiContiTable(CondiContiTable & ccTable1, ofstream & of1){ + of1 << "ccTable1.contiTableCounts: " << endl; + print_vec(ccTable1.contiTableCounts, of1); + + of1 << "ccTable1.contiTableAttrs: " << endl; + print_vec(ccTable1.contiTableAttrs, of1); + + of1 << "ccTable1.condiAttrs: " << endl; + print_vec(ccTable1.condiAttrs, of1); + + of1 << "ccTable1.condiVals: " << endl; + print_vec(ccTable1.condiVals, of1); + +} + + + +void print_vec(vector & vec1){ + for(int i = 0; i < (int) vec1.size(); i++){ + cout << vec1[i] << " "; + } + cout << endl; +} + + +void print_deque(deque & v1){ + for(int i = 0; i < (int) v1.size(); i++){ + cout << v1[i] << " "; + } + cout << endl; +} + +void print_CondiContiTable(CondiContiTable & ccTable1){ + cout << "ccTable1.contiTableCounts: " << endl; + print_vec(ccTable1.contiTableCounts); + + cout << "ccTable1.contiTableAttrs: " << endl; + print_vec(ccTable1.contiTableAttrs); + + cout << "ccTable1.condiAttrs: " << endl; + print_vec(ccTable1.condiAttrs); + + cout << "ccTable1.condiVals: " << endl; + print_vec(ccTable1.condiVals); + +} + + + +//HR: +// inquiryAttrs: {a_i_1, ..., a_i_n} +// ADNP: The current pointer to ADTreeNode +// givenAttrs, givenVals explicilty record the meaning of ADNP +// {b_i_1, ..., b_i_m} = {value_b_i_1, .. value_b_i_m} +// pre: +// inquiryAttrs is in the increasing order +// givenAttrs is in the increasing order +CondiContiTable MakeContab(deque inquiryAttrs, ADTreeNode * ADNP, + vector & givenAttrs, vector & givenVals, const vector & arities){ + +// cout << "\nMakeContab()" << endl; +// cout << "inquiryAttrs" << endl; +// print_deque(inquiryAttrs); +// cout << "givenAttrs" << endl; +// print_vec(givenAttrs); +// cout << "givenVals" << endl; +// print_vec(givenVals); + + + + + //base 1 + if(ADNP == NULL){ + CondiContiTable bottomTable; + + //base 1.1 + //if(inquiryAttrs.size() == 0){ + if(inquiryAttrs.empty()){ + bottomTable.contiTableCounts.push_back(0); + //bottomTable.contiTableAttrs is empty + bottomTable.condiAttrs = givenAttrs; + bottomTable.condiVals = givenVals; + + //cout << "//base 1.1: if(ADNP == NULL) and if(inquiryAttrs.empty())" << endl; + //print_CondiContiTable(bottomTable); + + } + //base 1.2 + //if(inquiryAttrs.size() > 0) + else{ + int rowNum = 1; + for(unsigned int i = 0; i < inquiryAttrs.size(); i++){ + //Note: I explicitly write inquiryAttrs[inquiryAttrs.size() - 1 - i ] instead of + // inquiryAttrs[i] because this is really the ordering to create the whole contiTable + int lastAttr = inquiryAttrs[inquiryAttrs.size() - 1 - i ]; + rowNum *= arities[lastAttr]; + + bottomTable.contiTableAttrs.push_back(lastAttr); + } + bottomTable.condiAttrs = givenAttrs; + bottomTable.condiVals = givenVals; + + for(int i = 0; i < rowNum; i++){ + bottomTable.contiTableCounts.push_back(0); + } + + //cout << "//base 1.2: if(ADNP == NULL) and if(!inquiryAttrs.empty())" << endl; + //print_CondiContiTable(bottomTable); + } + + + + return bottomTable; + + }// end if(ADNP == NULL) + + //base 2 + else if (inquiryAttrs.empty()){ + CondiContiTable bottomTable; + bottomTable.contiTableCounts.push_back(ADNP->count); + //bottomTable.contiTableAttrs is empty + bottomTable.condiAttrs = givenAttrs; + bottomTable.condiVals = givenVals; + + //cout << "//base 2: if(ADNP != NULL) and if(inquiryAttrs.empty())" << endl; + //print_CondiContiTable(bottomTable); + + return bottomTable; + } + else{ + int a_i_1 = inquiryAttrs.front(); + int VNIndex = a_i_1 - ADNP->startVaryNodeInd; + ADVaryNode * VNP = ADNP->ADVaryNodePs[VNIndex]; + int mcv = VNP->mcv; + int n_i_1 = arities[a_i_1]; + CondiContiTable CTs[n_i_1]; + + //delete the 1st element a_i_1 of inquiryAtts, will never use the 1st element a_i_1 again + //but does not allow the sub_program to delete more elements, so use deque inquiryAttrs instead of deque & inquiryAttrs + //vector & givenAttrs, vector & givenVals can be used because the sub_program will change but then restore it + inquiryAttrs.pop_front(); + + for(int k = 0; k < n_i_1; k++){ + if(k != mcv){ + ADTreeNode * ADNP_k = VNP -> ADTreeNodePs[k]; + + vector newGivenAttrs = givenAttrs; + newGivenAttrs.push_back(a_i_1); + vector newGivenVals = givenVals; + newGivenVals.push_back(k); + CTs[k] = MakeContab(inquiryAttrs, ADNP_k, newGivenAttrs, newGivenVals, arities); + + } + } + CondiContiTable sumCTs = MakeContab(inquiryAttrs, ADNP, givenAttrs, givenVals, arities); + MinusContab(sumCTs, CTs, n_i_1, mcv); + + CondiContiTable result; + ConcatContab(CTs, n_i_1, a_i_1, result); + + //cout << "//non-base: " << endl; + //print_CondiContiTable(result); + + return result; + } + +}// end CondiContiTable MakeContab() + + +void MinusContab(CondiContiTable & sumCTs, CondiContiTable * CTs, int n_i_1, int mcv){ + CTs[mcv].contiTableCounts = sumCTs.contiTableCounts; + CTs[mcv].contiTableAttrs = sumCTs.contiTableAttrs; + for(unsigned int i = 0; i < CTs[mcv].contiTableCounts.size(); i++){ + for(int k = 0; k < n_i_1; k++){ + if(k != mcv){ + CTs[mcv].contiTableCounts[i] -= CTs[k].contiTableCounts[i]; + } + } + } + //assume arities of each attr >=2 + if(mcv != 0){ + CTs[mcv].condiAttrs = CTs[0].condiAttrs; + CTs[mcv].condiVals = CTs[0].condiVals; + } + else{ + CTs[mcv].condiAttrs = CTs[1].condiAttrs; + CTs[mcv].condiVals = CTs[1].condiVals; + } + CTs[mcv].condiVals.pop_back(); + CTs[mcv].condiVals.push_back(mcv); + +} //end void MinusContab() + + +void ConcatContab(CondiContiTable * CTs, int n_i_1, int a_i_1, CondiContiTable & concatedTable){ + for(int i = 0; i < n_i_1; i++){ + //for each i, CTs[i].contiTableCounts.size() is the same + for(unsigned int j = 0; j < CTs[i].contiTableCounts.size(); j++){ + concatedTable.contiTableCounts.push_back(CTs[i].contiTableCounts[j]); + } + } + + //check: a_i_1 == CTs[0].condiAttrs.back(); + if(a_i_1 != CTs[0].condiAttrs.back()){ + cout << "Error: a_i_1 != CTs[0].condiAttrs.back()" << endl; + cout << "a_i_1 = " << a_i_1 << endl; + cout << "CTs[0].condiAttrs.back() = " << CTs[0].condiAttrs.back() << endl; + } + + //move a_i_1 from the end of CTs[0].condiAttrs to the end of concatedTable.contiTableAttrs + concatedTable.contiTableAttrs = CTs[0].contiTableAttrs; + concatedTable.contiTableAttrs.push_back(a_i_1); + + concatedTable.condiAttrs = CTs[0].condiAttrs; + concatedTable.condiAttrs.pop_back(); + + concatedTable.condiVals = CTs[0].condiVals; + concatedTable.condiVals.pop_back(); + +} //end void ContatContab() + + + + + + + diff --git a/sourcecodes/k-best/src/arcToNet.cc b/sourcecodes/k-best/src/arcToNet.cc new file mode 100644 index 00000000..b071baab --- /dev/null +++ b/sourcecodes/k-best/src/arcToNet.cc @@ -0,0 +1,233 @@ +#include +#include +#include +//#include "varpar.h" +#include +#include +#include +#include +#include +//#include "cfg.h" +#include +#include +//HR +//#define EPSILON 0.001 +#define EPSILON 0.00001 + + +//HR: Add for hash_set +#if __GNUC__ < 3 && __GNUC__ >= 2 && __GNUC_MINOR__ >= 95 +# include +# include +# define gnu_namespace std +#elif __GNUC__ >= 3 +# include +# if __GNUC_MINOR__ == 0 +# include +# define gnu_namespace std +# else +# include +# define gnu_namespace __gnu_cxx +# endif +#else +# include +# include +# define gnu_namespace std +#endif +using namespace gnu_namespace; +#include +// + + +using namespace std; + + +//HR: add from Engine.h in BayesNw +#define MARK 9e99 +#define LOG_ZERO -1e101 +//#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)); +//HR: add from UpdateHR.h in BayesNw to deal with the problem in sum = 0 in comp_post() +//call by reference to the original value +//the log sum result is stored in logA +//may change double into long double; no, use double +void logAddComp(double & logA, double & logB){ + if(logA == MARK){ + logA = logB; + } + else if (logA == LOG_ZERO){ + logA = logB; + } + else if(logB == MARK){ + logA = logA; + } + else if(logB == LOG_ZERO){ + logA = logA; + } + //Nither logA nor logB is MARK/LOG_ZERO + else{ + //1 + if(logA < 0 && logB < 0 ){ + logA = logA + log( 1 + exp(logB-logA) ); + } + //2 + else if(logA > 0 && logB > 0){ + //always make exp(.): . > 0 + //2.1 + if(logA >= logB){ + logA = -( -logA + log( 1 + exp(-logB+logA) ) ); + } + //2.2 + else{ + logA = -( -logB + log( 1 + exp(-logA+logB) ) ); + } + + } + //3 + else if(logA < 0 && logB > 0){ + //3.1 + if(logA + logB > 0){ + //logA = logA + log( 1 - exp(-logB-logA) ); + logA = -logB + log( exp(logA + logB) - 1 ); + + } + //3.2 + else if(logA + logB < 0){ + logA = -( logA + log( exp(-logB-logA) - 1) ); + } + //3.3 + else{ + logA = LOG_ZERO; + } + } + //4 + else if(logA > 0 && logB < 0){ + //4.1 + if(logA + logB > 0){ + logA = -logA + log( exp(logB+logA) - 1); + } + //4.2 + else if(logA + logB < 0){ + //logA = -( -logA + log(1 - exp(logB+logA)) ); + logA = -( logB + log( exp(-logA-logB) - 1) ); + + } + //4.3 + else{ + logA = LOG_ZERO; + } + } + + } + +} + + +//HR: Add +//Will be used for checking edges-based feasure for Engineer in BayesNW +void arcToNet(char * inFeasFileName, char* outFeasFileName){ + + cout << "inFeasFileName: " << inFeasFileName << endl; + + cout << "outFeasFileName: " << outFeasFileName << endl; + + int MAX_VAR_NUM = 32; + + int inFeatureEdges[MAX_VAR_NUM]; + for(int i = 0; i < MAX_VAR_NUM; i++){ + inFeatureEdges[i] = 0; + } + + string str; + str.assign(inFeasFileName); + + + FILE* fp = fopen(str.c_str(), "r"); + + int fromNode; + int toNode; + //int i=0; + while(!feof(fp)){ + char c; + fromNode = -1; + toNode = -1; + + fscanf(fp,"%d",&fromNode); + fscanf(fp,"%d",&toNode); + fscanf(fp,"%c",&c); + + cout << "fromNode = "<< fromNode << endl; + cout << "toNode = "<< toNode << endl; + + inFeatureEdges[toNode] += (1 << fromNode); + + + + } + fclose(fp); + int no_var = toNode; + + for(int i = 0; i < MAX_VAR_NUM; i++){ + cout << "inFeatureEdges[" << i << "] " << inFeatureEdges[i] << endl; + } + + + int outFeatureEdges[MAX_VAR_NUM]; + for(int i = 0; i < MAX_VAR_NUM; i++){ + outFeatureEdges[i] = 0; + } + + string str2; + str2.assign(outFeasFileName); + + + fp=fopen(str2.c_str(), "r"); + + //int fromNode; + //int toNode; + //int i=0; + while(!feof(fp)){ + char c; + fromNode = -1; + toNode = -1; + + fscanf(fp,"%d",&fromNode); + fscanf(fp,"%d",&toNode); + fscanf(fp,"%c",&c); + + cout << "fromNode = "<< fromNode << endl; + cout << "toNode = "<< toNode << endl; + + outFeatureEdges[toNode] += (1 << fromNode); + + + + } + fclose(fp); + //int no_var = toNode; + + + for(int i = 0; i < MAX_VAR_NUM; i++){ + cout << "outFeatureEdges[" << i << "] " << outFeatureEdges[i] << endl; + } + + + + +} + + + +//HR:F +int main(int argc, char* argv[]){ + + if (argc!=3){ + cout << "Please type the correct format: arcToNet inFeatureFileName outFeatureFileName\n" << endl; + return 1; + } + + arcToNet(argv[1], argv[2]); + +} + + + diff --git a/sourcecodes/k-best/src/arcs2dot.c b/sourcecodes/k-best/src/arcs2dot.c new file mode 100644 index 00000000..12f48158 --- /dev/null +++ b/sourcecodes/k-best/src/arcs2dot.c @@ -0,0 +1,60 @@ +#include "cfg.h" + +#include +#include +#include + +#define BUFSIZE (1024) + +void arcs2dot(char* vdfile, char* arcfile, char* dotfile) +{ + char buffer[BUFSIZE]; + FILE* dotf = strcmp(dotfile,"-") ? fopen(dotfile, "w") : stdout; + + fprintf(dotf, "digraph Bene {\n"); + fprintf(dotf, " /*@@@BEFORE_ALL@@@*/\n"); + + {/* read vdfile to get names of variables */ + + FILE* vdf = fopen(vdfile, "r"); + int v; + fprintf(dotf, " /*@@@BEFORE_NODES@@@*/\n"); + for(v=0; NULL != fgets(buffer, BUFSIZE, vdf); ++v){ + char* token = strtok(buffer, "\t\r\n"); + fprintf(dotf, " V%d [label=\"%s\" /*@@@NODE_ATTS@@@*/];\n", v, token); + } + fprintf(dotf, " /*@@@AFTER_NODES@@@*/\n"); + fclose(vdf); + } + + { + + int from, to; + FILE* arcf = strcmp(arcfile,"-") ? fopen(arcfile, "r") : stdin; + + fprintf(dotf, " /*@@@BEFORE_EDGES@@@*/\n"); + while(2 == fscanf(arcf, "%d %d", &from, &to)){ + fprintf(dotf, " V%d -> V%d /*@@@EDGE_ATTS@@@*/;\n", from, to); + } + fprintf(dotf, " /*@@@AFTER_EDGES@@@*/\n"); + + fclose(arcf); + } + + fprintf(dotf, " /*@@@AFTER_ALL@@@*/\n"); + fprintf(dotf, "}\n"); + fclose(dotf); +} + +int main(int argc, char* argv[]) +{ + + if (argc!=4) { + fprintf(stderr, "Usage: arcs2dot vdfile arcfile dotfile\n"); + return 1; + } + + arcs2dot(argv[1], argv[2], argv[3]); + + return 0; +} diff --git a/sourcecodes/k-best/src/buildk.sh b/sourcecodes/k-best/src/buildk.sh new file mode 100644 index 00000000..9bd06a3e --- /dev/null +++ b/sourcecodes/k-best/src/buildk.sh @@ -0,0 +1,23 @@ +#!/bin/bash + +if [ ${1:-UNIX} == "WIN" ]; then + CC="i586-mingw32msvc-gcc" + EXT=".exe" +else + CC=gcc +fi + +D=g++ + +#CFLAGS="-Wextra -g -ansi -pedantic" +#CFLAGS="-Wall -O3 -g -pg" +CFLAGS="-Wall -O3" + +$CC $CFLAGS -c -o files.o files.c +$CC $CFLAGS -c -o varpar.o varpar.c + +$CC $CFLAGS -o get_local_scores$EXT files.c reg.c ilogi.c ls_XIC.c ls_NML.c ls_BDe.c ls_LOO.c get_local_scores.c -lm +$CC $CFLAGS -o split_local_scores$EXT split_local_scores.c files.o +$CC $CFLAGS -o reverse_local_scores$EXT reverse_local_scores.c files.o +$D -o get_kbest_parents$EXT get_kbest_parents.cc +$D -o get_kbest_nets$EXT get_kbest_nets.cc diff --git a/sourcecodes/k-best/src/buildk_poster.sh b/sourcecodes/k-best/src/buildk_poster.sh new file mode 100644 index 00000000..a9964f8a --- /dev/null +++ b/sourcecodes/k-best/src/buildk_poster.sh @@ -0,0 +1,9 @@ +#!/bin/bash + + +D=g++ + + +$D -o get_kbest_parents$EXT get_kbest_parents.cc +$D -o get_kbest_nets$EXT get_kbest_nets.cc + diff --git a/sourcecodes/k-best/src/cfg.h b/sourcecodes/k-best/src/cfg.h new file mode 100644 index 00000000..15e3fe89 --- /dev/null +++ b/sourcecodes/k-best/src/cfg.h @@ -0,0 +1,9 @@ +#ifndef __CFG_H_ +#define __CFG_H_ + +#define MAX_NOF_VARS (32) +#define LARGEST_SET(NOF_VARS) ((NOF_VARS)==MAX_NOF_VARS?(varset_t)~0:(1U<<(NOF_VARS))-1) +typedef unsigned int varset_t; +typedef double score_t; + +#endif diff --git a/sourcecodes/k-best/src/clean.sh b/sourcecodes/k-best/src/clean.sh new file mode 100644 index 00000000..3db60eb5 --- /dev/null +++ b/sourcecodes/k-best/src/clean.sh @@ -0,0 +1,6 @@ +#!/bin/bash + +rm -f get_local_scores split_local_scores reverse_local_scores \ + get_best_parents get_best_sinks get_best_order get_best_net \ + score_net score_nets net2parents parents2arcs arcs2dot \ + *.exe *.o *~ diff --git a/sourcecodes/k-best/src/comp_pHat_D.cc b/sourcecodes/k-best/src/comp_pHat_D.cc new file mode 100644 index 00000000..fc789912 --- /dev/null +++ b/sourcecodes/k-best/src/comp_pHat_D.cc @@ -0,0 +1,213 @@ +#include +#include +#include +//#include "varpar.h" +#include +#include +#include +#include +#include +//#include "cfg.h" +#include +#include +//HR +//#define EPSILON 0.001 +#define EPSILON 0.00001 + + +//HR: Add for hash_set +#if __GNUC__ < 3 && __GNUC__ >= 2 && __GNUC_MINOR__ >= 95 +# include +# include +# define gnu_namespace std +#elif __GNUC__ >= 3 +# include +# if __GNUC_MINOR__ == 0 +# include +# define gnu_namespace std +# else +# include +# define gnu_namespace __gnu_cxx +# endif +#else +# include +# include +# define gnu_namespace std +#endif +using namespace gnu_namespace; +#include +// + + +using namespace std; + + +//HR: add from Engine.h in BayesNw +#define MARK 9e99 +#define LOG_ZERO -1e101 +//#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)); +//HR: add from UpdateHR.h in BayesNw to deal with the problem in sum = 0 in comp_post() +//call by reference to the original value +//the log sum result is stored in logA +//may change double into long double; no, use double +void logAddComp(double & logA, double & logB){ + if(logA == MARK){ + logA = logB; + } + else if (logA == LOG_ZERO){ + logA = logB; + } + else if(logB == MARK){ + logA = logA; + } + else if(logB == LOG_ZERO){ + logA = logA; + } + //Nither logA nor logB is MARK/LOG_ZERO + else{ + //1 + if(logA < 0 && logB < 0 ){ + logA = logA + log( 1 + exp(logB-logA) ); + } + //2 + else if(logA > 0 && logB > 0){ + //always make exp(.): . > 0 + //2.1 + if(logA >= logB){ + logA = -( -logA + log( 1 + exp(-logB+logA) ) ); + } + //2.2 + else{ + logA = -( -logB + log( 1 + exp(-logA+logB) ) ); + } + + } + //3 + else if(logA < 0 && logB > 0){ + //3.1 + if(logA + logB > 0){ + //logA = logA + log( 1 - exp(-logB-logA) ); + logA = -logB + log( exp(logA + logB) - 1 ); + + } + //3.2 + else if(logA + logB < 0){ + logA = -( logA + log( exp(-logB-logA) - 1) ); + } + //3.3 + else{ + logA = LOG_ZERO; + } + } + //4 + else if(logA > 0 && logB < 0){ + //4.1 + if(logA + logB > 0){ + logA = -logA + log( exp(logB+logA) - 1); + } + //4.2 + else if(logA + logB < 0){ + //logA = -( -logA + log(1 - exp(logB+logA)) ); + logA = -( logB + log( exp(-logA-logB) - 1) ); + + } + //4.3 + else{ + logA = LOG_ZERO; + } + } + + } + +} + + +//HR: Add +//Update from comp_post +void comp_pHat(){ + //double logSum = -450.491; //iris + double logSum = -9418.56; //tic + + //double logSum = -637.91; //zoo + + cerr << "logSum = " << logSum << endl; + + //double truelogP = -450.491; //iris + double truelogP = -9418.29; //tic + + //double truelogP = -619.502; //zoo + + cerr << "truelogP = " << truelogP << endl; + + int sampleSize = 1000; + + long double post[sampleSize]; + double relaLogScore[sampleSize]; + + double sumPost[sampleSize]; + + string str; + char str1[32]; + + FILE *fp = fopen("netpost" ,"r"); + int i = 0; + + + + + while(!feof(fp) && i < sampleSize){ + char c; + fscanf(fp,"%Lg",&post[i]); + //cout << " post[" << i << "] = " << post[i] << endl; + fscanf(fp,"%c",&c); + i++; + } + + for(int j = 0; j < sampleSize; j++){ + relaLogScore[j] = log(post[j]) + logSum; + //cout << " relaLogScore[" << i << "] = " << relaLogScore[j] << endl; + } + //cout << " relaLogScore[" << 1 << "] = " << relaLogScore[1] << endl; + //cout << " relaLogScore[" << 2 << "] = " << relaLogScore[2] << endl; + + for(int j = 0; j < sampleSize; j++){ + sumPost[j] = LOG_ZERO; + } + + for(int j = 0; j < sampleSize; j++){ + for(int k = 0; k <= j; k++){ + logAddComp(sumPost[j], relaLogScore[k]); + } + //cout << " sumPost[" << j << "] = " << sumPost[j] << endl; + } + + double DRatio[sampleSize]; + + + for(int j = 0; j < sampleSize; j++){ + DRatio[j] = exp( sumPost[j] - truelogP ) * 100; + //cout << "" << DRatio[j] << endl; + } + + double GRatio[sampleSize]; + for(int j = 0; j < sampleSize; j++){ + GRatio[j] = post[0]/post[j]; + //cout << "" << GRatio[j] << endl; + } + + cout << "\n\nDRatio\t GRatio" << endl; + for(int j = 0; j < sampleSize; j++){ + cout << "" << DRatio[j] << "\t " << GRatio[j] << endl; + } +} + + + +//HR:F +int main(int argc, char* argv[]){ + comp_pHat(); + +} + + + diff --git a/sourcecodes/k-best/src/data.sh b/sourcecodes/k-best/src/data.sh new file mode 100644 index 00000000..d775b92a --- /dev/null +++ b/sourcecodes/k-best/src/data.sh @@ -0,0 +1,27 @@ +#!/bin/bash + +binpath=`dirname $0` + +if [ $# -ne 4 ]; then + echo Usage: data2net.sh vdfile datafile score resultdir 1>&2 + exit 1 +fi + +vdfile=$1; shift +datafile=$1; shift +score=$1; shift +rdir=$1; shift + +mkdir -p $rdir +nof_vars=`cat $vdfile | wc -l` + +START=$(date +%s) +$binpath/get_best_parents $nof_vars ${rdir} +$binpath/get_best_sinks $nof_vars ${rdir} ${rdir}/sinks +$binpath/get_best_order $nof_vars ${rdir}/sinks ${rdir}/ord +$binpath/get_best_net $nof_vars ${rdir} ${rdir}/ord ${rdir}/net +$binpath/score_net ${rdir}/net ${rdir} +END=$(date +%s) +DIFF=$(( $END - $START )) +echo "Best network took $DIFF seconds" + diff --git a/sourcecodes/k-best/src/data2net.sh b/sourcecodes/k-best/src/data2net.sh new file mode 100644 index 00000000..3bc5d133 --- /dev/null +++ b/sourcecodes/k-best/src/data2net.sh @@ -0,0 +1,35 @@ +#!/bin/bash + +binpath=`dirname $0` + +if [ $# -ne 4 ]; then + echo Usage: data2net.sh vdfile datafile score resultdir 1>&2 + exit 1 +fi + +vdfile=$1; shift +datafile=$1; shift +score=$1; shift +rdir=$1; shift + +mkdir -p $rdir +nof_vars=`cat $vdfile | wc -l` + +START=$(date +%s) +$binpath/get_local_scores $vdfile $datafile $score 1 0 ${rdir}/res +$binpath/split_local_scores $nof_vars ${rdir} +$binpath/reverse_local_scores $nof_vars ${rdir} +END=$(date +%s) +DIFF=$(( $END - $START )) +echo "Local score took $DIFF seconds" + +START=$(date +%s) +$binpath/get_best_parents $nof_vars ${rdir} +$binpath/get_best_sinks $nof_vars ${rdir} ${rdir}/sinks +$binpath/get_best_order $nof_vars ${rdir}/sinks ${rdir}/ord +$binpath/get_best_net $nof_vars ${rdir} ${rdir}/ord ${rdir}/net +$binpath/score_net ${rdir}/net ${rdir} +END=$(date +%s) +DIFF=$(( $END - $START )) +echo "Best network took $DIFF seconds" + diff --git a/sourcecodes/k-best/src/data2netk.sh b/sourcecodes/k-best/src/data2netk.sh new file mode 100644 index 00000000..e20f12c9 --- /dev/null +++ b/sourcecodes/k-best/src/data2netk.sh @@ -0,0 +1,55 @@ +#!/bin/bash + +binpath=`dirname $0` + +if [ $# -ne 5 ]; then + echo Usage: data2net.sh vdfile datafile score resultdir k 1>&2 + exit 1 +fi + +vdfile=$1; shift +datafile=$1; shift +score=$1; shift +rdir=$1; shift +k=$1;shift + +mkdir -p $rdir +nof_vars=`cat $vdfile | wc -l` + +START=$(date +%s%N) +$binpath/get_local_scores $vdfile $datafile $score 1 0 ${rdir}/res +#HR +echo "get_local_scores is done." +$binpath/split_local_scores $nof_vars ${rdir} +#HR +echo "split_local_scores is done." +$binpath/reverse_local_scores $nof_vars ${rdir} +#HR +echo "reverse_local_scores is done." +END=$(date +%s%N) +#DIFF=$(( $END - $START )) +DIFF1=$(( $END - $START )) +echo "Local score took $DIFF1 seconds" + +#the above is the same as original bene + +START=$(date +%s%N) +$binpath/get_kbest_parents $nof_vars ${rdir} $k +#HR +END=$(date +%s%N) +#DIFF=$(( $END - $START )) +DIFF2=$(( $END - $START )) +echo "get_kbest_parents is done." +echo "get_kbest_parents took $DIFF2 seconds" + +START=$(date +%s%N) +$binpath/get_kbest_nets $nof_vars ${rdir} $k +#HR +END=$(date +%s%N) +#DIFF=$(( $END - $START )) +DIFF3=$(( $END - $START )) +echo "get_kbest_nets is done." +echo "get_kbest_nets took $DIFF3 seconds" + +DIFF=`expr $DIFF1 + $DIFF2 + $DIFF3` +echo "total process took $DIFF seconds" diff --git a/sourcecodes/k-best/src/data2netk_poster.sh b/sourcecodes/k-best/src/data2netk_poster.sh new file mode 100644 index 00000000..6b3a34ca --- /dev/null +++ b/sourcecodes/k-best/src/data2netk_poster.sh @@ -0,0 +1,56 @@ +#!/bin/bash + +binpath=`dirname $0` + +if [ $# -ne 4 ]; then + echo Usage: data2net.sh datafile resultdir k 1>&2 + exit 1 +fi + +datafile=$1; shift +#score=$1; shift +rdir=$1; shift +k=$1;shift +n=$1;shift + +nof_vars=`head -1 $datafile|wc -w` +echo "The number of variables is $nof_vars." + +#HR: +if [ "$n" -eq 0 ] +then + maxindegree=`expr $nof_vars - 1` +else + maxindegree=$n +fi + +#maxindegree=`expr $nof_vars - 1` +echo "maxindegree is $maxindegree." + +#HR: +nof_instances=`cat $datafile | wc -l` +#delete the header and type +nof_instances=`expr $nof_instances - 2` +echo "The number of instances is $nof_instances." + + +START=$(date +%s%N) +$binpath/get_kbest_parents $nof_vars ${rdir} $k +#HR +END=$(date +%s%N) +#DIFF=$(( $END - $START )) +DIFF2=$(( $END - $START )) +echo "get_kbest_parents is done." +echo "get_kbest_parents took $DIFF2 n-seconds." + +START=$(date +%s%N) +$binpath/get_kbest_nets $nof_vars ${rdir} $k +#HR +END=$(date +%s%N) +#DIFF=$(( $END - $START )) +DIFF3=$(( $END - $START )) +echo "get_kbest_nets is done." +echo "get_kbest_nets took $DIFF3 n-seconds." + +DIFF=`expr $DIFF1 + $DIFF2 + $DIFF3` +echo "The total process took $DIFF n-seconds." diff --git a/sourcecodes/k-best/src/data2netk_poster.sh.bk b/sourcecodes/k-best/src/data2netk_poster.sh.bk new file mode 100644 index 00000000..6b3a34ca --- /dev/null +++ b/sourcecodes/k-best/src/data2netk_poster.sh.bk @@ -0,0 +1,56 @@ +#!/bin/bash + +binpath=`dirname $0` + +if [ $# -ne 4 ]; then + echo Usage: data2net.sh datafile resultdir k 1>&2 + exit 1 +fi + +datafile=$1; shift +#score=$1; shift +rdir=$1; shift +k=$1;shift +n=$1;shift + +nof_vars=`head -1 $datafile|wc -w` +echo "The number of variables is $nof_vars." + +#HR: +if [ "$n" -eq 0 ] +then + maxindegree=`expr $nof_vars - 1` +else + maxindegree=$n +fi + +#maxindegree=`expr $nof_vars - 1` +echo "maxindegree is $maxindegree." + +#HR: +nof_instances=`cat $datafile | wc -l` +#delete the header and type +nof_instances=`expr $nof_instances - 2` +echo "The number of instances is $nof_instances." + + +START=$(date +%s%N) +$binpath/get_kbest_parents $nof_vars ${rdir} $k +#HR +END=$(date +%s%N) +#DIFF=$(( $END - $START )) +DIFF2=$(( $END - $START )) +echo "get_kbest_parents is done." +echo "get_kbest_parents took $DIFF2 n-seconds." + +START=$(date +%s%N) +$binpath/get_kbest_nets $nof_vars ${rdir} $k +#HR +END=$(date +%s%N) +#DIFF=$(( $END - $START )) +DIFF3=$(( $END - $START )) +echo "get_kbest_nets is done." +echo "get_kbest_nets took $DIFF3 n-seconds." + +DIFF=`expr $DIFF1 + $DIFF2 + $DIFF3` +echo "The total process took $DIFF n-seconds." diff --git a/sourcecodes/k-best/src/datak.sh b/sourcecodes/k-best/src/datak.sh new file mode 100644 index 00000000..8ae66439 --- /dev/null +++ b/sourcecodes/k-best/src/datak.sh @@ -0,0 +1,23 @@ +#!/bin/bash + +binpath=`dirname $0` + +if [ $# -ne 5 ]; then + echo Usage: data2net.sh vdfile datafile score resultdir k 1>&2 + exit 1 +fi + +vdfile=$1; shift +datafile=$1; shift +score=$1; shift +rdir=$1; shift +k=$1; shift + +mkdir -p $rdir +nof_vars=`cat $vdfile | wc -l` +START=$(date +%s) +$binpath/get_kbest_parents $nof_vars ${rdir} $k +$binpath/get_kbest_nets $nof_vars ${rdir} $k +END=$(date +%s) +DIFF=$(( $END - $START )) +echo "Top k took $DIFF seconds" diff --git a/sourcecodes/k-best/src/drawgraph.sh b/sourcecodes/k-best/src/drawgraph.sh new file mode 100644 index 00000000..a0f59c88 --- /dev/null +++ b/sourcecodes/k-best/src/drawgraph.sh @@ -0,0 +1,8 @@ +#!/bin/bash +k=$1 +rdir=$2 +for (( i = 0 ; i < k; i++ )) +do + ${rdir}/../../src/arcs2dot ${rdir}/../iris.vd ${rdir}/$arc$i ${rdir}/dot_$i +done + diff --git a/sourcecodes/k-best/src/files.c b/sourcecodes/k-best/src/files.c new file mode 100644 index 00000000..3227f63c --- /dev/null +++ b/sourcecodes/k-best/src/files.c @@ -0,0 +1,59 @@ +#include +#include +#include + +#include "files.h" + +int nof_digits(int n) +{ + int c=0; + do { + c += 1; + n /= 10; + } while(n); + return c; +} + +int nof_lines(char* filename) { + int nof_lines = 0; + FILE* f = fopen(filename,"r"); + int c = 0; + while(EOF != (c=fgetc(f))) + if(c=='\n') ++nof_lines; + fclose(f); + return nof_lines; +} + +char* create_fn(char* dirname, int i, char* ext) +{ char* fn = malloc((strlen(dirname)+1+nof_digits(i)+strlen(ext)+1) + *sizeof(char)); + sprintf(fn,"%s/%d%s", dirname,i,ext); + return fn; +} + +FILE* open_file(char* dirname, int i, char* ext, char* mode) +{ + FILE* f; + char* fn = create_fn(dirname, i, ext); + sprintf(fn,"%s/%d%s", dirname,i,ext); + f = fopen(fn, mode); + free(fn); + return f; +} + +FILE** open_files(int nof_vars, char* dirname, char* ext, char* mode) +{ + FILE** files = malloc(nof_vars * sizeof(FILE*)); + int i; + for(i=0; i + +int nof_lines(char* filename); +char* create_fn(char* dirname, int i, char* ext); +FILE* open_file(char* dirname, int i, char* ext, char* mode); +FILE** open_files(int nof_vars, char* dirname, char* ext, char* mode); +void free_files(int nof_vars, FILE** files); + +#endif diff --git a/sourcecodes/k-best/src/files.o b/sourcecodes/k-best/src/files.o new file mode 100644 index 00000000..33f376d3 Binary files /dev/null and b/sourcecodes/k-best/src/files.o differ diff --git a/sourcecodes/k-best/src/get_best_net.c b/sourcecodes/k-best/src/get_best_net.c new file mode 100644 index 00000000..5ed6ebfb --- /dev/null +++ b/sourcecodes/k-best/src/get_best_net.c @@ -0,0 +1,63 @@ +#include "cfg.h" + +#include +#include +#include + +#include "files.h" +#include "varpar.h" + +void get_best_net(int nof_vars, char* ordfile, FILE** bpsfiles, varset_t* net) +{ + FILE* ordf = fopen(ordfile,"r"); + varset_t parcands = 0; + int i; + for(i=0;i +#include +#include + +void get_best_order(int nof_vars, char* sinkfile, char* ord) +{ + FILE* sinkf = fopen(sinkfile, "rb"); + varset_t set = LARGEST_SET(nof_vars); + // printf("\n\n\nset = %u",set); + int i; + for(i=nof_vars-1; i>=0;--i){ + char sink; + // printf("\n\n\n\nset =%u",set); + fseek(sinkf, set, SEEK_SET); + + fread(&sink, sizeof(char), 1, sinkf); +// printf(" sink = %d",sink); + ord[i] = sink; + + set ^= 1U< +#include +#include + +#include "files.h" + + +void get_best_parents(int nof_vars, char* dirname) { +//printf("I'm here!!!"); + + varset_t nof_parsets = 1U<<(nof_vars-1); + int i; +// FILE* fp; + score_t* scores = malloc(nof_parsets*sizeof(score_t)); + varset_t* bsps = malloc(nof_parsets*sizeof(varset_t)); +// fp = fopen("/home/grad/lram/RA stuff/check/src/out.txt","w"); + + for(i=0;i= scores[ps]){ + scores[ps] = scores[subset]; + bsps[ps] = bsps[subset]; + } + } + } + } + } + + + { /* WRITE BEST PARENTS FOR i */ + + FILE* fout = open_file(dirname, i, ".bps", "wb"); + fwrite(bsps, sizeof(varset_t), nof_parsets, fout); + +//fprintf(fp,"\n\nBest parent sets for var %d",i); + //for(j=0;j +#include +#include + +#include "files.h" + +void get_best_sinks(int nof_vars, FILE** files, char* sinkfile) +{ + + FILE* sinkf = fopen(sinkfile, "wb"); +FILE *fp = fopen("/home/grad/lram/RA stuff/check/src/sinks.txt","w"); + +score_t* scores = malloc((1U< +#include +#include +//#include "varpar.h" +#include +#include +#include +#include +#include +#include "cfg.h" +#include +#include +//HR +//#define EPSILON 0.001 +#define EPSILON 0.00001 + + +//HR: Add for hash_set +#if __GNUC__ < 3 && __GNUC__ >= 2 && __GNUC_MINOR__ >= 95 +# include +# include +# define gnu_namespace std +#elif __GNUC__ >= 3 +# include +# if __GNUC_MINOR__ == 0 +# include +# define gnu_namespace std +# else +# include +# define gnu_namespace __gnu_cxx +# endif +#else +# include +# include +# define gnu_namespace std +#endif +using namespace gnu_namespace; +#include +// + + +using namespace std; + + +#define MARK 9e99 +#define LOG_ZERO -1e101 + +//HR: add from UpdateHR.h to deal with the problem in sum = 0 in comp_post() +//call by reference to the original value +//the log sum result is stored in logA + +#define epsilon 1.0e5 + +bool logSpecValEquals(double logX, double logSpecVal){ + return (logX >= logSpecVal - epsilon && logX <= logSpecVal + epsilon); +} + +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; + } + + + //Neither 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)c + + + +//HR: Note: The whole program supposes that max no. of vars = 32 and hardcord 32 + + +//Convert the parent set to the appropriate variable set +//HR:F +varset_t parset2varset(int v, varset_t set){ + varset_t sinkleton = 1U<(net_set s1,net_set s2); + friend bool same_obj(net_set s1, net_set s2, int nof_vars); + +}; + +//HR:F +/*================================================================*/ + +int operator<(net_set s1,net_set s2){ + if(fabs(s1.score-s2.score)(net_set s1,net_set s2){ + if(fabs(s1.score-s2.score)s2.score); +} + +int operator==(net_set s1,net_set s2){ + return(fabs(s1.score-s2.score)* v,net_set s, int k, int nof_vars){ + list::iterator it; + //HR: can improve its time + //first check whether v->size() > 0 && s < v->back() + //see Lav's comment + + + bool flag=false; + + + //HR: Refined + // + if( (v->size() > 0) && (s < v->back()) ){ + v->push_back(s); + flag=true; + } + else{ + // + + for(it=v->begin();it!=v->end();it++){ + if(same_obj(s,*it,nof_vars)){ +// cout<<"Same object so not inserted!!!!!!!!!!!!!"<(*it)){ + // cout<<"score of s:"<insert(it,s); + flag=true; + break; + } + + } + + } + + if(flag==false){ + v->push_back(s); + } + if(v->size()>k){ + v->pop_back(); + } + + //v->unique(); +} + +/*================================================================*/ +//class for bparents +//HR:F +//HR: (variable set, score) pair +class score_network{ + public: + score_t scores; + varset_t vset; + + friend int operator<(score_network s1,score_network s2); + friend int operator==(score_network s1,score_network s2); + friend int operator>(score_network s1,score_network s2); + friend bool same_obj(score_network s1, score_network s2); + +}; + + +//HR:F +int operator<(score_network s1,score_network s2){ + if(fabs(s1.scores-s2.scores)(score_network s1,score_network s2){ + if(fabs(s1.scores-s2.scores)s2.scores); + +} + +int operator==(score_network s1,score_network s2){ + return(fabs(s1.scores-s2.scores)* v,score_network s, int k) + { + list::iterator it; + + if( s < v->back()) + { + v->push_back(s); + } + else + { + for(it=v->begin();it!=v->end();it++) + { + // if(same_obj(s,*it)) + // { + // break; + // } +// else + if(s>(*it)) + { + v->insert(it,s); + break; + } + + } + } + + if(v->size()>k) + { + v->pop_back(); + } + + } + +/*================================================================*/ +//class for node +//HR:F +class node{ + public: + int ip; + int in; + score_t score; + + node(int x, int y, score_t z){ + ip=x; + in=y; + score=z; + } + + friend int operator<(node s1,node s2); + friend int operator==(node s1,node s2); + friend int operator>(node s1,node s2); + friend bool same_obj(node s1, node s2); +}; + +/*================================================================*/ + +int operator<(node s1,node s2){ + if(fabs(s1.score-s2.score)(node s1,node s2){ + + if(fabs(s1.score-s2.score)s2.score); +} + + +int operator==(node s1,node s2){ + return(fabs(s1.score-s2.score) v[],int size,int nof_vars){ + + list::iterator iter; + + for(int j=0;j* v){ + + vector::iterator iter; + + for(iter=v->begin();iter!=v->end();iter++){ + cout<<"bparents "<<(*iter).vset<<","<<(*iter).scores<<" ,, "; + } + +} + +void print_queue(list* v){ + + list::iterator iter; + cout<<"fringe"<begin();iter!=v->end();iter++){ + cout<<" ip "<<(*iter).ip<<", in "<<(*iter).in<<" score "<<(*iter).score; + } + +} + +void print_queue(vector* v,int nof_vars){ + + vector::iterator iter; + + for(iter=v->begin();iter!=v->end();iter++){ +// cout<<" "<<(*iter).bsps<<","<<(*iter).scores; + cout<<"bnets "<<" Score:"<<(*iter).score<<" net "<* v,int nof_vars){ + + list::iterator iter; + + for(iter=v->begin();iter!=v->end();iter++){ +// cout<<" "<<(*iter).bsps<<","<<(*iter).scores; + cout<<"kbnet "<<" Score:"<<(*iter).score<<" net "<* v,node s, int k){ + list::iterator it; + + //HR: test + //bool sameObj = false; + //bool inserted = false; + // + + //if( s < v->back()){ same + if( (v->size() == 0) || (s < v->back())){ + v->push_back(s); + } + else{ + for(it=v->begin();it!=v->end();it++){ + //comment it out will give the different answer. Because Lav's original is wrong, it will still insert the same object. It will eventually pop the correct end > k + // + if(same_obj(s, *it)){ + break; +// sameObj = true; +// cout << "s: "; +// print_node(s); +// cout << "*it: "; +// print_node(*it); +// cout << endl; +// print_queue(v); + + } + else{ + // + if(s>(*it)){ + v->insert(it,s); + //inserted = true; + break; + } + } + + } + + //HR: test +// if(sameObj == true && inserted == true){ +// cout << "Error: sameObj == true && inserted = true" << endl; + +// cout << "s: "; +// print_node(s); +// print_queue(v); +// } + } + + if(v->size()>k){ + //cout << "\n\nv->size()>k happens, v->pop_back()\n" << endl; + v->pop_back(); + } + // cout<<"Inserted:"< + + struct hash { + + hash hasher_ns; + + size_t operator()(const net_set* ns) const { + return hasher_ns((unsigned) roundl(- ns->score*10)); // wrong hash function?? + } + + }; //end of struct + +} //end of namespace + + + +struct eqNetScore{ + bool operator()(const net_set* ns1, const net_set* ns2) const{ + return same_obj(*ns1, *ns2); + } +}; + + + +typedef hash_set, eqNetScore> NetScoreHashSet; + +NetScoreHashSet knets; + + +//HR: Add +void print_nsHashSet(NetScoreHashSet nsHSet, int nof_vars){ + for(NetScoreHashSet::iterator it = nsHSet.begin(); it != nsHSet.end(); it++){ + print_netscore(**it, nof_vars); + } +} +// + + + +//=====================================================================================// + +/*================================================================*/ + +//HR:F +//HR: call gettopk(&bnets,&bparents,&kbnetscore[varset],k,sink,nof_vars); +void gettopk(vector* bnets, vector* bparents, list* kbnet, int k, varset_t sink,int nof_vars){ + int h; + int ipmax=bparents->size(); +// cout<<"ipmax: "<size(); + // cout<<"inmax: "< fringe; + + //HR: different from Algo4, here start from 0 vs. 1 + node n(0,0,(bparents->at(0)).scores + (bnets->at(0)).score); + insert_vec(&fringe,n,k); + // print_queue(&fringe); + // scanf("%d",&h); + + NetScoreHashSet::iterator iter1; + + while(!fringe.empty()){ + // cout<<"Fringe not empty"<size())at(in)).net[j]; + // cout<<"Net w j"<at(ip)).vset); + //cout<<"net w sinkleton"<back()).score)>EPSILON)&&(score>(kbnet->back().score))){ + // cout<<"score>kbnetback's score"<at(in)).net[j]; + // cout<<"Net w j"<at(ip)).vset; + //cout<<"net w sinkleton"<back(); + + //Lav + insert_vec(kbnet,*netwP,k,nof_vars); + //cout<<"After inserting"<at(ip+1)).scores + (bnets->at(in)).score); + //cout<<"after child"<at(ip)).scores + (bnets->at(in+1)).score); + + //may add hash_set generated check here + + insert_vec(&fringe,child,k); + //cout<<"Fringe after inserting"<* kbnet, int k,int nof_vars){ + list::iterator it; + //int j=0; + //each it is one of k best + for(it=kbnet->begin();it!=kbnet->end();it++){ + for(int p=0;p kbnetscore[],char* dirname,int k,int nof_vars) +{ + string str; + list::iterator it; + char str1[10]; + int j=0; + //HR: each it is the one of the k best + for(it=kbnetscore[(1U<<(nof_vars))-1].begin();it!=kbnetscore[(1U<<(nof_vars))-1].end();it++){ + str.assign(dirname); + str.append("/"); + sprintf(str1,"%d",j++); + str.append(str1); + str.append("net"); + FILE* netf=fopen(str.c_str(),"w"); + fprintf(netf,"%f \n",(*it).score); + // cout<<(*it).score< kbnetscore[],char* dirname,int k,int nof_vars){ +// string str; + list::iterator it; + + int j=0; + for(j=0;j<((1U<* kbnet, int k, char* dirname, int nof_vars){ + + //cout << "Start comp_postExact_compl()" << endl; + //long double Sum=0.0; + //double logSum = LOG_ZERO; + + list::iterator it; + + //HR: assume k < 5000 + long double postExact[5000]; + + string str; + char str1[32]; + + //cout << "dirname = " << dirname << endl; + str.assign(dirname); + str.append("/"); + str.append("exactPD.txt"); + //str.assign("exactPD.txt"); + + FILE *fp; + fp=fopen(str.c_str(),"r"); + double exactPD = 0.0; + + if(!feof(fp)){ + char c; + fscanf(fp, "%lf", &exactPD); + fscanf(fp,"%c",&c); + } + fclose(fp); + //printf("exactPD = %18.8f\n", exactPD); + + + str.assign(dirname); + str.append("/netpostExact"); + fp=fopen(str.c_str(),"w"); + + //for each it of k best networks +// for(it=kbnet->begin();it!=kbnet->end();it++){ +// //HR: change it +// //sum+=exp((*it).score); +// logAddComp(logSum, (*it).score); //HR +// } +// cout << "logSum: " << logSum << endl; //HR +// printf("logSum = %18.6f\n", logSum); +// // +// long double Sum = exp(logSum); //HR +// cout << "Sum: " << exp(logSum) << endl; //HR + // + + int i=0; + //for each it of k best networks, compute equation (4) + //each post[i] ( 0 <= i <= k - 1). repre post prob for i'th best network + for(it=kbnet->begin();it!=kbnet->end();it++){ + //post[i]=(exp((*it).score))/sum; + //double logScore = (*it).score; + //LOGMINUS_NEW(logScore, logSum) + postExact[i]=exp( (*it).score - exactPD ); // + fprintf(fp,"%Lg\n",postExact[i]); + i++; + } + fclose(fp); + + + //Start the update + + //define and init + long double postProbExact[nof_vars][nof_vars]; + for(int i = 0; i < nof_vars; i++){ + for (int j = 0; j < nof_vars; j++){ + postProbExact[i][j] = 0.0; + } + } + + int listSize = kbnet->size(); + int num1 = -1; + int num2 = -1; + + + //HR: for each j of k best networks, j >=1 + for(int j=0; j " << i << " \t" << postProbExact[j][i] << endl; + } + } + +} //end void comp_postExact_compl() + + + + +//HR: Add +//HR: Update from comp_post +//Compute each edge out of nof_vars * nof_vars, instead of just the best network +void comp_post_compl(list* kbnet, int k, char* dirname, int nof_vars){ + //long double Sum=0.0; + double logSum = LOG_ZERO; + list::iterator it; + //HR: assume k < 5000 + long double post[5000]; + + string str; + char str1[32]; + str.assign(dirname); + str.append("/netpost"); + FILE *fp=fopen(str.c_str(),"w"); + + //for each it of k best networks + for(it=kbnet->begin();it!=kbnet->end();it++){ + //HR: change it + //sum+=exp((*it).score); + logAddComp(logSum, (*it).score); //HR + } + //cout << "logSum: " << logSum << endl; //HR + //printf("logSum = %18.6f\n", logSum); + // + long double Sum = exp(logSum); //HR + //cout << "Sum: " << exp(logSum) << endl; //HR + // + + int i=0; + //for each it of k best networks, compute equation (4) + //each post[i] ( 0 <= i <= k - 1). repre post prob for i'th best network + for(it=kbnet->begin();it!=kbnet->end();it++){ + //post[i]=(exp((*it).score))/sum; + //double logScore = (*it).score; + //LOGMINUS_NEW(logScore, logSum) + post[i]=exp( (*it).score - logSum ); // + fprintf(fp,"%Lg\n",post[i]); + i++; + } + fclose(fp); + + + //Start the update + + //define and init + long double postProb[nof_vars][nof_vars]; + for(int i = 0; i < nof_vars; i++){ + for (int j = 0; j < nof_vars; j++){ + postProb[i][j] = 0.0; + } + } + + int listSize = kbnet->size(); + int num1 = -1; + int num2 = -1; + + + //HR: for each j of k best networks, j >=1 + for(int j=0; j " << i << " \t" << postProb[j][i] << endl; + } + } + +} + + + +//Computation of posterior probabilities - This function's call has been commented. Remove the commenting slashes if the posterior probability needs to be computed. +//HR: call comp_post(&kbnetscore[nof_combs-1],k,dirname); + +void comp_post(list* kbnet, int k, char* dirname){ + //long double Sum=0.0; + double logSum = LOG_ZERO; + list::iterator it; + long double post[5000]; + + string str; + char str1[32]; + str.assign(dirname); + str.append("/netpost"); + FILE *fp=fopen(str.c_str(),"w"); + + //for each it of k best networks + for(it=kbnet->begin();it!=kbnet->end();it++){ + //HR: change it + //sum+=exp((*it).score); + logAddComp(logSum, (*it).score); //HR + } + cout << "logSum: " << logSum << endl; //HR + // + long double Sum = exp(logSum); //HR + cout << "Sum: " << exp(logSum) << endl; //HR + // + + int i=0; + //for each it of k best networks, compute equation (4) + //each post[i] ( 0 <= i <= k - 1). repre post prob for i'th best network + for(it=kbnet->begin();it!=kbnet->end();it++){ + //post[i]=(exp((*it).score))/sum; + //double logScore = (*it).score; + //LOGMINUS_NEW(logScore, logSum) + post[i]=exp( (*it).score - logSum ); // + fprintf(fp,"%Lg\n",post[i]); + i++; + } + fclose(fp); + + + ifstream indata; + int n=0; + int num[50][2]; + int num1,num2; + // char f[50][20]; + long double postf[50]; + str.assign(dirname); + str.append("/"); + str.append("arc0"); + fp=fopen(str.c_str(),"r"); + // fp=fopen(str.c_str(),"r"); + i=0; + while(!feof(fp)){ + char c; + //HR: for parent + fscanf(fp,"%d",&num[i][0]); + //HR: for var + fscanf(fp,"%d",&num[i][1]); + fscanf(fp,"%c",&c); + + //indata>>num[i][1]; + //HR: + // cout<<"HIIIII"<=1 + for(j=1;j>num1; + // indata>>num2; + + //HR: + // cout< kbnetscore[nof_combs]" << endl; + //cerr << "before list kbnetscore[nof_combs]" << endl; + //HR: net_set is the (net, score) pair + //HR: have not yet set each list kbnetscore[W] the size k + list kbnetscore[nof_combs]; + //cerr << "after list kbnetscore[nof_combs]" << endl; + //cout << "after list kbnetscore[nof_combs]" << endl; + list::iterator it; + + vector bnets; + + //HR: score_network is (variable set, score) pair + score_network s; + + vector bparents; + + //HR: Add hash_set knets here; global may be better so that gettopk can access + //NetScoreHashSet knets; + // + //cout << "before net_set netw" << endl; + + net_set netw; + netw.score=0; + for(int j=0;j v[],int size,int nof_vars) + //cout << "print_queue(kbnetscore)" << endl; + //print_queue(kbnetscore,(1U< +#include +#include +#include "cfg.h" + +#include +#include +#include +#include +#define EPSILON 0.001 + + +#include "files.h" + +using namespace std; + +int h; + +//class for bparents +//HR:F +//HR: (variable set, score) pair +//This class is to store the best parent set and the corresponding score +class score_network{ + public: + score_t scores; + varset_t bsps; + + score_network(score_t a,varset_t b){ + scores=a; + bsps=b; + } + + score_network(){ + scores=bsps=0; + } + + friend int operator<(score_network s1,score_network s2); + friend int operator==(score_network s1,score_network s2); + friend int operator>(score_network s1,score_network s2); + friend bool same_obj(score_network s1, score_network s2); + + +}; + +//HR:F +/*================================================================*/ + +void print_queue(list* v){ + + list::iterator iter; + + for(iter=v->begin();iter!=v->end();iter++){ + cout<<" "<<(*iter).bsps<<","<<(*iter).scores; + } + +} + +//HR:F +/*================================================================*/ + +int operator<(score_network s1,score_network s2){ + if(fabs(s1.scores-s2.scores)(score_network s1,score_network s2){ + + if(fabs(s1.scores-s2.scores)s2.scores); + +} + +int operator==(score_network s1,score_network s2){ + return(fabs(s1.scores-s2.scores)* v,score_network s, int k){ + list::iterator it; + + if( s < v->back()){ + v->push_back(s); + } + else{ + for(it=v->begin();it!=v->end();it++){ + if(same_obj(s,*it)){ + break; + } + else if(s>(*it)){ + v->insert(it,s); + break; + } + + } + } + + if(v->size()>k){ + v->pop_back(); + } + +} + +/*================================================================*/ + +void get_best_parents(int nof_vars, char* dirname, int k){ + varset_t nof_parsets = 1U<<(nof_vars-1); + int i; +// cout< v[nof_parsets]; + list::iterator it; + score_network s; + string str; + char str1[40]; + + //HR: For each var + for(i=0;iscores), sizeof(double), 1, fout); + fwrite(&(it->bsps),sizeof(unsigned int),1,fout); + //print_queue(&v[j]); + //printf("\n"); + //fprintf(fout,"\n"); + } + + fwrite(&sentinel,sizeof(double),1,fout); + v[j].clear(); + + } + //printf("\n"); + fclose(fout); + + }// endfor(i=0;i /dev/null; then + ${rdir}/../../src/net2parents ${rdir}/${i}net ${rdir}/${i}parent \ + | ${rdir}/../../src/parents2arcs ${rdir}/${i}parent ${rdir}/${i}arc \ + | ${rdir}/../../src/arcs2dot ${rdir}/../iris.vd ${rdir}/${i}arc ${rdir}/dot_$i \ + | dot ${rdir}/dot_$i -Tps -o ${rdir}/dot_res_${i}.ps + echo See $d/resdir/iris.ps for a postscript picture of the net. +#fi +done \ No newline at end of file diff --git a/sourcecodes/k-best/src/ilogi.c b/sourcecodes/k-best/src/ilogi.c new file mode 100644 index 00000000..50a8e651 --- /dev/null +++ b/sourcecodes/k-best/src/ilogi.c @@ -0,0 +1,37 @@ +#include +#include +#include "cfg.h" +#include "ilogi.h" + +score_t* ilogi = NULL; +int len_ilogi = 0; + +static int refcount = 0; + +int ensure_ilogi(int maxi){ + if (maxi >= len_ilogi) + len_ilogi = maxi+1; + + ilogi = (score_t*) realloc(ilogi, len_ilogi * sizeof(score_t)); + ++refcount; + + { + int i; + + if (ilogi) + ilogi[0] = 0.0; + + for(i=1;i +#include +#include +#include +#include "get_local_scores.h" +#include "ilogi.h" +#include "ls_BDe.h" + +extern int* nof_vals; +extern int N; + +score_t ess; + +int* freq2mem = NULL; +int* freq2mem2 = NULL; +int nof_freq2s; + + +#define BIG_BDE_DATA (1<<10) + +extern double lgamma(double); /* to avoid the warning with ansi-flag */ + +/* BDe */ + +score_t big_bde_score(int i, varset_t psi, int nof_freqs){ + + int vc_v = nof_vals[i]; + score_t pcc = get_nof_cfgs(psi); + score_t ess_per_pcc = ess / pcc; + score_t ess_per_cc = ess_per_pcc / vc_v; + + int* freqp = freqmem; + int* end_freqp = freqp + nof_freqs * vc_v; + + score_t res = nof_freqs * lgamma(ess_per_pcc); + int nof_zeros = 0; + + memset(freq2mem, 0, nof_freq2s*sizeof(int)); /* could reset later */ + memset(freq2mem2, 0, nof_freq2s*sizeof(int)); /* could reset later */ + + for(;freqp < end_freqp; freqp += vc_v) { + int pcfreq = 0; + int v; + for(v = 0; v 0.0 ? 0.0 : res; +} + + +score_t bde_score(int i, varset_t psi, int nof_freqs){ +//HR +//printf ("bde score %d, %u %d",i, psi, nof_freqs); + + int vc_v = nof_vals[i]; + score_t pcc = get_nof_cfgs(psi); + score_t ess_per_pcc = ess / pcc; + score_t ess_per_cc = ess_per_pcc / vc_v; + + int* freqp = freqmem; + int* end_freqp = freqp + nof_freqs * vc_v; + + score_t res = nof_freqs * lgamma(ess_per_pcc); + int nof_zeros = 0; + + memset(freq2mem, 0, nof_freq2s*sizeof(int)); /* could reset later */ + memset(freq2mem2, 0, nof_freq2s*sizeof(int)); /* could reset later */ + + for(;freqp < end_freqp; freqp += vc_v) { + int pcfreq = 0; + int v; + for(v = 0; v 0.0 ? 0.0 : res; +} + + + +scorefun init_BDe_scorer(char* arg){ + + scorefun sf; + + ess = atof(arg); + + if (N +#include +#include +#include "get_local_scores.h" +#include "ls_LOO.h" + +extern int* nof_vals; +extern int N; + +static score_t ess; +static int splitess; + +static int* freq2mem = NULL; +static int* freq2mem2 = NULL; +static int nof_freq2s; + + +#define BIG_LOO_DATA (1<<10) + +/* LOO */ + +score_t big_loo_score(int i, varset_t psi, int nof_freqs){ + + int vc_v = nof_vals[i]; + score_t pcc = get_nof_cfgs(psi); + score_t ess_per_pcc = splitess ? (ess / pcc) : (ess*vc_v); + score_t ess_per_cc = splitess ? (ess_per_pcc / vc_v) : ess; + + int* freqp = freqmem; + int* end_freqp = freqp + nof_freqs * vc_v; + + score_t res = - log(N); + + memset(freq2mem, 0, nof_freq2s*sizeof(int)); /* could reset later */ + memset(freq2mem2, 0, nof_freq2s*sizeof(int)); /* could reset later */ + + for(;freqp < end_freqp; freqp += vc_v) { + int pcfreq = 0; + int v; + for(v = 0; v 0.0 ? 0.0 : res; +} + + +score_t loo_score(int i, varset_t psi, int nof_freqs){ + + int vc_v = nof_vals[i]; + score_t pcc = get_nof_cfgs(psi); + + score_t ess_per_pcc = splitess ? (ess / pcc) : (ess*vc_v); + score_t ess_per_cc = splitess ? (ess_per_pcc / vc_v) : ess; + + int* freqp = freqmem; + int* end_freqp = freqp + nof_freqs * vc_v; + + score_t res = -log(N); + + memset(freq2mem, 0, nof_freq2s*sizeof(int)); /* could reset later */ + memset(freq2mem2, 0, nof_freq2s*sizeof(int)); /* could reset later */ + + for(;freqp < end_freqp; freqp += vc_v) { + int pcfreq = 0; + int v; + for(v = 0; v 0.0 ? 0.0 : res; +} + + + +scorefun init_LOO_scorer(char* arg){ + + scorefun sf; + int arglen = strlen(arg); + + ess = atof(arg); + splitess = arg[arglen-1] == 'l'; + + if (N +#include "get_local_scores.h" +#include "ilogi.h" +#include "reg.h" +#include "ls_NML.h" + +extern int* nof_vals; +extern int N; + +#define BIG_NML_DATA (1<<16) + +/* NML = log(ML) - log(REGRET) */ + +score_t big_nml_score(int i, varset_t psi, int nof_freqs){ + + int vc_v = nof_vals[i]; + int* freqp = freqmem; + int* end_freqp = freqp + nof_freqs * vc_v; + + score_t res = 0.0; + + for(;freqp < end_freqp; freqp += vc_v) { + int pcfreq = 0; + int v; + for(v = 0; v +#include +#include "get_local_scores.h" +#include "ilogi.h" +#include "ls_XIC.h" + +extern int* nof_vals; +extern int N; + +static int bic = 0; + +#define BIG_XIC_DATA (1<<16) + +/* XIC = log(ML) - 0.5*nof_params*log(n) */ + +score_t big_xic_score(int i, varset_t psi, int nof_freqs){ + + int vc_v = nof_vals[i]; + score_t pcc = get_nof_cfgs(psi); + score_t nof_params = pcc * (vc_v-1); + int* freqp = freqmem; + int* end_freqp = freqp + nof_freqs * vc_v; + + score_t res = bic ? (-0.5 * nof_params * log(N)) : nof_params; + + for(;freqp < end_freqp; freqp += vc_v) { + int pcfreq = 0; + int v; + for(v = 0; v +#include +#include"math.h" + +#include"Arguments.h" +#include"Model.h" +#include"Engine.h" + +void welcome(){ + fprintf(stderr, + " ~~~ Welcome ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"); + fprintf(stderr, + " Poster: The Tool for Computing Posterior \n"); + fprintf(stderr, + " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"); + +} +void goodbye(){ + fprintf(stderr, + " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"); +} + + + +int main(int argc, char **argv){ + welcome(); + + Arguments::init(argc, argv); + + Model model; + + model.init(); + + //HR + //model.makeADTree(); + ////model.MakeADContiTable(); + //model.testHR_ADtree(); + //model.freeADTree(); + //return 1; + + //cerr << "Start Engine.h " << endl; + + Engine engine; + + engine.init(&model); + + engine.compute_edge_probabilities(); + + goodbye(); + return 1; +} diff --git a/sourcecodes/k-best/src/net2parents.c b/sourcecodes/k-best/src/net2parents.c new file mode 100644 index 00000000..f12e9f64 --- /dev/null +++ b/sourcecodes/k-best/src/net2parents.c @@ -0,0 +1,50 @@ +#include "cfg.h" + +#include +#include +#include + +void net2parents(char* netfile, char* parfile) +{ + int v; + varset_t varset; +//score_t s; +//FILE *fp=fopen("/home/grad/lram/RA stuff/try/src/kd.log","w+"); + + + FILE* netf = strcmp(netfile,"-") ? fopen(netfile, "r") : stdin; + FILE* parf = strcmp(parfile,"-") ? fopen(parfile, "w") : stdout; +//fprintf(fp, "great one was here"); +//fclose(fp); + +//fscanf(netf,"%f",&s); + for(v=0; 1 == fscanf(netf, "%u", &varset); ++v){ + int p; +//printf("Var %d", v); + + for(p=0; varset; ++p, varset>>=1){ + if(varset & 1) { + fprintf(parf, "%d", p); + if(varset>1)fprintf(parf, " "); + } + } +// printf("Var %d", v); + fprintf(parf,"\n"); + } +//printf("Here"); + fclose(parf); + fclose(netf); +} + +int main(int argc, char* argv[]) +{ + + if (argc!=3) { + fprintf(stderr, "Usage: net2parents netfile parentfile\n"); + return 1; + } + + net2parents(argv[1], argv[2]); + + return 0; +} diff --git a/sourcecodes/k-best/src/parents2arcs.c b/sourcecodes/k-best/src/parents2arcs.c new file mode 100644 index 00000000..81e8b3b1 --- /dev/null +++ b/sourcecodes/k-best/src/parents2arcs.c @@ -0,0 +1,37 @@ +#include "cfg.h" + +#include +#include +#include + +#define BUFSIZE (1024) + +void parents2arcs(char* parfile, char* arcfile) +{ + char buffer[BUFSIZE]; + int v; + FILE* parf = strcmp(parfile,"-") ? fopen(parfile, "r") : stdin; + FILE* arcf = strcmp(arcfile,"-") ? fopen(arcfile, "w") : stdout; + for(v=0; NULL != fgets(buffer, BUFSIZE, parf); ++v){ + /* could check that the last character in buffer is newline */ + char* token; + for(token=strtok(buffer, " \t\r\n"); token; token=strtok(NULL, " \t\r\n")){ + fprintf(arcf, "%d %d\n", atoi(token), v); + } + } + fclose(arcf); + fclose(parf); +} + +int main(int argc, char* argv[]) +{ + + if (argc!=3) { + fprintf(stderr, "Usage: parents2arcs parentfile arcfile\n"); + return 1; + } + + parents2arcs(argv[1], argv[2]); + + return 0; +} diff --git a/sourcecodes/k-best/src/reg.c b/sourcecodes/k-best/src/reg.c new file mode 100644 index 00000000..311bf59c --- /dev/null +++ b/sourcecodes/k-best/src/reg.c @@ -0,0 +1,320 @@ +#include +#include "reg.h" + +#define pi 3.1415926535897932384626433832795 + + +static score_t r2_1000[] = { +1.0, 2.0, 2.5, 2.88888888889, 3.21875, 3.5104, 3.77469135802, +4.01813870071, 4.24501800537, 4.45831574489, 4.66021568, +4.85237205074, 5.0360736751, 5.21234791295, 5.38202942438, +5.54580728515, 5.70425824707, 5.8578708208, 6.00706309899, +6.15219620096, 6.293584586, 6.43150408335, 6.56619822786, +6.69788331895, 6.82675250392, 6.95297910625, 7.07671936269, +7.19811469288, 7.31729359534, 7.43437324249, 7.54946083103, +7.66265473218, 7.77404547667, 7.88371660286, 7.99174539033, +8.09820349737, 8.2031575174, 8.30666946646, 8.408797212, +8.50959485133, 8.60911304685, 8.70739932387, 8.80449833604, +8.90045210255, 8.9953002208, 9.08908005739, 9.1818269202, +9.27357421386, 9.36435358028, 9.45419502625, 9.54312703934, +9.63117669346, 9.71836974515, 9.80473072166, 9.89028300153, +9.97504888851, 10.0590496796, 10.1423057273, 10.2248364976, +10.3066606229, 10.3877959509, 10.4682595908, 10.5480679547, +10.627236797, 10.7057812508, 10.7837158615, 10.8610546184, +10.937810984, 11.0139979213, 11.0896279195, 11.1647130175, +11.239264827, 11.3132945527, 11.3868130127, 11.4598306561, +11.5323575813, 11.6044035514, 11.6759780103, 11.7470900966, +11.8177486575, 11.8879622617, 11.9577392113, 12.0270875534, +12.0960150909, 12.1645293929, 12.2326378041, 12.3003474542, +12.3676652666, 12.4345979668, 12.5011520898, 12.5673339879, +12.6331498379, 12.6986056474, 12.7637072616, 12.8284603691, +12.8928705079, 12.9569430709, 13.0206833111, 13.0840963465, +13.1471871653, 13.2099606302, 13.2724214826, 13.3345743473, +13.3964237358, 13.4579740508, 13.5192295895, 13.5801945472, +13.6408730205, 13.7012690107, 13.7613864269, 13.8212290886, +13.8808007292, 13.9401049979, 13.9991454631, 14.0579256142, +14.1164488646, 14.1747185536, 14.2327379487, 14.2905102481, +14.3480385822, 14.405326016, 14.4623755507, 14.5191901259, +14.5757726211, 14.6321258572, 14.6882525988, 14.7441555552, +14.799837382, 14.8553006829, 14.9105480108, 14.9655818692, +15.0204047137, 15.0750189533, 15.1294269511, 15.1836310264, +15.2376334551, 15.291436471, 15.3450422672, 15.3984529967, +15.4516707737, 15.5046976745, 15.5575357382, 15.6101869683, +15.6626533327, 15.7149367651, 15.7670391659, 15.8189624026, +15.870708311, 15.9222786955, 15.9736753304, 16.0248999601, +16.0759543002, 16.1268400377, 16.1775588321, 16.2281123159, +16.2785020951, 16.3287297498, 16.3787968348, 16.4287048805, +16.4784553929, 16.5280498543, 16.5774897242, 16.6267764392, +16.6759114139, 16.7248960413, 16.7737316932, 16.8224197206, +16.8709614543, 16.9193582051, 16.9676112644, 17.0157219046, +17.0636913794, 17.1115209241, 17.1592117562, 17.2067650757, +17.2541820653, 17.3014638909, 17.3486117017, 17.3956266311, +17.4425097961, 17.4892622985, 17.5358852248, 17.5823796464, +17.6287466201, 17.6749871882, 17.721102379, 17.7670932069, +17.8129606726, 17.8587057635, 17.9043294539, 17.9498327054, +17.9952164666, 18.0404816741, 18.0856292522, 18.130660113, +18.1755751573, 18.2203752741, 18.2650613411, 18.309634225, +18.3540947814, 18.3984438554, 18.4426822814, 18.4868108834, +18.5308304754, 18.5747418613, 18.6185458352, 18.6622431813, +18.7058346748, 18.749321081, 18.7927031564, 18.8359816483, +18.8791572952, 18.9222308268, 18.9652029643, 19.0080744204, +19.0508458994, 19.0935180976, 19.136091703, 19.1785673961, +19.2209458492, 19.2632277271, 19.3054136872, 19.3475043792, +19.3895004457, 19.431402522, 19.4732112364, 19.5149272102, +19.5565510578, 19.598083387, 19.6395247987, 19.6808758874, +19.7221372412, 19.7633094418, 19.8043930645, 19.8453886786, +19.8862968473, 19.9271181279, 19.9678530715, 20.0085022239, +20.0490661246, 20.0895453079, 20.1299403023, 20.1702516311, +20.2104798118, 20.250625357, 20.2906887737, 20.3306705641, +20.3705712249, 20.4103912481, 20.4501311206, 20.4897913245, +20.5293723369, 20.5688746305, 20.608298673, 20.6476449276, +20.6869138529, 20.7261059032, 20.7652215281, 20.804261173, +20.843225279, 20.8821142829, 20.9209286172, 20.9596687105, +20.9983349871, 21.0369278674, 21.0754477678, 21.1138951007, +21.1522702748, 21.1905736947, 21.2288057616, 21.2669668728, +21.3050574217, 21.3430777984, 21.3810283893, 21.4189095773, +21.4567217417, 21.4944652585, 21.5321405003, 21.5697478361, +21.6072876319, 21.6447602504, 21.6821660507, 21.7195053891, +21.7567786187, 21.7939860892, 21.8311281476, 21.8682051376, +21.9052174, 21.9421652726, 21.9790490904, 22.0158691853, +22.0526258864, 22.0893195201, 22.1259504099, 22.1625188766, +22.1990252382, 22.2354698102, 22.271852905, 22.3081748329, +22.3444359014, 22.3806364152, 22.4167766768, 22.4528569861, +22.4888776404, 22.5248389346, 22.5607411612, 22.5965846105, +22.63236957, 22.6680963252, 22.7037651592, 22.7393763528, +22.7749301845, 22.8104269307, 22.8458668654, 22.8812502606, +22.9165773859, 22.9518485091, 22.9870638955, 23.0222238086, +23.0573285098, 23.0923782582, 23.1273733112, 23.1623139241, +23.1972003501, 23.2320328405, 23.2668116448, 23.3015370104, +23.3362091828, 23.3708284059, 23.4053949214, 23.4399089693, +23.4743707879, 23.5087806136, 23.5431386809, 23.5774452228, +23.6117004703, 23.6459046528, 23.6800579981, 23.7141607321, +23.7482130793, 23.7822152621, 23.8161675019, 23.8500700178, +23.8839230279, 23.9177267482, 23.9514813936, 23.985187177, +24.0188443102, 24.0524530031, 24.0860134643, 24.1195259009, +24.1529905184, 24.1864075209, 24.2197771112, 24.2530994905, +24.2863748586, 24.3196034138, 24.3527853533, 24.3859208726, +24.4190101661, 24.4520534266, 24.4850508458, 24.5180026138, +24.5509089198, 24.5837699513, 24.6165858946, 24.6493569349, +24.682083256, 24.7147650405, 24.7474024696, 24.7799957236, +24.8125449812, 24.8450504203, 24.8775122172, 24.9099305473, +24.9423055847, 24.9746375025, 25.0069264724, 25.0391726652, +25.0713762504, 25.1035373964, 25.1356562707, 25.1677330394, +25.1997678678, 25.2317609198, 25.2637123586, 25.295622346, +25.3274910431, 25.3593186096, 25.3911052044, 25.4228509853, +25.4545561092, 25.4862207318, 25.517845008, 25.5494290916, +25.5809731355, 25.6124772915, 25.6439417106, 25.6753665429, +25.7067519372, 25.7380980417, 25.7694050037, 25.8006729694, +25.8319020841, 25.8630924924, 25.8942443377, 25.9253577628, +25.9564329095, 25.9874699187, 26.0184689304, 26.0494300839, +26.0803535176, 26.111239369, 26.1420877748, 26.1728988708, +26.2036727921, 26.234409673, 26.2651096468, 26.2957728462, +26.326399403, 26.3569894483, 26.3875431124, 26.4180605247, +26.448541814, 26.4789871082, 26.5093965346, 26.5397702197, +26.5701082891, 26.6004108679, 26.6306780804, 26.66091005, +26.6911068996, 26.7212687513, 26.7513957266, 26.781487946, +26.8115455297, 26.841568597, 26.8715572664, 26.9015116561, +26.9314318831, 26.9613180642, 26.9911703153, 27.0209887518, +27.0507734882, 27.0805246386, 27.1102423164, 27.1399266342, +27.1695777042, 27.1991956378, 27.2287805458, 27.2583325386, +27.2878517257, 27.3173382162, 27.3467921183, 27.3762135401, +27.4056025886, 27.4349593705, 27.4642839918, 27.4935765581, +27.5228371742, 27.5520659444, 27.5812629725, 27.6104283616, +27.6395622145, 27.6686646331, 27.697735719, 27.7267755733, +27.7557842963, 27.7847619879, 27.8137087475, 27.842624674, +27.8715098657, 27.9003644204, 27.9291884353, 27.9579820072, +27.9867452325, 28.0154782067, 28.0441810252, 28.0728537828, +28.1014965736, 28.1301094915, 28.1586926296, 28.1872460809, +28.2157699376, 28.2442642915, 28.2727292341, 28.3011648561, +28.329571248, 28.3579484998, 28.3862967009, 28.4146159405, +28.442906307, 28.4711678886, 28.499400773, 28.5276050474, +28.5557807987, 28.5839281132, 28.6120470769, 28.6401377751, +28.6682002931, 28.6962347154, 28.7242411262, 28.7522196094, +28.7801702483, 28.808093126, 28.8359883249, 28.8638559271, +28.8916960146, 28.9195086685, 28.9472939699, 28.9750519992, +29.0027828366, 29.0304865619, 29.0581632545, 29.0858129932, +29.1134358568, 29.1410319234, 29.1686012709, 29.1961439767, +29.2236601179, 29.2511497712, 29.2786130131, 29.3060499193, +29.3334605657, 29.3608450274, 29.3882033793, 29.415535696, +29.4428420516, 29.4701225201, 29.4973771748, 29.524606089, +29.5518093353, 29.5789869864, 29.6061391142, 29.6332657906, +29.6603670869, 29.6874430744, 29.7144938237, 29.7415194054, +29.7685198894, 29.7954953457, 29.8224458438, 29.8493714526, +29.8762722411, 29.9031482778, 29.9299996309, 29.9568263682, +29.9836285574, 30.0104062657, 30.03715956, 30.0638885071, +30.0905931733, 30.1172736246, 30.1439299267, 30.1705621452, +30.1971703452, 30.2237545916, 30.2503149489, 30.2768514815, +30.3033642532, 30.3298533279, 30.3563187689, 30.3827606395, +30.4091790024, 30.4355739203, 30.4619454555, 30.48829367, +30.5146186255, 30.5409203836, 30.5671990054, 30.593454552, +30.619687084, 30.6458966617, 30.6720833454, 30.698247195, +30.7243882701, 30.75050663, 30.7766023339, 30.8026754406, +30.8287260087, 30.8547540966, 30.8807597623, 30.9067430638, +30.9327040585, 30.958642804, 30.9845593572, 31.010453775, +31.0363261141, 31.0621764309, 31.0880047814, 31.1138112216, +31.1395958072, 31.1653585935, 31.1910996358, 31.216818989, +31.242516708, 31.268192847, 31.2938474606, 31.3194806026, +31.345092327, 31.3706826872, 31.3962517367, 31.4217995287, +31.4473261161, 31.4728315515, 31.4983158875, 31.5237791764, +31.5492214701, 31.5746428207, 31.6000432796, 31.6254228984, +31.6507817282, 31.67611982, 31.7014372247, 31.7267339927, +31.7520101746, 31.7772658204, 31.8025009802, 31.8277157038, +31.8529100406, 31.87808404, 31.9032377513, 31.9283712234, +31.9534845051, 31.9785776449, 32.0036506913, 32.0287036925, +32.0537366963, 32.0787497507, 32.1037429033, 32.1287162015, +32.1536696925, 32.1786034234, 32.2035174411, 32.2284117922, +32.2532865233, 32.2781416806, 32.3029773103, 32.3277934584, +32.3525901706, 32.3773674925, 32.4021254696, 32.426864147, +32.4515835699, 32.4762837831, 32.5009648313, 32.5256267592, +32.5502696109, 32.5748934309, 32.599498263, 32.6240841512, +32.6486511391, 32.6731992702, 32.6977285879, 32.7222391355, +32.7467309558, 32.7712040919, 32.7956585863, 32.8200944816, +32.8445118202, 32.8689106443, 32.8932909959, 32.9176529169, +32.9419964492, 32.9663216342, 32.9906285133, 33.0149171279, +33.0391875191, 33.0634397278, 33.0876737947, 33.1118897607, +33.1360876661, 33.1602675513, 33.1844294565, 33.2085734218, +33.2326994871, 33.256807692, 33.2808980763, 33.3049706794, +33.3290255406, 33.353062699, 33.3770821938, 33.4010840637, +33.4250683475, 33.4490350839, 33.4729843113, 33.4969160679, +33.520830392, 33.5447273217, 33.5686068947, 33.5924691489, +33.616314122, 33.6401418513, 33.6639523743, 33.6877457282, +33.71152195, 33.7352810768, 33.7590231452, 33.7827481922, +33.8064562541, 33.8301473674, 33.8538215684, 33.8774788933, +33.9011193781, 33.9247430588, 33.948349971, 33.9719401506, +33.9955136329, 34.0190704534, 34.0426106474, 34.0661342501, +34.0896412964, 34.1131318213, 34.1366058595, 34.1600634458, +34.1835046147, 34.2069294006, 34.2303378377, 34.2537299604, +34.2771058026, 34.3004653983, 34.3238087813, 34.3471359853, +34.370447044, 34.3937419907, 34.4170208589, 34.4402836819, +34.4635304926, 34.4867613243, 34.5099762097, 34.5331751816, +34.5563582728, 34.5795255158, 34.6026769431, 34.625812587, +34.6489324797, 34.6720366535, 34.6951251402, 34.7181979718, +34.7412551801, 34.7642967968, 34.7873228536, 34.8103333817, +34.8333284128, 34.8563079779, 34.8792721083, 34.9022208351, +34.9251541892, 34.9480722013, 34.9709749024, 34.993862323, +35.0167344937, 35.0395914449, 35.0624332071, 35.0852598103, +35.1080712847, 35.1308676605, 35.1536489675, 35.1764152355, +35.1991664944, 35.2219027737, 35.2446241031, 35.2673305119, +35.2900220296, 35.3126986853, 35.3353605082, 35.3580075275, +35.3806397721, 35.4032572708, 35.4258600525, 35.4484481458, +35.4710215794, 35.4935803816, 35.5161245811, 35.538654206, +35.5611692847, 35.5836698452, 35.6061559156, 35.6286275239, +35.6510846979, 35.6735274653, 35.695955854, 35.7183698915, +35.7407696053, 35.7631550228, 35.7855261715, 35.8078830784, +35.8302257708, 35.8525542758, 35.8748686204, 35.8971688314, +35.9194549357, 35.94172696, 35.963984931, 35.9862288753, +36.0084588193, 36.0306747894, 36.0528768121, 36.0750649133, +36.0972391195, 36.1193994565, 36.1415459505, 36.1636786272, +36.1857975127, 36.2079026325, 36.2299940124, 36.252071678, +36.2741356547, 36.296185968, 36.3182226433, 36.3402457058, +36.3622551806, 36.3842510931, 36.4062334681, 36.4282023307, +36.4501577057, 36.472099618, 36.4940280922, 36.5159431531, +36.5378448252, 36.559733133, 36.581608101, 36.6034697535, +36.6253181149, 36.6471532092, 36.6689750607, 36.6907836935, +36.7125791315, 36.7343613986, 36.7561305187, 36.7778865156, +36.7996294129, 36.8213592344, 36.8430760035, 36.8647797438, +36.8864704787, 36.9081482315, 36.9298130256, 36.951464884, +36.9731038301, 36.9947298868, 37.0163430771, 37.0379434241, +37.0595309504, 37.0811056791, 37.1026676327, 37.124216834, +37.1457533055, 37.1672770699, 37.1887881495, 37.2102865668, +37.231772344, 37.2532455036, 37.2747060675, 37.2961540581, +37.3175894974, 37.3390124074, 37.36042281, 37.3818207271, +37.4032061805, 37.4245791921, 37.4459397834, 37.4672879761, +37.4886237918, 37.5099472519, 37.531258378, 37.5525571913, +37.5738437133, 37.5951179651, 37.616379968, 37.6376297431, +37.6588673114, 37.680092694, 37.7013059119, 37.7225069858, +37.7436959367, 37.7648727854, 37.7860375524, 37.8071902586, +37.8283309244, 37.8494595705, 37.8705762172, 37.8916808851, +37.9127735945, 37.9338543657, 37.9549232189, 37.9759801743, +37.9970252521, 38.0180584723, 38.039079855, 38.0600894201, +38.0810871876, 38.1020731772, 38.1230474088, 38.1440099022, +38.1649606769, 38.1858997527, 38.2068271491, 38.2277428857, +38.2486469818, 38.269539457, 38.2904203306, 38.3112896218, +38.3321473499, 38.3529935342, 38.3738281938, 38.3946513477, +38.4154630151, 38.4362632149, 38.457051966, 38.4778292874, +38.4985951977, 38.519349716, 38.5400928608, 38.5608246508, +38.5815451047, 38.602254241, 38.6229520783, 38.643638635, +38.6643139297, 38.6849779805, 38.7056308059, 38.7262724241, +38.7469028534, 38.767522112, 38.7881302179, 38.8087271892, +38.829313044, 38.8498878003, 38.870451476, 38.8910040889, +38.9115456569, 38.9320761978, 38.9525957293, 38.9731042692, +38.993601835, 39.0140884444, 39.0345641149, 39.055028864, +39.0754827091, 39.0959256678, 39.1163577573, 39.136778995, +39.1571893981, 39.1775889839, 39.1979777695, 39.2183557721, +39.2387230088, 39.2590794966, 39.2794252526, 39.2997602936, +39.3200846366, 39.3403982984, 39.3607012959, 39.3809936458, +39.401275365, 39.4215464701, 39.4418069776, 39.4620569043, +39.4822962667, 39.5025250812, 39.5227433645, 39.5429511329, +39.5631484027, 39.5833351904, 39.6035115122, 39.6236773845, +39.6438328234, 39.6639778451, 39.6841124657, 39.7042367014, +39.7243505682, 39.744454082, 39.7645472589, 39.7846301149, +39.8047026656, 39.8247649272, 39.8448169153, 39.8648586456, +39.884890134, 39.9049113962, 39.9249224477, 39.9449233041, +39.9649139811, 39.9848944942, 40.0048648588, 40.0248250904, +40.0447752044, 40.0647152162, 40.0846451411, 40.1045649944, +40.1244747914, 40.1443745472, 40.164264277, 40.1841439961, +40.2040137194, 40.2238734621, 40.2437232392, 40.2635630656, +40.2833929563, 40.3032129263 +}; + +static score_t reg2(int N){ + if (N<=1000) + return r2_1000[N]; + else + return exp(0.5*log(N*pi/2) + + sqrt(8/(9*N*pi)) + + (1.0/12 - 4/(9*pi))/N); +} + +score_t reg(int N, int K){ + + if (K == 1) + return 1.0 ; + else if (K == 2) + return reg2(N); + else { + int k = 0; + score_t rk_2 = reg(N,1); + score_t rk_1 = reg(N,2); + score_t rk = 0; + for (k=3; k<=K;++k){ + rk = rk_1 + rk_2 /(k-2)*N; + rk_2 = rk_1; + rk_1 = rk; + } + return rk; + } +} + +#ifdef MAIN +#include +#include + + int main(int argc, char* argv[]){ + int N, K; + + if (argc != 3) { + fprintf(stderr, "Usage: reg N K\n"); + exit(2); + } + + N = atoi(argv[1]); + K = atoi(argv[2]); + + if (K<=0){ + fprintf(stderr, "K has to be a positive integer.\n"); + exit(2); + } + + if (N<0){ + fprintf(stderr, "N has to be a non-negative integer.\n"); + exit(2); + } + + + printf("%f\n", reg(N,K)); + + return 0; + } + +#endif diff --git a/sourcecodes/k-best/src/reg.h b/sourcecodes/k-best/src/reg.h new file mode 100644 index 00000000..2422f2e5 --- /dev/null +++ b/sourcecodes/k-best/src/reg.h @@ -0,0 +1,7 @@ +#ifndef __REG_H__ +#define __REG_H__ +#include "cfg.h" + +extern score_t reg(int N, int K); + +#endif diff --git a/sourcecodes/k-best/src/score_net.c b/sourcecodes/k-best/src/score_net.c new file mode 100644 index 00000000..fccec93f --- /dev/null +++ b/sourcecodes/k-best/src/score_net.c @@ -0,0 +1,39 @@ +#include "cfg.h" + +#include +#include +#include + +#include "files.h" +#include "varpar.h" + +score_t score_net(char* netfile, char* dirprefix){ + score_t score = 0.0; + + int v; + varset_t varset; + FILE* netf = strcmp(netfile,"-") ? fopen(netfile, "rt") : stdin; + for(v=0; 1 == fscanf(netf, "%u", &varset); ++v){ + score_t vscore; + FILE* fin = open_file(dirprefix, v, "", "rb"); + fseek(fin, varset2parset(v, varset)*sizeof(score_t), SEEK_SET); + fread(&vscore, sizeof(score_t), 1, fin); + fclose(fin); + score += vscore; + } + fclose(netf); + return score; +} + +int main(int argc, char* argv[]) +{ + + if (argc!=3) { + fprintf(stderr, "Usage: score_net netfile dirname\n"); + return 1; + } + + printf("%.3f\n", score_net(argv[1], argv[2])); + + return 0; +} diff --git a/sourcecodes/k-best/src/score_nets.c b/sourcecodes/k-best/src/score_nets.c new file mode 100644 index 00000000..44ab3a27 --- /dev/null +++ b/sourcecodes/k-best/src/score_nets.c @@ -0,0 +1,54 @@ +#include "cfg.h" + +#include +#include +#include + +#include "files.h" +#include "varpar.h" + +score_t score_net(char* netfile, char* dirprefix, int nof_vars, FILE** fins){ + score_t score = 0.0; + + int v; + varset_t varset; + FILE* netf = strcmp(netfile,"-") ? fopen(netfile, "r") : stdin; + for(v=0; 1 == fscanf(netf, "%u", &varset); ++v){ + score_t vscore; + fseek(fins[v], varset2parset(v, varset)*sizeof(score_t), SEEK_SET); + fread(&vscore, sizeof(score_t), 1, fins[v]); + rewind(fins[v]); + score += vscore; + } + fclose(netf); + return score; +} + +void score_nets(char* netfiles, char* dirprefix){ + char netfile[1024]; + FILE* netfs = strcmp(netfiles,"-") ? fopen(netfiles, "r") : stdin; + int nof_vars = 0; + FILE** fins = NULL; + while(EOF != fscanf(netfs, "%s", netfile)){ + if (0 == nof_vars){ + nof_vars = nof_lines(netfile); + fins = open_files(nof_vars, dirprefix,"","r"); + } + printf("%.3f\n", score_net(netfile, dirprefix, nof_vars, fins)); + } + free_files(nof_vars,fins); + fclose(netfs); +} + +int main(int argc, char* argv[]) +{ + + if (argc!=3) { + fprintf(stderr, "Usage: score_nets netfiles dirname\n"); + return 1; + } + + score_nets(argv[1], argv[2]); + + return 0; +} diff --git a/sourcecodes/k-best/src/varpar.c b/sourcecodes/k-best/src/varpar.c new file mode 100644 index 00000000..584799c6 --- /dev/null +++ b/sourcecodes/k-best/src/varpar.c @@ -0,0 +1,19 @@ +#include "varpar.h" + +varset_t parset2varset(int v, varset_t set){ + varset_t sinkleton = 1U<>1) | low; +} + diff --git a/sourcecodes/k-best/src/varpar.h b/sourcecodes/k-best/src/varpar.h new file mode 100644 index 00000000..3f28b64b --- /dev/null +++ b/sourcecodes/k-best/src/varpar.h @@ -0,0 +1,9 @@ +#ifndef __VARPAR_H__ +#define __VARPAR_H__ + +#include "cfg.h" + +varset_t parset2varset(int v, varset_t set); +varset_t varset2parset(int v, varset_t set); + +#endif diff --git a/sourcecodes/k-best/src/varpar.o b/sourcecodes/k-best/src/varpar.o new file mode 100644 index 00000000..4fa1b165 Binary files /dev/null and b/sourcecodes/k-best/src/varpar.o differ diff --git a/sourcecodes/k-best/src/xtab.h b/sourcecodes/k-best/src/xtab.h new file mode 100644 index 00000000..51d424c1 --- /dev/null +++ b/sourcecodes/k-best/src/xtab.h @@ -0,0 +1,80 @@ +#ifndef _XTAB_H_ +#define _XTAB_H_ + +#include + +typedef struct xentry xentry; +typedef unsigned char uchar; +typedef unsigned int uint; + + +//Each xentry is one node in the linked list +struct xentry { + uchar* key; //Key value as a string + int* val; //value as int array + xentry* next; //pointer to the next xentry i.e. the next node in this linked list + uint h; /* h denotes hash key i.e to which row in the hash table it belongs to*/ +}; + +typedef struct xtab xtab; + +struct xtab { + int range; /* range of hash keys, for example 1024 */ //i.e. number of rows in hash table + xentry** ix; /* hashkey to xentry pointer mapping */ //pointer pointing to start of each row in the hash table + xentry* free; /* pointer to unused xentries */ + xentry* xentries; //Each entry in hash table +}; + +//Number of entries in the hash table are counted by (free-xentries) +#define xcount(t) ((t)->free - (t)->xentries) + +/*Reset memory to 0's from ix to the end of last xentry*/ + +# define xreset(t) \ +{\ + memset((t)->ix, 0, (t)->range * sizeof(xentry*));\ + memset((t)->xentries, 0, xcount(t) * sizeof(xentry));\ + (t)->free = (t)->xentries;\ +} + +xentry* xadd(xtab* t, uint h, uchar* key, int klen, int* new) { + + /* starting from ix[h], go through x:s and try to find key */ + + //Search in the linked list for that particular hash value and see if such a linked list has the particular key already + xentry* x = t->ix[h]; + xentry** p = t->ix + h; /* prev-pointer to chain the possibly new xentry */ + + //Compares the entry to be added with old entries. If the same key value does not already exist in the linked list, then just go to the next xentry and compare + //memcmp returns 0 if the strings are identical + for(; x && memcmp(x->key, key, klen); x = x->next) + p = &(x->next); + + //If the key value is not already present in the linked list for that hash value + *new = !x; + + if(*new){ + x = (t->free)++; //Free memory starts one block later because new element is inserted + x->h = h; //Set hash value + *p = x; /* link the new x */ + } + + return x; +} + +xtab* xcreate(int range, int nof_xentries){ + xtab* t = calloc(1,sizeof(xtab)); + t->range = range; + t->ix = calloc(range, sizeof(xentry*)); + t->xentries = calloc(nof_xentries, sizeof(xentry)); + t->free = t->xentries; + return t; +} + +void xdestroy(xtab* t) { + free(t->ix); + free(t->xentries); + free(t); +} + +#endif -- cgit 1.4.1