import os
import sys
import time
import math
import bisect
import numpy as np
import pandas as pd
from scipy import stats
from scipy.optimize import minimize
from scipy.optimize import curve_fit
from statsmodels.tsa import stattools
from scipy.interpolate import interp1d
from Base import Base
import lomb
class Amplitude(Base):
"""Half the difference between the maximum and the minimum magnitude"""
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
N = len(magnitude)
sorted_mag = np.sort(magnitude)
return (np.median(sorted_mag[-int(math.ceil(0.05 * N)):]) -
np.median(sorted_mag[0:int(math.ceil(0.05 * N))])) / 2.0
# return sorted_mag[10]
class Rcs(Base):
"""Range of cumulative sum"""
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
sigma = np.std(magnitude)
N = len(magnitude)
m = np.mean(magnitude)
s = np.cumsum(magnitude - m) * 1.0 / (N * sigma)
R = np.max(s) - np.min(s)
return R
class StetsonK(Base):
def __init__(self):
self.Data = ['magnitude', 'error']
def fit(self, data):
magnitude = data[0]
error = data[2]
mean_mag = (np.sum(magnitude/(error*error)) /
np.sum(1.0 / (error * error)))
N = len(magnitude)
sigmap = (np.sqrt(N * 1.0 / (N - 1)) *
(magnitude - mean_mag) / error)
K = (1 / np.sqrt(N * 1.0) *
np.sum(np.abs(sigmap)) / np.sqrt(np.sum(sigmap ** 2)))
return K
class Meanvariance(Base):
"""variability index"""
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
return np.std(magnitude) / np.mean(magnitude)
class Autocor_length(Base):
def __init__(self, lags=100):
self.Data = ['magnitude']
self.nlags = lags
def fit(self, data):
magnitude = data[0]
AC = stattools.acf(magnitude, nlags=self.nlags)
k = next((index for index, value in
enumerate(AC) if value < np.exp(-1)), None)
while k is None:
self.nlags = self.nlags + 100
AC = stattools.acf(magnitude, nlags=self.nlags)
k = next((index for index, value in
enumerate(AC) if value < np.exp(-1)), None)
return k
class SlottedA_length(Base):
def __init__(self, T=-99):
"""
lc: MACHO lightcurve in a pandas DataFrame
k: lag (default: 1)
T: tau (slot size in days. default: 4)
"""
self.Data = ['magnitude', 'time']
SlottedA_length.SAC = []
self.T = T
def slotted_autocorrelation(self, data, time, T, K,
second_round=False, K1=100):
slots = np.zeros((K, 1))
i = 1
# make time start from 0
time = time - np.min(time)
# subtract mean from mag values
m = np.mean(data)
data = data - m
prod = np.zeros((K, 1))
pairs = np.subtract.outer(time, time)
pairs[np.tril_indices_from(pairs)] = 10000000
ks = np.int64(np.floor(np.abs(pairs) / T + 0.5))
# We calculate the slotted autocorrelation for k=0 separately
idx = np.where(ks == 0)
prod[0] = ((sum(data ** 2) + sum(data[idx[0]] *
data[idx[1]])) / (len(idx[0]) + len(data)))
slots[0] = 0
# We calculate it for the rest of the ks
if second_round is False:
for k in np.arange(1, K):
idx = np.where(ks == k)
if len(idx[0]) != 0:
prod[k] = sum(data[idx[0]] * data[idx[1]]) / (len(idx[0]))
slots[i] = k
i = i + 1
else:
prod[k] = np.infty
else:
for k in np.arange(K1, K):
idx = np.where(ks == k)
if len(idx[0]) != 0:
prod[k] = sum(data[idx[0]] * data[idx[1]]) / (len(idx[0]))
slots[i - 1] = k
i = i + 1
else:
prod[k] = np.infty
np.trim_zeros(prod, trim='b')
slots = np.trim_zeros(slots, trim='b')
return prod / prod[0], np.int64(slots).flatten()
def fit(self, data):
magnitude = data[0]
time = data[1]
N = len(time)
if self.T == -99:
deltaT = time[1:] - time[:-1]
sorted_deltaT = np.sort(deltaT)
self.T = sorted_deltaT[int(N * 0.05)+1]
K = 100
[SAC, slots] = self.slotted_autocorrelation(magnitude, time, self.T, K)
# SlottedA_length.SAC = SAC
# SlottedA_length.slots = slots
SAC2 = SAC[slots]
SlottedA_length.autocor_vector = SAC2
k = next((index for index, value in
enumerate(SAC2) if value < np.exp(-1)), None)
while k is None:
K = K+K
if K > (np.max(time) - np.min(time)) / self.T:
break
else:
[SAC, slots] = self.slotted_autocorrelation(magnitude,
time, self.T, K,
second_round=True,
K1=K/2)
SAC2 = SAC[slots]
k = next((index for index, value in
enumerate(SAC2) if value < np.exp(-1)), None)
return slots[k] * self.T
def getAtt(self):
# return SlottedA_length.SAC, SlottedA_length.slots
return SlottedA_length.autocor_vector
class StetsonK_AC(SlottedA_length):
def __init__(self):
self.Data = ['magnitude', 'time', 'error']
def fit(self, data):
try:
a = StetsonK_AC()
# [autocor_vector, slots] = a.getAtt()
autocor_vector = a.getAtt()
# autocor_vector = autocor_vector[slots]
N_autocor = len(autocor_vector)
sigmap = (np.sqrt(N_autocor * 1.0 / (N_autocor - 1)) *
(autocor_vector - np.mean(autocor_vector)) /
np.std(autocor_vector))
K = (1 / np.sqrt(N_autocor * 1.0) *
np.sum(np.abs(sigmap)) / np.sqrt(np.sum(sigmap ** 2)))
return K
except:
print "error: please run SlottedA_length first to generate values for StetsonK_AC "
class StetsonL(Base):
def __init__(self):
self.Data = ['magnitude', 'time', 'error', 'magnitude2', 'error2']
def fit(self, data):
aligned_magnitude = data[4]
aligned_magnitude2 = data[5]
aligned_error = data[7]
aligned_error2 = data[8]
N = len(aligned_magnitude)
mean_mag = (np.sum(aligned_magnitude/(aligned_error*aligned_error)) /
np.sum(1.0 / (aligned_error * aligned_error)))
mean_mag2 = (np.sum(aligned_magnitude2/(aligned_error2*aligned_error2)) /
np.sum(1.0 / (aligned_error2 * aligned_error2)))
sigmap = (np.sqrt(N * 1.0 / (N - 1)) *
(aligned_magnitude[:N] - mean_mag) /
aligned_error)
sigmaq = (np.sqrt(N * 1.0 / (N - 1)) *
(aligned_magnitude2[:N] - mean_mag2) /
aligned_error2)
sigma_i = sigmap * sigmaq
J = (1.0 / len(sigma_i) *
np.sum(np.sign(sigma_i) * np.sqrt(np.abs(sigma_i))))
K = (1 / np.sqrt(N * 1.0) *
np.sum(np.abs(sigma_i)) / np.sqrt(np.sum(sigma_i ** 2)))
return J * K / 0.798
class Con(Base):
"""Index introduced for selection of variable starts from OGLE database.
To calculate Con, we counted the number of three consecutive measurements
that are out of 2sigma range, and normalized by N-2
Pavlos not happy
"""
def __init__(self, consecutiveStar=3):
self.Data = ['magnitude']
self.consecutiveStar = consecutiveStar
def fit(self, data):
magnitude = data[0]
N = len(magnitude)
if N < self.consecutiveStar:
return 0
sigma = np.std(magnitude)
m = np.mean(magnitude)
count = 0
for i in xrange(N - self.consecutiveStar + 1):
flag = 0
for j in xrange(self.consecutiveStar):
if(magnitude[i + j] > m + 2 * sigma or magnitude[i + j] < m - 2 * sigma):
flag = 1
else:
flag = 0
break
if flag:
count = count + 1
return count * 1.0 / (N - self.consecutiveStar + 1)
# class VariabilityIndex(Base):
# # Eta. Removed, it is not invariant to time sampling
# '''
# The index is the ratio of mean of the square of successive difference to
# the variance of data points
# '''
# def __init__(self):
# self.category='timeSeries'
# def fit(self, data):
# N = len(data)
# sigma2 = np.var(data)
# return 1.0/((N-1)*sigma2) * np.sum(np.power(data[1:] - data[:-1] , 2)
# )
class Color(Base):
"""Average color for each MACHO lightcurve
mean(B1) - mean(B2)
"""
def __init__(self):
self.Data = ['magnitude', 'time', 'magnitude2']
def fit(self, data):
magnitude = data[0]
magnitude2 = data[3]
return np.mean(magnitude) - np.mean(magnitude2)
# The categories of the following featurs should be revised
class Beyond1Std(Base):
"""Percentage of points beyond one st. dev. from the weighted
(by photometric errors) mean
"""
def __init__(self):
self.Data = ['magnitude', 'error']
def fit(self, data):
magnitude = data[0]
error = data[2]
n = len(magnitude)
weighted_mean = np.average(magnitude, weights=1 / error ** 2)
# Standard deviation with respect to the weighted mean
var = sum((magnitude - weighted_mean) ** 2)
std = np.sqrt((1.0 / (n - 1)) * var)
count = np.sum(np.logical_or(magnitude > weighted_mean + std,
magnitude < weighted_mean - std))
return float(count) / n
class SmallKurtosis(Base):
"""Small sample kurtosis of the magnitudes.
See http://www.xycoon.com/peakedness_small_sample_test_1.htm
"""
def __init__(self):
self.category = 'basic'
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
n = len(magnitude)
mean = np.mean(magnitude)
std = np.std(magnitude)
S = sum(((magnitude - mean) / std) ** 4)
c1 = float(n * (n + 1)) / ((n - 1) * (n - 2) * (n - 3))
c2 = float(3 * (n - 1) ** 2) / ((n - 2) * (n - 3))
return c1 * S - c2
class Std(Base):
"""Standard deviation of the magnitudes"""
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
return np.std(magnitude)
class Skew(Base):
"""Skewness of the magnitudes"""
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
return stats.skew(magnitude)
class StetsonJ(Base):
"""Stetson (1996) variability index, a robust standard deviation"""
def __init__(self):
self.Data = ['magnitude', 'time', 'error', 'magnitude2', 'error2']
def fit(self, data):
aligned_magnitude = data[4]
aligned_magnitude2 = data[5]
aligned_error = data[7]
aligned_error2 = data[8]
N = len(aligned_magnitude)
mean_mag = (np.sum(aligned_magnitude/(aligned_error*aligned_error)) /
np.sum(1.0 / (aligned_error * aligned_error)))
mean_mag2 = (np.sum(aligned_magnitude2 / (aligned_error2*aligned_error2)) /
np.sum(1.0 / (aligned_error2 * aligned_error2)))
sigmap = (np.sqrt(N * 1.0 / (N - 1)) *
(aligned_magnitude[:N] - mean_mag) /
aligned_error)
sigmaq = (np.sqrt(N * 1.0 / (N - 1)) *
(aligned_magnitude2[:N] - mean_mag2) /
aligned_error2)
sigma_i = sigmap * sigmaq
J = (1.0 / len(sigma_i) * np.sum(np.sign(sigma_i) *
np.sqrt(np.abs(sigma_i))))
return J
class MaxSlope(Base):
"""
Examining successive (time-sorted) magnitudes, the maximal first difference
(value of delta magnitude over delta time)
"""
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
magnitude = data[0]
time = data[1]
slope = np.abs(magnitude[1:] - magnitude[:-1]) / (time[1:] - time[:-1])
np.max(slope)
return np.max(slope)
class MedianAbsDev(Base):
def __init__(self):
self.category = 'basic'
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
median = np.median(magnitude)
devs = (abs(magnitude - median))
return np.median(devs)
class MedianBRP(Base):
"""Median buffer range percentage
Fraction (<= 1) of photometric points within amplitude/10
of the median magnitude
"""
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
median = np.median(magnitude)
amplitude = (np.max(magnitude) - np.min(magnitude)) / 10
n = len(magnitude)
count = np.sum(np.logical_and(magnitude < median + amplitude,
magnitude > median - amplitude))
return float(count) / n
class PairSlopeTrend(Base):
"""
Considering the last 30 (time-sorted) measurements of source magnitude,
the fraction of increasing first differences minus the fraction of
decreasing first differences.
"""
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
data_last = magnitude[-30:]
return (float(len(np.where(np.diff(data_last) > 0)[0]) -
len(np.where(np.diff(data_last) <= 0)[0])) / 30)
class FluxPercentileRatioMid20(Base):
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
sorted_data = np.sort(magnitude)
lc_length = len(sorted_data)
F_60_index = int(math.ceil(0.60 * lc_length))
F_40_index = int(math.ceil(0.40 * lc_length))
F_5_index = int(math.ceil(0.05 * lc_length))
F_95_index = int(math.ceil(0.95 * lc_length))
F_40_60 = sorted_data[F_60_index] - sorted_data[F_40_index]
F_5_95 = sorted_data[F_95_index] - sorted_data[F_5_index]
F_mid20 = F_40_60 / F_5_95
return F_mid20
class FluxPercentileRatioMid35(Base):
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
sorted_data = np.sort(magnitude)
lc_length = len(sorted_data)
F_325_index = int(math.ceil(0.325 * lc_length))
F_675_index = int(math.ceil(0.675 * lc_length))
F_5_index = int(math.ceil(0.05 * lc_length))
F_95_index = int(math.ceil(0.95 * lc_length))
F_325_675 = sorted_data[F_675_index] - sorted_data[F_325_index]
F_5_95 = sorted_data[F_95_index] - sorted_data[F_5_index]
F_mid35 = F_325_675 / F_5_95
return F_mid35
class FluxPercentileRatioMid50(Base):
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
sorted_data = np.sort(magnitude)
lc_length = len(sorted_data)
F_25_index = int(math.ceil(0.25 * lc_length))
F_75_index = int(math.ceil(0.75 * lc_length))
F_5_index = int(math.ceil(0.05 * lc_length))
F_95_index = int(math.ceil(0.95 * lc_length))
F_25_75 = sorted_data[F_75_index] - sorted_data[F_25_index]
F_5_95 = sorted_data[F_95_index] - sorted_data[F_5_index]
F_mid50 = F_25_75 / F_5_95
return F_mid50
class FluxPercentileRatioMid65(Base):
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
sorted_data = np.sort(magnitude)
lc_length = len(sorted_data)
F_175_index = int(math.ceil(0.175 * lc_length))
F_825_index = int(math.ceil(0.825 * lc_length))
F_5_index = int(math.ceil(0.05 * lc_length))
F_95_index = int(math.ceil(0.95 * lc_length))
F_175_825 = sorted_data[F_825_index] - sorted_data[F_175_index]
F_5_95 = sorted_data[F_95_index] - sorted_data[F_5_index]
F_mid65 = F_175_825 / F_5_95
return F_mid65
class FluxPercentileRatioMid80(Base):
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
sorted_data = np.sort(magnitude)
lc_length = len(sorted_data)
F_10_index = int(math.ceil(0.10 * lc_length))
F_90_index = int(math.ceil(0.90 * lc_length))
F_5_index = int(math.ceil(0.05 * lc_length))
F_95_index = int(math.ceil(0.95 * lc_length))
F_10_90 = sorted_data[F_90_index] - sorted_data[F_10_index]
F_5_95 = sorted_data[F_95_index] - sorted_data[F_5_index]
F_mid80 = F_10_90 / F_5_95
return F_mid80
class PercentDifferenceFluxPercentile(Base):
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
median_data = np.median(magnitude)
sorted_data = np.sort(magnitude)
lc_length = len(sorted_data)
F_5_index = int(math.ceil(0.05 * lc_length))
F_95_index = int(math.ceil(0.95 * lc_length))
F_5_95 = sorted_data[F_95_index] - sorted_data[F_5_index]
percent_difference = F_5_95 / median_data
return percent_difference
class PercentAmplitude(Base):
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
median_data = np.median(magnitude)
distance_median = np.abs(magnitude - median_data)
max_distance = np.max(distance_median)
percent_amplitude = max_distance / median_data
return percent_amplitude
class LinearTrend(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
magnitude = data[0]
time = data[1]
regression_slope = stats.linregress(time, magnitude)[0]
return regression_slope
class Eta_color(Base):
def __init__(self):
self.Data = ['magnitude', 'time', 'magnitude2']
def fit(self, data):
aligned_magnitude = data[4]
aligned_magnitude2 = data[5]
aligned_time = data[6]
N = len(aligned_magnitude)
B_Rdata = aligned_magnitude - aligned_magnitude2
w = 1.0 / np.power(aligned_time[1:] - aligned_time[:-1], 2)
w_mean = np.mean(w)
N = len(aligned_time)
sigma2 = np.var(B_Rdata)
S1 = sum(w * (B_Rdata[1:] - B_Rdata[:-1]) ** 2)
S2 = sum(w)
eta_B_R = (w_mean * np.power(aligned_time[N - 1] -
aligned_time[0], 2) * S1 / (sigma2 * S2 * N ** 2))
return eta_B_R
class Eta_e(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
magnitude = data[0]
time = data[1]
w = 1.0 / np.power(np.subtract(time[1:], time[:-1]), 2)
w_mean = np.mean(w)
N = len(time)
sigma2 = np.var(magnitude)
S1 = sum(w * (magnitude[1:] - magnitude[:-1]) ** 2)
S2 = sum(w)
eta_e = (w_mean * np.power(time[N - 1] -
time[0], 2) * S1 / (sigma2 * S2 * N ** 2))
return eta_e
class Mean(Base):
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
B_mean = np.mean(magnitude)
return B_mean
class Q31(Base):
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
return np.percentile(magnitude, 75) - np.percentile(magnitude, 25)
class Q31_color(Base):
def __init__(self):
self.Data = ['magnitude', 'time', 'magnitude2']
def fit(self, data):
aligned_magnitude = data[4]
aligned_magnitude2 = data[5]
N = len(aligned_magnitude)
b_r = aligned_magnitude[:N] - aligned_magnitude2[:N]
return np.percentile(b_r, 75) - np.percentile(b_r, 25)
class AndersonDarling(Base):
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = data[0]
ander = stats.anderson(magnitude)[0]
return 1 / (1.0 + np.exp(-10 * (ander - 0.3)))
class PeriodLS(Base):
def __init__(self, ofac=6.):
self.Data = ['magnitude', 'time']
self.ofac = ofac
def fit(self, data):
magnitude = data[0]
time = data[1]
global new_time
global prob
global period
fx, fy, nout, jmax, prob = lomb.fasper(time, magnitude, self.ofac, 100.)
period = fx[jmax]
T = 1.0 / period
new_time = np.mod(time, 2 * T) / (2 * T)
return T
class Period_fit(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return prob
except:
print "error: please run PeriodLS first to generate values for Period_fit"
class Psi_CS(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
magnitude = data[0]
time = data[1]
folded_data = magnitude[np.argsort(new_time)]
sigma = np.std(folded_data)
N = len(folded_data)
m = np.mean(folded_data)
s = np.cumsum(folded_data - m) * 1.0 / (N * sigma)
R = np.max(s) - np.min(s)
return R
except:
print "error: please run PeriodLS first to generate values for Psi_CS"
class Psi_eta(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
# folded_time = np.sort(new_time)
try:
magnitude = data[0]
folded_data = magnitude[np.argsort(new_time)]
# w = 1.0 / np.power(folded_time[1:]-folded_time[:-1] ,2)
# w_mean = np.mean(w)
# N = len(folded_time)
# sigma2=np.var(folded_data)
# S1 = sum(w*(folded_data[1:]-folded_data[:-1])**2)
# S2 = sum(w)
# Psi_eta = w_mean * np.power(folded_time[N-1]-folded_time[0],2) * S1 /
# (sigma2 * S2 * N**2)
N = len(folded_data)
sigma2 = np.var(folded_data)
Psi_eta = (1.0 / ((N - 1) * sigma2) *
np.sum(np.power(folded_data[1:] - folded_data[:-1], 2)))
return Psi_eta
except:
print "error: please run PeriodLS first to generate values for Psi_eta"
class CAR_sigma(Base):
def __init__(self):
self.Data = ['magnitude', 'time', 'error']
def CAR_Lik(self, parameters, t, x, error_vars):
sigma = parameters[0]
tau = parameters[1]
# b = parameters[1] #comment it to do 2 pars estimation
# tau = params(1,1);
# sigma = sqrt(2*var(x)/tau);
b = np.mean(x) / tau
epsilon = 1e-300
cte_neg = -np.infty
num_datos = np.size(x)
Omega = []
x_hat = []
a = []
x_ast = []
# Omega = np.zeros((num_datos,1))
# x_hat = np.zeros((num_datos,1))
# a = np.zeros((num_datos,1))
# x_ast = np.zeros((num_datos,1))
# Omega[0]=(tau*(sigma**2))/2.
# x_hat[0]=0.
# a[0]=0.
# x_ast[0]=x[0] - b*tau
Omega.append((tau * (sigma ** 2)) / 2.)
x_hat.append(0.)
a.append(0.)
x_ast.append(x[0] - b * tau)
loglik = 0.
for i in range(1, num_datos):
a_new = np.exp(-(t[i] - t[i - 1]) / tau)
x_ast.append(x[i] - b * tau)
x_hat.append(
a_new * x_hat[i - 1] +
(a_new * Omega[i - 1] / (Omega[i - 1] + error_vars[i - 1])) *
(x_ast[i - 1] - x_hat[i - 1]))
Omega.append(
Omega[0] * (1 - (a_new ** 2)) + ((a_new ** 2)) * Omega[i - 1] *
(1 - (Omega[i - 1] / (Omega[i - 1] + error_vars[i - 1]))))
# x_ast[i]=x[i] - b*tau
# x_hat[i]=a_new*x_hat[i-1] + (a_new*Omega[i-1]/(Omega[i-1] +
# error_vars[i-1]))*(x_ast[i-1]-x_hat[i-1])
# Omega[i]=Omega[0]*(1-(a_new**2)) + ((a_new**2))*Omega[i-1]*
# ( 1 - (Omega[i-1]/(Omega[i-1]+ error_vars[i-1])))
loglik_inter = np.log(
((2 * np.pi * (Omega[i] + error_vars[i])) ** -0.5) *
(np.exp(-0.5 * (((x_hat[i] - x_ast[i]) ** 2) /
(Omega[i] + error_vars[i]))) + epsilon))
loglik = loglik + loglik_inter
if(loglik <= cte_neg):
print('CAR lik se fue a inf')
return None
# the minus one is to perfor maximization using the minimize function
return -loglik
def calculateCAR(self, time, data, error):
x0 = [10, 0.5]
bnds = ((0, 100), (0, 100))
# res = minimize(self.CAR_Lik, x0, args=(LC[:,0],LC[:,1],LC[:,2]) ,
# method='nelder-mead',bounds = bnds)
res = minimize(self.CAR_Lik, x0, args=(time, data, error),
method='powell', bounds=bnds)
# options={'disp': True}
sigma = res.x[0]
global tau
tau = res.x[1]
return sigma
# def getAtt(self):
# return CAR_sigma.tau
def fit(self, data):
# LC = np.hstack((self.time , data.reshape((self.N,1)), self.error))
N = len(data[0])
magnitude = data[0].reshape((N, 1))
time = data[1].reshape((N, 1))
error = data[2].reshape((N, 1)) ** 2
a = self.calculateCAR(time, magnitude, error)
return a
class CAR_tau(Base):
def __init__(self):
self.Data = ['magnitude', 'time', 'error']
def fit(self, data):
try:
return tau
except:
print "error: please run CAR_sigma first to generate values for CAR_tau"
class CAR_mean(Base):
def __init__(self):
self.Data = ['magnitude', 'time', 'error']
def fit(self, data):
magnitude = data[0]
try:
return np.mean(magnitude) / tau
except:
print "error: please run CAR_sigma first to generate values for CAR_mean"
class Freq1_harmonics_amplitude_0(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
magnitude = data[0]
time = data[1]
time = time - np.min(time)
global A
global PH
global scaledPH
A = []
PH = []
scaledPH = []
def model(x, a, b, c, Freq):
return a*np.sin(2*np.pi*Freq*x)+b*np.cos(2*np.pi*Freq*x)+c
for i in range(3):
wk1, wk2, nout, jmax, prob = lomb.fasper(time, magnitude, 6., 100.)
fundamental_Freq = wk1[jmax]
# fit to a_i sin(2pi f_i t) + b_i cos(2 pi f_i t) + b_i,o
# a, b are the parameters we care about
# c is a constant offset
# f is the fundamental Frequency
def yfunc(Freq):
def func(x, a, b, c):
return a*np.sin(2*np.pi*Freq*x)+b*np.cos(2*np.pi*Freq*x)+c
return func
Atemp = []
PHtemp = []
popts = []
for j in range(4):
popt, pcov = curve_fit(yfunc((j+1)*fundamental_Freq), time, magnitude)
Atemp.append(np.sqrt(popt[0]**2+popt[1]**2))
PHtemp.append(np.arctan(popt[1] / popt[0]))
popts.append(popt)
A.append(Atemp)
PH.append(PHtemp)
for j in range(4):
magnitude = np.array(magnitude) - model(time, popts[j][0], popts[j][1], popts[j][2], (j+1)*fundamental_Freq)
for ph in PH:
scaledPH.append(np.array(ph) - ph[0])
return A[0][0]
class Freq1_harmonics_amplitude_1(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return A[0][1]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq1_harmonics_amplitude_2(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return A[0][2]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq1_harmonics_amplitude_3(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return A[0][3]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq2_harmonics_amplitude_0(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return A[1][0]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq2_harmonics_amplitude_1(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return A[1][1]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq2_harmonics_amplitude_2(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return A[1][2]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq2_harmonics_amplitude_3(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return A[1][3]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq3_harmonics_amplitude_0(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return A[2][0]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq3_harmonics_amplitude_1(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return A[2][1]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq3_harmonics_amplitude_2(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return A[2][2]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq3_harmonics_amplitude_3(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return A[2][3]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq1_harmonics_rel_phase_0(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return scaledPH[0][0]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq1_harmonics_rel_phase_1(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return scaledPH[0][1]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq1_harmonics_rel_phase_2(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return scaledPH[0][2]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq1_harmonics_rel_phase_3(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return scaledPH[0][3]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq2_harmonics_rel_phase_0(Base):
def __init__(self):
self.category = 'timeSeries'
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return scaledPH[1][0]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq2_harmonics_rel_phase_1(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return scaledPH[1][1]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq2_harmonics_rel_phase_2(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return scaledPH[1][2]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq2_harmonics_rel_phase_3(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return scaledPH[1][3]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq3_harmonics_rel_phase_0(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return scaledPH[2][0]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq3_harmonics_rel_phase_1(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return scaledPH[2][1]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq3_harmonics_rel_phase_2(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return scaledPH[2][2]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Freq3_harmonics_rel_phase_3(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return scaledPH[2][3]
except:
print "error: please run Freq1_harmonics_amplitude_0 first to generate values for all harmonics"
class Gskew(Base):
"""Median-based measure of the skew"""
def __init__(self):
self.Data = ['magnitude']
def fit(self, data):
magnitude = np.array(data[0])
median_mag = np.median(magnitude)
F_3_value = np.percentile(magnitude, 3)
F_97_value = np.percentile(magnitude, 97)
return (np.median(magnitude[magnitude <= F_3_value]) +
np.median(magnitude[magnitude >= F_97_value])
- 2*median_mag)
class StructureFunction_index_21(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
magnitude = data[0]
time = data[1]
global m_21
global m_31
global m_32
Nsf = 100
Np = 100
sf1 = np.zeros(Nsf)
sf2 = np.zeros(Nsf)
sf3 = np.zeros(Nsf)
f = interp1d(time, magnitude)
time_int = np.linspace(np.min(time), np.max(time), Np)
mag_int = f(time_int)
for tau in np.arange(1, Nsf):
sf1[tau-1] = np.mean(np.power(np.abs(mag_int[0:Np-tau] - mag_int[tau:Np]) , 1.0))
sf2[tau-1] = np.mean(np.abs(np.power(np.abs(mag_int[0:Np-tau] - mag_int[tau:Np]) , 2.0)))
sf3[tau-1] = np.mean(np.abs(np.power(np.abs(mag_int[0:Np-tau] - mag_int[tau:Np]) , 3.0)))
sf1_log = np.log10(np.trim_zeros(sf1))
sf2_log = np.log10(np.trim_zeros(sf2))
sf3_log = np.log10(np.trim_zeros(sf3))
m_21, b_21 = np.polyfit(sf1_log, sf2_log, 1)
m_31, b_31 = np.polyfit(sf1_log, sf3_log, 1)
m_32, b_32 = np.polyfit(sf2_log, sf3_log, 1)
return m_21
class StructureFunction_index_31(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return m_31
except:
print "error: please run StructureFunction_index_21 first to generate values for all Structure Function"
class StructureFunction_index_32(Base):
def __init__(self):
self.Data = ['magnitude', 'time']
def fit(self, data):
try:
return m_32
except:
print "error: please run StructureFunction_index_21 first to generate values for all Structure Function"
