Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
t = np.sort(np.random.rand(np.max(N)))
yerr = np.random.uniform(0.1, 0.2, len(t))
y = np.sin(t)
def numpy_compute(kernel, x, yerr):
K = kernel.get_value(x[:, None] - x[None, :])
K[np.diag_indices_from(K)] += yerr ** 2
return cho_factor(K)
def numpy_log_like(factor, y):
np.dot(y, cho_solve(factor, y)) + np.sum(np.log(np.diag(factor[0])))
for xi, j in enumerate(J):
kernel = terms.RealTerm(1.0, 0.1)
for k in range((2*j - 1) % 2):
kernel += terms.RealTerm(1.0, 0.1)
for k in range((2*j - 1) // 2):
kernel += terms.ComplexTerm(0.1, 2.0, 1.6)
coeffs = kernel.coefficients
assert 2*j == len(coeffs[0]) + 2*len(coeffs[2]), "Wrong number of terms"
if args.george:
george_kernel = None
for a, c in zip(*(coeffs[:2])):
k = CeleriteKernel(a=a, b=0.0, c=c, d=0.0)
george_kernel = k if george_kernel is None else george_kernel + k
for a, b, c, d in zip(*(coeffs[2:])):
k = CeleriteKernel(a=a, b=0.0, c=c, d=0.0)
george_kernel = k if george_kernel is None else george_kernel + k
solver = george.GP(george_kernel, solver=george.HODLRSolver)
elif args.carma:
arparams = np.random.randn(2*j)
class SHOTerm2(terms.SHOTerm):
parameter_names = ("Ampl", "Plife", "Prot")
def __init__(self, params):
Ampl, Plife, Prot = params
S0 = (Ampl*(Prot**2))/(2. * (np.pi**2) * Plife)
w0 = (2.0*np.pi)/Prot
Q = Plife(np.pi)/Prot
kernel = terms.SHOTerm(log_S0=np.log(S0), log_Q=np.log(Q), log_omega0=np.log(w0) )
return kernel
class Matern32(terms.Matern32Term):
parameter_names = ("sigma", "rho", "eps")
def __init__(self, params):
sigma, rho, eps = params
kernel = terms.Matern32Term(log_sigma=np.log(sigma), log_rho=np.log(rho), eps=eps)
return kernel
def call_gp(params):
log_sigma, log_rho, log_error_scale = params
if GP_CODE=='celerite':
kernel = terms.Matern32Term(log_sigma=log_sigma, log_rho=log_rho)
gp = celerite.GP(kernel, mean=MEAN, fit_mean=False) #log_white_noise=np.log(yerr),
gp.compute(xx, yerr=yyerr/err_norm*np.exp(log_error_scale))
return gp
elif GP_CODE=='george':
kernel = np.exp(log_sigma) * kernels.Matern32Kernel(log_rho)
gp = george.GP(kernel, mean=MEAN, fit_mean=False) #log_white_noise=np.log(yerr),
gp.compute(xx, yerr=yyerr/err_norm*np.exp(log_error_scale))
return gp
else:
raise ValueError('A bad thing happened.')
from celerite import plot_setup
plot_setup.setup()
# Set up the dimensions of the problem
N = 2**np.arange(6, 20)
times = np.empty((len(N), 3))
times[:] = np.nan
# Simulate a "dataset"
np.random.seed(42)
t = np.sort(np.random.rand(np.max(N)))
yerr = np.random.uniform(0.1, 0.2, len(t))
y = np.sin(t)
# Set up the GP model
kernel = terms.RealTerm(1.0, 0.1) + terms.ComplexTerm(0.1, 2.0, 1.6)
gp = GP(kernel)
for i, n in enumerate(N):
times[i, 0] = benchmark("gp.compute(t[:{0}], yerr[:{0}])".format(n),
"from __main__ import gp, t, yerr")
gp.compute(t[:n], yerr[:n])
times[i, 1] = benchmark("gp.log_likelihood(y[:{0}])".format(n),
"from __main__ import gp, y")
if n <= 4096:
times[i, 2] = benchmark("""
C = gp.get_matrix(t[:{0}])
C[np.diag_indices_from(C)] += yerr[:{0}]**2
cho_factor(C)
""".format(n), """