/
bgls.py
100 lines (68 loc) · 2.15 KB
/
bgls.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
# -*- coding: utf-8 -*-
#================================================================================
# Copyright (c) 2014-2015 Annelies Mortier, João Faria
# Distributed under the MIT License.
# (See accompanying file LICENSE or copy at http://opensource.org/licenses/MIT)
#================================================================================
import numpy as np
import sys
try:
import mpmath # https://code.google.com/p/mpmath/
except ImportError, e1:
try:
from sympy import mpmath # http://sympy.org/en/index.html
except ImportError, e2:
raise e2
finally:
raise e1
pi = np.pi
def bgls(t, y, err, plow=0.5, phigh=100, ofac=1):
n_steps = int(ofac*len(t)*(1./plow - 1./phigh))
f = np.linspace(1./phigh, 1./plow, n_steps)
omegas = 2. * pi * f
err2 = err * err
w = 1./err2
W = sum(w)
bigY = sum(w*y) # Eq. (10)
p = []
constants = []
exponents = []
for i, omega in enumerate(omegas):
theta = 0.5 * np.arctan2(sum(w*np.sin(2.*omega*t)), sum(w*np.cos(2.*omega*t)))
x = omega*t - theta
cosx = np.cos(x)
sinx = np.sin(x)
wcosx = w*cosx
wsinx = w*sinx
C = sum(wcosx)
S = sum(wsinx)
YCh = sum(y*wcosx)
YSh = sum(y*wsinx)
CCh = sum(wcosx*cosx)
SSh = sum(wsinx*sinx)
if (CCh != 0 and SSh != 0):
K = (C*C*SSh + S*S*CCh - W*CCh*SSh)/(2.*CCh*SSh)
L = (bigY*CCh*SSh - C*YCh*SSh - S*YSh*CCh)/(CCh*SSh)
M = (YCh*YCh*SSh + YSh*YSh*CCh)/(2.*CCh*SSh)
constants.append(1./np.sqrt(CCh*SSh*abs(K)))
elif (CCh == 0):
K = (S*S - W*SSh)/(2.*SSh)
L = (bigY*SSh - S*YSh)/(SSh)
M = (YSh*YSh)/(2.*SSh)
constants.append(1./np.sqrt(SSh*abs(K)))
elif (SSh == 0):
K = (C*C - W*CCh)/(2.*CCh)
L = (bigY*CCh - C*YCh)/(CCh)
M = (YCh*YCh)/(2.*CCh)
constants.append(1./np.sqrt(CCh*abs(K)))
if K > 0:
raise RuntimeError('K is positive. This should not happen.')
exponents.append(M - L*L/(4.*K))
constants = np.array(constants)
exponents = np.array(exponents)
logp = np.log10(constants) + (exponents * np.log10(np.exp(1.)))
p = [10**mpmath.mpf(x) for x in logp]
p = np.array(p) / max(p) # normalize
p[p < (sys.float_info.min * 10)] = 0
p = np.array([float(pp) for pp in p])
return 1./f, p