Default threshold bugs

.claude/settings.local.json

@@ -11,7 +11,9 @@
       "Bash(./HEC.FDA.TestingUtility.exe:*)",
       "Bash(./HEC.FDA.TestingUtility/bin/Debug/net9.0-windows/HEC.FDA.TestingUtility.exe:*)",
       "Bash(dir /s \"C:\\\\Programs\\\\Source\\\\HEC-FDA\\\\HEC.FDA.TestingUtility\\\\*.cs\")",
-      "Bash(findstr:*)"
+      "Bash(findstr:*)",
+      "Bash(dotnet run:*)",
+      "Bash(ls *.sln)"
     ],
     "deny": []
   }

HEC.FDA.Benchmarks/HEC.FDA.Benchmarks.csproj

@@ -11,6 +11,7 @@
 
 	<ItemGroup>
 		<ProjectReference Include="..\HEC.FDA.Statistics\HEC.FDA.Statistics.csproj" />
+		<ProjectReference Include="..\HEC.FDA.Model\HEC.FDA.Model.csproj" />
 	</ItemGroup>
 
 </Project>

HEC.FDA.Benchmarks/SumYsForGivenXBenchmarks.cs

@@ -0,0 +1,94 @@
+using BenchmarkDotNet.Attributes;
+using HEC.FDA.Model.paireddata;
+
+namespace HEC.FDA.Benchmarks;
+
+[MemoryDiagnoser]
+public class SumYsForGivenXBenchmarks
+{
+    private PairedData _subject = null!;
+    private PairedData _input = null!;
+
+    [Params("same_grid", "disjoint_grids", "input_sparse", "subject_sparse")]
+    public string Scenario { get; set; } = null!;
+
+    [GlobalSetup]
+    public void Setup()
+    {
+        switch (Scenario)
+        {
+            case "same_grid":
+            {
+                // Both curves have the same 200-point x-grid (typical stage-damage resolution)
+                double[] x = new double[200];
+                double[] ySubject = new double[200];
+                double[] yInput = new double[200];
+                for (int i = 0; i < 200; i++)
+                {
+                    x[i] = i * 0.1;
+                    ySubject[i] = i * 10.0;
+                    yInput[i] = i * 15.0;
+                }
+                _subject = new PairedData(x, ySubject);
+                _input = new PairedData((double[])x.Clone(), yInput);
+                break;
+            }
+            case "disjoint_grids":
+            {
+                // Two curves with completely different x-grids (worst case for union)
+                double[] xSubject = new double[200];
+                double[] ySubject = new double[200];
+                double[] xInput = new double[200];
+                double[] yInput = new double[200];
+                for (int i = 0; i < 200; i++)
+                {
+                    xSubject[i] = i * 0.2;       // 0.0, 0.2, 0.4, ...
+                    ySubject[i] = i * 10.0;
+                    xInput[i] = i * 0.2 + 0.1;   // 0.1, 0.3, 0.5, ...
+                    yInput[i] = i * 15.0;
+                }
+                _subject = new PairedData(xSubject, ySubject);
+                _input = new PairedData(xInput, yInput);
+                break;
+            }
+            case "input_sparse":
+            {
+                // Subject has 200 points, input has only 2 — the bug scenario
+                double[] xSubject = new double[200];
+                double[] ySubject = new double[200];
+                for (int i = 0; i < 200; i++)
+                {
+                    xSubject[i] = i * 0.1;
+                    ySubject[i] = i * 10.0;
+                }
+                double[] xInput = new double[] { 0.0, 19.9 };
+                double[] yInput = new double[] { 0.0, 2985.0 };
+                _subject = new PairedData(xSubject, ySubject);
+                _input = new PairedData(xInput, yInput);
+                break;
+            }
+            case "subject_sparse":
+            {
+                // Subject has only 2 points, input has 200
+                double[] xSubject = new double[] { 0.0, 19.9 };
+                double[] ySubject = new double[] { 0.0, 1990.0 };
+                double[] xInput = new double[200];
+                double[] yInput = new double[200];
+                for (int i = 0; i < 200; i++)
+                {
+                    xInput[i] = i * 0.1;
+                    yInput[i] = i * 15.0;
+                }
+                _subject = new PairedData(xSubject, ySubject);
+                _input = new PairedData(xInput, yInput);
+                break;
+            }
+        }
+    }
+
+    [Benchmark]
+    public PairedData SumYsForGivenX()
+    {
+        return _subject.SumYsForGivenX(_input);
+    }
+}

