changes

crates/domain/src/correlation/adjustment.rs

@@ -0,0 +1,76 @@
+use std::collections::BTreeMap;
+
+use super::{CorrelationStrategy, Family, PValue};
+
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Tested {
+    pub family: Family,
+    pub strategy: CorrelationStrategy,
+    pub p_value: PValue,
+}
+
+#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)]
+pub struct Adjustment {
+    false_discovery_rate: f64,
+}
+
+impl Adjustment {
+    pub fn controlling_false_discovery_at(false_discovery_rate: f64) -> Self {
+        Self {
+            false_discovery_rate,
+        }
+    }
+
+    pub fn holds_up_across(&self, tested: &[Tested]) -> Vec<bool> {
+        let mut groups: BTreeMap<(Family, CorrelationStrategy), Vec<usize>> = BTreeMap::new();
+        for (index, entry) in tested.iter().enumerate() {
+            groups
+                .entry((entry.family, entry.strategy))
+                .or_default()
+                .push(index);
+        }
+
+        let mut held = vec![false; tested.len()];
+
+        for places in groups.values() {
+            let p_values: Vec<PValue> = places.iter().map(|index| tested[*index].p_value).collect();
+
+            for (place, holds) in places.iter().zip(self.holds_up(&p_values)) {
+                held[*place] = holds;
+            }
+        }
+
+        held
+    }
+
+    pub fn holds_up(&self, p_values: &[PValue]) -> Vec<bool> {
+        let tested = p_values.len();
+
+        let mut by_size: Vec<usize> = (0..tested).collect();
+        by_size.sort_by(|left, right| p_values[*left].value().total_cmp(&p_values[*right].value()));
+
+        let largest_passing_rank = by_size
+            .iter()
+            .enumerate()
+            .filter(|(position, index)| {
+                p_values[**index].value() <= self.threshold_at(position + 1, tested)
+            })
+            .map(|(position, _)| position + 1)
+            .max();
+
+        let Some(rank) = largest_passing_rank else {
+            return vec![false; tested];
+        };
+
+        let mut held = vec![false; tested];
+        for index in &by_size[..rank] {
+            held[*index] = true;
+        }
+
+        held
+    }
+
+    fn threshold_at(&self, rank: usize, tested: usize) -> f64 {
+        rank as f64 / tested as f64 * self.false_discovery_rate
+    }
+}

crates/domain/src/correlation/agreement.rs

@@ -0,0 +1,29 @@
+use super::Coefficient;
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub struct Agreement {
+    agreeing: usize,
+    applicable: usize,
+}
+
+impl Agreement {
+    pub fn of(scores: &[Coefficient], applicable: usize) -> Self {
+        let positive = scores.iter().filter(|score| score.value() > 0.0).count();
+        let negative = scores.iter().filter(|score| score.value() < 0.0).count();
+
+        Self {
+            agreeing: positive.max(negative),
+            applicable,
+        }
+    }
+}
+
+impl Agreement {
+    pub fn agreeing(&self) -> usize {
+        self.agreeing
+    }
+
+    pub fn applicable(&self) -> usize {
+        self.applicable
+    }
+}

crates/domain/src/correlation/coefficient.rs

@@ -0,0 +1,29 @@
+use crate::errors::DomainError;
+
+const STRONGEST_NEGATIVE: f64 = -1.0;
+const STRONGEST_POSITIVE: f64 = 1.0;
+
+#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)]
+pub struct Coefficient(f64);
+
+impl Coefficient {
+    pub fn new(value: f64) -> Result<Self, DomainError> {
+        if !value.is_finite() {
+            return Err(DomainError::InvalidInput(
+                "a correlation coefficient must be a finite number".into(),
+            ));
+        }
+
+        if !(STRONGEST_NEGATIVE..=STRONGEST_POSITIVE).contains(&value) {
+            return Err(DomainError::InvalidInput(format!(
+                "a correlation coefficient must be between {STRONGEST_NEGATIVE} and {STRONGEST_POSITIVE}, got {value}"
+            )));
+        }
+
+        Ok(Self(value))
+    }
+
+    pub fn value(&self) -> f64 {
+        self.0
+    }
+}

