comparison toolboxes/FullBNT-1.0.7/bnt/potentials/Tables/marg_table.c @ 0:e9a9cd732c1e tip

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 /* marg_table.c ../potential/tables */
2
3
4 /******************************************/
5 /* 5 input & 1 output */
6 /* Big table */
7 /* Big domain */
8 /* Big sizes */
9 /* onto */
10 /* maximize, if missed, maximize=0 */
11 /* */
12 /* small table */
13 /******************************************/
14
15 #include "mex.h"
16
17 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
18 int i, j, count, NB, NS, siz_b, siz_s, ndim, temp, maximize;
19 int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *ssize;
20 double *pb, *ps, *bp, *sp, *pbd;
21
22
23 siz_b = mxGetNumberOfElements(prhs[1]);
24 siz_s = mxGetNumberOfElements(prhs[3]);
25 pb = mxGetPr(prhs[1]);
26 ps = mxGetPr(prhs[3]);
27
28 NB = mxGetNumberOfElements(prhs[0]);
29 bp = mxGetPr(prhs[0]);
30
31 pbd = mxGetPr(prhs[2]);
32
33 if(nrhs < 5) maximize = 0;
34 else maximize = (int)mxGetScalar(prhs[4]);
35
36 if(siz_s == 0){
37 plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
38 sp = mxGetPr(plhs[0]);
39 if(maximize){
40 for(i=0; i<NB; i++){
41 *sp = (*sp < bp[i])? bp[i] : *sp;
42 }
43 }
44 else{
45 for(i=0; i<NB; i++){
46 *sp += bp[i];
47 }
48 }
49 return;
50 }
51
52 mask = malloc(siz_s * sizeof(int));
53 ssize = malloc(siz_s * sizeof(int));
54 count = 0;
55 for(i=0; i<siz_s; i++){
56 for(j=0; j<siz_b; j++){
57 if(ps[i] == pb[j]){
58 mask[count] = j;
59 count++;
60 break;
61 }
62 }
63 }
64
65 ndim = siz_b;
66 sx = (int *)malloc(sizeof(int)*ndim);
67 sy = (int *)malloc(sizeof(int)*ndim);
68 for(i=0; i<ndim; i++){
69 sx[i] = (int)pbd[i];
70 sy[i] = 1;
71 }
72 for(i=0; i<siz_s; i++){
73 temp = mask[i];
74 sy[temp] = sx[temp];
75 ssize[i] = sx[temp];
76 }
77
78 NS = 1;
79 for(i=0; i<ndim; i++){
80 NS *= sy[i];
81 }
82
83 plhs[0] = mxCreateNumericArray(siz_s, ssize, mxDOUBLE_CLASS, mxREAL);
84 sp = mxGetPr(plhs[0]);
85
86 if(NS == 1){
87 if(maximize){
88 for(i=0; i<NB; i++){
89 *sp = (*sp < bp[i])? bp[i] : *sp;
90 }
91 }
92 else{
93 for(i=0; i<NB; i++){
94 *sp += bp[i];
95 }
96 }
97 free(mask);
98 free(sx);
99 free(sy);
100 free(ssize);
101 return;
102 }
103
104 if(NS == NB){
105 for(i=0; i<NB; i++) *sp++ = *bp++;
106 free(mask);
107 free(sx);
108 free(sy);
109 free(ssize);
110 return;
111 }
112
113 s = (int *)malloc(sizeof(int)*ndim);
114 *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
115 subs = (int *)malloc(sizeof(int)*ndim);
116 cpsy2 = (int *)malloc(sizeof(int)*ndim);
117 for(i = 0; i < ndim; i++){
118 subs[i] = 0;
119 s[i] = sx[i] - 1;
120 }
121
122 for(i = 0; i < ndim-1; i++){
123 cpsy[i+1] = cpsy[i]*sy[i]--;
124 cpsy2[i] = cpsy[i]*sy[i];
125 }
126 cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
127
128 if(maximize){
129 for(j=0; j<NB; j++){
130 *sp = (*sp < *bp)? *bp : *sp;
131 bp++;
132 for(i = 0; i < ndim; i++){
133 if(subs[i] == s[i]){
134 subs[i] = 0;
135 if(sy[i])
136 sp -= cpsy2[i];
137 }
138 else{
139 subs[i]++;
140 if(sy[i])
141 sp += cpsy[i];
142 break;
143 }
144 }
145 }
146 }
147 else{
148 for(j=0; j<NB; j++){
149 *sp += *bp++;
150 for(i = 0; i < ndim; i++){
151 if(subs[i] == s[i]){
152 subs[i] = 0;
153 if(sy[i])
154 sp -= cpsy2[i];
155 }
156 else{
157 subs[i]++;
158 if(sy[i])
159 sp += cpsy[i];
160 break;
161 }
162 }
163 }
164 }
165
166 free(sx);
167 free(sy);
168 free(s);
169 free(cpsy);
170 free(subs);
171 free(cpsy2);
172 free(mask);
173 free(ssize);
174 }
175