diff --git a/NAMESPACE b/NAMESPACE index 930ee8c..e8843bf 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,11 +15,14 @@ export(get_fx) export(get_price_scale_factor) export(get_prices) export(get_quantity_scale_factor) +export(get_scenario_name) export(get_slope) export(get_supply) -export(run_to_period) +export(print_xmldb) +export(run_period) export(set_data) export(set_prices) +export(set_scenario_name) export(set_slope) importFrom(Rcpp,cpp_object_initializer) importFrom(Rcpp,loadModule) diff --git a/R/gcamwrapper.R b/R/gcamwrapper.R index 3f30b4e..03eb31e 100644 --- a/R/gcamwrapper.R +++ b/R/gcamwrapper.R @@ -21,18 +21,34 @@ create_and_initialize <- function(configuration = "configuration.xml", workdir = #' been run will be kept track of and will not be run again. HOWEVER, #' we do not attempt to keep track of if those model periods are "dirty" #' such as if a user has called `set_data` in such a way that would invalidate -#' that solution. +#' that solution. Users can supply a call back function which will be called +#' after initCalc and before solving. Such a call back could be necessary if +#' if data a user wanted to set would have been overridden during initializations. #' @param gcam (gcam) An initialized GCAM instance #' @param period (integer) The GCAM model period to run up to or the #" `get_current_period` + 1 if \code{NULL} +#' @param post_init_calback (function) A call back function +#' @param ... Additional params to pass to the callback function #' @return GCAM instance #' @importFrom Rcpp cpp_object_initializer #' @export -run_to_period <- function(gcam, period = NULL) { +run_period <- function(gcam, period = NULL, post_init_calback = NULL, ...) { if(is.null(period)) { period <- get_current_period(gcam) + 1 } - gcam$run_to_period(period) + + if(is.null(post_init_calback)) { + # if we are not running a call back function just do the whole thing together + # which will generate less confusing log messages + gcam$run_period(period) + } else { + # otherwise we will need to partition the calls + # if we are not already at period the log messages may be confusing as it will + # show as running period-1 then a seperate running period + gcam$run_period_pre(period) + post_init_calback(gcam, ...) + gcam$run_period_post(period) + } invisible(gcam) } @@ -120,10 +136,10 @@ convert_period_to_year <- function(gcam, period) { ret } -#' Convert from a GCAM model year to model period +#' Convert from a GCAM model year to model period #' @param gcam (gcam) An initialized GCAM instance #' @param year (integer vector) The model year to convert to period -#' @return (integer vector) The corresponding model period +#' @return (integer vector) The corresponding model period #' @export convert_year_to_period <- function(gcam, year) { if(length(year) == 0) { @@ -139,6 +155,34 @@ convert_year_to_period <- function(gcam, year) { ret } +#' Write the full results to an XML Database +#' @details At the moment all of the options governing XML Database output are +#' locked into the values set in the \code{configuration} as well as the +#' \code{XMLDBDriver.properties} file. +#' @param gcam (gcam) An initialized GCAM instance +#' @param xmldb_location (string) Override the location to write the XMLDB if non-empty +#' @return (string) The XMLDB location, particuarly useful if relying on defaults from the config. +#' @export +print_xmldb <- function(gcam, xmldb_location = "") { + gcam$print_xmldb(xmldb_location) +} + +#' Get the scenario name +#' @param gcam (gcam) An initialized GCAM instance +#' @return (string) The name of the scenario +#' @export +get_scenario_name <- function(gcam) { + gcam$get_scenario_name() +} + +#' Set the scenario name +#' @param gcam (gcam) An initialized GCAM instance +#' @param name (string) The name of the scenario to set +#' @export +set_scenario_name <- function(gcam, name) { + gcam$set_scenario_name(name) +} + #' Create a solution debugging object #' @details Create a solution debugging object which can be used a single #' evaluation of the model and see how it affects prices, supplies, @@ -146,13 +190,15 @@ convert_year_to_period <- function(gcam, year) { #' @param gcam (gcam) An initialized GCAM instance #' @param period (integer) GCAM model period to create the debugger or if NULL #' use the last run model period. +#' @param market_filter (string) A \code{} string in the same +#' format as in the solver config XML. The default is "solvable". #' @return A SolutionDebugger object #' @export -create_solution_debugger <- function(gcam, period = NULL) { +create_solution_debugger <- function(gcam, period = NULL, market_filter = "solvable") { if(is.null(period)) { period <- get_current_period(gcam) } - gcam$create_solution_debugger(period) + gcam$create_solution_debugger(period, market_filter) } #' Gets the prices diff --git a/R/query_library.R b/R/query_library.R index f1ce5e4..3956aac 100644 --- a/R/query_library.R +++ b/R/query_library.R @@ -47,6 +47,7 @@ wrap_get_query <- function() { #' @return (string) The raw query from the query file with units set as an attribute #' if specified. #' @importFrom yaml read_yaml +#' @export get_query <- wrap_get_query() diff --git a/gcamwrapper/main.py b/gcamwrapper/main.py index 66378c1..78296ca 100644 --- a/gcamwrapper/main.py +++ b/gcamwrapper/main.py @@ -10,7 +10,7 @@ class Gcam(gcam_module.gcam): instance. """ - def run_to_period(self, period=None): + def run_period(self, period=None, post_init_calback=None): """ Run GCAM up to and including some model period. Model periods which have already been run will be kept track of and will not be run again. HOWEVER, we do not attempt to @@ -25,7 +25,18 @@ def run_to_period(self, period=None): if period is None: period = self.get_current_period() + 1 - super(Gcam, self).run_to_period(period) + + if post_init_calback is None: + # if we are not running a call back function just do the whole thing together + # which will generate less confusing log messages + super(Gcam, self).run_period(period) + else: + # otherwise we will need to partition the calls + # if we are not already at period the log messages may be confusing as it will + # show as running period-1 then a seperate running period + super(Gcam, self).run_period_pre(period) + post_init_calback(self) + super(Gcam, self).run_period_post(period) def get_data(self, query, *args, **kwargs): """Queries for arbitrary data from a running instance of GCAM. @@ -143,13 +154,32 @@ def convert_year_to_period(self, year): else: return list(map(lambda x: self.convert_year_to_period(x), year)) - def create_solution_debugger(self, period=None): + def print_xmldb(self, xmldb_location = ""): + """Write the full results to an XML Database + At the moment all of the options governing XML Database output are + locked into the values set in the `configuration` as well as the + `XMLDBDriver.properties` file. + """ + + return super(Gcam, self).print_xmldb(xmldb_location) + + def get_scenario_name(self): + return super(Gcam, self).get_scenario_name() + + def set_scenario_name(self, name): + super(Gcam, self).set_scenario_name(name) + + + def create_solution_debugger(self, period=None, market_filter="solvable"): """Create a solution debugging object which can be used a single evaluation of the model and see how it affects prices, supplies, and demands amongst other things.. :param period: GCAM model period to create the debugger :type period: integer + :param market_filter: A `` string in the same + format as in the solver config XML. The default is "solvable". + :type market_filter: str :returns: A SolutionDebugger object @@ -157,7 +187,7 @@ def create_solution_debugger(self, period=None): if period is None: period = self.get_current_period() - sd = super(Gcam, self).create_solution_debugger(period) + sd = super(Gcam, self).create_solution_debugger(period, market_filter) sd.__class__ = SolutionDebugger return sd diff --git a/inst/extdata/query_library.yml b/inst/extdata/query_library.yml index ed8f35b..d08f449 100644 --- a/inst/extdata/query_library.yml +++ b/inst/extdata/query_library.yml @@ -1,14 +1,112 @@ +resource: + production: ["world/region{region@name}/resource{resource@name}/subresource/technology{tech@name}/period/output/physical-output{year@year}", "EJ/km3/others"] + price: ["world/region{region@name}/resource{resource@name}/price{year@year}", "1975$/GJ or m3"] + tech_change: ["world/region{region@name}/resource{resource@name}/subresource{subresource@name}/techChange{year@year}", "%/yr"] + cal_price_adder: ["world/region{region@name}/resource{resource@name}/subresource{subresource@name}/price-adder{year@year}", "1975$/GJ"] + # max resource for renewable resources + max_available: ["world/region{region@name}/resource{resource@name}/subresource{subresource@name}/maxSubResource{year@year}", "EJ or km3"] + # params for smooth renewable resources + price_exponent: ["world/region{region@name}/resource{resource@name}/subresource{subresource@name}/price-exponent", "unitless"] + mid_price: ["world/region{region@name}/resource{resource@name}/subresource{subresource@name}/mid-price", "1975$/GJ"] + grade: + available: ["world/region{region@name}/resource{resource@name}/subresource{subresource@name}/grade{grade@name}/available", "EJ/km3/others"] + extraction_cost: ["world/region{region@name}/resource{resource@name}/subresource{subresource@name}/grade{grade@name}/extractioncost", "1975$/GJ or m3"] + # the extraction cost adjusted for tech change + adjusted_cost: ["world/region{region@name}/resource{resource@name}/subresource{subresource@name}/grade{grade@name}/totalcost{year@year}", "1975$/GJ or m3"] + electricity: - share_weight: ["world/region{region@name}/sector[NamedFilter,StringEquals,electricity]/subsector{subsector@name}/share-weight{year@year}", "unitless"] + generation: ["world/region{region@name}/sector[+NamedFilter,StringEquals,electricity]/subsector{subsector@name}/technology{tech@name}/period/output{output@name}/physical-output{year@year}", "EJ"] + generation_cooling: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^elec_]/subsector{subsector@name}/technology{tech@name}/period/output{output@name}/physical-output{year@year}", "EJ"] + generation_cooling_vintage: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^elec_]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/output{output@name}/physical-output{year@year}", "EJ"] + inputs_cooling: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^elec_]/subsector{subsector@name}/technology{tech@name}/period/input{input@name}/physical-demand{year@year}", "EJ or km3"] + subsector_share_weight: ["world/region{region@name}/sector[NamedFilter,StringEquals,electricity]/subsector{subsector@name}/share-weight{year@year}", "unitless"] + +refining: + production: ["world/region{region@name}/sector[+NamedFilter,StringEquals,refining]/subsector{subsector@name}/technology{tech@name}/period/output{output@name}/physical-output{year@year}", "EJ"] + production_vintage: ["world/region{region@name}/sector[+NamedFilter,StringEquals,refining]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/output{output@name}/physical-output{year@year}", "EJ"] + input: ["world/region{region@name}/sector[+NamedFilter,StringEquals,refining]/subsector{subsector@name}/technology{tech@name}/period/input{input@name}/physical-demand{year@year}", "EJ"] + sector_cost: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^refin]/price", "1975$/GJ"] + tech_cost: ["world/region{region@name}/sector[+NamedFilter,StringEquals,refining]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/cost", "1975$/GJ"] + subsector_share_weight: ["world/region{region@name}/sector[+NamedFilter,StringEquals,refining]/subsector{subsector@name}/share-weight{year@year}", "unitless"] + tech_share_weight: ["world/region{region@name}/sector[+NamedFilter,StringEquals,refining]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/share-weight", "unitless"] + input_coef: ["world/region{region@name}/sector[+NamedFilter,StringEquals,refining]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/input{input@name}/current-coef{year@year}", "unitless"] + non_energy_cost: ["world/region{region@name}/sector[+NamedFilter,StringEquals,refining]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/input{input@name}/adjusted-cost{year@year}", "1975$/GJ"] + +buildings: + floorspace: ["world/region{region@name}/consumer{consumer@name}//nodeInput[+NamedFilter,StringRegexMatches,building]/base-building-size{year@year}", "billion m2"] + service: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(resid|comm)]/subsector{subsector@name}/technology{tech@name}/period/output{output@name}/physical-output{year@year}", "EJ"] + input: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(resid|comm)]/subsector{subsector@name}/technology{tech@name}/period/input{input@name}/physical-demand{year@year}", "EJ"] + sector_cost: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(resid|comm)]/price", "1975$/GJ"] + tech_cost: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(resid|comm)]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/cost", "1975$/GJ"] + subsector_share_weight: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(resid|comm)]/subsector{subsector@name}/share-weight{year@year}", "unitless"] + tech_share_weight: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(resid|comm)]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/share-weight", "unitless"] + input_coef: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(resid|comm)]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/input{input@name}/current-coef{year@year}", "unitless"] + non_energy_cost: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(resid|comm)]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/input{input@name}/adjusted-cost{year@year}", "1975$/GJ"] + +industry: + service: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(N fertilizer|cement|^industrial)]/subsector{subsector@name}/technology{tech@name}/period/output{output@name}/physical-output{year@year}", "EJ"] + input: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(N fertilizer|cement|^industrial)]/subsector{subsector@name}/technology{tech@name}/period/input{input@name}/physical-demand{year@year}", "EJ"] + sector_cost: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(N fertilizer|cement|^industrial)]/price", "1975$/GJ"] + tech_cost: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(N fertilizer|cement|^industrial)]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/cost", "1975$/GJ"] + subsector_share_weight: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(N fertilizer|cement|^industrial)]/subsector{subsector@name}/share-weight{year@year}", "unitless"] + tech_share_weight: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(N fertilizer|cement|^industrial)]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/share-weight", "unitless"] + input_coef: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(N fertilizer|cement|^industrial)]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/input{input@name}/current-coef{year@year}", "unitless"] + non_energy_cost: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,(N fertilizer|cement|^industrial)]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/input{input@name}/adjusted-cost{year@year}", "1975$/GJ"] + +transportation: + service: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^trn_]/subsector{subsector@name}/technology{tech@name}/period/output{output@name}/physical-output{year@year}", "million pass|ton-km"] + input: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^trn_]/subsector{subsector@name}/technology{tech@name}/period/input{input@name}/physical-demand{year@year}", "EJ"] + sector_cost: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^trn_]/price", "1990$/pass|ton-km"] + tech_cost: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^trn_]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/cost", "1990$/pass|ton-km"] + subsector_share_weight: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^trn_]/subsector{subsector@name}/share-weight{year@year}", "unitless"] + tech_share_weight: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^trn_]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/share-weight", "unitless"] + input_coef: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^trn_]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/input{input@name}/current-coef{year@year}", "unitless"] + non_energy_cost: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^trn_]/subsector{subsector@name}/technology{tech@name}/period{vintage@year}/input{input@name}/adjusted-cost{year@year}", "1975$/GJ"] + + socioeconomic: labor_productivity: ["world/region{region@name}/GDP/laborproductivity{year@year}", "unitless"] population: ["world/region{region@name}/demographic/population{year@year}/totalPop", "thous ppl"] + gdp: ["world/region{region@name}/GDP/gdpValue{year@year}", "million 1990$"] + emissions: #TODO: useful for get data but not set - co2_emissions: ["world/region{region@name}//ghg[NamedFilter,StringEquals,CO2]/emissions{year@year}", "MTC"] + co2_emissions_agg: ["world/region{region@name}//ghg[NamedFilter,StringEquals,CO2]/emissions{year@year}", "MTC"] + ghg_emissions_agg: ["world/region{region@name}//ghg{ghg@name}/emissions{year@year}", "Tg or MTC or Gg"] + emissions_sector: ["world/region{region@name}/sector{sector@name}//technology{tech@name}/period/ghg{ghg@name}/emissions{year@year}", "Tg or MTC or Gg"] + emissions_resource: ["world/region{region@name}/resource{resource@name}//technology{tech@name}/period/ghg{ghg@name}/emissions{year@year}", "Tg or MTC or Gg"] + sequestration: ["world/region{region@name}/sector{sector@name}//technology{tech@name}/period/capture-component/sequestered-amount{year@year}", "MTC"] + + ag: - prices: ["world/region{region@name}/sector{sector@name}/price{year@year}", "1975$/kg"] + price: ["world/region{region@name}/sector{sector@name}/price{year@year}", "1975$/kg"] + production: ["world/region{region@name}/sector{sector@name}/subsector{subsector@name}/technology[+NamedFilter,StringRegexMatches,(IRR|RFD)_(hi|lo)$]/period/output{output@name}/physical-output{year@year}", "Mt"] + input: ["world/region{region@name}/sector{sector@name}/subsector{subsector@name}/technology[+NamedFilter,StringRegexMatches,(IRR|RFD)_(hi|lo)$]/period/input{input@name}/physical-input{year@year}", "Mt N or km3"] yield: ["world/region{region@name}/sector{sector@name}/subsector/technology{tech@name}/period{year@year}/yield", "Mt/thous km2"] + variable_cost: ["world/region{region@name}/sector{sector@name}/subsector/technology{tech@name}/period{year@year}/nonLandVariableCost", "1975$/Mt"] land: land_allocation: ["world/region{region@name}//child-nodes{leaf@name}/land-allocation{year@year}", "thous km2"] + # WARNING: this query takes a long time to run, may want to apply filters + luc_emissions: ["world/region{region@name}//child-nodes{leaf@name}/carbon-calc/land-use-change-emissions{year@year}", "MTC"] + above_carbon_density: ["world/region{region@name}//child-nodes{leaf@name}/carbon-calc/above-ground-carbon-density", "kg/m2"] + below_carbon_density: ["world/region{region@name}//child-nodes{leaf@name}/carbon-calc/below-ground-carbon-density", "kg/m2"] + mature_ag: ["world/region{region@name}//child-nodes{leaf@name}/carbon-calc/mature-age", "kg/m2"] + +water: + withdrawals: + supply: ["world/region{region@name}/resource{resource@name}/subresource[NamedFilter,StringRegexMatches,(runoff|groundwater)]/technology{tech@name}/period/output/physical-output{year@year}", "km3"] + supply_mapping_desal: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^water_td_.*_W$]/subsector{subsector@name}/technology{tech@name}/period/output/physical-output{year@year}", "km3"] + demand_mapping_sector: ["world/region{region@name}/sector/subsector/technology/period/input[+NamedFilter,StringRegexMatches,_W$]/physical-demand{year@year}", "km3"] + demand: ["world/region{region@name}/sector{sector@name}/subsector{subsector@name}/technology{tech@name}/period/input[+NamedFilter,StringRegexMatches,_W$]/physical-demand{year@year}", "km3"] + price: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^water_td_.*_W$]/price{year@year}", "1975$/m3"] + consumption: + supply_mapping_desal: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^water_td_.*_C$]/subsector{subsector@name}/technology{tech@name}/period/output/physical-output{year@year}", "km3"] + demand_mapping_sector: ["world/region{region@name}/sector/subsector/technology/period/input[+NamedFilter,StringRegexMatches,_C$]/physical-demand{year@year}", "km3"] + demand: ["world/region{region@name}/sector{sector@name}/subsector{subsector@name}/technology{tech@name}/period/input[+NamedFilter,StringRegexMatches,_C$]/physical-demand{year@year}", "km3"] + price: ["world/region{region@name}/sector[+NamedFilter,StringRegexMatches,^water_td_.*_C$]/price{year@year}", "1975$/m3"] + +market: + price: ["marketplace/market{market@name}/market-period{year@year}/price", "various"] + supply: ["marketplace/market{market@name}/market-period{year@year}/supply", "various"] + demand: ["marketplace/market{market@name}/market-period{year@year}/demand", "various"] diff --git a/inst/include/interp_interface.h b/inst/include/interp_interface.h index 0774f97..3b16cb2 100644 --- a/inst/include/interp_interface.h +++ b/inst/include/interp_interface.h @@ -4,6 +4,10 @@ #include #include +// use boost::iostreams to wrap Interp API for cout +#include +#include + #if __has_include("Rcpp.h") #define IS_INTERP_R #define STRICT_R_HEADERS @@ -12,6 +16,8 @@ #define IS_INTERP_PYTHON #include #include + +#include #else #error "Could not determine, or using an unknown interpreter. Only R and Python are currently supported." #endif @@ -27,6 +33,26 @@ class gcam_exception : public std::exception { }; #if defined(USING_R) + +/*! + * \brief Wrap Rprintf in a boost::iostreams so we can use a C++ interface. + * \details We couldn't just use Rcout directly because we need to add a flush + * for it to show up in a notebook in a timeley manner. + */ +class RStdoutSink { +public: + typedef char char_type; + typedef boost::iostreams::sink_tag category; + + inline std::streamsize write(const char* aString, std::streamsize aSize) { + Rprintf("%.*s", aSize, aString); + // using R_FlushConsole() doesn't work for whatever reason + Rcpp::Function flush("flush.console"); + flush(); + return aSize; + } +}; + namespace Interp { using Rcpp::stop; using Rcpp::warning; @@ -36,6 +62,11 @@ namespace Interp { using Rcpp::IntegerVector; using Rcpp::NumericMatrix; + inline std::ostream& getInterpCout() { + static boost::iostreams::stream sRCout; + return sRCout; + } + inline std::string extract(const Rcpp::String& aStr) { return aStr; } @@ -84,6 +115,26 @@ namespace Interp { } #elif defined(PY_VERSION_HEX) + +/*! + * \brief Wrap PySys_WriteStdout in a boost::iostreams so we can use a C++ interface. + * \details This also bring us more inline with Rcpp and therefore keep a consisten abstract + * Interp interface. + */ +class PySysStdoutSink { +public: + typedef char char_type; + typedef boost::iostreams::sink_tag category; + + inline std::streamsize write(const char* aString, std::streamsize aSize) { + const std::streamsize MAX_PY_SIZE = 1000; + std::streamsize actualSize = std::min(aSize, MAX_PY_SIZE); + PySys_WriteStdout((boost::format("%%.%1%s") % actualSize).str().c_str(), aString); + boost::python::call_method(PySys_GetObject("stdout"), "flush"); + return actualSize; + } +}; + namespace Interp { static void stop(const std::string& aMessage) { throw gcam_exception(aMessage); @@ -94,6 +145,11 @@ namespace Interp { namespace bp = boost::python; namespace bnp = boost::python::numpy; + inline std::ostream& getInterpCout() { + static boost::iostreams::stream sPyCout; + return sPyCout; + } + inline std::string extract(const bp::str& aStr) { return bp::extract(aStr); } diff --git a/inst/include/query_processor_base.h b/inst/include/query_processor_base.h index 2f096ea..aebc476 100644 --- a/inst/include/query_processor_base.h +++ b/inst/include/query_processor_base.h @@ -34,7 +34,7 @@ class QueryProcessorBase { void parseFilterString(const std::string& aFilterStr ); - FilterStep* parseFilterStepStr( const std::string& aFilterStepStr, int& aCol ); + FilterStep* parseFilterStepStr( const std::string& aFilterStepStr, int& aCol, const bool aIsLastStep ); AMatchesValue* createMatchesAny() const; diff --git a/inst/include/solution_debugger.h b/inst/include/solution_debugger.h index 6387251..e5cf1c6 100644 --- a/inst/include/solution_debugger.h +++ b/inst/include/solution_debugger.h @@ -12,7 +12,7 @@ class Marketplace; class SolutionDebugger { public: - static SolutionDebugger createInstance(const int aPeriod); + static SolutionDebugger createInstance(const int aPeriod, const std::string& aMarketFilterStr); SolutionDebugger(World *w, Marketplace *m, SolutionInfoSet &sisin, int per); @@ -57,6 +57,8 @@ class SolutionDebugger { } #endif + void resetScales(const Interp::NumericVector& aPriceScale, const Interp::NumericVector& aQuantityScale); + private: World* world; @@ -68,8 +70,6 @@ class SolutionDebugger { UBVECTOR x; UBVECTOR fx; Interp::StringVector marketNames; - Interp::NumericVector priceScaleFactor; - Interp::NumericVector quantityScaleFactor; }; #endif // __SOLUTION_DEBUGGER_H__ diff --git a/man/create_solution_debugger.Rd b/man/create_solution_debugger.Rd index 8dc80ca..69279c0 100644 --- a/man/create_solution_debugger.Rd +++ b/man/create_solution_debugger.Rd @@ -4,13 +4,16 @@ \alias{create_solution_debugger} \title{Create a solution debugging object} \usage{ -create_solution_debugger(gcam, period = NULL) +create_solution_debugger(gcam, period = NULL, market_filter = "solvable") } \arguments{ \item{gcam}{(gcam) An initialized GCAM instance} \item{period}{(integer) GCAM model period to create the debugger or if NULL use the last run model period.} + +\item{market_filter}{(string) A \code{} string in the same +format as in the solver config XML. The default is "solvable".} } \value{ A SolutionDebugger object diff --git a/man/get_scenario_name.Rd b/man/get_scenario_name.Rd new file mode 100644 index 0000000..8a0d0a0 --- /dev/null +++ b/man/get_scenario_name.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gcamwrapper.R +\name{get_scenario_name} +\alias{get_scenario_name} +\title{Get the scenario name} +\usage{ +get_scenario_name(gcam) +} +\arguments{ +\item{gcam}{(gcam) An initialized GCAM instance} +} +\value{ +(string) The name of the scenario +} +\description{ +Get the scenario name +} diff --git a/man/print_xmldb.Rd b/man/print_xmldb.Rd new file mode 100644 index 0000000..7c3e3a6 --- /dev/null +++ b/man/print_xmldb.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gcamwrapper.R +\name{print_xmldb} +\alias{print_xmldb} +\title{Write the full results to an XML Database} +\usage{ +print_xmldb(gcam, xmldb_location = "") +} +\arguments{ +\item{gcam}{(gcam) An initialized GCAM instance} + +\item{xmldb_location}{(string) Override the location to write the XMLDB if non-empty} +} +\value{ +(string) The XMLDB location, particuarly useful if relying on defaults from the config. +} +\description{ +Write the full results to an XML Database +} +\details{ +At the moment all of the options governing XML Database output are +locked into the values set in the \code{configuration} as well as the +\code{XMLDBDriver.properties} file. +} diff --git a/man/run_to_period.Rd b/man/run_period.Rd similarity index 59% rename from man/run_to_period.Rd rename to man/run_period.Rd index c99d807..5f8190a 100644 --- a/man/run_to_period.Rd +++ b/man/run_period.Rd @@ -1,15 +1,19 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/gcamwrapper.R -\name{run_to_period} -\alias{run_to_period} +\name{run_period} +\alias{run_period} \title{Run model period} \usage{ -run_to_period(gcam, period = NULL) +run_period(gcam, period = NULL, post_init_calback = NULL, ...) } \arguments{ \item{gcam}{(gcam) An initialized GCAM instance} \item{period}{(integer) The GCAM model period to run up to or the} + +\item{post_init_calback}{(function) A call back function} + +\item{...}{Additional params to pass to the callback function} } \value{ GCAM instance @@ -23,5 +27,7 @@ and including the given period. Model periods which have already been run will be kept track of and will not be run again. HOWEVER, we do not attempt to keep track of if those model periods are "dirty" such as if a user has called `set_data` in such a way that would invalidate -that solution. +that solution. Users can supply a call back function which will be called +after initCalc and before solving. Such a call back could be necessary if +if data a uses wanted to set would have been overridden during initializations. } diff --git a/man/set_scenario_name.Rd b/man/set_scenario_name.Rd new file mode 100644 index 0000000..46a3809 --- /dev/null +++ b/man/set_scenario_name.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gcamwrapper.R +\name{set_scenario_name} +\alias{set_scenario_name} +\title{Set the scenario name} +\usage{ +set_scenario_name(gcam, name) +} +\arguments{ +\item{gcam}{(gcam) An initialized GCAM instance} + +\item{name}{(string) The name of the scenario to set} +} +\description{ +Set the scenario name +} diff --git a/setup.py b/setup.py index 44f6b8e..56be278 100644 --- a/setup.py +++ b/setup.py @@ -17,16 +17,15 @@ def get_requirements(): GCAM_LIB = os.environ["GCAM_LIB"] BOOST_INCLUDE = os.environ["BOOST_INCLUDE"] BOOST_LIB = os.environ["BOOST_LIB"] -XERCES_INCLUDE = os.environ["XERCES_INCLUDE"] -XERCES_LIB = os.environ["XERCES_LIB"] TBB_INCLUDE = os.environ["TBB_INCLUDE"] TBB_LIB = os.environ["TBB_LIB"] EIGEN_INCLUDE = os.environ["EIGEN_INCLUDE"] JAVA_INCLUDE = os.environ["JAVA_INCLUDE"] +JAVA_PLATFORM_INCLUDE = JAVA_INCLUDE + '/' + os.uname()[0].lower() JAVA_LIB = os.environ["JAVA_LIB"] # the exact set of libs to link will be platform specific due to boost -gcam_libs = ['gcam', 'hector', 'jvm', 'xerces-c', 'tbb', 'tbbmalloc', 'tbbmalloc_proxy'] +gcam_libs = ['gcam', 'hector', 'jvm', 'tbb', 'tbbmalloc', 'tbbmalloc_proxy'] gcam_compile_args = [] gcam_link_args = [] if platform.system() == "Windows" : @@ -41,7 +40,7 @@ def get_requirements(): # ensure we use the correct c++ std gcam_compile_args += ['-std=c++14'] # add rpath info to find the dynamic linked libs - gcam_link_args += ['-Wl,-rpath,'+XERCES_LIB, '-Wl,-rpath,'+JAVA_LIB, '-Wl,-rpath,'+BOOST_LIB, '-Wl,-rpath,'+TBB_LIB] + gcam_link_args += ['-Wl,-rpath,'+JAVA_LIB, '-Wl,-rpath,'+BOOST_LIB, '-Wl,-rpath,'+TBB_LIB] gcam_module = Extension( 'gcam_module', @@ -50,8 +49,8 @@ def get_requirements(): 'src/query_processor_base.cpp', 'src/set_data_helper.cpp', 'src/get_data_helper.cpp'], - include_dirs=['inst/include', GCAM_INCLUDE, BOOST_INCLUDE, XERCES_INCLUDE, TBB_INCLUDE, EIGEN_INCLUDE, JAVA_INCLUDE], - library_dirs=[GCAM_LIB, BOOST_LIB, XERCES_LIB, TBB_LIB, JAVA_LIB], + include_dirs=['inst/include', GCAM_INCLUDE, BOOST_INCLUDE, TBB_INCLUDE, EIGEN_INCLUDE, JAVA_INCLUDE, JAVA_PLATFORM_INCLUDE], + library_dirs=[GCAM_LIB, BOOST_LIB, TBB_LIB, JAVA_LIB], libraries=gcam_libs, define_macros=[('NDEBUG', '1'), ('BOOST_DATE_TIME_NO_LIB', '1'), diff --git a/src/Makevars b/src/Makevars index fc178cb..4ae1426 100644 --- a/src/Makevars +++ b/src/Makevars @@ -5,8 +5,6 @@ CXX_STD = CXX14 #GCAM_LIB = /gcam/cvs/objects/build/linux #BOOST_INCLUDE=/usr/include #BOOST_LIB=/usr/lib -#XERCES_INCLUDE=/usr/include -#XERCES_LIB=/usr/lib/x86_64-linux-gnu #TBB_INCLUDE=/usr/include #TBB_LIB=/usr/lib/x86_64-linux-gnu #EIGEN_INCLUDE=/usr/include @@ -15,7 +13,8 @@ CXX_STD = CXX14 #JAVA_INCLUDE=${JAVA_HOME}/include #JAVA_LIB=${JAVA_HOME}/lib/server -PKG_CPPFLAGS = -I../inst/include/ -I$(GCAM_INCLUDE) -I$(BOOST_INCLUDE) -I$(JAVA_INCLUDE) -I$(JAVA_INCLUDE)/darwin -I$(XERCES_INCLUDE) -I$(TBB_INCLUDE) -I$(EIGEN_INCLUDE) -DNDEBUG +OSNAME_LOWERCASE := $(shell uname -s | tr '[:upper:]' '[:lower:]') +PKG_CPPFLAGS = -I../inst/include/ -I$(GCAM_INCLUDE) -I$(BOOST_INCLUDE) -I$(JAVA_INCLUDE) -I$(JAVA_INCLUDE)/$(OSNAME_LOWERCASE) -I$(TBB_INCLUDE) -I$(EIGEN_INCLUDE) -DNDEBUG -PKG_LIBS = -Wl,-rpath,$(XERCES_LIB) -Wl,-rpath,$(JAVA_LIB) -lgcam -L/usr/local/lib -L$(XERCES_LIB) -L$(GCAM_LIB) -L$(JAVA_LIB) -lxerces-c -ljvm -L$(BOOST_LIB) -Wl,-rpath,$(BOOST_LIB) -L$(TBB_LIB) -Wl,-rpath,$(TBB_LIB) -lhector -lboost_system -lboost_filesystem -ltbb -ltbbmalloc -ltbbmalloc_proxy -lm +PKG_LIBS = -Wl,-rpath,$(JAVA_LIB) -lgcam -L/usr/local/lib -L$(GCAM_LIB) -L$(JAVA_LIB) -ljvm -L$(BOOST_LIB) -Wl,-rpath,$(BOOST_LIB) -L$(TBB_LIB) -Wl,-rpath,$(TBB_LIB) -lhector -lboost_system -lboost_filesystem -ltbb -ltbbmalloc -ltbbmalloc_proxy -lm diff --git a/src/gcam.cpp b/src/gcam.cpp index 6481bde..25d13f0 100644 --- a/src/gcam.cpp +++ b/src/gcam.cpp @@ -12,8 +12,13 @@ #include "util/logger/include/logger_factory.h" #include "util/base/include/timer.h" #include "util/base/include/version.h" -#include "util/base/include/xml_helper.h" +#include "util/base/include/xml_parse_helper.h" #include "util/base/include/manage_state_variables.hpp" +#include "util/base/include/time_vector.h" +#include "reporting/include/xml_db_outputter.h" +#include "containers/include/imodel_feedback_calc.h" +#include "marketplace/include/marketplace.h" +#include "containers/include/world.h" #include @@ -29,7 +34,8 @@ Scenario* scenario; class gcam { public: - gcam(string aConfiguration, string aWorkDir):isInitialized(false), mCurrentPeriod(0) { + gcam(string aConfiguration, string aWorkDir):isInitialized(false), mCurrentPeriod(0), mIsMidPeriod(false) { + mCoutOrig = cout.rdbuf(Interp::getInterpCout().rdbuf()); try { boost::filesystem::current_path(boost::filesystem::path(aWorkDir)); } catch(...) { @@ -45,20 +51,104 @@ class gcam { scenario->mManageStateVars = 0; runner.reset(0); scenario = 0; + Configuration::reset(); + cout.rdbuf(mCoutOrig); } - void runToPeriod(const int aPeriod ) { + void runPeriod(const int aPeriod ) { + if(!isInitialized) { + Interp::stop("GCAM did not successfully initialize."); + } + Timer timer; + + bool success = runner->runScenarios(aPeriod, false, timer); + if(!success) { + Interp::warning("Failed to solve period "+util::toString(aPeriod)); + } + mCurrentPeriod = aPeriod; + } + + void runPeriodPre(const int aPeriod ) { if(!isInitialized) { Interp::stop("GCAM did not successfully initialize."); } Timer timer; - bool success = runner->runScenarios(aPeriod, false, timer); - mCurrentPeriod = aPeriod; - if(!success) { - Interp::warning("Failed to solve period "+util::toString(aPeriod)); + if(aPeriod > 0 && (mCurrentPeriod+1) < aPeriod) { + bool success = runner->runScenarios(aPeriod-1, false, timer); + if(!success) { + Interp::warning("Failed to solve period "+util::toString(aPeriod)); + } } + mCurrentPeriod = aPeriod; + + scenario->logPeriodBeginning( aPeriod ); + + // If this is period 0 initialize market price. + if( aPeriod == 0 ){ + scenario->mMarketplace->initPrices(); // initialize prices + } + + // Run the iteration of the model. + scenario->mMarketplace->nullSuppliesAndDemands( aPeriod ); // initialize market demand to null + scenario->mMarketplace->init_to_last( aPeriod ); // initialize to last period's info + scenario->mWorld->initCalc( aPeriod ); // call to initialize anything that won't change during calc + scenario->mMarketplace->assignMarketSerialNumbers( aPeriod ); // give the markets their serial numbers for this period. + + // Call any model feedback objects before we begin solving this period but after + // we are initialized and ready to go. + for( auto modelFeedback : scenario->mModelFeedbacks ) { + modelFeedback->calcFeedbacksBeforePeriod( scenario, scenario->mWorld->getClimateModel(), aPeriod ); + } + + // Set up the state data for the current period. + delete scenario->mManageStateVars; + scenario->mManageStateVars = new ManageStateVariables( aPeriod ); + + // Be sure to clear out any supplies and demands in the marketplace before making our + // initial call to world.calc. There may already be values in there if for instance + // they got set from a restart file. + scenario->mMarketplace->nullSuppliesAndDemands( aPeriod ); + + scenario->mWorld->calc( aPeriod ); // call to calculate initial supply and demand + mIsMidPeriod = true; } + + void runPeriodPost(const int aPeriod) { + bool success = scenario->solve( aPeriod ); // solution uses Bisect and NR routine to clear markets + + scenario->mWorld->postCalc( aPeriod ); + + // Mark that the period is now valid. + scenario->mIsValidPeriod[ aPeriod ] = true; + + // Run the climate model for this period (only if the solver is successful) + if( !success ) { + ILogger& climatelog = ILogger::getLogger( "climate-log" ); + climatelog.setLevel( ILogger::WARNING ); + climatelog << "Solver unsuccessful for period " << aPeriod << "." << endl; + } + scenario->mWorld->runClimateModel( aPeriod ); + + // Call any model feedbacks now that we are done solving the current period and + // the climate model has been run. + for( auto modelFeedback : scenario->mModelFeedbacks ) { + modelFeedback->calcFeedbacksAfterPeriod( scenario, scenario->mWorld->getClimateModel(), aPeriod ); + } + + scenario->logPeriodEnding( aPeriod ); + + // Write out the results for debugging. + /* + if( aPrintDebugging ){ + scenario->writeDebuggingFiles( aXMLDebugFile, aTabs, aPeriod ); + } + */ + + delete scenario->mManageStateVars; + scenario->mManageStateVars = 0; + mIsMidPeriod = false; + } void setData(const Interp::DataFrame& aData, const std::string& aHeader) { if(!isInitialized) { Interp::stop("GCAM did not successfully initialize."); @@ -74,10 +164,16 @@ class gcam { return helper.run(runner->getInternalScenario()); } - SolutionDebugger createSolutionDebugger(const int aPeriod) { - delete scenario->mManageStateVars; - scenario->mManageStateVars = new ManageStateVariables(aPeriod); - return SolutionDebugger::createInstance(aPeriod); + SolutionDebugger createSolutionDebugger(const int aPeriod, const std::string& aMarketFilterStr) { + int period = aPeriod; + if(!mIsMidPeriod) { + delete scenario->mManageStateVars; + scenario->mManageStateVars = new ManageStateVariables(aPeriod); + } else if(aPeriod != mCurrentPeriod) { + Interp::warning("Solution debugger can only be created for current period "+util::toString(mCurrentPeriod)+" when running feedbacks."); + period = mCurrentPeriod; + } + return SolutionDebugger::createInstance(period, aMarketFilterStr); } int getCurrentPeriod() const { @@ -100,11 +196,52 @@ class gcam { return scenario->getModeltime()->getyr_to_per(aYear); } + std::string printXMLDB(const std::string& aXMLDBName) const { + if(!isInitialized) { + Interp::stop("GCAM did not successfully initialize."); + } + const std::string XMLDB_KEY = "xmldb-location"; + Configuration* conf = Configuration::getInstance(); + bool origShouldWrite = conf->shouldWriteFile(XMLDB_KEY); + bool origAppendScnName = conf->shouldAppendScnToFile(XMLDB_KEY); + std::string origXMLDBName = conf->getFile(XMLDB_KEY, XMLDB_KEY, false); + // override the config to always write + conf->mShouldWriteFileMap[XMLDB_KEY] = true; + if(!aXMLDBName.empty()) { + // override the config to not adjust the XMLDB name given + conf->mShouldAppendScnFileMap[XMLDB_KEY] = false; + conf->fileMap[XMLDB_KEY] = aXMLDBName; + } + XMLDBOutputter xmldb; + scenario->accept(&xmldb, -1); + xmldb.finish(); + xmldb.finalizeAndClose(); + if(!aXMLDBName.empty()) { + // reset original settings + conf->mShouldWriteFileMap[XMLDB_KEY] = origShouldWrite; + conf->mShouldAppendScnFileMap[XMLDB_KEY] = origAppendScnName; + conf->fileMap[XMLDB_KEY] = origXMLDBName; + } + return !aXMLDBName.empty() ? aXMLDBName : + origAppendScnName ? origXMLDBName.append(scenario->getName()) : + origXMLDBName; + } + + std::string getScenarioName() const { + return scenario->getName(); + } + + void setScenarioName(const std::string& aName) { + scenario->setName(aName); + } + private: bool isInitialized; + std::streambuf* mCoutOrig; int mCurrentPeriod; + bool mIsMidPeriod; LoggerFactoryWrapper loggerFactoryWrapper; - auto_ptr runner; + unique_ptr runner; void initializeScenario(string configurationArg) { string loggerFactoryArg = "log_conf.xml"; @@ -118,7 +255,8 @@ class gcam { // Initialize the LoggerFactory - bool success = XMLHelper::parseXML( loggerFileName, &loggerFactoryWrapper ); + XMLParseHelper::initParser(); + bool success = XMLParseHelper::parseXML( loggerFileName, &loggerFactoryWrapper ); // Check if parsing succeeded. Non-zero return codes from main indicate // failure. @@ -139,7 +277,7 @@ class gcam { mainLog << "Configuration file: " << configurationFileName << endl; mainLog << "Parsing input files..." << endl; Configuration* conf = Configuration::getInstance(); - success = XMLHelper::parseXML( configurationFileName, conf ); + success = XMLParseHelper::parseXML( configurationFileName, conf ); // Check if parsing succeeded. Non-zero return codes from main indicate // failure. if( !success ){ @@ -163,7 +301,7 @@ class gcam { } // Cleanup Xerces. This should be encapsulated with an initializer object to ensure against leakage. - XMLHelper::cleanupParser(); + XMLParseHelper::cleanupParser(); timer.stop(); @@ -182,13 +320,18 @@ RCPP_MODULE(gcam_module) { .constructor("constructor") - .method("run_to_period", &gcam::runToPeriod, "run to model period") + .method("run_period", &gcam::runPeriod, "run to model period") + .method("run_period_pre", &gcam::runPeriodPre, "run to model period pre solve") + .method("run_period_post", &gcam::runPeriodPost, "run to model period solve and post") .method("set_data", &gcam::setData, "set data") .method("get_data", &gcam::getData, "get data") .method("create_solution_debugger", &gcam::createSolutionDebugger, "create solution debugger") .method("get_current_period", &gcam::getCurrentPeriod, "get the last run model period") .method("convert_period_to_year", &gcam::convertPeriodToYear, "convert a GCAM model period to year") .method("convert_year_to_period", &gcam::convertYearToPeriod, "convert a GCAM model year to model period") + .method("print_xmldb", &gcam::printXMLDB, "write the full XML Database results") + .method("set_scenario_name", &gcam::setScenarioName, "set scenario name") + .method("get_scenario_name", &gcam::getScenarioName, "get scenario name") ; Rcpp::class_("SolutionDebugger") @@ -205,6 +348,7 @@ RCPP_MODULE(gcam_module) { .method("calc_derivative", &SolutionDebugger::calcDerivative, "calcDerivative") .method("get_slope", &SolutionDebugger::getSlope, "getSlope") .method("set_slope", &SolutionDebugger::setSlope, "setSlope") + .method("reset_scales", &SolutionDebugger::resetScales, "resetScales") ; } #elif defined(IS_INTERP_PYTHON) @@ -214,13 +358,18 @@ BOOST_PYTHON_MODULE(gcam_module) { boost::python::numpy::initialize(); class_("gcam", init()) - .def("run_to_period", &gcam::runToPeriod, "run to model period") + .def("run_period", &gcam::runPeriod, "run to model period") + .def("run_period_pre", &gcam::runPeriodPre, "run to model period pre solve") + .def("run_period_post", &gcam::runPeriodPost, "run to model period solve and post") .def("set_data", &gcam::setData, "set data") .def("get_data", &gcam::getData, "get data") .def("create_solution_debugger", &gcam::createSolutionDebugger, "create solution debugger") .def("get_current_period", &gcam::getCurrentPeriod, "get the last run model period") .def("convert_period_to_year", &gcam::convertPeriodToYear, "convert a GCAM model period to year") .def("convert_year_to_period", &gcam::convertYearToPeriod, "convert a GCAM model year to model period") + .def("print_xmldb", &gcam::printXMLDB, "write the full XML Database results") + .def("set_scenario_name", &gcam::setScenarioName, "set scenario name") + .def("get_scenario_name", &gcam::getScenarioName, "get scenario name") ; to_python_converter >(); to_python_converter >(); @@ -240,6 +389,7 @@ BOOST_PYTHON_MODULE(gcam_module) { .def("calc_derivative", &SolutionDebugger::calcDerivative, "calcDerivative") .def("get_slope", &SolutionDebugger::getSlope, "getSlope") .def("set_slope", &SolutionDebugger::setSlope_wrap, "setSlope") + .def("reset_scales", &SolutionDebugger::resetScales, "resetScales") ; } #endif diff --git a/src/get_data_helper.cpp b/src/get_data_helper.cpp index 3d70b63..e02c2be 100644 --- a/src/get_data_helper.cpp +++ b/src/get_data_helper.cpp @@ -222,7 +222,7 @@ template<> void GetDataHelper::processData(std::map& aData) { // the user did not explicitly specify a year filter for this data but // it seems better to assume they did rather than implicitly aggregate - // accross model periods + // across model periods // so let's add the path tracking filter to record all years the first // time we see this for them if(!mHasYearInPath) { @@ -233,7 +233,7 @@ void GetDataHelper::processData(std::map& aData) { // update the year path tracker which would not have had to match thus far // as an ARRAY is not a CONTAINER (*mPathTracker.rbegin())->matchesInt(GetIndexAsYear::convertIterToYear(aData, iter)); - // deletegate to processData to take care of the rest + // delegate to processData to take care of the rest processData((*iter).second); } } @@ -241,7 +241,7 @@ template void GetDataHelper::vectorDataHelper(VecType& aDataVec) { // the user did not explicitly specify a year filter for this data but // it seems better to assume they did rather than implicitly aggregate - // accross model periods + // across model periods // so let's add the path tracking filter to record all years the first // time we see this for them if(!mHasYearInPath) { @@ -252,7 +252,7 @@ void GetDataHelper::vectorDataHelper(VecType& aDataVec) { // update the year path tracker which would not have had to match thus far // as an ARRAY is not a CONTAINER (*mPathTracker.rbegin())->matchesInt(GetIndexAsYear::convertIterToYear(aDataVec, iter)); - // deletegate to processData to take care of the rest + // delegate to processData to take care of the rest processData(*iter); } } diff --git a/src/query_processor_base.cpp b/src/query_processor_base.cpp index b151ac9..ccbb4e7 100644 --- a/src/query_processor_base.cpp +++ b/src/query_processor_base.cpp @@ -106,15 +106,18 @@ AMatchesValue* QueryProcessorBase::wrapPredicate(AMatchesValue* aToWrap, const s /*! * \brief Parses a single filter step. - * \details Follows the same basic processing as the standard GCAM Fusion + * \details Follows the same basic processing as the standard GCAM Fusion * implementation of parseFilterStepStr however adds some additional * support for the `+` extension to map filters to a get / set operation. * \param aFilterStepStr The unparsed filter step string. * \param aCol The current column number which should correspond to a DataFrame * which is useful should this be a Set Data operation. + * \param aIsLastStep A boolean flag to indicate if this is the last filter step. This + * is useful to know to be able to implement some implicit behavior with + * respect to adding a "year" column if we match a period vector of data. * \return A FilterStep which is the parsed representation of the given aFilterStepStr. */ -FilterStep* QueryProcessorBase::parseFilterStepStr( const std::string& aFilterStepStr, int& aCol ) { +FilterStep* QueryProcessorBase::parseFilterStepStr( const std::string& aFilterStepStr, int& aCol, const bool aIsLastStep ) { auto openBracketIter = std::find( aFilterStepStr.begin(), aFilterStepStr.end(), '[' ); if( openBracketIter == aFilterStepStr.end() ) { // no filter just the data name @@ -160,7 +163,7 @@ FilterStep* QueryProcessorBase::parseFilterStepStr( const std::string& aFilterSt } else if( filterOptions[ 0 ] == "YearFilter" ) { if(isRead) { - matcher = wrapPredicate(matcher, "year", true); + matcher = wrapPredicate(matcher, aIsLastStep ? "year" : dataName, true); } filterStep = new FilterStep( dataName, new YearFilter( matcher ) ); @@ -188,7 +191,7 @@ void QueryProcessorBase::parseFilterString(const std::string& aFilterStr ) { mFilterSteps.resize(filterStepsStr.size()); int col = 0; for( size_t i = 0; i < filterStepsStr.size(); ++i ) { - mFilterSteps[ i ] = parseFilterStepStr( filterStepsStr[ i ], col ); + mFilterSteps[ i ] = parseFilterStepStr( filterStepsStr[ i ], col, (i+1) == filterStepsStr.size() ); } mDataColName = mFilterSteps.back()->mDataName; } diff --git a/src/solution_debugger.cpp b/src/solution_debugger.cpp index 1e63084..892a66b 100644 --- a/src/solution_debugger.cpp +++ b/src/solution_debugger.cpp @@ -1,9 +1,13 @@ #include "solution_debugger.h" +#include +#include + #include "containers/include/world.h" #include "solution/util/include/solution_info.h" -#include "solution/util/include/solvable_solution_info_filter.h" +#include "solution/util/include/isolution_info_filter.h" +#include "solution/util/include/solution_info_filter_factory.h" #include "solution/util/include/solution_info_param_parser.h" #include "marketplace/include/marketplace.h" #include "solution/util/include/fdjac.hpp" @@ -11,11 +15,19 @@ using namespace Interp; -SolutionDebugger SolutionDebugger::createInstance(const int aPeriod) { +SolutionDebugger SolutionDebugger::createInstance(const int aPeriod, const std::string& aMarketFilterStr) { SolutionInfoSet solnInfoSet( scenario->getMarketplace() ); SolutionInfoParamParser solnParams; solnInfoSet.init( aPeriod, 0.001, 0.001, &solnParams ); + std::unique_ptr filter(SolutionInfoFilterFactory::createSolutionInfoFilterFromString(aMarketFilterStr)); + if(filter.get()) { + solnInfoSet.updateSolvable(filter.get()); + } + else { + Interp::stop("Could not parse info filter: " + aMarketFilterStr); + } + return SolutionDebugger(scenario->getWorld(), scenario->getMarketplace(), solnInfoSet, aPeriod); } @@ -28,19 +40,13 @@ SolutionDebugger::SolutionDebugger(World *w, Marketplace *m, SolutionInfoSet &si F(solnInfoSet,w,m,per,false), x(nsolv), fx(nsolv), - marketNames(createVector(nsolv)), - priceScaleFactor(createVector(nsolv)), - quantityScaleFactor(createVector(nsolv)) + marketNames(createVector(nsolv)) { std::vector smkts(solnInfoSet.getSolvableSet()); for(unsigned int i = 0; i < smkts.size(); ++i) { x[i] = smkts[i].getPrice(); marketNames[i] = smkts[i].getName().c_str(); - priceScaleFactor[i] = std::max(smkts[i].getForecastPrice(), 0.001); - quantityScaleFactor[i] = std::abs(smkts[i].getForecastDemand()); } - setVectorNames(priceScaleFactor, marketNames); - setVectorNames(quantityScaleFactor, marketNames); F.scaleInitInputs(x); F(x, fx); @@ -55,7 +61,7 @@ NumericVector SolutionDebugger::getPrices(const bool aScaled) { for(int i = 0; i < nsolv; ++i) { double val = x[i]; if(!aScaled) { - val *= priceScaleFactor[i]; + val *= F.mxscl[i]; } ret[i] = val; } @@ -75,7 +81,7 @@ NumericVector SolutionDebugger::getSupply(const bool aScaled) { for(int i = 0; i < nsolv; ++i) { double val = solnInfoSet.getSolvable(i).getSupply(); if(aScaled) { - val /= quantityScaleFactor[i]; + val /= 1.0/F.mfxscl[i]; } ret[i] = val; } @@ -88,7 +94,7 @@ NumericVector SolutionDebugger::getDemand(const bool aScaled) { for(int i = 0; i < nsolv; ++i) { double val = solnInfoSet.getSolvable(i).getDemand(); if(aScaled) { - val /= quantityScaleFactor[i]; + val /= 1.0/F.mfxscl[i]; } ret[i] = val; } @@ -97,11 +103,21 @@ NumericVector SolutionDebugger::getDemand(const bool aScaled) { } NumericVector SolutionDebugger::getPriceScaleFactor() { - return priceScaleFactor; + NumericVector ret(createVector(nsolv)); + for(int i = 0; i < nsolv; ++i) { + ret[i] = F.mxscl[i]; + } + setVectorNames(ret, marketNames); + return ret; } NumericVector SolutionDebugger::getQuantityScaleFactor() { - return quantityScaleFactor; + NumericVector ret(createVector(nsolv)); + for(int i = 0; i < nsolv; ++i) { + ret[i] = 1.0/F.mfxscl[i]; + } + setVectorNames(ret, marketNames); + return ret; } void SolutionDebugger::setPrices(const NumericVector& aPrices, const bool aScaled) { @@ -122,14 +138,18 @@ NumericVector SolutionDebugger::evaluate(const NumericVector& aPrices, const boo } setPrices(aPrices, aScaled); + std::unique_ptr resetState; if(aResetAfterCalc) { - scenario->getManageStateVariables()->setPartialDeriv(true); - F.partial(1); + resetState.reset(new double[scenario->getManageStateVariables()->mNumCollected]); + memcpy(resetState.get(), + scenario->getManageStateVariables()->mStateData[0], + (sizeof( double)) * scenario->getManageStateVariables()->mNumCollected); } F(x,fx); NumericVector fx_ret = getFX(); if(aResetAfterCalc) { - F.partial(-1); + memcpy(scenario->getManageStateVariables()->mStateData[0], resetState.get(), + (sizeof( double)) * scenario->getManageStateVariables()->mNumCollected); x = x_restore; fx = fx_restore; } @@ -139,7 +159,7 @@ NumericVector SolutionDebugger::evaluate(const NumericVector& aPrices, const boo NumericVector SolutionDebugger::evaluatePartial(const double aPrice, const int aIndex, const bool aScaled) { double x_restore = x[aIndex]; - x[aIndex] = aScaled ? aPrice : aPrice / priceScaleFactor[aIndex]; + x[aIndex] = aScaled ? aPrice : aPrice / F.mxscl[aIndex]; scenario->getManageStateVariables()->setPartialDeriv(true); F.partial(aIndex); UBVECTOR fx_restore = fx; @@ -153,8 +173,12 @@ NumericVector SolutionDebugger::evaluatePartial(const double aPrice, const int a } NumericMatrix SolutionDebugger::calcDerivative() { + std::list indicies; + for(int i = 0; i < nsolv; ++i) { + indicies.push_back(i); + } UBMATRIX jac(nsolv,nsolv); - fdjac(F, x, fx, jac, true); + fdjac(F, x, fx, jac, indicies, true); NumericMatrix jacRet = wrapMatrix(jac, nsolv); Interp::setMatrixNames(jacRet, marketNames); @@ -164,7 +188,7 @@ NumericMatrix SolutionDebugger::calcDerivative() { NumericVector SolutionDebugger::getSlope() { NumericVector slope(createVector(nsolv)); for(int i = 0; i < nsolv; ++i) { - //slope[i] = solnInfoSet.getSolvable(i).getCorrectionSlope(); + slope[i] = solnInfoSet.getSolvable(i).getCorrectionSlope(F.mxscl[i], 1.0/F.mfxscl[i]); } setVectorNames(slope, marketNames); return slope; @@ -178,3 +202,19 @@ void SolutionDebugger::setSlope(const NumericVector& aDX) { F.setSlope(dx); } +void SolutionDebugger::resetScales(const NumericVector& aPriceScale, + const NumericVector& aQuantityScale) +{ + // changing scales will invalidate F, x, and fx so we will have to recreate these + std::vector smkts(solnInfoSet.getSolvableSet()); + for(int i = 0; i < nsolv; ++i) { + x[i] = smkts[i].getPrice(); + smkts[i].setForecastPrice(aPriceScale[i]); + smkts[i].setForecastDemand(aQuantityScale[i]); + } + LogEDFun Fnew(solnInfoSet, world, marketplace, period, false); + //std::swap(F, Fnew); + F.scaleInitInputs(x); + F(x, fx); +} +