Skip to content

Target completeness caching #378

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

Draft
maxwell-zweig wants to merge 6 commits into
base: master
Choose a base branch
from
Draft

Target completeness caching #378

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
3ee22ce
target completeness caching
Jun 16, 2024
3bfba44
cache target completeness fix
Jun 16, 2024
2edf47b
blackening code
Jun 16, 2024
297e68e
changed integer completness filename variable name
Jun 20, 2024
1a7ab39
blackened code and removed a comment
Jun 20, 2024
762669e
removed test reports folder
Jun 20, 2024
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
164 changes: 108 additions & 56 deletions EXOSIMS/Completeness/GarrettCompleteness.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,9 @@ def __init__(self, order_of_quadrature=15, **specs):

self._outspec["order_of_quadrature"] = self.order_of_quadrature

# setting a base filename that isn't mutated by a completeness calculation
self.filename_intcompleteness = self.filename

def completeness_setup(self):
"""Preform any preliminary calculations needed for this flavor of completeness

Expand Down Expand Up @@ -206,20 +209,10 @@ def target_completeness(self, TL):
int_comp: 1D numpy array of completeness values for each target star

"""

OS = TL.OpticalSystem
if TL.calc_char_int_comp:
mode = list(
filter(lambda mode: "spec" in mode["inst"]["name"], OS.observingModes)
)[0]
else:
mode = list(filter(lambda mode: mode["detectionMode"], OS.observingModes))[
0
]

# To limit the amount of computation, we want to find the most common
# int_dMag value (typically the one the user sets in the input json since
# int_dMag is either the user input or the intCutoff_dMag).

vals, counts = np.unique(TL.int_dMag, return_counts=True)
self.mode_dMag = vals[np.argwhere(counts == np.max(counts))[0][0]]
mode_dMag_mask = TL.int_dMag == self.mode_dMag
Expand All @@ -234,60 +227,119 @@ def target_completeness(self, TL):
):
extstr += "%s: " % att + str(getattr(self.PlanetPopulation, att)) + " "
# include mode_dMag and intCutoff_dMag

extstr += (
"%s: " % "mode_dMag"
+ str(self.mode_dMag)
+ f"intCutoff_dMag: {TL.intCutoff_dMag}"
+ " "
)
ext = hashlib.md5(extstr.encode("utf-8")).hexdigest()
self.filename += ext
Cpath = os.path.join(self.cachedir, self.filename + ".acomp")

# calculate separations based on IWA
IWA = mode["IWA"]
OWA = mode["OWA"]
smin = (np.tan(IWA) * TL.dist).to("AU").value
if np.isinf(OWA):
smax = np.array([self.rmax] * len(smin))
else:
smax = (np.tan(OWA) * TL.dist).to("AU").value
smax[smax > self.rmax] = self.rmax

int_comp = np.zeros(smin.shape)
# calculate dMags based on maximum dMag
if self.PlanetPopulation.scaleOrbits:
L = np.where(TL.L > 0, TL.L, 1e-10) # take care of zero/negative values
smin = smin / np.sqrt(L)
smax = smax / np.sqrt(L)
dMag_vals = TL.int_dMag - 2.5 * np.log10(L)
separation_mask = smin < self.rmax
int_comp[separation_mask] = self.comp_s(
smin[separation_mask], smax[separation_mask], dMag_vals[separation_mask]
)
Cpath = os.path.join(
self.cachedir,
self.filename_intcompleteness
+ hashlib.md5("int_comp".encode("utf-8")).hexdigest()
+ ext
+ ".acomp",
)

if os.path.exists(Cpath):
self.vprint('Loading cached int completeness file from "%s".' % Cpath)
try:
with open(Cpath, "rb") as ff:
int_comp = pickle.load(ff)
except UnicodeDecodeError:
with open(Cpath, "rb") as ff:
int_comp = pickle.load(ff, encoding="latin1")
self.vprint("int Completeness loaded from cache.")
else:
# In this case we find where the mode dMag value is also in the
# separation range and use the vectorized integral since they have
# the same dMag value. Where the dMag values are not the mode we
# must use comp_s which is slower
dMag_vals = TL.int_dMag
separation_mask = smin < self.rmax
dist_s = self.genComp(Cpath, TL)
dist_sv = np.vectorize(dist_s.integral, otypes=[np.float64])
separation_mode_mask = separation_mask & mode_dMag_mask
separation_not_mode_mask = separation_mask & ~mode_dMag_mask
int_comp[separation_mode_mask] = dist_sv(
smin[separation_mode_mask], smax[separation_mode_mask]
)
int_comp[separation_not_mode_mask] = self.comp_s(
smin[separation_not_mode_mask],
smax[separation_not_mode_mask],
dMag_vals[separation_not_mode_mask],

self.vprint('Cached int completeness file not found at "%s".' % Cpath)
self.vprint("Beginning int_comp completeness calculations.")

OS = TL.OpticalSystem
if TL.calc_char_int_comp:
mode = list(
filter(
lambda mode: "spec" in mode["inst"]["name"], OS.observingModes
)
)[0]
else:
mode = list(
filter(lambda mode: mode["detectionMode"], OS.observingModes)
)[0]

# important PlanetPopulation attributes
atts = list(self.PlanetPopulation.__dict__)
extstr = ""
for att in sorted(atts, key=str.lower):
if (
not callable(getattr(self.PlanetPopulation, att))
and att != "PlanetPhysicalModel"
):
extstr += (
"%s: " % att + str(getattr(self.PlanetPopulation, att)) + " "
)
# include mode_dMag and intCutoff_dMag
extstr += (
"%s: " % "mode_dMag"
+ str(self.mode_dMag)
+ f"intCutoff_dMag: {TL.intCutoff_dMag}"
+ " "
)

# ensure that completeness values are between 0 and 1
int_comp = np.clip(int_comp, 0.0, 1.0)
ext = hashlib.md5(extstr.encode("utf-8")).hexdigest()
Cpath_aux = os.path.join(self.cachedir, self.filename + ext + ".acomp")
self.filename += ext

# calculate separations based on IWA
IWA = mode["IWA"]
OWA = mode["OWA"]
smin = (np.tan(IWA) * TL.dist).to("AU").value
if np.isinf(OWA):
smax = np.array([self.rmax] * len(smin))
else:
smax = (np.tan(OWA) * TL.dist).to("AU").value
smax[smax > self.rmax] = self.rmax

int_comp = np.zeros(smin.shape)
# calculate dMags based on maximum dMag
if self.PlanetPopulation.scaleOrbits:
L = np.where(TL.L > 0, TL.L, 1e-10) # take care of zero/negative values
smin = smin / np.sqrt(L)
smax = smax / np.sqrt(L)
dMag_vals = TL.int_dMag - 2.5 * np.log10(L)
separation_mask = smin < self.rmax
int_comp[separation_mask] = self.comp_s(
smin[separation_mask],
smax[separation_mask],
dMag_vals[separation_mask],
)
else:
# In this case we find where the mode dMag value is also in the
# separation range and use the vectorized integral since they have
# the same dMag value. Where the dMag values are not the mode we
# must use comp_s which is slower
dMag_vals = TL.int_dMag
separation_mask = smin < self.rmax
dist_s = self.genComp(Cpath_aux, TL)
dist_sv = np.vectorize(dist_s.integral, otypes=[np.float64])
separation_mode_mask = separation_mask & mode_dMag_mask
separation_not_mode_mask = separation_mask & ~mode_dMag_mask
int_comp[separation_mode_mask] = dist_sv(
smin[separation_mode_mask], smax[separation_mode_mask]
)
int_comp[separation_not_mode_mask] = self.comp_s(
smin[separation_not_mode_mask],
smax[separation_not_mode_mask],
dMag_vals[separation_not_mode_mask],
)

# ensure that completeness values are between 0 and 1
int_comp = np.clip(int_comp, 0.0, 1.0)
with open(Cpath, "wb") as ff:
pickle.dump(int_comp, ff)
self.vprint("Int Completeness data stored in %s" % Cpath)
return int_comp

def genComp(self, Cpath, TL):
Expand All @@ -308,7 +360,7 @@ def genComp(self, Cpath, TL):

if os.path.exists(Cpath):
# dist_s interpolant already exists for parameters
self.vprint("Loading cached completeness file from %s" % Cpath)
self.vprint("Loading cached completeness pdf file from %s" % Cpath)
try:
with open(Cpath, "rb") as ff:
H = pickle.load(ff)
Expand All @@ -334,7 +386,7 @@ def genComp(self, Cpath, TL):
H = {"dist_s": dist_s}
with open(Cpath, "wb") as ff:
pickle.dump(H, ff)
self.vprint("Completeness data stored in %s" % Cpath)
self.vprint("Completeness pdf data stored in %s" % Cpath)

return dist_s

Expand Down
Loading