Mercurial > hg > camir-aes2014

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