HEC.FDA.Model/compute/ImpactAreaScenarioSimulation.cs

@@ -151,7 +151,7 @@ private void SetupPerformanceThresholds(ConvergenceCriteria convergenceCriteria)
                 FrequencyStageCurves curves = GetFrequencyStageSample(computeIsDeterministic: true, 1);
                 //Floodplain stage because we're computing risk.
                 ComputeRiskFromStageFrequency(curves.FloodplainStage, 1, 1, computeIsDeterministic: true, _FailureStageDamageFunctions, ConsequenceType.Damage, false, true);
-                //default threshold represents a channel stage, because the stage damage functions are in channel stage. 
+                //default threshold represents a channel stage, because the stage damage functions are in channel stage.
                 Threshold defaultThreshold = ComputeDefaultThreshold(convergenceCriteria, damageFrequencyFunctions: _ImpactAreaScenarioResults.ConsequenceFrequencyFunctions.Select((c) => c.FrequencyCurve).ToList());
                 _ImpactAreaScenarioResults.PerformanceByThresholds.AddThreshold(defaultThreshold);
             }
@@ -712,9 +712,8 @@ private Threshold ComputeDefaultThreshold(ConvergenceCriteria convergenceCriteri
             }
             PairedData totalStageDamage = ComputeTotalStageDamage(_FailureStageDamageFunctions);
             PairedData totalFrequencyDamage = damageFrequencyFunctions[0];
-            for (int i = 0; i < damageFrequencyFunctions.Count; i++)
+            for (int i = 1; i < damageFrequencyFunctions.Count; i++)
             {
-                //Some unnecessary GC happening here. Not happening in the big parallel for. 10s to 100s , not hundreds of thousands.
                 totalFrequencyDamage = totalFrequencyDamage.SumYsForGivenX(damageFrequencyFunctions[i]);
             }
             double thresholdDamage = THRESHOLD_DAMAGE_PERCENT * totalFrequencyDamage.f(THRESHOLD_DAMAGE_RECURRENCE_INTERVAL);

HEC.FDA.Model/paireddata/PairedData.cs

@@ -196,45 +196,28 @@ public PairedData SumYsForGivenX(IPairedData inputPairedData)
             {
                 throw new ArgumentException("X values must be in increasing order.");
             }
-            double[] x = new double[inputPairedData.Xvals.Count];
-            double[] ySum = new double[inputPairedData.Yvals.Count];
+
+            // Build a sorted union of both x-grids
+            SortedSet<double> unionXSet = new SortedSet<double>();
+            for (int i = 0; i < Xvals.Count; i++)
+            {
+                unionXSet.Add(Xvals[i]);
+            }
             for (int i = 0; i < inputPairedData.Xvals.Count; i++)
             {
-                int index = Array.BinarySearch(_xVals, inputPairedData.Xvals[i]);
-                double yValueFromSubject;
-                if (index >= 0)
-                {
-                    //value matched exactly
-                    yValueFromSubject = Yvals[index];
-
-                }
-                else
-                {
-                    index = ~index;
-
-                    //if the x value from input is greater than all x values in subject
-                    if (index == Xvals.Count)
-                    {
-                        yValueFromSubject = Yvals[^1];
-                    }
-
-                    //if the x value from the input is less than all x values in subject
-                    else if (index == 0)
-                    {
-                        yValueFromSubject = Yvals[0];
-                    }
-                    else //x value from the input is within range of x values in subject, but does not match exactly
-                    { //Interpolate Y=mx+b
-                        double slope = (Yvals[index] - Yvals[index - 1]) / (Xvals[index] - Xvals[index - 1]);
-                        double intercept = Yvals[index] - slope * Xvals[index];
-                        yValueFromSubject = intercept + slope * inputPairedData.Xvals[i];
-                    }
-                }
-                ySum[i] = inputPairedData.Yvals[i] + yValueFromSubject;
-                x[i] = inputPairedData.Xvals[i];
+                unionXSet.Add(inputPairedData.Xvals[i]);
             }
