Index: trunk/sources/HeuristicLab.Problems.DataAnalysis/3.4/Implementation/Classification/ThresholdCalculators/NormalDistributionCutPointsThresholdCalculator.cs
===================================================================
--- trunk/sources/HeuristicLab.Problems.DataAnalysis/3.4/Implementation/Classification/ThresholdCalculators/NormalDistributionCutPointsThresholdCalculator.cs	(revision 8637)
+++ trunk/sources/HeuristicLab.Problems.DataAnalysis/3.4/Implementation/Classification/ThresholdCalculators/NormalDistributionCutPointsThresholdCalculator.cs	(revision 8638)
@@ -53,7 +53,6 @@
 
     public static void CalculateThresholds(IClassificationProblemData problemData, IEnumerable<double> estimatedValues, IEnumerable<double> targetClassValues, out double[] classValues, out double[] thresholds) {
-      double maxEstimatedValue = estimatedValues.Max();
-      double minEstimatedValue = estimatedValues.Min();
       var estimatedTargetValues = Enumerable.Zip(estimatedValues, targetClassValues, (e, t) => new { EstimatedValue = e, TargetValue = t }).ToList();
+      double estimatedValuesRange = estimatedValues.Range();
 
       Dictionary<double, double> classMean = new Dictionary<double, double>();
@@ -82,6 +81,13 @@
           // calculate all thresholds
           CalculateCutPoints(classMean[class0], classStdDev[class0], classMean[class1], classStdDev[class1], out x1, out x2);
-          if (!thresholdList.Any(x => x.IsAlmost(x1))) thresholdList.Add(x1);
-          if (!thresholdList.Any(x => x.IsAlmost(x2))) thresholdList.Add(x2);
+
+          // if the two cut points are too close (for instance because the stdDev=0)
+          // then move them by 0.1% of the range of estimated values
+          if (x1.IsAlmost(x2)) {
+            x1 -= 0.001 * estimatedValuesRange;
+            x2 += 0.001 * estimatedValuesRange;
+          }
+          if (!double.IsInfinity(x1) && !thresholdList.Any(x => x.IsAlmost(x1))) thresholdList.Add(x1);
+          if (!double.IsInfinity(x2) && !thresholdList.Any(x => x.IsAlmost(x2))) thresholdList.Add(x2);
         }
       }
@@ -92,26 +98,34 @@
       // all points in the partition are classified as the class with the maximal density in the parition
       List<double> classValuesList = new List<double>();
