From aa81ce7f68d137798126b6e425dd069c3b4fde3a Mon Sep 17 00:00:00 2001 From: Maurizio D'Addona Date: Tue, 28 Mar 2023 13:21:14 +0200 Subject: [PATCH] Fixed for python3 syntax --- FATS/Base.py | 7 +- FATS/Feature.py | 64 ++++-- FATS/FeatureFunctionLib.py | 170 ++++++++------ FATS/PreprocessLC.py | 54 +++-- FATS/alignLC.py | 9 +- FATS/featureFunction.py | 6 +- FATS/import_lc_cluster.py | 16 +- FATS/import_lightcurve.py | 16 +- FATS/lomb.py | 414 +++++++++++++++++----------------- FATS/test_library.py | 450 ++++++++++++++++++++++--------------- requirements.txt | 6 +- 11 files changed, 678 insertions(+), 534 deletions(-) diff --git a/FATS/Base.py b/FATS/Base.py index ea71912..555a9db 100644 --- a/FATS/Base.py +++ b/FATS/Base.py @@ -1,7 +1,10 @@ -import os,sys,time +import os, sys, time import numpy as np + + class Base: def __init__(self): - self.category='all' + self.category = 'all' + def fit(self, data): return self diff --git a/FATS/Feature.py b/FATS/Feature.py index 3fc1817..ac42c5b 100644 --- a/FATS/Feature.py +++ b/FATS/Feature.py @@ -7,7 +7,7 @@ import pandas as pd import matplotlib.pyplot as plt -import featureFunction +from . import featureFunction class FeatureSpace: @@ -53,21 +53,25 @@ def __init__(self, Data=None, featureList=None, excludeList=[], **kwargs): self.Data = Data if self.Data == 'all': - if featureList == None: + if featureList is None: - if excludeList == None: + if excludeList is None: for name, obj in inspect.getmembers(featureFunction): if inspect.isclass(obj) and name != 'Base': if obj.__module__.endswith('FeatureFunctionLib'): # if set(obj().Data).issubset(self.Data): - self.featureOrder.append((inspect.getsourcelines(obj)[-1:])[0]) + self.featureOrder.append( + (inspect.getsourcelines(obj)[-1:])[0] + ) self.featureList.append(name) else: for name, obj in inspect.getmembers(featureFunction): if inspect.isclass(obj) and name != 'Base' and not name in excludeList: if obj.__module__.endswith('FeatureFunctionLib'): # if set(obj().Data).issubset(self.Data): - self.featureOrder.append((inspect.getsourcelines(obj)[-1:])[0]) + self.featureOrder.append( + (inspect.getsourcelines(obj)[-1:])[0] + ) self.featureList.append(name) else: @@ -85,15 +89,23 @@ def __init__(self, Data=None, featureList=None, excludeList=[], **kwargs): if obj.__module__.endswith('FeatureFunctionLib'): if name in kwargs.keys(): if set(obj(kwargs[name]).Data).issubset(self.Data): - self.featureOrder.append((inspect.getsourcelines(obj)[-1:])[0]) + self.featureOrder.append( + (inspect.getsourcelines(obj)[-1:])[0] + ) self.featureList.append(name) else: if set(obj().Data).issubset(self.Data): - self.featureOrder.append((inspect.getsourcelines(obj)[-1:])[0]) + self.featureOrder.append( + (inspect.getsourcelines(obj)[-1:])[0] + ) self.featureList.append(name) else: - print "Warning: the feature", name, "could not be calculated because", obj().Data, "are needed." + print( + "Warning: the feature", name, + "could not be calculated because", + obj().Data, "are needed." + ) else: for feature in featureList: @@ -103,12 +115,19 @@ def __init__(self, Data=None, featureList=None, excludeList=[], **kwargs): if set(obj().Data).issubset(self.Data): self.featureList.append(name) else: - print "Warning: the feature", name, "could not be calculated because", obj().Data, "are needed." + print( + "Warning: the feature", name, + "could not be calculated because", + obj().Data, "are needed." + ) if self.featureOrder != []: self.sort = True self.featureOrder = np.argsort(self.featureOrder) - self.featureList = [self.featureList[i] for i in self.featureOrder] + self.featureList = [ + self.featureList[i] + for i in self.featureOrder + ] self.idx = np.argsort(self.featureList) else: @@ -120,20 +139,20 @@ def __init__(self, Data=None, featureList=None, excludeList=[], **kwargs): if item in kwargs.keys(): try: a = getattr(m, item)(kwargs[item]) - except: - print "error in feature " + item + except Exception: + print("error in feature", item) sys.exit(1) else: try: a = getattr(m, item)() - except: - print " could not find feature " + item + except Exception: + print(" could not find feature", item) # discuss -- should we exit? sys.exit(1) try: self.featureFunc.append(a.fit) - except: - print "could not initilize " + item + except Exception: + print("could not initilize", item) def calculateFeature(self, data): self._X = np.asarray(data) @@ -144,17 +163,22 @@ def calculateFeature(self, data): def result(self, method='array'): if method == 'array': - if self.sort == True: + if self.sort is True: return [np.asarray(self.__result)[i] for i in self.idx] else: return np.asarray(self.__result) elif method == 'dict': - if self.sort == True: - return dict(zip([self.featureList[i] for i in self.idx], [np.asarray(self.__result)[i] for i in self.idx])) + if self.sort is True: + return dict( + zip( + [self.featureList[i] for i in self.idx], + [np.asarray(self.__result)[i] for i in self.idx] + ) + ) else: return dict(zip(self.featureList, np.asarray(self.__result))) elif method == 'features': - if self.sort == True: + if self.sort is True: return [self.featureList[i] for i in self.idx] else: return self.featureList diff --git a/FATS/FeatureFunctionLib.py b/FATS/FeatureFunctionLib.py index 8ae0ff3..64a3efe 100644 --- a/FATS/FeatureFunctionLib.py +++ b/FATS/FeatureFunctionLib.py @@ -12,8 +12,13 @@ from statsmodels.tsa import stattools from scipy.interpolate import interp1d -from Base import Base -import lomb +from .Base import Base +from . import lomb + +ERR_MSG_HARM = "error: please run Freq1_harmonics_amplitude_0 first to " \ + "generate values for all harmonics" + +ERR_MSG_RUN = "error: please run {} first to generate values for {}" class Amplitude(Base): @@ -232,9 +237,8 @@ def fit(self, data): return K - except: - - print "error: please run SlottedA_length first to generate values for StetsonK_AC " + except Exception: + print(ERR_MSG_RUN.format('SlottedA_length', 'StetsonK_AC')) class StetsonL(Base): @@ -296,9 +300,9 @@ def fit(self, data): m = np.mean(magnitude) count = 0 - for i in xrange(N - self.consecutiveStar + 1): + for i in range(N - self.consecutiveStar + 1): flag = 0 - for j in xrange(self.consecutiveStar): + for j in range(self.consecutiveStar): if(magnitude[i + j] > m + 2 * sigma or magnitude[i + j] < m - 2 * sigma): flag = 1 else: @@ -804,7 +808,12 @@ def fit(self, data): global prob global period - fx, fy, nout, jmax, prob = lomb.fasper(time, magnitude, self.ofac, 100.) + 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) @@ -822,8 +831,8 @@ def fit(self, data): try: return prob - except: - print "error: please run PeriodLS first to generate values for Period_fit" + except Exception: + print(ERR_MSG_RUN.format('PeriodLS', 'Period_fit')) class Psi_CS(Base): @@ -845,8 +854,8 @@ def fit(self, data): R = np.max(s) - np.min(s) return R - except: - print "error: please run PeriodLS first to generate values for Psi_CS" + except Exception: + print(ERR_MSG_RUN.format('PeriodLS', 'Psi_CS')) class Psi_eta(Base): @@ -881,8 +890,8 @@ def fit(self, data): 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" + except Exception: + print(ERR_MSG_RUN.format('PeriodLS', 'Psi_eta')) class CAR_sigma(Base): @@ -952,7 +961,7 @@ def CAR_Lik(self, parameters, t, x, error_vars): loglik = loglik + loglik_inter - if(loglik <= cte_neg): + if loglik <= cte_neg: print('CAR lik se fue a inf') return None @@ -1000,8 +1009,8 @@ def fit(self, data): try: return tau - except: - print "error: please run CAR_sigma first to generate values for CAR_tau" + except Exception: + print(ERR_MSG_RUN.format('CAR_sigma', 'CAR_tau')) class CAR_mean(Base): @@ -1016,8 +1025,8 @@ def fit(self, data): try: return np.mean(magnitude) / tau - except: - print "error: please run CAR_sigma first to generate values for CAR_mean" + except Exception: + print(ERR_MSG_RUN.format('CAR_sigma', 'CAR_mean')) class Freq1_harmonics_amplitude_0(Base): @@ -1061,7 +1070,11 @@ def func(x, a, b, c): popts = [] for j in range(4): - popt, pcov = curve_fit(yfunc((j+1)*fundamental_Freq), time, magnitude) + 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) @@ -1070,7 +1083,14 @@ def func(x, a, b, c): 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) + magnitude = np.array(magnitude) + 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]) @@ -1087,8 +1107,8 @@ 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" + except Exception: + print(ERR_MSG_HARM) class Freq1_harmonics_amplitude_2(Base): @@ -1100,8 +1120,8 @@ 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" + except Exception: + print(ERR_MSG_HARM) class Freq1_harmonics_amplitude_3(Base): @@ -1112,8 +1132,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq2_harmonics_amplitude_0(Base): @@ -1124,8 +1144,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq2_harmonics_amplitude_1(Base): @@ -1136,8 +1156,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq2_harmonics_amplitude_2(Base): @@ -1148,8 +1168,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq2_harmonics_amplitude_3(Base): @@ -1160,8 +1180,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq3_harmonics_amplitude_0(Base): @@ -1172,8 +1192,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq3_harmonics_amplitude_1(Base): @@ -1184,8 +1204,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq3_harmonics_amplitude_2(Base): @@ -1196,8 +1216,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq3_harmonics_amplitude_3(Base): @@ -1208,8 +1228,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq1_harmonics_rel_phase_0(Base): @@ -1220,8 +1240,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq1_harmonics_rel_phase_1(Base): @@ -1232,8 +1252,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq1_harmonics_rel_phase_2(Base): @@ -1244,8 +1264,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq1_harmonics_rel_phase_3(Base): @@ -1256,8 +1276,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq2_harmonics_rel_phase_0(Base): @@ -1269,8 +1289,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq2_harmonics_rel_phase_1(Base): @@ -1281,8 +1301,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq2_harmonics_rel_phase_2(Base): @@ -1293,8 +1313,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq2_harmonics_rel_phase_3(Base): @@ -1305,8 +1325,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq3_harmonics_rel_phase_0(Base): @@ -1317,8 +1337,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq3_harmonics_rel_phase_1(Base): @@ -1329,8 +1349,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq3_harmonics_rel_phase_2(Base): @@ -1341,8 +1361,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Freq3_harmonics_rel_phase_3(Base): @@ -1353,8 +1373,8 @@ def __init__(self): 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" + except Exception: + print(ERR_MSG_HARM) class Gskew(Base): @@ -1420,8 +1440,12 @@ def __init__(self): def fit(self, data): try: return m_31 - except: - print "error: please run StructureFunction_index_21 first to generate values for all Structure Function" + except Exception: + print( + ERR_MSG_RUN.format( + 'StructureFunction_index_21', 'Structure Function' + ) + ) class StructureFunction_index_32(Base): @@ -1432,5 +1456,9 @@ def __init__(self): def fit(self, data): try: return m_32 - except: - print "error: please run StructureFunction_index_21 first to generate values for all Structure Function" + except Exception: + print( + ERR_MSG_RUN.format( + 'StructureFunction_index_21', 'Structure Function' + ) + ) diff --git a/FATS/PreprocessLC.py b/FATS/PreprocessLC.py index ab7fe3e..50a6c32 100644 --- a/FATS/PreprocessLC.py +++ b/FATS/PreprocessLC.py @@ -2,31 +2,29 @@ class Preprocess_LC: - def __init__(self, data, mjd, error): - - self.N = len(mjd) - self.m = np.mean(error) - self.mjd = mjd - self.data = data - self.error = error - - def Preprocess(self): - - mjd_out = [] - data_out = [] - error_out = [] - - for i in xrange(len(self.data)): - - if self.error[i] < (3 * self.m) and (np.absolute(self.data[i] - np.mean(self.data)) / np.std(self.data)) < 5 : - - mjd_out.append(self.mjd[i]) - data_out.append(self.data[i]) - error_out.append(self.error[i]) - - - data_out = np.asarray(data_out) - mjd_out = np.asarray(mjd_out) - error_out = np.asarray(error_out) - - return [data_out, mjd_out, error_out] + def __init__(self, data, mjd, error): + self.N = len(mjd) + self.m = np.mean(error) + self.mjd = mjd + self.data = data + self.error = error + + def Preprocess(self): + + mjd_out = [] + data_out = [] + error_out = [] + + for i in range(len(self.data)): + abs_over_std = np.absolute(self.data[i] - np.mean(self.data)) + abs_over_std = abs_over_std / np.std(self.data) + if self.error[i] < (3 * self.m) and abs_over_std < 5: + mjd_out.append(self.mjd[i]) + data_out.append(self.data[i]) + error_out.append(self.error[i]) + + data_out = np.asarray(data_out) + mjd_out = np.asarray(mjd_out) + error_out = np.asarray(error_out) + + return [data_out, mjd_out, error_out] diff --git a/FATS/alignLC.py b/FATS/alignLC.py index ee9c5ce..f6f06a9 100644 --- a/FATS/alignLC.py +++ b/FATS/alignLC.py @@ -13,7 +13,7 @@ def Align_LC(mjd, mjd2, data, data2, error, error2): new_data = np.copy(data) count = 0 - for index in xrange(len(data)): + for index in range(len(data)): where = np.where(mjd2 == mjd[index]) @@ -31,7 +31,6 @@ def Align_LC(mjd, mjd2, data, data2, error, error2): new_data2 = np.asarray(new_data2).flatten() new_error2 = np.asarray(new_error2).flatten() - else: new_data = [] @@ -41,7 +40,7 @@ def Align_LC(mjd, mjd2, data, data2, error, error2): new_error2 = np.copy(error2) new_data2 = np.copy(data2) count = 0 - for index in xrange(len(data2)): + for index in range(len(data2)): where = np.where(mjd == mjd2[index]) if np.array_equal(where[0], []) is False: @@ -56,7 +55,7 @@ def Align_LC(mjd, mjd2, data, data2, error, error2): new_data = np.asarray(new_data).flatten() new_mjd = np.asarray(new_mjd).flatten() - new_error = np.asarray(new_error).flatten() + new_error = np.asarray(new_error).flatten() return new_data, new_data2, new_mjd, new_error, new_error2 - #return new_mjd, new_data, new_error, new_mjd2, new_data2, new_error2 + # return new_mjd, new_data, new_error, new_mjd2, new_data2, new_error2 diff --git a/FATS/featureFunction.py b/FATS/featureFunction.py index bb34086..85a009f 100644 --- a/FATS/featureFunction.py +++ b/FATS/featureFunction.py @@ -1,6 +1,6 @@ -import os,sys,time +import os, sys, time import numpy as np import pandas as pd import matplotlib.pyplot as plt -import Base -from FeatureFunctionLib import * \ No newline at end of file +from . import Base +from .FeatureFunctionLib import * \ No newline at end of file diff --git a/FATS/import_lc_cluster.py b/FATS/import_lc_cluster.py index 45b1127..66a79d0 100644 --- a/FATS/import_lc_cluster.py +++ b/FATS/import_lc_cluster.py @@ -1,13 +1,12 @@ - -#from Feature import FeatureSpace +# from Feature import FeatureSpace import numpy as np -class ReadLC_MACHO: +class ReadLC_MACHO: - def __init__(self,lc): - self.content1=lc + def __init__(self, lc): + self.content1 = lc def ReadLC(self): @@ -15,12 +14,11 @@ def ReadLC(self): mjd = [] error = [] # Opening the blue band - #fid = open(self.id,'r') + # fid = open(self.id,'r') self.content1 = self.content1[3:] - - for i in xrange(len(self.content1)): + for i in range(len(self.content1)): if not self.content1[i]: break else: @@ -32,5 +30,3 @@ def ReadLC(self): # Opening the red band return [data, mjd, error] - - \ No newline at end of file diff --git a/FATS/import_lightcurve.py b/FATS/import_lightcurve.py index 5e8b287..7312985 100644 --- a/FATS/import_lightcurve.py +++ b/FATS/import_lightcurve.py @@ -1,13 +1,11 @@ - -#from Feature import FeatureSpace +# from Feature import FeatureSpace import numpy as np class ReadLC_MACHO: def __init__(self,id): - - self.id=id + self.id = id def ReadLC(self): @@ -15,25 +13,23 @@ def ReadLC(self): fid = open(self.id,'r') saltos_linea = 3 - delimiter = ' ' - for i in range(0,saltos_linea): + + for i in range(0, saltos_linea): fid.next() LC = [] for lines in fid: str_line = lines.strip().split() floats = map(float, str_line) - #numbers = (number for number in str_line.split()) + # numbers = (number for number in str_line.split()) LC.append(floats) LC = np.asarray(LC) - data = LC[:,1] + data = LC[:,1] error = LC[:,2] mjd = LC[:,0] # Opening the red band return [data, mjd, error] - - \ No newline at end of file diff --git a/FATS/lomb.py b/FATS/lomb.py index 7654481..a975b88 100644 --- a/FATS/lomb.py +++ b/FATS/lomb.py @@ -1,206 +1,208 @@ -#!/usr/bin/env python -""" Fast algorithm for spectral analysis of unevenly sampled data - -The Lomb-Scargle method performs spectral analysis on unevenly sampled -data and is known to be a powerful way to find, and test the -significance of, weak periodic signals. The method has previously been -thought to be 'slow', requiring of order 10(2)N(2) operations to analyze -N data points. We show that Fast Fourier Transforms (FFTs) can be used -in a novel way to make the computation of order 10(2)N log N. Despite -its use of the FFT, the algorithm is in no way equivalent to -conventional FFT periodogram analysis. - -Keywords: - DATA SAMPLING, FAST FOURIER TRANSFORMATIONS, - SPECTRUM ANALYSIS, SIGNAL PROCESSING - -Example: - > import numpy - > import lomb - > x = numpy.arange(10) - > y = numpy.sin(x) - > fx,fy, nout, jmax, prob = lomb.fasper(x,y, 6., 6.) - -Reference: - Press, W. H. & Rybicki, G. B. 1989 - ApJ vol. 338, p. 277-280. - Fast algorithm for spectral analysis of unevenly sampled data - bib code: 1989ApJ...338..277P - -""" -from numpy import * -from numpy.fft import * - -def __spread__(y, yy, n, x, m): - """ - Given an array yy(0:n-1), extirpolate (spread) a value y into - m actual array elements that best approximate the "fictional" - (i.e., possible noninteger) array element number x. The weights - used are coefficients of the Lagrange interpolating polynomial - Arguments: - y : - yy : - n : - x : - m : - Returns: - - """ - nfac=[0,1,1,2,6,24,120,720,5040,40320,362880] - if m > 10. : - print 'factorial table too small in spread' - return - - ix=long(x) - if x == float(ix): - yy[ix]=yy[ix]+y - else: - ilo = long(x-0.5*float(m)+1.0) - ilo = min( max( ilo , 1 ), n-m+1 ) - ihi = ilo+m-1 - nden = nfac[m] - fac=x-ilo - for j in range(ilo+1,ihi+1): fac = fac*(x-j) - yy[ihi] = yy[ihi] + y*fac/(nden*(x-ihi)) - for j in range(ihi-1,ilo-1,-1): - nden=(nden/(j+1-ilo))*(j-ihi) - yy[j] = yy[j] + y*fac/(nden*(x-j)) - -def fasper(x,y,ofac,hifac, MACC=4): - """ function fasper - Given abscissas x (which need not be equally spaced) and ordinates - y, and given a desired oversampling factor ofac (a typical value - being 4 or larger). this routine creates an array wk1 with a - sequence of nout increasing frequencies (not angular frequencies) - up to hifac times the "average" Nyquist frequency, and creates - an array wk2 with the values of the Lomb normalized periodogram at - those frequencies. The arrays x and y are not altered. This - routine also returns jmax such that wk2(jmax) is the maximum - element in wk2, and prob, an estimate of the significance of that - maximum against the hypothesis of random noise. A small value of prob - indicates that a significant periodic signal is present. - - Reference: - Press, W. H. & Rybicki, G. B. 1989 - ApJ vol. 338, p. 277-280. - Fast algorithm for spectral analysis of unevenly sampled data - (1989ApJ...338..277P) - - Arguments: - X : Abscissas array, (e.g. an array of times). - Y : Ordinates array, (e.g. corresponding counts). - Ofac : Oversampling factor. - Hifac : Hifac * "average" Nyquist frequency = highest frequency - for which values of the Lomb normalized periodogram will - be calculated. - - Returns: - Wk1 : An array of Lomb periodogram frequencies. - Wk2 : An array of corresponding values of the Lomb periodogram. - Nout : Wk1 & Wk2 dimensions (number of calculated frequencies) - Jmax : The array index corresponding to the MAX( Wk2 ). - Prob : False Alarm Probability of the largest Periodogram value - MACC : Number of interpolation points per 1/4 cycle - of highest frequency - - History: - 02/23/2009, v1.0, MF - Translation of IDL code (orig. Numerical recipies) - """ - #Check dimensions of input arrays - n = long(len(x)) - if n != len(y): - print 'Incompatible arrays.' - return - - nout = 0.5*ofac*hifac*n - if round(nout) != nout: - print("Warning: nout is not an integer and will be rounded down.") - nout = int(nout) - nfreqt = long(ofac*hifac*n*MACC) #Size the FFT as next power - nfreq = 64L # of 2 above nfreqt. - - while nfreq < nfreqt: - nfreq = 2*nfreq - - ndim = long(2*nfreq) - - #Compute the mean, variance - ave = y.mean() - ##sample variance because the divisor is N-1 - var = ((y-y.mean())**2).sum()/(len(y)-1) - # and range of the data. - xmin = x.min() - xmax = x.max() - xdif = xmax-xmin - - #extirpolate the data into the workspaces - wk1 = zeros(ndim, dtype='complex') - wk2 = zeros(ndim, dtype='complex') - - fac = ndim/(xdif*ofac) - fndim = ndim - ck = ((x-xmin)*fac) % fndim - ckk = (2.0*ck) % fndim - - for j in range(0L, n): - __spread__(y[j]-ave,wk1,ndim,ck[j],MACC) - __spread__(1.0,wk2,ndim,ckk[j],MACC) - - #Take the Fast Fourier Transforms - wk1 = ifft( wk1 )*len(wk1) - wk2 = ifft( wk2 )*len(wk1) - - wk1 = wk1[1:nout+1] - wk2 = wk2[1:nout+1] - rwk1 = wk1.real - iwk1 = wk1.imag - rwk2 = wk2.real - iwk2 = wk2.imag - - df = 1.0/(xdif*ofac) - - #Compute the Lomb value for each frequency - hypo2 = 2.0 * abs( wk2 ) - hc2wt = rwk2/hypo2 - hs2wt = iwk2/hypo2 - - cwt = sqrt(0.5+hc2wt) - swt = sign(hs2wt)*(sqrt(0.5-hc2wt)) - den = 0.5*n+hc2wt*rwk2+hs2wt*iwk2 - cterm = (cwt*rwk1+swt*iwk1)**2./den - sterm = (cwt*iwk1-swt*rwk1)**2./(n-den) - - wk1 = df*(arange(nout, dtype='float')+1.) - wk2 = (cterm+sterm)/(2.0*var) - pmax = wk2.max() - jmax = wk2.argmax() - - - #Significance estimation - #expy = exp(-wk2) - #effm = 2.0*(nout)/ofac - #sig = effm*expy - #ind = (sig > 0.01).nonzero() - #sig[ind] = 1.0-(1.0-expy[ind])**effm - - #Estimate significance of largest peak value - expy = exp(-pmax) - effm = 2.0*(nout)/ofac - prob = effm*expy - - if prob > 0.01: - prob = 1.0-(1.0-expy)**effm - - return wk1,wk2,nout,jmax,prob - -def getSignificance(wk1, wk2, nout, ofac): - """ returns the peak false alarm probabilities - Hence the lower is the probability and the more significant is the peak - """ - expy = exp(-wk2) - effm = 2.0*(nout)/ofac - sig = effm*expy - ind = (sig > 0.01).nonzero() - sig[ind] = 1.0-(1.0-expy[ind])**effm - return sig \ No newline at end of file +#!/usr/bin/env python +""" Fast algorithm for spectral analysis of unevenly sampled data + +The Lomb-Scargle method performs spectral analysis on unevenly sampled +data and is known to be a powerful way to find, and test the +significance of, weak periodic signals. The method has previously been +thought to be 'slow', requiring of order 10(2)N(2) operations to analyze +N data points. We show that Fast Fourier Transforms (FFTs) can be used +in a novel way to make the computation of order 10(2)N log N. Despite +its use of the FFT, the algorithm is in no way equivalent to +conventional FFT periodogram analysis. + +Keywords: + DATA SAMPLING, FAST FOURIER TRANSFORMATIONS, + SPECTRUM ANALYSIS, SIGNAL PROCESSING + +Example: + > import numpy + > import lomb + > x = numpy.arange(10) + > y = numpy.sin(x) + > fx,fy, nout, jmax, prob = lomb.fasper(x,y, 6., 6.) + +Reference: + Press, W. H. & Rybicki, G. B. 1989 + ApJ vol. 338, p. 277-280. + Fast algorithm for spectral analysis of unevenly sampled data + bib code: 1989ApJ...338..277P + +""" +import numpy as np +from numpy import fft + +def __spread__(y, yy, n, x, m): + """ + Given an array yy(0:n-1), extirpolate (spread) a value y into + m actual array elements that best approximate the "fictional" + (i.e., possible noninteger) array element number x. The weights + used are coefficients of the Lagrange interpolating polynomial + Arguments: + y : + yy : + n : + x : + m : + Returns: + + """ + nfac = [0, 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880] + if m > 10.: + print('factorial table too small in spread') + return + + ix = int(x) + if x == float(ix): + yy[ix] = yy[ix]+y + else: + ilo = int(x-0.5*float(m)+1.0) + ilo = min(max(ilo, 1), n-m+1) + ihi = ilo+m-1 + nden = nfac[m] + fac = x-ilo + for j in range(ilo+1, ihi+1): + fac = fac*(x-j) + yy[ihi] = yy[ihi] + y*fac/(nden*(x-ihi)) + for j in range(ihi-1, ilo-1, -1): + nden = (nden/(j+1-ilo))*(j-ihi) + yy[j] = yy[j] + y*fac/(nden*(x-j)) + +def fasper(x,y,ofac,hifac, MACC=4): + """ function fasper + Given abscissas x (which need not be equally spaced) and ordinates + y, and given a desired oversampling factor ofac (a typical value + being 4 or larger). this routine creates an array wk1 with a + sequence of nout increasing frequencies (not angular frequencies) + up to hifac times the "average" Nyquist frequency, and creates + an array wk2 with the values of the Lomb normalized periodogram at + those frequencies. The arrays x and y are not altered. This + routine also returns jmax such that wk2(jmax) is the maximum + element in wk2, and prob, an estimate of the significance of that + maximum against the hypothesis of random noise. A small value of prob + indicates that a significant periodic signal is present. + + Reference: + Press, W. H. & Rybicki, G. B. 1989 + ApJ vol. 338, p. 277-280. + Fast algorithm for spectral analysis of unevenly sampled data + (1989ApJ...338..277P) + + Arguments: + X : Abscissas array, (e.g. an array of times). + Y : Ordinates array, (e.g. corresponding counts). + Ofac : Oversampling factor. + Hifac : Hifac * "average" Nyquist frequency = highest frequency + for which values of the Lomb normalized periodogram will + be calculated. + + Returns: + Wk1 : An array of Lomb periodogram frequencies. + Wk2 : An array of corresponding values of the Lomb periodogram. + Nout : Wk1 & Wk2 dimensions (number of calculated frequencies) + Jmax : The array index corresponding to the MAX( Wk2 ). + Prob : False Alarm Probability of the largest Periodogram value + MACC : Number of interpolation points per 1/4 cycle + of highest frequency + + History: + 02/23/2009, v1.0, MF + Translation of IDL code (orig. Numerical recipies) + """ + # Check dimensions of input arrays + n = int(len(x)) + if n != len(y): + print("Incompatible arrays.") + return + + nout = 0.5*ofac*hifac*n + if round(nout) != nout: + print("Warning: nout is not an integer and will be rounded down.") + + nout = int(nout) + nfreqt = int(ofac*hifac*n*MACC) # Size the FFT as next power + nfreq = 64 # of 2 above nfreqt. + + while nfreq < nfreqt: + nfreq = 2*nfreq + + ndim = int(2*nfreq) + + # Compute the mean, variance + ave = y.mean() + # sample variance because the divisor is N-1 + var = ((y-y.mean())**2).sum()/(len(y)-1) + # and range of the data. + xmin = x.min() + xmax = x.max() + xdif = xmax-xmin + + # extirpolate the data into the workspaces + wk1 = np.zeros(ndim, dtype='complex') + wk2 = np.zeros(ndim, dtype='complex') + + fac = ndim/(xdif*ofac) + fndim = ndim + ck = ((x-xmin)*fac) % fndim + ckk = (2.0*ck) % fndim + + for j in range(0, n): + __spread__(y[j]-ave, wk1, ndim, ck[j], MACC) + __spread__(1.0, wk2, ndim, ckk[j], MACC) + + # Take the Fast Fourier Transforms + wk1 = fft.ifft(wk1)*len(wk1) + wk2 = fft.ifft(wk2)*len(wk1) + + wk1 = wk1[1:nout+1] + wk2 = wk2[1:nout+1] + rwk1 = wk1.real + iwk1 = wk1.imag + rwk2 = wk2.real + iwk2 = wk2.imag + + df = 1.0/(xdif*ofac) + + # Compute the Lomb value for each frequency + hypo2 = 2.0 * abs(wk2) + hc2wt = rwk2/hypo2 + hs2wt = iwk2/hypo2 + + cwt = np.sqrt(0.5+hc2wt) + swt = np.sign(hs2wt)*(np.sqrt(0.5-hc2wt)) + den = 0.5*n+hc2wt*rwk2+hs2wt*iwk2 + cterm = (cwt*rwk1+swt*iwk1)**2./den + sterm = (cwt*iwk1-swt*rwk1)**2./(n-den) + + wk1 = df*(np.arange(nout, dtype='float')+1.) + wk2 = (cterm+sterm)/(2.0*var) + pmax = wk2.max() + jmax = wk2.argmax() + + # Significance estimation + # expy = exp(-wk2) + # effm = 2.0*(nout)/ofac + # sig = effm*expy + # ind = (sig > 0.01).nonzero() + # sig[ind] = 1.0-(1.0-expy[ind])**effm + + # Estimate significance of largest peak value + expy = np.exp(-pmax) + effm = 2.0*(nout)/ofac + prob = effm*expy + + if prob > 0.01: + prob = 1.0-(1.0-expy)**effm + + return (wk1, wk2, nout, jmax, prob) + + +def getSignificance(wk1, wk2, nout, ofac): + """ returns the peak false alarm probabilities + Hence the lower is the probability and the more significant is the peak + """ + expy = np.exp(-wk2) + effm = 2.0*(nout)/ofac + sig = effm*expy + ind = (sig > 0.01).nonzero() + sig[ind] = 1.0-(1.0-expy[ind])**effm + return sig diff --git a/FATS/test_library.py b/FATS/test_library.py index 287c590..f6ee243 100644 --- a/FATS/test_library.py +++ b/FATS/test_library.py @@ -1,10 +1,10 @@ #!/usr/bin/env python -from Feature import FeatureSpace +from .Feature import FeatureSpace import numpy as np -from import_lc_cluster import ReadLC_MACHO -from PreprocessLC import Preprocess_LC -from alignLC import Align_LC +from .import_lc_cluster import ReadLC_MACHO +from .PreprocessLC import Preprocess_LC +from .alignLC import Align_LC import os.path import tarfile import sys @@ -13,297 +13,395 @@ @pytest.fixture def white_noise(): - data = np.random.normal(size=10000) - mjd=np.arange(10000) - error = np.random.normal(loc=0.01, scale =0.8, size=10000) - second_data = np.random.normal(size=10000) - mjd2=np.arange(10000) - error2 = np.random.normal(loc=0.01, scale =0.8, size=10000) - aligned_data = data - aligned_second_data = second_data - aligned_mjd = mjd - lc = np.array([data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd]) - return lc + data = np.random.normal(size=10000) + mjd = np.arange(10000) + error = np.random.normal(loc=0.01, scale=0.8, size=10000) + error2 = np.random.normal(loc=0.01, scale=0.8, size=10000) + second_data = np.random.normal(size=10000) + aligned_data = data + aligned_second_data = second_data + aligned_error = error + aligned_second_error = error2 + aligned_mjd = mjd + lc = np.array( + [ + data, mjd, error, + second_data, aligned_data, aligned_second_data, aligned_mjd, + aligned_error, aligned_second_error + ] + ) + return lc + @pytest.fixture def periodic_lc(): - N=100 - mjd_periodic = np.arange(N) - Period = 20 - cov = np.zeros([N,N]) - mean = np.zeros(N) - for i in np.arange(N): - for j in np.arange(N): - cov[i,j] = np.exp( -(np.sin( (np.pi/Period) *(i-j))**2)) - data_periodic=np.random.multivariate_normal(mean, cov) - lc = np.array([data_periodic, mjd_periodic]) - return lc + N = 100 + mjd_periodic = np.arange(N) + Period = 20 + cov = np.zeros([N, N]) + mean = np.zeros(N) + for i in np.arange(N): + for j in np.arange(N): + cov[i, j] = np.exp(-(np.sin((np.pi/Period) * (i-j))**2)) + data_periodic = np.random.multivariate_normal(mean, cov) + lc = np.array([data_periodic, mjd_periodic]) + return lc @pytest.fixture def uniform_lc(): - mjd_uniform=np.arange(1000000) - data_uniform=np.random.uniform(size=1000000) - lc = np.array([data_uniform, mjd_uniform]) - return lc + mjd_uniform = np.arange(1000000) + data_uniform = np.random.uniform(size=1000000) + lc = np.array([data_uniform, mjd_uniform]) + return lc @pytest.fixture def random_walk(): - N = 10000 - alpha = 1. - sigma = 0.5 - data_rw = np.zeros([N,1]) - data_rw[0] = 1 - time_rw = xrange(1, N) - for t in time_rw: - data_rw[t] = alpha * data_rw[t-1] + np.random.normal(loc=0.0, scale=sigma) - time_rw = np.array(range(0,N)) + 1 * np.random.uniform(size=N) - data_rw = data_rw.squeeze() - lc = np.array([data_rw, time_rw]) - return lc + N = 10000 + alpha = 1. + sigma = 0.5 + data_rw = np.zeros([N, 1]) + data_rw[0] = 1 + time_rw = range(1, N) + for t in time_rw: + data_rw[t] = alpha * data_rw[t-1] + np.random.normal( + loc=0.0, scale=sigma + ) + time_rw = np.array(range(0, N)) + 1 * np.random.uniform(size=N) + data_rw = data_rw.squeeze() + lc = np.array([data_rw, time_rw]) + return lc # def test_Amplitude(white_noise): -# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() +# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() -# a = FeatureSpace(featureList=['Amplitude']) -# a=a.calculateFeature(white_noise[0]) +# a = FeatureSpace(featureList=['Amplitude']) +# a = a.calculateFeature(white_noise[0]) -# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) +# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) # def test_Autocor(white_noise): -# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() +# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() -# a = FeatureSpace(featureList=['Autocor'] ) -# a=a.calculateFeature(white_noise[0]) +# a = FeatureSpace(featureList=['Autocor'] ) +# a = a.calculateFeature(white_noise[0]) -# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) +# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) # def test_Automean(white_noise): -# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() +# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() -# a = FeatureSpace(featureList=['Automean'] , Automean=[0,0]) -# a=a.calculateFeature(white_noise[0]) +# a = FeatureSpace(featureList=['Automean'] , Automean=[0,0]) +# a = a.calculateFeature(white_noise[0]) -# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) +# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) # def test_B_R(white_noise): -# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() +# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() -# a = FeatureSpace(featureList=['B_R'] , B_R=second_data) -# a=a.calculateFeature(white_noise[0]) +# a = FeatureSpace(featureList=['B_R'] , B_R=second_data) +# a = a.calculateFeature(white_noise[0]) -# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) +# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) def test_Beyond1Std(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['Beyond1Std']) - a=a.calculateFeature(white_noise) + a = FeatureSpace(featureList=['Beyond1Std']) + a = a.calculateFeature(white_noise) - assert(a.result(method='array') >= 0.30 and a.result(method='array') <= 0.40) + assert ( + a.result(method='array') >= 0.30 and + a.result(method='array') <= 0.40 + ) def test_Mean(white_noise): - a = FeatureSpace(featureList=['Mean']) - a=a.calculateFeature(white_noise) + a = FeatureSpace(featureList=['Mean']) + a = a.calculateFeature(white_noise) - assert(a.result(method='array') >= -0.1 and a.result(method='array') <= 0.1) + assert ( + a.result(method='array') >= -0.1 and + a.result(method='array') <= 0.1 + ) # def test_CAR(white_noise): -# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() +# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() -# a = FeatureSpace(featureList=['CAR_sigma', 'CAR_tau', 'CAR_tmean'] , CAR_sigma=[mjd, error]) -# a=a.calculateFeature(white_noise[0]) +# a = FeatureSpace(featureList=['CAR_sigma', 'CAR_tau', 'CAR_tmean'] , CAR_sigma=[mjd, error]) +# a = a.calculateFeature(white_noise[0]) -# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) +# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) def test_Con(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['Con'] , Con=1) - a=a.calculateFeature(white_noise) + a = FeatureSpace(featureList=['Con'], Con=1) + a = a.calculateFeature(white_noise) - assert(a.result(method='array') >= 0.04 and a.result(method='array') <= 0.05) + assert ( + a.result(method='array') >= 0.04 and + a.result(method='array') <= 0.05 + ) def test_Eta_color(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['Eta_color']) - a=a.calculateFeature(white_noise) + a = FeatureSpace(featureList=['Eta_color']) + a = a.calculateFeature(white_noise) - assert(a.result(method='array') >= 1.9 and a.result(method='array') <= 2.1) + assert ( + a.result(method='array') >= 1.9 and + a.result(method='array') <= 2.1 + ) def test_Eta_e(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['Eta_e']) - a=a.calculateFeature(white_noise) + a = FeatureSpace(featureList=['Eta_e']) + a = a.calculateFeature(white_noise) - assert(a.result(method='array') >= 1.9 and a.result(method='array') <= 2.1) + assert ( + a.result(method='array') >= 1.9 and + a.result(method='array') <= 2.1 + ) def test_FluxPercentile(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - - a = FeatureSpace(featureList=['FluxPercentileRatioMid20','FluxPercentileRatioMid35','FluxPercentileRatioMid50','FluxPercentileRatioMid65','FluxPercentileRatioMid80'] ) - a=a.calculateFeature(white_noise) - - assert(a.result(method='array')[0] >= 0.145 and a.result(method='array')[0] <= 0.160) - assert(a.result(method='array')[1] >= 0.260 and a.result(method='array')[1] <= 0.290) - assert(a.result(method='array')[2] >= 0.350 and a.result(method='array')[2] <= 0.450) - assert(a.result(method='array')[3] >= 0.540 and a.result(method='array')[3] <= 0.580) - assert(a.result(method='array')[4] >= 0.760 and a.result(method='array')[4] <= 0.800) - + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + + a = FeatureSpace(featureList=['FluxPercentileRatioMid20','FluxPercentileRatioMid35','FluxPercentileRatioMid50','FluxPercentileRatioMid65','FluxPercentileRatioMid80'] ) + a = a.calculateFeature(white_noise) + + assert ( + a.result(method='array')[0] >= 0.145 and + a.result(method='array')[0] <= 0.160 + ) + assert ( + a.result(method='array')[1] >= 0.260 and + a.result(method='array')[1] <= 0.290 + ) + assert ( + a.result(method='array')[2] >= 0.350 and + a.result(method='array')[2] <= 0.450 + ) + assert ( + a.result(method='array')[3] >= 0.540 and + a.result(method='array')[3] <= 0.580 + ) + assert ( + a.result(method='array')[4] >= 0.760 and + a.result(method='array')[4] <= 0.800 + ) def test_LinearTrend(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['LinearTrend']) - a=a.calculateFeature(white_noise) + a = FeatureSpace(featureList=['LinearTrend']) + a = a.calculateFeature(white_noise) - assert(a.result(method='array') >= -0.1 and a.result(method='array') <= 0.1) + assert ( + a.result(method='array') >= -0.1 and + a.result(method='array') <= 0.1 + ) # def test_MaxSlope(white_noise): -# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() +# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() -# a = FeatureSpace(featureList=['MaxSlope'] , MaxSlope=mjd) -# a=a.calculateFeature(white_noise[0]) +# a = FeatureSpace(featureList=['MaxSlope'] , MaxSlope=mjd) +# a = a.calculateFeature(white_noise[0]) -# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) +# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) def test_Meanvariance(uniform_lc): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['Meanvariance']) - a=a.calculateFeature(uniform_lc) + a = FeatureSpace(featureList=['Meanvariance']) + a = a.calculateFeature(uniform_lc) - assert(a.result(method='array') >= 0.575 and a.result(method='array') <= 0.580) + assert ( + a.result(method='array') >= 0.575 and + a.result(method='array') <= 0.580 + ) def test_MedianAbsDev(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['MedianAbsDev']) - a=a.calculateFeature(white_noise) + a = FeatureSpace(featureList=['MedianAbsDev']) + a = a.calculateFeature(white_noise) - assert(a.result(method='array') >= 0.630 and a.result(method='array') <= 0.700) + assert ( + a.result(method='array') >= 0.630 and + a.result(method='array') <= 0.700 + ) # def test_MedianBRP(white_noise): -# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() +# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() -# a = FeatureSpace(featureList=['MedianBRP'] , MaxSlope=mjd) -# a=a.calculateFeature(white_noise[0]) +# a = FeatureSpace(featureList=['MedianBRP'] , MaxSlope=mjd) +# a = a.calculateFeature(white_noise[0]) -# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) +# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) def test_PairSlopeTrend(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['PairSlopeTrend']) - a=a.calculateFeature(white_noise) + a = FeatureSpace(featureList=['PairSlopeTrend']) + a = a.calculateFeature(white_noise) - assert(a.result(method='array') >= -0.25 and a.result(method='array') <= 0.25) + assert ( + a.result(method='array') >= -0.25 and + a.result(method='array') <= 0.25 + ) # def test_PercentAmplitude(white_noise): -# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() +# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() -# a = FeatureSpace(featureList=['PercentAmplitude']) -# a=a.calculateFeature(white_noise[0]) +# a = FeatureSpace(featureList=['PercentAmplitude']) +# a = a.calculateFeature(white_noise[0]) -# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) +# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) # def test_PercentDifferenceFluxPercentile(white_noise): -# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() +# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() -# a = FeatureSpace(featureList=['PercentDifferenceFluxPercentile']) -# a=a.calculateFeature(white_noise[0]) +# a = FeatureSpace(featureList=['PercentDifferenceFluxPercentile']) +# a = a.calculateFeature(white_noise[0]) -# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) +# assert(a.result(method='array') >= 0.043 and a.result(method='array') <= 0.046) def test_Period_Psi(periodic_lc): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['PeriodLS', 'Period_fit','Psi_CS','Psi_eta']) - a=a.calculateFeature(periodic_lc) - # print a.result(method='array'), len(periodic_lc[0]) - assert(a.result(method='array')[0] >= 19 and a.result(method='array')[0] <= 21) + a = FeatureSpace(featureList=['PeriodLS', 'Period_fit','Psi_CS','Psi_eta']) + a = a.calculateFeature(periodic_lc) + # print a.result(method='array'), len(periodic_lc[0]) + assert ( + a.result(method='array')[0] >= 19 and + a.result(method='array')[0] <= 21 + ) def test_Q31(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['Q31']) - a=a.calculateFeature(white_noise) - assert(a.result(method='array') >= 1.30 and a.result(method='array') <= 1.38) + a = FeatureSpace(featureList=['Q31']) + a = a.calculateFeature(white_noise) + assert ( + a.result(method='array') >= 1.30 and + a.result(method='array') <= 1.38 + ) # def test_Q31B_R(white_noise): -# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() +# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() -# a = FeatureSpace(featureList=['Q31B_R'], Q31B_R = [aligned_second_data, aligned_data]) -# a=a.calculateFeature(white_noise[0]) +# a = FeatureSpace(featureList=['Q31B_R'], Q31B_R = [aligned_second_data, aligned_data]) +# a = a.calculateFeature(white_noise[0]) def test_Rcs(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['Rcs']) - a=a.calculateFeature(white_noise) - assert(a.result(method='array') >= 0 and a.result(method='array') <= 0.1) + a = FeatureSpace(featureList=['Rcs']) + a = a.calculateFeature(white_noise) + assert ( + a.result(method='array') >= 0 and + a.result(method='array') <= 0.1 + ) def test_Skew(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['Skew']) - a=a.calculateFeature(white_noise) - assert(a.result(method='array') >= -0.1 and a.result(method='array') <= 0.1) + a = FeatureSpace(featureList=['Skew']) + a = a.calculateFeature(white_noise) + assert ( + a.result(method='array') >= -0.1 and + a.result(method='array') <= 0.1 + ) # def test_SlottedA(white_noise): -# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() +# # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() -# a = FeatureSpace(featureList=['SlottedA'], SlottedA = [mjd, 1]) -# a=a.calculateFeature(white_noise[0]) +# a = FeatureSpace(featureList=['SlottedA'], SlottedA = [mjd, 1]) +# a = a.calculateFeature(white_noise[0]) def test_SmallKurtosis(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['SmallKurtosis']) - a=a.calculateFeature(white_noise) - assert(a.result(method='array') >= -0.2 and a.result(method='array') <= 0.2) + a = FeatureSpace(featureList=['SmallKurtosis']) + a = a.calculateFeature(white_noise) + assert ( + a.result(method='array') >= -0.2 and + a.result(method='array') <= 0.2 + ) def test_Std(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['Std']) - a=a.calculateFeature(white_noise) + a = FeatureSpace(featureList=['Std']) + a = a.calculateFeature(white_noise) - assert(a.result(method='array') >= 0.9 and a.result(method='array') <= 1.1) + assert ( + a.result(method='array') >= 0.9 and + a.result(method='array') <= 1.1 + ) def test_Stetson(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - - a = FeatureSpace(featureList=['SlottedA_length','StetsonK', 'StetsonK_AC', 'StetsonJ', 'StetsonL']) - a=a.calculateFeature(white_noise) - - assert(a.result(method='array')[1] >= 0.790 and a.result(method='array')[1] <= 0.85) - assert(a.result(method='array')[2] >= 0.20 and a.result(method='array')[2] <= 0.45) - assert(a.result(method='array')[3] >= -0.1 and a.result(method='array')[3] <= 0.1) - assert(a.result(method='array')[4] >= -0.1 and a.result(method='array')[4] <= 0.1) + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + + a = FeatureSpace(featureList=['SlottedA_length','StetsonK', 'StetsonK_AC', 'StetsonJ', 'StetsonL']) + a = a.calculateFeature(white_noise) + + assert ( + a.result(method='array')[1] >= 0.790 and + a.result(method='array')[1] <= 0.85 + ) + assert ( + a.result(method='array')[2] >= 0.20 and + a.result(method='array')[2] <= 0.45 + ) + assert ( + a.result(method='array')[3] >= -0.1 and + a.result(method='array')[3] <= 0.1 + ) + assert ( + a.result(method='array')[4] >= -0.1 and + a.result(method='array')[4] <= 0.1 + ) def test_Gskew(white_noise): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - a = FeatureSpace(featureList=['Gskew']) - a=a.calculateFeature(white_noise) - assert(a.result(method='array') >= -0.2 and a.result(method='array') <= 0.2) + a = FeatureSpace(featureList=['Gskew']) + a = a.calculateFeature(white_noise) + assert ( + a.result(method='array') >= -0.2 and + a.result(method='array') <= 0.2 + ) def test_StructureFunction(random_walk): - # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() - - a = FeatureSpace(featureList=['StructureFunction_index_21', 'StructureFunction_index_31', - 'StructureFunction_index_32']) - a = a.calculateFeature(random_walk) - assert(a.result(method='array')[0] >= 1.520 and a.result(method='array')[0] <= 2.067) - assert(a.result(method='array')[1] >= 1.821 and a.result(method='array')[1] <= 3.162) - assert(a.result(method='array')[2] >= 1.243 and a.result(method='array')[2] <= 1.562) + # data, mjd, error, second_data, aligned_data, aligned_second_data, aligned_mjd = white_noise() + + a = FeatureSpace( + featureList=[ + 'StructureFunction_index_21', + 'StructureFunction_index_31', + 'StructureFunction_index_32' + ] + ) + a = a.calculateFeature(random_walk) + assert ( + a.result(method='array')[0] >= 1.520 and + a.result(method='array')[0] <= 2.067 + ) + assert ( + a.result(method='array')[1] >= 1.821 and + a.result(method='array')[1] <= 3.162 + ) + assert ( + a.result(method='array')[2] >= 1.243 and + a.result(method='array')[2] <= 1.562 + ) diff --git a/requirements.txt b/requirements.txt index 6d2f179..0c0f71a 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,3 @@ -matplotlib==2.1.0 -pandas==0.13.1 -statsmodels==0.8.0 +matplotlib +pandas +statsmodels