File Coverage

pp-hmmviterbiq.c
Criterion Covered Total %
statement 0 125 0.0
branch 0 256 0.0
condition n/a
subroutine n/a
pod n/a
total 0 381 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-hmmviterbiq.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             static inline double logadd1(double x, double y) {
41             if (y-x > LOG_BIG) return y;
42             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             else if (x
45             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_hmmviterbiq_readdata(pdl_trans *__privtrans) {
78             pdl_error PDL_err = {0, NULL, 0};
79             #line 80 "pp-hmmviterbiq.c"
80 0 0         if (!__privtrans->broadcast.incs) return PDL->make_error(PDL_EUSERERROR, "Error in hmmviterbiq:" "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_deltaq = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,0);
92 0           register PDL_Indx __tinc0_psiq = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,6,0);
93 0           register PDL_Indx __tinc1_a = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,1);
94 0           register PDL_Indx __tinc1_b = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,1);
95 0           register PDL_Indx __tinc1_pi = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,1);
96 0           register PDL_Indx __tinc1_o = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,1);
97 0           register PDL_Indx __tinc1_oq = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,1);
98 0           register PDL_Indx __tinc1_deltaq = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,1);
99 0           register PDL_Indx __tinc1_psiq = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,6,1);
100             #define PDL_BROADCASTLOOP_START_hmmviterbiq_readdata PDL_BROADCASTLOOP_START( \
101             readdata, \
102             __privtrans->broadcast, \
103             __privtrans->vtable, \
104             a_datap += __offsp[0]; \
105             b_datap += __offsp[1]; \
106             pi_datap += __offsp[2]; \
107             o_datap += __offsp[3]; \
108             oq_datap += __offsp[4]; \
109             deltaq_datap += __offsp[5]; \
110             psiq_datap += __offsp[6]; \
111             , \
112             ( ,a_datap += __tinc1_a - __tinc0_a * __tdims0 \
113             ,b_datap += __tinc1_b - __tinc0_b * __tdims0 \
114             ,pi_datap += __tinc1_pi - __tinc0_pi * __tdims0 \
115             ,o_datap += __tinc1_o - __tinc0_o * __tdims0 \
116             ,oq_datap += __tinc1_oq - __tinc0_oq * __tdims0 \
117             ,deltaq_datap += __tinc1_deltaq - __tinc0_deltaq * __tdims0 \
118             ,psiq_datap += __tinc1_psiq - __tinc0_psiq * __tdims0 \
119             ), \
120             ( ,a_datap += __tinc0_a \
121             ,b_datap += __tinc0_b \
122             ,pi_datap += __tinc0_pi \
123             ,o_datap += __tinc0_o \
124             ,oq_datap += __tinc0_oq \
125             ,deltaq_datap += __tinc0_deltaq \
126             ,psiq_datap += __tinc0_psiq \
127             ) \
128             )
129             #define PDL_BROADCASTLOOP_END_hmmviterbiq_readdata PDL_BROADCASTLOOP_END( \
130             __privtrans->broadcast, \
131             a_datap -= __tinc1_a * __tdims1 + __offsp[0]; \
132             b_datap -= __tinc1_b * __tdims1 + __offsp[1]; \
133             pi_datap -= __tinc1_pi * __tdims1 + __offsp[2]; \
134             o_datap -= __tinc1_o * __tdims1 + __offsp[3]; \
135             oq_datap -= __tinc1_oq * __tdims1 + __offsp[4]; \
136             deltaq_datap -= __tinc1_deltaq * __tdims1 + __offsp[5]; \
137             psiq_datap -= __tinc1_psiq * __tdims1 + __offsp[6]; \
138             )
139 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;
140 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;
141 0           register PDL_Indx __inc_deltaq_Q = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,5,0)]; (void)__inc_deltaq_Q;register PDL_Indx __inc_deltaq_T = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,5,1)]; (void)__inc_deltaq_T;
142 0           register PDL_Indx __inc_o_T = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,3,0)]; (void)__inc_o_T;
143 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;
144 0           register PDL_Indx __inc_pi_N = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,2,0)]; (void)__inc_pi_N;
145 0           register PDL_Indx __inc_psiq_Q = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,6,0)]; (void)__inc_psiq_Q;register PDL_Indx __inc_psiq_T = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,6,1)]; (void)__inc_psiq_T;
146             #ifndef PDL_DECLARE_PARAMS_hmmviterbiq_1
147             #define PDL_DECLARE_PARAMS_hmmviterbiq_1(PDL_TYPE_OP,PDL_PPSYM_OP,PDL_TYPE_PARAM_psiq,PDL_PPSYM_PARAM_psiq) \
148             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, a, (__privtrans->pdls[0]), 1, PDL_PPSYM_OP) \
149             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, b, (__privtrans->pdls[1]), 1, PDL_PPSYM_OP) \
150             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, pi, (__privtrans->pdls[2]), 1, PDL_PPSYM_OP) \
151             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, o, (__privtrans->pdls[3]), 1, PDL_PPSYM_OP) \
152             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, oq, (__privtrans->pdls[4]), 1, PDL_PPSYM_OP) \
153             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, deltaq, (__privtrans->pdls[5]), 1, PDL_PPSYM_OP) \
154             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_psiq, psiq, (__privtrans->pdls[6]), 1, PDL_PPSYM_PARAM_psiq)
155             #endif
156             #define PDL_IF_BAD(t,f) f
157 0           switch (__privtrans->__datatype) { /* Start generic switch */
158 0           case PDL_F: {
159 0 0         PDL_DECLARE_PARAMS_hmmviterbiq_1(PDL_Float,F,PDL_Indx,N)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
160 0 0         PDL_BROADCASTLOOP_START_hmmviterbiq_readdata {
    0          
    0          
    0          
    0          
    0          
    0          
161             int qi,qj, i,j, t, o_t;
162             double deltaq_jt, deltaq_tmp;
163             int psiq_jt;
164              
165             /*-- Initialize: t=0: Loop: q_(0)=qi: state_(0)=oq(qi,0)=i --*/
166 0           o_t = (o_datap)[0+(__inc_o_T*(0))];
167 0 0         for (qi=0; qi<__privtrans->ind_sizes[2]; qi++) {
168 0           i = (oq_datap)[0+(__inc_oq_Q*(qi))+(__inc_oq_T*(0))];
169 0           (psiq_datap)[0+(__inc_psiq_Q*(qi))+(__inc_psiq_T*(0))] = 0;
170 0 0         (deltaq_datap)[0+(__inc_deltaq_Q*(qi))+(__inc_deltaq_T*(0))] = (i>=0 ? ((pi_datap)[0+(__inc_pi_N*(i))]+(b_datap)[0+(__inc_b_N*(i))+(__inc_b_M*(o_t))]) : (PDL_Float)LOG_ZERO);
171             }
172              
173             /*-- Loop: t>0: Loop: time==t --*/
174 0 0         for (t=1; t<__privtrans->ind_sizes[3]; t++) {
175 0           o_t = (o_datap)[0+(__inc_o_T*(t))];
176              
177             /*-- Loop: t>0: q_(t)=qj : state_(t)=oq(qj,t)=j --*/
178 0 0         for (qj=0; qj<__privtrans->ind_sizes[2]; qj++) {
179 0           j = (oq_datap)[0+(__inc_oq_Q*(qj))+(__inc_oq_T*(t))];
180 0           i = (oq_datap)[0+(__inc_oq_Q*(0))+(__inc_oq_T*(t-1))];
181 0           psiq_jt = 0;
182              
183 0 0         if (j >= 0 && i >=0) {
    0          
184 0           deltaq_jt = (deltaq_datap)[0+(__inc_deltaq_Q*(0))+(__inc_deltaq_T*(t-1))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))];
185             } else {
186 0           deltaq_jt = (deltaq_datap)[0+(__inc_deltaq_Q*(0))+(__inc_deltaq_T*(t-1))] + LOG_ZERO;
187             }
188              
189             /*-- Loop: t>0: q_(t-1)=qi : state_(t-1)=oq(qi,t)=i --*/
190 0 0         for (qi=1; qi<__privtrans->ind_sizes[2]; qi++) {
191 0           i = (oq_datap)[0+(__inc_oq_Q*(qi))+(__inc_oq_T*(t-1))];
192 0 0         if (j < 0 || i < 0) break;
    0          
193              
194 0           deltaq_tmp = (deltaq_datap)[0+(__inc_deltaq_Q*(qi))+(__inc_deltaq_T*(t-1))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))];
195              
196 0 0         if (deltaq_tmp > deltaq_jt) {
197 0           deltaq_jt = deltaq_tmp;
198 0           psiq_jt = qi;
199             }
200              
201             }
202             /*-- End Loop: t>0: q_(t-1)=qi : state_(t-1)=oq(qi,t)=i --*/
203              
204             /*-- Main: t>0: Store data for stateIndex,time=(qj,t) --*/
205 0 0         (deltaq_datap)[0+(__inc_deltaq_Q*(qj))+(__inc_deltaq_T*(t))] = deltaq_jt + (j>=0 ? (b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_t))] : LOG_ZERO);
206 0           (psiq_datap)[0+(__inc_psiq_Q*(qj))+(__inc_psiq_T*(t))] = psiq_jt;
207              
208             }
209             /*-- Loop: t>0: q_(t)=qj : state_(t)=oq(qj,t)=j --*/
210              
211             }
212             /*-- End Loop: t>0: Loop: time==t --*/
213 0 0         }PDL_BROADCASTLOOP_END_hmmviterbiq_readdata
    0          
