From 5e20030b8d7be9e96836c2f49a33211a12a0fb55 Mon Sep 17 00:00:00 2001 From: Kyle A Pearson Date: Wed, 15 Dec 2021 09:23:33 -0800 Subject: [PATCH 01/17] flux weighted centroid --- exotic/exotic.py | 60 ++++++++++++++++++++++++++++++++++++------------ 1 file changed, 45 insertions(+), 15 deletions(-) diff --git a/exotic/exotic.py b/exotic/exotic.py index 785d534a..a2346ade 100644 --- a/exotic/exotic.py +++ b/exotic/exotic.py @@ -90,7 +90,6 @@ from scipy.stats import mode from scipy.signal import savgol_filter from scipy.ndimage import binary_erosion -# from scipy.ndimage import binary_dilation, label, generic_filter from skimage.transform import SimilarityTransform # error handling for scraper from tenacity import retry, stop_after_delay @@ -915,12 +914,16 @@ def mesh_box(pos, box, maxx=0, maxy=0): # Method fits a 2D gaussian function that matches the star_psf to the star image and returns its pixel coordinates -def fit_centroid(data, pos, psf_function=gaussian_psf, box=10): +def fit_centroid(data, pos, psf_function=gaussian_psf, box=15, weightedcenter=True): # get sub field in image xv, yv = mesh_box(pos, box, maxx=data.shape[1], maxy=data.shape[0]) subarray = data[yv, xv] init = [np.nanmax(subarray) - np.nanmin(subarray), 1, 1, 0, np.nanmin(subarray)] + # compute flux weighted centroid in x and y + wx = np.sum(xv[0]*subarray.sum(0))/subarray.sum(0).sum() + wy = np.sum(yv[:,0]*subarray.sum(1))/subarray.sum(1).sum() + # lower bound: [xc, yc, amp, sigx, sigy, rotation, bg] lo = [pos[0] - box * 0.5, pos[1] - box * 0.5, 0, 0.5, 0.5, -np.pi / 4, np.nanmin(subarray) - 1] up = [pos[0] + box * 0.5, pos[1] + box * 0.5, 1e7, 20, 20, np.pi / 4, np.nanmax(subarray) + 1] @@ -936,6 +939,10 @@ def fcn2min(pars): res = least_squares(fcn2min, x0=[*pos, *init], jac='3-point', xtol=None, method='lm') + if weightedcenter: + res.x[0] = wx + res.x[1] = wy + return res.x @@ -1188,7 +1195,9 @@ def fit_lightcurve(times, tFlux, cFlux, airmass, ld, pDict): si = np.argsort(times) dt = np.mean(np.diff(np.sort(times))) ndt = int(25. / 24. / 60. / dt) * 2 + 1 - filtered_data = sigma_clip((tFlux / cFlux)[si], sigma=3, dt=ndt) + if ndt > len(times): + ndt = int(len(times)/4) * 2 + 1 + filtered_data = sigma_clip((tFlux / cFlux)[si], sigma=3, dt=max(5,ndt)) arrayFinalFlux = (tFlux / cFlux)[si][~filtered_data] f1 = tFlux[si][~filtered_data] sigf1 = f1 ** 0.5 @@ -1206,10 +1215,15 @@ def fit_lightcurve(times, tFlux, cFlux, airmass, ld, pDict): arrayAirmass) | np.less_equal(arrayFinalFlux, 0) | np.less_equal(arrayNormUnc, 0) nanmask = nanmask | np.isinf(arrayFinalFlux) | np.isinf(arrayNormUnc) | np.isinf(arrayTimes) | np.isinf( arrayAirmass) - arrayFinalFlux = arrayFinalFlux[~nanmask] - arrayNormUnc = arrayNormUnc[~nanmask] - arrayTimes = arrayTimes[~nanmask] - arrayAirmass = arrayAirmass[~nanmask] + + if np.sum(~nanmask) <= 1: + log_info('No data left after filtering', warn=True) + else: + arrayFinalFlux = arrayFinalFlux[~nanmask] + arrayNormUnc = arrayNormUnc[~nanmask] + arrayTimes = arrayTimes[~nanmask] + arrayAirmass = arrayAirmass[~nanmask] + # -----LM LIGHTCURVE FIT-------------------------------------- prior = { @@ -1247,6 +1261,9 @@ def fit_lightcurve(times, tFlux, cFlux, airmass, ld, pDict): 'a2': [-1, 1] } + if np.isnan(arrayTimes).any() or np.isnan(arrayFinalFlux).any() or np.isnan(arrayNormUnc).any(): + log_info("\nWarning: NANs in time, flux or error", warn=True) + myfit = lc_fitter( arrayTimes, arrayFinalFlux, @@ -1920,14 +1937,21 @@ def main(): # Calculate the proper timeseries uncertainties from the residuals of the out-of-transit data OOT = (bestlmfit.transit == 1) # find out-of-transit portion of the lightcurve - if len(OOT) == 0: # if user does not get any out-of-transit data, normalize by the max data instead - OOT = (bestlmfit.transit == np.nanmax(bestlmfit.transit)) - OOTscatter = np.std((bestlmfit.data / bestlmfit.airmass_model)[OOT]) # calculate the scatter in the data - goodNormUnc = OOTscatter * bestlmfit.airmass_model # scale this scatter back up by the airmass model and then adopt these as the uncertainties - # Normalize by OOT per AAVSO Database upload requirements - goodNormUnc = goodNormUnc / np.nanmedian(goodFluxes[OOT]) - goodFluxes = goodFluxes / np.nanmedian(goodFluxes[OOT]) + if sum(OOT) <= 1: + OOTscatter = np.std(bestlmfit.residuals) + goodNormUnc = OOTscatter * bestlmfit.airmass_model + goodNormUnc = goodNormUnc / np.nanmedian(goodFluxes) + goodFluxes = goodFluxes / np.nanmedian(goodFluxes) + else: + OOTscatter = np.std((bestlmfit.data / bestlmfit.airmass_model)[OOT]) # calculate the scatter in the data + goodNormUnc = OOTscatter * bestlmfit.airmass_model # scale this scatter back up by the airmass model and then adopt these as the uncertainties + goodNormUnc = goodNormUnc / np.nanmedian(goodFluxes[OOT]) + goodFluxes = goodFluxes / np.nanmedian(goodFluxes[OOT]) + + if np.isnan(goodFluxes).all(): + log_info("Error: No valid photometry data found.", error=True) + return goodTimes = goodTimes[si][gi] goodFluxes = goodFluxes[si][gi] @@ -1986,11 +2010,13 @@ def main(): goodTimes = timeConvert(goodTimes, exotic_infoDict['file_time'], pDict, exotic_infoDict) if exotic_infoDict['file_units'] != 'flux': + print("check flux convert") + import pdb; pdb.set_trace() goodFluxes, goodNormUnc = fluxConvert(goodFluxes, goodNormUnc, exotic_infoDict['file_units']) # for k in myfit.bounds.keys(): # print(f"{myfit.parameters[k]:.6f} +- {myfit.errors[k]}") - + if args.photometry: log_info("\nPhotometric Extraction Complete.") return @@ -2034,6 +2060,10 @@ def main(): 'a2': [-3, 3], } + if np.isnan(goodFluxes).all(): + log_info("Error: No valid photometry data found.", error=True) + return + # final light curve fit myfit = lc_fitter(goodTimes, goodFluxes, goodNormUnc, goodAirmasses, prior, mybounds, mode='ns') # myfit.dataerr *= np.sqrt(myfit.chi2 / myfit.data.shape[0]) # scale errorbars by sqrt(rchi2) From 22ef662edf24af2dab04b9c1f68bb3489ee7aa25 Mon Sep 17 00:00:00 2001 From: Kyle A Pearson Date: Wed, 15 Dec 2021 09:24:38 -0800 Subject: [PATCH 02/17] flux weighted centroid --- exotic/exotic.py | 1 + 1 file changed, 1 insertion(+) diff --git a/exotic/exotic.py b/exotic/exotic.py index a2346ade..aac3c844 100644 --- a/exotic/exotic.py +++ b/exotic/exotic.py @@ -939,6 +939,7 @@ def fcn2min(pars): res = least_squares(fcn2min, x0=[*pos, *init], jac='3-point', xtol=None, method='lm') + # override psf fit results with weighted centroid if weightedcenter: res.x[0] = wx res.x[1] = wy From 5eb2d4930cefeacb9d0c11025fe4f0dde41c6f3b Mon Sep 17 00:00:00 2001 From: Kyle A Pearson Date: Wed, 15 Dec 2021 09:25:38 -0800 Subject: [PATCH 03/17] flux weighted centroid --- exotic/exotic.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/exotic/exotic.py b/exotic/exotic.py index aac3c844..919d6711 100644 --- a/exotic/exotic.py +++ b/exotic/exotic.py @@ -2011,8 +2011,6 @@ def main(): goodTimes = timeConvert(goodTimes, exotic_infoDict['file_time'], pDict, exotic_infoDict) if exotic_infoDict['file_units'] != 'flux': - print("check flux convert") - import pdb; pdb.set_trace() goodFluxes, goodNormUnc = fluxConvert(goodFluxes, goodNormUnc, exotic_infoDict['file_units']) # for k in myfit.bounds.keys(): From 36962e9d8183d66703a562443d3e17bb73f1b1c8 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Fri, 17 Dec 2021 17:16:41 +0000 Subject: [PATCH 04/17] Update cython requirement from ~=0.29.25 to ~=0.29.26 Updates the requirements on [cython](https://github.com/cython/cython) to permit the latest version. - [Release notes](https://github.com/cython/cython/releases) - [Changelog](https://github.com/cython/cython/blob/master/CHANGES.rst) - [Commits](https://github.com/cython/cython/compare/0.29.25...0.29.26) --- updated-dependencies: - dependency-name: cython dependency-type: direct:production ... Signed-off-by: dependabot[bot] --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 948e030a..dd49c76a 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,7 +7,7 @@ astroalign~=2.4 astropy>=4.3 astroquery~=0.4 barycorrpy~=0.4 -cython~=0.29.25 ; platform_system != "Windows" +cython~=0.29.26 ; platform_system != "Windows" dynesty~=1.1 ; platform_system == "Windows" holoviews~=1.14 LDTk~=1.7 From e7ba0de9643f384a42faeb5a0ca861584475b2cb Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 20 Dec 2021 17:18:10 +0000 Subject: [PATCH 05/17] Update numpy requirement from ~=1.20 to ~=1.21 Updates the requirements on [numpy](https://github.com/numpy/numpy) to permit the latest version. - [Release notes](https://github.com/numpy/numpy/releases) - [Changelog](https://github.com/numpy/numpy/blob/main/doc/HOWTO_RELEASE.rst.txt) - [Commits](https://github.com/numpy/numpy/compare/v1.20.0...v1.21.5) --- updated-dependencies: - dependency-name: numpy dependency-type: direct:production ... Signed-off-by: dependabot[bot] --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 948e030a..9ff0affb 100644 --- a/requirements.txt +++ b/requirements.txt @@ -13,7 +13,7 @@ holoviews~=1.14 LDTk~=1.7 lmfit~=1.0 matplotlib>=3.4 -numpy~=1.20 +numpy~=1.21 pandas~=1.3 panel~=0.12 photutils>=0.7 From 6ff71bd1118c02457cdbc41dc45a93c02a2117eb Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 3 Jan 2022 17:19:19 +0000 Subject: [PATCH 06/17] Update requests requirement from ~=2.26 to ~=2.27 Updates the requirements on [requests](https://github.com/psf/requests) to permit the latest version. - [Release notes](https://github.com/psf/requests/releases) - [Changelog](https://github.com/psf/requests/blob/main/HISTORY.md) - [Commits](https://github.com/psf/requests/compare/v2.26.0...v2.27.0) --- updated-dependencies: - dependency-name: requests dependency-type: direct:production ... Signed-off-by: dependabot[bot] --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 948e030a..a4f33923 100644 --- a/requirements.txt +++ b/requirements.txt @@ -19,7 +19,7 @@ panel~=0.12 photutils>=0.7 python_dateutil~=2.8 pyvo~=1.1 -requests~=2.26 +requests~=2.27 scipy~=1.7 scikit-image~=0.18 tenacity~=8.0 From ffa483648f925392b3766d474c31983a93caf163 Mon Sep 17 00:00:00 2001 From: Kyle A Pearson Date: Wed, 5 Jan 2022 14:23:42 -0800 Subject: [PATCH 07/17] new image alignment fallback with smoothing --- exotic/exotic.py | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/exotic/exotic.py b/exotic/exotic.py index 919d6711..6fe49872 100644 --- a/exotic/exotic.py +++ b/exotic/exotic.py @@ -90,6 +90,7 @@ from scipy.stats import mode from scipy.signal import savgol_filter from scipy.ndimage import binary_erosion +from skimage.util import view_as_windows from skimage.transform import SimilarityTransform # error handling for scraper from tenacity import retry, stop_after_delay @@ -776,6 +777,22 @@ def transformation(image_data, file_name, roi=1): except Exception as ee: log_info(ee) + ws = 5 + # smooth image and try to align again + windows = view_as_windows(image_data[0], (ws,ws), step=1) + medimg = np.median(windows, axis=(2,3)) + + windows = view_as_windows(image_data[1], (ws,ws), step=1) + medimg1 = np.median(windows, axis=(2,3)) + + try: + results = aa.find_transform(medimg1[roiy, roix], medimg[roiy, roix]) + return results[0] + except Exception as ee: + log_info(ee) + + log_info(file_name) + for p in [99, 98, 95, 90]: for it in [2, 1, 0]: @@ -1700,7 +1717,10 @@ def main(): if i == 0: tform = SimilarityTransform(scale=1, rotation=0, translation=[0, 0]) else: - tform = transformation(np.array([imageData, firstImage]), fileName) + try: + tform = transformation(np.array([imageData, firstImage]), fileName) + except: + pass tx, ty = tform([exotic_UIprevTPX, exotic_UIprevTPY])[0] psf_data['target'][i] = fit_centroid(imageData, [tx, ty]) @@ -2011,6 +2031,8 @@ def main(): goodTimes = timeConvert(goodTimes, exotic_infoDict['file_time'], pDict, exotic_infoDict) if exotic_infoDict['file_units'] != 'flux': + print("check flux convert") + import pdb; pdb.set_trace() goodFluxes, goodNormUnc = fluxConvert(goodFluxes, goodNormUnc, exotic_infoDict['file_units']) # for k in myfit.bounds.keys(): From c359fccff3655cf48b3d209394ed1919b413b6d8 Mon Sep 17 00:00:00 2001 From: Kyle A Pearson Date: Wed, 5 Jan 2022 14:24:53 -0800 Subject: [PATCH 08/17] new image alignment fallback with smoothing --- exotic/exotic.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/exotic/exotic.py b/exotic/exotic.py index 6fe49872..74595575 100644 --- a/exotic/exotic.py +++ b/exotic/exotic.py @@ -1717,10 +1717,7 @@ def main(): if i == 0: tform = SimilarityTransform(scale=1, rotation=0, translation=[0, 0]) else: - try: - tform = transformation(np.array([imageData, firstImage]), fileName) - except: - pass + tform = transformation(np.array([imageData, firstImage]), fileName) tx, ty = tform([exotic_UIprevTPX, exotic_UIprevTPY])[0] psf_data['target'][i] = fit_centroid(imageData, [tx, ty]) From 1dc152b284287d046901398478a019885480b4e9 Mon Sep 17 00:00:00 2001 From: Tamim Fatahi Date: Thu, 6 Jan 2022 19:11:39 -0800 Subject: [PATCH 09/17] Clarifying input for image scale and removing numba temporarily. --- exotic/api/elca.py | 4 ++-- exotic/exotic.py | 5 ++--- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/exotic/api/elca.py b/exotic/api/elca.py index c0bbab23..ec71e8fc 100644 --- a/exotic/api/elca.py +++ b/exotic/api/elca.py @@ -42,7 +42,7 @@ # ########################################################################### # import copy -from numba import njit +# from numba import njit import numpy as np import matplotlib.pyplot as plt from scipy.optimize import least_squares @@ -288,7 +288,7 @@ def loglike(pars): model *= np.median(detrend) return -0.5 * np.sum(((self.data - model) / self.dataerr) ** 2) - @njit(fastmath=True) + # @njit(fastmath=True) def prior_transform(upars): # transform unit cube to prior volume return boundarray[:, 0] + bounddiff * upars diff --git a/exotic/exotic.py b/exotic/exotic.py index 74595575..afd6d65a 100644 --- a/exotic/exotic.py +++ b/exotic/exotic.py @@ -81,7 +81,7 @@ from matplotlib.animation import FuncAnimation # Pyplot imports import matplotlib.pyplot as plt -from numba import njit +# from numba import njit import numpy as np # photometry from photutils import CircularAperture @@ -835,7 +835,7 @@ def get_pixel_scale(wcs_header, header, pixel_init): image_scale = f"Image scale in arcsecs/pixel: {pixel_init}" else: log_info("Not able to find Image Scale in the Image Header.") - image_scale_num = user_input("Please enter Image Scale (e.g., 5 arcsec/pixel): ", type_=float) + image_scale_num = user_input("Please enter Image Scale (units in arcsec/pixel): ", type_=float) image_scale = f"Image scale in arcsecs/pixel: {image_scale_num}" return image_scale @@ -907,7 +907,6 @@ def find_target(target, hdufile, verbose=False): return pixcoord[0] -@njit def gaussian_psf(x, y, x0, y0, a, sigx, sigy, rot, b): rx = (x - x0) * np.cos(rot) - (y - y0) * np.sin(rot) ry = (x - x0) * np.sin(rot) + (y - y0) * np.cos(rot) From 2b0aeae5e0e96b265d491af3ec5e048dc9a0fad6 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Fri, 7 Jan 2022 03:13:18 +0000 Subject: [PATCH 10/17] Update pyvo requirement from ~=1.1 to ~=1.2 Updates the requirements on [pyvo](https://github.com/astropy/pyvo) to permit the latest version. - [Release notes](https://github.com/astropy/pyvo/releases) - [Changelog](https://github.com/astropy/pyvo/blob/main/CHANGES.rst) - [Commits](https://github.com/astropy/pyvo/compare/v1.1...v1.2) --- updated-dependencies: - dependency-name: pyvo dependency-type: direct:production ... Signed-off-by: dependabot[bot] --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index a4f33923..b1f2eb02 100644 --- a/requirements.txt +++ b/requirements.txt @@ -18,7 +18,7 @@ pandas~=1.3 panel~=0.12 photutils>=0.7 python_dateutil~=2.8 -pyvo~=1.1 +pyvo~=1.2 requests~=2.27 scipy~=1.7 scikit-image~=0.18 From b672d9c13ac1ae5bcaeb9b22cc79958ed9c6546e Mon Sep 17 00:00:00 2001 From: Tamim Fatahi Date: Thu, 6 Jan 2022 19:21:10 -0800 Subject: [PATCH 11/17] Modifying output for image scale entry. --- exotic/exotic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exotic/exotic.py b/exotic/exotic.py index afd6d65a..71e97d55 100644 --- a/exotic/exotic.py +++ b/exotic/exotic.py @@ -835,7 +835,7 @@ def get_pixel_scale(wcs_header, header, pixel_init): image_scale = f"Image scale in arcsecs/pixel: {pixel_init}" else: log_info("Not able to find Image Scale in the Image Header.") - image_scale_num = user_input("Please enter Image Scale (units in arcsec/pixel): ", type_=float) + image_scale_num = user_input("Please enter Image Scale (arcsec/pixel): ", type_=float) image_scale = f"Image scale in arcsecs/pixel: {image_scale_num}" return image_scale From c4eb4ee2d0b5ffebf5bcdeb4daafaff93c92049f Mon Sep 17 00:00:00 2001 From: Kyle A Pearson Date: Thu, 13 Jan 2022 13:58:02 -0800 Subject: [PATCH 12/17] global fitter --- examples/global_fitter_unit_test.py | 73 +++++++++++++ exotic/api/elca.py | 162 ++++++++++++++++++++++++++++ 2 files changed, 235 insertions(+) create mode 100644 examples/global_fitter_unit_test.py diff --git a/examples/global_fitter_unit_test.py b/examples/global_fitter_unit_test.py new file mode 100644 index 00000000..10f263bb --- /dev/null +++ b/examples/global_fitter_unit_test.py @@ -0,0 +1,73 @@ +import numpy as np +import matplotlib.pyplot as plt + +from exotic.api.elca import transit, glc_fitter + +if __name__ == "__main__": + + # simulate input data + epochs = np.random.choice(np.arange(100), 5, replace=False) + input_data = [] + + for i, epoch in enumerate(epochs): + + nobs = np.random.randint(50) + 100 + phase = np.linspace(-0.02-0.01*np.random.random(), 0.02+0.01*np.random.random(), nobs) + + prior = { + 'rprs':0.1, # Rp/Rs + 'ars':14.25, # a/Rs + 'per':3.5, # Period [day] + 'inc':87.5, # Inclination [deg] + 'u0': 1.349, 'u1': -0.709, # exotethys - limb darkening (nonlinear) + 'u2': 0.362, 'u3': -0.087, + 'ecc':0, # Eccentricity + 'omega':0, # Arg of periastron + 'tmid':1, # time of mid transit [day], + + 'a1':5000 + 2500*np.random.random(), # airmass coeffcients + 'a2':-0.25 + 0.1*np.random.random() + } + + time = prior['tmid'] + prior['per']*(phase+epoch) + stime = time-time[0] + alt = 90* np.cos(4*stime-np.pi/6) + airmass = 1./np.cos( np.deg2rad(90-alt)) + model = transit(time, prior)*prior['a1']*np.exp(prior['a2']*airmass) + flux = model*np.random.normal(1, np.mean(np.sqrt(model)/model)*0.25, model.shape) + ferr = flux**0.5 + + input_data.append({ + 'time':time, + 'flux':flux, + 'ferr':ferr, + 'airmass':airmass, + 'priors':prior + }) + + #plt.plot(time,flux,marker='o') + #plt.plot(time, model,ls='-') + #plt.show() + + # shared properties between light curves + global_bounds = { + 'per':[3.5-0.0001,3.5+0.0001], + 'tmid':[1-0.01,1+0.01], + 'ars':[14,14.5], + } + + # individual properties + local_bounds = { + 'rprs':[0,0.2], + 'a1':[0,1e4], + 'a2':[-0.5,0] + } + + print('epochs:',epochs) + myfit = glc_fitter(input_data, global_bounds, local_bounds, individual_fit=True, verbose=True) + + myfit.plot_triangle() + plt.show() + + myfit.plot_bestfit() + plt.show() \ No newline at end of file diff --git a/exotic/api/elca.py b/exotic/api/elca.py index c0bbab23..40c83727 100644 --- a/exotic/api/elca.py +++ b/exotic/api/elca.py @@ -518,6 +518,168 @@ def plot_triangle(self): return fig +# simultaneously fit multiple data sets with global and local parameters +class glc_fitter(lc_fitter): + # needed for lc_fitter + ns_type = 'ultranest' + + def __init__(self, input_data, global_bounds, local_bounds, individual_fit=False, verbose=False): + # keys for input_data: time, flux, ferr, airmass, priors all numpy arrays + self.data = copy.deepcopy(input_data) + self.global_bounds = global_bounds + self.local_bounds = local_bounds + self.individual_fit = individual_fit + self.verbose = verbose + + self.fit_nested() + + def fit_nested(self): + + # create bound arrays for generating samples + nobs = len(self.data) + gfreekeys = list(self.global_bounds.keys()) + + if isinstance(self.local_bounds, dict): + lfreekeys = list(self.local_bounds.keys()) + boundarray = np.vstack([ [self.global_bounds[k] for k in gfreekeys], [self.local_bounds[k] for k in lfreekeys]*nobs ]) + else: + # if list type + lfreekeys = list(self.local_bounds[0].keys()) + boundarray = [self.global_bounds[k] for k in gfreekeys] + for i in range(nobs): + boundarray.extend([self.local_bounds[i][k] for k in lfreekeys]) + boundarray = np.array(boundarray) + + # fit individual light curves to constrain priors + if self.individual_fit: + for i in range(nobs): + + print(f"Fitting individual light curve {i+1}/{nobs}") + mybounds = dict(**self.local_bounds, **self.global_bounds) + if 'per' in mybounds: del(mybounds['per']) + + myfit = lc_fitter( + self.data[i]['time'], + self.data[i]['flux'], + self.data[i]['ferr'], + self.data[i]['airmass'], + self.data[i]['priors'], + mybounds + ) + + self.data[i]['individual'] = myfit.parameters.copy() + + # update local priors + for j, key in enumerate(self.local_bounds.keys()): + + boundarray[j+i*len(self.local_bounds)+len(gfreekeys),0] = myfit.parameters[key] - 5*myfit.errors[key] + boundarray[j+i*len(self.local_bounds)+len(gfreekeys),1] = myfit.parameters[key] + 5*myfit.errors[key] + if key == 'rprs': + boundarray[j+i*len(self.local_bounds)+len(gfreekeys),0] = max(0,myfit.parameters[key] - 5*myfit.errors[key]) + + del(myfit) + + # transform unit cube to prior volume + bounddiff = np.diff(boundarray,1).reshape(-1) + def prior_transform(upars): + return (boundarray[:,0] + bounddiff*upars) + + def loglike(pars): + chi2 = 0 + + # for each light curve + for i in range(nobs): + + # global keys + for j, key in enumerate(gfreekeys): + self.data[i]['priors'][key] = pars[j] + + # local keys + for j, key in enumerate(lfreekeys): + self.data[i]['priors'][key] = pars[j+i*len(lfreekeys)+len(gfreekeys)] + + # compute model + model = transit(self.data[i]['time'], self.data[i]['priors']) + model *= np.exp(self.data[i]['priors']['a2']*self.data[i]['airmass']) + detrend = self.data[i]['flux']/model + model *= np.median(detrend) + + chi2 += np.sum( ((self.data[i]['flux']-model)/self.data[i]['ferr'])**2 ) + + # maximization metric for nested sampling + return -0.5*chi2 + + freekeys = []+gfreekeys + for n in range(nobs): + for k in lfreekeys: + freekeys.append(f"local_{n}_{k}") + + if self.verbose: + self.results = ReactiveNestedSampler(freekeys, loglike, prior_transform).run(max_ncalls=4e5) + else: + self.results = ReactiveNestedSampler(freekeys, loglike, prior_transform).run(max_ncalls=4e5, show_status=self.verbose, viz_callback=self.verbose) + + self.errors = {} + self.quantiles = {} + self.parameters = {} + + for i, key in enumerate(freekeys): + + self.parameters[key] = self.results['maximum_likelihood']['point'][i] + self.errors[key] = self.results['posterior']['stdev'][i] + self.quantiles[key] = [ + self.results['posterior']['errlo'][i], + self.results['posterior']['errup'][i]] + + for n in range(nobs): + for k in lfreekeys: + pkey = f"local_{n}_{k}" + self.data[n]['priors'][k] = self.parameters[pkey] + + # solve for a1 + model = transit(self.data[n]['time'], self.data[n]['priors']) + airmass = np.exp(self.data[n]['airmass']*self.data[n]['priors']['a2']) + detrend = self.data[n]['flux']/(model*airmass) + self.data[n]['priors']['a1'] = np.median(detrend) + + def plot_bestfit(self): + # TODO mosaic + nobs = len(self.data) + nrows = nobs//4+1 + fig,ax = plt.subplots(nrows, 4, figsize=(5+5*nrows,4*nrows)) + + # turn off all axes + for i in range(nrows*4): + ri = int(i/4) + ci = i%4 + if ax.ndim == 1: + ax[i].axis('off') + else: + ax[ri,ci].axis('off') + + # plot observations + for i in range(nobs): + ri = int(i/4) + ci = i%4 + model = transit(self.data[i]['time'], self.data[i]['priors']) + airmass = np.exp(self.data[i]['airmass']*self.data[i]['priors']['a2']) + detrend = self.data[i]['flux']/(model*airmass) + + if ax.ndim == 1: + ax[i].axis('on') + ax[i].errorbar(self.data[i]['time'], self.data[i]['flux']/airmass/detrend.mean(), yerr=self.data[i]['ferr']/airmass/detrend.mean(), ls='none', marker='o', color='black', alpha=0.5) + ax[i].plot(self.data[i]['time'], model, 'r-') + ax[i].set_xlabel("Time") + + else: + ax[ri,ci].axis('on') + ax[ri,ci].errorbar(self.data[i]['time'], self.data[i]['flux']/airmass/detrend.mean(), yerr=self.data[i]['ferr']/airmass/detrend.mean(), ls='none', marker='o', color='black', alpha=0.5) + ax[ri,ci].plot(self.data[i]['time'], model, 'r-') + ax[ri,ci].set_xlabel("Time") + plt.tight_layout() + return fig + + if __name__ == "__main__": prior = { From 4ef746f406a530faa04d9276e4fdfd58b9dd121d Mon Sep 17 00:00:00 2001 From: Kyle A Pearson Date: Thu, 13 Jan 2022 14:09:45 -0800 Subject: [PATCH 13/17] additional lc fit example --- examples/lc_fitter_unit_test.py | 75 +++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 examples/lc_fitter_unit_test.py diff --git a/examples/lc_fitter_unit_test.py b/examples/lc_fitter_unit_test.py new file mode 100644 index 00000000..7778a20e --- /dev/null +++ b/examples/lc_fitter_unit_test.py @@ -0,0 +1,75 @@ +from exotic.exotic import LimbDarkening +from exotic.api.elca import transit, lc_fitter +from ldtk.filters import create_tess +import matplotlib.pyplot as plt +import numpy as np + +if __name__ == "__main__": + prior = { + 'rprs': 0.02, # Rp/Rs + 'ars': 14.25, # a/Rs + 'per': 3.33, # Period [day] + 'inc': 88.5, # Inclination [deg] + 'u0': 0, 'u1': 0, 'u2': 0, 'u3': 0, # limb darkening (nonlinear) + 'ecc': 0.5, # Eccentricity + 'omega': 120, # Arg of periastron + 'tmid': 0.75, # Time of mid transit [day], + 'a1': 50, # Airmass coefficients + 'a2': 0., # trend = a1 * np.exp(a2 * airmass) + + 'teff':5000, + 'tefferr':50, + 'met': 0, + 'meterr': 0, + 'logg': 3.89, + 'loggerr': 0.01 + } + + # example generating LD coefficients + tessfilter = create_tess() + + ld_obj = LimbDarkening( + teff=prior['teff'], teffpos=prior['tefferr'], teffneg=prior['tefferr'], + met=prior['met'], metpos=prior['meterr'], metneg=prior['meterr'], + logg=prior['logg'], loggpos=prior['loggerr'], loggneg=prior['loggerr'], + wl_min=tessfilter.wl.min(), wl_max=tessfilter.wl.max(), filter_type="Clear") + + ld0, ld1, ld2, ld3, filt, wlmin, wlmax = ld_obj.nonlinear_ld() + + prior['u0'],prior['u1'],prior['u2'],prior['u3'] = [ld0[0], ld1[0], ld2[0], ld3[0]] + + time = np.linspace(0.7, 0.8, 1000) # [day] + + # simulate extinction from airmass + stime = time-time[0] + alt = 90 * np.cos(4*stime-np.pi/6) + #airmass = 1./np.cos(np.deg2rad(90-alt)) + airmass = np.zeros(time.shape[0]) + + # GENERATE NOISY DATA + data = transit(time, prior)*prior['a1']*np.exp(prior['a2']*airmass) + data += np.random.normal(0, prior['a1']*250e-6, len(time)) + dataerr = np.random.normal(300e-6, 50e-6, len(time)) + np.random.normal(300e-6, 50e-6, len(time)) + + # add bounds for free parameters only + mybounds = { + 'rprs': [0, 0.1], + 'tmid': [prior['tmid']-0.01, prior['tmid']+0.01], + 'ars': [13, 15], + #'a2': [0, 0.3] # uncomment if you want to fit for airmass + # never list 'a1' in bounds, it is perfectly correlated to exp(a2*airmass) + # and is solved for during the fit + } + + myfit = lc_fitter(time, data, dataerr, airmass, prior, mybounds, mode='ns') + + for k in myfit.bounds.keys(): + print(f"{myfit.parameters[k]:.6f} +- {myfit.errors[k]}") + + fig, axs = myfit.plot_bestfit() + plt.tight_layout() + plt.show() + + fig = myfit.plot_triangle() + plt.tight_layout() + plt.show() \ No newline at end of file From e37ea49d1e3d84dcf9cb29473f8d139485b67363 Mon Sep 17 00:00:00 2001 From: Ayush Nayak Date: Mon, 17 Jan 2022 16:38:18 -0800 Subject: [PATCH 14/17] Allow for window destroy events. This is just a quality of life improvement, allowing for the X button to work properly in GUI versions of EXOTIC. --- exotic/exotic_gui.py | 63 +++++++++++++++++++++++++++++--------------- 1 file changed, 42 insertions(+), 21 deletions(-) diff --git a/exotic/exotic_gui.py b/exotic/exotic_gui.py index 828b8d98..691b2de2 100644 --- a/exotic/exotic_gui.py +++ b/exotic/exotic_gui.py @@ -133,7 +133,9 @@ def main(): else: print(f"\nSUCCESS: Valid Python version {platform.python_version()} detected!\n") - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) + root.title(f"EXOTIC v{__version__}") tk.Label(root, @@ -165,7 +167,8 @@ def main(): root.mainloop() - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") reduction_opt = tk.IntVar() @@ -198,7 +201,8 @@ def main(): if reduction_opt.get() == 1: # First ask user how they want to enter the observing information - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) obsinfo = tk.StringVar() input_data = {} @@ -231,7 +235,8 @@ def main(): root.mainloop() if obsinfo.get() == "manual": - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") window_label = tk.Label(root, @@ -287,7 +292,8 @@ def save_input(): tk.mainloop() else: - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") # # "Directory with FITS files": "sample-data/HatP32Dec202017", @@ -306,7 +312,8 @@ def save_input(): tk.mainloop() if obsinfo.get() == "manual": - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") window_label = tk.Label(root, @@ -374,7 +381,8 @@ def save_input(): json.dump(new_inits, initsf, indent=4) print(f"{fname} saved!") - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") window_label = tk.Label(root, @@ -419,7 +427,8 @@ def save_input(): print("################################################\n\n") pass else: - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") fitsortext = tk.IntVar() @@ -441,6 +450,7 @@ def save_input(): padx=20, variable=fitsortext, value=2).pack(anchor=tk.W) + fitsortext.set(2) # Button for closing exit_button = tk.Button(root, text="Next", command=root.destroy) @@ -450,7 +460,8 @@ def save_input(): root.mainloop() if fitsortext.get() == 2: - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") i=0; j=0 @@ -528,7 +539,8 @@ def save_input(): root.mainloop() # First ask user how they want to enter the observing information - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) obsinfo = tk.StringVar() root.title(f"EXOTIC v{__version__}") @@ -560,7 +572,8 @@ def save_input(): if obsinfo.get() == "manual": - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") initparams = tk.IntVar() @@ -755,7 +768,8 @@ def save_input(): obsnotes_entry.grid(row=i, column=j + 1, sticky=tk.W, pady=2) i += 1 - # root = tk.Tk() + # root=tk.Tk() + # root.protocol("WM_DELETE_WINDOW", exit) # root.title(f"EXOTIC v{__version__}") # # window_label = tk.Label(root, @@ -833,7 +847,8 @@ def save_input(): try: if filteroptions.get() == "N/A": - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") i = 0 @@ -876,7 +891,8 @@ def save_input(): pass # then ask user how they want to enter the planetary information - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") planetparams = tk.StringVar() @@ -913,7 +929,8 @@ def save_input(): root.mainloop() if planetparams.get() == "manual": - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") initparams = tk.IntVar() @@ -1169,7 +1186,8 @@ def save_input(): tk.mainloop() elif planetparams.get() == 'nea': - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") initparams = tk.IntVar() @@ -1227,7 +1245,8 @@ def save_input(): tk.mainloop() if (planetparams.get() == "inits") or (obsinfo.get() == "inits"): - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") # # "Directory with FITS files": "sample-data/HatP32Dec202017", @@ -1246,7 +1265,8 @@ def save_inputs(): tk.mainloop() # if (obsinfo.get() == "manual"): - # root = tk.Tk() + # root=tk.Tk() + # root.protocol("WM_DELETE_WINDOW", exit) # root.title(f"EXOTIC v{__version__}") # # window_label = tk.Label(root, @@ -1285,12 +1305,12 @@ def save_inputs(): # # exit_button.pack(pady=20) # root.update() # root.mainloop() - # root.mainloop() # Create an initialization file if it does not already exist if (planetparams.get() != "inits") or (obsinfo.get() != "inits"): - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") initparams = tk.IntVar() @@ -1485,7 +1505,8 @@ def save_input(): json.dump(new_inits, initsf, indent=4) print(f"{fname} saved!") - root = tk.Tk() + root=tk.Tk() + root.protocol("WM_DELETE_WINDOW", exit) root.title(f"EXOTIC v{__version__}") window_label = tk.Label(root, From 6f0822a093dbd662f0d8083bef4be4e8f4664fca Mon Sep 17 00:00:00 2001 From: Kyle A Pearson Date: Wed, 19 Jan 2022 12:51:38 -0800 Subject: [PATCH 15/17] updates to global fitter --- examples/global_fitter_unit_test.py | 16 ++++---- exotic/api/elca.py | 57 +++++++++++++---------------- exotic/exotic.py | 2 + 3 files changed, 35 insertions(+), 40 deletions(-) diff --git a/examples/global_fitter_unit_test.py b/examples/global_fitter_unit_test.py index 10f263bb..26c4d7c3 100644 --- a/examples/global_fitter_unit_test.py +++ b/examples/global_fitter_unit_test.py @@ -8,6 +8,7 @@ # simulate input data epochs = np.random.choice(np.arange(100), 5, replace=False) input_data = [] + local_bounds = [] for i, epoch in enumerate(epochs): @@ -45,6 +46,12 @@ 'priors':prior }) + # individual properties + local_bounds.append({ + 'rprs':[0,0.2], + 'a2':[-0.5,0] + }) + #plt.plot(time,flux,marker='o') #plt.plot(time, model,ls='-') #plt.show() @@ -56,13 +63,6 @@ 'ars':[14,14.5], } - # individual properties - local_bounds = { - 'rprs':[0,0.2], - 'a1':[0,1e4], - 'a2':[-0.5,0] - } - print('epochs:',epochs) myfit = glc_fitter(input_data, global_bounds, local_bounds, individual_fit=True, verbose=True) @@ -70,4 +70,4 @@ plt.show() myfit.plot_bestfit() - plt.show() \ No newline at end of file + plt.show() diff --git a/exotic/api/elca.py b/exotic/api/elca.py index 40c83727..3d7f2e5f 100644 --- a/exotic/api/elca.py +++ b/exotic/api/elca.py @@ -539,23 +539,24 @@ def fit_nested(self): nobs = len(self.data) gfreekeys = list(self.global_bounds.keys()) - if isinstance(self.local_bounds, dict): - lfreekeys = list(self.local_bounds.keys()) - boundarray = np.vstack([ [self.global_bounds[k] for k in gfreekeys], [self.local_bounds[k] for k in lfreekeys]*nobs ]) - else: - # if list type - lfreekeys = list(self.local_bounds[0].keys()) - boundarray = [self.global_bounds[k] for k in gfreekeys] - for i in range(nobs): - boundarray.extend([self.local_bounds[i][k] for k in lfreekeys]) - boundarray = np.array(boundarray) + # if isinstance(self.local_bounds, dict): + # lfreekeys = list(self.local_bounds.keys()) + # boundarray = np.vstack([ [self.global_bounds[k] for k in gfreekeys], [self.local_bounds[k] for k in lfreekeys]*nobs ]) + # else: + # # if list type + lfreekeys = [] + boundarray = [self.global_bounds[k] for k in gfreekeys] + for i in range(nobs): + lfreekeys.append(list(self.local_bounds[i].keys())) + boundarray.extend([self.local_bounds[i][k] for k in lfreekeys[-1]]) + boundarray = np.array(boundarray) # fit individual light curves to constrain priors if self.individual_fit: for i in range(nobs): print(f"Fitting individual light curve {i+1}/{nobs}") - mybounds = dict(**self.local_bounds, **self.global_bounds) + mybounds = dict(**self.local_bounds[i], **self.global_bounds) if 'per' in mybounds: del(mybounds['per']) myfit = lc_fitter( @@ -568,14 +569,14 @@ def fit_nested(self): ) self.data[i]['individual'] = myfit.parameters.copy() - + ti = sum([len(self.local_bounds[k]) for k in range(i)]) # update local priors - for j, key in enumerate(self.local_bounds.keys()): + for j, key in enumerate(self.local_bounds[i].keys()): - boundarray[j+i*len(self.local_bounds)+len(gfreekeys),0] = myfit.parameters[key] - 5*myfit.errors[key] - boundarray[j+i*len(self.local_bounds)+len(gfreekeys),1] = myfit.parameters[key] + 5*myfit.errors[key] + boundarray[j+ti+len(gfreekeys),0] = myfit.parameters[key] - 5*myfit.errors[key] + boundarray[j+ti+len(gfreekeys),1] = myfit.parameters[key] + 5*myfit.errors[key] if key == 'rprs': - boundarray[j+i*len(self.local_bounds)+len(gfreekeys),0] = max(0,myfit.parameters[key] - 5*myfit.errors[key]) + boundarray[j+ti+len(gfreekeys),0] = max(0,myfit.parameters[key] - 5*myfit.errors[key]) del(myfit) @@ -595,8 +596,9 @@ def loglike(pars): self.data[i]['priors'][key] = pars[j] # local keys - for j, key in enumerate(lfreekeys): - self.data[i]['priors'][key] = pars[j+i*len(lfreekeys)+len(gfreekeys)] + ti = sum([len(self.local_bounds[k]) for k in range(i)]) + for j, key in enumerate(lfreekeys[i]): + self.data[i]['priors'][key] = pars[j+ti+len(gfreekeys)] # compute model model = transit(self.data[i]['time'], self.data[i]['priors']) @@ -611,7 +613,7 @@ def loglike(pars): freekeys = []+gfreekeys for n in range(nobs): - for k in lfreekeys: + for k in lfreekeys[n]: freekeys.append(f"local_{n}_{k}") if self.verbose: @@ -632,7 +634,7 @@ def loglike(pars): self.results['posterior']['errup'][i]] for n in range(nobs): - for k in lfreekeys: + for k in lfreekeys[n]: pkey = f"local_{n}_{k}" self.data[n]['priors'][k] = self.parameters[pkey] @@ -703,20 +705,11 @@ def plot_bestfit(self): } # example generating LD coefficients - from exotic.exotic import LimbDarkening - from ldtk.filters import create_tess - - tessfilter = create_tess() - - ld_obj = LimbDarkening( - teff=prior['teff'], teffpos=prior['tefferr'], teffneg=prior['tefferr'], - met=prior['met'], metpos=prior['meterr'], metneg=prior['meterr'], - logg=prior['logg'], loggpos=prior['loggerr'], loggneg=prior['loggerr'], - wl_min=tessfilter.wl.min(), wl_max=tessfilter.wl.max(), filter_type="Clear") + from pylightcurve import exotethys - ld0, ld1, ld2, ld3, filt, wlmin, wlmax = ld_obj.nonlinear_ld() + u0,u1,u2,u3 = exotethys(prior['logg'], prior['teff'], prior['met'], 'TESS', method='claret', stellar_model='phoenix') - prior['u0'],prior['u1'],prior['u2'],prior['u3'] = [ld0[0], ld1[0], ld2[0], ld3[0]] + prior['u0'],prior['u1'],prior['u2'],prior['u3'] = u0,u1,u2,u3 time = np.linspace(0.7, 0.8, 1000) # [day] diff --git a/exotic/exotic.py b/exotic/exotic.py index 74595575..f25a9472 100644 --- a/exotic/exotic.py +++ b/exotic/exotic.py @@ -1919,6 +1919,8 @@ def main(): log_info(f"Optimal Annulus: {np.round(minAnnulus, 2)}") bestaperture = str(abs(np.round(minAperture, 2))) log_info("*********************************************\n") + + import pdb; pdb.set_trace() # Take the BJD times from the image headers if "BJD_TDB" in image_header or "BJD" in image_header or "BJD_TBD" in image_header: From 0787266bfc07f6b9e85e80a4c858f1d57bccc418 Mon Sep 17 00:00:00 2001 From: Kyle A Pearson Date: Wed, 26 Jan 2022 10:58:56 -0800 Subject: [PATCH 16/17] global fitter fixes --- examples/global_fitter_unit_test.py | 12 ++- exotic/api/elca.py | 110 ++++++++++++++++++++++++---- exotic/exotic.py | 8 +- 3 files changed, 112 insertions(+), 18 deletions(-) diff --git a/examples/global_fitter_unit_test.py b/examples/global_fitter_unit_test.py index 26c4d7c3..b04cdfe9 100644 --- a/examples/global_fitter_unit_test.py +++ b/examples/global_fitter_unit_test.py @@ -6,7 +6,7 @@ if __name__ == "__main__": # simulate input data - epochs = np.random.choice(np.arange(100), 5, replace=False) + epochs = np.random.choice(np.arange(100), 3, replace=False) input_data = [] local_bounds = [] @@ -64,10 +64,16 @@ } print('epochs:',epochs) - myfit = glc_fitter(input_data, global_bounds, local_bounds, individual_fit=True, verbose=True) + myfit = glc_fitter(input_data, global_bounds, local_bounds, individual_fit=False, verbose=True) + + myfit.plot_bestfit() + plt.show() myfit.plot_triangle() plt.show() - myfit.plot_bestfit() + myfit.plot_bestfits() plt.show() + + + diff --git a/exotic/api/elca.py b/exotic/api/elca.py index 3d7f2e5f..728d0dbc 100644 --- a/exotic/api/elca.py +++ b/exotic/api/elca.py @@ -621,12 +621,11 @@ def loglike(pars): else: self.results = ReactiveNestedSampler(freekeys, loglike, prior_transform).run(max_ncalls=4e5, show_status=self.verbose, viz_callback=self.verbose) - self.errors = {} self.quantiles = {} - self.parameters = {} + self.errors = {} + self.parameters = self.data[0]['priors'].copy() for i, key in enumerate(freekeys): - self.parameters[key] = self.results['maximum_likelihood']['point'][i] self.errors[key] = self.results['posterior']['stdev'][i] self.quantiles[key] = [ @@ -634,20 +633,26 @@ def loglike(pars): self.results['posterior']['errup'][i]] for n in range(nobs): + self.data[n]['errors'] = {} for k in lfreekeys[n]: pkey = f"local_{n}_{k}" self.data[n]['priors'][k] = self.parameters[pkey] + self.data[n]['errors'][k] = self.errors[pkey] + + if k == 'rprs' and 'rprs' not in freekeys: + self.parameters[k] = self.data[n]['priors'][k] + self.errors[k] = self.data[n]['errors'][k] # solve for a1 model = transit(self.data[n]['time'], self.data[n]['priors']) airmass = np.exp(self.data[n]['airmass']*self.data[n]['priors']['a2']) detrend = self.data[n]['flux']/(model*airmass) self.data[n]['priors']['a1'] = np.median(detrend) + self.data[n]['residuals'] = self.data[n]['flux'] - model*airmass*self.data[n]['priors']['a1'] + self.data[n]['detrend'] = self.data[n]['flux']/(airmass*self.data[n]['priors']['a1']) - def plot_bestfit(self): - # TODO mosaic - nobs = len(self.data) - nrows = nobs//4+1 + def plot_bestfits(self): + nrows = len(self.data)//4+1 fig,ax = plt.subplots(nrows, 4, figsize=(5+5*nrows,4*nrows)) # turn off all axes @@ -659,8 +664,11 @@ def plot_bestfit(self): else: ax[ri,ci].axis('off') + markers = ['o','v','^','<','>','s','p','*','h','H','D','d','P','X'] + colors = ['black','blue','green','orange','purple','brown','pink','grey','magenta','cyan','yellow','lime'] + # plot observations - for i in range(nobs): + for i in range(len(self.data)): ri = int(i/4) ci = i%4 model = transit(self.data[i]['time'], self.data[i]['priors']) @@ -669,17 +677,93 @@ def plot_bestfit(self): if ax.ndim == 1: ax[i].axis('on') - ax[i].errorbar(self.data[i]['time'], self.data[i]['flux']/airmass/detrend.mean(), yerr=self.data[i]['ferr']/airmass/detrend.mean(), ls='none', marker='o', color='black', alpha=0.5) - ax[i].plot(self.data[i]['time'], model, 'r-') + ax[i].errorbar(self.data[i]['time'], self.data[i]['flux']/airmass/detrend.mean(), yerr=self.data[i]['ferr']/airmass/detrend.mean(), + ls='none', marker=markers[i], color=colors[i], alpha=0.5) + ax[i].plot(self.data[i]['time'], model, 'r-', zorder=2) ax[i].set_xlabel("Time") else: ax[ri,ci].axis('on') - ax[ri,ci].errorbar(self.data[i]['time'], self.data[i]['flux']/airmass/detrend.mean(), yerr=self.data[i]['ferr']/airmass/detrend.mean(), ls='none', marker='o', color='black', alpha=0.5) - ax[ri,ci].plot(self.data[i]['time'], model, 'r-') + ax[ri,ci].errorbar(self.data[i]['time'], self.data[i]['flux']/airmass/detrend.mean(), yerr=self.data[i]['ferr']/airmass/detrend.mean(), + ls='none', marker=markers[i], color=colors[i], alpha=0.5) + ax[ri,ci].plot(self.data[i]['time'], model, 'r-', zorder=2) ax[ri,ci].set_xlabel("Time") plt.tight_layout() - return fig + return fig + + def plot_bestfit(self, title="", bin_dt=30./(60*24)): + f = plt.figure(figsize=(9,6)) + f.subplots_adjust(top=0.92,bottom=0.09,left=0.14,right=0.98, hspace=0) + ax_lc = plt.subplot2grid((4,5), (0,0), colspan=5,rowspan=3) + ax_res = plt.subplot2grid((4,5), (3,0), colspan=5, rowspan=1) + axs = [ax_lc, ax_res] + + axs[0].set_title(title) + axs[0].set_ylabel("Relative Flux", fontsize=14) + axs[0].grid(True,ls='--') + + rprs2 = self.parameters['rprs']**2 + rprs2err = 2*self.parameters['rprs']*self.errors['rprs'] + lclabel1 = r"$R^{2}_{p}/R^{2}_{s}$ = %s $\pm$ %s" %( + str(round_to_2(rprs2, rprs2err)), + str(round_to_2(rprs2err)) + ) + + lclabel2 = r"$T_{mid}$ = %s $\pm$ %s BJD$_{TDB}$" %( + str(round_to_2(self.parameters['tmid'], self.errors.get('tmid',0))), + str(round_to_2(self.errors.get('tmid',0))) + ) + + lclabel = lclabel1 + "\n" + lclabel2 + minp = 1 + maxp = 0 + + min_std = 1 + markers = ['o','v','^','<','>','s','p','*','h','H','D','d','P','X'] + colors = ['black','blue','green','orange','purple','brown','pink','grey','magenta','cyan','yellow','lime'] + for n in range(len(self.data)): + phase = get_phase(self.data[n]['time'], self.parameters['per'], self.data[n]['priors']['tmid']) + si = np.argsort(phase) + bt2, br2, _ = time_bin(phase[si]*self.parameters['per'], self.data[n]['residuals'][si]/np.median(self.data[n]['flux'])*1e2, bin_dt) + + # plot data + axs[0].errorbar(phase, self.data[n]['detrend'], yerr=np.std(self.data[n]['residuals'])/np.median(self.data[n]['flux']), + ls='none', marker=markers[n], color=colors[n], zorder=1, alpha=0.2) + + # plot residuals + axs[1].plot(phase, self.data[n]['residuals']/np.median(self.data[n]['flux'])*1e2, color=colors[n], marker=markers[n], ls='none', + alpha=0.2, label=r'$\sigma$ = {:.2f} %'.format( np.std(self.data[n]['residuals']/np.median(self.data[n]['flux'])*1e2))) + + # plot binned data + bt2, bf2, bs = time_bin(phase[si]*self.data[n]['priors']['per'], self.data[n]['detrend'][si], bin_dt) + axs[0].errorbar(bt2/self.data[n]['priors']['per'],bf2,yerr=bs,alpha=1,zorder=2,color=colors[n],ls='none',marker=markers[n]) + + # replace min and max for upsampled lc model + minp = min(minp, min(phase)) + maxp = max(maxp, max(phase)) + min_std = min(min_std, np.std(self.data[n]['residuals']/np.median(self.data[n]['flux']))) + + # best fit model + self.time_upsample = np.linspace(minp*self.parameters['per']+self.parameters['tmid'], + maxp*self.parameters['per']+self.parameters['tmid'], 10000) + self.transit_upsample = transit(self.time_upsample, self.parameters) + self.phase_upsample = get_phase(self.time_upsample, self.parameters['per'], self.parameters['tmid']) + sii = np.argsort(self.phase_upsample) + axs[0].plot(self.phase_upsample[sii], self.transit_upsample[sii], 'r-', zorder=3, label=lclabel) + + axs[0].set_xlim([min(self.phase_upsample), max(self.phase_upsample)]) + axs[0].set_xlabel("Phase ", fontsize=14) + axs[0].set_ylim([1-self.parameters['rprs']**2-5*min_std, 1+5*min_std]) + axs[1].set_xlim([min(self.phase_upsample), max(self.phase_upsample)]) + axs[1].set_xlabel("Phase", fontsize=14) + axs[1].set_ylim([-5*min_std*1e2, 5*min_std*1e2]) + + axs[0].get_xaxis().set_visible(False) + axs[1].legend(loc='best') + axs[0].legend(loc='best') + axs[1].set_ylabel("Residuals [%]", fontsize=14) + axs[1].grid(True,ls='--',axis='y') + return f,axs if __name__ == "__main__": diff --git a/exotic/exotic.py b/exotic/exotic.py index f25a9472..343de3aa 100644 --- a/exotic/exotic.py +++ b/exotic/exotic.py @@ -1919,8 +1919,6 @@ def main(): log_info(f"Optimal Annulus: {np.round(minAnnulus, 2)}") bestaperture = str(abs(np.round(minAperture, 2))) log_info("*********************************************\n") - - import pdb; pdb.set_trace() # Take the BJD times from the image headers if "BJD_TDB" in image_header or "BJD" in image_header or "BJD_TBD" in image_header: @@ -2067,6 +2065,12 @@ def main(): upper = pDict['midT'] + 35 * pDict['midTUnc'] + np.floor(phase).max() * (pDict['pPer'] + 35 * pDict['pPerUnc']) lower = pDict['midT'] - 35 * pDict['midTUnc'] + np.floor(phase).max() * (pDict['pPer'] - 35 * pDict['pPerUnc']) + # clip bounds so they're within 1 orbit + if upper > prior['tmid'] + 0.5*prior['per']: + upper = prior['tmid'] + 0.5*prior['per'] + if lower < prior['tmid'] - 0.5*prior['per']: + lower = prior['tmid'] - 0.5*prior['per'] + if np.floor(phase).max() - np.floor(phase).min() == 0: log_info("Error: Estimated mid-transit not in observation range (check priors or observation time)", error=True) log_info(f"start:{np.min(goodTimes)}", error=True) From d5e67a99088562b5643fcfda9c969ad5f842f19d Mon Sep 17 00:00:00 2001 From: Rob Zellem Date: Wed, 26 Jan 2022 13:06:18 -0800 Subject: [PATCH 17/17] v1.10.0: updates to priors --- exotic/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exotic/version.py b/exotic/version.py index 38cf6dbe..fcfdf383 100644 --- a/exotic/version.py +++ b/exotic/version.py @@ -1 +1 @@ -__version__ = "1.9.1" +__version__ = "1.10.0"