/usr/lib/python3/dist-packages/astroML/time_series/ACF.py is in python3-astroml 0.3-6.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 | """
Auto-correlation functions
"""
from __future__ import division
import numpy as np
from scipy import fftpack
from .periodogram import lomb_scargle
def ACF_scargle(t, y, dy, n_omega=2 ** 10, omega_max=100):
"""Compute the Auto-correlation function via Scargle's method
Parameters
----------
t : array_like
times of observation. Assumed to be in increasing order.
y : array_like
values of each observation. Should be same shape as t
dy : float or array_like
errors in each observation.
n_omega : int (optional)
number of angular frequencies at which to evaluate the periodogram
default is 2^10
omega_max : float (optional)
maximum value of omega at which to evaluate the periodogram
default is 100
Returns
-------
ACF, t : ndarrays
The auto-correlation function and associated times
"""
t = np.asarray(t)
y = np.asarray(y)
if y.shape != t.shape:
raise ValueError("shapes of t and y must match")
dy = np.asarray(dy) * np.ones(y.shape)
d_omega = omega_max * 1. / (n_omega + 1)
omega = d_omega * np.arange(1, n_omega + 1)
# recall that P(omega = 0) = (chi^2(0) - chi^2(0)) / chi^2(0)
# = 0
# compute P and shifted full-frequency array
P = lomb_scargle(t, y, dy, omega,
generalized=True)
P = np.concatenate([[0], P, P[-2::-1]])
# compute PW, the power of the window function
PW = lomb_scargle(t, np.ones(len(t)), dy, omega,
generalized=False, subtract_mean=False)
PW = np.concatenate([[0], PW, PW[-2::-1]])
# compute the inverse fourier transform of P and PW
rho = fftpack.ifft(P).real
rhoW = fftpack.ifft(PW).real
ACF = fftpack.fftshift(rho / rhoW) / np.sqrt(2)
N = len(ACF)
dt = 2 * np.pi / N / (omega[1] - omega[0])
t = dt * (np.arange(N) - N // 2)
return ACF, t
def ACF_EK(t, y, dy, bins=20):
"""Auto-correlation function via the Edelson-Krolik method
Parameters
----------
t : array_like
times of observation. Assumed to be in increasing order.
y : array_like
values of each observation. Should be same shape as t
dy : float or array_like
errors in each observation.
bins : int or array_like (optional)
if integer, the number of bins to use in the analysis.
if array, the (nbins + 1) bin edges.
Default is bins=20.
Returns
-------
ACF : ndarray
The auto-correlation function and associated times
err : ndarray
the error in the ACF
bins : ndarray
bin edges used in computation
"""
t = np.asarray(t)
y = np.asarray(y)
if y.shape != t.shape:
raise ValueError("shapes of t and y must match")
if t.ndim != 1:
raise ValueError("t should be a 1-dimensional array")
dy = np.asarray(dy) * np.ones(y.shape)
# compute mean and standard deviation of y
w = 1. / dy / dy
w /= w.sum()
mu = np.dot(w, y)
sigma = np.std(y, ddof=1)
dy2 = dy[:, None]
dt = t - t[:, None]
UDCF = ((y - mu) * (y - mu)[:, None] /
np.sqrt((sigma ** 2 - dy ** 2) *
(sigma ** 2 - dy2 ** 2)))
# determine binning
bins = np.asarray(bins)
if bins.size == 1:
dt_min = dt.min()
dt_max = dt.max()
bins = np.linspace(dt_min, dt_max + 1E-10, bins + 1)
ACF = np.zeros(len(bins) - 1)
M = np.zeros(len(bins) - 1)
for i in range(len(bins) - 1):
flag = (dt >= bins[i]) & (dt < bins[i + 1])
M[i] = flag.sum()
ACF[i] = np.sum(UDCF[flag])
ACF /= M
return ACF, np.sqrt(2. / M), bins
|