1 /**************************************************************
2 *** RHmm package
3 ***
4 *** File: LogProb.cpp
5 ***
6 *** Author: Ollivier TARAMASCO <Ollivier.Taramasco@imag.fr>
7 *** Author: Sebastian BAUER <sebastian.bauer@charite.de>
8 ***
9 **************************************************************/
10
11 #include "LogProb.h"
12
eexp(const double theX)13 double eexp(const double theX)
14 {
15 if (theX <= LOGZERO)
16 return(0.0) ;
17 else
18 return(exp(theX)) ;
19 }
20
eln(const double theX)21 double eln(const double theX)
22 {
23 if (theX > 0.0)
24 return(log(theX)) ;
25 else
26 return(LOGZERO) ;
27 }
28
elnsum1(const double theX,const double theY)29 double elnsum1(const double theX, const double theY)
30 {
31 double myeLnX = eln(theX),
32 myeLnY = eln(theY) ;
33
34 if ( (myeLnX <= LOGZERO) || (myeLnY <= LOGZERO) )
35 { if (myeLnX <= LOGZERO)
36 return(myeLnY) ;
37 else
38 return(myeLnX) ;
39 }
40 else
41 { if (myeLnX > myeLnY)
42 return(myeLnX + eln(1.0+exp(myeLnY-myeLnX))) ;
43 else
44 return(myeLnY + eln(1.0+exp(myeLnX-myeLnY))) ;
45 }
46 }
47
elnsum(const double theeLnX,const double theeLnY)48 double elnsum(const double theeLnX, const double theeLnY)
49 {
50 // elnsum(eln(x), eln(y)) = eln(x+y) pour x, y > LOGZERO
51 // elnsum(LOGZERO, eln(y)) = eln(y)
52 // elnsum(eln(x), LOGZERO) = eln(x)
53 double myeLnX = MAX(theeLnX, theeLnY),
54 myeLnY = MIN(theeLnX, theeLnY) ;
55
56 if (myeLnY <= LOGZERO)
57 return(myeLnX) ;
58 else
59 return(myeLnX + eln(1.0+exp(myeLnY-myeLnX))) ;
60 }
61
62
elnproduct1(const double theX,const double theY)63 double elnproduct1(const double theX, const double theY)
64 {
65 double myeLnX = eln(theX),
66 myeLnY = eln(theY) ;
67
68 if ( (myeLnX <= LOGZERO) || (myeLnY <= LOGZERO) )
69 return(LOGZERO) ;
70 else
71 return(myeLnX + myeLnY) ;
72 }
73
elnproduct(const double theeLnX,const double theeLnY)74 double elnproduct(const double theeLnX, const double theeLnY)
75 // elnproduct(eln(x), eln(y)) = eln(x) + eln(y) pour x, y > 0
76 // elnproduct(LOGZERO, eln(y)) = elnproduct(eln(x), LOGZERO) = LOGZERO
77 {
78 if ( (theeLnX <= LOGZERO) || (theeLnY <= LOGZERO) )
79 return(LOGZERO) ;
80 else
81 return(theeLnX + theeLnY) ;
82 }
83