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
|