strange_loop/
lipschitz_loop.rs

1//! Lipschitz-continuous strange loops implementation
2//!
3//! This module provides mathematical guarantees for convergence through
4//! Lipschitz continuity constraints and implements various loop topologies.
5
6use crate::error::{LoopError, Result};
7use crate::types::{Vector3D, Matrix3D};
8use crate::types::NalgebraVec3;
9use nalgebra::Matrix3;
10use serde::{Deserialize, Serialize};
11use std::collections::VecDeque;
12
13/// Parameters for Lipschitz-continuous strange loops
14#[derive(Clone, Debug, Serialize, Deserialize, PartialEq)]
15pub struct LipschitzParams {
16    /// Lipschitz constant L (must be < 1 for contraction)
17    pub lipschitz_constant: f64,
18    /// Convergence tolerance
19    pub tolerance: f64,
20    /// Maximum iterations before giving up
21    pub max_iterations: usize,
22    /// Enable adaptive Lipschitz estimation
23    pub adaptive_estimation: bool,
24    /// Damping factor for numerical stability
25    pub damping: f64,
26}
27
28impl Default for LipschitzParams {
29    fn default() -> Self {
30        Self {
31            lipschitz_constant: 0.9,
32            tolerance: 1e-12,
33            max_iterations: 10_000,
34            adaptive_estimation: true,
35            damping: 0.99,
36        }
37    }
38}
39
40impl LipschitzParams {
41    /// Create parameters optimized for fast convergence
42    pub fn fast_convergence() -> Self {
43        Self {
44            lipschitz_constant: 0.5,
45            tolerance: 1e-9,
46            max_iterations: 1_000,
47            adaptive_estimation: true,
48            damping: 0.95,
49        }
50    }
51
52    /// Create parameters for high-precision convergence
53    pub fn high_precision() -> Self {
54        Self {
55            lipschitz_constant: 0.8,
56            tolerance: 1e-15,
57            max_iterations: 100_000,
58            adaptive_estimation: true,
59            damping: 0.999,
60        }
61    }
62
63    /// Validate parameters for mathematical soundness
64    pub fn validate(&self) -> Result<()> {
65        if self.lipschitz_constant <= 0.0 || self.lipschitz_constant >= 1.0 {
66            return Err(LoopError::invalid_policy(
67                "Lipschitz constant must be in (0, 1) for convergence"
68            ));
69        }
70        if self.tolerance <= 0.0 {
71            return Err(LoopError::invalid_policy("Tolerance must be positive"));
72        }
73        if self.max_iterations == 0 {
74            return Err(LoopError::invalid_policy("Max iterations must be positive"));
75        }
76        if self.damping <= 0.0 || self.damping > 1.0 {
77            return Err(LoopError::invalid_policy("Damping must be in (0, 1]"));
78        }
79        Ok(())
80    }
81
82    /// Estimate convergence rate
83    pub fn convergence_rate(&self) -> f64 {
84        -self.lipschitz_constant.ln()
85    }
86
87    /// Estimate iterations needed for convergence
88    pub fn estimated_iterations(&self, initial_distance: f64) -> usize {
89        if initial_distance <= self.tolerance {
90            return 0;
91        }
92        let rate = self.convergence_rate();
93        if rate <= 0.0 {
94            return self.max_iterations;
95        }
96        ((initial_distance / self.tolerance).ln() / rate).ceil() as usize
97    }
98}
99
100/// Types of strange loop topologies
101#[derive(Clone, Debug, Serialize, Deserialize, PartialEq, Eq)]
102pub enum LoopTopology {
103    /// Simple fixed-point iteration: x_{n+1} = f(x_n)
104    FixedPoint,
105    /// Newton-Raphson with Lipschitz constraints
106    Newton,
107    /// Secant method variant
108    Secant,
109    /// Accelerated fixed-point (Anderson acceleration)
110    Accelerated,
111    /// Conjugate gradient style
112    ConjugateGradient,
113    /// Custom topology with user-defined function
114    Custom,
115}
116
117/// Convergence result information
118#[derive(Clone, Debug, Serialize, Deserialize)]
119pub struct ConvergenceResult {
120    /// Whether convergence was achieved
121    pub converged: bool,
122    /// Number of iterations taken
123    pub iterations: usize,
124    /// Final residual/error
125    pub final_residual: f64,
126    /// Estimated Lipschitz constant
127    pub estimated_lipschitz: f64,
128    /// Convergence rate achieved
129    pub convergence_rate: f64,
130    /// Time to convergence in nanoseconds
131    pub convergence_time_ns: u128,
132    /// Trajectory of residuals
133    pub residual_history: Vec<f64>,
134}
135
136/// Lipschitz-continuous strange loop implementation
137pub struct LipschitzLoop {
138    params: LipschitzParams,
139    topology: LoopTopology,
140    state_history: VecDeque<NalgebraVec3>,
141    residual_history: VecDeque<f64>,
142    estimated_lipschitz: f64,
143    iteration_count: usize,
144}
145
146impl LipschitzLoop {
147    /// Create a new Lipschitz loop
148    pub fn new(params: LipschitzParams, topology: LoopTopology) -> Result<Self> {
149        params.validate()?;
150
151        Ok(Self {
152            params,
153            topology,
154            state_history: VecDeque::with_capacity(1000),
155            residual_history: VecDeque::with_capacity(1000),
156            estimated_lipschitz: 0.0,
157            iteration_count: 0,
158        })
159    }
160
161    /// Execute the strange loop until convergence
162    pub fn execute<F>(&mut self, mut f: F, initial_state: NalgebraVec3) -> Result<ConvergenceResult>
163    where
164        F: FnMut(NalgebraVec3) -> NalgebraVec3,
165    {
166        let start_time = std::time::Instant::now();
167        let mut current_state = initial_state;
168        let mut previous_state = initial_state;
169        self.state_history.clear();
170        self.residual_history.clear();
171        self.iteration_count = 0;
172
173        // Store initial state
174        self.state_history.push_back(current_state);
175
176        for iteration in 0..self.params.max_iterations {
177            self.iteration_count = iteration;
178
179            // Apply the loop function based on topology
180            let next_state = match self.topology {
181                LoopTopology::FixedPoint => self.fixed_point_step(&mut f, current_state)?,
182                LoopTopology::Newton => self.newton_step(&mut f, current_state, previous_state)?,
183                LoopTopology::Secant => self.secant_step(&mut f, current_state, previous_state)?,
184                LoopTopology::Accelerated => self.accelerated_step(&mut f, current_state)?,
185                LoopTopology::ConjugateGradient => self.conjugate_gradient_step(&mut f, current_state, previous_state)?,
186                LoopTopology::Custom => f(current_state), // Direct application for custom
187            };
188
189            // Calculate residual
190            let residual = (next_state - current_state).norm();
191            self.residual_history.push_back(residual);
192
193            // Update Lipschitz estimate if enabled
194            if self.params.adaptive_estimation && iteration > 0 {
195                self.update_lipschitz_estimate(current_state, next_state, previous_state);
196            }
197
198            // Check convergence
199            if residual < self.params.tolerance {
200                return Ok(ConvergenceResult {
201                    converged: true,
202                    iterations: iteration + 1,
203                    final_residual: residual,
204                    estimated_lipschitz: self.estimated_lipschitz,
205                    convergence_rate: self.calculate_convergence_rate(),
206                    convergence_time_ns: start_time.elapsed().as_nanos(),
207                    residual_history: self.residual_history.iter().copied().collect(),
208                });
209            }
210
211            // Check Lipschitz constraint violation
212            if self.estimated_lipschitz > self.params.lipschitz_constant && iteration > 2 {
213                return Err(LoopError::lipschitz_violation(
214                    self.estimated_lipschitz,
215                    self.params.lipschitz_constant,
216                ));
217            }
218
219            // Update state
220            previous_state = current_state;
221            current_state = next_state;
222            self.state_history.push_back(current_state);
223
224            // Maintain history size
225            if self.state_history.len() > 1000 {
226                self.state_history.pop_front();
227            }
228            if self.residual_history.len() > 1000 {
229                self.residual_history.pop_front();
230            }
231        }
232
233        // Did not converge
234        Err(LoopError::convergence_failure(self.params.max_iterations))
235    }
236
237    /// Fixed-point iteration step
238    fn fixed_point_step<F>(&self, f: &mut F, current: NalgebraVec3) -> Result<NalgebraVec3>
239    where
240        F: FnMut(NalgebraVec3) -> NalgebraVec3,
241    {
242        let new_state = f(current);
243        Ok(current + self.params.damping * (new_state - current))
244    }
245
246    /// Newton-Raphson style step with numerical differentiation
247    fn newton_step<F>(&self, f: &mut F, current: NalgebraVec3, _previous: NalgebraVec3) -> Result<NalgebraVec3>
248    where
249        F: FnMut(NalgebraVec3) -> NalgebraVec3,
250    {
251        let residual = f(current) - current;
252
253        // Numerical Jacobian estimation
254        let h = 1e-8;
255        let mut jacobian = Matrix3::zeros();
256
257        for i in 0..3 {
258            let mut perturbed = current;
259            perturbed[i] += h;
260            let f_perturbed = f(perturbed) - perturbed;
261
262            for j in 0..3 {
263                jacobian[(j, i)] = (f_perturbed[j] - residual[j]) / h;
264            }
265        }
266
267        // Add identity to make it (J - I)
268        for i in 0..3 {
269            jacobian[(i, i)] -= 1.0;
270        }
271
272        // Solve linear system: (J - I) * delta = -residual
273        let delta = match jacobian.lu().solve(&(-residual)) {
274            Some(solution) => solution,
275            None => {
276                // Fallback to steepest descent
277                -self.params.damping * residual
278            }
279        };
280
281        Ok(current + delta)
282    }
283
284    /// Secant method step
285    fn secant_step<F>(&self, f: &mut F, current: NalgebraVec3, previous: NalgebraVec3) -> Result<NalgebraVec3>
286    where
287        F: FnMut(NalgebraVec3) -> NalgebraVec3,
288    {
289        let f_current = f(current);
290        let f_previous = f(previous);
291
292        // Secant approximation
293        let denominator = f_current - f_previous;
294        if denominator.norm() < 1e-12 {
295            // Fallback to fixed-point
296            return self.fixed_point_step(f, current);
297        }
298
299        let secant_direction = (current - previous).component_div(&denominator);
300        let residual = f_current - current;
301
302        Ok(current - secant_direction.component_mul(&residual) * self.params.damping)
303    }
304
305    /// Anderson acceleration step
306    fn accelerated_step<F>(&self, f: &mut F, current: NalgebraVec3) -> Result<NalgebraVec3>
307    where
308        F: FnMut(NalgebraVec3) -> NalgebraVec3,
309    {
310        if self.state_history.len() < 3 {
311            return self.fixed_point_step(f, current);
312        }
313
314        let f_current = f(current);
315
316        // Simple Anderson acceleration with last two iterates
317        let m = self.state_history.len().min(3);
318        if m >= 2 {
319            let x_prev = self.state_history[self.state_history.len() - 2];
320            let f_prev = f(x_prev);
321
322            let df = f_current - f_prev;
323            let dx = current - x_prev;
324
325            if df.norm() > 1e-12 {
326                let alpha = dx.dot(&df) / df.norm_squared();
327                let alpha_clamped = alpha.clamp(0.0, 1.0);
328
329                return Ok(current + alpha_clamped * (f_current - current));
330            }
331        }
332
333        self.fixed_point_step(f, current)
334    }
335
336    /// Conjugate gradient style step
337    fn conjugate_gradient_step<F>(&self, f: &mut F, current: NalgebraVec3, previous: NalgebraVec3) -> Result<NalgebraVec3>
338    where
339        F: FnMut(NalgebraVec3) -> NalgebraVec3,
340    {
341        let gradient = f(current) - current;
342
343        if self.state_history.len() < 2 {
344            return Ok(current - self.params.damping * gradient);
345        }
346
347        // Previous gradient
348        let prev_gradient = f(previous) - previous;
349        let direction = if prev_gradient.norm() > 1e-12 {
350            let beta = gradient.norm_squared() / prev_gradient.norm_squared();
351            let prev_direction = current - previous;
352            gradient + beta * prev_direction
353        } else {
354            gradient
355        };
356
357        Ok(current - self.params.damping * direction)
358    }
359
360    /// Update Lipschitz constant estimate
361    fn update_lipschitz_estimate(&mut self, current: NalgebraVec3, next: NalgebraVec3, previous: NalgebraVec3) {
362        if (current - previous).norm() < 1e-12 {
363            return;
364        }
365
366        let output_distance = (next - current).norm();
367        let input_distance = (current - previous).norm();
368
369        let local_lipschitz = output_distance / input_distance;
370
371        // Exponentially weighted moving average
372        const ALPHA: f64 = 0.1;
373        if self.estimated_lipschitz == 0.0 {
374            self.estimated_lipschitz = local_lipschitz;
375        } else {
376            self.estimated_lipschitz = (1.0 - ALPHA) * self.estimated_lipschitz + ALPHA * local_lipschitz;
377        }
378    }
379
380    /// Calculate convergence rate from residual history
381    fn calculate_convergence_rate(&self) -> f64 {
382        if self.residual_history.len() < 2 {
383            return 0.0;
384        }
385
386        let n = self.residual_history.len();
387        let initial_residual = self.residual_history[0];
388        let final_residual = self.residual_history[n - 1];
389
390        if initial_residual <= 0.0 || final_residual <= 0.0 {
391            return 0.0;
392        }
393
394        -(final_residual / initial_residual).ln() / (n as f64)
395    }
396
397    /// Get current parameters
398    pub fn params(&self) -> &LipschitzParams {
399        &self.params
400    }
401
402    /// Update parameters
403    pub fn update_params(&mut self, params: LipschitzParams) -> Result<()> {
404        params.validate()?;
405        self.params = params;
406        Ok(())
407    }
408
409    /// Get current topology
410    pub fn topology(&self) -> &LoopTopology {
411        &self.topology
412    }
413
414    /// Change topology
415    pub fn set_topology(&mut self, topology: LoopTopology) {
416        self.topology = topology;
417    }
418
419    /// Get state history
420    pub fn state_history(&self) -> &VecDeque<NalgebraVec3> {
421        &self.state_history
422    }
423
424    /// Get residual history
425    pub fn residual_history(&self) -> &VecDeque<f64> {
426        &self.residual_history
427    }
428
429    /// Get current Lipschitz estimate
430    pub fn estimated_lipschitz(&self) -> f64 {
431        self.estimated_lipschitz
432    }
433
434    /// Reset the loop state
435    pub fn reset(&mut self) {
436        self.state_history.clear();
437        self.residual_history.clear();
438        self.estimated_lipschitz = 0.0;
439        self.iteration_count = 0;
440    }
441
442    /// Execute with custom convergence criteria
443    pub fn execute_with_criteria<F, C>(
444        &mut self,
445        mut f: F,
446        initial_state: NalgebraVec3,
447        mut convergence_check: C,
448    ) -> Result<ConvergenceResult>
449    where
450        F: FnMut(NalgebraVec3) -> NalgebraVec3,
451        C: FnMut(NalgebraVec3, NalgebraVec3, usize) -> bool,
452    {
453        let start_time = std::time::Instant::now();
454        let mut current_state = initial_state;
455        let mut previous_state = initial_state;
456        self.state_history.clear();
457        self.residual_history.clear();
458        self.iteration_count = 0;
459
460        self.state_history.push_back(current_state);
461
462        for iteration in 0..self.params.max_iterations {
463            self.iteration_count = iteration;
464
465            let next_state = match self.topology {
466                LoopTopology::FixedPoint => self.fixed_point_step(&mut f, current_state)?,
467                LoopTopology::Newton => self.newton_step(&mut f, current_state, previous_state)?,
468                LoopTopology::Secant => self.secant_step(&mut f, current_state, previous_state)?,
469                LoopTopology::Accelerated => self.accelerated_step(&mut f, current_state)?,
470                LoopTopology::ConjugateGradient => self.conjugate_gradient_step(&mut f, current_state, previous_state)?,
471                LoopTopology::Custom => f(current_state),
472            };
473
474            let residual = (next_state - current_state).norm();
475            self.residual_history.push_back(residual);
476
477            // Custom convergence check
478            if convergence_check(current_state, next_state, iteration) {
479                return Ok(ConvergenceResult {
480                    converged: true,
481                    iterations: iteration + 1,
482                    final_residual: residual,
483                    estimated_lipschitz: self.estimated_lipschitz,
484                    convergence_rate: self.calculate_convergence_rate(),
485                    convergence_time_ns: start_time.elapsed().as_nanos(),
486                    residual_history: self.residual_history.iter().copied().collect(),
487                });
488            }
489
490            if self.params.adaptive_estimation && iteration > 0 {
491                self.update_lipschitz_estimate(current_state, next_state, previous_state);
492            }
493
494            previous_state = current_state;
495            current_state = next_state;
496            self.state_history.push_back(current_state);
497
498            if self.state_history.len() > 1000 {
499                self.state_history.pop_front();
500            }
501            if self.residual_history.len() > 1000 {
502                self.residual_history.pop_front();
503            }
504        }
505
506        Err(LoopError::convergence_failure(self.params.max_iterations))
507    }
508
509    /// Analyze stability around a fixed point
510    pub fn analyze_stability<F>(&self, f: F, fixed_point: NalgebraVec3) -> Result<StabilityAnalysis>
511    where
512        F: Fn(NalgebraVec3) -> NalgebraVec3,
513    {
514        // Numerical Jacobian at fixed point
515        let h = 1e-8;
516        let mut jacobian = Matrix3::zeros();
517
518        for i in 0..3 {
519            let mut perturbed = fixed_point;
520            perturbed[i] += h;
521            let f_perturbed = f(perturbed);
522            let f_fixed = f(fixed_point);
523
524            for j in 0..3 {
525                jacobian[(j, i)] = (f_perturbed[j] - f_fixed[j]) / h;
526            }
527        }
528
529        // Calculate eigenvalues (simplified - real eigenvalues only)
530        let eigenvalues = jacobian.eigenvalues().unwrap_or_default();
531        let max_eigenvalue = eigenvalues.iter()
532            .map(|e| e.abs())
533            .fold(0.0f64, f64::max);
534
535        let stability = if max_eigenvalue < 1.0 {
536            StabilityType::Stable
537        } else if max_eigenvalue == 1.0 {
538            StabilityType::Marginal
539        } else {
540            StabilityType::Unstable
541        };
542
543        Ok(StabilityAnalysis {
544            stability,
545            max_eigenvalue,
546            spectral_radius: max_eigenvalue,
547            jacobian,
548            eigenvalues: eigenvalues.iter().map(|e| e.abs()).collect(),
549        })
550    }
551}
552
553/// Stability analysis result
554#[derive(Clone, Debug)]
555pub struct StabilityAnalysis {
556    /// Stability classification
557    pub stability: StabilityType,
558    /// Maximum eigenvalue magnitude
559    pub max_eigenvalue: f64,
560    /// Spectral radius
561    pub spectral_radius: f64,
562    /// Jacobian matrix at fixed point
563    pub jacobian: Matrix3D,
564    /// Eigenvalue magnitudes
565    pub eigenvalues: Vec<f64>,
566}
567
568/// Types of stability
569#[derive(Clone, Debug, PartialEq, Eq)]
570pub enum StabilityType {
571    /// All eigenvalues have magnitude < 1
572    Stable,
573    /// At least one eigenvalue has magnitude = 1
574    Marginal,
575    /// At least one eigenvalue has magnitude > 1
576    Unstable,
577}
578
579/// Factory for creating common loop functions
580pub struct LoopFunctionFactory;
581
582impl LoopFunctionFactory {
583    /// Create a simple scalar function mapped to 3D
584    pub fn scalar_function(target: f64, step_size: f64) -> impl Fn(NalgebraVec3) -> NalgebraVec3 {
585        move |x: NalgebraVec3| {
586            let gradient = NalgebraVec3::new(
587                x[0] - target,
588                x[1] - target,
589                x[2] - target,
590            );
591            x - step_size * gradient
592        }
593    }
594
595    /// Create a quadratic bowl function
596    pub fn quadratic_bowl(center: NalgebraVec3, curvature: f64) -> impl Fn(NalgebraVec3) -> NalgebraVec3 {
597        move |x: NalgebraVec3| {
598            let gradient = curvature * (x - center);
599            x - 0.1 * gradient
600        }
601    }
602
603    /// Create a Rosenbrock-like function
604    pub fn rosenbrock_like(a: f64, b: f64) -> impl Fn(NalgebraVec3) -> NalgebraVec3 {
605        move |x: NalgebraVec3| {
606            let grad_x = -2.0 * a * (1.0 - x[0]) - 4.0 * b * x[0] * (x[1] - x[0] * x[0]);
607            let grad_y = 2.0 * b * (x[1] - x[0] * x[0]);
608            let grad_z = 2.0 * (x[2] - 1.0);
609
610            NalgebraVec3::new(
611                x[0] - 0.001 * grad_x,
612                x[1] - 0.001 * grad_y,
613                x[2] - 0.001 * grad_z,
614            )
615        }
616    }
617
618    /// Create an attractor-like function
619    pub fn attractor_function(attractor_point: NalgebraVec3, strength: f64) -> impl Fn(NalgebraVec3) -> NalgebraVec3 {
620        move |x: NalgebraVec3| {
621            let direction = attractor_point - x;
622            x + strength * direction
623        }
624    }
625}
626
627#[cfg(test)]
628mod tests {
629    use super::*;
630    use approx::assert_relative_eq;
631
632    #[test]
633    fn test_lipschitz_params_validation() {
634        let params = LipschitzParams::default();
635        assert!(params.validate().is_ok());
636
637        let bad_params = LipschitzParams {
638            lipschitz_constant: 1.5, // > 1, invalid
639            ..params
640        };
641        assert!(bad_params.validate().is_err());
642    }
643
644    #[test]
645    fn test_convergence_rate_estimation() {
646        let params = LipschitzParams::default();
647        let rate = params.convergence_rate();
648        assert!(rate > 0.0);
649
650        let iterations = params.estimated_iterations(10.0);
651        assert!(iterations > 0);
652    }
653
654    #[test]
655    fn test_fixed_point_convergence() {
656        let params = LipschitzParams::fast_convergence();
657        let mut loop_solver = LipschitzLoop::new(params, LoopTopology::FixedPoint).unwrap();
658
659        // Simple contractive function: x' = 0.5 * x
660        let function = |x: NalgebraVec3| 0.5 * x;
661        let initial_state = NalgebraVec3::new(10.0, 10.0, 10.0);
662
663        let result = loop_solver.execute(function, initial_state).unwrap();
664
665        assert!(result.converged);
666        assert!(result.final_residual < 1e-9);
667        assert!(result.iterations > 0);
668    }
669
670    #[test]
671    fn test_newton_convergence() {
672        let params = LipschitzParams::default();
673        let mut loop_solver = LipschitzLoop::new(params, LoopTopology::Newton).unwrap();
674
675        // Quadratic function with known minimum
676        let target = NalgebraVec3::new(1.0, 2.0, 3.0);
677        let function = LoopFunctionFactory::quadratic_bowl(target, 0.1);
678        let initial_state = NalgebraVec3::new(5.0, 5.0, 5.0);
679
680        let result = loop_solver.execute(function, initial_state).unwrap();
681
682        assert!(result.converged);
683        assert!(result.estimated_lipschitz > 0.0);
684    }
685
686    #[test]
687    fn test_accelerated_convergence() {
688        let params = LipschitzParams::default();
689        let mut loop_solver = LipschitzLoop::new(params, LoopTopology::Accelerated).unwrap();
690
691        let function = LoopFunctionFactory::scalar_function(0.0, 0.1);
692        let initial_state = NalgebraVec3::new(5.0, 5.0, 5.0);
693
694        let result = loop_solver.execute(function, initial_state).unwrap();
695
696        assert!(result.converged);
697        assert!(!result.residual_history.is_empty());
698    }
699
700    #[test]
701    fn test_custom_convergence_criteria() {
702        let params = LipschitzParams::default();
703        let mut loop_solver = LipschitzLoop::new(params, LoopTopology::FixedPoint).unwrap();
704
705        let function = |x: NalgebraVec3| 0.9 * x;
706        let initial_state = NalgebraVec3::new(1.0, 1.0, 1.0);
707
708        // Custom criteria: stop when any component is < 0.1
709        let convergence_check = |_current: NalgebraVec3, next: NalgebraVec3, _iter: usize| {
710            next.iter().any(|&x| x.abs() < 0.1)
711        };
712
713        let result = loop_solver.execute_with_criteria(function, initial_state, convergence_check).unwrap();
714
715        assert!(result.converged);
716    }
717
718    #[test]
719    fn test_lipschitz_violation_detection() {
720        let params = LipschitzParams {
721            lipschitz_constant: 0.5, // Very strict
722            adaptive_estimation: true,
723            ..LipschitzParams::default()
724        };
725        let mut loop_solver = LipschitzLoop::new(params, LoopTopology::FixedPoint).unwrap();
726
727        // Expanding function (not contractive)
728        let function = |x: NalgebraVec3| 1.1 * x;
729        let initial_state = NalgebraVec3::new(1.0, 1.0, 1.0);
730
731        let result = loop_solver.execute(function, initial_state);
732
733        // Should detect Lipschitz violation
734        match result {
735            Err(LoopError::LipschitzViolation { .. }) => (),
736            _ => panic!("Expected Lipschitz violation"),
737        }
738    }
739
740    #[test]
741    fn test_stability_analysis() {
742        let params = LipschitzParams::default();
743        let loop_solver = LipschitzLoop::new(params, LoopTopology::FixedPoint).unwrap();
744
745        // Stable function: x' = 0.5 * x
746        let function = |x: NalgebraVec3| 0.5 * x;
747        let fixed_point = NalgebraVec3::zeros();
748
749        let analysis = loop_solver.analyze_stability(function, fixed_point).unwrap();
750
751        assert_eq!(analysis.stability, StabilityType::Stable);
752        assert!(analysis.max_eigenvalue < 1.0);
753    }
754
755    #[test]
756    fn test_loop_function_factory() {
757        let function = LoopFunctionFactory::scalar_function(0.0, 0.1);
758        let input = NalgebraVec3::new(1.0, 1.0, 1.0);
759        let output = function(input);
760
761        // Should move toward zero
762        assert!(output.norm() < input.norm());
763
764        let quadratic = LoopFunctionFactory::quadratic_bowl(NalgebraVec3::zeros(), 0.1);
765        let quad_output = quadratic(input);
766        assert!(quad_output.norm() < input.norm());
767    }
768
769    #[test]
770    fn test_rosenbrock_function() {
771        let function = LoopFunctionFactory::rosenbrock_like(1.0, 100.0);
772        let input = NalgebraVec3::new(0.0, 0.0, 0.0);
773        let output = function(input);
774
775        // Function should be well-defined
776        assert!(output.iter().all(|x| x.is_finite()));
777    }
778
779    #[test]
780    fn test_attractor_function() {
781        let attractor_point = NalgebraVec3::new(1.0, 2.0, 3.0);
782        let function = LoopFunctionFactory::attractor_function(attractor_point, 0.1);
783        let input = NalgebraVec3::zeros();
784        let output = function(input);
785
786        // Should move toward attractor
787        let distance_before = (input - attractor_point).norm();
788        let distance_after = (output - attractor_point).norm();
789        assert!(distance_after < distance_before);
790    }
791
792    #[test]
793    fn test_topology_switching() {
794        let params = LipschitzParams::default();
795        let mut loop_solver = LipschitzLoop::new(params, LoopTopology::FixedPoint).unwrap();
796
797        assert_eq!(*loop_solver.topology(), LoopTopology::FixedPoint);
798
799        loop_solver.set_topology(LoopTopology::Newton);
800        assert_eq!(*loop_solver.topology(), LoopTopology::Newton);
801    }
802
803    #[test]
804    fn test_parameter_updates() {
805        let initial_params = LipschitzParams::default();
806        let mut loop_solver = LipschitzLoop::new(initial_params.clone(), LoopTopology::FixedPoint).unwrap();
807
808        let new_params = LipschitzParams::fast_convergence();
809        loop_solver.update_params(new_params.clone()).unwrap();
810
811        assert_eq!(loop_solver.params().lipschitz_constant, new_params.lipschitz_constant);
812    }
813
814    #[test]
815    fn test_reset_functionality() {
816        let params = LipschitzParams::default();
817        let mut loop_solver = LipschitzLoop::new(params, LoopTopology::FixedPoint).unwrap();
818
819        // Run some iterations
820        let function = |x: NalgebraVec3| 0.9 * x;
821        let initial_state = NalgebraVec3::new(1.0, 1.0, 1.0);
822        let _ = loop_solver.execute(function, initial_state);
823
824        assert!(!loop_solver.state_history().is_empty());
825        assert!(!loop_solver.residual_history().is_empty());
826
827        loop_solver.reset();
828
829        assert!(loop_solver.state_history().is_empty());
830        assert!(loop_solver.residual_history().is_empty());
831        assert_eq!(loop_solver.estimated_lipschitz(), 0.0);
832    }
833}
strange_loop/lipschitz_loop.rs

strange_loop/
lipschitz_loop.rs
