perturbation.rs
1 #![allow(missing_docs)] 2 3 //! Functionality for constructing perturbations for testing. 4 // TODO: add support for curves 5 // TODO: see how this looks in arbtest 6 7 use proptest::{arbitrary::any, prop_oneof, strategy::Strategy}; 8 9 use crate::geom::Point; 10 11 #[derive(Clone, Copy, Debug)] 12 pub enum F64Perturbation { 13 /// Perturb by between -128 and 127 ulps. 14 Ulp(i8), 15 /// Perturb by a bounded additive amount. 16 Eps(f64), 17 } 18 19 // The Debug bound is because some of the proptest strategy builders want it. Having it here is 20 // easier than copying it everywhere. 21 pub trait FloatPerturbation: std::fmt::Debug { 22 type Float; 23 24 fn apply(&self, f: &Self::Float) -> Self::Float; 25 } 26 27 impl FloatPerturbation for F64Perturbation { 28 type Float = f64; 29 30 fn apply(&self, f: &f64) -> f64 { 31 match self { 32 F64Perturbation::Ulp(n) => { 33 let same_sign = (*n > 0) == (*f > 0.0); 34 let sign_bit = 1 << 63; 35 match f.classify() { 36 std::num::FpCategory::Nan => unreachable!(), 37 std::num::FpCategory::Infinite => *f, 38 std::num::FpCategory::Zero => { 39 let mut bits = n.unsigned_abs().into(); 40 if *n < 0 { 41 bits |= sign_bit; 42 } 43 f64::from_bits(bits) 44 } 45 std::num::FpCategory::Subnormal => { 46 let bits = f.abs().to_bits(); 47 let bits = if same_sign { 48 bits + u64::from(n.unsigned_abs()) 49 } else { 50 bits.abs_diff(u64::from(n.unsigned_abs())) 51 }; 52 f.signum() * f64::from_bits(bits) 53 } 54 std::num::FpCategory::Normal => { 55 let delta = if same_sign { 56 (*n as i64).abs() 57 } else { 58 -(*n as i64).abs() 59 }; 60 // Taking the absolute value sets the sign bit to zero, so the 61 // addition should never overflow. 62 let bits = f.abs().to_bits().checked_add_signed(delta).unwrap(); 63 let x = if bits & sign_bit != 0 { 64 f64::INFINITY 65 } else { 66 f64::from_bits(bits) 67 }; 68 f.signum() * x 69 } 70 } 71 } 72 F64Perturbation::Eps(x) => f + x, 73 } 74 } 75 } 76 77 #[derive(Clone, Copy, Debug)] 78 pub struct PointPerturbation<P: FloatPerturbation> { 79 pub x: P, 80 pub y: P, 81 } 82 83 impl<P: FloatPerturbation<Float = f64>> PointPerturbation<P> { 84 pub fn apply(&self, p: Point) -> Point { 85 Point { 86 x: self.x.apply(&p.x), 87 y: self.y.apply(&p.y), 88 } 89 } 90 } 91 92 #[derive(Clone, Debug)] 93 pub enum Perturbation<P: FloatPerturbation> { 94 Base { 95 idx: usize, 96 }, 97 Point { 98 perturbation: PointPerturbation<P>, 99 idx: usize, 100 next: Box<Perturbation<P>>, 101 }, 102 Subdivision { 103 // Between 0.0 and 1.0 104 t: P::Float, 105 idx: usize, 106 next: Box<Perturbation<P>>, 107 }, 108 Superimposition { 109 left: Box<Perturbation<P>>, 110 right: Box<Perturbation<P>>, 111 }, 112 } 113 114 pub fn f64_perturbation(eps: f64) -> impl Strategy<Value = F64Perturbation> + Clone { 115 prop_oneof![ 116 any::<i8>().prop_map(F64Perturbation::Ulp), 117 (-eps..=eps).prop_map(F64Perturbation::Eps) 118 ] 119 } 120 121 pub fn point_perturbation<P: FloatPerturbation>( 122 fp: impl Strategy<Value = P> + Clone, 123 ) -> impl Strategy<Value = PointPerturbation<P>> { 124 (fp.clone(), fp).prop_map(|(x, y)| PointPerturbation { x, y }) 125 } 126 127 pub fn perturbation( 128 fp: impl Strategy<Value = F64Perturbation> + Clone + 'static, 129 ) -> impl Strategy<Value = Perturbation<F64Perturbation>> { 130 let leaf = any::<usize>().prop_map(|idx| Perturbation::Base { idx }); 131 leaf.prop_recursive(3, 16, 8, move |inner| { 132 prop_oneof![ 133 ( 134 point_perturbation(fp.clone()), 135 any::<usize>(), 136 inner.clone() 137 ) 138 .prop_map(|(perturbation, idx, next)| { 139 Perturbation::Point { 140 perturbation, 141 idx, 142 next: Box::new(next), 143 } 144 }), 145 (0.0f64..1.0, any::<usize>(), inner.clone()).prop_map(|(t, idx, next)| { 146 Perturbation::Subdivision { 147 t, 148 idx, 149 next: Box::new(next), 150 } 151 }), 152 (inner.clone(), inner.clone()).prop_map(|(left, right)| { 153 Perturbation::Superimposition { 154 left: Box::new(left), 155 right: Box::new(right), 156 } 157 }) 158 ] 159 }) 160 } 161 162 fn index<T>(arr: &[T], idx: usize) -> &T { 163 &arr[idx % arr.len()] 164 } 165 166 fn index_mut<T>(arr: &mut [T], idx: usize) -> &mut T { 167 &mut arr[idx % arr.len()] 168 } 169 170 pub fn realize_perturbation( 171 base_cases: &[Vec<Point>], 172 pert: &Perturbation<F64Perturbation>, 173 ) -> Vec<Point> { 174 match pert { 175 Perturbation::Base { idx } => index(base_cases, *idx).to_owned(), 176 Perturbation::Point { 177 perturbation, 178 idx, 179 next, 180 } => { 181 let mut next = realize_perturbation(base_cases, next); 182 let p = index_mut(&mut next, *idx); 183 *p = perturbation.apply(*p); 184 next 185 } 186 Perturbation::Subdivision { t, idx, next } => { 187 let mut next = realize_perturbation(base_cases, next); 188 let idx = *idx % next.len(); 189 let p0 = *index(&next, idx); 190 let p1 = *index(&next, idx + 1); 191 next.insert(idx + 1, p0.affine(&p1, *t)); 192 next 193 } 194 Perturbation::Superimposition { left, right } => { 195 let mut next = realize_perturbation(base_cases, left); 196 next.extend(realize_perturbation(base_cases, right)); 197 next 198 } 199 } 200 }