use std::ops::Mul;
use rand::Rng;
pub type Point = [f64; 2];
#[derive(Debug, Clone, Copy)]
pub struct TMA {
pub matrix: [[f64; 2]; 2],
pub vector: Point,
pub probability: Option<f64>,
}
impl TMA {
pub fn new(matrix: [[f64; 2]; 2], vector: Point) -> Self {
TMA {
matrix,
vector,
probability: None,
}
}
pub fn identity() -> Self {
TMA {
matrix: [[1.0, 0.0], [0.0, 1.0]],
vector: [0.0, 0.0],
probability: None,
}
}
pub fn from_scale(s: f64) -> Self {
TMA::new([[s, 0.0], [0.0, s]], [0.0, 0.0])
}
pub fn from_translation(tx: f64, ty: f64) -> Self {
TMA::new([[1.0, 0.0], [0.0, 1.0]], [tx, ty])
}
pub fn from_rotation(theta: f64) -> Self {
let (sin_t, cos_t) = theta.sin_cos();
TMA::new([[cos_t, -sin_t], [sin_t, cos_t]], [0.0, 0.0])
}
pub fn with_probability(mut self, p: f64) -> Self {
self.probability = Some(p);
self
}
pub fn apply(&self, p: Point) -> Point {
let x = p[0];
let y = p[1];
let new_x = self.matrix[0][0] * x + self.matrix[0][1] * y + self.vector[0];
let new_y = self.matrix[1][0] * x + self.matrix[1][1] * y + self.vector[1];
[new_x, new_y]
}
pub fn compose(&self, other: &TMA) -> Self {
let m1 = self.matrix;
let m2 = other.matrix;
let new_matrix = [
[
m1[0][0] * m2[0][0] + m1[0][1] * m2[1][0],
m1[0][0] * m2[0][1] + m1[0][1] * m2[1][1],
],
[
m1[1][0] * m2[0][0] + m1[1][1] * m2[1][0],
m1[1][0] * m2[0][1] + m1[1][1] * m2[1][1],
],
];
let c1 = self.vector;
let c2 = other.vector;
let new_vector = [
m1[0][0] * c2[0] + m1[0][1] * c2[1] + c1[0],
m1[1][0] * c2[0] + m1[1][1] * c2[1] + c1[1],
];
TMA::new(new_matrix, new_vector)
}
}
impl Mul<TMA> for TMA {
type Output = TMA;
fn mul(self, rhs: TMA) -> Self::Output {
self.compose(&rhs)
}
}
impl Mul<Point> for TMA {
type Output = Point;
fn mul(self, rhs: Point) -> Self::Output {
self.apply(rhs)
}
}
pub struct IFS {
transformations: Vec<TMA>,
cumulative_probs: Vec<f64>,
}
impl IFS {
pub fn new(transformations: Vec<TMA>) -> Result<Self, &'static str> {
if transformations.is_empty() {
return Err("IFS cannot be empty.");
}
let mut cumulative_probs = Vec::with_capacity(transformations.len());
let mut total_prob = 0.0;
for tma in &transformations {
let prob = tma
.probability
.ok_or("All TMAs in an IFS must have a probability for stochastic generation.")?;
total_prob += prob;
cumulative_probs.push(total_prob);
}
if (total_prob - 1.0).abs() > 1e-6 {
return Err("Probabilities of all TMAs must sum to 1.0.");
}
Ok(IFS {
transformations,
cumulative_probs,
})
}
pub fn choose_transformation(&self) -> &TMA {
let mut rng = rand::thread_rng();
let r: f64 = rng.gen_range(0.0..1.0);
for (i, &cumulative_prob) in self.cumulative_probs.iter().enumerate() {
if r < cumulative_prob {
return &self.transformations[i];
}
}
self.transformations.last().unwrap()
}
pub fn run_chaos_game(
&self,
num_points: usize,
warmup_iterations: usize,
) -> Vec<(Point, usize)> {
let mut points = Vec::with_capacity(num_points);
let mut current_point: Point = [0.0, 0.0];
let mut chosen_index = 0;
let total_iterations = num_points + warmup_iterations;
for i in 0..total_iterations {
let mut rng = rand::thread_rng();
let r: f64 = rng.gen_range(0.0..1.0);
for (j, &cumulative_prob) in self.cumulative_probs.iter().enumerate() {
if r < cumulative_prob {
chosen_index = j;
break;
}
}
let tma = &self.transformations[chosen_index];
current_point = tma.apply(current_point);
if i >= warmup_iterations {
points.push((current_point, chosen_index));
}
}
points
}
}