picard/
density.rs

1// src/density.rs
2
3//! Density functions for ICA.
4//!
5//! The density function determines the non-linearity used in the ICA algorithm.
6//! Different densities are suited for different source distributions.
7
8use ndarray::{Array1, Array2};
9
10/// Trait for density functions used in ICA.
11///
12/// A density must provide methods for computing the log-likelihood,
13/// score function (derivative of log-likelihood), and score derivative.
14pub trait Density: Clone + Send + Sync {
15    /// Compute the log-likelihood for a 1D signal.
16    fn log_lik(&self, y: &Array1<f64>) -> Array1<f64>;
17
18    /// Compute the score function and its derivative for a 2D signal matrix.
19    ///
20    /// Returns `(score, score_derivative)` where both have the same shape as input.
21    fn score_and_der(&self, y: &Array2<f64>) -> (Array2<f64>, Array2<f64>);
22}
23
24/// Hyperbolic tangent density.
25///
26/// This is the default and most commonly used density. It works well for
27/// super-Gaussian sources (e.g., speech, sparse signals).
28///
29/// The log-likelihood is: `|y| + log(1 + exp(-2α|y|)) / α`
30#[derive(Clone, Debug)]
31pub struct Tanh {
32    /// Scaling parameter (default: 1.0).
33    pub alpha: f64,
34}
35
36impl Default for Tanh {
37    fn default() -> Self {
38        Self { alpha: 1.0 }
39    }
40}
41
42impl Tanh {
43    /// Create a new Tanh density with the given alpha parameter.
44    pub fn new(alpha: f64) -> Self {
45        Self { alpha }
46    }
47}
48
49impl Density for Tanh {
50    fn log_lik(&self, y: &Array1<f64>) -> Array1<f64> {
51        let alpha = self.alpha;
52        y.mapv(|v| {
53            let abs_y = v.abs();
54            abs_y + (1.0 + (-2.0 * alpha * abs_y).exp()).ln() / alpha
55        })
56    }
57
58    fn score_and_der(&self, y: &Array2<f64>) -> (Array2<f64>, Array2<f64>) {
59        let alpha = self.alpha;
60        let score = y.mapv(|v| (alpha * v).tanh());
61        let score_der = score.mapv(|s| alpha * (1.0 - s * s));
62        (score, score_der)
63    }
64}
65
66/// Exponential density.
67///
68/// Suited for super-Gaussian sources with heavy tails.
69///
70/// The log-likelihood is: `-exp(-αy²/2) / α`
71#[derive(Clone, Debug)]
72pub struct Exp {
73    /// Scaling parameter (default: 1.0).
74    pub alpha: f64,
75}
76
77impl Default for Exp {
78    fn default() -> Self {
79        Self { alpha: 1.0 }
80    }
81}
82
83impl Exp {
84    /// Create a new Exp density with the given alpha parameter.
85    pub fn new(alpha: f64) -> Self {
86        Self { alpha }
87    }
88}
89
90impl Density for Exp {
91    fn log_lik(&self, y: &Array1<f64>) -> Array1<f64> {
92        let a = self.alpha;
93        y.mapv(|v| -(-a * v * v / 2.0).exp() / a)
94    }
95
96    fn score_and_der(&self, y: &Array2<f64>) -> (Array2<f64>, Array2<f64>) {
97        let a = self.alpha;
98        let y_sq = y.mapv(|v| v * v);
99        let k = y_sq.mapv(|v| (-a / 2.0 * v).exp());
100        let score = y * &k;
101        let score_der = (1.0 - a * &y_sq) * k;
102        (score, score_der)
103    }
104}
105
106/// Cubic density.
107///
108/// Suited for sub-Gaussian sources (e.g., uniform distributions).
109///
110/// The log-likelihood is: `y⁴/4`
111#[derive(Clone, Debug, Default)]
112pub struct Cube;
113
114impl Cube {
115    /// Create a new Cube density.
116    pub fn new() -> Self {
117        Self
118    }
119}
120
121impl Density for Cube {
122    fn log_lik(&self, y: &Array1<f64>) -> Array1<f64> {
123        y.mapv(|v| v.powi(4) / 4.0)
124    }
125
126    fn score_and_der(&self, y: &Array2<f64>) -> (Array2<f64>, Array2<f64>) {
127        let score = y.mapv(|v| v.powi(3));
128        let score_der = y.mapv(|v| 3.0 * v * v);
129        (score, score_der)
130    }
131}
132
133/// Enumeration of built-in density types.
134///
135/// This allows specifying a density without type parameters.
136#[derive(Clone, Debug)]
137pub enum DensityType {
138    /// Hyperbolic tangent density.
139    Tanh(Tanh),
140    /// Exponential density.
141    Exp(Exp),
142    /// Cubic density.
143    Cube(Cube),
144}
145
146impl Default for DensityType {
147    fn default() -> Self {
148        DensityType::Tanh(Tanh::default())
149    }
150}
151
152impl DensityType {
153    /// Create a Tanh density with default parameters.
154    pub fn tanh() -> Self {
155        DensityType::Tanh(Tanh::default())
156    }
157
158    /// Create a Tanh density with custom alpha.
159    pub fn tanh_with_alpha(alpha: f64) -> Self {
160        DensityType::Tanh(Tanh::new(alpha))
161    }
162
163    /// Create an Exp density with default parameters.
164    pub fn exp() -> Self {
165        DensityType::Exp(Exp::default())
166    }
167
168    /// Create an Exp density with custom alpha.
169    pub fn exp_with_alpha(alpha: f64) -> Self {
170        DensityType::Exp(Exp::new(alpha))
171    }
172
173    /// Create a Cube density.
174    pub fn cube() -> Self {
175        DensityType::Cube(Cube::new())
176    }
177
178    /// Compute the log-likelihood.
179    pub fn log_lik(&self, y: &Array1<f64>) -> Array1<f64> {
180        match self {
181            DensityType::Tanh(d) => d.log_lik(y),
182            DensityType::Exp(d) => d.log_lik(y),
183            DensityType::Cube(d) => d.log_lik(y),
184        }
185    }
186
187    /// Compute the score function and its derivative.
188    pub fn score_and_der(&self, y: &Array2<f64>) -> (Array2<f64>, Array2<f64>) {
189        match self {
190            DensityType::Tanh(d) => d.score_and_der(y),
191            DensityType::Exp(d) => d.score_and_der(y),
192            DensityType::Cube(d) => d.score_and_der(y),
193        }
194    }
195}
picard/density.rs

picard/
density.rs
