From 24a0566fd4ed8668e720bc39ad253e5a1205f87d Mon Sep 17 00:00:00 2001 From: "joseph.lizier" Date: Fri, 14 Nov 2014 10:24:53 +0000 Subject: [PATCH] Patching bug in digamma function which is susceptible to race condition in multi-threaded applications. Also explicitly setting return array values in Kraskov calculators. --- ...tualInfoCalculatorMultiVariateKraskov.java | 8 +++- ...ualInfoCalculatorMultiVariateKraskov1.java | 16 +++++-- .../source/infodynamics/utils/MathsUtils.java | 43 +++++++++++++------ 3 files changed, 49 insertions(+), 18 deletions(-) diff --git a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java index 6783c49..8c4216a 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java @@ -513,7 +513,9 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov System.out.printf("Av = digamma(k)=%.3f + =%.3f = %.3f \n", MathsUtils.digamma(k), averageDiGammas, MathsUtils.digamma(k) + averageDiGammas); } - return new double[] { MathsUtils.digamma(k) + averageDiGammas }; + double[] result = new double[1]; + result[0] = MathsUtils.digamma(k) + averageDiGammas; + return result; } else { // Algorithm 2: // We also retrieve the sums of inverses for debugging purposes: @@ -531,7 +533,9 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov (2.0 / (double)k), averageMeasure, averageInverseCountInJointYZ, averageInverseCountInJointXZ); } - return new double[] { averageMeasure }; + double[] result = new double[1]; + result[0] = averageMeasure; + return result; } } } diff --git a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java index 22105d6..07e0994 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java +++ b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov1.java @@ -204,8 +204,10 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1 if (returnLocals) { localCondMi[t-startTimePoint] = digammaK - digammaNxzPlusOne - digammaNyzPlusOne + digammaNzPlusOne; if (debug) { - System.out.printf("t=%d, n_xz=%d, n_yz=%d, n_z=%d, local=%.4f\n", - t, n_xz, n_yz, n_z, localCondMi[t-startTimePoint]); + System.out.printf("t=%d, n_xz=%d, n_yz=%d, n_z=%d, local=%.4f," + + " digamma(n_xz+1)=%.5f, digamma(n_yz+1)=%.5f, digamma(n_z+1)=%.5f, \n", + t, n_xz, n_yz, n_z, localCondMi[t-startTimePoint], + digammaNxzPlusOne, digammaNyzPlusOne, digammaNzPlusOne); } } } @@ -221,6 +223,9 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1 System.out.println("\tz search: " + (conditionalTime/1000.0)); System.out.println("\tzx search: " + (conditionalXTime/1000.0)); System.out.println("\tzy search: " + (conditionalYTime/1000.0)); + System.out.printf("%d:%d -- Returning: %.4f, %.4f, %.4f, %.4f\n", + startTimePoint, (startTimePoint + numTimePoints), + sumDiGammas, sumNxz, sumNyz, sumNz); } // Select what to return: @@ -229,7 +234,12 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov1 } else { // Pad return array with two values, to allow compatibility in // return length with algorithm 2 - return new double[] {sumDiGammas, sumNxz, sumNyz, sumNz, 0, 0}; + double[] results = new double[6]; + results[0] = sumDiGammas; + results[1] = sumNxz; + results[2] = sumNyz; + results[3] = sumNz; + return results; } } } diff --git a/java/source/infodynamics/utils/MathsUtils.java b/java/source/infodynamics/utils/MathsUtils.java index 8f6b2b0..bb40844 100755 --- a/java/source/infodynamics/utils/MathsUtils.java +++ b/java/source/infodynamics/utils/MathsUtils.java @@ -259,10 +259,10 @@ public class MathsUtils { } /** - * Compute digamma(d). + *

Compute digamma(d).

* - * Stores previous calculations to speed up computation here, though some precision may - * be lost because we're adding in larger numbers first. + *

Stores previous calculations to speed up computation here, though some precision may + * be lost because we're adding in larger numbers first.

* * @param d * @return @@ -272,13 +272,15 @@ public class MathsUtils { if (d < 1) { return Double.NaN; } - if (highestDigammaArgCalced == 0) { - // We're not initialised yet: - storedDigammas[0] = Double.NaN; - storedDigammas[1] = -EULER_MASCHERONI_CONSTANT; - highestDigammaArgCalced = 1; - } - if (d <= highestDigammaArgCalced) { + // Need to take care with multi-threaded applications + // here -- another thread could be attempting to update + // highestDigammaArgCalced at the same time. + // Sample the value of highestDigammaArgCalced and + // run with it through the rest of the method, to avoid + // any race conditions if it is concurrently updated: + int highestDigammaArgCalcedAtStartForThisThread = + highestDigammaArgCalced; + if (d <= highestDigammaArgCalcedAtStartForThisThread) { // We've already calculated this one return storedDigammas[d]; } @@ -287,14 +289,29 @@ public class MathsUtils { // directly use commons.math: return Gamma.digamma(d); } + if (highestDigammaArgCalcedAtStartForThisThread == 0) { + // We're not initialised yet: + storedDigammas[0] = Double.NaN; + storedDigammas[1] = -EULER_MASCHERONI_CONSTANT; + highestDigammaArgCalcedAtStartForThisThread = 1; + } // Else we'll calculate it and update the storage: - double result = storedDigammas[highestDigammaArgCalced]; - for (int n = highestDigammaArgCalced + 1; n <= d; n++) { + double result = storedDigammas[highestDigammaArgCalcedAtStartForThisThread]; + for (int n = highestDigammaArgCalcedAtStartForThisThread + 1; n <= d; n++) { result += 1.0 / (double) (n-1); // n must be < NUM_STORED_DIGAMMAS by earlier if statement on d storedDigammas[n] = result; } - highestDigammaArgCalced = d; + // Finally, update highestDigammaArgCalced: + // We could synchronize on highestDigammaArgCalced, however + // this is costly; the following check may be compromised under + // a race condition, but the worst outcome is that we have to + // recalculate a few values -- this is probably faster than + // enforcing synchronisation + if (d > highestDigammaArgCalced) { + // We make the check in case another thread has already set this higher + highestDigammaArgCalced = d; + } return result; }