Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/inference/static/@jtree_sparse_inf_engine/collect_evidence.c @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 /* C mex for collect_evidence.c in @jtree_sparse_inf_engine directory */
wolffd@0 2 /* File enter_evidence.m in directory @jtree_sparse_inf_engine call it*/
wolffd@0 3
wolffd@0 4 /******************************************/
wolffd@0 5 /* collect_evidence has 3 input & 2 output*/
wolffd@0 6 /* engine */
wolffd@0 7 /* clpot */
wolffd@0 8 /* seppot */
wolffd@0 9 /* */
wolffd@0 10 /* clpot */
wolffd@0 11 /* seppot */
wolffd@0 12 /******************************************/
wolffd@0 13
wolffd@0 14 #include <math.h>
wolffd@0 15 #include <stdlib.h>
wolffd@0 16 #include "mex.h"
wolffd@0 17
wolffd@0 18 int compare(const void* src1, const void* src2){
wolffd@0 19 int i1 = *(int*)src1 ;
wolffd@0 20 int i2 = *(int*)src2 ;
wolffd@0 21 return i1-i2 ;
wolffd@0 22 }
wolffd@0 23
wolffd@0 24 void ind_subv(int index, const int *cumprod, int n, int *bsubv){
wolffd@0 25 int i;
wolffd@0 26
wolffd@0 27 for (i = n-1; i >= 0; i--) {
wolffd@0 28 bsubv[i] = ((int)floor(index / cumprod[i]));
wolffd@0 29 index = index % cumprod[i];
wolffd@0 30 }
wolffd@0 31 }
wolffd@0 32
wolffd@0 33 int subv_ind(const int n, const int *cumprod, const int *subv){
wolffd@0 34 int i, index=0;
wolffd@0 35
wolffd@0 36 for(i=0; i<n; i++){
wolffd@0 37 index += subv[i] * cumprod[i];
wolffd@0 38 }
wolffd@0 39 return index;
wolffd@0 40 }
wolffd@0 41
wolffd@0 42 void compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){
wolffd@0 43 int i, j;
wolffd@0 44 int *eff_cumprod, *subv, *diffsize, *diff_cumprod;
wolffd@0 45
wolffd@0 46 subv = malloc(diffdim * sizeof(int));
wolffd@0 47 eff_cumprod = malloc(diffdim * sizeof(int));
wolffd@0 48 diffsize = malloc(diffdim * sizeof(int));
wolffd@0 49 diff_cumprod = malloc(diffdim * sizeof(int));
wolffd@0 50 for(i=0; i<diffdim; i++){
wolffd@0 51 eff_cumprod[i] = bCumprod[dmask[i]];
wolffd@0 52 diffsize[i] = (int)pbSize[dmask[i]];
wolffd@0 53 }
wolffd@0 54 diff_cumprod[0] = 1;
wolffd@0 55 for(i=0; i<diffdim-1; i++){
wolffd@0 56 diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i];
wolffd@0 57 }
wolffd@0 58 for(i=0; i<ND; i++){
wolffd@0 59 ind_subv(i, diff_cumprod, diffdim, subv);
wolffd@0 60 weight[i] = 0;
wolffd@0 61 for(j=0; j<diffdim; j++){
wolffd@0 62 weight[i] += eff_cumprod[j] * subv[j];
wolffd@0 63 }
wolffd@0 64 }
wolffd@0 65 free(eff_cumprod);
wolffd@0 66 free(subv);
wolffd@0 67 free(diffsize);
wolffd@0 68 free(diff_cumprod);
wolffd@0 69 }
wolffd@0 70
wolffd@0 71 void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){
wolffd@0 72 double *ptr;
wolffd@0 73 void *newptr;
wolffd@0 74 int *ir, *jc;
wolffd@0 75 int nbytes;
wolffd@0 76
wolffd@0 77 if(new_nzmax == old_nzmax) return;
wolffd@0 78 nbytes = new_nzmax * sizeof(*ptr);
wolffd@0 79 ptr = mxGetPr(spArray);
wolffd@0 80 newptr = mxRealloc(ptr, nbytes);
wolffd@0 81 mxSetPr(spArray, newptr);
wolffd@0 82 nbytes = new_nzmax * sizeof(*ir);
wolffd@0 83 ir = mxGetIr(spArray);
wolffd@0 84 newptr = mxRealloc(ir, nbytes);
wolffd@0 85 mxSetIr(spArray, newptr);
wolffd@0 86 jc = mxGetJc(spArray);
wolffd@0 87 jc[0] = 0;
wolffd@0 88 jc[1] = new_nzmax;
wolffd@0 89 mxSetNzmax(spArray, new_nzmax);
wolffd@0 90 }
wolffd@0 91
wolffd@0 92 mxArray* convert_ill_table_to_sparse(const double *Table, const int *sequence, const int nzCounts, const int N){
wolffd@0 93 mxArray *spTable;
wolffd@0 94 int i, temp, *irs, *jcs, count=0;
wolffd@0 95 double *sr;
wolffd@0 96
wolffd@0 97 spTable = mxCreateSparse(N, 1, nzCounts, mxREAL);
wolffd@0 98 sr = mxGetPr(spTable);
wolffd@0 99 irs = mxGetIr(spTable);
wolffd@0 100 jcs = mxGetJc(spTable);
wolffd@0 101
wolffd@0 102 jcs[0] = 0;
wolffd@0 103 jcs[1] = nzCounts;
wolffd@0 104
wolffd@0 105 for(i=0; i<nzCounts; i++){
wolffd@0 106 irs[i] = sequence[count];
wolffd@0 107 count++;
wolffd@0 108 temp = sequence[count];
wolffd@0 109 sr[i] = Table[temp];
wolffd@0 110 count++;
wolffd@0 111 }
wolffd@0 112 return spTable;
wolffd@0 113 }
wolffd@0 114
wolffd@0 115 void multiply_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){
wolffd@0 116 int i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex, nzCounts=0;
wolffd@0 117 int *samemask, *diffmask, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *sequence, *weight;
wolffd@0 118 double *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *spr, *bpr;
wolffd@0 119 mxArray *pTemp, *pTemp1;
wolffd@0 120
wolffd@0 121 pTemp = mxGetField(bigPot, 0, "domain");
wolffd@0 122 pbDomain = mxGetPr(pTemp);
wolffd@0 123 bdim = mxGetNumberOfElements(pTemp);
wolffd@0 124 pTemp = mxGetField(smallPot, 0, "domain");
wolffd@0 125 psDomain = mxGetPr(pTemp);
wolffd@0 126 sdim = mxGetNumberOfElements(pTemp);
wolffd@0 127
wolffd@0 128 pTemp = mxGetField(bigPot, 0, "sizes");
wolffd@0 129 pbSize = mxGetPr(pTemp);
wolffd@0 130 pTemp = mxGetField(smallPot, 0, "sizes");
wolffd@0 131 psSize = mxGetPr(pTemp);
wolffd@0 132
wolffd@0 133 pTemp = mxGetField(smallPot, 0, "T");
wolffd@0 134 spr = mxGetPr(pTemp);
wolffd@0 135 sir = mxGetIr(pTemp);
wolffd@0 136 sjc = mxGetJc(pTemp);
wolffd@0 137 NZS = sjc[1];
wolffd@0 138
wolffd@0 139 NB = 1;
wolffd@0 140 for(i=0; i<bdim; i++){
wolffd@0 141 NB *= (int)pbSize[i];
wolffd@0 142 }
wolffd@0 143
wolffd@0 144 if(sdim == 0){
wolffd@0 145 pTemp = mxCreateSparse(NB, 1, NB, mxREAL);
wolffd@0 146 mxSetField(bigPot, 0, "T", pTemp);
wolffd@0 147 bpr = mxGetPr(pTemp);
wolffd@0 148 sir = mxGetIr(pTemp);
wolffd@0 149 sjc = mxGetJc(pTemp);
wolffd@0 150 sjc[0] = 0;
wolffd@0 151 sjc[1] = NB;
wolffd@0 152 for(i=0; i<NB; i++){
wolffd@0 153 bpr[i] = *spr;
wolffd@0 154 sir[i] = i;
wolffd@0 155 }
wolffd@0 156 return;
wolffd@0 157 }
wolffd@0 158
wolffd@0 159 NS = 1;
wolffd@0 160 for(i=0; i<sdim; i++){
wolffd@0 161 NS *= (int)psSize[i];
wolffd@0 162 }
wolffd@0 163 ND = NB / NS;
wolffd@0 164
wolffd@0 165 if(ND == 1){
wolffd@0 166 pTemp1 = mxGetField(smallPot, 0, "T");
wolffd@0 167 pTemp = mxDuplicateArray(pTemp1);
wolffd@0 168 mxSetField(bigPot, 0, "T", pTemp);
wolffd@0 169 return;
wolffd@0 170 }
wolffd@0 171
wolffd@0 172
wolffd@0 173 NZB = ND * NZS;
wolffd@0 174
wolffd@0 175 diffdim = bdim - sdim;
wolffd@0 176 sequence = malloc(NZB * 2 * sizeof(int));
wolffd@0 177 bigTable = malloc(NZB * sizeof(double));
wolffd@0 178 samemask = malloc(sdim * sizeof(int));
wolffd@0 179 diffmask = malloc(diffdim * sizeof(int));
wolffd@0 180 bCumprod = malloc(bdim * sizeof(int));
wolffd@0 181 sCumprod = malloc(sdim * sizeof(int));
wolffd@0 182 weight = malloc(ND * sizeof(int));
wolffd@0 183 ssubv = malloc(sdim * sizeof(int));
wolffd@0 184
wolffd@0 185 count = 0;
wolffd@0 186 count1 = 0;
wolffd@0 187 for(i=0; i<bdim; i++){
wolffd@0 188 match = 0;
wolffd@0 189 for(j=0; j<sdim; j++){
wolffd@0 190 if(pbDomain[i] == psDomain[j]){
wolffd@0 191 samemask[count] = i;
wolffd@0 192 match = 1;
wolffd@0 193 count++;
wolffd@0 194 break;
wolffd@0 195 }
wolffd@0 196 }
wolffd@0 197 if(match == 0){
wolffd@0 198 diffmask[count1] = i;
wolffd@0 199 count1++;
wolffd@0 200 }
wolffd@0 201 }
wolffd@0 202
wolffd@0 203 bCumprod[0] = 1;
wolffd@0 204 for(i=0; i<bdim-1; i++){
wolffd@0 205 bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
wolffd@0 206 }
wolffd@0 207 sCumprod[0] = 1;
wolffd@0 208 for(i=0; i<sdim-1; i++){
wolffd@0 209 sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
wolffd@0 210 }
wolffd@0 211
wolffd@0 212 count = 0;
wolffd@0 213 compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim);
wolffd@0 214 for(i=0; i<NZS; i++){
wolffd@0 215 sindex = sir[i];
wolffd@0 216 ind_subv(sindex, sCumprod, sdim, ssubv);
wolffd@0 217 temp = 0;
wolffd@0 218 for(j=0; j<sdim; j++){
wolffd@0 219 temp += ssubv[j] * bCumprod[samemask[j]];
wolffd@0 220 }
wolffd@0 221 for(j=0; j<ND; j++){
wolffd@0 222 bindex = weight[j] + temp;
wolffd@0 223 bigTable[nzCounts] = spr[i];
wolffd@0 224 sequence[count] = bindex;
wolffd@0 225 count++;
wolffd@0 226 sequence[count] = nzCounts;
wolffd@0 227 nzCounts++;
wolffd@0 228 count++;
wolffd@0 229 }
wolffd@0 230 }
wolffd@0 231
wolffd@0 232 pTemp = mxGetField(bigPot, 0, "T");
wolffd@0 233 if(pTemp)mxDestroyArray(pTemp);
wolffd@0 234 qsort(sequence, nzCounts, sizeof(int) * 2, compare);
wolffd@0 235 pTemp = convert_ill_table_to_sparse(bigTable, sequence, nzCounts, NB);
wolffd@0 236 mxSetField(bigPot, 0, "T", pTemp);
wolffd@0 237
wolffd@0 238 free(sequence);
wolffd@0 239 free(bigTable);
wolffd@0 240 free(samemask);
wolffd@0 241 free(diffmask);
wolffd@0 242 free(bCumprod);
wolffd@0 243 free(sCumprod);
wolffd@0 244 free(weight);
wolffd@0 245 free(ssubv);
wolffd@0 246 }
wolffd@0 247
wolffd@0 248 void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){
wolffd@0 249 int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0;
wolffd@0 250 int *mask, *result, *bir, *sir, *rir, *bjc, *sjc, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
wolffd@0 251 double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr;
wolffd@0 252 mxArray *pTemp, *pTemp1;
wolffd@0 253
wolffd@0 254 pTemp = mxGetField(bigPot, 0, "domain");
wolffd@0 255 pbDomain = mxGetPr(pTemp);
wolffd@0 256 bdim = mxGetNumberOfElements(pTemp);
wolffd@0 257 pTemp = mxGetField(smallPot, 0, "domain");
wolffd@0 258 psDomain = mxGetPr(pTemp);
wolffd@0 259 sdim = mxGetNumberOfElements(pTemp);
wolffd@0 260
wolffd@0 261 pTemp = mxGetField(bigPot, 0, "sizes");
wolffd@0 262 pbSize = mxGetPr(pTemp);
wolffd@0 263 pTemp = mxGetField(smallPot, 0, "sizes");
wolffd@0 264 psSize = mxGetPr(pTemp);
wolffd@0 265
wolffd@0 266 NB = 1;
wolffd@0 267 for(i=0; i<bdim; i++){
wolffd@0 268 NB *= (int)pbSize[i];
wolffd@0 269 }
wolffd@0 270
wolffd@0 271 pTemp = mxGetField(bigPot, 0, "T");
wolffd@0 272 bpr = mxGetPr(pTemp);
wolffd@0 273 bir = mxGetIr(pTemp);
wolffd@0 274 bjc = mxGetJc(pTemp);
wolffd@0 275 NZB = bjc[1];
wolffd@0 276
wolffd@0 277 pTemp = mxGetField(smallPot, 0, "T");
wolffd@0 278 spr = mxGetPr(pTemp);
wolffd@0 279 sir = mxGetIr(pTemp);
wolffd@0 280 sjc = mxGetJc(pTemp);
wolffd@0 281 NZS = sjc[1];
wolffd@0 282
wolffd@0 283 if(sdim == 0){
wolffd@0 284 for(i=0; i<NZB; i++){
wolffd@0 285 bpr[i] *= *spr;
wolffd@0 286 }
wolffd@0 287 return;
wolffd@0 288 }
wolffd@0 289
wolffd@0 290 pTemp1 = mxCreateSparse(NB, 1, NZB, mxREAL);
wolffd@0 291 rpr = mxGetPr(pTemp1);
wolffd@0 292 rir = mxGetIr(pTemp1);
wolffd@0 293 rjc = mxGetJc(pTemp1);
wolffd@0 294 rjc[0] = 0;
wolffd@0 295 rjc[1] = NZB;
wolffd@0 296
wolffd@0 297 mask = malloc(sdim * sizeof(int));
wolffd@0 298 bCumprod = malloc(bdim * sizeof(int));
wolffd@0 299 sCumprod = malloc(sdim * sizeof(int));
wolffd@0 300 bsubv = malloc(bdim * sizeof(int));
wolffd@0 301 ssubv = malloc(sdim * sizeof(int));
wolffd@0 302
wolffd@0 303 count = 0;
wolffd@0 304 for(i=0; i<sdim; i++){
wolffd@0 305 for(j=0; j<bdim; j++){
wolffd@0 306 if(psDomain[i] == pbDomain[j]){
wolffd@0 307 mask[count] = j;
wolffd@0 308 count++;
wolffd@0 309 break;
wolffd@0 310 }
wolffd@0 311 }
wolffd@0 312 }
wolffd@0 313
wolffd@0 314 bCumprod[0] = 1;
wolffd@0 315 for(i=0; i<bdim-1; i++){
wolffd@0 316 bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
wolffd@0 317 }
wolffd@0 318 sCumprod[0] = 1;
wolffd@0 319 for(i=0; i<sdim-1; i++){
wolffd@0 320 sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
wolffd@0 321 }
wolffd@0 322
wolffd@0 323 for(i=0; i<NZB; i++){
wolffd@0 324 bindex = bir[i];
wolffd@0 325 ind_subv(bindex, bCumprod, bdim, bsubv);
wolffd@0 326 for(j=0; j<sdim; j++){
wolffd@0 327 ssubv[j] = bsubv[mask[j]];
wolffd@0 328 }
wolffd@0 329 sindex = subv_ind(sdim, sCumprod, ssubv);
wolffd@0 330 result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare);
wolffd@0 331 if(result){
wolffd@0 332 position = result - sir;
wolffd@0 333 rpr[nzCounts] = bpr[i] * spr[position];
wolffd@0 334 rir[nzCounts] = bindex;
wolffd@0 335 nzCounts++;
wolffd@0 336 }
wolffd@0 337 }
wolffd@0 338
wolffd@0 339 pTemp = mxGetField(bigPot, 0, "T");
wolffd@0 340 if(pTemp)mxDestroyArray(pTemp);
wolffd@0 341 reset_nzmax(pTemp1, NZB, nzCounts);
wolffd@0 342 mxSetField(bigPot, 0, "T", pTemp1);
wolffd@0 343
wolffd@0 344 free(mask);
wolffd@0 345 free(bCumprod);
wolffd@0 346 free(sCumprod);
wolffd@0 347 free(bsubv);
wolffd@0 348 free(ssubv);
wolffd@0 349 }
wolffd@0 350
wolffd@0 351 mxArray* marginal_null_to_spPot(const mxArray *bigPot, const mxArray *sDomain, const int maximize){
wolffd@0 352 int i, j, count, bdim, sdim, NB, NS, ND;
wolffd@0 353 int *mask, *sir, *sjc;
wolffd@0 354 double *pbDomain, *psDomain, *pbSize, *psSize, *spr;
wolffd@0 355 mxArray *pTemp, *smallPot;
wolffd@0 356 const char *field_names[] = {"domain", "T", "sizes"};
wolffd@0 357
wolffd@0 358 pTemp = mxGetField(bigPot, 0, "domain");
wolffd@0 359 pbDomain = mxGetPr(pTemp);
wolffd@0 360 bdim = mxGetNumberOfElements(pTemp);
wolffd@0 361 psDomain = mxGetPr(sDomain);
wolffd@0 362 sdim = mxGetNumberOfElements(sDomain);
wolffd@0 363 pTemp = mxGetField(bigPot, 0, "sizes");
wolffd@0 364 pbSize = mxGetPr(pTemp);
wolffd@0 365
wolffd@0 366 smallPot = mxCreateStructMatrix(1, 1, 3, field_names);
wolffd@0 367 pTemp = mxDuplicateArray(sDomain);
wolffd@0 368 mxSetField(smallPot, 0, "domain", pTemp);
wolffd@0 369
wolffd@0 370 NB = 1;
wolffd@0 371 for(i=0; i<bdim; i++){
wolffd@0 372 NB *= (int)pbSize[i];
wolffd@0 373 }
wolffd@0 374
wolffd@0 375 if(sdim == 0){
wolffd@0 376 pTemp = mxCreateSparse(1, 1, 1, mxREAL);
wolffd@0 377 mxSetField(smallPot, 0, "T", pTemp);
wolffd@0 378 spr = mxGetPr(pTemp);
wolffd@0 379 sir = mxGetIr(pTemp);
wolffd@0 380 sjc = mxGetJc(pTemp);
wolffd@0 381 *spr = 0;
wolffd@0 382 *sir = 0;
wolffd@0 383 sjc[0] = 0;
wolffd@0 384 sjc[1] = 1;
wolffd@0 385 if(maximize) *spr = 1;
wolffd@0 386 else *spr = NB;
wolffd@0 387
wolffd@0 388 pTemp = mxCreateDoubleMatrix(1, 1, mxREAL);
wolffd@0 389 *mxGetPr(pTemp) = 1;
wolffd@0 390 mxSetField(smallPot, 0, "sizes", pTemp);
wolffd@0 391 return smallPot;
wolffd@0 392 }
wolffd@0 393
wolffd@0 394 mask = malloc(sdim * sizeof(int));
wolffd@0 395 count = 0;
wolffd@0 396 for(i=0; i<sdim; i++){
wolffd@0 397 for(j=0; j<bdim; j++){
wolffd@0 398 if(psDomain[i] == pbDomain[j]){
wolffd@0 399 mask[count] = j;
wolffd@0 400 count++;
wolffd@0 401 break;
wolffd@0 402 }
wolffd@0 403 }
wolffd@0 404 }
wolffd@0 405 pTemp = mxCreateDoubleMatrix(1, count, mxREAL);
wolffd@0 406 psSize = mxGetPr(pTemp);
wolffd@0 407 NS = 1;
wolffd@0 408 for(i=0; i<count; i++){
wolffd@0 409 psSize[i] = pbSize[mask[i]];
wolffd@0 410 NS *= (int)psSize[i];
wolffd@0 411 }
wolffd@0 412 mxSetField(smallPot, 0, "sizes", pTemp);
wolffd@0 413
wolffd@0 414 ND = NB / NS;
wolffd@0 415
wolffd@0 416 pTemp = mxCreateSparse(NS, 1, NS, mxREAL);
wolffd@0 417 mxSetField(smallPot, 0, "T", pTemp);
wolffd@0 418 spr = mxGetPr(pTemp);
wolffd@0 419 sir = mxGetIr(pTemp);
wolffd@0 420 sjc = mxGetJc(pTemp);
wolffd@0 421 if(maximize){
wolffd@0 422 for(i=0; i<NS; i++){
wolffd@0 423 spr[i] = 1;
wolffd@0 424 sir[i] = i;
wolffd@0 425 }
wolffd@0 426 }
wolffd@0 427 else{
wolffd@0 428 for(i=0; i<NS; i++){
wolffd@0 429 spr[i] = ND;
wolffd@0 430 sir[i] = i;
wolffd@0 431 }
wolffd@0 432 }
wolffd@0 433 sjc[0] = 0;
wolffd@0 434 sjc[1] = NS;
wolffd@0 435
wolffd@0 436 free(mask);
wolffd@0 437 return smallPot;
wolffd@0 438 }
wolffd@0 439
wolffd@0 440 mxArray* marginal_spPot_to_spPot(const mxArray *bigPot, const mxArray *sDomain, const int maximize){
wolffd@0 441 int i, j, count, bdim, sdim, NB, NS, NZB, position, bindex, sindex, nzCounts=0;
wolffd@0 442 int *mask, *sequence, *result, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
wolffd@0 443 double *sTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr;
wolffd@0 444 mxArray *pTemp, *smallPot;
wolffd@0 445 const char *field_names[] = {"domain", "T", "sizes"};
wolffd@0 446
wolffd@0 447 pTemp = mxGetField(bigPot, 0, "domain");
wolffd@0 448 pbDomain = mxGetPr(pTemp);
wolffd@0 449 bdim = mxGetNumberOfElements(pTemp);
wolffd@0 450 psDomain = mxGetPr(sDomain);
wolffd@0 451 sdim = mxGetNumberOfElements(sDomain);
wolffd@0 452 pTemp = mxGetField(bigPot, 0, "sizes");
wolffd@0 453 pbSize = mxGetPr(pTemp);
wolffd@0 454
wolffd@0 455 pTemp = mxGetField(bigPot, 0, "T");
wolffd@0 456 bpr = mxGetPr(pTemp);
wolffd@0 457 bir = mxGetIr(pTemp);
wolffd@0 458 bjc = mxGetJc(pTemp);
wolffd@0 459 NZB = bjc[1];
wolffd@0 460
wolffd@0 461 smallPot = mxCreateStructMatrix(1, 1, 3, field_names);
wolffd@0 462 pTemp = mxDuplicateArray(sDomain);
wolffd@0 463 mxSetField(smallPot, 0, "domain", pTemp);
wolffd@0 464
wolffd@0 465 if(sdim == 0){
wolffd@0 466 pTemp = mxCreateSparse(1, 1, 1, mxREAL);
wolffd@0 467 mxSetField(smallPot, 0, "T", pTemp);
wolffd@0 468 spr = mxGetPr(pTemp);
wolffd@0 469 bir = mxGetIr(pTemp);
wolffd@0 470 bjc = mxGetJc(pTemp);
wolffd@0 471 *spr = 0;
wolffd@0 472 *bir = 0;
wolffd@0 473 bjc[0] = 0;
wolffd@0 474 bjc[1] = 1;
wolffd@0 475 if(maximize){
wolffd@0 476 for(i=0; i<NZB; i++){
wolffd@0 477 *spr = (*spr < bpr[i])? bpr[i] : *spr;
wolffd@0 478 }
wolffd@0 479 }
wolffd@0 480 else{
wolffd@0 481 for(i=0; i<NZB; i++){
wolffd@0 482 *spr += bpr[i];
wolffd@0 483 }
wolffd@0 484 }
wolffd@0 485
wolffd@0 486 pTemp = mxCreateDoubleMatrix(1, 1, mxREAL);
wolffd@0 487 *mxGetPr(pTemp) = 1;
wolffd@0 488 mxSetField(smallPot, 0, "sizes", pTemp);
wolffd@0 489 return smallPot;
wolffd@0 490 }
wolffd@0 491
wolffd@0 492 NB = 1;
wolffd@0 493 for(i=0; i<bdim; i++){
wolffd@0 494 NB *= (int)pbSize[i];
wolffd@0 495 }
wolffd@0 496
wolffd@0 497 mask = malloc(sdim * sizeof(int));
wolffd@0 498 count = 0;
wolffd@0 499 for(i=0; i<sdim; i++){
wolffd@0 500 for(j=0; j<bdim; j++){
wolffd@0 501 if(psDomain[i] == pbDomain[j]){
wolffd@0 502 mask[count] = j;
wolffd@0 503 count++;
wolffd@0 504 break;
wolffd@0 505 }
wolffd@0 506 }
wolffd@0 507 }
wolffd@0 508 pTemp = mxCreateDoubleMatrix(1, count, mxREAL);
wolffd@0 509 psSize = mxGetPr(pTemp);
wolffd@0 510 NS = 1;
wolffd@0 511 for(i=0; i<count; i++){
wolffd@0 512 psSize[i] = pbSize[mask[i]];
wolffd@0 513 NS *= (int)psSize[i];
wolffd@0 514 }
wolffd@0 515 mxSetField(smallPot, 0, "sizes", pTemp);
wolffd@0 516
wolffd@0 517
wolffd@0 518 sTable = malloc(NZB * sizeof(double));
wolffd@0 519 sequence = malloc(NZB * 2 * sizeof(double));
wolffd@0 520 bCumprod = malloc(bdim * sizeof(int));
wolffd@0 521 sCumprod = malloc(sdim * sizeof(int));
wolffd@0 522 bsubv = malloc(bdim * sizeof(int));
wolffd@0 523 ssubv = malloc(sdim * sizeof(int));
wolffd@0 524
wolffd@0 525 for(i=0; i<NZB; i++)sTable[i] = 0;
wolffd@0 526
wolffd@0 527 bCumprod[0] = 1;
wolffd@0 528 for(i=0; i<bdim-1; i++){
wolffd@0 529 bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
wolffd@0 530 }
wolffd@0 531 sCumprod[0] = 1;
wolffd@0 532 for(i=0; i<sdim-1; i++){
wolffd@0 533 sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
wolffd@0 534 }
wolffd@0 535
wolffd@0 536 count = 0;
wolffd@0 537 for(i=0; i<NZB; i++){
wolffd@0 538 bindex = bir[i];
wolffd@0 539 ind_subv(bindex, bCumprod, bdim, bsubv);
wolffd@0 540 for(j=0; j<sdim; j++){
wolffd@0 541 ssubv[j] = bsubv[mask[j]];
wolffd@0 542 }
wolffd@0 543 sindex = subv_ind(sdim, sCumprod, ssubv);
wolffd@0 544 result = (int *) bsearch(&sindex, sequence, nzCounts, sizeof(int)*2, compare);
wolffd@0 545 if(result){
wolffd@0 546 position = (result - sequence) / 2;
wolffd@0 547 if(maximize)
wolffd@0 548 sTable[position] = (sTable[position] < bpr[i]) ? bpr[i] : sTable[position];
wolffd@0 549 else sTable[position] += bpr[i];
wolffd@0 550 }
wolffd@0 551 else {
wolffd@0 552 if(maximize)
wolffd@0 553 sTable[nzCounts] = (sTable[nzCounts] < bpr[i]) ? bpr[i] : sTable[nzCounts];
wolffd@0 554 else sTable[nzCounts] += bpr[i];
wolffd@0 555 sequence[count] = sindex;
wolffd@0 556 count++;
wolffd@0 557 sequence[count] = nzCounts;
wolffd@0 558 nzCounts++;
wolffd@0 559 count++;
wolffd@0 560 }
wolffd@0 561 }
wolffd@0 562
wolffd@0 563 qsort(sequence, nzCounts, sizeof(int) * 2, compare);
wolffd@0 564 pTemp = convert_ill_table_to_sparse(sTable, sequence, nzCounts, NS);
wolffd@0 565 mxSetField(smallPot, 0, "T", pTemp);
wolffd@0 566
wolffd@0 567 free(sTable);
wolffd@0 568 free(sequence);
wolffd@0 569 free(mask);
wolffd@0 570 free(bCumprod);
wolffd@0 571 free(sCumprod);
wolffd@0 572 free(bsubv);
wolffd@0 573 free(ssubv);
wolffd@0 574
wolffd@0 575 return smallPot;
wolffd@0 576 }
wolffd@0 577
wolffd@0 578
wolffd@0 579 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
wolffd@0 580 int i, n, p, np, pn, loop, loops, nCliques, temp, maximize;
wolffd@0 581 int *collect_order;
wolffd@0 582 double *pr, *pr1;
wolffd@0 583 mxArray *pTemp, *pTemp1, *pPostP, *pClpot, *pSeppot, *pSeparator;
wolffd@0 584
wolffd@0 585 pTemp = mxGetField(prhs[0], 0, "cliques");
wolffd@0 586 nCliques = mxGetNumberOfElements(pTemp);
wolffd@0 587 loops = nCliques - 1;
wolffd@0 588 pTemp = mxGetField(prhs[0], 0, "maximize");
wolffd@0 589 maximize = (int)mxGetScalar(pTemp);
wolffd@0 590 pSeparator = mxGetField(prhs[0], 0, "separator");
wolffd@0 591
wolffd@0 592 collect_order = malloc(2 * loops * sizeof(int));
wolffd@0 593
wolffd@0 594 pTemp = mxGetField(prhs[0], 0, "postorder");
wolffd@0 595 pr = mxGetPr(pTemp);
wolffd@0 596 pPostP = mxGetField(prhs[0], 0, "postorder_parents");
wolffd@0 597 for(i=0; i<loops; i++){
wolffd@0 598 temp = (int)pr[i] - 1;
wolffd@0 599 pTemp = mxGetCell(pPostP, temp);
wolffd@0 600 pr1 = mxGetPr(pTemp);
wolffd@0 601 collect_order[i] = (int)pr1[0] - 1;
wolffd@0 602 collect_order[i+loops] = temp;
wolffd@0 603 }
wolffd@0 604
wolffd@0 605 plhs[0] = mxDuplicateArray(prhs[1]);
wolffd@0 606 plhs[1] = mxDuplicateArray(prhs[2]);
wolffd@0 607
wolffd@0 608 for(loop=0; loop<loops; loop++){
wolffd@0 609 p = collect_order[loop];
wolffd@0 610 n = collect_order[loop+loops];
wolffd@0 611 np = p * nCliques + n;
wolffd@0 612 pn = n * nCliques + p;
wolffd@0 613 pClpot = mxGetCell(plhs[0], n);
wolffd@0 614 pTemp1 = mxGetField(pClpot, 0, "T");
wolffd@0 615 pTemp = mxGetCell(pSeparator, pn);
wolffd@0 616 if(pTemp1)
wolffd@0 617 pSeppot = marginal_spPot_to_spPot(pClpot, pTemp, maximize);
wolffd@0 618 else pSeppot = marginal_null_to_spPot(pClpot, pTemp, maximize);
wolffd@0 619 mxSetCell(plhs[1], pn, pSeppot);
wolffd@0 620
wolffd@0 621 pClpot = mxGetCell(plhs[0], p);
wolffd@0 622 pTemp1 = mxGetField(pClpot, 0, "T");
wolffd@0 623 if(pTemp1)
wolffd@0 624 multiply_spPot_by_spPot(pClpot, pSeppot);
wolffd@0 625 else multiply_null_by_spPot(pClpot, pSeppot);
wolffd@0 626 }
wolffd@0 627 free(collect_order);
wolffd@0 628 }
wolffd@0 629
wolffd@0 630
wolffd@0 631
wolffd@0 632
wolffd@0 633
wolffd@0 634