From 993cbd0c4e2becacf29188357c8b53d6605e584e Mon Sep 17 00:00:00 2001 From: leventov Date: Sun, 20 Sep 2015 01:39:18 +0300 Subject: [PATCH] Add ChronicleHashBuilder.nonTieredSegmentsPercentile() configuration and use proper math for segment sizing --- .../chronicle/hash/ChronicleHashBuilder.java | 37 +- .../impl/util/math/ContinuedFraction.java | 136 +++++ .../chronicle/hash/impl/util/math/Gamma.java | 509 ++++++++++++++++++ .../impl/util/math/PoissonDistribution.java | 162 ++++++ .../hash/impl/util/math/Precision.java | 98 ++++ .../hash/impl/util/math/package-info.java | 22 + .../chronicle/map/ChronicleMapBuilder.java | 76 ++- .../chronicle/set/ChronicleSetBuilder.java | 6 + 8 files changed, 1002 insertions(+), 44 deletions(-) create mode 100644 src/main/java/net/openhft/chronicle/hash/impl/util/math/ContinuedFraction.java create mode 100644 src/main/java/net/openhft/chronicle/hash/impl/util/math/Gamma.java create mode 100644 src/main/java/net/openhft/chronicle/hash/impl/util/math/PoissonDistribution.java create mode 100644 src/main/java/net/openhft/chronicle/hash/impl/util/math/Precision.java create mode 100644 src/main/java/net/openhft/chronicle/hash/impl/util/math/package-info.java diff --git a/src/main/java/net/openhft/chronicle/hash/ChronicleHashBuilder.java b/src/main/java/net/openhft/chronicle/hash/ChronicleHashBuilder.java index 4b8b50b64..d89d9350f 100755 --- a/src/main/java/net/openhft/chronicle/hash/ChronicleHashBuilder.java +++ b/src/main/java/net/openhft/chronicle/hash/ChronicleHashBuilder.java @@ -54,8 +54,8 @@ import java.util.concurrent.TimeUnit; * *

There are some "low-level" configurations in this builder, * that require deep understanding of the Chronicle implementation design to be properly used. - * Know what you do. These configurations are picked up strictly as-is, without extra round-ups, - * adjustments, etc. + * Know what you do. These configurations are applied as-is, without extra round-ups, adjustments, + * etc. * * @param the type of keys in hash containers, created by this builder * @param the container type, created by this builder, i. e. {@link ChronicleMap} or {@link @@ -201,8 +201,8 @@ public interface ChronicleHashBuilder, * .create();} * *

