File Coverage

pp-hmmfwq.c
Criterion Covered Total %
statement 0 105 0.0
branch 0 226 0.0
condition n/a
subroutine n/a
pod n/a
total 0 331 0.0


line stmt bran cond sub pod time code
1              
2             #line 453 "lib/PDL/PP.pm"
3             /*
4             * THIS FILE WAS GENERATED BY PDL::PP from HMM.pd! Do not modify!
5             */
6              
7             #define PDL_FREE_CODE(trans, destroy, comp_free_code, ntpriv_free_code) \
8             if (destroy) { \
9             comp_free_code \
10             } \
11             if ((trans)->dims_redone) { \
12             ntpriv_free_code \
13             }
14              
15             #include "EXTERN.h"
16             #include "perl.h"
17             #include "XSUB.h"
18             #include "pdl.h"
19             #include "pdlcore.h"
20             #define PDL PDL_HMM
21             extern Core* PDL; /* Structure hold core C functions */
22             #line 23 "pp-hmmfwq.c"
23              
24             #include
25              
26             /*#define DEBUG_ALPHA*/
27             /*#define DEBUG_BETA*/
28             /*#define DEBUG_VITERBI*/
29              
30              
31             /* logadd(x,y) = log(exp(x)+exp(y))
32             * + Code from Manning & Schütze (1997), Sec. 9.4, page 337
33             *
34             * LOG_BIG = log(1E31)
35             */
36             #define LOG_BIG 71.3801378828154
37             #define LOG_ZERO -1E+38
38             #define LOG_ONE 0
39             #define LOG_NONE 1
40 0           static inline double logadd1(double x, double y) {
41 0 0         if (y-x > LOG_BIG) return y;
42 0 0         else if (x-y > LOG_BIG) return x;
43             /*else return min(x,y) + log(exp(x-min(x,y)) + exp(y-min(x,y))); */
44 0 0         else if (x
45 0           else return y + log(exp(x-y) + 1);
46             }
47             static inline double logadd0(double x, double y) {
48             return log(exp(x)+exp(y));
49             }
50              
51             /* logdiff(x,y) = log(exp(x)-exp(y))
52             * + adapted from above
53             * + always returns positive (i.e. symmetric difference)
54             */
55             static inline double logdiff1(double x, double y) {
56             if (y-x > LOG_BIG) { return y; }
57             else if (x-y > LOG_BIG) { return x; }
58             /*else { return max(x,y) + log(exp(max(x,y)-max(x,y)) - exp(min(x,y)-max(x,y))); } */
59             /* = max(x,y) + log( 1 - exp(min(x,y)-max(x,y))); } */
60             else if (x>y) { return x + log( 1 - exp(y-x)); }
61             else { return y + log( 1 - exp(x-y)); }
62             }
63             static inline double logdiff0(double x, double y) {
64             return log(x>y ? (exp(x)-exp(y)) : (exp(y)-exp(x)));
65             }
66              
67             /*
68             #define logadd(x,y) logadd0(x,y)
69             #define logdiff(x,y) logdiff0(x,y)
70             */
71              
72             #define logadd(x,y) logadd1(x,y)
73             #define logdiff(x,y) logdiff1(x,y)
74              
75              
76             #line 1857 "lib/PDL/PP.pm"
77             pdl_error pdl_hmmfwq_readdata(pdl_trans *__privtrans) {
78             pdl_error PDL_err = {0, NULL, 0};
79             #line 80 "pp-hmmfwq.c"
80 0 0         if (!__privtrans->broadcast.incs) return PDL->make_error(PDL_EUSERERROR, "Error in hmmfwq:" "broadcast.incs NULL");
81             /* broadcastloop declarations */
82             int __brcloopval;
83             register PDL_Indx __tind0,__tind1; /* counters along dim */
84 0           register PDL_Indx __tnpdls = __privtrans->broadcast.npdls;
85             /* dims here are how many steps along those dims */
86 0           register PDL_Indx __tinc0_a = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,0);
87 0           register PDL_Indx __tinc0_b = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,0);
88 0           register PDL_Indx __tinc0_pi = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,0);
89 0           register PDL_Indx __tinc0_o = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,0);
90 0           register PDL_Indx __tinc0_oq = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,0);
91 0           register PDL_Indx __tinc0_alphaq = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,0);
92 0           register PDL_Indx __tinc1_a = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,1);
93 0           register PDL_Indx __tinc1_b = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,1);
94 0           register PDL_Indx __tinc1_pi = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,1);
95 0           register PDL_Indx __tinc1_o = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,1);
96 0           register PDL_Indx __tinc1_oq = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,1);
97 0           register PDL_Indx __tinc1_alphaq = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,1);
98             #define PDL_BROADCASTLOOP_START_hmmfwq_readdata PDL_BROADCASTLOOP_START( \
99             readdata, \
100             __privtrans->broadcast, \
101             __privtrans->vtable, \
102             a_datap += __offsp[0]; \
103             b_datap += __offsp[1]; \
104             pi_datap += __offsp[2]; \
105             o_datap += __offsp[3]; \
106             oq_datap += __offsp[4]; \
107             alphaq_datap += __offsp[5]; \
108             , \
109             ( ,a_datap += __tinc1_a - __tinc0_a * __tdims0 \
110             ,b_datap += __tinc1_b - __tinc0_b * __tdims0 \
111             ,pi_datap += __tinc1_pi - __tinc0_pi * __tdims0 \
112             ,o_datap += __tinc1_o - __tinc0_o * __tdims0 \
113             ,oq_datap += __tinc1_oq - __tinc0_oq * __tdims0 \
114             ,alphaq_datap += __tinc1_alphaq - __tinc0_alphaq * __tdims0 \
115             ), \
116             ( ,a_datap += __tinc0_a \
117             ,b_datap += __tinc0_b \
118             ,pi_datap += __tinc0_pi \
119             ,o_datap += __tinc0_o \
120             ,oq_datap += __tinc0_oq \
121             ,alphaq_datap += __tinc0_alphaq \
122             ) \
123             )
124             #define PDL_BROADCASTLOOP_END_hmmfwq_readdata PDL_BROADCASTLOOP_END( \
125             __privtrans->broadcast, \
126             a_datap -= __tinc1_a * __tdims1 + __offsp[0]; \
127             b_datap -= __tinc1_b * __tdims1 + __offsp[1]; \
128             pi_datap -= __tinc1_pi * __tdims1 + __offsp[2]; \
129             o_datap -= __tinc1_o * __tdims1 + __offsp[3]; \
130             oq_datap -= __tinc1_oq * __tdims1 + __offsp[4]; \
131             alphaq_datap -= __tinc1_alphaq * __tdims1 + __offsp[5]; \
132             )
133 0           register PDL_Indx __inc_a_N0 = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,0,0)]; (void)__inc_a_N0;register PDL_Indx __inc_a_N1 = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,0,1)]; (void)__inc_a_N1;
134 0           register PDL_Indx __inc_alphaq_Q = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,5,0)]; (void)__inc_alphaq_Q;register PDL_Indx __inc_alphaq_T = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,5,1)]; (void)__inc_alphaq_T;
135 0           register PDL_Indx __inc_b_N = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,1,0)]; (void)__inc_b_N;register PDL_Indx __inc_b_M = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,1,1)]; (void)__inc_b_M;
136 0           register PDL_Indx __inc_o_T = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,3,0)]; (void)__inc_o_T;
137 0           register PDL_Indx __inc_oq_Q = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,4,0)]; (void)__inc_oq_Q;register PDL_Indx __inc_oq_T = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,4,1)]; (void)__inc_oq_T;
138 0           register PDL_Indx __inc_pi_N = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,2,0)]; (void)__inc_pi_N;
139             #ifndef PDL_DECLARE_PARAMS_hmmfwq_1
140             #define PDL_DECLARE_PARAMS_hmmfwq_1(PDL_TYPE_OP,PDL_PPSYM_OP) \
141             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, a, (__privtrans->pdls[0]), 1, PDL_PPSYM_OP) \
142             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, b, (__privtrans->pdls[1]), 1, PDL_PPSYM_OP) \
143             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, pi, (__privtrans->pdls[2]), 1, PDL_PPSYM_OP) \
144             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, o, (__privtrans->pdls[3]), 1, PDL_PPSYM_OP) \
145             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, oq, (__privtrans->pdls[4]), 1, PDL_PPSYM_OP) \
146             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, alphaq, (__privtrans->pdls[5]), 1, PDL_PPSYM_OP)
147             #endif
148             #define PDL_IF_BAD(t,f) f
149 0           switch (__privtrans->__datatype) { /* Start generic switch */
150 0           case PDL_F: {
151 0 0         PDL_DECLARE_PARAMS_hmmfwq_1(PDL_Float,F)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
152 0 0         PDL_BROADCASTLOOP_START_hmmfwq_readdata {
    0          
    0          
    0          
    0          
    0          
    0          
153             /*-- Initialize: t==0 --*/
154 0           int i,j,t, o_tp1 = (o_datap)[0+(__inc_o_T*(0))];
155             int qi,qj;
156              
157 0 0         for (qi=0; qi < __privtrans->ind_sizes[2]; qi++) {
158 0           j = (oq_datap)[0+(__inc_oq_Q*(qi))+(__inc_oq_T*(0))];
159 0 0         (alphaq_datap)[0+(__inc_alphaq_Q*(qi))+(__inc_alphaq_T*(0))] = (j>=0 ? (pi_datap)[0+(__inc_pi_N*(j))] + (b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_tp1))] : (PDL_Float)LOG_ZERO);
160             }
161              
162             #ifdef DEBUG_ALPHA
163             printf("\n\n");
164             #endif
165              
166             /*-- Loop: time t>0 --*/
167 0 0         for (t=0; t < __privtrans->ind_sizes[3]-1; t++) {
168 0           o_tp1 = (o_datap)[0+(__inc_o_T*(t+1))];
169              
170             /*-- Loop: q_(t+1)=qj : state_(t+1)=oq(qj,t+1)=j --*/
171 0 0         for (qj=0; qj < __privtrans->ind_sizes[2]; qj++) {
172 0           PDL_Float alpha_j_tp1 = (PDL_Float)LOG_ZERO;
173 0           j = (oq_datap)[0+(__inc_oq_Q*(qj))+(__inc_oq_T*(t+1))];
174              
175             /*-- Loop: q_(t)=qi : state_(t)=oq(qi,t)=i --*/
176 0 0         for (qi=0; j>=0 && qi < __privtrans->ind_sizes[2]; qi++) {
    0          
177 0           i = (oq_datap)[0+(__inc_oq_Q*(qi))+(__inc_oq_T*(t))];
178 0 0         if (i < 0) break;
179              
180 0           alpha_j_tp1 = logadd( (alphaq_datap)[0+(__inc_alphaq_Q*(qi))+(__inc_alphaq_T*(t))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))], alpha_j_tp1 );
181              
182             #ifdef DEBUG_ALPHA
183             printf("qi=%u,i=%d,qj=%u,j=%d,t=%u,o=%d: alphaq(qi=%u,t=%u)=%.2e a(i=%d,j=%d)=%.2e b(j=%d,o=%d)=%.2e prod=%.2e sum=%.2e\n",
184             qi,i, qj,j, t,o_tp1,
185             i,t, exp((alphaq_datap)[0+(__inc_alphaq_Q*(qi))+(__inc_alphaq_T*(t))]),
186             i,j, ((i>=0 && j>=0) ? exp((a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))]) : 0),
187             j,o_tp1, ((i>=0 && j>=0) ? exp((b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_tp1))]) : 0),
188             ((i>=0 && j>=0) ? exp( (alphaq_datap)[0+(__inc_alphaq_Q*(qi))+(__inc_alphaq_T*(t))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))] ) : 0),
189             exp(alpha_j_tp1));
190             #endif
191             }
192             /*-- End Loop: q_(t)=qi : state_(t)=oq(qi,t)=i --*/
193              
194             /*-- Storage: alphaq(time=t+1, stateIndex=qj) --*/
195 0 0         if (j>=0) {
196 0           (alphaq_datap)[0+(__inc_alphaq_Q*(qj))+(__inc_alphaq_T*(t+1))] = alpha_j_tp1 + (b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_tp1))];
197             } else {
198 0           (alphaq_datap)[0+(__inc_alphaq_Q*(qj))+(__inc_alphaq_T*(t+1))] = (PDL_Float)LOG_ZERO;
199             }
200              
201             #ifdef DEBUG_ALPHA
202             printf("----> alphaq(qj=%u [j=%d], t=%u)=%.2E\n", qj,j,t+1, exp((alphaq_datap)[0+(__inc_alphaq_Q*(qj))+(__inc_alphaq_T*(t+1))]));
203             #endif
204             }
205             /*-- End Loop: q_(t+1)=qj : state_(t+1)=oq(qj,t+1)=j --*/
206              
207             #ifdef DEBUG_ALPHA
208             printf("\n\n");
209             #endif
210             }
211             /*-- End Loop: time t>0 --*/
212 0 0         }PDL_BROADCASTLOOP_END_hmmfwq_readdata
    0          
