In [1]:

    

#Load a specific lightcurve file for testing purposes
import json
import sys
import numpy as np
import emcee
import george
import matplotlib.pyplot as plt
sys.path.append('../classification')
import os
filename = '../gen_lightcurves/gp_smoothed/SN2005el_gpsmoothed.json'
with open(filename, 'r') as f:
     file_data = json.load(f)
print(file_data.keys())


    

    




dict_keys(['B_', 'B__CSP', 'B_kait3', 'H_PAIRITEL', 'H__CSP', 'I_kait3', 'J_PAIRITEL', 'J__CSP', 'Ks_PAIRITEL', 'R_kait3', 'U_', 'V_', 'V__CSP', 'V_kait3', 'Y__CSP', 'g__CSP', 'i__CSP', 'r__CSP', 'u__CSP'])

In [ ]:

    

import bandMap
lcurve_mapped = bandMap.remapBands(file_data)
print(lcurve_mapped.keys())


    

In [ ]:

    

for filt in file_data:
    mjd = np.array(file_data[filt]['mjd'])
    mag = np.array(file_data[filt]['mag'])
    magerr = np.array(file_data[filt]['dmag'])
    modelmjd = file_data[filt]['modeldate']
    modelmag = file_data[filt]['modelmag']
    
    plt.errorbar(mjd,mag,yerr=magerr)
    plt.plot(modelmjd, modelmag)
    plt.title(filt)
    plt.show()


    

In [2]:

    

print(file_data['B_'].keys())
filter_name = 'I_kait3'
t = np.array(file_data[filter_name]['mjd'])
y = np.array(file_data[filter_name]['mag'])
err = np.array(file_data[filter_name]['dmag'])
plt.plot(t,y, 'ko')
plt.show()


    

    




dict_keys(['bsplinemag', 'dmag', 'goodstatus', 'kernel', 'mag', 'mjd', 'modeldate', 'modelerr', 'modelmag', 'modelmag_sub', 'type'])






    

In [11]:

    

from george import kernels
kernel = kernels.ConstantKernel(1*10**-3) * kernels.ExpSquaredKernel(100.) * kernels.DotProductKernel()
#kernel = kernels.ExpSquaredKernel(100.) * kernels.DotProductKernel()


    

In [12]:

    

import scipy.optimize as op

gp = george.GP(kernel, mean=np.mean(y))

# Define the objective function (negative log-likelihood in this case).
def nll(p):
    # Update the kernel parameters and compute the likelihood.
    gp.kernel[:] = p
    ll = gp.lnlikelihood(y, quiet=True)

    # The scipy optimizer doesn't play well with infinities.
    return -ll if np.isfinite(ll) else 1e25

# And the gradient of the objective function.
def grad_nll(p):
    # Update the kernel parameters and compute the likelihood.
    gp.kernel[:] = p
    return -gp.grad_lnlikelihood(y, quiet=True)

# You need to compute the GP once before starting the optimization.
gp.compute(t, err)

# Print the initial ln-likelihood.
print(gp.lnlikelihood(y))

# Run the optimization routine.
p0 = gp.kernel.vector
results = op.minimize(nll, p0, jac=grad_nll)

# Update the kernel and print the final log-likelihood.
gp.kernel[:] = results.x
print(gp.lnlikelihood(y))


    

    




-159.010392971
-76.1852862535

In [13]:

    

x = np.linspace(min(t), max(t), 100)
mu, cov = gp.predict(y,x)
std = np.sqrt(np.diag(cov))
plt.plot(x,mu)
plt.errorbar(t,y,yerr=err)
plt.show()


    

In [6]:

    

kernel.pars


    

    
Out[6]:





array([  1.00000000e-03,   1.00000000e+02])

In [7]:

    

gp = george.GP(kernel, mean=np.mean(y))
gp.compute(t, err)
print(gp.lnlikelihood(y))
print(gp.grad_lnlikelihood(y))

def lnprob(p):
    # Trivial improper prior: uniform in the log.
    if np.any((-10 > p) + (p > 10)):
        return -np.inf
    lnprior = 0.0

    # Update the kernel and compute the lnlikelihood.
    kernel.pars = np.exp(p)
    return lnprior + gp.lnlikelihood(y, quiet=True)
# You need to compute the GP once before starting. Then the sample list
# will be saved.

