qualia_core_db/specialized_libs/computational_economics/
markov.rs

1//! Discrete-time Markov chain core.
2//!
3//! Implements transition-matrix validation, stationary-distribution power
4//! iteration, deterministic seeded simulation, and mean first-passage-time
5//! solves. Part of the QualiaDB computational economics library (plan §5.3 /
6//! P3).
7//!
8//! # Allocation class: `HotZeroHeap`
9//!
10//! Every public kernel operates on caller-owned slices and fixed-capacity
11//! stack arrays (`[0.0f64; MAX_STATES]`). No `Vec`, `String`, or `Box` is
12//! constructed on the hot path. The only heap traffic is whatever the caller
13//! chose to allocate for the input/output buffers.
14//!
15//! # Assumptions
16//!
17//! `stationary_distribution_into` assumes the chain is **ergodic** (aperiodic
18//! and irreducible). Power iteration converges to the unique stationary
19//! distribution only under that assumption; otherwise the kernel returns
20//! `MarkovError::NonConverged` once the iteration budget is exhausted.
21//!
22//! `mean_first_passage_time_into` likewise assumes the target state is
23//! reachable from every other state (irreducibility); unreachable states
24//! produce divergent hitting times and surface as `NonConverged`.
25
26use super::error::{EconConvergence, EconStatus};
27
28/// Maximum number of states supported by the stack-array kernels.
29pub const MAX_STATES: usize = 32;
30
31/// Tolerance used when checking that a transition-matrix row sums to 1.0.
32const ROW_SUM_TOLERANCE: f64 = 1e-9;
33
34/// Markov-chain kernel error.
35#[derive(Debug, Clone, Copy, PartialEq, Eq)]
36pub enum MarkovError {
37    /// Bad dimensions, zero states, or `n > MAX_STATES`.
38    InvalidInput,
39    /// A transition matrix row does not sum to 1.0, or contains a negative /
40    /// non-finite entry.
41    InvalidTransitionMatrix,
42    /// A state index was out of range for the supplied matrix.
43    InvalidState,
44    /// A caller-owned output buffer was too small for the request.
45    BufferTooSmall,
46    /// A non-finite value (NaN / infinity) appeared during iteration.
47    NonFinite,
48    /// The iterative solver did not converge within the iteration budget.
49    NonConverged,
50}
51
52impl MarkovError {
53    /// Map to the ABI-stable `EconStatus` used by the shared error vocabulary.
54    pub fn to_status(self) -> EconStatus {
55        match self {
56            MarkovError::InvalidInput => EconStatus::InvalidInput,
57            MarkovError::InvalidTransitionMatrix => EconStatus::InvalidInput,
58            MarkovError::InvalidState => EconStatus::InvalidInput,
59            MarkovError::BufferTooSmall => EconStatus::BufferTooSmall,
60            MarkovError::NonFinite => EconStatus::NonFinite,
61            MarkovError::NonConverged => EconStatus::MaxIterations,
62        }
63    }
64}
65
66impl From<MarkovError> for EconStatus {
67    #[inline]
68    fn from(err: MarkovError) -> Self {
69        err.to_status()
70    }
71}
72
73// ---------------------------------------------------------------------------
74// Deterministic RNG — local SplitMix64 (reimplemented to avoid cross-module
75// coupling with `domains::financial::economics::stochastic`).
76// ---------------------------------------------------------------------------
77
78/// SplitMix64 bit-mixing PRNG. Deterministic for a given seed; used by the
79/// seeded simulation kernel so that identical seeds reproduce identical paths.
80#[derive(Debug, Clone, Copy)]
81struct SplitMix64 {
82    state: u64,
83}
84
85impl SplitMix64 {
86    #[inline]
87    fn new(seed: u64) -> Self {
88        Self { state: seed }
89    }
90
91    #[inline]
92    fn next_u64(&mut self) -> u64 {
93        self.state = self.state.wrapping_add(0x9E37_79B9_7F4A_7C15);
94        let mut z = self.state;
95        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
96        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
97        z ^ (z >> 31)
98    }
99
100    /// Uniform draw on the open interval (0, 1) with 53 random mantissa bits.
101    #[inline]
102    fn unit_open(&mut self) -> f64 {
103        let bits = self.next_u64() >> 11;
104        ((bits as f64) + 0.5) * (1.0 / ((1u64 << 53) as f64))
105    }
106}
107
108// ---------------------------------------------------------------------------
109// Helpers
110// ---------------------------------------------------------------------------
111
112/// Validate the dimension and buffer length for an `n`-state matrix.
113fn check_dim(n: usize) -> Result<(), MarkovError> {
114    if n == 0 || n > MAX_STATES {
115        return Err(MarkovError::InvalidInput);
116    }
117    Ok(())
118}
119
120/// Row-major index into a flat `n x n` transition matrix.
121#[inline]
122fn idx(row: usize, col: usize, n: usize) -> usize {
123    row * n + col
124}
125
126// ---------------------------------------------------------------------------
127// Public API
128// ---------------------------------------------------------------------------
129
130/// Validate that `p` is a row-stochastic `n x n` transition matrix.
131///
132/// Each row must sum to `1.0` within `ROW_SUM_TOLERANCE`, and every entry must
133/// be finite and non-negative. `p` must hold at least `n * n` elements.
134pub fn validate_transition_matrix(p: &[f64], n: usize) -> Result<(), MarkovError> {
135    check_dim(n)?;
136    if p.len() < n * n {
137        return Err(MarkovError::InvalidInput);
138    }
139    for row in 0..n {
140        let mut sum = 0.0f64;
141        for col in 0..n {
142            let entry = p[idx(row, col, n)];
143            if !entry.is_finite() {
144                return Err(MarkovError::InvalidTransitionMatrix);
145            }
146            if entry < 0.0 {
147                return Err(MarkovError::InvalidTransitionMatrix);
148            }
149            sum += entry;
150        }
151        if (sum - 1.0).abs() > ROW_SUM_TOLERANCE {
152            return Err(MarkovError::InvalidTransitionMatrix);
153        }
154    }
155    Ok(())
156}
157
158/// Look up `P[from][to]` with bounds checking.
159pub fn transition_probability(
160    p: &[f64],
161    n: usize,
162    from: usize,
163    to: usize,
164) -> Result<f64, MarkovError> {
165    check_dim(n)?;
166    if p.len() < n * n {
167        return Err(MarkovError::InvalidInput);
168    }
169    if from >= n || to >= n {
170        return Err(MarkovError::InvalidState);
171    }
172    Ok(p[idx(from, to, n)])
173}
174
175/// Expected holding time for `state`: `1 / (1 - P[state][state])`.
176///
177/// Returns `MarkovError::InvalidState` for an out-of-range state and
178/// `MarkovError::NonFinite` when the self-loop probability is `1.0` (the
179/// state is absorbing, so the holding time is infinite).
180pub fn expected_holding_time(p: &[f64], n: usize, state: usize) -> Result<f64, MarkovError> {
181    check_dim(n)?;
182    if p.len() < n * n {
183        return Err(MarkovError::InvalidInput);
184    }
185    if state >= n {
186        return Err(MarkovError::InvalidState);
187    }
188    let self_loop = p[idx(state, state, n)];
189    if !self_loop.is_finite() {
190        return Err(MarkovError::NonFinite);
191    }
192    let denom = 1.0 - self_loop;
193    if denom <= 0.0 {
194        return Err(MarkovError::NonFinite);
195    }
196    Ok(1.0 / denom)
197}
198
199/// Compute the stationary distribution of an ergodic Markov chain via power
200/// iteration: `pi_{t+1} = pi_t * P`.
201///
202/// Writes the stationary distribution into `out[..n]`. The initial guess is
203/// the uniform distribution. Convergence is declared when the infinity-norm
204/// of `pi_{t+1} - pi_t` falls below `tolerance`. Returns an `EconConvergence`
205/// report; `status` is `Converged` on success and `MaxIterations` (mapped from
206/// `MarkovError::NonConverged`) when the budget is exhausted.
207///
208/// Assumes the chain is ergodic (aperiodic + irreducible). Non-ergodic chains
209/// may fail to converge, in which case `NonConverged` is returned.
210pub fn stationary_distribution_into(
211    p: &[f64],
212    n: usize,
213    max_iterations: u32,
214    tolerance: f64,
215    out: &mut [f64],
216) -> Result<EconConvergence, MarkovError> {
217    check_dim(n)?;
218    if p.len() < n * n {
219        return Err(MarkovError::InvalidInput);
220    }
221    if out.len() < n {
222        return Err(MarkovError::BufferTooSmall);
223    }
224    if max_iterations == 0 || !tolerance.is_finite() || tolerance <= 0.0 {
225        return Err(MarkovError::InvalidInput);
226    }
227    validate_transition_matrix(p, n)?;
228
229    // Stack scratch arrays — HotZeroHeap.
230    let mut current = [0.0f64; MAX_STATES];
231    let mut next = [0.0f64; MAX_STATES];
232
233    let uniform = 1.0 / n as f64;
234    for i in 0..n {
235        current[i] = uniform;
236    }
237
238    let mut iter: u32 = 0;
239    let mut residual = f64::INFINITY;
240    let mut converged = false;
241
242    while iter < max_iterations {
243        // next[j] = sum_i current[i] * P[i][j]
244        for j in 0..n {
245            next[j] = 0.0;
246        }
247        for i in 0..n {
248            let ci = current[i];
249            if ci == 0.0 {
250                continue;
251            }
252            let base = i * n;
253            for j in 0..n {
254                next[j] += ci * p[base + j];
255            }
256        }
257
258        // Residual: infinity norm of (next - current).
259        let mut diff = 0.0f64;
260        for j in 0..n {
261            if !next[j].is_finite() {
262                return Err(MarkovError::NonFinite);
263            }
264            let d = (next[j] - current[j]).abs();
265            if d > diff {
266                diff = d;
267            }
268        }
269        residual = diff;
270
271        // Swap current and next.
272        for j in 0..n {
273            current[j] = next[j];
274        }
275
276        iter += 1;
277        if residual <= tolerance {
278            converged = true;
279            break;
280        }
281    }
282
283    if !converged {
284        // Still copy the best estimate into the caller buffer.
285        for j in 0..n {
286            out[j] = current[j];
287        }
288        return Err(MarkovError::NonConverged);
289    }
290
291    for j in 0..n {
292        out[j] = current[j];
293    }
294    Ok(EconConvergence::converged(iter, residual))
295}
296
297/// Deterministic seeded simulation of a Markov chain.
298///
299/// Starting from `initial_state`, draws `steps` successive states using a
300/// local SplitMix64 RNG seeded with `seed` and inverse-CDF sampling per step.
301/// Writes the state indices (including the initial state at index 0) into
302/// `out[..=steps]` — i.e. `out` must hold at least `steps + 1` elements.
303/// Returns the number of indices written (`steps + 1`).
304pub fn simulate_chain_into(
305    p: &[f64],
306    n: usize,
307    initial_state: usize,
308    steps: usize,
309    seed: u64,
310    out: &mut [usize],
311) -> Result<usize, MarkovError> {
312    check_dim(n)?;
313    if p.len() < n * n {
314        return Err(MarkovError::InvalidInput);
315    }
316    if initial_state >= n {
317        return Err(MarkovError::InvalidState);
318    }
319    if out.len() < steps + 1 {
320        return Err(MarkovError::BufferTooSmall);
321    }
322    validate_transition_matrix(p, n)?;
323
324    let mut rng = SplitMix64::new(seed);
325    let mut state = initial_state;
326    out[0] = state;
327
328    for step in 0..steps {
329        let u = rng.unit_open();
330        let base = state * n;
331        let mut acc = 0.0f64;
332        let mut next_state = n - 1; // fallback for floating-point tail
333        for col in 0..n {
334            acc += p[base + col];
335            if u < acc {
336                next_state = col;
337                break;
338            }
339        }
340        state = next_state;
341        out[step + 1] = state;
342    }
343    Ok(steps + 1)
344}
345
346/// Iterative solve for the mean first-passage time to `target`.
347///
348/// For an ergodic chain, `m_i` is the expected number of steps to first reach
349/// `target` starting from state `i`. `m_target = 0`. For `i != target`:
350///
351/// ```text
352/// m_i = 1 + sum_{j != target} P[i][j] * m_j
353/// ```
354///
355/// Solved by fixed-point iteration. Writes results into `out[..n]` with
356/// `out[target] = 0`. Returns an `EconConvergence` report.
357pub fn mean_first_passage_time_into(
358    p: &[f64],
359    n: usize,
360    target: usize,
361    max_iterations: u32,
362    tolerance: f64,
363    out: &mut [f64],
364) -> Result<EconConvergence, MarkovError> {
365    check_dim(n)?;
366    if p.len() < n * n {
367        return Err(MarkovError::InvalidInput);
368    }
369    if target >= n {
370        return Err(MarkovError::InvalidState);
371    }
372    if out.len() < n {
373        return Err(MarkovError::BufferTooSmall);
374    }
375    if max_iterations == 0 || !tolerance.is_finite() || tolerance <= 0.0 {
376        return Err(MarkovError::InvalidInput);
377    }
378    validate_transition_matrix(p, n)?;
379
380    let mut current = [0.0f64; MAX_STATES];
381    let mut next = [0.0f64; MAX_STATES];
382
383    // Initial guess: zero everywhere (target stays zero).
384    for i in 0..n {
385        current[i] = 0.0;
386    }
387
388    let mut iter: u32 = 0;
389    let mut residual = f64::INFINITY;
390    let mut converged = false;
391
392    while iter < max_iterations {
393        let mut diff = 0.0f64;
394        for i in 0..n {
395            if i == target {
396                next[i] = 0.0;
397                continue;
398            }
399            let base = i * n;
400            let mut s = 1.0f64;
401            for j in 0..n {
402                if j == target {
403                    continue;
404                }
405                s += p[base + j] * current[j];
406            }
407            if !s.is_finite() {
408                return Err(MarkovError::NonFinite);
409            }
410            next[i] = s;
411            let d = (s - current[i]).abs();
412            if d > diff {
413                diff = d;
414            }
415        }
416        residual = diff;
417
418        for i in 0..n {
419            current[i] = next[i];
420        }
421
422        iter += 1;
423        if residual <= tolerance {
424            converged = true;
425            break;
426        }
427    }
428
429    for i in 0..n {
430        out[i] = current[i];
431    }
432
433    if !converged {
434        return Err(MarkovError::NonConverged);
435    }
436    Ok(EconConvergence::converged(iter, residual))
437}
438
439// ---------------------------------------------------------------------------
440// Tests
441// ---------------------------------------------------------------------------
442
443#[cfg(test)]
444mod tests {
445    use super::*;
446
447    /// 2-state symmetric chain: stationary distribution is [0.5, 0.5].
448    #[test]
449    fn symmetric_two_state_stationary() {
450        // [[0.5, 0.5], [0.5, 0.5]]
451        let p = [0.5, 0.5, 0.5, 0.5];
452        let mut out = [0.0f64; MAX_STATES];
453        let conv = stationary_distribution_into(&p, 2, 10_000, 1e-12, &mut out).unwrap();
454        assert_eq!(conv.status, EconStatus::Converged);
455        assert!((out[0] - 0.5).abs() < 1e-9, "out[0]={}", out[0]);
456        assert!((out[1] - 0.5).abs() < 1e-9, "out[1]={}", out[1]);
457    }
458
459    /// 2-state asymmetric chain [[0.9,0.1],[0.5,0.5]] -> stationary [5/6, 1/6].
460    #[test]
461    fn asymmetric_two_state_stationary() {
462        let p = [0.9, 0.1, 0.5, 0.5];
463        let mut out = [0.0f64; MAX_STATES];
464        let conv = stationary_distribution_into(&p, 2, 100_000, 1e-12, &mut out).unwrap();
465        assert_eq!(conv.status, EconStatus::Converged);
466        assert!((out[0] - 5.0 / 6.0).abs() < 1e-6, "out[0]={}", out[0]);
467        assert!((out[1] - 1.0 / 6.0).abs() < 1e-6, "out[1]={}", out[1]);
468    }
469
470    /// Validation rejects a row that does not sum to 1.
471    #[test]
472    fn validation_rejects_row_not_summing_to_one() {
473        let p = [0.9, 0.2, 0.5, 0.5]; // row 0 sums to 1.1
474        let err = validate_transition_matrix(&p, 2).unwrap_err();
475        assert_eq!(err, MarkovError::InvalidTransitionMatrix);
476    }
477
478    /// Validation rejects a negative entry.
479    #[test]
480    fn validation_rejects_negative_entry() {
481        let p = [1.2, -0.2, 0.5, 0.5];
482        let err = validate_transition_matrix(&p, 2).unwrap_err();
483        assert_eq!(err, MarkovError::InvalidTransitionMatrix);
484    }
485
486    /// Validation rejects a NaN entry.
487    #[test]
488    fn validation_rejects_nan_entry() {
489        let p = [f64::NAN, 1.0, 0.5, 0.5];
490        let err = validate_transition_matrix(&p, 2).unwrap_err();
491        assert_eq!(err, MarkovError::InvalidTransitionMatrix);
492    }
493
494    /// Validation rejects an infinity entry.
495    #[test]
496    fn validation_rejects_infinity_entry() {
497        let p = [f64::INFINITY, 0.0, 0.5, 0.5];
498        let err = validate_transition_matrix(&p, 2).unwrap_err();
499        assert_eq!(err, MarkovError::InvalidTransitionMatrix);
500    }
501
502    /// Validation accepts a well-formed matrix.
503    #[test]
504    fn validation_accepts_valid_matrix() {
505        let p = [0.9, 0.1, 0.5, 0.5];
506        assert!(validate_transition_matrix(&p, 2).is_ok());
507    }
508
509    /// Validation rejects zero states and over-capacity states.
510    #[test]
511    fn validation_rejects_bad_dimensions() {
512        let p = [0.0; 4];
513        assert_eq!(
514            validate_transition_matrix(&p, 0).unwrap_err(),
515            MarkovError::InvalidInput
516        );
517        assert_eq!(
518            validate_transition_matrix(&p, MAX_STATES + 1).unwrap_err(),
519            MarkovError::InvalidInput
520        );
521    }
522
523    /// Validation rejects a buffer too small for the matrix.
524    #[test]
525    fn validation_rejects_undersized_slice() {
526        let p = [0.9, 0.1]; // only 2 elements for a 2x2
527        assert_eq!(
528            validate_transition_matrix(&p, 2).unwrap_err(),
529            MarkovError::InvalidInput
530        );
531    }
532
533    /// Same seed reproduces the same path.
534    #[test]
535    fn simulation_is_reproducible() {
536        let p = [0.9, 0.1, 0.5, 0.5];
537        let steps = 200;
538        let mut a = [0usize; 256];
539        let mut b = [0usize; 256];
540        let na = simulate_chain_into(&p, 2, 0, steps, 42, &mut a).unwrap();
541        let nb = simulate_chain_into(&p, 2, 0, steps, 42, &mut b).unwrap();
542        assert_eq!(na, nb);
543        assert_eq!(na, steps + 1);
544        for k in 0..na {
545            assert_eq!(a[k], b[k], "path divergence at {}", k);
546        }
547    }
548
549    /// Different seeds produce (almost certainly) different paths.
550    #[test]
551    fn different_seeds_diverge() {
552        let p = [0.5, 0.5, 0.5, 0.5];
553        let steps = 100;
554        let mut a = [0usize; 256];
555        let mut b = [0usize; 256];
556        simulate_chain_into(&p, 2, 0, steps, 1, &mut a).unwrap();
557        simulate_chain_into(&p, 2, 0, steps, 2, &mut b).unwrap();
558        // At least one step should differ.
559        let mut any_diff = false;
560        for k in 0..=steps {
561            if a[k] != b[k] {
562                any_diff = true;
563                break;
564            }
565        }
566        assert!(any_diff, "two different seeds produced identical paths");
567    }
568
569    /// Simulation respects transition probabilities statistically.
570    #[test]
571    fn simulation_respects_transition_probabilities() {
572        // [[0.9, 0.1], [0.5, 0.5]]
573        let p = [0.9, 0.1, 0.5, 0.5];
574        let steps = 200_000;
575        let mut path = vec![0usize; steps + 1];
576        simulate_chain_into(&p, 2, 0, steps, 12345, &mut path).unwrap();
577
578        // Count transitions out of state 0.
579        let mut from0_total = 0usize;
580        let mut from0_to1 = 0usize;
581        let mut from1_total = 0usize;
582        let mut from1_to0 = 0usize;
583        for k in 0..steps {
584            let s = path[k];
585            let ns = path[k + 1];
586            if s == 0 {
587                from0_total += 1;
588                if ns == 1 {
589                    from0_to1 += 1;
590                }
591            } else {
592                from1_total += 1;
593                if ns == 0 {
594                    from1_to0 += 1;
595                }
596            }
597        }
598        let p01 = from0_to1 as f64 / from0_total as f64;
599        let p10 = from1_to0 as f64 / from1_total as f64;
600        assert!((p01 - 0.1).abs() < 0.01, "empirical p01={}", p01);
601        assert!((p10 - 0.5).abs() < 0.01, "empirical p10={}", p10);
602    }
603
604    /// Holding time for a state with p=0.9 self-loop is 10.
605    #[test]
606    fn holding_time_self_loop() {
607        let p = [0.9, 0.1, 0.5, 0.5];
608        let h = expected_holding_time(&p, 2, 0).unwrap();
609        assert!((h - 10.0).abs() < 1e-9, "holding time={}", h);
610    }
611
612    /// Holding time for an absorbing state (self-loop = 1.0) errors.
613    #[test]
614    fn holding_time_absorbing_errors() {
615        let p = [1.0, 0.0, 0.5, 0.5];
616        let err = expected_holding_time(&p, 2, 0).unwrap_err();
617        assert_eq!(err, MarkovError::NonFinite);
618    }
619
620    /// Mean first-passage time on a small chain.
621    ///
622    /// Chain: 0 -> 1 -> 2 (deterministic), plus self-loops to make it ergodic.
623    /// Use [[0.5, 0.5, 0], [0, 0.5, 0.5], [0, 0, 1.0]] is not ergodic (2
624    /// absorbing). Instead use a 3-state chain where target=2.
625    #[test]
626    fn mean_first_passage_time_small_chain() {
627        // [[0.5, 0.5, 0.0],
628        //  [0.0, 0.5, 0.5],
629        //  [0.1, 0.1, 0.8]]
630        let p = [0.5, 0.5, 0.0, 0.0, 0.5, 0.5, 0.1, 0.1, 0.8];
631        let mut out = [0.0f64; MAX_STATES];
632        let conv = mean_first_passage_time_into(&p, 3, 2, 100_000, 1e-12, &mut out).unwrap();
633        assert_eq!(conv.status, EconStatus::Converged);
634        assert!(out[2].abs() < 1e-9, "target mfp must be 0, got {}", out[2]);
635        // m_0 and m_1 must be positive and finite.
636        assert!(out[0] > 0.0 && out[0].is_finite(), "m_0={}", out[0]);
637        assert!(out[1] > 0.0 && out[1].is_finite(), "m_1={}", out[1]);
638        // m_1 < m_0 since state 1 is closer to target 2.
639        assert!(out[1] < out[0], "m_1={} should be < m_0={}", out[1], out[0]);
640    }
641
642    /// Mean first-passage time on a 2-state chain with a closed-form check.
643    ///
644    /// [[0.9, 0.1], [0.5, 0.5]], target = 1.
645    /// m_1 = 0. m_0 = 1 + 0.9 * m_0 => m_0 = 1 / 0.1 = 10.
646    #[test]
647    fn mean_first_passage_time_two_state_closed_form() {
648        let p = [0.9, 0.1, 0.5, 0.5];
649        let mut out = [0.0f64; MAX_STATES];
650        let conv = mean_first_passage_time_into(&p, 2, 1, 100_000, 1e-12, &mut out).unwrap();
651        assert_eq!(conv.status, EconStatus::Converged);
652        assert!(out[1].abs() < 1e-9);
653        assert!((out[0] - 10.0).abs() < 1e-6, "m_0={}", out[0]);
654    }
655
656    /// Buffer-too-small error for stationary distribution.
657    #[test]
658    fn stationary_buffer_too_small() {
659        let p = [0.5, 0.5, 0.5, 0.5];
660        let mut out = [0.0f64; 1]; // need 2
661        let err = stationary_distribution_into(&p, 2, 100, 1e-9, &mut out).unwrap_err();
662        assert_eq!(err, MarkovError::BufferTooSmall);
663    }
664
665    /// Buffer-too-small error for simulation.
666    #[test]
667    fn simulation_buffer_too_small() {
668        let p = [0.5, 0.5, 0.5, 0.5];
669        let mut out = [0usize; 5]; // need steps+1 = 11
670        let err = simulate_chain_into(&p, 2, 0, 10, 1, &mut out).unwrap_err();
671        assert_eq!(err, MarkovError::BufferTooSmall);
672    }
673
674    /// Buffer-too-small error for mean first-passage time.
675    #[test]
676    fn mfp_buffer_too_small() {
677        let p = [0.5, 0.5, 0.5, 0.5];
678        let mut out = [0.0f64; 1];
679        let err = mean_first_passage_time_into(&p, 2, 0, 100, 1e-9, &mut out).unwrap_err();
680        assert_eq!(err, MarkovError::BufferTooSmall);
681    }
682
683    /// Invalid-state error for transition_probability.
684    #[test]
685    fn transition_probability_invalid_state() {
686        let p = [0.5, 0.5, 0.5, 0.5];
687        assert_eq!(
688            transition_probability(&p, 2, 2, 0).unwrap_err(),
689            MarkovError::InvalidState
690        );
691        assert_eq!(
692            transition_probability(&p, 2, 0, 2).unwrap_err(),
693            MarkovError::InvalidState
694        );
695    }
696
697    /// Invalid-state error for expected_holding_time.
698    #[test]
699    fn holding_time_invalid_state() {
700        let p = [0.5, 0.5, 0.5, 0.5];
701        assert_eq!(
702            expected_holding_time(&p, 2, 5).unwrap_err(),
703            MarkovError::InvalidState
704        );
705    }
706
707    /// Invalid-state error for simulation initial state.
708    #[test]
709    fn simulation_invalid_initial_state() {
710        let p = [0.5, 0.5, 0.5, 0.5];
711        let mut out = [0usize; 16];
712        assert_eq!(
713            simulate_chain_into(&p, 2, 5, 10, 1, &mut out).unwrap_err(),
714            MarkovError::InvalidState
715        );
716    }
717
718    /// Invalid-state error for mean first-passage time target.
719    #[test]
720    fn mfp_invalid_target() {
721        let p = [0.5, 0.5, 0.5, 0.5];
722        let mut out = [0.0f64; MAX_STATES];
723        assert_eq!(
724            mean_first_passage_time_into(&p, 2, 5, 100, 1e-9, &mut out).unwrap_err(),
725            MarkovError::InvalidState
726        );
727    }
728
729    /// transition_probability returns the correct entry.
730    #[test]
731    fn transition_probability_lookup() {
732        let p = [0.9, 0.1, 0.5, 0.5];
733        assert!((transition_probability(&p, 2, 0, 1).unwrap() - 0.1).abs() < 1e-12);
734        assert!((transition_probability(&p, 2, 1, 0).unwrap() - 0.5).abs() < 1e-12);
735    }
736
737    /// Stationary distribution of a 3-state chain sums to 1.
738    #[test]
739    fn stationary_distribution_normalizes() {
740        // Ergodic 3-state chain.
741        let p = [0.2, 0.6, 0.2, 0.3, 0.4, 0.3, 0.1, 0.2, 0.7];
742        let mut out = [0.0f64; MAX_STATES];
743        let conv = stationary_distribution_into(&p, 3, 100_000, 1e-12, &mut out).unwrap();
744        assert_eq!(conv.status, EconStatus::Converged);
745        let sum = out[0] + out[1] + out[2];
746        assert!((sum - 1.0).abs() < 1e-9, "sum={}", sum);
747        // Verify pi * P = pi.
748        for j in 0..3 {
749            let mut pip = 0.0;
750            for i in 0..3 {
751                pip += out[i] * p[i * 3 + j];
752            }
753            assert!((pip - out[j]).abs() < 1e-6, "pi*P[{}] != pi[{}]", j, j);
754        }
755    }
756
757    /// Non-convergence surfaces when the iteration budget is too small.
758    #[test]
759    fn stationary_non_converged_on_tiny_budget() {
760        // Near-identity but asymmetric so uniform start is not stationary.
761        let p = [0.999, 0.001, 0.1, 0.9];
762        let mut out = [0.0f64; MAX_STATES];
763        let err = stationary_distribution_into(&p, 2, 1, 1e-15, &mut out).unwrap_err();
764        assert_eq!(err, MarkovError::NonConverged);
765    }
766
767    /// Mean first-passage time non-convergence on a tiny budget.
768    #[test]
769    fn mfp_non_converged_on_tiny_budget() {
770        let p = [0.9, 0.1, 0.5, 0.5];
771        let mut out = [0.0f64; MAX_STATES];
772        let err = mean_first_passage_time_into(&p, 2, 1, 1, 1e-15, &mut out).unwrap_err();
773        assert_eq!(err, MarkovError::NonConverged);
774    }
775
776    /// Zero-step simulation writes only the initial state.
777    #[test]
778    fn simulation_zero_steps_writes_initial() {
779        let p = [0.5, 0.5, 0.5, 0.5];
780        let mut out = [0usize; 4];
781        let n = simulate_chain_into(&p, 2, 1, 0, 7, &mut out).unwrap();
782        assert_eq!(n, 1);
783        assert_eq!(out[0], 1);
784    }
785
786    /// Error-to-status mapping covers every variant.
787    #[test]
788    fn error_to_status_mapping() {
789        assert_eq!(
790            MarkovError::InvalidInput.to_status(),
791            EconStatus::InvalidInput
792        );
793        assert_eq!(
794            MarkovError::InvalidTransitionMatrix.to_status(),
795            EconStatus::InvalidInput
796        );
797        assert_eq!(
798            MarkovError::InvalidState.to_status(),
799            EconStatus::InvalidInput
800        );
801        assert_eq!(
802            MarkovError::BufferTooSmall.to_status(),
803            EconStatus::BufferTooSmall
804        );
805        assert_eq!(MarkovError::NonFinite.to_status(), EconStatus::NonFinite);
806        assert_eq!(
807            MarkovError::NonConverged.to_status(),
808            EconStatus::MaxIterations
809        );
810    }
811
812    /// `MAX_STATES` is 32.
813    #[test]
814    fn max_states_is_32() {
815        assert_eq!(MAX_STATES, 32);
816    }
817
818    /// A 4-state ring chain has a uniform stationary distribution.
819    #[test]
820    fn ring_chain_uniform_stationary() {
821        // 0->1->2->3->0, deterministic ring (periodic). Add small self-loops
822        // to make it aperiodic / ergodic.
823        let p = [
824            0.1, 0.9, 0.0, 0.0, //
825            0.0, 0.1, 0.9, 0.0, //
826            0.0, 0.0, 0.1, 0.9, //
827            0.9, 0.0, 0.0, 0.1, //
828        ];
829        let mut out = [0.0f64; MAX_STATES];
830        let conv = stationary_distribution_into(&p, 4, 100_000, 1e-12, &mut out).unwrap();
831        assert_eq!(conv.status, EconStatus::Converged);
832        for i in 0..4 {
833            assert!((out[i] - 0.25).abs() < 1e-6, "out[{}]={}", i, out[i]);
834        }
835    }
836}
qualia_core_db/specialized_libs/computational_economics/markov.rs

qualia_core_db/specialized_libs/computational_economics/
markov.rs
