Pull request #218: Plp/bugfix/060121

Merge in JGCRI/gcam-core from plp/bugfix/060121 to master

* commit '64ab9a0040a16aab53a3a059f84ba54d3b92db82': (31 commits)
  Update links to slides and videos
  Add GCAMv5.1 publication to README
  Add changes that attempts to detect if Broyden is stuck in a local minimum and allows it to try to jump out if so.  Otherwise it can not becuase we only allow each step to strictly make the solution better.
  Bump version to 5.4
  One last rebuild of package data and fix some errors / warnings after CEDS update
  Remove TODO in comments now that we added SolutionInfo::getSolutionFloor
  Address comments on zchunk_L112.ceds_ghg_en_R_S_T_Y.R
  Address PR comments including initialization of params in GHGPolicy and use solution floor instead of arbitrary hard coding it in the Preconditioner.  To do this we have to make the solution floor available in SolutionInfo.
  Fix bug where fgas emissions unit only got set from some years and not all (applicable) years resulting in inconsistent units reported to the database by year
  Update zchunk_L112.ceds_ghg_en_R_S_T_Y.R
  Update initial tax guesses for 2.6 target finder files for 5.4
  Fix primary energy rewrites for traded*CCS
  Clean up a couple more precursor warnings
  Some fixes to primary energy and food demand queries
  Take case of deprecation warning about `group_by` / `add`
  Update the compiled .jar files and the modelinterface submodule to reflect the bugfix changes
  Some solution changes that got lost in the shuffle which were meant to be in with the last set of solution fixes (the CMP even describes that it is!).  The only change of consequence is to set the scale of demands to 1 if the market is solved below the solution floor so that Broyden doesn't waste iterations on such markets.  In theory this should help with variability in solution performance with parallel turned on.
  Fix bug with negative emissions budget + linked policies + constraint which was causing GHG constraints to have poor solution behavior
  Fix GHGPolicy so that taxes and constraints do not use the sentinal value -1 to indicate not set and use the Value class in this case.
  Update electricity cost files / calculations to ensure that all technologies receive a variable O&M cost input (even if it is zero).  This adds an OM-var cost for CSP & CSP_storage, which will (modestly) impact results.
  ...

README.md