crates/domain/src/correlation/correlation_input.rs

@@ -0,0 +1,48 @@
+use crate::activity::ActivityId;
+use crate::metric::MetricKind;
+
+use super::{Family, SeriesShape};
+
+const MOON_PHASE: &str = "moonPhase";
+const ACTIVITY: &str = "activity";
+const CYCLE_PROGRESS: &str = "cycleProgress";
+const TEMPERATURE: &str = "temperature";
+
+#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Hash)]
+pub enum CorrelationInput {
+    Metric(MetricKind),
+    MoonPhase,
+    CycleProgress,
+    Temperature,
+    Activity(ActivityId),
+}
+
+impl CorrelationInput {
+    pub fn name(&self) -> &'static str {
+        match self {
+            Self::Metric(kind) => kind.name(),
+            Self::MoonPhase => MOON_PHASE,
+            Self::CycleProgress => CYCLE_PROGRESS,
+            Self::Temperature => TEMPERATURE,
+            Self::Activity(_) => ACTIVITY,
+        }
+    }
+
+    pub fn family(&self) -> Family {
+        match self {
+            Self::Metric(_) | Self::MoonPhase | Self::CycleProgress | Self::Temperature => {
+                Family::Measurements
+            }
+            Self::Activity(_) => Family::Activities,
+        }
+    }
+
+    pub fn series(&self) -> SeriesShape {
+        match self {
+            Self::Metric(_) | Self::MoonPhase | Self::CycleProgress | Self::Temperature => {
+                SeriesShape::Continuous
+            }
+            Self::Activity(_) => SeriesShape::Presence,
+        }
+    }
+}

crates/domain/src/correlation/correlation_strategy.rs

@@ -0,0 +1,113 @@
+use super::PValue;
+use super::kendall::kendall_tau_b;
+use super::mean_difference::mean_difference;
+use super::pearson::pearson;
+use super::ranks::average_ranks;
+use super::welch::welch_standard_score;
+use super::{Coefficient, CorrelationInput, Observation, SeriesShape};
+
+const PEARSON: &str = "pearson";
+const SPEARMAN: &str = "spearman";
+const KENDALL: &str = "kendall";
+const MEAN_DIFFERENCE: &str = "meanDifference";
+const FEWEST_OBSERVATIONS: usize = 2;
+const FEWEST_FOR_A_P_VALUE: usize = 4;
+const SPEARMAN_VARIANCE_INFLATION: f64 = 1.06;
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
+pub enum CorrelationStrategy {
+    Pearson,
+    Spearman,
+    Kendall,
+    MeanDifference,
+}
+
+impl CorrelationStrategy {
+    pub const ALL: [CorrelationStrategy; 4] = [
+        Self::Pearson,
+        Self::Spearman,
+        Self::Kendall,
+        Self::MeanDifference,
+    ];
+
+    pub fn name(&self) -> &'static str {
+        match self {
+            Self::Pearson => PEARSON,
+            Self::Spearman => SPEARMAN,
+            Self::Kendall => KENDALL,
+            Self::MeanDifference => MEAN_DIFFERENCE,
+        }
+    }
+
+    pub fn scores(&self) -> SeriesShape {
+        match self {
+            Self::Pearson | Self::Spearman | Self::Kendall => SeriesShape::Continuous,
+            Self::MeanDifference => SeriesShape::Presence,
+        }
+    }
+
+    pub fn can_score(&self, input: &CorrelationInput) -> bool {
+        self.scores() == input.series()
+    }
+
+    pub fn significance(&self, observations: &[Observation]) -> Option<PValue> {
+        if observations.len() < FEWEST_FOR_A_P_VALUE {
+            return None;
+        }
+
+        self.standard_score(observations)
+            .map(PValue::from_standard_score)
+    }
+
+    fn standard_score(&self, observations: &[Observation]) -> Option<f64> {
+        let (values, moods) = series(observations);
+        let count = observations.len() as f64;
+
+        match self {
+            Self::Pearson => {
+                let scored = pearson(&values, &moods)?;
+                Some(fisher(scored, count, 1.0))
+            }
+            Self::Spearman => {
+                let scored = pearson(&average_ranks(&values), &average_ranks(&moods))?;
+                Some(fisher(scored, count, SPEARMAN_VARIANCE_INFLATION))
+            }
+            Self::Kendall => {
+                let scored = kendall_tau_b(&values, &moods)?;
+                Some(scored * (9.0 * count * (count - 1.0) / (2.0 * (2.0 * count + 5.0))).sqrt())
+            }
+            Self::MeanDifference => welch_standard_score(&values, &moods),
+        }
+    }
+
+    pub fn score(&self, observations: &[Observation]) -> Option<Coefficient> {
+        if observations.len() < FEWEST_OBSERVATIONS {
+            return None;
+        }
+
+        let (values, moods) = series(observations);
+
+        let scored = match self {
+            Self::Pearson => pearson(&values, &moods),
+            Self::Spearman => pearson(&average_ranks(&values), &average_ranks(&moods)),
+            Self::Kendall => kendall_tau_b(&values, &moods),
+            Self::MeanDifference => mean_difference(&values, &moods),
+        }?;
+
+        Coefficient::new(scored.clamp(-1.0, 1.0)).ok()
+    }
+}
+
+fn series(observations: &[Observation]) -> (Vec<f64>, Vec<f64>) {
+    let values = observations.iter().map(Observation::value).collect();
+    let moods = observations
+        .iter()
+        .map(|observation| observation.day_mood().value())
+        .collect();
+
+    (values, moods)
+}
+
+fn fisher(scored: f64, count: f64, variance_inflation: f64) -> f64 {
+    scored.atanh() * ((count - 3.0) / variance_inflation).sqrt()
+}

