Patching bug in digamma function which is susceptible to race condition in multi-threaded applications. Also explicitly setting return array values in Kraskov calculators.

java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java

@@ -513,7 +513,9 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov
 					System.out.printf("Av = digamma(k)=%.3f + <digammas>=%.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;
 			}			
 		}
 	}

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;
 		}		
 	}
 }

java/source/infodynamics/utils/MathsUtils.java

@@ -259,10 +259,10 @@ public class MathsUtils {
 	}
 
 	/**
-	 * Compute digamma(d).
+	 * <p>Compute digamma(d).</p>
 	 * 
-	 * Stores previous calculations to speed up computation here, though some precision may
-	 *  be lost because we're adding in larger numbers first.
+	 * <p>Stores previous calculations to speed up computation here, though some precision may
+	 *  be lost because we're adding in larger numbers first.</p>
 	 * 
 	 * @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;
 	}