4993920423

Build # 4993920423

Build Type

github

Committed by GitHub

Commit Message

Merge efd26c040 into d1de71ffa

Pull Request Pull Request #34: Preparation for JOSS submission

Run Details

184 of 184 new or added lines in 28 files covered. (100.0%)

1071 of 2324 relevant lines covered (46.08%)

0.46 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/coolest/api/analysis.py

__author__ = 'aymgal'

import numpy as np
from astropy.coordinates import SkyCoord
import warnings

from coolest.api.composable_models import *
from coolest.api import util


class Analysis(object):
    """Handles computation of model-independent quantities 
    and other analysis computations.

    Parameters
    ----------
    coolest_object : COOLEST
        COOLEST instance
    coolest_directory : str
        Directory which contains the COOLEST template and other data files
    supersampling : int, optional
        Supersampling factor (relative to the instrument pixel size)
        that defines the grid on which computations are performed, by default 1
    """

    def __init__(self, coolest_object, coolest_directory, supersampling=1):
        self.coolest = coolest_object
        self.coolest_dir = coolest_directory
        base_coordinates = util.get_coordinates(self.coolest)
        if supersampling > 1:
            self.coordinates = base_coordinates.create_new_coordinates(pixel_scale_factor=1./supersampling)
        else:
            self.coordinates = base_coordinates

    def effective_einstein_radius(self, center=None, initial_guess=1, initial_delta_pix=10, 
                                  n_iter=5, return_accuracy=False, **kwargs_selection):
        """Calculates Einstein radius for a kappa grid starting from an initial guess with large step size and zeroing in from there.
        Uses the grid from the create_kappa_image which is built from the coolest file.

        Parameters
        ----------
        center : (float, float), optional
            (x, y)-coordinates of the center from which to calculate Einstein radius; if None, use the value from create_kappa_image, by default None
        initial_guess : int, optional
            initial guess for Einstein radius, by default 1
        initial_delta_pix : int, optional
            initial step size before shrinking in future iterations, by default 10
        n_iter : int, optional
            number of iterations, by default 5
        return_accuracy : bool, optional
            if True, return estimate of accuracy as well as R_Ein value, by default False

        Returns
        -------
        float, or (float, float) if return_accuracy is True
            Effective Einstein radius

        Raises
        ------
        Warning
            If the algorithm is running for more than 100 loops.
        """
        if kwargs_selection is None:
            kwargs_selection = {}
        mass_model = ComposableMassModel(self.coolest, self.coolest_dir, **kwargs_selection)

        # get an image of the convergence
        x, y = self.coordinates.pixel_coordinates
        kappa_image = mass_model.evaluate_convergence(x, y)

        # select a center
        if center is None:
            center_x, center_y = mass_model.estimate_center()
        else:
            center_x, center_y = center

        def loop_until_overshoot(r_Ein, delta, direction, runningtotal, total_pix):
            """
            this subfunction iteratively adjusts the mask radius by delta either inward or outward until the sign flips on mean_kappa-area
            """
            loopcount=0
            if r_Ein==np.nan: #return nan if given nan (needed when I put this in a loop below)
                return np.nan, np.nan, 0, runningtotal, total_pix
            while (runningtotal-total_pix)*direction>0:
                if loopcount > 100:
                    raise Warning('Stuck in very long (possibly infinite) loop')
                    break
                if direction > 0:
                    mask=only_between_r1r2(r_Ein, r_Ein+delta,r_grid)
                else:
                    mask=only_between_r1r2(r_Ein-delta, r_Ein, r_grid)
                kappa_annulus=np.sum(kappa_image[mask])
                added_pix=np.sum(mask)
                runningtotal=runningtotal+kappa_annulus*direction
                total_pix=total_pix+added_pix*direction
                r_Ein=r_Ein+delta*direction

                if r_Ein < min_guess:
                    #Either we overshot into the center region between pixels (rare), or kappa is subcritical.
                    #check starting from zero
                    mask=only_between_r1r2(0,min_guess,r_grid)
                    runningtotal=np.sum(kappa_image[mask])
                    total_pix=np.sum(mask)
                    if runningtotal-total_pix < 0:
                        print('WARNING: kappa is sub-critical, Einstein radius undefined.')
                        return np.nan, np.nan, 0, runningtotal, total_pix
                loopcount+=1
            return r_Ein, delta, direction, runningtotal, total_pix

        def only_between_r1r2(r1, r2, r_grid):
            if r1 > r2:
                raise ValueError('r2 must be greater than r1')
            return (r_grid >= r1) & (r_grid < r2)

        #some initialization: if there is an even number of pixels we don't want to be stuck in the center between 4 pixels.
        grid_res=np.abs(x[0,0]-x[0,1])
        initial_delta=grid_res*initial_delta_pix #default inital step size is 10 pixels
        min_guess=grid_res*np.sqrt(2)+1e-6
        initial_guess=max(initial_guess,min_guess) #initial guess must include at least one pixel

        r_grid=np.sqrt((x-center_x)**2+(y-center_y)**2)
        mask=only_between_r1r2(0,initial_guess,r_grid)
        runningtotal=np.sum(kappa_image[mask])
        total_pix=np.sum(mask)
        if runningtotal > total_pix: #move outward
            direction=1
        elif runningtotal < total_pix: #move inward
            direction=-1
        else:
            return initial_guess
        r_Ein=initial_guess
        delta=initial_delta

        for n in range(n_iter):
            #overshoots, turn around and backtrack at higher precision
            r_Ein, delta, direction, runningtotal, total_pix = loop_until_overshoot(r_Ein, delta, direction, runningtotal, total_pix)
            direction=direction*-1
            delta=delta/2
        accuracy=grid_res/2 #after testing, accuracy is about grid_res/2
        if np.isnan(r_Ein):
            accuracy=np.nan

        if return_accuracy:
            return r_Ein, accuracy
        else:
            return r_Ein

    def kappa_1d_profile(self, center=None, r_vec=np.linspace(0, 10, 100),**kwargs_selection):
        """Calculates 1D profile using the kappa grid.

        Parameters
        ----------
        center : (float, float), optional
            (x, y)-coordinates of the center from which to calculate; if None, use the value from , by default None
        r_vec : _type_, optional
            range of radii over which to calculate the 1D profile, by default np.linspace(0, 10, 100)

        Returns
        -------
        (array, array)
            kappa values and associated radius values
        """
        if kwargs_selection is None:
            kwargs_selection = {}
        mass_model = ComposableMassModel(self.coolest, self.coolest_dir, **kwargs_selection)
        x, y = self.coordinates.pixel_coordinates
        kappa_image = mass_model.evaluate_convergence(x, y)
        if center is None:
            center_x, center_y = mass_model.estimate_center()
        else:
            center_x, center_y = center
        r_grid=np.sqrt((x-center_x)**2+(y-center_y)**2)
        kappa_profile=np.zeros_like(r_vec)
        for r_i,r in enumerate(r_vec):
            if r_i==0:
                r_min=0
            else:
                r_min=r_vec[r_i-1]
            in_radial_bin=(r_grid > r_min) & (r_grid <= r)
            kappa_profile[r_i]=np.mean(kappa_image[in_radial_bin])
        return kappa_profile, r_vec

    def effective_radial_slope(self, r_eval=None, center=None, r_vec=np.linspace(0, 10, 100),**kwargs_selection):
        """Numerically calculates slope of the kappa profile. Because this is defined on a grid, it is not as accurate or robust as
        an analytical calculation. 

        Parameters
        ----------
        r_eval : float, optional
            radius at which to return a single value of the slope (e.g. Einstein radius). If None (default), returns slope for all values in r_vec, by default None
        center : (float, float), optional
            (x, y)-coordinates of the center from which to calculate; if None, use the value from create_kappa_image, by default None
        r_vec : array-like, optional
            range of radii over which to calculate the 1D profile, by default np.linspace(0, 10, 100)

        Returns
        -------
        float
            Effective slope
        """
        kappa_profile, r_vec =self.kappa_1d_profile(center=center, r_vec=r_vec, **kwargs_selection)
        rise=np.log10(kappa_profile[:-1])-np.log10(kappa_profile[1:])
        run=np.log10(r_vec[:-1])-np.log10(r_vec[1:])
        slope=np.append(0,rise/run)#add a zero at r=0 so that slope as same size as r_vec
        if r_eval==None:
            return slope 
        else:
            closest_r = self.find_nearest(r_vec,r_eval) #just takes closest r. Could rebuild it to interpolate.
            return slope[r_vec==closest_r]
            
    def effective_radius_light(self, outer_radius=10, center=None, coordinates=None,
                               initial_guess=1, initial_delta_pix=10, 
                               n_iter=5, **kwargs_selection):
        """        


        Parameters
        ----------
        outer_radius : int, optional
            outer limit of integration within which half the light is calculated to estimate the effective radius, by default 10
        center : (float, float), optional
            (x, y)-coordinates of the center from which to calculate Einstein radius; if None, use the value from create_kappa_image, by default None
        coordinates : Coordinates, optional
            Instance of a Coordinates object to be used for the computation.
            If None, will use an instance based on the Instrument, by default None
        initial_guess : int, optional
            initial guess for effective radius, by default 1
        initial_delta_pix : int, optional
            initial step size before shrinking in future iterations, by default 10
        n_iter : int, optional
            number of iterations, by default 5

        Returns
        -------
        float
            Effective radius

        Raises
        ------
        Warning
            If integration loop exceeds outer bound before convergence.
        """
        if kwargs_selection is None:
            kwargs_selection = {}
        light_model = ComposableLightModel(self.coolest, self.coolest_dir, **kwargs_selection)
        # get an image of the convergence
        if coordinates is None:
            x, y = self.coordinates.pixel_coordinates
        else:
            x, y = coordinates.pixel_coordinates
        light_image = light_model.evaluate_surface_brightness(x, y)
        light_image[np.isnan(light_image)] = 0.

        # select a center
        if center is None:
            center_x, center_y = light_model.estimate_center()
        else:
            center_x, center_y = center
        
        #if limit of integration exceeds FoV, raise warning
        x_FoV=self.coolest.observation.pixels.field_of_view_x
        y_FoV=self.coolest.observation.pixels.field_of_view_y
        out_of_FoV=False
        if center_x - outer_radius < x_FoV[0] or center_x + outer_radius > x_FoV[1]:
            out_of_FoV=True
        if center_y - outer_radius < y_FoV[0] or center_y + outer_radius > y_FoV[1]:
            out_of_FoV=True
        if out_of_FoV is True:
            warnings.warn("Warning: Outer limit of integration exceeds FoV; effective radius may not be accurate.")
            
        #initialize
        grid_res=np.abs(x[0,0]-x[0,1])
        initial_delta=grid_res*initial_delta_pix #default inital step size is 10 pixels
        r_grid=np.sqrt((x-center_x)**2+(y-center_y)**2)
        total_light=np.sum(light_image[r_grid<outer_radius])
        cumulative_light=np.sum(light_image[r_grid<initial_guess])
        if cumulative_light < total_light/2: #move outward
            direction=1
        elif cumulative_light > total_light/2: #move inward
            direction=-1
        else:
            return initial_guess
        r_eff=initial_guess
        delta=initial_delta
        loopcount=0
        
        for n in range(n_iter): #overshoots, turn around and backtrack at higher precision
            while (total_light/2-cumulative_light)*direction>0: 
                if loopcount > 100:
                    raise Warning('Stuck in very long (possibly infinite) loop')
                    break
                r_eff=r_eff+delta*direction
                cumulative_light=np.sum(light_image[r_grid<r_eff])
                loopcount+=1
            direction=direction*-1
            delta=delta/2
            
        return r_eff

    
    def find_nearest(self, array, value):
        """subfunction to find nearest closest element in array to value"""
        array = np.asarray(array)
        idx = (np.abs(array - value)).argmin()
        return array[idx]
    

