forked from mmaus96/Lens_Modeling_Auto
-
Notifications
You must be signed in to change notification settings - Fork 0
/
auto_modeling_functions_1p0.py
executable file
·1421 lines (1086 loc) · 61.7 KB
/
auto_modeling_functions_1p0.py
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
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
import os, sys, psutil
import PIL
import astropy.io.fits as pyfits
import astropy.io.ascii as ascii
import scipy
import pandas
from scipy.ndimage.filters import gaussian_filter as gauss1D
from scipy import optimize
from pylab import figure, cm
from matplotlib.colors import LogNorm
from lenstronomy.Workflow.fitting_sequence import FittingSequence
from lenstronomy.Data.imaging_data import ImageData
from lenstronomy.Data.psf import PSF
from lenstronomy.Util.util import make_grid_with_coordtransform
from lenstronomy.LensModel.lens_param import LensParam
from lenstronomy.LightModel.light_param import LightParam
from copy import deepcopy
from scipy.ndimage import gaussian_laplace
from skimage.feature import peak_local_max
from scipy.ndimage import gaussian_laplace
from matplotlib.patches import Circle
from matplotlib.colors import SymLogNorm
# from lenstronomy.Util.mask_util import mask_sphere
from lenstronomy.Util.mask_util import mask_center_2d
from lenstronomy.Util.param_util import ellipticity2phi_q
from astropy.table import QTable
from astropy.stats import sigma_clipped_stats
from scipy.spatial import distance
def calcBackgroundRMS(image_data):
'''
Calculate background root-mean-square (rms) of each image.
:image_data: List of images (dim = 3)
Returns:
:background_rms: list of floats. Background rms value of each image in image_data.
'''
#Equation for root-mean-square.
def rms(array):
return np.sqrt(np.mean(np.square(array)))
background_rms = []
for im in image_data:
#take 5x5 pixel slices from each corner of image to use for rms computation
corners = np.zeros([4,5,5])
means = np.zeros(4)
rms_array = np.zeros(4)
corners[0] = im[0:5,0:5]
corners[1] = im[-5:,0:5]
corners[2] = im[0:5,-5:]
corners[3] = im[-5:,-5:]
#mean flux of each corner
for i in range(len(corners)):
means[i] = np.mean(corners[i])
#rms of each corner
for i in range(len(corners)):
rms_array[i] = rms(corners[i])
# print('Mean values: ', means, '\n',
# 'rms values: ',rms_array)
#Exclude corners with average flux > twice the minimum flux
rms_good = [rms_array[means == means.min()]]
for x in means[means != means.min()]:
if (x <= means.min() + 1.0 * np.abs(means.min())):
rms_good.append(rms_array[means == x])
# print('good rms: ',rms_good, '\n',
# 'background_rms: ',np.mean(rms_good))
background_rms.append(np.mean(rms_good)) #array of background rms values of each image
return background_rms
def prepareData(lens_info,lens_data,psf_data):
'''
Takes raw data, psf, noise map, etc and formats them as lists of dictionaries, "kwargs_data" and "kwargs_psf", which can then
be passed to Lenstronomy.
:lens_info: list of dictionaries (length = # of bands). Contains info about data in each band, including pixel scale
(deltaPix in arcsec/pixel), psf upsample factor, data array length (numPix), exposure time, noise map, and background
rms.
:lens_data: raw image data arrays with shape (# of bands, numPix, numPix) (e.g. [data_band1, data_band2, ...], where each
data_band is a 2D array of the image in that band).
:psf_data: psf data arrays in the form [psf_band1, psf_band2, ...]
Returns:
:kwargs_data: list of dictionaries with length = # bands, containing data and info.
:kwargs_psf: list of dictionaries with length = # bands, containing psf data and info.
'''
kwargs_data = []
kwargs_psf = []
#loop though bands
for i,x in enumerate(lens_info):
#define image coordinate system/grid using image dimensions (numPix) and pixel scale (deltaPix)
_,_,ra_at_xy_0,dec_at_xy_0,_,_,transform_pix2angle,_ = make_grid_with_coordtransform(numPix = x['numPix'],
deltapix = x['deltaPix'])
#make kwargs_data dictionary for band[i] with image data, coordinate grid, and info from lens_info
kwargs_data.append({'image_data': lens_data[i],
'background_rms': x['background_rms'],
'exposure_time': x['exposure_time'],
'noise_map': x['noise_map'],
'ra_at_xy_0':ra_at_xy_0,
'dec_at_xy_0': dec_at_xy_0,
'transform_pix2angle': transform_pix2angle})
#make kwargs_psf for band[i] with psf data and relevant info
kwargs_psf.append({'psf_type': x['psf_type'],
'kernel_point_source': psf_data[i],
'point_source_supersampling_factor': x['psf_upsample_factor']})
return kwargs_data, kwargs_psf
def removekeys(d, keys):
'''
Create new dictionary, r, out of dictionary d with keys removed.
'''
r = dict(d)
for k in keys:
del r[k]
return r
def openFITS(filepath):
'''
Opens FITS files and reads data and headers. Supports datacubes with multiple bands in one FITS file.
'''
data = []
hdr = []
with pyfits.open(filepath) as file:
for f in file:
data.append(deepcopy(f.data))
hdr.append(deepcopy(f.header))
del f.data
file.close()
return data, hdr
def printMemory(location):
'''
Print total memory in Megabytes currently being used by script.
:location: string. Just so you know where in the modeling process this amount of memory is currently being used.
'''
print('$$$$$$$$$$$$$$$$$$$ Memory Usage ({}) $$$$$$$$$$$$$$$$$$$$$$$$$$'.format(location))
process = psutil.Process(os.getpid())
print('{} Megabytes'.format(float(process.memory_info().rss) / 2**(20))) # in Megabytes
print('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$')
def optParams(kwargs_init,opt_params,model_kwarg_names):
'''
Function for setting up init and fixed dictionaries for optimizing specific model parameters and fixing the rest. Can be
used for making a chain of PSO fits in pre-optimization scheme as desired in a less messy way. Used in optimize_dynamic.py
script.
:kwargs_init: initial parameter values for modeling. Same format as kwargs_result after a fit.
:opt_params: dictionary with 'kwargs_lens','kwargs_source', and 'kwargs_lens_light' as keys. Values are arrays with names
of parameters to optimize in the next fit.
:model_kwarg_names: Dictionary. All parameter names in profiles in source, lens, and lens_light model lists.
Returns:
:args_init: dictionary of initialization parameters
:args_fixed: dictionary of fixed parametes
'''
fixed_lens = []
kwargs_lens_init = []
for i in range(len(model_kwarg_names['kwargs_lens'])): # num iterations correspond to num profiles in lens_model_list
kwargs_lens_init.append(kwargs_init['kwargs_lens'][i])
#fix all params (and their values) that are not being optimized
fixed_lens.append(removekeys(kwargs_init['kwargs_lens'][i],opt_params['kwargs_lens'][i]))
fixed_source = []
kwargs_source_init = []
for i in range(len(model_kwarg_names['kwargs_source'])):
kwargs_source_init.append(kwargs_init['kwargs_source'][i])
fixed_source.append(removekeys(kwargs_init['kwargs_source'][i],opt_params['kwargs_source'][i]))
fixed_lens_light = []
kwargs_lens_light_init = []
for i in range(len(model_kwarg_names['kwargs_lens_light'])):
kwargs_lens_light_init.append(kwargs_init['kwargs_lens_light'][i])
fixed_lens_light.append(removekeys(kwargs_init['kwargs_lens_light'][i],opt_params['kwargs_lens_light'][i]))
args_init = {'kwargs_lens': kwargs_lens_init,
'kwargs_source': kwargs_source_init,
'kwargs_lens_light': kwargs_lens_light_init}
args_fixed = {'kwargs_lens': fixed_lens,
'kwargs_source': fixed_source,
'kwargs_lens_light': fixed_lens_light}
return args_init, args_fixed
def get_kwarg_names(lens_model_list,source_model_list,lens_light_model_list,kwargs_fixed = None):
'''
Function to get all parameter names for specific profiles used in lens, source, and lens_light model_lists.
:lens_model_list:list of strings. Profiles used for modeling lens mass
:source_model_list:list of strings. Profiles used for modeling source light
:lens_light_model_list:list of strings. Profiles used for modeling lens light
:kwargs_fixed: dictionary corresponding to fixed parameters in lens, source, and lens light profiles.
Returns:
:model_kwarg_names: dictionary with parameter names of lens, source, and lens_light profiles
'''
if kwargs_fixed == None:
kwargs_fixed = {'kwargs_lens': [], 'kwargs_source': [] , 'kwargs_lens_light': []}
for i in range(len(lens_model_list)): kwargs_fixed['kwargs_lens'].append({})
for j in range(len(source_model_list)): kwargs_fixed['kwargs_source'].append({})
for k in range(len(lens_light_model_list)): kwargs_fixed['kwargs_lens_light'].append({})
lens_model = LensParam(lens_model_list,kwargs_fixed['kwargs_lens'])
lens_params_list = lens_model._param_name_list
lens_params=np.array([np.array(xi) for xi in lens_params_list])
source_model = LightParam(source_model_list,kwargs_fixed['kwargs_source'])
source_params_list = source_model._param_name_list
source_params = np.array([np.array(xi) for xi in source_params_list])
lens_light_model = LightParam(lens_light_model_list,kwargs_fixed['kwargs_lens_light'])
lens_light_params_list = lens_light_model._param_name_list
lens_light_params = np.array([np.array(xi) for xi in lens_light_params_list])
model_kwarg_names = {'kwargs_lens': np.array(lens_params),
'kwargs_source': np.array(source_params),
'kwargs_lens_light': np.array(lens_light_params)}
return model_kwarg_names
def prior_phi_q_gaussian(kwargs_list, prior_list):
"""
Function for placing gaussian priors on aspect ratio q by first converting e1 and e2 ellipticity parameters (used in
Lenstronomy modeling) to q and phi.
:param kwargs_list: keyword argument list
:param prior_list: prior list
:return: logL
"""
logL = 0
if not kwargs_list:
pass #So nothing crashes if there is not lens light or source light model
else:
for i in range(len(prior_list)):
index, param_name, value, sigma = prior_list[i]
if (('e1' in kwargs_list[index]) and ('e2' in kwargs_list[index])):
model_e1 = kwargs_list[index]['e1']
model_e2 = kwargs_list[index]['e2']
model_vals = {}
model_vals['phi'], model_vals['q'] = ellipticity2phi_q(model_e1,model_e2)
dist = (model_vals[param_name] - value) ** 2 / sigma ** 2 / 2
logL -= np.sum(dist)
# print('prior: {} \n model value: {} \n mean value: {} \n sigma: {}'.format(param_name,
# model_vals[param_name],
# value,sigma))
else:
pass
return logL
def join_param_between_bands(kwargs_list, prior_list):
"""
Gaussian prior for joining parameter in one band to that of another band.
:param kwargs_list: keyword argument list
:param prior_list: prior list
:return: logL
"""
logL = 0
if not kwargs_list:
pass
else:
for i in range(len(prior_list)):
mean_index,index,param,sigma = prior_list[i]
model_val = kwargs_list[index][param] #model parameter on which to place prior
mean_val = kwargs_list[mean_index][param] #mean value around which gaussian is centered
dist = (model_val - mean_val)**2 / sigma ** 2 / 2
logL -= np.sum(dist)
# print('prior: {} \n model value: {} \n mean value: {} \n sigma: {}'.format(param,model_val,mean_val,sigma))
return logL
def join_lens_with_light_loose(kwargs_lens,kwargs_lens_light,prior_list):
"""
Gaussian Prior for joining parameter in lens mass to parameter in lens light
:param kwargs_lens (lens_light): keyword argument list of lens (lens_light) profiles
:param prior_list: prior list
:return: logL
"""
logL = 0
if not kwargs_lens:
pass
elif not kwargs_lens_light:
pass
else:
for i in range(len(prior_list)):
light_index, lens_index, param, sigma = prior_list[i]
model_val = kwargs_lens[lens_index][param]
mean_val = kwargs_lens_light[light_index][param]
dist = (model_val - mean_val)**2 / sigma ** 2 / 2
logL -= np.sum(dist)
# print('prior: {}_lens \n lens_val: {} \n light value: {} \n sigma: {}'.format(param,model_val,mean_val,sigma))
return logL
def custom_logL_addition(kwargs_lens=None, kwargs_source=None, kwargs_lens_light=None, kwargs_ps=None, kwargs_special=None,
kwargs_extinction=None):
'''
custom_logL_addition function. Joins e1,e2 parameters between bands using gaussian priors.
Joins e1,e2 parameters of lens mass to lens light profile with gaussian priors.
Puts gaussian priors on q for source and lens light sersic profiles
:kwargs_*: kwargs lists
:return: logL
'''
#make list of indices in kwargs_source corresponding to sersic profile
if not kwargs_source:
sersic_indices_source = []
else:
sersic_indices_source = [i for i,x in enumerate(kwargs_source) if 'R_sersic' in x]
#make list of indices in kwargs_lens_light corresponding to sersic profile
if not kwargs_lens_light:
sersic_indices_lens_light = []
else:
sersic_indices_lens_light = [i for i,x in enumerate(kwargs_lens_light) if 'R_sersic' in x]
#prior list for source gaussian prior on q
prior_list_source = []
for i,x in enumerate(sersic_indices_source):
prior_list_source.append([x,'q',0.8,0.1]) #[index,pram,mean,sigma]
#prior list for joining source ellipticity between bands
join_list_source = []
for i in range(len(sersic_indices_source) - 1):
#gaussian centered on PREVIOUS band:
join_list_source.append([sersic_indices_source[i],sersic_indices_source[i+1],'e1',0.01]) #[mean_index,index,param,sigma]
join_list_source.append([sersic_indices_source[i],sersic_indices_source[i+1],'e2',0.01])
#gaussian centered on FIRST band:
# join_list_source.append([0,sersic_indices_source[i+1],'e1',0.01])
# join_list_source.append([0,sersic_indices_source[i+1],'e2',0.01])
#prior list for lens light gaussian prior on q
prior_list_lens_light = []
for i,x in enumerate(sersic_indices_lens_light):
prior_list_lens_light.append([x,'q',0.8,0.1]) #[mean_index,index,param,sigma]
#prior list for joining lens light ellipticity between bands
join_list_lens_light = []
for i in range(len(sersic_indices_lens_light) - 1):
#gaussian centered on PREVIOUS band: #[mean_index,index,param,sigma]
join_list_lens_light.append([sersic_indices_lens_light[i],sersic_indices_lens_light[i+1],'e1',0.01])
join_list_lens_light.append([sersic_indices_lens_light[i],sersic_indices_lens_light[i+1],'e2',0.01])
#gaussian centered on FIRST band:
# join_list_lens_light.append([0,sersic_indices_lens_light[i+1],'e1',0.01])
# join_list_lens_light.append([0,sersic_indices_lens_light[i+1],'e2',0.01])
#prior list for joining lens mass ellipticity with lens light ellipticity using gaussian prior
if not kwargs_lens:
join_list_lens_with_light = []
else:
join_list_lens_with_light = [[sersic_indices_lens_light[-1],0,'e1',0.01], #[light_index, lens_index, param, sigma]
[sersic_indices_lens_light[-1],0,'e2',0.01]]
logL = 0
# print('Prior on q (lens_light)')
logL += prior_phi_q_gaussian(kwargs_lens_light, prior_list_lens_light)
# print('Prior on q (source)')
logL += prior_phi_q_gaussian(kwargs_source, prior_list_source)
# print('Loose join of lens with light e1 and e2')
logL += join_lens_with_light_loose(kwargs_lens,kwargs_lens_light,join_list_lens_with_light)
# print('Loose join of e1 and e2 between source bands')
logL += join_param_between_bands(kwargs_source, join_list_source)
# print('Loose join of e1 and e2 between lens_light bands')
logL += join_param_between_bands(kwargs_lens_light, join_list_lens_light)
return logL
def prepareFit(kwargs_data, kwargs_psf, lens_model_list,
source_model_list, lens_light_model_list, image_mask_list = None):
'''
Function for preparing fit settings in appropriate kwargs to pass into Lenstronomy.
:kwargs_data: list of dicts; image data and info
:kwargs_psf:list of dicts; psf data and info
:lens_model_list: list of strings; lens model profiles for a single band
:source_model_list: list of strings; source model profiles for a single band
:lens_light_model_list: list of strings; lens light model profiles for a single band
:image_mask_list: List of masks. One for each band. Masks are 2D arrays of ints, with values 0 or 1
Returns:
:kwargs_likelihood: dictionary with info for likelihood computations (mask list and priors)
:kwargs_model: dictionary for model lists, indices of profiles in model lists to be used for each band
:kwargs_data_joint: dictionary containing multi_band_list and multi band type
:multi_band_list: list with length = # bands. Each item (i) is another list with i'th kwargs_data, kwargs_psf, and
kwargs_numerics for that band.
:kwargs_constraints: dictionary for constraints (e.g. joining parameters between profiles)
'''
kwargs_likelihood = {'source_marg': False,
'image_likelihood_mask_list': image_mask_list,
'prior_lens_light': [],
# 'prior_lens_light_lognormal': [],
'prior_source': [],
# 'check_positive_flux': True,
'custom_logL_addition': custom_logL_addition,
}
multi_source_model_list = []
multi_lens_light_model_list = []
#model lists with profiles repeated for number of bands
for i in range(len(kwargs_data)):
multi_source_model_list.extend(deepcopy(source_model_list))
multi_lens_light_model_list.extend(deepcopy(lens_light_model_list))
kwargs_model = {'lens_model_list': deepcopy(lens_model_list),
'source_light_model_list': multi_source_model_list,
'lens_light_model_list': multi_lens_light_model_list,
'index_lens_model_list': [[i for i in range(len(lens_model_list))] for x in kwargs_data],
'index_source_light_model_list': [[i+j*len(source_model_list) for i in range(len(source_model_list))]
for j in range(len(kwargs_data))],
'index_lens_light_model_list': [[i+j*len(lens_light_model_list) for i in
range(len(lens_light_model_list))] for j in range(len(kwargs_data))]}
#indices in multi_source_model_list corresponding to sersic profile:
SERSIC_indices_source = [i for i,x in enumerate(kwargs_model['source_light_model_list']) if x == 'SERSIC_ELLIPSE']
#indices in multi_source_model_list corresponding to shapelets profile:
SHAPELETS_indices_source = [i for i,x in enumerate(kwargs_model['source_light_model_list']) if x == 'SHAPELETS']
#list of indices of source profiles
indices_source = [i for i in range(len(kwargs_model['source_light_model_list']))]
#indices in multi_lens_light_model_list corresponding to sersic profile:
SERSIC_indices_lens_light = [i for i,x in enumerate(kwargs_model['lens_light_model_list']) if x == 'SERSIC_ELLIPSE']
kwargs_constraints = {'joint_source_with_source': [],
'joint_lens_light_with_lens_light': [],
'joint_lens_with_light': []}
print('Check 1')
#Join lens_light centroids between bands
for i in range(len(SERSIC_indices_lens_light) - 1):
kwargs_constraints['joint_lens_light_with_lens_light'].append([SERSIC_indices_lens_light[0],
SERSIC_indices_lens_light[i+1],
# ['center_x','center_y','e1','e2']
['center_x','center_y']
])
print('Check 2')
#Join source centroids between bands
for i in range(len(indices_source) - 1):
kwargs_constraints['joint_source_with_source'].append([indices_source[0],indices_source[i+1],['center_x','center_y']])
#Join source
# for i in range(len(SERSIC_indices_source) - 1):
# kwargs_constraints['joint_source_with_source'].append([SERSIC_indices_source[i],SERSIC_indices_source[i+1],['e1','e2']])
for i,j in enumerate(SERSIC_indices_source):
kwargs_likelihood['prior_source'].append([j,'R_sersic', 1.0,2.0])
kwargs_likelihood['prior_source'].append([j,'n_sersic', 3.0,1.0])
# # kwargs_likelihood['prior_source'].append([i,'e1', 0.,0.15])
# # kwargs_likelihood['prior_source'].append([i,'e2', 0.,0.15])
for i,j in enumerate(SHAPELETS_indices_source):
kwargs_likelihood['prior_source'].append([j,'beta', 0.01,0.05])
for i in range(len(multi_lens_light_model_list)):
if len(lens_model_list) != 0:
kwargs_constraints['joint_lens_with_light'].append([i,0,['center_x','center_y']])
# kwargs_constraints['joint_lens_with_light'].append([i,0,['center_x','center_y','e1','e2']])
kwargs_likelihood['prior_lens_light'].append([i,'R_sersic', 1.0,2.0])
kwargs_likelihood['prior_lens_light'].append([i,'n_sersic', 3.0,1.0])
# # kwargs_likelihood['prior_lens_light'].append([i,'e1', 0.,0.15])
# # kwargs_likelihood['prior_lens_light'].append([i,'e2', 0.,0.15])
multi_band_list = []
for i in range(len(kwargs_data)):
kwargs_numerics = {'supersampling_factor': 2, # each pixel gets super-sampled (in each axis direction)
'supersampling_convolution': False}
if kwargs_psf[i]['point_source_supersampling_factor'] != 1:
kwargs_numerics['supersampling_convolution'] = True
kwargs_numerics['supersampling_factor'] = kwargs_psf[i]['point_source_supersampling_factor']
multi_band_list.append([kwargs_data[i], kwargs_psf[i], kwargs_numerics]) # if you have multiple bands to be modeled simultaneously,you can append them to the mutli_band_list
kwargs_data_joint = {'multi_band_list': multi_band_list,
'multi_band_type': 'multi-linear' # 'multi-linear': every imaging band has independent solutions of the surface brightness, 'joint-linear': there is one joint solution of the linear coefficients demanded across the bands.
}
return kwargs_likelihood, kwargs_model, kwargs_data_joint, multi_band_list, kwargs_constraints
def runFit(fitting_kwargs_list, kwargs_params, kwargs_likelihood, kwargs_model, kwargs_data_joint, kwargs_constraints = {}):
'''
Function for running a fit.
:fitting_kwargs_list: list. Type of fit (PSO or MCMC), and fitting parameters
:kwargs_params: Dictionary. Parameters (init, fixed, sigma, upper/lower bounds) for lens, lens light, and source light
profiles
:kwargs_likelihood: Dictionary with info for likelihood computations (mask list and priors)
:kwargs_model: dictionary for model lists, indices of profiles in model lists to be used for each band
:kwargs_data_joint: dictionary containing multi_band_list and multi band type
:kwargs_constraints: dictionary for constraints (e.g. joining parameters between profiles)
Returns:
:chain_list: fitting results
:kwargs_result: optimized model parameter values
'''
#set up fitting sequence
fitting_seq = FittingSequence(kwargs_data_joint, kwargs_model, kwargs_constraints, kwargs_likelihood, kwargs_params)
#print fit settings before starting fit
fitting_seq.param_class.print_setting()
#run fit
chain_list = fitting_seq.fit_sequence(fitting_kwargs_list)
kwargs_result = fitting_seq.best_fit()
return chain_list, kwargs_result
def find_components_old(image,deltaPix,lens_rad_arcsec = 5.0,lens_rad_ratio = None, gal_rad_ratio = 0.1,min_size_arcsec=0.3,thresh=0.4, show_locations=False):
"""
Detect multiple components in a pixelated model of the source and return their coordinates.
min_size_arcsec: minimum size of features to be considered as individual components
thresh: threshold all values below `tresh` times the max of the image after LoG filtering
"""
# convert minimum component size in pixel units
min_size = int(min_size_arcsec / deltaPix)
#Convert lens radius and central galaxy radius to pixels
if lens_rad_ratio == None:
lens_rad = int(lens_rad_arcsec / deltaPix)
else: lens_rad = int(len(image) * lens_rad_ratio)
gal_rad = int(len(image) * gal_rad_ratio)
# downscale source image to data resolution (for speed + easier for converting to data units)
#down = image_util.re_size(image, factor=supersampling_factor_source)
# apply laplacian of gaussian (LoG) filter to enhance maxima
filtered = - gaussian_laplace(deepcopy(image), sigma = min_size, mode='constant', cval=0.)
# print(filtered.min(),filtered.max(),filtered.min() + thresh * np.abs(filtered.min()))
# assume all value below max*threshold can not be maxima, so put all to zero
# filtered[filtered < thresh*filtered.max()] = 0.
# assume all value below min*threshold can not be maxima, so put all to zero
filtered[filtered < filtered.min() + thresh * np.abs(filtered.min())] = 0.
if show_locations:
plt.figure(figsize = (8,8))
plt.subplot(1,2,1)
plt.imshow(image, origin='lower', norm=SymLogNorm(5))
plt.title('Image')
plt.subplot(1,2,2)
plt.imshow(filtered, origin='lower', norm=SymLogNorm(5))
plt.title('Filtered Image')
plt.show()
# find coordinates of local maxima
#print(int(0.5 * min_size))
max_idx_2d_small = peak_local_max(filtered, min_distance=0)
max_idx_2d_large = peak_local_max(filtered, min_distance=1)
x_list_small, y_list_small = max_idx_2d_small[:, 1], max_idx_2d_small[:, 0]
x_list_large, y_list_large = max_idx_2d_large[:, 1], max_idx_2d_large[:, 0]
im_center_x, im_center_y = len(image) / 2., len(image) / 2.
R = np.sqrt((x_list_large - im_center_x)**2 + (y_list_large - im_center_y)**2)
new_center_x, new_center_y = x_list_large[R < gal_rad], y_list_large[R < gal_rad]
if (len(new_center_x) > 1) and (len(x_list_large[R == R.min()]) ==1 ):
new_center_x, new_center_y = x_list_large[R == R.min()], y_list_large[R == R.min()]
elif (len(new_center_x) > 1) and (len(x_list_large[R == R.min()]) > 1 ):
new_center_x, new_center_y = im_center_x, im_center_y
elif len(new_center_x) == 0:
new_center_x, new_center_y = im_center_x, im_center_y
R_small = np.sqrt((x_list_small - new_center_x)**2 + (y_list_small - new_center_y)**2)
R_large = np.sqrt((x_list_large - new_center_x)**2 + (y_list_large - new_center_y)**2)
x_sats, y_sats = x_list_small[R_small > lens_rad], y_list_small[R_small > lens_rad]
# show maxima on image for debug
if show_locations:
fig = plt.figure(figsize=(4, 4))
#plt.imshow(image, origin='lower', cmap=cmap_flux, norm=LogNorm(1e-2))
plt.imshow(image, origin='lower', norm=SymLogNorm(5))
for i in range(len(x_sats)):
plt.scatter([x_sats[i]], [y_sats[i]], c='red', s=60, marker='+')
# plt.annotate(i+1, (x_list[i], y_list[i]), color='black')
# for i in range(len(x_mask)):
# plt.scatter([x_mask[i]], [y_mask[i]], c='red', s=100, marker='*')
# plt.annotate(i+1, (x_mask[i], y_mask[i]), color='red')
plt.scatter(new_center_x, new_center_y,c='red', s=100, marker='*')
draw_lens_circle = Circle((new_center_x, new_center_y),lens_rad ,fill=False)
draw_gal_circle = Circle((new_center_x, new_center_y),gal_rad, fill = False)
plt.gcf().gca().add_artist(draw_lens_circle)
plt.gcf().gca().add_artist(draw_gal_circle)
plt.title('Detected Components')
plt.text(1, 1, "detected components", color='red')
fig.axes[0].get_xaxis().set_visible(True); fig.axes[0].get_yaxis().set_visible(True)
plt.show()
return (x_sats, y_sats), (new_center_x, new_center_y)
def mask_for_sat_old(image,deltaPix,lens_rad_arcsec = 5.0,lens_rad_ratio = 0.4, show_plot = False):
gal_rad_ratio = 0.1
min_size_arcsec= 0.7
thresh = 1.4
numPix = len(image)
(x_sats, y_sats), (new_center_x, new_center_y) = find_components(image, deltaPix,lens_rad_arcsec = lens_rad_arcsec,
gal_rad_ratio = gal_rad_ratio,
min_size_arcsec= min_size_arcsec,
lens_rad_ratio = lens_rad_ratio,
thresh=thresh, show_locations=show_plot)
satellites = []
mask_final = np.ones([numPix,numPix])
for i in range(len(x_sats)):
satellites.append(np.zeros([numPix,numPix]))
for j in range(len(image)):
# satellites[i][j] = mask_sphere(j,np.linspace(0,numPix - 1,numPix),y_sats[i],x_sats[i],2)
satellites[i][j] = mask_center_2d(x_sats[i], y_sats[i], 2, np.linspace(0,numPix - 1,numPix), j)
# mask_final[satellites[i] ==1] = 0
mask_final[satellites[i] ==0] = 0
#mask_list.append(mask_final)
#mask_final = np.ones([numPix,numPix])
#mask_final[mask_both == 1] = 0.
return mask_final
def mask_for_lens_gal(image,deltaPix, gal_rad_ratio = 0.1, show_plot = False):
lens_rad_arcsec = 3.0
min_size_arcsec= 0.7
thresh = 1.4
lens_rad_ratio = 0.4
gal_rad = int(len(image) * gal_rad_ratio)
numPix = len(image)
(x_sats, y_sats), (new_center_x, new_center_y) = find_components(image, deltaPix,lens_rad_arcsec = lens_rad_arcsec,
gal_rad_ratio = gal_rad_ratio,
min_size_arcsec= min_size_arcsec,
lens_rad_ratio = lens_rad_ratio,
thresh=thresh, show_locations=show_plot)
# print(new_center_x, new_center_y)
mask = np.zeros([numPix,numPix])
for i in range(len(image)):
# mask[i] = mask_sphere(i,np.linspace(0,numPix - 1,numPix),new_center_y,new_center_x,gal_rad)
lens_gal[i] = mask_center_2d(new_center_x,new_center_y, gal_rad, np.linspace(0,numPix - 1,numPix), i)
mask[lens_gal[i] == 0] = 1
return mask
def find_lens_gal(image,deltaPix,show_plot=False,title=None):
corners = np.zeros([4,5,5])
corners[0] = image[0:5,0:5]
corners[1] = image[-5:,0:5]
corners[2] = image[0:5,-5:]
corners[3] = image[-5:,-5:]
means, medians, stds = [],[],[]
for c in corners:
mn,med,s = sigma_clipped_stats(c,sigma=3.0)
means.append(mn)
medians.append(med)
stds.append(s)
mean_bg = np.mean(means)
median_bg = np.mean(medians)
std_bg = np.mean(stds)
# print('means: {}, median: {}, std: {}'.format(means,median_bg,std_bg))
thresh = std_bg / 4
min_size_arcsec = 0.7
min_size = int(min_size_arcsec / deltaPix)
LoG = - gaussian_laplace(deepcopy(image), sigma = min_size, mode='constant', cval=0.)
filtered = deepcopy(LoG)
corners = np.zeros([4,5,5])
corners[0] = LoG[0:5,0:5]
corners[1] = LoG[-5:,0:5]
corners[2] = LoG[0:5,-5:]
corners[3] = LoG[-5:,-5:]
means = []
for c in corners:
mn,med,s = sigma_clipped_stats(c,sigma=3.0)
means.append(mn)
means = np.array(means)
means_std = np.std(means)
means_good = means[(means >= means.mean() - 1.0 * means_std) & (means <= means.mean() + 1.0 * means_std)]
mean_bg = np.mean(means_good)
filtered[filtered < mean_bg + thresh] = 0.
max_idx_2d_large = peak_local_max(filtered, min_distance=1)
x_list_large, y_list_large = max_idx_2d_large[:, 1], max_idx_2d_large[:, 0]
gal_rad = 5.
im_center_x, im_center_y = len(image) / 2., len(image) / 2.
R = np.sqrt((x_list_large - im_center_x)**2 + (y_list_large - im_center_y)**2)
new_center_x, new_center_y = deepcopy(x_list_large[R < gal_rad]), deepcopy(y_list_large[R < gal_rad])
if (len(new_center_x) == 1):
print('c1:')
print(new_center_x[0], new_center_y[0])
elif (len(new_center_x) > 1):
if (len(x_list_large[R == R.min()]) ==1 ):
new_center_x, new_center_y = x_list_large[R == R.min()], y_list_large[R == R.min()]
print('c2:')
print(new_center_x[0], new_center_y[0])
elif (len(x_list_large[R == R.min()]) > 1):
new_center_x, new_center_y = x_list_large[R == R.min()][0], y_list_large[R == R.min()][0]
print('c3:')
print(new_center_x[0], new_center_y[0])
elif (len(new_center_x) == 0):
gal_rad = 10
new_center_x, new_center_y = x_list_large[R < gal_rad], y_list_large[R < gal_rad]
if (len(new_center_x) > 1):
if (len(x_list_large[R == R.min()]) ==1 ):
new_center_x, new_center_y = x_list_large[R == R.min()], y_list_large[R == R.min()]
print('c4:')
print(new_center_x[0], new_center_y[0])
elif (len(x_list_large[R == R.min()]) > 1):
new_center_x, new_center_y = x_list_large[R == R.min()][0], y_list_large[R == R.min()][0]
print('c5:')
print(new_center_x[0], new_center_y[0])
elif (len(new_center_x) == 0):
new_center_x, new_center_y = im_center_x, im_center_y
print('c6:')
print(new_center_x[0], new_center_y[0])
if show_plot:
plt.figure(figsize=(8,6))
plt.imshow(image,origin='lower',norm=SymLogNorm(5))
circle = Circle((new_center_x, new_center_y),gal_rad, fill = False)
plt.gcf().gca().add_artist(circle)
plt.scatter(new_center_x, new_center_y,c='red', s=100, marker='*')
plt.text(1, 1, "loc = ({},{})".format(new_center_x[0], new_center_y[0]),fontsize=30,fontweight='bold',color='red')
plt.axis('off')
plt.title(title)
return new_center_x[0], new_center_y[0]
def find_components(image,deltaPix,lens_rad_arcsec = 6.0,lens_rad_ratio = None,
center_x = None,center_y = None, gal_rad_ratio = 0.1,
min_size_arcsec=0.7,thresh=0.5, many_sources = True,
show_locations=False, title = None):
"""
Detect multiple components in a pixelated model of the source and return their coordinates.
min_size_arcsec: minimum size of features to be considered as individual components
thresh: threshold all values below `tresh` times the max of the image after LoG filtering
"""
# convert minimum component size in pixel units
min_size = int(min_size_arcsec / deltaPix)
#Convert lens radius and central galaxy radius to pixels
if lens_rad_ratio == None:
lens_rad = int(lens_rad_arcsec / deltaPix)
else: lens_rad = int(len(image) * lens_rad_ratio)
gal_rad = int(len(image) * gal_rad_ratio)
# im2[im2 < im2.min() + 10.*thresh] = 0.
# downscale source image to data resolution (for speed + easier for converting to data units)
#down = image_util.re_size(image, factor=supersampling_factor_source)
# apply laplacian of gaussian (LoG) filter to enhance maxima
LoG = - gaussian_laplace(deepcopy(image), sigma = min_size, mode='constant', cval=0.)
# LoG = - gaussian_laplace(deepcopy(im2), sigma = 2., mode='constant', cval=0.)
filtered = deepcopy(LoG)
# print(LoG.min(),LoG.max(),np.abs(LoG.min()) + thresh )
# print(type(filtered))
corners = np.zeros([4,5,5])
corners[0] = LoG[0:5,0:5]
corners[1] = LoG[-5:,0:5]
corners[2] = LoG[0:5,-5:]
corners[3] = LoG[-5:,-5:]
means = []
for c in corners:
mn,med,s = sigma_clipped_stats(c,sigma=3.0)
means.append(mn)
means = np.array(means)
means_std = np.std(means)
means_good = means[(means >= means.mean() - 1.0 * means_std) & (means <= means.mean() + 1.0 * means_std)]
mean_bg = np.mean(means_good)
# print('LoG means: {}, Log means std: {}, Log means good: {}, LoG avg mean: {}'.format(means,means_std,means_good,mean_bg))
# print('min: {}, max: {}, cut: {}'.format(LoG.min(),LoG.max(),mean_bg + thresh))
# print(LoG.min(),LoG.max(),filtered.min() + thresh)
# assume all value below max*threshold can not be maxima, so put all to zero
# filtered[filtered < thresh*filtered.max()] = 0.
# assume all value below min*threshold can not be maxima, so put all to zero
# filtered[filtered < filtered.min() + thresh * np.abs(filtered.min())] = 0.
filtered[filtered < mean_bg + thresh] = 0.
# find coordinates of local maxima
#print(int(0.5 * min_size))
max_idx_2d_small = peak_local_max(filtered, min_distance=0)
max_idx_2d_large = peak_local_max(filtered, min_distance=1)
x_list_small, y_list_small = max_idx_2d_small[:, 1], max_idx_2d_small[:, 0]
x_list_large, y_list_large = max_idx_2d_large[:, 1], max_idx_2d_large[:, 0]
im_center_x, im_center_y = len(image) / 2., len(image) / 2.
if (center_x == None) & (center_y == None):
R = np.sqrt((x_list_large - im_center_x)**2 + (y_list_large - im_center_y)**2)
new_center_x, new_center_y = x_list_large[R < gal_rad], y_list_large[R < gal_rad]
if (len(new_center_x) > 1) and (len(x_list_large[R == R.min()]) ==1 ):
new_center_x, new_center_y = x_list_large[R == R.min()], y_list_large[R == R.min()]
elif (len(new_center_x) > 1) and (len(x_list_large[R == R.min()]) > 1 ):
new_center_x, new_center_y = im_center_x, im_center_y
elif len(new_center_x) == 0:
new_center_x, new_center_y = im_center_x, im_center_y
else:
new_center_x, new_center_y = center_x,center_y
print(new_center_x,new_center_y)
R_small = np.sqrt((x_list_small - new_center_x)**2 + (y_list_small - new_center_y)**2)
R_large = np.sqrt((x_list_large - new_center_x)**2 + (y_list_large - new_center_y)**2)
x_sats, y_sats = x_list_small[R_small > lens_rad], y_list_small[R_small > lens_rad]
if many_sources:
x_lens, y_lens = deepcopy(x_list_small), deepcopy(y_list_small)
else:
x_lens, y_lens = deepcopy(x_list_large), deepcopy(y_list_large)
# x_lens, y_lens = x_list_small[R_small <= lens_rad], y_list_small[R_small <= lens_rad]
if (len(x_lens) == 0) & (len(y_lens) == 0):
x_lens = [0,15]
y_lens = [0,15]
sources = QTable([x_lens, y_lens],names={'x_local_peak','y_local_peak'})
# show maxima on image for debug
if show_locations:
# fig = plt.figure(figsize=(4, 4))
#plt.imshow(image, origin='lower', cmap=cmap_flux, norm=LogNorm(1e-2))
f, axes = plt.subplots(1, 5, figsize=(20,5), sharex=False, sharey=False)
# plt.figure(figsize = (8,8))
# plt.subplot(1,2,1)
axes[0].imshow(image, origin='lower', norm=SymLogNorm(5))
axes[0].set_title('Image')
axes[0].set_axis_off()
axes[1].imshow(LoG, origin='lower', norm=SymLogNorm(5))
axes[1].set_title('LoG Filtered Image')
axes[1].set_axis_off()
# plt.subplot(1,2,2)
axes[2].imshow(filtered, origin='lower', norm=SymLogNorm(5))
axes[2].set_title('Final Filtered Image')
axes[2].set_axis_off()
axes[3].imshow(image, origin='lower', norm=SymLogNorm(5))
for i in range(len(x_lens)):
axes[3].scatter([x_lens[i]], [y_lens[i]], c='red', s=60, marker='+')
axes[3].set_title('lens detections')
axes[3].set_axis_off()
axes[4].imshow(image, origin='lower', norm=SymLogNorm(5))
for i in range(len(x_sats)):
axes[4].scatter([x_sats[i]], [y_sats[i]], c='red', s=60, marker='+')
# plt.annotate(i+1, (x_list[i], y_list[i]), color='black')
# for i in range(len(x_mask)):
# plt.scatter([x_mask[i]], [y_mask[i]], c='red', s=100, marker='*')
# plt.annotate(i+1, (x_mask[i], y_mask[i]), color='red')
axes[4].scatter(new_center_x, new_center_y,c='red', s=100, marker='*')
draw_lens_circle = Circle((new_center_x, new_center_y),lens_rad ,fill=False)
draw_gal_circle = Circle((new_center_x, new_center_y),gal_rad, fill = False)
# plt.gcf().gca().add_artist(draw_lens_circle)
# plt.gcf().gca().add_artist(draw_gal_circle)
axes[4].add_patch(draw_lens_circle)
axes[4].add_patch(draw_gal_circle)
axes[4].set_title('Detected Components: \n Center = ({:.1f},{:.1f}) \n r = {:.3f}'.format(new_center_x, new_center_y,lens_rad_arcsec))
axes[4].text(1, 1, "detected components", color='red')
axes[4].set_axis_off()
if title != None:
f.suptitle(title, fontsize = 15)
# plt.show()
return (x_sats, y_sats), (new_center_x, new_center_y), sources
def estimate_radius_stat(image,deltaPix,center_x = None,center_y = None,show_plot=False, name = None):
corners = np.zeros([4,5,5])
corners[0] = image[0:5,0:5]
corners[1] = image[-5:,0:5]
corners[2] = image[0:5,-5:]
corners[3] = image[-5:,-5:]
means, medians, stds = [],[],[]