213 0           } break;
214 0           case PDL_D: {
215 0 0         PDL_DECLARE_PARAMS_hmmfwq_1(PDL_Double,D)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
216 0 0         PDL_BROADCASTLOOP_START_hmmfwq_readdata {
    0          
    0          
    0          
    0          
    0          
    0          
217             /*-- Initialize: t==0 --*/
218 0           int i,j,t, o_tp1 = (o_datap)[0+(__inc_o_T*(0))];
219             int qi,qj;
220              
221 0 0         for (qi=0; qi < __privtrans->ind_sizes[2]; qi++) {
222 0           j = (oq_datap)[0+(__inc_oq_Q*(qi))+(__inc_oq_T*(0))];
223 0 0         (alphaq_datap)[0+(__inc_alphaq_Q*(qi))+(__inc_alphaq_T*(0))] = (j>=0 ? (pi_datap)[0+(__inc_pi_N*(j))] + (b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_tp1))] : (PDL_Double)LOG_ZERO);
224             }
225              
226             #ifdef DEBUG_ALPHA
227             printf("\n\n");
228             #endif
229              
230             /*-- Loop: time t>0 --*/
231 0 0         for (t=0; t < __privtrans->ind_sizes[3]-1; t++) {
232 0           o_tp1 = (o_datap)[0+(__inc_o_T*(t+1))];
233              
234             /*-- Loop: q_(t+1)=qj : state_(t+1)=oq(qj,t+1)=j --*/
235 0 0         for (qj=0; qj < __privtrans->ind_sizes[2]; qj++) {
236 0           PDL_Double alpha_j_tp1 = (PDL_Double)LOG_ZERO;
237 0           j = (oq_datap)[0+(__inc_oq_Q*(qj))+(__inc_oq_T*(t+1))];
238              
239             /*-- Loop: q_(t)=qi : state_(t)=oq(qi,t)=i --*/
240 0 0         for (qi=0; j>=0 && qi < __privtrans->ind_sizes[2]; qi++) {
    0          
241 0           i = (oq_datap)[0+(__inc_oq_Q*(qi))+(__inc_oq_T*(t))];
242 0 0         if (i < 0) break;
243              
244 0           alpha_j_tp1 = logadd( (alphaq_datap)[0+(__inc_alphaq_Q*(qi))+(__inc_alphaq_T*(t))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))], alpha_j_tp1 );
245              
246             #ifdef DEBUG_ALPHA
247             printf("qi=%u,i=%d,qj=%u,j=%d,t=%u,o=%d: alphaq(qi=%u,t=%u)=%.2e a(i=%d,j=%d)=%.2e b(j=%d,o=%d)=%.2e prod=%.2e sum=%.2e\n",
248             qi,i, qj,j, t,o_tp1,
249             i,t, exp((alphaq_datap)[0+(__inc_alphaq_Q*(qi))+(__inc_alphaq_T*(t))]),
250             i,j, ((i>=0 && j>=0) ? exp((a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))]) : 0),
251             j,o_tp1, ((i>=0 && j>=0) ? exp((b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_tp1))]) : 0),
252             ((i>=0 && j>=0) ? exp( (alphaq_datap)[0+(__inc_alphaq_Q*(qi))+(__inc_alphaq_T*(t))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))] ) : 0),
253             exp(alpha_j_tp1));
254             #endif
255             }
256             /*-- End Loop: q_(t)=qi : state_(t)=oq(qi,t)=i --*/
257              
258             /*-- Storage: alphaq(time=t+1, stateIndex=qj) --*/
259 0 0         if (j>=0) {
260 0           (alphaq_datap)[0+(__inc_alphaq_Q*(qj))+(__inc_alphaq_T*(t+1))] = alpha_j_tp1 + (b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_tp1))];
261             } else {
262 0           (alphaq_datap)[0+(__inc_alphaq_Q*(qj))+(__inc_alphaq_T*(t+1))] = (PDL_Double)LOG_ZERO;
263             }
264              
265             #ifdef DEBUG_ALPHA
266             printf("----> alphaq(qj=%u [j=%d], t=%u)=%.2E\n", qj,j,t+1, exp((alphaq_datap)[0+(__inc_alphaq_Q*(qj))+(__inc_alphaq_T*(t+1))]));
267             #endif
268             }
269             /*-- End Loop: q_(t+1)=qj : state_(t+1)=oq(qj,t+1)=j --*/
270              
271             #ifdef DEBUG_ALPHA
272             printf("\n\n");
273             #endif
274             }
275             /*-- End Loop: time t>0 --*/
276 0 0         }PDL_BROADCASTLOOP_END_hmmfwq_readdata
    0          