214 0           } break;
215 0           case PDL_D: {
216 0 0         PDL_DECLARE_PARAMS_hmmviterbiq_1(PDL_Double,D,PDL_Indx,N)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
217 0 0         PDL_BROADCASTLOOP_START_hmmviterbiq_readdata {
    0          
    0          
    0          
    0          
    0          
    0          
218             int qi,qj, i,j, t, o_t;
219             double deltaq_jt, deltaq_tmp;
220             int psiq_jt;
221              
222             /*-- Initialize: t=0: Loop: q_(0)=qi: state_(0)=oq(qi,0)=i --*/
223 0           o_t = (o_datap)[0+(__inc_o_T*(0))];
224 0 0         for (qi=0; qi<__privtrans->ind_sizes[2]; qi++) {
225 0           i = (oq_datap)[0+(__inc_oq_Q*(qi))+(__inc_oq_T*(0))];
226 0           (psiq_datap)[0+(__inc_psiq_Q*(qi))+(__inc_psiq_T*(0))] = 0;
227 0 0         (deltaq_datap)[0+(__inc_deltaq_Q*(qi))+(__inc_deltaq_T*(0))] = (i>=0 ? ((pi_datap)[0+(__inc_pi_N*(i))]+(b_datap)[0+(__inc_b_N*(i))+(__inc_b_M*(o_t))]) : (PDL_Double)LOG_ZERO);
228             }
229              
230             /*-- Loop: t>0: Loop: time==t --*/
231 0 0         for (t=1; t<__privtrans->ind_sizes[3]; t++) {
232 0           o_t = (o_datap)[0+(__inc_o_T*(t))];
233              
234             /*-- Loop: t>0: q_(t)=qj : state_(t)=oq(qj,t)=j --*/
235 0 0         for (qj=0; qj<__privtrans->ind_sizes[2]; qj++) {
236 0           j = (oq_datap)[0+(__inc_oq_Q*(qj))+(__inc_oq_T*(t))];
237 0           i = (oq_datap)[0+(__inc_oq_Q*(0))+(__inc_oq_T*(t-1))];
238 0           psiq_jt = 0;
239              
240 0 0         if (j >= 0 && i >=0) {
    0          
241 0           deltaq_jt = (deltaq_datap)[0+(__inc_deltaq_Q*(0))+(__inc_deltaq_T*(t-1))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))];
242             } else {
243 0           deltaq_jt = (deltaq_datap)[0+(__inc_deltaq_Q*(0))+(__inc_deltaq_T*(t-1))] + LOG_ZERO;
244             }
245              
246             /*-- Loop: t>0: q_(t-1)=qi : state_(t-1)=oq(qi,t)=i --*/
247 0 0         for (qi=1; qi<__privtrans->ind_sizes[2]; qi++) {
248 0           i = (oq_datap)[0+(__inc_oq_Q*(qi))+(__inc_oq_T*(t-1))];
249 0 0         if (j < 0 || i < 0) break;
    0          
250              
251 0           deltaq_tmp = (deltaq_datap)[0+(__inc_deltaq_Q*(qi))+(__inc_deltaq_T*(t-1))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))];
252              
253 0 0         if (deltaq_tmp > deltaq_jt) {
254 0           deltaq_jt = deltaq_tmp;
255 0           psiq_jt = qi;
256             }
257              
258             }
259             /*-- End Loop: t>0: q_(t-1)=qi : state_(t-1)=oq(qi,t)=i --*/
260              
261             /*-- Main: t>0: Store data for stateIndex,time=(qj,t) --*/
262 0 0         (deltaq_datap)[0+(__inc_deltaq_Q*(qj))+(__inc_deltaq_T*(t))] = deltaq_jt + (j>=0 ? (b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_t))] : LOG_ZERO);
263 0           (psiq_datap)[0+(__inc_psiq_Q*(qj))+(__inc_psiq_T*(t))] = psiq_jt;
264              
265             }
266             /*-- Loop: t>0: q_(t)=qj : state_(t)=oq(qj,t)=j --*/
267              
268             }
269             /*-- End Loop: t>0: Loop: time==t --*/
270 0 0         }PDL_BROADCASTLOOP_END_hmmviterbiq_readdata
    0          
