organicomplex 0.7.0

/*
 * OrganiComplex
 *
 * Interactive complex-valued cellular automaton on 2D and 3D grids in search
 * of that stuff - emergence, open-endedness, organicity etc.
 * 
 * https://sunkware.org/organicomplex
 * 
 * mediator@sunkware.org
 * 
 * Copyright (c) 2026 Sunkware
 * 
 * OrganiComplex is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * OrganiComplex is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with OrganiComplex. If not, see <https://www.gnu.org/licenses/>.
 */

use rand::prelude::*;
use rand_xoshiro::Xoshiro256PlusPlus;
use rayon::iter::{IndexedParallelIterator, IntoParallelIterator, ParallelIterator};

use std::{
    fs::File,
    io::{
        BufReader,
        BufWriter,
        Read,
        Write
    }
};

use crate::base::{
    Complex,
    COMPLEX_ZERO,
    ReadBin,
    Real,
    REAL_UNIT,
    WriteBin
};

use super::neighbourhood::Neighbourhood;

const DEF_HEIGHT_BINLOG: usize = 0;
const DEF_WIDTH_BINLOG: usize = 9;
const DEF_LENGTH_BINLOG: usize = 9;

#[derive(Clone, Copy, PartialEq, Debug)]
pub struct Cell {
    // pub neighbours: Vec<&'a Cell<'a>>,
    pub amph: Complex,

    pub amphunc: Complex
}

pub struct Field {
    pub cells: Vec<Vec<Vec<Cell>>>,

    pub step: i128,

    pub next_cells: Vec<Vec<Vec<Cell>>>,
    pub energy: f64
}

fn termfunc(a: Complex) -> Complex {
	// a
    a.sqr()
    // let sqr = a.sqr();
    // sqr * (COMPLEX_UNIT + sqr * (COMPLEX_UNIT + sqr))
    // a.exp()
}

fn energy(cells: &Vec<Vec<Vec<Cell>>>) -> f64 {
    let mut energy: f64 = 0.0;
    for layer in cells {
        for row in layer {
            for cell in row {
                energy += cell.amph.absqr().0;
            }
        }
    }
    energy
}

impl<W: Write> WriteBin<&Cell> for W {
    fn write_bin(&mut self, cell: &Cell) -> Result<(), String> {
        self.write_bin(& cell.amph)?;
        Ok(())
    }
}

impl<R: Read> ReadBin<Cell> for R {
    fn read_bin(&mut self) -> Result<Cell, String> {
        let amph: Complex = self.read_bin()?;
        let amphunc = termfunc(amph);
        Ok(Cell{amph, amphunc}) // types are probably derived from here...
    }
}

impl<W: Write> WriteBin<&Field> for W {
    fn write_bin(&mut self, field: &Field) -> Result<(), String> {
        self.write_bin(field.height())?; // z-size
        self.write_bin(field.width())?; // y-size
        self.write_bin(field.length())?; // x-size
        let cells = & field.cells;
        for layer in cells {
            for row in layer {
                for cell in row {
                    self.write_bin(cell)?;
                }            
            }
        }
        self.write_bin(field.step)?;
        Ok(())
    }
}

impl<R: Read> ReadBin<Field> for R {
    fn read_bin(&mut self) -> Result<Field, String> {
        let height: usize = self.read_bin()?;
        let width: usize = self.read_bin()?;
        let length: usize = self.read_bin()?;

        let mut cells = Vec::<Vec<Vec<Cell>>>::with_capacity(height);
        let mut next_cells = Vec::<Vec<Vec<Cell>>>::with_capacity(height);

        for _ in 0..height {
            let mut layer = Vec::<Vec<Cell>>::with_capacity(width);
            let mut next_layer = Vec::<Vec<Cell>>::with_capacity(width);
            for _ in 0..width {
                let mut row = Vec::<Cell>::with_capacity(length);
                let mut next_row = Vec::<Cell>::with_capacity(length);
                for _ in 0..length {
                    let cell = self.read_bin()?;
                    let next_cell = Cell::new_default();
                    row.push(cell);
                    next_row.push(next_cell);
                }
                layer.push(row);
                next_layer.push(next_row);
            }
            cells.push(layer);
            next_cells.push(next_layer);
        }

        let step = self.read_bin()?;

        let energy = energy(&cells);

        Ok(Field{cells, step, next_cells, energy})
    }
}

