File Coverage

pp-logdiff.c

Criterion	Covered	Total	%
statement	34	39	87.1
branch	51	124	41.1
condition			n/a
subroutine			n/a
pod			n/a
total	85	163	52.1

line	stmt	bran	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-logdiff.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	7		static inline double logdiff1(double x, double y) {
56	7	100	if (y-x > LOG_BIG) { return y; }
57	6	100	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	5	100	else if (x>y) { return x + log( 1 - exp(y-x)); }
61	2		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_logdiff_readdata(pdl_trans *__privtrans) {
78			pdl_error PDL_err = {0, NULL, 0};
79			#line 80 "pp-logdiff.c"
80	7	50	if (!__privtrans->broadcast.incs) return PDL->make_error(PDL_EUSERERROR, "Error in logdiff:" "broadcast.incs NULL");
81			/* broadcastloop declarations */
82			int __brcloopval;
83			register PDL_Indx __tind0,__tind1; /* counters along dim */
84	7		register PDL_Indx __tnpdls = __privtrans->broadcast.npdls;
85			/* dims here are how many steps along those dims */
86	7		register PDL_Indx __tinc0_a = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,0);
87	7		register PDL_Indx __tinc0_b = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,0);
88	7		register PDL_Indx __tinc0_c = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,0);
89	7		register PDL_Indx __tinc1_a = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,1);
90	7		register PDL_Indx __tinc1_b = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,1);
91	7		register PDL_Indx __tinc1_c = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,1);
92			#define PDL_BROADCASTLOOP_START_logdiff_readdata PDL_BROADCASTLOOP_START( \
93			readdata, \
94			__privtrans->broadcast, \
95			__privtrans->vtable, \
96			a_datap += __offsp[0]; \
97			b_datap += __offsp[1]; \
98			c_datap += __offsp[2]; \
99			, \
100			( ,a_datap += __tinc1_a - __tinc0_a * __tdims0 \
101			,b_datap += __tinc1_b - __tinc0_b * __tdims0 \
102			,c_datap += __tinc1_c - __tinc0_c * __tdims0 \
103			), \
104			( ,a_datap += __tinc0_a \
105			,b_datap += __tinc0_b \
106			,c_datap += __tinc0_c \
107			) \
108			)
109			#define PDL_BROADCASTLOOP_END_logdiff_readdata PDL_BROADCASTLOOP_END( \
110			__privtrans->broadcast, \
111			a_datap -= __tinc1_a * __tdims1 + __offsp[0]; \
112			b_datap -= __tinc1_b * __tdims1 + __offsp[1]; \
113			c_datap -= __tinc1_c * __tdims1 + __offsp[2]; \
114			)
115			#ifndef PDL_DECLARE_PARAMS_logdiff_1
116			#define PDL_DECLARE_PARAMS_logdiff_1(PDL_TYPE_OP,PDL_PPSYM_OP) \
117			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, a, (__privtrans->pdls[0]), 1, PDL_PPSYM_OP) \
118			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, b, (__privtrans->pdls[1]), 1, PDL_PPSYM_OP) \
119			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, c, (__privtrans->pdls[2]), 1, PDL_PPSYM_OP)
120			#endif
121			#define PDL_IF_BAD(t,f) f
122	7		switch (__privtrans->__datatype) { /* Start generic switch */
123	0		case PDL_F: {
124	0	0	PDL_DECLARE_PARAMS_logdiff_1(PDL_Float,F)
		0
		0
		0
		0
		0
		0
		0
		0
125	0	0	PDL_BROADCASTLOOP_START_logdiff_readdata {(c_datap)[0] = logdiff((a_datap)[0],(b_datap)[0]);}PDL_BROADCASTLOOP_END_logdiff_readdata
		0
		0
		0
		0
		0
		0
		0
		0
126	0		} break;
127	6		case PDL_D: {
128	6	50	PDL_DECLARE_PARAMS_logdiff_1(PDL_Double,D)
		50
		50
		50
		50
		50
		50
		50
		50
129	24	50	PDL_BROADCASTLOOP_START_logdiff_readdata {(c_datap)[0] = logdiff((a_datap)[0],(b_datap)[0]);}PDL_BROADCASTLOOP_END_logdiff_readdata
		50
		50
		50
		50
		100
		100
		50
		50
130	6		} break;
131	1		case PDL_LD: {
132	1	50	PDL_DECLARE_PARAMS_logdiff_1(PDL_LDouble,E)
		50
		50
		50
		50
		50
		50
		50
		50
133	4	50	PDL_BROADCASTLOOP_START_logdiff_readdata {(c_datap)[0] = logdiff((a_datap)[0],(b_datap)[0]);}PDL_BROADCASTLOOP_END_logdiff_readdata
		50
		50
		50
		50
		100
		100
		50
		50
134	1		} break;
135	0		default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in logdiff: unhandled datatype(%d), only handles (FDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
136			}
137			#undef PDL_IF_BAD
138	7		return PDL_err;
139			}
140
141			static pdl_datatypes pdl_logdiff_vtable_gentypes[] = { PDL_F, PDL_D, PDL_LD, -1 };
142			static PDL_Indx pdl_logdiff_vtable_realdims[] = { 0, 0, 0 };
143			static char *pdl_logdiff_vtable_parnames[] = { "a","b","c" };
144			static short pdl_logdiff_vtable_parflags[] = {
145			0,
146			0,
147			PDL_PARAM_ISCREAT\|PDL_PARAM_ISOUT\|PDL_PARAM_ISWRITE
148			};
149			static pdl_datatypes pdl_logdiff_vtable_partypes[] = { -1, -1, -1 };
150			static PDL_Indx pdl_logdiff_vtable_realdims_starts[] = { 0, 0, 0 };
151			static PDL_Indx pdl_logdiff_vtable_realdims_ind_ids[] = { 0 };
152			static char *pdl_logdiff_vtable_indnames[] = { "" };
153			pdl_transvtable pdl_logdiff_vtable = {
154			PDL_TRANS_DO_BROADCAST, 0, pdl_logdiff_vtable_gentypes, 2, 3, NULL /CORE21/,
155			pdl_logdiff_vtable_realdims, pdl_logdiff_vtable_parnames,
156			pdl_logdiff_vtable_parflags, pdl_logdiff_vtable_partypes,
157			pdl_logdiff_vtable_realdims_starts, pdl_logdiff_vtable_realdims_ind_ids, 0,
158			0, pdl_logdiff_vtable_indnames,
159			NULL, pdl_logdiff_readdata, NULL,
160			NULL,
161			0,"PDL::HMM::logdiff"
162			};
163
164
165	7		pdl_error pdl_run_logdiff(pdl a,pdl b,pdl *c) {
166	7		pdl_error PDL_err = {0, NULL, 0};
167	7	50	if (!PDL) return (pdl_error){PDL_EFATAL, "PDL core struct is NULL, can't continue",0};
168	7		pdl_trans *__privtrans = PDL->create_trans(&pdl_logdiff_vtable);
169	7	50	if (!__privtrans) return PDL->make_error_simple(PDL_EFATAL, "Couldn't create trans");
170	7		__privtrans->pdls[0] = a;
171	7		__privtrans->pdls[1] = b;
172	7		__privtrans->pdls[2] = c;
173	7	50	PDL_RETERROR(PDL_err, PDL->type_coerce(__privtrans));
174	7	50	PDL_RETERROR(PDL_err, PDL->make_trans_mutual(__privtrans));
175	7		return PDL_err;
176			}