LAL  7.1.7.1-ab2cc12
LALMarcumQTest.py
Go to the documentation of this file.
1 # Copyright (C) 2014 Kipp Cannon
2 #
3 # This program is free software; you can redistribute it and/or modify it
4 # under the terms of the GNU General Public License as published by the
5 # Free Software Foundation; either version 3 of the License, or (at your
6 # option) any later version.
7 #
8 # This program is distributed in the hope that it will be useful, but
9 # WITHOUT ANY WARRANTY; without even the implied warranty of
10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11 # Public License for more details.
12 #
13 # You should have received a copy of the GNU General Public License along
14 # with this program; if not, write to the Free Software Foundation, Inc.,
15 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
16 
17 
18 #
19 # =============================================================================
20 #
21 # Preamble
22 #
23 # =============================================================================
24 #
25 
26 
27 import random
28 import lal
29 
30 
31 #
32 # =============================================================================
33 #
34 # Test Data
35 #
36 # =============================================================================
37 #
38 
39 
40 #
41 # test data copied from testmarq.f90 in reference implementation by Gil et
42 # al. published as part of
43 #
44 # A. Gil, J. Segura, and N. M. Temme. Algorithm 939: Computation of the
45 # Marcum Q-Function. ACM Transactions on Mathematical Software (TOMS),
46 # Volume 40 Issue 3, April 2014, Article No. 20. arXiv:1311.0681
47 #
48 
49 mu = [1.0, 3.0, 4.0, 6.0, 8.0, 10.0, 20.0, 22.0, 25.0, 27.0, 30.0, 32.0, 40.0, 50.0, 200.0, 350.0, 570.0, 1000.0]
50 
51 x = [0.3, 2.0, 8.0, 25.0, 13.0, 45.0, 47.0, 100.0, 85.0, 120.0, 130.0, 140.0, 30.0, 40.0, 0.01, 100.0, 1.0, 0.08]
52 
53 y = [0.01, 0.1, 50.0, 10.0, 15.0, 25.0, 30.0, 150.0, 60.0, 205.0, 90.0, 100.0, 120.0, 150.0, 190.0, 320.0, 480.0, 799.0]
54 
55 p = [.7382308441994e-02, .2199222796783e-04, .9999999768807, .1746869995977e-03, .1483130042637, .1748328323235e-03, .1340769184710e-04, .9646591215441, .1783991673043e-04, .9994542406431, .1220231636641e-05, .1757487653307e-05, .9999894753719, .9999968347378, .2431297758708, .3851423018735e-09, .2984493152360e-04, .4191472999694e-11]
56 
57 q = [.9926176915580, .9999780077720, .2311934913546e-07, .9998253130004, .8516869957363, .9998251671677, .9999865923082, .3534087845586e-01, .9999821600833, .5457593568564e-03, .9999987797684, .9999982425123, .1052462813144e-04, .3165262228904e-05, .7568702241292, .9999999996149, .9999701550685, .9999999999958]
58 
59 for mu, x, y, p, q in zip(mu, x, y, p, q):
60  lal_q = lal.MarcumQmodified(mu, x, y)
61  rel_err_q = abs(q - lal_q) / q
62  print("Q(%g,%g,%g) = %.16g\t%.16g\t%.16g" % (mu, x, y, q, lal_q, rel_err_q))
63  assert rel_err_q < 1e-12
64 
65 
66 #
67 # =============================================================================
68 #
69 # Extra Tests
70 #
71 # =============================================================================
72 #
73 
74 
75 #
76 # exercise specific code paths
77 #
78 
79 
80 # series expansion in section 3, x < 30
81 print()
82 print(lal.MarcumQmodified(1., 18, 8.))
83 
84 # asymptotic expansion in section 4.1, xsi > 30 && M*M < 2 * xsi
85 # for x < y
86 print()
87 print(lal.MarcumQmodified(1., 32., 40.5))
88 # for x == y
89 print(lal.MarcumQmodified(1., 32., 32.))
90 # for y > x
91 print(lal.MarcumQmodified(1., 32., 8.))
92 
93 # recurrence relation in (14), f1 < y < f2 && M < 135.
94 print()
95 print(lal.MarcumQmodified(50., 32., 80.))
96 
97 # asymptotic expansion in section 4.2, f1 < y < f2 && M >= 135.
98 #print
99 #print lal.MarcumQmodified(150., 32., 180.)
100 
101 # quadrature integration in section 5
102 print()
103 print(lal.MarcumQmodified(50., 32., 98.))
104 
105 # more tests
106 print()
107 print(lal.MarcumQmodified(1., 31.99999, 32.00001), lal.MarcumQmodified(1., 32.00001, 31.99999))
108 print(lal.MarcumQmodified(1., 31.9999999999, 32.0000000001), lal.MarcumQmodified(1., 32.0000000001, 31.9999999999))
109 print(lal.MarcumQmodified(1., 32., 32.), lal.MarcumQmodified(1., 32., 32.))
110 
111 
112 #
113 # =============================================================================
114 #
115 # Random Input
116 #
117 # =============================================================================
118 #
119 
120 
121 for i in range(100000):
122  M, x, y = random.uniform(1, 10000), random.uniform(0, 10000), random.uniform(0, 10000)
123  if x + M - (4.*x + 2.*M)**.5 < y < x + M + (4.*x + 2.*M)**.5 and M >= 135:
124  continue
125  print("trying XLALMarcumQmodified(%.17g, %.17g, %.17g)" % (M, x, y))
126  lal.MarcumQmodified(M, x, y)