LSST Applications 24.1.6,g063fba187b+56b85ce14a,g0f08755f38+df8a265115,g12f32b3c4e+891a09f10d,g1524ad2192+7a5d7b3fbd,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g28da252d5a+07cb1400be,g2bbee38e9b+ae03bbfc84,g2bc492864f+ae03bbfc84,g3156d2b45e+6e55a43351,g347aa1857d+ae03bbfc84,g35bb328faa+a8ce1bb630,g3a166c0a6a+ae03bbfc84,g3e281a1b8c+c5dd892a6c,g414038480c+6b9177ef31,g41af890bb2+8f257c4c0b,g781aacb6e4+a8ce1bb630,g7af13505b9+7137b3b17d,g80478fca09+6df6903293,g82479be7b0+091ce1d07f,g858d7b2824+df8a265115,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,g9726552aa6+414189b318,ga5288a1d22+4a2bca08d7,gacef1a1666+c9a8ff65f4,gb58c049af0+d64f4d3760,gbcfae0f0a0+de1d42d831,gc28159a63d+ae03bbfc84,gcf0d15dbbd+72117bf34e,gda6a2b7d83+72117bf34e,gdaeeff99f8+1711a396fd,ge500cccec5+c8c9c9af63,ge79ae78c31+ae03bbfc84,gf0baf85859+c1f95f4921,gfa517265be+df8a265115,gfa999e8aa5+17cd334064,gfb92a5be7c+df8a265115
LSST Data Management Base Package
Loading...
Searching...
No Matches
test_model.py
Go to the documentation of this file.
1import lsst.gauss2d as g2d
2import lsst.gauss2d.fit as g2f
3import numpy as np
4import pytest
5
6# TODO: Investigate why clang/osx builds have non-zero values in DM-45308
7# GCC/linux are exactly 0. Hopefully it's not something like ffast-math
8max_diff_ll = 1e-25
9
10
11@pytest.fixture(scope="module")
13 return tuple(g2f.Channel(x) for x in "rgb")
14
15
16@pytest.fixture(scope="module")
17def data(channels):
18 n_x, n_y = 5, 5
19 return g2f.DataD(
20 [
21 g2f.ObservationD(
22 g2d.ImageD(np.zeros((n_y, n_x))),
23 g2d.ImageD(np.full((n_y, n_x), 1e3)),
24 g2d.ImageB(np.ones((n_y, n_x))),
25 channel,
26 )
27 for channel in channels
28 ]
29 )
30
31
32@pytest.fixture(scope="module")
33def psfmodels(data):
34 psfmodels = tuple(
36 [
38 g2f.GaussianParametricEllipse(1.0, 1.0, 0.0),
39 None,
40 g2f.LinearIntegralModel([(g2f.Channel.NONE, g2f.IntegralParameterD(1.0))]),
41 )
42 ]
43 )
44 for _ in range(len(data))
45 )
46 for psfmodel in psfmodels:
47 for param in psfmodel.parameters():
48 param.fixed = True
49 return psfmodels
50
51
52@pytest.fixture(scope="module")
55
56
57@pytest.fixture(scope="module")
58def sources(channels, interptype_sersic_default):
59 last = None
60 n_sources, n_components = 2, 2
61 sources = [None] * n_sources
62 n_comp_max = n_components - 1
63 log10 = g2f.Log10TransformD()
64 for i in range(n_sources):
65 components = [None] * n_components
66 for c in range(n_components):
67 is_last = c == n_comp_max
69 [
70 (
71 channel,
73 (is_last == 1) or 0.25 * (1 + idx_channel), fixed=is_last
74 ),
75 )
76 for idx_channel, channel in enumerate(channels)
77 ],
79 [(channel, g2f.IntegralParameterD(0.5 + 0.5 * (i + 1))) for channel in channels]
80 )
81 if (c == 0)
82 else last,
83 is_last,
84 )
85 components[c] = g2f.GaussianComponent(
86 centroid=None,
87 ellipse=g2f.GaussianParametricEllipse(1.0, 1.0, 0.0),
88 integral=last,
89 )
90 components.append(
92 centroid=None,
94 size_x=g2f.ReffXParameterD(2.6, transform=log10, fixed=False),
95 size_y=g2f.ReffYParameterD(3.5, transform=log10, fixed=False),
97 -0.1,
98 fixed=False,
99 transform=g2f.LogitLimitedTransformD(
100 limits=g2f.LimitsD(min=-0.999, max=0.999, name="ref_logit_sersic[0.5, 6.0]"),
101 ),
102 ),
103 ),
105 [
106 (
107 channel,
108 g2f.IntegralParameterD(1.0, transform=log10, label=f"Sersic {channel}-band"),
109 )
110 for channel in channels
111 ]
112 ),
113 # Linear interpolation fails at exact knot values
114 # Adding a small offset solves the problem
116 value=1.0 + 5e-4 * (interptype_sersic_default == g2f.InterpType.linear),
117 fixed=False,
118 transform=g2f.LogitLimitedTransformD(
119 limits=g2f.LimitsD(min=0.5, max=6.0, name="ref_logit_sersic[0.5, 6.0]"),
120 ),
121 ),
122 ),
123 )
124 sources[i] = g2f.Source(components)
125
126 return sources
127
128
129@pytest.fixture(scope="module")
130def model(data, psfmodels, sources):
131 model = g2f.ModelD(data, list(psfmodels), list(sources))
132 with pytest.raises(RuntimeError):
133 model.evaluate()
134 model.setup_evaluators()
135 model.evaluate()
136 result = model.evaluate()
137 for value in result:
138 assert value == 0
139 for data, output in zip(model.data, model.outputs):
140 data.image.data.flat = output.data.flat
141 return model
142
143
144def test_data(channels, data):
145 assert len(data) == len(channels)
146
147
148def test_model(channels, model):
149 assert str(model) != ""
150 gaussians = model.gaussians(channels[0])
151 assert len(gaussians) == 12
152 params = model.parameters()
153 assert len(params) == 80
154 assert len(set(params)) == 68
155
156
158 # Force is needed in case other tests use a jacobian evaluator
159 # with different parameters
160 model.setup_evaluators(g2f.EvaluatorMode.jacobian, force=True)
161 result = np.array(model.evaluate())
162 # TODO: Investigate why clang/osx builds have non-zero values in DM-45308
163 # GCC/linux are exactly 0. Hopefully it's not something like ffast-math
164 assert (np.abs(result) <= max_diff_ll).all()
165 # TODO: Investigate why this rtol needs to be so high
166 errors = model.verify_jacobian(atol=1e-3, rtol=5e-3, max_diff_ll=max_diff_ll)
167 # All of the IntegralParameters got double-counted - see DM-40674
168 # ... but also, linear interpolators will fail the Sersic index params
169 assert len(errors) == 6
170
171
173 model.setup_evaluators(g2f.EvaluatorMode.loglike_grad, force=True)
174 result = np.array(model.evaluate())
175 # TODO: Investigate why clang/osx builds have non-zero values in DM-45308
176 assert (np.abs(result) <= max_diff_ll).all()
177 dloglike = np.array(model.compute_loglike_grad(True))
178 assert (np.abs(dloglike) < max_diff_ll).all()
179
180
182 options = g2f.HessianOptions(return_negative=False)
183 params_fixed = set()
184
185 for idx in range(2):
186 hessians = {
187 is_transformed: model.compute_hessian(transformed=is_transformed, options=options)
188 for is_transformed in (False, True)
189 }
190 # Note that since there is no noise in the image and none of the
191 # param values have changed, the Hessian terms are not all negative -
192 # the sign corrections are based on the sign of loglike_grad,
193 # and those are all exactly zero
194 assert all(np.isfinite(hessian.data).all() for hessian in hessians.values())
195 if idx == 0:
196 for param in set(model.parameters()):
197 # TODO: Remove this when DM-40674 is fixed
198 if isinstance(param, g2f.IntegralParameterD) and param.free:
199 param.free = False
200 params_fixed.add(param)
201
202 for is_transformed in (False, True):
203 # can set print=True for debugging, but the output is too verbose for
204 # pytest output even on a failure
205 hessian_from_jac = model.compute_hessian(transformed=is_transformed, options=None).data
206 # TODO: investigate these unreasonably large differences after DM-40674
207 assert np.all(np.isclose(hessians[is_transformed].data, hessian_from_jac, rtol=1e-3, atol=1))
208
209 for param in params_fixed:
210 param.free = True
char * data
Definition BaseRecord.cc:61
An observational channel, usually representing some range of wavelengths of light.
Definition channel.h:29
An IntegralModel that returns a Parameter-dependent fraction of the flux of another IntegralModel.
A Component consisting of a 2D Gaussian.
A Parameter-based implementation of lsst::gauss2d::EllipseData and ParametricEllipse.
An IntegralModel with a single linear IntegralParameterD per Channel.
A Gaussian mixture model of a point spread function.
Definition psfmodel.h:26
A Gaussian mixture approximation to a Sersic profile Component.
A SersicIndexParameter for a Gaussian mixture Component.
A ParametericEllipse with effective radius Parameters.
An association of physically-related Component instances.
Definition source.h:22
interptype_sersic_default()
Definition test_model.py:53
sources(channels, interptype_sersic_default)
Definition test_model.py:58
test_model_eval_loglike_grad(model)
test_model_hessian(model)
psfmodels(data)
Definition test_model.py:33
model(data, psfmodels, sources)
test_model_eval_jacobian(model)
test_data(channels, data)