impl Cell {
    fn new(amph: Complex) -> Self {
        let amphunc = termfunc(amph);
        Self{amph, amphunc}
    }

    fn new_default() -> Self {
        Self::new(COMPLEX_ZERO)
    }

    fn reset_random(&mut self, irng: &mut impl Rng) {
        self.amph = Complex::new_random(1.0, irng);
        self.amphunc = termfunc(self.amph);
    }

    fn deim(&mut self) {
        self.amph = self.amph.deim();
        self.amphunc = termfunc(self.amph);
    }

    fn set(&mut self, v: &Complex) {
        self.amph = *v;
        self.amphunc = termfunc(self.amph);
    }

    #[allow(dead_code)]
    pub fn nullify(&mut self) {
        self.set(&COMPLEX_ZERO);
    }

}

impl Field {
    pub fn height(&self) -> usize {
        self.cells.len()
    }

    pub fn width(&self) -> usize {
        self.cells[0].len()
    }

    pub fn length(&self) -> usize {
        self.cells[0][0].len()
    }

    #[allow(dead_code)]
    pub fn volume(&self) -> usize {
        self.height() * self.width() * self.length()
    }

    pub fn new(height: usize, width: usize, length: usize) -> Self {
        let cells: Vec<Vec<Vec<Cell>>> = (0..height).map(|_| 
            (0..width).map(|_|
                (0..length).map(|_|
                    Cell::new_default()).collect()).collect()).collect();
        let step = 0;

        let next_cells = cells.clone();
        let energy = energy(&cells);

        Self{cells, next_cells, step, energy}
    }

    pub fn new_default() -> Self {
        Self::new(
            1 << DEF_HEIGHT_BINLOG,
            1 << DEF_WIDTH_BINLOG,
            1 << DEF_LENGTH_BINLOG
        )
    }

    pub fn get(&mut self, u: usize, v: usize, w: usize) -> Result<Complex, String> {
        if (u < self.height()) && (v < self.width()) && (w < self.length()) {
            Ok(self.cells[u][v][w].amph)
        } else {
            Err(String::from("cell outside of field"))
        }
    }

    pub fn set(&mut self, u0: usize, v0: usize, w0: usize, vw_radius: usize, s: &Complex) {
        if u0 < self.height() {
            let (v0, w0, vw_radius) = (v0 as isize, w0 as isize, vw_radius as isize);
            let (width_binmask, length_binmask) = ((self.width() - 1) as isize, (self.length() - 1) as isize);
            for dv in -vw_radius..=vw_radius {
                for dw in -vw_radius..=vw_radius {
                    self.cells[u0][((v0 + dv) & width_binmask) as usize][((w0 + dw) & length_binmask) as usize].set(s);
                }
            }
            self.calc_energy();
        }
    }

    pub fn set_outside(&mut self, u0: usize, v0: usize, w0: usize, vw_radius: usize, s: &Complex) {
        if u0 < self.height() {
            let (v0, w0, vw_radius) = (v0 as isize, w0 as isize, vw_radius as isize);
            let (width, length) = (self.width(), self.length());
            for v in 0..width {
                for w in 0..length {
                    if (((v as isize) - v0).abs() > vw_radius) || (((w as isize) - w0).abs() > vw_radius) {
                        self.cells[u0][v][w].set(s);
                    }
                }
            }
            self.calc_energy();
        }
    }

    pub fn fill(&mut self, s: &Complex) {
        let cells = &mut self.cells;
        for layer in cells {
            for row in layer {
                for cell in row {
                    cell.set(s);
                }
            }
        }
        self.step = 0;
        self.calc_energy();
    }

