Coverage for tests / test_truncatedGaussian.py: 16%
141 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-06 08:47 +0000
1#
2# LSST Data Management System
3#
4# Copyright 2008-2016 AURA/LSST.
5#
6# This product includes software developed by the
7# LSST Project (http://www.lsst.org/).
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the LSST License Statement and
20# the GNU General Public License along with this program. If not,
21# see <https://www.lsstcorp.org/LegalNotices/>.
22#
23import unittest
24import numpy
25try:
26 import scipy.integrate
27 import scipy.stats
28 import scipy.special
29except ImportError:
30 scipy = None
32import lsst.utils.logging
33import lsst.utils.tests
34import lsst.meas.modelfit
37if False:
38 lsst.utils.logging.trace_set_at("lsst.meas.modelfit.integrals", 5)
39 lsst.utils.logging.trace_set_at("lsst.meas.modelfit.TruncatedGaussian", 5)
42class TruncatedGaussianTestCase(lsst.utils.tests.TestCase):
44 def setUp(self):
45 numpy.random.seed(500)
47 def check1d(self, mu, hessian, tg):
48 evaluator = tg.evaluate()
49 logEvaluator = tg.evaluateLog()
50 dist = scipy.stats.norm(loc=mu[0], scale=hessian[0, 0]**-0.5)
51 self.assertFloatsAlmostEqual(1.0 - dist.cdf(0.0), tg.getUntruncatedFraction())
52 eps = 1E-7
53 if numpy.all(mu >= 0.0):
54 self.assertFloatsAlmostEqual(logEvaluator(mu), tg.getLogPeakAmplitude())
55 self.assertGreater(logEvaluator(mu+eps), tg.getLogPeakAmplitude())
56 self.assertGreater(logEvaluator(mu-eps), tg.getLogPeakAmplitude())
57 peak = tg.maximize()
58 self.assertGreater(evaluator(peak), 0.0)
59 self.assertLess(evaluator(peak+eps), evaluator(peak))
60 self.assertLess(evaluator(peak-eps), evaluator(peak))
62 def altLogEval(x):
63 if numpy.any(x < 0):
64 return float("inf")
65 return tg.getLogPeakAmplitude() + 0.5*hessian[0, 0]*(x-mu[0])**2
66 for alpha in (numpy.random.randn(10, 1) * hessian[0, 0]**-0.5 + mu[0]):
67 x1 = logEvaluator(alpha)
68 x2 = altLogEval(alpha[0])
69 if numpy.isfinite(x1) and numpy.isfinite(x2):
70 self.assertFloatsAlmostEqual(x1, x2, rtol=1E-14)
71 else:
72 self.assertEqual(x1, x2)
73 integral, check = self.integrate1d(tg)
74 self.assertLess(check, 1E-7)
75 self.assertFloatsAlmostEqual(integral, numpy.exp(-tg.getLogIntegral()), atol=check)
77 def check2d(self, mu, hessian, tg, isDegenerate=False):
78 evaluator = tg.evaluate()
79 logEvaluator = tg.evaluateLog()
80 unit1 = numpy.array([1.0, 0.0])
81 unit2 = numpy.array([0.0, 1.0])
82 eps = 1E-7
83 if numpy.all(mu >= 0.0):
84 self.assertFloatsAlmostEqual(logEvaluator(mu), tg.getLogPeakAmplitude())
85 self.assertGreater(logEvaluator(mu + unit1*eps), tg.getLogPeakAmplitude())
86 self.assertGreater(logEvaluator(mu - unit1*eps), tg.getLogPeakAmplitude())
87 self.assertGreater(logEvaluator(mu + unit2*eps), tg.getLogPeakAmplitude())
88 self.assertGreater(logEvaluator(mu - unit2*eps), tg.getLogPeakAmplitude())
89 peak = tg.maximize()
90 self.assertGreater(evaluator(peak), 0.0)
91 self.assertLess(evaluator(peak + unit1*eps) / evaluator(peak), 1.0)
92 self.assertLess(evaluator(peak - unit1*eps) / evaluator(peak), 1.0)
93 self.assertLess(evaluator(peak + unit2*eps) / evaluator(peak), 1.0)
94 self.assertLess(evaluator(peak - unit2*eps) / evaluator(peak), 1.0)
96 def altLogEval(a):
97 if numpy.any(a < 0):
98 return float("inf")
99 return tg.getLogPeakAmplitude() + 0.5*numpy.dot(numpy.dot(hessian, a - mu).transpose(), a - mu)
100 for alpha in (numpy.random.randn(10, 2) * hessian.diagonal()**-0.5 + mu):
101 x1 = logEvaluator(alpha)
102 x2 = altLogEval(alpha)
103 if numpy.isfinite(x1) and numpy.isfinite(x2):
104 self.assertFloatsAlmostEqual(x1, x2, rtol=1E-14)
105 else:
106 self.assertEqual(x1, x2)
107 integral, check = self.integrate2d(tg)
108 self.assertLess(check, 1E-7)
109 self.assertFloatsAlmostEqual(integral, numpy.exp(-tg.getLogIntegral()), atol=check)
111 def integrate1d(self, tg):
112 evaluator = tg.evaluate()
114 def func(x):
115 return evaluator(numpy.array([x]))
116 return scipy.integrate.quad(func, 0.0, numpy.inf)
118 def integrate2d(self, tg):
119 evaluator = tg.evaluate()
121 def func(x, y):
122 return evaluator(numpy.array([x, y]))
123 return scipy.integrate.dblquad(func, 0.0, numpy.inf, lambda x: 0.0, lambda x: numpy.inf)
125 @unittest.skipIf(scipy is None, "Test requires SciPy")
126 def test1d(self):
127 for i in range(5):
128 sigma = (numpy.random.randn(1, 1)**2 + 1)*5
129 mu = (numpy.random.randn(1))*3
130 q0 = float(numpy.random.randn())
131 hessian = numpy.linalg.inv(sigma)
132 gradient = -numpy.dot(hessian, mu)
133 tg1 = lsst.meas.modelfit.TruncatedGaussian.fromStandardParameters(mu, sigma)
134 tg2 = lsst.meas.modelfit.TruncatedGaussian.fromSeriesParameters(q0, gradient, hessian)
135 self.assertEqual(tg1.getLogIntegral(), 0.0)
136 self.assertFloatsAlmostEqual(tg1.getLogPeakAmplitude(),
137 (0.5*numpy.log(numpy.linalg.det(2.0*numpy.pi*sigma))
138 + numpy.log(tg1.getUntruncatedFraction())),
139 rtol=1E-13)
140 self.assertFloatsAlmostEqual(tg2.getLogPeakAmplitude(),
141 q0 + 0.5*numpy.dot(mu, gradient),
142 rtol=1E-13)
143 self.check1d(mu, hessian, tg1)
144 self.check1d(mu, hessian, tg2)
146 @unittest.skipIf(scipy is None, "Test requires SciPy")
147 def test2d(self):
148 for i in range(5):
149 x = numpy.linspace(-1, 1, 5)
150 model = numpy.zeros((x.size, 2), dtype=float)
151 model[:, 0] = x
152 model[:, 1] = x**2 + x
153 data = numpy.random.randn(x.size) + model[:, 0]*0.9 + model[:, 1]*1.1
154 q0 = 0.5*float(numpy.dot(data, data))
155 gradient = -numpy.dot(model.transpose(), data)
156 hessian = numpy.dot(model.transpose(), model)
157 sigma = numpy.linalg.inv(hessian)
158 self.assertFloatsAlmostEqual(numpy.linalg.inv(sigma), hessian, rtol=1E-15, atol=1E-15)
159 mu = -numpy.dot(sigma, gradient)
160 tg1 = lsst.meas.modelfit.TruncatedGaussian.fromStandardParameters(mu, sigma)
161 self.assertFloatsAlmostEqual(tg1.getLogPeakAmplitude(),
162 (numpy.log(tg1.getUntruncatedFraction())
163 + 0.5*numpy.log(numpy.linalg.det(2.0*numpy.pi*sigma))),
164 rtol=1E-13)
165 self.assertEqual(tg1.getLogIntegral(), 0.0)
166 self.check2d(mu, hessian, tg1)
167 tg2 = lsst.meas.modelfit.TruncatedGaussian.fromSeriesParameters(q0, gradient, hessian)
168 self.assertFloatsAlmostEqual(tg2.getLogPeakAmplitude(),
169 q0+0.5*numpy.dot(mu, gradient),
170 rtol=1E-13)
171 self.check2d(mu, hessian, tg2)
173 @unittest.skipIf(scipy is None, "Test requires SciPy")
174 def testDegenerate(self):
175 for i in range(5):
176 x = numpy.linspace(-1, 1, 5)
177 model = numpy.zeros((x.size, 2), dtype=float)
178 model[:, 0] = x
179 model[:, 1] = 2*x
180 data = numpy.random.randn(x.size) + model[:, 0]*0.9 + model[:, 1]*1.1
181 q0 = 0.5*float(numpy.dot(data, data))
182 gradient = -numpy.dot(model.transpose(), data)
183 hessian = numpy.dot(model.transpose(), model)
184 mu, _, _, _ = numpy.linalg.lstsq(model, data, rcond=None)
185 tg = lsst.meas.modelfit.TruncatedGaussian.fromSeriesParameters(q0, gradient, hessian)
186 self.assertFloatsAlmostEqual(tg.getLogPeakAmplitude(),
187 q0+0.5*numpy.dot(mu, gradient),
188 rtol=1E-13)
189 self.check2d(mu, hessian, tg, isDegenerate=True)
192class TestMemory(lsst.utils.tests.MemoryTestCase):
193 pass
196def setup_module(module):
197 lsst.utils.tests.init()
200if __name__ == "__main__": 200 ↛ 201line 200 didn't jump to line 201 because the condition on line 200 was never true
201 lsst.utils.tests.init()
202 unittest.main()