Skip to content

Instantly share code, notes, and snippets.

@arodland
Last active January 12, 2022 21:08
Show Gist options
  • Save arodland/53f228bf631504aa2e02ff8be2c57304 to your computer and use it in GitHub Desktop.
Save arodland/53f228bf631504aa2e02ff8be2c57304 to your computer and use it in GitHub Desktop.
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment