# -*- coding: utf-8 -*-
# This file as well as the whole tsfresh package are licenced under the MIT licence (see the LICENCE.txt)
# Maximilian Christ (maximilianchrist.com), Blue Yonder Gmbh, 2016
# Thanks to Andreas W. Kempa-Liehr for providing this snippet
import logging
import numpy as np
import pandas as pd
_logger = logging.getLogger(__name__)
[docs]
class velocity:
"""
Simulates the velocity of a dissipative soliton (kind of self organized particle) [6]_.
The equilibrium velocity without noise R=0 for
$\tau>1.0/\\kappa_3$ is $\\kappa_3 \\sqrt{(tau - 1.0/\\kappa_3)/Q}.
Before the drift-bifurcation $\tau \\le 1.0/\\kappa_3$ the velocity is zero.
References
----------
.. [6] Andreas Kempa-Liehr (2013, p. 159-170)
Dynamics of Dissipative Soliton
Dissipative Solitons in Reaction Diffusion Systems.
Springer: Berlin
>>> ds = velocity(tau=3.5) # Dissipative soliton with equilibrium velocity 1.5e-3
>>> print(ds.label) # Discriminating before or beyond Drift-Bifurcation
1
# Equilibrium velocity
>>> print(ds.deterministic)
0.0015191090506254991
# Simulated velocity as a time series with 20000 time steps being disturbed by Gaussian white noise
>>> v = ds.simulate(20000)
"""
def __init__(self, tau=3.8, kappa_3=0.3, Q=1950.0, R=3e-4, delta_t=0.05, seed=None):
"""
:param tau: Bifurcation parameter determining the intrinsic velocity of the dissipative soliton,
which is zero for tau<=1.0/kappa_3 and np.sqrt(kappa_3**3/Q * (tau - 1.0/kappa_3)) otherwise
:type tau: float
:param kappa_3: Inverse bifurcation point.
:type kappa_3:
:param Q: Shape parameter of dissipative soliton
:type Q: float
:param R: Noise amplitude
:type R: float
:param delta_t: temporal discretization
:type delta_t: float
"""
# done: add start seed
self.delta_t = delta_t
self.kappa_3 = kappa_3
self.Q = Q
self.tau = tau
self.a = self.delta_t * kappa_3**2 * (tau - 1.0 / kappa_3)
self.b = self.delta_t * Q / kappa_3
self.label = int(tau > 1.0 / kappa_3)
self.c = np.sqrt(self.delta_t) * R
self.delta_t = self.delta_t
if seed is not None:
np.random.seed(seed)
if tau <= 1.0 / kappa_3:
self.deterministic = 0.0
else:
self.deterministic = kappa_3**1.5 * np.sqrt((tau - 1.0 / kappa_3) / Q)
def __call__(self, v):
"""
returns deterministic dynamic = acceleration (without noise)
:param v: initial velocity vector
:rtype v: ndarray
:return: velocity vector of next time step
:return type: ndarray
"""
return v * (1.0 + self.a - self.b * np.dot(v, v))
[docs]
def simulate(self, N, v0=np.zeros(2)):
"""
:param N: number of time steps
:type N: int
:param v0: initial velocity vector
:type v0: ndarray
:return: time series of velocity vectors with shape (N, v0.shape[0])
:rtype: ndarray
"""
v = [v0] # first value is initial condition
n = N - 1 # Because we are returning the initial condition,
# only (N-1) time steps are computed
gamma = np.random.randn(n, v0.size)
for i in range(n):
next_v = self.__call__(v[i]) + self.c * gamma[i]
v.append(next_v)
v_vec = np.array(v)
return v_vec
[docs]
def sample_tau(n=10, kappa_3=0.3, ratio=0.5, rel_increase=0.15):
"""
Return list of control parameters
:param n: number of samples
:type n: int
:param kappa_3: inverse bifurcation point
:type kappa_3: float
:param ratio: ratio (default 0.5) of samples before and beyond drift-bifurcation
:type ratio: float
:param rel_increase: relative increase from bifurcation point
:type rel_increase: float
:return: tau. List of sampled bifurcation parameter
:rtype tau: list
"""
assert ratio > 0 and ratio <= 1
assert kappa_3 > 0
assert rel_increase > 0 and rel_increase <= 1
tau_c = 1.0 / kappa_3
tau_max = tau_c * (1.0 + rel_increase)
tau = tau_c + (tau_max - tau_c) * (np.random.rand(n) - ratio)
return tau.tolist()
[docs]
def load_driftbif(n, length, m=2, classification=True, kappa_3=0.3, seed=False):
"""
Simulates n time-series with length time steps each for the m-dimensional velocity of a dissipative soliton
classification=True:
target 0 means tau<=1/0.3, Dissipative Soliton with Brownian motion (purely noise driven)
target 1 means tau> 1/0.3, Dissipative Soliton with Active Brownian motion (intrinsiv velocity with overlaid noise)
classification=False:
target is bifurcation parameter tau
:param n: number of samples
:type n: int
:param length: length of the time series
:type length: int
:param m: number of spatial dimensions (default m=2) the dissipative soliton is propagating in
:type m: int
:param classification: distinguish between classification (default True) and regression target
:type classification: bool
:param kappa_3: inverse bifurcation parameter (default 0.3)
:type kappa_3: float
:param seed: random seed (default False)
:type seed: float
:return: X, y. Time series container and target vector
:rtype X: pandas.DataFrame
:rtype y: pandas.DataFrame
"""
# todo: add ratio of classes
if m > 2:
logging.warning(
"You set the dimension parameter for the dissipative soliton to m={}, however it is only"
"properly defined for m=1 or m=2.".format(m)
)
id = np.repeat(range(n), length * m)
dimensions = list(np.repeat(range(m), length)) * n
labels = list()
values = list()
ls_tau = sample_tau(n, kappa_3=kappa_3)
for i, tau in enumerate(ls_tau):
ds = velocity(tau=tau, kappa_3=kappa_3, seed=seed)
if classification:
labels.append(ds.label)
else:
labels.append(ds.tau)
values.append(ds.simulate(length, v0=np.zeros(m)).transpose().flatten())
time = np.stack([ds.delta_t * np.arange(length)] * n * m).flatten()
df = pd.DataFrame(
{
"id": id,
"time": time,
"value": np.stack(values).flatten(),
"dimension": dimensions,
}
)
y = pd.Series(labels)
y.index = range(n)
return df, y