From 5859ab0eaaae070c0900021123d210dcff9131a1 Mon Sep 17 00:00:00 2001 From: Kyle Pearson Date: Fri, 13 Oct 2023 15:45:17 -0700 Subject: [PATCH] cleaned up ephemeris fitting code --- exotic/api/ephemeris.py | 515 ++++++----------------------- exotic/api/nested_linear_fitter.py | 270 +-------------- 2 files changed, 105 insertions(+), 680 deletions(-) diff --git a/exotic/api/ephemeris.py b/exotic/api/ephemeris.py index 208a2af8..f9bff13b 100644 --- a/exotic/api/ephemeris.py +++ b/exotic/api/ephemeris.py @@ -925,363 +925,117 @@ def plot_triangle(self): return fig - -def empty_data(N): # orbit parameters in Nbody simulation - return { - 'x':np.zeros(N), - 'y':np.zeros(N), - 'z':np.zeros(N), - 'P':np.zeros(N), - 'a':np.zeros(N), - 'e':np.zeros(N), - 'inc':np.zeros(N), - 'Omega':np.zeros(N), - 'omega':np.zeros(N), - 'M':np.zeros(N), - } - -def integrate(sim, objects, Ndays, Noutputs): - # ps is now an array of pointers and will change as the simulation runs - ps = sim.particles - - times = np.linspace(0., Ndays, Noutputs) # units of days - - pdata = [empty_data(Noutputs) for i in range(len(objects)-1)] - star = {'x':np.zeros(Noutputs),'y':np.zeros(Noutputs),'z':np.zeros(Noutputs) } - - # run simulation time steps - for i,time in enumerate(times): - sim.integrate(time) - - # record data - for k in star.keys(): - star[k][i] = getattr(ps[0], k) - - for j in range(1,len(objects)): - for k in pdata[j-1].keys(): - pdata[j-1][k][i] = getattr(ps[j],k) - - sim_data = ({ - 'pdata':pdata, - 'star':star, - 'times':times, - 'objects':objects, - 'dt': Noutputs/(Ndays*24*60*60) # conversion factor to get seconds for RV semi-amp - }) - - return sim_data - -def generate(objects, Ndays=None, Noutputs=None): - # create rebound simulation - # for object parameters see: - # https://rebound.readthedocs.io/en/latest/_modules/rebound/particle.html - sim = rebound.Simulation() - - # sim.integrator = "tes" - # sim.ri_tes.dq_max = 1e-2 - # sim.ri_tes.recti_per_orbit = 1.61803398875 - # sim.ri_tes.epsilon = 1e-6 - - # sim.integrator = "whfast" - # sim.ri_whfast.corrector = 5 - # sim.ri_whfast.safe_mode = 0 - sim.units = ('day', 'AU', 'Msun') - for i in range(len(objects)): - sim.add( **objects[i] ) - sim.move_to_com() - - if Ndays and Noutputs: - return integrate(sim, objects, Ndays, Noutputs) - else: - return sim - -def find_zero(t1,dx1, t2,dx2): - # find zero with linear interpolation - m = (dx2-dx1)/(t2-t1) - T0 = -dx1/m + t1 - return T0 - -def transit_times(xp,xs,times): - # check for sign change in position difference - dx = xp-xs - tt = [] - for i in range(1,len(dx)): - if dx[i-1] >= 0 and dx[i] <= 0: - tt.append( find_zero(times[i-1],dx[i-1], times[i],dx[i]) ) - return np.array(tt) - -def TTV(epochs, tt): - N = len(epochs) - A = np.vstack([np.ones(N), epochs]).T - b, m = np.linalg.lstsq(A, tt, rcond=None)[0] - ttv = (tt-m*np.array(epochs)-b) - return [ttv,m,b] - -class nbody_fitter(): - - def __init__(self, data, prior=None, bounds=None, verbose=True): - self.data = data - self.bounds = bounds - self.prior = prior - self.verbose = verbose - self.fit_nested() - - def fit_nested(self): - - # set up some arrays for mapping sampler output - freekeys = [] - boundarray = [] - for i,planet in enumerate(self.bounds): - for bound in planet: - freekeys.append(f"{i}_{bound}") - boundarray.append(planet[bound]) - - # find min and max time for simulation - min_time = np.min(self.data[1]['Tc']) - max_time = np.max(self.data[1]['Tc']) - # todo extend for multiplanet systems - - sim_time = (max_time-min_time)*1.05 - # TODO extend for multiplanet systems - Tc_norm = self.data[1]['Tc'] - min_time # normalize the data to the first observation - self.orbit = np.rint(Tc_norm / self.prior[1]['P']).astype(int) # number of orbits since first observation (rounded to nearest integer) - - # numpify - boundarray = np.array(boundarray) - bounddiff = np.diff(boundarray,1).reshape(-1) - - # create queue and save simulation to - def loglike(pars): - chi2 = 0 - - # set parameters - for i,par in enumerate(pars): - idx,key = freekeys[i].split('_') - idx = int(idx) - if key == 'tmid': - continue - # this dict goes to REBOUND and needs specific keys - self.prior[idx][key] = par - - # take difference between data and simulation - # run N-body simulation - sim_data = generate(self.prior, sim_time, int(sim_time*24)) # uses linspace behind the scenes - - # json structure with analytics of interest from the simulation - # ttv_data = analyze(sim_data) # slow - - sim_shift = 0 - # loop over planets, check for data - for i,planet in enumerate(self.prior): - if self.data[i]: - # compute transit times from N-body simulation - Tc_sim = transit_times( sim_data['pdata'][i-1]['x'], sim_data['star']['x'], sim_data['times'] ) - - # derive an offset in time from the first planet - if i-1==0: - sim_shift = Tc_sim.min() - - # shift the first mid-transit in simulation to observation - Tc_sim -= sim_shift # first orbit has tmid at 0 - - # scale Tc_sim to data - try: - residual = self.data[i]['Tc'] - Tc_sim[self.orbit] - except: - #import pdb; pdb.set_trace - #ttv,m,b = TTV(np.arange(Tc_sim.shape[0]), Tc_sim) - #ttv1,m1,b1 = TTV(self.orbit, self.data[i]['Tc']) - # could be unstable orbit or not enough data - # switch to average error and clip by max epoch? - print(self.prior) - chi2 += -1e6 - continue - - Tc_sim += residual.mean() - - # take difference between data and simulation - try: - chi2 += -0.5*np.sum(((self.data[i]['Tc'] - Tc_sim[self.orbit])/self.data[i]['Tc_err'])**2) - except: - chi2 += -1e6 - print(self.prior) - # usually unstable orbit - - return chi2 - - - def prior_transform(upars): - return (boundarray[:,0] + bounddiff*upars) - - if self.verbose: - self.results = ReactiveNestedSampler(freekeys, loglike, prior_transform).run(max_ncalls=1e5) - else: - self.results = ReactiveNestedSampler(freekeys, loglike, prior_transform).run(max_ncalls=1e5, show_status=self.verbose, viz_callback=self.verbose) - - self.errors = {} - self.quantiles = {} - self.parameters = copy.deepcopy(self.prior) - - #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]] - - -#def main(): if __name__ == "__main__": Tc = np.array([ # measured mid-transit times #Exoplanet Watch transit times array -2457025.7867,2457347.7655,2457694.8263,2458112.8495,2458492.6454,2459532.7834,2459604.81036,2459604.8137,2459614.6365,2459616.8209, -2459651.7466,2459903.8588,2459914.7783,2459915.8684,2459925.6921,2459939.8793,2459949.7047,2459959.5249,2459962.8037,2459973.7129, -2459974.798,2459986.8057,2459994.4489,2460009.7312,2458867.01587,2459501.1295, - -#ExoClock transit times array -2454508.977,2454515.525,2454801.477,2454840.769,2455123.447,2455140.91,2455147.459,2455147.459,2455148.551,2455148.552,2455158.372, -2455158.373,2455159.464,2455159.467,2455160.556,2455163.831,2455172.562,2455192.205,2455209.669,2455210.762,2455215.129,2455230.407, -2455238.046,2455254.419,2455265.331,2455494.53,2455494.531,2455509.81,2455510.902,2455530.547,2455542.553,2455548.011,2455565.472, -2455566.563,2455575.296,2455575.297,2455578.569,2455589.483,2455590.576,2455598.216,2455600.397,2455600.398,2455600.398,2455601.49, -2455601.49,2455601.49,2455603.673,2455612.404,2455623.318,2455623.319,2455624.41,2455624.41,2455663.702,2455842.693,2455842.694, -2455876.528,2455887.442,2455887.442,2455887.442,2455888.533,2455888.534,2455888.534,2455889.625,2455890.716,2455901.63,2455903.814, -2455903.814,2455920.185,2455923.459,2455926.733,2455939.829,2455945.287,2455946.378,2455946.379,2455947.47,2455948.561,2455951.835, -2455952.927,2455959.475,2455960.567,2455970.39,2455971.481,2455981.304,2455982.395,2455983.487,2455984.578,2455985.67,2455996.584, -2456000.95,2456005.315,2456006.406,2456014.046,2456175.577,2456245.427,2456249.794,2456273.805,2456282.536,2456284.719,2456297.816, -2456302.182,2456305.455,2456319.644,2456328.376,2456329.467,2456604.505,2456605.596,2456606.688,2456607.779,2456629.607,2456630.699, -2456654.71,2456659.076,2456662.35,2456663.441,2456664.533,2456674.356,2456677.63,2456688.544,2456694.002,2456703.824,2456711.464, -2456719.104,2456721.287,2456722.378,2456986.502,2457010.513,2457012.696,2457033.434,2457045.438,2457046.53,2457059.627,2457060.718, -2457067.267,2457068.358,2457103.284,2457345.579,2457368.5,2457378.321,2457390.327,2457391.418,2457426.343,2457426.344,2457427.435, -2457451.446,2457474.364,2457486.372,2457671.913,2457691.559,2457703.564,2457706.838,2457726.484,2457727.575,2457760.318,2457765.775, -2457766.866,2457772.324,2457773.415,2457776.689,2457786.512,2457788.695,2457800.7,2457808.34,2457809.432,2457809.434,2457810.523, -2457832.348,2457832.351,2458026.624,2458050.635,2458060.459,2458072.462,2458073.555,2458074.647,2458077.921,2458123.76,2458124.852, -2458132.491,2458134.675,2458136.858,2458155.41,2458155.412,2458156.503,2458159.778,2458166.326,2458178.331,2458214.347,2458411.895, -2458467.557,2458471.923,2458478.472,2458479.564,2458490.477,2458494.843,2458501.39,2458501.392,2458506.848,2458536.315,2458537.408, -2458871.382,2458880.113,2458883.384,2458893.209,2458903.032,2459088.575,2459148.602,2459156.242,2459157.334,2459190.076,2459194.441, -2459194.443,2459196.623,2459197.717,2459204.264, - -#ETD mid transit times array -2454508.976854,2454836.403393,2454840.768603,2454840.771193,2454908.437992,2454931.358182,2455164.923956,2455164.924316,2455172.561816,2455172.562786, -2455198.756735,2455230.406730,2455254.418870,2455256.596033,2455265.334053,2455506.533285,2455509.808915,2455519.632074,2455532.729964,2455541.461084, -2455576.387242,2455577.475172,2455577.476782,2455577.477112,2455578.568422,2455593.850192,2455600.398282,2455601.489211,2455603.671161,2455612.400801, -2455612.401321,2455615.675911,2455647.328760,2455857.973783,2455857.973783,2455867.797832,2455873.253612,2455877.618482,2455877.619162,2455889.623302, -2455899.448441,2455899.451251,2455900.539781,2455904.905271,2455905.995241,2455912.544561,2455912.544561,2455913.636131,2455914.726741,2455922.370121, -2455928.913321,2455936.552560,2455936.561520,2455937.645430,2455938.738730,2455946.383030,2455947.468480,2455949.650900,2455949.650900,2455950.742510, -2455957.287600,2455957.294400,2455957.296780,2455957.296780,2455961.659130,2455962.749930,2455970.392369,2455982.394919,2455993.309869,2455993.313119, -2455994.398719,2455994.399359,2456003.129689,2456003.130169,2456005.315328,2456222.509695,2456246.516835,2456249.793845,2456269.425195,2456271.623265, -2456296.724384,2456304.364534,2456304.365804,2456304.367454,2456304.367784,2456329.469294,2456364.392544,2456584.860173,2456595.772813,2456596.866383, -2456604.503073,2456608.869564,2456616.511354,2456626.333324,2456627.429034,2456628.515004,2456628.517074,2456630.699754,2456639.430344,2456639.432894, -2456640.522554,2456641.610274,2456650.341584,2456663.443614,2456688.544734,2456722.372535,2456722.382615,2456734.383505,2456734.388675,2456746.385515, -2456746.391885,2456927.546138,2456950.485119,2456960.306869,2456963.578439,2456963.578959,2456963.582479,2456986.501140,2456998.508180,2456998.508180, -2457032.333311,2457033.433561,2457069.448912,2457080.367572,2457080.367572,2457092.368353,2457344.486723,2457344.488693,2457345.578813,2457345.579833, -2457356.493733,2457357.585013,2457357.586493,2457368.498964,2457368.499654,2457369.589474,2457379.414424,2457379.414424,2457414.334726,2457426.343356, -2457451.447237,2457474.364078,2457474.364348,2457486.371248,2457668.637786,2457726.484138,2457749.405049,2457750.494949,2457751.585649,2457751.586939, -2457757.043281,2457758.134851,2457760.317471,2457760.318171,2457761.408951,2457770.140361,2457772.327372,2457773.415162,2457774.506522,2457809.431273, -2457809.435043,2457832.350874,2458087.745843,2458096.474134,2458109.571734,2458131.400315,2458396.615162,2458396.617562,2458411.896052,2458455.554593, -2458457.734253,2458479.568743,2458489.383454,2458489.387804,2458490.477904,2458490.482824,2458501.391674,2458501.392014,2458501.396544,2458513.402254, -2458537.407634,2458538.507734,2458837.550156,2458848.461686,2458859.375956,2458860.465876,2458860.469616,2458871.379136,2458871.382556,2458872.470796, -2458883.387165,2458883.389555,2458884.478175,2458884.478425,2458884.479455,2458906.306925,2458930.332215,2459147.510251,2459171.522921,2459172.608741, -2459172.617851,2459183.531910,2459193.352480,2459194.442340,2459194.443610,2459196.623950,2459197.719250,2459220.632699,2459242.461499,2459265.384258, -2459265.384318,2459505.495609,2459553.514627,2459553.515137,2459564.427746,2459565.520376,2459566.615006,2459576.438096]) + 2457025.7867,2457347.7655,2457694.8263,2458112.8495,2458492.6454,2459532.7834,2459604.81036,2459604.8137,2459614.6365,2459616.8209, + 2459651.7466,2459903.8588,2459914.7783,2459915.8684,2459925.6921,2459939.8793,2459949.7047,2459959.5249,2459962.8037,2459973.7129, + 2459974.798,2459986.8057,2459994.4489,2460009.7312,2458867.01587,2459501.1295, + + #ExoClock transit times array + 2454508.977,2454515.525,2454801.477,2454840.769,2455123.447,2455140.91,2455147.459,2455147.459,2455148.551,2455148.552,2455158.372, + 2455158.373,2455159.464,2455159.467,2455160.556,2455163.831,2455172.562,2455192.205,2455209.669,2455210.762,2455215.129,2455230.407, + 2455238.046,2455254.419,2455265.331,2455494.53,2455494.531,2455509.81,2455510.902,2455530.547,2455542.553,2455548.011,2455565.472, + 2455566.563,2455575.296,2455575.297,2455578.569,2455589.483,2455590.576,2455598.216,2455600.397,2455600.398,2455600.398,2455601.49, + 2455601.49,2455601.49,2455603.673,2455612.404,2455623.318,2455623.319,2455624.41,2455624.41,2455663.702,2455842.693,2455842.694, + 2455876.528,2455887.442,2455887.442,2455887.442,2455888.533,2455888.534,2455888.534,2455889.625,2455890.716,2455901.63,2455903.814, + 2455903.814,2455920.185,2455923.459,2455926.733,2455939.829,2455945.287,2455946.378,2455946.379,2455947.47,2455948.561,2455951.835, + 2455952.927,2455959.475,2455960.567,2455970.39,2455971.481,2455981.304,2455982.395,2455983.487,2455984.578,2455985.67,2455996.584, + 2456000.95,2456005.315,2456006.406,2456014.046,2456175.577,2456245.427,2456249.794,2456273.805,2456282.536,2456284.719,2456297.816, + 2456302.182,2456305.455,2456319.644,2456328.376,2456329.467,2456604.505,2456605.596,2456606.688,2456607.779,2456629.607,2456630.699, + 2456654.71,2456659.076,2456662.35,2456663.441,2456664.533,2456674.356,2456677.63,2456688.544,2456694.002,2456703.824,2456711.464, + 2456719.104,2456721.287,2456722.378,2456986.502,2457010.513,2457012.696,2457033.434,2457045.438,2457046.53,2457059.627,2457060.718, + 2457067.267,2457068.358,2457103.284,2457345.579,2457368.5,2457378.321,2457390.327,2457391.418,2457426.343,2457426.344,2457427.435, + 2457451.446,2457474.364,2457486.372,2457671.913,2457691.559,2457703.564,2457706.838,2457726.484,2457727.575,2457760.318,2457765.775, + 2457766.866,2457772.324,2457773.415,2457776.689,2457786.512,2457788.695,2457800.7,2457808.34,2457809.432,2457809.434,2457810.523, + 2457832.348,2457832.351,2458026.624,2458050.635,2458060.459,2458072.462,2458073.555,2458074.647,2458077.921,2458123.76,2458124.852, + 2458132.491,2458134.675,2458136.858,2458155.41,2458155.412,2458156.503,2458159.778,2458166.326,2458178.331,2458214.347,2458411.895, + 2458467.557,2458471.923,2458478.472,2458479.564,2458490.477,2458494.843,2458501.39,2458501.392,2458506.848,2458536.315,2458537.408, + 2458871.382,2458880.113,2458883.384,2458893.209,2458903.032,2459088.575,2459148.602,2459156.242,2459157.334,2459190.076,2459194.441, + 2459194.443,2459196.623,2459197.717,2459204.264, + + #ETD mid transit times array + 2454508.976854,2454836.403393,2454840.768603,2454840.771193,2454908.437992,2454931.358182,2455164.923956,2455164.924316,2455172.561816,2455172.562786, + 2455198.756735,2455230.406730,2455254.418870,2455256.596033,2455265.334053,2455506.533285,2455509.808915,2455519.632074,2455532.729964,2455541.461084, + 2455576.387242,2455577.475172,2455577.476782,2455577.477112,2455578.568422,2455593.850192,2455600.398282,2455601.489211,2455603.671161,2455612.400801, + 2455612.401321,2455615.675911,2455647.328760,2455857.973783,2455857.973783,2455867.797832,2455873.253612,2455877.618482,2455877.619162,2455889.623302, + 2455899.448441,2455899.451251,2455900.539781,2455904.905271,2455905.995241,2455912.544561,2455912.544561,2455913.636131,2455914.726741,2455922.370121, + 2455928.913321,2455936.552560,2455936.561520,2455937.645430,2455938.738730,2455946.383030,2455947.468480,2455949.650900,2455949.650900,2455950.742510, + 2455957.287600,2455957.294400,2455957.296780,2455957.296780,2455961.659130,2455962.749930,2455970.392369,2455982.394919,2455993.309869,2455993.313119, + 2455994.398719,2455994.399359,2456003.129689,2456003.130169,2456005.315328,2456222.509695,2456246.516835,2456249.793845,2456269.425195,2456271.623265, + 2456296.724384,2456304.364534,2456304.365804,2456304.367454,2456304.367784,2456329.469294,2456364.392544,2456584.860173,2456595.772813,2456596.866383, + 2456604.503073,2456608.869564,2456616.511354,2456626.333324,2456627.429034,2456628.515004,2456628.517074,2456630.699754,2456639.430344,2456639.432894, + 2456640.522554,2456641.610274,2456650.341584,2456663.443614,2456688.544734,2456722.372535,2456722.382615,2456734.383505,2456734.388675,2456746.385515, + 2456746.391885,2456927.546138,2456950.485119,2456960.306869,2456963.578439,2456963.578959,2456963.582479,2456986.501140,2456998.508180,2456998.508180, + 2457032.333311,2457033.433561,2457069.448912,2457080.367572,2457080.367572,2457092.368353,2457344.486723,2457344.488693,2457345.578813,2457345.579833, + 2457356.493733,2457357.585013,2457357.586493,2457368.498964,2457368.499654,2457369.589474,2457379.414424,2457379.414424,2457414.334726,2457426.343356, + 2457451.447237,2457474.364078,2457474.364348,2457486.371248,2457668.637786,2457726.484138,2457749.405049,2457750.494949,2457751.585649,2457751.586939, + 2457757.043281,2457758.134851,2457760.317471,2457760.318171,2457761.408951,2457770.140361,2457772.327372,2457773.415162,2457774.506522,2457809.431273, + 2457809.435043,2457832.350874,2458087.745843,2458096.474134,2458109.571734,2458131.400315,2458396.615162,2458396.617562,2458411.896052,2458455.554593, + 2458457.734253,2458479.568743,2458489.383454,2458489.387804,2458490.477904,2458490.482824,2458501.391674,2458501.392014,2458501.396544,2458513.402254, + 2458537.407634,2458538.507734,2458837.550156,2458848.461686,2458859.375956,2458860.465876,2458860.469616,2458871.379136,2458871.382556,2458872.470796, + 2458883.387165,2458883.389555,2458884.478175,2458884.478425,2458884.479455,2458906.306925,2458930.332215,2459147.510251,2459171.522921,2459172.608741, + 2459172.617851,2459183.531910,2459193.352480,2459194.442340,2459194.443610,2459196.623950,2459197.719250,2459220.632699,2459242.461499,2459265.384258, + 2459265.384318,2459505.495609,2459553.514627,2459553.515137,2459564.427746,2459565.520376,2459566.615006,2459576.438096]) Tc_error = np.array([ 0.0039,0.0058,0.0054,0.0022,0.0065,0.0024,0.00091,0.00097,0.0022,0.0018, -0.002,0.0021,0.0065,0.0036,0.0071,0.0031,0.0044,0.0021,0.0028,0.0023, -0.0039,0.0071,0.0015,0.0055,0.00092,0.00092, - -#ExoClock error array -0.0002,0.00013,0.00041,0.00047,0.0007,0.000417,0.0009,0.00042,0.001,0.0013,0.00057, -0.0013,0.00095,0.00061,0.0007,0.000324,0.00044,0.0005,0.000463,0.000405,0.00082,0.00011, -0.00066,0.00014,0.0016,0.00048,0.00063,0.00037,0.000313,0.00077,0.00029,0.00024,0.0011, -0.00022,0.00063,0.00079,0.0008,0.00094,0.001,0.00036,0.00083,0.00026,0.00078,0.000289, -0.00017,0.00039,0.001,0.00088,0.0009,0.00054,0.00086,0.0014,0.00063,0.00093,0.0013, -0.0003,0.00038,0.00062,0.00075,0.00059,0.0003,0.00038,0.0009,0.00023,0.00083,0.00093, -0.000324,0.00046,0.0002,0.00092,0.0019,0.00079,0.00036,0.00034,0.00022,0.00028,0.00011, -0.0001,0.00018,0.0004,0.00029,0.00029,0.00094,0.00047,0.00029,0.000324,0.000417,0.00037, -0.0004,0.00038,0.0004,0.00023,0.00033,0.00033,0.000394,0.000301,0.0003,0.000301,0.000301, -0.00046,0.00026,0.000382,0.00027,0.00029,0.0002,0.0003,0.00034,0.000706,0.00019,0.00043, -0.000336,0.00034,0.00019,0.00019,0.00032,0.00028,0.000324,0.00041,0.00029,0.00029,0.00026, -0.00034,0.00034,0.00046,0.00043,0.00039,0.000486,0.0005,0.00049,0.00049,0.000347,0.000359, -0.00022,0.00021,0.0003,0.00042,0.0004,0.0013,0.00034,0.00033,0.00055,0.0006,0.00023,0.00021, -0.0007,0.0013,0.00035,0.00025,0.00034,0.00037,0.00028,0.00023,0.0006,0.00028,0.00039, -0.00024,0.00022,0.00029,0.00026,0.00048,0.00032,0.0004,0.00018,0.0009,0.00021,0.0006, -0.0006,0.00056,0.00023,0.0003,0.0003,0.00022,0.00034,0.00028,0.00027,0.00035,0.00031, -0.00032,0.00033,0.0005,0.00031,0.00032,0.00091,0.00034,0.00038,0.0017,0.0004,0.0005, -0.00026,0.0006,0.0006,0.0008,0.0003,0.0009,0.0003,0.00044,0.0008,0.0007,0.0009, -0.0003,0.0007,0.0005,0.0003,0.0013,0.0007,0.0003,0.0004,0.0003,0.0012,0.0006, -0.0005,0.0013,0.0004, - -#ETD error array -0.000200,0.000600,0.000470,0.001000,0.001000,0.000980,0.001490,0.001290,0.000440,0.000140, -0.001410,0.000110,0.000140,0.000590,0.001290,0.001120,0.000870,0.000700,0.000350,0.000500, -0.001350,0.000380,0.000690,0.000850,0.000820,0.000470,0.000420,0.001090,0.000940,0.000830, -0.000780,0.000540,0.001730,0.000710,0.000710,0.000750,0.001260,0.000920,0.001290,0.000730, -0.000630,0.000570,0.000380,0.001030,0.001350,0.000570,0.000570,0.000780,0.000470,0.000900, -0.000730,0.000910,0.001230,0.001190,0.000720,0.000770,0.001020,0.000590,0.000590,0.000660, -0.001100,0.001040,0.000570,0.000570,0.001070,0.001320,0.000860,0.001160,0.000600,0.000760, -0.000680,0.000760,0.000630,0.000600,0.000440,0.000810,0.000740,0.000670,0.000900,0.000550, -0.000520,0.001460,0.000890,0.001560,0.000580,0.001640,0.001170,0.000510,0.000960,0.000510, -0.000920,0.000710,0.000900,0.000510,0.001050,0.000970,0.000880,0.000440,0.000740,0.000680, -0.000820,0.000800,0.001040,0.000760,0.000540,0.000890,0.001080,0.001120,0.000730,0.001390, -0.001410,0.001640,0.000740,0.000570,0.000930,0.000890,0.000610,0.000510,0.001450,0.001450, -0.001130,0.000460,0.000510,0.001190,0.001190,0.000600,0.000650,0.001070,0.001090,0.000490, -0.000610,0.000520,0.000430,0.000580,0.000460,0.000340,0.000890,0.000890,0.001310,0.000650, -0.000860,0.000770,0.000550,0.001350,0.000820,0.000550,0.000840,0.000530,0.000940,0.000720, -0.000370,0.000520,0.000670,0.001030,0.000630,0.000330,0.000990,0.000550,0.000620,0.000720, -0.000940,0.000580,0.000350,0.000570,0.001380,0.000640,0.001180,0.000700,0.000700,0.000710, -0.000550,0.000580,0.000680,0.001030,0.000860,0.000850,0.000300,0.000760,0.000680,0.000820, -0.000610,0.001450,0.000730,0.000700,0.001350,0.000820,0.000670,0.000940,0.000490,0.001130, -0.000540,0.000540,0.000700,0.000790,0.000840,0.000600,0.000520,0.000730,0.000640,0.001020, -0.000780,0.000610,0.001330,0.000770,0.000610,0.000520,0.001130,0.001130,0.000530,0.000780, -0.000420,0.001250,0.000380,0.000720,0.000860,0.000470,0.000950,0.000540]) + 0.002,0.0021,0.0065,0.0036,0.0071,0.0031,0.0044,0.0021,0.0028,0.0023, + 0.0039,0.0071,0.0015,0.0055,0.00092,0.00092, + + #ExoClock error array + 0.0002,0.00013,0.00041,0.00047,0.0007,0.000417,0.0009,0.00042,0.001,0.0013,0.00057, + 0.0013,0.00095,0.00061,0.0007,0.000324,0.00044,0.0005,0.000463,0.000405,0.00082,0.00011, + 0.00066,0.00014,0.0016,0.00048,0.00063,0.00037,0.000313,0.00077,0.00029,0.00024,0.0011, + 0.00022,0.00063,0.00079,0.0008,0.00094,0.001,0.00036,0.00083,0.00026,0.00078,0.000289, + 0.00017,0.00039,0.001,0.00088,0.0009,0.00054,0.00086,0.0014,0.00063,0.00093,0.0013, + 0.0003,0.00038,0.00062,0.00075,0.00059,0.0003,0.00038,0.0009,0.00023,0.00083,0.00093, + 0.000324,0.00046,0.0002,0.00092,0.0019,0.00079,0.00036,0.00034,0.00022,0.00028,0.00011, + 0.0001,0.00018,0.0004,0.00029,0.00029,0.00094,0.00047,0.00029,0.000324,0.000417,0.00037, + 0.0004,0.00038,0.0004,0.00023,0.00033,0.00033,0.000394,0.000301,0.0003,0.000301,0.000301, + 0.00046,0.00026,0.000382,0.00027,0.00029,0.0002,0.0003,0.00034,0.000706,0.00019,0.00043, + 0.000336,0.00034,0.00019,0.00019,0.00032,0.00028,0.000324,0.00041,0.00029,0.00029,0.00026, + 0.00034,0.00034,0.00046,0.00043,0.00039,0.000486,0.0005,0.00049,0.00049,0.000347,0.000359, + 0.00022,0.00021,0.0003,0.00042,0.0004,0.0013,0.00034,0.00033,0.00055,0.0006,0.00023,0.00021, + 0.0007,0.0013,0.00035,0.00025,0.00034,0.00037,0.00028,0.00023,0.0006,0.00028,0.00039, + 0.00024,0.00022,0.00029,0.00026,0.00048,0.00032,0.0004,0.00018,0.0009,0.00021,0.0006, + 0.0006,0.00056,0.00023,0.0003,0.0003,0.00022,0.00034,0.00028,0.00027,0.00035,0.00031, + 0.00032,0.00033,0.0005,0.00031,0.00032,0.00091,0.00034,0.00038,0.0017,0.0004,0.0005, + 0.00026,0.0006,0.0006,0.0008,0.0003,0.0009,0.0003,0.00044,0.0008,0.0007,0.0009, + 0.0003,0.0007,0.0005,0.0003,0.0013,0.0007,0.0003,0.0004,0.0003,0.0012,0.0006, + 0.0005,0.0013,0.0004, + + #ETD error array + 0.000200,0.000600,0.000470,0.001000,0.001000,0.000980,0.001490,0.001290,0.000440,0.000140, + 0.001410,0.000110,0.000140,0.000590,0.001290,0.001120,0.000870,0.000700,0.000350,0.000500, + 0.001350,0.000380,0.000690,0.000850,0.000820,0.000470,0.000420,0.001090,0.000940,0.000830, + 0.000780,0.000540,0.001730,0.000710,0.000710,0.000750,0.001260,0.000920,0.001290,0.000730, + 0.000630,0.000570,0.000380,0.001030,0.001350,0.000570,0.000570,0.000780,0.000470,0.000900, + 0.000730,0.000910,0.001230,0.001190,0.000720,0.000770,0.001020,0.000590,0.000590,0.000660, + 0.001100,0.001040,0.000570,0.000570,0.001070,0.001320,0.000860,0.001160,0.000600,0.000760, + 0.000680,0.000760,0.000630,0.000600,0.000440,0.000810,0.000740,0.000670,0.000900,0.000550, + 0.000520,0.001460,0.000890,0.001560,0.000580,0.001640,0.001170,0.000510,0.000960,0.000510, + 0.000920,0.000710,0.000900,0.000510,0.001050,0.000970,0.000880,0.000440,0.000740,0.000680, + 0.000820,0.000800,0.001040,0.000760,0.000540,0.000890,0.001080,0.001120,0.000730,0.001390, + 0.001410,0.001640,0.000740,0.000570,0.000930,0.000890,0.000610,0.000510,0.001450,0.001450, + 0.001130,0.000460,0.000510,0.001190,0.001190,0.000600,0.000650,0.001070,0.001090,0.000490, + 0.000610,0.000520,0.000430,0.000580,0.000460,0.000340,0.000890,0.000890,0.001310,0.000650, + 0.000860,0.000770,0.000550,0.001350,0.000820,0.000550,0.000840,0.000530,0.000940,0.000720, + 0.000370,0.000520,0.000670,0.001030,0.000630,0.000330,0.000990,0.000550,0.000620,0.000720, + 0.000940,0.000580,0.000350,0.000570,0.001380,0.000640,0.001180,0.000700,0.000700,0.000710, + 0.000550,0.000580,0.000680,0.001030,0.000860,0.000850,0.000300,0.000760,0.000680,0.000820, + 0.000610,0.001450,0.000730,0.000700,0.001350,0.000820,0.000670,0.000940,0.000490,0.001130, + 0.000540,0.000540,0.000700,0.000790,0.000840,0.000600,0.000520,0.000730,0.000640,0.001020, + 0.000780,0.000610,0.001330,0.000770,0.000610,0.000520,0.001130,0.001130,0.000530,0.000780, + 0.000420,0.001250,0.000380,0.000720,0.000860,0.000470,0.000950,0.000540]) labels = np.full(len(Tc), 'EpW') - # # Tc = np.array([ # measured mid-transit times - # # 2459150.837905, 2459524.851045, - # # 2459546.613126, 2459565.643663, 2459584.690470, 2459584.686476, - # # 2459909.739104, 2459957.337739, 2459169.880602, 2458416.424861, - # # 2458428.664794, 2459145.400124, 2458430.025044, 2459164.440695, - # # 2458415.064846, 2458435.465269, 2458412.344658, 2459150.840070, - # # 2459160.360691, 2458431.384897, 2459146.760284, 2459154.920516, - # # 2458417.784945, 2459161.720466, 2459167.160761, 2458427.304939, - # # 2458413.705181, 2459163.080605, 2459153.560229, 2459168.520861, - # # 2458425.945068, 2459148.120269, 2458434.105274, 2458432.745423, - # # 2459152.200558, 2459165.800918, 2459159.000645, 2459149.480289, - # # 2455870.450027, 2456663.347570, 2457420.884390, 2457658.888900]) - - # # Tc_error = np.array([ - # # 0.001674, 0.000715, - # # 0.000758, 0.001560, 0.000371, 0.001651, - # # 0.000955, 0.000176, 0.000154, 0.000357, - # # 0.000381, 0.000141, 0.000362, 0.000149, - # # 0.000336, 0.000368, 0.000379, 0.000153, - # # 0.000153, 0.000349, 0.000149, 0.000146, - # # 0.000385, 0.000146, 0.000153, 0.000360, - # # 0.000356, 0.000147, 0.000147, 0.000146, - # # 0.000363, 0.000142, 0.000357, 0.000368, - # # 0.000160, 0.000160, 0.000151, 0.000160, - # # 0.000140, 0.000120, 0.000800, 0.000140]) - - # # labels for a legend - # labels = np.array(['TESS', 'TESS', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar', - # 'TESS', 'EpW', 'ExoClock', 'Unistellar']) P = 1.0914203 # orbital period for your target Tc_norm = Tc - Tc.min() # normalize the data to the first observation - # print(Tc_norm) orbit = np.rint(Tc_norm / P) # number of orbits since first observation (rounded to nearest integer) - # print(orbit) # make a n x 2 matrix with 1's in the first column and values of orbit in the second A = np.vstack([np.ones(len(Tc)), orbit]).T @@ -1317,17 +1071,17 @@ def prior_transform(upars): # min and max values to search between for fitting bounds = { - 'm': [P - 0.1, P + 0.1], # orbital period - 'b': [intercept - 0.1, intercept + 0.1] # mid-transit time + 'P': [slope - 0.1, slope + 0.1], # orbital period + 'T0': [intercept - 0.1, intercept + 0.1] # mid-transit time } # used to plot red overlay in O-C figure prior = { - 'm': [slope, slope_std_dev], # value from WLS (replace with literature value) - 'b': [intercept, intercept_std_dev] # value from WLS (replace with literature value) + 'P': [slope, slope_std_dev], # value from WLS (replace with literature value) + 'T0': [intercept, intercept_std_dev] # value from WLS (replace with literature value) } - lf = linear_fitter(Tc, Tc_error, bounds, prior=prior, labels=labels) + lf = ephemeris_fitter(Tc, Tc_error, bounds, prior=prior, labels=labels) lf.plot_triangle() plt.subplots_adjust(top=0.9, hspace=0.2, wspace=0.2) @@ -1363,7 +1117,7 @@ def prior_transform(upars): 'dPdN': [0, 0.0001] # best guess } - nlf = non_linear_fitter(Tc, Tc_error, bounds, prior=prior, labels=labels) + nlf = decay_fitter(Tc, Tc_error, bounds, prior=prior, labels=labels) nlf.plot_triangle() plt.subplots_adjust(top=0.9, hspace=0.2, wspace=0.2) @@ -1382,65 +1136,4 @@ def prior_transform(upars): for key in nlf.parameters: print(f"Parameter {key} = {nlf.parameters[key]:.2e} +- {nlf.errors[key]:.2e}") - # TODO extend plot, BIC Values - - # should this be here? - dude() - # mearth = u.M_earth.to(u.kg).value - msun = u.M_sun.to(u.kg) - mjup = u.M_jup.to(u.kg) - - # Parameters for N-body Retrieval - # read more about adding particles: https://rebound.readthedocs.io/en/latest/addingparticles/ - nbody_prior = [ - {'m':1.12}, # star - - {'m':0.28*mjup/msun, - 'P':lf.parameters['P'], # inner planet - 'inc':3.14159/2, - 'e':0, - 'omega':0}, - - {'m':0.988*mjup/msun, - 'P':lf.parameters['P']*3, # outer planet - 'inc':3.14159/2, - 'e':0, - 'omega':0 }, - ] - - - # specify data to fit - data = [ - {}, # data for star (e.g. RV) - {'Tc':Tc, 'Tc_err':Tc_error}, # data for inner planet (e.g. Mid-transit times) - {} # data for outer planet (e.g. Mid-transit times) - ] - - - # set up where to look for a solution - bounds = [ - - {}, # no bounds on star parameters - - { # bounds for inner planet - 'P': [nbody_prior[1]['P']-0.025, nbody_prior[1]['P']+0.025], # based on solution from linear fit - # mid-transit is automatically adjusted in the fitter class - }, - - { # bounds for outer planet - 'P':[lf.parameters['P']*1.75, lf.parameters['P']*4.25], # orbital period [day] - 'm':[0.1*mjup/msun,1.5*mjup/msun], # mass [msun] - #'e':[0,0.05], # eccentricity - #'omega':[-20*np.pi/180, 20*np.pi/180] # arg of periastron [radians] - } - ] - - # run the fitter - nfit = nbody_fitter(data, nbody_prior, bounds) - # todo how could I phase fold the data to just shorter nbody simulations? - - print(nfit.parameters) - print(nfit.errors) - - -# main() + # TODO BIC Values \ No newline at end of file diff --git a/exotic/api/nested_linear_fitter.py b/exotic/api/nested_linear_fitter.py index dda4fa99..229d4c8a 100644 --- a/exotic/api/nested_linear_fitter.py +++ b/exotic/api/nested_linear_fitter.py @@ -48,7 +48,6 @@ from itertools import cycle import matplotlib.pyplot as plt import numpy as np -import rebound import statsmodels.api as sm try: from ultranest import ReactiveNestedSampler @@ -914,210 +913,6 @@ def plot_triangle(self): return fig -def empty_data(N): # orbit parameters in Nbody simulation - return { - 'x':np.zeros(N), - 'y':np.zeros(N), - 'z':np.zeros(N), - 'P':np.zeros(N), - 'a':np.zeros(N), - 'e':np.zeros(N), - 'inc':np.zeros(N), - 'Omega':np.zeros(N), - 'omega':np.zeros(N), - 'M':np.zeros(N), - } - -def integrate(sim, objects, Ndays, Noutputs): - # ps is now an array of pointers and will change as the simulation runs - ps = sim.particles - - times = np.linspace(0., Ndays, Noutputs) # units of days - - pdata = [empty_data(Noutputs) for i in range(len(objects)-1)] - star = {'x':np.zeros(Noutputs),'y':np.zeros(Noutputs),'z':np.zeros(Noutputs) } - - # run simulation time steps - for i,time in enumerate(times): - sim.integrate(time) - - # record data - for k in star.keys(): - star[k][i] = getattr(ps[0], k) - - for j in range(1,len(objects)): - for k in pdata[j-1].keys(): - pdata[j-1][k][i] = getattr(ps[j],k) - - sim_data = ({ - 'pdata':pdata, - 'star':star, - 'times':times, - 'objects':objects, - 'dt': Noutputs/(Ndays*24*60*60) # conversion factor to get seconds for RV semi-amp - }) - - return sim_data - -def generate(objects, Ndays=None, Noutputs=None): - # create rebound simulation - # for object parameters see: - # https://rebound.readthedocs.io/en/latest/_modules/rebound/particle.html - sim = rebound.Simulation() - - # sim.integrator = "tes" - # sim.ri_tes.dq_max = 1e-2 - # sim.ri_tes.recti_per_orbit = 1.61803398875 - # sim.ri_tes.epsilon = 1e-6 - - # sim.integrator = "whfast" - # sim.ri_whfast.corrector = 5 - # sim.ri_whfast.safe_mode = 0 - sim.units = ('day', 'AU', 'Msun') - for i in range(len(objects)): - sim.add( **objects[i] ) - sim.move_to_com() - - if Ndays and Noutputs: - return integrate(sim, objects, Ndays, Noutputs) - else: - return sim - -def find_zero(t1,dx1, t2,dx2): - # find zero with linear interpolation - m = (dx2-dx1)/(t2-t1) - t0 = -dx1/m + t1 - return t0 - -def transit_times(xp,xs,times): - # check for sign change in position difference - dx = xp-xs - tt = [] - for i in range(1,len(dx)): - if dx[i-1] >= 0 and dx[i] <= 0: - tt.append( find_zero(times[i-1],dx[i-1], times[i],dx[i]) ) - return np.array(tt) - -def TTV(epochs, tt): - N = len(epochs) - A = np.vstack([np.ones(N), epochs]).T - b, m = np.linalg.lstsq(A, tt, rcond=None)[0] - ttv = (tt-m*np.array(epochs)-b) - return [ttv,m,b] - -class nbody_fitter(): - - def __init__(self, data, prior=None, bounds=None, verbose=True): - self.data = data - self.bounds = bounds - self.prior = prior - self.verbose = verbose - self.fit_nested() - - def fit_nested(self): - - # set up some arrays for mapping sampler output - freekeys = [] - boundarray = [] - for i,planet in enumerate(self.bounds): - for bound in planet: - freekeys.append(f"{i}_{bound}") - boundarray.append(planet[bound]) - - # find min and max time for simulation - min_time = np.min(self.data[1]['Tc']) - max_time = np.max(self.data[1]['Tc']) - # todo extend for multiplanet systems - - sim_time = (max_time-min_time)*1.05 - # TODO extend for multiplanet systems - Tc_norm = self.data[1]['Tc'] - min_time # normalize the data to the first observation - self.orbit = np.rint(Tc_norm / self.prior[1]['P']).astype(int) # number of orbits since first observation (rounded to nearest integer) - - # numpify - boundarray = np.array(boundarray) - bounddiff = np.diff(boundarray,1).reshape(-1) - - # create queue and save simulation to - def loglike(pars): - chi2 = 0 - - # set parameters - for i,par in enumerate(pars): - idx,key = freekeys[i].split('_') - idx = int(idx) - if key == 'tmid': - continue - # this dict goes to REBOUND and needs specific keys - self.prior[idx][key] = par - - # take difference between data and simulation - # run N-body simulation - sim_data = generate(self.prior, sim_time, int(sim_time*24)) # uses linspace behind the scenes - - # json structure with analytics of interest from the simulation - # ttv_data = analyze(sim_data) # slow - - sim_shift = 0 - # loop over planets, check for data - for i,planet in enumerate(self.prior): - if self.data[i]: - # compute transit times from N-body simulation - Tc_sim = transit_times( sim_data['pdata'][i-1]['x'], sim_data['star']['x'], sim_data['times'] ) - - # derive an offset in time from the first planet - if i-1==0: - sim_shift = Tc_sim.min() - - # shift the first mid-transit in simulation to observation - Tc_sim -= sim_shift # first orbit has tmid at 0 - - # scale Tc_sim to data - try: - residual = self.data[i]['Tc'] - Tc_sim[self.orbit] - except: - #import pdb; pdb.set_trace - #ttv,m,b = TTV(np.arange(Tc_sim.shape[0]), Tc_sim) - #ttv1,m1,b1 = TTV(self.orbit, self.data[i]['Tc']) - # could be unstable orbit or not enough data - # switch to average error and clip by max epoch? - print(self.prior) - chi2 += -1e6 - continue - - Tc_sim += residual.mean() - - # take difference between data and simulation - try: - chi2 += -0.5*np.sum(((self.data[i]['Tc'] - Tc_sim[self.orbit])/self.data[i]['Tc_err'])**2) - except: - chi2 += -1e6 - print(self.prior) - # usually unstable orbit - - return chi2 - - - def prior_transform(upars): - return (boundarray[:,0] + bounddiff*upars) - - if self.verbose: - self.results = ReactiveNestedSampler(freekeys, loglike, prior_transform).run(max_ncalls=1e5) - else: - self.results = ReactiveNestedSampler(freekeys, loglike, prior_transform).run(max_ncalls=1e5, show_status=self.verbose, viz_callback=self.verbose) - - self.errors = {} - self.quantiles = {} - self.parameters = copy.deepcopy(self.prior) - - #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]] - - #def main(): if __name__ == "__main__": @@ -1368,67 +1163,4 @@ def prior_transform(upars): # for each free key print parameter and error for key in nlf.parameters: - print(f"Parameter {key} = {nlf.parameters[key]:.2e} +- {nlf.errors[key]:.2e}") - - # TODO extend plot, BIC Values - - # should this be here? - dude() - # mearth = u.M_earth.to(u.kg).value - msun = u.M_sun.to(u.kg) - mjup = u.M_jup.to(u.kg) - - # Parameters for N-body Retrieval - # read more about adding particles: https://rebound.readthedocs.io/en/latest/addingparticles/ - nbody_prior = [ - {'m':1.12}, # star - - {'m':0.28*mjup/msun, - 'P':lf.parameters['P'], # inner planet - 'inc':3.14159/2, - 'e':0, - 'omega':0}, - - {'m':0.988*mjup/msun, - 'P':lf.parameters['P']*3, # outer planet - 'inc':3.14159/2, - 'e':0, - 'omega':0 }, - ] - - - # specify data to fit - data = [ - {}, # data for star (e.g. RV) - {'Tc':Tc, 'Tc_err':Tc_error}, # data for inner planet (e.g. Mid-transit times) - {} # data for outer planet (e.g. Mid-transit times) - ] - - - # set up where to look for a solution - bounds = [ - - {}, # no bounds on star parameters - - { # bounds for inner planet - 'P': [nbody_prior[1]['P']-0.025, nbody_prior[1]['P']+0.025], # based on solution from linear fit - # mid-transit is automatically adjusted in the fitter class - }, - - { # bounds for outer planet - 'P':[lf.parameters['P']*1.75, lf.parameters['P']*4.25], # orbital period [day] - 'm':[0.1*mjup/msun,1.5*mjup/msun], # mass [msun] - #'e':[0,0.05], # eccentricity - #'omega':[-20*np.pi/180, 20*np.pi/180] # arg of periastron [radians] - } - ] - - # run the fitter - nfit = nbody_fitter(data, nbody_prior, bounds) - # todo how could I phase fold the data to just shorter nbody simulations? - - print(nfit.parameters) - print(nfit.errors) - - -# main() + print(f"Parameter {key} = {nlf.parameters[key]:.2e} +- {nlf.errors[key]:.2e}") \ No newline at end of file