Merge pull request #951 from rzellem/develop

v1.10.0 - updates to priors

Author: rzellem

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()
+
+
+

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()