gamlss-formula 0.2.1

Experimental optional formula and builder layer for typed GAMLSS models
Documentation 
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

use std::ops::Range;

use gamlss_core::{Penalty, SegmentPenalty};
use gamlss_spline::{PreparedCyclicDifferencePenalty, PreparedDifferencePenalty};

use crate::FittedTerm;

/// Formula-local segment penalty representation.
#[derive(Debug, Clone, PartialEq, Default)]
pub struct FormulaPenalty {
    spline_segments: Vec<SegmentPenalty<SegmentPenaltyKind>>,
}

#[derive(Debug, Clone, PartialEq)]
enum SegmentPenaltyKind {
    Difference(PreparedDifferencePenalty),
    Cyclic(PreparedCyclicDifferencePenalty),
}

impl Penalty for SegmentPenaltyKind {
    fn value(&self, beta: &[f64]) -> f64 {
        match self {
            Self::Difference(penalty) => penalty.value(beta),
            Self::Cyclic(penalty) => penalty.value(beta),
        }
    }

    fn add_gradient(&self, beta: &[f64], grad: &mut [f64]) {
        match self {
            Self::Difference(penalty) => penalty.add_gradient(beta, grad),
            Self::Cyclic(penalty) => penalty.add_gradient(beta, grad),
        }
    }
}

impl FormulaPenalty {
    pub(crate) fn add_spline(&mut self, range: Range<usize>, lambda: f64, order: usize) {
        self.spline_segments.push(SegmentPenalty::new(
            range,
            SegmentPenaltyKind::Difference(PreparedDifferencePenalty::new(lambda, order)),
        ));
    }

    pub(crate) fn add_cyclic_spline(&mut self, range: Range<usize>, lambda: f64, order: usize) {
        self.spline_segments.push(SegmentPenalty::new(
            range,
            SegmentPenaltyKind::Cyclic(PreparedCyclicDifferencePenalty::new(lambda, order)),
        ));
    }
}

impl Penalty for FormulaPenalty {
    fn value(&self, beta: &[f64]) -> f64 {
        self.spline_segments
            .iter()
            .map(|segment| segment.value(beta))
            .sum()
    }

    fn add_gradient(&self, beta: &[f64], grad: &mut [f64]) {
        debug_assert_eq!(beta.len(), grad.len());

        for segment in &self.spline_segments {
            segment.add_gradient(beta, grad);
        }
    }
}

pub(crate) fn prediction_penalty(terms: &[FittedTerm]) -> FormulaPenalty {
    let mut penalty = FormulaPenalty::default();
    for term in terms {
        match term {
            FittedTerm::PSpline {
                range,
                lambda,
                penalty_order,
                ..
            } => penalty.add_spline(range.clone(), *lambda, *penalty_order),
            FittedTerm::CyclicPSpline {
                range,
                lambda,
                penalty_order,
                ..
            } => penalty.add_cyclic_spline(range.clone(), *lambda, *penalty_order),
            _ => {}
        }
    }
    penalty
}