diff --git a/wntr/epanet/io.py b/wntr/epanet/io.py index fea4e3651..8add03604 100644 --- a/wntr/epanet/io.py +++ b/wntr/epanet/io.py @@ -2806,11 +2806,26 @@ def read(self, filename, convergence_error=False, darcy_weisbach=False, convert= self.results.link["headloss"] = pd.DataFrame(data=headloss, columns=linknames, index=reporttimes) status = np.array(df['linkstatus']) + if self.convert_status: - status[status <= 2] = 0 - status[status == 3] = 1 - status[status >= 5] = 1 - status[status == 4] = 2 + """ + EPANET status codes + 0 = closed (max head exceeded) + 1 = temporarily closed + 2 = closed + 3 = open + 4 = active (partially open) + 5 = open (max flow exceeded) + 6 = open (flow setting not met) + 7 = open (press setting not met) + """ + # 0 = 0, treat closed (max head exceeded) pump as closed + # 1 = 1, treat temporarily closed pipe as open + status[status == 2] = 0 # 2 = 0, closed + status[status == 3] = 1 # 3 = 1, open + status[status == 4] = 2 # 4 = 2, active + status[status >= 5] = 1 # 5,6,7 = 1, treat valve open under different conditions as open + self.results.link['status'] = pd.DataFrame(data=status, columns=linknames, index=reporttimes) setting = np.array(df['linksetting']) diff --git a/wntr/network/controls.py b/wntr/network/controls.py index 7838acd08..bd2783652 100644 --- a/wntr/network/controls.py +++ b/wntr/network/controls.py @@ -253,6 +253,25 @@ def __init__(self): def _reset(self): pass + def _shift(self, value): + """ + Shift any SimTimeConditions within larger condition rules by value seconds (backward). + + I.e., if a control is scheduled at simulation time 7200 and you shift by 3600, the new + control threshold with be sim time 3600. + + Parameters + ---------- + value : float + seconds to subtract from threshold + + Returns + ------- + bool + is this still a valid control? + """ + return True + @abc.abstractmethod def requires(self): """ @@ -585,6 +604,13 @@ def _compare(self, other): return False return True + def _shift(self, value): + self._threshold -= value + if self._threshold >= 0: + return True + self._threshold = 0 + return False + @property def name(self): if not self._repeat: @@ -996,6 +1022,11 @@ def _reset(self): self._condition_1._reset() self._condition_2._reset() + def _shift(self, value): + success1 = self._condition_1._shift(value) + success2 = self._condition_2._shift(value) + return success1 or success2 + def _compare(self, other): """ Parameters @@ -1062,6 +1093,11 @@ def _reset(self): self._condition_1._reset() self._condition_2._reset() + def _shift(self, value): + success1 = self._condition_1._shift(value) + success2 = self._condition_2._shift(value) + return success1 or success2 + def _compare(self, other): """ Parameters @@ -2012,6 +2048,9 @@ def epanet_control_type(self): control_type: _ControlType """ return self._control_type + + def _shift(self, step): + return self._condition._shift(step) def requires(self): req = self._condition.requires() @@ -2177,6 +2216,9 @@ def __init__(self, condition, then_action: BaseControlAction, priority=ControlPr else: self._control_type = _ControlType.postsolve + def _shift(self, step): + return self._condition._shift(step) + @classmethod def _time_control(cls, wnm, run_at_time, time_flag, daily_flag, control_action, name=None): """ diff --git a/wntr/network/model.py b/wntr/network/model.py index 26aaeffc1..f8d0bad0b 100644 --- a/wntr/network/model.py +++ b/wntr/network/model.py @@ -1552,7 +1552,139 @@ def reset_initial_values(self): for name, control in self.controls(): control._reset() + def set_initial_conditions(self, results, ts=None, remove_controls=True, warn=False): + """ + Set the initial conditions of the network based on prior simulation results. + + Parameters + ---------- + results : SimulationResults + Results from a prior simulation + ts : int, optional + The time value (in seconds) from the results to use to select initial conditions, + by default None (which will use the final values) + remove_controls : bool, optional + If a rule or control has a SimTimeCondition that now occurs prior to simulation start, remove + the control, by default True. + warn : bool + Send a warning to the logger that the rule has been deleted, by default False. + When False, information is sent to the logger at the `info` level. + + Returns + ------- + list + Control names that have been, when `remove_controls is True`, + or need to be, when `remove_controls is False`, + removed from the water network model + + + Raises + ------ + NameError + If both `ts` and `idx` are passed in + IndexError + If `ts` is passed, but no such time exists in the results + ValueError + If the time selected is not a multiple of the pattern timestep + + + """ + if ts is None: + end_time = results.node['demand'].index[-1] + else: + ts = int(ts) + if ts in results.node['demand'].index: + end_time = ts + else: + raise IndexError('There is no time "{}" in the results'.format(ts)) + + # if end_time / self.options.time.pattern_timestep != end_time // self.options.time.pattern_timestep: + # raise ValueError('You must give a time step that is a multiple of the pattern_timestep ({})'.format(self.options.time.pattern_timestep)) + + self.sim_time = 0.0 + self._prev_sim_time = None + + for name, node in self.nodes(Junction): + node._head = None + node._demand = None + node._pressure = None + try: node.initial_quality = float(results.node['quality'].loc[end_time, name]) + except KeyError: pass + node._leak_demand = None + node._leak_status = False + node._is_isolated = False + + for name, node in self.nodes(Tank): + node._head = None + node._demand = None + node._pressure = None + node.init_level = float(results.node['head'].loc[end_time, name] - node.elevation) + try: node.initial_quality = float(results.node['quality'].loc[end_time, name]) + except KeyError: pass + node._prev_head = node.head + node._leak_demand = None + node._leak_status = False + node._is_isolated = False + + for name, node in self.nodes(Reservoir): + node._head = None + node._demand = None + node._pressure = None + try: node.initial_quality = float(results.node['quality'].loc[end_time, name]) + except KeyError: pass + node._leak_demand = None + node._is_isolated = False + + for name, link in self.links(Pipe): + link.initial_status = results.link['status'].loc[end_time, name] + try: link.initial_setting = results.link['setting'].loc[end_time, name] + except KeyError: link.initial_setting = link.setting + link._user_status = link.initial_status + link._internal_status = LinkStatus.Active + link._is_isolated = False + link._flow = None + link._prev_setting = None + + for name, link in self.links(Pump): + link.initial_status = results.link['status'].loc[end_time, name] + try: link.initial_setting = results.link['setting'].loc[end_time, name] + except KeyError: link.initial_setting = link.setting + link._user_status = link.initial_status + link._setting = link.initial_setting + link._internal_status = LinkStatus.Active + link._is_isolated = False + link._flow = None + if isinstance(link, PowerPump): + link.power = link._base_power + link._prev_setting = None + + for name, link in self.links(Valve): + link.initial_status = results.link['status'].loc[end_time, name] + try: link.initial_setting = results.link['setting'].loc[end_time, name] + except KeyError: link.initial_setting = link.setting + # print(name, link.initial_status, link.initial_setting) + link._user_status = link.initial_status + link._setting = link.initial_setting + link._internal_status = LinkStatus.Active + link._is_isolated = False + link._flow = None + link._prev_setting = None + + to_delete = [] + for name, control in self.controls(): + control._reset() + still_good = control._shift(end_time) + if not still_good: + to_delete.append(name) + + for name in to_delete: + msg = 'Rule {} {} removed from the network'.format(name, 'has been' if remove_controls else 'needs to be') + if warn: logger.warning(msg) + else: logger.info(msg) + if remove_controls: + self.remove_control(name) + return to_delete class PatternRegistry(Registry): """A registry for patterns.""" diff --git a/wntr/sim/results.py b/wntr/sim/results.py index 9c9127654..d39289d5e 100644 --- a/wntr/sim/results.py +++ b/wntr/sim/results.py @@ -2,22 +2,277 @@ import datetime import enum +import numpy as np +import pandas as pd class ResultsStatus(enum.IntEnum): converged = 1 error = 0 - -class SimulationResults(object): +class SimulationResults: """ Water network simulation results class. - """ - def __init__(self): + A small number of mathematical and statistical functions are also provided. + These functions are applied to all dataframes within the results object (or + between two results objects) by name, elementwise. + + Assuming ``A`` and ``B`` are both results objects that have the same time + indices for the results and which describe the same water network physical + model (i.e., have the same nodes and links), then the following functions + are defined: + ================== ============================================================ + Example function Description + ------------------ ------------------------------------------------------------ + ``C = A + B`` Add the values from A and B for each property + ``C = A - B`` Subtract the property values in B from A + ``C = A / B`` Divide the property values in A by the values in B + ``C = A / n`` Divide the property values in A by n [int]; + note that this only makes sense if calculating an average + ``C = pow(A, p)`` Raise the property values in A to the p-th power; + note the syntax ``C = pow(A, p, mod)`` can also be used + ``C = abs(A)`` Take the absolute value of the property values in A + ``C = -A`` Take the negative of all property values in A + ``C = +A`` Take the positive of all property values in A + ``C = A.max()`` Get the maximum value for each property for node/link + ``C = A.min()`` Get the minimum value for each property for node/link + ``C = A.sum()`` Take the sum of each property across time for each node/link + ================== ============================================================ + + As an example, to calculate the relative difference between the results of two + simulations, one could do: ``rel_dif = abs(A - B) / A`` (warning - this will operate + on link statuses as well, which may result in meaningless results for that + parameter). + + """ + _data_attributes = ["link", "node"] + + def __init__(self): # Simulation time series self.timestamp = str(datetime.datetime.now()) self.network_name = None - self.link = None - self.node = None + for attr in self._data_attributes: + setattr(self, attr, dict()) + + def __add__(self, other): + if not isinstance(other, SimulationResults): + raise ValueError(f"operating on a results object requires both be SimulationResults") + new = SimulationResults() + new.network_name = "{}[{}] + {}[{}]".format( + self.network_name, self.timestamp, other.network_name, other.timestamp + ) + + for attr in new._data_attributes: + for key in getattr(self, attr).keys(): + if key in getattr(other, attr).keys(): + self_dict = getattr(self, attr) + other_dict = getattr(other, attr) + getattr(new, attr)[key] = self_dict[key] + other_dict[key] + return new + + + def __sub__(self, other): + if not isinstance(other, SimulationResults): + raise ValueError(f"operating on a results object requires both be SimulationResults") + new = SimulationResults() + new.network_name = "{}[{}] - {}[{}]".format( + self.network_name, self.timestamp, other.network_name, other.timestamp + ) + for attr in new._data_attributes: + for key in getattr(self, attr).keys(): + if key in getattr(other, attr).keys(): + self_dict = getattr(self, attr) + other_dict = getattr(other, attr) + getattr(new, attr)[key] = self_dict[key] - other_dict[key] + return new + + def __abs__(self): + new = SimulationResults() + new.network_name = "|{}[{}]|".format(self.network_name, self.timestamp) + + for attr in new._data_attributes: + for key in getattr(self, attr).keys(): + self_dict = getattr(self, attr) + new_dict = getattr(new, attr) + new_dict[key] = abs(self_dict[key]) + return new + + def __neg__(self): + new = SimulationResults() + new.network_name = "-{}[{}]".format(self.network_name, self.timestamp) + + for attr in new._data_attributes: + for key in getattr(self, attr).keys(): + self_dict = getattr(self, attr) + new_dict = getattr(new, attr) + new_dict[key] = -self_dict[key] + return new + + def min(self): + """Min operates on each axis separately, therefore it needs to be a function. + The built-in ``min`` function will not work.""" + new = SimulationResults() + new.network_name = "min({}[{}])".format(self.network_name, self.timestamp) + for attr in new._data_attributes: + for key in getattr(self, attr).keys(): + self_dict = getattr(self, attr) + new_dict = getattr(new, attr) + new_dict[key] = self_dict[key].min(axis=0) + return new + + def max(self): + """Max operates on each axis separately, therefore it needs to be a function. + The built-in ``max`` function will not work.""" + new = SimulationResults() + new.network_name = "max({}[{}])".format(self.network_name, self.timestamp) + + for attr in new._data_attributes: + for key in getattr(self, attr).keys(): + self_dict = getattr(self, attr) + new_dict = getattr(new, attr) + new_dict[key] = self_dict[key].max(axis=0) + return new + + def sum(self): + """Sum across time for each node/link for each property.""" + new = SimulationResults() + new.network_name = "sum({}[{}])".format(self.network_name, self.timestamp) + + for attr in new._data_attributes: + for key in getattr(self, attr).keys(): + self_dict = getattr(self, attr) + new_dict = getattr(new, attr) + new_dict[key] = self_dict[key].sum(axis=0) + return new + + def __pos__(self): + new = SimulationResults() + new.network_name = "+{}[{}]".format(self.network_name, self.timestamp) + for attr in new._data_attributes: + for key in getattr(self, attr).keys(): + self_dict = getattr(self, attr) + new_dict = getattr(new, attr) + new_dict[key] = +self_dict[key] + return new + + def __truediv__(self, other): + new = SimulationResults() + + if isinstance(other, SimulationResults): + new.network_name = "{}[{}] / {}[{}]".format( + self.network_name, self.timestamp, other.network_name, other.timestamp + ) + for attr in new._data_attributes: + for key in getattr(self, attr).keys(): + if key in getattr(other, attr).keys(): + self_dict = getattr(self, attr) + other_dict = getattr(other, attr) + getattr(new, attr)[key] = self_dict[key] / other_dict[key] + return new + + + elif isinstance(other, int): + new.network_name = "{}[{}] / {}".format(self.network_name, self.timestamp, other) + for key in self.link.keys(): + new.link[key] = self.link[key] / other + for key in self.node.keys(): + new.node[key] = self.node[key] / other + + for attr in new._data_attributes: + for key in getattr(self, attr).keys(): + self_dict = getattr(self, attr) + getattr(new, attr)[key] = self_dict[key] / other + return new + + else: + raise ValueError(f"using / on a results object requires the divisor to be SimulationResults or int") + + def __pow__(self, exp, mod=None): + new = SimulationResults() + new.network_name = "{}[{}] ** {}".format(self.network_name, self.timestamp, exp) + for attr in new._data_attributes: + for key in getattr(self, attr).keys(): + self_dict = getattr(self, attr) + new_dict = getattr(new, attr) + new_dict[key] = pow(self_dict[key], exp, mod) + return new + + def _adjust_time(self, ts: int): + """ + Adjust the time index for the results object by `ts`. + + Parameters + ---------- + ts : int + The number of seconds by which to adjust the result dataframe index + + Returns + ------- + self: SimulationResults + """ + ts = int(ts) + + for attr in self._data_attributes: + for key in getattr(self, attr).keys(): + self_dict = getattr(self, attr) + self_dict[key].index += ts + return self + + def append(self, other): + """ + Combine two results objects into a single, new result object. + If the times overlap, then the results from the `other` object will take precedence + over the values in the calling object. I.e., given ``A.append(B)``, + where ``A`` and ``B`` + are both `SimluationResults`, any results from ``A`` that relate to times equal to or + greater than the starting time of results in ``B`` will be dropped. + + .. warning:: + + This operations will be performed "in-place" and will change ``A`` + + + Parameters + ---------- + other : SimulationResults + Results objects from a different, and subsequent, simulation. + + Returns + ------- + self : SimulationResults + + Raises + ------ + ValueError + if `other` is the wrong type + + """ + if not isinstance(other, SimulationResults): + raise ValueError(f"operating on a results object requires both be SimulationResults") + #NOTE: Below two lines assume that results object has the node attribute and "head" key + + for attr in self._data_attributes: + self_dict = getattr(self, attr) + other_dict = getattr(other, attr) + + # query for any dataframe to get shape of df in each results object + self_attr_dataframe = list(self_dict.values())[0] + other_attr_dataframe = list(other_dict.values())[0] + + + all_keys = list(self_dict.keys()) + list(other_dict.keys()) + for key in all_keys: + self_df = self_dict[key] + other_df = other_dict[key] + overlap = self_df.index.intersection(other_df.index) + if key in self_dict.keys() and key in other_dict.keys(): + self_dict[key] = pd.concat([self_df[~self_df.index.isin(overlap)], other_df]) + elif key not in other_dict.keys(): + temp = other_attr_dataframe * np.nan + self_dict[key] = pd.concat([self_df[~self_df.index.isin(overlap)], temp]) + elif key not in self_dict.keys(): + temp = self_attr_dataframe * np.nan + self_dict[key] = pd.concat([temp[~temp.index.isin(overlap)], other_df]) + return self diff --git a/wntr/tests/test_multiple_simulations.py b/wntr/tests/test_multiple_simulations.py deleted file mode 100644 index 2bb344eaf..000000000 --- a/wntr/tests/test_multiple_simulations.py +++ /dev/null @@ -1,259 +0,0 @@ -# These tests run a demand-driven simulation with both WNTR and Epanet and compare the results for the example networks -import pickle -import unittest -from os.path import abspath, dirname, join - -import pandas as pd - -testdir = dirname(abspath(str(__file__))) -test_datadir = join(testdir, "networks_for_testing") -ex_datadir = join(testdir, "..", "..", "examples", "networks") - - -class TestResetInitialValues(unittest.TestCase): - @classmethod - def setUpClass(self): - import wntr - - self.wntr = wntr - - inp_file = join(ex_datadir, "Net3.inp") - self.wn = self.wntr.network.WaterNetworkModel(inp_file) - self.wn.options.time.hydraulic_timestep = 3600 - self.wn.options.time.duration = 24 * 3600 - - sim = self.wntr.sim.WNTRSimulator(self.wn) - self.res1 = sim.run_sim(solver_options={"TOL": 1e-8}) - - self.wn.reset_initial_values() - self.res2 = sim.run_sim(solver_options={"TOL": 1e-8}) - - @classmethod - def tearDownClass(self): - pass - - def test_link_flowrate(self): - for link_name, link in self.wn.links(): - for t in self.res1.time: - self.assertAlmostEqual( - self.res1.link["flowrate"].at[t, link_name], - self.res2.link["flowrate"].at[t, link_name], - 7, - ) - - def test_link_velocity(self): - for link_name, link in self.wn.links(): - for t in self.res1.link["velocity"].index: - self.assertAlmostEqual( - self.res1.link["velocity"].at[t, link_name], - self.res2.link["velocity"].at[t, link_name], - 7, - ) - - def test_node_demand(self): - for node_name, node in self.wn.nodes(): - for t in self.res1.node["demand"].index: - self.assertAlmostEqual( - self.res1.node["demand"].at[t, node_name], - self.res2.node["demand"].at[t, node_name], - 7, - ) - - def test_node_head(self): - for node_name, node in self.wn.nodes(): - for t in self.res1.node["head"].index: - self.assertAlmostEqual( - self.res1.node["head"].at[t, node_name], - self.res2.node["head"].at[t, node_name], - 7, - ) - - def test_node_pressure(self): - for node_name, node in self.wn.nodes(): - for t in self.res1.node["pressure"].index: - self.assertAlmostEqual( - self.res1.node["pressure"].at[t, node_name], - self.res2.node["pressure"].at[t, node_name], - 7, - ) - - -class TestStopStartSim(unittest.TestCase): - @classmethod - def setUpClass(self): - import wntr - - self.wntr = wntr - - inp_file = join(ex_datadir, "Net3.inp") - - parser = self.wntr.epanet.InpFile() - self.wn = parser.read(inp_file) - self.wn.options.time.hydraulic_timestep = 3600 - self.wn.options.time.duration = 24 * 3600 - sim = self.wntr.sim.WNTRSimulator(self.wn) - self.res1 = sim.run_sim(solver_options={"TOL": 1e-8}) - - parser = self.wntr.epanet.InpFile() - self.wn = parser.read(inp_file) - self.wn.options.time.hydraulic_timestep = 3600 - self.wn.options.time.duration = 10 * 3600 - sim = self.wntr.sim.WNTRSimulator(self.wn) - self.res2 = sim.run_sim(solver_options={'TOL':1e-8}) - self.wn.options.time.duration = 24*3600 - self.res3 = sim.run_sim(solver_options={'TOL':1e-8}) - - node_res = {} - link_res = {} - for key in self.res2.node.keys(): - node_res[key] = pd.concat([self.res2.node[key],self.res3.node[key]],axis=0) - for key in self.res2.link.keys(): - link_res[key] = pd.concat([self.res2.link[key],self.res3.link[key]],axis=0) - self.res2.node = node_res - self.res2.link = link_res - - @classmethod - def tearDownClass(self): - pass - - def test_link_flowrate(self): - for link_name, link in self.wn.links(): - for t in self.res1.link["flowrate"].index: - self.assertAlmostEqual( - self.res1.link["flowrate"].at[t, link_name], - self.res2.link["flowrate"].at[t, link_name], - 7, - ) - - def test_link_velocity(self): - for link_name, link in self.wn.links(): - for t in self.res1.link["velocity"].index: - self.assertAlmostEqual( - self.res1.link["velocity"].at[t, link_name], - self.res2.link["velocity"].at[t, link_name], - 7, - ) - - def test_node_demand(self): - for node_name, node in self.wn.nodes(): - for t in self.res1.node["demand"].index: - self.assertAlmostEqual( - self.res1.node["demand"].at[t, node_name], - self.res2.node["demand"].at[t, node_name], - 7, - ) - - def test_node_head(self): - for node_name, node in self.wn.nodes(): - for t in self.res1.node["head"].index: - self.assertAlmostEqual( - self.res1.node["head"].at[t, node_name], - self.res2.node["head"].at[t, node_name], - 7, - ) - - def test_node_pressure(self): - for node_name, node in self.wn.nodes(): - for t in self.res1.node["pressure"].index: - self.assertAlmostEqual( - self.res1.node["pressure"].at[t, node_name], - self.res2.node["pressure"].at[t, node_name], - 7, - ) - - -class TestPickle(unittest.TestCase): - @classmethod - def setUpClass(self): - import wntr - - self.wntr = wntr - - inp_file = join(ex_datadir, "Net3.inp") - - parser = self.wntr.epanet.InpFile() - self.wn = parser.read(inp_file) - self.wn.options.time.hydraulic_timestep = 3600 - self.wn.options.time.duration = 24 * 3600 - sim = self.wntr.sim.WNTRSimulator(self.wn) - self.res1 = sim.run_sim(solver_options={"TOL": 1e-8}) - - parser = self.wntr.epanet.InpFile() - self.wn = parser.read(inp_file) - self.wn.options.time.hydraulic_timestep = 3600 - self.wn.options.time.duration = 10 * 3600 - sim = self.wntr.sim.WNTRSimulator(self.wn) - self.res2 = sim.run_sim(solver_options={"TOL": 1e-8}) - f = open("temp.pickle", "wb") - pickle.dump(self.wn, f) - f.close() - f = open("temp.pickle", "rb") - wn2 = pickle.load(f) - f.close() - #wn2.set_initial_conditions(self.res2) - wn2.options.time.duration = 24*3600 - sim = self.wntr.sim.WNTRSimulator(wn2) - self.res3 = sim.run_sim(solver_options={'TOL':1e-8}) - - node_res = {} - link_res = {} - for key in self.res2.node.keys(): - node_res[key] = pd.concat([self.res2.node[key],self.res3.node[key]],axis=0) - for key in self.res2.link.keys(): - link_res[key] = pd.concat([self.res2.link[key],self.res3.link[key]],axis=0) - - self.res2.node = node_res - self.res2.link = link_res - - @classmethod - def tearDownClass(self): - pass - - def test_link_flowrate(self): - for link_name, link in self.wn.links(): - for t in self.res1.link["flowrate"].index: - self.assertAlmostEqual( - self.res1.link["flowrate"].at[t, link_name], - self.res2.link["flowrate"].at[t, link_name], - 7, - ) - - def test_link_velocity(self): - for link_name, link in self.wn.links(): - for t in self.res1.link["velocity"].index: - self.assertAlmostEqual( - self.res1.link["velocity"].at[t, link_name], - self.res2.link["velocity"].at[t, link_name], - 7, - ) - - def test_node_demand(self): - for node_name, node in self.wn.nodes(): - for t in self.res1.node["demand"].index: - self.assertAlmostEqual( - self.res1.node["demand"].at[t, node_name], - self.res2.node["demand"].at[t, node_name], - 7, - ) - - def test_node_head(self): - for node_name, node in self.wn.nodes(): - for t in self.res1.node["head"].index: - self.assertAlmostEqual( - self.res1.node["head"].at[t, node_name], - self.res2.node["head"].at[t, node_name], - 7, - ) - - def test_node_pressure(self): - for node_name, node in self.wn.nodes(): - for t in self.res1.node["pressure"].index: - self.assertAlmostEqual( - self.res1.node["pressure"].at[t, node_name], - self.res2.node["pressure"].at[t, node_name], - 7, - ) - - -if __name__ == "__main__": - unittest.main() diff --git a/wntr/tests/test_pickle_reset.py b/wntr/tests/test_pickle_reset.py new file mode 100644 index 000000000..5d04b5975 --- /dev/null +++ b/wntr/tests/test_pickle_reset.py @@ -0,0 +1,90 @@ +# Test the use of pickling and then unpickling a water network object for resetting values +# Also tests pickling and unpickling in the middle of a simulation +import pickle +import unittest +from os.path import abspath, dirname, join + +import pandas as pd +import wntr + +testdir = dirname(abspath(str(__file__))) +test_datadir = join(testdir, "networks_for_testing") +ex_datadir = join(testdir, "..", "..", "examples", "networks") + +tolerances = { + "flowrate": 1.0e-5, # 10 mL/s + "velocity": 1.0e-2, # 0.01 m/s + "demand": 1.0e-5, # 10 mL/s + "head": 1.0e-2, # 0.01 m + "pressure": 1.0e-2, # 0.01 m H2O + "headloss": 1.0e-2, # 0.01 + "status": 0.5, # i.e., 0 since this is integer + "setting": 1.0e-2, # 0.01 +} + +class TestPickle(unittest.TestCase): + @classmethod + def setUpClass(self): + + inp_file = join(ex_datadir, "Net3.inp") + + self.wn = wntr.network.WaterNetworkModel(inp_file) + self.wn.options.time.hydraulic_timestep = 3600 + self.wn.options.time.duration = 24 * 3600 + + # pickle WN before running sims + f = open("temp.pickle", "wb") + pickle.dump(self.wn, f) + f.close() + + + # run sim with initial wn + sim = wntr.sim.WNTRSimulator(self.wn) + self.res1 = sim.run_sim(solver_options={"TOL": 1e-8}) + + # load pickled inital wn + f = open("temp.pickle", "rb") + self.wn2 = pickle.load(f) + f.close() + + self.wn2.options.time.hydraulic_timestep = 3600 + self.wn2.options.time.duration = 10 * 3600 + sim = wntr.sim.WNTRSimulator(self.wn2) + self.res2 = sim.run_sim(solver_options={"TOL": 1e-8}) + self.wn2.set_initial_conditions(self.res2, 36000) + self.wn2.options.time.pattern_start = self.wn2.options.time.pattern_start + 10 * 3600 + + # pickle wn and reload and set conditions from previous sim + f = open("temp.pickle", "wb") + pickle.dump(self.wn2, f) + f.close() + f = open("temp.pickle", "rb") + self.wn2 = pickle.load(f) + f.close() + + self.wn2.options.time.duration = 14*3600 + sim = wntr.sim.WNTRSimulator(self.wn2) + self.res3 = sim.run_sim(solver_options={'TOL':1e-8}) + self.res3._adjust_time(10*3600) + + self.res2.append(self.res3) + self.res4 = abs(self.res1 - self.res2).max() + + @classmethod + def tearDownClass(self): + pass + + def test_nodes(self): + node_keys = ["demand", "head", "pressure"] + for key in node_keys: + max_res_diff = self.res4.node[key].max() + self.assertLess(max_res_diff, tolerances[key]) + + def test_links(self): + link_keys = ["flowrate", "velocity", "status", "setting"] + for key in link_keys: + max_res_diff = self.res4.link[key].max() + self.assertLess(max_res_diff, tolerances[key]) + +if __name__ == "__main__": + unittest.main() diff --git a/wntr/tests/test_reset_initial_values.py b/wntr/tests/test_reset_initial_values.py new file mode 100644 index 000000000..d142201a7 --- /dev/null +++ b/wntr/tests/test_reset_initial_values.py @@ -0,0 +1,61 @@ +# These tests run a demand-driven simulation with both WNTR and Epanet and compare the results for the example networks +import pickle +import unittest +from os.path import abspath, dirname, join + +import pandas as pd + +testdir = dirname(abspath(str(__file__))) +test_datadir = join(testdir, "networks_for_testing") +ex_datadir = join(testdir, "..", "..", "examples", "networks") + +tolerances = { + "flowrate": 1.0e-5, # 10 mL/s + "velocity": 1.0e-2, # 0.01 m/s + "demand": 1.0e-5, # 10 mL/s + "head": 1.0e-2, # 0.01 m + "pressure": 1.0e-2, # 0.01 m H2O + "headloss": 1.0e-2, # 0.01 + "status": 0.5, # i.e., 0 since this is integer + "setting": 1.0e-2, # 0.01 +} + +class TestResetInitialValues(unittest.TestCase): + @classmethod + def setUpClass(self): + import wntr + + self.wntr = wntr + self.tol = 1e-8 + + inp_file = join(ex_datadir, "Net3.inp") + self.wn = self.wntr.network.WaterNetworkModel(inp_file) + self.wn.options.time.hydraulic_timestep = 3600 + self.wn.options.time.duration = 24 * 3600 + + sim = self.wntr.sim.WNTRSimulator(self.wn) + self.res1 = sim.run_sim(solver_options={"TOL": 1e-8}) + + self.wn.reset_initial_values() + self.res2 = sim.run_sim(solver_options={"TOL": 1e-8}) + + self.res4 = abs(self.res1 - self.res2).max() + + @classmethod + def tearDownClass(self): + pass + + def test_nodes(self): + node_keys = ["demand", "head", "pressure"] + for key in node_keys: + max_res_diff = self.res4.node[key].max() + self.assertLess(max_res_diff, tolerances[key]) + + def test_links(self): + link_keys = ["flowrate", "velocity", "status", "setting"] + for key in link_keys: + max_res_diff = self.res4.link[key].max() + self.assertLess(max_res_diff, tolerances[key]) + +if __name__ == "__main__": + unittest.main() diff --git a/wntr/tests/test_results.py b/wntr/tests/test_results.py new file mode 100644 index 000000000..5ea5b9f73 --- /dev/null +++ b/wntr/tests/test_results.py @@ -0,0 +1,133 @@ +# These tests run a demand-driven simulation with both WNTR and Epanet and compare the results for the example networks +import copy +import unittest +from os.path import abspath, dirname, join +from pandas.testing import assert_frame_equal, assert_series_equal, assert_index_equal + +import pandas as pd + +testdir = dirname(abspath(str(__file__))) +test_datadir = join(testdir, "networks_for_testing") +ex_datadir = join(testdir, "..", "..", "examples", "networks") + +tolerances = { + "flowrate": 1.0e-5, # 10 mL/s + "velocity": 1.0e-2, # 0.01 m/s + "demand": 1.0e-5, # 10 mL/s + "head": 1.0e-2, # 0.01 m + "pressure": 1.0e-2, # 0.01 m H2O + "headloss": 1.0e-2, # 0.01 + "status": 0.5, # i.e., 0 since this is integer + "setting": 1.0e-2, # 0.01 +} + +class TestResultsOperations(unittest.TestCase): + @classmethod + def setUpClass(self): + import wntr + + self.wntr = wntr + self.tol = 1e-8 + self.ts = 3600 + + inp_file = join(ex_datadir, "Net3.inp") + self.wn = self.wntr.network.WaterNetworkModel(inp_file) + self.wn.options.time.hydraulic_timestep = 3600 + self.wn.options.time.duration = 4 * 3600 + + sim = self.wntr.sim.WNTRSimulator(self.wn) + self.res = sim.run_sim(solver_options={"TOL": self.tol}) + + @classmethod + def tearDownClass(self): + pass + + def test_adjust_time(self): + res_copy = copy.deepcopy(self.res) + res_copy._adjust_time(self.ts) + + for attr in self.res._data_attributes: + for key in getattr(self.res, attr).keys(): + assert_index_equal(getattr(self.res, attr)[key].index + self.ts, getattr(res_copy, attr)[key].index) + + def test_append(self): + res_copy1 = copy.deepcopy(self.res) + res_copy2 = copy.deepcopy(self.res) + res_copy2._adjust_time(self.wn.options.time.duration) + res_copy1.append(res_copy2) + assert res_copy1.node["head"].shape == (9,97) + assert (res_copy1.node["head"].loc[0] == res_copy1.node["head"].loc[14400]).all() + + def test_arithmetic(self): + df = self.res.node["head"] + + # addition + temp_res = self.res + self.res + temp_df = temp_res.node["head"] + test_df = df + df + assert_frame_equal(temp_df, test_df) + + # subtraction + temp_res = self.res - self.res + temp_df = temp_res.node["head"] + test_df = df - df + assert_frame_equal(temp_df, test_df) + + # division + temp_res = self.res / self.res + temp_df = temp_res.node["head"] + test_df = df / df + assert_frame_equal(temp_df, test_df) + + # int division + temp_res = self.res / 2 + temp_df = temp_res.node["head"] + test_df = df / 2 + assert_frame_equal(temp_df, test_df) + + # power + temp_res = pow(self.res, 1/2) + temp_df = temp_res.node["head"] + test_df = pow(df, 1/2) + assert_frame_equal(temp_df, test_df) + + # abs + temp_res = abs(self.res) + temp_df = temp_res.node["head"] + test_df = abs(df) + assert_frame_equal(temp_df, test_df) + + # neg + temp_res = -self.res + temp_df = temp_res.node["head"] + test_df = -df + assert_frame_equal(temp_df, test_df) + + # pos + temp_res = +self.res + temp_df = temp_res.node["head"] + test_df = +df + assert_frame_equal(temp_df, test_df) + + # max + temp_res = self.res.max() + temp_df = temp_res.node["head"] + test_df = df.max(axis=0) + assert_series_equal(temp_df, test_df) + + # min + temp_res = self.res.min() + temp_df = temp_res.node["head"] + test_df = df.min(axis=0) + assert_series_equal(temp_df, test_df) + + # sum + temp_res = self.res.sum() + temp_df = temp_res.node["head"] + test_df = df.sum(axis=0) + assert_series_equal(temp_df, test_df) + + + +if __name__ == "__main__": + unittest.main() diff --git a/wntr/tests/test_stepwise_sim.py b/wntr/tests/test_stepwise_sim.py new file mode 100644 index 000000000..b1fa165dc --- /dev/null +++ b/wntr/tests/test_stepwise_sim.py @@ -0,0 +1,82 @@ +import unittest +from os.path import abspath, dirname, join +import wntr +import pandas as pd + +testdir = dirname(abspath(str(__file__))) +test_datadir = join(testdir, "networks_for_testing") +ex_datadir = join(testdir, "..", "..", "examples", "networks") + +tolerances = { + "flowrate": 1.0e-5, # 10 mL/s + "velocity": 1.0e-2, # 0.01 m/s + "demand": 1.0e-4, # 10 mL/s + "head": 1.0e-2, # 0.01 m + "pressure": 1.0e-2, # 0.01 m H2O + "headloss": 1.0e-2, # 0.01 + "status": 0.5, # i.e., 0 since this is integer + "setting": 1.0e-2, # 0.01 +} + + +class TestStepwiseSim(unittest.TestCase): + @classmethod + def setUpClass(self): + + inp_file = join(ex_datadir, "Net3.inp") + + # straight 24 hour simulation + wn = wntr.network.WaterNetworkModel(inp_file) + wn.options.time.hydraulic_timestep = 3600 + wn.options.time.duration = 24 * 3600 + sim = wntr.sim.EpanetSimulator(wn) + continuous_res = sim.run_sim() + + # 24 hour simulation done in 24 1-hour chunks + wn = wntr.network.WaterNetworkModel(inp_file) + wn.options.time.hydraulic_timestep = 3600 + wn.options.time.duration = 1 * 3600 + + for i in range(0,24): + # run simulation for ith step + sim = wntr.sim.EpanetSimulator(wn) + i_res = sim.run_sim() + + # update wn with ith results + wn.set_initial_conditions(i_res) + wn.options.time.pattern_start = wn.options.time.pattern_start + (1 * 3600) + + + # adjust time of ith results + start_time = i * 3600 + i_res._adjust_time(start_time) + + # concatenate step results + if i == 0: + stepwise_res = i_res + else: + stepwise_res.append(i_res, overwrite=True) + + self.diff_res = abs(continuous_res - stepwise_res).max() + + def test_nodes(self): + node_keys = ["demand", "head", "pressure"] + for key in node_keys: + max_res_diff = self.diff_res.node[key].max() + self.assertLess(max_res_diff, tolerances[key]) + + def test_links(self): + link_keys = ["flowrate", "velocity", "status", "setting", "headloss"] + for key in link_keys: + max_res_diff = self.diff_res.link[key].max() + self.assertLess(max_res_diff, tolerances[key]) + + + +# self.diff_res.node["pressure"].plot() + +# self.diff_res.link["status"].plot() + + +if __name__ == "__main__": + unittest.main() \ No newline at end of file diff --git a/wntr/tests/test_stop_start_sim.py b/wntr/tests/test_stop_start_sim.py new file mode 100644 index 000000000..d303dcf66 --- /dev/null +++ b/wntr/tests/test_stop_start_sim.py @@ -0,0 +1,113 @@ +import unittest +from os.path import abspath, dirname, join +import wntr + + +testdir = dirname(abspath(str(__file__))) +test_datadir = join(testdir, "networks_for_testing") +ex_datadir = join(testdir, "..", "..", "examples", "networks") + + +tolerances = { + "flowrate": 1.0e-5, # 10 mL/s + "velocity": 1.0e-2, # 0.01 m/s + "demand": 1.0e-4, # 10 mL/s + "head": 1.0e-2, # 0.01 m + "pressure": 1.0e-2, # 0.01 m H2O + "headloss": 1.0e-2, # 0.01 + "status": 0.5, # i.e., 0 since this is integer + "setting": 1.0e-2, # 0.01 +} + +class TestStopStartSim(unittest.TestCase): + @classmethod + def setUpClass(self): + inp_file = join(ex_datadir, "Net3.inp") + simulator = wntr.sim.WNTRSimulator + + # straight 24 hour simulation + self.wn = wntr.network.WaterNetworkModel(inp_file) + self.wn.options.time.hydraulic_timestep = 3600 + self.wn.options.time.duration = 24 * 3600 + sim = simulator(self.wn) + self.res1 = sim.run_sim(solver_options={"TOL": 1e-8}) + + # 24 hour simulation done in one 10 hour chunk and one 14 hour chunk + self.wn = wntr.network.WaterNetworkModel(inp_file) + self.wn.options.time.hydraulic_timestep = 3600 + self.wn.options.time.duration = 10 * 3600 + sim = simulator(self.wn) + self.res2 = sim.run_sim(solver_options={'TOL':1e-8}) + self.wn.set_initial_conditions(self.res2) + self.wn.options.time.pattern_start = self.wn.options.time.pattern_start + 10 * 3600 + self.wn.options.time.duration = 14 * 3600 + sim = simulator(self.wn) + self.res3 = sim.run_sim(solver_options={'TOL':1e-8}) + self.res3._adjust_time(10*3600) + self.res2.append(self.res3) + self.res4 = abs(self.res1 - self.res2).max() + + @classmethod + def tearDownClass(self): + pass + + def test_nodes(self): + node_keys = ["demand", "head", "pressure"] + for key in node_keys: + max_res_diff = self.res4.node[key].max() + self.assertLess(max_res_diff, tolerances[key]) + + def test_links(self): + link_keys = ["flowrate", "velocity", "status", "setting"] + for key in link_keys: + max_res_diff = self.res4.link[key].max() + self.assertLess(max_res_diff, tolerances[key]) + + +class TestStopStartEpanetSim(unittest.TestCase): + @classmethod + def setUpClass(self): + inp_file = join(ex_datadir, "Net3.inp") + + # straight 24 hour simulation + self.wn = wntr.network.WaterNetworkModel(inp_file) + self.wn.options.time.hydraulic_timestep = 3600 + self.wn.options.time.duration = 24 * 3600 + sim = wntr.sim.EpanetSimulator(self.wn) + self.res1 = sim.run_sim() + + # 24 hour simulation done in one 10 hour chunk and one 14 hour chunk + self.wn = wntr.network.WaterNetworkModel(inp_file) + self.wn.options.time.hydraulic_timestep = 3600 + self.wn.options.time.duration = 10 * 3600 + sim = wntr.sim.EpanetSimulator(self.wn) + self.res2 = sim.run_sim() + self.wn.set_initial_conditions(self.res2) + self.wn.options.time.pattern_start = self.wn.options.time.pattern_start + 10 * 3600 + self.wn.options.time.duration = 14 * 3600 + sim = wntr.sim.EpanetSimulator(self.wn) + self.res3 = sim.run_sim() + self.res3._adjust_time(10*3600) + + self.res2.append(self.res3) + self.res4 = abs(self.res1 - self.res2).max() + + @classmethod + def tearDownClass(self): + pass + + def test_nodes(self): + node_keys = ["demand", "head", "pressure"] + for key in node_keys: + max_res_diff = self.res4.node[key].max() + self.assertLess(max_res_diff, tolerances[key]) + + def test_links(self): + link_keys = ["flowrate", "velocity", "status", "setting", "headloss"] + for key in link_keys: + max_res_diff = self.res4.link[key].max() + self.assertLess(max_res_diff, tolerances[key]) + + +if __name__ == "__main__": + unittest.main() \ No newline at end of file