    pub fn nullify(&mut self, u0: usize, v0: usize, w0: usize, vw_radius: usize) {
        self.set(u0, v0, w0, vw_radius, &COMPLEX_ZERO);
    }

    pub fn nullify_outside(&mut self, u0: usize, v0: usize, w0: usize, vw_radius: usize) {
        self.set_outside(u0, v0, w0, vw_radius, &COMPLEX_ZERO);
    }

    pub fn nullify_all(&mut self) {
        self.fill(&COMPLEX_ZERO);
    }

    pub fn reset_random(&mut self, irng: &mut impl Rng) {
        let cells = &mut self.cells;
        for layer in cells {
            for row in layer {
                for cell in row {
                    cell.reset_random(irng);
                }
            }
        }
        self.step = 0;
        self.calc_energy();
    }

    pub fn deim(&mut self) {
        let cells = &mut self.cells;
        for layer in cells {
            for row in layer {
                for cell in row {
                    cell.deim();
                }
            }
        }
        self.calc_energy();
    }

    pub fn add_embryo(&mut self, u0: usize, v0: usize, w0: usize, mult: f64) {
        let v0 = v0 as isize;
        let w0 = w0 as isize;
        let width_binmask = (self.width() - 1) as isize;
        let length_binmask = (self.length() - 1) as isize;

        let half_size: isize = 6;
        let radius = mult / 8.0;

        for dy in -half_size..=half_size {
            let angle = (dy as f64) / ((half_size + 1) as f64) * std::f64::consts::PI;
            let re = radius * (1.0 + angle.cos());
            let im = radius * angle.sin();
            let s = Complex::from((re, im));
            let v = ((v0 + dy) & width_binmask) as usize;
            for dx in 0..(half_size << 1) {
                let w = ((w0 + dx) & length_binmask) as usize;
                self.set(u0, v, w, 0, &s);
            }
            for dx in (-dy - half_size)..0 {
                let v = ((v0 + dy + dx) & width_binmask) as usize;
                let w = ((w0 + dx) & length_binmask) as usize;
                self.set(u0, v, w, 0, &s);
            }
        }
    }

    pub fn calc_energy(&mut self) {
        self.energy = energy(& self.cells);
    }

