opendp/transformations/lipschitz_mul/
mod.rs

1use dashu::integer::IBig;
2use opendp_derive::bootstrap;
3
4use crate::{
5    core::{Function, StabilityMap, Transformation},
6    domains::AtomDomain,
7    error::Fallible,
8    metrics::AbsoluteDistance,
9    traits::Float,
10};
11
12#[cfg(feature = "ffi")]
13mod ffi;
14
15#[cfg(test)]
16mod test;
17
18#[bootstrap(
19    features("contrib"),
20    arguments(constant(c_type = "void *")),
21    generics(TA(suppress)),
22    derived_types(TA = "$get_atom(get_carrier_type(input_domain))")
23)]
24/// Make a transformation that multiplies an aggregate by a constant.
25///
26/// The bounds clamp the input, in order to bound the increase in sensitivity from float rounding.
27///
28/// # Arguments
29/// * `input_domain` - The domain of the input.
30/// * `input_metric` - The metric of the input.
31/// * `constant` - The constant to multiply aggregates by.
32/// * `bounds` - Tuple of inclusive lower and upper bounds.
33///
34/// # Generics
35/// * `TA` - Atomic type.
36pub fn make_lipschitz_float_mul<TA: Float>(
37    input_domain: AtomDomain<TA>,
38    input_metric: AbsoluteDistance<TA>,
39    constant: TA,
40    bounds: (TA, TA),
41) -> Fallible<
42    Transformation<AtomDomain<TA>, AbsoluteDistance<TA>, AtomDomain<TA>, AbsoluteDistance<TA>>,
43> {
44    if input_domain.nan() {
45        return fallible!(MakeTransformation, "input_domain may not contain NaN.");
46    }
47
48    let _2 = TA::exact_int_cast(2)?;
49    let (lower, upper) = bounds;
50
51    // Float arithmetic is computed with effectively infinite precision, and rounded to the nearest float
52    // Consider * to be the float multiplication operation
53    // d_out =  max_{v~v'} |f(v) - f(v')|
54    //       <= max_{v~v'} |v * c - v' * c| where * is the float multiplication operation
55    //       =  max_{v~v'} |(vc + e) - (v'c + e')| where e is the float error
56    //       <= max_{v~v'} |(v - v')c + 2e|
57    //       <= max_{v~v'} |v - v'||c| + 2|e| by triangle inequality
58    //       =  d_in|c| + 2 max_{v~v'} |e|
59    //       =  d_in|c| + 2 (ulp(w) / 2) where w = max(|L|, U).inf_mul(|c|), the greatest magnitude output
60    //       =  d_in|c| + ulp(w)
61    //       =  d_in|c| + 2^(exp - k) where exp = exponent(w) and k is the number of bits in the mantissa
62
63    // max(|L|, U)
64    let input_mag = lower.alerting_abs()?.total_max(upper)?;
65
66    // w = max(|L|, U) * |c|
67    let output_mag = input_mag.inf_mul(&constant.alerting_abs()?)?;
68
69    // retrieve the unbiased exponent from the largest output magnitude float w
70    let max_exponent: IBig = output_mag.raw_exponent().into();
71    let max_unbiased_exponent = max_exponent - TA::EXPONENT_BIAS.into();
72
73    // greatest possible error is the ulp of the greatest possible output
74    let output_ulp = _2.inf_powi(max_unbiased_exponent - TA::MANTISSA_BITS.into())?;
75
76    Transformation::new(
77        input_domain,
78        input_metric.clone(),
79        AtomDomain::new_non_nan(),
80        input_metric,
81        Function::new_fallible(move |arg: &TA| {
82            Ok(arg.total_clamp(lower, upper)?.saturating_mul(&constant))
83        }),
84        StabilityMap::new_fallible(move |d_in| {
85            constant.alerting_abs()?.inf_mul(d_in)?.inf_add(&output_ulp)
86        }),
87    )
88}
opendp/transformations/lipschitz_mul/mod.rs

opendp/transformations/lipschitz_mul/
mod.rs