-            return new PairedData(x, ySum);
 
+            double[] unionX = new double[unionXSet.Count];
+            double[] unionY = new double[unionXSet.Count];
+            int idx = 0;
+            foreach (double x in unionXSet)
+            {
+                unionX[idx] = x;
+                unionY[idx] = f(x) + inputPairedData.f(x);
+                idx++;
+            }
+            return new PairedData(unionX, unionY);
         }
 
 

HEC.FDA.ModelTest/unittests/DefaultThresholdShould.cs

@@ -0,0 +1,152 @@
+using HEC.FDA.Model.compute;
+using HEC.FDA.Model.metrics;
+using HEC.FDA.Model.paireddata;
+using Statistics;
+using Statistics.Distributions;
+using System.Collections.Generic;
+using System.Linq;
+using System.Threading;
+using Xunit;
+
+namespace HEC.FDA.ModelTest.unittests
+{
+    [Trait("RunsOn", "Remote")]
+    public class DefaultThresholdShould
+    {
+        static double[] Flows = { 0, 100000 };
+        static double[] Stages = { 0, 15, 30 };
+        static string damCat = "residential";
+        static string assetCat = "Structure";
+        static int id = 1;
+
+        private static CurveMetaData MakeMetaData(string damageCategory = "residential", string assetCategory = "Structure")
+        {
+            return new CurveMetaData("x", "y", "name", damageCategory, assetCategory);
+        }
+
+        private static UncertainPairedData MakeFlowStage()
+        {
+            IDistribution[] stages = new IDistribution[2];
+            for (int i = 0; i < 2; i++)
+            {
+                stages[i] = IDistributionFactory.FactoryUniform(0, 30 * i, 10);
+            }
+            return new UncertainPairedData(Flows, stages, MakeMetaData());
+        }
+
+        private static UncertainPairedData MakeStageDamage(string damageCategory = "residential", string assetCategory = "Structure")
+        {
+            IDistribution[] damages = new IDistribution[3]
+            {
+                new Uniform(0, 0, 10),
+                new Uniform(0, 600000, 10),
+                new Uniform(0, 600000, 10)
+            };
+            return new UncertainPairedData(Stages, damages, MakeMetaData(damageCategory, assetCategory));
+        }
+
+        /// <summary>
+        /// Regression test: adding a second identical damage category should not change the default threshold stage.
+        /// Before the fix, the first damage-frequency curve was double-counted, inflating the threshold.
+        /// With N identical categories the threshold damage scales by N, but the stage-damage also scales by N,
+        /// so the inverse lookup should return the same stage.
+        /// </summary>
+        [Fact]
+        public void ProduceSameThresholdStageRegardlessOfDamageCategoryCount()
+        {
+            ConvergenceCriteria convergenceCriteria = new ConvergenceCriteria(minIterations: 1, maxIterations: 1);
+            ContinuousDistribution flow_frequency = new Uniform(0, 100000, 1000);
+            UncertainPairedData flow_stage = MakeFlowStage();
+
+            // Single damage category
+            List<UncertainPairedData> singleCategory = new List<UncertainPairedData> { MakeStageDamage() };
+            ImpactAreaScenarioSimulation simSingle = ImpactAreaScenarioSimulation.Builder(id)
+                .WithFlowFrequency(flow_frequency)
+                .WithFlowStage(flow_stage)
+                .WithStageDamages(singleCategory)
+                .Build();
+            ImpactAreaScenarioResults resultsSingle = simSingle.Compute(convergenceCriteria, new CancellationToken(), computeIsDeterministic: true);
+
+            // Two damage categories (different category names so they're treated separately)
+            List<UncertainPairedData> twoCategories = new List<UncertainPairedData>
+            {
+                MakeStageDamage("residential", "Structure"),
+                MakeStageDamage("commercial", "Structure")
+            };
+            ImpactAreaScenarioSimulation simTwo = ImpactAreaScenarioSimulation.Builder(id)
+                .WithFlowFrequency(flow_frequency)
+                .WithFlowStage(flow_stage)
+                .WithStageDamages(twoCategories)
+                .Build();
+            ImpactAreaScenarioResults resultsTwo = simTwo.Compute(convergenceCriteria, new CancellationToken(), computeIsDeterministic: true);
+
+            Threshold defaultSingle = resultsSingle.PerformanceByThresholds.ListOfThresholds
+                .First(t => t.ThresholdID == 0);
+            Threshold defaultTwo = resultsTwo.PerformanceByThresholds.ListOfThresholds
+                .First(t => t.ThresholdID == 0);
+
+            // Both should produce the same threshold stage because adding identical categories
+            // scales both numerator (frequency-damage) and denominator (stage-damage) equally.
+            Assert.Equal(defaultSingle.ThresholdValue, defaultTwo.ThresholdValue, precision: 2);
+        }
+
+        /// <summary>
+        /// Verifies that a computed default threshold exists, has the correct type, and has a reasonable stage value
+        /// (between the min and max stages in the stage-damage function).
+        /// </summary>
+        [Fact]
+        public void BeCreatedAutomaticallyWhenNoThresholdProvided()
+        {
+            ConvergenceCriteria convergenceCriteria = new ConvergenceCriteria(minIterations: 1, maxIterations: 1);
+            ContinuousDistribution flow_frequency = new Uniform(0, 100000, 1000);
+            UncertainPairedData flow_stage = MakeFlowStage();
+            List<UncertainPairedData> stageDamageList = new List<UncertainPairedData> { MakeStageDamage() };
+
+            ImpactAreaScenarioSimulation simulation = ImpactAreaScenarioSimulation.Builder(id)
+                .WithFlowFrequency(flow_frequency)
+                .WithFlowStage(flow_stage)
+                .WithStageDamages(stageDamageList)
+                .Build();
+
+            ImpactAreaScenarioResults results = simulation.Compute(convergenceCriteria, new CancellationToken(), computeIsDeterministic: true);
+
+            Threshold defaultThreshold = results.PerformanceByThresholds.ListOfThresholds
+                .FirstOrDefault(t => t.ThresholdID == 0);
+
+            Assert.NotNull(defaultThreshold);
+            Assert.Equal(ThresholdEnum.DefaultExteriorStage, defaultThreshold.ThresholdType);
+            // The threshold stage should fall within the range of the stage-damage function
+            Assert.True(defaultThreshold.ThresholdValue >= Stages[0], "Default threshold stage should be >= minimum stage");
+            Assert.True(defaultThreshold.ThresholdValue <= Stages[^1], "Default threshold stage should be <= maximum stage");
+        }
+
+        /// <summary>
+        /// When a user provides an explicit default threshold (ID=0), the compute should use that
+        /// value and not override it with a calculated one.
+        /// </summary>
+        [Fact]
+        public void NotOverrideUserProvidedDefaultThreshold()
+        {
+            double userStage = 20.0;
+            ConvergenceCriteria convergenceCriteria = new ConvergenceCriteria(minIterations: 1, maxIterations: 1);
+            ContinuousDistribution flow_frequency = new Uniform(0, 100000, 1000);
+            UncertainPairedData flow_stage = MakeFlowStage();
+            List<UncertainPairedData> stageDamageList = new List<UncertainPairedData> { MakeStageDamage() };
+
+            Threshold userThreshold = new Threshold(0, convergenceCriteria, ThresholdEnum.AdditionalExteriorStage, userStage);
+            ImpactAreaScenarioSimulation simulation = ImpactAreaScenarioSimulation.Builder(id)
+                .WithFlowFrequency(flow_frequency)
+                .WithFlowStage(flow_stage)
+                .WithStageDamages(stageDamageList)
+                .WithAdditionalThreshold(userThreshold)
+                .Build();
+
+            ImpactAreaScenarioResults results = simulation.Compute(convergenceCriteria, new CancellationToken(), computeIsDeterministic: true);
+
+            Threshold defaultThreshold = results.PerformanceByThresholds.ListOfThresholds
+                .First(t => t.ThresholdID == 0);
+
+            Assert.Equal(userStage, defaultThreshold.ThresholdValue);
+        }
+    }
+}

