diff --git a/examples/global_fitter_unit_test.py b/examples/global_fitter_unit_test.py new file mode 100644 index 00000000..b04cdfe9 --- /dev/null +++ b/examples/global_fitter_unit_test.py @@ -0,0 +1,79 @@ +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), 3, replace=False) + input_data = [] + local_bounds = [] + + 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 + }) + + # 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() + + # 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], + } + + print('epochs:',epochs) + 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_bestfits() + plt.show() + + + 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 diff --git a/exotic/api/elca.py b/exotic/api/elca.py index c0bbab23..c6685742 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 @@ -518,6 +518,254 @@ 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 = [] + 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[i], **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() + ti = sum([len(self.local_bounds[k]) for k in range(i)]) + # update local priors + for j, key in enumerate(self.local_bounds[i].keys()): + + 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+ti+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 + 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']) + 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[n]: + 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.quantiles = {} + 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] = [ + self.results['posterior']['errlo'][i], + 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_bestfits(self): + nrows = len(self.data)//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') + + 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(len(self.data)): + 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=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=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 + + 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__": prior = { @@ -541,20 +789,11 @@ def plot_triangle(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 785d534a..346adb6c 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 @@ -90,7 +90,7 @@ 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.util import view_as_windows from skimage.transform import SimilarityTransform # error handling for scraper from tenacity import retry, stop_after_delay @@ -777,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]: @@ -819,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 (arcsec/pixel): ", type_=float) image_scale = f"Image scale in arcsecs/pixel: {image_scale_num}" return image_scale @@ -891,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) @@ -915,12 +930,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 +955,11 @@ 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 + return res.x @@ -1188,7 +1212,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 +1232,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 +1278,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 +1954,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 +2027,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 @@ -2021,6 +2064,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) @@ -2034,6 +2083,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) 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, 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" diff --git a/requirements.txt b/requirements.txt index 948e030a..9aef81b8 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,19 +7,19 @@ 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 lmfit~=1.0 matplotlib>=3.4 -numpy~=1.20 +numpy~=1.21 pandas~=1.3 panel~=0.12 photutils>=0.7 python_dateutil~=2.8 -pyvo~=1.1 -requests~=2.26 +pyvo~=1.2 +requests~=2.27 scipy~=1.7 scikit-image~=0.18 tenacity~=8.0