271 0           } break;
272 0           case PDL_LD: {
273 0 0         PDL_DECLARE_PARAMS_hmmviterbiq_1(PDL_LDouble,E,PDL_Indx,N)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
274 0 0         PDL_BROADCASTLOOP_START_hmmviterbiq_readdata {
    0          
    0          
    0          
    0          
    0          
    0          
275             int qi,qj, i,j, t, o_t;
276             double deltaq_jt, deltaq_tmp;
277             int psiq_jt;
278              
279             /*-- Initialize: t=0: Loop: q_(0)=qi: state_(0)=oq(qi,0)=i --*/
280 0           o_t = (o_datap)[0+(__inc_o_T*(0))];
281 0 0         for (qi=0; qi<__privtrans->ind_sizes[2]; qi++) {
282 0           i = (oq_datap)[0+(__inc_oq_Q*(qi))+(__inc_oq_T*(0))];
283 0           (psiq_datap)[0+(__inc_psiq_Q*(qi))+(__inc_psiq_T*(0))] = 0;
284 0 0         (deltaq_datap)[0+(__inc_deltaq_Q*(qi))+(__inc_deltaq_T*(0))] = (i>=0 ? ((pi_datap)[0+(__inc_pi_N*(i))]+(b_datap)[0+(__inc_b_N*(i))+(__inc_b_M*(o_t))]) : (PDL_LDouble)LOG_ZERO);
285             }
286              
287             /*-- Loop: t>0: Loop: time==t --*/
288 0 0         for (t=1; t<__privtrans->ind_sizes[3]; t++) {
289 0           o_t = (o_datap)[0+(__inc_o_T*(t))];
290              
291             /*-- Loop: t>0: q_(t)=qj : state_(t)=oq(qj,t)=j --*/
292 0 0         for (qj=0; qj<__privtrans->ind_sizes[2]; qj++) {
293 0           j = (oq_datap)[0+(__inc_oq_Q*(qj))+(__inc_oq_T*(t))];
294 0           i = (oq_datap)[0+(__inc_oq_Q*(0))+(__inc_oq_T*(t-1))];
295 0           psiq_jt = 0;
296              
297 0 0         if (j >= 0 && i >=0) {
    0          
298 0           deltaq_jt = (deltaq_datap)[0+(__inc_deltaq_Q*(0))+(__inc_deltaq_T*(t-1))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))];
299             } else {
300 0           deltaq_jt = (deltaq_datap)[0+(__inc_deltaq_Q*(0))+(__inc_deltaq_T*(t-1))] + LOG_ZERO;
301             }
302              
303             /*-- Loop: t>0: q_(t-1)=qi : state_(t-1)=oq(qi,t)=i --*/
304 0 0         for (qi=1; qi<__privtrans->ind_sizes[2]; qi++) {
305 0           i = (oq_datap)[0+(__inc_oq_Q*(qi))+(__inc_oq_T*(t-1))];
306 0 0         if (j < 0 || i < 0) break;
    0          
307              
308 0           deltaq_tmp = (deltaq_datap)[0+(__inc_deltaq_Q*(qi))+(__inc_deltaq_T*(t-1))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))];
309              
310 0 0         if (deltaq_tmp > deltaq_jt) {
311 0           deltaq_jt = deltaq_tmp;
312 0           psiq_jt = qi;
313             }
314              
315             }
316             /*-- End Loop: t>0: q_(t-1)=qi : state_(t-1)=oq(qi,t)=i --*/
317              
318             /*-- Main: t>0: Store data for stateIndex,time=(qj,t) --*/
319 0 0         (deltaq_datap)[0+(__inc_deltaq_Q*(qj))+(__inc_deltaq_T*(t))] = deltaq_jt + (j>=0 ? (b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_t))] : LOG_ZERO);
320 0           (psiq_datap)[0+(__inc_psiq_Q*(qj))+(__inc_psiq_T*(t))] = psiq_jt;
321              
322             }
323             /*-- Loop: t>0: q_(t)=qj : state_(t)=oq(qj,t)=j --*/
324              
325             }
326             /*-- End Loop: t>0: Loop: time==t --*/
327 0 0         }PDL_BROADCASTLOOP_END_hmmviterbiq_readdata
    0          