HEC.FDA.ModelTest/unittests/GraphicalUncertaintyPairedDataTests.cs

@@ -138,5 +138,28 @@ public void SamplePairedDataShould(double[] inputProbabilities, double[] inputSt
 
         }
 
+        [Theory]
+        [InlineData(new double[] { .99, .5, .1, .02, .01, .002 }, new double[] { 500, 2000, 34900, 66900, 86000, 146000 }, 5, false)]
+        [InlineData(new double[] { .99, .95, .90, .85, .8, .75, .7, .65, .6, .55, .5, .45, .4, .35, .3, .25, .2, .15, .1, .05, .02, .01, .005, .0025 },
+            new double[] { 6.6, 7.4, 8.55, 9.95, 11.5, 12.7, 13.85, 14.7, 15.8, 16.7, 17.5, 18.25, 19, 19.7, 20.3, 21.1, 21.95, 23, 24.2, 25.7, 27.4, 28.4, 29.1, 29.4 },
+            20, true)]
+        [InlineData(new double[] { 0.5, 0.2, 0.1, 0.04, 0.02, 0.01, 0.005, 0.002 }, new double[] { -2, -1.9, 1.93, 1.98, 2.02, 2.04, 2.18, 2.3 }, 40, true)]
+        [InlineData(new double[] { 0.999, 0.5, 0.2, 0.1, 0.04, 0.02, 0.01, 0.004, 0.002 }, new double[] { 80, 82, 84, 84.5, 84.8, 85, 86, 88, 90 }, 50, true)]
+        public void DeterministicSamplingShouldReturnInputRelationship(double[] exceedanceProbabilities, double[] flowOrStageValues, int erl, bool usingStagesNotFlows)
+        {
+            GraphicalUncertainPairedData graphical = new GraphicalUncertainPairedData(exceedanceProbabilities, flowOrStageValues, erl, curveMetaData, usingStagesNotFlows);
+            PairedData sampled = graphical.SamplePairedData(1, computeIsDeterministic: true);
+
+            for (int i = 0; i < exceedanceProbabilities.Length; i++)
+            {
+                double nonExceedanceProbability = 1 - exceedanceProbabilities[i];
+                double actual = sampled.f(nonExceedanceProbability);
+                double expected = flowOrStageValues[i];
+                double tolerance = usingStagesNotFlows ? 0.01 : Math.Max(1, Math.Abs(expected) * 0.01);
+                Assert.True(Math.Abs(actual - expected) < tolerance,
+                    $"At exceedance probability {exceedanceProbabilities[i]}: expected {expected}, got {actual}");
+            }
+        }
+
     }
 }

HEC.FDA.ModelTest/unittests/PairedDataShould.cs

@@ -74,8 +74,9 @@ public void SumYsForGivenX_Test()
             double[] subjectY = new double[5] { 20, 30, 40, 50, 60 };
             PairedData subjectPairedData = new PairedData(subjectX, subjectY);
 
-            double[] expectedX = new double[6] { 1, 2, 3, 4.5, 5, 6 };
-            double[] expectedY = new double[6] { 30, 40, 60, 90, 100, 120 };
+            // Union of x-grids: {1, 2, 3, 4, 4.5, 5, 6}
+            double[] expectedX = new double[7] { 1, 2, 3, 4, 4.5, 5, 6 };
+            double[] expectedY = new double[7] { 30, 40, 60, 80, 90, 100, 120 };
             PairedData expected = new PairedData(expectedX, expectedY);
 
             IPairedData actual = subjectPairedData.SumYsForGivenX(inputPairedData);