    pub fn update(&mut self, neighbourhood: &Neighbourhood, amplification: f64, noise: f64, synchronicity: f64, irng: &mut impl Rng) -> Result<(), String> {
        let width = self.width() as usize;

        let (height_binmask, width_binmask, length_binmask) = (
            (self.height() - 1) as isize,
            (width - 1) as isize,
            (self.length() - 1) as isize
        );

        let scalemult = Real::from(amplification) * Real::from(neighbourhood.perpetuator());

        let cells = & self.cells;
        let next_cells = &mut self.next_cells;
        
        let neighbourhood: Vec<(isize, isize, isize)> = neighbourhood.0.iter().map(|&(du, dv, dw)| (du, dv, dw)).collect();

        let mut seeds = vec![0u64; width as usize];

        for (u, next_layer) in next_cells.into_iter().enumerate() {
            let u = u as isize;
            irng.fill(&mut seeds[..]);
            next_layer.into_par_iter().enumerate().for_each(|(v, next_row)| {
                let mut frng = Xoshiro256PlusPlus::seed_from_u64(seeds[v]);
                let v = v as isize;
                for (w, next_cell) in next_row.into_iter().enumerate() {
                    let w = w as isize;

                    let this_cell = & cells[u as usize][v as usize][w as usize];

                    if frng.random_bool(synchronicity) {
                        let this_amph = this_cell.amph;

                        let mut nextamph = COMPLEX_ZERO;
/*
                        // +50% to speed for such Cartesian neighbourhoods
                        let du = 0isize;
                        for dv in -2isize..=2 {
                            for dw in -2isize..=2 {
*/
                        for &(du, dv, dw) in &neighbourhood {
                            let neigh_cell = & cells[((u + du) & height_binmask) as usize][((v + dv) & width_binmask) as usize][((w + dw) & length_binmask) as usize];
                            nextamph += neigh_cell.amphunc;
                            // nextamph += wght * neigh_cell.amph;
                        }
/*
                            }
                        }
*/
                        // nextamph *= nextamph;
                        // nextamph *= this_amph;
                        nextamph = nextamph.conj();

                        // Conway's "Game of Life"
                        /*
                        (1 - c) * s * (s - 1) * (s - 2) * (4 - s) * (5 - s) * (6 - s) * (7 - s) * (8 - s) / 720 +
                        c * (
                            s * (s - 1) * (3 - s) * (4 - s) * (5 - s) * (6 - s) * (7 - s) * (8 - s) / 1440 +
                            s * (s - 1) * (s - 2) * (4 - s) * (5 - s) * (6 - s) * (7 - s) * (8 - s) / 720
                        ) = s * (s - 1) * (4 - s) * (5 - s) * (6 - s) * (7 - s) * (8 - s) / 1440 * (
                            2 * (1 - c) * (s - 2) +
                            c * (
                                3 - s +
                                2 * (s - 2)
                            )
                        ) = s * (s - 1) * (4 - s) * (5 - s) * (6 - s) * (7 - s) * (8 - s) / 1440 * (
                            2 * (1 - c) * (s - 2) + c * (s - 1)
                        ) = s * (s - 1) * (4 - s) * (5 - s) * (6 - s) * (7 - s) * (8 - s) / 1440 * (
                            2s - 2cs - 4 + 4c + cs - c
                        ) = s * (s - 1) * (s - 4) * (s - 5) * (s - 6) * (s - 7) * (s - 8) / 1440 * (cs - 3c - 2s + 4)
                        */
                        
                        nextamph += this_amph;
                        
                        let noise = Complex::new_random(noise, &mut frng);
                        nextamph += noise;

                        let clampabs = nextamph.abs().min(REAL_UNIT);

                        let normult = scalemult * (REAL_UNIT - clampabs);

                        // let abs = nextamph.abs();
                        // let absmod = abs.fract();
                        // let absmod_nl = absmod.exp2() - 1.0;

                        // let normult = absmod * (1.0 - absmod) / abs.max(1e-16);
                        // let normult = (std::f64::consts::PI * absmod).sin() / abs.max(1e-64);
                        // let normult = (1.0 + absmod * (1.0 - absmod)).ln() / abs.max(1e-16);

                        // let normult = (absmod * (1.0 - absmod)).powi(2) / abs.max(1e-16);
                        // let normult = (0.05 * (absmod * 50.0).sin() + absmod * (1.0 - absmod))/ abs.max(1e-16);

                        // let normult = 1.0 / (1.0 + abs);
                        // let normult = 2.0 / std::f64::consts::PI * abs.atan() / abs.max(1e-16);

                        // let normult = scalemult * Real::from(normult);

                        nextamph *= normult;

                        next_cell.amph = nextamph;
                        next_cell.amphunc = termfunc(nextamph);

                    } else {
                        *next_cell = *this_cell;
                    }
                }
            });
        }

        self.cells.swap_with_slice(self.next_cells.as_mut_slice());

        self.step += 1;

        self.calc_energy();

        Ok(())
    }

    pub fn load(filepath: &str) -> Result<Self, String> {
        // "Binary deserialization"... assumes that structures are fixed. Cf. save()
        let mut reader = BufReader::new(File::open(filepath).map_err(|e| e.to_string())?);
        reader.read_bin()
    }

    pub fn save(&self, filepath: &str) -> Result<(), String> {
        // "Binary serialization"... assumes that structures are fixed. Cf. load()
        let mut writer = BufWriter::new(File::create(filepath).map_err(|e| e.to_string())?);
        writer.write_bin(self)?;
        writer.flush().map_err(|e| e.to_string())?;
        Ok(())
    }

}