Mercurial > hg > camir-aes2014

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