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
|