Commit 492052d7 authored by SirBlueRabbit's avatar SirBlueRabbit
Browse files

refactoring

parent 92212040
......@@ -2,7 +2,7 @@ mod triangulation;
use clap::Parser;
use serde_derive::Deserialize;
use std::{error::Error, fs, io::Write, path::Path, time};
use triangulation::{observable::Observable, MHProbs, Model, StepCount, Triangulation};
use triangulation::{observable::Observable, MHProbs, StepCount, Triangulation};
fn main() -> Result<(), Box<dyn Error>> {
std::env::set_var("RUST_BACKTRACE", "1");
......@@ -118,6 +118,21 @@ struct Config {
measurement: Measurement,
}
#[derive(Clone, Copy, Debug, Deserialize)]
pub struct Model {
pub target_volume: usize,
eps: f32,
weights: Weights,
kappa_0: f32,
kappa_3: f32,
}
#[derive(Clone, Copy, Debug, Deserialize)]
pub struct Weights {
shard: usize,
stars: usize,
}
#[derive(Clone, Debug, Deserialize)]
struct MarkovChain {
thermalisation: usize,
......
use super::MHProbs;
use super::{collections::Label, HalfEdge, Step, Triangulation, Weights};
use super::{collections::Label, HalfEdge, Model, Step, Triangulation, Weights};
impl Triangulation {
pub fn choose_step(&self, weights: Weights, mh_probs: &MHProbs) -> Step {
......@@ -129,6 +128,112 @@ impl Triangulation {
}
}
#[derive(Clone, Debug)]
pub struct MHProbs {
shard_select_grow: Box<[f32]>,
shard_accept_grow: Box<[f32]>,
shard_accept_shrink: Box<[f32]>,
stars_select_grow: Box<[f32]>,
stars_accept_grow: Box<[f32]>,
stars_accept_shrink: Box<[f32]>,
}
impl MHProbs {
pub fn new(model: &Model) -> MHProbs {
let shard = mh_probs(model, 2, 1);
let stars = mh_probs(model, 1, 0);
MHProbs {
shard_select_grow: shard.0,
shard_accept_grow: shard.1,
shard_accept_shrink: shard.2,
stars_select_grow: stars.0,
stars_accept_grow: stars.1,
stars_accept_shrink: stars.2,
}
}
}
fn mh_probs(model: &Model, dn3: usize, dn0: usize) -> (Box<[f32]>, Box<[f32]>, Box<[f32]>) {
let ratio = ratio(model, dn3, dn0);
let norm = norm(&ratio, model.target_volume, dn3);
let select_grow = select_grow(&ratio, &norm, model.target_volume, dn3);
let (accept_grow, accept_shrink) = accept(&norm, model.target_volume, dn3);
(select_grow, accept_grow, accept_shrink)
}
fn ratio(model: &Model, dn3: usize, dn0: usize) -> Box<[f32]> {
let n = model.target_volume;
let eps = model.eps;
let k0 = model.kappa_0;
let k3 = model.kappa_3;
(0..=(2 * n))
.into_iter()
.map(|n3| {
let prefactor = n3 as f32 / ((n3 + dn3) as f32);
let potential = dn3 as isize * (dn3 as isize + 2 * (n3 as isize - n as isize));
let exponent = -k3 * (dn3 as f32) + k0 * (dn0 as f32) - eps * (potential as f32);
prefactor * exponent.exp()
})
.collect()
}
fn norm(ratio: &[f32], n: usize, dn3: usize) -> Box<[f32]> {
(0..=(2 * n))
.into_iter()
.map(|n3| {
if n3 < dn3 {
std::f32::NAN
} else {
1_f32.min(ratio[n3]) + 1_f32.min(1.0 / ratio[n3 - dn3])
}
})
.collect()
}
fn select_grow(ratio: &[f32], norm: &[f32], n: usize, dn3: usize) -> Box<[f32]> {
(0..=(2 * n))
.into_iter()
.map(|n3| {
if n3 < 1 + dn3 {
1_f32
} else if n3 + dn3 > 2 * n {
0_f32
} else {
1_f32.min(ratio[n3]) / norm[n3]
}
})
.collect()
}
fn accept(norm: &[f32], n: usize, dn3: usize) -> (Box<[f32]>, Box<[f32]>) {
let accept_grow = (0..=(2 * n))
.into_iter()
.map(|n3| {
if n3 < 1 + dn3 {
1_f32
} else if n3 + dn3 > 2 * n {
0_f32
} else {
1_f32.min(norm[n3] / norm[n3 + dn3])
}
})
.collect();
let accept_shrink = (0..(norm.len()))
.into_iter()
.map(|n3| {
if n3 < 1 + dn3 {
0_f32
} else if n3 + dn3 > 2 * n {
1_f32
} else {
1_f32.min(norm[n3] / norm[n3 - dn3])
}
})
.collect();
(accept_grow, accept_shrink)
}
impl Weights {
fn cumulative(&self) -> [usize; 2] {
[self.shard, self.shard + self.stars]
......
......@@ -4,6 +4,12 @@ use std::{
ops::{Index, IndexMut},
};
#[derive(Debug, Clone, Copy)]
enum Element<T: Copy> {
Object(T),
Hole(usize),
}
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Label<T> {
pub value: usize,
......@@ -25,12 +31,6 @@ impl<T: Copy> fmt::Display for Label<T> {
}
}
#[derive(Debug, Clone, Copy)]
enum Element<T: Copy> {
Object(T),
Hole(usize),
}
#[derive(Debug, Clone)]
pub struct Pool<T: Copy> {
elements: Box<[Element<T>]>,
......
......@@ -3,10 +3,10 @@ mod collections;
mod do_step;
pub mod observable;
use std::fmt;
use super::{Model, Weights};
pub use choose_step::MHProbs;
use collections::{Bag, Label, Pool};
use serde_derive::Deserialize;
use std::fmt;
const NEXT: [usize; 12] = [1, 2, 0, 4, 5, 3, 7, 8, 6, 10, 11, 9];
const ADJ_INT: [usize; 12] = [3, 7, 11, 0, 10, 8, 9, 1, 5, 6, 4, 2];
......@@ -26,57 +26,6 @@ const VERTICES: [[usize; 2]; 12] = [
[0, 2],
];
#[derive(Clone, Copy, Debug, Deserialize)]
pub struct Model {
pub target_volume: usize,
eps: f32,
weights: Weights,
kappa_0: f32,
kappa_3: f32,
}
#[derive(Clone, Copy, Debug, Deserialize)]
pub struct Weights {
shard: usize,
stars: usize,
}
#[derive(Clone, Debug)]
pub struct MHProbs {
shard_select_grow: Box<[f32]>,
shard_accept_grow: Box<[f32]>,
shard_accept_shrink: Box<[f32]>,
stars_select_grow: Box<[f32]>,
stars_accept_grow: Box<[f32]>,
stars_accept_shrink: Box<[f32]>,
}
#[derive(Debug, Clone)]
pub struct Triangulation {
half_edges: Pool<HalfEdge>,
half_edge_bag: Bag<HalfEdge>,
vertices: Pool<Vertex>,
}
#[derive(Copy, Clone, Debug, PartialEq)]
pub struct HalfEdge {
adj_ext: Label<HalfEdge>,
vertex_tail: Label<Vertex>,
vertex_head: Label<Vertex>,
}
#[derive(Copy, Clone, Debug, PartialEq)]
pub struct Vertex;
#[derive(Clone, Copy, Debug, Default)]
pub struct StepCount {
move02: usize,
move20: usize,
move23: usize,
move32: usize,
invalid: usize,
rejected: usize,
}
#[derive(Clone, Copy, Debug, PartialEq)]
pub enum Step {
Move02(Label<HalfEdge>),
......@@ -87,144 +36,11 @@ pub enum Step {
Rejected,
}
impl fmt::Display for StepCount {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
writeln!(f, "move02: {}", self.move02)?;
writeln!(f, "move20: {}", self.move20)?;
writeln!(f, "move23: {}", self.move23)?;
writeln!(f, "move32: {}", self.move32)?;
writeln!(f, "invalid: {}", self.invalid)?;
writeln!(f, "rejected: {}", self.rejected)?;
writeln!(f, "total: {}", self.total())?;
Ok(())
}
}
impl StepCount {
fn total(&self) -> usize {
self.move02 + self.move20 + self.move23 + self.move32 + self.invalid + self.rejected
}
}
impl MHProbs {
pub fn new(model: &Model) -> MHProbs {
let shard = model.mh_probs(2, 1);
let stars = model.mh_probs(1, 0);
MHProbs {
shard_select_grow: shard.0,
shard_accept_grow: shard.1,
shard_accept_shrink: shard.2,
stars_select_grow: stars.0,
stars_accept_grow: stars.1,
stars_accept_shrink: stars.2,
}
}
}
impl Model {
#[cfg(test)]
pub fn new(
target_volume: usize,
eps: f32,
kappa_0: f32,
kappa_3: f32,
shard: usize,
stars: usize,
) -> Self {
let weights = Weights { shard, stars };
Model {
target_volume,
eps,
kappa_0,
kappa_3,
weights,
}
}
fn mh_probs(&self, dn3: usize, dn0: usize) -> (Box<[f32]>, Box<[f32]>, Box<[f32]>) {
let ratio = self.ratio(dn3, dn0);
let norm = self.norm(&ratio, dn3);
let select_grow = self.select_grow(&ratio, &norm, dn3);
let (accept_grow, accept_shrink) = self.accept(&norm, dn3);
(select_grow, accept_grow, accept_shrink)
}
fn ratio(&self, dn3: usize, dn0: usize) -> Box<[f32]> {
let n = self.target_volume;
let eps = self.eps;
let k0 = self.kappa_0;
let k3 = self.kappa_3;
(0..=(2 * n))
.into_iter()
.map(|n3| {
let prefactor = n3 as f32 / ((n3 + dn3) as f32);
let potential = dn3 as isize * (dn3 as isize + 2 * (n3 as isize - n as isize));
let exponent = -k3 * (dn3 as f32) + k0 * (dn0 as f32) - eps * (potential as f32);
prefactor * exponent.exp()
})
.collect()
}
fn norm(&self, ratio: &[f32], dn3: usize) -> Box<[f32]> {
let n = self.target_volume;
(0..=(2 * n))
.into_iter()
.map(|n3| {
if n3 < dn3 {
std::f32::NAN
} else {
1_f32.min(ratio[n3]) + 1_f32.min(1.0 / ratio[n3 - dn3])
}
})
.collect()
}
fn select_grow(&self, ratio: &[f32], norm: &[f32], dn3: usize) -> Box<[f32]> {
let n = self.target_volume;
(0..=(2 * n))
.into_iter()
.map(|n3| {
if n3 < 1 + dn3 {
1_f32
} else if n3 + dn3 > 2 * n {
0_f32
} else {
1_f32.min(ratio[n3]) / norm[n3]
}
})
.collect()
}
fn accept(&self, norm: &[f32], dn3: usize) -> (Box<[f32]>, Box<[f32]>) {
let n = self.target_volume;
let accept_grow = (0..=(2 * n))
.into_iter()
.map(|n3| {
if n3 < 1 + dn3 {
1_f32
} else if n3 + dn3 > 2 * n {
0_f32
} else {
1_f32.min(norm[n3] / norm[n3 + dn3])
}
})
.collect();
let accept_shrink = (0..(norm.len()))
.into_iter()
.map(|n3| {
if n3 < 1 + dn3 {
0_f32
} else if n3 + dn3 > 2 * n {
1_f32
} else {
1_f32.min(norm[n3] / norm[n3 - dn3])
}
})
.collect();
(accept_grow, accept_shrink)
}
#[derive(Debug, Clone)]
pub struct Triangulation {
half_edges: Pool<HalfEdge>,
half_edge_bag: Bag<HalfEdge>,
vertices: Pool<Vertex>,
}
impl Triangulation {
......@@ -288,6 +104,22 @@ impl Triangulation {
}
}
impl fmt::Display for Triangulation {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
writeln!(f, "half edges: \n{}", self.half_edges)?;
writeln!(f, "half edge bag: {}", self.half_edge_bag)?;
writeln!(f, "vertices: \n{}", self.vertices)?;
Ok(())
}
}
#[derive(Copy, Clone, Debug, PartialEq)]
pub struct HalfEdge {
adj_ext: Label<HalfEdge>,
vertex_tail: Label<Vertex>,
vertex_head: Label<Vertex>,
}
impl Label<HalfEdge> {
pub fn next(&self) -> Label<HalfEdge> {
let val = self.value;
......@@ -352,15 +184,6 @@ impl Pool<HalfEdge> {
}
}
impl fmt::Display for Triangulation {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
writeln!(f, "half edges: \n{}", self.half_edges)?;
writeln!(f, "half edge bag: {}", self.half_edge_bag)?;
writeln!(f, "vertices: \n{}", self.vertices)?;
Ok(())
}
}
impl fmt::Display for HalfEdge {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(
......@@ -371,12 +194,45 @@ impl fmt::Display for HalfEdge {
}
}
#[derive(Copy, Clone, Debug, PartialEq)]
pub struct Vertex;
impl fmt::Display for Vertex {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "")
}
}
#[derive(Clone, Copy, Debug, Default)]
pub struct StepCount {
move02: usize,
move20: usize,
move23: usize,
move32: usize,
invalid: usize,
rejected: usize,
}
impl StepCount {
fn total(&self) -> usize {
self.move02 + self.move20 + self.move23 + self.move32 + self.invalid + self.rejected
}
}
impl fmt::Display for StepCount {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
writeln!(f, "move02: {}", self.move02)?;
writeln!(f, "move20: {}", self.move20)?;
writeln!(f, "move23: {}", self.move23)?;
writeln!(f, "move32: {}", self.move32)?;
writeln!(f, "invalid: {}", self.invalid)?;
writeln!(f, "rejected: {}", self.rejected)?;
writeln!(f, "total: {}", self.total())?;
Ok(())
}
}
#[cfg(test)]
mod tests {
use super::*;
......@@ -510,4 +366,24 @@ mod tests {
assert_eq!(half_edge.vertex_tail, adj_ext.vertex_head);
assert_eq!(half_edge.vertex_head, adj_ext.vertex_tail);
}
impl Model {
pub fn new(
target_volume: usize,
eps: f32,
kappa_0: f32,
kappa_3: f32,
shard: usize,
stars: usize,
) -> Self {
let weights = Weights { shard, stars };
Model {
target_volume,
eps,
kappa_0,
kappa_3,
weights,
}
}
}
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment