Skip to content

Commit a481765

Browse files
timstaleytimstaley
committed
BUG: Fix return-value ordering of gaussian2d_jacobian.
Add unit-test to verify analytical vs numerical derivative of gaussian2d.
1 parent 239b22f commit a481765

2 files changed

Lines changed: 148 additions & 65 deletions

File tree

‎src/fastimgproto/sourcefind/fit.py‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -275,5 +275,8 @@ def gaussian2d_jac(x, y, x_centre, y_centre, amplitude, x_stddev, y_stddev, thet
275275
dg_dtheta = g * (-(da_dtheta * xdiff2 +
276276
db_dtheta * xdiff * ydiff +
277277
dc_dtheta * ydiff2))
278-
return [dg_dA, dg_dx_mean, dg_dy_mean, dg_dx_stddev, dg_dy_stddev,
278+
279+
return [dg_dx_mean, dg_dy_mean,
280+
dg_dA,
281+
dg_dx_stddev, dg_dy_stddev,
279282
dg_dtheta]

‎tests/test_sourcefind/test_fitting.py‎

Lines changed: 144 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22

33
import inspect
44
import logging
5+
from functools import partial
56

67
import attr
78
import numpy as np
@@ -13,7 +14,11 @@
1314
check_single_source_extraction_successful,
1415
generate_random_source_params,
1516
)
16-
from fastimgproto.sourcefind.fit import gaussian2d, gaussian2d_jac
17+
from fastimgproto.sourcefind.fit import (
18+
Gaussian2dParams,
19+
gaussian2d,
20+
gaussian2d_jac,
21+
)
1722
from fastimgproto.sourcefind.image import SourceFindImage
1823

1924
logger = logging.getLogger(__name__)
@@ -22,66 +27,141 @@
2227
def test_gauss2d_func_and_jacobian():
2328
assert (inspect.getargspec(gaussian2d).args ==
2429
inspect.getargspec(gaussian2d_jac).args)
25-
scipy.optimize.check_grad(gaussian2d, gaussian2d_jac)
26-
27-
28-
29-
def test_basic_source_detection_and_fitting():
30-
ydim = 64
31-
xdim = 32
32-
image_shape = (ydim, xdim)
33-
seed = 123456
34-
35-
base_x = 18
36-
base_y = 34
37-
n_sources = 100
38-
positive_sources = generate_random_source_params(n_sources=n_sources,
39-
base_x=base_x,
40-
base_y=base_y,
41-
amplitude_range=(4.5, 42.),
42-
semiminor_range=(0.4, 2.),
43-
axis_ratio_range=(1., 5.),
44-
seed=123456,
45-
)
46-
n_islands = 0
47-
fits = []
48-
for src in positive_sources:
49-
img = np.zeros(image_shape)
50-
add_gaussian2d_to_image(src, img)
51-
detection_thresh = 4.
52-
sfimg = SourceFindImage(img, detection_n_sigma=detection_thresh,
53-
analysis_n_sigma=3.,
54-
rms_est=1.,
55-
find_negative_sources=True)
56-
57-
check_single_source_extraction_successful(src, sfimg)
58-
if sfimg.islands:
59-
n_islands += 1
60-
assert len(sfimg.islands) == 1
61-
lsq_fit = sfimg.fit_gaussian_2d(sfimg.islands[0], verbose=1)
62-
fits.append(lsq_fit)
63-
else:
64-
fits.append(None)
65-
66-
n_completed_fits = sum(1 for f in fits if f) # Count where f is not False
67-
logger.debug(
68-
"{} of {} island-fits completed".format(n_completed_fits, n_islands))
69-
70-
success = np.zeros_like(positive_sources, dtype=int)
71-
for idx, lsq_fit in enumerate(fits):
72-
if lsq_fit is None:
73-
success[idx] = 1 # Trivial no-island case. Count as OK.
74-
elif (positive_sources[idx].comparable_params ==
75-
approx(lsq_fit.comparable_params,rel=1e-2)):
76-
success[idx] = 1
77-
78-
n_successful = success.sum()
79-
logger.debug("{} of {} sources fitted accurately".format(
80-
n_successful, len(positive_sources)
81-
))
82-
# assert n_successful == len(positive_sources )
83-
84-
85-
bad_idx = np.where(success == 0)
86-
bad_truth = np.array([attr.astuple(s) for s in positive_sources])[bad_idx]
87-
bad_fits = np.array([attr.astuple(s) for s in fits])[bad_idx]
30+
31+
profile = Gaussian2dParams(x_centre=4,
32+
y_centre=9,
33+
amplitude=10.,
34+
semimajor=1.5,
35+
semiminor=1.2,
36+
theta=np.pi / 4.
37+
)
38+
# Generate some random positions to test
39+
rstate = np.random.RandomState(seed=42)
40+
n_samples = 10
41+
42+
pars = attr.astuple(profile)
43+
n_tests_run = 0
44+
for x in rstate.uniform(profile.x_centre - profile.semimajor * 5.,
45+
profile.x_centre + profile.semimajor * 5.,
46+
size=n_samples
47+
):
48+
for y in rstate.uniform(profile.y_centre - profile.semimajor * 5.,
49+
profile.y_centre + profile.semimajor * 5.,
50+
size=n_samples
51+
):
52+
# Need to partially bind the functions for each x/y location
53+
def located_g2d(x_centre, y_centre, amplitude, x_stddev, y_stddev,
54+
theta):
55+
return gaussian2d(x, y, x_centre, y_centre, amplitude,
56+
x_stddev, y_stddev,
57+
theta)
58+
59+
def located_g2d_jac(x_centre, y_centre, amplitude, x_stddev,
60+
y_stddev,
61+
theta):
62+
return gaussian2d_jac(x, y, x_centre, y_centre, amplitude,
63+
x_stddev, y_stddev,
64+
theta)
65+
66+
def wrapped_gaussian2d(pars):
67+
assert len(pars) == 6
68+
return located_g2d(*pars)
69+
70+
def wrapped_gaussian2d_jac(pars):
71+
assert len(pars) == 6
72+
return located_g2d_jac(*pars)
73+
74+
n_tests_run += 1
75+
logger.debug("Gradient check {} of {}".format(
76+
n_tests_run, n_samples ** 2))
77+
logger.debug("Checking gradient at x/y position {} {}".format(x, y))
78+
79+
eps = np.sqrt(np.finfo(float).eps) # default for check_grad
80+
# eps = 1e-12
81+
82+
numerical = scipy.optimize.approx_fprime(pars, wrapped_gaussian2d,
83+
epsilon=eps, )
84+
analytic = wrapped_gaussian2d_jac(pars)
85+
check_grad_result = scipy.optimize.check_grad(
86+
wrapped_gaussian2d,
87+
wrapped_gaussian2d_jac,
88+
x0=pars,
89+
epsilon=eps,
90+
)
91+
logger.debug("Analytic {}".format(analytic))
92+
logger.debug("Numerical {}".format(numerical))
93+
94+
logger.debug("Diffs: manual {}, check grad {}".format(
95+
np.abs(analytic - numerical),
96+
check_grad_result))
97+
logger.debug("")
98+
99+
np.testing.assert_allclose(analytic,numerical,
100+
atol=1e-7)
101+
assert check_grad_result == approx(0., abs=2e-6)
102+
103+
def test_basic_source_detection_and_fitting():
104+
ydim = 64
105+
xdim = 32
106+
image_shape = (ydim, xdim)
107+
seed = 123456
108+
109+
base_x = 18
110+
base_y = 34
111+
n_sources = 10
112+
positive_sources = generate_random_source_params(n_sources=n_sources,
113+
base_x=base_x,
114+
base_y=base_y,
115+
amplitude_range=(
116+
4.5, 42.),
117+
semiminor_range=(
118+
0.4, 2.),
119+
axis_ratio_range=(
120+
1., 5.),
121+
seed=123456,
122+
)
123+
n_islands = 0
124+
fits = []
125+
for src in positive_sources:
126+
img = np.zeros(image_shape)
127+
add_gaussian2d_to_image(src, img)
128+
detection_thresh = 4.
129+
sfimg = SourceFindImage(img, detection_n_sigma=detection_thresh,
130+
analysis_n_sigma=3.,
131+
rms_est=1.,
132+
find_negative_sources=True)
133+
134+
check_single_source_extraction_successful(src, sfimg)
135+
if sfimg.islands:
136+
n_islands += 1
137+
assert len(sfimg.islands) == 1
138+
lsq_fit = sfimg.fit_gaussian_2d(sfimg.islands[0], verbose=1)
139+
fits.append(lsq_fit)
140+
else:
141+
fits.append(None)
142+
143+
n_completed_fits = sum(
144+
1 for f in fits if f) # Count where f is not False
145+
logger.debug(
146+
"{} of {} island-fits completed".format(n_completed_fits,
147+
n_islands))
148+
149+
success = np.zeros_like(positive_sources, dtype=int)
150+
for idx, lsq_fit in enumerate(fits):
151+
if lsq_fit is None:
152+
success[idx] = 1 # Trivial no-island case. Count as OK.
153+
elif (positive_sources[idx].comparable_params ==
154+
approx(lsq_fit.comparable_params, rel=1e-2)):
155+
success[idx] = 1
156+
157+
n_successful = success.sum()
158+
logger.debug("{} of {} sources fitted accurately".format(
159+
n_successful, len(positive_sources)
160+
))
161+
# assert n_successful == len(positive_sources )
162+
163+
164+
bad_idx = np.where(success == 0)
165+
bad_truth = np.array([attr.astuple(s) for s in positive_sources])[
166+
bad_idx]
167+
bad_fits = np.array([attr.astuple(s) for s in fits])[bad_idx]

0 commit comments

Comments
 (0)