diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine')
23 files changed, 4760 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries new file mode 100644 index 00000000..cc6f3f5b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries @@ -0,0 +1,12 @@ +/cliques_from_engine.m/1.1.1.1/Wed May 29 15:59:56 2002// +/clq_containing_nodes.m/1.1.1.1/Wed May 29 15:59:56 2002// +/collect_evidence.c/1.1.1.1/Wed May 29 15:59:56 2002// +/distribute_evidence.c/1.1.1.1/Wed May 29 15:59:56 2002// +/enter_evidence.m/1.1.1.1/Wed May 29 15:59:56 2002// +/enter_soft_evidence.m/1.1.1.1/Wed May 29 15:59:56 2002// +/init_pot.c/1.1.1.1/Wed May 29 15:59:56 2002// +/jtree_sparse_inf_engine.m/1.1.1.1/Sat Jan 18 22:11:32 2003// +/marginal_family.m/1.1.1.1/Wed May 29 15:59:56 2002// +/marginal_nodes.m/1.1.1.1/Wed May 29 15:59:56 2002// +/set_fields.m/1.1.1.1/Wed May 29 15:59:56 2002// +D diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries.Log b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries.Log new file mode 100644 index 00000000..61d96f3f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Entries.Log @@ -0,0 +1 @@ +A D/old//// diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Repository new file mode 100644 index 00000000..ccd02123 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/inference/static/@jtree_sparse_inf_engine diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Root b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/cliques_from_engine.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/cliques_from_engine.m new file mode 100644 index 00000000..cd9d871d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/cliques_from_engine.m @@ -0,0 +1,5 @@ +function cliques = cliques_from_engine(engine) +% CLIQUES_FROM_ENGINE Return the cliques stored inside the inf. engine (jtree) +% cliques = cliques_from_engine(engine) + +cliques = engine.cliques; diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/clq_containing_nodes.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/clq_containing_nodes.m new file mode 100644 index 00000000..8904fa49 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/clq_containing_nodes.m @@ -0,0 +1,24 @@ +function c = clq_containing_nodes(engine, nodes, fam) +% CLQ_CONTAINING_NODES Find the lightest clique (if any) that contains the set of nodes +% c = clq_containing_nodes(engine, nodes, family) +% +% If the optional 'family' argument is specified, it means nodes = family(nodes(end)). +% (This is useful since clq_ass_to_node is not accessible to outsiders.) +% Returns c=-1 if there is no such clique. + +if nargin < 3, fam = 0; else fam = 1; end + +if length(nodes)==1 + c = engine.clq_ass_to_node(nodes(1)); +%elseif fam +% c = engine.clq_ass_to_node(nodes(end)); +else + B = engine.cliques_bitv; + w = engine.clique_weight; + clqs = find(all(B(:,nodes), 2)); % all selected columns must be 1 + if isempty(clqs) + c = -1; + else + c = clqs(argmin(w(clqs))); + end +end diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/collect_evidence.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/collect_evidence.c new file mode 100644 index 00000000..8480c701 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/collect_evidence.c @@ -0,0 +1,634 @@ +/* C mex for collect_evidence.c in @jtree_sparse_inf_engine directory */ +/* File enter_evidence.m in directory @jtree_sparse_inf_engine call it*/ + +/******************************************/ +/* collect_evidence has 3 input & 2 output*/ +/* engine */ +/* clpot */ +/* seppot */ +/* */ +/* clpot */ +/* seppot */ +/******************************************/ + +#include <math.h> +#include <stdlib.h> +#include "mex.h" + +int compare(const void* src1, const void* src2){ + int i1 = *(int*)src1 ; + int i2 = *(int*)src2 ; + return i1-i2 ; +} + +void ind_subv(int index, const int *cumprod, int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i<n; i++){ + index += subv[i] * cumprod[i]; + } + return index; +} + +void compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){ + int i, j; + int *eff_cumprod, *subv, *diffsize, *diff_cumprod; + + subv = malloc(diffdim * sizeof(int)); + eff_cumprod = malloc(diffdim * sizeof(int)); + diffsize = malloc(diffdim * sizeof(int)); + diff_cumprod = malloc(diffdim * sizeof(int)); + for(i=0; i<diffdim; i++){ + eff_cumprod[i] = bCumprod[dmask[i]]; + diffsize[i] = (int)pbSize[dmask[i]]; + } + diff_cumprod[0] = 1; + for(i=0; i<diffdim-1; i++){ + diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i]; + } + for(i=0; i<ND; i++){ + ind_subv(i, diff_cumprod, diffdim, subv); + weight[i] = 0; + for(j=0; j<diffdim; j++){ + weight[i] += eff_cumprod[j] * subv[j]; + } + } + free(eff_cumprod); + free(subv); + free(diffsize); + free(diff_cumprod); +} + +void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){ + double *ptr; + void *newptr; + int *ir, *jc; + int nbytes; + + if(new_nzmax == old_nzmax) return; + nbytes = new_nzmax * sizeof(*ptr); + ptr = mxGetPr(spArray); + newptr = mxRealloc(ptr, nbytes); + mxSetPr(spArray, newptr); + nbytes = new_nzmax * sizeof(*ir); + ir = mxGetIr(spArray); + newptr = mxRealloc(ir, nbytes); + mxSetIr(spArray, newptr); + jc = mxGetJc(spArray); + jc[0] = 0; + jc[1] = new_nzmax; + mxSetNzmax(spArray, new_nzmax); +} + +mxArray* convert_ill_table_to_sparse(const double *Table, const int *sequence, const int nzCounts, const int N){ + mxArray *spTable; + int i, temp, *irs, *jcs, count=0; + double *sr; + + spTable = mxCreateSparse(N, 1, nzCounts, mxREAL); + sr = mxGetPr(spTable); + irs = mxGetIr(spTable); + jcs = mxGetJc(spTable); + + jcs[0] = 0; + jcs[1] = nzCounts; + + for(i=0; i<nzCounts; i++){ + irs[i] = sequence[count]; + count++; + temp = sequence[count]; + sr[i] = Table[temp]; + count++; + } + return spTable; +} + +void multiply_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex, nzCounts=0; + int *samemask, *diffmask, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *sequence, *weight; + double *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *spr, *bpr; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + if(sdim == 0){ + pTemp = mxCreateSparse(NB, 1, NB, mxREAL); + mxSetField(bigPot, 0, "T", pTemp); + bpr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + sjc[0] = 0; + sjc[1] = NB; + for(i=0; i<NB; i++){ + bpr[i] = *spr; + sir[i] = i; + } + return; + } + + NS = 1; + for(i=0; i<sdim; i++){ + NS *= (int)psSize[i]; + } + ND = NB / NS; + + if(ND == 1){ + pTemp1 = mxGetField(smallPot, 0, "T"); + pTemp = mxDuplicateArray(pTemp1); + mxSetField(bigPot, 0, "T", pTemp); + return; + } + + + NZB = ND * NZS; + + diffdim = bdim - sdim; + sequence = malloc(NZB * 2 * sizeof(int)); + bigTable = malloc(NZB * sizeof(double)); + samemask = malloc(sdim * sizeof(int)); + diffmask = malloc(diffdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + weight = malloc(ND * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + count1 = 0; + for(i=0; i<bdim; i++){ + match = 0; + for(j=0; j<sdim; j++){ + if(pbDomain[i] == psDomain[j]){ + samemask[count] = i; + match = 1; + count++; + break; + } + } + if(match == 0){ + diffmask[count1] = i; + count1++; + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + count = 0; + compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim); + for(i=0; i<NZS; i++){ + sindex = sir[i]; + ind_subv(sindex, sCumprod, sdim, ssubv); + temp = 0; + for(j=0; j<sdim; j++){ + temp += ssubv[j] * bCumprod[samemask[j]]; + } + for(j=0; j<ND; j++){ + bindex = weight[j] + temp; + bigTable[nzCounts] = spr[i]; + sequence[count] = bindex; + count++; + sequence[count] = nzCounts; + nzCounts++; + count++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + qsort(sequence, nzCounts, sizeof(int) * 2, compare); + pTemp = convert_ill_table_to_sparse(bigTable, sequence, nzCounts, NB); + mxSetField(bigPot, 0, "T", pTemp); + + free(sequence); + free(bigTable); + free(samemask); + free(diffmask); + free(bCumprod); + free(sCumprod); + free(weight); + free(ssubv); +} + +void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0; + int *mask, *result, *bir, *sir, *rir, *bjc, *sjc, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + if(sdim == 0){ + for(i=0; i<NZB; i++){ + bpr[i] *= *spr; + } + return; + } + + pTemp1 = mxCreateSparse(NB, 1, NZB, mxREAL); + rpr = mxGetPr(pTemp1); + rir = mxGetIr(pTemp1); + rjc = mxGetJc(pTemp1); + rjc[0] = 0; + rjc[1] = NZB; + + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); + if(result){ + position = result - sir; + rpr[nzCounts] = bpr[i] * spr[position]; + rir[nzCounts] = bindex; + nzCounts++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + reset_nzmax(pTemp1, NZB, nzCounts); + mxSetField(bigPot, 0, "T", pTemp1); + + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + +mxArray* marginal_null_to_spPot(const mxArray *bigPot, const mxArray *sDomain, const int maximize){ + int i, j, count, bdim, sdim, NB, NS, ND; + int *mask, *sir, *sjc; + double *pbDomain, *psDomain, *pbSize, *psSize, *spr; + mxArray *pTemp, *smallPot; + const char *field_names[] = {"domain", "T", "sizes"}; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + psDomain = mxGetPr(sDomain); + sdim = mxGetNumberOfElements(sDomain); + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + + smallPot = mxCreateStructMatrix(1, 1, 3, field_names); + pTemp = mxDuplicateArray(sDomain); + mxSetField(smallPot, 0, "domain", pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + if(sdim == 0){ + pTemp = mxCreateSparse(1, 1, 1, mxREAL); + mxSetField(smallPot, 0, "T", pTemp); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + *spr = 0; + *sir = 0; + sjc[0] = 0; + sjc[1] = 1; + if(maximize) *spr = 1; + else *spr = NB; + + pTemp = mxCreateDoubleMatrix(1, 1, mxREAL); + *mxGetPr(pTemp) = 1; + mxSetField(smallPot, 0, "sizes", pTemp); + return smallPot; + } + + mask = malloc(sdim * sizeof(int)); + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + pTemp = mxCreateDoubleMatrix(1, count, mxREAL); + psSize = mxGetPr(pTemp); + NS = 1; + for(i=0; i<count; i++){ + psSize[i] = pbSize[mask[i]]; + NS *= (int)psSize[i]; + } + mxSetField(smallPot, 0, "sizes", pTemp); + + ND = NB / NS; + + pTemp = mxCreateSparse(NS, 1, NS, mxREAL); + mxSetField(smallPot, 0, "T", pTemp); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + if(maximize){ + for(i=0; i<NS; i++){ + spr[i] = 1; + sir[i] = i; + } + } + else{ + for(i=0; i<NS; i++){ + spr[i] = ND; + sir[i] = i; + } + } + sjc[0] = 0; + sjc[1] = NS; + + free(mask); + return smallPot; +} + +mxArray* marginal_spPot_to_spPot(const mxArray *bigPot, const mxArray *sDomain, const int maximize){ + int i, j, count, bdim, sdim, NB, NS, NZB, position, bindex, sindex, nzCounts=0; + int *mask, *sequence, *result, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *sTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr; + mxArray *pTemp, *smallPot; + const char *field_names[] = {"domain", "T", "sizes"}; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + psDomain = mxGetPr(sDomain); + sdim = mxGetNumberOfElements(sDomain); + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + smallPot = mxCreateStructMatrix(1, 1, 3, field_names); + pTemp = mxDuplicateArray(sDomain); + mxSetField(smallPot, 0, "domain", pTemp); + + if(sdim == 0){ + pTemp = mxCreateSparse(1, 1, 1, mxREAL); + mxSetField(smallPot, 0, "T", pTemp); + spr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + *spr = 0; + *bir = 0; + bjc[0] = 0; + bjc[1] = 1; + if(maximize){ + for(i=0; i<NZB; i++){ + *spr = (*spr < bpr[i])? bpr[i] : *spr; + } + } + else{ + for(i=0; i<NZB; i++){ + *spr += bpr[i]; + } + } + + pTemp = mxCreateDoubleMatrix(1, 1, mxREAL); + *mxGetPr(pTemp) = 1; + mxSetField(smallPot, 0, "sizes", pTemp); + return smallPot; + } + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + mask = malloc(sdim * sizeof(int)); + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + pTemp = mxCreateDoubleMatrix(1, count, mxREAL); + psSize = mxGetPr(pTemp); + NS = 1; + for(i=0; i<count; i++){ + psSize[i] = pbSize[mask[i]]; + NS *= (int)psSize[i]; + } + mxSetField(smallPot, 0, "sizes", pTemp); + + + sTable = malloc(NZB * sizeof(double)); + sequence = malloc(NZB * 2 * sizeof(double)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + for(i=0; i<NZB; i++)sTable[i] = 0; + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + count = 0; + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sequence, nzCounts, sizeof(int)*2, compare); + if(result){ + position = (result - sequence) / 2; + if(maximize) + sTable[position] = (sTable[position] < bpr[i]) ? bpr[i] : sTable[position]; + else sTable[position] += bpr[i]; + } + else { + if(maximize) + sTable[nzCounts] = (sTable[nzCounts] < bpr[i]) ? bpr[i] : sTable[nzCounts]; + else sTable[nzCounts] += bpr[i]; + sequence[count] = sindex; + count++; + sequence[count] = nzCounts; + nzCounts++; + count++; + } + } + + qsort(sequence, nzCounts, sizeof(int) * 2, compare); + pTemp = convert_ill_table_to_sparse(sTable, sequence, nzCounts, NS); + mxSetField(smallPot, 0, "T", pTemp); + + free(sTable); + free(sequence); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); + + return smallPot; +} + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, n, p, np, pn, loop, loops, nCliques, temp, maximize; + int *collect_order; + double *pr, *pr1; + mxArray *pTemp, *pTemp1, *pPostP, *pClpot, *pSeppot, *pSeparator; + + pTemp = mxGetField(prhs[0], 0, "cliques"); + nCliques = mxGetNumberOfElements(pTemp); + loops = nCliques - 1; + pTemp = mxGetField(prhs[0], 0, "maximize"); + maximize = (int)mxGetScalar(pTemp); + pSeparator = mxGetField(prhs[0], 0, "separator"); + + collect_order = malloc(2 * loops * sizeof(int)); + + pTemp = mxGetField(prhs[0], 0, "postorder"); + pr = mxGetPr(pTemp); + pPostP = mxGetField(prhs[0], 0, "postorder_parents"); + for(i=0; i<loops; i++){ + temp = (int)pr[i] - 1; + pTemp = mxGetCell(pPostP, temp); + pr1 = mxGetPr(pTemp); + collect_order[i] = (int)pr1[0] - 1; + collect_order[i+loops] = temp; + } + + plhs[0] = mxDuplicateArray(prhs[1]); + plhs[1] = mxDuplicateArray(prhs[2]); + + for(loop=0; loop<loops; loop++){ + p = collect_order[loop]; + n = collect_order[loop+loops]; + np = p * nCliques + n; + pn = n * nCliques + p; + pClpot = mxGetCell(plhs[0], n); + pTemp1 = mxGetField(pClpot, 0, "T"); + pTemp = mxGetCell(pSeparator, pn); + if(pTemp1) + pSeppot = marginal_spPot_to_spPot(pClpot, pTemp, maximize); + else pSeppot = marginal_null_to_spPot(pClpot, pTemp, maximize); + mxSetCell(plhs[1], pn, pSeppot); + + pClpot = mxGetCell(plhs[0], p); + pTemp1 = mxGetField(pClpot, 0, "T"); + if(pTemp1) + multiply_spPot_by_spPot(pClpot, pSeppot); + else multiply_null_by_spPot(pClpot, pSeppot); + } + free(collect_order); +} + + + + + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/distribute_evidence.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/distribute_evidence.c new file mode 100644 index 00000000..8147c403 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/distribute_evidence.c @@ -0,0 +1,618 @@ +/* C mex for distribute_evidence.c in @jtree_sparse_inf_engine directory*/ +/* File enter_evidence.m in directory @jtree_sparse_inf_engine call it */ + +/*********************************************/ +/* distribute_evidence has 3 input & 2 output*/ +/* engine */ +/* clpot */ +/* seppot */ +/* */ +/* clpot */ +/* seppot */ +/*********************************************/ + +#include <math.h> +#include <stdlib.h> +#include "mex.h" + +int compare(const void* src1, const void* src2){ + int i1 = *(int*)src1 ; + int i2 = *(int*)src2 ; + return i1-i2 ; +} + +void ind_subv(int index, const int *cumprod, int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i<n; i++){ + index += subv[i] * cumprod[i]; + } + return index; +} + +void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){ + double *ptr; + void *newptr; + int *ir, *jc; + int nbytes; + + if(new_nzmax == old_nzmax) return; + nbytes = new_nzmax * sizeof(*ptr); + ptr = mxGetPr(spArray); + newptr = mxRealloc(ptr, nbytes); + mxSetPr(spArray, newptr); + nbytes = new_nzmax * sizeof(*ir); + ir = mxGetIr(spArray); + newptr = mxRealloc(ir, nbytes); + mxSetIr(spArray, newptr); + jc = mxGetJc(spArray); + jc[0] = 0; + jc[1] = new_nzmax; + mxSetNzmax(spArray, new_nzmax); +} + +void compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){ + int i, j; + int *eff_cumprod, *subv, *diffsize, *diff_cumprod; + + subv = malloc(diffdim * sizeof(int)); + eff_cumprod = malloc(diffdim * sizeof(int)); + diffsize = malloc(diffdim * sizeof(int)); + diff_cumprod = malloc(diffdim * sizeof(int)); + for(i=0; i<diffdim; i++){ + eff_cumprod[i] = bCumprod[dmask[i]]; + diffsize[i] = (int)pbSize[dmask[i]]; + } + diff_cumprod[0] = 1; + for(i=0; i<diffdim-1; i++){ + diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i]; + } + for(i=0; i<ND; i++){ + ind_subv(i, diff_cumprod, diffdim, subv); + weight[i] = 0; + for(j=0; j<diffdim; j++){ + weight[i] += eff_cumprod[j] * subv[j]; + } + } + free(eff_cumprod); + free(subv); + free(diffsize); + free(diff_cumprod); +} + +mxArray* convert_ill_table_to_sparse(const double *Table, const int *sequence, const int nzCounts, const int N){ + mxArray *spTable; + int i, temp, *irs, *jcs, count=0; + double *sr; + + spTable = mxCreateSparse(N, 1, nzCounts, mxREAL); + sr = mxGetPr(spTable); + irs = mxGetIr(spTable); + jcs = mxGetJc(spTable); + + jcs[0] = 0; + jcs[1] = nzCounts; + + for(i=0; i<nzCounts; i++){ + irs[i] = sequence[count]; + count++; + temp = sequence[count]; + sr[i] = Table[temp]; + count++; + } + return spTable; +} + +void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0; + int *mask, *result, *bir, *sir, *rir, *bjc, *sjc, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + if(sdim == 0){ + for(i=0; i<NZB; i++){ + bpr[i] *= *spr; + } + return; + } + + pTemp1 = mxCreateSparse(NB, 1, NZB, mxREAL); + rpr = mxGetPr(pTemp1); + rir = mxGetIr(pTemp1); + rjc = mxGetJc(pTemp1); + rjc[0] = 0; + rjc[1] = NZB; + + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); + if(result){ + position = result - sir; + rpr[nzCounts] = bpr[i] * spr[position]; + rir[nzCounts] = bindex; + nzCounts++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + reset_nzmax(pTemp1, NZB, nzCounts); + mxSetField(bigPot, 0, "T", pTemp1); + + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + +void marginal_spPot_to_spPot(const mxArray *bigPot, mxArray *smallPot, const int maximize){ + int i, j, count, bdim, sdim, NB, NS, NZB, position, bindex, sindex, nzCounts=0; + int *mask, *sequence, *result, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *sTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr; + mxArray *pTemp; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + NS = 1; + for(i=0; i<sdim; i++){ + NS *= (int)psSize[i]; + } + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + if(sdim == 0){ + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + *spr = 0; + if(maximize){ + for(i=0; i<NZB; i++){ + *spr = (*spr < bpr[i])? bpr[i] : *spr; + } + } + else{ + for(i=0; i<NZB; i++){ + *spr += bpr[i]; + } + } + return; + } + + mask = malloc(sdim * sizeof(int)); + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + + sTable = malloc(NZB * sizeof(double)); + sequence = malloc(NZB * 2 * sizeof(double)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + for(i=0; i<NZB; i++){ + sTable[i] = 0; + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + count = 0; + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sequence, nzCounts, sizeof(int)*2, compare); + if(result){ + position = (result - sequence) / 2; + if(maximize) + sTable[position] = (sTable[position] < bpr[i]) ? bpr[i] : sTable[position]; + else sTable[position] += bpr[i]; + } + else { + if(maximize) + sTable[nzCounts] = (sTable[nzCounts] < bpr[i]) ? bpr[i] : sTable[nzCounts]; + else sTable[nzCounts] += bpr[i]; + sequence[count] = sindex; + count++; + sequence[count] = nzCounts; + nzCounts++; + count++; + } + } + + pTemp = mxGetField(smallPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + qsort(sequence, nzCounts, sizeof(int) * 2, compare); + pTemp = convert_ill_table_to_sparse(sTable, sequence, nzCounts, NS); + mxSetField(smallPot, 0, "T", pTemp); + + free(sTable); + free(sequence); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + +void divide_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex; + int *samemask, *diffmask, *rir, *rjc, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *weight; + double *pbDomain, *psDomain, *pbSize, *psSize, *rpr, *spr, value; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + if(sdim == 0){ + pTemp1 = mxGetField(bigPot, 0, "T"); + if(pTemp1)mxDestroyArray(pTemp1); + pTemp = mxCreateSparse(NB, 1, NB, mxREAL); + mxSetField(bigPot, 0, "T", pTemp); + rpr = mxGetPr(pTemp); + rir = mxGetIr(pTemp); + rjc = mxGetJc(pTemp); + rjc[0] = 0; + rjc[1] = NB; + value = *spr; + if(value == 0) value = 1; + for(i=0; i<NB; i++){ + rpr[i] = 1 / value; + rir[i] = i; + } + return; + } + + NS = 1; + for(i=0; i<sdim; i++){ + NS *= (int)psSize[i]; + } + ND = NB / NS; + + pTemp = mxCreateSparse(NB, 1, NB, mxREAL); + rpr = mxGetPr(pTemp); + rir = mxGetIr(pTemp); + rjc = mxGetJc(pTemp); + rjc[0] = 0; + rjc[1] = NB; + for(i=0; i<NB; i++){ + rpr[i] = 1; + rir[i] = i; + } + + NZB = ND * NZS; + + diffdim = bdim - sdim; + samemask = malloc(sdim * sizeof(int)); + diffmask = malloc(diffdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + weight = malloc(ND * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + count1 = 0; + for(i=0; i<bdim; i++){ + match = 0; + for(j=0; j<sdim; j++){ + if(pbDomain[i] == psDomain[j]){ + samemask[count] = i; + match = 1; + count++; + break; + } + } + if(match == 0){ + diffmask[count1] = i; + count1++; + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + count = 0; + compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim); + for(i=0; i<NZS; i++){ + sindex = sir[i]; + ind_subv(sindex, sCumprod, sdim, ssubv); + temp = 0; + for(j=0; j<sdim; j++){ + temp += ssubv[j] * bCumprod[samemask[j]]; + } + for(j=0; j<ND; j++){ + bindex = weight[j] + temp; + rpr[bindex] = 1 / (spr[i]); + } + } + + pTemp1 = mxGetField(bigPot, 0, "T"); + if(pTemp1)mxDestroyArray(pTemp1); + mxSetField(bigPot, 0, "T", pTemp); + + free(samemask); + free(diffmask); + free(bCumprod); + free(sCumprod); + free(weight); + free(ssubv); +} + +void divide_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex; + int *mask, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp1 = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp1); + bir = mxGetIr(pTemp1); + bjc = mxGetJc(pTemp1); + NZB = bjc[1]; + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + if(sdim == 0){ + value = *spr; + if(value == 0)value = 1; + for(i=0; i<NZB; i++){ + bpr[i] /= value; + } + return; + } + + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); + if(result){ + position = result - sir; + bpr[i] /= spr[position]; + } + } + + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, loop, loops, nCliques, temp, count, parent, child, maximize, *distribute_order; + double *pr, *pr1; + mxArray *pTemp, *pPreCh, *pClpot, *pSeppot; + + pTemp = mxGetField(prhs[0], 0, "cliques"); + nCliques = mxGetNumberOfElements(pTemp); + loops = nCliques - 1; + pTemp = mxGetField(prhs[0], 0, "maximize"); + maximize = (int)mxGetScalar(pTemp); + + distribute_order = malloc(2 * loops * sizeof(int)); + pTemp = mxGetField(prhs[0], 0, "preorder"); + pr = mxGetPr(pTemp); + pPreCh = mxGetField(prhs[0], 0, "preorder_children"); + count = 0; + for(i=0; i<nCliques; i++){ + temp = (int)pr[i] - 1; + pTemp = mxGetCell(pPreCh, temp); + pr1 = mxGetPr(pTemp); + loop = mxGetNumberOfElements(pTemp); + for(j=0; j<loop; j++){ + distribute_order[count] = temp; + distribute_order[count + loops] = (int)pr1[j] - 1; + count++; + } + } + + plhs[0] = mxDuplicateArray(prhs[1]); + plhs[1] = mxDuplicateArray(prhs[2]); + + for(loop=0; loop<loops; loop++){ + parent = distribute_order[loop]; + child = distribute_order[loop+loops]; + i = nCliques * child + parent; + pClpot = mxGetCell(plhs[0], child); + pTemp = mxGetField(pClpot, 0, "T"); + pSeppot = mxGetCell(plhs[1], i); + if(pTemp){ + if(mxIsEmpty(pTemp)) + divide_null_by_spPot(pClpot, pSeppot); + else + divide_spPot_by_spPot(pClpot, pSeppot); + } + else divide_null_by_spPot(pClpot, pSeppot); + + pClpot = mxGetCell(plhs[0], parent); + marginal_spPot_to_spPot(pClpot, pSeppot, maximize); + mxSetCell(plhs[1], i, pSeppot); + + pClpot = mxGetCell(plhs[0], child); + multiply_spPot_by_spPot(pClpot, pSeppot); + } + free(distribute_order); +} diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_evidence.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_evidence.m new file mode 100644 index 00000000..86041be2 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_evidence.m @@ -0,0 +1,100 @@ +function [engine, loglik] = enter_evidence(engine, evidence, varargin) +% ENTER_EVIDENCE Add the specified evidence to the network (jtree) +% [engine, loglik] = enter_evidence(engine, evidence, ...) +% +% evidence{i} = [] if X(i) is hidden, and otherwise contains its observed value (scalar or column vector). +% +% The following optional arguments can be specified in the form of name/value pairs: +% [default value in brackets] +% +% maximize - if 1, does max-product instead of sum-product [0] +% soft - a cell array of soft/virtual evidence; +% soft{i} is a prob. distrib. over i's values, or [] [ cell(1,N) ] +% +% e.g., engine = enter_evidence(engine, ev, 'soft', soft_ev) +% +% For backwards compatibility with BNT2, you can also specify the parameters in the following order +% engine = enter_evidence(engine, ev, soft_ev) + +bnet = bnet_from_engine(engine); +ns = bnet.node_sizes(:); +N = length(bnet.dag); + +engine.evidence = evidence; % store this for marginal_nodes with add_ev option + +% set default params +exclude = []; +soft_evidence = cell(1,N); +maximize = 0; + +% parse optional params +args = varargin; +nargs = length(args); +if nargs > 0 + if iscell(args{1}) + soft_evidence = args{1}; + else + for i=1:2:nargs + switch args{i}, + case 'soft', soft_evidence = args{i+1}; + case 'maximize', maximize = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end + end + end +end + +engine.maximize = maximize; + +onodes = find(~isemptycell(evidence)); +hnodes = find(isemptycell(evidence)); +pot_type = determine_pot_type(bnet, onodes); + if strcmp(pot_type, 'cg') + check_for_cd_arcs(onodes, bnet.cnodes, bnet.dag); +end + +hard_nodes = 1:N; +soft_nodes = find(~isemptycell(soft_evidence)); +S = length(soft_nodes); +if S > 0 + assert(pot_type == 'd'); + assert(mysubset(soft_nodes, bnet.dnodes)); +end + +% Evaluate CPDs with evidence, and convert to potentials +pot = cell(1, N+S); +for n=1:N + fam = family(bnet.dag, n); + e = bnet.equiv_class(n); + if isempty(bnet.CPD{e}) + error(['must define CPD ' num2str(e)]) + else + pot{n} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence); + end +end + +for i=1:S + n = soft_nodes(i); + pot{N+i} = dpot(n, ns(n), soft_evidence{n}); +end +clqs = engine.clq_ass_to_node([hard_nodes soft_nodes]); + +[clpot, seppot] = init_pot(engine, clqs, pot, pot_type, onodes); +[clpot, seppot] = collect_evidence(engine, clpot, seppot); +[clpot, seppot] = distribute_evidence(engine, clpot, seppot); +C = length(clpot); +ll = zeros(1, C); +for i=1:C + domain = clpot{i}.domain; + sizes = clpot{i}.sizes; + T = clpot{i}.T; + clpot{i} = dpot(domain, sizes, T); +end + +for i=1:C + [clpot{i}, ll(i)] = normalize_pot(clpot{i}); +end +loglik = ll(1); % we can extract the likelihood from any clique + +engine.clpot = clpot; diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_soft_evidence.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_soft_evidence.m new file mode 100644 index 00000000..59671415 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/enter_soft_evidence.m @@ -0,0 +1,19 @@ +function [clpot, loglik] = enter_soft_evidence(engine, clique, potential, onodes, pot_type) +% ENTER_SOFT_EVIDENCE Add the specified potentials to the network (jtree) +% [clpot, loglik] = enter_soft_evidence(engine, clique, potential, onodes, pot_type, maximize) +% +% We multiply potential{i} onto clique(i) before propagating. +% We return all the modified clique potentials. + +[clpot, seppot] = init_pot(engine, clique, potential, pot_type, onodes); +[clpot, seppot] = collect_evidence(engine, clpot, seppot); +[clpot, seppot] = distribute_evidence(engine, clpot, seppot); + +C = length(clpot); +ll = zeros(1, C); +for i=1:C + [clpot{i}, ll(i)] = normalize_pot(clpot{i}); +end +loglik = ll(1); % we can extract the likelihood from any clique + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/init_pot.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/init_pot.c new file mode 100644 index 00000000..86e09eae --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/init_pot.c @@ -0,0 +1,624 @@ +/* C mex init_pot for in @jtree_sparse_inf_engine directory */ +/* The file enter_evidence.m in directory @jtree_sparse_inf_engine call it*/ + +/**************************************/ +/* init_pot.c has 5 input & 2 output */ +/* engine */ +/* clqs */ +/* pots */ +/* pot_type */ +/* onodes */ +/* */ +/* clpot */ +/* seppot */ +/**************************************/ +#include <math.h> +#include <stdlib.h> +#include "mex.h" + +int compare(const void* src1, const void* src2){ + int i1 = *(int*)src1 ; + int i2 = *(int*)src2 ; + return i1-i2 ; +} + +void ind_subv(int index, const int *cumprod, int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i<n; i++){ + index += subv[i] * cumprod[i]; + } + return index; +} + +void compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){ + int i, j; + int *eff_cumprod, *subv, *diffsize, *diff_cumprod; + + subv = malloc(diffdim * sizeof(int)); + eff_cumprod = malloc(diffdim * sizeof(int)); + diffsize = malloc(diffdim * sizeof(int)); + diff_cumprod = malloc(diffdim * sizeof(int)); + for(i=0; i<diffdim; i++){ + eff_cumprod[i] = bCumprod[dmask[i]]; + diffsize[i] = (int)pbSize[dmask[i]]; + } + diff_cumprod[0] = 1; + for(i=0; i<diffdim-1; i++){ + diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i]; + } + for(i=0; i<ND; i++){ + ind_subv(i, diff_cumprod, diffdim, subv); + weight[i] = 0; + for(j=0; j<diffdim; j++){ + weight[i] += eff_cumprod[j] * subv[j]; + } + } + free(eff_cumprod); + free(subv); + free(diffsize); + free(diff_cumprod); +} + +void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){ + double *ptr; + void *newptr; + int *ir, *jc; + int nbytes; + + if(new_nzmax == old_nzmax) return; + nbytes = new_nzmax * sizeof(*ptr); + ptr = mxGetPr(spArray); + newptr = mxRealloc(ptr, nbytes); + mxSetPr(spArray, newptr); + nbytes = new_nzmax * sizeof(*ir); + ir = mxGetIr(spArray); + newptr = mxRealloc(ir, nbytes); + mxSetIr(spArray, newptr); + jc = mxGetJc(spArray); + jc[0] = 0; + jc[1] = new_nzmax; + mxSetNzmax(spArray, new_nzmax); +} + +mxArray* convert_ill_table_to_sparse(const double *bigTable, const int *sequence, const int nzCounts, const int NB){ + mxArray *spTable; + int i, temp, *irs, *jcs, count=0; + double *sr; + + spTable = mxCreateSparse(NB, 1, nzCounts, mxREAL); + sr = mxGetPr(spTable); + irs = mxGetIr(spTable); + jcs = mxGetJc(spTable); + + jcs[0] = 0; + jcs[1] = nzCounts; + + for(i=0; i<nzCounts; i++){ + irs[i] = sequence[count]; + count++; + temp = sequence[count]; + sr[i] = bigTable[temp]; + count++; + } + return spTable; +} + +void multiply_null_by_fuPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, NB, NS, siz_b, siz_s, ndim, nzCounts=0; + int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *bir, *bjc; + double *pbDomain, *psDomain, *pbSize, *psSize, *spr, *bpr, value; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + siz_b = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + siz_s = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<siz_b; i++){ + NB *= (int)pbSize[i]; + } + NS = 1; + for(i=0; i<siz_s; i++){ + NS *= (int)psSize[i]; + } + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + + pTemp1 = mxCreateSparse(NB, 1, NB, mxREAL); + bpr = mxGetPr(pTemp1); + bir = mxGetIr(pTemp1); + bjc = mxGetJc(pTemp1); + bjc[0] = 0; + bjc[1] = NB; + + if(NS == 1){ + value = *spr; + for(i=0; i<NB; i++){ + bpr[i] = value; + bir[i] = i; + } + nzCounts = NB; + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + reset_nzmax(pTemp1, NB, nzCounts); + mxSetField(bigPot, 0, "T", pTemp1); + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++){ + if(spr[i] != 0){ + bpr[nzCounts] = spr[i]; + bir[nzCounts] = i; + nzCounts++; + } + } + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + reset_nzmax(pTemp1, NB, nzCounts); + mxSetField(bigPot, 0, "T", pTemp1); + return; + } + + mask = malloc(siz_s * sizeof(int)); + count = 0; + for(i=0; i<siz_s; i++){ + for(j=0; j<siz_b; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + ndim = siz_b; + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + sx[i] = (int)pbSize[i]; + sy[i] = 1; + } + for(i=0; i<count; i++){ + sy[mask[i]] = sx[mask[i]]; + } + + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i = 0; i < ndim; i++){ + subs[i] = 0; + s[i] = sx[i] - 1; + } + + for(i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + for(j=0; j<NB; j++){ + if(*spr != 0){ + bpr[nzCounts] = *spr; + bir[nzCounts] = j; + nzCounts++; + } + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + spr -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + spr += cpsy[i]; + break; + } + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + reset_nzmax(pTemp1, NB, nzCounts); + mxSetField(bigPot, 0, "T", pTemp1); + + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); + free(mask); +} + +void multiply_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex, nzCounts=0; + int *samemask, *diffmask, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *sequence, *weight; + double *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *spr; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + NS = 1; + for(i=0; i<sdim; i++){ + NS *= (int)psSize[i]; + } + ND = NB / NS; + + if(ND == 1){ + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + pTemp1 = mxGetField(smallPot, 0, "T"); + pTemp = mxDuplicateArray(pTemp1); + mxSetField(bigPot, 0, "T", pTemp); + return; + } + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + NZB = ND * NZS; + + diffdim = bdim - sdim; + sequence = malloc(NZB * 2 * sizeof(int)); + bigTable = malloc(NZB * sizeof(double)); + samemask = malloc(sdim * sizeof(int)); + diffmask = malloc(diffdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + weight = malloc(ND * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + count1 = 0; + for(i=0; i<bdim; i++){ + match = 0; + for(j=0; j<sdim; j++){ + if(pbDomain[i] == psDomain[j]){ + samemask[count] = i; + match = 1; + count++; + break; + } + } + if(match == 0){ + diffmask[count1] = i; + count1++; + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + count = 0; + compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim); + for(i=0; i<NZS; i++){ + sindex = sir[i]; + ind_subv(sindex, sCumprod, sdim, ssubv); + temp = 0; + for(j=0; j<sdim; j++){ + temp += ssubv[j] * bCumprod[samemask[j]]; + } + for(j=0; j<ND; j++){ + bindex = weight[j] + temp; + bigTable[nzCounts] = spr[i]; + sequence[count] = bindex; + count++; + sequence[count] = nzCounts; + nzCounts++; + count++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + qsort(sequence, nzCounts, sizeof(int) * 2, compare); + pTemp = convert_ill_table_to_sparse(bigTable, sequence, nzCounts, NB); + mxSetField(bigPot, 0, "T", pTemp); + + free(sequence); + free(bigTable); + free(samemask); + free(diffmask); + free(bCumprod); + free(sCumprod); + free(weight); + free(ssubv); +} + +void multiply_spPot_by_fuPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, bdim, sdim, NB, NZB, bindex, sindex, nzCounts=0; + int *mask, *bir, *bjc, *rir, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr, value; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + + pTemp1 = mxCreateSparse(NB, 1, NZB, mxREAL); + rpr = mxGetPr(pTemp1); + rir = mxGetIr(pTemp1); + rjc = mxGetJc(pTemp1); + rjc[0] = 0; + rjc[1] = NZB; + + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + value = spr[sindex]; + if(value != 0){ + rpr[nzCounts] = bpr[i] * value; + rir[nzCounts] = bindex; + nzCounts++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + reset_nzmax(pTemp1, NZB, nzCounts); + mxSetField(bigPot, 0, "T", pTemp1); + + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + +void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0; + int *mask, *result, *bir, *sir, *rir, *bjc, *sjc, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + pTemp1 = mxCreateSparse(NB, 1, NZB, mxREAL); + rpr = mxGetPr(pTemp1); + rir = mxGetIr(pTemp1); + rjc = mxGetJc(pTemp1); + rjc[0] = 0; + rjc[1] = NZB; + + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); + if(result){ + position = result - sir; + rpr[nzCounts] = bpr[i] * spr[position]; + rir[nzCounts] = bindex; + nzCounts++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + reset_nzmax(pTemp1, NZB, nzCounts); + mxSetField(bigPot, 0, "T", pTemp1); + + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, c, loop, nNodes, nCliques, ndomain, ns_num, nOnodes, dims[2]; + double *pClqs, *pr, *pt, *pSize, *eff_ns; + mxArray *pTemp, *pTemp1, *pStruct, *pCliques, *pBigpot, *pSmallpot; + const char *field_names[] = {"domain", "T", "sizes"}; + + nNodes = mxGetNumberOfElements(prhs[1]); + pCliques = mxGetField(prhs[0], 0, "cliques"); + nCliques = mxGetNumberOfElements(pCliques); + pTemp = mxGetField(prhs[0], 0, "actual_node_sizes"); + ns_num = mxGetNumberOfElements(pTemp); + pSize = mxGetPr(pTemp); + + eff_ns = (double *)malloc(ns_num * sizeof(double)); + for(i=0; i<ns_num; i++) eff_ns[i] = pSize[i]; + nOnodes = mxGetNumberOfElements(prhs[4]); + pr = mxGetPr(prhs[4]); + for(i=0; i<nOnodes; i++) eff_ns[(int)pr[i] - 1] = 1; + + plhs[0] = mxCreateCellArray(1, &nCliques); + for(i=0; i<nCliques; i++){ + pStruct = mxCreateStructMatrix(1, 1, 3, field_names); + mxSetCell(plhs[0], i, pStruct); + pTemp = mxGetCell(pCliques, i); + ndomain = mxGetNumberOfElements(pTemp); + pt = mxGetPr(pTemp); + pTemp1 = mxDuplicateArray(pTemp); + mxSetField(pStruct, 0, "domain", pTemp1); + + pTemp = mxCreateDoubleMatrix(1, ndomain, mxREAL); + mxSetField(pStruct, 0, "sizes", pTemp); + pr = mxGetPr(pTemp); + for(j=0; j<ndomain; j++){ + pr[j] = eff_ns[(int)pt[j]-1]; + } + } + + pClqs = mxGetPr(prhs[1]); + for(loop=0; loop<nNodes; loop++){ + c = (int)pClqs[loop] - 1; + pSmallpot = mxGetCell(prhs[2], loop); + pTemp = mxGetField(pSmallpot, 0, "T"); + pBigpot = mxGetCell(plhs[0], c); + pTemp1 = mxGetField(pBigpot, 0, "T"); + if(pTemp1){ + if(mxIsSparse(pTemp)) + multiply_spPot_by_spPot(pBigpot, pSmallpot); + else multiply_spPot_by_fuPot(pBigpot, pSmallpot); + } + else{ + if(mxIsSparse(pTemp)) + multiply_null_by_spPot(pBigpot, pSmallpot); + else multiply_null_by_fuPot(pBigpot, pSmallpot); + } + } + + free(eff_ns); + dims[0] = nCliques; + dims[1] = nCliques; + plhs[1] = mxCreateCellArray(2, dims); +} + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/jtree_sparse_inf_engine.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/jtree_sparse_inf_engine.m new file mode 100644 index 00000000..49dcd69f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/jtree_sparse_inf_engine.m @@ -0,0 +1,126 @@ +function engine = jtree_sparse_inf_engine(bnet, varargin) +% JTREE_SPARSE_INF_ENGINE Junction tree inference engine when CPTs and Potentials are sparse +% engine = jtree_sparse_inf_engine(bnet, ...) +% It differs from jtree_inf_engine with all CPTs and potentials are 1D sparse arrays. +% +% The following optional arguments can be specified in the form of name/value pairs: +% [default value in brackets] +% +% clusters - a cell array of sets of nodes we want to ensure are in the same clique (in addition to families) [ {} ] +% root - the root of the junction tree will be a clique that contains this set of nodes [N] +% stages - stages{t} is a set of nodes we want to eliminate before stages{t+1}, ... [ {1:N} ] +% +% e.g., engine = jtree_inf_engine(bnet, 'maximize', 1); +% +% For more details on the junction tree algorithm, see +% - "Probabilistic networks and expert systems", Cowell, Dawid, Lauritzen and Spiegelhalter, Springer, 1999 +% - "Inference in Belief Networks: A procedural guide", C. Huang and A. Darwiche, +% Intl. J. Approximate Reasoning, 15(3):225-263, 1996. + + +% set default params +N = length(bnet.dag); +clusters = {}; +root = N; +stages = { 1:N }; + +if nargin >= 2 + args = varargin; + nargs = length(args); + if ~isstr(args{1}) + error('the interface to jtree has changed; now, onodes is not allowed and all optional params must be passed by name') + end + for i=1:2:nargs + switch args{i}, + case 'clusters', clusters = args{i+1}; + case 'root', root = args{i+1}; + case 'stages', stages = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end + end +end + +engine = init_fields; +engine = class(engine, 'jtree_sparse_inf_engine', inf_engine(bnet)); + +onodes = bnet.observed; +%[engine.jtree, dummy, engine.cliques, B, w] = dag_to_jtree(bnet, onodes, stages, clusters); + +porder = determine_elim_constraints(bnet, onodes); +strong = ~isempty(porder); +ns = bnet.node_sizes(:); +ns(onodes) = 1; % observed nodes have only 1 possible value +[engine.jtree, root2, engine.cliques, B, w] = ... + graph_to_jtree(moralize(bnet.dag), ns, porder, stages, clusters); + +engine.cliques_bitv = B; +engine.clique_weight = w; +C = length(engine.cliques); +engine.clpot = cell(1,C); + +% Compute the separators between connected cliques. +[is,js] = find(engine.jtree > 0); +engine.separator = cell(C,C); +for k=1:length(is) + i = is(k); j = js(k); + engine.separator{i,j} = find(B(i,:) & B(j,:)); % intersect(cliques{i}, cliques{j}); +end + +% A node can be a member of many cliques, but is assigned to exactly one, to avoid +% double-counting its CPD. We assign node i to clique c if c is the "lightest" clique that +% contains i's family, so it can accomodate its CPD. + +engine.clq_ass_to_node = zeros(1, N); +for i=1:N + %c = clq_containing_nodes(engine, family(bnet.dag, i)); + clqs_containing_family = find(all(B(:,family(bnet.dag, i)), 2)); % all selected columns must be 1 + c = clqs_containing_family(argmin(w(clqs_containing_family))); + engine.clq_ass_to_node(i) = c; +end + +% Make the jtree rooted, so there is a fixed message passing order. +engine.root_clq = clq_containing_nodes(engine, root); +if engine.root_clq <= 0 + error(['no clique contains ' num2str(root)]); +end + +[engine.jtree, engine.preorder, engine.postorder] = mk_rooted_tree(engine.jtree, engine.root_clq); + +% collect +engine.postorder_parents = cell(1,length(engine.postorder)); +for n=engine.postorder(:)' + engine.postorder_parents{n} = parents(engine.jtree, n); +end +% distribute +engine.preorder_children = cell(1,length(engine.preorder)); +for n=engine.preorder(:)' + engine.preorder_children{n} = children(engine.jtree, n); +end + +ns = bnet.node_sizes; +engine.actual_node_sizes = ns; + + +%%%%%%%% + +function engine = init_fields() + +engine.jtree = []; +engine.cliques = []; +engine.separator = []; +engine.cliques_bitv = []; +engine.clique_weight = []; +engine.clpot = []; +engine.clq_ass_to_node = []; +engine.root_clq = []; +engine.preorder = []; +engine.postorder = []; +engine.preorder_children = []; +engine.postorder_parents = []; +engine.maximize = []; +engine.evidence = []; +engine.actual_node_sizes = []; + + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_family.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_family.m new file mode 100644 index 00000000..eff60ca2 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_family.m @@ -0,0 +1,11 @@ +function marginal = marginal_family(engine, i, add_ev) +% MARGINAL_FAMILY Compute the marginal on the specified family (jtree) +% marginal = marginal_family(engine, i) + +if nargin < 3, add_ev = 0; end +assert(~add_ev); + +bnet = bnet_from_engine(engine); +fam = family(bnet.dag, i); +c = engine.clq_ass_to_node(i); +marginal = pot_to_marginal(marginalize_pot(engine.clpot{c}, fam)); diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_nodes.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_nodes.m new file mode 100644 index 00000000..6413172c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/marginal_nodes.m @@ -0,0 +1,22 @@ +function marginal = marginal_nodes(engine, query, add_ev) +% MARGINAL_NODES Compute the marginal on the specified query nodes (jtree) +% marginal = marginal_nodes(engine, query, add_ev) +% +% 'query' must be a subset of some clique; an error will be raised if not. +% add_ev is an optional argument; if 1, we will "inflate" the marginal of observed nodes +% to their original size, adding 0s to the positions which contradict the evidence + +if nargin < 3, add_ev = 0; end + +c = clq_containing_nodes(engine, query); +if c == -1 + error(['no clique contains ' num2str(query)]); +end +marginal = pot_to_marginal(marginalize_pot(engine.clpot{c}, query, engine.maximize)); + +if add_ev + bnet = bnet_from_engine(engine); + %marginal = add_ev_to_dmarginal(marginal, engine.evidence, bnet.node_sizes); + marginal = add_evidence_to_gmarginal(marginal, engine.evidence, bnet.node_sizes, bnet.cnodes); +end + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Entries new file mode 100644 index 00000000..f74fd729 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Entries @@ -0,0 +1,6 @@ +/collect_evidence.c/1.1.1.1/Wed May 29 15:59:56 2002// +/distribute_evidence.c/1.1.1.1/Wed May 29 15:59:56 2002// +/init_pot.c/1.1.1.1/Wed May 29 15:59:56 2002// +/init_pot1.c/1.1.1.1/Wed May 29 15:59:56 2002// +/init_pot1.m/1.1.1.1/Wed May 29 15:59:56 2002// +D diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Repository new file mode 100644 index 00000000..eb323e83 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/inference/static/@jtree_sparse_inf_engine/old diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Root b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/collect_evidence.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/collect_evidence.c new file mode 100644 index 00000000..3e6d35c7 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/collect_evidence.c @@ -0,0 +1,635 @@ +/* C mex for collect_evidence.c in @jtree_sparse_inf_engine directory */ +/* File enter_evidence.m in directory @jtree_sparse_inf_engine call it*/ + +/******************************************/ +/* collect_evidence has 3 input & 2 output*/ +/* engine */ +/* clpot */ +/* seppot */ +/* */ +/* clpot */ +/* seppot */ +/******************************************/ + +#include <math.h> +#include <search.h> +#include "mex.h" + +int compare(const void* src1, const void* src2){ + int i1 = *(int*)src1 ; + int i2 = *(int*)src2 ; + return i1-i2 ; +} + +void ind_subv(int index, const int *cumprod, int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i<n; i++){ + index += subv[i] * cumprod[i]; + } + return index; +} + +void compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){ + int i, j; + int *eff_cumprod, *subv, *diffsize, *diff_cumprod; + + subv = malloc(diffdim * sizeof(int)); + eff_cumprod = malloc(diffdim * sizeof(int)); + diffsize = malloc(diffdim * sizeof(int)); + diff_cumprod = malloc(diffdim * sizeof(int)); + for(i=0; i<diffdim; i++){ + eff_cumprod[i] = bCumprod[dmask[i]]; + diffsize[i] = (int)pbSize[dmask[i]]; + } + diff_cumprod[0] = 1; + for(i=0; i<diffdim-1; i++){ + diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i]; + } + for(i=0; i<ND; i++){ + ind_subv(i, diff_cumprod, diffdim, subv); + weight[i] = 0; + for(j=0; j<diffdim; j++){ + weight[i] += eff_cumprod[j] * subv[j]; + } + } + free(eff_cumprod); + free(subv); + free(diffsize); + free(diff_cumprod); +} + +mxArray* convert_table_to_sparse(const double *bT, const int *index, const int nzCounts, const int N){ + mxArray *spTable; + int i, *irs, *jcs; + double *sr; + + spTable = mxCreateSparse(N, 1, nzCounts, mxREAL); + sr = mxGetPr(spTable); + irs = mxGetIr(spTable); + jcs = mxGetJc(spTable); + + jcs[0] = 0; + jcs[1] = nzCounts; + + for(i=0; i<nzCounts; i++){ + sr[i] = bT[i]; + irs[i] = index[i]; + } + return spTable; +} + +mxArray* convert_ill_table_to_sparse(const double *Table, const int *sequence, const int nzCounts, const int N){ + mxArray *spTable; + int i, temp, *irs, *jcs, count=0; + double *sr; + + spTable = mxCreateSparse(N, 1, nzCounts, mxREAL); + sr = mxGetPr(spTable); + irs = mxGetIr(spTable); + jcs = mxGetJc(spTable); + + jcs[0] = 0; + jcs[1] = nzCounts; + + for(i=0; i<nzCounts; i++){ + irs[i] = sequence[count]; + count++; + temp = sequence[count]; + sr[i] = Table[temp]; + count++; + } + return spTable; +} + +void multiply_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex, nzCounts=0; + int *samemask, *diffmask, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *sequence, *weight; + double *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *spr, *bpr; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + if(sdim == 0){ + pTemp = mxCreateSparse(NB, 1, NB, mxREAL); + mxSetField(bigPot, 0, "T", pTemp); + bpr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + sjc[0] = 0; + sjc[1] = NB; + for(i=0; i<NB; i++){ + bpr[i] = *spr; + sir[i] = i; + } + return; + } + + NS = 1; + for(i=0; i<sdim; i++){ + NS *= (int)psSize[i]; + } + ND = NB / NS; + + if(ND == 1){ + pTemp1 = mxGetField(smallPot, 0, "T"); + pTemp = mxDuplicateArray(pTemp1); + mxSetField(bigPot, 0, "T", pTemp); + return; + } + + + NZB = ND * NZS; + + diffdim = bdim - sdim; + sequence = malloc(NZB * 2 * sizeof(int)); + bigTable = malloc(NZB * sizeof(double)); + samemask = malloc(sdim * sizeof(int)); + diffmask = malloc(diffdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + weight = malloc(ND * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + count1 = 0; + for(i=0; i<bdim; i++){ + match = 0; + for(j=0; j<sdim; j++){ + if(pbDomain[i] == psDomain[j]){ + samemask[count] = i; + match = 1; + count++; + break; + } + } + if(match == 0){ + diffmask[count1] = i; + count1++; + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + count = 0; + compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim); + for(i=0; i<NZS; i++){ + sindex = sir[i]; + ind_subv(sindex, sCumprod, sdim, ssubv); + temp = 0; + for(j=0; j<sdim; j++){ + temp += ssubv[j] * bCumprod[samemask[j]]; + } + for(j=0; j<ND; j++){ + bindex = weight[j] + temp; + bigTable[nzCounts] = spr[i]; + sequence[count] = bindex; + count++; + sequence[count] = nzCounts; + nzCounts++; + count++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + qsort(sequence, nzCounts, sizeof(int) * 2, compare); + pTemp = convert_ill_table_to_sparse(bigTable, sequence, nzCounts, NB); + mxSetField(bigPot, 0, "T", pTemp); + + free(sequence); + free(bigTable); + free(samemask); + free(diffmask); + free(bCumprod); + free(sCumprod); + free(weight); + free(ssubv); +} + +void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0; + int *mask, *index, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value; + mxArray *pTemp; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + if(sdim == 0){ + for(i=0; i<NZB; i++){ + bpr[i] *= *spr; + } + return; + } + + bigTable = malloc(NZB * sizeof(double)); + index = malloc(NZB * sizeof(double)); + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + for(i=0; i<NZB; i++){ + bigTable[i] = 0; + } + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + value = bpr[i]; + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); + if(result){ + position = result - sir; + value *= spr[position]; + bigTable[nzCounts] = value; + index[nzCounts] = bindex; + nzCounts++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB); + mxSetField(bigPot, 0, "T", pTemp); + + free(bigTable); + free(index); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + +mxArray* marginal_null_to_spPot(const mxArray *bigPot, const mxArray *sDomain, const int maximize){ + int i, j, count, bdim, sdim, NB, NS, ND; + int *mask, *sir, *sjc; + double *pbDomain, *psDomain, *pbSize, *psSize, *spr; + mxArray *pTemp, *smallPot; + const char *field_names[] = {"domain", "T", "sizes"}; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + psDomain = mxGetPr(sDomain); + sdim = mxGetNumberOfElements(sDomain); + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + + smallPot = mxCreateStructMatrix(1, 1, 3, field_names); + pTemp = mxDuplicateArray(sDomain); + mxSetField(smallPot, 0, "domain", pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + if(sdim == 0){ + pTemp = mxCreateSparse(1, 1, 1, mxREAL); + mxSetField(smallPot, 0, "T", pTemp); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + *spr = 0; + *sir = 0; + sjc[0] = 0; + sjc[1] = 1; + if(maximize) *spr = 1; + else *spr = NB; + + pTemp = mxCreateDoubleMatrix(1, 1, mxREAL); + *mxGetPr(pTemp) = 1; + mxSetField(smallPot, 0, "sizes", pTemp); + return smallPot; + } + + mask = malloc(sdim * sizeof(int)); + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + pTemp = mxCreateDoubleMatrix(1, count, mxREAL); + psSize = mxGetPr(pTemp); + NS = 1; + for(i=0; i<count; i++){ + psSize[i] = pbSize[mask[i]]; + NS *= (int)psSize[i]; + } + mxSetField(smallPot, 0, "sizes", pTemp); + + ND = NB / NS; + + pTemp = mxCreateSparse(NS, 1, NS, mxREAL); + mxSetField(smallPot, 0, "T", pTemp); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + if(maximize){ + for(i=0; i<NS; i++){ + spr[i] = 1; + sir[i] = i; + } + } + else{ + for(i=0; i<NS; i++){ + spr[i] = ND; + sir[i] = i; + } + } + sjc[0] = 0; + sjc[1] = NS; + + free(mask); + return smallPot; +} + +mxArray* marginal_spPot_to_spPot(const mxArray *bigPot, const mxArray *sDomain, const int maximize){ + int i, j, count, bdim, sdim, NB, NS, NZB, position, bindex, sindex, nzCounts=0; + int *mask, *sequence, *result, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *sTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr; + mxArray *pTemp, *smallPot; + const char *field_names[] = {"domain", "T", "sizes"}; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + psDomain = mxGetPr(sDomain); + sdim = mxGetNumberOfElements(sDomain); + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + smallPot = mxCreateStructMatrix(1, 1, 3, field_names); + pTemp = mxDuplicateArray(sDomain); + mxSetField(smallPot, 0, "domain", pTemp); + + if(sdim == 0){ + pTemp = mxCreateSparse(1, 1, 1, mxREAL); + mxSetField(smallPot, 0, "T", pTemp); + spr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + *spr = 0; + *bir = 0; + bjc[0] = 0; + bjc[1] = 1; + if(maximize){ + for(i=0; i<NZB; i++){ + *spr = (*spr < bpr[i])? bpr[i] : *spr; + } + } + else{ + for(i=0; i<NZB; i++){ + *spr += bpr[i]; + } + } + + pTemp = mxCreateDoubleMatrix(1, 1, mxREAL); + *mxGetPr(pTemp) = 1; + mxSetField(smallPot, 0, "sizes", pTemp); + return smallPot; + } + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + mask = malloc(sdim * sizeof(int)); + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + pTemp = mxCreateDoubleMatrix(1, count, mxREAL); + psSize = mxGetPr(pTemp); + NS = 1; + for(i=0; i<count; i++){ + psSize[i] = pbSize[mask[i]]; + NS *= (int)psSize[i]; + } + mxSetField(smallPot, 0, "sizes", pTemp); + + + sTable = malloc(NZB * sizeof(double)); + sequence = malloc(NZB * 2 * sizeof(double)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + for(i=0; i<NZB; i++)sTable[i] = 0; + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + count = 0; + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sequence, nzCounts, sizeof(int)*2, compare); + if(result){ + position = (result - sequence) / 2; + if(maximize) + sTable[position] = (sTable[position] < bpr[i]) ? bpr[i] : sTable[position]; + else sTable[position] += bpr[i]; + } + else { + if(maximize) + sTable[nzCounts] = (sTable[nzCounts] < bpr[i]) ? bpr[i] : sTable[nzCounts]; + else sTable[nzCounts] += bpr[i]; + sequence[count] = sindex; + count++; + sequence[count] = nzCounts; + nzCounts++; + count++; + } + } + + qsort(sequence, nzCounts, sizeof(int) * 2, compare); + pTemp = convert_ill_table_to_sparse(sTable, sequence, nzCounts, NS); + mxSetField(smallPot, 0, "T", pTemp); + + free(sTable); + free(sequence); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); + + return smallPot; +} + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, n, p, np, pn, loop, loops, nCliques, temp, maximize; + int *collect_order; + double *pr, *pr1; + mxArray *pTemp, *pTemp1, *pPostP, *pClpot, *pSeppot, *pSeparator; + + pTemp = mxGetField(prhs[0], 0, "cliques"); + nCliques = mxGetNumberOfElements(pTemp); + loops = nCliques - 1; + pTemp = mxGetField(prhs[0], 0, "maximize"); + maximize = (int)mxGetScalar(pTemp); + pSeparator = mxGetField(prhs[0], 0, "separator"); + + collect_order = malloc(2 * loops * sizeof(int)); + + pTemp = mxGetField(prhs[0], 0, "postorder"); + pr = mxGetPr(pTemp); + pPostP = mxGetField(prhs[0], 0, "postorder_parents"); + for(i=0; i<loops; i++){ + temp = (int)pr[i] - 1; + pTemp = mxGetCell(pPostP, temp); + pr1 = mxGetPr(pTemp); + collect_order[i] = (int)pr1[0] - 1; + collect_order[i+loops] = temp; + } + + plhs[0] = mxDuplicateArray(prhs[1]); + plhs[1] = mxDuplicateArray(prhs[2]); + + for(loop=0; loop<loops; loop++){ + p = collect_order[loop]; + n = collect_order[loop+loops]; + np = p * nCliques + n; + pn = n * nCliques + p; + pClpot = mxGetCell(plhs[0], n); + pTemp1 = mxGetField(pClpot, 0, "T"); + pTemp = mxGetCell(pSeparator, pn); + if(pTemp1) + pSeppot = marginal_spPot_to_spPot(pClpot, pTemp, maximize); + else pSeppot = marginal_null_to_spPot(pClpot, pTemp, maximize); + mxSetCell(plhs[1], pn, pSeppot); + + pClpot = mxGetCell(plhs[0], p); + pTemp1 = mxGetField(pClpot, 0, "T"); + if(pTemp1) + multiply_spPot_by_spPot(pClpot, pSeppot); + else multiply_null_by_spPot(pClpot, pSeppot); + } + free(collect_order); +} + + + + + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/distribute_evidence.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/distribute_evidence.c new file mode 100644 index 00000000..3d8ec66b --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/distribute_evidence.c @@ -0,0 +1,613 @@ +/* C mex for distribute_evidence.c in @jtree_sparse_inf_engine directory*/ +/* File enter_evidence.m in directory @jtree_sparse_inf_engine call it */ + +/*********************************************/ +/* distribute_evidence has 3 input & 2 output*/ +/* engine */ +/* clpot */ +/* seppot */ +/* */ +/* clpot */ +/* seppot */ +/*********************************************/ + +#include "mex.h" + +#include <math.h> +#include <search.h> +#include "mex.h" + +int compare(const void* src1, const void* src2){ + int i1 = *(int*)src1 ; + int i2 = *(int*)src2 ; + return i1-i2 ; +} + +void ind_subv(int index, const int *cumprod, int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i<n; i++){ + index += subv[i] * cumprod[i]; + } + return index; +} + +void compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){ + int i, j; + int *eff_cumprod, *subv, *diffsize, *diff_cumprod; + + subv = malloc(diffdim * sizeof(int)); + eff_cumprod = malloc(diffdim * sizeof(int)); + diffsize = malloc(diffdim * sizeof(int)); + diff_cumprod = malloc(diffdim * sizeof(int)); + for(i=0; i<diffdim; i++){ + eff_cumprod[i] = bCumprod[dmask[i]]; + diffsize[i] = (int)pbSize[dmask[i]]; + } + diff_cumprod[0] = 1; + for(i=0; i<diffdim-1; i++){ + diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i]; + } + for(i=0; i<ND; i++){ + ind_subv(i, diff_cumprod, diffdim, subv); + weight[i] = 0; + for(j=0; j<diffdim; j++){ + weight[i] += eff_cumprod[j] * subv[j]; + } + } + free(eff_cumprod); + free(subv); + free(diffsize); + free(diff_cumprod); +} + +mxArray* convert_table_to_sparse(const double *bT, const int *index, const int nzCounts, const int N){ + mxArray *spTable; + int i, *irs, *jcs; + double *sr; + + spTable = mxCreateSparse(N, 1, nzCounts, mxREAL); + sr = mxGetPr(spTable); + irs = mxGetIr(spTable); + jcs = mxGetJc(spTable); + + jcs[0] = 0; + jcs[1] = nzCounts; + + for(i=0; i<nzCounts; i++){ + sr[i] = bT[i]; + irs[i] = index[i]; + } + return spTable; +} + +mxArray* convert_ill_table_to_sparse(const double *Table, const int *sequence, const int nzCounts, const int N){ + mxArray *spTable; + int i, temp, *irs, *jcs, count=0; + double *sr; + + spTable = mxCreateSparse(N, 1, nzCounts, mxREAL); + sr = mxGetPr(spTable); + irs = mxGetIr(spTable); + jcs = mxGetJc(spTable); + + jcs[0] = 0; + jcs[1] = nzCounts; + + for(i=0; i<nzCounts; i++){ + irs[i] = sequence[count]; + count++; + temp = sequence[count]; + sr[i] = Table[temp]; + count++; + } + return spTable; +} + +void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0; + int *mask, *index, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value; + mxArray *pTemp; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + if(sdim == 0){ + for(i=0; i<NZB; i++){ + bpr[i] *= *spr; + } + return; + } + + bigTable = malloc(NZB * sizeof(double)); + index = malloc(NZB * sizeof(double)); + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + value = bpr[i]; + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); + if(result){ + position = result - sir; + value *= spr[position]; + bigTable[nzCounts] = value; + index[nzCounts] = bindex; + nzCounts++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB); + mxSetField(bigPot, 0, "T", pTemp); + + free(bigTable); + free(index); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + +void marginal_spPot_to_spPot(const mxArray *bigPot, mxArray *smallPot, const int maximize){ + int i, j, count, bdim, sdim, NB, NS, NZB, position, bindex, sindex, nzCounts=0; + int *mask, *sequence, *result, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *sTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr; + mxArray *pTemp; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + NS = 1; + for(i=0; i<sdim; i++){ + NS *= (int)psSize[i]; + } + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + if(sdim == 0){ + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + *spr = 0; + if(maximize){ + for(i=0; i<NZB; i++){ + *spr = (*spr < bpr[i])? bpr[i] : *spr; + } + } + else{ + for(i=0; i<NZB; i++){ + *spr += bpr[i]; + } + } + return; + } + + mask = malloc(sdim * sizeof(int)); + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + + sTable = malloc(NZB * sizeof(double)); + sequence = malloc(NZB * 2 * sizeof(double)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + for(i=0; i<NZB; i++){ + sTable[i] = 0; + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + count = 0; + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sequence, nzCounts, sizeof(int)*2, compare); + if(result){ + position = (result - sequence) / 2; + if(maximize) + sTable[position] = (sTable[position] < bpr[i]) ? bpr[i] : sTable[position]; + else sTable[position] += bpr[i]; + } + else { + if(maximize) + sTable[nzCounts] = (sTable[nzCounts] < bpr[i]) ? bpr[i] : sTable[nzCounts]; + else sTable[nzCounts] += bpr[i]; + sequence[count] = sindex; + count++; + sequence[count] = nzCounts; + nzCounts++; + count++; + } + } + + pTemp = mxGetField(smallPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + qsort(sequence, nzCounts, sizeof(int) * 2, compare); + pTemp = convert_ill_table_to_sparse(sTable, sequence, nzCounts, NS); + mxSetField(smallPot, 0, "T", pTemp); + + free(sTable); + free(sequence); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + +void divide_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex; + int *samemask, *diffmask, *rir, *rjc, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *weight; + double *pbDomain, *psDomain, *pbSize, *psSize, *rpr, *spr, value; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + if(sdim == 0){ + pTemp = mxCreateSparse(NB, 1, NB, mxREAL); + mxSetField(bigPot, 0, "T", pTemp); + rpr = mxGetPr(pTemp); + rir = mxGetIr(pTemp); + rjc = mxGetJc(pTemp); + rjc[0] = 0; + rjc[1] = NB; + value = *spr; + if(value == 0) value = 1; + for(i=0; i<NB; i++){ + rpr[i] = 1 / value; + rir[i] = i; + } + return; + } + + NS = 1; + for(i=0; i<sdim; i++){ + NS *= (int)psSize[i]; + } + ND = NB / NS; + + + pTemp = mxCreateSparse(NB, 1, NB, mxREAL); + rpr = mxGetPr(pTemp); + rir = mxGetIr(pTemp); + rjc = mxGetJc(pTemp); + rjc[0] = 0; + rjc[1] = NB; + for(i=0; i<NB; i++){ + rpr[i] = 1; + rir[i] = i; + } + + NZB = ND * NZS; + + diffdim = bdim - sdim; + samemask = malloc(sdim * sizeof(int)); + diffmask = malloc(diffdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + weight = malloc(ND * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + count1 = 0; + for(i=0; i<bdim; i++){ + match = 0; + for(j=0; j<sdim; j++){ + if(pbDomain[i] == psDomain[j]){ + samemask[count] = i; + match = 1; + count++; + break; + } + } + if(match == 0){ + diffmask[count1] = i; + count1++; + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + count = 0; + compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim); + for(i=0; i<NZS; i++){ + sindex = sir[i]; + ind_subv(sindex, sCumprod, sdim, ssubv); + temp = 0; + for(j=0; j<sdim; j++){ + temp += ssubv[j] * bCumprod[samemask[j]]; + } + for(j=0; j<ND; j++){ + bindex = weight[j] + temp; + rpr[bindex] = 1 / (spr[i]); + } + } + + pTemp1 = mxGetField(bigPot, 0, "T"); + if(pTemp1)mxDestroyArray(pTemp1); + mxSetField(bigPot, 0, "T", pTemp); + + free(samemask); + free(diffmask); + free(bCumprod); + free(sCumprod); + free(weight); + free(ssubv); +} + +void divide_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex; + int *mask, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp1 = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp1); + bir = mxGetIr(pTemp1); + bjc = mxGetJc(pTemp1); + NZB = bjc[1]; + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + if(sdim == 0){ + value = *spr; + if(value == 0)value = 1; + for(i=0; i<NZB; i++){ + bpr[i] /= value; + } + return; + } + + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); + if(result){ + position = result - sir; + bpr[i] /= spr[position]; + } + } + + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, loop, loops, nCliques, temp, count, parent, child, maximize, *distribute_order; + double *pr, *pr1; + mxArray *pTemp, *pPreCh, *pClpot, *pSeppot; + + pTemp = mxGetField(prhs[0], 0, "cliques"); + nCliques = mxGetNumberOfElements(pTemp); + loops = nCliques - 1; + pTemp = mxGetField(prhs[0], 0, "maximize"); + maximize = (int)mxGetScalar(pTemp); + + distribute_order = malloc(2 * loops * sizeof(int)); + pTemp = mxGetField(prhs[0], 0, "preorder"); + pr = mxGetPr(pTemp); + pPreCh = mxGetField(prhs[0], 0, "preorder_children"); + count = 0; + for(i=0; i<nCliques; i++){ + temp = (int)pr[i] - 1; + pTemp = mxGetCell(pPreCh, temp); + pr1 = mxGetPr(pTemp); + loop = mxGetNumberOfElements(pTemp); + for(j=0; j<loop; j++){ + distribute_order[count] = temp; + distribute_order[count + loops] = (int)pr1[j] - 1; + count++; + } + } + + plhs[0] = mxDuplicateArray(prhs[1]); + plhs[1] = mxDuplicateArray(prhs[2]); + + for(loop=0; loop<loops; loop++){ + parent = distribute_order[loop]; + child = distribute_order[loop+loops]; + i = nCliques * child + parent; + pClpot = mxGetCell(plhs[0], child); + pTemp = mxGetField(pClpot, 0, "T"); + pSeppot = mxGetCell(plhs[1], i); + if(pTemp) + divide_spPot_by_spPot(pClpot, pSeppot); + else divide_null_by_spPot(pClpot, pSeppot); + + pClpot = mxGetCell(plhs[0], parent); + marginal_spPot_to_spPot(pClpot, pSeppot, maximize); + mxSetCell(plhs[1], i, pSeppot); + + pClpot = mxGetCell(plhs[0], child); + multiply_spPot_by_spPot(pClpot, pSeppot); + } + free(distribute_order); +} diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot.c new file mode 100644 index 00000000..5d0ed8a3 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot.c @@ -0,0 +1,637 @@ +/* C mex init_pot for in @jtree_sparse_inf_engine directory */ +/* The file enter_evidence.m in directory @jtree_sparse_inf_engine call it*/ + +/**************************************/ +/* init_pot.c has 6 input & 2 output */ +/* engine */ +/* clqs */ +/* pots */ +/* pot_type */ +/* onodes */ +/* ndx */ +/* */ +/* clpot */ +/* seppot */ +/**************************************/ +#include <math.h> +#include <search.h> +#include "mex.h" + +int compare(const void* src1, const void* src2){ + int i1 = *(int*)src1 ; + int i2 = *(int*)src2 ; + return i1-i2 ; +} + +void ind_subv(int index, const int *cumprod, int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i<n; i++){ + index += subv[i] * cumprod[i]; + } + return index; +} + +void compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){ + int i, j; + int *eff_cumprod, *subv, *diffsize, *diff_cumprod; + + subv = malloc(diffdim * sizeof(int)); + eff_cumprod = malloc(diffdim * sizeof(int)); + diffsize = malloc(diffdim * sizeof(int)); + diff_cumprod = malloc(diffdim * sizeof(int)); + for(i=0; i<diffdim; i++){ + eff_cumprod[i] = bCumprod[dmask[i]]; + diffsize[i] = (int)pbSize[dmask[i]]; + } + diff_cumprod[0] = 1; + for(i=0; i<diffdim-1; i++){ + diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i]; + } + for(i=0; i<ND; i++){ + ind_subv(i, diff_cumprod, diffdim, subv); + weight[i] = 0; + for(j=0; j<diffdim; j++){ + weight[i] += eff_cumprod[j] * subv[j]; + } + } + free(eff_cumprod); + free(subv); + free(diffsize); + free(diff_cumprod); +} + +mxArray* convert_to_sparse(const double *table, const int NB, const int counts){ + mxArray *spTable; + int i, k, *ir, *jc; + double *sr; + + spTable = mxCreateSparse(NB, 1, counts, mxREAL); + sr = mxGetPr(spTable); + ir = mxGetIr(spTable); + jc = mxGetJc(spTable); + + k = 0; + jc[0] = 0; + jc[1] = counts; + for(i=0; i<NB; i++){ + if(table[i] != 0.0){ + sr[k] = table[i]; + ir[k] = i; + k++; + } + } + + return spTable; +} + +mxArray* convert_table_to_sparse(const double *bT, const int *index, const int nzCounts, const int NB){ + mxArray *spTable; + int i, *irs, *jcs; + double *sr; + + spTable = mxCreateSparse(NB, 1, nzCounts, mxREAL); + sr = mxGetPr(spTable); + irs = mxGetIr(spTable); + jcs = mxGetJc(spTable); + + jcs[0] = 0; + jcs[1] = nzCounts; + + for(i=0; i<nzCounts; i++){ + sr[i] = bT[i]; + irs[i] = index[i]; + } + return spTable; +} + +mxArray* convert_ill_table_to_sparse(const double *bigTable, const int *sequence, const int nzCounts, const int NB){ + mxArray *spTable; + int i, temp, *irs, *jcs, count=0; + double *sr; + + spTable = mxCreateSparse(NB, 1, nzCounts, mxREAL); + sr = mxGetPr(spTable); + irs = mxGetIr(spTable); + jcs = mxGetJc(spTable); + + jcs[0] = 0; + jcs[1] = nzCounts; + + for(i=0; i<nzCounts; i++){ + irs[i] = sequence[count]; + count++; + temp = sequence[count]; + sr[i] = bigTable[temp]; + count++; + } + return spTable; +} + +void multiply_null_by_fuPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, NB, NS, siz_b, siz_s, ndim, nzCounts=0; + int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *jc; + double *pbDomain, *psDomain, *pbSize, *psSize, *bTable, *sTable, value; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + siz_b = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + siz_s = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<siz_b; i++){ + NB *= (int)pbSize[i]; + } + NS = 1; + for(i=0; i<siz_s; i++){ + NS *= (int)psSize[i]; + } + + pTemp = mxGetField(smallPot, 0, "T"); + sTable = mxGetPr(pTemp); + bTable = malloc(NB * sizeof(double)); + for(i=0; i<NB; i++){ + bTable[i] = 0; + } + + if(NS == 1){ + value = *sTable; + for(i=0; i<NB; i++){ + bTable[i] = value; + } + nzCounts = NB; + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + pTemp = convert_to_sparse(bTable, NB, NB); + mxSetField(bigPot, 0, "T", pTemp); + free(bTable); + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++){ + bTable[i] = sTable[i]; + if(sTable[i] != 0) nzCounts++; + } + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + pTemp = convert_to_sparse(bTable, NB, nzCounts); + mxSetField(bigPot, 0, "T", pTemp); + free(bTable); + return; + } + + mask = malloc(siz_s * sizeof(int)); + count = 0; + for(i=0; i<siz_s; i++){ + for(j=0; j<siz_b; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + ndim = siz_b; + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + sx[i] = (int)pbSize[i]; + sy[i] = 1; + } + for(i=0; i<count; i++){ + sy[mask[i]] = sx[mask[i]]; + } + + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i = 0; i < ndim; i++){ + subs[i] = 0; + s[i] = sx[i] - 1; + } + + for(i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + for(j=0; j<NB; j++){ + bTable[j] = *sTable; + if(*sTable != 0.0) nzCounts++; + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + sTable -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + sTable += cpsy[i]; + break; + } + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + pTemp = convert_to_sparse(bTable, NB, nzCounts); + mxSetField(bigPot, 0, "T", pTemp); + pTemp1 = mxGetField(bigPot, 0, "T"); + jc = mxGetJc(pTemp1); + + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); + free(mask); + free(bTable); +} + +void multiply_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex, nzCounts=0; + int *samemask, *diffmask, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *sequence, *weight; + double *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *spr; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + NS = 1; + for(i=0; i<sdim; i++){ + NS *= (int)psSize[i]; + } + ND = NB / NS; + + if(ND == 1){ + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + pTemp1 = mxGetField(smallPot, 0, "T"); + pTemp = mxDuplicateArray(pTemp1); + mxSetField(bigPot, 0, "T", pTemp); + return; + } + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + NZB = ND * NZS; + + diffdim = bdim - sdim; + sequence = malloc(NZB * 2 * sizeof(int)); + bigTable = malloc(NZB * sizeof(double)); + samemask = malloc(sdim * sizeof(int)); + diffmask = malloc(diffdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + weight = malloc(ND * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + count1 = 0; + for(i=0; i<bdim; i++){ + match = 0; + for(j=0; j<sdim; j++){ + if(pbDomain[i] == psDomain[j]){ + samemask[count] = i; + match = 1; + count++; + break; + } + } + if(match == 0){ + diffmask[count1] = i; + count1++; + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + count = 0; + compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim); + for(i=0; i<NZS; i++){ + sindex = sir[i]; + ind_subv(sindex, sCumprod, sdim, ssubv); + temp = 0; + for(j=0; j<sdim; j++){ + temp += ssubv[j] * bCumprod[samemask[j]]; + } + for(j=0; j<ND; j++){ + bindex = weight[j] + temp; + bigTable[nzCounts] = spr[i]; + sequence[count] = bindex; + count++; + sequence[count] = nzCounts; + nzCounts++; + count++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + qsort(sequence, nzCounts, sizeof(int) * 2, compare); + pTemp = convert_ill_table_to_sparse(bigTable, sequence, nzCounts, NB); + mxSetField(bigPot, 0, "T", pTemp); + + free(sequence); + free(bigTable); + free(samemask); + free(diffmask); + free(bCumprod); + free(sCumprod); + free(weight); + free(ssubv); +} + +void multiply_spPot_by_fuPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, bdim, sdim, NB, NZB, bindex, sindex, nzCounts=0; + int *mask, *index, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value; + mxArray *pTemp; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + + bigTable = malloc(NZB * sizeof(double)); + index = malloc(NZB * sizeof(double)); + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + for(i=0; i<NZB; i++){ + bigTable[i] = 0; + } + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + value = spr[sindex]; + if(value != 0){ + bigTable[nzCounts] = bpr[i] * value; + index[nzCounts] = bindex; + nzCounts++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB); + mxSetField(bigPot, 0, "T", pTemp); + + free(bigTable); + free(index); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + +void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0; + int *mask, *index, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value; + mxArray *pTemp; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + bigTable = malloc(NZB * sizeof(double)); + index = malloc(NZB * sizeof(double)); + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + for(i=0; i<NZB; i++){ + bigTable[i] = 0; + } + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + value = bpr[i]; + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); + if(result){ + position = result - sir; + value *= spr[position]; + bigTable[nzCounts] = value; + index[nzCounts] = bindex; + nzCounts++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB); + mxSetField(bigPot, 0, "T", pTemp); + + free(bigTable); + free(index); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, c, loop, nNodes, nCliques, ndomain, dims[2]; + double *pClqs, *pr, *pt, *pSize; + mxArray *pTemp, *pTemp1, *pStruct, *pCliques, *pBigpot, *pSmallpot; + const char *field_names[] = {"domain", "T", "sizes"}; + + nNodes = mxGetNumberOfElements(prhs[1]); + pCliques = mxGetField(prhs[0], 0, "cliques"); + nCliques = mxGetNumberOfElements(pCliques); + pTemp = mxGetField(prhs[0], 0, "eff_node_sizes"); + pSize = mxGetPr(pTemp); + + plhs[0] = mxCreateCellArray(1, &nCliques); + for(i=0; i<nCliques; i++){ + pStruct = mxCreateStructMatrix(1, 1, 3, field_names); + mxSetCell(plhs[0], i, pStruct); + pTemp = mxGetCell(pCliques, i); + ndomain = mxGetNumberOfElements(pTemp); + pt = mxGetPr(pTemp); + pTemp1 = mxDuplicateArray(pTemp); + mxSetField(pStruct, 0, "domain", pTemp1); + + pTemp = mxCreateDoubleMatrix(1, ndomain, mxREAL); + mxSetField(pStruct, 0, "sizes", pTemp); + pr = mxGetPr(pTemp); + for(j=0; j<ndomain; j++){ + pr[j] = pSize[(int)pt[j]-1]; + } + } + + pClqs = mxGetPr(prhs[1]); + for(loop=0; loop<nNodes; loop++){ + c = (int)pClqs[loop] - 1; + pSmallpot = mxGetCell(prhs[2], loop); + pTemp = mxGetField(pSmallpot, 0, "T"); + pBigpot = mxGetCell(plhs[0], c); + pTemp1 = mxGetField(pBigpot, 0, "T"); + if(pTemp1){ + if(mxIsSparse(pTemp)) + multiply_spPot_by_spPot(pBigpot, pSmallpot); + else multiply_spPot_by_fuPot(pBigpot, pSmallpot); + } + else{ + if(mxIsSparse(pTemp)) + multiply_null_by_spPot(pBigpot, pSmallpot); + else multiply_null_by_fuPot(pBigpot, pSmallpot); + } + } + + dims[0] = nCliques; + dims[1] = nCliques; + plhs[1] = mxCreateCellArray(2, dims); +} + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.c b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.c new file mode 100644 index 00000000..b3a6a66d --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.c @@ -0,0 +1,636 @@ +/* C mex init_pot for in @jtree_sparse_inf_engine directory */ +/* The file enter_evidence.m in directory @jtree_sparse_inf_engine call it*/ + +/**************************************/ +/* init_pot.c has 6 input & 2 output */ +/* engine */ +/* clqs */ +/* pots */ +/* pot_type */ +/* onodes */ +/* ndx */ +/* */ +/* clpot */ +/* seppot */ +/**************************************/ +#include <math.h> +#include <search.h> +#include "mex.h" + +int compare(const void* src1, const void* src2){ + int i1 = *(int*)src1 ; + int i2 = *(int*)src2 ; + return i1-i2 ; +} + +void ind_subv(int index, const int *cumprod, int n, int *bsubv){ + int i; + + for (i = n-1; i >= 0; i--) { + bsubv[i] = ((int)floor(index / cumprod[i])); + index = index % cumprod[i]; + } +} + +int subv_ind(const int n, const int *cumprod, const int *subv){ + int i, index=0; + + for(i=0; i<n; i++){ + index += subv[i] * cumprod[i]; + } + return index; +} + +void compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){ + int i, j; + int *eff_cumprod, *subv, *diffsize, *diff_cumprod; + + subv = malloc(diffdim * sizeof(int)); + eff_cumprod = malloc(diffdim * sizeof(int)); + diffsize = malloc(diffdim * sizeof(int)); + diff_cumprod = malloc(diffdim * sizeof(int)); + for(i=0; i<diffdim; i++){ + eff_cumprod[i] = bCumprod[dmask[i]]; + diffsize[i] = (int)pbSize[dmask[i]]; + } + diff_cumprod[0] = 1; + for(i=0; i<diffdim-1; i++){ + diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i]; + } + for(i=0; i<ND; i++){ + ind_subv(i, diff_cumprod, diffdim, subv); + weight[i] = 0; + for(j=0; j<diffdim; j++){ + weight[i] += eff_cumprod[j] * subv[j]; + } + } + free(eff_cumprod); + free(subv); + free(diffsize); + free(diff_cumprod); +} + +void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){ + double *ptr; + void *newptr; + int *ir, *jc; + int nbytes; + + if(new_nzmax == old_nzmax) return; + nbytes = new_nzmax * sizeof(*ptr); + ptr = mxGetPr(spArray); + newptr = mxRealloc(ptr, nbytes); + mxSetPr(spArray, newptr); + nbytes = new_nzmax * sizeof(*ir); + ir = mxGetIr(spArray); + newptr = mxRealloc(ir, nbytes); + mxSetIr(spArray, newptr); + jc = mxGetJc(spArray); + jc[0] = 0; + jc[1] = new_nzmax; + mxSetNzmax(spArray, new_nzmax); +} + +mxArray* convert_table_to_sparse(const double *bT, const int *index, const int nzCounts, const int NB){ + mxArray *spTable; + int i, *irs, *jcs; + double *sr; + + spTable = mxCreateSparse(NB, 1, nzCounts, mxREAL); + sr = mxGetPr(spTable); + irs = mxGetIr(spTable); + jcs = mxGetJc(spTable); + + jcs[0] = 0; + jcs[1] = nzCounts; + + for(i=0; i<nzCounts; i++){ + sr[i] = bT[i]; + irs[i] = index[i]; + } + return spTable; +} + +mxArray* convert_ill_table_to_sparse(const double *bigTable, const int *sequence, const int nzCounts, const int NB){ + mxArray *spTable; + int i, temp, *irs, *jcs, count=0; + double *sr; + + spTable = mxCreateSparse(NB, 1, nzCounts, mxREAL); + sr = mxGetPr(spTable); + irs = mxGetIr(spTable); + jcs = mxGetJc(spTable); + + jcs[0] = 0; + jcs[1] = nzCounts; + + for(i=0; i<nzCounts; i++){ + irs[i] = sequence[count]; + count++; + temp = sequence[count]; + sr[i] = bigTable[temp]; + count++; + } + return spTable; +} + +void multiply_null_by_fuPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, NB, NS, siz_b, siz_s, ndim, nzCounts=0; + int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *bir, *bjc; + double *pbDomain, *psDomain, *pbSize, *psSize, *spr, *bpr, value; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + siz_b = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + siz_s = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<siz_b; i++){ + NB *= (int)pbSize[i]; + } + NS = 1; + for(i=0; i<siz_s; i++){ + NS *= (int)psSize[i]; + } + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + + pTemp1 = mxCreateSparse(NB, 1, NB, mxREAL); + bpr = mxGetPr(pTemp1); + bir = mxGetIr(pTemp1); + bjc = mxGetJc(pTemp1); + bjc[0] = 0; + bjc[1] = NB; + + if(NS == 1){ + value = *spr; + for(i=0; i<NB; i++){ + bpr[i] = value; + bir[i] = i; + } + nzCounts = NB; + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + reset_nzmax(pTemp1, NB, nzCounts); + mxSetField(bigPot, 0, "T", pTemp1); + return; + } + + if(NS == NB){ + for(i=0; i<NB; i++){ + if(spr[i] != 0){ + bpr[nzCounts] = spr[i]; + bir[nzCounts] = i; + nzCounts++; + } + } + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + reset_nzmax(pTemp1, NB, nzCounts); + mxSetField(bigPot, 0, "T", pTemp1); + return; + } + + mask = malloc(siz_s * sizeof(int)); + count = 0; + for(i=0; i<siz_s; i++){ + for(j=0; j<siz_b; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + ndim = siz_b; + sx = (int *)malloc(sizeof(int)*ndim); + sy = (int *)malloc(sizeof(int)*ndim); + for(i=0; i<ndim; i++){ + sx[i] = (int)pbSize[i]; + sy[i] = 1; + } + for(i=0; i<count; i++){ + sy[mask[i]] = sx[mask[i]]; + } + + s = (int *)malloc(sizeof(int)*ndim); + *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1; + subs = (int *)malloc(sizeof(int)*ndim); + cpsy2 = (int *)malloc(sizeof(int)*ndim); + for(i = 0; i < ndim; i++){ + subs[i] = 0; + s[i] = sx[i] - 1; + } + + for(i = 0; i < ndim-1; i++){ + cpsy[i+1] = cpsy[i]*sy[i]--; + cpsy2[i] = cpsy[i]*sy[i]; + } + cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]); + + for(j=0; j<NB; j++){ + if(*spr != 0){ + bpr[nzCounts] = *spr; + bir[nzCounts] = j; + nzCounts++; + } + for(i = 0; i < ndim; i++){ + if(subs[i] == s[i]){ + subs[i] = 0; + if(sy[i]) + spr -= cpsy2[i]; + } + else{ + subs[i]++; + if(sy[i]) + spr += cpsy[i]; + break; + } + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + reset_nzmax(pTemp1, NB, nzCounts); + mxSetField(bigPot, 0, "T", pTemp1); + + free(sx); + free(sy); + free(s); + free(cpsy); + free(subs); + free(cpsy2); + free(mask); +} + +void multiply_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex, nzCounts=0; + int *samemask, *diffmask, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *sequence, *weight; + double *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *spr; + mxArray *pTemp, *pTemp1; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + NS = 1; + for(i=0; i<sdim; i++){ + NS *= (int)psSize[i]; + } + ND = NB / NS; + + if(ND == 1){ + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + pTemp1 = mxGetField(smallPot, 0, "T"); + pTemp = mxDuplicateArray(pTemp1); + mxSetField(bigPot, 0, "T", pTemp); + return; + } + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + NZB = ND * NZS; + + diffdim = bdim - sdim; + sequence = malloc(NZB * 2 * sizeof(int)); + bigTable = malloc(NZB * sizeof(double)); + samemask = malloc(sdim * sizeof(int)); + diffmask = malloc(diffdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + weight = malloc(ND * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + count1 = 0; + for(i=0; i<bdim; i++){ + match = 0; + for(j=0; j<sdim; j++){ + if(pbDomain[i] == psDomain[j]){ + samemask[count] = i; + match = 1; + count++; + break; + } + } + if(match == 0){ + diffmask[count1] = i; + count1++; + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + count = 0; + compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim); + for(i=0; i<NZS; i++){ + sindex = sir[i]; + ind_subv(sindex, sCumprod, sdim, ssubv); + temp = 0; + for(j=0; j<sdim; j++){ + temp += ssubv[j] * bCumprod[samemask[j]]; + } + for(j=0; j<ND; j++){ + bindex = weight[j] + temp; + bigTable[nzCounts] = spr[i]; + sequence[count] = bindex; + count++; + sequence[count] = nzCounts; + nzCounts++; + count++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + qsort(sequence, nzCounts, sizeof(int) * 2, compare); + pTemp = convert_ill_table_to_sparse(bigTable, sequence, nzCounts, NB); + mxSetField(bigPot, 0, "T", pTemp); + + free(sequence); + free(bigTable); + free(samemask); + free(diffmask); + free(bCumprod); + free(sCumprod); + free(weight); + free(ssubv); +} + +void multiply_spPot_by_fuPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, bdim, sdim, NB, NZB, bindex, sindex, nzCounts=0; + int *mask, *index, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value; + mxArray *pTemp; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + + bigTable = malloc(NZB * sizeof(double)); + index = malloc(NZB * sizeof(double)); + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + value = spr[sindex]; + if(value != 0){ + bigTable[nzCounts] = bpr[i] * value; + index[nzCounts] = bindex; + nzCounts++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB); + mxSetField(bigPot, 0, "T", pTemp); + + free(bigTable); + free(index); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + +void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){ + int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0; + int *mask, *index, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv; + double *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value; + mxArray *pTemp; + + pTemp = mxGetField(bigPot, 0, "domain"); + pbDomain = mxGetPr(pTemp); + bdim = mxGetNumberOfElements(pTemp); + pTemp = mxGetField(smallPot, 0, "domain"); + psDomain = mxGetPr(pTemp); + sdim = mxGetNumberOfElements(pTemp); + + pTemp = mxGetField(bigPot, 0, "sizes"); + pbSize = mxGetPr(pTemp); + pTemp = mxGetField(smallPot, 0, "sizes"); + psSize = mxGetPr(pTemp); + + NB = 1; + for(i=0; i<bdim; i++){ + NB *= (int)pbSize[i]; + } + + pTemp = mxGetField(bigPot, 0, "T"); + bpr = mxGetPr(pTemp); + bir = mxGetIr(pTemp); + bjc = mxGetJc(pTemp); + NZB = bjc[1]; + + pTemp = mxGetField(smallPot, 0, "T"); + spr = mxGetPr(pTemp); + sir = mxGetIr(pTemp); + sjc = mxGetJc(pTemp); + NZS = sjc[1]; + + bigTable = malloc(NZB * sizeof(double)); + index = malloc(NZB * sizeof(double)); + mask = malloc(sdim * sizeof(int)); + bCumprod = malloc(bdim * sizeof(int)); + sCumprod = malloc(sdim * sizeof(int)); + bsubv = malloc(bdim * sizeof(int)); + ssubv = malloc(sdim * sizeof(int)); + + for(i=0; i<NZB; i++){ + bigTable[i] = 0; + } + count = 0; + for(i=0; i<sdim; i++){ + for(j=0; j<bdim; j++){ + if(psDomain[i] == pbDomain[j]){ + mask[count] = j; + count++; + break; + } + } + } + + bCumprod[0] = 1; + for(i=0; i<bdim-1; i++){ + bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; + } + sCumprod[0] = 1; + for(i=0; i<sdim-1; i++){ + sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; + } + + for(i=0; i<NZB; i++){ + value = bpr[i]; + bindex = bir[i]; + ind_subv(bindex, bCumprod, bdim, bsubv); + for(j=0; j<sdim; j++){ + ssubv[j] = bsubv[mask[j]]; + } + sindex = subv_ind(sdim, sCumprod, ssubv); + result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); + if(result){ + position = result - sir; + value *= spr[position]; + bigTable[nzCounts] = value; + index[nzCounts] = bindex; + nzCounts++; + } + } + + pTemp = mxGetField(bigPot, 0, "T"); + if(pTemp)mxDestroyArray(pTemp); + pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB); + mxSetField(bigPot, 0, "T", pTemp); + + free(bigTable); + free(index); + free(mask); + free(bCumprod); + free(sCumprod); + free(bsubv); + free(ssubv); +} + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int i, j, c, loop, nNodes, nCliques, ndomain, dims[2]; + double *pClqs, *pr, *pt, *pSize; + mxArray *pTemp, *pTemp1, *pStruct, *pCliques, *pBigpot, *pSmallpot; + const char *field_names[] = {"domain", "T", "sizes"}; + + nNodes = mxGetNumberOfElements(prhs[1]); + pCliques = mxGetField(prhs[0], 0, "cliques"); + nCliques = mxGetNumberOfElements(pCliques); + pTemp = mxGetField(prhs[0], 0, "eff_node_sizes"); + pSize = mxGetPr(pTemp); + + plhs[0] = mxCreateCellArray(1, &nCliques); + for(i=0; i<nCliques; i++){ + pStruct = mxCreateStructMatrix(1, 1, 3, field_names); + mxSetCell(plhs[0], i, pStruct); + pTemp = mxGetCell(pCliques, i); + ndomain = mxGetNumberOfElements(pTemp); + pt = mxGetPr(pTemp); + pTemp1 = mxDuplicateArray(pTemp); + mxSetField(pStruct, 0, "domain", pTemp1); + + pTemp = mxCreateDoubleMatrix(1, ndomain, mxREAL); + mxSetField(pStruct, 0, "sizes", pTemp); + pr = mxGetPr(pTemp); + for(j=0; j<ndomain; j++){ + pr[j] = pSize[(int)pt[j]-1]; + } + } + + pClqs = mxGetPr(prhs[1]); + for(loop=0; loop<nNodes; loop++){ + c = (int)pClqs[loop] - 1; + pSmallpot = mxGetCell(prhs[2], loop); + pTemp = mxGetField(pSmallpot, 0, "T"); + pBigpot = mxGetCell(plhs[0], c); + pTemp1 = mxGetField(pBigpot, 0, "T"); + if(pTemp1){ + if(mxIsSparse(pTemp)) + multiply_spPot_by_spPot(pBigpot, pSmallpot); + else multiply_spPot_by_fuPot(pBigpot, pSmallpot); + } + else{ + if(mxIsSparse(pTemp)) + multiply_null_by_spPot(pBigpot, pSmallpot); + else multiply_null_by_fuPot(pBigpot, pSmallpot); + } + } + + dims[0] = nCliques; + dims[1] = nCliques; + plhs[1] = mxCreateCellArray(2, dims); +} + + diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.m new file mode 100644 index 00000000..857e6266 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/old/init_pot1.m @@ -0,0 +1,20 @@ +function [clpot, seppot] = init_pot(engine, clqs, pots, pot_type, onodes, ndx) +% INIT_POT Initialise potentials with evidence (jtree_inf) +% function [clpot, seppot] = init_pot(engine, clqs, pots, pot_type, onodes) + +cliques = engine.cliques; +bnet = bnet_from_engine(engine); +% Set the clique potentials to all 1s +C = length(cliques); +clpot = cell(1,C); +for i=1:C + clpot{i} = mk_initial_pot(pot_type, cliques{i}, bnet.node_sizes(:), bnet.cnodes(:), onodes); +end + +% Multiply on specified potentials +for i=1:length(clqs) + c = clqs(i); + clpot{c} = multiply_by_pot(clpot{c}, pots{i}); +end + +seppot = cell(C,C); % implicitely initialized to 1 diff --git a/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/set_fields.m b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/set_fields.m new file mode 100644 index 00000000..e75cfa45 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/static/@jtree_sparse_inf_engine/set_fields.m @@ -0,0 +1,13 @@ +function engine = set_fields(engine, varargin) +% SET_FIELDS Set the fields for a generic engine +% engine = set_fields(engine, name/value pairs) +% +% e.g., engine = set_fields(engine, 'maximize', 1) + +args = varargin; +nargs = length(args); +for i=1:2:nargs + switch args{i}, + case 'maximize', engine.maximize = args{i+1}; + end +end |