Mercurial > hg > camir-aes2014
diff toolboxes/FullBNT-1.0.7/bnt/inference/static/@jtree_sparse_inf_engine/collect_evidence.c @ 0:e9a9cd732c1e tip
first hg version after svn
author | wolffd |
---|---|
date | Tue, 10 Feb 2015 15:05:51 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/FullBNT-1.0.7/bnt/inference/static/@jtree_sparse_inf_engine/collect_evidence.c Tue Feb 10 15:05:51 2015 +0000 @@ -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); +} + + + + + +