Skip to content

Commit

Permalink
Merge 12b401e into 1214b26
Browse files Browse the repository at this point in the history
  • Loading branch information
@gvernard
gvernard committed Jul 20, 2023
2 parents 1214b26 + 12b401e commit de31ab3
Show file tree
Hide file tree
Showing 12 changed files with 550 additions and 106 deletions.
197 changes: 111 additions & 86 deletions coolest/api/analysis.py
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
__author__ = 'aymgal', 'mattgomer'
__author__ = 'aymgal', 'mattgomer', 'gvernard'

import numpy as np
from astropy.coordinates import SkyCoord
Expand All @@ -11,6 +11,11 @@
class Analysis(object):
"""Handles computation of model-independent quantities
and other analysis computations.
NOTE: Except for methods that do have a `coordinates` keyword argument,
the grid used to performed the computations will always be the one corresponding
to the instrument / observation field-of-view, with resolution controlled by the
`supersampling` keyword argument below.
Parameters
----------
Expand All @@ -20,7 +25,7 @@ class Analysis(object):
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
that defines the grid on which computations are performed, by default 1.
"""

def __init__(self, coolest_object, coolest_directory, supersampling=1):
Expand Down Expand Up @@ -212,94 +217,12 @@ def effective_radial_slope(self, r_eval=None, center=None, r_vec=np.linspace(0,
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 two_point_correlation(self,Nbins=100,rmax=None,coordinates=None,**kwargs_selection):
"""
The two point correlation function can be obtained from the covariance matrix of an image and the distances between its pixels.
By binning the covariance matrix entries in distance (or radial) bins, one can obtain the 1D correlation function.
There are two ways to obtain the covariance matrix:
1) it is equivalent to the inverse Fourier transform of the power spectrum, and
2) by calculating explicitly the covariance between any two pixels
Here we use the first way.
Parameters
----------
Nbins : int, optional
The number of radial bins to use for converting the 2D covariance matrix into a 1D correlation function.
rmax : float, optional
A value for the maximum extent of the radial bins. If none is given then it is equal to half the diagonal of the provided image.
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
Returns
-------
(array, array, array)
The location, value, and uncertainty of the 1D bins of the two-point correlation function.
The location (radius/distance) is in the same units as the coordinates.
"""
if kwargs_selection is None:
kwargs_selection = {}
light_model = ComposableLightModel(self.coolest, self.coolest_dir, **kwargs_selection)
if coordinates is None:
x, y = self.coordinates.pixel_coordinates
else:
x, y = coordinates.pixel_coordinates
if rmax is None:
rmax = math.hypot(x[0,0]-x[-1,-1],y[0,0]-y[-1,-1])/2.0
light_image = light_model.evaluate_surface_brightness(x, y)
light_image[np.isnan(light_image)] = 0.


# Fourier transform image
fouriertf = np.fft.fft2(light_image,norm="ortho")
# Power spectrum (the square of the signal)
absval2 = fouriertf.real**2 + fouriertf.imag**2
# Covariance matrix (the inverse fourier transform of the power spectrum)
complex_cov = np.fft.fftshift(np.fft.ifft2(absval2,norm="ortho"))
cov = complex_cov.real

## Write the covariance matrix
#newhdu = fits.PrimaryHDU(corr)
#newhdu.writeto('matrix_strue.fits',overwrite=True)


# Bin the 2D covariance matrix into radial bins
rmin = 0.0
dr = (rmax-rmin)/Nbins
bins = np.arange(rmin,rmax,dr)
vals = []
for i in range(0,len(bins)):
vals[i] = []

# Ni = Nj = the total number of pixels in the image
Ni = cov.shape[1]
Nj = cov.shape[0]
for i in range(0,Ni):
for j in range(0,Nj):
r = math.hypot((j-Nx/2.0)*dpix,(i-Ny/2.0)*dpix)
if r < rmax and i != j:
index = int(math.floor(r/dr))
vals[index].append(cov[i][j])

means = np.zeros(len(bins))
sdevs = np.zeros(len(bins))
for i in range(0,len(bins)):
if len(vals[i]) > 0:
means[i] = np.mean(vals[i])
sdevs[i] = np.std(vals[i])

return bins,means,sdevs




def effective_radius_light(self, outer_radius=10, center=None, coordinates=None,
initial_guess=1, initial_delta_pix=10,
n_iter=10, **kwargs_selection):
"""
"""Computes the effective radius of the 2D surface brightness profile,
based on a definition similar to the half-light radius.
Parameters
----------
Expand Down Expand Up @@ -395,3 +318,105 @@ def find_nearest(self, array, value):
idx = (np.abs(array - value)).argmin()
return array[idx]


def two_point_correlation(self, Nbins=100, rmax=None, normalize=False,
use_profile_coordinates=True, coordinates=None,
**kwargs_selection):
"""
The two point correlation function can be obtained from the covariance matrix of an image and the distances between its pixels.
By binning the covariance matrix entries in distance (or radial) bins, one can obtain the 1D correlation function.
There are two ways to obtain the covariance matrix:
1) it is equivalent to the inverse Fourier transform of the power spectrum, and
2) by calculating explicitly the covariance between any two pixels
Here we use the first way.
Parameters
----------
Nbins : int, optional
The number of radial bins to use for converting the 2D covariance matrix into a 1D correlation function.
rmax : float, optional
A value for the maximum extent of the radial bins. If none is given then it is equal to half the diagonal of the provided image.
normalize : bool, optional
Normalize the given image by its maximum. Defualt is false.
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
use_profile_coordinates : bool, optional
If True and `coordinates=None`, uses the coordinates attached to the light profile, if available. Default is True.
Returns
-------
(array, array, array)
The location, value, and uncertainty of the 1D bins of the two-point correlation function.
The location (radius/distance) is in the same units as the coordinates.
"""
if kwargs_selection is None:
kwargs_selection = {}
if coordinates is None:
coordinates = self.coordinates

light_model = ComposableLightModel(self.coolest, self.coolest_dir, **kwargs_selection)

if use_profile_coordinates is True and coordinates is None:
light_image, _, coordinates = light_model.surface_brightness(return_extra=True)
if coordinates is None:
# can be known if e.g. the underlying light profile is not pixelated
raise ValueError("Light profile does not have any coordinates grid attached to it.")
else:
if coordinates is None:
coordinates = self.coordinates
x, y = coordinates.pixel_coordinates
light_image = light_model.evaluate_surface_brightness(x, y)

light_image = np.nan_to_num(light_image, nan=0.)
extent = coordinates.extent
dpix = coordinates.pixel_size

if rmax is None:
rmax = math.hypot(extent[0]-extent[1],extent[2]-extent[3])/2.0

if normalize:
max_image = np.amax(light_image)
light_image = np.divide(light_image,max_image)

# Fourier transform image
fouriertf = np.fft.fft2(light_image,norm="ortho")
# Power spectrum (the square of the signal)
absval2 = fouriertf.real**2 + fouriertf.imag**2
# Covariance matrix (the inverse fourier transform of the power spectrum)
complex_cov = np.fft.fftshift(np.fft.ifft2(absval2,norm="ortho"))
cov = complex_cov.real

## Write the covariance matrix
#newhdu = fits.PrimaryHDU(corr)
#newhdu.writeto('matrix_strue.fits',overwrite=True)


# Bin the 2D covariance matrix into radial bins
rmin = 0.0
dr = (rmax-rmin)/Nbins
bins = np.arange(rmin,rmax,dr)
vals = []
for i in range(0,len(bins)):
vals.append( [] )

# Ni = Nj = the total number of pixels in the image
Ni = cov.shape[1]
Nj = cov.shape[0]
for i in range(0,Ni):
for j in range(0,Nj):
r = math.hypot((j-Nj/2.0)*dpix,(i-Ni/2.0)*dpix)
if r < rmax and i != j:
index = int(math.floor(r/dr))
vals[index].append(cov[i][j])

means = np.zeros(len(bins))
sdevs = np.zeros(len(bins))
for i in range(0,len(bins)):
if len(vals[i]) > 0:
means[i] = np.mean(vals[i])
sdevs[i] = np.std(vals[i])

return bins,means,sdevs


12 changes: 7 additions & 5 deletions coolest/template/json.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,14 +104,14 @@ def load(self, verbose=False):
COOLEST object that corresponds to the JSON template
"""
try:
content = self.load_jsonpickle()
instance = self.load_jsonpickle()
assert isinstance(instance, COOLEST)
except Exception as e:
if verbose is True:
print(f"Failed reading '{self._api_suffix}' template with jsonpickle, "
f"now trying with the pure json template (original error: {e})")
content = self.load_simple()
assert isinstance(content, COOLEST)
return content
instance = self.load_simple()
return instance

def load_simple(self, as_object=True):
"""Read the JSON template file and build up the corresponding COOLEST object.
Expand Down Expand Up @@ -173,7 +173,7 @@ def _json_to_coolest(self, json_content):
# LENSING ENTITIES {GALAXY, MASSFIELDS}
lensing_entities = self._setup_lensing_entities(c['lensing_entities'])

# CSOMOLOGY
# COSMOLOGY
cosmology = self._setup_cosmology(c['cosmology'])

# COORDINATES
Expand Down Expand Up @@ -206,6 +206,8 @@ def _validate_global(coolest):
"""Performs consistency checks regarding some key properties of the COOLEST object.
For instance, it checks that the pixel size of both the observation and
the instrument are consistent.
The checks performed here are those that cannot be handled by individual
class constructors called during instantiation of the COOLEST object.
Parameters
----------
Expand Down
2 changes: 1 addition & 1 deletion docs/index.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ caption: Examples
maxdepth: 2
---
Codes compliant with COOLEST <examples/existing_interfaces>
Python notebooks <notebooks/index>
Tutorial notebooks <notebooks/index>
```

```{toctree}
Expand Down
16 changes: 2 additions & 14 deletions docs/notebooks/04-plot_lens_models.ipynb
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"id": "1a14f34a",
"metadata": {},
Expand Down Expand Up @@ -55,7 +54,6 @@
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "1106d6cc",
"metadata": {},
Expand All @@ -64,7 +62,6 @@
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "e7727434",
"metadata": {},
Expand Down Expand Up @@ -94,7 +91,6 @@
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "31a6ff1e",
"metadata": {},
Expand Down Expand Up @@ -124,7 +120,6 @@
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "c8b2eff6",
"metadata": {},
Expand Down Expand Up @@ -156,7 +151,6 @@
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "09222291",
"metadata": {},
Expand Down Expand Up @@ -192,7 +186,6 @@
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "d8945179",
"metadata": {},
Expand All @@ -213,7 +206,6 @@
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "b59e0074",
"metadata": {},
Expand Down Expand Up @@ -260,7 +252,6 @@
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "7bb72250",
"metadata": {},
Expand Down Expand Up @@ -317,7 +308,6 @@
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "6bf3b797",
"metadata": {},
Expand All @@ -326,7 +316,6 @@
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "5604ef50",
"metadata": {},
Expand Down Expand Up @@ -357,7 +346,6 @@
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "4045dbbd",
"metadata": {},
Expand Down Expand Up @@ -491,7 +479,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "coolest38",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
Expand All @@ -505,7 +493,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.17"
"version": "3.8.8"
}
},
"nbformat": 4,
Expand Down
Loading

0 comments on commit de31ab3

Please sign in to comment.