Pull request #306: Enl/bugfix/misc 6p1 fixes

Merge in JGCRI/gcam-core from enl/bugfix/misc_6p1_fixes to master

* commit '596d50c7d2c137e82f9f709f496c19bef901575b':
  Some final things to do before next release: * Update target finder 2.6 initial guesses * Update calibrate macro TFP * Update to version.h to say 7.0
  Rebuild package documentation one last time
  Up max iteration count to adjust for derivative calc inflation
  Fix missing metadata warning on negative emissions tibble
  Avoid warnings about CO2_LUC not related to an Activity
  Pull in Hector Xcode project fix
  A set of solution tweaks
  Add prune_empty_ag.xml to config and batch files
  Clean up precursor warnings
  Data system changes to add prune_empty_ag.xml
  small cleanup: - fix units in rsrc_unconv_oil_prod_bbld.csv (million barrels of oil per day, not billions) - forcing_target_ files:  fix comment related to hotelling rate (updated to 3% in GCAM v6.0, but the comment says <!-- path-discount-rate | default: 0.05 | The hotelling rate -->)
  revise LB162 to make sure no negative AgProdChange
  Ensure emissions unit tag is always copied forward for consistent reporting
  Fix emission control bugs

cvs/objects/containers/source/market_dependency_finder.cpp

