Revert "Normalized Mean Model and PSF metrics"#70

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
christinahedges merged 1 commit into from
Dec 8, 2022
Merged

Revert "Normalized Mean Model and PSF metrics" #70

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
185 changes: 2 additions & 183 deletions src/psfmachine/machine.py
View file
Open in desktop
Original file line numberDiff line numberDiff line change
Expand Up@@ -5,7 +5,6 @@
import numpy as np
import pandas as pd
from scipy import sparse
from scipy import stats
import astropy.units as u
from tqdm import tqdm
import matplotlib.pyplot as plt
Expand DownExpand Up@@ -1057,11 +1056,8 @@ def build_shape_model(

self.psf_w = psf_w
self.psf_w_err = psf_w_err
self.normalized_shape_model = False

# We then build the same design matrix for all pixels with flux
# this non-normalized mean model is temporary and used to re-create a better
# `source_mask`
self._get_mean_model()
# remove background pixels and recreate mean model
self._update_source_mask_remove_bkg_pixels(
Expand DownExpand Up@@ -1140,8 +1136,8 @@ def _update_source_mask_remove_bkg_pixels(self, flux_cut_off=1, frame_index="mea
# (self.mean_model.max(axis=1) < 1)
# )

# create the final normalized mean model!
self._get_normalized_mean_model()
# Recreate mean model!
self._get_mean_model()

def _get_mean_model(self):
"""Convenience function to make the scene model"""
Expand All@@ -1163,183 +1159,6 @@ def _get_mean_model(self):
mean_model.eliminate_zeros()
self.mean_model = mean_model

def _get_normalized_mean_model(self, npoints=300, plot=False):
"""Renomarlize shape model to sum 1"""

# create a high resolution polar grid
r = self.source_mask.multiply(self.r).data
phi_hd = np.linspace(-np.pi, np.pi, npoints)
r_hd = np.linspace(0, r.max(), npoints)
phi_hd, r_hd = np.meshgrid(phi_hd, r_hd)

# high res DM
Ap = _make_A_polar(
phi_hd.ravel(),
r_hd.ravel(),
rmin=self.rmin,
rmax=self.rmax,
cut_r=self.cut_r,
n_r_knots=self.n_r_knots,
n_phi_knots=self.n_phi_knots,
)
# evaluate the high res model
mean_model_hd = Ap.dot(self.psf_w)
mean_model_hd[~np.isfinite(mean_model_hd)] = np.nan
mean_model_hd = mean_model_hd.reshape(phi_hd.shape)

# mask out datapoint that don't contribuite to the psf
mean_model_hd_ma = mean_model_hd.copy()
mask = mean_model_hd > -3
mean_model_hd_ma[~mask] = -np.inf
mask &= ~((r_hd > 14) & (np.gradient(mean_model_hd_ma, axis=0) > 0))
mean_model_hd_ma[~mask] = -np.inf

# double integral using trapezoidal rule
integral = np.trapz(
np.trapz(10 ** mean_model_hd_ma, r_hd[:, 0], axis=0),
phi_hd[0, :],
axis=0,
)
# renormalize weights and build new shape model
if not self.normalized_shape_model:
self.psf_w *= np.log10(integral)
self.normalized_shape_model = True
self._get_mean_model()

if plot:
fig, ax = plt.subplots(1, 2, figsize=(9, 5))
im = ax[0].scatter(
phi_hd.ravel(),
r_hd.ravel(),
c=mean_model_hd_ma.ravel(),
vmin=-3,
vmax=-1,
s=1,
label=r"$\int = $" + f"{integral:.4f}",
)
im = ax[1].scatter(
r_hd.ravel() * np.cos(phi_hd.ravel()),
r_hd.ravel() * np.sin(phi_hd.ravel()),
c=mean_model_hd_ma.ravel(),
vmin=-3,
vmax=-1,
s=1,
)
ax[0].legend()
fig.colorbar(im, ax=ax, location="bottom")
plt.show()

def get_psf_metrics(self, npoints_per_pixel=10):
"""
Computes three metrics for the PSF model:
source_psf_fraction: the amount of PSF in the data. Tells how much of a
sources is used to estimate the PSF, values are in between [0, 1].
perturbed_ratio_mean: the ratio between the mean model and perturbed model
for each source. Usefull to find when the time model affects the
mean value of the light curve.
perturbed_std: the standard deviation of the perturbed model for each
source. USeful to find when the time model introduces variability in the
light curve.

If npoints_per_pixel > 0, it creates high npoints_per_pixel shape models for each source by
dividing each pixels into a grid of [npoints_per_pixel x npoints_per_pixel]. This provides
a better estimate of `source_psf_fraction`.

Parameters
----------
npoints_per_pixel : int
Value in which each pixel axis is split to increase npoints_per_pixel. Default is
0 for no subpixel npoints_per_pixel.

"""
if npoints_per_pixel > 0:
# find from which observation (TPF) a sources comes
obs_per_pixel = self.source_mask.multiply(self.pix2obs).tocsr()
tpf_idx = []
for k in range(self.source_mask.shape[0]):
pix = obs_per_pixel[k].data
mode = stats.mode(pix)[0]
if len(mode) > 0:
tpf_idx.append(mode[0])
else:
tpf_idx.append(
[x for x, ss in enumerate(self.tpf_meta["sources"]) if k in ss][
0
]
)
tpf_idx = np.array(tpf_idx)

# get the pix coord for each source, we know how to increase resolution in
# the pixel space but not in WCS
row = self.source_mask.multiply(self.row).tocsr()
col = self.source_mask.multiply(self.column).tocsr()
mean_model_hd_sum = []
# iterating per sources avoids creating a new super large `source_mask`
# with high resolution, which a priori is hard
for k in range(self.nsources):
# find row, col combo for each source
row_ = row[k].data
col_ = col[k].data
colhd, rowhd = [], []
# pixels are divided into `resolution` - 1 subpixels
for c, r in zip(col_, row_):
x = np.linspace(c - 0.5, c + 0.5, npoints_per_pixel + 1)
y = np.linspace(r - 0.5, r + 0.5, npoints_per_pixel + 1)
x, y = np.meshgrid(x, y)
colhd.extend(x[:, :-1].ravel())
rowhd.extend(y[:-1].ravel())
colhd = np.array(colhd)
rowhd = np.array(rowhd)
# convert to ra, dec beacuse machine shape model works in sky coord
rahd, dechd = self.tpfs[tpf_idx[k]].wcs.wcs_pix2world(
colhd - self.tpfs[tpf_idx[k]].column,
rowhd - self.tpfs[tpf_idx[k]].row,
0,
)
drahd = rahd - self.sources["ra"][k]
ddechd = dechd - self.sources["dec"][k]
drahd = drahd * (u.deg)
ddechd = ddechd * (u.deg)
rhd = np.hypot(drahd, ddechd).to("arcsec").value
phihd = np.arctan2(ddechd, drahd).value
# create a high resolution DM
Ap = _make_A_polar(
phihd.ravel(),
rhd.ravel(),
rmin=self.rmin,
rmax=self.rmax,
cut_r=self.cut_r,
n_r_knots=self.n_r_knots,
n_phi_knots=self.n_phi_knots,
)
# evaluate the HD model
modelhd = 10 ** Ap.dot(self.psf_w)
# compute the model sum for source, how much of the source is in data
mean_model_hd_sum.append(
np.trapz(modelhd, dx=1 / npoints_per_pixel ** 2)
)

# get normalized psf fraction metric
self.source_psf_fraction = np.array(
mean_model_hd_sum
) # / np.nanmax(mean_model_hd_sum)
else:
self.source_psf_fraction = np.array(self.mean_model.sum(axis=1)).ravel()

# time model metrics
if hasattr(self, "P"):
perturbed_lcs = np.vstack(
[
np.array(self.perturbed_model(time_index=k).sum(axis=1)).ravel()
for k in range(self.time.shape[0])
]
)
self.perturbed_ratio_mean = (
np.nanmean(perturbed_lcs, axis=0)
/ np.array(self.mean_model.sum(axis=1)).ravel()
)
self.perturbed_std = np.nanstd(perturbed_lcs, axis=0)

def plot_shape_model(self, radius=20, frame_index="mean", bin_data=False):
"""
Diagnostic plot of shape model.
Expand Down
6 changes: 0 additions & 6 deletions src/psfmachine/tpf.py
View file
Open in desktop
Original file line numberDiff line numberDiff line change
Expand Up@@ -689,8 +689,6 @@ def load_shape_model(self, input=None, plot=False):
self.rmax = hdu[1].header["rmax"]
self.cut_r = hdu[1].header["cut_r"]
self.psf_w = hdu[1].data["psf_w"]
# read from header if weights come from a normalized model.
self.normalized_shape_model = bool(hdu[1].header.get("normalized"))
del hdu

# create mean model, but PRF shapes from FFI are in pixels! and TPFMachine
Expand DownExpand Up@@ -750,10 +748,6 @@ def save_shape_model(self, output=None):
)
# spline degree is hardcoded in `_make_A_polar` implementation.
table.header["spln_deg"] = (3, "Degree of the spline basis")
table.header["normalized"] = (
int(self.normalized_shape_model),
"Normalized weights",
)

table.writeto(output, checksum=True, overwrite=True)

Expand Down
, 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Add copy buttons to all 
 blocks\n(function() {\n function addCopyButtons() {\n document.querySelectorAll('pre code').forEach(function(codeBlock) {\n if (codeBlock.parentElement.hasAttribute('data-copy-added')) return;\n codeBlock.parentElement.setAttribute('data-copy-added', 'true');\n \n var btn = document.createElement('button');\n btn.textContent = 'Copy';\n btn.style.cssText = 'position:absolute;top:4px;right:4px;padding:2px 8px;font-size:11px;background:#4ecdc4;border:none;border-radius:4px;color:#1a1a2e;cursor:pointer;opacity:0.7;transition:opacity 0.2s;';\n btn.onmouseover = function() { this.style.opacity = '1'; };\n btn.onmouseout = function() { this.style.opacity = '0.7'; };\n btn.onclick = function() {\n navigator.clipboard.writeText(codeBlock.textContent).then(function() {\n btn.textContent = 'Copied!';\n setTimeout(function() { btn.textContent = 'Copy'; }, 1500);\n });\n };\n codeBlock.parentElement.style.position = 'relative';\n codeBlock.parentElement.appendChild(btn);\n });\n }\n \n addCopyButtons();\n \n // Re-run on dynamic content\n var observer = new MutationObserver(addCopyButtons);\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Add Copy Buttons to Code Blocks");
}
} catch(__e) { console.warn('[Userscript:Add Copy Buttons to Code Blocks]', __e); }
})();
(function(){
try {
var __m = "github.com";
var __re = new RegExp('^' + "github\\.com" + '

Skip to content
Revert "Normalized Mean Model and PSF metrics"#70

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

christinahedges
merged 1 commit into

from

Dec 8, 2022
 Merged


Revert "Normalized Mean Model and PSF metrics"
#70
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


1 commit


Select commit


File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


185 changes: 2 additions & 183 deletions
src/psfmachine/machine.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -5,7 +5,6 @@
import numpy as np
import pandas as pd
from scipy import sparse
from scipy import stats
import astropy.units as u
from tqdm import tqdm
import matplotlib.pyplot as plt
Expand DownExpand Up@@ -1057,11 +1056,8 @@ def build_shape_model(

 self.psf_w = psf_w
 self.psf_w_err = psf_w_err
 self.normalized_shape_model = False

 # We then build the same design matrix for all pixels with flux
 # this non-normalized mean model is temporary and used to re-create a better
 # `source_mask`
 self._get_mean_model()
 # remove background pixels and recreate mean model
 self._update_source_mask_remove_bkg_pixels(
Expand DownExpand Up@@ -1140,8 +1136,8 @@ def _update_source_mask_remove_bkg_pixels(self, flux_cut_off=1, frame_index="mea
 # (self.mean_model.max(axis=1) < 1)
 # )

 # create the final normalized mean model!
 self._get_normalized_mean_model()
 # Recreate mean model!
 self._get_mean_model()

 def _get_mean_model(self):
 """Convenience function to make the scene model"""
Expand All@@ -1163,183 +1159,6 @@ def _get_mean_model(self):
 mean_model.eliminate_zeros()
 self.mean_model = mean_model

 def _get_normalized_mean_model(self, npoints=300, plot=False):
 """Renomarlize shape model to sum 1"""

 # create a high resolution polar grid
 r = self.source_mask.multiply(self.r).data
 phi_hd = np.linspace(-np.pi, np.pi, npoints)
 r_hd = np.linspace(0, r.max(), npoints)
 phi_hd, r_hd = np.meshgrid(phi_hd, r_hd)

 # high res DM
 Ap = _make_A_polar(
 phi_hd.ravel(),
 r_hd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the high res model
 mean_model_hd = Ap.dot(self.psf_w)
 mean_model_hd[~np.isfinite(mean_model_hd)] = np.nan
 mean_model_hd = mean_model_hd.reshape(phi_hd.shape)

 # mask out datapoint that don't contribuite to the psf
 mean_model_hd_ma = mean_model_hd.copy()
 mask = mean_model_hd > -3
 mean_model_hd_ma[~mask] = -np.inf
 mask &= ~((r_hd > 14) & (np.gradient(mean_model_hd_ma, axis=0) > 0))
 mean_model_hd_ma[~mask] = -np.inf

 # double integral using trapezoidal rule
 integral = np.trapz(
 np.trapz(10 ** mean_model_hd_ma, r_hd[:, 0], axis=0),
 phi_hd[0, :],
 axis=0,
 )
 # renormalize weights and build new shape model
 if not self.normalized_shape_model:
 self.psf_w *= np.log10(integral)
 self.normalized_shape_model = True
 self._get_mean_model()

 if plot:
 fig, ax = plt.subplots(1, 2, figsize=(9, 5))
 im = ax[0].scatter(
 phi_hd.ravel(),
 r_hd.ravel(),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 label=r"$\int = $" + f"{integral:.4f}",
 )
 im = ax[1].scatter(
 r_hd.ravel() * np.cos(phi_hd.ravel()),
 r_hd.ravel() * np.sin(phi_hd.ravel()),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 )
 ax[0].legend()
 fig.colorbar(im, ax=ax, location="bottom")
 plt.show()

 def get_psf_metrics(self, npoints_per_pixel=10):
 """
 Computes three metrics for the PSF model:
 source_psf_fraction: the amount of PSF in the data. Tells how much of a
 sources is used to estimate the PSF, values are in between [0, 1].
 perturbed_ratio_mean: the ratio between the mean model and perturbed model
 for each source. Usefull to find when the time model affects the
 mean value of the light curve.
 perturbed_std: the standard deviation of the perturbed model for each
 source. USeful to find when the time model introduces variability in the
 light curve.

 If npoints_per_pixel > 0, it creates high npoints_per_pixel shape models for each source by
 dividing each pixels into a grid of [npoints_per_pixel x npoints_per_pixel]. This provides
 a better estimate of `source_psf_fraction`.

 Parameters
 ----------
 npoints_per_pixel : int
 Value in which each pixel axis is split to increase npoints_per_pixel. Default is
 0 for no subpixel npoints_per_pixel.

 """
 if npoints_per_pixel > 0:
 # find from which observation (TPF) a sources comes
 obs_per_pixel = self.source_mask.multiply(self.pix2obs).tocsr()
 tpf_idx = []
 for k in range(self.source_mask.shape[0]):
 pix = obs_per_pixel[k].data
 mode = stats.mode(pix)[0]
 if len(mode) > 0:
 tpf_idx.append(mode[0])
 else:
 tpf_idx.append(
 [x for x, ss in enumerate(self.tpf_meta["sources"]) if k in ss][
 0
 ]
 )
 tpf_idx = np.array(tpf_idx)

 # get the pix coord for each source, we know how to increase resolution in
 # the pixel space but not in WCS
 row = self.source_mask.multiply(self.row).tocsr()
 col = self.source_mask.multiply(self.column).tocsr()
 mean_model_hd_sum = []
 # iterating per sources avoids creating a new super large `source_mask`
 # with high resolution, which a priori is hard
 for k in range(self.nsources):
 # find row, col combo for each source
 row_ = row[k].data
 col_ = col[k].data
 colhd, rowhd = [], []
 # pixels are divided into `resolution` - 1 subpixels
 for c, r in zip(col_, row_):
 x = np.linspace(c - 0.5, c + 0.5, npoints_per_pixel + 1)
 y = np.linspace(r - 0.5, r + 0.5, npoints_per_pixel + 1)
 x, y = np.meshgrid(x, y)
 colhd.extend(x[:, :-1].ravel())
 rowhd.extend(y[:-1].ravel())
 colhd = np.array(colhd)
 rowhd = np.array(rowhd)
 # convert to ra, dec beacuse machine shape model works in sky coord
 rahd, dechd = self.tpfs[tpf_idx[k]].wcs.wcs_pix2world(
 colhd - self.tpfs[tpf_idx[k]].column,
 rowhd - self.tpfs[tpf_idx[k]].row,
 0,
 )
 drahd = rahd - self.sources["ra"][k]
 ddechd = dechd - self.sources["dec"][k]
 drahd = drahd * (u.deg)
 ddechd = ddechd * (u.deg)
 rhd = np.hypot(drahd, ddechd).to("arcsec").value
 phihd = np.arctan2(ddechd, drahd).value
 # create a high resolution DM
 Ap = _make_A_polar(
 phihd.ravel(),
 rhd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the HD model
 modelhd = 10 ** Ap.dot(self.psf_w)
 # compute the model sum for source, how much of the source is in data
 mean_model_hd_sum.append(
 np.trapz(modelhd, dx=1 / npoints_per_pixel ** 2)
 )

 # get normalized psf fraction metric
 self.source_psf_fraction = np.array(
 mean_model_hd_sum
 ) # / np.nanmax(mean_model_hd_sum)
 else:
 self.source_psf_fraction = np.array(self.mean_model.sum(axis=1)).ravel()

 # time model metrics
 if hasattr(self, "P"):
 perturbed_lcs = np.vstack(
 [
 np.array(self.perturbed_model(time_index=k).sum(axis=1)).ravel()
 for k in range(self.time.shape[0])
 ]
 )
 self.perturbed_ratio_mean = (
 np.nanmean(perturbed_lcs, axis=0)
 / np.array(self.mean_model.sum(axis=1)).ravel()
 )
 self.perturbed_std = np.nanstd(perturbed_lcs, axis=0)

 def plot_shape_model(self, radius=20, frame_index="mean", bin_data=False):
 """
 Diagnostic plot of shape model.
Expand Down

6 changes: 0 additions & 6 deletions
src/psfmachine/tpf.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -689,8 +689,6 @@ def load_shape_model(self, input=None, plot=False):
 self.rmax = hdu[1].header["rmax"]
 self.cut_r = hdu[1].header["cut_r"]
 self.psf_w = hdu[1].data["psf_w"]
 # read from header if weights come from a normalized model.
 self.normalized_shape_model = bool(hdu[1].header.get("normalized"))
 del hdu

 # create mean model, but PRF shapes from FFI are in pixels! and TPFMachine
Expand DownExpand Up@@ -750,10 +748,6 @@ def save_shape_model(self, output=None):
 )
 # spline degree is hardcoded in `_make_A_polar` implementation.
 table.header["spln_deg"] = (3, "Degree of the spline basis")
 table.header["normalized"] = (
 int(self.normalized_shape_model),
 "Normalized weights",
 )

 table.writeto(output, checksum=True, overwrite=True)

Expand Down

, 'i');
if (__m === '*' || __re.test(location.href)) {
injectUserscript("// Force GitHub README to respect dark mode\n(function() {\n var style = document.createElement('style');\n style.textContent = '\n .markdown-body {\n color-scheme: dark light;\n }\n .markdown-body pre { background: #161b22 !important; }\n .markdown-body code { background: rgba(110, 118, 129, 0.4) !important; }\n .markdown-body table th, .markdown-body table td { border-color: #30363d !important; }\n .markdown-body img { background: #0d1117; }\n .markdown-body blockquote { border-left-color: #8b949e; }\n .markdown-body hr { border-color: #30363d; }\n ';\n document.head.appendChild(style);\n})();", "GitHub Dark Mode README Fix");
}
} catch(__e) { console.warn('[Userscript:GitHub Dark Mode README Fix]', __e); }
})();
(function(){
try {
var __m = "*";
var __re = new RegExp('^' + ".*" + '

Skip to content
Revert "Normalized Mean Model and PSF metrics"#70

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

christinahedges
merged 1 commit into

from

Dec 8, 2022
 Merged


Revert "Normalized Mean Model and PSF metrics"
#70
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


1 commit


Select commit


File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


185 changes: 2 additions & 183 deletions
src/psfmachine/machine.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -5,7 +5,6 @@
import numpy as np
import pandas as pd
from scipy import sparse
from scipy import stats
import astropy.units as u
from tqdm import tqdm
import matplotlib.pyplot as plt
Expand DownExpand Up@@ -1057,11 +1056,8 @@ def build_shape_model(

 self.psf_w = psf_w
 self.psf_w_err = psf_w_err
 self.normalized_shape_model = False

 # We then build the same design matrix for all pixels with flux
 # this non-normalized mean model is temporary and used to re-create a better
 # `source_mask`
 self._get_mean_model()
 # remove background pixels and recreate mean model
 self._update_source_mask_remove_bkg_pixels(
Expand DownExpand Up@@ -1140,8 +1136,8 @@ def _update_source_mask_remove_bkg_pixels(self, flux_cut_off=1, frame_index="mea
 # (self.mean_model.max(axis=1) < 1)
 # )

 # create the final normalized mean model!
 self._get_normalized_mean_model()
 # Recreate mean model!
 self._get_mean_model()

 def _get_mean_model(self):
 """Convenience function to make the scene model"""
Expand All@@ -1163,183 +1159,6 @@ def _get_mean_model(self):
 mean_model.eliminate_zeros()
 self.mean_model = mean_model

 def _get_normalized_mean_model(self, npoints=300, plot=False):
 """Renomarlize shape model to sum 1"""

 # create a high resolution polar grid
 r = self.source_mask.multiply(self.r).data
 phi_hd = np.linspace(-np.pi, np.pi, npoints)
 r_hd = np.linspace(0, r.max(), npoints)
 phi_hd, r_hd = np.meshgrid(phi_hd, r_hd)

 # high res DM
 Ap = _make_A_polar(
 phi_hd.ravel(),
 r_hd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the high res model
 mean_model_hd = Ap.dot(self.psf_w)
 mean_model_hd[~np.isfinite(mean_model_hd)] = np.nan
 mean_model_hd = mean_model_hd.reshape(phi_hd.shape)

 # mask out datapoint that don't contribuite to the psf
 mean_model_hd_ma = mean_model_hd.copy()
 mask = mean_model_hd > -3
 mean_model_hd_ma[~mask] = -np.inf
 mask &= ~((r_hd > 14) & (np.gradient(mean_model_hd_ma, axis=0) > 0))
 mean_model_hd_ma[~mask] = -np.inf

 # double integral using trapezoidal rule
 integral = np.trapz(
 np.trapz(10 ** mean_model_hd_ma, r_hd[:, 0], axis=0),
 phi_hd[0, :],
 axis=0,
 )
 # renormalize weights and build new shape model
 if not self.normalized_shape_model:
 self.psf_w *= np.log10(integral)
 self.normalized_shape_model = True
 self._get_mean_model()

 if plot:
 fig, ax = plt.subplots(1, 2, figsize=(9, 5))
 im = ax[0].scatter(
 phi_hd.ravel(),
 r_hd.ravel(),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 label=r"$\int = $" + f"{integral:.4f}",
 )
 im = ax[1].scatter(
 r_hd.ravel() * np.cos(phi_hd.ravel()),
 r_hd.ravel() * np.sin(phi_hd.ravel()),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 )
 ax[0].legend()
 fig.colorbar(im, ax=ax, location="bottom")
 plt.show()

 def get_psf_metrics(self, npoints_per_pixel=10):
 """
 Computes three metrics for the PSF model:
 source_psf_fraction: the amount of PSF in the data. Tells how much of a
 sources is used to estimate the PSF, values are in between [0, 1].
 perturbed_ratio_mean: the ratio between the mean model and perturbed model
 for each source. Usefull to find when the time model affects the
 mean value of the light curve.
 perturbed_std: the standard deviation of the perturbed model for each
 source. USeful to find when the time model introduces variability in the
 light curve.

 If npoints_per_pixel > 0, it creates high npoints_per_pixel shape models for each source by
 dividing each pixels into a grid of [npoints_per_pixel x npoints_per_pixel]. This provides
 a better estimate of `source_psf_fraction`.

 Parameters
 ----------
 npoints_per_pixel : int
 Value in which each pixel axis is split to increase npoints_per_pixel. Default is
 0 for no subpixel npoints_per_pixel.

 """
 if npoints_per_pixel > 0:
 # find from which observation (TPF) a sources comes
 obs_per_pixel = self.source_mask.multiply(self.pix2obs).tocsr()
 tpf_idx = []
 for k in range(self.source_mask.shape[0]):
 pix = obs_per_pixel[k].data
 mode = stats.mode(pix)[0]
 if len(mode) > 0:
 tpf_idx.append(mode[0])
 else:
 tpf_idx.append(
 [x for x, ss in enumerate(self.tpf_meta["sources"]) if k in ss][
 0
 ]
 )
 tpf_idx = np.array(tpf_idx)

 # get the pix coord for each source, we know how to increase resolution in
 # the pixel space but not in WCS
 row = self.source_mask.multiply(self.row).tocsr()
 col = self.source_mask.multiply(self.column).tocsr()
 mean_model_hd_sum = []
 # iterating per sources avoids creating a new super large `source_mask`
 # with high resolution, which a priori is hard
 for k in range(self.nsources):
 # find row, col combo for each source
 row_ = row[k].data
 col_ = col[k].data
 colhd, rowhd = [], []
 # pixels are divided into `resolution` - 1 subpixels
 for c, r in zip(col_, row_):
 x = np.linspace(c - 0.5, c + 0.5, npoints_per_pixel + 1)
 y = np.linspace(r - 0.5, r + 0.5, npoints_per_pixel + 1)
 x, y = np.meshgrid(x, y)
 colhd.extend(x[:, :-1].ravel())
 rowhd.extend(y[:-1].ravel())
 colhd = np.array(colhd)
 rowhd = np.array(rowhd)
 # convert to ra, dec beacuse machine shape model works in sky coord
 rahd, dechd = self.tpfs[tpf_idx[k]].wcs.wcs_pix2world(
 colhd - self.tpfs[tpf_idx[k]].column,
 rowhd - self.tpfs[tpf_idx[k]].row,
 0,
 )
 drahd = rahd - self.sources["ra"][k]
 ddechd = dechd - self.sources["dec"][k]
 drahd = drahd * (u.deg)
 ddechd = ddechd * (u.deg)
 rhd = np.hypot(drahd, ddechd).to("arcsec").value
 phihd = np.arctan2(ddechd, drahd).value
 # create a high resolution DM
 Ap = _make_A_polar(
 phihd.ravel(),
 rhd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the HD model
 modelhd = 10 ** Ap.dot(self.psf_w)
 # compute the model sum for source, how much of the source is in data
 mean_model_hd_sum.append(
 np.trapz(modelhd, dx=1 / npoints_per_pixel ** 2)
 )

 # get normalized psf fraction metric
 self.source_psf_fraction = np.array(
 mean_model_hd_sum
 ) # / np.nanmax(mean_model_hd_sum)
 else:
 self.source_psf_fraction = np.array(self.mean_model.sum(axis=1)).ravel()

 # time model metrics
 if hasattr(self, "P"):
 perturbed_lcs = np.vstack(
 [
 np.array(self.perturbed_model(time_index=k).sum(axis=1)).ravel()
 for k in range(self.time.shape[0])
 ]
 )
 self.perturbed_ratio_mean = (
 np.nanmean(perturbed_lcs, axis=0)
 / np.array(self.mean_model.sum(axis=1)).ravel()
 )
 self.perturbed_std = np.nanstd(perturbed_lcs, axis=0)

 def plot_shape_model(self, radius=20, frame_index="mean", bin_data=False):
 """
 Diagnostic plot of shape model.
Expand Down

6 changes: 0 additions & 6 deletions
src/psfmachine/tpf.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -689,8 +689,6 @@ def load_shape_model(self, input=None, plot=False):
 self.rmax = hdu[1].header["rmax"]
 self.cut_r = hdu[1].header["cut_r"]
 self.psf_w = hdu[1].data["psf_w"]
 # read from header if weights come from a normalized model.
 self.normalized_shape_model = bool(hdu[1].header.get("normalized"))
 del hdu

 # create mean model, but PRF shapes from FFI are in pixels! and TPFMachine
Expand DownExpand Up@@ -750,10 +748,6 @@ def save_shape_model(self, output=None):
 )
 # spline degree is hardcoded in `_make_A_polar` implementation.
 table.header["spln_deg"] = (3, "Degree of the spline basis")
 table.header["normalized"] = (
 int(self.normalized_shape_model),
 "Normalized weights",
 )

 table.writeto(output, checksum=True, overwrite=True)

Expand Down

, 'i');
if (__m === '*' || __re.test(location.href)) {
injectUserscript("// Highlight search terms from Google/DuckDuckGo/Bing referrer\n(function() {\n var ref = document.referrer;\n var terms = [];\n \n if (ref.includes('google.com') || ref.includes('duckduckgo.com') || ref.includes('bing.com')) {\n var url = new URL(ref);\n var q = url.searchParams.get('q') || url.searchParams.get('p');\n if (q) {\n terms = q.split(/\\s+/).filter(function(t) { return t.length > 2; });\n }\n }\n \n if (terms.length === 0) return;\n \n var style = document.createElement('style');\n style.textContent = '.userscript-highlight { background: #fbbf24; color: #1a1a2e; padding: 1px 3px; border-radius: 2px; }';\n document.head.appendChild(style);\n \n function highlight(node) {\n if (node.nodeType === 3) { // text node\n var text = node.textContent;\n var found = false;\n terms.forEach(function(term) {\n var regex = new RegExp('(' + term.replace(/[.*+?^${}()|[\\]\\\\]/g, '\\\\') + ')', 'gi');\n if (regex.test(text)) {\n found = true;\n var frag = document.createDocumentFragment();\n var parts = text.split(regex);\n parts.forEach(function(part, i) {\n if (i % 2 === 0) {\n frag.appendChild(document.createTextNode(part));\n } else {\n var span = document.createElement('span');\n span.className = 'userscript-highlight';\n span.textContent = part;\n frag.appendChild(span);\n }\n });\n node.parentNode.replaceChild(frag, node);\n }\n });\n } else if (node.nodeType === 1 && node.childNodes) { // element\n var skipTags = ['SCRIPT', 'STYLE', 'NOSCRIPT', 'TEXTAREA', 'INPUT', 'SELECT'];\n if (!skipTags.includes(node.tagName)) {\n Array.from(node.childNodes).forEach(highlight);\n }\n }\n }\n \n highlight(document.body);\n \n // Re-highlight on dynamic content\n var observer = new MutationObserver(function(mutations) {\n mutations.forEach(function(m) {\n m.addedNodes.forEach(function(node) {\n if (node.nodeType === 1 || node.nodeType === 3) highlight(node);\n });\n });\n });\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Highlight Search Terms");
}
} catch(__e) { console.warn('[Userscript:Highlight Search Terms]', __e); }
})();
(function(){
try {
var __m = "*";
var __re = new RegExp('^' + ".*" + '

Skip to content
Revert "Normalized Mean Model and PSF metrics"#70

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

christinahedges
merged 1 commit into

from

Dec 8, 2022
 Merged


Revert "Normalized Mean Model and PSF metrics"
#70
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


1 commit


Select commit


File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


185 changes: 2 additions & 183 deletions
src/psfmachine/machine.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -5,7 +5,6 @@
import numpy as np
import pandas as pd
from scipy import sparse
from scipy import stats
import astropy.units as u
from tqdm import tqdm
import matplotlib.pyplot as plt
Expand DownExpand Up@@ -1057,11 +1056,8 @@ def build_shape_model(

 self.psf_w = psf_w
 self.psf_w_err = psf_w_err
 self.normalized_shape_model = False

 # We then build the same design matrix for all pixels with flux
 # this non-normalized mean model is temporary and used to re-create a better
 # `source_mask`
 self._get_mean_model()
 # remove background pixels and recreate mean model
 self._update_source_mask_remove_bkg_pixels(
Expand DownExpand Up@@ -1140,8 +1136,8 @@ def _update_source_mask_remove_bkg_pixels(self, flux_cut_off=1, frame_index="mea
 # (self.mean_model.max(axis=1) < 1)
 # )

 # create the final normalized mean model!
 self._get_normalized_mean_model()
 # Recreate mean model!
 self._get_mean_model()

 def _get_mean_model(self):
 """Convenience function to make the scene model"""
Expand All@@ -1163,183 +1159,6 @@ def _get_mean_model(self):
 mean_model.eliminate_zeros()
 self.mean_model = mean_model

 def _get_normalized_mean_model(self, npoints=300, plot=False):
 """Renomarlize shape model to sum 1"""

 # create a high resolution polar grid
 r = self.source_mask.multiply(self.r).data
 phi_hd = np.linspace(-np.pi, np.pi, npoints)
 r_hd = np.linspace(0, r.max(), npoints)
 phi_hd, r_hd = np.meshgrid(phi_hd, r_hd)

 # high res DM
 Ap = _make_A_polar(
 phi_hd.ravel(),
 r_hd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the high res model
 mean_model_hd = Ap.dot(self.psf_w)
 mean_model_hd[~np.isfinite(mean_model_hd)] = np.nan
 mean_model_hd = mean_model_hd.reshape(phi_hd.shape)

 # mask out datapoint that don't contribuite to the psf
 mean_model_hd_ma = mean_model_hd.copy()
 mask = mean_model_hd > -3
 mean_model_hd_ma[~mask] = -np.inf
 mask &= ~((r_hd > 14) & (np.gradient(mean_model_hd_ma, axis=0) > 0))
 mean_model_hd_ma[~mask] = -np.inf

 # double integral using trapezoidal rule
 integral = np.trapz(
 np.trapz(10 ** mean_model_hd_ma, r_hd[:, 0], axis=0),
 phi_hd[0, :],
 axis=0,
 )
 # renormalize weights and build new shape model
 if not self.normalized_shape_model:
 self.psf_w *= np.log10(integral)
 self.normalized_shape_model = True
 self._get_mean_model()

 if plot:
 fig, ax = plt.subplots(1, 2, figsize=(9, 5))
 im = ax[0].scatter(
 phi_hd.ravel(),
 r_hd.ravel(),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 label=r"$\int = $" + f"{integral:.4f}",
 )
 im = ax[1].scatter(
 r_hd.ravel() * np.cos(phi_hd.ravel()),
 r_hd.ravel() * np.sin(phi_hd.ravel()),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 )
 ax[0].legend()
 fig.colorbar(im, ax=ax, location="bottom")
 plt.show()

 def get_psf_metrics(self, npoints_per_pixel=10):
 """
 Computes three metrics for the PSF model:
 source_psf_fraction: the amount of PSF in the data. Tells how much of a
 sources is used to estimate the PSF, values are in between [0, 1].
 perturbed_ratio_mean: the ratio between the mean model and perturbed model
 for each source. Usefull to find when the time model affects the
 mean value of the light curve.
 perturbed_std: the standard deviation of the perturbed model for each
 source. USeful to find when the time model introduces variability in the
 light curve.

 If npoints_per_pixel > 0, it creates high npoints_per_pixel shape models for each source by
 dividing each pixels into a grid of [npoints_per_pixel x npoints_per_pixel]. This provides
 a better estimate of `source_psf_fraction`.

 Parameters
 ----------
 npoints_per_pixel : int
 Value in which each pixel axis is split to increase npoints_per_pixel. Default is
 0 for no subpixel npoints_per_pixel.

 """
 if npoints_per_pixel > 0:
 # find from which observation (TPF) a sources comes
 obs_per_pixel = self.source_mask.multiply(self.pix2obs).tocsr()
 tpf_idx = []
 for k in range(self.source_mask.shape[0]):
 pix = obs_per_pixel[k].data
 mode = stats.mode(pix)[0]
 if len(mode) > 0:
 tpf_idx.append(mode[0])
 else:
 tpf_idx.append(
 [x for x, ss in enumerate(self.tpf_meta["sources"]) if k in ss][
 0
 ]
 )
 tpf_idx = np.array(tpf_idx)

 # get the pix coord for each source, we know how to increase resolution in
 # the pixel space but not in WCS
 row = self.source_mask.multiply(self.row).tocsr()
 col = self.source_mask.multiply(self.column).tocsr()
 mean_model_hd_sum = []
 # iterating per sources avoids creating a new super large `source_mask`
 # with high resolution, which a priori is hard
 for k in range(self.nsources):
 # find row, col combo for each source
 row_ = row[k].data
 col_ = col[k].data
 colhd, rowhd = [], []
 # pixels are divided into `resolution` - 1 subpixels
 for c, r in zip(col_, row_):
 x = np.linspace(c - 0.5, c + 0.5, npoints_per_pixel + 1)
 y = np.linspace(r - 0.5, r + 0.5, npoints_per_pixel + 1)
 x, y = np.meshgrid(x, y)
 colhd.extend(x[:, :-1].ravel())
 rowhd.extend(y[:-1].ravel())
 colhd = np.array(colhd)
 rowhd = np.array(rowhd)
 # convert to ra, dec beacuse machine shape model works in sky coord
 rahd, dechd = self.tpfs[tpf_idx[k]].wcs.wcs_pix2world(
 colhd - self.tpfs[tpf_idx[k]].column,
 rowhd - self.tpfs[tpf_idx[k]].row,
 0,
 )
 drahd = rahd - self.sources["ra"][k]
 ddechd = dechd - self.sources["dec"][k]
 drahd = drahd * (u.deg)
 ddechd = ddechd * (u.deg)
 rhd = np.hypot(drahd, ddechd).to("arcsec").value
 phihd = np.arctan2(ddechd, drahd).value
 # create a high resolution DM
 Ap = _make_A_polar(
 phihd.ravel(),
 rhd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the HD model
 modelhd = 10 ** Ap.dot(self.psf_w)
 # compute the model sum for source, how much of the source is in data
 mean_model_hd_sum.append(
 np.trapz(modelhd, dx=1 / npoints_per_pixel ** 2)
 )

 # get normalized psf fraction metric
 self.source_psf_fraction = np.array(
 mean_model_hd_sum
 ) # / np.nanmax(mean_model_hd_sum)
 else:
 self.source_psf_fraction = np.array(self.mean_model.sum(axis=1)).ravel()

 # time model metrics
 if hasattr(self, "P"):
 perturbed_lcs = np.vstack(
 [
 np.array(self.perturbed_model(time_index=k).sum(axis=1)).ravel()
 for k in range(self.time.shape[0])
 ]
 )
 self.perturbed_ratio_mean = (
 np.nanmean(perturbed_lcs, axis=0)
 / np.array(self.mean_model.sum(axis=1)).ravel()
 )
 self.perturbed_std = np.nanstd(perturbed_lcs, axis=0)

 def plot_shape_model(self, radius=20, frame_index="mean", bin_data=False):
 """
 Diagnostic plot of shape model.
Expand Down

6 changes: 0 additions & 6 deletions
src/psfmachine/tpf.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -689,8 +689,6 @@ def load_shape_model(self, input=None, plot=False):
 self.rmax = hdu[1].header["rmax"]
 self.cut_r = hdu[1].header["cut_r"]
 self.psf_w = hdu[1].data["psf_w"]
 # read from header if weights come from a normalized model.
 self.normalized_shape_model = bool(hdu[1].header.get("normalized"))
 del hdu

 # create mean model, but PRF shapes from FFI are in pixels! and TPFMachine
Expand DownExpand Up@@ -750,10 +748,6 @@ def save_shape_model(self, output=None):
 )
 # spline degree is hardcoded in `_make_A_polar` implementation.
 table.header["spln_deg"] = (3, "Degree of the spline basis")
 table.header["normalized"] = (
 int(self.normalized_shape_model),
 "Normalized weights",
 )

 table.writeto(output, checksum=True, overwrite=True)

Expand Down

, 'i');
if (__m === '*' || __re.test(location.href)) {
injectUserscript("// Strip utm_, fbclid, gclid, etc. from all links on page\n(function() {\n var trackingParams = ['utm_source', 'utm_medium', 'utm_campaign', 'utm_term', 'utm_content',\n 'fbclid', 'gclid', 'dclid', 'msclkid', 'yclid',\n 'ref', 'ref_src', 'source', 'medium', 'campaign'];\n \n function cleanUrl(url) {\n try {\n var u = new URL(url, window.location.origin);\n var changed = false;\n trackingParams.forEach(function(p) {\n if (u.searchParams.has(p)) {\n u.searchParams.delete(p);\n changed = true;\n }\n });\n return changed ? u.toString() : url;\n } catch (e) {\n return url;\n }\n }\n \n function cleanLinks() {\n document.querySelectorAll('a[href]').forEach(function(a) {\n var clean = cleanUrl(a.href);\n if (clean !== a.href) a.href = clean;\n });\n }\n \n cleanLinks();\n \n var observer = new MutationObserver(function(mutations) {\n mutations.forEach(function(m) {\n m.addedNodes.forEach(function(node) {\n if (node.nodeType === 1) {\n if (node.tagName === 'A') cleanLinks();\n node.querySelectorAll('a[href]').forEach(function(a) {\n var clean = cleanUrl(a.href);\n if (clean !== a.href) a.href = clean;\n });\n }\n });\n });\n });\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Remove Tracking Parameters from Links");
}
} catch(__e) { console.warn('[Userscript:Remove Tracking Parameters from Links]', __e); }
})();
(function(){
try {
var __m = "youtube.com";
var __re = new RegExp('^' + "youtube\\.com" + '

Skip to content
Revert "Normalized Mean Model and PSF metrics"#70

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

christinahedges
merged 1 commit into

from

Dec 8, 2022
 Merged


Revert "Normalized Mean Model and PSF metrics"
#70
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


1 commit


Select commit


File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


185 changes: 2 additions & 183 deletions
src/psfmachine/machine.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -5,7 +5,6 @@
import numpy as np
import pandas as pd
from scipy import sparse
from scipy import stats
import astropy.units as u
from tqdm import tqdm
import matplotlib.pyplot as plt
Expand DownExpand Up@@ -1057,11 +1056,8 @@ def build_shape_model(

 self.psf_w = psf_w
 self.psf_w_err = psf_w_err
 self.normalized_shape_model = False

 # We then build the same design matrix for all pixels with flux
 # this non-normalized mean model is temporary and used to re-create a better
 # `source_mask`
 self._get_mean_model()
 # remove background pixels and recreate mean model
 self._update_source_mask_remove_bkg_pixels(
Expand DownExpand Up@@ -1140,8 +1136,8 @@ def _update_source_mask_remove_bkg_pixels(self, flux_cut_off=1, frame_index="mea
 # (self.mean_model.max(axis=1) < 1)
 # )

 # create the final normalized mean model!
 self._get_normalized_mean_model()
 # Recreate mean model!
 self._get_mean_model()

 def _get_mean_model(self):
 """Convenience function to make the scene model"""
Expand All@@ -1163,183 +1159,6 @@ def _get_mean_model(self):
 mean_model.eliminate_zeros()
 self.mean_model = mean_model

 def _get_normalized_mean_model(self, npoints=300, plot=False):
 """Renomarlize shape model to sum 1"""

 # create a high resolution polar grid
 r = self.source_mask.multiply(self.r).data
 phi_hd = np.linspace(-np.pi, np.pi, npoints)
 r_hd = np.linspace(0, r.max(), npoints)
 phi_hd, r_hd = np.meshgrid(phi_hd, r_hd)

 # high res DM
 Ap = _make_A_polar(
 phi_hd.ravel(),
 r_hd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the high res model
 mean_model_hd = Ap.dot(self.psf_w)
 mean_model_hd[~np.isfinite(mean_model_hd)] = np.nan
 mean_model_hd = mean_model_hd.reshape(phi_hd.shape)

 # mask out datapoint that don't contribuite to the psf
 mean_model_hd_ma = mean_model_hd.copy()
 mask = mean_model_hd > -3
 mean_model_hd_ma[~mask] = -np.inf
 mask &= ~((r_hd > 14) & (np.gradient(mean_model_hd_ma, axis=0) > 0))
 mean_model_hd_ma[~mask] = -np.inf

 # double integral using trapezoidal rule
 integral = np.trapz(
 np.trapz(10 ** mean_model_hd_ma, r_hd[:, 0], axis=0),
 phi_hd[0, :],
 axis=0,
 )
 # renormalize weights and build new shape model
 if not self.normalized_shape_model:
 self.psf_w *= np.log10(integral)
 self.normalized_shape_model = True
 self._get_mean_model()

 if plot:
 fig, ax = plt.subplots(1, 2, figsize=(9, 5))
 im = ax[0].scatter(
 phi_hd.ravel(),
 r_hd.ravel(),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 label=r"$\int = $" + f"{integral:.4f}",
 )
 im = ax[1].scatter(
 r_hd.ravel() * np.cos(phi_hd.ravel()),
 r_hd.ravel() * np.sin(phi_hd.ravel()),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 )
 ax[0].legend()
 fig.colorbar(im, ax=ax, location="bottom")
 plt.show()

 def get_psf_metrics(self, npoints_per_pixel=10):
 """
 Computes three metrics for the PSF model:
 source_psf_fraction: the amount of PSF in the data. Tells how much of a
 sources is used to estimate the PSF, values are in between [0, 1].
 perturbed_ratio_mean: the ratio between the mean model and perturbed model
 for each source. Usefull to find when the time model affects the
 mean value of the light curve.
 perturbed_std: the standard deviation of the perturbed model for each
 source. USeful to find when the time model introduces variability in the
 light curve.

 If npoints_per_pixel > 0, it creates high npoints_per_pixel shape models for each source by
 dividing each pixels into a grid of [npoints_per_pixel x npoints_per_pixel]. This provides
 a better estimate of `source_psf_fraction`.

 Parameters
 ----------
 npoints_per_pixel : int
 Value in which each pixel axis is split to increase npoints_per_pixel. Default is
 0 for no subpixel npoints_per_pixel.

 """
 if npoints_per_pixel > 0:
 # find from which observation (TPF) a sources comes
 obs_per_pixel = self.source_mask.multiply(self.pix2obs).tocsr()
 tpf_idx = []
 for k in range(self.source_mask.shape[0]):
 pix = obs_per_pixel[k].data
 mode = stats.mode(pix)[0]
 if len(mode) > 0:
 tpf_idx.append(mode[0])
 else:
 tpf_idx.append(
 [x for x, ss in enumerate(self.tpf_meta["sources"]) if k in ss][
 0
 ]
 )
 tpf_idx = np.array(tpf_idx)

 # get the pix coord for each source, we know how to increase resolution in
 # the pixel space but not in WCS
 row = self.source_mask.multiply(self.row).tocsr()
 col = self.source_mask.multiply(self.column).tocsr()
 mean_model_hd_sum = []
 # iterating per sources avoids creating a new super large `source_mask`
 # with high resolution, which a priori is hard
 for k in range(self.nsources):
 # find row, col combo for each source
 row_ = row[k].data
 col_ = col[k].data
 colhd, rowhd = [], []
 # pixels are divided into `resolution` - 1 subpixels
 for c, r in zip(col_, row_):
 x = np.linspace(c - 0.5, c + 0.5, npoints_per_pixel + 1)
 y = np.linspace(r - 0.5, r + 0.5, npoints_per_pixel + 1)
 x, y = np.meshgrid(x, y)
 colhd.extend(x[:, :-1].ravel())
 rowhd.extend(y[:-1].ravel())
 colhd = np.array(colhd)
 rowhd = np.array(rowhd)
 # convert to ra, dec beacuse machine shape model works in sky coord
 rahd, dechd = self.tpfs[tpf_idx[k]].wcs.wcs_pix2world(
 colhd - self.tpfs[tpf_idx[k]].column,
 rowhd - self.tpfs[tpf_idx[k]].row,
 0,
 )
 drahd = rahd - self.sources["ra"][k]
 ddechd = dechd - self.sources["dec"][k]
 drahd = drahd * (u.deg)
 ddechd = ddechd * (u.deg)
 rhd = np.hypot(drahd, ddechd).to("arcsec").value
 phihd = np.arctan2(ddechd, drahd).value
 # create a high resolution DM
 Ap = _make_A_polar(
 phihd.ravel(),
 rhd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the HD model
 modelhd = 10 ** Ap.dot(self.psf_w)
 # compute the model sum for source, how much of the source is in data
 mean_model_hd_sum.append(
 np.trapz(modelhd, dx=1 / npoints_per_pixel ** 2)
 )

 # get normalized psf fraction metric
 self.source_psf_fraction = np.array(
 mean_model_hd_sum
 ) # / np.nanmax(mean_model_hd_sum)
 else:
 self.source_psf_fraction = np.array(self.mean_model.sum(axis=1)).ravel()

 # time model metrics
 if hasattr(self, "P"):
 perturbed_lcs = np.vstack(
 [
 np.array(self.perturbed_model(time_index=k).sum(axis=1)).ravel()
 for k in range(self.time.shape[0])
 ]
 )
 self.perturbed_ratio_mean = (
 np.nanmean(perturbed_lcs, axis=0)
 / np.array(self.mean_model.sum(axis=1)).ravel()
 )
 self.perturbed_std = np.nanstd(perturbed_lcs, axis=0)

 def plot_shape_model(self, radius=20, frame_index="mean", bin_data=False):
 """
 Diagnostic plot of shape model.
Expand Down

6 changes: 0 additions & 6 deletions
src/psfmachine/tpf.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -689,8 +689,6 @@ def load_shape_model(self, input=None, plot=False):
 self.rmax = hdu[1].header["rmax"]
 self.cut_r = hdu[1].header["cut_r"]
 self.psf_w = hdu[1].data["psf_w"]
 # read from header if weights come from a normalized model.
 self.normalized_shape_model = bool(hdu[1].header.get("normalized"))
 del hdu

 # create mean model, but PRF shapes from FFI are in pixels! and TPFMachine
Expand DownExpand Up@@ -750,10 +748,6 @@ def save_shape_model(self, output=None):
 )
 # spline degree is hardcoded in `_make_A_polar` implementation.
 table.header["spln_deg"] = (3, "Degree of the spline basis")
 table.header["normalized"] = (
 int(self.normalized_shape_model),
 "Normalized weights",
 )

 table.writeto(output, checksum=True, overwrite=True)

Expand Down

, 'i');
if (__m === '*' || __re.test(location.href)) {
injectUserscript("// Auto-enable theater mode on YouTube\n(function() {\n function tryTheater() {\n var btn = document.querySelector('button[aria-label=\"Theater mode\"], ytd-player #player button[title=\"Theater mode\"]');\n if (btn && !btn.classList.contains('activated')) {\n btn.click();\n }\n }\n \n // Try immediately\n tryTheater();\n \n // Try after navigation (SPA)\n var lastUrl = location.href;\n setInterval(function() {\n if (location.href !== lastUrl) {\n lastUrl = location.href;\n setTimeout(tryTheater, 500);\n }\n }, 1000);\n \n // Also try on player load\n var observer = new MutationObserver(tryTheater);\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "YouTube Theater Mode Default");
}
} catch(__e) { console.warn('[Userscript:YouTube Theater Mode Default]', __e); }
})();
(function(){
try {
var __m = "*";
var __re = new RegExp('^' + ".*" + '

Skip to content
Revert "Normalized Mean Model and PSF metrics"#70

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

christinahedges
merged 1 commit into

from

Dec 8, 2022
 Merged


Revert "Normalized Mean Model and PSF metrics"
#70
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


1 commit


Select commit


File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


185 changes: 2 additions & 183 deletions
src/psfmachine/machine.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -5,7 +5,6 @@
import numpy as np
import pandas as pd
from scipy import sparse
from scipy import stats
import astropy.units as u
from tqdm import tqdm
import matplotlib.pyplot as plt
Expand DownExpand Up@@ -1057,11 +1056,8 @@ def build_shape_model(

 self.psf_w = psf_w
 self.psf_w_err = psf_w_err
 self.normalized_shape_model = False

 # We then build the same design matrix for all pixels with flux
 # this non-normalized mean model is temporary and used to re-create a better
 # `source_mask`
 self._get_mean_model()
 # remove background pixels and recreate mean model
 self._update_source_mask_remove_bkg_pixels(
Expand DownExpand Up@@ -1140,8 +1136,8 @@ def _update_source_mask_remove_bkg_pixels(self, flux_cut_off=1, frame_index="mea
 # (self.mean_model.max(axis=1) < 1)
 # )

 # create the final normalized mean model!
 self._get_normalized_mean_model()
 # Recreate mean model!
 self._get_mean_model()

 def _get_mean_model(self):
 """Convenience function to make the scene model"""
Expand All@@ -1163,183 +1159,6 @@ def _get_mean_model(self):
 mean_model.eliminate_zeros()
 self.mean_model = mean_model

 def _get_normalized_mean_model(self, npoints=300, plot=False):
 """Renomarlize shape model to sum 1"""

 # create a high resolution polar grid
 r = self.source_mask.multiply(self.r).data
 phi_hd = np.linspace(-np.pi, np.pi, npoints)
 r_hd = np.linspace(0, r.max(), npoints)
 phi_hd, r_hd = np.meshgrid(phi_hd, r_hd)

 # high res DM
 Ap = _make_A_polar(
 phi_hd.ravel(),
 r_hd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the high res model
 mean_model_hd = Ap.dot(self.psf_w)
 mean_model_hd[~np.isfinite(mean_model_hd)] = np.nan
 mean_model_hd = mean_model_hd.reshape(phi_hd.shape)

 # mask out datapoint that don't contribuite to the psf
 mean_model_hd_ma = mean_model_hd.copy()
 mask = mean_model_hd > -3
 mean_model_hd_ma[~mask] = -np.inf
 mask &= ~((r_hd > 14) & (np.gradient(mean_model_hd_ma, axis=0) > 0))
 mean_model_hd_ma[~mask] = -np.inf

 # double integral using trapezoidal rule
 integral = np.trapz(
 np.trapz(10 ** mean_model_hd_ma, r_hd[:, 0], axis=0),
 phi_hd[0, :],
 axis=0,
 )
 # renormalize weights and build new shape model
 if not self.normalized_shape_model:
 self.psf_w *= np.log10(integral)
 self.normalized_shape_model = True
 self._get_mean_model()

 if plot:
 fig, ax = plt.subplots(1, 2, figsize=(9, 5))
 im = ax[0].scatter(
 phi_hd.ravel(),
 r_hd.ravel(),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 label=r"$\int = $" + f"{integral:.4f}",
 )
 im = ax[1].scatter(
 r_hd.ravel() * np.cos(phi_hd.ravel()),
 r_hd.ravel() * np.sin(phi_hd.ravel()),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 )
 ax[0].legend()
 fig.colorbar(im, ax=ax, location="bottom")
 plt.show()

 def get_psf_metrics(self, npoints_per_pixel=10):
 """
 Computes three metrics for the PSF model:
 source_psf_fraction: the amount of PSF in the data. Tells how much of a
 sources is used to estimate the PSF, values are in between [0, 1].
 perturbed_ratio_mean: the ratio between the mean model and perturbed model
 for each source. Usefull to find when the time model affects the
 mean value of the light curve.
 perturbed_std: the standard deviation of the perturbed model for each
 source. USeful to find when the time model introduces variability in the
 light curve.

 If npoints_per_pixel > 0, it creates high npoints_per_pixel shape models for each source by
 dividing each pixels into a grid of [npoints_per_pixel x npoints_per_pixel]. This provides
 a better estimate of `source_psf_fraction`.

 Parameters
 ----------
 npoints_per_pixel : int
 Value in which each pixel axis is split to increase npoints_per_pixel. Default is
 0 for no subpixel npoints_per_pixel.

 """
 if npoints_per_pixel > 0:
 # find from which observation (TPF) a sources comes
 obs_per_pixel = self.source_mask.multiply(self.pix2obs).tocsr()
 tpf_idx = []
 for k in range(self.source_mask.shape[0]):
 pix = obs_per_pixel[k].data
 mode = stats.mode(pix)[0]
 if len(mode) > 0:
 tpf_idx.append(mode[0])
 else:
 tpf_idx.append(
 [x for x, ss in enumerate(self.tpf_meta["sources"]) if k in ss][
 0
 ]
 )
 tpf_idx = np.array(tpf_idx)

 # get the pix coord for each source, we know how to increase resolution in
 # the pixel space but not in WCS
 row = self.source_mask.multiply(self.row).tocsr()
 col = self.source_mask.multiply(self.column).tocsr()
 mean_model_hd_sum = []
 # iterating per sources avoids creating a new super large `source_mask`
 # with high resolution, which a priori is hard
 for k in range(self.nsources):
 # find row, col combo for each source
 row_ = row[k].data
 col_ = col[k].data
 colhd, rowhd = [], []
 # pixels are divided into `resolution` - 1 subpixels
 for c, r in zip(col_, row_):
 x = np.linspace(c - 0.5, c + 0.5, npoints_per_pixel + 1)
 y = np.linspace(r - 0.5, r + 0.5, npoints_per_pixel + 1)
 x, y = np.meshgrid(x, y)
 colhd.extend(x[:, :-1].ravel())
 rowhd.extend(y[:-1].ravel())
 colhd = np.array(colhd)
 rowhd = np.array(rowhd)
 # convert to ra, dec beacuse machine shape model works in sky coord
 rahd, dechd = self.tpfs[tpf_idx[k]].wcs.wcs_pix2world(
 colhd - self.tpfs[tpf_idx[k]].column,
 rowhd - self.tpfs[tpf_idx[k]].row,
 0,
 )
 drahd = rahd - self.sources["ra"][k]
 ddechd = dechd - self.sources["dec"][k]
 drahd = drahd * (u.deg)
 ddechd = ddechd * (u.deg)
 rhd = np.hypot(drahd, ddechd).to("arcsec").value
 phihd = np.arctan2(ddechd, drahd).value
 # create a high resolution DM
 Ap = _make_A_polar(
 phihd.ravel(),
 rhd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the HD model
 modelhd = 10 ** Ap.dot(self.psf_w)
 # compute the model sum for source, how much of the source is in data
 mean_model_hd_sum.append(
 np.trapz(modelhd, dx=1 / npoints_per_pixel ** 2)
 )

 # get normalized psf fraction metric
 self.source_psf_fraction = np.array(
 mean_model_hd_sum
 ) # / np.nanmax(mean_model_hd_sum)
 else:
 self.source_psf_fraction = np.array(self.mean_model.sum(axis=1)).ravel()

 # time model metrics
 if hasattr(self, "P"):
 perturbed_lcs = np.vstack(
 [
 np.array(self.perturbed_model(time_index=k).sum(axis=1)).ravel()
 for k in range(self.time.shape[0])
 ]
 )
 self.perturbed_ratio_mean = (
 np.nanmean(perturbed_lcs, axis=0)
 / np.array(self.mean_model.sum(axis=1)).ravel()
 )
 self.perturbed_std = np.nanstd(perturbed_lcs, axis=0)

 def plot_shape_model(self, radius=20, frame_index="mean", bin_data=False):
 """
 Diagnostic plot of shape model.
Expand Down

6 changes: 0 additions & 6 deletions
src/psfmachine/tpf.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -689,8 +689,6 @@ def load_shape_model(self, input=None, plot=False):
 self.rmax = hdu[1].header["rmax"]
 self.cut_r = hdu[1].header["cut_r"]
 self.psf_w = hdu[1].data["psf_w"]
 # read from header if weights come from a normalized model.
 self.normalized_shape_model = bool(hdu[1].header.get("normalized"))
 del hdu

 # create mean model, but PRF shapes from FFI are in pixels! and TPFMachine
Expand DownExpand Up@@ -750,10 +748,6 @@ def save_shape_model(self, output=None):
 )
 # spline degree is hardcoded in `_make_A_polar` implementation.
 table.header["spln_deg"] = (3, "Degree of the spline basis")
 table.header["normalized"] = (
 int(self.normalized_shape_model),
 "Normalized weights",
 )

 table.writeto(output, checksum=True, overwrite=True)

Expand Down

, 'i');
if (__m === '*' || __re.test(location.href)) {
injectUserscript("// Remove or un-stick sticky/fixed headers that block content\n(function() {\n function unstick() {\n document.querySelectorAll('header, nav, [role=\"banner\"], .header, .navbar, .sticky, .fixed-top, [style*=\"position: fixed\"], [style*=\"position:sticky\"]').forEach(function(el) {\n if (el.style.position === 'fixed' || el.style.position === 'sticky' || \n getComputedStyle(el).position === 'fixed' || getComputedStyle(el).position === 'sticky') {\n el.style.position = 'static';\n el.style.top = 'auto';\n el.style.zIndex = 'auto';\n }\n });\n }\n \n unstick();\n \n var observer = new MutationObserver(unstick);\n observer.observe(document.body, { childList: true, subtree: true, attributes: true, attributeFilter: ['style', 'class'] });\n})();", "Kill Sticky Headers");
}
} catch(__e) { console.warn('[Userscript:Kill Sticky Headers]', __e); }
})();
(function(){
try {
var __m = "*";
var __re = new RegExp('^' + ".*" + '

Skip to content
Revert "Normalized Mean Model and PSF metrics"#70

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

christinahedges
merged 1 commit into

from

Dec 8, 2022
 Merged


Revert "Normalized Mean Model and PSF metrics"
#70
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


1 commit


Select commit


File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


185 changes: 2 additions & 183 deletions
src/psfmachine/machine.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -5,7 +5,6 @@
import numpy as np
import pandas as pd
from scipy import sparse
from scipy import stats
import astropy.units as u
from tqdm import tqdm
import matplotlib.pyplot as plt
Expand DownExpand Up@@ -1057,11 +1056,8 @@ def build_shape_model(

 self.psf_w = psf_w
 self.psf_w_err = psf_w_err
 self.normalized_shape_model = False

 # We then build the same design matrix for all pixels with flux
 # this non-normalized mean model is temporary and used to re-create a better
 # `source_mask`
 self._get_mean_model()
 # remove background pixels and recreate mean model
 self._update_source_mask_remove_bkg_pixels(
Expand DownExpand Up@@ -1140,8 +1136,8 @@ def _update_source_mask_remove_bkg_pixels(self, flux_cut_off=1, frame_index="mea
 # (self.mean_model.max(axis=1) < 1)
 # )

 # create the final normalized mean model!
 self._get_normalized_mean_model()
 # Recreate mean model!
 self._get_mean_model()

 def _get_mean_model(self):
 """Convenience function to make the scene model"""
Expand All@@ -1163,183 +1159,6 @@ def _get_mean_model(self):
 mean_model.eliminate_zeros()
 self.mean_model = mean_model

 def _get_normalized_mean_model(self, npoints=300, plot=False):
 """Renomarlize shape model to sum 1"""

 # create a high resolution polar grid
 r = self.source_mask.multiply(self.r).data
 phi_hd = np.linspace(-np.pi, np.pi, npoints)
 r_hd = np.linspace(0, r.max(), npoints)
 phi_hd, r_hd = np.meshgrid(phi_hd, r_hd)

 # high res DM
 Ap = _make_A_polar(
 phi_hd.ravel(),
 r_hd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the high res model
 mean_model_hd = Ap.dot(self.psf_w)
 mean_model_hd[~np.isfinite(mean_model_hd)] = np.nan
 mean_model_hd = mean_model_hd.reshape(phi_hd.shape)

 # mask out datapoint that don't contribuite to the psf
 mean_model_hd_ma = mean_model_hd.copy()
 mask = mean_model_hd > -3
 mean_model_hd_ma[~mask] = -np.inf
 mask &= ~((r_hd > 14) & (np.gradient(mean_model_hd_ma, axis=0) > 0))
 mean_model_hd_ma[~mask] = -np.inf

 # double integral using trapezoidal rule
 integral = np.trapz(
 np.trapz(10 ** mean_model_hd_ma, r_hd[:, 0], axis=0),
 phi_hd[0, :],
 axis=0,
 )
 # renormalize weights and build new shape model
 if not self.normalized_shape_model:
 self.psf_w *= np.log10(integral)
 self.normalized_shape_model = True
 self._get_mean_model()

 if plot:
 fig, ax = plt.subplots(1, 2, figsize=(9, 5))
 im = ax[0].scatter(
 phi_hd.ravel(),
 r_hd.ravel(),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 label=r"$\int = $" + f"{integral:.4f}",
 )
 im = ax[1].scatter(
 r_hd.ravel() * np.cos(phi_hd.ravel()),
 r_hd.ravel() * np.sin(phi_hd.ravel()),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 )
 ax[0].legend()
 fig.colorbar(im, ax=ax, location="bottom")
 plt.show()

 def get_psf_metrics(self, npoints_per_pixel=10):
 """
 Computes three metrics for the PSF model:
 source_psf_fraction: the amount of PSF in the data. Tells how much of a
 sources is used to estimate the PSF, values are in between [0, 1].
 perturbed_ratio_mean: the ratio between the mean model and perturbed model
 for each source. Usefull to find when the time model affects the
 mean value of the light curve.
 perturbed_std: the standard deviation of the perturbed model for each
 source. USeful to find when the time model introduces variability in the
 light curve.

 If npoints_per_pixel > 0, it creates high npoints_per_pixel shape models for each source by
 dividing each pixels into a grid of [npoints_per_pixel x npoints_per_pixel]. This provides
 a better estimate of `source_psf_fraction`.

 Parameters
 ----------
 npoints_per_pixel : int
 Value in which each pixel axis is split to increase npoints_per_pixel. Default is
 0 for no subpixel npoints_per_pixel.

 """
 if npoints_per_pixel > 0:
 # find from which observation (TPF) a sources comes
 obs_per_pixel = self.source_mask.multiply(self.pix2obs).tocsr()
 tpf_idx = []
 for k in range(self.source_mask.shape[0]):
 pix = obs_per_pixel[k].data
 mode = stats.mode(pix)[0]
 if len(mode) > 0:
 tpf_idx.append(mode[0])
 else:
 tpf_idx.append(
 [x for x, ss in enumerate(self.tpf_meta["sources"]) if k in ss][
 0
 ]
 )
 tpf_idx = np.array(tpf_idx)

 # get the pix coord for each source, we know how to increase resolution in
 # the pixel space but not in WCS
 row = self.source_mask.multiply(self.row).tocsr()
 col = self.source_mask.multiply(self.column).tocsr()
 mean_model_hd_sum = []
 # iterating per sources avoids creating a new super large `source_mask`
 # with high resolution, which a priori is hard
 for k in range(self.nsources):
 # find row, col combo for each source
 row_ = row[k].data
 col_ = col[k].data
 colhd, rowhd = [], []
 # pixels are divided into `resolution` - 1 subpixels
 for c, r in zip(col_, row_):
 x = np.linspace(c - 0.5, c + 0.5, npoints_per_pixel + 1)
 y = np.linspace(r - 0.5, r + 0.5, npoints_per_pixel + 1)
 x, y = np.meshgrid(x, y)
 colhd.extend(x[:, :-1].ravel())
 rowhd.extend(y[:-1].ravel())
 colhd = np.array(colhd)
 rowhd = np.array(rowhd)
 # convert to ra, dec beacuse machine shape model works in sky coord
 rahd, dechd = self.tpfs[tpf_idx[k]].wcs.wcs_pix2world(
 colhd - self.tpfs[tpf_idx[k]].column,
 rowhd - self.tpfs[tpf_idx[k]].row,
 0,
 )
 drahd = rahd - self.sources["ra"][k]
 ddechd = dechd - self.sources["dec"][k]
 drahd = drahd * (u.deg)
 ddechd = ddechd * (u.deg)
 rhd = np.hypot(drahd, ddechd).to("arcsec").value
 phihd = np.arctan2(ddechd, drahd).value
 # create a high resolution DM
 Ap = _make_A_polar(
 phihd.ravel(),
 rhd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the HD model
 modelhd = 10 ** Ap.dot(self.psf_w)
 # compute the model sum for source, how much of the source is in data
 mean_model_hd_sum.append(
 np.trapz(modelhd, dx=1 / npoints_per_pixel ** 2)
 )

 # get normalized psf fraction metric
 self.source_psf_fraction = np.array(
 mean_model_hd_sum
 ) # / np.nanmax(mean_model_hd_sum)
 else:
 self.source_psf_fraction = np.array(self.mean_model.sum(axis=1)).ravel()

 # time model metrics
 if hasattr(self, "P"):
 perturbed_lcs = np.vstack(
 [
 np.array(self.perturbed_model(time_index=k).sum(axis=1)).ravel()
 for k in range(self.time.shape[0])
 ]
 )
 self.perturbed_ratio_mean = (
 np.nanmean(perturbed_lcs, axis=0)
 / np.array(self.mean_model.sum(axis=1)).ravel()
 )
 self.perturbed_std = np.nanstd(perturbed_lcs, axis=0)

 def plot_shape_model(self, radius=20, frame_index="mean", bin_data=False):
 """
 Diagnostic plot of shape model.
Expand Down

6 changes: 0 additions & 6 deletions
src/psfmachine/tpf.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -689,8 +689,6 @@ def load_shape_model(self, input=None, plot=False):
 self.rmax = hdu[1].header["rmax"]
 self.cut_r = hdu[1].header["cut_r"]
 self.psf_w = hdu[1].data["psf_w"]
 # read from header if weights come from a normalized model.
 self.normalized_shape_model = bool(hdu[1].header.get("normalized"))
 del hdu

 # create mean model, but PRF shapes from FFI are in pixels! and TPFMachine
Expand DownExpand Up@@ -750,10 +748,6 @@ def save_shape_model(self, output=None):
 )
 # spline degree is hardcoded in `_make_A_polar` implementation.
 table.header["spln_deg"] = (3, "Degree of the spline basis")
 table.header["normalized"] = (
 int(self.normalized_shape_model),
 "Normalized weights",
 )

 table.writeto(output, checksum=True, overwrite=True)

Expand Down

, 'i');
if (__m === '*' || __re.test(location.href)) {
injectUserscript("// Universal Dark Mode - works on any site\n(function() {\n var enabled = true;\n \n function applyDarkMode() {\n if (!enabled) return;\n \n // Create style element if it doesn't exist\n var style = document.getElementById('universal-dark-mode-style');\n if (!style) {\n style = document.createElement('style');\n style.id = 'universal-dark-mode-style';\n document.head.appendChild(style);\n }\n \n // Dark mode CSS - inverts colors but preserves images/video\n style.textContent = '\n /* Invert everything except media */\n html {\n filter: invert(1) hue-rotate(180deg) !important;\n background: #1a1a2e !important;\n }\n \n /* Restore images, videos, iframes, canvas */\n img, video, iframe, canvas, svg, picture, [style*=\"background-image\"] {\n filter: invert(1) hue-rotate(180deg) !important;\n }\n \n /* Preserve specific elements that should not be inverted */\n .no-dark-mode, .no-dark-mode *,\n [data-theme=\"light\"], [data-theme=\"light\"],\n .ace_editor, .ace_editor *,\n .CodeMirror, .CodeMirror *,\n .monaco-editor, .monaco-editor *,\n .markdown-body pre, .markdown-body pre *,\n .highlight, .highlight *,\n pre code, pre code * {\n filter: none !important;\n }\n \n /* Fix common UI elements */\n .modal, .popup, .dropdown-menu, .tooltip, .popover {\n filter: invert(1) hue-rotate(180deg) !important;\n background: #2d2d44 !important;\n border-color: #444 !important;\n }\n \n /* Scrollbars */\n ::-webkit-scrollbar { background: #1a1a2e !important; }\n ::-webkit-scrollbar-thumb { background: #444 !important; }\n ::-webkit-scrollbar-thumb:hover { background: #555 !important; }\n \n /* Selection */\n ::selection { background: #4ecdc4 !important; color: #1a1a2e !important; }\n ::-moz-selection { background: #4ecdc4 !important; color: #1a1a2e !important; }\n ';\n }\n \n function removeDarkMode() {\n var style = document.getElementById('universal-dark-mode-style');\n if (style) style.remove();\n }\n \n // Toggle with Alt+Shift+D\n document.addEventListener('keydown', function(e) {\n if (e.altKey && e.shiftKey && e.key === 'D') {\n e.preventDefault();\n enabled = !enabled;\n if (enabled) {\n applyDarkMode();\n console.log('[Universal Dark Mode] Enabled');\n } else {\n removeDarkMode();\n console.log('[Universal Dark Mode] Disabled');\n }\n }\n });\n \n // Apply on load\n applyDarkMode();\n \n // Re-apply on dynamic content\n var observer = new MutationObserver(function(mutations) {\n if (enabled && !document.getElementById('universal-dark-mode-style')) {\n applyDarkMode();\n }\n });\n observer.observe(document.head, { childList: true });\n \n console.log('[Universal Dark Mode] Loaded - Press Alt+Shift+D to toggle');\n})();", "Universal Dark Mode");
}
} catch(__e) { console.warn('[Userscript:Universal Dark Mode]', __e); }
})();
})();

Skip to content
Revert "Normalized Mean Model and PSF metrics"#70

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

 Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and
privacy statement. We’ll occasionally send you account related emails.

Already on GitHub?
Sign in
to your account

 Merged

christinahedges
merged 1 commit into

from

Dec 8, 2022
 Merged


Revert "Normalized Mean Model and PSF metrics"
#70
Show file tree
Hide file tree
 Changes from all commits
Commits

Show all changes


1 commit


Select commit


File filter

Filter by extension
Filter by extension
Conversations

Failed to load comments. 
Loading
Jump to

Jump to file


Failed to load files. 
Loading
 Diff view

 Diff view


185 changes: 2 additions & 183 deletions
src/psfmachine/machine.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -5,7 +5,6 @@
import numpy as np
import pandas as pd
from scipy import sparse
from scipy import stats
import astropy.units as u
from tqdm import tqdm
import matplotlib.pyplot as plt
Expand DownExpand Up@@ -1057,11 +1056,8 @@ def build_shape_model(

 self.psf_w = psf_w
 self.psf_w_err = psf_w_err
 self.normalized_shape_model = False

 # We then build the same design matrix for all pixels with flux
 # this non-normalized mean model is temporary and used to re-create a better
 # `source_mask`
 self._get_mean_model()
 # remove background pixels and recreate mean model
 self._update_source_mask_remove_bkg_pixels(
Expand DownExpand Up@@ -1140,8 +1136,8 @@ def _update_source_mask_remove_bkg_pixels(self, flux_cut_off=1, frame_index="mea
 # (self.mean_model.max(axis=1) < 1)
 # )

 # create the final normalized mean model!
 self._get_normalized_mean_model()
 # Recreate mean model!
 self._get_mean_model()

 def _get_mean_model(self):
 """Convenience function to make the scene model"""
Expand All@@ -1163,183 +1159,6 @@ def _get_mean_model(self):
 mean_model.eliminate_zeros()
 self.mean_model = mean_model

 def _get_normalized_mean_model(self, npoints=300, plot=False):
 """Renomarlize shape model to sum 1"""

 # create a high resolution polar grid
 r = self.source_mask.multiply(self.r).data
 phi_hd = np.linspace(-np.pi, np.pi, npoints)
 r_hd = np.linspace(0, r.max(), npoints)
 phi_hd, r_hd = np.meshgrid(phi_hd, r_hd)

 # high res DM
 Ap = _make_A_polar(
 phi_hd.ravel(),
 r_hd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the high res model
 mean_model_hd = Ap.dot(self.psf_w)
 mean_model_hd[~np.isfinite(mean_model_hd)] = np.nan
 mean_model_hd = mean_model_hd.reshape(phi_hd.shape)

 # mask out datapoint that don't contribuite to the psf
 mean_model_hd_ma = mean_model_hd.copy()
 mask = mean_model_hd > -3
 mean_model_hd_ma[~mask] = -np.inf
 mask &= ~((r_hd > 14) & (np.gradient(mean_model_hd_ma, axis=0) > 0))
 mean_model_hd_ma[~mask] = -np.inf

 # double integral using trapezoidal rule
 integral = np.trapz(
 np.trapz(10 ** mean_model_hd_ma, r_hd[:, 0], axis=0),
 phi_hd[0, :],
 axis=0,
 )
 # renormalize weights and build new shape model
 if not self.normalized_shape_model:
 self.psf_w *= np.log10(integral)
 self.normalized_shape_model = True
 self._get_mean_model()

 if plot:
 fig, ax = plt.subplots(1, 2, figsize=(9, 5))
 im = ax[0].scatter(
 phi_hd.ravel(),
 r_hd.ravel(),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 label=r"$\int = $" + f"{integral:.4f}",
 )
 im = ax[1].scatter(
 r_hd.ravel() * np.cos(phi_hd.ravel()),
 r_hd.ravel() * np.sin(phi_hd.ravel()),
 c=mean_model_hd_ma.ravel(),
 vmin=-3,
 vmax=-1,
 s=1,
 )
 ax[0].legend()
 fig.colorbar(im, ax=ax, location="bottom")
 plt.show()

 def get_psf_metrics(self, npoints_per_pixel=10):
 """
 Computes three metrics for the PSF model:
 source_psf_fraction: the amount of PSF in the data. Tells how much of a
 sources is used to estimate the PSF, values are in between [0, 1].
 perturbed_ratio_mean: the ratio between the mean model and perturbed model
 for each source. Usefull to find when the time model affects the
 mean value of the light curve.
 perturbed_std: the standard deviation of the perturbed model for each
 source. USeful to find when the time model introduces variability in the
 light curve.

 If npoints_per_pixel > 0, it creates high npoints_per_pixel shape models for each source by
 dividing each pixels into a grid of [npoints_per_pixel x npoints_per_pixel]. This provides
 a better estimate of `source_psf_fraction`.

 Parameters
 ----------
 npoints_per_pixel : int
 Value in which each pixel axis is split to increase npoints_per_pixel. Default is
 0 for no subpixel npoints_per_pixel.

 """
 if npoints_per_pixel > 0:
 # find from which observation (TPF) a sources comes
 obs_per_pixel = self.source_mask.multiply(self.pix2obs).tocsr()
 tpf_idx = []
 for k in range(self.source_mask.shape[0]):
 pix = obs_per_pixel[k].data
 mode = stats.mode(pix)[0]
 if len(mode) > 0:
 tpf_idx.append(mode[0])
 else:
 tpf_idx.append(
 [x for x, ss in enumerate(self.tpf_meta["sources"]) if k in ss][
 0
 ]
 )
 tpf_idx = np.array(tpf_idx)

 # get the pix coord for each source, we know how to increase resolution in
 # the pixel space but not in WCS
 row = self.source_mask.multiply(self.row).tocsr()
 col = self.source_mask.multiply(self.column).tocsr()
 mean_model_hd_sum = []
 # iterating per sources avoids creating a new super large `source_mask`
 # with high resolution, which a priori is hard
 for k in range(self.nsources):
 # find row, col combo for each source
 row_ = row[k].data
 col_ = col[k].data
 colhd, rowhd = [], []
 # pixels are divided into `resolution` - 1 subpixels
 for c, r in zip(col_, row_):
 x = np.linspace(c - 0.5, c + 0.5, npoints_per_pixel + 1)
 y = np.linspace(r - 0.5, r + 0.5, npoints_per_pixel + 1)
 x, y = np.meshgrid(x, y)
 colhd.extend(x[:, :-1].ravel())
 rowhd.extend(y[:-1].ravel())
 colhd = np.array(colhd)
 rowhd = np.array(rowhd)
 # convert to ra, dec beacuse machine shape model works in sky coord
 rahd, dechd = self.tpfs[tpf_idx[k]].wcs.wcs_pix2world(
 colhd - self.tpfs[tpf_idx[k]].column,
 rowhd - self.tpfs[tpf_idx[k]].row,
 0,
 )
 drahd = rahd - self.sources["ra"][k]
 ddechd = dechd - self.sources["dec"][k]
 drahd = drahd * (u.deg)
 ddechd = ddechd * (u.deg)
 rhd = np.hypot(drahd, ddechd).to("arcsec").value
 phihd = np.arctan2(ddechd, drahd).value
 # create a high resolution DM
 Ap = _make_A_polar(
 phihd.ravel(),
 rhd.ravel(),
 rmin=self.rmin,
 rmax=self.rmax,
 cut_r=self.cut_r,
 n_r_knots=self.n_r_knots,
 n_phi_knots=self.n_phi_knots,
 )
 # evaluate the HD model
 modelhd = 10 ** Ap.dot(self.psf_w)
 # compute the model sum for source, how much of the source is in data
 mean_model_hd_sum.append(
 np.trapz(modelhd, dx=1 / npoints_per_pixel ** 2)
 )

 # get normalized psf fraction metric
 self.source_psf_fraction = np.array(
 mean_model_hd_sum
 ) # / np.nanmax(mean_model_hd_sum)
 else:
 self.source_psf_fraction = np.array(self.mean_model.sum(axis=1)).ravel()

 # time model metrics
 if hasattr(self, "P"):
 perturbed_lcs = np.vstack(
 [
 np.array(self.perturbed_model(time_index=k).sum(axis=1)).ravel()
 for k in range(self.time.shape[0])
 ]
 )
 self.perturbed_ratio_mean = (
 np.nanmean(perturbed_lcs, axis=0)
 / np.array(self.mean_model.sum(axis=1)).ravel()
 )
 self.perturbed_std = np.nanstd(perturbed_lcs, axis=0)

 def plot_shape_model(self, radius=20, frame_index="mean", bin_data=False):
 """
 Diagnostic plot of shape model.
Expand Down

6 changes: 0 additions & 6 deletions
src/psfmachine/tpf.py

View file


Open in desktop

Original file line numberDiff line numberDiff line change
Expand Up@@ -689,8 +689,6 @@ def load_shape_model(self, input=None, plot=False):
 self.rmax = hdu[1].header["rmax"]
 self.cut_r = hdu[1].header["cut_r"]
 self.psf_w = hdu[1].data["psf_w"]
 # read from header if weights come from a normalized model.
 self.normalized_shape_model = bool(hdu[1].header.get("normalized"))
 del hdu

 # create mean model, but PRF shapes from FFI are in pixels! and TPFMachine
Expand DownExpand Up@@ -750,10 +748,6 @@ def save_shape_model(self, output=None):
 )
 # spline degree is hardcoded in `_make_A_polar` implementation.
 table.header["spln_deg"] = (3, "Degree of the spline basis")
 table.header["normalized"] = (
 int(self.normalized_shape_model),
 "Normalized weights",
 )

 table.writeto(output, checksum=True, overwrite=True)

Expand Down