diff --git a/src/DemandExpressionHandling.cpp b/src/DemandExpressionHandling.cpp index d0848a89..abfa0665 100644 --- a/src/DemandExpressionHandling.cpp +++ b/src/DemandExpressionHandling.cpp @@ -10,7 +10,7 @@ const int MAX_TERMS_PER_GROUP=10; ////////////////////////////////////////////////////////////////// /// constructor and destructor of expression term structure // -expressionTerm::expressionTerm() +expressionTerm::expressionTerm() { type =TERM_UNKNOWN; mult =1.0; @@ -62,7 +62,7 @@ dv_constraint::dv_constraint() penalty_over=1.0; slack_ind1=DOESNT_EXIST; slack_ind2=DOESNT_EXIST; - + penalty_value=0; nConditions=0; @@ -70,7 +70,7 @@ dv_constraint::dv_constraint() } dv_constraint::~dv_constraint() { - delete pExpression; + delete pExpression; delete [] pConditions; nConditions=0; } void dv_constraint::AddCondition(dv_condition* pCond) @@ -148,7 +148,7 @@ bool CDemandOptimizer::ConvertToExpressionTerm(const string s, expressionTerm* t size_t i1=s.find("["); size_t i2=s.find("]"); size_t i3=s.find("("); - + if ((i1 != NPOS) && (i2 != NPOS) && (i3==NPOS)) { //array NOT nested within function string contents=s.substr(i1+1,i2-i1); cout<<" BRACKET CONTENTS: ["<DV_ind = GetDVColumnInd(DV_QOUTRES, _aResIndices[p]); - } + term->DV_ind = GetDVColumnInd(DV_QOUTRES, _aResIndices[p]); + } else if (_aSBIndices[p] != DOESNT_EXIST) { - term->DV_ind=GetDVColumnInd(DV_QOUT,_aSBIndices[p]); + term->DV_ind=GetDVColumnInd(DV_QOUT,_aSBIndices[p]); } else { - //invalid p - likely not enabled + //invalid p - likely not enabled } - } + } else if (s[1]=='D') { term->DV_ind=GetDVColumnInd(DV_DELIVERY,d); - } + } else if (s[1]=='h'){ - term->DV_ind=GetDVColumnInd(DV_STAGE,_aResIndices[p]); - } + term->DV_ind=GetDVColumnInd(DV_STAGE,_aResIndices[p]); + } term->type=TERM_DV; - } + } //---------------------------------------------------------------------- else if (s.find("@ts(") != NPOS) //time series (e.g., @ts(my_time_series,n) { @@ -235,7 +235,7 @@ bool CDemandOptimizer::ConvertToExpressionTerm(const string s, expressionTerm* t warn="ConvertToExpressionTerm: missing end parentheses in @ts expression"+warnstring; ExitGracefully(warn.c_str(), BAD_DATA_WARN); } - if ((is != NPOS) && (ie != NPOS)) + if ((is != NPOS) && (ie != NPOS)) { bool found=false; name = s.substr(is+4,ie-is); @@ -246,7 +246,7 @@ bool CDemandOptimizer::ConvertToExpressionTerm(const string s, expressionTerm* t term->pTS = _pUserTimeSeries[i]; found=true; } - } + } term->timeshift=index; term->type =TERM_TS; if (!found) { @@ -271,7 +271,7 @@ bool CDemandOptimizer::ConvertToExpressionTerm(const string s, expressionTerm* t warn="ConvertToExpressionTerm: missing end paretheses in @lookup expression"+warnstring; ExitGracefully(warn.c_str(), BAD_DATA_WARN); } - if ((is != NPOS) && (ie != NPOS)) + if ((is != NPOS) && (ie != NPOS)) { bool found=false; name = s.substr(is+8,ie-is); @@ -282,7 +282,7 @@ bool CDemandOptimizer::ConvertToExpressionTerm(const string s, expressionTerm* t term->pLT = _pUserLookupTables[i]; found=true; } - } + } term->nested_exp1 =x_in; term->type =TERM_LT; if (!found) { @@ -307,7 +307,7 @@ bool CDemandOptimizer::ConvertToExpressionTerm(const string s, expressionTerm* t warn="ConvertToExpressionTerm: missing end paretheses in @max expression"+warnstring; ExitGracefully(warn.c_str(), BAD_DATA_WARN); } - if ((is != NPOS) && (ie != NPOS)) + if ((is != NPOS) && (ie != NPOS)) { x_in = s.substr(is+5,ie-is); y_in = s.substr(ie+1,ip-ie); @@ -332,7 +332,7 @@ bool CDemandOptimizer::ConvertToExpressionTerm(const string s, expressionTerm* t warn="ConvertToExpressionTerm: missing end paretheses in @min expression"+warnstring; ExitGracefully(warn.c_str(), BAD_DATA_WARN); } - if ((is != NPOS) && (ie != NPOS)) + if ((is != NPOS) && (ie != NPOS)) { x_in = s.substr(is+5,ie-is); y_in = s.substr(ie+1,ip-ie); @@ -357,12 +357,12 @@ bool CDemandOptimizer::ConvertToExpressionTerm(const string s, expressionTerm* t warn="ConvertToExpressionTerm: missing end parentheses in @convert expression"+warnstring; ExitGracefully(warn.c_str(), BAD_DATA_WARN); } - if ((is != NPOS) && (ie != NPOS)) + if ((is != NPOS) && (ie != NPOS)) { bool found=false; x_in = s.substr(is+8,ie-is); units = s.substr(ie+1,ip-ie); - + if (units == "ACREFTD_TO_CMS") {term->value=1.0/ACREFTD_PER_CMS; found=true;} if (units == "CMS_TO_ACREFTD") {term->value=ACREFTD_PER_CMS; found=true;} if (units == "INCHES_TO_MM" ) {term->value=MM_PER_INCH; found=true;} @@ -398,7 +398,7 @@ bool CDemandOptimizer::ConvertToExpressionTerm(const string s, expressionTerm* t } } //---------------------------------------------------------------------- - else if (GetNamedConstant(s)!=RAV_BLANK_DATA) // named constant + else if (GetNamedConstant(s)!=RAV_BLANK_DATA) // named constant { term->type=TERM_CONST; term->value=GetNamedConstant(s); @@ -431,7 +431,7 @@ string expTypeToString(termtype &typ){ } return "?"; } -void SummarizeExpression(const char **s, const int Len, expressionStruct* exp) +void SummarizeExpression(const char **s, const int Len, expressionStruct* exp) { cout<<"*"<')){ - type[i] = EXP_OP; - } + type[i] = EXP_OP; + } if (s[i-1][0]=='/'){type[i]=EXP_INV; } if ((s[i][0]=='=') || (s[i][0]=='~')){ - comp=COMPARE_IS_EQUAL; rhs_ind=i; + comp=COMPARE_IS_EQUAL; rhs_ind=i; } - if (s[i][0]=='<'){//works for <=, too - comp=COMPARE_LESSTHAN; rhs_ind=i; + if (s[i][0]=='<'){//works for <=, too + comp=COMPARE_LESSTHAN; rhs_ind=i; } if (s[i][0]=='>'){ //works for >=, too - comp=COMPARE_GREATERTHAN; rhs_ind=i; + comp=COMPARE_GREATERTHAN; rhs_ind=i; } } if (comp == COMPARE_BETWEEN) { @@ -492,9 +492,9 @@ expressionStruct *CDemandOptimizer::ParseExpression(const char **s, int k=0; int DVcount=0; bool valid; - for (int i = 1; i < Len; i++) + for (int i = 1; i < Len; i++) { - if (type[i] != EXP_OP) + if (type[i] != EXP_OP) { terms[j][k] = new expressionTerm(); valid=ConvertToExpressionTerm(to_string(s[i]),terms[j][k],lineno,filename); @@ -511,8 +511,8 @@ expressionStruct *CDemandOptimizer::ParseExpression(const char **s, if (type[i] == EXP_INV) { terms[j][k]->reciprocal=true; } - //Handle nested expressions - if (terms[j][k]->nested_exp1 != "") + //Handle nested expressions + if (terms[j][k]->nested_exp1 != "") { terms[j][k] = new expressionTerm(); valid=ConvertToExpressionTerm(terms[j][k]->nested_exp1, terms[j][k + 1],lineno, filename); @@ -522,7 +522,7 @@ expressionStruct *CDemandOptimizer::ParseExpression(const char **s, //TODO - THROW WARNING IF DV IS IN NESTED EXP k++; } - if (terms[j][k]->nested_exp2 != "") + if (terms[j][k]->nested_exp2 != "") { terms[j][k] = new expressionTerm(); valid=ConvertToExpressionTerm(terms[j][k]->nested_exp2, terms[j][k + 1],lineno, filename); @@ -534,13 +534,13 @@ expressionStruct *CDemandOptimizer::ParseExpression(const char **s, } k++; } - else if (i!=1){ //operator - expression product term finished and expression doesn't start with + or - + else if (i!=1){ //operator - expression product term finished and expression doesn't start with + or - termspergrp[j]=k; j++; k=0; DVcount=0; } } - //copy to new structure + //copy to new structure int nGroups=j; expressionStruct *tmp; tmp=new expressionStruct(); @@ -561,13 +561,13 @@ expressionStruct *CDemandOptimizer::ParseExpression(const char **s, } void CDemandOptimizer::CheckConditions(const int ii, const time_struct &tt) -{ +{ double dv_value; int p; dv_constraint *pC=_pConstraints[ii]; _pConstraints[ii]->conditions_satisfied=false; - //Check if conditionals are satisfied - for (int j = 0; j < pC->nConditions; j++) + //Check if conditionals are satisfied + for (int j = 0; j < pC->nConditions; j++) { dv_value=0.0; if (pC->pConditions[j]->dv_name == "DAY_OF_YEAR") { @@ -576,7 +576,7 @@ void CDemandOptimizer::CheckConditions(const int ii, const time_struct &tt) else if (pC->pConditions[j]->dv_name == "MONTH") { dv_value=(double)(tt.month); } - else if (pC->pConditions[j]->dv_name[0] == '!') { //decision variable + else if (pC->pConditions[j]->dv_name[0] == '!') { //decision variable char tmp =pC->pConditions[j]->dv_name[1]; string tmp2=pC->pConditions[j]->dv_name.substr(2); long ind=s_to_l(tmp2.c_str()); @@ -595,7 +595,7 @@ void CDemandOptimizer::CheckConditions(const int ii, const time_struct &tt) //int d=GetDemandIndex(tmp2); //int p=_aDemandBasin[d]; //dv_value = _pModel->GetSubBasinByID(p)->GetDemandDelivery(); - //TODO + //TODO } //TODO : handle user specified DVs } @@ -614,15 +614,15 @@ void CDemandOptimizer::CheckConditions(const int ii, const time_struct &tt) _pConstraints[ii]->conditions_satisfied=true; } ////////////////////////////////////////////////////////////////// -/// adds constraint ii to LP solve problem statement +/// adds constraint ii to LP solve problem statement /// \params ii [in] - index of constraint in _pConstraints array /// \param pLinProg [out] - pointer to valid lpsolve structure to be modified -/// \param tt [in] - time structure -/// \param *col_ind [in] - empty array (with memory reserved) for storing column indices +/// \param tt [in] - time structure +/// \param *col_ind [in] - empty array (with memory reserved) for storing column indices /// \param *row_val [in] - empty array (with memory reserved) for storing row values // #ifdef _LPSOLVE_ -void CDemandOptimizer::AddConstraintToLP(const int ii, lp_lib::lprec* pLinProg, const time_struct &tt, int *col_ind, double *row_val) const +void CDemandOptimizer::AddConstraintToLP(const int ii, lp_lib::lprec* pLinProg, const time_struct &tt, int *col_ind, double *row_val) const { dv_constraint *pC=_pConstraints[ii]; @@ -632,11 +632,11 @@ void CDemandOptimizer::AddConstraintToLP(const int ii, lp_lib::lprec* pLinProg, double RHS; bool group_has_dv; int DV_ind; - + int nn=(int)((tt.model_time+TIME_CORRECTION)/1.0);//Options.timestep; //TMP DEBUG - + expressionStruct *pE= pC->pExpression; - + RHS=0.0; for (int j = 0; j < pE->nGroups; j++) { @@ -654,7 +654,7 @@ void CDemandOptimizer::AddConstraintToLP(const int ii, lp_lib::lprec* pLinProg, if (pE->pTerms[j][k]->reciprocal == true) { //TODO: SHOULD CHECK FOR DIVIDE BY ZERO, BUT DO WHAT? coeff /= EvaluateTerm(pE->pTerms[j], k, nn)/(pE->pTerms[j][k]->mult); - } + } else { coeff *= (pE->pTerms[j][k]->mult) * EvaluateTerm(pE->pTerms[j], k, nn); } @@ -692,9 +692,9 @@ void CDemandOptimizer::AddConstraintToLP(const int ii, lp_lib::lprec* pLinProg, retval = lp_lib::add_constraintex(pLinProg,i,row_val,col_ind,constr_type,RHS); ExitGracefullyIf(retval==0,"AddConstraintToLP::Error adding user-specified constraint/goal",RUNTIME_ERR); } -#endif +#endif ////////////////////////////////////////////////////////////////// -/// evaluates term in constraint/goal expression to numerical value (unless term is DV) +/// evaluates term in constraint/goal expression to numerical value (unless term is DV) /// \params pTerms [in] - pointer to array of terms in expression group (e.g., (A*B*C(D,E) stored as [A,B,C,D,E]) /// \param k [in] index of term to be evaluated (e.g., 2 would correspond to C) /// \param nn [in] current timestep index @@ -707,8 +707,8 @@ double CDemandOptimizer::EvaluateTerm(expressionTerm **pTerms,const int k, const int kk; int p; const expressionTerm *pT=pTerms[k]; - - if (pT->type == TERM_DV) + + if (pT->type == TERM_DV) { //do nothing - can't evaluate return -1; @@ -718,22 +718,22 @@ double CDemandOptimizer::EvaluateTerm(expressionTerm **pTerms,const int k, const nshift=pT->timeshift; return pT->pTS->GetValue(nn+nshift); } - else if (pT->type == TERM_LT) + else if (pT->type == TERM_LT) { kk=pT->nested_ind1; x=EvaluateTerm(pTerms,kk,nn); - return pT->pLT->GetValue(x); //decision variable can't be in lookup table + return pT->pLT->GetValue(x); //decision variable can't be in lookup table } - else if (pT->type == TERM_CONST) + else if (pT->type == TERM_CONST) { return pT->value; } else if (pT->type == TERM_HISTORY) //e.g., !Q100[-3] - { + { nshift=pT->timeshift; char tmp =pT->origexp[1]; string tmp2=pT->origexp.substr(2); - p=_pModel->GetSubBasinIndex(s_to_l(tmp2.c_str())); + p=_pModel->GetSubBasinIndex(s_to_l(tmp2.c_str())); if (tmp=='Q'){return _aQhist[p][nshift]; } if (tmp=='h'){return _ahhist[p][nshift]; } if (tmp=='D'){return _aDhist[p][nshift]; } @@ -747,10 +747,10 @@ double CDemandOptimizer::EvaluateTerm(expressionTerm **pTerms,const int k, const x=EvaluateTerm(pTerms,pT->nested_ind1,nn); y=EvaluateTerm(pTerms,pT->nested_ind2,nn); return min(x,y); - } + } else if (pT->type == TERM_CONVERT) { x=EvaluateTerm(pTerms,pT->nested_ind1,nn); return x*pT->value; - } + } return 0; -} \ No newline at end of file +} diff --git a/src/DemandOptimization.cpp b/src/DemandOptimization.cpp index 02323f45..14b87b91 100644 --- a/src/DemandOptimization.cpp +++ b/src/DemandOptimization.cpp @@ -11,7 +11,7 @@ ////////////////////////////////////////////////////////////////// /// \brief Implementation of the Demand optimization constructor -/// \details Creates empty instance of the demand optimization class +/// \details Creates empty instance of the demand optimization class // CDemandOptimizer::CDemandOptimizer(CModel *pMod) { @@ -38,7 +38,7 @@ CDemandOptimizer::CDemandOptimizer(CModel *pMod) _nUserConstants=0; _aUserConstNames=NULL; _aUserConstants=NULL; - + _nUserTimeSeries=0; _pUserTimeSeries=NULL; @@ -50,7 +50,7 @@ CDemandOptimizer::CDemandOptimizer(CModel *pMod) _aDhist=NULL; _ahhist=NULL; - _aCumDelDate=NULL; + _aCumDelDate=NULL; _aCumDelivery=NULL; _nConstraints=0; @@ -60,12 +60,12 @@ CDemandOptimizer::CDemandOptimizer(CModel *pMod) ////////////////////////////////////////////////////////////////// /// \brief Implementation of the Demand optimization destructor // -CDemandOptimizer::~CDemandOptimizer() +CDemandOptimizer::~CDemandOptimizer() { - for (int i=0;i<_nDecisionVars; i++){delete _pDecisionVars[i]; }delete [] _pDecisionVars; - for (int i=0;i<_nUserTimeSeries; i++){delete _pUserTimeSeries[i]; }delete [] _pUserTimeSeries; + for (int i=0;i<_nDecisionVars; i++){delete _pDecisionVars[i]; }delete [] _pDecisionVars; + for (int i=0;i<_nUserTimeSeries; i++){delete _pUserTimeSeries[i]; }delete [] _pUserTimeSeries; for (int i=0;i<_nUserLookupTables;i++){delete _pUserLookupTables[i];}delete [] _pUserLookupTables; - for (int p = 0; p < _nEnabledSubBasins;p++){ + for (int p = 0; p < _nEnabledSubBasins;p++){ if (_aQhist!=NULL){ delete [] _aQhist[p]; delete [] _ahhist[p]; @@ -87,19 +87,19 @@ CDemandOptimizer::~CDemandOptimizer() ////////////////////////////////////////////////////////////////// /// \brief gets demand index d from name /// \params demand_tag [in] - demand ID value (e.g., '123') or Alias (e.g., 'FarmerBob' ) -WITHOUT !D or !D. start -/// \returns demand index d or DOESNT_EXIST if invalid dname +/// \returns demand index d or DOESNT_EXIST if invalid dname /// must be called after _aDemandIDs is populated in CDemandOptimizer::Initialize() // -int CDemandOptimizer::GetDemandIndexFromName(const string demand_tag) const -{ +int CDemandOptimizer::GetDemandIndexFromName(const string demand_tag) const +{ for (int d = 0; d < _nDemands; d++) { if (s_to_i(demand_tag.c_str()) == _aDemandIDs[d]) {return d;} if (demand_tag == _aDemandAliases[d]) {return d;} - } + } return DOESNT_EXIST; } ////////////////////////////////////////////////////////////////// -/// \brief sets history length +/// \brief sets history length /// \params n [in] - number of timesteps in history // void CDemandOptimizer::SetHistoryLength(const int n) @@ -107,11 +107,11 @@ void CDemandOptimizer::SetHistoryLength(const int n) _nHistoryItems=n; } ////////////////////////////////////////////////////////////////// -/// \brief date at which cumulative delivery is rebooted to zero (usually winter date) +/// \brief date at which cumulative delivery is rebooted to zero (usually winter date) /// \params julian_date [in] reboot date /// \params demandID [in] demand identifier (e.g., 'D123') // -void CDemandOptimizer::SetCumulativeDate(const int julian_date, const string demandID) +void CDemandOptimizer::SetCumulativeDate(const int julian_date, const string demandID) { int d; d=GetDemandIndexFromName(demandID); @@ -123,10 +123,10 @@ void CDemandOptimizer::SetCumulativeDate(const int julian_date, const string dem } } ////////////////////////////////////////////////////////////////// -/// \brief adds decision variable structure to _pDecisionVars member aarray -/// \params pDV [in] - user-specified decision variable to be added +/// \brief adds decision variable structure to _pDecisionVars member aarray +/// \params pDV [in] - user-specified decision variable to be added // -void CDemandOptimizer::AddDecisionVar(const decision_var* pDV) +void CDemandOptimizer::AddDecisionVar(const decision_var* pDV) { if (!DynArrayAppend((void**&)(_pDecisionVars),(void*)(pDV),_nDecisionVars)){ ExitGracefully("CDemandOptimizer::AddDecisionVar: adding NULL DV",BAD_DATA);} @@ -136,9 +136,9 @@ void CDemandOptimizer::AddDecisionVar(const decision_var* pDV) } } ////////////////////////////////////////////////////////////////// -/// \brief sets bounds for user specified decision variable +/// \brief sets bounds for user specified decision variable // -void CDemandOptimizer::SetDecisionVarBounds(const string name, const double& min, const double& max) +void CDemandOptimizer::SetDecisionVarBounds(const string name, const double& min, const double& max) { for (int i = 0; i < _nDecisionVars; i++) { if (_pDecisionVars[i]->name==name){ @@ -151,9 +151,9 @@ void CDemandOptimizer::SetDecisionVarBounds(const string name, const double& min WriteWarning(warn.c_str(),_pModel->GetOptStruct()->noisy); } ////////////////////////////////////////////////////////////////// -/// \brief adds user constants +/// \brief adds user constants // -void CDemandOptimizer::AddUserConstant(const string name, const double& val) +void CDemandOptimizer::AddUserConstant(const string name, const double& val) { string *tmp = new string[_nUserConstants+1]; double *tmp2 = new double[_nUserConstants+1]; @@ -170,7 +170,7 @@ void CDemandOptimizer::AddUserConstant(const string name, const double& val) _aUserConstants=tmp2; } ////////////////////////////////////////////////////////////////// -/// \brief adds user time series +/// \brief adds user time series void CDemandOptimizer::AddUserTimeSeries(const CTimeSeries *pTS) { @@ -178,7 +178,7 @@ void CDemandOptimizer::AddUserTimeSeries(const CTimeSeries *pTS) ExitGracefully("CDemandOptimizer::AddUserTimeSeries: adding NULL time series",BAD_DATA);} } ////////////////////////////////////////////////////////////////// -/// \brief adds user lookup table +/// \brief adds user lookup table // void CDemandOptimizer::AddUserLookupTable(const CLookupTable *pLUT) { @@ -200,12 +200,12 @@ void CDemandOptimizer::SetDemandPenalty(const string dname, const double& pen) { } } ////////////////////////////////////////////////////////////////// -/// \brief adds decision variable constraint OR goal to _pConstraints member array +/// \brief adds decision variable constraint OR goal to _pConstraints member array /// \params name [in] - constraint name /// \params exp [in] - expression structure /// \param soft_constraint [in] - true if this is a goal, not a constraint -/// -dv_constraint *CDemandOptimizer::AddConstraint(string name, expressionStruct *exp, bool soft_constraint) +/// +dv_constraint *CDemandOptimizer::AddConstraint(string name, expressionStruct *exp, bool soft_constraint) { dv_constraint *pC=new dv_constraint(); pC->name=name; @@ -213,17 +213,17 @@ dv_constraint *CDemandOptimizer::AddConstraint(string name, expressionStruct *ex pC->pExpression=exp; if (!DynArrayAppend((void**&)(_pConstraints),(void*)(pC),_nConstraints)){ - ExitGracefully("CDemandOptimizer::AddConstraint: adding NULL constraint",BAD_DATA);} + ExitGracefully("CDemandOptimizer::AddConstraint: adding NULL constraint",BAD_DATA);} return _pConstraints[_nConstraints-1]; } ////////////////////////////////////////////////////////////////// /// \brief Initializes Demand optimization instance -/// \notes to be called after .rvh file read and subbasin network initialization, but prior to reading .rvm file -/// \params pModel [in] - pointer to model -/// \params Options [in] - model options structure +/// \notes to be called after .rvh file read and subbasin network initialization, but prior to reading .rvm file +/// \params pModel [in] - pointer to model +/// \params Options [in] - model options structure // -void CDemandOptimizer::Initialize(CModel* pModel, const optStruct& Options) +void CDemandOptimizer::Initialize(CModel* pModel, const optStruct& Options) { #ifndef _LPSOLVE_ ExitGracefully("Demand optimization requires compilation with lpsolve library.",RUNTIME_ERR); @@ -232,7 +232,7 @@ void CDemandOptimizer::Initialize(CModel* pModel, const optStruct& Options) string name; CSubBasin *pSB; decision_var *pDV; - + //populate _aSBindices, _aResIndices //------------------------------------------------------------------ _nEnabledSubBasins=0; @@ -268,7 +268,7 @@ void CDemandOptimizer::Initialize(CModel* pModel, const optStruct& Options) AddDecisionVar(pDV); SB_count++; } - } + } int res_count=0; for (int pp=0;ppGetNumSubBasins();pp++) { @@ -336,7 +336,7 @@ void CDemandOptimizer::Initialize(CModel* pModel, const optStruct& Options) p=pModel->GetOrderedSubBasinIndex(pp); pSB=pModel->GetSubBasin(p); if (pSB->IsEnabled()){ - for (int i = 0; i < pSB->GetNumWaterDemands();i++) + for (int i = 0; i < pSB->GetNumWaterDemands();i++) { ID = pSB->GetWaterDemandID (i);//e.g., 132 name = pSB->GetWaterDemandName(i); @@ -350,13 +350,13 @@ void CDemandOptimizer::Initialize(CModel* pModel, const optStruct& Options) } } } - //TODO: Add code to get demand IDs for reservoir extractions - no slack vars for this + //TODO: Add code to get demand IDs for reservoir extractions - no slack vars for this /*for (p = 0; p < pModel->GetNumSubBasins(); p++) { pSB=pModel->GetSubBasin(p); if (pSB->IsEnabled()){ if (pModel->GetSubBasin(p)->GetReservoir()!=NULL){ - for (int i = 0; i < pModel->GetSubBasin(p)->GetReservoir()->GetNumWaterDemands();i++) + for (int i = 0; i < pModel->GetSubBasin(p)->GetReservoir()->GetNumWaterDemands();i++) { ID = pSB->GetReservoir()->GetWaterDemandID (i);//e.g., 132 name = pSB->GetReservoir()->GetWaterDemandName(i); @@ -380,7 +380,7 @@ void CDemandOptimizer::Initialize(CModel* pModel, const optStruct& Options) _aDelivery [d]=0; _aDemandPenalties[d]=1.0; //Default demand penalty (make accessible parameter?) _aCumDelivery [d]=0.0; - _aCumDelDate [d]=0; //Jan 1 + _aCumDelDate [d]=0; //Jan 1 } //print summary @@ -396,7 +396,7 @@ void CDemandOptimizer::Initialize(CModel* pModel, const optStruct& Options) cout<<" # Demands : "<<_nDemands<GetSubBasin(p); SBID= pSB->GetID(); if (pSB->IsEnabled()){ - for (int i = 0; i < pSB->GetNumWaterDemands();i++) + for (int i = 0; i < pSB->GetNumWaterDemands();i++) { - cout << " " << d << ": " << _aDemandIDs[d] << " (alias: "<<_aDemandAliases[d]<<") in basin "<< SBID <(rvn_floor(1.0/Options.timestep+REAL_SMALL));//default = 1 day - } + } _aDhist = NULL; _aQhist = new double *[_nEnabledSubBasins]; _ahhist = new double *[_nEnabledSubBasins]; @@ -449,10 +449,10 @@ void CDemandOptimizer::InitializePostRVMRead(CModel* pModel, const optStruct& Op } } - // create slack variables for environmental min flow goal + // create slack variables for environmental min flow goal //------------------------------------------------------------------ int p; - _nSlackVars = 0; + _nSlackVars = 0; decision_var *pDV; CSubBasin *pSB; for (int pp=0;ppGetNumSubBasins();pp++) @@ -460,9 +460,9 @@ void CDemandOptimizer::InitializePostRVMRead(CModel* pModel, const optStruct& Op p=pModel->GetOrderedSubBasinIndex(pp); pSB=pModel->GetSubBasin(p); if (pSB->IsEnabled()){ - //created for all basins with environmental min flow or non-zero unusable flow pct + //created for all basins with environmental min flow or non-zero unusable flow pct if ((pSB->GetUnusableFlowPercentage()>0) || (pSB->HasEnviroMinFlow())){ - pDV = new decision_var("E" + to_string(pModel->GetSubBasin(p)->GetID()), p, DV_SLACK,_nSlackVars); + pDV = new decision_var("E" + to_string(pModel->GetSubBasin(p)->GetID()), p, DV_SLACK,_nSlackVars); AddDecisionVar(pDV); _nSlackVars++; } @@ -473,14 +473,14 @@ void CDemandOptimizer::InitializePostRVMRead(CModel* pModel, const optStruct& Op //------------------------------------------------------------------ // ISSUE : THIS IS TOTAL POSSIBLE NUMBER OF SLACK VARIABLES; MAY NOT BE IN PLAY DUE TO CONSTRAINTS - for disabled constraints, just set penalty to zero>?? // - for (int j = 0; j < _nConstraints; j++) + for (int j = 0; j < _nConstraints; j++) { if (_pConstraints[j]->is_goal) { if (_pConstraints[j]->pExpression->compare == COMPARE_IS_EQUAL) { _pConstraints[j]->slack_ind1=_nSlackVars; _pConstraints[j]->slack_ind2=_nSlackVars+1; _nSlackVars+=2; - } + } else { _pConstraints[j]->slack_ind1=_nSlackVars; _nSlackVars++; @@ -506,9 +506,9 @@ bool CDemandOptimizer::UserTimeSeriesExists(string TSname) const } return false; } -void TokenizeString(string instring, char **s, int Len) +void TokenizeString(string instring, char **s, int Len) { - // Returns first token + // Returns first token static char str[256]; strcpy(str,instring.c_str()); char *token = strtok(str, " "); @@ -525,9 +525,9 @@ void TokenizeString(string instring, char **s, int Len) Len=i; } ////////////////////////////////////////////////////////////////// -// adds reservoir constraints as 'user-specified' constraints -// these are done as if the corresponding constraints were read in as expressions from the .rvm file -// +// adds reservoir constraints as 'user-specified' constraints +// these are done as if the corresponding constraints were read in as expressions from the .rvm file +// void CDemandOptimizer::AddReservoirConstraints() { int p; @@ -553,8 +553,8 @@ void CDemandOptimizer::AddReservoirConstraints() //------------------------------------------------------------------------ TSname="_ResQmin_" + SBIDs; if (UserTimeSeriesExists(TSname)) - { - expString="!Q" + SBIDs+ " < @ts(" + TSname + ",0)"; + { + expString="!Q" + SBIDs+ " < @ts(" + TSname + ",0)"; TokenizeString(expString, s, Len); exp = ParseExpression((const char**)(s), Len, 0, "internal"); pConst=AddConstraint(TSname, exp, true); @@ -566,8 +566,8 @@ void CDemandOptimizer::AddReservoirConstraints() //------------------------------------------------------------------------ TSname="_ResQmax_" + SBIDs; if (UserTimeSeriesExists(TSname)) - { - expString="!Q" + SBIDs+ " > @ts(" + TSname + ",0)"; + { + expString="!Q" + SBIDs+ " > @ts(" + TSname + ",0)"; TokenizeString(expString, s, Len); exp = ParseExpression((const char**)(s), Len, 0, "internal"); pConst=AddConstraint(TSname, exp, true); @@ -578,8 +578,8 @@ void CDemandOptimizer::AddReservoirConstraints() //------------------------------------------------------------------------ TSname="_ResFlow_" + SBIDs; if (UserTimeSeriesExists(TSname)) - { - expString="!Q" + SBIDs+ " = @ts(" + TSname + ",0)"; + { + expString="!Q" + SBIDs+ " = @ts(" + TSname + ",0)"; TokenizeString(expString, s, Len); exp = ParseExpression((const char**)(s), Len, 0, "internal"); pConst=AddConstraint(TSname, exp, true); @@ -590,57 +590,57 @@ void CDemandOptimizer::AddReservoirConstraints() //------------------------------------------------------------------------ TSname="_MaxQDelta_" + SBIDs; if (UserTimeSeriesExists(TSname)) - { + { expString = "!Q" + SBIDs + " < @ts(" + TSname + ",0) * 86400 +!Q"+SBIDs+"[-1]"; //TBD - check if this is !Q[0] or !Q[-1] - //expString = "!q" + SBIDs + " < @ts(" + TSname + ",0)"; + //expString = "!q" + SBIDs + " < @ts(" + TSname + ",0)"; TokenizeString(expString, s, Len); exp = ParseExpression((const char**)(s), Len, 0, "internal"); pConst=AddConstraint(TSname, exp, true); pConst->penalty_over=4.0; } - // Max Q Decrease + // Max Q Decrease } } - + } ////////////////////////////////////////////////////////////////// // indexing for DV vector -// 0 .. nSB-1 - outflow of stream reach p +// 0 .. nSB-1 - outflow of stream reach p // nSB .. nSB+nRes-1 - reservoir outflows -// nSB+nRes.. nSB+2*Res-1 - reservoir stages -// nSB+2*nRes .. nSB+2*nRes+nDem-1 - satisfied demand -// nSB+2*nRes+nDem .. nSB+2*nRes+nDem+nDV - user-specified decision variables +// nSB+nRes.. nSB+2*Res-1 - reservoir stages +// nSB+2*nRes .. nSB+2*nRes+nDem-1 - satisfied demand +// nSB+2*nRes+nDem .. nSB+2*nRes+nDem+nDV - user-specified decision variables // remaining - slack vars - [#enviro min flow+ # user <=/>= goals + 2* # user == goals + all reservoir goals ] -// -int CDemandOptimizer::GetDVColumnInd(const dv_type typ, const int counter) const +// +int CDemandOptimizer::GetDVColumnInd(const dv_type typ, const int counter) const { - if (typ==DV_QOUT ){return counter+1;} + if (typ==DV_QOUT ){return counter+1;} else if (typ==DV_QOUTRES ){return _nEnabledSubBasins+counter+1;} else if (typ==DV_STAGE ){return _nEnabledSubBasins+_nReservoirs+counter+1;} - else if (typ==DV_DELIVERY){return _nEnabledSubBasins+2*_nReservoirs+counter+1;} + else if (typ==DV_DELIVERY){return _nEnabledSubBasins+2*_nReservoirs+counter+1;} else if (typ==DV_USER ){return _nEnabledSubBasins+2*_nReservoirs+_nDemands+counter+1; } else if (typ==DV_SLACK ){return _nEnabledSubBasins+2*_nReservoirs+_nDemands+_nUserDecisionVars+counter+1;} - + return 0; } ////////////////////////////////////////////////////////////////// -/// \brief Solves demand optimization problem -/// \notes to be called every time step in lieu of routing mass balance -/// \params pModel [in] - pointer to model -/// \params Options [in] - model options structure +/// \brief Solves demand optimization problem +/// \notes to be called every time step in lieu of routing mass balance +/// \params pModel [in] - pointer to model +/// \params Options [in] - model options structure /// \params aSBrunoff [in] - array of current amount of runoff released to each subbasin (m3) [size:nSubBasins] -/// \params t [in] - local model time +/// \params t [in] - local model time // -void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Options, const double *aSBrunoff,const time_struct &tt ) +void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Options, const double *aSBrunoff,const time_struct &tt ) { #ifdef _LPSOLVE_ - + int d,i,p,s; double RHS; int retval; - double demand_penalty_sum=0; //sum (P * Dstar) to be added to global obj function such that min is zero + double demand_penalty_sum=0; //sum (P * Dstar) to be added to global obj function such that min is zero double tstep=Options.timestep; CSubBasin *pSB; @@ -648,17 +648,17 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio double t=tt.model_time; //int nn=tt.nn+1;//end of timestep int nn=static_cast((t+TIME_CORRECTION)/Options.timestep)+1;//end-of timestep index - + int *col_ind=new int [_nDecisionVars]; //index of column to insert value in current row (1:nDV, not zero-indexed) - double *row_val=new double [_nDecisionVars]; //values of row[col_ind] + double *row_val=new double [_nDecisionVars]; //values of row[col_ind] // instantiate linear programming solver // ---------------------------------------------------------------- - lp_lib::lprec *pLinProg; + lp_lib::lprec *pLinProg; pLinProg=lp_lib::make_lp(0,_nDecisionVars); if (pLinProg==NULL){ExitGracefully("Error in SolveDemandProblem(): couldn construct new linear programming model",RUNTIME_ERR);} - + char name[200]; strcpy(name,"RavenLPSolve"); lp_lib::set_lp_name(pLinProg,name); @@ -668,7 +668,7 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio strcpy(name,_pDecisionVars[i]->name.c_str()); lp_lib::set_col_name(pLinProg, i + 1, name); } - + // Set upper bounds of delivery (D) to demand (D*) (preferred to adding constraint) // ---------------------------------------------------------------- d=0; @@ -685,7 +685,7 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio } } } - // Set upper/lower bounds of reservoir stage + // Set upper/lower bounds of reservoir stage // ---------------------------------------------------------------- int res_count=0; double min_stage,max_stage; @@ -694,7 +694,7 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio p =pModel->GetOrderedSubBasinIndex(pp); pSB=pModel->GetSubBasin(p); if (pSB->IsEnabled()){ - if (pSB->GetReservoir() != NULL) + if (pSB->GetReservoir() != NULL) { min_stage=pSB->GetReservoir()->GetMinStage(nn-1); retval=lp_lib::set_lowbo(pLinProg,GetDVColumnInd(DV_STAGE,res_count), min_stage); //TMP DEBUG @@ -710,10 +710,10 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio } } } - // Set bounds of user-specified decision variables + // Set bounds of user-specified decision variables // ---------------------------------------------------------------- for (int i = 0; i < _nDecisionVars; i++) { - + if (_pDecisionVars[i]->max < ALMOST_INF * 0.99) { retval=lp_lib::set_upbo(pLinProg,GetDVColumnInd(DV_USER,i), _pDecisionVars[i]->max); ExitGracefullyIf(retval!=1,"SolveDemandProblem::Error adding decision variable upper bound",RUNTIME_ERR); @@ -757,27 +757,27 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio } } } - //ADD RESERVOIR EXTRACTION HERE - - // enviro min flow goal penalties + //ADD RESERVOIR EXTRACTION HERE + + // enviro min flow goal penalties // ---------------------------------------------------------------- for (int pp = 0; ppGetNumSubBasins(); pp++) { p =pModel->GetOrderedSubBasinIndex(pp); pSB=pModel->GetSubBasin(p); if (pSB->IsEnabled()){ - if ((pSB->GetNumWaterDemands()>0) && + if ((pSB->GetNumWaterDemands()>0) && ((pSB->GetUnusableFlowPercentage()>0) || (pSB->GetEnviroMinFlow(t)>0))) { col_ind[i]=GetDVColumnInd(DV_SLACK,s); - row_val[i]=10.0; //Is this a fair penalty?- needs to be stronger than demand penalty. make a parameter? + row_val[i]=10.0; //Is this a fair penalty?- needs to be stronger than demand penalty. make a parameter? i++; s++; } } } - - // user-specified goal penalties + + // user-specified goal penalties // ---------------------------------------------------------------- - for (int j = 0; j < _nConstraints; j++) + for (int j = 0; j < _nConstraints; j++) { CheckConditions(j,tt); @@ -791,8 +791,8 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio col_ind[i]=GetDVColumnInd(DV_SLACK,s); row_val[i]=_pConstraints[j]->penalty_under; i++; s++; - } - else // only one slack var (>= or <=) + } + else // only one slack var (>= or <=) { col_ind[i]=GetDVColumnInd(DV_SLACK,s); row_val[i]=_pConstraints[j]->penalty_over; @@ -800,14 +800,14 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio } } } - + retval = lp_lib::set_obj_fnex(pLinProg,i, row_val, col_ind); ExitGracefullyIf(retval==0,"SolveDemandProblem: Error specifying objective function",RUNTIME_ERR); - + // ---------------------------------------------------------------- // add mass balance constraints // ---------------------------------------------------------------- - + // river routing mass balance constraints (WORKING) // ---------------------------------------------------------------- d=0; //counter for demands @@ -841,13 +841,13 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio } i=0; - + col_ind[i]=GetDVColumnInd(DV_QOUT,_aSBIndices[p]); row_val[i]=1.0; i++; int id; - for (int in = 0; in < nInlets; in++) + for (int in = 0; in < nInlets; in++) { dvtype=DV_QOUT; id=_aSBIndices[p_in[in]]; if (pModel->GetSubBasin(p_in[in])->GetReservoir() != NULL) { @@ -855,7 +855,7 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio } //ExitGracefullyIf(id==DOESNT_EXIST,"SolveDemandProblem: invalid index",RUNTIME_ERR); col_ind[i]=GetDVColumnInd(dvtype,id); - row_val[i]=-U_0; + row_val[i]=-U_0; i++; } for (int dd = 0; ddGetNumWaterDemands();dd++) @@ -881,7 +881,7 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio ExitGracefullyIf(retval==0,"SolveDemandProblem::Error adding mass balance constraint",RUNTIME_ERR); } } - + // Lake mass balance constraints (WORKING) //---------------------------------------------------------------- CReservoir *pRes; @@ -902,9 +902,9 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio //dV/dh*(hnew-hold)/dt=(Qin-Qout-ET+P-resextract-seep) // A/dt * h -1.0 * Qin + 1.0 * Qout - 1.0*resextract = P-ET-seep +A/dt* hold // must fix area at A^n-1 and stage used for flows at h^n-1 - + h_old=pRes->GetResStage(); - Adt =pRes->GetSurfaceArea()/Options.timestep/SEC_PER_DAY; + Adt =pRes->GetSurfaceArea()/Options.timestep/SEC_PER_DAY; k =pRes->GetHRUIndex(); ET =pModel->GetHydroUnit(k)->GetForcingFunctions()->OW_PET / SEC_PER_DAY / MM_PER_METER; //average for timestep, in m/s ET *=pModel->GetHydroUnit(k)->GetSurfaceProps()->lake_PET_corr; @@ -918,7 +918,7 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio i++; col_ind[i]=GetDVColumnInd(DV_QOUTRES,_aResIndices[p]); - row_val[i]=+1.0; + row_val[i]=+1.0; i++; col_ind[i]=GetDVColumnInd(DV_STAGE ,_aResIndices[p]); @@ -932,7 +932,7 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio //i++; d++; //} - RHS=precip-ET-seepage+Adt*h_old; + RHS=precip-ET-seepage+Adt*h_old; retval = lp_lib::add_constraintex(pLinProg,i,row_val,col_ind,ROWTYPE_EQ,RHS); ExitGracefullyIf(retval==0,"SolveDemandProblem::Error adding mass balance constraint",RUNTIME_ERR); @@ -940,32 +940,32 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio // Second constraint equation (relation between Qout and h): // ------------------------------------------------------------------------ // linearized storage -discharge curve - // Qout = Qoutold + dQ/dh * (h-hold) - // 1.0 * Qout - dQ/dh * h = Qoutold - dQ/dh * hold - + // Qout = Qoutold + dQ/dh * (h-hold) + // 1.0 * Qout - dQ/dh * h = Qoutold - dQ/dh * hold + dQdh =pRes->GetStageDischargeDerivative(h_old,nn); Qout =pRes->GetOutflowRate(); i = 0; col_ind[i]=GetDVColumnInd(DV_QOUTRES,_aResIndices[p]); - row_val[i]=1.0; + row_val[i]=1.0; i++; col_ind[i]=GetDVColumnInd(DV_STAGE ,_aResIndices[p]); - row_val[i]=-dQdh; + row_val[i]=-dQdh; i++; //TODO: THIS SHOULD BE A GOAL INSTEAD. MAY BE OVERRIDDEN BY OTHER CONSTRAINTS (YES!) /*col_ind[i] = GetDVColumnInd(DV_SLACK, s); - row_val[i]=-1.0; + row_val[i]=-1.0; i++; s++; col_ind[i]=GetDVColumnInd(DV_SLACK ,s); - row_val[i]=+1.0; + row_val[i]=+1.0; i++; s++;*/ RHS=Qout-dQdh*h_old; - + retval = lp_lib::add_constraintex(pLinProg,i,row_val,col_ind,ROWTYPE_EQ,RHS); ExitGracefullyIf(retval==0,"SolveDemandProblem::Error adding mass balance constraint",RUNTIME_ERR); @@ -973,7 +973,7 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio } } - // demand goals + // demand goals // ---------------------------------------------------------------- d=0; double minQ,ufp; @@ -983,19 +983,19 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio { p =pModel->GetOrderedSubBasinIndex(pp); pSB=pModel->GetSubBasin(p); - + if (pSB->IsEnabled()) - { - // sum(D) <= (1-ufp)*((Qout+sum(D))-Qmin) + { + // sum(D) <= (1-ufp)*((Qout+sum(D))-Qmin) // (1-ufp) * Q - (ufp) * sum(D) + (1-ufp) * S >= (1-ufp) * Qmin - // important- must have higher penalty than demand undersatisfaction penalty, else enviro min flow could be violated by demand - + // important- must have higher penalty than demand undersatisfaction penalty, else enviro min flow could be violated by demand + ufp =pSB->GetUnusableFlowPercentage(); minQ=pSB->GetEnviroMinFlow(t); i=0; col_ind[i]=GetDVColumnInd(DV_QOUT,_aSBIndices[p]); - row_val[i] = (1.0-ufp); + row_val[i] = (1.0-ufp); i++; for (int dd = 0; ddGetNumWaterDemands();dd++) @@ -1006,7 +1006,7 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio } if (pSB->GetNumWaterDemands()>0){ col_ind[i]=GetDVColumnInd(DV_SLACK,slack_index); - row_val[i]= (1.0-ufp); + row_val[i]= (1.0-ufp); i++; slack_index++; } RHS = (1.0-ufp)*minQ; @@ -1019,16 +1019,16 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio } retval = lp_lib::add_constraintex(pLinProg,i,row_val,col_ind,ROWTYPE_GE,RHS); ExitGracefullyIf(retval==0,"SolveDemandProblem::Error adding environmental flow goal",RUNTIME_ERR); - } + } } } - // ADD RESERVOIR DEMAND HERE + // ADD RESERVOIR DEMAND HERE - // user-specified constraints / goals + // user-specified constraints / goals // ---------------------------------------------------------------- - for (int i = 0; i < _nConstraints; i++) + for (int i = 0; i < _nConstraints; i++) { if (_pConstraints[i]->conditions_satisfied) { @@ -1041,14 +1041,14 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio // ---------------------------------------------------------------- // SOLVE OPTIMIZATION PROBLEM WITH LP_SOLVE // ---------------------------------------------------------------- - lp_lib::set_minim(pLinProg); + lp_lib::set_minim(pLinProg); lp_lib::set_verbose(pLinProg,IMPORTANT); //lp_lib::write_LP(pLinProg); - + //retval=lp_lib::set_scaling(lp, CURTISREIDSCALE); //May wish to look for scaling improvements //ExitGracefullyIf(retval!=0,"SolveDemandProblem: unable to set scaling algorithm",RUNTIME_ERR); - //retval=lp_lib::set_BFP(pLinProg, "bfp_LUSOL"); //TEST THIS FOR SPEED + //retval=lp_lib::set_BFP(pLinProg, "bfp_LUSOL"); //TEST THIS FOR SPEED //ExitGracefullyIf(retval!=0,"SolveDemandProblem: unable to set BFP",RUNTIME_ERR); //lp_lib::print_debugdump(pLinProg,"lp_solve_dump.txt"); @@ -1057,10 +1057,10 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio retval = lp_lib::solve(pLinProg); if (retval!=OPTIMAL){cout<<"code: "<value=soln[i]; } for (int j=0;jp_index; value = _pDecisionVars[i]->value; typ = _pDecisionVars[i]->dvar_type; - + if (typ == DV_QOUT) { //pModel->GetSubBasin(p)->SetQout(value); - } + } else if (typ == DV_STAGE) { //double h_last=pModel->GetSubBasin(p)->GetReservoir()->GetResStage(); //pModel->GetSubBasin(p)->GetReservoir()->SetReservoirStage(value,h_last); - } + } else if (typ == DV_QOUTRES) { - //this is byproduct of stage determination - not passed back to model - } + //this is byproduct of stage determination - not passed back to model + } else if (typ == DV_DELIVERY) { - d = _pDecisionVars[i]->loc_index; + d = _pDecisionVars[i]->loc_index; _aDelivery [d]=value; //cout<<"delivered["<name<GetSubBasin(p)->SetDelivery(d,value); //this can be handled within CSubBasin::ApplyIrrigationDemand instead - } - + //pModel->GetSubBasin(p)->SetDelivery(d,value); //this can be handled within CSubBasin::ApplyIrrigationDemand instead + } + else if (typ == DV_SLACK) { s = _pDecisionVars[i]->loc_index; _aSlackValues[s]=value; @@ -1129,10 +1129,10 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio } else if (typ == DV_USER) { - //doesn't influence hydrologic model + //doesn't influence hydrologic model } } - + delete [] col_ind; delete [] row_val; @@ -1149,7 +1149,7 @@ void CDemandOptimizer::WriteOutputFileHeaders(const optStruct &Options) { CSubBasin *pSB; - // DemandOptimization.csv + // DemandOptimization.csv //-------------------------------------------------------- string tmpFilename=FilenamePrepare("DemandOptimization.csv",Options); _DEMANDOPT.open(tmpFilename.c_str()); @@ -1166,14 +1166,14 @@ void CDemandOptimizer::WriteOutputFileHeaders(const optStruct &Options) } _DEMANDOPT<GetSubBasin(_pDecisionVars[i]->p_index); if ( ((pSB!=NULL) && (pSB->IsGauged())) || (pSB==NULL) ){ //gauged basin or not associated with basin @@ -1208,7 +1208,7 @@ void CDemandOptimizer::WriteMinorOutput(const optStruct &Options,const time_stru _DEMANDOPT<<","<<_aDelivery[0]; _DEMANDOPT<GetNumSubBasins();pp++){ p = _pModel->GetOrderedSubBasinIndex(pp); pSB=_pModel->GetSubBasin(p); - for (int ii = 0; ii < pSB->GetNumWaterDemands(); ii++) + for (int ii = 0; ii < pSB->GetNumWaterDemands(); ii++) { if (include_pen){mult=_aDemandPenalties[d];} _GOALSAT << "," << mult*(pSB->GetWaterDemand(ii,tt.model_time)-_aDelivery[d]); d++; } } - int k=_nDemands; //first goals are due to environmental flow constraints + int k=_nDemands; //first goals are due to environmental flow constraints for (int i = 0; i < _nConstraints; i++) { if (_pConstraints[i]->is_goal) { if (_pConstraints[i]->pExpression->compare != COMPARE_IS_EQUAL) { @@ -1235,7 +1235,7 @@ void CDemandOptimizer::WriteMinorOutput(const optStruct &Options,const time_stru pen = mult*_aSlackValues[k]; _GOALSAT<<","<penalty_over;} pen =mult*_aSlackValues[k]; @@ -1257,4 +1257,4 @@ void CDemandOptimizer::CloseOutputStreams() { if (_DEMANDOPT.is_open()){_DEMANDOPT.close();} if ( _GOALSAT.is_open()){ _GOALSAT.close();} -} \ No newline at end of file +} diff --git a/src/DemandOptimization.h b/src/DemandOptimization.h index 43a4c1b0..535c76f6 100644 --- a/src/DemandOptimization.h +++ b/src/DemandOptimization.h @@ -15,12 +15,12 @@ namespace lp_lib { #include "../lib/lp_solve/lp_lib.h" } -#endif +#endif /////////////////////////////////////////////////////////////////// -/// \brief Data abstraction for general lookup table y(x) +/// \brief Data abstraction for general lookup table y(x) // -class CLookupTable +class CLookupTable { private: string _name; @@ -38,79 +38,79 @@ class CLookupTable }; /////////////////////////////////////////////////////////////////// -/// \brief different expression types +/// \brief different expression types // -enum exptype +enum exptype { EXP, //< default expression EXP_OP, //< operator - EXP_INV //< inverse operator + EXP_INV //< inverse operator }; /////////////////////////////////////////////////////////////////// -/// \brief different expression term types +/// \brief different expression term types // enum termtype { TERM_DV, //< decision variable - TERM_TS, //< time series @ts(name,n) + TERM_TS, //< time series @ts(name,n) TERM_LT, //< lookup table @lookup(x) - TERM_CONST, //< constant - TERM_HISTORY, //< bracketed - !Q123[-2] - TERM_MAX, //< @max(x,y) + TERM_CONST, //< constant + TERM_HISTORY, //< bracketed - !Q123[-2] + TERM_MAX, //< @max(x,y) TERM_MIN, //< @min(x,y) TERM_CONVERT, //< @convert(x,units) - TERM_UNKNOWN //< unknown + TERM_UNKNOWN //< unknown }; /////////////////////////////////////////////////////////////////// -/// \brief decision variable types +/// \brief decision variable types // -enum dv_type +enum dv_type { DV_QOUT, //< outflow from reach - DV_QOUTRES, //< outflow from reservoir + DV_QOUTRES, //< outflow from reservoir DV_STAGE, //< reservoir stage - DV_DELIVERY,//< delivery of water demand + DV_DELIVERY,//< delivery of water demand DV_SLACK, //< slack variable for goal satisfaction - DV_USER //< user specified decision variable + DV_USER //< user specified decision variable }; // data structures used by CDemandOptimizer class // ------------------------------------------------------------------- ////////////////////////////////////////////////////////////////// -/// expression term +/// expression term /// individual term in expression -// -struct expressionTerm +// +struct expressionTerm { - termtype type; //< type of expression + termtype type; //< type of expression double mult; //< multiplier of expression (+/- 1, depending upon operator and location in exp) - bool reciprocal; //< true if term is in denominator - + bool reciprocal; //< true if term is in denominator + double value; //< constant value or conversion multiplier CTimeSeries *pTS; //< pointer to time series, if this is a named time series CLookupTable *pLT; //< pointer to lookup table, if this is a named time series - bool is_nested; //< true if nested within (i.e., an argument to) another term - int timeshift; //< for time series (+ or -) or lookback value (+) - int DV_ind; //< index of decision variable + bool is_nested; //< true if nested within (i.e., an argument to) another term + int timeshift; //< for time series (+ or -) or lookback value (+) + int DV_ind; //< index of decision variable int nested_ind1; //< index k of first argument (e.g., for lookup table with term entry) int nested_ind2; //< index k of second argument (e.g., for min/max functions) string nested_exp1; //< contents of first argument to function - can be expression string nested_exp2; //< contents of second argument to function - - string origexp; //< original string expression - expressionTerm(); + string origexp; //< original string expression + + expressionTerm(); }; ////////////////////////////////////////////////////////////////// /// expression structure -/// abstraction of (A*B*C)+(D*E)+(-F)+(G*H) <= 0 +/// abstraction of (A*B*C)+(D*E)+(-F)+(G*H) <= 0 /// parenthetical collections are groups of terms -example has 4 groups with [3,2,1,2] Terms per group // -struct expressionStruct //full expression -{ +struct expressionStruct //full expression +{ expressionTerm ***pTerms; //< 2D irregular array of pointers to expression terms size:[nGroups][nExpPerGrp[j]] - int nGroups; //< total number of terms groups in expression + int nGroups; //< total number of terms groups in expression int *nTermsPerGrp;//< number of terms per group [size: nGroups] comparison compare; //< comparison operator (==, <, >) @@ -120,34 +120,34 @@ struct expressionStruct //full expression ////////////////////////////////////////////////////////////////// -// decision variable +// decision variable // struct decision_var { - string name; //< decision variable names: Qxxxx or Dxxxx where xxxx is SBID + string name; //< decision variable names: Qxxxx or Dxxxx where xxxx is SBID dv_type dvar_type; //< decision variable type: e.g., DV_QOUT or DV_DELIVERY - int p_index; //< raw subbasin index p (not SBID) of decision variable (or DOESNT_EXIST if not linked to SB) + int p_index; //< raw subbasin index p (not SBID) of decision variable (or DOESNT_EXIST if not linked to SB) int dem_index; //< demand index in subbasin p (or DOESNT_EXIST if type not DV_DELIVERY) (i..nDemands in SB p_index, usually <=1) - int loc_index; //< local index (rescount or subbasincount or demand count) + int loc_index; //< local index (rescount or subbasincount or demand count) double value; //< solution value for decision variable - double min; //< minimum bound (default=0) + double min; //< minimum bound (default=0) double max; //< maximum bound (default unbounded) decision_var(string nam, int p, dv_type typ, int loc_ind) { - name=nam; - p_index=p; - loc_index=loc_ind; + name=nam; + p_index=p; + loc_index=loc_ind; value=0.0;min=0.0;max=ALMOST_INF;dvar_type=typ;dem_index=DOESNT_EXIST; } }; ////////////////////////////////////////////////////////////////// -// decision variable constraint condition +// decision variable constraint condition // struct dv_condition { string dv_name; //< decision variable name (e.g., Q1023) or "MONTH" - double value; //< conditional value + double value; //< conditional value double value2; //< second conditional (if COMPARE_BETWEEN) comparison compare; //> comparison operator, e.g., COMPARE_IS_EQUAL @@ -163,9 +163,9 @@ struct dv_condition // struct dv_constraint { - string name; //< constraint name - expressionStruct *pExpression; //< constraint expression - + string name; //< constraint name + expressionStruct *pExpression; //< constraint expression + bool is_goal; //< true if constraint is soft (goal rather than constraint) int priority; //< priority (default==1, for goals only) double penalty_under; //< penalty if under specified value (for goals only) @@ -173,18 +173,18 @@ struct dv_constraint int slack_ind1; //< slack index (0.._nSlackVars) of under/over slack index for goal, or DOESNT_EXIST if constraint int slack_ind2; //< slack index (0.._nSlackVars) of over slack index for target goal, or DOESNT_EXIST if constraint - double penalty_value; //< (from solution) penalty incurred by not satisfying goal (or zero for constraint) + double penalty_value; //< (from solution) penalty incurred by not satisfying goal (or zero for constraint) - int nConditions; //< number of conditional statments - dv_condition **pConditions; //< array of pointers to conditional statements - bool conditions_satisfied; + int nConditions; //< number of conditional statments + dv_condition **pConditions; //< array of pointers to conditional statements + bool conditions_satisfied; dv_constraint(); ~dv_constraint(); void AddCondition(dv_condition *pCondition); }; /////////////////////////////////////////////////////////////////// -/// \brief Data abstraction for demand optimization +/// \brief Data abstraction for demand optimization // class CDemandOptimizer { @@ -195,37 +195,37 @@ class CDemandOptimizer int _nDecisionVars; //< total number of decision variables considered decision_var **_pDecisionVars; //< array of pointers to decision variable structures [size:_nDecisionVars] - int _nConstraints; //< number of user-defined enforced constraints/goals in management model + int _nConstraints; //< number of user-defined enforced constraints/goals in management model dv_constraint **_pConstraints; //< array of pointers to user-defined enforced constraints/goals in management model - int _nEnabledSubBasins; //< number of enabled subbasins in model + int _nEnabledSubBasins; //< number of enabled subbasins in model int *_aSBIndices; //< local index of enabled subbasins (0:_nEnabledSubBasins) [size: _nSubBasins] - int _nReservoirs; //< local storage of number of simulated lakes/reservoirs - int *_aResIndices; //< storage of enabled reservoir indices (0:_nReservoirs or DOESNT_EXIST) [size:_nSubBasins] + int _nReservoirs; //< local storage of number of simulated lakes/reservoirs + int *_aResIndices; //< storage of enabled reservoir indices (0:_nReservoirs or DOESNT_EXIST) [size:_nSubBasins] int _nDemands; //< local storage of number of demand locations (:IrrigationDemand/:WaterDemand + :ReservoirExtraction time series) int *_aDemandIDs; //< local storage of demand IDs [size:_nDemands] - string *_aDemandAliases; //< list of demand aliases [size: _nDemands] + string *_aDemandAliases; //< list of demand aliases [size: _nDemands] double *_aDemandPenalties; //< array of priority weights [size: _nDemands] double *_aDelivery; //< array of delivered water for each demand [m3/s] [size: _nDemands] - double *_aCumDelivery; //< array of cumulative deliveries of demand since _aCumDivDate of this year [m3] [size: _nDemands] + double *_aCumDelivery; //< array of cumulative deliveries of demand since _aCumDivDate of this year [m3] [size: _nDemands] int *_aCumDelDate; //< julian date to calculate cumulative deliveries from {default: Jan 1)[size: _nDemands] - - //int _nDemandGroups; //< number of demand groups - //CDemandGroup **_pDemandGroups; //< array of pointers to demand groups - + + //int _nDemandGroups; //< number of demand groups + //CDemandGroup **_pDemandGroups; //< array of pointers to demand groups + double *_aSlackValues; //< array of slack variable values [size: _nSlackVars] - int _nSlackVars; //< number of slack variables + int _nSlackVars; //< number of slack variables - int _nUserDecisionVars; //< number of user-specified decision variables + int _nUserDecisionVars; //< number of user-specified decision variables int _nUserConstants; //< number of user-specified named constants string *_aUserConstNames; //< array of names of user-specified constants double *_aUserConstants; //< array of values of user-specified constants int _nUserTimeSeries; //< number of user variable time series CTimeSeries **_pUserTimeSeries; //< array of pointers to user variable time series - int _nUserLookupTables; //< number of user variable lookup tables - CLookupTable **_pUserLookupTables; //< array of pointers to user variable lookup tables + int _nUserLookupTables; //< number of user variable lookup tables + CLookupTable **_pUserLookupTables; //< array of pointers to user variable lookup tables int _nHistoryItems; //< size of flow/demand/stage history that needs to be stored (in time steps) (from :LookbackDuration) double **_aQhist; //< history of subbasin discharge [size: nSBs * _nHistoryItems] @@ -240,10 +240,10 @@ class CDemandOptimizer #ifdef _LPSOLVE_ void AddConstraintToLP(const int i, lp_lib::lprec *pLinProg, const time_struct &tt,int *col_ind, double *row_val) const; -#endif +#endif double EvaluateTerm(expressionTerm **pTerms,const int k, const int nn) const; - void CheckConditions(const int ii, const time_struct &tt); + void CheckConditions(const int ii, const time_struct &tt); bool UserTimeSeriesExists(string TSname) const; void AddReservoirConstraints(); @@ -253,20 +253,20 @@ class CDemandOptimizer ~CDemandOptimizer(); int GetDemandIndexFromName(const string dname) const; - double GetNamedConstant (const string s) const; + double GetNamedConstant (const string s) const; double GetTotalWaterDemandDelivery(const int p) const; int GetNumUserDVs () const; void SetHistoryLength (const int n); void SetCumulativeDate (const int julian_date, const string demandID); - + void AddDecisionVar (const decision_var *pDV); void SetDecisionVarBounds (const string name, const double &min, const double &max); dv_constraint *AddConstraint (const string name, expressionStruct *exp, const bool soft_constraint); void AddUserConstant (const string name, const double &val); void AddUserTimeSeries (const CTimeSeries *pTS); void AddUserLookupTable (const CLookupTable *pLUT); - + //void MultiplyDemand (const string dname, const double &mult); //void MultiplyGroupDemand (const string groupname, const double &mult); void SetDemandPenalty (const string dname, const double &pen); @@ -275,11 +275,11 @@ class CDemandOptimizer expressionStruct *ParseExpression(const char **s, const int Len, const int lineno, const string filename) const; void Initialize (CModel *pModel, const optStruct &Options); - void InitializePostRVMRead (CModel *pModel, const optStruct &Options); + void InitializePostRVMRead (CModel *pModel, const optStruct &Options); void SolveDemandProblem (CModel *pModel, const optStruct &Options, const double *aSBrunoff, const time_struct &tt); void WriteOutputFileHeaders(const optStruct &Options); void WriteMinorOutput (const optStruct &Options,const time_struct &tt); void CloseOutputStreams (); }; -#endif \ No newline at end of file +#endif diff --git a/src/ForcingGrid.cpp b/src/ForcingGrid.cpp index 8025aa21..dec5f7eb 100644 --- a/src/ForcingGrid.cpp +++ b/src/ForcingGrid.cpp @@ -165,14 +165,14 @@ CForcingGrid::CForcingGrid( const CForcingGrid &grid ) _CellIDToIdx = new int [_nCells]; ExitGracefullyIf(_CellIDToIdx==NULL,"CForcingGrid::Copy Constructor(8)",OUT_OF_MEMORY); for(int c=0; c<_nCells; c++) { // loop over non-zero weighted cells - _CellIDToIdx[c]=grid._CellIDToIdx[c]; + _CellIDToIdx[c]=grid._CellIDToIdx[c]; } _IdxNonZeroGridCells = NULL; _IdxNonZeroGridCells = new int [_nNonZeroWeightedGridCells]; ExitGracefullyIf(_IdxNonZeroGridCells==NULL,"CForcingGrid::Copy Constructor(3)",OUT_OF_MEMORY); for (int ic=0; ic<_nNonZeroWeightedGridCells; ic++) { // loop over non-zero weighted cells - _IdxNonZeroGridCells[ic]=grid._IdxNonZeroGridCells[ic]; + _IdxNonZeroGridCells[ic]=grid._IdxNonZeroGridCells[ic]; } _aLatitude=NULL;_aLongitude=NULL;_aElevation=NULL;_aStationIDs=NULL; diff --git a/src/Model.cpp b/src/Model.cpp index 66a1124a..b7da2715 100644 --- a/src/Model.cpp +++ b/src/Model.cpp @@ -1615,7 +1615,7 @@ void CModel::SetEnsembleMode(CEnsemble *pEnsemble) { _pEnsemble=pEnsemble; } ////////////////////////////////////////////////////////////////// -/// \brief Sets demand optimizer +/// \brief Sets demand optimizer /// /// \param *pDO [in] (assumed valid) pointer to optimizer instance // diff --git a/src/Model.h b/src/Model.h index 5b1988e8..941ba014 100644 --- a/src/Model.h +++ b/src/Model.h @@ -111,7 +111,7 @@ class CModel: public CModelABC CGroundwaterModel *_pGWModel; ///< pointer to corresponding groundwater model CTransportModel *_pTransModel; ///< pointer to corresponding transport model CEnsemble *_pEnsemble; ///< pointer to model ensemble - CDemandOptimizer *_pDO; ///< pointer to demand optimizer (NULL w/o optimization) + CDemandOptimizer *_pDO; ///< pointer to demand optimizer (NULL w/o optimization) //For Diagnostics Calculation CTimeSeriesABC **_pObservedTS; ///< array of pointers of observation time series [size: _nObservedTS] diff --git a/src/ParseHRUFile.cpp b/src/ParseHRUFile.cpp index 25ba915a..772bb997 100644 --- a/src/ParseHRUFile.cpp +++ b/src/ParseHRUFile.cpp @@ -830,7 +830,7 @@ bool ParseHRUPropsFile(CModel *&pModel, const optStruct &Options, bool terrain_r int iter=0; for(int p=0;pGetNumSubBasins();p++) { - if ((pModel->IsSubBasinUpstream(pModel->GetSubBasin(p)->GetID(),SBID)) || + if ((pModel->IsSubBasinUpstream(pModel->GetSubBasin(p)->GetID(),SBID)) || (pModel->GetSubBasin(p)->GetID() == SBID) ) { pSBGroup->AddSubbasin(pModel->GetSubBasin(p)); advice=advice+to_string(pModel->GetSubBasin(p)->GetID())+" "; diff --git a/src/ParseLib.cpp b/src/ParseLib.cpp index 9f7e67ef..b8cb00f7 100644 --- a/src/ParseLib.cpp +++ b/src/ParseLib.cpp @@ -54,7 +54,7 @@ string CParser::Peek() string firstword = ""; if (Len > 0) {firstword=s[0];} if (!eof){ - INPUT->seekg(place ,std::ios_base::beg); // Return to position before peeked line + INPUT->seekg(place ,std::ios_base::beg); // Return to position before peeked line } return firstword; } diff --git a/src/ParseManagementFile.cpp b/src/ParseManagementFile.cpp index 69a8268e..fa70eb8c 100644 --- a/src/ParseManagementFile.cpp +++ b/src/ParseManagementFile.cpp @@ -7,11 +7,11 @@ Copyright (c) 2008-2024 the Raven Development Team #include "ParseLib.h" #include "DemandOptimization.h" -void SummarizeExpression(const char **s, const int Len, expressionStruct* exp); //defined in DemandExpressionHandling.cpp +void SummarizeExpression(const char **s, const int Len, expressionStruct* exp); //defined in DemandExpressionHandling.cpp ////////////////////////////////////////////////////////////////// /// \brief Parses Demand Management file -/// \details model.rvm: input file that defines management optimization problem and solution options +/// \details model.rvm: input file that defines management optimization problem and solution options /// /// \param *&pModel [out] Reference to model object /// \param &Options [out] Global model options information @@ -159,9 +159,9 @@ bool ParseManagementFile(CModel *&pModel,const optStruct &Options) break; } case(11): //---------------------------------------------- - {/*:DemandGroup [groupname] + {/*:DemandGroup [groupname] [demand1] [demand2] ... [demandN] - :EndDemandGroup + :EndDemandGroup */ if(Options.noisy) { cout <<"Demand Group"<AddDemandGroup(s[1]); @@ -208,20 +208,20 @@ bool ParseManagementFile(CModel *&pModel,const optStruct &Options) case(21): //---------------------------------------------- { /*:DefineDecisionVariable [name] = [expressionRHS] */ if(Options.noisy) { cout <<"Define Decision Variable"<GetNumUserDVs()); - + pDO->AddDecisionVar(pDV); pExp=pDO->ParseExpression((const char**)(s),Len,pp->GetLineNumber(),pp->GetFilename()); if (Options.noisy){ - SummarizeExpression((const char**)(s),Len,pExp); + SummarizeExpression((const char**)(s),Len,pExp); } if (pExp!=NULL){ pConst = pDO->AddConstraint(s[1], pExp, false); - } + } else { string warn ="Invalid expression in :DefineDecisionVariable command at line " + pp->GetLineNumber(); WriteWarning(warn.c_str(),Options.noisy); @@ -235,22 +235,22 @@ bool ParseManagementFile(CModel *&pModel,const optStruct &Options) break; } case(23): //---------------------------------------------- - {/*:ManagementConstraint [name] + {/*:ManagementConstraint [name] :Expresion [expression] :Conditional [condition] :Conditional [condition] :EndManagementConstraint or - :ManagementGoal [name] + :ManagementGoal [name] :Expresion [expression] :Conditional [condition] :Conditional [condition] - :Penalty [value] {value2} + :Penalty [value] {value2} :EndManagementGoal */ if(Options.noisy) { cout <<"Management Constraint or Management Goal"<Peek(); if (firstword == ":Expression") { @@ -262,7 +262,7 @@ bool ParseManagementFile(CModel *&pModel,const optStruct &Options) if(Options.noisy) { cout << "-->reading line " << pp->GetLineNumber() << ": "; } if (Len == 0) { if(Options.noisy) { cout << "#" << endl; } }//Do nothing else if(IsComment(s[0],Len)) { if(Options.noisy) { cout << "#" << endl; } } - else if(!strcmp(s[0], ":Expression")) + else if(!strcmp(s[0], ":Expression")) { pExp=pDO->ParseExpression((const char**)(s),Len,pp->GetLineNumber(),pp->GetFilename()); if (pExp!=NULL){ @@ -273,9 +273,9 @@ bool ParseManagementFile(CModel *&pModel,const optStruct &Options) WriteWarning(warn.c_str(),Options.noisy); } } - else if (!strcmp(s[0], ":Condition")) + else if (!strcmp(s[0], ":Condition")) { - //TODO: handle more complex, expression-based conditions + //TODO: handle more complex, expression-based conditions if (pConst!=NULL){ dv_condition *pCond = new dv_condition(); pCond->dv_name=s[1]; @@ -287,11 +287,11 @@ bool ParseManagementFile(CModel *&pModel,const optStruct &Options) pCond->value2 = s_to_d(s[4]); } pConst->AddCondition(pCond); - } + } else{ ExitGracefully("ParseManagementFile: :Condition statement must appear after valid :Expression in :ManagementConstraint command",BAD_DATA_WARN); } - } + } else if (!strcmp(s[0], ":Penalty")) { pConst->penalty_under = s_to_d(s[1]); pConst->penalty_over = s_to_d(s[1]); @@ -300,14 +300,14 @@ bool ParseManagementFile(CModel *&pModel,const optStruct &Options) } } else if (!strcmp(s[0], ":Priority")) { - pConst->priority = s_to_i(s[1]); //for later + pConst->priority = s_to_i(s[1]); //for later } else if (!strcmp(s[0], ":EndManagementGoal")) { break; } else if (!strcmp(s[0], ":EndManagementConstraint")) { break; - } + } else { WriteWarning("ParseManagementFile: Unrecognized command in :ManagementConstraint command block",Options.noisy); } @@ -319,16 +319,16 @@ bool ParseManagementFile(CModel *&pModel,const optStruct &Options) break; } case(30): //---------------------------------------------- - {/*:LookupTable [name] - N - x1 y1 + {/*:LookupTable [name] + N + x1 y1 x2 y2 ... - xN yN + xN yN :EndLookupTable */ if(Options.noisy) { cout <<"LookupTable"<Tokenize(s,Len); @@ -336,13 +336,13 @@ bool ParseManagementFile(CModel *&pModel,const optStruct &Options) double *aX = new double[N]; double *aY = new double[N]; - for(int i = 0; i < N; i++) + for(int i = 0; i < N; i++) { pp->Tokenize(s,Len); if(IsComment(s[0],Len)) { i--; } else { if(Len>=2) { - aX[i] = s_to_d(s[0]); + aX[i] = s_to_d(s[0]); aY[i] = s_to_d(s[1]); } else { @@ -407,4 +407,4 @@ bool ParseManagementFile(CModel *&pModel,const optStruct &Options) pp=NULL; return true; -} \ No newline at end of file +} diff --git a/src/Raven.vcxproj b/src/Raven.vcxproj index 71269e09..54b19ce0 100644 --- a/src/Raven.vcxproj +++ b/src/Raven.vcxproj @@ -1030,4 +1030,4 @@ - \ No newline at end of file + diff --git a/src/Raven.vcxproj.filters b/src/Raven.vcxproj.filters index bd9e5bbc..4c04646f 100644 --- a/src/Raven.vcxproj.filters +++ b/src/Raven.vcxproj.filters @@ -688,4 +688,4 @@ Resource Files - \ No newline at end of file + diff --git a/src/Raven_BMI.cpp b/src/Raven_BMI.cpp index 9b0bdefe..db868cf6 100644 --- a/src/Raven_BMI.cpp +++ b/src/Raven_BMI.cpp @@ -10,13 +10,28 @@ const int GRID_SUBBASIN=1; const int GRID_HRU =0; +// constants used in the config file (cfg_file) parsing +const std::string CFG_FILE_CLIARGS_TAG = "cli_args"; +const std::string CFG_FILE_INPVARS_TAG = "input_vars"; +const std::string CFG_FILE_OUTVARS_TAG = "output_vars"; +const std::string CFG_FILE_IOVAR_FORCE_TAG = "forcing"; +const std::string CFG_FILE_IOVAR_HRUST_TAG = "hru_state"; +const std::string CFG_FILE_IOVAR_SUBST_TAG = "subbasin_state"; + ////////////////////////////////////////////////////////////////// /// \brief RavenBMI class constructor and destructor // CRavenBMI::CRavenBMI() { - //not sure how this will work with global variable defined in RavenMain. pModel=NULL; + + // vectors storing input/output variable names and ids are filled in the Initialize function + input_var_names = std::vector(); + input_var_ids = std::vector(); + output_var_names= std::vector(); + output_var_ids = std::vector(); + output_var_layer_index = std::vector(); + output_var_type = std::vector(); } CRavenBMI::~CRavenBMI() {} @@ -27,7 +42,7 @@ CRavenBMI::~CRavenBMI() {} /// /// \param line [in] string to be split /// \return vector of char* containing the split string -std::vector CRavenBMI::SplitLine(std::string line) +std::vector CRavenBMI::_SplitLineByWhitespace(std::string line) { std::vector args; std::istringstream iss(line); @@ -46,38 +61,47 @@ std::vector CRavenBMI::SplitLine(std::string line) ////////////////////////////////////////////////////////////////// -/// @brief Sets values to Option struct based on command line arguments +/// \brief Splits a string by colon into a vector of char* /// -/// @param config_key Argument key of the yaml config file -/// @param config_value Value of the argument -/// @return void - sets values in Options struct -void CRavenBMI::ProcessConfigFileArgument(std::string config_key, std::string config_value) +/// \param line [in] string to be split +/// \return vector of char* containing the split string +std::vector CRavenBMI::_SplitLineByColon(std::string line) { - if (config_key == "rvi_file") { - Options.rvi_filename = config_value; - } else if (config_key == "rvp_file") { - Options.rvp_filename = config_value; - } else if (config_key == "rvh_file") { - Options.rvh_filename = config_value; - } else if (config_key == "rvt_file") { - Options.rvt_filename = config_value; - } else if (config_key == "rvc_file") { - Options.rvc_filename = config_value; - } else if (config_key == "rvg_file") { - Options.rvg_filename = config_value; - } else if (config_key == "rve_file") { - Options.rve_filename = config_value; - } else if (config_key == "rvl_file") { - Options.rvl_filename = config_value; - } else if (config_key == "output_directory") { - Options.output_dir = config_value; - } else if (config_key == "main_output_directory") { - Options.main_output_dir = config_value; - } else { - // TODO: include all remaining options - cout << "WARNING: Unknown key '" << config_key << "' in Raven config file. Ignoring it." << endl; + int config_str_ini, config_str_mid; + std::vector args; + + config_str_ini = 0; + config_str_mid = (int)line.find(":", config_str_ini); + + // simple case: no colon found - return a list of one element (the whole string) + if (config_str_mid == -1) { + args.push_back(line); + return args; } - return; + + // first string is the key, add it to the list + args.push_back(line.substr(config_str_ini, config_str_mid - config_str_ini)); + + // even if there is nothing after the colon, add an empty string to the list to communicate that + // the line HAD a final colon + args.push_back(line.substr(config_str_mid + 1)); + + return args; +} + + +////////////////////////////////////////////////////////////////// +/// \brief Checks if a subbasin variable is valid +/// +/// \param var_name [in] name of the variable +/// \return true if the variable is valid, false otherwise +bool CRavenBMI::_IsValidSubBasinStateVariable(std::string var_name) +{ + if ((var_name == "OUTFLOW") || (var_name == "STREAMFLOW") || + (var_name == "RESERVOIR_STAGE")) { + return true; + } + return false; } @@ -86,7 +110,7 @@ void CRavenBMI::ProcessConfigFileArgument(std::string config_key, std::string co /// /// \param config_file [in] name of configuration file /// \return void - sets values in Options struct -void CRavenBMI::ReadConfigFile(std::string config_file) +void CRavenBMI::_ReadConfigFile(std::string config_file) { ifstream CONFIG; int config_str_ini, config_str_end; // used to parse the lines of the config file @@ -94,14 +118,18 @@ void CRavenBMI::ReadConfigFile(std::string config_file) bool cli_args = false; // used to check if command line arguments were used std::vector args; // used to parse the command line arguments char** argv; + int config_value_id; + std::vector line_split_by_colon; + bool listing_inp_vars = false; // flag to check if the config file is listing input variables + bool listing_out_vars = false; // flag to check if the config file is listing output variables CONFIG.open(config_file); if (CONFIG.fail()) { - cout << "Cannot find configuration file " << config_file << endl; + throw std::logic_error("Cannot find configuration file " + config_file); return; } - // read line by line and parse + // read and parse line by line for( std::string line; getline( CONFIG, line ); ) { // skip blank lines and comments @@ -109,43 +137,120 @@ void CRavenBMI::ReadConfigFile(std::string config_file) continue; } - // parse line and update options struct accordingly - config_str_end = -1; - config_key = ""; - config_value = ""; - do { - // find ini and end of substring - config_str_ini = config_str_end + 1; - config_str_end = (int)line.find(":", config_str_ini); - - // extract substring and trim leading and trailing whitespaces - config_value = line.substr(config_str_ini, config_str_end - config_str_ini); - config_value.erase(0, config_value.find_first_not_of(" \t")); - config_value.erase(config_value.find_last_not_of(" \t") + 1); - - // set config_key if it is the first substring - if (config_key == "") { - config_key = config_value; - config_value = ""; + // all lines of the config file must have one colon + line_split_by_colon = _SplitLineByColon(line); + if (line_split_by_colon.size() != 2) { + throw std::logic_error("WARNING: Invalid line in Raven config file: '" + line + "'"); + } + + // just to make the code more readable + config_key = line_split_by_colon[0]; + config_value = line_split_by_colon[1]; + + // remove leading and trailing whitespaces + config_key.erase(0, config_key.find_first_not_of(" \t")); + config_key.erase(config_key.find_last_not_of(" \t") + 1); + config_value.erase(0, config_value.find_first_not_of(" \t")); + config_value.erase(config_value.find_last_not_of(" \t") + 1); + + if (config_key == CFG_FILE_CLIARGS_TAG) { + + cli_args = true; + // "Raven.exe" is a dummy argument to mimic the command line arguments + args = _SplitLineByWhitespace("Raven.exe " + config_value); + argv = args.data(); + ProcessExecutableArguments((int)(args.size())-1, argv, Options); + listing_inp_vars = false; + listing_out_vars = false; + continue; + + } else if (config_key == CFG_FILE_INPVARS_TAG) { + + listing_inp_vars = true; + listing_out_vars = false; + continue; + + } else if (config_key == CFG_FILE_OUTVARS_TAG) { + + listing_out_vars = true; + listing_inp_vars = false; + continue; + + } else if (listing_inp_vars) { + + // all lines after the "input_vars:" line must be in the form "- [VAR_NAME]: [var_type]" + if (config_key.substr(0, 2) != "- ") { + throw std::logic_error("WARNING: Invalid input var line in Raven config file: " + line); + } + + // remove the '- ' from the beginning of the line + config_key = config_key.substr(2); + + if (config_value == CFG_FILE_IOVAR_FORCE_TAG) { + // case: regular - forcing as input + config_value_id = GetForcingTypeFromString(config_key); + if (config_value_id == F_UNRECOGNIZED) { + throw std::logic_error("WARNING: Invalid forcing type in Raven config file: '" + line + "'"); continue; } + input_var_names.push_back(config_key); + input_var_ids.push_back(config_value_id); + continue; + } else if (config_value == CFG_FILE_IOVAR_HRUST_TAG) { + // case: data assimilation - HRU state as input + // TODO: implement + throw std::logic_error("WARNING: HRU state variables as input are not supported in the BMI interface yet."); + continue; + } else if (config_value == CFG_FILE_IOVAR_SUBST_TAG) { + // case: data assimilation - sub-basin state as input + // TODO: implement + throw std::logic_error("WARNING: Sub-basin state variables as input are not supported in the BMI interface yet."); + continue; + } + throw std::logic_error("WARNING: Invalid input var line in Raven config file: " + line); + continue; + + } else if (listing_out_vars) { + + // all lines after the "input_vars:" line must be in the form "- [VAR_NAME]: [var_type]" + if (config_key.substr(0, 2) != "- ") { + throw std::logic_error("WARNING: Invalid input var line in Raven config file: " + line ); + continue; + } + + // remove the '- ' from the beginning of the line + config_key = config_key.substr(2); - // set options struct accordingly - if (config_key == "cli_args") { - cli_args = true; - // "Raven.exe" is a dummy argument to mimic the command line arguments - args = SplitLine("Raven.exe " + config_value); - argv = args.data(); - ProcessExecutableArguments((int)(args.size())-1, argv, Options); - break; + if (config_value == CFG_FILE_IOVAR_HRUST_TAG) { + // case: regular - state as output (evaluated after the model initialized) + output_var_names.push_back(config_key); + output_var_type.push_back(CFG_FILE_IOVAR_HRUST_TAG); + continue; + + } else if (config_value == CFG_FILE_IOVAR_SUBST_TAG) { + // case: regular - special cases as output + if (_IsValidSubBasinStateVariable(config_key)) { + output_var_names.push_back(config_key); + output_var_type.push_back(CFG_FILE_IOVAR_SUBST_TAG); } else { - ProcessConfigFileArgument(config_key, config_value); + throw std::logic_error("WARNING: Invalid sub-basin state variable in Raven config file: '" + line + "'"); } + continue; - } while (config_str_end != -1); + } else if (config_value == CFG_FILE_IOVAR_FORCE_TAG) { + // case: forcing as output? does it make sense? + throw std::logic_error("WARNING: Forcing variables as output are not supported in Raven's BMI interface."); + continue; - // if a "cli_args:" line was found then we can stop parsing - if (cli_args) { break; } + } + + // if we got here, the line is invalid + throw std::logic_error("WARNING: Invalid output var line in Raven config file: " + line); + continue; + } + + // if we got here, the line is invalid + throw std::logic_error("WARNING: Invalid line in Raven config file: " + line); } CONFIG.close(); @@ -163,7 +268,7 @@ void CRavenBMI::Initialize(std::string config_file) { //NOTE: ENSEMBLE MODE NOT SUPPORTED WITH BMI - ReadConfigFile(config_file); + _ReadConfigFile(config_file); PrepareOutputdirectory(Options); @@ -177,8 +282,7 @@ void CRavenBMI::Initialize(std::string config_file) WARNINGS.close(); Options.in_bmi_mode = true; // flag needed to ignore some arguments of the .rvi file - Options.rvt_filename = ""; // just a dummy filename to avoid errors - Options.duration = ALMOST_INF; // "infinity": will run as long as "Update()" is called + Options.rvt_filename = ""; // just a dummy filename to avoid errors //Read input files, create model, set model options if (!ParseInputFiles(pModel, Options)){ @@ -192,8 +296,33 @@ void CRavenBMI::Initialize(std::string config_file) pModel->SummarizeToScreen (Options); pModel->GetEnsemble()->Initialize (pModel,Options); + // pModel->Initialize() sets the duration to 365 (days) by default, but in BMI mode the duration + // of the simulation is governed by the framwork consuming the BMI interface, so we set the + // duration to ALMOST_INF(inite) to avoid the model stopping by itself + Options.duration = ALMOST_INF; + CheckForErrorWarnings(false, pModel); + // all the output variables must be checked + int config_out_var_layer_index; + int config_value_id; + for (int i = 0; i < output_var_names.size(); i++) { + if (output_var_type[i] == CFG_FILE_IOVAR_SUBST_TAG) { + output_var_ids.push_back(-1); + output_var_layer_index.push_back(-1); + continue; + } + if (output_var_type[i] == CFG_FILE_IOVAR_HRUST_TAG) { + config_value_id = pModel->GetStateVarInfo()->StringToSVType(output_var_names[i], + config_out_var_layer_index, + true); + output_var_ids.push_back(config_value_id); + output_var_layer_index.push_back(config_out_var_layer_index); + continue; + } + throw std::logic_error("WARNING: Output variable '" + output_var_names[i] + "' has an invalid type '" + output_var_type[i] + "'."); + } + PrepareOutputdirectory(Options); //adds new output folders, if needed pModel->WriteOutputFileHeaders(Options); @@ -206,8 +335,8 @@ void CRavenBMI::Initialize(std::string config_file) pModel->UpdateHRUForcingFunctions (Options,tt); pModel->UpdateDiagnostics (Options,tt); pModel->WriteMinorOutput (Options,tt); - } + ////////////////////////////////////////////////////////////////// /// \brief run simulation for a single time step /// @@ -238,6 +367,7 @@ void CRavenBMI::Update() //if ((Options.use_stopfile) && (CheckForStopfile(step,tt))) { break; } //step++; } + ////////////////////////////////////////////////////////////////// /// \brief run simulation for multiple time steps, ending at time 'time' /// @@ -250,6 +380,7 @@ void CRavenBMI::UpdateUntil(double time) Update(); } } + ////////////////////////////////////////////////////////////////// /// \brief called after simulation is done // @@ -275,44 +406,41 @@ std::string CRavenBMI::GetComponentName() { return "Raven Hydrological Modelling Framework "+Options.version; } + ////////////////////////////////////////////////////////////////// /// \brief returns number of accessible input datasets /// \return number of accessible input datasets // int CRavenBMI::GetInputItemCount() { - return 2; + return (input_var_names.size()); } + ////////////////////////////////////////////////////////////////// /// \brief returns array of accessible input dataset names /// \return array of accessible input dataset names // std::vector CRavenBMI::GetInputVarNames() { - vector names; - names.push_back("precipitation"); - names.push_back("temp_ave"); - return names; + return(input_var_names); } + ////////////////////////////////////////////////////////////////// /// \brief returns number of accessible output datasets /// \return number of accessible output datasets // int CRavenBMI::GetOutputItemCount() { - return 3; + return (output_var_names.size()); } + ////////////////////////////////////////////////////////////////// /// \brief returns array of accessible output dataset names /// \return array of accessible output dataset names // std::vector CRavenBMI::GetOutputVarNames() { - vector names; - names.push_back("streamflow"); - names.push_back("soil[0]"); - names.push_back("snow"); - return names; + return(output_var_names); } //------------------------------------------------------------------ @@ -325,17 +453,24 @@ std::vector CRavenBMI::GetOutputVarNames() // int CRavenBMI::GetVarGrid(std::string name) { - //input variables - if (name=="precipitation"){return GRID_HRU;} - else if (name=="temp_ave" ){return GRID_HRU;} - //output variables - if (name=="streamflow" ){return GRID_SUBBASIN;} - else if (name=="soil[0]" ){return GRID_HRU;} - else if (name=="snow" ){return GRID_HRU;} + // subbasin states, surprisingly, follow the subbasin grid + if (_IsValidSubBasinStateVariable(name)) { + return(GRID_SUBBASIN); + } - return 0; + // forcings are applies to each HRU, so they follow the HRU grid + forcing_type Ftype = GetForcingTypeFromString(name); + if (Ftype != F_UNRECOGNIZED) { + return(GRID_HRU); // is this a good assumption? + } + + // if not a subbasin state nor a forcing, it MUST be a hru state var (or the code will just break) + int state_var_layer_idx; + sv_type SType = pModel->GetStateVarInfo()->StringToSVType(name, state_var_layer_idx, true); + return(GRID_HRU); } + ////////////////////////////////////////////////////////////////// /// \brief returns units of input or output variable as string /// \param name [in] - name of variable @@ -343,17 +478,29 @@ int CRavenBMI::GetVarGrid(std::string name) // std::string CRavenBMI::GetVarUnits(std::string name) { - //input variables - if (name=="precipitation"){return "mm/d";} - else if (name=="temp_ave" ){return "C";} + int state_var_layer_idx; + forcing_type Ftype; + sv_type SType; + + // "OUTFLOW" is just another name for "STREAMFLOW", but "OUTFLOW" does not + // have a sv_type, so it needs to be treated as a special case + if (name == "OUTFLOW") { + SType = pModel->GetStateVarInfo()->StringToSVType("STREAMFLOW", state_var_layer_idx, + true); + return(CStateVariable::GetStateVarUnits(SType)); + } - //output variables - if (name=="streamflow" ){return "m3/s";} - else if (name=="soil[0]" ){return "mm"; } - else if (name=="snow" ){return "mm"; } + // the variable can be a forcing, so we need to search for it + Ftype = GetForcingTypeFromString(name); + if (Ftype != F_UNRECOGNIZED) { + return(GetForcingTypeUnits(Ftype)); + } - return ""; + // if not a special case nor a forcing, it MUST be a state var (or the code will just break) + SType = pModel->GetStateVarInfo()->StringToSVType(name, state_var_layer_idx, true); + return(CStateVariable::GetStateVarUnits(SType)); } + ////////////////////////////////////////////////////////////////// /// \brief returns type of input or output variable as string /// \param name [in] - name of variable @@ -363,6 +510,7 @@ std::string CRavenBMI::GetVarType(std::string name) { return "double"; } + ////////////////////////////////////////////////////////////////// /// \brief returns size of input or output variable, in bytes /// \param name [in] - name of variable @@ -372,6 +520,7 @@ int CRavenBMI::GetVarItemsize(std::string name) { return sizeof(double); } + ////////////////////////////////////////////////////////////////// /// \brief returns total size of input or output variable array, in bytes /// \param name [in] - name of variable @@ -381,6 +530,7 @@ int CRavenBMI::GetVarNbytes(std::string name) { return GetVarItemsize(name)*GetGridSize(GetVarGrid(name)); } + ////////////////////////////////////////////////////////////////// /// \brief returns variable location on unstructured grid /// \param name [in] - name of variable @@ -390,6 +540,7 @@ std::string CRavenBMI::GetVarLocation(std::string name) { return "node"; } + ////////////////////////////////////////////////////////////////// /// \brief accessors for model basic elements // @@ -414,10 +565,11 @@ void CRavenBMI::GetValue(std::string name, void* dest) int k,p; bool is_HRU_SV=false; int iSV=DOESNT_EXIST; + sv_type Stype= (sv_type)DOESNT_EXIST; + int Slayer=DOESNT_EXIST; - //output variables - if (name=="streamflow") - { + if ((name == "OUTFLOW") || (name =="STREAMFLOW")) { + // special case: outflow = streamflow (subbasin state var) out=new double [pModel->GetNumSubBasins()]; for (p = 0; p < pModel->GetNumSubBasins(); p++) { @@ -430,13 +582,46 @@ void CRavenBMI::GetValue(std::string name, void* dest) } memcpy(dest,out,pModel->GetNumSubBasins()*sizeof(double)); delete [] out; + return; } - else if (name == "soil[0]") {is_HRU_SV=true; iSV=pModel->GetStateVarIndex(SOIL,0);} - else if (name == "snow") {is_HRU_SV=true; iSV=pModel->GetStateVarIndex(SNOW); } - else { - throw std::logic_error("RavenBMI.GetValue: variable '" + name + "' not covered by this function."); + else if (name == "RESERVOIR_STAGE") { + // special case: reservoir stage is a subbasin state var + out=new double [pModel->GetNumSubBasins()]; + for (p = 0; p < pModel->GetNumSubBasins(); p++) { + CSubBasin* pBasin = pModel->GetSubBasin(p); + if (pBasin->GetReservoir() != NULL) { + out[p] = pBasin->GetReservoir()->GetResStage(); + } else { + out[p] = 0.0; // or maybe '-1' to indicate that there is no reservoir? + } + } + memcpy(dest,out,pModel->GetNumSubBasins()*sizeof(double)); + delete [] out; return; } + else { + // search in the output variable names + for (int i = 0; i < output_var_names.size(); i++) { + if (output_var_names[i] != name) { + continue; + } + if (output_var_type[i] != CFG_FILE_IOVAR_HRUST_TAG) { + throw std::logic_error("CRavenBMI.GetValue: variable '" + name + "' is not an HRU state variable."); + return; + } + + Stype = (sv_type)output_var_ids[i]; + Slayer = output_var_layer_index[i]; + is_HRU_SV=true; + if (Slayer > -1) { + iSV = pModel->GetStateVarIndex(Stype, Slayer); + } else { + iSV = pModel->GetStateVarIndex(Stype); + } + + break; + } + } if ((is_HRU_SV) && (iSV!=DOESNT_EXIST)) { @@ -446,9 +631,14 @@ void CRavenBMI::GetValue(std::string name, void* dest) } memcpy(dest,out,pModel->GetNumSubBasins()*sizeof(double)); delete [] out; + return; } + // if it got here... dude, I have a bad news for you... + throw std::logic_error("RavenBMI.GetValue: variable '" + name + "' not covered by this function."); + return; } + ////////////////////////////////////////////////////////////////// /// \brief returns array of variable values for variable with supplied name at subset o flocations /// \param name [in] - name of variable @@ -463,16 +653,56 @@ void CRavenBMI::GetValueAtIndices(std::string name, void* dest, int* inds, int c bool is_HRU_SV=false; int iSV=DOESNT_EXIST; int i,k,p; + sv_type Stype=(sv_type)DOESNT_EXIST; + int Slayer=DOESNT_EXIST; - if (name=="streamflow") - { + if ((name == "OUTFLOW") || (name =="STREAMFLOW")) { + // special case: outflow = streamflow (subbasin state var) + for (i = 0; i GetSubBasin(p)->GetIntegratedOutflow(Options.timestep)/(Options.timestep*SEC_PER_DAY); + } + else{ + out[p]=pModel->GetSubBasin(p)->GetOutflowRate(); + } + } + } + else if (name == "RESERVOIR_STAGE") { + // special case: reservoir stage is a subbasin state var for (i = 0; i GetSubBasin(p)->GetIntegratedOutflow(Options.timestep); + CSubBasin* pBasin = pModel->GetSubBasin(p); + if (pBasin->GetReservoir() != NULL) { + out[p] = pBasin->GetReservoir()->GetResStage(); + } else { + out[p] = 0.0; // or maybe '-1' to indicate that there is no reservoir? + } + } + } + else { + // search in the output variable names + for (int i = 0; i < output_var_names.size(); i++) { + if (output_var_names[i] != name) { + continue; + } + if (output_var_type[i] != CFG_FILE_IOVAR_HRUST_TAG) { + throw std::logic_error("CRavenBMI.GetValue: variable '" + name + "' is not an HRU state variable."); + return; + } + + Stype = (sv_type)output_var_ids[i]; + Slayer = output_var_layer_index[i]; + is_HRU_SV=true; + if (Slayer > -1) { + iSV = pModel->GetStateVarIndex(Stype, Slayer); + } else { + iSV = pModel->GetStateVarIndex(Stype); + } + + break; } } - else if (name == "soil[0]") {is_HRU_SV=true; iSV=pModel->GetStateVarIndex(SOIL,0);} - else if (name == "snow") {is_HRU_SV=true; iSV=pModel->GetStateVarIndex(SNOW); } if ((is_HRU_SV) && (iSV!=DOESNT_EXIST)) { @@ -517,9 +747,14 @@ void CRavenBMI::SetValue(std::string name, void* src) forcing_type Ftype; bool is_forcing=false; - if (name=="precipitation"){is_forcing=true; Ftype=F_PRECIP; } - else if (name=="temp_ave") {is_forcing=true; Ftype=F_TEMP_AVE;} + Ftype = GetForcingTypeFromString(name); + if (Ftype != F_UNRECOGNIZED) { + is_forcing = true; + } else { + throw std::logic_error("CRavenBMI.SetValue: only accepts setting forcing values (for now)."); + } + // all HRUs receive the same forcing if (is_forcing){ for (int k=0;kGetNumHRUs();k++){ pModel->GetHydroUnit(k)->SetHRUForcing(Ftype, input[k]); @@ -540,9 +775,14 @@ void CRavenBMI::SetValueAtIndices(std::string name, int* inds, int count, void* forcing_type Ftype; bool is_forcing=false; - if (name=="precipitation"){is_forcing=true; Ftype=F_PRECIP;} - else if (name=="temp_ave") {is_forcing=true; Ftype=F_TEMP_AVE;} + Ftype = GetForcingTypeFromString(name); + if (Ftype != F_UNRECOGNIZED) { + is_forcing = true; + } else { + throw std::logic_error("WARNING: CRavenBMI.SetValueAtIndices: only accepts forcings now. Tried to set value for variable '"+name+"'"); + } + // all HRUs receive the same forcing if (is_forcing){ for (int i=0;i input_var_names; + std::vector input_var_ids; + std::vector output_var_names; + std::vector output_var_ids; + std::vector output_var_layer_index; + std::vector output_var_type; + // Internal functions for reading the YAML config file. - void ReadConfigFile(std::string config_file); - std::vector SplitLine(std::string line); - void ProcessConfigFileArgument(std::string config_key, std::string config_value); + void _ReadConfigFile(std::string config_file); + std::vector _SplitLineByWhitespace(std::string line); + std::vector _SplitLineByColon(std::string line); + bool _IsValidSubBasinStateVariable(std::string var_name); public: CRavenBMI(); diff --git a/src/Reservoir.cpp b/src/Reservoir.cpp index 717a562f..b3227cc1 100644 --- a/src/Reservoir.cpp +++ b/src/Reservoir.cpp @@ -401,14 +401,14 @@ double CReservoir::GetOldStage () const { return _stage_last; } ////////////////////////////////////////////////////////////////// /// \returns min stage [m] // -double CReservoir::GetMinStage (const int nn) const { +double CReservoir::GetMinStage (const int nn) const { if(_pMinStageTS !=NULL){ return _pMinStageTS->GetSampledValue(nn);} return _min_stage; } ////////////////////////////////////////////////////////////////// /// \returns max stage [m] // -double CReservoir::GetMaxStage (const int nn) const { +double CReservoir::GetMaxStage (const int nn) const { if(_pMaxStageTS !=NULL){ return _pMaxStageTS->GetSampledValue(nn);} return ALMOST_INF; } diff --git a/src/Solvers.cpp b/src/Solvers.cpp index 7de5ae36..f8ae54bb 100644 --- a/src/Solvers.cpp +++ b/src/Solvers.cpp @@ -541,11 +541,11 @@ void MassEnergyBalance( CModel *pModel, pBasin->UpdateSubBasin(tt,Options); } - // Management optimization - determines optimal demand delivery/reservoir outflows + // Management optimization - determines optimal demand delivery/reservoir outflows // ---------------------------------------------------------------------------------------- if (Options.management_optimization) - { - pModel->GetDemandOptimizer()->SolveDemandProblem(pModel, Options, aRouted, tt); + { + pModel->GetDemandOptimizer()->SolveDemandProblem(pModel, Options, aRouted, tt); } // Route water over timestep diff --git a/src/StandardOutput.cpp b/src/StandardOutput.cpp index f94c3d9c..057a06ed 100644 --- a/src/StandardOutput.cpp +++ b/src/StandardOutput.cpp @@ -140,7 +140,7 @@ void CModel::CloseOutputStreams() if ( _LEVELS.is_open()){ _LEVELS.close();} if (_pDO!=NULL){_pDO->CloseOutputStreams();} - + #ifdef _RVNETCDF_ int retval; // error value for NetCDF routines diff --git a/src/SubBasin.h b/src/SubBasin.h index 13351b21..c7492478 100644 --- a/src/SubBasin.h +++ b/src/SubBasin.h @@ -240,7 +240,7 @@ class CSubBasin double GetDownstreamIrrDemand (const double &t) const; //[m3/s] cumulative downstream irrigation demand, including from this subbasin double GetDemandDelivery () const; //[m3/s] instantaneous delivery rate Qirr double GetEnviroMinFlow (const double &t) const; //[m3/s] environmental minimum flow target from downstream outlet - bool HasEnviroMinFlow () const; // true if basin has enviro min flow time series + bool HasEnviroMinFlow () const; // true if basin has enviro min flow time series bool HasIrrigationDemand () const; // true if basin has specified irrigation demand int GetDiversionTargetIndex (const int i) const; // returns subbasin index p of diversion i