@@ -537,7 +537,9 @@ void MarketDependencyFinder::createOrdering() {
     // activities in the model as it may be an indication of misconfiguration.
     ILogger& depLog = ILogger::getLogger( "dependency_finder_log" );
     for( CItemIterator it = mDependencyItems.begin(); it != mDependencyItems.end(); ++it ) {
-        if( !(*it)->mHasIncomingDependency && (*it)->mDependentList.empty() ) {
+        // less than ideal but our heuristics do not properly deal with link CO2 markets
+        // so explicitly avoid warning about CO2_LUC (potentially others here)
+        if( !(*it)->mHasIncomingDependency && (*it)->mDependentList.empty() && (*it)->mName != "CO2_LUC" ) {
             depLog.setLevel( ILogger::SEVERE );
             depLog << (*it)->mName << " in " << (*it)->mLocatedInRegion << " is not related to any other activities." << endl;
         }

cvs/objects/emissions/source/nonco2_emissions.cpp

@@ -120,6 +120,8 @@ void NonCO2Emissions::copyGHGParameters( const AGHG* aPrevGHG ){
         abort();
     }
 
+    mEmissionsUnit = prevComplexGHG->mEmissionsUnit;
+
     if( !mEmissionsDriver ) {
         mEmissionsDriver = prevComplexGHG->mEmissionsDriver->clone();
     }

cvs/objects/marketplace/source/demand_market.cpp

@@ -123,5 +123,5 @@ bool DemandMarket::shouldSolveNR() const {
     // A demand market is solving an equality constraint, so it should
     // be solved even if it has negative values in it (edfun.cpp makes
     // special allowances for price, demand, and trial value markets)
-    return shouldSolve();
+    return shouldSolve() && fabs(mPrice) > util::getSmallNumber() && fabs(mDemand) > util::getSmallNumber();
 }

cvs/objects/marketplace/source/price_market.cpp

@@ -141,5 +141,5 @@ bool PriceMarket::shouldSolveNR() const {
     // A price market is solving an equality constraint, so we should
     // try to solve it even when price, supply, and demand are all
     // negative.
-    return shouldSolve();
+    return shouldSolve() && fabs(mPrice) > util::getSmallNumber() && fabs(mDemand) > util::getSmallNumber();
 }

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

@@ -401,7 +401,7 @@ int LogBroyden::bsolve(VecFVec &F, UBVECTOR &x, UBVECTOR &fx,
     const int TRACK_NUM_PAST_F_VALUES = 4;
     std::queue<double> past_f_values;
     past_f_values.push(f0);
-    const double LARGE_DIAG_THRESHOLD = 100.0;
+    const double LARGE_DIAG_THRESHOLD = 1000000.0;
   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
@@ -497,7 +497,7 @@ int LogBroyden::bsolve(VecFVec &F, UBVECTOR &x, UBVECTOR &fx,
       // start with partial pivot LU decomposition as it is so fast to execute and
       // see if we need to fall back to an alternative if it doesn't "perform" well
       Eigen::PartialPivLU<UBMATRIX> luPartialPiv(B);
-      dx = luPartialPiv.solve(-1.0 * fx);
+      dx = luPartialPiv.solve(-fx);
       double dxmag = sqrt(dx.dot(dx));
       if(luPartialPiv.determinant() == 0 || !util::isValidNumber(dxmag)) {
           // singular or badly messed up Jacobian, going to have to use SVD
@@ -512,7 +512,7 @@ int LogBroyden::bsolve(VecFVec &F, UBVECTOR &x, UBVECTOR &fx,
           // without that we should set this to zero (i.e. only suppress truly singular cols)
           const double small_threshold = 0;
           svdSolver.setThreshold(small_threshold);
-          dx = svdSolver.solve(-1.0 * fx);
+          dx = svdSolver.solve(-fx);
           dxmag = sqrt(dx.dot(dx));
           solverLog << " new dxmag: " << dxmag << std::endl;
       }
@@ -558,7 +558,7 @@ int LogBroyden::bsolve(VecFVec &F, UBVECTOR &x, UBVECTOR &fx,
       // 1) B is not a descent direction.  If this is the first
       // failure starting from this x value, try a finite difference
       // jacobian
-      if(!lsfail) {
+      /*if(!lsfail) {
           // linesearch will have left off on some other price vector thus we
           // could have bad state data from which we calculate derivatives or
           // check if we we can return due to a relaxed convergence test
@@ -590,7 +590,7 @@ int LogBroyden::bsolve(VecFVec &F, UBVECTOR &x, UBVECTOR &fx,
           past_f_values.push(f0);
         // start the next iteration *without* updating x
         continue;
-      }
+      }*/
 
       // 2) The descent only continues for a very short distance
       // (roughly TOL*x0).  Maybe we're really close to a solution and

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

@@ -287,7 +287,7 @@ SolverComponent::ReturnCode Preconditioner::solve( SolutionInfoSet& aSolutionSet
                         chg = true;
                         ++nchg;
                     } 
-                    else if (oldprice > ub &&
+                    /*else if (oldprice > ub &&
                              oldsply > olddmnd) {
                       // price is above the top of the supply curve,
                       // and there little demand.  This is not
@@ -317,7 +317,7 @@ SolverComponent::ReturnCode Preconditioner::solve( SolutionInfoSet& aSolutionSet
                       solvable[i].setPrice(newprice);
                       chg = true;
                       ++nchg;
-                    }
+                    }*/
                     else if(fd < mFTOL && oldsply >= mFTOL && olddmnd >= mFTOL) {
                         // reset a small demand scale when it looked like a market was "off"
                         // but the supplies and demands are now in a normal range
@@ -354,6 +354,17 @@ SolverComponent::ReturnCode Preconditioner::solve( SolutionInfoSet& aSolutionSet
                         chg = true;
                         ++nchg;
                     }
+                    else if(abs(fp) != abs(fd)) {
+                        // not sure why the price and quantity scales divereged but they should
+                        // be the same so reset them to be
+                        double newScale = abs(fp) > util::getSmallNumber() ? fp : fd;
+                        fp = newScale;
+                        fd = newScale;
+                        solvable[i].setForecastPrice(newScale);
+                        solvable[i].setForecastDemand(newScale);
+                        chg = true;
+                        ++nchg;
+                    }
                     break; 
                 case IMarketType::PRICE:
                     // price markets are solving a consistency

cvs/objects/util/base/include/version.h

@@ -41,10 +41,10 @@
  * NOTE: ADD 1 TO LATEST SUBVERSION REVISION NUMBER
  */
 //! The latest SVN revision number for identification of the build.
-#define __REVISION_NUMBER__ "gcam-v6.0"
+#define __REVISION_NUMBER__ "gcam-v7.0"
 /*****************************************************************************/
 
 //! GCAM model version.
-#define __ObjECTS_VER__ "6.0"
+#define __ObjECTS_VER__ "7.0"
 
 #endif // _VERSION_H_ 

exe/batch_SSP_REF.xml

@@ -137,6 +137,7 @@
 			<!-- WASTEWATER TREATMENT FRACTIONS -->
 			<Value name="wwtrt">../input/gcamdata/xml/EFW_WWtrt_coefs_SSP4.xml</Value>
 
+            <Value name = "prune_ag">../input/gcamdata/xml/prune_empty_ag.xml</Value>
 		</FileSet>
 
 		<FileSet name="SSP5">
@@ -174,6 +175,7 @@
 			<!-- WASTEWATER TREATMENT FRACTIONS -->
 			<Value name="wwtrt">../input/gcamdata/xml/EFW_WWtrt_coefs_SSP5.xml</Value>
 
+            <Value name = "prune_ag">../input/gcamdata/xml/prune_empty_ag.xml</Value>
 		</FileSet>
 
 		<FileSet name="SSP1">
@@ -229,6 +231,7 @@
 			<!-- WASTEWATER TREATMENT FRACTIONS -->
 			<Value name="wwtrt">../input/gcamdata/xml/EFW_WWtrt_coefs_SSP1.xml</Value>
 
+            <Value name = "prune_ag">../input/gcamdata/xml/prune_empty_ag.xml</Value>
 		</FileSet>
 
 		<FileSet name="SSP2">
@@ -264,6 +267,7 @@
 			<!-- WASTEWATER TREATMENT FRACTIONS -->
 			<Value name="wwtrt">../input/gcamdata/xml/EFW_WWtrt_coefs_SSP2.xml</Value>
 
+            <Value name = "prune_ag">../input/gcamdata/xml/prune_empty_ag.xml</Value>
 		</FileSet>
 
 		<FileSet name="SSP3">
@@ -321,6 +325,7 @@
 			<!-- WASTEWATER TREATMENT FRACTIONS -->
 			<Value name="wwtrt">../input/gcamdata/xml/EFW_WWtrt_coefs_SSP3.xml</Value>
 
+            <Value name = "prune_ag">../input/gcamdata/xml/prune_empty_ag.xml</Value>
 		</FileSet>		
 	</ComponentSet>
 

exe/batch_SSP_SPA1.xml

@@ -145,6 +145,7 @@
 			<Value name="link">../input/policy/global_uct_spa1.xml</Value>
             <Value name = "bio_neg_emiss_budget">../input/gcamdata/xml/negative_emissions_budget.xml</Value>
             <Value name = "bio_externality">../input/gcamdata/xml/bio_externality.xml</Value>
+            <Value name = "prune_ag">../input/gcamdata/xml/prune_empty_ag.xml</Value>
 		</FileSet>
 
 

exe/batch_SSP_SPA23.xml

@@ -125,7 +125,7 @@
 			<Value name="bio_tax">../input/policy/global_uct_spa2.xml</Value>
             <Value name = "bio_neg_emiss_budget">../input/gcamdata/xml/negative_emissions_budget.xml</Value>
             <Value name = "bio_externality">../input/gcamdata/xml/bio_externality.xml</Value>
-			
+            <Value name = "prune_ag">../input/gcamdata/xml/prune_empty_ag.xml</Value>
 		</FileSet>
 		<FileSet name="SSP3">
 			<!-- SOCIOECONOMICS -->
@@ -189,6 +189,7 @@
 			<Value name="link">../input/policy/2040_target_finder_phasein.xml</Value>
             <Value name = "bio_neg_emiss_budget">../input/gcamdata/xml/negative_emissions_budget_SSP3.xml</Value>
             <Value name = "bio_externality">../input/gcamdata/xml/bio_externality.xml</Value>
+            <Value name = "prune_ag">../input/gcamdata/xml/prune_empty_ag.xml</Value>
 		</FileSet>	
 	</ComponentSet>
 	<runner-set name="policy-target-runner">

exe/batch_SSP_SPA4.xml

@@ -147,6 +147,7 @@
 			<Value name="bio_trade">../input/gcamdata/xml/ssp4_ag_bio_trade.xml</Value>
             <Value name = "bio_neg_emiss_budget">../input/gcamdata/xml/negative_emissions_budget.xml</Value>
             <Value name = "bio_externality">../input/gcamdata/xml/bio_externality.xml</Value>
+            <Value name = "prune_ag">../input/gcamdata/xml/prune_empty_ag.xml</Value>
 		</FileSet>
 
 	</ComponentSet>	

exe/batch_SSP_SPA5.xml

@@ -127,6 +127,7 @@
 			<Value name="link">../input/policy/global_uct_spa5.xml</Value>
             <Value name = "bio_neg_emiss_budget">../input/gcamdata/xml/negative_emissions_budget.xml</Value>
             <Value name = "bio_externality">../input/gcamdata/xml/bio_externality.xml</Value>
+            <Value name = "prune_ag">../input/gcamdata/xml/prune_empty_ag.xml</Value>
 		</FileSet>
 	</ComponentSet>
 	<runner-set name="policy-target-runner">

exe/configuration_policy.xml

@@ -105,6 +105,8 @@
 
 		<Value name = "macro">../input/gcamdata/xml/socioeconomics_macro.xml</Value>
 
+        <Value name = "prune_ag">../input/gcamdata/xml/prune_empty_ag.xml</Value>
+
         <!-- set up the near-term / long-term policy inputs -->
 		<Value name = "long-term-co2">../input/policy/carbon_tax_0.xml</Value>
 		<Value name = "near-term-co2">../input/policy/spa14_tax.xml</Value>

exe/configuration_ref.xml

@@ -102,6 +102,8 @@
 		<Value name = "bio_externality">../input/gcamdata/xml/bio_externality.xml</Value>
         <Value name = "gas_trade">../input/gcamdata/xml/gas_trade.xml</Value>
 
+        <Value name = "prune_ag">../input/gcamdata/xml/prune_empty_ag.xml</Value>
+
 		<Value name = "solver">../input/solution/cal_broyden_config.xml</Value>
 
 	</ScenarioComponents>

exe/configuration_usa.xml

@@ -103,6 +103,8 @@
         <Value name = "wind_update">../input/gcamdata/xml/onshore_wind.xml</Value>
         <Value name = "gas_trade">../input/gcamdata/xml/gas_trade.xml</Value>
 
+        <Value name = "prune_ag">../input/gcamdata/xml/prune_empty_ag.xml</Value>
+
 		<Value name = "socio_usa">../input/gcamdata/xml/socioeconomics_USA.xml</Value>
 		<Value name = "resource_usa">../input/gcamdata/xml/resources_USA.xml</Value>
 		<Value name = "Cstorage_usa">../input/gcamdata/xml/Cstorage_USA.xml</Value>

input/gcamdata/R/zaglu_L100.FAO_SUA_PrimaryEquivalent.R

@@ -714,38 +714,34 @@ module_aglu_L100.FAO_SUA_PrimaryEquivalent <- function(command, ...) {
                                     "Food", "Feed", "Other uses") %in%
                                     c(.DF %>% distinct(element) %>% pull)))
       # 0. Check NA
-      if (isTRUE(.DF %>% filter(is.na(value)) %>% nrow == 0)) {
-        message("Good! No NA in value") } else{
-          warning("NA in value")
-        }
+      if (.DF %>% filter(is.na(value)) %>% nrow() > 0) {
+        warning("NA values in SUA Balance")
+      }
 
 
       # 1. Positive value except stock variation and residues
-      if (isTRUE(.DF %>% filter(!element %in% c("Stock Variation", "Other uses")) %>%
+      if (isFALSE(.DF %>% filter(!element %in% c("Stock Variation", "Other uses")) %>%
                  summarise(min = min(value, na.rm = T)) %>% pull(min) >= -0.001)) {
-        message("Good! Signs checked") } else{
-          warning("Negative values in key elements (not including stock variation and other uses)")
-        }
+        warning("Negative values in key elements (not including stock variation and other uses)")
+      }
 
       # 2. Trade balance in all year and items
-      if (isTRUE(.DF %>% filter(element %in% c("Import", "Export")) %>%
+      if (isFALSE(.DF %>% filter(element %in% c("Import", "Export")) %>%
                  group_by(year, GCAM_commodity, element) %>%
                  summarise(value = sum(value), .groups = "drop") %>%
                  spread(element, value) %>% filter(abs(Import - Export) > 0.0001) %>% nrow() == 0)) {
-        message("Good! Gross trade in balance") } else{
-          warning("Gross trade imbalance")
-        }
+        warning("Gross trade imbalance")
+      }
 
       # 3. SUA balance check
-      if (isTRUE(.DF %>%
+      if (isFALSE(.DF %>%
                  spread(element, value) %>%
                  mutate(`Regional supply` = Production + `Import`,
                         `Regional demand` = `Export` + Feed + Food  + `Other uses`,
                         bal = abs(`Regional supply` -  `Regional demand`)) %>%
                  filter(bal > 0.0001) %>% nrow() == 0)) {
-        message("Good! Regional supply = Regional demand + Residuals") } else{
-          warning("Regional supply != Regional demand + Residuals")
-        }
+        warning("Regional supply != Regional demand + Residuals")
+      }
 
       # 4. Balanced in all dimensions
       assertthat::assert_that(.DF %>% nrow() ==

input/gcamdata/R/zaglu_L100.GTAP_downscale_ctry.R

@@ -25,8 +25,6 @@ module_aglu_L100.GTAP_downscale_ctry <- function(command, ...) {
       FILE = "aglu/FAO/FAO_ag_items_PRODSTAT",
       "L100.LDS_value_milUSD",
       "L100.LDS_ag_prod_t",
-      FILE = "common/iso_GCAM_regID",
-      FILE = "common/GCAM_region_names",
       FILE = "socioeconomics/GTAP/GCAM_GTAP_region_mapping",
       FILE = "socioeconomics/GTAP/GTAP_sector_aggregation_mapping",
       OPTIONAL_FILE = "socioeconomics/GTAP/GTAPv10_basedata_VKB_SAVE_VDEP",

input/gcamdata/R/zaglu_L162.ag_prodchange_R_C_Y_GLU_irr.R

@@ -112,6 +112,11 @@ module_aglu_L162.ag_prodchange_R_C_Y_GLU_irr <- function(command, ...) {
       mutate(Irr_Rfd = "RFD") %>%
       rename(yield_kgHa = Yield_kgHa_rainfed) %>%
       bind_rows(L162.ag_irrYield_kgHa_Rcrs_Ccrs_Y) %>%
+      group_by(CROSIT_ctry, CROSIT_crop, Irr_Rfd) %>%
+      mutate(tag1 = ifelse(yield_kgHa[year == 2030] < yield_kgHa[year == 2005], 1, 0), # if 2030 < 2005, then AgProdChange1 = 0
+             yield_kgHa = ifelse(tag1 == 1 & year == 2030, yield_kgHa[year == 2005], yield_kgHa),
+             tag2 = ifelse(yield_kgHa[year == 2050] < yield_kgHa[year == 2030], 1, 0), # if 2050 < 2030, then AgProdChange2 = AgProdChange1
+             yield_kgHa = ifelse(tag2 == 1 & year == 2050, yield_kgHa[year == 2030] + 4*(yield_kgHa[year == 2030] - yield_kgHa[year == 2005])/5, yield_kgHa)) %>%
       # add the missing aglu.SPEC_AG_PROD_YEARS and interpolate the yields
       complete(year = c(year, aglu.SPEC_AG_PROD_YEARS) ,
                       CROSIT_ctry, CROSIT_crop, Irr_Rfd) %>%

input/gcamdata/R/zaglu_L2072.ag_water_irr_mgmt.R

@@ -26,7 +26,6 @@ module_aglu_L2072.ag_water_irr_mgmt <- function(command, ...) {
   MODULE_INPUTS <-
     c(FILE = "common/GCAM_region_names",
       FILE = "water/basin_to_country_mapping",
-      "L2012.AgSupplySector",
       "L161.ag_irrProd_Mt_R_C_Y_GLU",
       "L161.ag_rfdProd_Mt_R_C_Y_GLU",
       "L165.BlueIrr_m3kg_R_C_GLU",