Coverage for tests / test_softenedLinearPrior.py: 16%
127 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-01 08:35 +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
26import lsst.utils.tests
27import lsst.utils.logging
28import lsst.meas.modelfit
30try:
31 import scipy.integrate
32except ImportError:
33 scipy = None
35lsst.utils.logging.trace_set_at("lsst.meas.modelfit.SoftenedLinearPrior", 5)
38class SoftenedLinearPriorTestCase(lsst.utils.tests.TestCase):
40 NUM_DIFF_STEP = 1E-3
42 def setUp(self):
43 # a prior with broad ramps and non-zero slope; broad ramps makes evaluating numerical
44 # derivatives easier, and we want to do that to check the analytic ones
45 numpy.random.seed(500)
46 ctrl = lsst.meas.modelfit.SoftenedLinearPrior.Control()
47 ctrl.logRadiusMinOuter = ctrl.logRadiusMinInner - 2.0
48 ctrl.logRadiusMaxOuter = ctrl.logRadiusMaxInner + 2.0
49 ctrl.ellipticityMaxOuter = ctrl.ellipticityMaxInner + 2.0
50 ctrl.logRadiusMinMaxRatio = 2.0
51 self.prior = lsst.meas.modelfit.SoftenedLinearPrior(ctrl)
52 self.amplitudes = numpy.array([1.0], dtype=lsst.meas.modelfit.Scalar)
54 def tearDown(self):
55 del self.prior
56 del self.amplitudes
58 def evaluatePrior(self, e1, e2, r):
59 b = numpy.broadcast(e1, e2, r)
60 p = numpy.zeros(b.shape, dtype=lsst.meas.modelfit.Scalar)
61 for i, (e1i, e2i, ri) in enumerate(b):
62 p.flat[i] = self.prior.evaluate(numpy.array([e1i, e2i, ri]), self.amplitudes)
63 return p
65 def checkDerivatives(self, e1, e2, r):
66 nonlinear = numpy.array([e1, e2, r], dtype=lsst.meas.modelfit.Scalar)
67 amplitudeGradient = numpy.zeros(1, dtype=lsst.meas.modelfit.Scalar)
68 amplitudeHessian = numpy.zeros((1, 1), dtype=lsst.meas.modelfit.Scalar)
69 crossHessian = numpy.zeros((3, 1), dtype=lsst.meas.modelfit.Scalar)
70 nonlinearGradient = numpy.zeros(3, dtype=lsst.meas.modelfit.Scalar)
71 nonlinearHessian = numpy.zeros((3, 3), dtype=lsst.meas.modelfit.Scalar)
72 self.prior.evaluateDerivatives(nonlinear, self.amplitudes,
73 nonlinearGradient, amplitudeGradient,
74 nonlinearHessian, amplitudeHessian,
75 crossHessian)
76 p = self.prior.evaluate(nonlinear, self.amplitudes)
77 for i in range(3):
78 nonlinearA = nonlinear.copy()
79 nonlinearB = nonlinear.copy()
80 nonlinearA[i] -= self.NUM_DIFF_STEP
81 nonlinearB[i] += self.NUM_DIFF_STEP
82 pA = self.prior.evaluate(nonlinearA, self.amplitudes)
83 pB = self.prior.evaluate(nonlinearB, self.amplitudes)
84 dp = (pB - pA) / (2*self.NUM_DIFF_STEP)
85 self.assertFloatsAlmostEqual(nonlinearGradient[i], dp, rtol=1E-3, atol=1E-8)
86 d2p = (pA + pB - 2*p) / self.NUM_DIFF_STEP**2
87 self.assertFloatsAlmostEqual(nonlinearHessian[i, i], d2p, rtol=1E-3, atol=1E-8)
88 for j in range(i+1, 3):
89 nonlinearAA = nonlinearA.copy()
90 nonlinearAB = nonlinearA.copy()
91 nonlinearBA = nonlinearB.copy()
92 nonlinearBB = nonlinearB.copy()
93 nonlinearAA[j] -= self.NUM_DIFF_STEP
94 nonlinearAB[j] += self.NUM_DIFF_STEP
95 nonlinearBA[j] -= self.NUM_DIFF_STEP
96 nonlinearBB[j] += self.NUM_DIFF_STEP
97 pAA = self.prior.evaluate(nonlinearAA, self.amplitudes)
98 pAB = self.prior.evaluate(nonlinearAB, self.amplitudes)
99 pBA = self.prior.evaluate(nonlinearBA, self.amplitudes)
100 pBB = self.prior.evaluate(nonlinearBB, self.amplitudes)
101 d2p = (pBB - pAB - pBA + pAA) / (4*self.NUM_DIFF_STEP**2)
102 self.assertFloatsAlmostEqual(nonlinearHessian[i, j], d2p, rtol=1E-3, atol=1E-8)
104 def testDerivatives(self):
105 """Test that evaluateDerivatives() returns results similar to finite-differences
106 on evaluate().
107 """
108 ctrl = self.prior.getControl()
109 # a single |e| value for each ellipticity zone
110 ellipticityPoints = numpy.array([0.5*ctrl.ellipticityMaxInner,
111 0.5*(ctrl.ellipticityMaxInner + ctrl.ellipticityMaxOuter)])
112 # a single ln(radius) value for each logRadius zone
113 logRadiusPoints = numpy.array([0.5*(ctrl.logRadiusMinOuter + ctrl.logRadiusMinInner),
114 0.5*(ctrl.logRadiusMinInner + ctrl.logRadiusMaxInner),
115 0.5*(ctrl.logRadiusMaxInner + ctrl.logRadiusMaxOuter)])
116 # a range of position angles
117 thetaPoints = numpy.linspace(0.0, numpy.pi, 5)
118 for theta in thetaPoints:
119 for e in ellipticityPoints:
120 e1 = e*numpy.cos(2.0*theta)
121 e2 = e*numpy.sin(2.0*theta)
122 for r in logRadiusPoints:
123 self.checkDerivatives(e1, e2, r)
125 @unittest.skipIf(scipy is None, "could not import scipy")
126 def testIntegral(self):
127 """Test that the prior is properly normalized.
129 Normally, this test has a very low bar, because it's expensive to compute a high-quality
130 numerical integral to compare with. Even so, the scipy integrator does better than it
131 thinks it does, and we use that smaller tolerance for the test. That means this test
132 could fail if something about the scipy integrator changes, because we're not telling it
133 that it has to get as close as it currently is (because doing so would take way too long).
135 If this class is ever changed, we should do at least one of this test with the tolerances
136 tightened.
137 """
138 ctrl = self.prior.getControl()
139 integral, absErr = scipy.integrate.tplquad(
140 self.evaluatePrior,
141 ctrl.logRadiusMinOuter, ctrl.logRadiusMaxOuter,
142 lambda logR: -ctrl.ellipticityMaxOuter,
143 lambda logR: ctrl.ellipticityMaxOuter,
144 lambda logR, e2: -(ctrl.ellipticityMaxOuter**2 - e2**2)**0.5,
145 lambda logR, e2: (ctrl.ellipticityMaxOuter**2 - e2**2)**0.5,
146 epsabs=1.0,
147 epsrel=1.0,
148 )
149 self.assertFloatsAlmostEqual(integral, 1.0, atol=0.01)
151 def testEllipticityDistribution(self):
152 """Test that the ellipticity distribution is constant in the inner region,
153 mononotically decreasing in the ramp, and zero in the outer region, according
154 to evaluate().
155 """
156 ctrl = self.prior.getControl()
157 # a range of |e| values in each ellipticity zone
158 eInnerPoints = numpy.linspace(0.0, ctrl.ellipticityMaxInner, 5)
159 eRampPoints = numpy.linspace(ctrl.ellipticityMaxInner, ctrl.ellipticityMaxOuter, 5)
160 eOuterPoints = numpy.linspace(ctrl.ellipticityMaxOuter, ctrl.ellipticityMaxOuter + 5.0, 5)
161 # a range of position angles
162 thetaPoints = numpy.linspace(0.0, numpy.pi, 5)
163 # a single ln(radius) value for each logRadius zone
164 logRadiusPoints = numpy.array([0.5*(ctrl.logRadiusMinOuter + ctrl.logRadiusMinInner),
165 0.5*(ctrl.logRadiusMinInner + ctrl.logRadiusMaxInner),
166 0.5*(ctrl.logRadiusMaxInner + ctrl.logRadiusMaxOuter)])
167 for logRadius in logRadiusPoints:
168 for theta in thetaPoints:
169 # All inner points should have the same value
170 pInner = self.evaluatePrior(eInnerPoints*numpy.cos(2*theta),
171 eInnerPoints*numpy.sin(2*theta),
172 logRadius)
173 self.assertFloatsAlmostEqual(pInner.mean(), pInner)
175 # Each ramp point should be greater than the next one
176 pRamp = self.evaluatePrior(eRampPoints*numpy.cos(2*theta),
177 eRampPoints*numpy.sin(2*theta),
178 logRadius)
179 numpy.testing.assert_array_less(pRamp[1:], pRamp[:-1])
181 # Each outer point should be zero
182 pOuter = self.evaluatePrior(eOuterPoints*numpy.cos(2*theta),
183 eOuterPoints*numpy.sin(2*theta),
184 logRadius)
185 self.assertFloatsAlmostEqual(pOuter, 0.0)
187 def testLogRadiusDistribution(self):
188 """Test that the ln(radius) distribution is constant in the inner region,
189 mononotically decreasing in the ramps, and zero in the outer regions, according
190 to evaluate().
191 """
192 ctrl = self.prior.getControl()
193 # a range of ln(radius) values in each logRadius zone
194 rLowerOuterPoints = numpy.linspace(ctrl.logRadiusMinOuter - 2.0, ctrl.logRadiusMinOuter, 5)
195 rLowerRampPoints = numpy.linspace(ctrl.logRadiusMinOuter, ctrl.logRadiusMinInner, 5)
196 rInnerPoints = numpy.linspace(ctrl.logRadiusMinInner, ctrl.logRadiusMaxInner, 5)
197 rUpperRampPoints = numpy.linspace(ctrl.logRadiusMaxInner, ctrl.logRadiusMaxOuter, 5)
198 rUpperOuterPoints = numpy.linspace(ctrl.logRadiusMaxOuter, ctrl.logRadiusMaxOuter + 2.0, 5)
200 # a range of position angles
201 thetaPoints = numpy.linspace(0.0, numpy.pi, 5)
202 # a single |e| value for each ellipticity zone
203 ellipticityPoints = numpy.array([0.5*ctrl.ellipticityMaxInner,
204 0.5*(ctrl.ellipticityMaxInner + ctrl.ellipticityMaxOuter)])
205 for ellipticity in ellipticityPoints:
206 for theta in thetaPoints:
207 e1 = ellipticity*numpy.cos(2*theta)
208 e2 = ellipticity*numpy.sin(2*theta)
209 # Outer points should be zero
210 pLowerOuter = self.evaluatePrior(e1, e2, rLowerOuterPoints)
211 self.assertFloatsAlmostEqual(pLowerOuter, 0.0)
212 # Each ramp point should be less than the next one
213 pLowerRamp = self.evaluatePrior(e1, e2, rLowerRampPoints)
214 numpy.testing.assert_array_less(pLowerRamp[:-1], pLowerRamp[1:])
215 # All adjacent inner points should have the same distance between them (constant slope)
216 pInner = self.evaluatePrior(e1, e2, rInnerPoints)
217 diffs = pInner[1:] - pInner[:-1]
218 self.assertFloatsAlmostEqual(diffs.mean(), diffs)
219 # Each ramp point should be greater than the next one
220 pUpperRamp = self.evaluatePrior(e1, e2, rUpperRampPoints)
221 numpy.testing.assert_array_less(pUpperRamp[1:], pUpperRamp[:-1])
222 # Outer points should be zero
223 pUpperOuter = self.evaluatePrior(e1, e2, rUpperOuterPoints)
224 self.assertFloatsAlmostEqual(pUpperOuter, 0.0)
227class TestMemory(lsst.utils.tests.MemoryTestCase):
228 pass
231def setup_module(module):
232 lsst.utils.tests.init()
235if __name__ == "__main__": 235 ↛ 236line 235 didn't jump to line 236 because the condition on line 235 was never true
236 lsst.utils.tests.init()
237 unittest.main()