diff options
| -rw-r--r-- | src/likelihood.c | 4 | ||||
| -rw-r--r-- | src/misc.c | 119 | ||||
| -rw-r--r-- | src/misc.h | 2 | ||||
| -rw-r--r-- | src/test.c | 23 |
4 files changed, 146 insertions, 2 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 6ddd576..76364fd 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -83,7 +83,7 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m * different transit times across the face of the PMT, it seems better to * convolve first which is what we do here. In addition, the problem is not * analytically tractable if you do things the other way around. */ - return log(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)); + return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)); } static double gsl_muon_time(double x, void *params) @@ -414,7 +414,7 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl logp[j] = log(pq(ev->pmt_hits[i].qhs,j)) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], 1, &ts[i], tmean, 1.5); if (j == 1 || logp[j] > max_logp) max_logp = logp[j]; - if (logp[j] - max_logp < MIN_RATIO*log(10)) { + if (logp[j] - max_logp < MIN_RATIO*ln(10)) { j++; break; } @@ -6,6 +6,113 @@ static struct { int n; double f; +} ln_table[LN_MAX + 1] = { + {0,-INFINITY}, + {1,0}, + {2,0.69314718055994529}, + {3,1.0986122886681098}, + {4,1.3862943611198906}, + {5,1.6094379124341003}, + {6,1.791759469228055}, + {7,1.9459101490553132}, + {8,2.0794415416798357}, + {9,2.1972245773362196}, + {10,2.3025850929940459}, + {11,2.3978952727983707}, + {12,2.4849066497880004}, + {13,2.5649493574615367}, + {14,2.6390573296152584}, + {15,2.7080502011022101}, + {16,2.7725887222397811}, + {17,2.8332133440562162}, + {18,2.8903717578961645}, + {19,2.9444389791664403}, + {20,2.9957322735539909}, + {21,3.044522437723423}, + {22,3.0910424533583161}, + {23,3.1354942159291497}, + {24,3.1780538303479458}, + {25,3.2188758248682006}, + {26,3.2580965380214821}, + {27,3.2958368660043291}, + {28,3.3322045101752038}, + {29,3.3672958299864741}, + {30,3.4011973816621555}, + {31,3.4339872044851463}, + {32,3.4657359027997265}, + {33,3.4965075614664802}, + {34,3.5263605246161616}, + {35,3.5553480614894135}, + {36,3.5835189384561099}, + {37,3.6109179126442243}, + {38,3.6375861597263857}, + {39,3.6635616461296463}, + {40,3.6888794541139363}, + {41,3.713572066704308}, + {42,3.7376696182833684}, + {43,3.7612001156935624}, + {44,3.784189633918261}, + {45,3.8066624897703196}, + {46,3.8286413964890951}, + {47,3.8501476017100584}, + {48,3.8712010109078911}, + {49,3.8918202981106265}, + {50,3.912023005428146}, + {51,3.9318256327243257}, + {52,3.9512437185814275}, + {53,3.970291913552122}, + {54,3.9889840465642745}, + {55,4.0073331852324712}, + {56,4.0253516907351496}, + {57,4.0430512678345503}, + {58,4.0604430105464191}, + {59,4.0775374439057197}, + {60,4.0943445622221004}, + {61,4.1108738641733114}, + {62,4.1271343850450917}, + {63,4.1431347263915326}, + {64,4.1588830833596715}, + {65,4.1743872698956368}, + {66,4.1896547420264252}, + {67,4.2046926193909657}, + {68,4.219507705176107}, + {69,4.2341065045972597}, + {70,4.2484952420493594}, + {71,4.2626798770413155}, + {72,4.2766661190160553}, + {73,4.290459441148391}, + {74,4.3040650932041702}, + {75,4.3174881135363101}, + {76,4.3307333402863311}, + {77,4.3438054218536841}, + {78,4.3567088266895917}, + {79,4.3694478524670215}, + {80,4.3820266346738812}, + {81,4.3944491546724391}, + {82,4.4067192472642533}, + {83,4.4188406077965983}, + {84,4.4308167988433134}, + {85,4.4426512564903167}, + {86,4.4543472962535073}, + {87,4.4659081186545837}, + {88,4.4773368144782069}, + {89,4.4886363697321396}, + {90,4.499809670330265}, + {91,4.5108595065168497}, + {92,4.5217885770490405}, + {93,4.5325994931532563}, + {94,4.5432947822700038}, + {95,4.5538768916005408}, + {96,4.5643481914678361}, + {97,4.5747109785033828}, + {98,4.5849674786705723}, + {99,4.5951198501345898}, + {100,4.6051701859880918}, +}; + +static struct { + int n; + double f; } ln_fact_table[LNFACT_MAX + 1] = { {0,0}, {1,0}, @@ -110,6 +217,18 @@ static struct { {100,363.73937555556347}, }; +double ln(unsigned int n) +{ + /* Returns the logarithm of n. + * + * Uses a lookup table to return results for n < 100. */ + + if (n <= LN_MAX) + return ln_table[n].f; + + return log(n); +} + double lnfact(unsigned int n) { /* Returns the logarithm of n!. @@ -3,8 +3,10 @@ #include <stdlib.h> /* for size_t */ +#define LN_MAX 100 #define LNFACT_MAX 100 +double ln(unsigned int n); double lnfact(unsigned int n); double kahan_sum(double *x, size_t n); double interp1d(double x, double *xp, double *yp, size_t n); @@ -489,6 +489,28 @@ err: return 1; } +int test_ln(char *err) +{ + /* Tests the lnfact() function. */ + size_t i; + double x, expected; + + for (i = 1; i < 1000; i++) { + x = ln(i); + expected = log(i); + + if (!isclose(x, expected, 1e-9, 1e-9)) { + sprintf(err, "ln(%zu) returned %.5g, but expected %.5g", i, x, expected); + goto err; + } + } + + return 0; + +err: + return 1; +} + struct tests { testFunction *test; char *name; @@ -505,6 +527,7 @@ struct tests { {test_interp1d, "test_interp1d"}, {test_kahan_sum, "test_kahan_sum"}, {test_lnfact, "test_lnfact"}, + {test_ln, "test_ln"}, }; int main(int argc, char **argv) |