This is a low-level configuration. The configured number - * of bytes is strictly used as-is, without anything like round-up to the multiple of 8 or - * 16, or any other adjustment. + * of bytes is used as-is, without anything like round-up to the multiple of 8 or 16, or any + * other adjustment. * * @param actualChunkSize the "chunk size" in bytes * @return this builder back @@ -312,6 +312,35 @@ public interface ChronicleHashBuilder, */ B allowSegmentTiering(boolean allowSegmentTiering); + /** + * Configures probabilistic fraction of segments, which shouldn't become tiered, if Chronicle + * Hash size is {@link #entries(long)}, assuming hash code distribution of the keys, inserted + * into configured Chronicle Hash, is good. + * + *

The last caveat means that the configured percentile and affects segment size relying on + * Poisson distribution law, if inserted entries (keys) fall into all segments randomly. If + * e. g. the keys, inserted into the Chronicle Hash, are purposely selected to collide by + * a certain range of hash code bits, so that they all fall into the same segment (a DOS + * attacker might do this), this segment is obviously going to be tiered. + * + *

This configuration affects the actual number of segments, if {@link #entries(long)} and + * {@link #entriesPerSegment(long)} or {@link #actualChunksPerSegment(long)} are configured. It + * affects the actual number of entries/chunks per segment, if {@link #entries(long)} and + * {@link #actualSegments(int)} are configured. If all 4 configurations, mentioned in this + * paragraph, are specified, {@code nonTieredSegmentsPercentile} is irrelevant. + * + *

Default value is 0.99999, i. e. if hash code distribution of the keys is good, only one + * segment of 100K is tiered on average. If your segment size is small and you want to improve + * memory footprint of Chronicle Hash (probably compromising latency percentiles), you might + * want to configure more "relaxed" value, e. g. 0.99. + * + * @param nonTieredSegmentsPercentile Fraction of segments which shouldn't be tiered + * @return this builder back + * @throws IllegalArgumentException if {@code nonTieredSegmentsPercentile} is out of (0.5, 1.0) + * range (both bounds are excluded) + */ + B nonTieredSegmentsPercentile(double nonTieredSegmentsPercentile); + /** * Configures the actual maximum number entries, that could be inserted into any single segment * of the hash containers, created by this builder. Configuring both the actual number of diff --git a/src/main/java/net/openhft/chronicle/hash/impl/util/math/ContinuedFraction.java b/src/main/java/net/openhft/chronicle/hash/impl/util/math/ContinuedFraction.java new file mode 100644 index 000000000..ba94ba4eb --- /dev/null +++ b/src/main/java/net/openhft/chronicle/hash/impl/util/math/ContinuedFraction.java @@ -0,0 +1,136 @@ +/* + * Copyright (C) 2015 higherfrequencytrading.com + * Copyright 2001-2015 The Apache Software Foundation + * Copyright 2010-2012 CS Systèmes d'Information + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +package net.openhft.chronicle.hash.impl.util.math; + +/** + * Provides a generic means to evaluate continued fractions. Subclasses simply + * provided the a and b coefficients to evaluate the continued fraction. + * + *

+ * References: + *

+ *

+ * + */ +abstract class ContinuedFraction { + + /** + * Access the n-th a coefficient of the continued fraction. Since a can be + * a function of the evaluation point, x, that is passed in as well. + * @param n the coefficient index to retrieve. + * @param x the evaluation point. + * @return the n-th a coefficient. + */ + protected abstract double getA(int n, double x); + + /** + * Access the n-th b coefficient of the continued fraction. Since b can be + * a function of the evaluation point, x, that is passed in as well. + * @param n the coefficient index to retrieve. + * @param x the evaluation point. + * @return the n-th b coefficient. + */ + protected abstract double getB(int n, double x); + + /** + * Evaluates the continued fraction at the value x. + *

+ * The implementation of this method is based on the modified Lentz algorithm as described + * on page 18 ff. in: + *

+ * Note: the implementation uses the terms ai and bi as defined in + * Continued Fraction @ MathWorld. + *

+ * + * @param x the evaluation point. + * @param epsilon maximum error allowed. + * @param maxIterations maximum number of convergents + * @return the value of the continued fraction evaluated at x. + * @throws IllegalStateException if the algorithm fails to converge. + * @throws IllegalStateException if maximal number of iterations is reached + */ + public double evaluate(double x, double epsilon, int maxIterations) { + final double small = 1e-50; + double hPrev = getA(0, x); + + // use the value of small as epsilon criteria for zero checks + if (Precision.equals(hPrev, 0.0, small)) { + hPrev = small; + } + + int n = 1; + double dPrev = 0.0; + double cPrev = hPrev; + double hN = hPrev; + + while (n < maxIterations) { + final double a = getA(n, x); + final double b = getB(n, x); + + double dN = a + b * dPrev; + if (Precision.equals(dN, 0.0, small)) { + dN = small; + } + double cN = a + b / cPrev; + if (Precision.equals(cN, 0.0, small)) { + cN = small; + } + + dN = 1 / dN; + final double deltaN = cN * dN; + hN = hPrev * deltaN; + + if (Double.isInfinite(hN)) { + throw new IllegalStateException( + "Continued fraction convergents diverged to +/- infinity for value " + x); + } + if (Double.isNaN(hN)) { + throw new IllegalStateException( + "Continued fraction diverged to NaN for value " + x); + } + + if (Math.abs(deltaN - 1.0) < epsilon) { + break; + } + + dPrev = dN; + cPrev = cN; + hPrev = hN; + n++; + } + + if (n >= maxIterations) { + throw new IllegalStateException( + "Continued fraction convergents failed to converge (in less than " + + maxIterations + " iterations) for value " + x); + } + + return hN; + } + +} diff --git a/src/main/java/net/openhft/chronicle/hash/impl/util/math/Gamma.java b/src/main/java/net/openhft/chronicle/hash/impl/util/math/Gamma.java new file mode 100644 index 000000000..a9ca5fad4 --- /dev/null +++ b/src/main/java/net/openhft/chronicle/hash/impl/util/math/Gamma.java @@ -0,0 +1,509 @@ +/* + * Copyright (C) 2015 higherfrequencytrading.com + * Copyright 2001-2015 The Apache Software Foundation + * Copyright 2010-2012 CS Systèmes d'Information + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +package net.openhft.chronicle.hash.impl.util.math; + + +/** + *

+ * This is a utility class that provides computation methods related to the + * Γ (Gamma) family of functions. + *

+ *

+ * Implementation of {@link #invGamma1pm1(double)} and + * {@link #logGamma1p(double)} is based on the algorithms described in + *

+ * and implemented in the + * NSWC Library of Mathematical Functions, + * available + * here. + * This library is "approved for public release", and the + * Copyright guidance + * indicates that unless otherwise stated in the code, all FORTRAN functions in + * this library are license free. Since no such notice appears in the code these + * functions can safely be ported to Commons-Math. + *

+ * + */ +class Gamma { + + /** + * The value of the {@code g} constant in the Lanczos approximation, see + * {@link #lanczos(double)}. + * @since 3.1 + */ + public static final double LANCZOS_G = 607.0 / 128.0; + + /** Lanczos coefficients */ + private static final double[] LANCZOS = { + 0.99999999999999709182, + 57.156235665862923517, + -59.597960355475491248, + 14.136097974741747174, + -0.49191381609762019978, + .33994649984811888699e-4, + .46523628927048575665e-4, + -.98374475304879564677e-4, + .15808870322491248884e-3, + -.21026444172410488319e-3, + .21743961811521264320e-3, + -.16431810653676389022e-3, + .84418223983852743293e-4, + -.26190838401581408670e-4, + .36899182659531622704e-5, + }; + + /** Avoid repeated computation of log of 2 PI in logGamma */ + private static final double HALF_LOG_2_PI = 0.5 * Math.log(2.0 * Math.PI); + + /* + * Constants for the computation of double invGamma1pm1(double). + * Copied from DGAM1 in the NSWC library. + */ + + /** The constant {@code A0} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_A0 = .611609510448141581788E-08; + + /** The constant {@code A1} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_A1 = .624730830116465516210E-08; + + /** The constant {@code B1} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B1 = .203610414066806987300E+00; + + /** The constant {@code B2} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B2 = .266205348428949217746E-01; + + /** The constant {@code B3} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B3 = .493944979382446875238E-03; + + /** The constant {@code B4} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B4 = -.851419432440314906588E-05; + + /** The constant {@code B5} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B5 = -.643045481779353022248E-05; + + /** The constant {@code B6} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B6 = .992641840672773722196E-06; + + /** The constant {@code B7} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B7 = -.607761895722825260739E-07; + + /** The constant {@code B8} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_B8 = .195755836614639731882E-09; + + /** The constant {@code P0} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P0 = .6116095104481415817861E-08; + + /** The constant {@code P1} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P1 = .6871674113067198736152E-08; + + /** The constant {@code P2} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P2 = .6820161668496170657918E-09; + + /** The constant {@code P3} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P3 = .4686843322948848031080E-10; + + /** The constant {@code P4} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P4 = .1572833027710446286995E-11; + + /** The constant {@code P5} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P5 = -.1249441572276366213222E-12; + + /** The constant {@code P6} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_P6 = .4343529937408594255178E-14; + + /** The constant {@code Q1} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_Q1 = .3056961078365221025009E+00; + + /** The constant {@code Q2} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_Q2 = .5464213086042296536016E-01; + + /** The constant {@code Q3} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_Q3 = .4956830093825887312020E-02; + + /** The constant {@code Q4} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_Q4 = .2692369466186361192876E-03; + + /** The constant {@code C} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C = -.422784335098467139393487909917598E+00; + + /** The constant {@code C0} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C0 = .577215664901532860606512090082402E+00; + + /** The constant {@code C1} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C1 = -.655878071520253881077019515145390E+00; + + /** The constant {@code C2} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C2 = -.420026350340952355290039348754298E-01; + + /** The constant {@code C3} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C3 = .166538611382291489501700795102105E+00; + + /** The constant {@code C4} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C4 = -.421977345555443367482083012891874E-01; + + /** The constant {@code C5} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C5 = -.962197152787697356211492167234820E-02; + + /** The constant {@code C6} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C6 = .721894324666309954239501034044657E-02; + + /** The constant {@code C7} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C7 = -.116516759185906511211397108401839E-02; + + /** The constant {@code C8} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C8 = -.215241674114950972815729963053648E-03; + + /** The constant {@code C9} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C9 = .128050282388116186153198626328164E-03; + + /** The constant {@code C10} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C10 = -.201348547807882386556893914210218E-04; + + /** The constant {@code C11} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C11 = -.125049348214267065734535947383309E-05; + + /** The constant {@code C12} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C12 = .113302723198169588237412962033074E-05; + + /** The constant {@code C13} defined in {@code DGAM1}. */ + private static final double INV_GAMMA1P_M1_C13 = -.205633841697760710345015413002057E-06; + + /** + * Default constructor. Prohibit instantiation. + */ + private Gamma() {} + + /** + *

+ * Returns the value of log Γ(x) for x > 0. + *

+ *

+ * For x ≤ 8, the implementation is based on the double precision + * implementation in the NSWC Library of Mathematics Subroutines, + * {@code DGAMLN}. For x > 8, the implementation is based on + *

+ * + * + * @param x Argument. + * @return the value of {@code log(Gamma(x))}, {@code Double.NaN} if + * {@code x <= 0.0}. + */ + public static double logGamma(double x) { + double ret; + + if (Double.isNaN(x) || (x <= 0.0)) { + ret = Double.NaN; + } else if (x < 0.5) { + return logGamma1p(x) - Math.log(x); + } else if (x <= 2.5) { + return logGamma1p((x - 0.5) - 0.5); + } else if (x <= 8.0) { + final int n = (int) Math.floor(x - 1.5); + double prod = 1.0; + for (int i = 1; i <= n; i++) { + prod *= x - i; + } + return logGamma1p(x - (n + 1)) + Math.log(prod); + } else { + double sum = lanczos(x); + double tmp = x + LANCZOS_G + .5; + ret = ((x + .5) * Math.log(tmp)) - tmp + + HALF_LOG_2_PI + Math.log(sum / x); + } + + return ret; + } + + /** + * Returns the regularized gamma function P(a, x). + * + * The implementation of this method is based on: + * + * + * @param a the a parameter. + * @param x the value. + * @param epsilon When the absolute value of the nth item in the + * series is less than epsilon the approximation ceases to calculate + * further elements in the series. + * @param maxIterations Maximum number of "iterations" to complete. + * @return the regularized gamma function P(a, x) + * @throws IllegalStateException if the algorithm fails to converge. + */ + public static double regularizedGammaP(double a, + double x, + double epsilon, + int maxIterations) { + double ret; + + if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) { + ret = Double.NaN; + } else if (x == 0.0) { + ret = 0.0; + } else if (x >= a + 1) { + // use regularizedGammaQ because it should converge faster in this + // case. + ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations); + } else { + // calculate series + double n = 0.0; // current element index + double an = 1.0 / a; // n-th element in the series + double sum = an; // partial sum + while (Math.abs(an/sum) > epsilon && + n < maxIterations && + sum < Double.POSITIVE_INFINITY) { + // compute next element in the series + n += 1.0; + an *= x / (a + n); + + // update partial sum + sum += an; + } + if (n >= maxIterations) { + throw new IllegalStateException("max count exceeded: " + maxIterations); + } else if (Double.isInfinite(sum)) { + ret = 1.0; + } else { + ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * sum; + } + } + + return ret; + } + + /** + * Returns the regularized gamma function Q(a, x) = 1 - P(a, x). + * + * The implementation of this method is based on: + * + * + * @param a the a parameter. + * @param x the value. + * @param epsilon When the absolute value of the nth item in the + * series is less than epsilon the approximation ceases to calculate + * further elements in the series. + * @param maxIterations Maximum number of "iterations" to complete. + * @return the regularized gamma function P(a, x) + * @throws IllegalStateException if the algorithm fails to converge. + */ + public static double regularizedGammaQ(final double a, + double x, + double epsilon, + int maxIterations) { + double ret; + + if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) { + ret = Double.NaN; + } else if (x == 0.0) { + ret = 1.0; + } else if (x < a + 1.0) { + // use regularizedGammaP because it should converge faster in this + // case. + ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations); + } else { + // create continued fraction + ContinuedFraction cf = new ContinuedFraction() { + + @Override + protected double getA(int n, double x) { + return ((2.0 * n) + 1.0) - a + x; + } + + @Override + protected double getB(int n, double x) { + return n * (a - n); + } + }; + + ret = 1.0 / cf.evaluate(x, epsilon, maxIterations); + ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * ret; + } + + return ret; + } + + /** + *

+ * Returns the Lanczos approximation used to compute the gamma function. + * The Lanczos approximation is related to the Gamma function by the + * following equation + *

+ * {@code gamma(x) = sqrt(2 * pi) / x * (x + g + 0.5) ^ (x + 0.5) + * * exp(-x - g - 0.5) * lanczos(x)}, + *
+ * where {@code g} is the Lanczos constant. + *

+ * + * @param x Argument. + * @return The Lanczos approximation. + * @see Lanczos Approximation + * equations (1) through (5), and Paul Godfrey's + * Note on the computation + * of the convergent Lanczos complex Gamma approximation + * @since 3.1 + */ + public static double lanczos(final double x) { + double sum = 0.0; + for (int i = LANCZOS.length - 1; i > 0; --i) { + sum += LANCZOS[i] / (x + i); + } + return sum + LANCZOS[0]; + } + + /** + * Returns the value of 1 / Γ(1 + x) - 1 for -0.5 ≤ x ≤ + * 1.5. This implementation is based on the double precision + * implementation in the NSWC Library of Mathematics Subroutines, + * {@code DGAM1}. + * + * @param x Argument. + * @return The value of {@code 1.0 / Gamma(1.0 + x) - 1.0}. + * @since 3.1 + */ + public static double invGamma1pm1(final double x) { + assert x >= -0.5; + assert x <= 1.5; + + final double ret; + final double t = x <= 0.5 ? x : (x - 0.5) - 0.5; + if (t < 0.0) { + final double a = INV_GAMMA1P_M1_A0 + t * INV_GAMMA1P_M1_A1; + double b = INV_GAMMA1P_M1_B8; + b = INV_GAMMA1P_M1_B7 + t * b; + b = INV_GAMMA1P_M1_B6 + t * b; + b = INV_GAMMA1P_M1_B5 + t * b; + b = INV_GAMMA1P_M1_B4 + t * b; + b = INV_GAMMA1P_M1_B3 + t * b; + b = INV_GAMMA1P_M1_B2 + t * b; + b = INV_GAMMA1P_M1_B1 + t * b; + b = 1.0 + t * b; + + double c = INV_GAMMA1P_M1_C13 + t * (a / b); + c = INV_GAMMA1P_M1_C12 + t * c; + c = INV_GAMMA1P_M1_C11 + t * c; + c = INV_GAMMA1P_M1_C10 + t * c; + c = INV_GAMMA1P_M1_C9 + t * c; + c = INV_GAMMA1P_M1_C8 + t * c; + c = INV_GAMMA1P_M1_C7 + t * c; + c = INV_GAMMA1P_M1_C6 + t * c; + c = INV_GAMMA1P_M1_C5 + t * c; + c = INV_GAMMA1P_M1_C4 + t * c; + c = INV_GAMMA1P_M1_C3 + t * c; + c = INV_GAMMA1P_M1_C2 + t * c; + c = INV_GAMMA1P_M1_C1 + t * c; + c = INV_GAMMA1P_M1_C + t * c; + if (x > 0.5) { + ret = t * c / x; + } else { + ret = x * ((c + 0.5) + 0.5); + } + } else { + double p = INV_GAMMA1P_M1_P6; + p = INV_GAMMA1P_M1_P5 + t * p; + p = INV_GAMMA1P_M1_P4 + t * p; + p = INV_GAMMA1P_M1_P3 + t * p; + p = INV_GAMMA1P_M1_P2 + t * p; + p = INV_GAMMA1P_M1_P1 + t * p; + p = INV_GAMMA1P_M1_P0 + t * p; + + double q = INV_GAMMA1P_M1_Q4; + q = INV_GAMMA1P_M1_Q3 + t * q; + q = INV_GAMMA1P_M1_Q2 + t * q; + q = INV_GAMMA1P_M1_Q1 + t * q; + q = 1.0 + t * q; + + double c = INV_GAMMA1P_M1_C13 + (p / q) * t; + c = INV_GAMMA1P_M1_C12 + t * c; + c = INV_GAMMA1P_M1_C11 + t * c; + c = INV_GAMMA1P_M1_C10 + t * c; + c = INV_GAMMA1P_M1_C9 + t * c; + c = INV_GAMMA1P_M1_C8 + t * c; + c = INV_GAMMA1P_M1_C7 + t * c; + c = INV_GAMMA1P_M1_C6 + t * c; + c = INV_GAMMA1P_M1_C5 + t * c; + c = INV_GAMMA1P_M1_C4 + t * c; + c = INV_GAMMA1P_M1_C3 + t * c; + c = INV_GAMMA1P_M1_C2 + t * c; + c = INV_GAMMA1P_M1_C1 + t * c; + c = INV_GAMMA1P_M1_C0 + t * c; + + if (x > 0.5) { + ret = (t / x) * ((c - 0.5) - 0.5); + } else { + ret = x * c; + } + } + + return ret; + } + + /** + * Returns the value of log Γ(1 + x) for -0.5 ≤ x ≤ 1.5. + * This implementation is based on the double precision implementation in + * the NSWC Library of Mathematics Subroutines, {@code DGMLN1}. + * + * @param x Argument. + * @return The value of {@code log(Gamma(1 + x))}. + * @since 3.1 + */ + public static double logGamma1p(final double x) { + assert x >= -0.5; + assert x <= 1.5; + return -Math.log1p(invGamma1pm1(x)); + } + +} + diff --git a/src/main/java/net/openhft/chronicle/hash/impl/util/math/PoissonDistribution.java b/src/main/java/net/openhft/chronicle/hash/impl/util/math/PoissonDistribution.java new file mode 100644 index 000000000..f6944bad4 --- /dev/null +++ b/src/main/java/net/openhft/chronicle/hash/impl/util/math/PoissonDistribution.java @@ -0,0 +1,162 @@ +/* + * Copyright (C) 2015 higherfrequencytrading.com + * Copyright 2001-2015 The Apache Software Foundation + * Copyright 2010-2012 CS Systèmes d'Information + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +package net.openhft.chronicle.hash.impl.util.math; + + +public class PoissonDistribution { + + private static final double EPSILON = 1e-12; + /** + * Poisson distribution is used to estimate segment fillings. Segments are not bound with + * Integer.MAX_VALUE, but it's not clear if algorithms from Commons Math could work with such + * big values as Long.MAX_VALUES + */ + private static final long UPPER_BOUND = 1L << 36; + private static final int MAX_ITERATIONS = 10000000; + + public static double cumulativeProbability(double mean, long x) { + if (x < 0) { + return 0; + } + if (x >= UPPER_BOUND) { + return 1; + } + return Gamma.regularizedGammaQ((double) x + 1, mean, EPSILON, MAX_ITERATIONS); + } + + public static long inverseCumulativeProbability(double mean, double p) { + checkProbability(p); + + long lower = 0; + if (p == 0.0) { + return lower; + } + lower -= 1; // this ensures cumulativeProbability(lower) < p, which + // is important for the solving step + + long upper = UPPER_BOUND; + if (p == 1.0) { + return upper; + } + + // use the one-sided Chebyshev inequality to narrow the bracket + // cf. AbstractRealDistribution.inverseCumulativeProbability(double) + final double mu = mean; + // in Poisson distribution, variance == mean + double variance = mean; + final double sigma = Math.sqrt(variance); + final boolean chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) || + Double.isInfinite(sigma) || Double.isNaN(sigma) || sigma == 0.0); + if (chebyshevApplies) { + double k = Math.sqrt((1.0 - p) / p); + double tmp = mu - k * sigma; + if (tmp > lower) { + lower = ((int) Math.ceil(tmp)) - 1; + } + k = 1.0 / k; + tmp = mu + k * sigma; + if (tmp < upper) { + upper = ((int) Math.ceil(tmp)) - 1; + } + } + + return solveInverseCumulativeProbability(mean, p, lower, upper); + } + + private static void checkProbability(double p) { + if (p < 0.0 || p > 1.0) { + throw new IllegalArgumentException("probability should be in [0.0, 1.0] bounds, " + p + + " given"); + } + } + + /** + * This is a utility function used by {@link + * #inverseCumulativeProbability(double, double)}. It assumes {@code 0 < p < 1} and + * that the inverse cumulative probability lies in the bracket {@code + * (lower, upper]}. The implementation does simple bisection to find the + * smallest {@code p}-quantile inf{x in Z | P(X<=x) >= p}. + * + * @param p the cumulative probability + * @param lower a value satisfying {@code cumulativeProbability(lower) < p} + * @param upper a value satisfying {@code p <= cumulativeProbability(upper)} + * @return the smallest {@code p}-quantile of this distribution + */ + private static long solveInverseCumulativeProbability(double mean, final double p, + long lower, long upper) { + + while (lower + 1 < upper) { + long xm = (lower + upper) / 2; + if (xm < lower || xm > upper) { + /* + * Overflow. + * There will never be an overflow in both calculation methods + * for xm at the same time + */ + xm = lower + (upper - lower) / 2; + } + + double pm = checkedCumulativeProbability(mean, xm); + if (pm >= p) { + upper = xm; + } else { + lower = xm; + } + } + return upper; + } + + public static double meanByCumulativeProbabilityAndValue(double p, long x, double precision) { + checkProbability(p); + assert x > 0 && x < UPPER_BOUND; + + double lower = 0; + double upper = UPPER_BOUND; + while (lower + precision < upper) { + double m = (lower + upper) / 2; + double pm = checkedCumulativeProbability(m, x); + if (pm >= p) { + upper = m; + } else { + lower = m; + } + } + return lower; + } + + /** + * Computes the cumulative probability function and checks for {@code NaN} + * values returned. Throws {@code MathInternalError} if the value is + * {@code NaN}. Rethrows any exception encountered evaluating the cumulative + * probability function. Throws {@code MathInternalError} if the cumulative + * probability function returns {@code NaN}. + * + * @param argument input value + * @return the cumulative probability + * @throws AssertionError if the cumulative probability is {@code NaN} + */ + private static double checkedCumulativeProbability(double mean, long argument) { + double result = cumulativeProbability(mean, argument); + if (Double.isNaN(result)) { + throw new AssertionError("Discrete cumulative probability function returned NaN " + + "for argument " + argument); + } + return result; + } +} diff --git a/src/main/java/net/openhft/chronicle/hash/impl/util/math/Precision.java b/src/main/java/net/openhft/chronicle/hash/impl/util/math/Precision.java new file mode 100644 index 000000000..6e8f2aef3 --- /dev/null +++ b/src/main/java/net/openhft/chronicle/hash/impl/util/math/Precision.java @@ -0,0 +1,98 @@ +/* + * Copyright (C) 2015 higherfrequencytrading.com + * Copyright 2001-2015 The Apache Software Foundation + * Copyright 2010-2012 CS Systèmes d'Information + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +package net.openhft.chronicle.hash.impl.util.math; + +class Precision { + + /** Offset to order signed double numbers lexicographically. */ + private static final long SGN_MASK = 0x8000000000000000L; + /** Positive zero bits. */ + private static final long POSITIVE_ZERO_DOUBLE_BITS = Double.doubleToRawLongBits(+0.0); + /** Negative zero bits. */ + private static final long NEGATIVE_ZERO_DOUBLE_BITS = Double.doubleToRawLongBits(-0.0); + + /** + * Returns {@code true} if there is no double value strictly between the + * arguments or the difference between them is within the range of allowed + * error (inclusive). + * + * @param x First value. + * @param y Second value. + * @param eps Amount of allowed absolute error. + * @return {@code true} if the values are two adjacent floating point + * numbers or they are within range of each other. + */ + public static boolean equals(double x, double y, double eps) { + return equals(x, y, 1) || Math.abs(y - x) <= eps; + } + + /** + * Returns true if both arguments are equal or within the range of allowed + * error (inclusive). + *

+ * Two float numbers are considered equal if there are {@code (maxUlps - 1)} + * (or fewer) floating point numbers between them, i.e. two adjacent + * floating point numbers are considered equal. + *

+ *

+ * Adapted from + * Bruce Dawson + *

+ * + * @param x first value + * @param y second value + * @param maxUlps {@code (maxUlps - 1)} is the number of floating point + * values between {@code x} and {@code y}. + * @return {@code true} if there are fewer than {@code maxUlps} floating + * point values between {@code x} and {@code y}. + */ + public static boolean equals(final double x, final double y, final int maxUlps) { + + final long xInt = Double.doubleToRawLongBits(x); + final long yInt = Double.doubleToRawLongBits(y); + + final boolean isEqual; + if (((xInt ^ yInt) & SGN_MASK) == 0l) { + // number have same sign, there is no risk of overflow + isEqual = Math.abs(xInt - yInt) <= maxUlps; + } else { + // number have opposite signs, take care of overflow + final long deltaPlus; + final long deltaMinus; + if (xInt < yInt) { + deltaPlus = yInt - POSITIVE_ZERO_DOUBLE_BITS; + deltaMinus = xInt - NEGATIVE_ZERO_DOUBLE_BITS; + } else { + deltaPlus = xInt - POSITIVE_ZERO_DOUBLE_BITS; + deltaMinus = yInt - NEGATIVE_ZERO_DOUBLE_BITS; + } + + if (deltaPlus > maxUlps) { + isEqual = false; + } else { + isEqual = deltaMinus <= (maxUlps - deltaPlus); + } + + } + + return isEqual && !Double.isNaN(x) && !Double.isNaN(y); + + } +} diff --git a/src/main/java/net/openhft/chronicle/hash/impl/util/math/package-info.java b/src/main/java/net/openhft/chronicle/hash/impl/util/math/package-info.java new file mode 100644 index 000000000..632e6bf3e --- /dev/null +++ b/src/main/java/net/openhft/chronicle/hash/impl/util/math/package-info.java @@ -0,0 +1,22 @@ +/* + * Copyright (C) 2015 higherfrequencytrading.com + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + */ + +/** + * All classes in this package are copied from Apache Commons Math 3.5 + * with minor simplifications, removing unused methods and replacing usages + * of {@code FastMath} with vanilla {@link java.lang.Math} + */ +package net.openhft.chronicle.hash.impl.util.math; \ No newline at end of file diff --git a/src/main/java/net/openhft/chronicle/map/ChronicleMapBuilder.java b/src/main/java/net/openhft/chronicle/map/ChronicleMapBuilder.java index 86a473c08..c27296116 100755 --- a/src/main/java/net/openhft/chronicle/map/ChronicleMapBuilder.java +++ b/src/main/java/net/openhft/chronicle/map/ChronicleMapBuilder.java @@ -19,6 +19,7 @@ package net.openhft.chronicle.map; import net.openhft.chronicle.hash.ChronicleHashBuilder; import net.openhft.chronicle.hash.ChronicleHashBuilderPrivateAPI; import net.openhft.chronicle.hash.ChronicleHashInstanceBuilder; +import net.openhft.chronicle.hash.impl.util.math.PoissonDistribution; import net.openhft.chronicle.hash.replication.*; import net.openhft.chronicle.hash.serialization.*; import net.openhft.chronicle.hash.serialization.internal.MetaBytesWriter; @@ -159,6 +160,7 @@ public final class ChronicleMapBuilder implements private int metaDataBytes = 0; private double maxBloatFactor = 1.0; private boolean allowSegmentTiering = true; + private double nonTieredSegmentsPercentile = 0.99999; private boolean putReturnsNull = false; private boolean removeReturnsNull = false; @@ -221,13 +223,15 @@ public final class ChronicleMapBuilder implements return roundUp; } - private static void checkSegments(int segments) { - if (segments <= 0 || segments > MAX_SEGMENTS) + private static void checkSegments(long segments) { + if (segments <= 0) { throw new IllegalArgumentException("segments should be positive, " + segments + " given"); - if (segments > MAX_SEGMENTS) + } + if (segments > MAX_SEGMENTS) { throw new IllegalArgumentException("Max segments is " + MAX_SEGMENTS + ", " + segments + " given"); + } } private static long divideUpper(long dividend, long divisor) { @@ -750,8 +754,9 @@ public final class ChronicleMapBuilder implements entriesPerSegment = this.entriesPerSegment; } else { int actualSegments = actualSegments(replicated); - long totalEntries = totalEntriesIfPoorDistribution(actualSegments); - entriesPerSegment = divideUpper(totalEntries, actualSegments); + double averageEntriesPerSegment = entries() * 1.0 / actualSegments; + entriesPerSegment = PoissonDistribution.inverseCumulativeProbability( + averageEntriesPerSegment, nonTieredSegmentsPercentile); } if (actualChunksPerSegment > 0) return entriesPerSegment; @@ -809,28 +814,6 @@ public final class ChronicleMapBuilder implements return round(entriesPerSegment * averageChunksPerEntry(replicated)); } - private long totalEntriesIfPoorDistribution(long segments) { - if (segments == 1) - return entries(); - long entries = entries(); - // check if the number of entries is small compared with the number of segments. - long s3 = Math.min(8, segments) * Math.min(32, segments) * Math.min(128, segments); - if (entries * 4 <= s3) - entries *= 1.8; - else if (entries <= s3) - entries *= 1.45; - else if (entries <= s3 * 4) - entries *= 1.33; - else if (entries <= s3 * 8) - entries *= 1.22; - else if (entries <= s3 * 16) - entries *= 1.15; - else - entries *= 1.1; - return Math.min(segments * entries(), - entries + 4 * segments + 8); - } - @Override public ChronicleMapBuilder minSegments(int minSegments) { checkSegments(minSegments); @@ -892,16 +875,7 @@ public final class ChronicleMapBuilder implements if (actualSegments > 0) return actualSegments; if (entriesPerSegment > 0) { - long segments = 1; - for (int i = 0; i < 3; i++) { // try to converge - long totalEntries = totalEntriesIfPoorDistribution((int) segments); - segments = divideUpper(totalEntries, entriesPerSegment); - if (segments > MAX_SEGMENTS) - throw new IllegalStateException(); - } - if (minSegments > 0) - segments = Math.max(minSegments, segments); - return (int) segments; + return (int) segmentsGivenEntriesPerSegmentFixed(entriesPerSegment, replicated); } // iteratively try to fit 3..8 bytes per hash lookup slot. Trying to apply small slot // sizes (=> segment sizes, because slot size depends on segment size) not only because @@ -971,13 +945,23 @@ public final class ChronicleMapBuilder implements } private long trySegments(long entriesPerSegment, int maxSegments, boolean replicated) { - long firstSegmentsApproximation = Math.max(minSegments(), entries() / entriesPerSegment); - long adjustedEntries = totalEntriesIfPoorDistribution(firstSegmentsApproximation); - long segments = divideUpper(adjustedEntries, entriesPerSegment); + long segments = segmentsGivenEntriesPerSegmentFixed(entriesPerSegment, replicated); segments = Maths.nextPower2(Math.max(segments, minSegments()), 1L); return segments <= maxSegments ? segments : -segments; } + private long segmentsGivenEntriesPerSegmentFixed(long entriesPerSegment, boolean replicated) { + double precision = 1.0 / averageChunksPerEntry(replicated); + double entriesPerSegmentShouldBe = + PoissonDistribution.meanByCumulativeProbabilityAndValue( + nonTieredSegmentsPercentile, entriesPerSegment, precision); + long segments = ((long) (entries() / entriesPerSegmentShouldBe)) + 1; + checkSegments(segments); + if (minSegments > 0) + segments = Math.max(minSegments, segments); + return segments; + } + int segmentHeaderSize(boolean replicated) { int segments = actualSegments(replicated); @@ -1040,6 +1024,18 @@ public final class ChronicleMapBuilder implements return this; } + @Override + public ChronicleMapBuilder nonTieredSegmentsPercentile( + double nonTieredSegmentsPercentile) { + if (Double.isNaN(nonTieredSegmentsPercentile) || + 0.5 <= nonTieredSegmentsPercentile || nonTieredSegmentsPercentile >= 1.0) { + throw new IllegalArgumentException("nonTieredSegmentsPercentile should be in (0.5, " + + "1.0) range, " + nonTieredSegmentsPercentile + " is given"); + } + this.nonTieredSegmentsPercentile = nonTieredSegmentsPercentile; + return this; + } + long maxExtraTiers(boolean replicated) { if (!allowSegmentTiering) return 0; diff --git a/src/main/java/net/openhft/chronicle/set/ChronicleSetBuilder.java b/src/main/java/net/openhft/chronicle/set/ChronicleSetBuilder.java index 750e340dd..cb6f19f16 100755 --- a/src/main/java/net/openhft/chronicle/set/ChronicleSetBuilder.java +++ b/src/main/java/net/openhft/chronicle/set/ChronicleSetBuilder.java @@ -177,6 +177,12 @@ public final class ChronicleSetBuilder return this; } + @Override + public ChronicleSetBuilder nonTieredSegmentsPercentile(double nonTieredSegmentsPercentile) { + chronicleMapBuilder.nonTieredSegmentsPercentile(nonTieredSegmentsPercentile); + return this; + } + @Override public String toString() { return " ChronicleSetBuilder{" +