diff --git a/include/bmi/Bmi_Adapter.hpp b/include/bmi/Bmi_Adapter.hpp index 54c3f1b50b..5a0cc455a8 100644 --- a/include/bmi/Bmi_Adapter.hpp +++ b/include/bmi/Bmi_Adapter.hpp @@ -32,9 +32,10 @@ namespace models { { // This replicates a lot of Initialize, but it's necessary to be able to do it separately to support // "initializing" on construction, given using Initialize requires use of virtual functions + errno = 0; if (!utils::FileChecker::file_is_readable(this->bmi_init_config)) { init_exception_msg = "Cannot create and initialize " + this->model_name + " using unreadable file '" - + this->bmi_init_config + "'"; + + this->bmi_init_config + "'. Error: "+std::strerror(errno); throw std::runtime_error(init_exception_msg); } } @@ -125,6 +126,7 @@ namespace models { void Initialize() { // If there was previous init attempt but w/ failure exception, throw runtime error and include previous // message + errno = 0; if (model_initialized && !init_exception_msg.empty()) { throw std::runtime_error( "Previous " + model_name + " init attempt had exception: \n\t" + init_exception_msg); @@ -135,7 +137,7 @@ namespace models { } else if (!utils::FileChecker::file_is_readable(bmi_init_config)) { init_exception_msg = "Cannot initialize " + model_name + " using unreadable file '" - + bmi_init_config + "'"; + + bmi_init_config + "'. Error: "+std::strerror(errno);; throw std::runtime_error(init_exception_msg); } else { diff --git a/include/core/DomainLayer.hpp b/include/core/DomainLayer.hpp new file mode 100644 index 0000000000..a98d3de3d8 --- /dev/null +++ b/include/core/DomainLayer.hpp @@ -0,0 +1,50 @@ +#ifndef __NGEN_DOMAIN_LAYER__ +#define __NGEN_DOMAIN_LAYER__ + +#include "Catchment_Formulation.hpp" +#include "Layer.hpp" + +namespace ngen +{ + class DomainLayer : public Layer + { + public: + // DomainLayer() = delete; + DomainLayer( + const LayerDescription& desc, + const Simulation_Time& s_t, + feature_type& features, + long idx, + std::shared_ptr formulation): + //description(desc), features(features), output_time_index(0), simulation_time(s_t), formulation(formulation) + Layer(desc, s_t, features, idx), formulation(formulation) + { + formulation->write_output("Time Step,""Time,"+formulation->get_output_header_line(",")+"\n"); + } + + /*** + * @brief Run one simulation timestep for each model in this layer + */ + void update_models() override{ + formulation->get_response(output_time_index, simulation_time.get_output_interval_seconds()); + std::string current_timestamp = simulation_time.get_timestamp(output_time_index); + std::string output = std::to_string(output_time_index)+","+current_timestamp+","+ + formulation->get_output_line_for_timestep(output_time_index)+"\n"; + formulation->write_output(output); + ++output_time_index; + if ( output_time_index < simulation_time.get_total_output_times() ) + { + simulation_time.advance_timestep(); + } + + } + + private: + // LayerDescription description; + // Simulation_Time simulation_time; + // long output_time_index; + std::shared_ptr formulation; + }; +} + +#endif \ No newline at end of file diff --git a/include/core/HY_Features.hpp b/include/core/HY_Features.hpp index 7f1333c5fe..07ad407f6a 100644 --- a/include/core/HY_Features.hpp +++ b/include/core/HY_Features.hpp @@ -79,6 +79,7 @@ namespace hy_features { * @param id * @return std::shared_ptr */ + // std::shared_ptr catchment_at(std::string id) std::shared_ptr catchment_at(std::string id) { if( _catchments.find(id) != _catchments.end() ) diff --git a/include/core/Layer.hpp b/include/core/Layer.hpp index 5aed53ac28..92e8c1528e 100644 --- a/include/core/Layer.hpp +++ b/include/core/Layer.hpp @@ -40,6 +40,19 @@ namespace ngen } + Layer( + const LayerDescription& desc, + const Simulation_Time& s_t, + feature_type& f, + long idx) : + description(desc), + simulation_time(s_t), + features(f), + output_time_index(idx) + { + + } + virtual ~Layer() {} /*** diff --git a/include/core/LayerData.hpp b/include/core/LayerData.hpp index 9820402ba5..47f2208633 100644 --- a/include/core/LayerData.hpp +++ b/include/core/LayerData.hpp @@ -51,6 +51,17 @@ namespace ngen void put_layer(const LayerDescription& desc, int id) { + // check to see if this layer was allready defined + if ( exists(id) ) + { + std::string message = "A layer with id = "; + message += std::to_string(id); + message += " was defined more than once"; + + std::runtime_error r_error(message); + + throw r_error; + } stored_layers[id] = desc; keys.push_back(id); } diff --git a/include/core/catchment/HY_Catchment.hpp b/include/core/catchment/HY_Catchment.hpp index 5b0329ec32..d86ce5d162 100644 --- a/include/core/catchment/HY_Catchment.hpp +++ b/include/core/catchment/HY_Catchment.hpp @@ -7,6 +7,7 @@ #include "HY_CatchmentRealization.hpp" #include "HY_HydroFeature.hpp" +#include "Formulation.hpp" class HY_HydroNexus; @@ -21,6 +22,7 @@ class HY_Catchment : public HY_HydroFeature public: HY_Catchment(); + // HY_Catchment(std::string id, Nexuses inflows, Nexuses outflows, std::shared_ptr realization, long lyr = 0): HY_Catchment(std::string id, Nexuses inflows, Nexuses outflows, std::shared_ptr realization, long lyr = 0): id(std::move(id)), inflows(std::move(inflows)), @@ -34,6 +36,7 @@ class HY_Catchment : public HY_HydroFeature layer(lyr){} virtual ~HY_Catchment(); const Nexuses& get_outflow_nexuses(){ return outflows; } + // std::shared_ptr realization; std::shared_ptr realization; /*! \brief get the hydrofabric layer of this catchment diff --git a/include/core/catchment/HY_CatchmentRealization.hpp b/include/core/catchment/HY_CatchmentRealization.hpp index ac0962b484..45bafebcc2 100644 --- a/include/core/catchment/HY_CatchmentRealization.hpp +++ b/include/core/catchment/HY_CatchmentRealization.hpp @@ -5,6 +5,7 @@ #include #include #include "GenericDataProvider.hpp" +#include "Formulation.hpp" using std::shared_ptr; @@ -15,6 +16,7 @@ class HY_Catchment; typedef long time_step_t; //TODO template //TODO template +// class HY_CatchmentRealization: public realization::Formulation class HY_CatchmentRealization { public: diff --git a/include/core/network.hpp b/include/core/network.hpp index f0db312e92..dd4b637111 100644 --- a/include/core/network.hpp +++ b/include/core/network.hpp @@ -136,7 +136,48 @@ namespace network { * */ using IndexPair = std::pair< NetworkIndexT::const_iterator, NetworkIndexT::const_iterator>; +struct detect_loops : public boost::dfs_visitor<> +{ + using colormap = std::map; + colormap vertex_coloring; + + using edgeColorMap = std::map; + edgeColorMap edge_coloring; + using vertex_t = Graph::vertex_descriptor; + + template + void tree_edge(Edge e, const Graph& g) { + //std::cout << "tree_edge: " << boost::source(e, g) << " --> " << boost::target(e, g) << std::endl; + + //edgeVisited.push(e); + if (vertexVisited.empty()) { + vertexVisited.push(boost::source(e, g)); + } + vertexVisited.push(boost::target(e, g)); + } + template + void back_edge(Edge e, const Graph& g) { + std::cout << source(e, g) + << " -- " + << target(e, g) << "\n"; + std::cout<< get(boost::vertex_name, g)[ source(e, g) ] << + " -- " << get(boost::vertex_name, g)[ target(e,g) ]<<"\n"; + vertex_t v2; + std::cout << "Cycle end= " << boost::target(e, g) << std::endl; + while ( vertexVisited.top() != boost::target(e, g) ) + { + //std::cout << " Cycle middle=" << vertexVisited.top() << std::endl; + v2 = vertexVisited.top(); + vertexVisited.pop(); + } + std::cout << "Cycle starting= " << vertexVisited.top() << std::endl; + vertexVisited.push(v2); + cycles.push_back(std::make_pair(source(e,g), target(e,g))); + } + std::stack vertexVisited; + std::vector> cycles; +}; /** * @brief A lightweight, graph based index of hydrologic features. * diff --git a/include/forcing/DataProviderSelectors.hpp b/include/forcing/DataProviderSelectors.hpp index c0186b5ba8..e64cf99eaa 100644 --- a/include/forcing/DataProviderSelectors.hpp +++ b/include/forcing/DataProviderSelectors.hpp @@ -9,20 +9,15 @@ * */ -class CatchmentAggrDataSelector -{ - public: +class DataSelector{ - /** - * @brief Construct a new Catchment Aggregate Data Selector object with default values - * - */ - CatchmentAggrDataSelector() : + public: + + DataSelector(): variable_name(), init_time(0), duration_s(1), - output_units(), - id_str() + output_units() {} /** @@ -34,8 +29,7 @@ class CatchmentAggrDataSelector * @param units The units to output the result in * @param id the id of the associated catchment */ - CatchmentAggrDataSelector(std::string id, std::string var, time_t start, long dur, std::string units) : - id_str(id), + DataSelector(std::string var, time_t start, long dur, std::string units) : variable_name(var), init_time(start), duration_s(dur), @@ -98,6 +92,41 @@ class CatchmentAggrDataSelector */ void set_output_units(std::string units) { output_units = units; } + private: + + std::string variable_name; //!< The variable name that should be queried + time_t init_time; //!< The inital time to query the requested variable + long duration_s; //!< The duration of the query + std::string output_units; //!< required units for the result to be return in +}; + +class CatchmentAggrDataSelector: public DataSelector +{ + public: + + /** + * @brief Construct a new Catchment Aggregate Data Selector object with default values + * + */ + CatchmentAggrDataSelector() : + DataSelector(), + id_str() + {} + + /** + * @brief Construct a new Data Selector object with inital values + * + * @param var The variable to be accessed by the selector + * @param start THe start time for this selector + * @param dur The duration for this selector + * @param units The units to output the result in + * @param id the id of the associated catchment + */ + CatchmentAggrDataSelector(std::string id, std::string var, time_t start, long dur, std::string units) : + DataSelector(var, start, dur, units), + id_str(id) + {} + /** * @brief Get the id string for this NetCDF Data Selector * @@ -114,10 +143,6 @@ class CatchmentAggrDataSelector private: - std::string variable_name; //!< The variable name that should be queried - time_t init_time; //!< The inital time to query the requested variable - long duration_s; //!< The duration of the query - std::string output_units; //!< required units for the result to be return in std::string id_str; //< the catchment to access data for }; diff --git a/include/forcing/selectors/GridSelectors.hpp b/include/forcing/selectors/GridSelectors.hpp new file mode 100644 index 0000000000..0c114d235f --- /dev/null +++ b/include/forcing/selectors/GridSelectors.hpp @@ -0,0 +1,211 @@ +#ifndef NGEN_GRID_SELECTORS_HPP +#define NGEN_GRID_SELECTORS_HPP + +#include +#include + +/** + * @brief DataProvider selector for regular gridded data + * + */ + +class GridSelector : DataSelector +{ + public: + + /** + * @brief Construct a new Grid Selector object with default values + * + */ + GridSelector() : + DataSelector(), + bounds(), + indicies() + {} + + /** + * @brief Construct a new Data Selector object with inital values + * + * @param var The variable to be accessed by the selector + * @param start THe start time for this selector + * @param dur The duration for this selector + * @param units The units to output the result in + * @param id the id of the associated catchment + */ + GridSelector(std::string id, std::string var, time_t start, long dur, std::string units, std::vector bounds, std::vector indicies) : + DataSelector(var, start, dur, units), + bounds(bounds), + indicies(indicies) + {} + + //TODO write get_bounds() get_indicies functions + + private: + + std::string variable_name; //!< The variable name that should be queried + time_t init_time; //!< The inital time to query the requested variable + long duration_s; //!< The duration of the query + std::string output_units; //!< required units for the result to be return in + std::string id_str; //< the catchment to access data for + std::vector bounds; //!< Bounding box for selection from regular grid + std::vector indicies; //!< Indicies for selection from regular grid. + // Assumed flatten indices along z,y,x +}; + +/* FIXME defined in DataProviderSelectors.hpp, redefine in this unit??? +#if !defined(__APPLE__) && !defined(__llvm__) && defined(__GNUC__) && __GNUC__ < 8 +// Enable a workaround for old GCC versions (but not Clang or Apple's fake-GCC that's actually Clang) +#define NGEN_SELECTOR_CAST +#endif +*/ +/* FIXME is this needed for this selector? + #if defined(NGEN_SELECTOR_CAST) + operator const CatchmentAggrDataSelector&() const { return *this; } + #endif +*/ +/** + * @brief This a data selector intended for use with CSV data + * + */ + +#endif + + +/** + * @brief DataProvider selector for regular gridded data + * + */ + +class RegularGridSelector +{ + public: + + /** + * @brief Construct a new Grid Selector object with default values + * + */ + RegularGridSelector() : + variable_name(), + init_time(0), + duration_s(1), + output_units(), + bounds(), + indicies() + {} + + /** + * @brief Construct a new Data Selector object with inital values + * + * @param var The variable to be accessed by the selector + * @param start THe start time for this selector + * @param dur The duration for this selector + * @param units The units to output the result in + * @param id the id of the associated catchment + */ + CatchmentAggrDataSelector(std::string id, std::string var, time_t start, long dur, std::string units) : + id_str(id), + variable_name(var), + init_time(start), + duration_s(dur), + output_units(units) + {} + + /** + * @brief Get the variable name for this selector + * + * @return std::string + */ + std::string get_variable_name() const { return variable_name; } + + /** + * @brief Get the initial time for this selector + * + * @return time_t + */ + time_t get_init_time() const { return init_time; } + + /** + * @brief Get the duration in seconds that is requested by this selector + * + * @return long + */ + long get_duration_secs() const { return duration_s; } + + /** + * @brief Get the output units that is requested by this selector + * + * @return std::string + */ + std::string get_output_units() const { return output_units; } + + /** + * @brief Set the variable name for this selector + * + * @param var The name of the variable to access + */ + void set_variable_name(std::string var) { variable_name = var; } + + /** + * @brief Set the init time for this selector + * + * @param start The inital time for this query + */ + void set_init_time(time_t start) { init_time = start; } + + /** + * @brief Set the duration in secs for the query represented by this selector + * + * @param dur the duration of the query in seconds + */ + void set_duration_secs(long dur) { duration_s = dur; } + + /** + * @brief Set the output units for the data to be returned with this selector + * + * @param units The units of the output + */ + void set_output_units(std::string units) { output_units = units; } + + /** + * @brief Get the id string for this NetCDF Data Selector + * + * @return std::string + */ + std::string get_id() const { return id_str; } + + /** + * @brief Set the id string for this NetCDF Data Selector + * + * @param s + */ + void set_id(std::string s) { id_str = s; } + + private: + + std::string variable_name; //!< The variable name that should be queried + time_t init_time; //!< The inital time to query the requested variable + long duration_s; //!< The duration of the query + std::string output_units; //!< required units for the result to be return in + std::string id_str; //< the catchment to access data for + std::vector bounds; //!< Bounding box for selection from regular grid + std::vector indicies; //!< Indicies for selection from regular grid. + // Assumed flatten indices along z,y,x +}; + +/* FIXME defined in DataProviderSelectors.hpp, redefine in this unit??? +#if !defined(__APPLE__) && !defined(__llvm__) && defined(__GNUC__) && __GNUC__ < 8 +// Enable a workaround for old GCC versions (but not Clang or Apple's fake-GCC that's actually Clang) +#define NGEN_SELECTOR_CAST +#endif +*/ +/* FIXME is this needed for this selector? + #if defined(NGEN_SELECTOR_CAST) + operator const CatchmentAggrDataSelector&() const { return *this; } + #endif +*/ +/** + * @brief This a data selector intended for use with CSV data + * + */ + +#endif diff --git a/include/geojson/JSONProperty.hpp b/include/geojson/JSONProperty.hpp index e5a7ab69ba..43fda70509 100644 --- a/include/geojson/JSONProperty.hpp +++ b/include/geojson/JSONProperty.hpp @@ -5,7 +5,7 @@ #include #include #include - +#include #include #include #include @@ -388,6 +388,10 @@ namespace geojson { } } + JSONProperty(const std::string& value_key, const JSONProperty&original):JSONProperty(original){ + key = value_key; + } + /** * A basic destructor */ @@ -419,6 +423,46 @@ namespace geojson { data = Object( &values ); } + static void print_property(const geojson::JSONProperty& p, int tab=0){ + + switch( p.get_type() ){ + case geojson::PropertyType::String: + std::cout< +#include +#include "Formulation.hpp" +#include "GenericDataProvider.hpp" +#include "GM_Object.hpp" +#include "FileStreamHandler.hpp" + +namespace realization { + + class Layer_Formulation : public Formulation { + public: + Layer_Formulation(std::string id, std::shared_ptr forcing, utils::StreamHandler output_stream) + : Formulation(id), forcing(forcing), output(output_stream) { + + }; + + Layer_Formulation(std::string id) : Formulation(id){ + + }; + + /** + * Get a header line appropriate for a file made up of entries from this type's implementation of + * ``get_output_line_for_timestep``. + * + * Note that like the output generating function, this line does not include anything for time step. + * + * A default implementation is provided for inheritors of this type, which includes only "Total Discharge." + * + * @return An appropriate header line for this type. + */ + std::string get_output_header_line(std::string delimiter) override { + return "Total Discharge"; + } + + /** + * Execute the backing model formulation for the given time step, where it is of the specified size, and + * return the response output. + * + * Any inputs and additional parameters must be made available as instance members. + * + * Types should clearly document the details of their particular response output. + * + * @param t_index The index of the time step for which to run model calculations. + * @param d_delta_s The duration, in seconds, of the time step for which to run model calculations. + * @return The response output of the model for this time step. + */ + double get_response(time_step_t t_index, time_step_t t_delta) override = 0; + + const std::vector& get_required_parameters() override = 0; + + void create_formulation(boost::property_tree::ptree &config, geojson::PropertyMap *global = nullptr) override = 0; + void create_formulation(geojson::PropertyMap properties) override = 0; + virtual ~Layer_Formulation(){}; + + protected: + polygon_t bounds; + std::shared_ptr forcing; + utils::StreamHandler output; + + private: + std::string cat_id; + }; +} +#endif // LAYER_FORMULATION_H diff --git a/include/realizations/catchment/Bmi_Module_Formulation.hpp b/include/realizations/catchment/Bmi_Module_Formulation.hpp index a9b63a7866..d7de0e9a9d 100644 --- a/include/realizations/catchment/Bmi_Module_Formulation.hpp +++ b/include/realizations/catchment/Bmi_Module_Formulation.hpp @@ -898,7 +898,7 @@ namespace realization { if(values.size() == 1){ //FIXME this isn't generic broadcasting, but works for scalar implementations #ifndef NGEN_QUIET - std::cerr << "WARN: broadcasting variable '" << var_name << "' from scalar to expected array\n"; + //std::cerr << "WARN: broadcasting variable '" << var_name << "' from scalar to expected array\n"; #endif values.resize(numItems, values[0]); } else if (values.size() != numItems) { diff --git a/include/realizations/catchment/Catchment_Formulation.hpp b/include/realizations/catchment/Catchment_Formulation.hpp index cce4a4cf2b..ccc048ee15 100644 --- a/include/realizations/catchment/Catchment_Formulation.hpp +++ b/include/realizations/catchment/Catchment_Formulation.hpp @@ -9,9 +9,11 @@ namespace realization { + // class Catchment_Formulation : public Formulation { class Catchment_Formulation : public Formulation, public HY_CatchmentArea { public: Catchment_Formulation(std::string id, std::shared_ptr forcing, utils::StreamHandler output_stream) + // : Formulation(id) { : Formulation(id), HY_CatchmentArea(forcing, output_stream) { // Assume the catchment ID is equal to or embedded in the formulation `id` size_t idx = id.find("."); @@ -92,11 +94,11 @@ namespace realization { virtual ~Catchment_Formulation(){}; protected: - std::string get_catchment_id() override { + std::string get_catchment_id() { return this->cat_id; } - void set_catchment_id(std::string cat_id) override { + void set_catchment_id(std::string cat_id) { this->cat_id = cat_id; } diff --git a/include/realizations/catchment/Formulation.hpp b/include/realizations/catchment/Formulation.hpp index 161252e16f..6395f6de23 100644 --- a/include/realizations/catchment/Formulation.hpp +++ b/include/realizations/catchment/Formulation.hpp @@ -8,7 +8,9 @@ #include #include "JSONProperty.hpp" - +#include "StreamHandler.hpp" +#include "FileStreamHandler.hpp" +#include "GenericDataProvider.hpp" #include #include @@ -95,6 +97,8 @@ namespace realization { virtual void create_formulation(boost::property_tree::ptree &config, geojson::PropertyMap *global = nullptr) = 0; virtual void create_formulation(geojson::PropertyMap properties) = 0; + //void set_output_stream(std::string file_path){output = utils::FileStreamHandler(file_path.c_str());} + //void write_output(std::string out){ output<& get_required_parameters() = 0; @@ -145,6 +149,8 @@ namespace realization { } std::string id; + // std::shared_ptr forcing; + // utils::StreamHandler output; }; } diff --git a/include/realizations/catchment/Formulation_Constructors.hpp b/include/realizations/catchment/Formulation_Constructors.hpp index a5cf348f04..0eb22ec3f2 100644 --- a/include/realizations/catchment/Formulation_Constructors.hpp +++ b/include/realizations/catchment/Formulation_Constructors.hpp @@ -45,6 +45,14 @@ namespace realization { #endif // ACTIVATE_PYTHON }; + static std::string valid_formulation_keys(){ + std::string keys = ""; + for(const auto& kv : formulations){ + keys.append(kv.first+" "); + } + return keys; + } + static bool formulation_exists(std::string formulation_type) { return formulations.count(formulation_type) > 0; } @@ -55,6 +63,7 @@ namespace realization { forcing_params &forcing_config, utils::StreamHandler output_stream ) { + std::cout<<"FORMULATION: "< fp; if (forcing_config.provider == "CsvPerFeature" || forcing_config.provider == ""){ @@ -83,13 +92,16 @@ namespace realization { return node.first; } }*/ + std::cout<<"Getting optional name\n"; boost::optional key = tree.get_optional("name"); + std::cout<<"Got name??? "< #include -#include +#include //TODO check if needed #include #include -#include "JSONProperty.hpp" +#include "JSONProperty.hpp" //TODO check if needed here #include "features/Features.hpp" #include #include "Formulation_Constructors.hpp" -#include "Simulation_Time.hpp" -#include "routing/Routing_Params.h" -#include "LayerData.hpp" +#include "Simulation_Time.hpp" //TODO check if needed here +#include "routing/Routing_Params.h" //TODO check if needed here +#include "LayerData.hpp" //TODO check if needed here +#include "../config.hpp" + +//TODO consider a ngen::forcing namespace +class ProviderManager{ + public: + ProviderManager(){}; + + private: + +}; namespace realization { class Formulation_Manager { @@ -55,79 +65,24 @@ namespace realization { //TODO seperate the parsing of configuration options like time //and routing and other non feature specific tasks from this main function //which has to iterate the entire hydrofabric. + auto possible_global_config = tree.get_child_optional("global"); if (possible_global_config) { - this->global_formulation_tree = *possible_global_config; - - //get forcing info - for (auto &forcing_parameter : (*possible_global_config).get_child("forcing")) { - this->global_forcing.emplace( - forcing_parameter.first, - geojson::JSONProperty(forcing_parameter.first, forcing_parameter.second) - ); - } - - //get first empty key under formulations (corresponds to first json array element) - auto formulation = (*possible_global_config).get_child("formulations.."); - - for (std::pair global_setting : formulation.get_child("params")) { - this->global_formulation_parameters.emplace( - global_setting.first, - geojson::JSONProperty(global_setting.first, global_setting.second) - ); - } + global_config = realization::config::Formulation_Config(*possible_global_config); } - /** - * Read simulation time from configuration file - * /// \todo TODO: Separate input_interval from output_interval - */ auto possible_simulation_time = tree.get_child_optional("time"); if (!possible_simulation_time) { throw std::runtime_error("ERROR: No simulation time period defined."); } - - geojson::JSONProperty simulation_time_parameters("time", *possible_simulation_time); - - std::vector missing_simulation_time_parameters; - - if (!simulation_time_parameters.has_key("start_time")) { - missing_simulation_time_parameters.push_back("start_time"); - } - - if (!simulation_time_parameters.has_key("end_time")) { - missing_simulation_time_parameters.push_back("end_time"); - } - - if (!simulation_time_parameters.has_key("output_interval")) { - missing_simulation_time_parameters.push_back("output_interval"); - } - - if (missing_simulation_time_parameters.size() > 0) { - std::string message = "ERROR: A simulation time parameter cannot be created; the following parameters are missing: "; - - for (int missing_parameter_index = 0; missing_parameter_index < missing_simulation_time_parameters.size(); missing_parameter_index++) { - message += missing_simulation_time_parameters[missing_parameter_index]; - - if (missing_parameter_index < missing_simulation_time_parameters.size() - 1) { - message += ", "; - } - } - - throw std::runtime_error(message); - } - - simulation_time_params simulation_time_config( - simulation_time_parameters.at("start_time").as_string(), - simulation_time_parameters.at("end_time").as_string(), - simulation_time_parameters.at("output_interval").as_natural_number() - ); - + config::Time time = config::Time(*possible_simulation_time); + auto simulation_time_config = time.make_params(); /** * Call constructor to construct a Simulation_Time object - */ + */ + this->Simulation_Time_Object = std::make_shared(simulation_time_config); /** @@ -136,56 +91,62 @@ namespace realization { // try to get the json node auto layers_json_array = tree.get_child_optional("layers"); - + config::Layer layer; + // layer description struct + ngen::LayerDescription layer_desc; + layer_desc = layer.get_descriptor(); + // add the default surface layer to storage + layer_storage.put_layer(layer_desc, layer_desc.id); // check to see if the node existed - if (!layers_json_array) { - // layer description struct - ngen::LayerDescription layer_desc; - - // extract and store layer data from the json - layer_desc.name = "surface layer"; - layer_desc.id = 0; - layer_desc.time_step = 3600; - layer_desc.time_step_units = "s"; - - // add the layer to storage - layer_storage.put_layer(layer_desc, layer_desc.id); - } - else - { + // if (!layers_json_array) { + // layer_desc = layer.get_descriptor(); + // // add the layer to storage + // layer_storage.put_layer(layer_desc, layer_desc.id); + // } + // else + if(layers_json_array){ + std::cout<<"Parsing Layers as Array???\n"; + for (std::pair layer_config : *layers_json_array) { - // layer description struct - ngen::LayerDescription layer_desc; - - // extract and store layer data from the json - layer_desc.name = layer_config.second.get("name"); - layer_desc.id = layer_config.second.get("id"); - layer_desc.time_step = layer_config.second.get("time_step"); - boost::optional layer_units = layer_config.second.get_optional("time_step_units"); - if (*layer_units == "") layer_units = "s"; - layer_desc.time_step_units = *layer_units; - - // check to see if this layer was allready defined - if (layer_storage.exists(layer_desc.id) ) - { - std::string message = "A layer with id = "; - message += std::to_string(layer_desc.id); - message += " was defined more than once"; - - std::runtime_error r_error(message); - - throw r_error; - } - + layer = config::Layer(layer_config.second); + layer_desc = layer.get_descriptor(); + std::cout<<"key: "<set_output_stream(get_output_root() + layer_desc.name + "_layer_"+std::to_string(layer_desc.id) + ".csv"); + domain_formulations.emplace( + layer_desc.id, + construct_formulation_from_tree(simulation_time_config, + "layer-"+std::to_string(layer_desc.id), + layer_config.second, //FIXME this isn't even used in the func anymore + layer.formulation, + output_stream + ) + ); + domain_formulations.at(layer_desc.id)->set_output_stream(get_output_root() + layer_desc.name + "_layer_"+std::to_string(layer_desc.id) + ".csv"); + } + + //TODO for each layer, create deferred providers for use by other layers + //VERY SIMILAR TO NESTED MODULE INIT } } + //TODO use the set of layer providers as input for catchments to lookup from + /** * Read routing configurations from configuration file */ @@ -195,11 +156,7 @@ namespace realization { //Since it is possible to build NGEN without routing support, if we see it in the config //but it isn't enabled in the build, we should at least put up a warning #ifdef NGEN_ROUTING_ACTIVE - geojson::JSONProperty routing_parameters("routing", *possible_routing_configs); - - this->routing_config = std::make_shared( - routing_parameters.at("t_route_config_file_with_path").as_string() - ); + this->routing_config = (config::Routing(*possible_routing_configs)).params; using_routing = true; #else using_routing = false; @@ -207,7 +164,7 @@ namespace realization { <<", but routing support isn't enabled. No routing will occur."<get_feature(catchment_index); - for (auto &formulation: *formulations) { + // for (auto &formulation: catchment_formulation.formulation_tree) { // Handle single-bmi - decltype(auto) model_params = formulation.second.get_child_optional("params.model_params"); - if (model_params) { - parse_external_model_params(*model_params, catchment_feature); - } - + // decltype(auto) model_params = formulation.second.get_child_optional("params.model_params"); + // if (model_params) { + // parse_external_model_params(*model_params, catchment_feature); + // } + // if( catchment_formulation.formulation.parameters.count("model_params") > 0){ + // parse_external_model_params(catchment_formulation.formulation.parameters, catchment_feature); + // } // Handle multi-bmi // FIXME: this will not handle doubly nested multi-BMI configs, // might need a recursive helper here? - decltype(auto) nested_modules = formulation.second.get_child_optional("params.modules"); - if (nested_modules) { - for (decltype(auto) nested_formulation : *nested_modules) { - decltype(auto) nested_model_params = nested_formulation.second.get_child_optional("params.model_params"); - if (nested_model_params) { - parse_external_model_params(*nested_model_params, catchment_feature); - } - } - } - + // decltype(auto) nested_modules = formulation.second.get_child_optional("params.modules"); + // if (nested_modules) { + // for (decltype(auto) nested_formulation : *nested_modules) { + // decltype(auto) nested_model_params = nested_formulation.second.get_child_optional("params.model_params"); + // if (nested_model_params) { + // parse_external_model_params(*nested_model_params, catchment_feature); + // } + // } + // } + //for(auto& formulation : catchment_formulation.formulation.nested){ + //auto params = formulation.get_values(); + // parse_external_model_params(formulation.parameters, catchment_feature); + //} + std::cout<<"READ: Link External\n"; + catchment_formulation.formulation.link_external(catchment_feature); this->add_formulation( this->construct_formulation_from_tree( simulation_time_config, catchment_config.first, catchment_config.second, - formulation.second, + catchment_formulation, output_stream ) ); - break; //only construct one for now FIXME - } //end for formulaitons + // break; //only construct one for now FIXME + // } //end for formulaitons }//end for catchments @@ -288,6 +251,14 @@ namespace realization { return this->formulations.at(id); } + virtual std::shared_ptr get_domain_formulation(long id) const { + return this->domain_formulations.at(id); + } + + virtual bool has_domain_formulation(long id) const { + return this->domain_formulations.count( id ) > 0; + } + virtual bool contains(std::string identifier) const { return this->formulations.count(identifier) > 0; } @@ -371,31 +342,27 @@ namespace realization { simulation_time_params &simulation_time_config, std::string identifier, boost::property_tree::ptree &tree, - const boost::property_tree::ptree &formulation, + //const boost::property_tree::ptree &formulation, + const realization::config::Formulation_Config& catchment_formulation, utils::StreamHandler output_stream ) { - auto params = formulation.get_child("params"); - std::string formulation_type_key; - try { - formulation_type_key = get_formulation_key(formulation); - } - catch(std::exception& e) { - throw std::runtime_error("Catchment " + identifier + " failed initialization: " + e.what()); + //Formulation_Config catchment_formulation(tree); + std::cout<<"FS: "< missing_parameters; - if (!forcing_parameters.has_key("path")) { + if (!catchment_formulation.forcing.has_key("path")) { missing_parameters.push_back("path"); } @@ -413,71 +380,93 @@ namespace realization { throw std::runtime_error(message); } - geojson::PropertyMap local_forcing; - for (auto &forcing_parameter : *possible_forcing) { - local_forcing.emplace( - forcing_parameter.first, - geojson::JSONProperty(forcing_parameter.first, forcing_parameter.second) - ); - } - - forcing_params forcing_config = this->get_forcing_params(local_forcing, identifier, simulation_time_config); - - std::shared_ptr constructed_formulation = construct_formulation(formulation_type_key, identifier, forcing_config, output_stream); + forcing_params forcing_config = this->get_forcing_params(catchment_formulation.forcing.parameters, identifier, simulation_time_config); + std::cout<<"HERERE: "< constructed_formulation = construct_formulation(catchment_formulation.formulation.type, identifier, forcing_config, output_stream); //, geometry); - constructed_formulation->create_formulation(formulation_config, &global_formulation_parameters); + + constructed_formulation->create_formulation(catchment_formulation.formulation.parameters); + std::cout<<"CREATED FROM TREE\n"; return constructed_formulation; } std::shared_ptr construct_missing_formulation(geojson::Feature& feature, utils::StreamHandler output_stream, simulation_time_params &simulation_time_config){ const std::string identifier = feature->get_id(); - std::string formulation_type_key = get_formulation_key(global_formulation_tree.get_child("formulations..")); - - forcing_params forcing_config = this->get_forcing_params(this->global_forcing, identifier, simulation_time_config); - - std::shared_ptr missing_formulation = construct_formulation(formulation_type_key, identifier, forcing_config, output_stream); + //std::string formulation_type_key = get_formulation_key(global_config.formulation_tree); + std::cout<<"GETTING GLOBAL FORCING PARAMS\n"; + auto fubar = global_config.formulation_tree.get_child("forcing");//(global_config.forcing.parameters) + geojson::PropertyMap parameters; + for (auto &forcing_parameter : fubar) { + parameters.emplace( + forcing_parameter.first, + geojson::JSONProperty(forcing_parameter.first, forcing_parameter.second) + ); + } + forcing_params forcing_config = this->get_forcing_params(parameters, identifier, simulation_time_config); + std::cout<<"DOING STUFF\n"; + std::shared_ptr missing_formulation = construct_formulation(global_config.formulation.type, identifier, forcing_config, output_stream); // Need to work with a copy, since it is altered in-place - geojson::PropertyMap global_properties_copy = global_formulation_parameters; - Catchment_Formulation::config_pattern_substitution(global_properties_copy, + //geojson::PropertyMap global_properties_copy = global_config.formulation.parameters; + + // FIXME something isn't getting parsed here correct run the nom_Grid example realization and you can see + // only sloth is created/called!!!!!! + + realization::config::Formulation_Config global_copy = global_config; + Catchment_Formulation::config_pattern_substitution(global_copy.formulation.parameters, BMI_REALIZATION_CFG_PARAM_REQ__INIT_CONFIG, "{{id}}", identifier); - + //for(auto& formulation : global_copy.formulation.nested){ + // Catchment_Formulation::config_pattern_substitution(formulation.parameters, + // BMI_REALIZATION_CFG_PARAM_REQ__INIT_CONFIG, "{{id}}", + // identifier); + //} // parse any external model parameters in this config - // handle single-bmi - if (global_properties_copy.count("model_params") > 0) { - decltype(auto) model_params = global_properties_copy.at("model_params"); - geojson::PropertyMap model_params_copy = model_params.get_values(); - parse_external_model_params(model_params_copy, feature); - global_properties_copy.at("model_params") = geojson::JSONProperty("model_params", model_params_copy); - } + // // handle single-bmi + // if (global_properties_copy.count("model_params") > 0) { + // decltype(auto) model_params = global_properties_copy.at("model_params"); + // geojson::PropertyMap model_params_copy = model_params.get_values(); + // parse_external_model_params(model_params_copy, feature); + // global_properties_copy.at("model_params") = geojson::JSONProperty("model_params", model_params_copy); + // } + + // // handle multi-bmi + // // FIXME: this seems inefficient -- is there a better way? + // if (global_properties_copy.count("modules") > 0) { + // decltype(auto) nested_modules = global_properties_copy.at("modules").as_list(); + // for (auto& bmi_module : nested_modules) { + // geojson::PropertyMap module_def = bmi_module.get_values(); + // geojson::PropertyMap module_params = module_def.at("params").get_values(); + // if (module_params.count("model_params") > 0) { + // decltype(auto) model_params = module_params.at("model_params"); + // geojson::PropertyMap model_params_copy = model_params.get_values(); + // parse_external_model_params(model_params_copy, feature); + // module_params.at("model_params") = geojson::JSONProperty("model_params", model_params_copy); + // } + // module_def.at("params") = geojson::JSONProperty("params", module_params); + // bmi_module = geojson::JSONProperty("", module_def); + // } + // global_properties_copy.at("modules") = geojson::JSONProperty("", nested_modules); + // } + + // if( global_copy.formulation.parameters.count("model_params") > 0){ + // parse_external_model_params(global_copy.formulation.parameters, feature); + // } + // for(auto& formulation : global_copy.formulation.nested){ + // //auto params = formulation.get_values(); + // parse_external_model_params(formulation.parameters, feature); + // } + // geojson::JSONProperty::print_property(global_config.formulation.parameters.at("modules")); + global_config.formulation.link_external(feature); + // geojson::JSONProperty::print_property(global_config.formulation.parameters.at("modules")); + missing_formulation->create_formulation(global_config.formulation.parameters); - // handle multi-bmi - // FIXME: this seems inefficient -- is there a better way? - if (global_properties_copy.count("modules") > 0) { - decltype(auto) nested_modules = global_properties_copy.at("modules").as_list(); - for (auto& bmi_module : nested_modules) { - geojson::PropertyMap module_def = bmi_module.get_values(); - geojson::PropertyMap module_params = module_def.at("params").get_values(); - if (module_params.count("model_params") > 0) { - decltype(auto) model_params = module_params.at("model_params"); - geojson::PropertyMap model_params_copy = model_params.get_values(); - parse_external_model_params(model_params_copy, feature); - module_params.at("model_params") = geojson::JSONProperty("model_params", model_params_copy); - } - module_def.at("params") = geojson::JSONProperty("params", module_params); - bmi_module = geojson::JSONProperty("", module_def); - } - global_properties_copy.at("modules") = geojson::JSONProperty("", nested_modules); - } - - missing_formulation->create_formulation(global_properties_copy); return missing_formulation; } - forcing_params get_forcing_params(geojson::PropertyMap &forcing_prop_map, std::string identifier, simulation_time_params &simulation_time_config) { - std::string path; + forcing_params get_forcing_params(const geojson::PropertyMap &forcing_prop_map, std::string identifier, simulation_time_params &simulation_time_config) { + std::string path = ""; if(forcing_prop_map.count("path") != 0){ path = forcing_prop_map.at("path").as_string(); } @@ -505,7 +494,7 @@ namespace realization { } std::string filepattern = forcing_prop_map.at("file_pattern").as_string(); - + std::cout<<"FSFSFS: "<> formulations; + //Store global layer formulation pointers + std::map > domain_formulations; + std::shared_ptr routing_config; bool using_routing = false; diff --git a/include/realizations/config.hpp b/include/realizations/config.hpp new file mode 100644 index 0000000000..9e09b92d45 --- /dev/null +++ b/include/realizations/config.hpp @@ -0,0 +1,270 @@ +#ifndef NGEN_REALIZATION_CONFIG_H +#define NGEN_REALIZATION_CONFIG_H + +#include +#include + +#include "JSONProperty.hpp" +#include "Simulation_Time.hpp" +#include "routing/Routing_Params.h" +#include "LayerData.hpp" + + +#include + +namespace realization{ + namespace config{ + + static const std::string ROUTING_CONFIG_KEY = "t_route_config_file_with_path"; + struct Routing{ + std::shared_ptr params; + Routing(const boost::property_tree::ptree& tree){ + params = std::make_shared(tree.get(ROUTING_CONFIG_KEY, "")); + } + }; + + struct Time{ + std::string start_time; + std::string end_time; + unsigned int output_interval; + + Time(const boost::property_tree::ptree& tree){ + start_time = tree.get("start_time", std::string()); + end_time = tree.get("end_time", std::string()); + output_interval = tree.get("output_interval", 0); + } + + simulation_time_params make_params(){ + std::vector missing_simulation_time_parameters; + + if (start_time.empty()){ + missing_simulation_time_parameters.push_back("start_time"); + } + + if (end_time.empty()) { + missing_simulation_time_parameters.push_back("end_time"); + } + + if (output_interval == 0) { + missing_simulation_time_parameters.push_back("output_interval"); + } + + if (missing_simulation_time_parameters.size() > 0) { + std::string message = "ERROR: A simulation time parameter cannot be created; the following parameters are missing or invalid: "; + + for (int missing_parameter_index = 0; missing_parameter_index < missing_simulation_time_parameters.size(); missing_parameter_index++) { + message += missing_simulation_time_parameters[missing_parameter_index]; + + if (missing_parameter_index < missing_simulation_time_parameters.size() - 1) { + message += ", "; + } + } + + throw std::runtime_error(message); + } + return simulation_time_params( + start_time, + end_time, + output_interval + ); + } + }; + + + struct Formulation_Config{ + + struct Formulation{ + std::string type; + geojson::PropertyMap parameters; + std::vector nested; + Formulation():type(std::string()), parameters(geojson::PropertyMap()){} + Formulation(std::string& type, geojson::PropertyMap params):type(std::move(type)), parameters(params){} + Formulation(const boost::property_tree::ptree& tree){ + type = tree.get("name"); + for (std::pair setting : tree.get_child("params")) { + //if( setting.first == "modules" ) continue; //FUCK STUFF if skipping modules, we somehow + //end up with "modules" as an object and it fucks up calls to as_list() later with a boost assertion error + //that has something to do with an empty invariant... + //std::cout<<"FORMULATION: KEY -> "< tmp; + for(auto& n : nested){ + if(n.parameters.count("model_params")){ + n.link_external(feature); + } + geojson::PropertyMap map = {}; + map.emplace("name", geojson::JSONProperty("name", n.type)); + map.emplace("params", geojson::JSONProperty("", n.parameters)); + // std::cout<<"TESTING "<has_property(param_name)) { + auto catchment_attribute = feature->get_property(param_name); + + // Use param.first in the `.emplace` calls instead of param_name since + // the expected name is given by the key of the model_params values. + switch (catchment_attribute.get_type()) { + case geojson::PropertyType::List: + case geojson::PropertyType::Object: + // TODO: Should list/object values be passed to model parameters? + // Typically, feature properties *should* be scalars. + std::cerr << "WARNING: property type " << static_cast(catchment_attribute.get_type()) << " not allowed as model parameter. " + << "Must be one of: Natural (int), Real (double), Boolean, or String" << '\n'; + break; + default: + //std::cout<<"NO CLUE: "<parameters.emplace( + forcing_parameter.first, + geojson::JSONProperty(forcing_parameter.first, forcing_parameter.second) + ); + } + } + bool has_key(const std::string& key) const{ + return parameters.count(key) > 0; + } + }; + + Formulation_Config(){}; + Formulation_Config(const boost::property_tree::ptree& tree):formulation_tree(tree){ + + auto possible_forcing = tree.get_child_optional("forcing"); + + if (possible_forcing) { + forcing = Forcing(*possible_forcing); + } + //get first empty key under formulations (corresponds to first json array element) + auto possible_formulation_tree = tree.get_child_optional("formulations.."); + if(possible_formulation_tree){ + formulation = Formulation(*possible_formulation_tree); + } + } + + bool has_formulation(){ + return !(formulation.type.empty() || formulation.parameters.empty()); + } + + Formulation formulation; + boost::property_tree::ptree formulation_tree; + Forcing forcing; + }; + + struct Layer{ + Layer():descriptor( {"surface layer", "s", 0, 3600 } ){}; + Layer(const boost::property_tree::ptree& tree):formulation(tree){ + auto name = tree.get_optional("name"); + if(!name) missing_keys.push_back("name"); + auto unit = tree.get("time_step_units", "s"); + + auto id = tree.get_optional("id"); + if(!id) missing_keys.push_back("id"); + + auto ts = tree.get_optional("time_step"); + //TODO is time_step required? e.g. need to add to missing_keys? + auto tmp = tree.get_optional("domain"); + if(tmp) domain = *tmp; + if(missing_keys.empty()){ + descriptor = {*name, unit, *id, *ts};// ngen::LayerDescription( ); + } + + std::cout<<"BUILDING LAYER FORMULATION: "< missing_keys; + + }; + };//end namespace config +}//end namespace realization +#endif //NGEN_REALIZATION_CONFIG_H \ No newline at end of file diff --git a/include/utilities/FileStreamHandler.hpp b/include/utilities/FileStreamHandler.hpp index 1e4108fbaa..1cdff9e064 100644 --- a/include/utilities/FileStreamHandler.hpp +++ b/include/utilities/FileStreamHandler.hpp @@ -17,6 +17,34 @@ namespace utils virtual ~FileStreamHandler(){} }; + // class FileStreamHandler : public StreamHandler + // { + // //FIXME!!!!!!!!! hack to prevent too many open files---also makes the output's empty :confused: + // public: + // FileStreamHandler(const char* path) : StreamHandler(), path(path) + // { + // stream = std::make_shared(); + // stream->open(path, std::ios::trunc); //clear any existing data + // stream->close(); //avoid too many open files + // output_stream = stream; + // } + // virtual ~FileStreamHandler(){ + // if ( stream->is_open() ) stream->close(); + // } + + // template std::ostream& operator<<(const DataType& val) + // { + // if( !stream->is_open() ) + // stream->open(path, std::ios_base::app); + // put(val); + // stream->close(); + // return *output_stream; + // } + // private: + // std::string path; + // std::shared_ptr stream; + // }; + } diff --git a/src/NGen.cpp b/src/NGen.cpp index 55a8f6e472..dbffbc2585 100644 --- a/src/NGen.cpp +++ b/src/NGen.cpp @@ -55,6 +55,7 @@ int mpi_num_procs; #include #include +#include std::unordered_map nexus_outfiles; @@ -440,9 +441,7 @@ int main(int argc, char *argv[]) { // get the keys for the existing layers std::vector& keys = layer_meta_data.get_keys(); - // get the converted time steps for layers - double min_time_step; - double max_time_step; + //FIXME refactor the layer building to avoid this mess std::vector time_steps; for(int i = 0; i < keys.size(); ++i) { @@ -451,10 +450,6 @@ int main(int argc, char *argv[]) { time_steps.push_back(c_value); } - std::vector output_time_index; - output_time_index.resize(keys.size()); - std::fill(output_time_index.begin(), output_time_index.end(),0); - // now create the layer objects // first make sure that the layer are listed in decreasing order @@ -470,15 +465,25 @@ int main(int argc, char *argv[]) { // make a new simulation time object with a different output interval Simulation_Time sim_time(*manager->Simulation_Time_Object, time_steps[i]); - - for ( std::string id : features.catchments(keys[i]) ) { cat_ids.push_back(id); } - if (keys[i] != 0 ) - { - layers[i] = std::make_shared(desc, cat_ids, sim_time, features, catchment_collection, 0); + std::cout<<"FUBAR: "<has_domain_formulation(keys[i])){ + //create a domain wide layer + auto formulation = manager->get_domain_formulation(keys[i]); + layers[i] = std::make_shared(desc, sim_time, features, 0, formulation); + std::cout<<"CREATING DOMAIN FORMULATION: "<(desc, cat_ids, sim_time, features, catchment_collection, 0, nexus_subset_ids, nexus_outfiles); + else{ + std::cout<<"FUCKING OTHER LAYER: "<(desc, cat_ids, sim_time, features, catchment_collection, 0); + } + else + { + layers[i] = std::make_shared(desc, cat_ids, sim_time, features, catchment_collection, 0, nexus_subset_ids, nexus_outfiles); + } } } @@ -513,7 +518,8 @@ int main(int argc, char *argv[]) { // only advance if you would not pass the master next time and the previous layer next time if ( layer_next_time <= next_time && layer_next_time <= prev_layer_time) { - layer->update_models(); + if(count%100==0) std::cout<<"Updating layer: "<get_name()<<"\n"; + layer->update_models(); //assume update_models() calls time->advance_timestep() prev_layer_time = layer_next_time; } else diff --git a/src/core/SurfaceLayer.cpp b/src/core/SurfaceLayer.cpp index 3eec688fa8..e044d55292 100644 --- a/src/core/SurfaceLayer.cpp +++ b/src/core/SurfaceLayer.cpp @@ -31,7 +31,8 @@ void ngen::SurfaceLayer::update_models() } else { //This is a terminal node, SHOULDN'T be remote, so ID shouldn't matter too much - cat_id = "terminal"; + cat_id = "terminal-fs"; + //continue; } //std::cerr << "Requesting water from nexus, id = " << id << " at time = " < forcing) :Formulation("NaN"), forcing(forcing) { } HY_CatchmentRealization::HY_CatchmentRealization(std::shared_ptr forcing) : forcing(forcing) { } HY_CatchmentRealization::~HY_CatchmentRealization() diff --git a/src/core/network.cpp b/src/core/network.cpp index 38580ef600..2c8ea43606 100644 --- a/src/core/network.cpp +++ b/src/core/network.cpp @@ -3,6 +3,8 @@ #include #include #include +#include +#include using namespace network; @@ -95,7 +97,12 @@ void Network::init_indicies(){ //std::cout<<"TW: "<<*it<graph, vis, make_assoc_property_map(vis.vertex_coloring), make_assoc_property_map(vis.edge_coloring)); + for(auto p : vis.cycles){ + std::cout<graph, std::back_inserter(this->topo_order), boost::vertex_index_map(get(boost::vertex_index, this->graph))); } diff --git a/src/core/nexus/HY_PointHydroNexus.cpp b/src/core/nexus/HY_PointHydroNexus.cpp index 3497e3b19b..7c6e70ec51 100644 --- a/src/core/nexus/HY_PointHydroNexus.cpp +++ b/src/core/nexus/HY_PointHydroNexus.cpp @@ -62,7 +62,11 @@ double HY_PointHydroNexus::get_downstream_flow(std::string catchment_id, time_st { // there are no recorded flows for this time. // throw exception - + std::cout<<"Cat id: "< JSONProperty::as_list() const { std::vector copy; - + std::cout<<"TRYING AS LIST\n"; if (type == PropertyType::List) { for( auto & val : value_list){ copy.push_back(JSONProperty(val)); @@ -61,6 +61,7 @@ std::vector JSONProperty::as_list() const { return copy; } else if (type != PropertyType::Object) { + std::cout<<"as_list object\n"; copy.push_back(JSONProperty(*this)); return copy; } diff --git a/src/realizations/catchment/Bmi_Multi_Formulation.cpp b/src/realizations/catchment/Bmi_Multi_Formulation.cpp index 04f8889f33..abdec4b177 100644 --- a/src/realizations/catchment/Bmi_Multi_Formulation.cpp +++ b/src/realizations/catchment/Bmi_Multi_Formulation.cpp @@ -8,6 +8,7 @@ using namespace realization; void Bmi_Multi_Formulation::create_multi_formulation(geojson::PropertyMap properties, bool needs_param_validation) { + std::cout<<"CREATING MULTI\n"; if (needs_param_validation) { validate_parameters(properties); } @@ -33,7 +34,9 @@ void Bmi_Multi_Formulation::create_multi_formulation(geojson::PropertyMap proper // TODO: go back and set this up properly in required params collection auto nested_module_configs_it = properties.find(BMI_REALIZATION_CFG_PARAM_REQ__MODULES); + std::cout<<"GOT SMOE VARS\n"; std::vector nested_module_configs = nested_module_configs_it->second.as_list(); + std::cout<<"FS LIST\n"; // By default, have the "primary" module be the last primary_module_index = nested_module_configs.size() - 1; modules = std::vector(nested_module_configs.size()); @@ -43,6 +46,7 @@ void Bmi_Multi_Formulation::create_multi_formulation(geojson::PropertyMap proper /* ************************ Begin outer loop: "for sub_formulations_list" ************************ */ for (size_t i = 0; i < nested_module_configs.size(); ++i) { geojson::JSONProperty formulation_config = nested_module_configs[i]; + std::cout<<"I DONT KNOW ANYMORE\n"; std::string type_name = formulation_config.at("name").as_string(); std::string identifier = get_catchment_id() + "." + std::to_string(i); nested_module_ptr module = nullptr; diff --git a/test/realizations/Formulation_Manager_Test.cpp b/test/realizations/Formulation_Manager_Test.cpp index c062c4717c..38d3dbc039 100644 --- a/test/realizations/Formulation_Manager_Test.cpp +++ b/test/realizations/Formulation_Manager_Test.cpp @@ -602,6 +602,7 @@ const std::string EXAMPLE_5_a = " \"init_config\": \"\"," " \"allow_exceed_end_time\": true," " \"main_output_variable\": \"OUTPUT_VAR_4\"," +" \"uses_forcing_file\": false," " \"modules\": [" " {" " \"name\": \"bmi_c++\"," @@ -920,6 +921,7 @@ TEST_F(Formulation_Manager_Test, read_external_attributes) { } for (auto& expect : expected) { + std::cout<<"CHECKING "<fabric, catchment_output); - + std::cout<<"HERERERERE\n"; ASSERT_EQ(manager.get_size(), 3); check_formulation_values(manager, "cat-67", { 1.70352, 10.0 }); check_formulation_values(manager, "cat-52", { 3.14159, 15.0 }); diff --git a/utilities/plotting/forcing_animation.py b/utilities/plotting/forcing_animation.py new file mode 100755 index 0000000000..4beafbb8af --- /dev/null +++ b/utilities/plotting/forcing_animation.py @@ -0,0 +1,81 @@ +#!/usr/bin/env python +import pandas as pd +import numpy as np +import geopandas as gpd +import matplotlib.pyplot as plt +from matplotlib.animation import FuncAnimation, PillowWriter, FFMpegWriter +from matplotlib.colors import ListedColormap +from itertools import cycle +import shapely as shp +import seaborn as sns +from pathlib import Path +import xarray as xr +from functools import partial + +import os +catchments = gpd.read_file("gauge_01073000.gpkg", layer="divides", include_fields=['divide_id']) +forcing = xr.open_dataset("data/Long_Range_Hydrofabric_Precip.nc", decode_cf=False) + +forecast="Short Range" #TODO add product (CFS, GFS ect...) +forecast="Long Range" #TODO add product (CFS, GFS ect...) + +forcing.time.attrs['units'] = forcing.time.attrs['reference_time'] +forcing = xr.decode_cf(forcing) +forcing = forcing.set_coords("catchment_ids") + +cmap = sns.color_palette('flare', as_cmap=True) + +e_color = 'none' +vmin = forcing['RAINRATE'].min().item() +vmax = forcing['RAINRATE'].max().item() + +def plot_frame(ax, i): + + t = forcing['RAINRATE'][i:i+1] + #print(t) + #each dataset is a slice of time + ts = t.time.dt.strftime("%Y-%m-%d %H:%M:%S")[0].item() + + df = t.to_dataframe().droplevel(0) + df['catchment_ids'] = df['catchment_ids'].apply(lambda x: "cat-"+str(x)) + #merge on catchment ids + rdf = catchments.merge(df, left_on="divide_id", right_on="catchment_ids") + #rdf.plot(ax=ax, column="RAINRATE", cmap=cmap, legend=False, categorical=False, edgecolor=e_color) + ax.set_title(forecast+"\n"+ts+"\nPrecipitation", fontsize=20) + #FIXME don't plot new, update old!!! + ax.collections[0].set_array(rdf['RAINRATE'].values) + #ax.redraw_in_frame() + print("Plotting time ", i) + return ax + +fig = plt.figure(figsize=(19.20, 10.80)) +ax = fig.gca() +#plt.suptitle("Short Range") +cax = fig.add_axes([0.9, 0.1, 0.03, 0.8]) + +sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=vmin, vmax=vmax)) +# fake up the array of the scalar mappable. Urgh... +sm._A = [] +cbar = fig.colorbar(sm, cax=cax) +cbar.set_label(forcing['RAINRATE'].attrs['units']) +#ax = fig.gca() +ax.set_axis_off() +ax.title.set_size(20) + +ax = catchments.plot(ax=ax, legend=False, edgecolor=e_color, facecolor='white') +ax.collections[0].set_cmap(cmap) +ax.collections[0].set_clim(vmin, vmax) + +print("NUM FRAMES") +print(forcing.coords['time'].size) +#plot_frame(ax, 224) +ani = FuncAnimation(fig, partial(plot_frame, ax), frames=forcing.coords['time'].size, interval=100) +#plt.show() +#ani.save("forcings_"+forecast.replace(" ","_")+".gif", dpi=1200, writer=PillowWriter(fps=20)) +plt.tight_layout() +plt.subplots_adjust(top=0.85) +ani.save("forcings_"+forecast.replace(" ","_")+".mp4", writer=FFMpegWriter(fps=60) ) +plt.close() +#plt.savefig("test.png", dpi=1200, bbox_inches='tight', +# transparent=True, +# pad_inches=0) \ No newline at end of file