Mercurial > hg > camir-aes2014

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