crates/domain/src/correlation/family.rs

@@ -0,0 +1,5 @@
+#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
+pub enum Family {
+    Measurements,
+    Activities,
+}

crates/domain/src/correlation/kendall.rs

@@ -0,0 +1,63 @@
+pub fn kendall_tau_b(left: &[f64], right: &[f64]) -> Option<f64> {
+    let mut concordant = 0i64;
+    let mut discordant = 0i64;
+
+    for first in 0..left.len() {
+        for second in (first + 1)..left.len() {
+            let agreement =
+                direction(left[first] - left[second]) * direction(right[first] - right[second]);
+
+            match agreement {
+                1 => concordant += 1,
+                -1 => discordant += 1,
+                _ => continue,
+            }
+        }
+    }
+
+    let pairs = pair_count(left.len());
+    let denominator = ((pairs - tied_pairs(left)) * (pairs - tied_pairs(right))).sqrt();
+
+    if denominator == 0.0 {
+        return None;
+    }
+
+    Some((concordant - discordant) as f64 / denominator)
+}
+
+fn direction(difference: f64) -> i64 {
+    if difference > 0.0 {
+        return 1;
+    }
+    if difference < 0.0 {
+        return -1;
+    }
+
+    0
+}
+
+fn pair_count(len: usize) -> f64 {
+    let len = len as f64;
+
+    len * (len - 1.0) / 2.0
+}
+
+fn tied_pairs(values: &[f64]) -> f64 {
+    let mut sorted = values.to_vec();
+    sorted.sort_by(|left, right| left.total_cmp(right));
+
+    let mut tied = 0.0;
+    let mut start = 0;
+
+    while start < sorted.len() {
+        let mut end = start;
+        while end + 1 < sorted.len() && sorted[end + 1] == sorted[start] {
+            end += 1;
+        }
+
+        tied += pair_count(end - start + 1);
+        start = end + 1;
+    }
+
+    tied
+}

crates/domain/src/correlation/mean_difference.rs

