File Coverage

pp-hmmviterbi.c
Criterion Covered Total %
statement 58 101 57.4
branch 44 202 21.7
condition n/a
subroutine n/a
pod n/a
total 102 303 33.6


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-hmmviterbi.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_hmmviterbi_readdata(pdl_trans *__privtrans) {
78             pdl_error PDL_err = {0, NULL, 0};
79             #line 80 "pp-hmmviterbi.c"
80 2           register PDL_Indx __N_size = __privtrans->ind_sizes[1];
81 2 50         if (!__privtrans->broadcast.incs) return PDL->make_error(PDL_EUSERERROR, "Error in hmmviterbi:" "broadcast.incs NULL");
82             /* broadcastloop declarations */
83             int __brcloopval;
84             register PDL_Indx __tind0,__tind1; /* counters along dim */
85 2           register PDL_Indx __tnpdls = __privtrans->broadcast.npdls;
86             /* dims here are how many steps along those dims */
87 2           register PDL_Indx __tinc0_a = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,0);
88 2           register PDL_Indx __tinc0_b = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,0);
89 2           register PDL_Indx __tinc0_pi = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,0);
90 2           register PDL_Indx __tinc0_o = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,0);
91 2           register PDL_Indx __tinc0_delta = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,0);
92 2           register PDL_Indx __tinc0_psi = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,0);
93 2           register PDL_Indx __tinc1_a = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,1);
94 2           register PDL_Indx __tinc1_b = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,1);
95 2           register PDL_Indx __tinc1_pi = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,1);
96 2           register PDL_Indx __tinc1_o = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,1);
97 2           register PDL_Indx __tinc1_delta = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,1);
98 2           register PDL_Indx __tinc1_psi = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,1);
99             #define PDL_BROADCASTLOOP_START_hmmviterbi_readdata PDL_BROADCASTLOOP_START( \
100             readdata, \
101             __privtrans->broadcast, \
102             __privtrans->vtable, \
103             a_datap += __offsp[0]; \
104             b_datap += __offsp[1]; \
105             pi_datap += __offsp[2]; \
106             o_datap += __offsp[3]; \
107             delta_datap += __offsp[4]; \
108             psi_datap += __offsp[5]; \
109             , \
110             ( ,a_datap += __tinc1_a - __tinc0_a * __tdims0 \
111             ,b_datap += __tinc1_b - __tinc0_b * __tdims0 \
112             ,pi_datap += __tinc1_pi - __tinc0_pi * __tdims0 \
113             ,o_datap += __tinc1_o - __tinc0_o * __tdims0 \
114             ,delta_datap += __tinc1_delta - __tinc0_delta * __tdims0 \
115             ,psi_datap += __tinc1_psi - __tinc0_psi * __tdims0 \
116             ), \
117             ( ,a_datap += __tinc0_a \
118             ,b_datap += __tinc0_b \
119             ,pi_datap += __tinc0_pi \
120             ,o_datap += __tinc0_o \
121             ,delta_datap += __tinc0_delta \
122             ,psi_datap += __tinc0_psi \
123             ) \
124             )
125             #define PDL_BROADCASTLOOP_END_hmmviterbi_readdata PDL_BROADCASTLOOP_END( \
126             __privtrans->broadcast, \
127             a_datap -= __tinc1_a * __tdims1 + __offsp[0]; \
128             b_datap -= __tinc1_b * __tdims1 + __offsp[1]; \
129             pi_datap -= __tinc1_pi * __tdims1 + __offsp[2]; \
130             o_datap -= __tinc1_o * __tdims1 + __offsp[3]; \
131             delta_datap -= __tinc1_delta * __tdims1 + __offsp[4]; \
132             psi_datap -= __tinc1_psi * __tdims1 + __offsp[5]; \
133             )
134 2           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;
135 2           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 2           register PDL_Indx __inc_delta_N = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,4,0)]; (void)__inc_delta_N;register PDL_Indx __inc_delta_T = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,4,1)]; (void)__inc_delta_T;
137 2           register PDL_Indx __inc_o_T = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,3,0)]; (void)__inc_o_T;
138 2           register PDL_Indx __inc_pi_N = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,2,0)]; (void)__inc_pi_N;
139 2           register PDL_Indx __inc_psi_N = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,5,0)]; (void)__inc_psi_N;register PDL_Indx __inc_psi_T = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,5,1)]; (void)__inc_psi_T;
140             #ifndef PDL_DECLARE_PARAMS_hmmviterbi_1
141             #define PDL_DECLARE_PARAMS_hmmviterbi_1(PDL_TYPE_OP,PDL_PPSYM_OP,PDL_TYPE_PARAM_psi,PDL_PPSYM_PARAM_psi) \
142             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, a, (__privtrans->pdls[0]), 1, PDL_PPSYM_OP) \
143             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, b, (__privtrans->pdls[1]), 1, PDL_PPSYM_OP) \
144             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, pi, (__privtrans->pdls[2]), 1, PDL_PPSYM_OP) \
145             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, o, (__privtrans->pdls[3]), 1, PDL_PPSYM_OP) \
146             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, delta, (__privtrans->pdls[4]), 1, PDL_PPSYM_OP) \
147             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_psi, psi, (__privtrans->pdls[5]), 1, PDL_PPSYM_PARAM_psi)
148             #endif
149             #define PDL_IF_BAD(t,f) f
150 2           switch (__privtrans->__datatype) { /* Start generic switch */
151 0           case PDL_F: {
152 0 0         PDL_DECLARE_PARAMS_hmmviterbi_1(PDL_Float,F,PDL_Long,L)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
153 0 0         PDL_BROADCASTLOOP_START_hmmviterbi_readdata {
    0          
    0          
    0          
    0          
    0          
    0          
154             int i,j, t, o_t;
155             double delta_jt, delta_tmp;
156             int psi_jt;
157              
158             /*-- Initialize: t==0: Loop: state_0==N --*/
159 0           o_t = (o_datap)[0+(__inc_o_T*(0))];
160 0 0         {/* Open N */ PDL_EXPAND2(register PDL_Indx N=0, __N_stop=(__N_size)); for(; N<__N_stop; N+=1) {
161 0           (delta_datap)[0+(__inc_delta_N*(N))+(__inc_delta_T*(0))] = (pi_datap)[0+(__inc_pi_N*(N))] + (b_datap)[0+(__inc_b_N*(N))+(__inc_b_M*(o_t))];
162 0           (psi_datap)[0+(__inc_psi_N*(N))+(__inc_psi_T*(0))] = 0;
163             #ifdef DEBUG_VITERBI
164             printf("t=0,j=%d,o_t=%d: delta(t=0,j=%d)=%.2e psi(t=0,j=%d)=%.0g b(j=%d,o=%d)=%.2e\n",
165             N,o_t,
166             N, exp((delta_datap)[0+(__inc_delta_N*(N))+(__inc_delta_T*(0))]),
167             N, (psi_datap)[0+(__inc_psi_N*(N))+(__inc_psi_T*(0))],
168             N,o_t, exp((b_datap)[0+(__inc_b_N*(N))+(__inc_b_M*(o_t))]));
169             #endif
170             }} /* Close N */
171              
172             #ifdef DEBUG_VITERBI
173             printf("\n");
174             #endif
175              
176             /*-- Main: t>0: Loop: time==t --*/
177 0 0         for (t=1; t<__privtrans->ind_sizes[2]; t++) {
178 0           o_t = (o_datap)[0+(__inc_o_T*(t))];
179              
180             /*-- Main: t>0: Loop: state_t==j --*/
181 0 0         for (j=0; j<__privtrans->ind_sizes[1]; j++) {
182 0           psi_jt = 0;
183 0           delta_jt = (delta_datap)[0+(__inc_delta_N*(0))+(__inc_delta_T*(t-1))] + (a_datap)[0+(__inc_a_N0*(0))+(__inc_a_N1*(j))];
184              
185             /*-- Main: t>0: Loop: state_(t-1)==i --*/
186 0 0         for (i=1; i<__privtrans->ind_sizes[1]; i++) {
187 0           delta_tmp = (delta_datap)[0+(__inc_delta_N*(i))+(__inc_delta_T*(t-1))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))];
188              
189 0 0         if (delta_tmp > delta_jt) {
190 0           delta_jt = delta_tmp;
191 0           psi_jt = i;
192             #ifdef DEBUG_VITERBI
193             printf("+");
194             #endif
195             }
196              
197             #ifdef DEBUG_VITERBI
198             printf("t=%d,i=%d,j=%d,o_t=%d: deltaX(i=%d,t=%d)=%.2e psi(j=%d,t=%d)=%.0g delta(j=%d,t=%d)=%.2e\n",
199             t,i,j,o_t,
200             i,t, exp(delta_tmp),
201             j,t, psi_jt,
202             j,t, exp(delta_jt));
203             #endif
204             }
205              
206             /*-- Main: t>0: Store data for state,time=(j,t) --*/
207 0           (delta_datap)[0+(__inc_delta_N*(j))+(__inc_delta_T*(t))] = delta_jt + (b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_t))];
208 0           (psi_datap)[0+(__inc_psi_N*(j))+(__inc_psi_T*(t))] = psi_jt;
209              
210             #ifdef DEBUG_VITERBI
211             printf("\n---> t=%d: b(j=%d,o=%d)=%.2e delta(j=%d,t=%d)=%.2e psi(j=%d,t=%d)=%.0g\n\n",
212             t, j,o_t, exp((b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_t))]),
213             j,t, exp((delta_datap)[0+(__inc_delta_N*(j))+(__inc_delta_T*(t))]),
214             j,t, (psi_datap)[0+(__inc_psi_N*(j))+(__inc_psi_T*(t))]);
215             #endif
216             }
217             #ifdef DEBUG_VITERBI
218             printf("\n");
219             #endif
220             }
221 0 0         }PDL_BROADCASTLOOP_END_hmmviterbi_readdata
    0          