277 0           } break;
278 0           case PDL_LD: {
279 0 0         PDL_DECLARE_PARAMS_hmmfwq_1(PDL_LDouble,E)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
280 0 0         PDL_BROADCASTLOOP_START_hmmfwq_readdata {
    0          
    0          
    0          
    0          
    0          
    0          
281             /*-- Initialize: t==0 --*/
282 0           int i,j,t, o_tp1 = (o_datap)[0+(__inc_o_T*(0))];
283             int qi,qj;
284              
285 0 0         for (qi=0; qi < __privtrans->ind_sizes[2]; qi++) {
286 0           j = (oq_datap)[0+(__inc_oq_Q*(qi))+(__inc_oq_T*(0))];
287 0 0         (alphaq_datap)[0+(__inc_alphaq_Q*(qi))+(__inc_alphaq_T*(0))] = (j>=0 ? (pi_datap)[0+(__inc_pi_N*(j))] + (b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_tp1))] : (PDL_LDouble)LOG_ZERO);
288             }
289              
290             #ifdef DEBUG_ALPHA
291             printf("\n\n");
292             #endif
293              
294             /*-- Loop: time t>0 --*/
295 0 0         for (t=0; t < __privtrans->ind_sizes[3]-1; t++) {
296 0           o_tp1 = (o_datap)[0+(__inc_o_T*(t+1))];
297              
298             /*-- Loop: q_(t+1)=qj : state_(t+1)=oq(qj,t+1)=j --*/
299 0 0         for (qj=0; qj < __privtrans->ind_sizes[2]; qj++) {
300 0           PDL_LDouble alpha_j_tp1 = (PDL_LDouble)LOG_ZERO;
301 0           j = (oq_datap)[0+(__inc_oq_Q*(qj))+(__inc_oq_T*(t+1))];
302              
303             /*-- Loop: q_(t)=qi : state_(t)=oq(qi,t)=i --*/
304 0 0         for (qi=0; j>=0 && qi < __privtrans->ind_sizes[2]; qi++) {
    0          
305 0           i = (oq_datap)[0+(__inc_oq_Q*(qi))+(__inc_oq_T*(t))];
306 0 0         if (i < 0) break;
307              
308 0           alpha_j_tp1 = logadd( (alphaq_datap)[0+(__inc_alphaq_Q*(qi))+(__inc_alphaq_T*(t))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))], alpha_j_tp1 );
309              
310             #ifdef DEBUG_ALPHA
311             printf("qi=%u,i=%d,qj=%u,j=%d,t=%u,o=%d: alphaq(qi=%u,t=%u)=%.2e a(i=%d,j=%d)=%.2e b(j=%d,o=%d)=%.2e prod=%.2e sum=%.2e\n",
312             qi,i, qj,j, t,o_tp1,
313             i,t, exp((alphaq_datap)[0+(__inc_alphaq_Q*(qi))+(__inc_alphaq_T*(t))]),
314             i,j, ((i>=0 && j>=0) ? exp((a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))]) : 0),
315             j,o_tp1, ((i>=0 && j>=0) ? exp((b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_tp1))]) : 0),
316             ((i>=0 && j>=0) ? exp( (alphaq_datap)[0+(__inc_alphaq_Q*(qi))+(__inc_alphaq_T*(t))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))] ) : 0),
317             exp(alpha_j_tp1));
318             #endif
319             }
320             /*-- End Loop: q_(t)=qi : state_(t)=oq(qi,t)=i --*/
321              
322             /*-- Storage: alphaq(time=t+1, stateIndex=qj) --*/
323 0 0         if (j>=0) {
324 0           (alphaq_datap)[0+(__inc_alphaq_Q*(qj))+(__inc_alphaq_T*(t+1))] = alpha_j_tp1 + (b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_tp1))];
325             } else {
326 0           (alphaq_datap)[0+(__inc_alphaq_Q*(qj))+(__inc_alphaq_T*(t+1))] = (PDL_LDouble)LOG_ZERO;
327             }
328              
329             #ifdef DEBUG_ALPHA
330             printf("----> alphaq(qj=%u [j=%d], t=%u)=%.2E\n", qj,j,t+1, exp((alphaq_datap)[0+(__inc_alphaq_Q*(qj))+(__inc_alphaq_T*(t+1))]));
331             #endif
332             }
333             /*-- End Loop: q_(t+1)=qj : state_(t+1)=oq(qj,t+1)=j --*/
334              
335             #ifdef DEBUG_ALPHA
336             printf("\n\n");
337             #endif
338             }
339             /*-- End Loop: time t>0 --*/
340 0 0         }PDL_BROADCASTLOOP_END_hmmfwq_readdata
    0          
