How to use the celerite.terms function in celerite

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), """