222 0           } break;
223 2           case PDL_D: {
224 2 50         PDL_DECLARE_PARAMS_hmmviterbi_1(PDL_Double,D,PDL_Long,L)
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
225 8 50         PDL_BROADCASTLOOP_START_hmmviterbi_readdata {
    50          
    50          
    50          
    50          
    100          
    100          
226             int i,j, t, o_t;
227             double delta_jt, delta_tmp;
228             int psi_jt;
229              
230             /*-- Initialize: t==0: Loop: state_0==N --*/
231 2           o_t = (o_datap)[0+(__inc_o_T*(0))];
232 6 100         {/* Open N */ PDL_EXPAND2(register PDL_Indx N=0, __N_stop=(__N_size)); for(; N<__N_stop; N+=1) {
233 4           (delta_datap)[0+(__inc_delta_N*(N))+(__inc_delta_T*(0))] = (pi_datap)[0+(__inc_pi_N*(N))] + (b_datap)[0+(__inc_b_N*(N))+(__inc_b_M*(o_t))];
234 4           (psi_datap)[0+(__inc_psi_N*(N))+(__inc_psi_T*(0))] = 0;
235             #ifdef DEBUG_VITERBI
236             printf("t=0,j=%d,o_t=%d: delta(t=0,j=%d)=%.2e psi(t=0,j=%d)=%.0g b(j=%d,o=%d)=%.2e\n",
237             N,o_t,
238             N, exp((delta_datap)[0+(__inc_delta_N*(N))+(__inc_delta_T*(0))]),
239             N, (psi_datap)[0+(__inc_psi_N*(N))+(__inc_psi_T*(0))],
240             N,o_t, exp((b_datap)[0+(__inc_b_N*(N))+(__inc_b_M*(o_t))]));
241             #endif
242             }} /* Close N */
243              
244             #ifdef DEBUG_VITERBI
245             printf("\n");
246             #endif
247              
248             /*-- Main: t>0: Loop: time==t --*/
249 5 100         for (t=1; t<__privtrans->ind_sizes[2]; t++) {
250 3           o_t = (o_datap)[0+(__inc_o_T*(t))];
251              
252             /*-- Main: t>0: Loop: state_t==j --*/
253 9 100         for (j=0; j<__privtrans->ind_sizes[1]; j++) {
254 6           psi_jt = 0;
255 6           delta_jt = (delta_datap)[0+(__inc_delta_N*(0))+(__inc_delta_T*(t-1))] + (a_datap)[0+(__inc_a_N0*(0))+(__inc_a_N1*(j))];
256              
257             /*-- Main: t>0: Loop: state_(t-1)==i --*/
258 12 100         for (i=1; i<__privtrans->ind_sizes[1]; i++) {
259 6           delta_tmp = (delta_datap)[0+(__inc_delta_N*(i))+(__inc_delta_T*(t-1))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))];
260              
261 6 100         if (delta_tmp > delta_jt) {
262 1           delta_jt = delta_tmp;
263 1           psi_jt = i;
264             #ifdef DEBUG_VITERBI
265             printf("+");
266             #endif
267             }
268              
269             #ifdef DEBUG_VITERBI
270             printf("t=%d,i=%d,j=%d,o_t=%d: deltaX(i=%d,t=%d)=%.2e psi(j=%d,t=%d)=%.0g delta(j=%d,t=%d)=%.2e\n",
271             t,i,j,o_t,
272             i,t, exp(delta_tmp),
273             j,t, psi_jt,
274             j,t, exp(delta_jt));
275             #endif
276             }
277              
278             /*-- Main: t>0: Store data for state,time=(j,t) --*/
279 6           (delta_datap)[0+(__inc_delta_N*(j))+(__inc_delta_T*(t))] = delta_jt + (b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_t))];
280 6           (psi_datap)[0+(__inc_psi_N*(j))+(__inc_psi_T*(t))] = psi_jt;
281              
282             #ifdef DEBUG_VITERBI
283             printf("\n---> t=%d: b(j=%d,o=%d)=%.2e delta(j=%d,t=%d)=%.2e psi(j=%d,t=%d)=%.0g\n\n",
284             t, j,o_t, exp((b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_t))]),
285             j,t, exp((delta_datap)[0+(__inc_delta_N*(j))+(__inc_delta_T*(t))]),
286             j,t, (psi_datap)[0+(__inc_psi_N*(j))+(__inc_psi_T*(t))]);
287             #endif
288             }
289             #ifdef DEBUG_VITERBI
290             printf("\n");
291             #endif
292             }
293 2 50         }PDL_BROADCASTLOOP_END_hmmviterbi_readdata
    50          
