src/algorithms/privacy/poisson.cc

#include "src/algorithms/privacy/poisson.h"

#include <fenv.h>

#include <cmath>
#include <limits>

#include "src/algorithms/random/distributions.h"

namespace cobalt {

std::vector<uint64_t> ApplyPoissonNoise(const std::vector<uint64_t>& indices, uint64_t max_index,
                                        PoissonParameter lambda,
                                        SecureBitGeneratorInterface<uint32_t>* gen) {
  // Set rounding mode to rounding upwards in order to prevent accidentally using a
  // lower-than-expected noise parameter.
  int rounding_mode = fegetround();
  fesetround(FE_UPWARD);
  double max_index_plus_1 = static_cast<double>(max_index) + 1;
  double lambda_times_num_index = static_cast<double>(max_index_plus_1 * lambda.get());
  // Restore rounding mode.
  fesetround(rounding_mode);

  // To minimize the number of calls to |gen|, we draw the total number of observations
  // then, we distribute those additional observations over the indices. This is more
  // efficient than drawing from a Poisson distribution once per index when lambda < 1.
  // We expect lambda << 1 in practice.
  PoissonDistribution poisson(gen, PoissonParameter(lambda_times_num_index));
  DiscreteUniformDistribution uniform(gen, 0, max_index);

  std::vector<uint64_t> with_noise(indices.begin(), indices.end());

  uint32_t added_ones = poisson.Sample();
  for (uint32_t i = 0; i < added_ones; i++) {
    with_noise.push_back(uniform.Sample());
  }

  return with_noise;
}

double DebiasPoissonCount(double raw_count, uint64_t num_contributions, PoissonParameter lambda) {
  return (raw_count - lambda.get() * static_cast<double>(num_contributions));
}

}  // namespace cobalt