compb-dla-model/src/tools/boxcount.rs
2023-03-17 18:58:03 +00:00

136 lines
3.9 KiB
Rust

use std::fs::File;
use std::os::unix::fs::symlink;
use bevy::tasks::ParallelSlice;
use crate::system::spaces::square_grid::{Grid2D, Grid3D};
use itertools::{Itertools, MinMaxResult};
use clap::Parser;
use serde::Serialize;
use crate::BoxCountCli;
use crate::cli::cli::OutputFormat;
use crate::system::{GriddedPosition, Position};
use crate::system::model::HistoryLine;
use crate::tools::read;
fn bb(data: &Vec<Grid2D>) -> ((i32, i32), (i32, i32)) {
let x = data
.iter().minmax_by(|a, b| a.x.cmp(&b.x));
let y = data
.iter().minmax_by(|a, b| a.y.cmp(&b.y));
match (x, y) {
(MinMaxResult::MinMax(min_x, max_x), MinMaxResult::MinMax(min_y, max_y)) => {
((min_x.x, min_y.y), (max_x.x, max_y.y))
}
_ => panic!("Cannot determine bounding box")
}
}
fn box_count_2d(data: &Vec<Grid2D>, box_number: u32) -> (f64, usize) {
let ((x_min, y_min), (x_max, y_max)) = bb(data);
let x_range = (x_max - x_min) as f64;
let y_range = (y_max - y_min) as f64;
let w: f64 = x_range / (box_number as f64);
let boxes_occupied = data.iter()
.map(|Grid2D { x, y }| [((x - x_min) as f64 / w) as i32, ((y - y_min) as f64 / w) as i32])
.unique()
.count();
(w, boxes_occupied)
}
fn box_count_3d(data: &Vec<Grid3D>, size: u32) -> usize {
let n = data.len();
let x_min = data
.iter()
.min_by(|Grid3D { x: x1, .. }, Grid3D { x: x2, .. }| x1.cmp(x2))
.unwrap().x;
let x_max = data
.iter()
.max_by(|Grid3D { x: x1, .. }, Grid3D { x: x2, .. }| x1.cmp(x2))
.unwrap().x;
let y_min = data
.iter()
.min_by(|Grid3D { y: v1, .. }, Grid3D { y: v2, .. }| v1.cmp(v2))
.unwrap().y;
let y_max = data
.iter()
.max_by(|Grid3D { y: v1, .. }, Grid3D { y: v2, .. }| v1.cmp(v2))
.unwrap().y;
let z_min = data
.iter()
.min_by(|Grid3D { z: v1, .. }, Grid3D { z: v2, .. }| v1.cmp(v2))
.unwrap().y;
let z_max = data
.iter()
.max_by(|Grid3D { z: v1, .. }, Grid3D { z: v2, .. }| v1.cmp(v2))
.unwrap().y;
let x_range = (x_max - x_min) as f64;
let y_range = (y_max - y_min) as f64;
let z_range = (z_max - z_min) as f64;
let w: f64 = x_range / (size as f64);
let grid_points = data.iter()
.map(|Grid3D { x, y, z }| [
((x - x_min) as f64 / w) as i32,
((y - y_min) as f64 / w) as i32,
((z - z_min) as f64 / w) as i32,
])
.collect::<Vec<_>>();
return grid_points.iter()
.unique()
.count();
}
fn box_count_nd<const N: usize>(data: &Vec<[f32; N]>, size: u32) -> usize {
let ranges = (0..N).map(|n|
match data.iter()
.minmax_by(|a, b| a[n].total_cmp(&b[n])) {
MinMaxResult::NoElements => panic!("No data"),
MinMaxResult::OneElement(_) => panic!("Needs more than one point to compute boxcount"),
MinMaxResult::MinMax(min, max) => [min[n], min[n]],
}).collect::<Vec<_>>();
let w: f32 = (ranges[0][1] - ranges[0][0]) / (size as f32);
return data.iter()
.map(|point| -> [i32; N] { std::array::from_fn(|n| ((point[n] - ranges[n][0]) / w) as i32) })
.unique()
.count();
}
#[derive(Serialize)]
struct FDRow {
w: f64,
n_occupied: usize,
}
pub(crate) fn main(cli: &BoxCountCli) {
let particles: Vec<Grid2D> = read(&cli.path, cli.format);
let n_particles = dbg!(particles.len());
let box_side_counts = 1..500;
let mut writer = csv::Writer::from_path(&cli.output)
.expect("Unable to create csv");
box_side_counts
.map(|box_side_count| box_count_2d(&particles, box_side_count))
// .filter(|(w, n_occupied)| *n_occupied < n_particles) // Remove saturated values
.map(|(w, n_occupied)| FDRow { w, n_occupied })
.for_each(|row| writer.serialize(row).expect("Failed to write row"));
writer.flush().unwrap();
}