wolffd@0
|
1 #include "mex.h"
|
wolffd@0
|
2
|
wolffd@0
|
3 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray
|
wolffd@0
|
4 *prhs[])
|
wolffd@0
|
5 {
|
wolffd@0
|
6 double *pn, *pi, *pj, *pm, *y, *ecElts, *pcpt, *famElts, *strideElts,
|
wolffd@0
|
7 *ev, *nsElts;
|
wolffd@0
|
8 int i, k, j, m, n;
|
wolffd@0
|
9 mxArray *ec, *cpt, *fam, *ns;
|
wolffd@0
|
10 int c1, famSize, nsj;
|
wolffd@0
|
11 int strideStride, startInd, stride, pos, numNodes;
|
wolffd@0
|
12
|
wolffd@0
|
13 const int BNET = 0;
|
wolffd@0
|
14 const int STATE = 1;
|
wolffd@0
|
15 const int STRIDES = 6;
|
wolffd@0
|
16 const int FAMILIES = 7;
|
wolffd@0
|
17 const int CPT = 8;
|
wolffd@0
|
18
|
wolffd@0
|
19 pn = mxGetPr(prhs[3]);
|
wolffd@0
|
20 n = (int) pn[0];
|
wolffd@0
|
21 pi = mxGetPr(prhs[2]);
|
wolffd@0
|
22 i = (int) pi[0];
|
wolffd@0
|
23 pj = mxGetPr(prhs[4]);
|
wolffd@0
|
24 j = (int) pj[0];
|
wolffd@0
|
25 pm = mxGetPr(prhs[5]);
|
wolffd@0
|
26 m = (int) pm[0];
|
wolffd@0
|
27 ev = mxGetPr(prhs[STATE]);
|
wolffd@0
|
28 ns = mxGetField (prhs[BNET], 0, "node_sizes");
|
wolffd@0
|
29 nsElts = mxGetPr (ns);
|
wolffd@0
|
30 numNodes = mxGetM(ns);
|
wolffd@0
|
31
|
wolffd@0
|
32 strideStride = mxGetM(prhs[STRIDES]);
|
wolffd@0
|
33 strideElts = mxGetPr(prhs[STRIDES]);
|
wolffd@0
|
34
|
wolffd@0
|
35
|
wolffd@0
|
36
|
wolffd@0
|
37 /* Treat the case n = 1 separately */
|
wolffd@0
|
38 if (pn[0] == 1) {
|
wolffd@0
|
39
|
wolffd@0
|
40 /* Get the appropriate CPT */
|
wolffd@0
|
41 ec = mxGetField (prhs[BNET], 0, "eclass1");
|
wolffd@0
|
42 ecElts = mxGetPr(ec);
|
wolffd@0
|
43 k = (int) ecElts[i-1];
|
wolffd@0
|
44 cpt = mxGetCell (prhs[8], k-1);
|
wolffd@0
|
45 pcpt = mxGetPr(cpt);
|
wolffd@0
|
46
|
wolffd@0
|
47 nsj = (int) nsElts[j-1];
|
wolffd@0
|
48
|
wolffd@0
|
49 /* Get the correct family vector */
|
wolffd@0
|
50 /* (Note : MEX is painful) */
|
wolffd@0
|
51 fam = mxGetCell (prhs[FAMILIES], i - 1);
|
wolffd@0
|
52 famSize = mxGetNumberOfElements(fam);
|
wolffd@0
|
53 famElts = mxGetPr(fam);
|
wolffd@0
|
54
|
wolffd@0
|
55
|
wolffd@0
|
56 /* Figure out starting position and stride */
|
wolffd@0
|
57 startInd = 0;
|
wolffd@0
|
58 for (c1 = 0, pos = k-1; c1 < famSize; c1++, pos+=strideStride) {
|
wolffd@0
|
59 if (famElts[c1] != j) {
|
wolffd@0
|
60 startInd += strideElts[pos]*(ev[(int)famElts[c1]-1]-1);
|
wolffd@0
|
61 }
|
wolffd@0
|
62 else {
|
wolffd@0
|
63 stride = strideElts[pos];
|
wolffd@0
|
64 }
|
wolffd@0
|
65 }
|
wolffd@0
|
66
|
wolffd@0
|
67 plhs[0] = mxCreateDoubleMatrix (1, nsj, mxREAL);
|
wolffd@0
|
68 y = mxGetPr(plhs[0]);
|
wolffd@0
|
69 for (c1 = 0, pos = startInd; c1 < nsj; c1++, pos+=stride) {
|
wolffd@0
|
70 y[c1] = pcpt[pos];
|
wolffd@0
|
71 }
|
wolffd@0
|
72 }
|
wolffd@0
|
73
|
wolffd@0
|
74 /* Handle the case n > 1 */
|
wolffd@0
|
75 else {
|
wolffd@0
|
76
|
wolffd@0
|
77 /* Get the appropriate CPT */
|
wolffd@0
|
78 ec = mxGetField (prhs[BNET], 0, "eclass2");
|
wolffd@0
|
79 ecElts = mxGetPr(ec);
|
wolffd@0
|
80 k = (int) ecElts[i-1];
|
wolffd@0
|
81 cpt = mxGetCell (prhs[8], k-1);
|
wolffd@0
|
82 pcpt = mxGetPr(cpt);
|
wolffd@0
|
83
|
wolffd@0
|
84 /* Figure out size of slice */
|
wolffd@0
|
85 if (m == 1) {
|
wolffd@0
|
86 nsj = (int) nsElts[j-1];
|
wolffd@0
|
87 }
|
wolffd@0
|
88 else {
|
wolffd@0
|
89 nsj = (int) nsElts[j-1+numNodes];
|
wolffd@0
|
90 }
|
wolffd@0
|
91
|
wolffd@0
|
92 /* Figure out family */
|
wolffd@0
|
93 fam = mxGetCell (prhs[FAMILIES], i - 1 + numNodes);
|
wolffd@0
|
94 famSize = mxGetNumberOfElements(fam);
|
wolffd@0
|
95 famElts = mxGetPr(fam);
|
wolffd@0
|
96
|
wolffd@0
|
97 startInd = 0;
|
wolffd@0
|
98 for (c1 = 0, pos = k-1; c1 < famSize; c1++, pos+=strideStride) {
|
wolffd@0
|
99 int f = (int) famElts[c1];
|
wolffd@0
|
100
|
wolffd@0
|
101 if (((f == j+numNodes) && (m == n)) || ((f == j) && (m ==
|
wolffd@0
|
102 n-1))) {
|
wolffd@0
|
103 stride = strideElts[pos];
|
wolffd@0
|
104 }
|
wolffd@0
|
105 else {
|
wolffd@0
|
106 startInd += strideElts[pos] * (ev[f-1+((n-2)*numNodes)]-1);
|
wolffd@0
|
107 }
|
wolffd@0
|
108 }
|
wolffd@0
|
109
|
wolffd@0
|
110 plhs[0] = mxCreateDoubleMatrix(1,nsj, mxREAL);
|
wolffd@0
|
111 y = mxGetPr(plhs[0]);
|
wolffd@0
|
112 for (c1 = 0, pos = startInd; c1 < nsj; c1++, pos+=stride) {
|
wolffd@0
|
113 y[c1] = pcpt[pos];
|
wolffd@0
|
114 }
|
wolffd@0
|
115 }
|
wolffd@0
|
116 }
|