1	__author__ = 'aymgal'	×
2
3	import numpy as np	×
4	from astropy.coordinates import SkyCoord	×
5	import warnings	×
6
7	from coolest.api.composable_models import *	×
8	from coolest.api import util	×
9
10
11	class Analysis(object):	×
12	"""Handles computation of model-independent quantities
13	and other analysis computations.
14
15	Parameters
16	----------
17	coolest_object : COOLEST
18	COOLEST instance
19	coolest_directory : str
20	Directory which contains the COOLEST template and other data files
21	supersampling : int, optional
22	Supersampling factor (relative to the instrument pixel size)
23	that defines the grid on which computations are performed, by default 1
24	"""
25
26	def __init__(self, coolest_object, coolest_directory, supersampling=1):	×
27	self.coolest = coolest_object	×
28	self.coolest_dir = coolest_directory	×
29	base_coordinates = util.get_coordinates(self.coolest)	×
30	if supersampling > 1:	×
31	self.coordinates = base_coordinates.create_new_coordinates(pixel_scale_factor=1./supersampling)	×
32	else:
33	self.coordinates = base_coordinates	×
34
35	def effective_einstein_radius(self, center=None, initial_guess=1, initial_delta_pix=10,	×
36	n_iter=5, return_accuracy=False, **kwargs_selection):
37	"""Calculates Einstein radius for a kappa grid starting from an initial guess with large step size and zeroing in from there.
38	Uses the grid from the create_kappa_image which is built from the coolest file.
39
40	Parameters
41	----------
42	center : (float, float), optional
43	(x, y)-coordinates of the center from which to calculate Einstein radius; if None, use the value from create_kappa_image, by default None
44	initial_guess : int, optional
45	initial guess for Einstein radius, by default 1
46	initial_delta_pix : int, optional
47	initial step size before shrinking in future iterations, by default 10
48	n_iter : int, optional
49	number of iterations, by default 5
50	return_accuracy : bool, optional
51	if True, return estimate of accuracy as well as R_Ein value, by default False
52
53	Returns
54	-------
55	float, or (float, float) if return_accuracy is True
56	Effective Einstein radius
57
58	Raises
59	------
60	Warning
61	If the algorithm is running for more than 100 loops.
62	"""
63	if kwargs_selection is None:	×
64	kwargs_selection = {}	×
65	mass_model = ComposableMassModel(self.coolest, self.coolest_dir, **kwargs_selection)	×
66
67	# get an image of the convergence
68	x, y = self.coordinates.pixel_coordinates	×
69	kappa_image = mass_model.evaluate_convergence(x, y)	×
70
71	# select a center
72	if center is None:	×
73	center_x, center_y = mass_model.estimate_center()	×
74	else:
75	center_x, center_y = center	×
76
77	def loop_until_overshoot(r_Ein, delta, direction, runningtotal, total_pix):	×
78	"""
79	this subfunction iteratively adjusts the mask radius by delta either inward or outward until the sign flips on mean_kappa-area
80	"""
81	loopcount=0	×
82	if r_Ein==np.nan: #return nan if given nan (needed when I put this in a loop below)	×
83	return np.nan, np.nan, 0, runningtotal, total_pix	×
84	while (runningtotal-total_pix)*direction>0:	×
85	if loopcount > 100:	×
86	raise Warning('Stuck in very long (possibly infinite) loop')	×
87	break	×
88	if direction > 0:	×
89	mask=only_between_r1r2(r_Ein, r_Ein+delta,r_grid)	×
90	else:
91	mask=only_between_r1r2(r_Ein-delta, r_Ein, r_grid)	×
92	kappa_annulus=np.sum(kappa_image[mask])	×
93	added_pix=np.sum(mask)	×
94	runningtotal=runningtotal+kappa_annulus*direction	×
95	total_pix=total_pix+added_pix*direction	×
96	r_Ein=r_Ein+delta*direction	×
97
98	if r_Ein < min_guess:	×
99	#Either we overshot into the center region between pixels (rare), or kappa is subcritical.
100	#check starting from zero
101	mask=only_between_r1r2(0,min_guess,r_grid)	×
102	runningtotal=np.sum(kappa_image[mask])	×
103	total_pix=np.sum(mask)	×
104	if runningtotal-total_pix < 0:	×
105	print('WARNING: kappa is sub-critical, Einstein radius undefined.')	×
106	return np.nan, np.nan, 0, runningtotal, total_pix	×
107	loopcount+=1	×
108	return r_Ein, delta, direction, runningtotal, total_pix	×
109
110	def only_between_r1r2(r1, r2, r_grid):	×
111	if r1 > r2:	×
112	raise ValueError('r2 must be greater than r1')	×
113	return (r_grid >= r1) & (r_grid < r2)	×
114
115	#some initialization: if there is an even number of pixels we don't want to be stuck in the center between 4 pixels.
116	grid_res=np.abs(x[0,0]-x[0,1])	×
117	initial_delta=grid_res*initial_delta_pix #default inital step size is 10 pixels	×
118	min_guess=grid_res*np.sqrt(2)+1e-6	×
119	initial_guess=max(initial_guess,min_guess) #initial guess must include at least one pixel	×
120
121	r_grid=np.sqrt((x-center_x)2+(y-center_y)2)	×
122	mask=only_between_r1r2(0,initial_guess,r_grid)	×
123	runningtotal=np.sum(kappa_image[mask])	×
124	total_pix=np.sum(mask)	×
125	if runningtotal > total_pix: #move outward	×
126	direction=1	×
127	elif runningtotal < total_pix: #move inward	×
128	direction=-1	×
129	else:
130	return initial_guess	×
131	r_Ein=initial_guess	×
132	delta=initial_delta	×
133
134	for n in range(n_iter):	×
135	#overshoots, turn around and backtrack at higher precision
136	r_Ein, delta, direction, runningtotal, total_pix = loop_until_overshoot(r_Ein, delta, direction, runningtotal, total_pix)	×
137	direction=direction*-1	×
138	delta=delta/2	×
139	accuracy=grid_res/2 #after testing, accuracy is about grid_res/2	×
140	if np.isnan(r_Ein):	×
141	accuracy=np.nan	×
142
143	if return_accuracy:	×
144	return r_Ein, accuracy	×
145	else:
146	return r_Ein	×
147
148	def kappa_1d_profile(self, center=None, r_vec=np.linspace(0, 10, 100),**kwargs_selection):	×
149	"""Calculates 1D profile using the kappa grid.
150
151	Parameters
152	----------
153	center : (float, float), optional
154	(x, y)-coordinates of the center from which to calculate; if None, use the value from , by default None
155	r_vec : _type_, optional
156	range of radii over which to calculate the 1D profile, by default np.linspace(0, 10, 100)
157
158	Returns
159	-------
160	(array, array)
161	kappa values and associated radius values
162	"""
163	if kwargs_selection is None:	×
164	kwargs_selection = {}	×
165	mass_model = ComposableMassModel(self.coolest, self.coolest_dir, **kwargs_selection)	×
166	x, y = self.coordinates.pixel_coordinates	×
167	kappa_image = mass_model.evaluate_convergence(x, y)	×
168	if center is None:	×
169	center_x, center_y = mass_model.estimate_center()	×
170	else:
171	center_x, center_y = center	×
172	r_grid=np.sqrt((x-center_x)2+(y-center_y)2)	×
173	kappa_profile=np.zeros_like(r_vec)	×
174	for r_i,r in enumerate(r_vec):	×
175	if r_i==0:	×
176	r_min=0	×
177	else:
178	r_min=r_vec[r_i-1]	×
179	in_radial_bin=(r_grid > r_min) & (r_grid <= r)	×
180	kappa_profile[r_i]=np.mean(kappa_image[in_radial_bin])	×
181	return kappa_profile, r_vec	×
182
183	def effective_radial_slope(self, r_eval=None, center=None, r_vec=np.linspace(0, 10, 100),**kwargs_selection):	×
184	"""Numerically calculates slope of the kappa profile. Because this is defined on a grid, it is not as accurate or robust as
185	an analytical calculation.
186
187	Parameters
188	----------
189	r_eval : float, optional
190	radius at which to return a single value of the slope (e.g. Einstein radius). If None (default), returns slope for all values in r_vec, by default None
191	center : (float, float), optional
192	(x, y)-coordinates of the center from which to calculate; if None, use the value from create_kappa_image, by default None
193	r_vec : array-like, optional
194	range of radii over which to calculate the 1D profile, by default np.linspace(0, 10, 100)
195
196	Returns
197	-------
198	float
199	Effective slope
200	"""
201	kappa_profile, r_vec =self.kappa_1d_profile(center=center, r_vec=r_vec, **kwargs_selection)	×
202	rise=np.log10(kappa_profile[:-1])-np.log10(kappa_profile[1:])	×
203	run=np.log10(r_vec[:-1])-np.log10(r_vec[1:])	×
204	slope=np.append(0,rise/run)#add a zero at r=0 so that slope as same size as r_vec	×
205	if r_eval==None:	×
206	return slope	×
207	else:
208	closest_r = self.find_nearest(r_vec,r_eval) #just takes closest r. Could rebuild it to interpolate.	×
209	return slope[r_vec==closest_r]	×
210
211	def effective_radius_light(self, outer_radius=10, center=None, coordinates=None,	×
212	initial_guess=1, initial_delta_pix=10,
213	n_iter=5, **kwargs_selection):
214	"""
215
216
217	Parameters
218	----------
219	outer_radius : int, optional
220	outer limit of integration within which half the light is calculated to estimate the effective radius, by default 10
221	center : (float, float), optional
222	(x, y)-coordinates of the center from which to calculate Einstein radius; if None, use the value from create_kappa_image, by default None
223	coordinates : Coordinates, optional
224	Instance of a Coordinates object to be used for the computation.
225	If None, will use an instance based on the Instrument, by default None
226	initial_guess : int, optional
227	initial guess for effective radius, by default 1
228	initial_delta_pix : int, optional
229	initial step size before shrinking in future iterations, by default 10
230	n_iter : int, optional
231	number of iterations, by default 5
232
233	Returns
234	-------
235	float
236	Effective radius
237
238	Raises
239	------
240	Warning
241	If integration loop exceeds outer bound before convergence.
242	"""
243	if kwargs_selection is None:	×
244	kwargs_selection = {}	×
245	light_model = ComposableLightModel(self.coolest, self.coolest_dir, **kwargs_selection)	×
246	# get an image of the convergence
247	if coordinates is None:	×
248	x, y = self.coordinates.pixel_coordinates	×
249	else:
250	x, y = coordinates.pixel_coordinates	×
251	light_image = light_model.evaluate_surface_brightness(x, y)	×
252	light_image[np.isnan(light_image)] = 0.	×
253
254	# select a center
255	if center is None:	×
256	center_x, center_y = light_model.estimate_center()	×
257	else:
258	center_x, center_y = center	×
259
260	#if limit of integration exceeds FoV, raise warning
261	x_FoV=self.coolest.observation.pixels.field_of_view_x	×
262	y_FoV=self.coolest.observation.pixels.field_of_view_y	×
263	out_of_FoV=False	×
264	if center_x - outer_radius < x_FoV[0] or center_x + outer_radius > x_FoV[1]:	×
265	out_of_FoV=True	×
266	if center_y - outer_radius < y_FoV[0] or center_y + outer_radius > y_FoV[1]:	×
267	out_of_FoV=True	×
268	if out_of_FoV is True:	×
269	warnings.warn("Warning: Outer limit of integration exceeds FoV; effective radius may not be accurate.")	×
270
271	#initialize
272	grid_res=np.abs(x[0,0]-x[0,1])	×
273	initial_delta=grid_res*initial_delta_pix #default inital step size is 10 pixels	×
274	r_grid=np.sqrt((x-center_x)2+(y-center_y)2)	×
275	total_light=np.sum(light_image[r_grid<outer_radius])	×
276	cumulative_light=np.sum(light_image[r_grid<initial_guess])	×
277	if cumulative_light < total_light/2: #move outward	×
278	direction=1	×
279	elif cumulative_light > total_light/2: #move inward	×
280	direction=-1	×
281	else:
282	return initial_guess	×
283	r_eff=initial_guess	×
284	delta=initial_delta	×
285	loopcount=0	×
286
287	for n in range(n_iter): #overshoots, turn around and backtrack at higher precision	×
288	while (total_light/2-cumulative_light)*direction>0:	×
289	if loopcount > 100:	×
290	raise Warning('Stuck in very long (possibly infinite) loop')	×
291	break	×
292	r_eff=r_eff+delta*direction	×
293	cumulative_light=np.sum(light_image[r_grid<r_eff])	×
294	loopcount+=1	×
295	direction=direction*-1	×
296	delta=delta/2	×
297
298	return r_eff	×
299
300
301	def find_nearest(self, array, value):	×
302	"""subfunction to find nearest closest element in array to value"""
303	array = np.asarray(array)	×
304	idx = (np.abs(array - value)).argmin()	×
305	return array[idx]	×
306

aymgal / COOLEST / 4993920423

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous