v501 - ECCC Assimilation; Non-linear LP Solver; LP Diagnostics [BENCHMARKED]

Implementing ECCC Assimilation strategy (Assimilate.cpp)
-CReservoir::AdjustStage() (Reservoir.cpp)
-Options.assim_method (ParseInput.cpp;RavenInclude.h)
-:AssimilationMethod command

Non-linear management optimization.cpp
-@pow() function in management optimization (DemandExpressionHandling.cpp)
-new member data, _pNonLinVars, _maxIterations, _iterTolerance, _relaxCoeff, assoc parsing commands (DemandOptimization.cpp;ParseManagementFile.cpp)
-new TERM_ITER, TERM_POW (Expression.h/DemandExpressionHandling.cpp)
-adapt AddConstraintToLP
-properly handle blanks in nested TS statements (DemandExpressionHandling.cpp)

Better LP_SOLVE infeasible soln handling/diagnostics (DemandOptimization.cpp)

retain UPSTREAM_OF_INCLUSIVE for backward compat (ParseHRUFile.cpp)

improved warnings in ParseInitialConditionFile.cpp

Minor changes triggered by strict debugging process:
-GWRiverConnection.cpp/.h
-EnergyTransport.cpp
-SubBasin.cpp

Renamed count to _count in CustomOutput.cpp/.h
Renamed SolveDemandProblem to SolveManagementProblem() (DemandOptimization.cpp;Solvers.cpp)
Minor renaming in IsotopeTransport.cpp

Author: James Craig

src/Assimilate.cpp

@@ -77,7 +77,7 @@ void CModel::InitializeDataAssimilation(const optStruct &Options)
 
 }
 /////////////////////////////////////////////////////////////////
-/// \brief Overrides flows with observations at gauges, or
+/// \brief Overrides flows with observations at gauges and propagates adjustments upstream of gauges 
 /// \param p [in] global subbasin index
 /// \param Options [in] current model options structure
 /// \param tt [in] current model time structure
@@ -111,8 +111,15 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim
 
   // actually scale flows
   //---------------------------------------------------------------------
-  double mass_added=_pSubBasins[p]->ScaleAllFlows(_aDAscale[p]/_aDAscale_last[p],_aDAoverride[p],Options.timestep,tt.model_time);
-  //double mass_added=_pSubBasins[p]->AdjustAllFlows(_aDAadjust[p],_aDAoverride[p],Options.timestep,tt.model_time);
+  double mass_added;
+  if (Options.assim_method==DA_RAVEN_DEFAULT){
+    mass_added=_pSubBasins[p]->ScaleAllFlows(_aDAscale[p]/_aDAscale_last[p],_aDAoverride[p],Options.timestep,tt.model_time);
+  }
+  else if (Options.assim_method==DA_ECCC) {
+    mass_added=_pSubBasins[p]->AdjustAllFlows(_aDAQadjust[p],_aDAoverride[p],Options.timestep,tt.model_time);
+  }
+
+  // 
   if(mass_added>0.0){_CumulInput +=mass_added/(_WatershedArea*M2_PER_KM2)*MM_PER_METER;}
   else              {_CumulOutput-=mass_added/(_WatershedArea*M2_PER_KM2)*MM_PER_METER;}
 
@@ -170,7 +177,7 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
             _aDAoverride [p]=true;
             _aDAobsQ     [p]=Qobs;
             _aDADrainSum [p]=0.0; //??? maybe doesnt matter
-            if (pdown != DOESNT_EXIST) {
+            if (pdown != DOESNT_EXIST) {      
               _aDADrainSum [pdown]+=_pSubBasins[p]->GetDrainageArea(); //DOES THIS HANDLE NESTING RIGHT?
             }
           }
@@ -182,17 +189,17 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
             _aDAlength   [p]=0.0;
             _aDAoverride [p]=false;
             _aDAobsQ     [p]=0.0;
-            if (pdown != DOESNT_EXIST) {
+            if (pdown != DOESNT_EXIST) {        
               _aDADrainSum[pdown] += _aDADrainSum[p];
             }
           }
           ObsExists=true;
           break; //avoids duplicate observations
         }
       }
-    }
+    } 
     else {
-      if (pdown != DOESNT_EXIST) {
+      if (pdown != DOESNT_EXIST) {        
         _aDADrainSum[pdown] += _aDADrainSum[p];
       }
     }
@@ -247,6 +254,8 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
     else {
       ECCCwt=1.0;
     }
-    _aDAQadjust[p] = _aDAQadjust[p]*ECCCwt;
+    if (!_aDAoverride[p]) { //no scaling for stations being overridden, only upstream
+      _aDAQadjust[p] = _aDAQadjust[p]*ECCCwt;
+    }
   }
 }

src/CustomOutput.cpp

@@ -69,7 +69,7 @@ CCustomOutput::CCustomOutput( const diagnostic    variable,
   _hist_max     =10;
   _nBins        =1; //Arbitrary default
 
-  count         =0;
+  _count        =0;
 
   kk_only       =kk;
 
@@ -870,19 +870,14 @@ void CCustomOutput::WriteCustomOutput(const time_struct &tt,
       else if (_spaceAgg==BY_SELECT_HRUS){val=pModel->GetHRUGroup (kk_only)->GetHRU(k)->GetCumulFluxBet(_svind,_svind2);}
     }
 
-    if (k==0){count++;}//increment number of data items stored
+    if (k==0){_count++;}//increment number of data items stored
 
     //---Update diagnostics--------------------------------------------------
     //-----------------------------------------------------------------------
     if      (_aggstat==AGG_AVERAGE)
     {
-      // \todo[funct] - should handle pointwise variables (e.g., state vars) differently from periodwise variables (e.g., forcings)
-      if      (_var == VAR_STATE_VAR){
-        data[k][0]=(double)(count-1)/(double)(count)*data[k][0]+val/count; //NOT CURRENTLY STRICTLY VALID
-      }
-      else {
-        data[k][0]=(double)(count-1)/(double)(count)*data[k][0]+val/count;
-      }
+      //incrementally adding to average - 1/N*((N-1)(old mean)+(val))
+      data[k][0]=(double)(_count-1)/(double)(_count)*data[k][0]+val/(double)(_count);
     }
     else if (_aggstat==AGG_MAXIMUM)
     {
@@ -907,7 +902,7 @@ void CCustomOutput::WriteCustomOutput(const time_struct &tt,
              (_aggstat==AGG_HISTOGRAM))
     {
       //populate data
-      data [k][count-1] = val;
+      data [k][_count-1] = val;
     }
     else
     {
@@ -933,31 +928,31 @@ void CCustomOutput::WriteCustomOutput(const time_struct &tt,
         else if(_aggstat==AGG_MEDIAN)
         {
           double Q1,Q2,Q3;
-          quickSort(data[k],0,count-1);
-          GetQuartiles(data[k],count,Q1,Q2,Q3);
+          quickSort(data[k],0,_count-1);
+          GetQuartiles(data[k],_count,Q1,Q2,Q3);
           _CUSTOM<<FormatDouble(Q2)<<sep;
         }
         else if(_aggstat==AGG_QUARTILES) //find lower quartile, median, then upper quartile
         {
           double Q1,Q2,Q3;
-          quickSort(data[k],0,count-1);
-          GetQuartiles(data[k],count,Q1,Q2,Q3);
+          quickSort(data[k],0,_count-1);
+          GetQuartiles(data[k],_count,Q1,Q2,Q3);
           _CUSTOM<<FormatDouble(Q1)<<sep<<FormatDouble(Q2)<<sep<<FormatDouble(Q3)<<sep;
         }
         else if(_aggstat==AGG_95CI)
         {
-          quickSort(data[k],0,count-1);//take floor and ceiling of lower and upper intervals to be conservative
-          _CUSTOM  << FormatDouble(data[k][(int)floor((double)(count-1)*0.025)])<<sep;
-          _CUSTOM  << FormatDouble(data[k][(int)ceil((double)(count-1)*0.975)])<<sep;
+          quickSort(data[k],0,_count-1);//take floor and ceiling of lower and upper intervals to be conservative
+          _CUSTOM  << FormatDouble(data[k][(int)floor((double)(_count-1)*0.025)])<<sep;
+          _CUSTOM  << FormatDouble(data[k][(int)ceil((double)(_count-1)*0.975)])<<sep;
         }
         else if(_aggstat==AGG_HISTOGRAM)
         {
-          quickSort(data[k],0,count-1);
+          quickSort(data[k],0,_count-1);
           double binsize = (_hist_max-_hist_min)/_nBins;
           int bincount[MAX_HISTOGRAM_BINS];
 
           for(int bin=0;bin<_nBins;bin++){ bincount[bin]=0; }
-          for(int a=0;a<count;a++)
+          for(int a=0;a<_count;a++)
           {
             for(int bin=0;bin<_nBins;bin++){
               if((data[k][a]>=(_hist_min+bin*binsize)) && (data[k][a]<(_hist_min+(bin+1)*binsize))){ bincount[bin]++; }
@@ -982,8 +977,8 @@ void CCustomOutput::WriteCustomOutput(const time_struct &tt,
         else if(_aggstat==AGG_MEDIAN)
         {
           double Q1,Q2,Q3;
-          quickSort(data[k],0,count-1);
-          GetQuartiles(data[k],count,Q1,Q2,Q3);
+          quickSort(data[k],0,_count-1);
+          GetQuartiles(data[k],_count,Q1,Q2,Q3);
           out=Q2;
         }
         output[k]=out;
@@ -1011,7 +1006,7 @@ void CCustomOutput::WriteCustomOutput(const time_struct &tt,
         }
       }
       //...more stats here
-      if (k==num_data-1){count=0;}//don't reboot count until last HRU/basin is done
+      if (k==num_data-1){_count=0;}//don't reboot count until last HRU/basin is done
     }//end if (reset)
   }//end for (k=0;...
 
@@ -1218,4 +1213,4 @@ CCustomOutput *CCustomOutput::ParseCustomOutputCommand(char *s[MAXINPUTITEMS], c
   }
 
   return pCustom;
-}
+}

src/CustomOutput.h

@@ -1,6 +1,6 @@
 /*----------------------------------------------------------------
   Raven Library Source Code
-  Copyright (c) 2008-2023 the Raven Development Team
+  Copyright (c) 2008-2025 the Raven Development Team
   ----------------------------------------------------------------
   Custom output generation
   ----------------------------------------------------------------*/
@@ -96,7 +96,7 @@ class CCustomOutput
                             //(e.g., =2 if max and min are both tracked)
   int         _time_index;  ///< index tracking current output line (e.g., 3=3 years/months/days passed, dependent upon _timeAgg
 
-  int          count;       ///< counts accumulated data (# of timesteps since last output dump)
+  int          _count;      ///< counts accumulated data (# of timesteps since last output dump)
 
   int          kk_only;     ///< index of HRUGroup for which output is generated when spaceAgg==BY_SELECT_HRUS