<commit_before>#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import division, print_function

import emcee
import fitsio
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

from plot_setup import setup, SQUARE_FIGSIZE, COLORS

import genrp
from genrp import kernels

np.random.seed(42)
setup()

# Define the custom kernel
class RotationKernel(kernels.Kernel):
    parameter_names = ("log_amp", "log_timescale", "log_period", "log_factor")

    @property
    def p_real(self):
        return 1

    @property
    def p_complex(self):
        return 1

    @property
    def alpha_real(self):
        f = np.exp(self.log_factor)
        return np.array([np.exp(self.log_amp) * (1.0 + f) / (2.0 + f)])

    @property
    def beta_real(self):
        return np.array([np.exp(-self.log_timescale)])

    @property
    def alpha_complex_real(self):
        f = np.exp(self.log_factor)
        return np.array([np.exp(self.log_amp) / (2.0 + f)])

    @property
    def alpha_complex_imag(self):
        return np.array([0.0])

    @property
    def beta_complex_real(self):
        return np.array([np.exp(-self.log_timescale)])

    @property
    def beta_complex_imag(self):
        return np.array([2*np.pi*np.exp(-self.log_period)])

# Load the data
data = fitsio.read("data/kplr001430163-2013011073258_llc.fits")
# data = fitsio.read("data/kplr001430163-2010355172524_llc.fits")
m = data["SAP_QUALITY"] == 0
m &= np.isfinite(data["TIME"]) & np.isfinite(data["PDCSAP_FLUX"])
t = np.ascontiguousarray(data["TIME"][m], dtype=np.float64)
t -= np.min(t)
y = np.ascontiguousarray(data["PDCSAP_FLUX"][m], dtype=np.float64)
yerr = np.ascontiguousarray(data["PDCSAP_FLUX_ERR"][m], dtype=np.float64)

# Normalize the data
mean = np.median(y)
y = (y / mean - 1.0) * 1e3
yerr *= 1e3 / mean

# Set up the GP model
kernel = RotationKernel(
    np.log(np.var(y)), np.log(10), np.log(2.0), np.log(1.0)
)
gp = genrp.GP(kernel, mean=np.median(y))
gp.compute(t, yerr)

# Define the model
def neg_log_like(params, y, gp):
    gp.set_parameter_vector(params)
    return -gp.log_likelihood(y)

# Optimize with random restarts
initial_params = gp.get_parameter_vector()
bounds = [
    np.log(np.var(y) * np.array([0.01, 100])),
    np.log([np.max(np.diff(t)), (t.max() - t.min())]),
    np.log([3*np.median(np.diff(t)), 0.5*(t.max() - t.min())]),
    [-5.0, np.log(5.0)],
]
best = (np.inf, initial_params)
for i in range(10):<commit_msg>updating rotation example to new terms api [ci skip]<commit_after>#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import division, print_function

import emcee
import fitsio
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

from plot_setup import setup, SQUARE_FIGSIZE, COLORS

import genrp
from genrp import terms

np.random.seed(42)
setup()

# Define the custom kernel
class RotationTerm(terms.Term):
    parameter_names = ("log_amp", "log_timescale", "log_period", "log_factor")

    def get_real_coefficients(self):
        f = np.exp(self.log_factor)
        return (
            np.exp(self.log_amp) * (1.0 + f) / (2.0 + f),
            np.exp(-self.log_timescale),
        )

    def get_complex_coefficients(self):
        f = np.exp(self.log_factor)
        return (
            np.exp(self.log_amp) / (2.0 + f),
            0.0,
            np.exp(-self.log_timescale),
            2*np.pi*np.exp(-self.log_period),
        )

# Load the data
data = fitsio.read("data/kplr001430163-2013011073258_llc.fits")
m = data["SAP_QUALITY"] == 0
m &= np.isfinite(data["TIME"]) & np.isfinite(data["PDCSAP_FLUX"])
t = np.ascontiguousarray(data["TIME"][m], dtype=np.float64)
t -= np.min(t)
y = np.ascontiguousarray(data["PDCSAP_FLUX"][m], dtype=np.float64)
yerr = np.ascontiguousarray(data["PDCSAP_FLUX_ERR"][m], dtype=np.float64)

# Normalize the data
mean = np.median(y)
y = (y / mean - 1.0) * 1e3
yerr *= 1e3 / mean

# Set up the GP model
kernel = RotationTerm(
    np.log(np.var(y)), np.log(10), np.log(2.0), np.log(1.0)
)
gp = genrp.GP(kernel, mean=np.median(y))
gp.compute(t, yerr)

# Define the model
def neg_log_like(params, y, gp):
    gp.set_parameter_vector(params)
    return -gp.log_likelihood(y)

# Optimize with random restarts
initial_params = gp.get_parameter_vector()
bounds = [
    np.log(np.var(y) * np.array([0.01, 100])),
    np.log([np.max(np.diff(t)), (t.max() - t.min())]),
    np.log([3*np.median(np.diff(t)), 0.5*(t.max() - t.min())]),
    [-5.0, np.log(5.0)],
]
best = (np.inf, initial_params)
for i in range(10):
{'fraction_non_alphanumeric': 0.1646180860403863, 'fraction_numerical': 0.03511852502194908, 'mean_word_length': 3.507418397626113, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3363342', 'n_tokens_mistral': 1882, 'n_tokens_neox': 1770, 'n_words': 353}