Skip to content

Instantly share code, notes, and snippets.

@sebcrozet
Created August 21, 2013 22:04
Show Gist options
  • Save sebcrozet/6300848 to your computer and use it in GitHub Desktop.
Save sebcrozet/6300848 to your computer and use it in GitHub Desktop.
A benchmark showing exhibiting a 80% slowdown after the LLVM update.
extern mod extra;
use std::util;
use std::vec;
#[deriving(Eq, Clone)]
struct Vec3d {
x: f64,
y: f64,
z: f64
}
impl Vec3d {
pub fn new(x: f64, y: f64, z: f64) -> Vec3d {
Vec3d {
x: x,
y: y,
z: z
}
}
pub fn sub_dot(&self, sub: &Vec3d, dot: &Vec3d) -> f64 {
(self.x - sub.x) * dot.x +
(self.y - sub.y) * dot.y +
(self.x - sub.z) * dot.z
}
}
impl Add<Vec3d, Vec3d> for Vec3d {
fn add(&self, other: &Vec3d) -> Vec3d {
Vec3d::new(self.x + other.x, self.y + other.y, self.z + other.z)
}
}
impl Mul<f64, Vec3d> for Vec3d {
fn mul(&self, n: &f64) -> Vec3d {
Vec3d::new(self.x * *n, self.y * *n, self.z * *n)
}
}
impl Div<f64, Vec3d> for Vec3d {
fn div(&self, n: &f64) -> Vec3d {
Vec3d::new(self.x / *n, self.y / *n, self.z / *n)
}
}
/// Simplex using the Johnson subalgorithm to compute the projection of the origin on the simplex.
#[deriving(Eq)]
struct JohnsonSimplex {
priv points: ~[Vec3d],
priv exchange_points: ~[Vec3d],
priv determinants: ~[f64]
}
static permutation_list_0: [uint, ..1] = [0];
static offset_0: [uint, ..2] = [0, 1];
static sub_determinants_0: [uint, ..1] = [0];
static permutation_list_1: [uint, ..6] = [0, 1, 1, 0, 0, 1];
static offset_1: [uint, ..3] = [0, 2, 6];
static sub_determinants_1: [uint, ..6] = [3, 2, 1, 3, 0, 2];
static permutation_list_2: [uint, ..24] = [0, 1, 2, 1, 0, 0, 1, 2, 0, 0, 2, 2, 1, 1, 2, 2, 1, 0, 1,
2, 0, 0, 2, 1];
static offset_2: [uint, ..4] = [0, 3, 15, 24];
static sub_determinants_2: [uint, ..24] = [11, 10, 9, 8, 11, 7, 10, 6, 11, 5, 9, 4, 10, 3, 9, 2, 8,
7, 1, 6, 5, 0, 4, 3];
static permutation_list_3: [uint, ..80] = [0, 1, 2, 3, 1, 0, 0, 1, 2, 0, 0, 2, 3, 0, 0, 3, 2, 1, 1,
2, 3, 1, 1, 3, 3, 2, 2, 3, 2, 1, 0, 1, 2, 0, 0, 2, 1, 3, 1, 0, 1, 3, 0, 0, 3, 1, 3, 2, 0, 2, 3, 0,
0, 3, 2, 3, 2, 1, 2, 3, 1, 1, 3, 2, 3, 2, 1, 0, 2, 3, 1, 0, 1, 3, 2, 0, 0, 3, 2, 1];
static offset_3: [uint, ..5] = [0, 4, 28, 64, 80];
static sub_determinants_3: [uint, ..80] = [31, 30, 29, 28, 27, 31, 26, 30, 25, 31, 24, 29, 23, 31,
22, 28, 21, 30, 20, 29, 19, 30, 18, 28, 17, 29, 16, 28, 15, 27, 26, 14, 25, 24, 13, 21, 20, 12, 27,
26, 11, 23,22, 10, 19, 18, 9, 25, 24, 8, 23, 22, 7, 17, 16, 6, 21, 20, 5, 19, 18, 4, 17, 16, 3, 15,
14, 13, 2, 12, 11, 10, 1, 9, 8, 7, 0, 6, 5, 4];
impl JohnsonSimplex {
/// Creates a new, empty, johnson simplex.
pub fn new() -> JohnsonSimplex {
let _dim = 3;
JohnsonSimplex {
points: vec::with_capacity(_dim + 1),
exchange_points: vec::with_capacity(_dim + 1),
determinants: vec::from_elem(32, 0.0f64),
}
}
}
impl JohnsonSimplex {
fn do_project_origin(&mut self, reduce: bool) -> Vec3d {
if self.points.is_empty() {
fail!("Cannot project the origin on an empty simplex.")
}
if self.points.len() == 1 {
return self.points[0].clone();
}
let max_num_pts = self.points.len();
let permutation_list;
let offset;
let sub_determinants;
if self.points.len() == 1 {
permutation_list = permutation_list_0.slice_to(1);
offset = offset_0.slice_to(2);
sub_determinants = sub_determinants_0.slice_to(1);
}
else if self.points.len() == 2 {
permutation_list = permutation_list_1.slice_to(6);
offset = offset_1.slice_to(3);
sub_determinants = sub_determinants_1.slice_to(6);
}
else if self.points.len() == 3 {
permutation_list = permutation_list_2.slice_to(24);
offset = offset_2.slice_to(4);
sub_determinants = sub_determinants_2.slice_to(24);
}
else {
permutation_list = permutation_list_3.slice_to(80);
offset = offset_3.slice_to(5);
sub_determinants = sub_determinants_3.slice_to(80);
}
let mut curr_num_pts = 1u;
let mut curr = max_num_pts;
for c in self.determinants.mut_slice_from(32 - max_num_pts).mut_iter() {
*c = 1.0;
}
/*
* first loop: compute all the determinants
*/
for &end in offset.slice_from(2).iter() {
// for each sub-simplex ...
while (curr != end) {
unsafe {
let mut determinant = 0.0;
let kpt = self.points.unsafe_get(permutation_list.unsafe_get(curr + 1u)).clone();
let jpt = self.points.unsafe_get(permutation_list.unsafe_get(curr)).clone();
// ... with curr_num_pts points ...
for i in range(curr + 1, curr + 1 + curr_num_pts) {
// ... compute its determinant.
let i_pid = permutation_list.unsafe_get(i);
let sub_determinant = self.determinants.unsafe_get(
sub_determinants.unsafe_get(i)).clone();
let delta = sub_determinant * kpt.sub_dot(&jpt, &self.points.unsafe_get(i_pid));
determinant = determinant + delta;
}
self.determinants.unsafe_set(sub_determinants.unsafe_get(curr), determinant);
curr = curr + curr_num_pts + 1;
}
}
curr_num_pts = curr_num_pts + 1;
}
/*
* second loop: find the subsimplex containing the projection
*/
let mut offsets_iter = offset.rev_iter();
offsets_iter.next(); // skip the first offset
for &end in offsets_iter {
// for each sub-simplex ...
while (curr != end) {
let mut foundit = true;
// ... with curr_num_pts points permutations ...
for i in range(0u, curr_num_pts) {
unsafe {
// ... see if its determinant is positive
let det_id = curr - (i + 1) * curr_num_pts;
let det = self.determinants.unsafe_get(sub_determinants.unsafe_get(det_id));
if det > 0.0 {
// invalidate the children determinant
if curr_num_pts > 1 {
let subdetid = sub_determinants.unsafe_get(det_id + 1);
if self.determinants.unsafe_get(subdetid) > 0.0 {
self.determinants.unsafe_set(subdetid, Bounded::max_value::<f64>())
}
}
// dont concider this sub-simplex if it has been invalidated by its
// parent(s)
if det == Bounded::max_value::<f64>() {
foundit = false
}
}
else {
// we found a negative determinant: no projection possible here
foundit = false
}
}
}
if foundit {
// we found a projection!
// re-run the same iteration but, this time, compute the projection
let mut total_det = 0.0f64;
let mut proj = Vec3d::new(0.0f64, 0.0, 0.0);
unsafe {
for i in range(0u, curr_num_pts) { // FIXME: change this when decreasing loops are implemented
// ... see if its determinant is positive
let id = curr - (i + 1) * curr_num_pts;
let det = self.determinants
.unsafe_get(sub_determinants.unsafe_get(id));
total_det = total_det + det;
proj = proj +
self.points.unsafe_get(permutation_list
.unsafe_get(id)) * det;
}
if reduce {
// we need to reduce the simplex
for i in range(0u, curr_num_pts) {
let id = curr - (i + 1) * curr_num_pts;
self.exchange_points.push(
self.points.unsafe_get(
permutation_list.unsafe_get(id)));
}
util::swap(&mut self.exchange_points, &mut self.points);
self.exchange_points.clear();
}
}
return proj / total_det;
}
curr = curr - curr_num_pts * curr_num_pts;
}
curr_num_pts = curr_num_pts - 1;
}
Vec3d::new(0.0, 0.0, 0.0)
}
#[inline]
fn reset(&mut self, pt: Vec3d) {
self.points.clear();
self.points.push(pt);
}
#[inline]
fn add_point(&mut self, pt: Vec3d) {
self.points.push(pt);
assert!(self.points.len() <= 3 + 1);
}
#[inline]
fn project_origin_and_reduce(&mut self) -> Vec3d {
self.do_project_origin(true)
}
}
#[cfg(test)]
mod test {
use super::*;
use extra::test::BenchHarness;
#[bench]
fn bench_johnson_simplex(bh: &mut BenchHarness) {
let a = Vec3d::new(-0.5, -0.5, -0.5);
let b = Vec3d::new(0.0, 0.5, 0.0);
let c = Vec3d::new(0.5, -0.5, -0.5);
let d = Vec3d::new(0.0, -0.5, -0.5);
do bh.iter {
let mut spl = JohnsonSimplex::new();
do 100.times {
spl.reset(a);
spl.add_point(b);
spl.add_point(c);
spl.add_point(d);
spl.project_origin_and_reduce();
}
}
}
}
fn main() {
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment