/usr/lib/python3/dist-packages/astroML/time_series/generate.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 | import numpy as np
from ..utils import check_random_state
def generate_power_law(N, dt, beta, generate_complex=False, random_state=None):
"""Generate a power-law light curve
This uses the method from Timmer & Koenig [1]_
Parameters
----------
N : integer
Number of equal-spaced time steps to generate
dt : float
Spacing between time-steps
beta : float
Power-law index. The spectrum will be (1 / f)^beta
generate_complex : boolean (optional)
if True, generate a complex time series rather than a real time series
random_state : None, int, or np.random.RandomState instance (optional)
random seed or random number generator
Returns
-------
x : ndarray
the length-N
References
----------
.. [1] Timmer, J. & Koenig, M. On Generating Power Law Noise. A&A 300:707
"""
random_state = check_random_state(random_state)
dt = float(dt)
N = int(N)
Npos = int(N / 2)
Nneg = int((N - 1) / 2)
domega = (2 * np.pi / dt / N)
if generate_complex:
omega = domega * np.fft.ifftshift(np.arange(N) - int(N / 2))
else:
omega = domega * np.arange(Npos + 1)
x_fft = np.zeros(len(omega), dtype=complex)
x_fft.real[1:] = random_state.normal(0, 1, len(omega) - 1)
x_fft.imag[1:] = random_state.normal(0, 1, len(omega) - 1)
x_fft[1:] *= (1. / omega[1:]) ** (0.5 * beta)
x_fft[1:] *= (1. / np.sqrt(2))
# by symmetry, the Nyquist frequency is real if x is real
if (not generate_complex) and (N % 2 == 0):
x_fft.imag[-1] = 0
if generate_complex:
x = np.fft.ifft(x_fft)
else:
x = np.fft.irfft(x_fft, N)
return x
def generate_damped_RW(t_rest, tau=300., z=2.0,
xmean=0, SFinf=0.3, random_state=None):
"""Generate a damped random walk light curve
This uses a damped random walk model to generate a light curve similar
to that of a QSO [1]_.
Parameters
----------
t_rest : array_like
rest-frame time. Should be in increasing order
tau : float
relaxation time
z : float
redshift
xmean : float (optional)
mean value of random walk; default=0
SFinf : float (optional
Structure function at infinity; default=0.3
random_state : None, int, or np.random.RandomState instance (optional)
random seed or random number generator
Returns
-------
x : ndarray
the sampled values corresponding to times t_rest
Notes
-----
The differential equation is (with t = time/tau):
dX = -X(t) * dt + sigma * sqrt(tau) * e(t) * sqrt(dt) + b * tau * dt
where e(t) is white noise with zero mean and unit variance, and
Xmean = b * tau
SFinf = sigma * sqrt(tau / 2)
so
dX(t) = -X(t) * dt + sqrt(2) * SFint * e(t) * sqrt(dt) + Xmean * dt
References
----------
.. [1] Kelly, B., Bechtold, J. & Siemiginowska, A. (2009)
Are the Variations in Quasar Optical Flux Driven by Thermal
Fluctuations? ApJ 698:895 (2009)
"""
# Xmean = b * tau
# SFinf = sigma * sqrt(tau / 2)
t_rest = np.atleast_1d(t_rest)
if t_rest.ndim != 1:
raise ValueError('t_rest should be a 1D array')
random_state = check_random_state(random_state)
N = len(t_rest)
t_obs = t_rest * (1. + z) / tau
x = np.zeros(N)
x[0] = random_state.normal(xmean, SFinf)
E = random_state.normal(0, 1, N)
for i in range(1, N):
dt = t_obs[i] - t_obs[i - 1]
x[i] = (x[i - 1]
- dt * (x[i - 1] - xmean)
+ np.sqrt(2) * SFinf * E[i] * np.sqrt(dt))
return x
|