0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 #include <math.h>
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 #include <stdio.h>
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 #include <stdlib.h>
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 #include <string.h>
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 #include <mex.h>
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 #include "matrix.h"
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 #include "version.h"
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 #include "ltpda_ssmsim.h"
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 #define DEBUG 0
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 /*
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 * A mex file to propagate an input signal with the given SS model
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 *
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 * M Hewitson 19-08-10
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 *
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 * $Id: ltpda_ssmsim.c,v 1.4 2010/08/23 11:56:09 ingo Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 */
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 /*
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 * function [y,lx] = ltpda_ssmsim(lastX, A.', Coutputs.', Cstates.', Baos.', Daos.', input);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 */
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 {
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 /* outputs */
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 double *y, *yptr;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 /* inputs */
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31 double *input, *iptr;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 double *Daos, *Dptr;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 double *Baos, *Bptr;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 double *Cstates;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35 double *Coutputs, *Coptr;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 double *A, *Aptr;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 double *lastX;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 double *SSini;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 double *tmpX;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 mwSize Ninputs, Nsamples, Nstates, Nstatesout, Noutputs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 int kk,jj,ll;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43 int mm,nn;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 int ki;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 /* parse input functions */
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 if( (nrhs == 0) || (nlhs == 0) )
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 {
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 print_usage(VERSION);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52 }
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54 if( (nrhs == 7) && (nlhs == 2) )/* let's go */
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 {
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 /*----------------- set inputs*/
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57 SSini = mxGetPr(prhs[0]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 A = mxGetPr(prhs[1]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 Coutputs = mxGetPr(prhs[2]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60 Cstates = mxGetPr(prhs[3]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 Baos = mxGetPr(prhs[4]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 Daos = mxGetPr(prhs[5]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63 input = mxGetPr(prhs[6]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65 Ninputs = mxGetM(prhs[4]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 Nsamples = mxGetN(prhs[6]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 Nstates = mxGetM(prhs[1]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 Nstatesout = mxGetM(prhs[3]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 Noutputs = mxGetN(prhs[2]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 #if DEBUG
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72 mexPrintf("Ninputs: %d\n", Ninputs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 mexPrintf("Nsamples: %d\n", Nsamples);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 mexPrintf("Nstates: %d\n", Nstates);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 mexPrintf("Nstatesout: %d\n", Nstatesout);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 mexPrintf("Noutputs: %d\n", Noutputs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 mexPrintf("N Coutputs: %d\n", mxGetNumberOfElements(prhs[2]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 mexPrintf("input: %dx%d\n", mxGetM(prhs[6]), mxGetN(prhs[6]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81 mexPrintf("D: %dx%d\n", mxGetN(prhs[5]), mxGetM(prhs[5]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 #endif
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 /* output y */
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85 plhs[0] = mxCreateDoubleMatrix(Noutputs, Nsamples, mxREAL);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 y = mxGetPr(plhs[0]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 /* output state vector*/
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 plhs[1] = mxCreateDoubleMatrix(Nstates, 1, mxREAL);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90 lastX = mxGetPr(plhs[1]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92 tmpX = (double*)calloc(Nstates, sizeof(double));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 memcpy(lastX, SSini, Nstates*sizeof(double));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95 /* do the business */
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96 yptr = &(y[0]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97 for (kk=0; kk<Nsamples; kk++) {
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 ki = kk*Ninputs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101 /* observation equation */
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102 Coptr = &(Coutputs[0]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103 Dptr = &(Daos[0]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
104 for (jj=0; jj<Noutputs; jj++) {
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
105 *yptr = 0.0;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
106 for (ll=0; ll<Nstates; ll++) {
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
107 *yptr += *Coptr * lastX[ll] ;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
108 Coptr++;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
109 }
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
110 iptr = &(input[ki]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
111 for (ll=0; ll<Ninputs; ll++) {
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
112 *yptr += *Dptr * (*iptr);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
113 Dptr++;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
114 iptr++;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
115 }
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
116 yptr++;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
117 }
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
118
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
119 /* state propagation */
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
120 memcpy(tmpX, lastX, Nstates*sizeof(double));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
121 Bptr = &(Baos[0]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
122 Aptr = &(A[0]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
123 for (jj=0; jj<Nstates; jj++) {
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
124 lastX[jj] = 0;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
125 for (ll=0; ll<Nstates; ll++) {
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
126 lastX[jj] += *Aptr * tmpX[ll];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
127 Aptr++;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
128 }
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
129 iptr = &(input[ki]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
130 for (ll=0; ll<Ninputs; ll++) {
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
131 lastX[jj] += *Bptr * (*iptr);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
132 Bptr++;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
133 iptr++;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
134 }
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
135 }
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
136
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
137
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
138 } /* end sample loop */
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
139
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
140
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
141 free(tmpX);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
142
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
143
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
144 }
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
145 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
146 {
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
147 print_usage(VERSION);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
148 }
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
149 }
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
150
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
151 void print_usage(char *version)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
152 {
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
153 mexPrintf("ltpda_ssmsim version %s\n", version);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
154 mexPrintf(" usage: [y,lx] = ltpda_ssmsim(lastX, A.', Coutputs.', Cstates.', Baos.', Daos.', input);");
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
155 mexErrMsgTxt("### incorrect usage");
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
156 }
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
157 |