@@ -0,0 +1,33 @@
+use crate::entry::Mood;
+
+const PRESENT: f64 = 1.0;
+
+pub fn mean_difference(values: &[f64], moods: &[f64]) -> Option<f64> {
+    let present: Vec<f64> = pick(values, moods, |value| value == PRESENT);
+    let absent: Vec<f64> = pick(values, moods, |value| value != PRESENT);
+
+    let difference = mean(&present)? - mean(&absent)?;
+
+    Some(difference / mood_span())
+}
+
+fn pick(values: &[f64], moods: &[f64], wanted: impl Fn(f64) -> bool) -> Vec<f64> {
+    values
+        .iter()
+        .zip(moods)
+        .filter(|(value, _)| wanted(**value))
+        .map(|(_, mood)| *mood)
+        .collect()
+}
+
+fn mean(moods: &[f64]) -> Option<f64> {
+    if moods.is_empty() {
+        return None;
+    }
+
+    Some(moods.iter().sum::<f64>() / moods.len() as f64)
+}
+
+fn mood_span() -> f64 {
+    f64::from(Mood::Rad.value() - Mood::Awful.value())
+}

crates/domain/src/correlation/mod.rs

@@ -0,0 +1,25 @@
+mod adjustment;
+mod agreement;
+mod coefficient;
+mod correlation_input;
+mod correlation_strategy;
+mod family;
+mod kendall;
+mod mean_difference;
+mod normal;
+mod observation;
+mod p_value;
+mod pearson;
+mod ranks;
+mod series_shape;
+mod welch;
+
+pub use adjustment::{Adjustment, Tested};
+pub use agreement::Agreement;
+pub use coefficient::Coefficient;
+pub use correlation_input::CorrelationInput;
+pub use correlation_strategy::CorrelationStrategy;
+pub use family::Family;
+pub use observation::Observation;
+pub use p_value::PValue;
+pub use series_shape::SeriesShape;

crates/domain/src/correlation/normal.rs

@@ -0,0 +1,23 @@
+const ERF_A1: f64 = 0.254_829_592;
+const ERF_A2: f64 = -0.284_496_736;
+const ERF_A3: f64 = 1.421_413_741;
+const ERF_A4: f64 = -1.453_152_027;
+const ERF_A5: f64 = 1.061_405_429;
+const ERF_P: f64 = 0.327_591_1;
+
+pub fn two_sided_tail(standard_score: f64) -> f64 {
+    let beyond = 1.0 - error_function(standard_score.abs() / std::f64::consts::SQRT_2);
+
+    beyond.clamp(0.0, 1.0)
+}
+
+fn error_function(x: f64) -> f64 {
+    if x < 0.0 {
+        return -error_function(-x);
+    }
+
+    let t = 1.0 / (1.0 + ERF_P * x);
+    let series = t * (ERF_A1 + t * (ERF_A2 + t * (ERF_A3 + t * (ERF_A4 + t * ERF_A5))));
+
+    1.0 - series * (-x * x).exp()
+}

crates/domain/src/correlation/observation.rs

@@ -0,0 +1,21 @@
+use crate::entry::DayMood;
+
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Observation {
+    value: f64,
+    day_mood: DayMood,
+}
+
+impl Observation {
+    pub fn new(value: f64, day_mood: DayMood) -> Self {
+        Self { value, day_mood }
+    }
+
+    pub fn value(&self) -> f64 {
+        self.value
+    }
+
+    pub fn day_mood(&self) -> DayMood {
+        self.day_mood
+    }
+}

crates/domain/src/correlation/p_value.rs

@@ -0,0 +1,29 @@
+use crate::errors::DomainError;
+
+use super::normal::two_sided_tail;
+
+const CERTAIN: f64 = 1.0;
+const IMPOSSIBLE: f64 = 0.0;
+
+#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)]
+pub struct PValue(f64);
+
+impl PValue {
+    pub fn new(value: f64) -> Result<Self, DomainError> {
+        if !value.is_finite() || !(IMPOSSIBLE..=CERTAIN).contains(&value) {
+            return Err(DomainError::InvalidInput(format!(
+                "a p-value must be between {IMPOSSIBLE} and {CERTAIN}, got {value}"
+            )));
+        }
+
+        Ok(Self(value))
+    }
+
+    pub fn from_standard_score(standard_score: f64) -> Self {
+        Self(two_sided_tail(standard_score))
+    }
+
+    pub fn value(&self) -> f64 {
+        self.0
+    }
+}