294 2           } break;
295 0           case PDL_LD: {
296 0 0         PDL_DECLARE_PARAMS_hmmviterbi_1(PDL_LDouble,E,PDL_Long,L)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
297 0 0         PDL_BROADCASTLOOP_START_hmmviterbi_readdata {
    0          
    0          
    0          
    0          
    0          
    0          
298             int i,j, t, o_t;
299             double delta_jt, delta_tmp;
300             int psi_jt;
301              
302             /*-- Initialize: t==0: Loop: state_0==N --*/
303 0           o_t = (o_datap)[0+(__inc_o_T*(0))];
304 0 0         {/* Open N */ PDL_EXPAND2(register PDL_Indx N=0, __N_stop=(__N_size)); for(; N<__N_stop; N+=1) {
305 0           (delta_datap)[0+(__inc_delta_N*(N))+(__inc_delta_T*(0))] = (pi_datap)[0+(__inc_pi_N*(N))] + (b_datap)[0+(__inc_b_N*(N))+(__inc_b_M*(o_t))];
306 0           (psi_datap)[0+(__inc_psi_N*(N))+(__inc_psi_T*(0))] = 0;
307             #ifdef DEBUG_VITERBI
308             printf("t=0,j=%d,o_t=%d: delta(t=0,j=%d)=%.2e psi(t=0,j=%d)=%.0g b(j=%d,o=%d)=%.2e\n",
309             N,o_t,
310             N, exp((delta_datap)[0+(__inc_delta_N*(N))+(__inc_delta_T*(0))]),
311             N, (psi_datap)[0+(__inc_psi_N*(N))+(__inc_psi_T*(0))],
312             N,o_t, exp((b_datap)[0+(__inc_b_N*(N))+(__inc_b_M*(o_t))]));
313             #endif
314             }} /* Close N */
315              
316             #ifdef DEBUG_VITERBI
317             printf("\n");
318             #endif
319              
320             /*-- Main: t>0: Loop: time==t --*/
321 0 0         for (t=1; t<__privtrans->ind_sizes[2]; t++) {
322 0           o_t = (o_datap)[0+(__inc_o_T*(t))];
323              
324             /*-- Main: t>0: Loop: state_t==j --*/
325 0 0         for (j=0; j<__privtrans->ind_sizes[1]; j++) {
326 0           psi_jt = 0;
327 0           delta_jt = (delta_datap)[0+(__inc_delta_N*(0))+(__inc_delta_T*(t-1))] + (a_datap)[0+(__inc_a_N0*(0))+(__inc_a_N1*(j))];
328              
329             /*-- Main: t>0: Loop: state_(t-1)==i --*/
330 0 0         for (i=1; i<__privtrans->ind_sizes[1]; i++) {
331 0           delta_tmp = (delta_datap)[0+(__inc_delta_N*(i))+(__inc_delta_T*(t-1))] + (a_datap)[0+(__inc_a_N0*(i))+(__inc_a_N1*(j))];
332              
333 0 0         if (delta_tmp > delta_jt) {
334 0           delta_jt = delta_tmp;
335 0           psi_jt = i;
336             #ifdef DEBUG_VITERBI
337             printf("+");
338             #endif
339             }
340              
341             #ifdef DEBUG_VITERBI
342             printf("t=%d,i=%d,j=%d,o_t=%d: deltaX(i=%d,t=%d)=%.2e psi(j=%d,t=%d)=%.0g delta(j=%d,t=%d)=%.2e\n",
343             t,i,j,o_t,
344             i,t, exp(delta_tmp),
345             j,t, psi_jt,
346             j,t, exp(delta_jt));
347             #endif
348             }
349              
350             /*-- Main: t>0: Store data for state,time=(j,t) --*/
351 0           (delta_datap)[0+(__inc_delta_N*(j))+(__inc_delta_T*(t))] = delta_jt + (b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_t))];
352 0           (psi_datap)[0+(__inc_psi_N*(j))+(__inc_psi_T*(t))] = psi_jt;
353              
354             #ifdef DEBUG_VITERBI
355             printf("\n---> t=%d: b(j=%d,o=%d)=%.2e delta(j=%d,t=%d)=%.2e psi(j=%d,t=%d)=%.0g\n\n",
356             t, j,o_t, exp((b_datap)[0+(__inc_b_N*(j))+(__inc_b_M*(o_t))]),
357             j,t, exp((delta_datap)[0+(__inc_delta_N*(j))+(__inc_delta_T*(t))]),
358             j,t, (psi_datap)[0+(__inc_psi_N*(j))+(__inc_psi_T*(t))]);
359             #endif
360             }
361             #ifdef DEBUG_VITERBI
362             printf("\n");
363             #endif
364             }
365 0 0         }PDL_BROADCASTLOOP_END_hmmviterbi_readdata
    0          
