use crate::evaluator::*;
use crate::extractor::FeatureExtractor;
use crate::periodogram;
use crate::periodogram::{AverageNyquistFreq, NyquistFreq, PeriodogramPower, PeriodogramPowerFft};
use crate::statistics::Statistics;
use std::iter;
use std::marker::PhantomData;
#[derive(Clone, Debug)]
struct PeriodogramPeaks<T: Float> {
info: EvaluatorInfo,
peaks: usize,
names: Vec<String>,
phantom: PhantomData<T>,
}
impl<T> PeriodogramPeaks<T>
where
T: Float,
{
fn new(peaks: usize) -> Self {
assert!(peaks > 0, "Number of peaks should be at least one");
Self {
info: EvaluatorInfo {
size: 2 * peaks,
min_ts_length: 1,
t_required: true,
m_required: true,
w_required: false,
sorting_required: true,
},
peaks,
names: (0..peaks)
.flat_map(|i| vec![format!("period_{}", i), format!("period_s_to_n_{}", i)])
.collect(),
phantom: PhantomData {},
}
}
}
impl<T> Default for PeriodogramPeaks<T>
where
T: Float,
{
fn default() -> Self {
Self::new(1)
}
}
impl<T> FeatureEvaluator<T> for PeriodogramPeaks<T>
where
T: Float,
{
fn eval(&self, ts: &mut TimeSeries<T>) -> Result<Vec<T>, EvaluatorError> {
self.check_ts_length(ts)?;
Ok(ts
.m
.sample
.peak_indices_reverse_sorted()
.iter()
.flat_map(|&i| {
iter::once(Periodogram::period(ts.t.sample[i]))
.chain(iter::once(ts.m.signal_to_noise(ts.m.sample[i])))
})
.chain(iter::repeat(T::zero()))
.take(2 * self.peaks)
.collect())
}
fn get_info(&self) -> &EvaluatorInfo {
&self.info
}
fn get_names(&self) -> Vec<&str> {
self.names.iter().map(|name| name.as_str()).collect()
}
}
#[derive(Clone, Debug)]
pub struct Periodogram<T: Float> {
info: EvaluatorInfo,
resolution: f32,
max_freq_factor: f32,
nyquist: Box<dyn NyquistFreq<T>>,
feature_extractor: FeatureExtractor<T>,
feature_names: Vec<String>,
periodogram_algorithm: fn() -> Box<dyn PeriodogramPower<T>>,
}
impl<T> Periodogram<T>
where
T: Float,
{
pub fn new(peaks: usize) -> Self {
let peaks = PeriodogramPeaks::<T>::new(peaks);
let peak_names = peaks.names.clone();
let peaks_size_hint = peaks.size_hint();
let peaks_min_ts_length = peaks.min_ts_length();
Self {
info: EvaluatorInfo {
size: peaks_size_hint,
min_ts_length: usize::max(peaks_min_ts_length, 2),
t_required: true,
m_required: true,
w_required: false,
sorting_required: true,
},
resolution: 10.0,
max_freq_factor: 1.0,
nyquist: Box::new(AverageNyquistFreq),
feature_extractor: feat_extr!(peaks),
feature_names: peak_names,
periodogram_algorithm: || Box::new(PeriodogramPowerFft),
}
}
pub fn set_freq_resolution(&mut self, resolution: f32) -> &mut Self {
self.resolution = resolution;
self
}
pub fn set_max_freq_factor(&mut self, max_freq_factor: f32) -> &mut Self {
self.max_freq_factor = max_freq_factor;
self
}
pub fn set_nyquist(&mut self, nyquist: Box<dyn NyquistFreq<T>>) -> &mut Self {
self.nyquist = nyquist;
self
}
pub fn add_feature(&mut self, feature: Box<dyn FeatureEvaluator<T>>) -> &mut Self {
self.info.size += feature.size_hint();
self.feature_names.extend(
feature
.get_names()
.iter()
.map(|name| "periodogram_".to_owned() + name),
);
self.feature_extractor.add_feature(feature);
self
}
pub fn set_periodogram_algorithm(
&mut self,
periodogram_power: fn() -> Box<dyn PeriodogramPower<T>>,
) -> &mut Self {
self.periodogram_algorithm = periodogram_power;
self
}
pub fn init_thread_local_fft_plan(n: &[usize]) {
periodogram::Periodogram::<T>::init_thread_local_fft_plans(n);
}
fn periodogram(&self, ts: &mut TimeSeries<T>) -> periodogram::Periodogram<T> {
periodogram::Periodogram::from_t(
(self.periodogram_algorithm)(),
ts.t.sample,
self.resolution,
self.max_freq_factor,
&self.nyquist,
)
}
pub fn power(&self, ts: &mut TimeSeries<T>) -> Vec<T> {
self.periodogram(ts).power(ts)
}
pub fn freq_power(&self, ts: &mut TimeSeries<T>) -> (Vec<T>, Vec<T>) {
let p = self.periodogram(ts);
let power = p.power(ts);
let freq: Vec<_> = (0..power.len()).map(|i| p.freq(i)).collect();
(freq, power)
}
fn transform_ts(&self, ts: &mut TimeSeries<T>) -> Result<TMWVectors<T>, EvaluatorError> {
self.check_ts_length(ts)?;
let (freq, power) = self.freq_power(ts);
Ok(TMWVectors {
t: freq,
m: power,
w: None,
})
}
fn period(omega: T) -> T {
T::two() * T::PI() / omega
}
}
impl<T> Default for Periodogram<T>
where
T: Float,
{
fn default() -> Self {
Self::new(1)
}
}
impl<T> FeatureEvaluator<T> for Periodogram<T>
where
T: Float,
{
transformer_eval!();
}
#[cfg(test)]
#[allow(clippy::unreadable_literal)]
#[allow(clippy::excessive_precision)]
mod tests {
use super::*;
use crate::features::amplitude::Amplitude;
use crate::periodogram::{PeriodogramPowerDirect, QuantileNyquistFreq};
use crate::tests::*;
eval_info_test!(periodogram_info_1, {
let mut periodogram = Periodogram::default();
periodogram.set_periodogram_algorithm(|| Box::new(PeriodogramPowerDirect {}));
periodogram
});
eval_info_test!(periodogram_info_2, {
let mut periodogram = Periodogram::new(5);
periodogram.set_periodogram_algorithm(|| Box::new(PeriodogramPowerDirect {}));
periodogram
});
eval_info_test!(periodogram_info_3, {
let mut periodogram = Periodogram::default();
periodogram.add_feature(Box::new(Amplitude::default()));
periodogram.set_periodogram_algorithm(|| Box::new(PeriodogramPowerDirect {}));
periodogram
});
#[test]
fn periodogram_plateau() {
let fe = FeatureExtractor::new(vec![Box::new(Periodogram::default())]);
let x = linspace(0.0_f32, 1.0, 100);
let y = [0.0_f32; 100];
let mut ts = TimeSeries::new(&x[..], &y[..], None);
let desired = vec![0.0, 0.0];
let actual = fe.eval(&mut ts).unwrap();
assert_eq!(desired, actual);
}
#[test]
fn periodogram_evenly_sinus() {
let fe = FeatureExtractor::new(vec![Box::new(Periodogram::default())]);
let mut rng = StdRng::seed_from_u64(0);
let period = 0.17;
let x = linspace(0.0_f32, 1.0, 101);
let y: Vec<_> = x
.iter()
.map(|&x| {
3.0 * f32::sin(2.0 * std::f32::consts::PI / period * x + 0.5)
+ 4.0
+ 0.01 * rng.gen::<f32>()
})
.collect();
let mut ts = TimeSeries::new(&x[..], &y[..], None);
let desired = [period];
let actual = [fe.eval(&mut ts).unwrap()[0]];
all_close(&desired[..], &actual[..], 5e-3);
}
#[test]
fn periodogram_unevenly_sinus() {
let fe = FeatureExtractor::new(vec![Box::new(Periodogram::default())]);
let period = 0.17;
let mut rng = StdRng::seed_from_u64(0);
let mut x: Vec<f32> = (0..100).map(|_| rng.gen()).collect();
x[..].sort_unstable_by(|a, b| a.partial_cmp(b).unwrap());
let y: Vec<_> = x
.iter()
.map(|&x| 3.0 * f32::sin(2.0 * std::f32::consts::PI / period * x + 0.5) + 4.0)
.collect();
let mut ts = TimeSeries::new(&x[..], &y[..], None);
let desired = [period];
let actual = [fe.eval(&mut ts).unwrap()[0]];
all_close(&desired[..], &actual[..], 5e-3);
}
#[test]
fn periodogram_one_peak_vs_two_peaks() {
let fe = FeatureExtractor::new(vec![
Box::new(Periodogram::new(1)),
Box::new(Periodogram::new(2)),
]);
let period = 0.17;
let mut rng = StdRng::seed_from_u64(0);
let mut x: Vec<f32> = (0..100).map(|_| rng.gen()).collect();
x[..].sort_unstable_by(|a, b| a.partial_cmp(b).unwrap());
let y: Vec<_> = x
.iter()
.map(|&x| 3.0 * f32::sin(2.0 * std::f32::consts::PI / period * x + 0.5) + 4.0)
.collect();
let mut ts = TimeSeries::new(&x[..], &y[..], None);
let features = fe.eval(&mut ts).unwrap();
all_close(
&[features[0], features[1]],
&[features[2], features[3]],
1e-6,
);
}
#[test]
fn periodogram_unevenly_sinus_cosine() {
let fe = FeatureExtractor::new(vec![Box::new(Periodogram::new(2))]);
let period1 = 0.0753;
let period2 = 0.45;
let mut rng = StdRng::seed_from_u64(0);
let mut x: Vec<f32> = (0..1000).map(|_| rng.gen()).collect();
x[..].sort_unstable_by(|a, b| a.partial_cmp(b).unwrap());
let y: Vec<_> = x
.iter()
.map(|&x| {
3.0 * f32::sin(2.0 * std::f32::consts::PI / period1 * x + 0.5)
+ -5.0 * f32::cos(2.0 * std::f32::consts::PI / period2 * x + 0.5)
+ 4.0
})
.collect();
let mut ts = TimeSeries::new(&x[..], &y[..], None);
let desired = [period2, period1];
let features = fe.eval(&mut ts).unwrap();
let actual = [features[0], features[2]];
all_close(&desired[..], &actual[..], 1e-2);
assert!(features[1] > features[3]);
}
#[test]
fn periodogram_unevenly_sinus_cosine_noised() {
let fe = FeatureExtractor::new(vec![Box::new(Periodogram::new(2))]);
let period1 = 0.0753;
let period2 = 0.46;
let mut rng = StdRng::seed_from_u64(0);
let mut x: Vec<f32> = (0..1000).map(|_| rng.gen()).collect();
x[..].sort_unstable_by(|a, b| a.partial_cmp(b).unwrap());
let y: Vec<_> = x
.iter()
.map(|&x| {
3.0 * f32::sin(2.0 * std::f32::consts::PI / period1 * x + 0.5)
+ -5.0 * f32::cos(2.0 * std::f32::consts::PI / period2 * x + 0.5)
+ 10.0 * rng.gen::<f32>()
+ 4.0
})
.collect();
let mut ts = TimeSeries::new(&x[..], &y[..], None);
let desired = [period2, period1];
let features = fe.eval(&mut ts).unwrap();
let actual = [features[0], features[2]];
all_close(&desired[..], &actual[..], 1e-2);
assert!(features[1] > features[3]);
}
#[test]
fn periodogram_different_time_scales() {
let mut periodogram = Periodogram::new(2);
periodogram
.set_nyquist(Box::new(QuantileNyquistFreq { quantile: 0.05 }))
.set_freq_resolution(10.0)
.set_max_freq_factor(1.0)
.set_periodogram_algorithm(|| Box::new(PeriodogramPowerFft));
let fe = FeatureExtractor::new(vec![Box::new(periodogram)]);
let period1 = 0.01;
let period2 = 1.0;
let n = 100;
let mut x = linspace(0.0, 0.1, n);
x.append(&mut linspace(1.0, 10.0, n));
let y: Vec<_> = x
.iter()
.map(|&x| {
3.0 * f32::sin(2.0 * std::f32::consts::PI / period1 * x + 0.5)
+ -5.0 * f32::cos(2.0 * std::f32::consts::PI / period2 * x + 0.5)
+ 4.0
})
.collect();
let mut ts = TimeSeries::new(&x, &y, None);
let features = fe.eval(&mut ts).unwrap();
assert!(f32::abs(features[0] - period2) / period2 < 1.0 / n as f32);
assert!(f32::abs(features[2] - period1) / period1 < 1.0 / n as f32);
assert!(features[1] > features[3]);
}
}