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
|