sparta/extract_phase_curve.py at master · ideasrule/sparta · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
import numpy as np
import sys
import matplotlib
#matplotlib.use("Agg")
import matplotlib.pyplot as plt
import batman
from algorithms import reject_beginning, bin_data, get_mad, \
    get_data_pickle, print_stats
import emcee
import argparse
from configparser import ConfigParser
from emcee_methods import lnprob as lnprob
from emcee_methods import get_batman_params, run_emcee
import time
import corner

def correct_lc(start_wave, end_wave, fluxes, errors, bjds, y, x, t0, per, rp, a, inc,
               limb_dark_coeffs, fp, C1, D1, C2, D2, output_file_prefix="chain",
               nwalkers=100, burn_in_runs=100, production_runs=1000, extra_phase_terms=False, output_txt="white_light_result.txt"):
    print("Median is", np.median(fluxes))
    fluxes /= np.median(fluxes)

    #We initialize transit and eclipse models so that lnprob doesn't have to.
    #Speeds up lnprob dramatically.
    initial_batman_params = get_batman_params(t0, per, rp, a, inc, limb_dark_coeffs) #, limb_dark_law="quadratic")
    transit_model = batman.TransitModel(initial_batman_params, bjds)
    eclipse_model = batman.TransitModel(initial_batman_params, bjds, transittype='secondary')

    #First guess for PLD corrections based on PCA components
    error_factor = 1.1 #Initial guess; slightly larger than 1 is usually good
    print("Error factor", error_factor)
    b = a*np.cos(inc*np.pi/180)
    if extra_phase_terms:
        initial_params = np.array([0, 0, fp, C1, D1, C2, D2, rp, a, b, error_factor, 1, 0, 0.1, 0, 0, 0])
        labels = ["delta_t0", "delta_te", "Fp", "C1", "D1", "C2", "D2", "rp", "a", "b", "error", "Fstar", "Aramp", "tau", "cy", "cx", "m"]
        rp_index = 7
    else:
        initial_params = np.array([0, 0, fp, C1, D1, rp, error_factor, 1, 0, 0.1, 0, 0, 0])
        labels = ["delta_t0", "delta_te", "Fp", "C1", "D1", "rp", "error", "Fstar", "Aramp", "tau", "cy", "cx", "m"]
        rp_index = 5

    #All arguments, aside from the parameters, that will be passed to lnprob
    w = 2*np.pi/per
    lnprob_args = (initial_batman_params, transit_model, eclipse_model, bjds, fluxes, errors, y, x, t0, extra_phase_terms)

    _, chain, lnprobs = run_emcee(lnprob, lnprob_args, initial_params, nwalkers, output_file_prefix, burn_in_runs, production_runs)
    length = len(chain)
    chain = chain[int(length/2):]
    lnprobs = lnprobs[int(length/2):]
    best_step = chain[np.argmax(lnprobs)]
    best_lnprob = lnprobs[np.argmax(lnprobs)]

    print("Best step", best_step)
    _, residuals = lnprob(best_step, *lnprob_args, plot_result=True, return_residuals=True)

    A = np.sqrt(chain[:,3]**2 + chain[:,4]**2)
    phi = np.arctan2(chain[:,4], chain[:,3]) * 180 / np.pi

    for i in range(chain.shape[1]):
        print_stats(chain[:,i], labels[i])

    print_stats(A, "A")
    print_stats(phi, "phi")

    plt.figure()
    corner.corner(chain, range=[0.99] * chain.shape[1], labels=labels)
    plt.show()

    night_Fp = chain[:,2] - 2*chain[:,3]

    with open(output_txt, "w") as f:
        f.write("#min_wavelength max_wavelength A_med A_lower_err A_upper_err phi_med phi_lower_err phi_upper_err Fp_med Fp_lower_err Fp_upper_err RpRs_med RpRs_lower_err RpRs_upper_err night_Fp_med night_Fp_lower_err night_Fp_upper_err lnprob\n")
        f.write("{} {} ".format(start_wave, end_wave))
        for var in (A, phi, chain[:,2], chain[:,rp_index], night_Fp):
            f.write("{} {} {} ".format(np.median(var), np.median(var) - np.percentile(var, 16), np.percentile(var, 84) - np.median(var)))
        f.write("{}".format(best_lnprob))
        f.write("\n")
    return chain, lnprobs

parser = argparse.ArgumentParser(description="Extracts phase curve and transit information from light curves")
parser.add_argument("config_file", help="Contains transit, eclipse, and phase curve parameters")
parser.add_argument("start_wave", type=float)
parser.add_argument("end_wave", type=float)
parser.add_argument("-b", "--bin-size", type=int, default=1, help="Bin size to use on data")
parser.add_argument("--burn-in-runs", type=int, default=2000, help="Number of burn in runs")
parser.add_argument("--production-runs", type=int, default=1000, help="Number of production runs")
parser.add_argument("--num-walkers", type=int, default=100, help="Number of walkers in the ensemble sampler")
parser.add_argument("-o", "--output", type=str, default="chain", help="Directory to store the chain and lnprob arrays")
parser.add_argument("--extra-phase-terms", action="store_true", help="Include C2 and D2 in phase curve fit")
parser.add_argument("-e", "--exclude", type=int, default=525, help="How many integrations to trim at the beginning")

args = parser.parse_args()

bjds, fluxes, flux_errors, wavelengths, y, x = get_data_pickle(args.start_wave, args.end_wave, args.exclude)

bin_size = args.bin_size
binned_fluxes = bin_data(fluxes, bin_size)
factor = np.median(binned_fluxes)
binned_fluxes /= factor
binned_errors = np.sqrt(bin_data(flux_errors**2, bin_size) / bin_size) / factor
binned_bjds = bin_data(bjds, bin_size)
binned_x = bin_data(x, bin_size)
binned_y = bin_data(y, bin_size)

delta_t = binned_bjds - np.median(binned_bjds)
coeffs = np.polyfit(delta_t, binned_x, 1)
smoothed_binned_x = np.polyval(coeffs, delta_t)

binned_y = bin_data(y, bin_size)
coeffs = np.polyfit(delta_t, binned_y, 1)
smoothed_binned_y = np.polyval(coeffs, delta_t)


#plt.scatter(binned_bjds, binned_fluxes, s=0.5, alpha=0.5)
#plt.show()

print("Num points", len(binned_fluxes))
#get values from configuration file
default_section_name = "DEFAULT"
config = ConfigParser()
config.read(args.config_file)
items = dict(config.items(default_section_name))

chain, lnprobs = correct_lc(args.start_wave/1000, args.end_wave/1000, binned_fluxes, binned_errors,
                            binned_bjds, binned_y - smoothed_binned_y, binned_x - smoothed_binned_x,
                            float(items["t0"]), float(items["per"]),
                            float(items["rp"]), float(items["a"]), float(items["inc"]), eval(items["limb_dark_coeffs"]),
                            float(items["fp"]), float(items["c1"]), float(items["d1"]), float(items["c2"]), float(items["d2"]),
                            args.output, args.num_walkers, args.burn_in_runs, args.production_runs, args.extra_phase_terms)
np.save("white_chain.npy", chain)
np.save("white_lnprobs.npy", lnprobs)