use super::super::{
ComplexNumberSpace, Domain, DspVec, ErrorReason, NumberSpace, RealNumberSpace, ScalarResult,
ToSlice, Vector,
};
use crate::array_to_complex;
use arrayvec::ArrayVec;
use crate::multicore_support::*;
use crate::numbers::*;
use crate::simd_extensions::*;
#[repr(C)]
#[derive(Copy, Clone, PartialEq, Debug)]
pub struct Statistics<T> {
pub sum: T,
pub count: usize,
pub average: T,
pub rms: T,
pub min: T,
pub min_index: usize,
pub max: T,
pub max_index: usize,
}
pub const STATS_VEC_CAPACTIY: usize = 16;
pub type StatsVec<T> = ArrayVec<[T; STATS_VEC_CAPACTIY]>;
pub trait StatisticsOps<T> {
type Result;
fn statistics(&self) -> Self::Result;
}
pub trait StatisticsSplitOps<T> {
type Result;
fn statistics_split(&self, len: usize) -> ScalarResult<Self::Result>;
}
pub trait SumOps<T>: Sized
where
T: Sized,
{
fn sum(&self) -> T;
fn sum_sq(&self) -> T;
}
pub trait Stats<T>: Sized {
fn empty() -> Self;
fn empty_vec(len: usize) -> StatsVec<Self>;
fn invalid() -> Self;
fn merge(stats: &[Self]) -> Self;
fn merge_cols(stats: &[StatsVec<Self>]) -> StatsVec<Self>;
fn add(&mut self, elem: T, index: usize);
}
macro_rules! impl_common_stats {
() => {
fn merge_cols(stats: &[StatsVec<Self>]) -> StatsVec<Self> {
if stats.is_empty() {
return StatsVec::new();
}
let len = stats[0].len();
let mut results = StatsVec::new();
for i in 0..len {
let mut reordered = StatsVec::new();
for s in stats.iter() {
reordered.push(s[i]);
}
let merged = Statistics::merge(&reordered);
results.push(merged);
}
results
}
fn empty_vec(len: usize) -> StatsVec<Self> {
let mut results = StatsVec::new();
for _ in 0..len {
results.push(Statistics::empty());
}
results
}
};
}
impl<T> Stats<T> for Statistics<T>
where
T: RealNumber,
{
fn empty() -> Self {
Statistics {
sum: T::zero(),
count: 0,
average: T::zero(),
min: T::infinity(),
max: T::neg_infinity(),
rms: T::zero(),
min_index: 0,
max_index: 0,
}
}
fn invalid() -> Self {
Statistics {
sum: T::zero(),
count: 0,
average: T::nan(),
min: T::nan(),
max: T::nan(),
rms: T::nan(),
min_index: 0,
max_index: 0,
}
}
fn merge(stats: &[Statistics<T>]) -> Statistics<T> {
if stats.is_empty() {
return Statistics::<T>::invalid();
}
let mut sum = T::zero();
let mut max = stats[0].max;
let mut min = stats[0].min;
let mut max_index = stats[0].max_index;
let mut min_index = stats[0].min_index;
let mut sum_squared = T::zero();
let mut len = 0;
for stat in stats {
sum = sum + stat.sum;
len += stat.count;
sum_squared = sum_squared + stat.rms;
if stat.max > max {
max = stat.max;
max_index = stat.max_index;
} else if stat.min < min {
min = stat.min;
min_index = stat.min_index;
}
}
Statistics {
sum,
count: len,
average: sum / (T::from(len).unwrap()),
min,
max,
rms: (sum_squared / (T::from(len).unwrap())).sqrt(),
min_index,
max_index,
}
}
impl_common_stats!();
#[inline]
fn add(&mut self, elem: T, index: usize) {
self.sum = self.sum + elem;
self.count += 1;
self.rms = self.rms + elem * elem;
if elem > self.max {
self.max = elem;
self.max_index = index;
}
if elem < self.min {
self.min = elem;
self.min_index = index;
}
}
}
impl<T> Stats<Complex<T>> for Statistics<Complex<T>>
where
T: RealNumber,
{
fn empty() -> Self {
Statistics {
sum: Complex::<T>::new(T::zero(), T::zero()),
count: 0,
average: Complex::<T>::new(T::zero(), T::zero()),
min: Complex::<T>::new(T::infinity(), T::infinity()),
max: Complex::<T>::new(T::zero(), T::zero()),
rms: Complex::<T>::new(T::zero(), T::zero()),
min_index: 0,
max_index: 0,
}
}
fn invalid() -> Self {
Statistics {
sum: Complex::<T>::new(T::zero(), T::zero()),
count: 0,
average: Complex::<T>::new(T::nan(), T::nan()),
min: Complex::<T>::new(T::nan(), T::nan()),
max: Complex::<T>::new(T::nan(), T::nan()),
rms: Complex::<T>::new(T::nan(), T::nan()),
min_index: 0,
max_index: 0,
}
}
fn merge(stats: &[Statistics<Complex<T>>]) -> Statistics<Complex<T>> {
if stats.is_empty() {
return Statistics::<Complex<T>>::invalid();
}
let mut sum = Complex::<T>::new(T::zero(), T::zero());
let mut max = stats[0].max;
let mut min = stats[0].min;
let mut count = 0;
let mut max_index = stats[0].max_index;
let mut min_index = stats[0].min_index;
let mut max_norm = max.norm();
let mut min_norm = min.norm();
let mut sum_squared = Complex::<T>::new(T::zero(), T::zero());
for stat in stats {
sum = sum + stat.sum;
count += stat.count;
sum_squared = sum_squared + stat.rms;
if stat.max.norm() > max_norm {
max = stat.max;
max_norm = max.norm();
max_index = stat.max_index;
} else if stat.min.norm() < min_norm {
min = stat.min;
min_norm = min.norm();
min_index = stat.min_index;
}
}
Statistics {
sum,
count,
average: sum / (T::from(count).unwrap()),
min,
max,
rms: (sum_squared / (T::from(count).unwrap())).sqrt(),
min_index,
max_index,
}
}
impl_common_stats!();
#[inline]
fn add(&mut self, elem: Complex<T>, index: usize) {
self.sum = self.sum + elem;
self.count += 1;
self.rms = self.rms + elem * elem;
if elem.norm() > self.max.norm() {
self.max = elem;
self.max_index = index;
}
if elem.norm() < self.min.norm() {
self.min = elem;
self.min_index = index;
}
}
}
impl<S, T, N, D> StatisticsOps<T> for DspVec<S, T, N, D>
where
S: ToSlice<T>,
T: RealNumber,
N: RealNumberSpace,
D: Domain,
{
type Result = Statistics<T>;
fn statistics(&self) -> Statistics<T> {
let data_length = self.len();
let array = self.data.to_slice();
let chunks = Chunk::get_chunked_results(
Complexity::Small,
&self.multicore_settings,
&array[0..data_length],
1,
(),
|array, range, _arg| {
let mut stats = Statistics::empty();
let mut j = range.start;
for num in array {
stats.add(*num, j);
j += 1;
}
stats
},
);
Statistics::merge(&chunks[..])
}
}
impl<S, T, N, D> StatisticsSplitOps<T> for DspVec<S, T, N, D>
where
S: ToSlice<T>,
T: RealNumber,
N: RealNumberSpace,
D: Domain,
{
type Result = StatsVec<Statistics<T>>;
fn statistics_split(&self, len: usize) -> ScalarResult<StatsVec<Statistics<T>>> {
if len == 0 {
return Ok(StatsVec::new());
}
if len > STATS_VEC_CAPACTIY {
return Err(ErrorReason::InvalidArgumentLength);
}
let data_length = self.len();
let array = self.data.to_slice();
let chunks = Chunk::get_chunked_results(
Complexity::Small,
&self.multicore_settings,
&array[0..data_length],
1,
len,
|array, range, len| {
let mut results = Statistics::empty_vec(len);
let mut j = range.start;
for num in array {
let stats = &mut results[j % len];
stats.add(*num, j / len);
j += 1;
}
results
},
);
Ok(Statistics::merge_cols(&chunks[..]))
}
}
impl<S, T, N, D> DspVec<S, T, N, D>
where
S: ToSlice<T>,
T: RealNumber,
N: NumberSpace,
D: Domain,
{
fn sum_real<Reg: SimdGeneric<T>>(&self, _: RegType<Reg>) -> T {
let data_length = self.len();
let array = self.data.to_slice();
let partition = Reg::calc_data_alignment_reqs(&array[0..data_length]);
let sum = {
let chunks = Chunk::get_chunked_results(
Complexity::Small,
&self.multicore_settings,
partition.center(array),
Reg::LEN,
(),
move |array, _, _| {
let array = Reg::array_to_regs(array);
let mut sum = Reg::splat(T::zero());
for reg in array {
sum = sum + *reg;
}
sum
},
);
(&chunks[..])
.iter()
.map(|v| v.sum_real())
.fold(T::zero(), |a, b| a + b)
};
sum + partition
.edge_iter(array)
.fold(T::zero(), |sum, x| sum + *x)
}
fn sum_sq_real<Reg: SimdGeneric<T>>(&self, _: RegType<Reg>) -> T {
let data_length = self.len();
let array = self.data.to_slice();
let partition = Reg::calc_data_alignment_reqs(&array[0..data_length]);
let sum = {
let chunks = Chunk::get_chunked_results(
Complexity::Small,
&self.multicore_settings,
partition.center(array),
Reg::LEN,
(),
move |array, _, _| {
let array = Reg::array_to_regs(array);
let mut sum = Reg::splat(T::zero());
for reg in array {
sum = sum + *reg * *reg;
}
sum
},
);
(&chunks[..])
.iter()
.map(|v| v.sum_real())
.fold(T::zero(), |a, b| a + b)
};
sum + partition
.edge_iter(array)
.fold(T::zero(), |sum, x| sum + *x * *x)
}
fn sum_complex<Reg: SimdGeneric<T>>(&self, _: RegType<Reg>) -> Complex<T> {
let data_length = self.len();
let array = self.data.to_slice();
let partition = Reg::calc_data_alignment_reqs(&array[0..data_length]);
let sum = {
let chunks = Chunk::get_chunked_results(
Complexity::Small,
&self.multicore_settings,
partition.center(array),
Reg::LEN,
(),
move |array, _, _| {
let array = Reg::array_to_regs(array);
let mut sum = Reg::splat(T::zero());
for reg in array {
sum = sum + *reg;
}
sum
},
);
(&chunks[..])
.iter()
.map(|v| v.sum_complex())
.fold(Complex::<T>::new(T::zero(), T::zero()), |acc, x| acc + x)
};
sum + partition
.cedge_iter(array_to_complex(array))
.fold(Complex::<T>::zero(), |sum, x| sum + *x)
}
fn sum_sq_complex<Reg: SimdGeneric<T>>(&self, _: RegType<Reg>) -> Complex<T> {
let data_length = self.len();
let array = self.data.to_slice();
let partition = Reg::calc_data_alignment_reqs(&array[0..data_length]);
let sum = {
let chunks = Chunk::get_chunked_results(
Complexity::Small,
&self.multicore_settings,
partition.center(array),
Reg::LEN,
(),
move |array, _, _| {
let array = Reg::array_to_regs(array);
let mut sum = Reg::splat(T::zero());
for reg in array {
sum = sum + reg.mul_complex(*reg);
}
sum
},
);
(&chunks[..])
.iter()
.map(|v| v.sum_complex())
.fold(Complex::<T>::new(T::zero(), T::zero()), |acc, x| acc + x)
};
sum + partition
.cedge_iter(array_to_complex(array))
.fold(Complex::<T>::zero(), |sum, x| sum + *x * *x)
}
}
impl<S, T, N, D> SumOps<T> for DspVec<S, T, N, D>
where
S: ToSlice<T>,
T: RealNumber,
N: RealNumberSpace,
D: Domain,
{
fn sum(&self) -> T {
sel_reg!(self.sum_real::<T>())
}
fn sum_sq(&self) -> T {
sel_reg!(self.sum_sq_real::<T>())
}
}
impl<S, T, N, D> StatisticsOps<Complex<T>> for DspVec<S, T, N, D>
where
S: ToSlice<T>,
T: RealNumber,
N: ComplexNumberSpace,
D: Domain,
{
type Result = Statistics<Complex<T>>;
fn statistics(&self) -> Statistics<Complex<T>> {
let data_length = self.len();
let array = self.data.to_slice();
let chunks = Chunk::get_chunked_results(
Complexity::Small,
&self.multicore_settings,
&array[0..data_length],
2,
(),
|array, range, _arg| {
let mut stat = Statistics::<Complex<T>>::empty();
let mut j = range.start / 2;
let array = array_to_complex(array);
for num in array {
stat.add(*num, j);
j += 1;
}
stat
},
);
Statistics::merge(&chunks[..])
}
}
impl<S, T, N, D> StatisticsSplitOps<Complex<T>> for DspVec<S, T, N, D>
where
S: ToSlice<T>,
T: RealNumber,
N: ComplexNumberSpace,
D: Domain,
{
type Result = StatsVec<Statistics<Complex<T>>>;
fn statistics_split(&self, len: usize) -> ScalarResult<StatsVec<Statistics<Complex<T>>>> {
if len == 0 {
return Ok(StatsVec::new());
}
if len > STATS_VEC_CAPACTIY {
return Err(ErrorReason::InvalidArgumentLength);
}
let data_length = self.len();
let array = self.data.to_slice();
let chunks = Chunk::get_chunked_results(
Complexity::Small,
&self.multicore_settings,
&array[0..data_length],
2,
len,
|array, range, len| {
let mut results = Statistics::<Complex<T>>::empty_vec(len);
let mut j = range.start / 2;
let array = array_to_complex(array);
for num in array {
let stat = &mut results[j % len];
stat.add(*num, j / len);
j += 1;
}
results
},
);
Ok(Statistics::merge_cols(&chunks[..]))
}
}
impl<S, T, N, D> SumOps<Complex<T>> for DspVec<S, T, N, D>
where
S: ToSlice<T>,
T: RealNumber,
N: ComplexNumberSpace,
D: Domain,
{
fn sum(&self) -> Complex<T> {
sel_reg!(self.sum_complex::<T>())
}
fn sum_sq(&self) -> Complex<T> {
sel_reg!(self.sum_sq_complex::<T>())
}
}