Mercurial > hg > camir-aes2014
view 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 source
/* 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); }