optimize_esfi.py

import os
import urllib.request, json
import pandas as pd
import numpy as np
import random
import george
from george.kernels import ExpSquaredKernel, ExpSine2Kernel, Matern32Kernel, ConstantKernel
from kernel import esfi_kernel, esfi_without_daily

import scipy.optimize as op

iter, fev, grev, best = 0, 0, 0, 9999

i = 0
tm = []
sfi = []

def get_data(url):
    with urllib.request.urlopen(url) as res:
        data = json.loads(res.read().decode())

    return data

data = get_data('http://localhost:%s/mixscale_essn.json?series=6h&var=sfi&points=8000' % (os.getenv('HISTORY_PORT')))
tm = [ x[0] / 86400. for x in data ]
sfi = [ x[1] for x in data ]

first_tm = tm[0]
last_tm = tm[len(tm)-1]
span = last_tm - first_tm

tm = np.array(tm)
sfi = np.array(sfi)
mean_sfi = np.mean(sfi)
sfi -= mean_sfi

gp = george.GP(esfi_kernel, white_noise=np.log(0.1**2), fit_white_noise=True)
gp.compute(tm)

def loss(p):
    ret = 0
    gp.set_parameter_vector(p)
    ll = gp.log_likelihood(sfi, quiet=True)
    return -ll if np.isfinite(ll) else 1e25

def nll(p):
    global fev
    global best
    fev = fev + 1

    lsum = loss(p)

    if lsum < best:
        best = lsum
    return lsum

def grad(p):
    gp.set_parameter_vector(p)
    return -gp.grad_log_likelihood(sfi, quiet=True)
    return ret

def grad_nll(p):
    global grev
    grev = grev + 1

    return grad(p)


def cb(p):
    global iter
    iter = iter + 1
    print("# iter=", iter, "fev=", fev, "grev=", grev, "best=", best, "p=", repr(p))

p0 = gp.get_parameter_vector()
print("# Init: ", p0)

bounds = gp.get_parameter_bounds()
print("# Bounds: ", bounds)

opt_result = op.minimize(nll, p0, jac=grad_nll, method='L-BFGS-B', callback=cb, options={'maxiter': 100}, bounds=bounds)
print("# RESULT ", repr(opt_result.x))

gp.set_parameter_vector(opt_result.x)
print(esfi_kernel)