Skip to content

Commit

Permalink
Merge pull request #1152 from HydrologicEngineeringCenter/1151-not-en…
Browse files Browse the repository at this point in the history 
…ough-compute-points-in-stage-damage-at-bottom

We need more compute points at the bottom for weird situations
  • Loading branch information
@Brennan1994
Brennan1994 authored Oct 1, 2024
2 parents af9be60 + 4e8fb8c commit ed3319c
Showing 1 changed file with 64 additions and 46 deletions.
110 changes: 64 additions & 46 deletions HEC.FDA.Model/stageDamage/ImpactAreaStageDamage.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,14 @@ public class ImpactAreaStageDamage : PropertyValidationHelper, IDontImplementVal
#region Hard Coded Compute Settings
private const double MIN_PROBABILITY = 0.0001;
private const double MAX_PROBABILITY = 0.9999;
//this setting is used to identify the quantity of compute points
//depth-percent damage functions typically defined at half-foot intervals
//0.25ft target more than preserves the level of information content
private const double FEET_PER_COORDINATE = 0.25;
//require at least four coordinates of stages at which we'll calculate damage outside the events in the hydraulic profiles
private const int MINIMUM_EXTRAPOLATION_COORDINATES = 4;
//require at least two coordinates of stages at which we'll calculate damage betweem the events in the hydraulic profiles
private const int MINIMUM_INTERPOLATION_COORDINATES = 2;
private readonly ConvergenceCriteria _ConvergenceCriteria = new(minIterations: 1000, maxIterations: 5000);
#endregion

Expand All @@ -40,15 +48,16 @@ public class ImpactAreaStageDamage : PropertyValidationHelper, IDontImplementVal
private double _MinStageForArea;
private double _MaxStageForArea;

//we have a way to figure out how many additional stages for which
//to calculate damage (we start with about 8 points)
//these are the number of stages interpolated at the top of the stage damage function
private int _MoreInterpolationPoints;
//these are the number of stages interpolated at the bottom and in the middle of the function
private int _LessInterpolationPoints;

//these are the number of stages extrapolated at the top of the stage-damage function
private int _TopExtrapolationPoints;
//these are the number of stages interpolated between profiles
private int _CentralInterpolationPoints;
//these are the number of stages extrapolated at the bottom of the stage-damage functiom
private int _BottomExtrapolationPoints;

private readonly string _HydraulicParentDirectory;
private readonly PairedData _StageFrequency;
private PairedData _StageFrequency;
private double[] _StagesAtIndexLocation;
#endregion

Expand Down Expand Up @@ -80,37 +89,46 @@ public ImpactAreaStageDamage(int impactAreaID, Inventory inventory, HydraulicDat
Inventory = inventory.GetInventoryTrimmedToImpactArea(impactAreaID);
}
_HydraulicDataset = hydraulicDataset;
SetMinAndMaxStage();
SetCoordinateQuantity();
_StageFrequency = CreateStageFrequency();
EstablishAggregationStages();
}
#endregion

#region Methods
//larger ranges need more points to preserve information content
private void SetCoordinateQuantity()
private void EstablishAggregationStages()
{
//depth-percent damage functions typically defined at half-foot intervals
//this preserves the level of information content
double feetPerCoordinate = 0.5;
double range = _MaxStageForArea - _MinStageForArea;
int setsOfCoordinatesBetweenProfiles = 10;
int coordinateQuantity = Convert.ToInt32(Math.Ceiling((range / feetPerCoordinate) / setsOfCoordinatesBetweenProfiles));

//require at least two coordinates to interpolate and extrapolate
if (coordinateQuantity < 4)
{
coordinateQuantity = 4;
}
_MoreInterpolationPoints = coordinateQuantity * 4;
_LessInterpolationPoints = coordinateQuantity;
_StageFrequency = IdentifyCentralStageFrequencyAtIndexLocation();
IdentifyMinAndMaxStageWithUncertainty();
SetCoordinateQuantity();
}

private void SetCoordinateQuantity()
{
//set bottom coordinate quantity
double stageAtAEPofMostFrequentHydraulicsProfile = _StageFrequency.f(_HydraulicDataset.HydraulicProfiles.First().Probability);
double rangeOfStagesAtBottom = stageAtAEPofMostFrequentHydraulicsProfile - _MinStageForArea;
_BottomExtrapolationPoints = Convert.ToInt32(Math.Ceiling(rangeOfStagesAtBottom/FEET_PER_COORDINATE));
//require at least 4 coordinates at the bottom
if (_BottomExtrapolationPoints < MINIMUM_EXTRAPOLATION_COORDINATES) { _BottomExtrapolationPoints = MINIMUM_EXTRAPOLATION_COORDINATES; }

//set middle coordinate quantity
double stageAtAEPofLeastFrequentHydraulicsProfile = _StageFrequency.f(_HydraulicDataset.HydraulicProfiles.Last().Probability);
double middleRange = stageAtAEPofLeastFrequentHydraulicsProfile - stageAtAEPofMostFrequentHydraulicsProfile;
//based upon the range of stages between most frequent and least frequent AEP in hydraulics and the number of intervals to be interpolated
_CentralInterpolationPoints = Convert.ToInt32(Math.Ceiling((middleRange / FEET_PER_COORDINATE) / (_HydraulicDataset.HydraulicProfiles.Count() - 1)));
//require at least two coordinates between profiles
if (_CentralInterpolationPoints < MINIMUM_INTERPOLATION_COORDINATES) { _CentralInterpolationPoints = MINIMUM_INTERPOLATION_COORDINATES; }

//set top coordinate quantity
double rangeOfStagesAtTop = _MaxStageForArea - stageAtAEPofLeastFrequentHydraulicsProfile;
_TopExtrapolationPoints = Convert.ToInt32(Math.Ceiling(rangeOfStagesAtTop/FEET_PER_COORDINATE));
//require at least four coordinates at the top
if ( _TopExtrapolationPoints < MINIMUM_EXTRAPOLATION_COORDINATES) { _TopExtrapolationPoints = MINIMUM_EXTRAPOLATION_COORDINATES; }
}