-      for (int i = 0; i < thresholdList.Count; i++) {
-        double m;
-        if (double.IsNegativeInfinity(thresholdList[i])) {
-          m = thresholdList[i + 1] - 1.0; // smaller than the smalles non-infinity threshold
-        } else if (i == thresholdList.Count - 1) {
-          // last threshold
-          m = thresholdList[i] + 1.0; // larger than the last threshold
-        } else {
-          m = thresholdList[i] + (thresholdList[i + 1] - thresholdList[i]) / 2.0; // middle of partition
+      if (thresholdList.Count == 1) {
+        // this happens if there are no thresholds (distributions for all classes are exactly the same)
+        // -> all samples should be classified as the first class
+        classValuesList.Add(originalClasses[0]);
+      } else {
+        // at least one reasonable threshold ...
+        // find the most likely class for the points between thresholds m
+        for (int i = 0; i < thresholdList.Count; i++) {
+          double m;
+          if (double.IsNegativeInfinity(thresholdList[i])) {
+            m = thresholdList[i + 1] - 1.0; // smaller than the smallest non-infinity threshold
+          } else if (i == thresholdList.Count - 1) {
+            // last threshold
+            m = thresholdList[i] + 1.0; // larger than the last threshold
+          } else {
+            m = thresholdList[i] + (thresholdList[i + 1] - thresholdList[i]) / 2.0; // middle of partition
+          }
+
+          // determine class with maximal probability density in m
+          double maxDensity = LogNormalDensity(m, classMean[originalClasses[0]], classStdDev[originalClasses[0]]);
+          double maxDensityClassValue = originalClasses[0];
+          foreach (var classValue in originalClasses.Skip(1)) {
+            double density = LogNormalDensity(m, classMean[classValue], classStdDev[classValue]);
+            if (density > maxDensity) {
+              maxDensity = density;
+              maxDensityClassValue = classValue;
+            }
+          }
+          classValuesList.Add(maxDensityClassValue);
         }
-
-        // determine class with maximal probability density in m
-        double maxDensity = double.MinValue;
-        double maxDensityClassValue = -1;
-        foreach (var classValue in originalClasses) {
-          double density = LogNormalDensity(m, classMean[classValue], classStdDev[classValue]);
-          if (density > maxDensity) {
-            maxDensity = density;
-            maxDensityClassValue = classValue;
-          }
-        }
-        classValuesList.Add(maxDensityClassValue);
       }
 
@@ -130,5 +144,5 @@
       filteredClassValues.Add(classValuesList[0]);
       for (int i = 0; i < classValuesList.Count - 1; i++) {
-        if (classValuesList[i] != classValuesList[i + 1]) {
+        if (!classValuesList[i].IsAlmost(classValuesList[i + 1])) {
           filteredThresholds.Add(thresholdList[i + 1]);
           filteredClassValues.Add(classValuesList[i + 1]);
@@ -140,11 +154,36 @@
 
     private static double LogNormalDensity(double x, double mu, double sigma) {
+      if (sigma.IsAlmost(0.0))
+        if (mu.IsAlmost(x)) return 0.0; // (log(1))
+        else return double.NegativeInfinity;
       return -0.5 * Math.Log(2.0 * Math.PI * sigma * sigma) - ((x - mu) * (x - mu)) / (2.0 * sigma * sigma);
     }
 
     private static void CalculateCutPoints(double m1, double s1, double m2, double s2, out double x1, out double x2) {
-      double a = (s1 * s1 - s2 * s2);
-      x1 = -(-m2 * s1 * s1 + m1 * s2 * s2 + Math.Sqrt(s1 * s1 * s2 * s2 * ((m1 - m2) * (m1 - m2) + 2.0 * (-s1 * s1 + s2 * s2) * Math.Log(s2 / s1)))) / a;
-      x2 = (m2 * s1 * s1 - m1 * s2 * s2 + Math.Sqrt(s1 * s1 * s2 * s2 * ((m1 - m2) * (m1 - m2) + 2.0 * (-s1 * s1 + s2 * s2) * Math.Log(s2 / s1)))) / a;
+      if (s1.IsAlmost(s2)) {
+        if (m1.IsAlmost(m2)) {
+          x1 = double.NegativeInfinity;
+          x2 = double.NegativeInfinity;
+        } else {
+          x1 = (m1 + m2) / 2;
+          x2 = double.NegativeInfinity;
+        }
+      } else if (s1.IsAlmost(0.0)) {
+        x1 = m1;
+        x2 = m1;
+      } else if (s2.IsAlmost(0.0)) {
+        x1 = m2;
+        x2 = m2;
+      } else {
+        // scale s1 and s2 for numeric stability
+        s2 = s2 / s1;
+        s1 = 1.0;
+        double a = (s1 + s2) * (s1 - s2);
+        double g = Math.Sqrt(s1 * s1 * s2 * s2 * ((m1 - m2) * (m1 - m2) + 2.0 * (s1 * s1 + s2 * s2) * Math.Log(s2 / s1)));
+        double m1s2 = m1 * s2 * s2;
+        double m2s1 = m2 * s1 * s1;
+        x1 = -(-m2s1 + m1s2 + g) / a;
+        x2 = (m2s1 - m1s2 + g) / a;
+      }
     }
   }
Index: trunk/sources/HeuristicLab.Tests/HeuristicLab.Problems.DataAnalysis-3.4/ThresholdCalculatorsTest.cs
===================================================================
--- trunk/sources/HeuristicLab.Tests/HeuristicLab.Problems.DataAnalysis-3.4/ThresholdCalculatorsTest.cs	(revision 8638)
+++ trunk/sources/HeuristicLab.Tests/HeuristicLab.Problems.DataAnalysis-3.4/ThresholdCalculatorsTest.cs	(revision 8638)
@@ -0,0 +1,122 @@
+#region License Information
+/* HeuristicLab
+ * Copyright (C) 2002-2012 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
+ *
+ * This file is part of HeuristicLab.
+ *
+ * HeuristicLab is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * HeuristicLab 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 General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
+ */
+#endregion
+
+using HeuristicLab.Problems.DataAnalysis;
+using Microsoft.VisualStudio.TestTools.UnitTesting;
+namespace HeuristicLab.Problems.DataAnalysis_34.Tests {
+
+  [TestClass()]
+  public class ThresholdCalculatorsTest {
+    [TestMethod]
+    public void NormalDistributionCutPointsThresholdCalculatorTest() {
+
+      {
+        // simple two-class case
+        double[] estimatedValues = new double[] { 1.0, 0.99, 1.01, 2.0, 1.99, 2.01 };
+        double[] targetClassValues = new double[] { 0.0, 0.0, 0.0, 1.0, 1.0, 1.0 };
+        double[] classValues;
+        double[] thresholds;
+        NormalDistributionCutPointsThresholdCalculator.CalculateThresholds(null, estimatedValues, targetClassValues,
+                                                                           out classValues, out thresholds);
+
+        var expectedClassValues = new double[] { 0.0, 1.0 };
+        var expectedTresholds = new double[] { double.NegativeInfinity, 1.5 };
+
+        AssertEqual(expectedClassValues, classValues);
+        AssertEqual(expectedTresholds, thresholds);
+      }
+
+      {
+        // switched classes two-class case
+        double[] estimatedValues = new double[] { 1.0, 0.99, 1.01, 2.0, 1.99, 2.01 };
+        double[] targetClassValues = new double[] { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
+        double[] classValues;
+        double[] thresholds;
+        NormalDistributionCutPointsThresholdCalculator.CalculateThresholds(null, estimatedValues, targetClassValues,
+                                                                           out classValues, out thresholds);
+
+        var expectedClassValues = new double[] { 1.0, 0.0 };
+        var expectedTresholds = new double[] { double.NegativeInfinity, 1.5 };
+
+        AssertEqual(expectedClassValues, classValues);
+        AssertEqual(expectedTresholds, thresholds);
+      }
+
+      {
+        // three-class case with permutated estimated values
+        double[] estimatedValues = new double[] { 1.0, 0.99, 1.01, 2.0, 1.99, 2.01, -1.0, -0.99, -1.01 };
+        double[] targetClassValues = new double[] { 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0 };
+        double[] classValues;
+        double[] thresholds;
+        NormalDistributionCutPointsThresholdCalculator.CalculateThresholds(null, estimatedValues, targetClassValues,
+                                                                           out classValues, out thresholds);
+
+        var expectedClassValues = new double[] { 1.0, 2.0, 0.0 };
+        var expectedTresholds = new double[] { double.NegativeInfinity, 0.0, 1.5 };
+
+        AssertEqual(expectedClassValues, classValues);
+        AssertEqual(expectedTresholds, thresholds);
+      }
+
+      {
+        // constant output values for all classes
+        double[] estimatedValues = new double[] { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 };
+        double[] targetClassValues = new double[] { 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0 };
+        double[] classValues;
+        double[] thresholds;
+        NormalDistributionCutPointsThresholdCalculator.CalculateThresholds(null, estimatedValues, targetClassValues,
+                                                                           out classValues, out thresholds);
+
+        var expectedClassValues = new double[] { 0.0 };
+        var expectedTresholds = new double[] { double.NegativeInfinity };
+
+        AssertEqual(expectedClassValues, classValues);
+        AssertEqual(expectedTresholds, thresholds);
+      }
+
+      {
+        // constant output values for two of three classes
+        double[] estimatedValues = new double[] { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -0.99, -1.01 };
+        double[] targetClassValues = new double[] { 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0 };
+        double[] classValues;
+        double[] thresholds;
+        NormalDistributionCutPointsThresholdCalculator.CalculateThresholds(null, estimatedValues, targetClassValues,
+                                                                           out classValues, out thresholds);
+
+
+        var expectedClassValues = new double[] { 1.0, 0.0, 1.0 };
+        double range = 1.0 + 1.01;
+        var expectedTresholds = new double[] { double.NegativeInfinity, 1.0 - 0.001 * range, 1.0 + 0.001 * range };
+
+        AssertEqual(expectedClassValues, classValues);
+        AssertEqual(expectedTresholds, thresholds);
+      }
+
+    }
+
+
+    private static void AssertEqual(double[] expected, double[] actual) {
+      Assert.AreEqual(expected.Length, actual.Length);
+      for (int i = 0; i < expected.Length; i++)
+        Assert.AreEqual(expected[i], actual[i]);
+    }
+  }
+}
Index: trunk/sources/HeuristicLab.Tests/HeuristicLab.Tests.csproj
===================================================================
--- trunk/sources/HeuristicLab.Tests/HeuristicLab.Tests.csproj	(revision 8637)
+++ trunk/sources/HeuristicLab.Tests/HeuristicLab.Tests.csproj	(revision 8638)
@@ -357,4 +357,5 @@
     <Compile Include="HeuristicLab.PluginInfraStructure-3.3\InstallationManagerTest.cs" />
     <Compile Include="HeuristicLab.PluginInfraStructure-3.3\TypeDiscoveryTest.cs" />
+    <Compile Include="HeuristicLab.Problems.DataAnalysis-3.4\ThresholdCalculatorsTest.cs" />
     <Compile Include="HeuristicLab.Problems.DataAnalysis-3.4\OnlineCalculatorPerformanceTest.cs" />
     <Compile Include="HeuristicLab.Problems.DataAnalysis-3.4\StatisticCalculatorsTest.cs" />
