From 1e01c6c5895b45d57b6eb4ae5b87d885c98225b0 Mon Sep 17 00:00:00 2001 From: hposborn Date: Wed, 3 Jul 2024 11:07:44 +0200 Subject: [PATCH] readding a bunch of dropped files --- MonoTools/data/tables/tess_lc_locations.csv | 140 +++--- MonoTools/fit.py | 452 ++++++++++++-------- MonoTools/lightcurve.py | 1 - __init__.py | 0 pyproject.toml | 0 setup.py | 2 +- 6 files changed, 357 insertions(+), 238 deletions(-) mode change 100644 => 100755 __init__.py mode change 100644 => 100755 pyproject.toml diff --git a/MonoTools/data/tables/tess_lc_locations.csv b/MonoTools/data/tables/tess_lc_locations.csv index 72ac4c3..344d303 100755 --- a/MonoTools/data/tables/tess_lc_locations.csv +++ b/MonoTools/data/tables/tess_lc_locations.csv @@ -1,64 +1,78 @@ ,date,runid -1,2018206045859,120 -2,2018234235059,121 -3,2018263035959,123 -4,2018292075959,124 -5,2018319095959,125 -6,2018349182500,126 -7,2019006130736,131 -8,2019032160000,136 -9,2019058134432,139 -10,2019085135100,140 -11,2019112060037,143 -12,2019140104343,144 -13,2019169103026,146 -14,2019198215352,150 -15,2019226182529,151 -16,2019253231442,152 -17,2019279210107,161 -18,2019306063752,162 -19,2019331140908,164 -20,2019357164649,165 -21,2020020091053,167 -22,2020049080258,174 -23,2020078014623,177 -24,2020106103520,180 -25,2020133194932,182 -26,2020160202036,188 -27,2020186164531,189 -28,2020212050318,190 -29,2020238165205,193 -30,2020266004630,195 -31,2020294194027,198 -32,2020324010417,200 -33,2020351194500,203 -34,2021014023720,204 -35,2021039152502,205 -36,2021065132309,207 -37,2021091135823,208 -38,2021118034608,209 -39,2021146024351,210 -40,2021175071901,211 -41,2021204101404,212 -42,2021232031932,213 -43,2021258175143,214 -44,2021284114741,215 -45,2021310001228,216 -46,2021336043614,217 -47,2021364111932,218 -48,2022027120115,219 -49,2022057073128,221 -50,2022085151738,222 -51,2022112184951,223 -52,2022138205153,224 -53,2022164095748,226 -54,2022190063128,227 -55,2022217014003,242 -56,2022244194134,243 -57,2022273165103,245 -58,2022302161335,247 -59,2022330142927,248 -60,2022357055054,249 -61,2023018032328,250 -62,2023043185947,254 -63,2023069172124,255 +1,2018206045859.0,120.0 +2,2018234235059.0,121.0 +3,2018263035959.0,123.0 +4,2018292075959.0,124.0 +5,2018319095959.0,125.0 +6,2018349182500.0,126.0 +7,2019006130736.0,131.0 +8,2019032160000.0,136.0 +9,2019058134432.0,139.0 +10,2019085135100.0,140.0 +11,2019112060037.0,143.0 +12,2019140104343.0,144.0 +13,2019169103026.0,146.0 +14,2019198215352.0,150.0 +15,2019226182529.0,151.0 +16,2019253231442.0,152.0 +17,2019279210107.0,161.0 +18,2019306063752.0,162.0 +19,2019331140908.0,164.0 +20,2019357164649.0,165.0 +21,2020020091053.0,167.0 +22,2020049080258.0,174.0 +23,2020078014623.0,177.0 +24,2020106103520.0,180.0 +25,2020133194932.0,182.0 +26,2020160202036.0,188.0 +27,2020186164531.0,189.0 +28,2020212050318.0,190.0 +29,2020238165205.0,193.0 +30,2020266004630.0,195.0 +31,2020294194027.0,198.0 +32,2020324010417.0,200.0 +33,2020351194500.0,203.0 +34,2021014023720.0,204.0 +35,2021039152502.0,205.0 +36,2021065132309.0,207.0 +37,2021091135823.0,208.0 +38,2021118034608.0,209.0 +39,2021146024351.0,210.0 +40,2021175071901.0,211.0 +41,2021204101404.0,212.0 +42,2021232031932.0,213.0 +43,2021258175143.0,214.0 +44,2021284114741.0,215.0 +45,2021310001228.0,216.0 +46,2021336043614.0,217.0 +47,2021364111932.0,218.0 +48,2022027120115.0,219.0 +49,2022057073128.0,221.0 +50,2022085151738.0,222.0 +51,2022112184951.0,223.0 +52,2022138205153.0,224.0 +53,2022164095748.0,226.0 +54,2022190063128.0,227.0 +55,2022217014003.0,242.0 +56,2022244194134.0,243.0 +57,2022273165103.0,245.0 +58,2022302161335.0,247.0 +59,2022330142927.0,248.0 +60,2022357055054.0,249.0 +61,2023018032328.0,250.0 +62,2023043185947.0,254.0 +63,2023069172124.0,255.0 +64,2023096110322.0,257.0 +65,2023124020739.0,259.0 +66,2023153011303.0,260.0 +67,2023181235917.0,261.0 +68,2023209231226.0,262.0 +69,2023237165326.0,264.0 +70,2023263165758.0,265.0 +71,2023289093419.0,266.0 +72,2023315124025.0,267.0 +73,2023341045131.0,268.0 +74,2024003055635.0,269.0 +75,2024030031500.0,270.0 +76,2024058030222.0,271.0 +77,2024085201119.0,272.0 diff --git a/MonoTools/fit.py b/MonoTools/fit.py index e415877..2395d97 100755 --- a/MonoTools/fit.py +++ b/MonoTools/fit.py @@ -139,6 +139,8 @@ def __init__(self, ID, mission, lc=None, rvs=None, planets=None, overwrite=False 'pred_all':False, # Do we predict all time array, or only a cut-down version? 'use_multinest':False, # use_multinest - bool - currently not supported 'use_pymc3':True, # use_pymc3 - bool + 'bin_all':False, # bin_all - bool - Bin all points to 10mins (to speed up certain) + 'bin_all_size':10/1440., # bin_all_size - float - Bin size if binning all points in minutes (default to 10mins) 'bin_oot':True, # bin_oot - bool - Bin points outside the cut_distance to 30mins 'model_t03_ttv':False, # model_t03_ttv - bool - Whether to model the third transit as a seperate parameter, otherwise it is constrained with the other params 'timing_sigma':0.1, # timing_sigma - float - the sigma (as a function of transit duration) to use when setting transit times. Default=0.1*t_D @@ -350,7 +352,7 @@ def add_rvplanet(self, pl_dic, name,**kwargs): self.rvplanets[name]=pl_dic - def add_multi(self, pl_dic, name,**kwargs): + def add_multi(self, pl_dic, name, update_per=False, **kwargs): """Adds a transiting planet with multiple, consecutive transits to planet properties dict Args: @@ -378,6 +380,14 @@ def add_multi(self, pl_dic, name,**kwargs): pl_dic['b']=np.clip((1+pl_dic['ror'])**2 - (pl_dic['tdur']*86400)**2 * \ ((3*pl_dic['period']*86400) / (np.pi**2*6.67e-11*rho_S*1410))**(-2/3), 0.01,2.0)**0.5 + if update_per: + if not hasattr(self.lc,'flux_flat'): + try: + self.lc.flatten() + except: + setattr(self.lc,'flux_flat',self.lc.flux[:]) + pl_dic['period']=tools.update_period_w_tls(self.lc.time[self.lc.mask], self.lc.flux_flat[self.lc.mask], pl_dic['period']) + phase=(self.lc.time-pl_dic['tcen']-0.5*pl_dic['period'])%pl_dic['period']-0.5*pl_dic['period'] self.lc.near_trans[name] = abs(phase)match_trans_thresh*tdur and abs((tcen_3-tcen-per*0.5)%per-per*0.5)>match_trans_thresh*tdur: # and (np.sum(trans2&intr)<0.75*days_in_known_transits[1] or np.sum(trans3&intr)<0.75*days_in_known_transits[2]): #Either second or third transit does not match with this period... Adding zero to list. - check_pers_ix+=[False] + #check_pers_ix+=[False]# + days_in_tr=np.sum([float(self.lc.cadence[ncad].split('_')[1])/86400 for ncad in np.arange(len(self.lc.cadence))[self.lc.mask][intr]]) + check_pers_ix+=[days_in_tr<(1.0+coverage_thresh)*np.sum(days_in_known_transits)] #print(per,"NO",abs((tcen_2-tcen-per*0.5)%per-per*0.5),match_trans_thresh*tdur,abs((tcen_3-tcen-per*0.5)%per-per*0.5),match_trans_thresh*tdur) else: #Here we need to add up the cadences in transit (and not simply count the points) to check coverage: @@ -653,13 +676,13 @@ def compute_period_aliases(self,pl_dic,dur=0.5,**kwargs): #Also need to check that the implied periods match the third period check_pers_ix = self.CheckPeriodsHaveGaps(pl_dic['period']/check_pers_ints,pl_dic['tdur'],pl_dic['tcen'],tcen_2=pl_dic['tcen_2'],tcen_3=pl_dic['tcen_3'],**kwargs) else: - check_pers_ix = self.CheckPeriodsHaveGaps(pl_dic['period']/check_pers_ints,pl_dic['tdur'],pl_dic['tcen'],tcen_2=pl_dic['tcen_2'],**kwargs) pl_dic['period_int_aliases']=check_pers_ints[check_pers_ix] - if pl_dic['period_int_aliases']==[]: + if len(pl_dic['period_int_aliases'])==0: print("problem in computing Duotransit aliases") else: + print(pl_dic) pl_dic['period_aliases']=pl_dic['period']/pl_dic['period_int_aliases'] pl_dic['P_min']=np.min(pl_dic['period_aliases']) return pl_dic @@ -763,6 +786,7 @@ def add_trio(self, pl_dic,name,maxint=24,**kwargs): if 'period_err' not in pl_dic or not np.isfinite(pl_dic['period_err']): pl_dic['period_err'] = 0.1666*pl_dic['tdur'] #Calculating P_min and the integer steps + print(pl_dic) pl_dic=self.compute_period_aliases(pl_dic,**kwargs) pl_dic['npers']=len(pl_dic['period_int_aliases']) @@ -885,6 +909,7 @@ def init_starpars(self,Rstar=None,Teff=None,logg=None,FeH=0.0,rhostar=None,Mstar rho_MR=[Mstar[0]/self.Rstar[0]**3] rho_MR+=[(Mstar[0]+Mstar[1])/(self.Rstar[0]-abs(self.Rstar[1]))**3/rho_MR[0]-1.0, 1.0-(Mstar[0]-abs(Mstar[2]))/(self.Rstar[0]+self.Rstar[2])**3/rho_MR[0]] + print(rho_MR) #Weighted sums of two avenues to density: rhostar=[rho_logg[0]*(rho_MR[1]+rho_MR[2])/(rho_logg[1]+rho_logg[2]+rho_MR[1]+rho_MR[2])+ rho_MR[0]*(rho_logg[1]+rho_logg[2])/(rho_logg[1]+rho_logg[2]+rho_MR[1]+rho_MR[2])] @@ -1030,11 +1055,36 @@ def init_lc(self, oot_binsize=1/48, **kwargs): """Initialise light curve. This can be done either after or before model initialisation. This function creates transit maskes and flattens/bins the light curve in ways to avoid influencing the transit. """ - # + + #Replacing all timeseries with binned timeseries: + if self.bin_all: + print("BIN ALL") + print(np.min(np.diff(self.lc.time))) + if hasattr(self,'in_trans'): + delattr(self, 'in_trans') + self.lc.remove_binned_arrs() + print([ts for ts in self.lc.timeseries if type(getattr(self.lc,ts))==np.ndarray and type(getattr(self.lc,ts)[0]) in (float,np.float64)]) + print(type(self.lc.bg_flux)) + self.lc.bin(timeseries=[ts for ts in self.lc.timeseries if type(getattr(self.lc,ts))==np.ndarray and type(getattr(self.lc,ts)[0]) in (float,np.float32,np.float64)], + use_masked=False, binsize=self.bin_all_size) + + for key in self.lc.timeseries: + if 'bin' not in key and "bin_"+key in self.lc.timeseries: + setattr(self.lc, key, getattr(self.lc, "bin_"+key)) + print(key) + elif 'bin' not in key: + print(key,"not binned? Removing.") + if 'mask' in key: + setattr(self.lc,key,np.tile(True,len(self.lc.bin_time))) + self.lc.remove_binned_arrs() + self.lc.make_mask(overwrite=True) + print("BIN ALL") + print(np.min(np.diff(self.lc.time))) + #step=0.133*np.min([self.init_soln['tdur_'+pl] for pl in self.planets]) if hasattr(self,'init_soln') else 0.133*np.min([self.planets[pl]['tdur'] for pl in self.planets]) #win=6.5*np.max([self.init_soln['tdur_'+pl] for pl in self.planets]) if hasattr(self,'init_soln') else 6.5*np.max([self.planets[pl]['tdur'] for pl in self.planets]) - + self.lc.near_trans = {'all':np.tile(False, len(self.lc.time))} self.lc.in_trans = {'all':np.tile(False, len(self.lc.time))} for pl in self.planets: @@ -1043,9 +1093,13 @@ def init_lc(self, oot_binsize=1/48, **kwargs): p = self.init_soln['per_'+pl] if hasattr(self,'init_soln') else self.planets[pl]['period'] phase=(self.lc.time-t0-0.5*p)%p-0.5*p elif pl in self.trios: - t0= self.init_soln['t0_'+pl] if hasattr(self,'init_soln') else self.planets[pl]['tcen_2'] p=np.max(self.init_soln['per_'+pl]) if hasattr(self,'init_soln') else np.max(self.planets[pl]['period_aliases']) - phase=(self.lc.time-t0-0.5*p)%p-0.5*p + if self.model_t03_ttv: + t0s= [self.init_soln['t0_'+pl],self.init_soln['t0_2_'+pl],self.init_soln['t0_3_'+pl]] if hasattr(self,'init_soln') else [self.planets[pl]['tcen'],self.planets[pl]['tcen_2'],self.planets[pl]['tcen_3']] + + else: + t0= self.init_soln['t0_'+pl] if hasattr(self,'init_soln') else self.planets[pl]['tcen_2'] + phase=self.make_phase(self.lc.time,t0s,p) elif pl in self.duos: t0= self.init_soln['t0_'+pl] if hasattr(self,'init_soln') else self.planets[pl]['tcen'] p=abs(self.init_soln['t0_2_'+pl]-self.init_soln['t0_'+pl]) if hasattr(self,'init_soln') else abs(self.planets[pl]['tcen_2']-self.planets[pl]['tcen']) @@ -1059,13 +1113,14 @@ def init_lc(self, oot_binsize=1/48, **kwargs): self.lc.in_trans[pl] = abs(phase)0 or not self.use_GP: if self.bin_oot: #Creating a pseudo-binned dataset where out-of-transit LC is binned to 30mins but near-transit is not. oot_binsize=1/12 if self.mission.lower()=='kepler' else oot_binsize + print(self.lc.near_trans['all']) self.lc.OOTbin(near_transit_mask=self.lc.near_trans['all'],use_flat=(not self.use_GP or self.fit_no_flatten), binsize=oot_binsize) self.model_time=self.lc.ootbin_time[:] @@ -1176,10 +1231,7 @@ def init_model(self, overwrite=False, **kwargs): self.model_tele_index[:,3:])) if self.use_GP: - self.gp={} - if self.train_GP and not hasattr(self,'gp_init_trace'): - self.GP_training() - + self.init_GP() ###################################### # Initialising sampling models: ###################################### @@ -1189,7 +1241,7 @@ def init_model(self, overwrite=False, **kwargs): elif self.use_multinest: self.run_multinest(**kwargs) - def GP_training(self,n_draws=900,max_len_lc=25000,use_binned=False): + def init_GP(self, n_draws=900, max_len_lc=25000, use_binned=False): """Function to train GPs on out-of-transit photometry Args: @@ -1197,10 +1249,14 @@ def GP_training(self,n_draws=900,max_len_lc=25000,use_binned=False): max_len_lc (int, optional): Maximum length of lightcurve to use (limiting to a few 1000 helps with compute time). Defaults to 25000. uselc (bool, optional): Specify lightcurve to use. Defaults to None, which takes the `mod.lc` light curve. """ - print("initialising and training the GP") + + self.gp={} + + print("initialising the GP") prefix='bin_' if use_binned else '' mask=(~self.lc.in_trans['all'])&self.lc.mask + if len(self.lc.time[(~self.lc.in_trans['all'])&self.lc.mask])>max_len_lc: mask[mask]*=np.arange(0,np.sum(mask),1)3: #Including mutual inclination prior (for high-order multi systems only) av_incl = pm.Deterministic("av_incl",tt.mean([incls[pl] for pl in self.planets])) sd_incl = pm.Deterministic("sd_incl",tt.std([incls[pl] for pl in self.planets])) mut_incl_prior = pm.Potential("mut_incl_prior", tt.exp(-0.5* ((sd_incl - self.mutual_incl_sigma)/self.mutual_incl_sigma)**2)) + phot_mean=pm.Normal("phot_mean",mu=np.median(self.model_flux), sd=2*np.std(self.model_flux)) elif self.local_spline: self.spline_params={} from patsy import dmatrix @@ -1950,6 +2014,7 @@ def gen_lc(i_orbit, i_rpl, n_pl, mask=None,prefix='',make_deterministic=False,pr #print(self.lc['tele_index'][mask,0].astype(bool),len(self.lc['tele_index'][mask,0]),cadmask[mask],len(cadmask[mask])) miss=cad.lower().split('_')[0] cad_index+=[cadmask] + #Have three transits and ttvs - need to modify the input time vector such that if miss=='ts': #Taking the "telescope" index, and adding those points with the matching cadences to the cadmask trans_pred+=[xo.LimbDarkLightCurve(u_star_tess).get_light_curve( @@ -2499,9 +2564,13 @@ def create_orbit(pl, Rs, rho_S, pers, t0s, bs, n_marg=1, eccs=None, omegas=None) initvars3+=[rv_logs2,Ms] if self.use_GP: - initvars3+=[logs2, phot_power, phot_w0, phot_mean] - if self.periodic_kernel is not None: - initvars3+=[periodic_power, periodic_w0, periodic_logQ] + if hasattr(self,'rotation_kernel') and self.rotation_kernel is not None: + initvars3+=[rotation_logamp, rotation_period, rotation_logdeltaQ, rotation_logQ0, rotation_mix] + else: + initvars3+=[logs2, phot_power, phot_w0, phot_mean] + if hasattr(self,'periodic_kernel') and self.periodic_kernel is not None: + initvars3+=[periodic_power, periodic_w0, periodic_logQ] + else: if self.local_spline: initvars3+=[self.spline_params['splines_'+pl+'_'+str(int(n))] for pl in self.planets for n in range(3) if 'splines_'+pl+'_'+str(int(n)) in self.spline_params] @@ -2657,7 +2726,10 @@ def init_gp_to_plot(self, n_samp=7, max_gp_len=12000, interp=True, newgp=False, nchunks=int(np.ceil(2+np.log10(len(self.lc.time)))) timechunks=np.percentile(np.hstack((np.min(self.lc.time)-0.25,self.lc.time,np.max(self.lc.time)-0.25)), np.linspace(0,100,nchunks+1)) #Splitting using a percentile, this way every time point is in the self.lc.time array - min_dist_to_lc=np.hstack([np.min(np.abs(self.lc.time[(self.lc.time>timechunks[tc])&(self.lc.time<=timechunks[tc+1]),None]-self.model_time[None,(self.model_time>timechunks[tc])&(self.model_time<=timechunks[tc+1])]),axis=1) for tc in range(nchunks)]) + if np.all([np.any((self.lc.time>timechunks[tc])&(self.lc.time<=timechunks[tc+1])) for tc in range(nchunks)]): + min_dist_to_lc=np.hstack([np.min(abs(self.lc.time[(self.lc.time>timechunks[tc])&(self.lc.time<=timechunks[tc+1]),None]-self.model_time[None,(self.model_time>timechunks[tc])&(self.model_time<=timechunks[tc+1])]),axis=1) for tc in range(nchunks)]) + else: + print("NO TIME HERE?",timechunks[tc]) self.gp_to_plot['gp_sd'] = np.tile(np.nanmedian(abs(np.diff(self.lc.bin_flux))),len(self.lc.time))*(np.clip(86400/1800*min_dist_to_lc,1.0,25)**0.33) elif newgp: for key in self.lc_regions['limits']: @@ -2705,7 +2777,7 @@ def init_gp_to_plot(self, n_samp=7, max_gp_len=12000, interp=True, newgp=False, elif n_samp>1: assert hasattr(self,'trace') if interp: - assert self.bin_oot + #assert self.bin_oot stacktime=np.hstack((self.lc.time[0]-1,self.model_time,self.lc.time[-1]+1)) preds=[] for i in np.random.choice(len(self.trace['phot_mean']),int(np.clip(10*n_samp,1,len(self.trace['phot_mean']))),replace=False): @@ -4078,7 +4150,7 @@ def Plot(self, interactive=False, n_samp=None, overwrite=False, interp=True, new spl+self.trans_to_plot['all']['allpl']['-1sig'][ix], spl+self.trans_to_plot['all']['allpl']['+1sig'][ix], alpha=0.3, color="C0",zorder=11, rasterized=raster) - if not self.fit_no_flatten: + if not self.fit_no_flatten and not self.fit_GP: f_alls[key].plot(self.lc.time[ix], spl+self.trans_to_plot['all']['allpl']['med'][ix], color="C0", label="transit fit", linewidth=2.5,alpha=0.5,zorder=12, rasterized=raster) @@ -4128,7 +4200,7 @@ def Plot(self, interactive=False, n_samp=None, overwrite=False, interp=True, new setattr(self.lc,'phase',{}) for n,pl in enumerate(self.planets): if hasattr(self,'trace'): - t0=np.nanmedian(self.trace['t0_'+pl]) + t0s=[np.nanmedian(self.trace['t0_'+pl]), np.nanmedian(self.trace['t0_2_'+pl]), np.nanmedian(self.trace['t0_3_'+pl])] if pl in self.trios else [np.nanmedian(self.trace['t0_'+pl])] if pl in self.multis or pl in self.rvplanets: per=np.nanmedian(self.trace['per_'+pl]) elif pl in self.duos+self.trios: @@ -4140,7 +4212,7 @@ def Plot(self, interactive=False, n_samp=None, overwrite=False, interp=True, new elif 'tdur_'+pl+'[0]' in self.init_soln: binsize=np.nanmedian(self.trace['tdur_'+pl+'[0]'])/n_intrans_bins elif hasattr(self,'init_soln'): - t0=self.init_soln['t0_'+pl] + t0s=[self.init_soln['t0_'+pl], self.init_soln['t0_2_'+pl], self.init_soln['t0_3_'+pl]] if pl in self.trios else [self.init_soln['t0_'+pl]] if pl in self.multis or pl in self.rvplanets: per=self.init_soln['per_'+pl] elif pl in self.duos+self.trios: @@ -4151,30 +4223,30 @@ def Plot(self, interactive=False, n_samp=None, overwrite=False, interp=True, new binsize=self.init_soln['tdur_'+pl]/n_intrans_bins elif 'tdur_'+pl+'[0]' in self.init_soln: binsize=self.init_soln['tdur_'+pl+'[0]']/n_intrans_bins - self.lc.phase[pl]=(self.lc.time-t0-0.5*per)%per-0.5*per + self.lc.phase[pl]=self.make_phase(self.lc.time,t0s,per) for nkey,key in enumerate(self.lc_regions): - if pl in self.multis or pl in self.rvplanets or pl in self.duos: + if pl in self.multis or pl in self.rvplanets or pl in self.duos or pl in self.trios: #print(key,pl,per,t0) - n_p_sta_end=np.array([np.floor((self.lc_regions[key]['start']-t0)/per),np.ceil((self.lc_regions[key]['end']-t0)/per)]) + n_p_sta_end=np.array([np.floor((self.lc_regions[key]['start']-np.min(t0s))/per),np.ceil((self.lc_regions[key]['end']-np.max(t0s))/per)]) #Adding ticks for the position of each planet below the data: trans_range=np.arange(n_p_sta_end[0],n_p_sta_end[1],1.0) if interactive: - f_alls[key].scatter(t0+trans_range*per,np.tile(self.lc_regions[key]['minmax_global'][0]+0.2*resid_sd+(0.5*resid_sd*n/len(self.planets)), + f_alls[key].scatter(np.min(t0s)+trans_range*per,np.tile(self.lc_regions[key]['minmax_global'][0]+0.2*resid_sd+(0.5*resid_sd*n/len(self.planets)), int(len(trans_range))), marker="triangle", s=12.5, fill_color=self.pal[4-n], alpha=0.85) else: - f_alls[key].scatter(t0+trans_range*per, + f_alls[key].scatter(np.min(t0s)+trans_range*per, np.tile(self.lc_regions[key]['minmax_global'][0]+0.2*resid_sd+(resid_sd*n/len(self.planets)), int(len(trans_range))), marker="^", s=12.5, color=self.pal[4-n], alpha=0.85) elif pl in self.monos: - if (t0>self.lc_regions[key]['start'])&(t0self.lc_regions[key]['start'])&(t0s[0]self.lc_regions[key]['start'])&(t01e-5",-4:"p>1e-4",-3:"p>0.1%",-2:"p>1%",-1:"p>10%",0:"p>100%"} plot_pers=self.duos+self.monos+self.trios ymin = 0-15*float(int(ylog)) if ymin is None else ymin - if len(plot_pers)>0: - plt.figure(figsize=(8.5,4.2)) + #plt.figure(figsize=(4.5+1.5*np.sqrt(len(self.monos+self.duos+self.trios)),4.2)) + plt.figure(figsize=(3.2+1.5*np.sqrt(len(self.monos+self.duos+self.trios)),3.5)) for npl, pl in enumerate(plot_pers): plt.subplot(1,len(plot_pers),npl+1) if pl in self.duos+self.trios: @@ -4394,7 +4487,7 @@ def PlotPeriods(self, plot_loc=None, ylog=True, xlog=True, nbins=25, cols+=[ncol] plt.plot(np.tile(pers[n],2), [ymin,probs[n]], - linewidth=5.0,color=pal[6+ncol],alpha=0.7,label=coldic[ncol]) + linewidth=5.0,color=pal[6+ncol],alpha=0.7,label="$"+coldic[ncol].replace('%','\%')+"$") else: plt.plot(np.tile(pers[n],2), [ymin,probs[n]], @@ -4405,8 +4498,10 @@ def PlotPeriods(self, plot_loc=None, ylog=True, xlog=True, nbins=25, plt.legend() if xlog: plt.xscale('log') - plt.xticks([20,30,40,60,80,100,150,200,250,300,350,400,450,500,600,800,1000,1500,2000,2500,3000,3500], - np.array([20,30,40,60,80,100,150,200,250,300,350,400,450,500,600,800,1000,1500,2000,2500,3000,3500]).astype(str)) + plt.xticks([20,30,40,60,80,100,150,200,300,400,600,800,1000,1500,2000,2500,3000,3500], + np.array([20,30,40,60,80,100,150,200,300,400,600,800,1000,1500,2000,2500,3000,3500]).astype(str)) + #plt.xticks([20,30,40,60,80,100,150,200,250,300,350,400,450,500,600,800,1000,1500,2000,2500,3000,3500], + # np.array([20,30,40,60,80,100,150,200,250,300,350,400,450,500,600,800,1000,1500,2000,2500,3000,3500]).astype(str)) #plt.xticklabels([20,40,60,80,100,150,200,250]) plt.ylabel("$\log_{10}{p}$") plt.xlabel("Period [d]") @@ -4420,12 +4515,12 @@ def PlotPeriods(self, plot_loc=None, ylog=True, xlog=True, nbins=25, pmax = np.nanmax(self.trace['per_'+pl].ravel()) if pmax is None else pmax pmin = 0.5*np.min(self.planets[pl]['per_gaps']['gap_starts']) if pmin is None else pmin cols=[] - for ngap in np.arange(self.planets[pl]['ngaps']): + for ngap in np.arange(self.planets[pl]['ngaps'])[np.argsort(np.nanmedian(self.trace['logprob_marg_'+pl],axis=0))]: if self.planets[pl]['per_gaps']['gap_starts'][ngap]pmin: - bins=np.arange(np.floor(self.planets[pl]['per_gaps']['gap_starts'][ngap]), - np.clip(np.ceil(self.planets[pl]['per_gaps']['gap_ends'][ngap])+1.0,0.0,pmax), + bins=np.arange(np.floor(self.planets[pl]['per_gaps']['gap_starts'][ngap])-0.5, + np.clip(np.ceil(self.planets[pl]['per_gaps']['gap_ends'][ngap])+0.5,0.0,pmax), 1.0) - + print(bins) ncol=int(np.floor(np.clip(np.nanmedian(self.trace['logprob_marg_'+pl][:,ngap])-total_av_prob,-6,0))) #print(self.planets[pl]['per_gaps']['gap_starts'][ngap], # ncol,np.nanmedian(self.trace['logprob_marg_'+pl][:,ngap])-total_av_prob) @@ -4462,7 +4557,7 @@ def PlotPeriods(self, plot_loc=None, ylog=True, xlog=True, nbins=25, #plt.xlim(60,80) plt.xlabel("Period [d]") plt.legend(title="Average prob") - + plt.tight_layout() if plot_loc is None: plt.savefig(self.savenames[0]+'_period_dists.pdf') else: @@ -4745,7 +4840,7 @@ def PredictFutureTransits(self, time_start=None, time_end=None, time_dur=180, in time_start = Time(datetime.now()) if time_end is None: time_end = Time(datetime.now())+time_dur*u.day - + if type(time_start)==Time: print("time range",time_start.isot,"->",time_end.isot) time_start = time_start.jd - self.lc.jd_base @@ -4760,6 +4855,7 @@ def PredictFutureTransits(self, time_start=None, time_end=None, time_dur=180, in assert type(time_start) in [int, np.int64, float, np.float64] print("time range",Time(time_start+self.lc.jd_base,format='jd').isot, "->",Time(time_end+self.lc.jd_base,format='jd').isot) + if check_TESS: sect_start_ends=self.CheckTESS() @@ -4875,6 +4971,12 @@ def CheckTESS(self,**kwargs): return np.column_stack((midtimes-0.5*sectdiffs[future_sect_ix]+0.2,midtimes+0.5*sectdiffs[future_sect_ix]-0.2)) #Now we have sector start & end times, let's check which future transit will be TESS observed: + def CheopsRMS(self, Gmag, tdur): + #RMS polynomial fits for 3 hour durations: + rms_brightfit = np.array([ 2.49847572, -6.41232409]) + rms_faintfit = np.array([ 30.2599025 , -256.41381477]) + rms = np.max([np.polyval(rms_faintfit,Gmag),np.polyval(rms_brightfit,Gmag)]) + return rms/np.sqrt(tdur/0.125) def MakeCheopsOR(self, DR2ID=None, pl=None, min_eff=45, oot_min_orbits=1.0, timing_sigma=3, t_start=None, t_end=None, Texp=None, max_orbits=14, min_pretrans_orbits=0.5, min_intrans_orbits=None, orbits_flex=1.4, observe_sigma=2, @@ -5012,6 +5114,10 @@ def MakeCheopsOR(self, DR2ID=None, pl=None, min_eff=45, oot_min_orbits=1.0, timi if max_ORsobserve_threshold)>max_ORs and ipl not in self.multis: observe_threshold=np.sort(allprobs)[::-1][max_ORs] + depth=1e6*np.nanmedian(self.trace['ror_'+ipl])**2 + print("SNR for whole transit is: ",depth/self.CheopsRMS(gaiainfo['phot_g_mean_mag'], np.nanmedian(self.trace['tdur_'+ipl]))) + print("SNR for single orbit in/egress is: ",depth/self.CheopsRMS(gaiainfo['phot_g_mean_mag'], 0.5*98/1440)) + prio_1_prob_threshold = np.ceil(np.sum(allprobs>observe_threshold)*prio_1_threshold) prio_3_prob_threshold = np.ceil(np.sum(allprobs>observe_threshold)*(1-prio_3_threshold)) print(allprobs,observe_threshold,allpers[allprobs>observe_threshold]) diff --git a/MonoTools/lightcurve.py b/MonoTools/lightcurve.py index eddcf51..7bfd86a 100755 --- a/MonoTools/lightcurve.py +++ b/MonoTools/lightcurve.py @@ -1161,7 +1161,6 @@ def get_tess_lc(self,sector,search=['all'],use_fast=True,use_eleanor=False,**kwa Returns: lightcurve.lc: TESS lightcurve """ - print(sector,search) #['all','spoc_20','spoc_120','spoc_1800','qlp_1800','eleanor_1800'] #use_ppt=True, coords=None, use_qlp=None, use_eleanor=None, data_loc=None, search_fast=False, **kwargs): diff --git a/__init__.py b/__init__.py old mode 100644 new mode 100755 diff --git a/pyproject.toml b/pyproject.toml old mode 100644 new mode 100755 diff --git a/setup.py b/setup.py index 8f25f5c..ff73153 100755 --- a/setup.py +++ b/setup.py @@ -39,7 +39,7 @@ packages=setuptools.find_packages(), package_data={'MonoTools': extrafiles}, install_requires=['matplotlib', - 'numpy==1.22', + 'numpy==1.21', 'pandas', 'scipy', 'astropy',