@@ -53,13 +53,12 @@ mitigation policy.
 * [GCAM Documentation](http://jgcri.github.io/gcam-doc/)
 * [Getting Started with GCAM](http://jgcri.github.io/gcam-doc/user-guide.html)
 * [GCAM Community](http://www.globalchange.umd.edu/models/gcam/gcam-community/)
-* [2019 Tutorial Slides](http://www.globalchange.umd.edu/data/gcam/GCAM_tutorial_2019.pdf)
-* GCAM lectures from the UMD Water and Climate class
-  * [Lecture 3: Introduction to GCAM](https://www.youtube.com/watch?v=xRF9lFwtMr0)
-  * [Lecture 4: GCAM Tutorial](https://www.youtube.com/watch?v=S7vAShH-dbs)
+* [GCAM Videos and Tutorial Slides](https://gcims.pnnl.gov/community)
 
 ## Selected Publications
 
+Calvin, K., Patel, P., Clarke, L., Asrar, G., Bond-Lamberty, B., Cui, R. Y., Di Vittorio, A., Dorheim, K., Edmonds, J., Hartin, C., Hejazi, M., Horowitz, R., Iyer, G., Kyle, P., Kim, S., Link, R., McJeon, H., Smith, S. J., Snyder, A., Waldhoff, S., and Wise, M.: GCAM v5.1: representing the linkages between energy, water, land, climate, and economic systems, Geosci. Model Dev., 12, 677–698, https://doi.org/10.5194/gmd-12-677-2019, 2019.
+
 Edmonds, J., and J. Reilly (1985)Global Energy: Assessing the Future (Oxford University Press, New York) pp.317.
 
 Edmonds, J., M. Wise, H. Pitcher, R. Richels, T. Wigley, and C. MacCracken. (1997) “An Integrated Assessment of Climate Change and the Accelerated Introduction of Advanced Energy Technologies”, Mitigation and Adaptation Strategies for Global Change, 1, pp. 311-39

cvs/objects/java/source/Makefile

@@ -1,6 +1,6 @@
 # Java compiler options to target specific JVM versions to allow
 # backwards compatibility.
-JAVA_TARGET = -target 1.7 -source 1.7
+JAVA_TARGET = -target 1.8 -source 1.8
 
 install: XMLDBDriver.jar
 	cp XMLDBDriver.jar ../../../../exe

cvs/objects/java/source/WildcardExpandingClassLoader.java

@@ -62,7 +62,7 @@ public class WildcardExpandingClassLoader extends URLClassLoader {
      * A pattern that looks for a "prefered" BaseX version if we find
      * duplicate versions of the library.
      */
-    final static String PREFERED_BASEX_VER = "^.*[Bb][Aa][Ss][Ee][Xx]-9.5.0.*$";
+    final static String PREFERED_BASEX_VER = "^.*[Bb][Aa][Ss][Ee][Xx]-9.5.[0-9].*$";
 
     /**
      * Constructor which will expand any wildcard jar specifications then call

cvs/objects/policy/include/policy_ghg.h

@@ -52,6 +52,7 @@
 
 #include "util/base/include/inamed.h"
 #include "util/base/include/time_vector.h"
+#include "util/base/include/value.h"
 #include "util/base/include/data_definition_util.h"
 
 // Need to forward declare the subclasses as well.
@@ -94,10 +95,10 @@ class GHGPolicy: public INamed, private boost::noncopyable {
         DEFINE_VARIABLE( SIMPLE, "market", mMarket, std::string ),
 
         //! Emissions constraint by year(tgC or MTC)
-        DEFINE_VARIABLE( ARRAY, "constraint", mConstraint, objects::PeriodVector<double> ),
+        DEFINE_VARIABLE( ARRAY, "constraint", mConstraint, objects::PeriodVector<Value> ),
 
         //! Fixed tax on Emissions by year($/TC)
-        DEFINE_VARIABLE( ARRAY, "fixedTax", mFixedTax, objects::PeriodVector<double> )
+        DEFINE_VARIABLE( ARRAY, "fixedTax", mFixedTax, objects::PeriodVector<Value> )
     )
 
     void copy( const GHGPolicy& aOther );

cvs/objects/policy/source/linked_ghg_policy.cpp

@@ -218,7 +218,7 @@ void LinkedGHGPolicy::completeInit( const string& aRegionName ) {
     // Add a dependency between this mand the linked market and include a dummy
     // activity to ensure the dependecy finder has an activity to traverse
     MarketDependencyFinder* depFinder = marketplace->getDependencyFinder();
-    depFinder->addDependency( mPolicyName, aRegionName, mLinkedPolicyName, aRegionName );
+    depFinder->addDependency( mPolicyName, aRegionName, mLinkedPolicyName, aRegionName, false );
     depFinder->resolveActivityToDependency( aRegionName, mPolicyName,
             new DummyActivity(), new DummyActivity() );
 

cvs/objects/policy/source/policy_ghg.cpp

@@ -62,9 +62,7 @@ using namespace xercesc;
 extern Scenario* scenario;
 
 /*! \brief Default constructor. */
-GHGPolicy::GHGPolicy():
-mConstraint( -1 ),
-mFixedTax( -1 )
+GHGPolicy::GHGPolicy()
 {
 }
 
@@ -73,24 +71,21 @@ mFixedTax( -1 )
 *        constraint.
 */
 GHGPolicy::GHGPolicy( const string aName, const string aMarket ):
-mConstraint( -1 ),
-mFixedTax( -1 )
+mName(aName),
+mMarket(aMarket)
 {
-    mName = aName;
-    mMarket = aMarket;
 }
 
 /*! \brief Constructor used when explicitly constructing a fixed tax.
 */
 GHGPolicy::GHGPolicy( const string aName, const string aMarket,
-                      const vector<double>& aTaxes )
+                      const vector<double>& aTaxes ):
+mName(aName),
+mMarket(aMarket)
 {
     // Ensure that the taxes vector passed in is the right size.
     assert( aTaxes.size() == mConstraint.size() );
 
-    mName = aName;
-    mMarket = aMarket;
-    mConstraint.assign( mConstraint.size(), -1 );
     std::copy( aTaxes.begin(), aTaxes.end(), mFixedTax.begin() );
 }
 
@@ -165,10 +160,10 @@ void GHGPolicy::XMLParse( const DOMNode* node ){
             mMarket = XMLHelper<string>::getValue( curr ); // should be only one market
         }
         else if( nodeName == "constraint" ){
-            XMLHelper<double>::insertValueIntoVector( curr, mConstraint, modeltime );
+            XMLHelper<Value>::insertValueIntoVector( curr, mConstraint, modeltime );
         }
         else if( nodeName == "fixedTax" ){
-            XMLHelper<double>::insertValueIntoVector( curr, mFixedTax, modeltime );
+            XMLHelper<Value>::insertValueIntoVector( curr, mFixedTax, modeltime );
         }
         else {
             ILogger& mainLog = ILogger::getLogger( "main_log" );
@@ -219,9 +214,9 @@ void GHGPolicy::completeInit( const string& aRegionName ) {
 
     // check for missing periods in which case interpolate
     for( int i = 1; i < modeltime->getmaxper(); ++i ) {
-        if( mFixedTax[ i ] == -1 && mFixedTax[ i - 1 ] != -1 ) {
+        if( !mFixedTax[ i ].isInited() && mFixedTax[ i - 1 ].isInited() ) {
             int j;
-            for( j = i + 1; j < modeltime->getmaxper() && mFixedTax[ j ] == -1; ++j ) {
+            for( j = i + 1; j < modeltime->getmaxper() && !mFixedTax[ j ].isInited(); ++j ) {
             }
             if( j < modeltime->getmaxper() ) {
                 mFixedTax[ i ] = util::linearInterpolateY( modeltime->getper_to_yr( i ),
@@ -231,9 +226,9 @@ void GHGPolicy::completeInit( const string& aRegionName ) {
                                                            mFixedTax[ j ] );
             }
         }
-        if( mConstraint[ i ] == -1 && mConstraint[ i - 1 ] != -1 ) {
+        if( !mConstraint[ i ].isInited() && mConstraint[ i - 1 ].isInited() ) {
             int j;
-            for( j = i + 1; j < modeltime->getmaxper() && mConstraint[ j ] == -1; ++j ) {
+            for( j = i + 1; j < modeltime->getmaxper() && !mConstraint[ j ].isInited(); ++j ) {
             }
             if( j < modeltime->getmaxper() ) {
                 mConstraint[ i ] = util::linearInterpolateY( modeltime->getper_to_yr( i ),
@@ -250,11 +245,13 @@ void GHGPolicy::completeInit( const string& aRegionName ) {
     // If it is a constraint, add the constraint to the market and set the 
     // market to solve.
     for( unsigned int i = 0; i < modeltime->getmaxper(); ++i ){
-        if( mFixedTax[ i ] != -1 ){
+        if( mFixedTax[ i ].isInited() ){
             marketplace->unsetMarketToSolve( mName, aRegionName, i );
             marketplace->setPrice( mName, aRegionName, mFixedTax[ i ], i );
         }
-        else if( mConstraint[ i ] != -1 ){
+        // if both a fixed tax and constraint are provided the constraint takes
+        // precedence and the fixed tax is used as initial guess
+        if( mConstraint[ i ].isInited() ){
             marketplace->setMarketToSolve( mName, aRegionName, i );
             // Adding the difference between the constraint for this period
             // and the current supply because addToSupply adds to the current
@@ -265,7 +262,7 @@ void GHGPolicy::completeInit( const string& aRegionName ) {
 
         // GHG policies must have a price >= 0.  It may be the case that the constraint is
         // non-binding at a zero price in which case the solver can use this information to
-        // make a supply currection to still ensure equality.
+        // make a supply correction to still ensure equality.
         SectorUtils::setSupplyBehaviorBounds( mName, aRegionName, 0, util::getLargeNumber(), i );
     }
 }

cvs/objects/reporting/source/xml_db_outputter.cpp

@@ -123,7 +123,7 @@
 #include <string>
 #include <sstream>
 
-#include <boost/math/tr1.hpp>
+#include <cmath>
 
 #include "reporting/include/xml_db_outputter.h"
 
@@ -666,7 +666,7 @@ void XMLDBOutputter::startVisitSector( const Sector* aSector, const int aPeriod
     const Modeltime* modeltime = scenario->getModeltime();
     for( int i = 0; i < modeltime->getmaxper(); ++i ){
         double currCost = aSector->getPrice( mGDP, i );
-        if( !boost::math::isnan( currCost ) ) {
+        if( !std::isnan( currCost ) ) {
             writeItem( "cost", mCurrentPriceUnit, currCost, i );
         }
     }
@@ -713,7 +713,7 @@ void XMLDBOutputter::startVisitSubsector( const Subsector* aSubsector,
     }
     for( int i = 0; i < modeltime->getmaxper(); ++i ){
         double currValue = aSubsector->getPrice( mGDP, i );
-        if( !objects::isEqual<double>( currValue, 0.0 ) && !boost::math::isnan( currValue ) ) {
+        if( !objects::isEqual<double>( currValue, 0.0 ) && !std::isnan( currValue ) ) {
             writeItem( "cost", mCurrentPriceUnit, currValue, i );
         }
     }
@@ -832,7 +832,7 @@ void XMLDBOutputter::startVisitTechnology( const Technology* aTechnology, const
         // "year" attribute in addition to the technology "year" attribute.
         if( mCurrentTechnology->mProductionState[ curr ] && mCurrentTechnology->mProductionState[ curr ]->isNewInvestment() ){
             double currValue = aTechnology->getCost( curr );
-            if( !objects::isEqual<double>( currValue, 0.0 ) && !boost::math::isnan( currValue ) ) {
+            if( !objects::isEqual<double>( currValue, 0.0 ) && !std::isnan( currValue ) ) {
                 writeItemToBuffer( aTechnology->getCost( curr ), "cost", 
                     *childBuffer, mTabs.get(), curr, mCurrentPriceUnit );
             }

cvs/objects/sectors/source/subsector.cpp

@@ -46,7 +46,7 @@
 #include <algorithm>
 #include <xercesc/dom/DOMNode.hpp>
 #include <xercesc/dom/DOMNodeList.hpp>
-#include <boost/math/tr1.hpp>
+#include <cmath>
 
 #include "util/base/include/configuration.h"
 #include "sectors/include/subsector.h"
@@ -472,7 +472,7 @@ void Subsector::calcCost( const int aPeriod ){
 double Subsector::calcShare( const IDiscreteChoice* aChoiceFn, const GDP* aGDP, const int aPeriod ) const {
     double subsectorPrice = getPrice( aGDP, aPeriod );
 
-    if( boost::math::isnan( subsectorPrice ) ) {
+    if( std::isnan( subsectorPrice ) ) {
         // Check for a NaN sentinel value.  If we find it, set the
         // subsector's share to zero.
         return -numeric_limits<double>::infinity();

cvs/objects/solution/solvers/source/logbroyden.cpp

@@ -40,6 +40,7 @@
 
 #include "util/base/include/definitions.h"
 #include <string>
+#include <queue>
 #include <algorithm>
 #include <iomanip>
 #include <math.h>
@@ -413,6 +414,11 @@ int LogBroyden::bsolve(VecFVec &F, UBVECTOR &x, UBVECTOR &fx,
 
   // We calculate a scalar value for F(x) using F(x) * F(x)
   double f0 = fx.dot(fx); // already have a value of F on input, so no need to call fnorm yet
+    // Keep track of several past f(x) values which we will use to determine if
+    // progress has stalled out or not.
+    const int TRACK_NUM_PAST_F_VALUES = 4;
+    std::queue<double> past_f_values;
+    past_f_values.push(f0);
   if(f0 < FTINY) {
     // Guard against F=0 since it can cause a NaN in our solver.  This
     // is a more stringent test than our regular convergence test
@@ -535,7 +541,21 @@ int LogBroyden::bsolve(VecFVec &F, UBVECTOR &x, UBVECTOR &fx,
     // dx now holds the newton step.  Execute the line search along
     // that direction.
     double fnew;
-    int lserr = linesearch(F,x,f0,gx,dx, xnew,fnew, neval, &solverLog);
+      // if we are making adequate progress we only allow linesearch to accept
+      // steps that are make f(x) smaller
+      double fxIncr = 0.0;
+      solverLog << "Past f size: " << past_f_values.size();
+      solverLog << " values F: " << past_f_values.front() << ", B: " << past_f_values.back() << std::endl;
+      if(past_f_values.size() == TRACK_NUM_PAST_F_VALUES && past_f_values.back()*1.1 > past_f_values.front() ) {
+          solverLog << "Taking a chance to try to jump out of local min.\n";
+          // very little progress, might be stuck in a local minima
+          // let's take a chance and take a step out of our comfort zone
+          // by accepting a step that makes f(x) at most 1000 time worse
+          fxIncr = 1000.0;
+          past_f_values = std::queue<double>();
+      }
+      UBVECTOR fxnew(fx.size());
+    int lserr = linesearch(F,x,f0,gx,dx, xnew,fnew, fxnew, fxIncr, neval, &solverLog);
 
     if(lserr != 0) {
       // line search failed.  There are a couple of things that could
@@ -600,16 +620,17 @@ int LogBroyden::bsolve(VecFVec &F, UBVECTOR &x, UBVECTOR &fx,
       // reset the line search fail flag
       lsfail = false;
     }
+      // keep track of the last TRACK_NUM_PAST_F_VALUES f(x) values only
+      if(past_f_values.size() == TRACK_NUM_PAST_F_VALUES) {
+          past_f_values.pop();
+      }
+      past_f_values.push(fnew);
 
     UBVECTOR xstep(xnew-x);    // step in x eventually taken
     double lambda = fabs(dx[0]) > 0.0 ? xstep[0] / dx[0] : 0.0;
     solverLog << "################Return from linesearch\nfold= " << f0 << "\tfnew= " << fnew
               << "\tlambda= " << lambda << "\n";
 
-    UBVECTOR fxnew(fx.size());
-      // TODO: if we return the last fx from linesearch we wouldn't have to recalculate
-      // it here
-      F(xnew, fxnew);
     //solverLog << "\nxnew: " << xnew << "\nfxnew: " << fxnew << "\n";
     UBVECTOR fxstep(fxnew -fx); // change in F( x ).  We will need this for the secant update
 
@@ -652,9 +673,7 @@ int LogBroyden::bsolve(VecFVec &F, UBVECTOR &x, UBVECTOR &fx,
         fxstep -= B * xstep;
       fxstep /= dx2;
         B += fxstep * xstep.transpose();
-        if(iter >0) {
       ageB++;                // increment the age of B
-        }
     }
     else {
       // Progress using the Broyden formula is anemic.  This usually

cvs/objects/solution/solvers/source/preconditioner.cpp

@@ -221,7 +221,7 @@ SolverComponent::ReturnCode Preconditioner::solve( SolutionInfoSet& aSolutionSet
             double newprice = oldprice;
             bool chg = false;
             double lb,ub;       // only used for normal markets, but need to be declared up here.
-            bool isSolved = solvable[i].getRelativeED() < mFTOL || fabs(solvable[i].getED()) < mFTOL;
+            bool isSolved = solvable[i].getRelativeED() < mFTOL || fabs(solvable[i].getED()) < mFTOL || solvable[i].isSolved();
 
             if(pass > 0) {
                 // If this market is close to solved update the "forecast" price and demand which
@@ -243,10 +243,14 @@ SolverComponent::ReturnCode Preconditioner::solve( SolutionInfoSet& aSolutionSet
                     fp = fd = newScale;
                 }
                 else if(isSolved && !(solvable[i].getType() == IMarketType::TAX || solvable[i].getType() == IMarketType::SUBSIDY)) {
+                    double solutionFloor = solvable[i].getSolutionFloor();
+                    // Check if the market is solved to the solution floor in which case don't let
+                    // NR make a big deal out of relative differences by reseting the demand scale to 1
+                    double demandScale = abs(olddmnd) < solutionFloor || abs(oldsply) < solutionFloor ? 1.0 : olddmnd;
                     solvable[i].setForecastPrice(oldprice);
-                    solvable[i].setForecastDemand(olddmnd);
+                    solvable[i].setForecastDemand(demandScale);
                     fp = oldprice;
-                    fd = olddmnd;
+                    fd = demandScale;
                 }
                 if(fd == 0.0) {
                     if(olddmnd > 0.0) {
@@ -344,7 +348,7 @@ SolverComponent::ReturnCode Preconditioner::solve( SolutionInfoSet& aSolutionSet
                     lb = solvable[i].getLowerBoundSupplyPrice();
                     ub = solvable[i].getUpperBoundSupplyPrice();
                     if(!isSolved && oldprice < lb) {
-                        newprice = 0.001;
+                        newprice = lb + 0.001;
                         solvable[i].setPrice(newprice);
                         chg = true;
                         ++nchg;

cvs/objects/solution/util/include/linesearch.hpp

@@ -61,7 +61,9 @@
  * \param[in] g0: Value of grad-f(x0)
  * \param[in] dx: Proposed step from the root finder
  * \param[out] x: Final step determined by the search
- * \param[out]fx: Value of f(x) 
+ * \param[out] fx: Value of f(x)
+ * \param[in] fxVec: The value of F(x)
+ * \param[in] fxIncr: The maximum increase of f(x) we will allow to consider linesearch a success.
  * \param[inout]neval: number of function evaluations. The subroutine
  * adds whatever value is passed in, allowing the caller to keep a
  * running total.
@@ -71,7 +73,7 @@
 int linesearch(VecFVec &f, const UBVECTOR &x0,
                double f0, const UBVECTOR &g0,
                const UBVECTOR &dx, UBVECTOR &x,
-               double &fx, int &neval, std::ostream *solverlog = 0)
+               double &fx, UBVECTOR& fxVec, const double fxIncr, int &neval, std::ostream *solverlog = 0)
 {
   const double lseps = 1.0e-7;   // part of the definition of "sufficient" decrease
   const double TOLX = 1.0e-6;    // tolerance for x values
@@ -85,7 +87,6 @@ int linesearch(VecFVec &f, const UBVECTOR &x0,
   double g0dx=g0.dot(dx); // initial rate of decrease, df/dlambda
   double lambda = 1.0;           // start with full step
   double maxval = 0.0;
-  UBVECTOR fxVec(x0.size());
 
   if(g0dx >= 0) {
     if(solverlog)
@@ -122,7 +123,7 @@ int linesearch(VecFVec &f, const UBVECTOR &x0,
     // it specifies a minimum rate of decrease.  lseps is set fairly
     // small, so we're biased toward accepting steps unless their rate
     // of decrease is painfully slow
-    if(fx <= f0 + lseps*lambda*g0dx)
+    if(fx <= f0 + lseps*lambda*g0dx + fxIncr)
       // SUCCESS
       return 0;
 

cvs/objects/solution/util/include/solution_info.h

@@ -105,6 +105,7 @@ class SolutionInfo {
     double getED() const;
     double getEDLeft() const;
     double getEDRight() const;
+    double getSolutionFloor() const;
     void expandBracket( const double aAdjFactor );
     double getRelativeED() const;
     bool isWithinTolerance() const;

cvs/objects/solution/util/source/solution_info.cpp

@@ -227,6 +227,10 @@ double SolutionInfo::getEDRight() const {
     return EDR;
 }
 
+double SolutionInfo::getSolutionFloor() const {
+    return mSolutionFloor;
+}
+
 /*! \brief Return the name of the SolutionInfo object.
 * \author Josh Lurz
 * \return The name of the market the SolutionInfo is connected to.