crates/domain/src/correlation/pearson.rs

@@ -0,0 +1,20 @@
+pub fn pearson(left: &[f64], right: &[f64]) -> Option<f64> {
+    let count = left.len() as f64;
+    let left_mean = left.iter().sum::<f64>() / count;
+    let right_mean = right.iter().sum::<f64>() / count;
+
+    let covariance: f64 = left
+        .iter()
+        .zip(right)
+        .map(|(l, r)| (l - left_mean) * (r - right_mean))
+        .sum();
+
+    let left_variance: f64 = left.iter().map(|l| (l - left_mean).powi(2)).sum();
+    let right_variance: f64 = right.iter().map(|r| (r - right_mean).powi(2)).sum();
+
+    if left_variance == 0.0 || right_variance == 0.0 {
+        return None;
+    }
+
+    Some(covariance / (left_variance * right_variance).sqrt())
+}

crates/domain/src/correlation/ranks.rs

@@ -0,0 +1,30 @@
+pub fn average_ranks(values: &[f64]) -> Vec<f64> {
+    let mut order: Vec<usize> = (0..values.len()).collect();
+    order.sort_by(|left, right| values[*left].total_cmp(&values[*right]));
+
+    let mut ranks = vec![0.0; values.len()];
+    let mut start = 0;
+
+    while start < order.len() {
+        let mut end = start;
+        while end + 1 < order.len() && values[order[end + 1]] == values[order[start]] {
+            end += 1;
+        }
+
+        let shared = shared_rank(start, end);
+        for position in &order[start..=end] {
+            ranks[*position] = shared;
+        }
+
+        start = end + 1;
+    }
+
+    ranks
+}
+
+fn shared_rank(start: usize, end: usize) -> f64 {
+    let first = start as f64 + 1.0;
+    let last = end as f64 + 1.0;
+
+    (first + last) / 2.0
+}

crates/domain/src/correlation/series_shape.rs

@@ -0,0 +1,5 @@
+#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
+pub enum SeriesShape {
+    Continuous,
+    Presence,
+}

crates/domain/src/correlation/welch.rs

@@ -0,0 +1,51 @@
+const PRESENT: f64 = 1.0;
+const FEWEST_FOR_A_SPREAD: usize = 2;
+
+pub fn welch_standard_score(values: &[f64], moods: &[f64]) -> Option<f64> {
+    let present = group(values, moods, |value| value == PRESENT);
+    let absent = group(values, moods, |value| value != PRESENT);
+
+    if present.len() < FEWEST_FOR_A_SPREAD || absent.len() < FEWEST_FOR_A_SPREAD {
+        return None;
+    }
+
+    let spread = (variance(&present)? / present.len() as f64
+        + variance(&absent)? / absent.len() as f64)
+        .sqrt();
+
+    if spread == 0.0 {
+        return None;
+    }
+
+    Some((mean(&present)? - mean(&absent)?) / spread)
+}
+
+fn group(values: &[f64], moods: &[f64], wanted: impl Fn(f64) -> bool) -> Vec<f64> {
+    values
+        .iter()
+        .zip(moods)
+        .filter(|(value, _)| wanted(**value))
+        .map(|(_, mood)| *mood)
+        .collect()
+}
+
+fn mean(moods: &[f64]) -> Option<f64> {
+    if moods.is_empty() {
+        return None;
+    }
+
+    Some(moods.iter().sum::<f64>() / moods.len() as f64)
+}
+
+fn variance(moods: &[f64]) -> Option<f64> {
+    let average = mean(moods)?;
+    let degrees_of_freedom = moods.len() as f64 - 1.0;
+
+    Some(
+        moods
+            .iter()
+            .map(|mood| (mood - average).powi(2))
+            .sum::<f64>()
+            / degrees_of_freedom,
+    )
+}