# Set up the sampler.
nwalkers, ndim = 36, len(kernel)
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob)

# Initialize the walkers.
p0 = [np.log(kernel.pars) + 1e-4 * np.random.randn(ndim)
      for i in range(nwalkers)]

print("Running burn-in")
p0, _, _ = sampler.run_mcmc(p0, 2000)

print("Running production chain")
sampler.run_mcmc(p0, 2000)


    

    




-91.0271154876
[ -9.41100448  45.55883867]
Running burn-in
Running production chain






    
Out[7]:





(array([[-9.85354098,  5.85154889],
        [-9.68128942,  5.81284183],
        [-9.94465188,  5.67354078],
        [-9.86822772,  5.87683923],
        [-9.88359671,  5.95511675],
        [-9.99874192,  5.71609669],
        [-9.95765161,  5.87066769],
        [-9.62664767,  5.96913593],
        [-9.90069598,  5.95875506],
        [-9.95989114,  5.57405727],
        [-9.95202784,  5.85492569],
        [-9.71421877,  5.75505987],
        [-9.97271626,  5.9008392 ],
        [-9.78184753,  5.92124371],
        [-9.25298183,  6.09209703],
        [-9.89405351,  6.03361663],
        [-9.87915766,  5.96556995],
        [-9.99744347,  5.71632428],
        [-9.91225281,  5.96458277],
        [-9.92162343,  5.90276465],
        [-9.85977183,  5.76984065],
        [-9.2381504 ,  5.90752193],
        [-9.82680942,  5.88855358],
        [-9.63890607,  5.59971144],
        [-9.68975533,  5.93891808],
        [-9.99876437,  5.78716909],
        [-9.85281704,  5.90428807],
        [-9.84070529,  5.69102542],
        [-9.73720292,  5.55036946],
        [-9.98844071,  5.76972759],
        [-9.840558  ,  5.83154688],
        [-9.72780855,  5.74233812],
        [-9.9120763 ,  6.09215482],
        [-9.74852173,  5.67452675],
        [-9.59379378,  5.85074837],
        [-9.69038775,  5.90469119]]),
 array([-39.39731764, -40.28041181, -40.08134042, -39.3292441 ,
        -39.59772337, -39.38609303, -38.97212846, -40.51400019,
        -39.57649112, -41.33224416, -38.98839684, -40.49820364,
        -38.98463656, -39.74575421, -42.91439153, -40.59507911,
        -39.70199305, -39.39075761, -39.59588384, -39.17630686,
        -39.67931234, -41.90741492, -39.50236857, -42.82473465,
        -40.14168005, -38.93931268, -39.43371935, -40.44688545,
        -43.05225771, -39.06988013, -39.49223977, -40.53601098,
        -41.83542742, -41.13512658, -40.51931849, -40.0516312 ]),
 ('MT19937', array([1046097296, 2648881685, 1233255482, 3264740253, 1913203028,
         2596676463, 2286397863,  310909077, 1537929949, 1356794061,
          437583856, 2752422509, 1897149676,   72028463, 2059935134,
         1080185334, 2789140688, 2220186527, 2841877375, 1014418599,
         3655696615, 2797054365, 3975629184, 1558655195, 3285859947,
         4223106420, 1290864562, 4217364073,  151492094, 2171049631,
         1896562182, 3579452787, 1102722289,   81392530, 3171727746,
         3581104289, 2854481076,  539362043, 1409524868, 3807931139,
         2470652721, 2554031803, 1938890518, 3241401642, 1746217406,
         3654872962, 4261379267, 2963362250, 3999465736,  945116028,
         2376283716, 1187453562, 3155976788,  730843374, 1361442044,
          464962950, 3780241470, 3158097649, 3293630761, 3556898841,
         3962235078, 3920264626,  740026288, 3562858237, 3975427521,
         2651061783, 2116455659, 2106007003, 3353172417,  694049559,
         2677471872, 3633587586, 3292622421, 2165646271, 2527314462,
         3707649180, 1086550758, 3919136151,  887629642,  776252688,
         1703437324, 1293276800, 2421035620, 1808114921, 2962210779,
          396833365, 2692874044, 3665124184, 1323490100, 2126697859,
         3357370577, 2888757300, 2830820617, 3738848617, 1956729753,
         2112333857, 3455710230, 1306613208,  570334966, 2068111282,
           91936260,  731735682,  953928587, 2459314239, 4052618752,
          233306106, 3662028466, 2023527789, 2389284351,   38308699,
         3020589523,  523614388, 2054455155, 1249741315, 1415142284,
          783276859, 2084868811,  516047716, 3656533840, 1655319309,
          627373902, 3021309919, 2024326590, 4077953382,  999663508,
         2040333918, 3693786754, 3559712727, 3667807986, 3207834296,
          607690846,  701123691,  183543890,  623219767,  473653949,
           36015729, 2322601498, 1073229363, 1907200730, 4230216187,
         1987747316, 2356035236, 2389329802, 1626659860,  250993580,
         2463387353, 2462543744, 3659614859, 4226434025, 1185152730,
         2211419166,  925917206, 2050751187, 3674762372, 1518948292,
          318044299, 2831229985, 2287960170, 1050222422, 1592743007,
          412708966, 1598213707, 1457209533,  129807265,  898404551,
         2507193043, 2541216730, 3977333481,  549687916, 4190265694,
         2446993860,  472058198, 4261521776, 3999959919, 1437744260,
         3157662673,  863773807, 3476211996, 2387366511, 2241953636,
         2990628638, 1760651582,  446265766, 1375457148,  149952335,
         1606551934,  768163822, 4251069235, 4074489244, 3745538390,
         1962181778, 2957420538, 1832074687, 1787854097, 1837123463,
         2973254627, 3641349655, 4221686555,  591067809, 4185681931,
          594747948, 3262582553,  413487646, 1943219359, 4110061384,
         2998515124, 1110730161,  715185694,  779109334, 2618536230,
          238435782, 3054691208, 4021561696, 2500148026, 2034200488,
         1134199390,  538596186, 2547993023,   63908374, 4167281310,
         2819340214,   16282229, 3709797977, 2630654952, 3695487049,
         1611389299,  611334735, 1598496365, 3256240117,  489471236,
         1432760204, 3335012207, 1440247910, 1581656418,  165472530,
         2911971687, 1638434673, 3230827009, 3543111350, 2247705330,
         1833025778, 1382694985, 2460448383,  491389498, 2588876029,
          365255171, 4137690881,  935869953, 1886250805, 2659231915,
         2324158958, 4107483919, 1249410166,  135879263,  100528232,
         2666105314, 1821545535, 2309848254, 3365331604, 1974490397,
         3655893032,  342050877, 1885047831,  528114548, 3266079719,
         1411875297,  404480340, 3617793673, 2441474688,  710425427,
         3285118247, 4239566905, 4016523291, 2905263461, 2514218511,
         2831203225, 2089040420,  537881159, 1654324632, 2942187194,
          736871581, 3290826499, 3827182863, 4183048147,  562158307,
         2379080845, 1137162757, 3881446188, 1751170003, 2765499177,
         1978424596, 4234536295,  585629012, 2873088534, 2718599423,
         1956048671,  765175744, 1421904541,  759090627,  868821590,
         2817311188,  986582187, 1505094876, 2616718049, 3532490763,
          287851533, 2148797188,  125668768, 1936883537,  303458563,
          383821088, 2802752377, 3398912265, 4239382754, 4107228461,
         2130948476, 2746817372, 1566487848, 3596861718, 3080458825,
         3402526413, 1123093596, 3004889795,  362435289, 1192696924,
         2259600409,   64339244,  361000807, 1925276286, 3490614985,
          740150462, 1830614750, 2836789929, 2365758879, 3574314397,
         3386055357, 4212473260, 2703414901,  846375880, 1799582455,
         2862982333, 1323965391, 3020494623, 3586799673,  536820327,
         1790288142, 2721171931, 3083505617,  122508174, 3254383169,
         1368469008, 3881096209, 4274514912,  217829909, 4180989800,
          536628096,  862065238, 1176585353, 2431298445, 1017468744,
          945846019, 2659326170, 1220694890,  657643237, 3215259456,
          532729999, 2292042399,  364061632, 1006777135, 3554148313,
         2016766404, 2171733701, 3582782679, 2475101265,  191435917,
         4012447224,  224833171, 4140042556, 2923609843, 2298130359,
         1435151234, 2187509886, 3479330983, 3547475846, 2421342298,
         1888682106, 2548840072, 4194590925, 2318161722, 3309251608,
         3660401757, 1460346669, 3195370790, 4026904172, 3067890921,
         3858362192, 3089271327, 3162049636, 4119986556, 1531950767,
         3169189185, 4091467946, 1835690472,  438308475, 1332353613,
         2810865803,  108528302, 2794011628,  165477012, 4159671419,
          312768103,  665879456,  557288463,  787519214,  796890443,
         4156564226, 4035350792, 2903265845, 4236407674,  986650467,
           16328182, 3005296620, 2595336200, 1273815480,  497595370,
         2816681890, 3349004891, 1749721188,  187856339, 3711846399,
          151241239, 3944143846, 1608803709, 2306050287, 3468552041,
         2592388550, 3130722053,  336807262, 2764661312, 3554666510,
          370077960, 3340910110, 2305337451, 1659404597, 3030487124,
          476556645,  158044683, 4216998647, 3127920093, 1719767143,
          138702977, 2556800631, 2730795618, 1048519012,  402174927,
         1068276451, 3005819505, 3320638464, 1156721523, 3567475181,
         2930594270,  450283476, 3480241656, 1503773128, 2085395249,
         3502653117, 1315429254, 3110814615, 4101668088,   31974587,
         3014245347,  572847665,  525973889, 2639360939, 1993802977,
         4067157484, 1664079318,  884227186, 1026562264, 2747534219,
         1513913480, 1040068124,  631456287,  199551533, 2711199931,
         3717195973, 1970196465, 2050908199, 2203318041,  996362102,
          450432128, 3229432280, 3889597818,  325131663, 1870472916,
         1977394448, 2620412085, 2585184385, 1714431877, 2978283590,
         1406099199, 1670761954, 1293625379,  212647755, 3259582863,
         4273181408, 1465178421, 3397402007, 1449194566, 3709918046,
         1607793205, 2554805932, 1129090100,  126993944, 1226329885,
         2275863811, 4258761521, 1499010402,  109146897, 2334551562,
          346244996, 2378709696, 3140106935, 4157881780, 3700679192,
         1997806425,   56236169,  923025383, 4120318880, 1612718619,
         3056201217, 1796019171, 1095181692, 3780951999, 3034459741,
         1587421976, 3016195884, 2352871429, 3946176827,   70386237,
         2394961327,  577820680, 2256761086,  349933455, 1137188868,
         1326280836, 3085669327,  441011629, 1459177862, 1854881946,
         1779855348,   78886160,   41640942, 2946961331,  353052432,
         2574571448,  980401057,  781308701, 1358071203, 1176153018,
         1026559546, 2562154780,  171464852,   54617817, 1480868366,
         3458176028, 1792433476,  830923862,  113086598, 2577667667,
         1718499997, 3638592214, 1128986265,   71355407, 1312060944,
         1962272847, 3542916749,   63993424, 3399824415, 1106111042,
          268209257, 1833221414, 2770741122, 3127147288, 3646111612,
         2640465173, 1982680836, 2868322076, 3240259067, 1125431802,
         2066726469,  867719965, 1200915936, 1195658101,  668637092,
         3129223704, 3672293694, 2207706438,  692066602, 1702298574,
         3400656876, 2035082172, 1276951799, 1374362514, 3834676748,
         3919682589, 3453292145, 2688945667, 1268004914, 1794678597,
         1042807928, 3686635487,  134443170, 2744772364, 1593111624,
         4106380967, 1526945894, 2310356961,  203893331, 3458258010,
           25742950, 3194188814,  480397720, 2802818263], dtype=uint32), 461, 0, 0.0))

In [8]:

    

x = np.linspace(min(t),max(t), 100)
chain = sampler.chain
for i in range(50):
    # Choose a random walker and step.
    w = np.random.randint(chain.shape[0])
    n = np.random.randint(2000, chain.shape[1])
    gp.kernel.pars = np.exp(chain[w, n])

    # Plot a single sample.
    plt.plot(x, gp.sample_conditional(y, x), "k", alpha=0.3)
plt.show()


    

    




/home/levicivita/.local/lib/python3.6/site-packages/george/utils.py:30: RuntimeWarning: covariance is not positive-semidefinite.
  samples = np.random.multivariate_normal(mean, matrix, N)






    

In [ ]: