v505: fix assimilation time lag; @dLookUp()

repair to Assimilate to remove time lag (Assimilate.cpp)
new @dLookUp() function (DemandExpressionHandling.cpp;Expression.h)

clean up:
deleted duplicative TermTypeToSstring() function (DemandExpressionHandling.cpp)
use :SeasonalCanopyHeight in Raven .rvp template (ParsePropertyFile.cpp)

Author: James Craig

src/Assimilate.cpp

@@ -90,17 +90,16 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim
   //---------------------------------------------------------------------
   if(_aDAoverride[p])
   {
-    double Qobs,Qmod;
+    double Qobs,Qmod,Qmodlast;
     double alpha = _pGlobalParams->GetParams()->assimilation_fact;
 
     Qobs = _aDAobsQ[p];
-    //Option A: mean flow
-    //Qmod = _pSubBasins[p]->GetIntegratedOutflow(Options.timestep)/(Options.timestep*SEC_PER_DAY);
-    //Option B: instantaneous flow
     Qmod = _pSubBasins[p]->GetOutflowRate();
+    Qmodlast= _pSubBasins[p]->GetLastOutflowRate();
     if(Qmod>PRETTY_SMALL) {
       _aDAscale  [p]=1.0+alpha*((Qobs-Qmod)/Qmod); //if alpha = 1, Q=Qobs in observation basin
-      _aDAQadjust[p]=alpha*(Qobs-Qmod);
+      //_aDAQadjust[p]=alpha*(Qobs-Qmod);//Option A: instantaneous flow
+      _aDAQadjust[p]=alpha*(2.0*Qobs-Qmodlast-Qmod);//Option B: mean flow
     }
     else {
       _aDAscale  [p]=1.0;
@@ -166,7 +165,7 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
       {
         if(IsContinuousFlowObs2(_pObservedTS[i],_pSubBasins[p]->GetID()))//flow observation is available and linked to this subbasin
         {
-          Qobs = _pObservedTS[i]->GetSampledValue(nn); //end of timestep flow
+          Qobs = _pObservedTS[i]->GetSampledValue(nn); //mean timestep flow
 
           //bool fakeblank=((tt.model_time>30) && (tt.model_time<40)) || ((tt.model_time>45) && (tt.model_time<47));//TMP DEBUG
           //if (fakeblank){Qobs=RAV_BLANK_DATA;}

src/DemandExpressionHandling.cpp

@@ -279,25 +279,6 @@ int CDemandOptimizer::GetNumUserDVs() const
 //////////////////////////////////////////////////////////////////
 /// \brief enum to string routines
 //
-string TermTypeToString(termtype t)
-{
-  if (t==TERM_DV      ){return "TERM_DV"; }
-  if (t==TERM_TS      ){return "TERM_TS"; }
-  if (t==TERM_LT      ){return "TERM_LT"; }
-  if (t==TERM_HRU     ){return "TERM_HRU";}
-  if (t==TERM_SB      ){return "TERM_SB"; }
-  if (t==TERM_CONST   ){return "TERM_CONST"; }
-  if (t==TERM_WORKFLOW){return "TERM_WORKFLOW"; }
-  if (t==TERM_HISTORY ){return "TERM_HISTORY"; }
-  if (t==TERM_ITER    ){return "TERM_ITER"; }   
-  if (t==TERM_MAX     ){return "TERM_MAX"; }
-  if (t==TERM_MIN     ){return "TERM_MIN"; }
-  if (t==TERM_POW     ){return "TERM_POW"; }
-  if (t==TERM_CUMUL   ){return "TERM_CUMUL";}
-  if (t==TERM_CUMUL_TS){return "TERM_CUMUL_TS";}
-  if (t==TERM_UNKNOWN ){return "TERM_UNKNOWN"; }
-  return "TERM_UNKNOWN";
-}
 string DVTypeToString(dv_type t)
 {
   if (t==DV_QOUT    ){return "DV_QOUT";     }
@@ -321,6 +302,7 @@ string expTypeToString(termtype &typ){
     case(TERM_DV):      return "TERM_DV"; break;
     case(TERM_TS):      return "TERM_TS"; break;
     case(TERM_LT):      return "TERM_LT"; break;
+    case(TERM_DLT):     return "TERM_DLT"; break;
     case(TERM_CONST):   return "TERM_CONST"; break;
     case(TERM_HRU):     return "TERM_HRU"; break;
     case(TERM_SB):      return "TERM_SB"; break;
@@ -333,7 +315,7 @@ string expTypeToString(termtype &typ){
     case(TERM_CUMUL):   return "TERM_CUMUL"; break;
     case(TERM_UNKNOWN): return "TERM_UNKNOWN"; break;
   }
-  return "?";
+  return "TERM_UNKNOWN";
 }
 //////////////////////////////////////////////////////////////////
 /// \brief parses individual expression string and converts to expression term structure
@@ -644,6 +626,43 @@ bool CDemandOptimizer::ConvertToExpressionTerm(const string s, expressionTerm* t
       }
     }
   }
+  //----------------------------------------------------------------------
+  else if (s.substr(0, 9) == "@dlookup(")//derivative of lookup table (e.g., @dlookup(my_table,EXPRESSION))
+  {
+    string name;
+    string x_in;
+    size_t is = s.find("@dlookup(");
+    size_t ie = s.find(",",is);
+    size_t ip = s.find_last_of(")");
+    if (ie == NPOS) {
+      warn="ConvertToExpressionTerm: missing comma in @dlookup expression"+warnstring;
+      ExitGracefully(warn.c_str(),BAD_DATA_WARN);
+    }
+    if (ip == NPOS) {
+      warn="ConvertToExpressionTerm: missing end parentheses in @dlookup expression"+warnstring;
+      ExitGracefully(warn.c_str(), BAD_DATA_WARN);
+    }
+    if ((is != NPOS) && (ie != NPOS))
+    {
+      bool found=false;
+      name = s.substr(is+9,ie-(is+9));
+      x_in = s.substr(ie+1,ip-(ie+1));
+      term->pLT=NULL;
+      for (int i = 0; i < _nUserLookupTables; i++) {
+        if (StringToUppercase(_pUserLookupTables[i]->GetName()) == StringToUppercase(name)) {
+          term->pLT = _pUserLookupTables[i];
+          found=true;
+        }
+      }
+      term->nested_exp1 =x_in;
+      term->type        =TERM_DLT;
+      if (!found) {
+        warn="ConvertToExpression: unrecognized lookup table name in @dlookup command"+warnstring;
+        ExitGracefully(warn.c_str(), BAD_DATA_WARN);
+        return false;
+      }
+    }
+  }
  //----------------------------------------------------------------------
  else if (s.substr(0, 9) == "@HRU_var(") // HRU state var (e.g., @HRU_var(SNOW,[id])
   {
@@ -922,7 +941,7 @@ expressionStruct *CDemandOptimizer::ParseExpression(const char **s,
       if (!valid){return NULL; }
 
       //TMP DEBUG
-      //cout<<"   TERM["<<k<<" in "<< j << "]: " << s[i] << " : type = " << TermTypeToString(terms[j][k]->type) << " index : " << terms[j][k]->DV_ind << endl;
+      //cout<<"   TERM["<<k<<" in "<< j << "]: " << s[i] << " : type = " << expTypeToString(terms[j][k]->type) << " index : " << terms[j][k]->DV_ind << endl;
 
       if ((i>rhs_ind)      && (k==0)){terms[j][k]->mult*=-1.0; } //only multiply leading term by -1 to move to LHS
       if ((s[i-1][0]=='-') && (k==0)){terms[j][k]->mult*=-1.0; } //and reverse if preceded with minus sign
@@ -1492,6 +1511,12 @@ double CDemandOptimizer::EvaluateTerm(expressionTerm **pTerms,const int k, const
     x=EvaluateTerm(pTerms,kk,t);
     return pT->pLT->GetValue(x); //decision variable can't be in lookup table
   }
+  else if (pT->type == TERM_DLT)
+  {
+    kk=pT->nested_ind1;
+    x=EvaluateTerm(pTerms,kk,t);
+    return pT->pLT->GetSlope(x); //decision variable can't be in lookup table
+  }
   else if (pT->type == TERM_HRU)
   {
     int i=pT->SV_index;
@@ -1533,7 +1558,7 @@ double CDemandOptimizer::EvaluateTerm(expressionTerm **pTerms,const int k, const
     else if (tmp=='h'){return _ahhist[_aSBIndices[p]][nshift]; }
     else if (tmp=='D'){return _aDhist[_aSBIndices[p]][nshift]; }
     else if (tmp=='I'){return _aIhist[_aSBIndices[p]][nshift]; }
-    else if (tmp=='B'){return _pModel->GetSubBasin(p)->GetSpecifiedInflow(t+pT->timeshift*1.0); } //ASSUMES DAILY TIMESTEP!
+    else if (tmp=='B'){return _pModel->GetSubBasin(p)->GetSpecifiedInflow(t+pT->timeshift*1.0); } //ASSUMES DAILY TIMESTEP! \todo[funct] should be Options.timestep
     else if (tmp=='E'){return _pModel->GetSubBasin(p)->GetEnviroMinFlow  (t+pT->timeshift*1.0); } //ASSUMES DAILY TIMESTEP!
     else {
       ExitGracefully("CDemandOptimizer::EvaluateTerm: Invalid history variable ",BAD_DATA);

src/DemandOptimization.cpp

@@ -11,7 +11,7 @@
 
 void SummarizeExpression(const char **s, const int Len, expressionStruct* exp); //defined in DemandExpressionHandling.cpp
 string DVTypeToString(dv_type t);
-string TermTypeToString(termtype t);
+string expTypeToString(termtype &t);
 
 //////////////////////////////////////////////////////////////////
 /// \brief Implementation of the Demand optimization constructor
@@ -880,7 +880,7 @@ void CDemandOptimizer::InitializePostRVMRead(CModel* pModel, const optStruct& Op
     expressionTerm* term=new expressionTerm();
     if (ConvertToExpressionTerm(_pNonLinVars[i]->target, term, 0, "internal")) { 
       _pNonLinVars[i]->DV_index=term->DV_ind-1; //using 0 indexing for dv rateher than lp 1 indexing
-      //cout<<" NONLINEAR TARGET FOUND : "<<term->DV_ind-1<<" "<<term->origexp<<" "<<TermTypeToString(term->type)<<endl;
+      //cout<<" NONLINEAR TARGET FOUND : "<<term->DV_ind-1<<" "<<term->origexp<<" "<<expTypeToString(term->type)<<endl;
     }
     else {
       string warn="ManagementOptimization:Initialize cannot convert target '"+_pNonLinVars[i]->target+"' of non-linear variable to decision variable";
@@ -2092,7 +2092,8 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
         if (fabs(_aSolverResiduals[j])>REAL_SMALL)
         cout<<"   -"<<_aSolverRowNames[j] << " " << _aSolverResiduals[j] << endl;
       }
-
+      
+      //lp_lib::print_debugdump(pLinProg, dumpfile.c_str()); 
       WriteLPSubMatrix(pLinProg,"overconstrained_lp_matrix.csv",Options);
       ExitGracefully("SolveManagementProblem: non-optimal solution found. Problem is over-constrained. Remove or adjust management constraints.",RUNTIME_ERR);
     }

src/Expression.h

@@ -1,3 +1,11 @@
+/*----------------------------------------------------------------
+  Raven Library Source Code
+  Copyright (c) 2008-2025 the Raven Development Team
+  ----------------------------------------------------------------
+  enum termtype
+  expressionTerm struct
+  expressionStruct structure
+  ----------------------------------------------------------------*/
 #include "RavenInclude.h"
 #include "TimeSeries.h"
 #include "LookupTable.h"
@@ -11,6 +19,7 @@ enum termtype
   TERM_DV,        //< decision variable !Q123 or named
   TERM_TS,        //< time series @ts(name,n)
   TERM_LT,        //< lookup table @lookup(x)
+  TERM_DLT,       //< lookup table derivative @dlookup(x)
   TERM_HRU,       //< state variable @HRU_var(SNOW,2345)
   TERM_SB,        //< state variable @SB_var(SNOW,234)
   TERM_CONST,     //< constant

src/ParsePropertyFile.cpp

@@ -2059,18 +2059,18 @@ void  CreateRVPTemplate(string *aP,class_type *aPC,int &nP,const optStruct &Opti
     for(int ii=0;ii<nP;ii++)
     {
       if((aPC[ii]==CLASS_VEGETATION) && (aP[ii]!="RELATIVE_LAI")){
-        TEMPLATE<<":SeasonalRelativeLAI"<<endl;
+        TEMPLATE<<":SeasonalCanopyLAI"<<endl;
         TEMPLATE<<"  *VEGET_1*, *J*,*F*,*M*,*A*,*M*,*J*,*J*,*A*,*S*,*O*,*N*,*D*"<<endl;
-        TEMPLATE<<":EndSeasonalRelativeLAI"<<endl;
+        TEMPLATE<<":EndSeasonalCanopyLAI"<<endl;
         break;
       }
     }
     for(int ii=0;ii<nP;ii++)
     {
       if((aPC[ii]==CLASS_VEGETATION) && (aP[ii]!="RELATIVE_HT")){
-        TEMPLATE<<":SeasonalRelativeHeight"<<endl;
+        TEMPLATE<<":SeasonalCanopyHeight"<<endl;
         TEMPLATE<<"  *VEGET_1*, *J*,*F*,*M*,*A*,*M*,*J*,*J*,*A*,*S*,*O*,*N*,*D*"<<endl;
-        TEMPLATE<<":EndSeasonalRelativeHeight"<<endl;
+        TEMPLATE<<":EndSeasonalCanopyHeight"<<endl;
         break;
       }
     }