366 0           } break;
367 0           default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in hmmviterbi: unhandled datatype(%d), only handles (FDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
368             }
369             #undef PDL_IF_BAD
370 2           return PDL_err;
371             }
372              
373             static pdl_datatypes pdl_hmmviterbi_vtable_gentypes[] = { PDL_F, PDL_D, PDL_LD, -1 };
374             static PDL_Indx pdl_hmmviterbi_vtable_realdims[] = { 2, 2, 1, 1, 2, 2 };
375             static char *pdl_hmmviterbi_vtable_parnames[] = { "a","b","pi","o","delta","psi" };
376             static short pdl_hmmviterbi_vtable_parflags[] = {
377             0,
378             0,
379             0,
380             0,
381             PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISWRITE,
382             PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISTYPED|PDL_PARAM_ISWRITE
383             };
384             static pdl_datatypes pdl_hmmviterbi_vtable_partypes[] = { -1, -1, -1, -1, -1, PDL_L };
385             static PDL_Indx pdl_hmmviterbi_vtable_realdims_starts[] = { 0, 2, 4, 5, 6, 8 };
386             static PDL_Indx pdl_hmmviterbi_vtable_realdims_ind_ids[] = { 1, 1, 1, 0, 1, 2, 1, 2, 1, 2 };
387             static char *pdl_hmmviterbi_vtable_indnames[] = { "M","N","T" };
388             pdl_transvtable pdl_hmmviterbi_vtable = {
389             PDL_TRANS_DO_BROADCAST, 0, pdl_hmmviterbi_vtable_gentypes, 4, 6, NULL /*CORE21*/,
390             pdl_hmmviterbi_vtable_realdims, pdl_hmmviterbi_vtable_parnames,
391             pdl_hmmviterbi_vtable_parflags, pdl_hmmviterbi_vtable_partypes,
392             pdl_hmmviterbi_vtable_realdims_starts, pdl_hmmviterbi_vtable_realdims_ind_ids, 10,
393             3, pdl_hmmviterbi_vtable_indnames,
394             NULL, pdl_hmmviterbi_readdata, NULL,
395             NULL,
396             0,"PDL::HMM::hmmviterbi"
397             };
398              
399              
400 2           pdl_error pdl_run_hmmviterbi(pdl *a,pdl *b,pdl *pi,pdl *o,pdl *delta,pdl *psi) {
401 2           pdl_error PDL_err = {0, NULL, 0};
402 2 50         if (!PDL) return (pdl_error){PDL_EFATAL, "PDL core struct is NULL, can't continue",0};
403 2           pdl_trans *__privtrans = PDL->create_trans(&pdl_hmmviterbi_vtable);
404 2 50         if (!__privtrans) return PDL->make_error_simple(PDL_EFATAL, "Couldn't create trans");
405 2           __privtrans->pdls[0] = a;
406 2           __privtrans->pdls[1] = b;
407 2           __privtrans->pdls[2] = pi;
408 2           __privtrans->pdls[3] = o;
409 2           __privtrans->pdls[4] = delta;
410 2           __privtrans->pdls[5] = psi;
411 2 50         PDL_RETERROR(PDL_err, PDL->type_coerce(__privtrans));
412 2 50         PDL_RETERROR(PDL_err, PDL->make_trans_mutual(__privtrans));
413 2           return PDL_err;
414             }