src/bin/config_parser/src/privacy/poisson.go - cobalt - Git at Google

 // Copyright 2023 The Fuchsia Authors. All rights reserved.
 // Use of this source code is governed by a BSD-style license that can be
 // found in the LICENSE file.

 package privacy

 import (
 	"fmt"
 	"math"

 	"gonum.org/v1/gonum/mathext"
 )

 type poissonDist struct {
 	lambda float64
 }

 func newPoissonDist(lambda float64) (*poissonDist, error) {
 	if lambda < 0 {
 		return nil, fmt.Errorf("Poisson parameter must be greater than 0.")
 	}
 	return &poissonDist{lambda: lambda}, nil
 }

 // Computes the probability of an output less than or equal to |x|.
 func (d poissonDist) cdf(x float64) float64 {
 	if x < 0 {
 		return 0
 	}

 	return mathext.GammaIncRegComp(math.Floor(x+1), d.lambda)
 }

 // Computes the probability of an output greater than or equal to |x|.
 //
 // This is equivalent to 1 - d.cdf(x), but may be more numerically stable.
 func (d poissonDist) sf(x float64) float64 {
 	if x < 0 {
 		return 1
 	}

 	return mathext.GammaIncReg(math.Floor(x+1), d.lambda)
 }

 // Computes the probability of an output of |x| for a poisson distribution with parameter |lambda|.
 func (d poissonDist) pmf(x float64) float64 {
 	if x < 0 || math.Floor(x) != x {
 		return 0
 	}

 	// Since the pmf for the Poisson distribution includes division by a factorial
 	// it is hard to compute directly. This is why we compute the logarithm of the
 	// pmf instead.
 	lg, _ := math.Lgamma(math.Floor(x) + 1.0) // Computes log(x!)
 	lpmf := x*math.Log(d.lambda) - d.lambda - lg
 	return math.Exp(lpmf)
 }

 // Assuming fn is true for integers in the range [0, n] and false for all
 // integers greater than n. This function will return n.
 func findBoundary(fn func(x float64) bool) float64 {
 	c1, c2 := 0.0, 1.0

 	for {
 		if !fn(c2) {
 			break
 		}

 		c1 = c2
 		c2 = c2 * 2
 	}

 	for c1 != c2 {
 		mid := math.Ceil((c1 + c2) / 2)

 		if fn(mid) {
 			c1 = mid
 		} else {
 			c2 = mid - 1
 		}
 	}

 	return c1
 }

 // Assuming fn(x) is true for lowGuess <= x < n (for some n) and false for
 // x >= n, this function will return some value y such that fn(y) is false and
 // n - precision < y < n + precision.
 func findBoundaryFloat(fn func(x float64) bool, lowGuess float64, highGuess float64, precision float64) float64 {
 	l := lowGuess
 	h := highGuess

 	for {
 		if !fn(h) {
 			break
 		}

 		l = h
 		h = h * 2
 	}

 	for (h - l) > precision {
 		mid := l + (h-l)/2
 		if fn(mid) {
 			l = mid
 		} else {
 			h = mid
 		}
 	}

 	return h
 }

 // Creates a privacy loss distribution relative to one observation being added.
 //
 // The Poisson mechanism for computing a scalar-valued function f outputs
 // the sum of the true value of the function and a noise drawn from the Poisson
 // distribution. Recall that the Poisson distribution with parameter p has
 // probability density function exp(-p) * p^x / x! at x for any non-negative
 // integer x. We only consider f with sensitivity 1 for classes defined here.
 //
 // The privacy loss distribution of the Poisson mechanism w.r.t. adding one is
 // equivalent to the privacy loss distribution between the Poisson distribution
 // and the same distribution shifted by +1. Specifically, the privacy loss
 // distribution is generated as follows: first pick x according to the Poisson
 // distribution and, let the privacy loss be ln(pmf(x) / pmf(x - 1)), which
 // is equal to ln(p / x).
 func newPoissonAddPld(lambda, discretization, truncation float64) *PrivacyLossDistribution {
 	dist, _ := newPoissonDist(lambda)
 	lower := findBoundary(func(x float64) bool { return dist.cdf(x-1) <= truncation/2 })
 	if lower < 1.0 {
 		lower = 1.0
 	}

 	upper := 1 + findBoundary(func(x float64) bool { return (dist.sf(x)) > truncation/2 })

 	pmf := make(ProbabilityMassFunction)

 	// Upper truncated tail
 	privacyLoss := math.Log(lambda / upper)
 	roundedPrivacyLoss := int(math.Ceil(privacyLoss / discretization))
 	pmf[roundedPrivacyLoss] = dist.sf(upper)

 	for x := lower; x <= upper; x++ {
 		privacyLoss := math.Log(lambda / x)
 		roundedPrivacyLoss := int(math.Ceil(privacyLoss / discretization))
 		val, ok := pmf[roundedPrivacyLoss]
 		if !ok {
 			val = 0.0
 		}
 		pmf[roundedPrivacyLoss] = val + dist.pmf(x)

 	}

 	// Lower truncated tail
 	infiniteMass := dist.cdf(lower - 1)

 	return &PrivacyLossDistribution{
 		discretization:        discretization,
 		infiniteMass:          infiniteMass,
 		pmf:                   pmf,
 		pessimisticEstimation: true,
 	}
 }

 // Creates a privacy loss distribution relative to one observation being removed.
 //
 // The Poisson mechanism for computing a scalar-valued function f simply outputs
 // the sum of the true value of the function and a noise drawn from the Poisson
 // distribution. Recall that the Poisson distribution with parameter p has
 // probability density function exp(-p) * p^x / x! at x for any non-negative
 // integer x. We only consider f with sensitivity 1 for classes defined here.
 //
 // The privacy loss distribution of the Poisson mechanism w.r.t. removing one is
 // equivalent to the privacy loss distribution between the Poisson distribution
 // and the same distribution shifted by -1. Specifically, the privacy loss
 // distribution is generated as follows: first pick x according to the Poisson
 // distribution. Then, let the privacy loss be ln(pmf(x) / pmf(x + 1)), which
 // is equal to ln((x + 1) / p).
 func newPoissonRemovePld(lambda, discretization, truncation float64) *PrivacyLossDistribution {
 	dist, _ := newPoissonDist(lambda)
 	lower := findBoundary(func(x float64) bool { return dist.cdf(x-2) <= truncation/2 })
 	upper := 1 + findBoundary(func(x float64) bool { return (dist.sf(x - 1)) > truncation/2 })

 	pmf := make(ProbabilityMassFunction)

 	// Lower truncated tail
 	privacyLoss := math.Log(lower / lambda)
 	roundedPrivacyLoss := int(math.Ceil(privacyLoss / discretization))
 	pmf[roundedPrivacyLoss] = dist.cdf(lower - 2)

 	for x := lower; x <= upper; x++ {
 		privacyLoss := math.Log(x / lambda)
 		roundedPrivacyLoss := int(math.Ceil(privacyLoss / discretization))
 		val, ok := pmf[roundedPrivacyLoss]
 		if !ok {
 			val = 0.0
 		}
 		pmf[roundedPrivacyLoss] = val + dist.pmf(x-1)
 	}

 	// Upper truncated tail
 	infiniteMass := dist.sf(upper - 1)

 	return &PrivacyLossDistribution{
 		discretization:        discretization,
 		infiniteMass:          infiniteMass,
 		pmf:                   pmf,
 		pessimisticEstimation: true,
 	}
 }

 // Creates a privacy loss distribution using the specified parameter for the Poisson mechanism.
 func newPoissonPld(lambda float64, sparsity uint64, discretization float64, truncation float64) (*PrivacyLossDistribution, error) {
 	add := newPoissonAddPld(lambda, discretization, truncation)
 	rm := newPoissonRemovePld(lambda, discretization, truncation)
 	pld, err := composeHeterogeneousPLD(add, rm, truncation)

 	if err != nil {
 		return nil, err
 	}

 	return composeHomogeneousPLD(pld, sparsity, truncation)
 }

 type PoissonParameterCalculator struct {
 	discretization float64
 	truncation     float64
 	precision      float64
 }

 // Builds a new PoissonParameterCalculator with default discretization and truncation parameters.
 func NewPoissonParameterCalculator() PoissonParameterCalculator {
 	return PoissonParameterCalculator{
 		discretization: 1e-5,
 		truncation:     1e-15,
 		precision:      1e-5,
 	}
 }

 // Checks if the provided lambda for the specified sparsity provides privacy at the specified epsilon and delta..
 func (c PoissonParameterCalculator) IsPrivate(lambda float64, sparsity uint64, epsilon float64, delta float64) bool {
 	if lambda == 0 {
 		return false
 	}

 	// Delta >= 1 implies no privacy is added.
 	if delta >= 1.0 {
 		return true
 	}

 	// Computes a lower bound on the infinite mass.
 	// This check allows Cobalt to handle very small lambda values.
 	infiniteMassLowerBound := -math.Expm1(float64(sparsity) * math.Log(-math.Expm1(-lambda)))
 	if infiniteMassLowerBound > delta {
 		return false
 	}

 	pld, err := newPoissonPld(lambda, sparsity, c.discretization, c.truncation)
 	if err != nil {
 		return false
 	}

 	return pld.getDeltaForPrivacyLossDistribution(epsilon) < delta
 }

 func (c PoissonParameterCalculator) GetBestPoissonMean(sparsity uint64, epsilon float64, delta float64) float64 {
 	fn := func(x float64) bool {
 		return !c.IsPrivate(x, sparsity, epsilon, delta)
 	}

 	return findBoundaryFloat(fn, 0.0, 1.0, c.precision)
 }
	// Copyright 2023 The Fuchsia Authors. All rights reserved.
	// Use of this source code is governed by a BSD-style license that can be
	// found in the LICENSE file.

	package privacy

	import (
	"fmt"
	"math"

	"gonum.org/v1/gonum/mathext"
	)

	type poissonDist struct {
	lambda float64
	}

	func newPoissonDist(lambda float64) (*poissonDist, error) {
	if lambda < 0 {
	return nil, fmt.Errorf("Poisson parameter must be greater than 0.")
	}
	return &poissonDist{lambda: lambda}, nil
	}

	// Computes the probability of an output less than or equal to \|x\|.
	func (d poissonDist) cdf(x float64) float64 {
	if x < 0 {
	return 0
	}

	return mathext.GammaIncRegComp(math.Floor(x+1), d.lambda)
	}

	// Computes the probability of an output greater than or equal to \|x\|.
	//
	// This is equivalent to 1 - d.cdf(x), but may be more numerically stable.
	func (d poissonDist) sf(x float64) float64 {
	if x < 0 {
	return 1
	}

	return mathext.GammaIncReg(math.Floor(x+1), d.lambda)
	}

	// Computes the probability of an output of \|x\| for a poisson distribution with parameter \|lambda\|.
	func (d poissonDist) pmf(x float64) float64 {
	if x < 0 \|\| math.Floor(x) != x {
	return 0
	}

	// Since the pmf for the Poisson distribution includes division by a factorial
	// it is hard to compute directly. This is why we compute the logarithm of the
	// pmf instead.
	lg, _ := math.Lgamma(math.Floor(x) + 1.0) // Computes log(x!)
	lpmf := x*math.Log(d.lambda) - d.lambda - lg
	return math.Exp(lpmf)
	}

	// Assuming fn is true for integers in the range [0, n] and false for all
	// integers greater than n. This function will return n.
	func findBoundary(fn func(x float64) bool) float64 {
	c1, c2 := 0.0, 1.0

	for {
	if !fn(c2) {
	break
	}

	c1 = c2
	c2 = c2 * 2
	}

	for c1 != c2 {
	mid := math.Ceil((c1 + c2) / 2)

	if fn(mid) {
	c1 = mid
	} else {
	c2 = mid - 1
	}
	}

	return c1
	}

	// Assuming fn(x) is true for lowGuess <= x < n (for some n) and false for
	// x >= n, this function will return some value y such that fn(y) is false and
	// n - precision < y < n + precision.
	func findBoundaryFloat(fn func(x float64) bool, lowGuess float64, highGuess float64, precision float64) float64 {
	l := lowGuess
	h := highGuess

	for {
	if !fn(h) {
	break
	}

	l = h
	h = h * 2
	}

	for (h - l) > precision {
	mid := l + (h-l)/2
	if fn(mid) {
	l = mid
	} else {
	h = mid
	}
	}

	return h
	}

	// Creates a privacy loss distribution relative to one observation being added.
	//
	// The Poisson mechanism for computing a scalar-valued function f outputs
	// the sum of the true value of the function and a noise drawn from the Poisson
	// distribution. Recall that the Poisson distribution with parameter p has
	// probability density function exp(-p) * p^x / x! at x for any non-negative
	// integer x. We only consider f with sensitivity 1 for classes defined here.
	//
	// The privacy loss distribution of the Poisson mechanism w.r.t. adding one is
	// equivalent to the privacy loss distribution between the Poisson distribution
	// and the same distribution shifted by +1. Specifically, the privacy loss
	// distribution is generated as follows: first pick x according to the Poisson
	// distribution and, let the privacy loss be ln(pmf(x) / pmf(x - 1)), which
	// is equal to ln(p / x).
	func newPoissonAddPld(lambda, discretization, truncation float64) *PrivacyLossDistribution {
	dist, _ := newPoissonDist(lambda)
	lower := findBoundary(func(x float64) bool { return dist.cdf(x-1) <= truncation/2 })
	if lower < 1.0 {
	lower = 1.0
	}

	upper := 1 + findBoundary(func(x float64) bool { return (dist.sf(x)) > truncation/2 })

	pmf := make(ProbabilityMassFunction)

	// Upper truncated tail
	privacyLoss := math.Log(lambda / upper)
	roundedPrivacyLoss := int(math.Ceil(privacyLoss / discretization))
	pmf[roundedPrivacyLoss] = dist.sf(upper)

	for x := lower; x <= upper; x++ {
	privacyLoss := math.Log(lambda / x)
	roundedPrivacyLoss := int(math.Ceil(privacyLoss / discretization))
	val, ok := pmf[roundedPrivacyLoss]
	if !ok {
	val = 0.0
	}
	pmf[roundedPrivacyLoss] = val + dist.pmf(x)

	}

	// Lower truncated tail
	infiniteMass := dist.cdf(lower - 1)

	return &PrivacyLossDistribution{
	discretization: discretization,
	infiniteMass: infiniteMass,
	pmf: pmf,
	pessimisticEstimation: true,
	}
	}

	// Creates a privacy loss distribution relative to one observation being removed.
	//
	// The Poisson mechanism for computing a scalar-valued function f simply outputs
	// the sum of the true value of the function and a noise drawn from the Poisson
	// distribution. Recall that the Poisson distribution with parameter p has
	// probability density function exp(-p) * p^x / x! at x for any non-negative
	// integer x. We only consider f with sensitivity 1 for classes defined here.
	//
	// The privacy loss distribution of the Poisson mechanism w.r.t. removing one is
	// equivalent to the privacy loss distribution between the Poisson distribution
	// and the same distribution shifted by -1. Specifically, the privacy loss
	// distribution is generated as follows: first pick x according to the Poisson
	// distribution. Then, let the privacy loss be ln(pmf(x) / pmf(x + 1)), which
	// is equal to ln((x + 1) / p).
	func newPoissonRemovePld(lambda, discretization, truncation float64) *PrivacyLossDistribution {
	dist, _ := newPoissonDist(lambda)
	lower := findBoundary(func(x float64) bool { return dist.cdf(x-2) <= truncation/2 })
	upper := 1 + findBoundary(func(x float64) bool { return (dist.sf(x - 1)) > truncation/2 })

	pmf := make(ProbabilityMassFunction)

	// Lower truncated tail
	privacyLoss := math.Log(lower / lambda)
	roundedPrivacyLoss := int(math.Ceil(privacyLoss / discretization))
	pmf[roundedPrivacyLoss] = dist.cdf(lower - 2)

	for x := lower; x <= upper; x++ {
	privacyLoss := math.Log(x / lambda)
	roundedPrivacyLoss := int(math.Ceil(privacyLoss / discretization))
	val, ok := pmf[roundedPrivacyLoss]
	if !ok {
	val = 0.0
	}
	pmf[roundedPrivacyLoss] = val + dist.pmf(x-1)
	}

	// Upper truncated tail
	infiniteMass := dist.sf(upper - 1)

	return &PrivacyLossDistribution{
	discretization: discretization,
	infiniteMass: infiniteMass,
	pmf: pmf,
	pessimisticEstimation: true,
	}
	}

	// Creates a privacy loss distribution using the specified parameter for the Poisson mechanism.
	func newPoissonPld(lambda float64, sparsity uint64, discretization float64, truncation float64) (*PrivacyLossDistribution, error) {
	add := newPoissonAddPld(lambda, discretization, truncation)
	rm := newPoissonRemovePld(lambda, discretization, truncation)
	pld, err := composeHeterogeneousPLD(add, rm, truncation)

	if err != nil {
	return nil, err
	}

	return composeHomogeneousPLD(pld, sparsity, truncation)
	}

	type PoissonParameterCalculator struct {
	discretization float64
	truncation float64
	precision float64
	}

	// Builds a new PoissonParameterCalculator with default discretization and truncation parameters.
	func NewPoissonParameterCalculator() PoissonParameterCalculator {
	return PoissonParameterCalculator{
	discretization: 1e-5,
	truncation: 1e-15,
	precision: 1e-5,
	}
	}

	// Checks if the provided lambda for the specified sparsity provides privacy at the specified epsilon and delta..
	func (c PoissonParameterCalculator) IsPrivate(lambda float64, sparsity uint64, epsilon float64, delta float64) bool {
	if lambda == 0 {
	return false
	}

	// Delta >= 1 implies no privacy is added.
	if delta >= 1.0 {
	return true
	}

	// Computes a lower bound on the infinite mass.
	// This check allows Cobalt to handle very small lambda values.
	infiniteMassLowerBound := -math.Expm1(float64(sparsity) * math.Log(-math.Expm1(-lambda)))
	if infiniteMassLowerBound > delta {
	return false
	}

	pld, err := newPoissonPld(lambda, sparsity, c.discretization, c.truncation)
	if err != nil {
	return false
	}

	return pld.getDeltaForPrivacyLossDistribution(epsilon) < delta
	}

	func (c PoissonParameterCalculator) GetBestPoissonMean(sparsity uint64, epsilon float64, delta float64) float64 {
	fn := func(x float64) bool {
	return !c.IsPrivate(x, sparsity, epsilon, delta)
	}

	return findBoundaryFloat(fn, 0.0, 1.0, c.precision)
	}