341 0           } break;
342 0           default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in hmmfwq: unhandled datatype(%d), only handles (FDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
343             }
344             #undef PDL_IF_BAD
345 0           return PDL_err;
346             }
347              
348             static pdl_datatypes pdl_hmmfwq_vtable_gentypes[] = { PDL_F, PDL_D, PDL_LD, -1 };
349             static PDL_Indx pdl_hmmfwq_vtable_realdims[] = { 2, 2, 1, 1, 2, 2 };
350             static char *pdl_hmmfwq_vtable_parnames[] = { "a","b","pi","o","oq","alphaq" };
351             static short pdl_hmmfwq_vtable_parflags[] = {
352             0,
353             0,
354             0,
355             0,
356             0,
357             PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISWRITE
358             };
359             static pdl_datatypes pdl_hmmfwq_vtable_partypes[] = { -1, -1, -1, -1, -1, -1 };
360             static PDL_Indx pdl_hmmfwq_vtable_realdims_starts[] = { 0, 2, 4, 5, 6, 8 };
361             static PDL_Indx pdl_hmmfwq_vtable_realdims_ind_ids[] = { 1, 1, 1, 0, 1, 3, 2, 3, 2, 3 };
362             static char *pdl_hmmfwq_vtable_indnames[] = { "M","N","Q","T" };
363             pdl_transvtable pdl_hmmfwq_vtable = {
364             PDL_TRANS_DO_BROADCAST, 0, pdl_hmmfwq_vtable_gentypes, 5, 6, NULL /*CORE21*/,
365             pdl_hmmfwq_vtable_realdims, pdl_hmmfwq_vtable_parnames,
366             pdl_hmmfwq_vtable_parflags, pdl_hmmfwq_vtable_partypes,
367             pdl_hmmfwq_vtable_realdims_starts, pdl_hmmfwq_vtable_realdims_ind_ids, 10,
368             4, pdl_hmmfwq_vtable_indnames,
369             NULL, pdl_hmmfwq_readdata, NULL,
370             NULL,
371             0,"PDL::HMM::hmmfwq"
372             };
373              
374              
375 0           pdl_error pdl_run_hmmfwq(pdl *a,pdl *b,pdl *pi,pdl *o,pdl *oq,pdl *alphaq) {
376 0           pdl_error PDL_err = {0, NULL, 0};
377 0 0         if (!PDL) return (pdl_error){PDL_EFATAL, "PDL core struct is NULL, can't continue",0};
378 0           pdl_trans *__privtrans = PDL->create_trans(&pdl_hmmfwq_vtable);
379 0 0         if (!__privtrans) return PDL->make_error_simple(PDL_EFATAL, "Couldn't create trans");
380 0           __privtrans->pdls[0] = a;
381 0           __privtrans->pdls[1] = b;
382 0           __privtrans->pdls[2] = pi;
383 0           __privtrans->pdls[3] = o;
384 0           __privtrans->pdls[4] = oq;
385 0           __privtrans->pdls[5] = alphaq;
386 0 0         PDL_RETERROR(PDL_err, PDL->type_coerce(__privtrans));
387 0 0         PDL_RETERROR(PDL_err, PDL->make_trans_mutual(__privtrans));
388 0           return PDL_err;
389             }