/// <summary>
/// This method is used to identify the minimum stage at the index location and the maximum stage at the index location for which we will calculate damage
/// </summary>
private void SetMinAndMaxStage()
private void IdentifyMinAndMaxStageWithUncertainty()
{
if (_AnalyticalFlowFrequency != null)
{
Expand Down Expand Up @@ -188,7 +206,7 @@ private void SetMinAndMaxStage()
/// This method grabs the input summary relationships and generates the median stage frequency function
/// The frequencies in the function are used to align the aggregation stages to the stages at the structures
/// </summary>
private PairedData CreateStageFrequency()
private PairedData IdentifyCentralStageFrequencyAtIndexLocation()
{
if (_AnalyticalFlowFrequency != null)
{
Expand Down Expand Up @@ -414,14 +432,14 @@ private List<StudyAreaConsequencesBinned> CreateConsequenceDistributionResults(s
private double[] ComputeStagesAtIndexLocation(List<double> profileProbabilities)
{
//less stages at the bottom, more stages at the top, less stages in between
int quantityStages = _LessInterpolationPoints + _MoreInterpolationPoints + (profileProbabilities.Count - 1) * _LessInterpolationPoints;
int quantityStages = _BottomExtrapolationPoints + _TopExtrapolationPoints + (profileProbabilities.Count - 1) * _CentralInterpolationPoints;
double[] stages = new double[quantityStages];
//extrapolate lower stages
double stageAtProbabilityOfLowestProfile = _StageFrequency.f(1 - profileProbabilities.Max());
float indexStationLowerStageDelta = (float)(stageAtProbabilityOfLowestProfile - _MinStageForArea);
float interval = indexStationLowerStageDelta / _LessInterpolationPoints;
float interval = indexStationLowerStageDelta / _BottomExtrapolationPoints;
int stageIndex = 0;
for (int i = 0; i < _LessInterpolationPoints + 1; i++)
for (int i = 0; i < _BottomExtrapolationPoints + 1; i++)
{
stages[i] = (_MinStageForArea + i * interval);
stageIndex++;
Expand All @@ -434,12 +452,12 @@ private double[] ComputeStagesAtIndexLocation(List<double> profileProbabilities)
double previousProbability = profileProbabilities[i - 1];
double currentProbability = profileProbabilities[i];

for (int j = 0; j < _LessInterpolationPoints; j++)
for (int j = 0; j < _CentralInterpolationPoints; j++)
{
double previousStageAtIndexLocation = _StageFrequency.f(1 - previousProbability);
double currentStageAtIndexLocation = _StageFrequency.f(1 - currentProbability);
double stageDeltaAtIndexLocation = currentStageAtIndexLocation - previousStageAtIndexLocation;
double intervalAtIndexLocation = stageDeltaAtIndexLocation / _LessInterpolationPoints;
double intervalAtIndexLocation = stageDeltaAtIndexLocation / _CentralInterpolationPoints;
double stageAtIndexLocation = previousStageAtIndexLocation + intervalAtIndexLocation * (j + 1);
stages[stageIndex] = stageAtIndexLocation;
stageIndex++;
Expand All @@ -449,11 +467,11 @@ private double[] ComputeStagesAtIndexLocation(List<double> profileProbabilities)
//extrapolate upper stages
double stageAtProbabilityOfHighestProfile = _StageFrequency.f(1 - profileProbabilities.Min());
float indexStationUpperStageDelta = (float)(_MaxStageForArea - stageAtProbabilityOfHighestProfile);
float upperInterval = indexStationUpperStageDelta / _MoreInterpolationPoints;
for (int i = 1; i < _MoreInterpolationPoints; i++)
float upperInterval = indexStationUpperStageDelta / _TopExtrapolationPoints;
for (int i = 1; i < _TopExtrapolationPoints; i++)
{

stages[stageIndex] = (_MaxStageForArea - upperInterval * (_MoreInterpolationPoints - i));
stages[stageIndex] = (_MaxStageForArea - upperInterval * (_TopExtrapolationPoints - i));
stageIndex++;
}

Expand Down Expand Up @@ -484,10 +502,10 @@ private void ComputeLowerStageDamage(ref List<StudyAreaConsequencesBinned> paral

float interval = CalculateLowerIncrementOfStages(profileProbabilities);
List<float[]> stagesAtAllStructuresAllEvents = new();
for (int stageIndex = 0; stageIndex < _LessInterpolationPoints + 1; stageIndex++)
for (int stageIndex = 0; stageIndex < _BottomExtrapolationPoints + 1; stageIndex++)
{
//for each stage, add the consequenceResult to the consequenceResultArray in the correct place
float[] WSEsParallelToIndexLocation = ExtrapolateFromBelowStagesAtIndexLocation(inventoryAndWaterCoupled.Item2[0], interval, stageIndex, _LessInterpolationPoints);
float[] WSEsParallelToIndexLocation = ExtrapolateFromBelowStagesAtIndexLocation(inventoryAndWaterCoupled.Item2[0], interval, stageIndex, _BottomExtrapolationPoints);
//Can we modify the below to push more of the calculation into the parallelization. So, instead of passing in a float[] wses, it might be a float[][].
stagesAtAllStructuresAllEvents.Add(WSEsParallelToIndexLocation);
}
Expand All @@ -507,7 +525,7 @@ private float CalculateLowerIncrementOfStages(List<double> profileProbabilities)
//the delta is the difference between the min stage at the index location and the stage at the index location for the lowest profile
float indexStationLowerStageDelta = (float)(stageAtProbabilityOfLowestProfile - _MinStageForArea);
//this interval defines the interval in stages by which we'll compute damage
float interval = indexStationLowerStageDelta / _LessInterpolationPoints;
float interval = indexStationLowerStageDelta / _BottomExtrapolationPoints;
//Collect damage for first part of function up to and including the stages at the lowest profile
return interval;
}
Expand All @@ -528,11 +546,11 @@ public static float[] ExtrapolateFromBelowStagesAtIndexLocation(float[] WSEsAtLo
private void ComputeMiddleStageDamage(ref List<StudyAreaConsequencesBinned> parallelConsequenceResultCollection, string damageCategory, List<DeterministicOccupancyType> deterministicOccTypes, (Inventory, List<float[]>) inventoryAndWaterCoupled, List<double> profileProbabilities, int iterationIndex)
{
int numProfiles = profileProbabilities.Count;
int stageIndex = _LessInterpolationPoints + 1;
int stageIndex = _BottomExtrapolationPoints + 1;
for (int profileIndex = 1; profileIndex < numProfiles; profileIndex++)
{
InterpolateBetweenProfiles(ref parallelConsequenceResultCollection, deterministicOccTypes, inventoryAndWaterCoupled.Item2[profileIndex - 1], inventoryAndWaterCoupled.Item2[profileIndex], damageCategory, inventoryAndWaterCoupled.Item1, stageIndex, iterationIndex);
stageIndex += _LessInterpolationPoints;
stageIndex += _CentralInterpolationPoints;
}
}

Expand All @@ -541,7 +559,7 @@ private void InterpolateBetweenProfiles(ref List<StudyAreaConsequencesBinned> pa
{
float[] intervalsAtStructures = CalculateIntervals(previousHydraulicProfile, currentHydraulicProfile);
List<float[]> stagesAllStructuresAllStages = new();
for (int interpolatorIndex = 0; interpolatorIndex < _LessInterpolationPoints; interpolatorIndex++)
for (int interpolatorIndex = 0; interpolatorIndex < _CentralInterpolationPoints; interpolatorIndex++)
{
float[] stages = CalculateIncrementOfStages(previousHydraulicProfile, intervalsAtStructures, interpolatorIndex + 1);
stagesAllStructuresAllStages.Add(stages);
Expand Down Expand Up @@ -570,7 +588,7 @@ private float[] CalculateIntervals(float[] previousStagesAtStructures, float[] c
float[] intervals = new float[previousStagesAtStructures.Length];
for (int j = 0; j < previousStagesAtStructures.Length; j++)
{
intervals[j] = (currentStagesAtStructures[j] - previousStagesAtStructures[j]) / _LessInterpolationPoints;
intervals[j] = (currentStagesAtStructures[j] - previousStagesAtStructures[j]) / _CentralInterpolationPoints;
}
return intervals;
}
Expand All @@ -580,13 +598,13 @@ private float[] CalculateIntervals(float[] previousStagesAtStructures, float[] c
private void ComputeUpperStageDamage(ref List<StudyAreaConsequencesBinned> parallelConsequenceResultCollection, string damageCategory, List<DeterministicOccupancyType> deterministicOccTypes, (Inventory, List<float[]>) inventoryAndWaterCoupled, List<double> profileProbabilities, int iterationIndex)
{
//the probability of a profile is an EXCEEDANCE probability but in the model we use NONEXCEEDANCE PROBABILITY
int stageIndex = _LessInterpolationPoints + _LessInterpolationPoints * (profileProbabilities.Count - 1);
int stageIndex = _BottomExtrapolationPoints + _CentralInterpolationPoints * (profileProbabilities.Count - 1);
double stageAtProbabilityOfHighestProfile = _StageFrequency.f(1 - profileProbabilities.Min());
float indexStationUpperStageDelta = (float)(_MaxStageForArea - stageAtProbabilityOfHighestProfile);
float upperInterval = indexStationUpperStageDelta / _MoreInterpolationPoints;
float upperInterval = indexStationUpperStageDelta / _TopExtrapolationPoints;

List<float[]> stagesAllStructuresAllEvents = new();
for (int extrapolatorIndex = 1; extrapolatorIndex < _MoreInterpolationPoints; extrapolatorIndex++)
for (int extrapolatorIndex = 1; extrapolatorIndex < _TopExtrapolationPoints; extrapolatorIndex++)
{
float[] WSEsParallelToIndexLocation = ExtrapolateFromAboveAtIndexLocation(inventoryAndWaterCoupled.Item2[^1], upperInterval, extrapolatorIndex);
stagesAllStructuresAllEvents.Add(WSEsParallelToIndexLocation);
Expand Down

0 comments on commit ed3319c

Please sign in to comment.