328 0           } break;
329 0           default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in hmmviterbiq: unhandled datatype(%d), only handles (FDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
330             }
331             #undef PDL_IF_BAD
332 0           return PDL_err;
333             }
334              
335             static pdl_datatypes pdl_hmmviterbiq_vtable_gentypes[] = { PDL_F, PDL_D, PDL_LD, -1 };
336             static PDL_Indx pdl_hmmviterbiq_vtable_realdims[] = { 2, 2, 1, 1, 2, 2, 2 };
337             static char *pdl_hmmviterbiq_vtable_parnames[] = { "a","b","pi","o","oq","deltaq","psiq" };
338             static short pdl_hmmviterbiq_vtable_parflags[] = {
339             0,
340             0,
341             0,
342             0,
343             0,
344             PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISWRITE,
345             PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISTYPED|PDL_PARAM_ISWRITE
346             };
347             static pdl_datatypes pdl_hmmviterbiq_vtable_partypes[] = { -1, -1, -1, -1, -1, -1, PDL_IND };
348             static PDL_Indx pdl_hmmviterbiq_vtable_realdims_starts[] = { 0, 2, 4, 5, 6, 8, 10 };
349             static PDL_Indx pdl_hmmviterbiq_vtable_realdims_ind_ids[] = { 1, 1, 1, 0, 1, 3, 2, 3, 2, 3, 2, 3 };
350             static char *pdl_hmmviterbiq_vtable_indnames[] = { "M","N","Q","T" };
351             pdl_transvtable pdl_hmmviterbiq_vtable = {
352             PDL_TRANS_DO_BROADCAST, 0, pdl_hmmviterbiq_vtable_gentypes, 5, 7, NULL /*CORE21*/,
353             pdl_hmmviterbiq_vtable_realdims, pdl_hmmviterbiq_vtable_parnames,
354             pdl_hmmviterbiq_vtable_parflags, pdl_hmmviterbiq_vtable_partypes,
355             pdl_hmmviterbiq_vtable_realdims_starts, pdl_hmmviterbiq_vtable_realdims_ind_ids, 12,
356             4, pdl_hmmviterbiq_vtable_indnames,
357             NULL, pdl_hmmviterbiq_readdata, NULL,
358             NULL,
359             0,"PDL::HMM::hmmviterbiq"
360             };
361              
362              
363 0           pdl_error pdl_run_hmmviterbiq(pdl *a,pdl *b,pdl *pi,pdl *o,pdl *oq,pdl *deltaq,pdl *psiq) {
364 0           pdl_error PDL_err = {0, NULL, 0};
365 0 0         if (!PDL) return (pdl_error){PDL_EFATAL, "PDL core struct is NULL, can't continue",0};
366 0           pdl_trans *__privtrans = PDL->create_trans(&pdl_hmmviterbiq_vtable);
367 0 0         if (!__privtrans) return PDL->make_error_simple(PDL_EFATAL, "Couldn't create trans");
368 0           __privtrans->pdls[0] = a;
369 0           __privtrans->pdls[1] = b;
370 0           __privtrans->pdls[2] = pi;
371 0           __privtrans->pdls[3] = o;
372 0           __privtrans->pdls[4] = oq;
373 0           __privtrans->pdls[5] = deltaq;
374 0           __privtrans->pdls[6] = psiq;
375 0 0         PDL_RETERROR(PDL_err, PDL->type_coerce(__privtrans));
376 0 0         PDL_RETERROR(PDL_err, PDL->make_trans_mutual(__privtrans));
377 0           return PDL_err;
378             }