Mercurial > hg > camir-aes2014
diff toolboxes/FullBNT-1.0.7/bnt/inference/static/@jtree_sparse_inf_engine/init_pot.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/init_pot.c Tue Feb 10 15:05:51 2015 +0000 @@ -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); +} + +