1 files changed, 253 insertions, 0 deletions
diff --git a/src/lib/math/prob_distr.h b/src/lib/math/prob_distr.h
new file mode 100644
index 0000000000..a036073b93
--- /dev/null
+++ b/src/lib/math/prob_distr.h
@@ -0,0 +1,253 @@
+
+/**
+ * \file prob_distr.h
+ *
+ * \brief Header for prob_distr.c
+ **/
+
+#ifndef TOR_PROB_DISTR_H
+#define TOR_PROB_DISTR_H
+
+#include "lib/cc/compat_compiler.h"
+#include "lib/cc/torint.h"
+#include "lib/testsupport/testsupport.h"
+
+/**
+ * Container for distribution parameters for sampling, CDF, &c.
+ */
+struct dist_t {
+  const struct dist_ops_t *ops;
+};
+
+/**
+ * Untyped initializer element for struct dist_t using the specified
+ * struct dist_ops_t pointer.  Don't actually use this directly -- use
+ * the type-specific macro built out of DIST_BASE_TYPED below -- but if
+ * you did use this directly, it would be something like:
+ *
+ *   struct weibull mydist = {
+ *     DIST_BASE(&weibull_ops),
+ *     .lambda = ...,
+ *     .k = ...,
+ *   };
+ *
+ * Note there is NO COMPILER FEEDBACK if you accidentally do something
+ * like
+ *
+ *   struct geometric mydist = {
+ *     DIST_BASE(&weibull_ops),
+ *     ...
+ *   };
+ */
+#define DIST_BASE(OPS)  { .ops = (OPS) }
+
+/** A compile-time type-checking macro for use with DIST_BASE_TYPED.
+ *
+ *  This macro works by checking that &OBJ is a pointer type that is the same
+ *  type (except for qualifiers) as (const TYPE *)&OBJ. It's a C constraint
+ *  violation (which requires a diagnostic) if two pointers are different types
+ *  and are subtracted. The sizeof() forces compile-time evaluation, and the
+ *  multiplication by zero is to discard the result of the sizeof() from the
+ *  expression.
+ *
+ *  We define this conditionally to suppress false positives from
+ *  Coverity, which gets confused by the sizeof business.
+ */
+#ifdef __COVERITY__
+#define TYPE_CHECK_OBJ(OPS, OBJ, TYPE) 0
+#else
+#define TYPE_CHECK_OBJ(OPS, OBJ, TYPE) \
+  (0*sizeof(&(OBJ) - (const TYPE *)&(OBJ)))
+#endif /* defined(__COVERITY__) */
+
+/**
+* Typed initializer element for struct dist_t using the specified struct
+* dist_ops_t pointer.  Don't actually use this directly -- use a
+* type-specific macro built out of it -- but if you did use this
+* directly, it would be something like:
+*
+*     struct weibull mydist = {
+*       DIST_BASE_TYPED(&weibull_ops, mydist, struct weibull_t),
+*       .lambda = ...,
+*       .k = ...,
+*     };
+*
+* If you want to define a distribution type, define a canonical set of
+* operations and define a type-specific initializer element like so:
+*
+*     struct foo_t {
+*       struct dist_t base;
+*       int omega;
+*       double tau;
+*       double phi;
+*     };
+*
+*     struct dist_ops_t foo_ops = ...;
+*
+*     #define FOO(OBJ) DIST_BASE_TYPED(&foo_ops, OBJ, struct foo_t)
+*
+* Then users can do:
+*
+*     struct foo_t mydist = {
+*       FOO(mydist),
+*       .omega = ...,
+*       .tau = ...,
+*       .phi = ...,
+*     };
+*
+* If you accidentally write
+*
+*     struct bar_t mydist = {
+*       FOO(mydist),
+*       ...
+*     };
+*
+* then the compiler will report a type mismatch in the sizeof
+* expression, which otherwise evaporates at runtime.
+*/
+#define DIST_BASE_TYPED(OPS, OBJ, TYPE)                         \
+  DIST_BASE((OPS) + TYPE_CHECK_OBJ(OPS,OBJ,TYPE))
+
+/**
+ * Generic operations on distributions.  These simply defer to the
+ * corresponding dist_ops_t function.  In the parlance of C++, these call
+ * virtual member functions.
+ */
+const char *dist_name(const struct dist_t *);
+double dist_sample(const struct dist_t *);
+double dist_cdf(const struct dist_t *, double x);
+double dist_sf(const struct dist_t *, double x);
+double dist_icdf(const struct dist_t *, double p);
+double dist_isf(const struct dist_t *, double p);
+
+/**
+ * Set of operations on a potentially parametric family of
+ * distributions.  In the parlance of C++, this would be called a
+ * `vtable' and the members are virtual member functions.
+ */
+struct dist_ops_t {
+  const char *name;
+  double (*sample)(const struct dist_t *);
+  double (*cdf)(const struct dist_t *, double x);
+  double (*sf)(const struct dist_t *, double x);
+  double (*icdf)(const struct dist_t *, double p);
+  double (*isf)(const struct dist_t *, double p);
+};
+
+/* Geometric distribution on positive number of trials before first success */
+
+struct geometric_t {
+  struct dist_t base;
+  double p; /* success probability */
+};
+
+extern const struct dist_ops_t geometric_ops;
+
+#define GEOMETRIC(OBJ)                                      \
+  DIST_BASE_TYPED(&geometric_ops, OBJ, struct geometric_t)
+
+/* Pareto distribution */
+
+struct genpareto_t {
+  struct dist_t base;
+  double mu;
+  double sigma;
+  double xi;
+};
+
+extern const struct dist_ops_t genpareto_ops;
+
+#define GENPARETO(OBJ)                                      \
+  DIST_BASE_TYPED(&genpareto_ops, OBJ, struct genpareto_t)
+
+/* Weibull distribution */
+
+struct weibull_t {
+  struct dist_t base;
+  double lambda;
+  double k;
+};
+
+extern const struct dist_ops_t weibull_ops;
+
+#define WEIBULL(OBJ)                                    \
+  DIST_BASE_TYPED(&weibull_ops, OBJ, struct weibull_t)
+
+/* Log-logistic distribution */
+
+struct log_logistic_t {
+  struct dist_t base;
+  double alpha;
+  double beta;
+};
+
+extern const struct dist_ops_t log_logistic_ops;
+
+#define LOG_LOGISTIC(OBJ)                                       \
+  DIST_BASE_TYPED(&log_logistic_ops, OBJ, struct log_logistic_t)
+
+/* Logistic distribution */
+
+struct logistic_t {
+  struct dist_t base;
+  double mu;
+  double sigma;
+};
+
+extern const struct dist_ops_t logistic_ops;
+
+#define LOGISTIC(OBJ)                                   \
+  DIST_BASE_TYPED(&logistic_ops, OBJ, struct logistic_t)
+
+/* Uniform distribution */
+
+struct uniform_t {
+  struct dist_t base;
+  double a;
+  double b;
+};
+
+extern const struct dist_ops_t uniform_ops;
+
+#define UNIFORM(OBJ)                                    \
+  DIST_BASE_TYPED(&uniform_ops, OBJ, struct uniform_t)
+
+/** Only by unittests */
+
+#ifdef PROB_DISTR_PRIVATE
+
+STATIC double logithalf(double p0);
+STATIC double logit(double p);
+
+STATIC double random_uniform_01(void);
+
+STATIC double logistic(double x);
+STATIC double cdf_logistic(double x, double mu, double sigma);
+STATIC double sf_logistic(double x, double mu, double sigma);
+STATIC double icdf_logistic(double p, double mu, double sigma);
+STATIC double isf_logistic(double p, double mu, double sigma);
+STATIC double sample_logistic(uint32_t s, double t, double p0);
+
+STATIC double cdf_log_logistic(double x, double alpha, double beta);
+STATIC double sf_log_logistic(double x, double alpha, double beta);
+STATIC double icdf_log_logistic(double p, double alpha, double beta);
+STATIC double isf_log_logistic(double p, double alpha, double beta);
+STATIC double sample_log_logistic(uint32_t s, double p0);
+
+STATIC double cdf_weibull(double x, double lambda, double k);
+STATIC double sf_weibull(double x, double lambda, double k);
+STATIC double icdf_weibull(double p, double lambda, double k);
+STATIC double isf_weibull(double p, double lambda, double k);
+STATIC double sample_weibull(uint32_t s, double p0, double lambda, double k);
+
+STATIC double sample_uniform_interval(double p0, double a, double b);
+
+STATIC double cdf_genpareto(double x, double mu, double sigma, double xi);
+STATIC double sf_genpareto(double x, double mu, double sigma, double xi);
+STATIC double icdf_genpareto(double p, double mu, double sigma, double xi);
+STATIC double isf_genpareto(double p, double mu, double sigma, double xi);
+STATIC double sample_genpareto(uint32_t s, double p0, double xi);
+
+#endif /* defined(PROB_DISTR_PRIVATE) */
+
+#endif /* !defined(TOR_PROB_DISTR_H) */