-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathtodo_pointsource.txt
389 lines (314 loc) · 15.3 KB
/
todo_pointsource.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
[X] Simple program to test 2D interpolation
-- Read in PSF image
-- Test interpolation shifting by whole pixels
-- Test interpolation shifting by 0.5 pixels
-- Test interpolation shifting by fractional pixels
[X] Simple derivation of FunctionObject to do point-source interpolation (PointSource)
-- assume external code generates PsfInterpolator instance
-- needs extra setup function where we pass in pointer to PsfInterpolator instance
[X] Make stub files, modify SConstruct; test compilation
[X] Add generation of instance to spline2dtest_main.cpp
[X] Write & test code for passing in pointer to PsfInterpolator instance
[X] Make alternate version of MakeShiftedImage2() in spline2dtest_main.cpp
that uses PointSource
[X] Test spline2dtest using PointSource
[X] Add point-source function to makeimage
-- ModelObject::CreateModelImage
-- ModelObject::GetSingleFunctionImage
-- [X]ModelObject::FindTotalFluxes -- OK, bcs we can call TotalFlux() method
[X] Add PointSource to add_functions.cpp, etc.
[X] Add PointSource setup to ModelObject
[X] Add normalization of localPsfPixels in ModelObject
[X] Look for places within ModelObject where we have to handle PointSource
presence/absence
[X] Exclude PointSource from general, pre-convolution image construction
[X] Handle case where there are *no* non-PointSource functions
[X] Add generation of PointSource contributions, after convolution step
[X] Rewrite GetSingleFunctionImage to handle PointSource output
-- rewrite code to do loop over functions, with
if func.is.PointSource:
generate image w/o PSF convolution
if not func.is.PointSource
generate image normally, w/ PSF convolution
[X] Test integration of PointSource image: does flux add up?
[X] Test imfit
[X] Test compilation and standard tests
[X] Test fitting with PointSource
$ ./imfit pointsource_test_200x200.fits -c config_imfit_pointsource.dat --psf=tests/psf_moffat_35.fits
CHI-SQUARE = 40535.202040 (39997 DOF)
INITIAL CHI^2 = 130759.860769
NPAR = 3
NFREE = 3
NPEGGED = 0
NITER = 7
NFEV = 25
Reduced Chi^2 = 1.013456
AIC = 40541.202640, BIC = 40566.991944
X0 120.4926 # +/- 0.0041
Y0 120.1928 # +/- 0.0041
FUNCTION PointSource
I_tot 100319 # +/- 326.27
time = 0.331s
$ ./imfit pointsource_test_200x200.fits -c config_imfit_fixedgauss.dat --psf=tests/psf_moffat_35.fits
CHI-SQUARE = 41947.858460 (39997 DOF)
INITIAL CHI^2 = 134104.030538
NPAR = 6
NFREE = 3
NPEGGED = 0
NITER = 21
NFEV = 84
Reduced Chi^2 = 1.048775
AIC = 41953.859060, BIC = 41979.648364
X0 120.4981 # +/- 0.0011
Y0 120.4066 # +/- 0.0016
FUNCTION Gaussian
PA 0 # +/- 0
ell 0 # +/- 0
I_0 1.587e+06 # +/- 5201.2
sigma 0.1 # +/- 0
time = 1.1s
# Using bicubic interpolation
$ ./imfit n3368i_star_cutout.fits -c config_imfit_n3368i_star.dat --psf=tests/psf_moffat_35.fits --save-params=bestfit_params_n3368i_star_bicubic.dat --save-residual=resid_bicubic.fits
Reduced Chi^2 = 7.989425
AIC = 79876.279987, BIC = 79897.908607
X0 49.6518 # +/- 0.0016
Y0 50.8098 # +/- 0.0016
FUNCTION PointSource
I_tot 145002 # +/- 177.07 counts
# Using lanczos2 interpolation
$ ./imfit n3368i_star_cutout.fits -c config_imfit_n3368i_star.dat --psf=tests/psf_moffat_35.fits --save-params=bestfit_params_n3368i_star_lanczos2.dat --save-residual=resid_lanczos2.fits
Reduced Chi^2 = 7.955269
AIC = 79534.822164, BIC = 79556.450784
X0 49.6637 # +/- 0.0017
Y0 50.8169 # +/- 0.0016
FUNCTION PointSource
I_tot 142029 # +/- 174.27 counts
Lanczos2 is formally a slightly better fit (DeltaAIC = -341)
[ ] Test imfit-mcmc
TESTING STAR ON HST IMAGE
Fit using FlatSky or TiltedPlane + PointSource
WFC3-IR images of CBS galaxies
ic2051 -- 381,468
n3887 -- 185,1840
n4440 -- 260,1504
n4612 -- 1773,1978 -- superbright (probably sat. w/in inner 3 pix or so)
n4984 -- 838,1823
1973,230 -- fainter, cleaner surrounds
n5363 -- 603,93
n6744 -- 483,954
n7177 -- 2168,523 (interesting as case with two faint neighbors)
134,1421 -- fainter, also has faint neighbors
* Bright star in IC2051 image:
$ ./imfit ic2051f160w_star_cutout.fits -c config_imfit_ic2051f160w_star.dat --psf=/Users/erwin/Beleriand/data/composite_bulges_sample/PSFs/psf_f160w_grizli_extracted.fits --save-params=bestfit_params_ic2051f160w_star_bicubic.dat --save-residual=resid_ic2051_bicubic.fits --save-model=model_ic2051_bicubic.fits
Reduced Chi^2 = 7.890131
AIC = 177504.381365, BIC = 177536.464669
X0 75.4919 # +/- 0.0008
Y0 74.6383 # +/- 0.0008
FUNCTION PointSource
OPTIONAL_PARAMS_START
method bicubic
OPTIONAL_PARAMS_END
I_tot 5.39718e+06 # +/- 2359.8 counts
FUNCTION FlatSky
I_sky 15.0931 # +/- 0.062134 counts/pixel
PROBLEM: outer spikes of star clearly extend outside our PSF image! (PSF image is too small...)
Visible PSF image in model is 60x57 pixels in size.
$ ./imfit ic2051f160w_star_cutout.fits -c config_imfit_ic2051f160w_star.dat --psf=/Users/erwin/Beleriand/data/composite_bulges_sample/PSFs/psf_f160w_grizli_extracted.fits --save-params=bestfit_params_ic2051f160w_star_lanczos2.dat --save-residual=resid_ic2051_lanczos2.fits --save-model=model_ic2051_lanczos2.fits
Reduced Chi^2 = 15.083530
AIC = 339327.094547, BIC = 339359.177851
X0 75.6407 # +/- 0.0009
Y0 74.4861 # +/- 0.0010
FUNCTION PointSource
OPTIONAL_PARAMS_START
method lanczos2
OPTIONAL_PARAMS_END
I_tot 5.05563e+06 # +/- 2253.9 counts
FUNCTION FlatSky
I_sky 15.231 # +/- 0.062133 counts/pixel
PROBLEM: outer spikes of star clearly extend outside our PSF image! (PSF image is too small...)
Visible PSF image in model is 56x60 pixels in size.
Bicubic is formally a better fit; residuals are better for bicubic, though not
perfect...
Actual PSF image has 60x56 as its visible size
(note that brightest pixel is at 31,31, with image size = 61,61)
* Fainter star in NGC4984 image:
$ ./imfit n4984f160w_star_cutout.fits -c config_imfit_n4984f160w_star.dat --psf=/Users/erwin/Beleriand/data/composite_bulges_sample/PSFs/psf_f160w_grizli_extracted.fits --save-params=bestfit_params_n4984f160w_star_bicubic.dat --save-residual=resid_n4984_bicubic.fits --save-model=model_n4984_bicubic.fits
Reduced Chi^2 = 6.760717
AIC = 67588.134101, BIC = 67616.971460
X0 50.6168 # +/- 0.0009
Y0 51.2036 # +/- 0.0008
FUNCTION PointSource
OPTIONAL_PARAMS_START
method bicubic
OPTIONAL_PARAMS_END
I_tot 4.82231e+06 # +/- 2267 counts
FUNCTION FlatSky
I_sky 177.199 # +/- 0.14453 counts/pixel
$ ./imfit n4984f160w_star_cutout.fits -c config_imfit_n4984f160w_star.dat --psf=/Users/erwin/Beleriand/data/composite_bulges_sample/PSFs/psf_f160w_grizli_extracted.fits --save-params=bestfit_params_n4984f160w_star_lanczos2.dat --save-residual=resid_n4984_lanczos2.fits --save-model=model_n4984_lanczos2.fits
Reduced Chi^2 = 19.099567
AIC = 190927.277633, BIC = 190956.114993
X0 50.7570 # +/- 0.0009
Y0 51.0169 # +/- 0.0007
FUNCTION PointSource
OPTIONAL_PARAMS_START
method lanczos2
OPTIONAL_PARAMS_END
I_tot 4.64736e+06 # +/- 2228.1 counts
FUNCTION FlatSky
I_sky 177.971 # +/- 0.14451 counts/pixel
$ ./imfit n4984f160w_star_cutout.fits -c config_imfit_n4984f160w_star.dat --psf=/Users/erwin/Beleriand/data/composite_bulges_sample/PSFs/psf_f160w_grizli_extracted.fits --save-params=bestfit_params_n4984f160w_star_lanczos3.dat --save-residual=resid_n4984_lanczos3.fits --save-model=model_n4984_lanczos3.fits
Reduced Chi^2 = 18.632666
AIC = 186260.133296, BIC = 186288.970656
X0 50.7394 # +/- 0.0009
Y0 51.0467 # +/- 0.0009
FUNCTION PointSource
OPTIONAL_PARAMS_START
method lanczos3
OPTIONAL_PARAMS_END
I_tot 4.71331e+06 # +/- 2246.6 counts
FUNCTION FlatSky
I_sky 177.939 # +/- 0.14451 counts/pixel
Let ModelObject handle supplying PSF images to PointSource functions?
-- ModelObject generates psfInterpolator object, passes pointer to that
into individual PointSource functions
Alternate interpolation
-- Lanczos looks simple enough to implement directly
-- Make PointSource have an optional parameter which specifies e.g.
bicubic vs Lanczos2 vs Lanczos3?
-- Lanczos3 seems to be more popular than Lanczos2?
-- zero-pad PointSource copy of PSF if we're using Lanczos
[X] Look for Python code implementing Lanczos interpolation (e.g., for
checking our code against "proper" output)
-- /Users/erwin/python/lanczos_test.py
[X] Write Python code to implement 1D Lanczos interpolation
[X] Write Python code to implement 2D Lanczos interpolation
[X] Investigate issue of normalization of Lanczos kernel
PROBLEM: Currently, the bicubic interpolation is hard-coded into
model_object.cpp and oversampled_region.cpp, which instantiate
PsfInterpolator_bicubic objects and then pass these to the FunctionObject
instances when the latter are PointSource objects.
ModelObject::SetupPsfInterpolation
OversampledRegion::AddPSFVector
(this saves a tiny amount of memory in that we don't need multiple
copies of the PsfInterpolator objects)
In ModelObject::AddFunction( FunctionObject *newFunctionObj_ptr )
// handle optional case of PointSource function
if (newFunctionObj_ptr->IsPointSource()) {
if (! psfInterpolator_allocated) {
// START NEW
string interpolationType = newFunctionObj_ptr->InterpolationType()
result = SetupPsfInterpolation(interpolationType);
// END NEW
if (result < 0)
return -1;
}
newFunctionObj_ptr->AddPsfInterpolator(psfInterpolator);
pointSourcesPresent = true;
}
[ ] Modify ModelObject::SetupPsfInterpolation to take string input
specifying interpolation type (default = "bicubic-spline")
[ ] Check if input is valid (length > 0, has one of allowed values);
return error if not
[X] Add InterpolationType method to FunctionObject
[X] Add stub method (returns "") to FunctionObject
[X] Add overridden method to PointSource
[ ] Add in Lanczsos2 interpolator
[X] Write PsfInterpolator subclass to do Lanczos2 interpolation
[X] Write unit tests for PsfInterpolator class
POSSIBLE SOLUTION: Give PointSource a method that tells caller what type
of interpolation has been requested? Then caller (e.g. ModelObject) knows
what class of PsfInterpolator to instantiate and pass to PointSource
Future options:
[] How to handle PSF oversampling regions?
-- want to make sure we don't accidentally overwrite a PointSource
with an oversampled region
-- [-]PointSource with two PSF images?
Optional boolean-setting method for FunctionObject: doing oversampled
region
-- Ensure internal flag is reset to false when Setup() is called;
new method should be called *after* Setup if we're doing oversampling
-- Problem: each oversampling region has its own (possibly unique) oversampled
PSF, so this is *not* a good idea
-- Conceptually, a PointSource is an image function like any other (e.g. Sersic),
so the user shouldn't have to specify anything unusual (*except* the actual
oversampled PSF images, which apply to regions)
-- Note that we could in principle have multiple PointSource objects within
the same oversampled-PSF region
-- Currently, ModelObject::AddFunction passes in a pointer to the main
localPsfPixels array to each new PointSource object (using their AddPsfData method).
-- So we need a way to pass in an oversampled PSF vector and tell a
PointSource object to interpolate using that instead of the main PSF
-- How do we handle point sources with centers outside oversampling region
but wings extending in to it? (We can still use low-res PointSource objects,
since it will just involve finer sampling with the interpolator)
Each oversampling region gets its own oversampled PSF (and scale). The
question is: how do we tell a given PointSource instance which oversampled
PSF to use?
-- Very general possibility: PointSource can accept pointer to psfInterpolator
object via special method (after Setup is called), so each OversampledRegion
object can keep its own psfInterpolator object and pass it to the PointSource
image function
-- PointSource::AddPsfInterpolator
-- Do this in ModelObject, too
-- This is in several ways simpler and more efficient
-- one (correct) interpolator for each oversampled region
-- no need to delete old interpolator and create new one each time
we call xxx from within ModelObject
-- Awkward: how do we allow different interpolators (e.g., GSL bicubic
vs Lanczos)?
-- Modify things to make a command-line option specifying the
interpolation mode? (Applied to all PointSource objects)
-- default is bicubic; later we can add things like
"--point-source-interpolation=lanczos3"
-- BUT: how do we take care of overall I_tot scaling?
OK: psfInterpolator object should carry its own total-counts factor
[X] Modify ModelImage to hold internal PsfInterpolator object
-- this will do interpolation for the main/global PSF
[X] Add PsfInterpolator data member (& allocation flag)
[X] Add initialization to constructor
[X] Add deallocation to destructor
[X] Run unit & regression tests!
[X] Modify ModelObject::AddFunction to remove code passing PSF pixel
data to PointSource function objects and replace it with code passing in
pointer to PsfInterpolator data member (via PointSource::AddPsfInterpolator)
[X] Modify ModelObject::CreateModelImage to add passing-in of PsfInterpolator
pointer to PointSource objects *if* oversampled regions exist (because if
they do, we'll be re-assigning the PointSource internal pointer within
the oversampling)
[X] Modify PsfInterpolator_bicubic to catch case of PSFs with x_size or y_size
< 4 (minimum size for GSL bicubic interplation)
[X] Add test to unit tests to catch the above
[X] Run unit tests!
[X] Modify ModelObject::GetSingleFunctionImage as for CreateModelImage
-- need to put this in after PSF convolution!
[X] Run unit & regression tests!
[X] Modify OversampledRegion to hold internal PsfInterpolator object
[X] Add PsfInterpolator data member (& allocation flag)
[X] Add initialization to constructor
[X] Add deallocation to destructor
[X] Run unit & regression tests!
[X] Modify OversampledRegion::AddPSFVector to create PsfInterpolator object
[X] Modify OversampledRegion::ComputeRegionAndDownsample to assign
PsfInterpolator object to PointSource function objects (if they exist)
[X] Modify OversampledRegion::ComputeRegionAndDownsample to handle
PointSource computations after regular image computation & PSF convolution
[X] Run regression tests!
[ ] Make simple config file & reference image for testing
-- One PointSource object in an oversampled region
[ ] Make sample config file with one PointSource instance, located
in what will be the oversampled region
[ ] Make reference image
[ ] make oversampled image with point source
[ ] downsample 2nd image to same resolution as base image
-- Two PointSource objects, one in an oversampled region
[ ] Make sample config file with two PointSource instances, one located
in what will be the oversampled region
[ ] Make reference image
[ ] make base image with (non-oversampled) point source
[ ] make oversampled image with point source
[ ] downsample 2nd image to same resolution as base image, then
add to base image
[X] How does this interact with weird radio-interferometric PSFs (cf. Corentin
Schreiber)? How do we scale those?
-- OK, Corentin says that as long as we don't try to normalize the PSF image,
the central peak flux can be interpreted as the total flux.