Created
August 21, 2013 22:04
-
-
Save sebcrozet/6300848 to your computer and use it in GitHub Desktop.
A benchmark showing exhibiting a 80% slowdown after the LLVM update.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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