Repository: SDXorg/pysd

Issue:
Warning for user when non-constant functions are modified by 'params'
Say the user wants to set the initial conditions for the simulation, but they accidentally use the 'params' keyword. This will overwrite the stock calling functions with the values they meant to use as initial conditions, and they will get very unusual behavior.

We could check for this, and give a warning.


Discussion:
pep8speaks: Hello @enekomartinmartinez! Thanks for opening this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:







* In the file [`pysd/py_backend/model.py`](https://github.com/SDXorg/pysd/blob/30bda49028168d092bb7a14315d0bd9cf08cf7d8/pysd/py_backend/model.py):

> [Line 504:80](https://github.com/SDXorg/pysd/blob/30bda49028168d092bb7a14315d0bd9cf08cf7d8/pysd/py_backend/model.py#L504): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (131 > 79 characters)
> [Line 580:80](https://github.com/SDXorg/pysd/blob/30bda49028168d092bb7a14315d0bd9cf08cf7d8/pysd/py_backend/model.py#L580): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (131 > 79 characters)
> [Line 995:80](https://github.com/SDXorg/pysd/blob/30bda49028168d092bb7a14315d0bd9cf08cf7d8/pysd/py_backend/model.py#L995): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (118 > 79 characters)
> [Line 1004:80](https://github.com/SDXorg/pysd/blob/30bda49028168d092bb7a14315d0bd9cf08cf7d8/pysd/py_backend/model.py#L1004): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (144 > 79 characters)
> [Line 1433:80](https://github.com/SDXorg/pysd/blob/30bda49028168d092bb7a14315d0bd9cf08cf7d8/pysd/py_backend/model.py#L1433): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (107 > 79 characters)
> [Line 1552:80](https://github.com/SDXorg/pysd/blob/30bda49028168d092bb7a14315d0bd9cf08cf7d8/pysd/py_backend/model.py#L1552): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (131 > 79 characters)
> [Line 1638:80](https://github.com/SDXorg/pysd/blob/30bda49028168d092bb7a14315d0bd9cf08cf7d8/pysd/py_backend/model.py#L1638): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (130 > 79 characters)
> [Line 1795:80](https://github.com/SDXorg/pysd/blob/30bda49028168d092bb7a14315d0bd9cf08cf7d8/pysd/py_backend/model.py#L1795): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (120 > 79 characters)

























Patch:
diff --git a/docs/advanced_usage.rst b/docs/advanced_usage.rst
index 0c8f775..c583b19 100644
--- a/docs/advanced_usage.rst
+++ b/docs/advanced_usage.rst
@@ -145,6 +145,8 @@ A submodel of a translated model can be run as a standalone model. This can be d
 .. automethod:: pysd.py_backend.model.Model.select_submodel
    :noindex:
 
+.. note::
+   This method will mutate the original model. If you want to have a copy of the model with the selected variables/modules you can use :py:data:`inplace=False` argument.
 
 In order to preview the needed exogenous variables, the :py:meth:`.get_dependencies` method can be used:
 
diff --git a/docs/getting_started.rst b/docs/getting_started.rst
index b7ec317..b164612 100644
--- a/docs/getting_started.rst
+++ b/docs/getting_started.rst
@@ -272,3 +272,11 @@ We can easily access the current value of a model component using square bracket
 If you try to get the current values of a lookup variable, the previous method will fail, as lookup variables take arguments. However, it is possible to get the full series of a lookup or data object with :py:meth:`.get_series_data` method::
 
    >>> model.get_series_data('Growth lookup')
+
+
+Copying a model
+---------------
+Sometimes, you may want to run several versions of a model. For this purpose, copying an already-loaded model to make changes while keeping an untouched one is useful. The :py:meth:`.copy` method will help do that; it will load a new model from the translated file and apply to it the same changes that have been applied to the original model (modifying components, selecting submodels, etc.). You can also load a copy of the source model (without applying) any change setting the argument :py:data:`reload=True`.
+
+.. warning::
+   The copy function will load a new model from the file and apply the same changes to it. If any of these changes have replaced a variable with a function that references other variables in the model, the copy will not work properly since the function will still reference the variables in the original model, in which case the function should be redefined.
diff --git a/docs/python_api/model_class.rst b/docs/python_api/model_class.rst
index 55f918d..3da4238 100644
--- a/docs/python_api/model_class.rst
+++ b/docs/python_api/model_class.rst
@@ -11,3 +11,4 @@ Macro class
 -----------
 .. autoclass:: pysd.py_backend.model.Macro
     :members:
+    :exclude-members: export
diff --git a/docs/whats_new.rst b/docs/whats_new.rst
index f0506ce..4f52480 100644
--- a/docs/whats_new.rst
+++ b/docs/whats_new.rst
@@ -1,5 +1,41 @@
 What's New
 ==========
+v3.13.0 (2023/12/25)
+--------------------
+New Features
+~~~~~~~~~~~~
+- Include new method :py:meth:`pysd.py_backend.model.Model.copy` which allows copying a model (:issue:`131`). (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+- :py:meth:`pysd.py_backend.model.Model.select_submodel` now  takes an optional argument `inplace` when set to :py:data:`False` it will return a modified copy of the model instead of modifying the original model (:issue:`131`). (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+- :py:meth:`pysd.py_backend.model.Model.export` will now save also time component information if changed (e.g. final time, time step...). (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+
+Breaking changes
+~~~~~~~~~~~~~~~~
+
+Deprecations
+~~~~~~~~~~~~
+
+Bug fixes
+~~~~~~~~~
+- Set the pointer of :py:class:`pysd.py_backend.statefuls.DelayFixed` to 0 during initialization (:issue:`427`). (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+- :py:meth:`pysd.py_backend.model.Model.export` now works with Macros. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+
+Documentation
+~~~~~~~~~~~~~
+- Improve documentation of methods in :py:class:`pysd.py_backend.model.Model` and :py:class:`pysd.py_backend.model.Macro` includying cross-references and rewrite the one from :py:meth:`pysd.py_backend.model.Macro.set_components`. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+- Include documentation about the new method :py:meth:`pysd.py_backend.model.Model.copy` and update documentation from :py:meth:`pysd.py_backend.model.Model.select_submodel`. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+
+Performance
+~~~~~~~~~~~
+- Improved performace of :py:class:`pysd.py_backend.output.DataFrameHandler` by creating the dataframe at the end of the run (:issue:`374` and :issue:`330`). (`@easyas314159 <https://github.com/easyas314159>`_ and `@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+
+Internal Changes
+~~~~~~~~~~~~~~~~
+- Move old :py:meth:`pysd.py_backend.model.Macro.set_components` to :py:meth:`pysd.py_backend.model.Macro._set_components`, and create new method with the same name without the `new` argument. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+- Move old :py:meth:`pysd.py_backend.model.Macro.set_stateful` to :py:meth:`pysd.py_backend.model.Macro._set_stateful`. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+- Make integration tests filter only specific warnings. (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+- Include warnings in :py:meth:`pysd.py_backend.model.Macro.set_components` when changing the behaviour of the component (:issue:`58`). (`@enekomartinmartinez <https://github.com/enekomartinmartinez>`_)
+
+
 v3.12.0 (2023/10/02)
 --------------------
 New Features
diff --git a/pysd/_version.py b/pysd/_version.py
index d1a7f1e..62ee17b 100644
--- a/pysd/_version.py
+++ b/pysd/_version.py
@@ -1,1 +1,1 @@
-__version__ = "3.12.0"
+__version__ = "3.13.0"
diff --git a/pysd/builders/python/python_model_builder.py b/pysd/builders/python/python_model_builder.py
index ce43283..019a69c 100644
--- a/pysd/builders/python/python_model_builder.py
+++ b/pysd/builders/python/python_model_builder.py
@@ -663,7 +663,7 @@ class ElementBuilder:
         if ", " in self.type:
             warn(
                 f"Variable '{self.name}' is defined with different types:"
-                f" '{self.type}'. This may cause bugs when trying to"
+                f" '{self.type}'. This may cause bugs when trying to "
                 "change its value or applying other methods from the "
                 "pysd.py_backend.model.Model class. Running the model "
                 "without modifying this variable should not cause any "
@@ -672,7 +672,7 @@ class ElementBuilder:
         elif ", " in self.subtype:
             warn(
                 f"Variable '{self.name}' is defined with different subtypes:"
-                f" '{self.subtype}'. This may cause bugs when trying to"
+                f" '{self.subtype}'. This may cause bugs when trying to "
                 "change its value or applying other methods from the "
                 "pysd.py_backend.model.Model class. Running the model "
                 "without modifying this variable should not cause any "
diff --git a/pysd/py_backend/components.py b/pysd/py_backend/components.py
index 140ac9b..8eb79d9 100644
--- a/pysd/py_backend/components.py
+++ b/pysd/py_backend/components.py
@@ -7,6 +7,7 @@ import os
 import random
 import inspect
 import importlib.util
+from copy import deepcopy
 
 import numpy as np
 
@@ -129,13 +130,45 @@ class Time(object):
         self._time = None
         self.stage = None
         self.return_timestamps = None
+        self._next_return = None
+        self._control_vars_tracker = {}
 
     def __call__(self):
         return self._time
 
+    def export(self):
+        """Exports time values to a dictionary."""
+        return {
+            "control_vars": self._get_control_vars(),
+            "stage": self.stage,
+            "_time": self._time,
+            "return_timestamps": self.return_timestamps,
+            "_next_return": self._next_return
+        }
+
+    def _get_control_vars(self):
+        """
+        Make control vars changes exportable.
+        """
+        out = {}
+        for cvar, value in self._control_vars_tracker.items():
+            if callable(value):
+                out[cvar] = value()
+            else:
+                out[cvar] = value
+        return out
+
+    def _set_time(self, time_dict):
+        """Copy values from other Time object, used by Model.copy"""
+        self.set_control_vars(**time_dict['control_vars'])
+        for key, value in time_dict.items():
+            if key == 'control_vars':
+                continue
+            setattr(self, key, value)
+
     def set_control_vars(self, **kwargs):
         """
-        Set the control variables valies
+        Set the control variables values
 
         Parameters
         ----------
@@ -149,6 +182,20 @@ class Time(object):
             saveper: float, callable or None
                 Saveper.
 
+        """
+        # filter None values
+        kwargs = {
+                key: value for key, value in kwargs.items()
+                if value is not None
+        }
+        # track changes
+        self._control_vars_tracker.update(kwargs)
+        self._set_control_vars(**kwargs)
+
+    def _set_control_vars(self, **kwargs):
+        """
+        Set the control variables values. Private version to be used
+        to avoid tracking changes.
         """
         def _convert_value(value):
             # this function is necessary to avoid copying the pointer in the
@@ -159,8 +206,7 @@ class Time(object):
                 return lambda: value
 
         for key, value in kwargs.items():
-            if value is not None:
-                setattr(self, key, _convert_value(value))
+            setattr(self, key, _convert_value(value))
 
         if "initial_time" in kwargs:
             self._initial_time = self.initial_time()
@@ -184,16 +230,16 @@ class Time(object):
 
         if self.return_timestamps is not None:
             # this allows managing float precision error
-            if self.next_return is None:
+            if self._next_return is None:
                 return False
-            if np.isclose(self._time, self.next_return, prec):
+            if np.isclose(self._time, self._next_return, prec):
                 self._update_next_return()
                 return True
             else:
-                while self.next_return is not None\
-                      and self._time > self.next_return:
+                while self._next_return is not None\
+                      and self._time > self._next_return:
                     warn(
-                        f"The returning time stamp '{self.next_return}' "
+                        f"The returning time stamp '{self._next_return}' "
                         "seems to not be a multiple of the time step. "
                         "This value will not be saved in the output. "
                         "Please, modify the returning timestamps or the "
@@ -218,12 +264,12 @@ class Time(object):
            and len(return_timestamps) > 0:
             self.return_timestamps = list(return_timestamps)
             self.return_timestamps.sort(reverse=True)
-            self.next_return = self.return_timestamps.pop()
+            self._next_return = self.return_timestamps.pop()
         elif isinstance(return_timestamps, (float, int)):
-            self.next_return = return_timestamps
+            self._next_return = return_timestamps
             self.return_timestamps = []
         else:
-            self.next_return = None
+            self._next_return = None
             self.return_timestamps = None
 
     def update(self, value):
@@ -233,9 +279,9 @@ class Time(object):
     def _update_next_return(self):
         """ Update the next_return value """
         if self.return_timestamps:
-            self.next_return = self.return_timestamps.pop()
+            self._next_return = self.return_timestamps.pop()
         else:
-            self.next_return = None
+            self._next_return = None
 
     def reset(self):
         """ Reset time value to the initial """
diff --git a/pysd/py_backend/model.py b/pysd/py_backend/model.py
index 4f52b33..3f69a18 100644
--- a/pysd/py_backend/model.py
+++ b/pysd/py_backend/model.py
@@ -10,6 +10,7 @@ import warnings
 import inspect
 import pickle
 from pathlib import Path
+from copy import deepcopy
 from typing import Union
 
 import numpy as np
@@ -65,6 +66,10 @@ class Macro(DynamicStateful):
     py_name: str or None
         The name of the Macro object. Default is None.
 
+    See also
+    --------
+    :class:`pysd.py_backend.model.Model`
+
     """
     def __init__(self, py_model_file, params=None, return_func=None,
                  time=None, time_initialization=None, data_files=None,
@@ -81,6 +86,11 @@ class Macro(DynamicStateful):
         self.lookups_loaded = False
         # Functions with constant cache
         self._constant_funcs = set()
+        # Attributes that are set later
+        self.stateful_initial_dependencies = None
+        self.initialize_order = None
+        self.cache_type = None
+        self._components_setter_tracker = {}
         # Load model/macro from file and save in components
         self.components = Components(str(py_model_file), self.set_components)
 
@@ -111,7 +121,7 @@ class Macro(DynamicStateful):
             # add params to namespace
             self._namespace.update(self.components._components._params)
             # create new components with the params
-            self.set_components(params, new=True)
+            self._set_components(params, new=True)
             # update dependencies
             for param in params:
                 self._dependencies[
@@ -204,7 +214,7 @@ class Macro(DynamicStateful):
 
     def clean_caches(self):
         """
-        Clean the cahce of the object and the macros objects that it
+        Clean the cache of the object and the macros objects that it
         contains
         """
         self.cache.clean()
@@ -478,59 +488,6 @@ class Macro(DynamicStateful):
         [component.update(val) for component, val
          in zip(self._dynamicstateful_elements, new_value)]
 
-    def export(self, file_name):
-        """
-        Export stateful values to pickle file.
-
-        Parameters
-        ----------
-        file_name: str or pathlib.Path
-          Name of the file to export the values.
-
-        """
-        warnings.warn(
-            "\nCompatibility of exported states could be broken between"
-            " different versions of PySD or xarray, current versions:\n"
-            f"\tPySD {__version__}\n\txarray {xr.__version__}\n"
-        )
-        stateful_elements = {
-            name: element.export()
-            for name, element in self._stateful_elements.items()
-        }
-
-        with open(file_name, 'wb') as file:
-            pickle.dump(
-                (self.time(),
-                 stateful_elements,
-                 {'pysd': __version__, 'xarray': xr.__version__}
-                 ), file)
-
-    def import_pickle(self, file_name):
-        """
-        Import stateful values from pickle file.
-
-        Parameters
-        ----------
-        file_name: str or pathlib.Path
-          Name of the file to import the values from.
-
-        """
-        with open(file_name, 'rb') as file:
-            time, stateful_dict, metadata = pickle.load(file)
-
-        if __version__ != metadata['pysd']\
-           or xr.__version__ != metadata['xarray']:  # pragma: no cover
-            warnings.warn(
-                "\nCompatibility of exported states could be broken between"
-                " different versions of PySD or xarray. Current versions:\n"
-                f"\tPySD {__version__}\n\txarray {xr.__version__}\n"
-                "Loaded versions:\n"
-                f"\tPySD {metadata['pysd']}\n\txarray {metadata['xarray']}\n"
-                )
-
-        self.set_stateful(stateful_dict)
-        self.time.set_control_vars(initial_time=time)
-
     def initialize_external_data(self, externals=None):
 
         """
@@ -543,6 +500,9 @@ class Macro(DynamicStateful):
         To get the full performance gain of loading the externals from a netCDF
         file, the model should be loaded with initialize=False first.
 
+        Examples of usage are available at
+        `Advanced Usage <https://pysd.readthedocs.io/en/master/advanced_usage.html#initializing-external-data-from-netcdf-file>`__.
+
         Parameters
         ----------
         externals: str or pathlib.Path (optional)
@@ -552,6 +512,10 @@ class Macro(DynamicStateful):
         -------
         None
 
+        See also
+        --------
+        :func:`pysd.py_backend.model.Macro.serialize_externals`
+
         """
 
         if not externals:
@@ -612,6 +576,9 @@ class Macro(DynamicStateful):
         Names of variables should be those in the model (python safe,
         without the _ext_type_ string in front).
 
+        Examples of usage are available at
+        `Advanced Usage <https://pysd.readthedocs.io/en/master/advanced_usage.html#initializing-external-data-from-netcdf-file>`__.
+
         Parameters
         ----------
         export_path: str or pathlib.Path (optional)
@@ -636,6 +603,10 @@ class Macro(DynamicStateful):
         -------
         None
 
+        See also
+        --------
+        :func:`pysd.py_backend.model.Macro.initialize_external_data`
+
         """
         data = {}
         metadata = {}
@@ -726,7 +697,9 @@ class Macro(DynamicStateful):
         all ExtData objects to have many nans when stored in a xarray Dataset.
         It does the same for the "lookup_dim" of all ExtLookup objects.
 
-        NOTE: Though subscripts can be read from Excel, they are hardcoded
+        Note
+        ----
+        Though subscripts can be read from Excel, they are hardcoded
         during the model building process. Therefore they will not be
         serialized.
 
@@ -835,6 +808,10 @@ class Macro(DynamicStateful):
         >>> model.get_args('birth_rate')
         >>> model.get_args('Birth Rate')
 
+        See also
+        --------
+        :func:`pysd.py_backend.model.Macro.get_coords`
+
         """
         if isinstance(param, str):
             func_name = utils.get_key_and_value_by_insensitive_key_or_value(
@@ -875,6 +852,10 @@ class Macro(DynamicStateful):
         >>> model.get_coords('birth_rate')
         >>> model.get_coords('Birth Rate')
 
+        See also
+        --------
+        :func:`pysd.py_backend.model.Macro.get_args`
+
         """
         if isinstance(param, str):
             func_name = utils.get_key_and_value_by_insensitive_key_or_value(
@@ -926,6 +907,10 @@ class Macro(DynamicStateful):
         ----
         It will crash if the model component takes arguments.
 
+        See also
+        --------
+        :func:`pysd.py_backend.model.Macro.get_series_data`
+
         """
         func_name = utils.get_key_and_value_by_insensitive_key_or_value(
             param,
@@ -981,25 +966,78 @@ class Macro(DynamicStateful):
                 "Trying to get the values of a constant variable. "
                 "'model.get_series_data' only works lookups/data objects.\n\n")
 
-    def set_components(self, params, new=False):
-        """ Set the value of exogenous model elements.
-        Element values can be passed as keyword=value pairs in the
+    def set_components(self, params):
+        """
+        Set the value of exogenous model elements.
+        Element values should be passed with a dictionary in the
         function call. Values can be numeric type or pandas Series.
         Series will be interpolated by integrator.
 
+        Parameters
+        ----------
+        params: dict
+            Dictionary with the name of the elements to modify and the
+            value that they would take. If the passed value is a
+            :class:`float` or a :class:`xarray.DataArray` (must be
+            compatible with the dimensions of the variable). In this
+            case, the variable will have a constant value, returning
+            the past value and broadcasting to all dimensions, if
+            necessary. If a :class:`pandas.Series` is passed, the
+            variable will be of type data and will use the time of the
+            model to interpolate the result having as reference the
+            indexes of the series. In this case, the series 'data' can
+            also be a :class:`float` or a :class:`xarray.DataArray`,
+            as with constant values. In the case of the target using
+            a :class:`pysd.py_backend.lookups.Lookup` object, it will
+            modify the object values to use the original arguments when
+            being call. More detailed information and examples of usage
+            are available at
+            `Getting Started <https://pysd.readthedocs.io/en/master/getting_started.html#setting-parameter-values>`__.
+
+            To write more complex relationships, which may or may not
+            include other model variables, a callable, e.g. a function,
+            can be passed that takes the same arguments as the original
+            function and returns a :class:`float` or a
+            :class:`xarray.DataArray` with exactly the same dimensions
+            as the original function. More detailed information and
+            examples of usage are available at
+            `Advanced Usage <https://pysd.readthedocs.io/en/master/advanced_usage.html#replacing-model-components-with-more-complex-objects>`__.
+
+
+        Note
+        ----
+        This function is to modify the value or equations of the
+        variables, it won't work properly with Stateful objects, e.g.
+        Integ, DelayFixed... In order to modify them it will be
+        necessary to do it manually, if you have other inputs it is
+        recommended to modify these ones. To change their initial
+        value :func:`pysd.py_backend.model.Model.set_initial_condition`
+        method could be used.
+
+
         Examples
         --------
         >>> model.set_components({'birth_rate': 10})
         >>> model.set_components({'Birth Rate': 10})
 
-        >>> br = pandas.Series(index=range(30), values=np.sin(range(30))
+        >>> br = pandas.Series(index=range(30), data=np.sin(range(30))
         >>> model.set_components({'birth_rate': br})
 
+        See also
+        --------
+        :func:`pysd.py_backend.model.Model.set_initial_condition`
+        :func:`pysd.py_backend.model.Macro.get_coords`
+        :func:`pysd.py_backend.model.Macro.get_args`
+
+        """
+        self._components_setter_tracker.update(params)
+        self._set_components(params, new=False)
+
+    def _set_components(self, params, new):
+        """
+        Set the value of exogenous model elements, giving the option to
+        set new components (used in Macros).
         """
-        # TODO: allow the params argument to take a pandas dataframe, where
-        # column names are variable names. However some variables may be
-        # constant or have no values for some index. This should be processed.
-        # TODO: make this compatible with loading outputs from other files
 
         for key, value in params.items():
             func_name = utils.get_key_and_value_by_insensitive_key_or_value(
@@ -1037,6 +1075,10 @@ class Macro(DynamicStateful):
                     obj
                 ).set_values(value)
 
+                if not isinstance(value, pd.Series):
+                    warnings.warn(
+                        "Replacing interpolation data with constant values.")
+
                 # Update dependencies
                 if func_type == "Data":
                     if isinstance(value, pd.Series):
@@ -1048,11 +1090,29 @@ class Macro(DynamicStateful):
 
                 continue
 
+            if func_type == "Stateful":
+                warnings.warn(
+                    "Replacing the value of Stateful variable with "
+                    "an expression. To set initial conditions use "
+                    "`set_initial_condition` instead..."
+                )
+
             if isinstance(value, pd.Series):
+                if func_type == "Constant":
+                    warnings.warn(
+                        "Replacing a constant value with a "
+                        "time-dependent value. The value will be "
+                        "interpolated over time."
+                    )
                 new_function, deps = self._timeseries_component(
                     value, dims)
                 self._dependencies[func_name] = deps
             elif callable(value):
+                if func_type == "Constant":
+                    warnings.warn(
+                        "Replacing a constant value with a callable. "
+                        "The value may not be constant anymore."
+                    )
                 new_function = value
                 # Using step cache adding time as dependency
                 # TODO it would be better if we can parse the content
@@ -1060,6 +1120,8 @@ class Macro(DynamicStateful):
                 self._dependencies[func_name] = {"time": 1}
 
             else:
+                if func_type != "Constant":
+                    warnings.warn("Replacing a variable by a constant value.")
                 new_function = self._constant_component(value, dims)
                 self._dependencies[func_name] = {}
 
@@ -1102,9 +1164,9 @@ class Macro(DynamicStateful):
 
         else:
             # the interpolation will be time dependent
-            return lambda:\
-                np.interp(self.time(), series.index, series.values),\
-                {'time': 1}
+            return lambda: np.interp(
+                self.time(), series.index, series.values
+                ), {'time': 1}
 
     def _constant_component(self, value, dims):
         """ Internal function for creating a constant model element """
@@ -1115,21 +1177,26 @@ class Macro(DynamicStateful):
         else:
             return lambda: value
 
-    def set_initial_value(self, t, initial_value):
-        """ Set the system initial value.
+    def set_initial_value(self, time, initial_value):
+        """
+        Set the system initial value.
 
         Parameters
         ----------
-        t : numeric
-            The system time
+        time : float or int
+            The system intial time to be set.
 
         initial_value : dict
             A (possibly partial) dictionary of the system initial values.
             The keys to this dictionary may be either pysafe names or
-            original model file names
+            original model file names.
+
+        See also
+        --------
+        :func:`pysd.py_backend.model.Model.set_initial_condition`
 
         """
-        self.time.set_control_vars(initial_time=t)
+        self.time.set_control_vars(initial_time=time)
         stateful_name = "_NONE"
         modified_statefuls = set()
 
@@ -1209,7 +1276,14 @@ class Macro(DynamicStateful):
 
         return elements_to_initialize
 
-    def set_stateful(self, stateful_dict):
+    def export(self):
+        """Exports stateful values to a dictionary."""
+        return {
+            name: element.export()
+            for name, element in self._stateful_elements.items()
+        }
+
+    def _set_stateful(self, stateful_dict):
         """
         Set stateful values.
 
@@ -1220,8 +1294,14 @@ class Macro(DynamicStateful):
 
         """
         for element, attrs in stateful_dict.items():
-            for attr, value in attrs.items():
-                setattr(getattr(self.components, element), attr, value)
+            component = getattr(self.components, element)
+            if hasattr(component, '_set_stateful'):
+                component._set_stateful(attrs)
+            else:
+                [
+                    setattr(component, attr, value)
+                    for attr, value in attrs.items()
+                ]
 
     def _build_doc(self):
         """
@@ -1301,23 +1381,42 @@ class Model(Macro):
         the missing values, this option may cause the integration to
         fail, but it may be used to check the quality of the data.
 
+    See also
+    --------
+    :class:`pysd.py_backend.model.Macro`
+
     """
     def __init__(self, py_model_file, data_files, initialize, missing_values):
         """ Sets up the Python objects """
         super().__init__(py_model_file, None, None, Time(),
                          data_files=data_files)
-        self.time.stage = 'Load'
-        self.time.set_control_vars(**self.components._control_vars)
         self.data_files = data_files
         self.missing_values = missing_values
+        # set time component
+        self.time.stage = 'Load'
+        # set control var privately to do not change it when copying
+        self.time._set_control_vars(**self.components._control_vars)
+        # Attributes that are set later
         self.progress = None
         self.output = None
+        self.capture_elements = None
+        self.return_addresses = None
+        self._stepper_mode = None
+        self._submodel_tracker = {}
 
         if initialize:
             self.initialize()
 
     def initialize(self):
-        """ Initializes the simulation model """
+        """
+        Initializes the simulation model.
+
+        See also
+        --------
+        :func:`pysd.py_backend.model.Macro.initialize`
+        :func:`pysd.py_backend.model.Model.reload`
+
+        """
         self.time.stage = 'Initialization'
         External.missing = self.missing_values
         super().initialize()
@@ -1328,8 +1427,10 @@ class Model(Macro):
             cache_output=True, output_file=None):
         """
         Simulate the model's behavior over time.
-        Return a pandas dataframe with timestamps as rows,
-        model elements as columns.
+        Return a pandas dataframe with timestamps as rows and model
+        elements as columns. More detailed information and examples
+        of usage are available at
+        `Getting Started <https://pysd.readthedocs.io/en/master/getting_started.html#running-the-model>`__.
 
         Parameters
         ----------
@@ -1338,6 +1439,8 @@ class Model(Macro):
             Values are numeric or pandas Series.
             Numeric values represent constants over the model integration.
             Timeseries will be interpolated to give time-varying input.
+            For more information, check the documentation of
+            :func:`pysd.py_backend.model.Macro.set_components`.
 
         return_timestamps: list, numeric, ndarray (1D) (optional)
             Timestamps in model execution at which to return state information.
@@ -1351,13 +1454,15 @@ class Model(Macro):
 
         initial_condition: str or (float, dict) (optional)
             The starting time, and the state of the system (the values of
-            all the stocks) at that starting time. 'original' or 'o'uses
+            all the stocks) at that starting time. 'original' or 'o' uses
             model-file specified initial condition. 'current' or 'c' uses
             the state of the model after the previous execution. Other str
             objects, loads initial conditions from the pickle file with the
             given name.(float, dict) tuple lets the user specify a starting
             time (float) and (possibly partial) dictionary of initial values
             for stock (stateful) objects. Default is 'original'.
+            For more information, check the documentation of
+            :func:`pysd.py_backend.model.Model.set_initial_condition`
 
         final_time: float or None
             Final time of the simulation. If float, the given value will be
@@ -1412,11 +1517,11 @@ class Model(Macro):
         >>> model.run(return_timestamps=np.linspace(1, 10, 20))
         >>> model.run(output_file="results.nc")
 
-
-        See Also
+        See also
         --------
-        pysd.set_components : handles setting model parameters
-        pysd.set_initial_condition : handles setting initial conditions
+        :func:`pysd.py_backend.model.Macro.set_components`
+        :func:`pysd.py_backend.model.Model.set_initial_condition`
+        :func:`pysd.py_backend.model.Model.reload`
 
         """
         self._stepper_mode = False
@@ -1442,7 +1547,9 @@ class Model(Macro):
                     initial_condition='original', final_time=None,
                     time_step=None, saveper=None, cache_output=True):
         """
-        Configure the model stepping behavior.
+        Configure the model stepping behavior. Examples of usage are
+        available at
+        `Advanced Usage <https://pysd.readthedocs.io/en/master/advanced_usage.html#running-models-one-or-more-step-s-at-a-time>`__.
 
         Parameters
         ----------
@@ -1505,6 +1612,10 @@ class Model(Macro):
            recommended to activate this feature, if time step << saveper
            it is recommended to deactivate it. Default is True.
 
+        See also
+        --------
+        :func:`pysd.py_backend.model.Model.step`
+
         """
         self.output = output_obj
 
@@ -1523,6 +1634,8 @@ class Model(Macro):
         Run a model step. Updates model variables first (optional), and then
         runs any number of model steps. To collect the outputs after one or
         more steps, use the collect method of the ModelOutput class.
+        Examples of usage are available at
+        `Advanced Usage <https://pysd.readthedocs.io/en/master/advanced_usage.html#running-models-one-or-more-step-s-at-a-time>`__
 
         Parameters
         ----------
@@ -1538,9 +1651,13 @@ class Model(Macro):
         -------
         None
 
+        See also
+        --------
+        :func:`pysd.py_backend.model.Model.set_stepper`
+
         """
         # TODO warn the user if we exceeded the final_time??
-        self.set_components(step_vars)
+        self._set_components(step_vars, new=False)
 
         for _ in range(num_steps):
             self._integrate_step()
@@ -1554,7 +1671,8 @@ class Model(Macro):
         Internal method to set all simulation config parameters. Arguments
         to this function are those of the run and set_stepper methods.
         """
-
+        # set control var at the beginning in case they are needed to
+        # initialize any object
         self._set_control_vars(return_timestamps, final_time, time_step,
                                saveper)
 
@@ -1568,7 +1686,11 @@ class Model(Macro):
         # update cache types after setting params
         self._assign_cache_type()
 
+        # set initial conditions
         self.set_initial_condition(initial_condition)
+        # set control vars again in case a pickle has been used
+        self._set_control_vars(return_timestamps, final_time, time_step,
+                               saveper)
 
         # progressbar only makes sense when not running step by step
         if not self._stepper_mode:
@@ -1659,16 +1781,18 @@ class Model(Macro):
             if self.time.return_timestamps:
                 final_time = self.time.return_timestamps[0]
             else:
-                final_time = self.time.next_return
+                final_time = self.time._next_return
 
         self.time.set_control_vars(
             final_time=final_time, time_step=time_step, saveper=saveper)
 
-    def select_submodel(self, vars=[], modules=[], exogenous_components={}):
+    def select_submodel(self, vars=[], modules=[], exogenous_components={},
+                        inplace=True):
         """
         Select a submodel from the original model. After selecting a submodel
         only the necessary stateful objects for integrating this submodel will
-        be computed.
+        be computed. Examples of usage are available at
+        `Advanced Usage <https://pysd.readthedocs.io/en/master/advanced_usage.html#selecting-and-running-a-submodel>`__.
 
         Parameters
         ----------
@@ -1691,12 +1815,19 @@ class Model(Macro):
             set_components method. By default it is an empty dict and
             the needed exogenous components will be set to a numpy.nan value.
 
+        inplace: bool (optional)
+            If True it will modify current object and will return None.
+            If False it will create a copy of the model and return it
+            keeping the original model unchange. Default is True.
+
         Returns
         -------
-        None
+        None or pysd.py_backend.model.Model
+            If inplace=False it will return a modified copy of the
+            original model.
 
-        Notes
-        -----
+        Note
+        ----
         modules can be only passed when the model has been split in
         different files during translation.
 
@@ -1728,7 +1859,23 @@ class Model(Macro):
         ...         "initial_value_stock3": 5})
         UserWarning: Selecting submodel, to run the full model again use model.reload()
 
+        See also
+        --------
+        :func:`pysd.py_backend.model.Model.get_vars_in_module`
+        :func:`pysd.py_backend.model.Model.get_dependencies`
+
         """
+        if inplace:
+            self._select_submodel(vars, modules, exogenous_components)
+        else:
+            return self.copy()._select_submodel(
+                vars, modules, exogenous_components)
+
+    def _select_submodel(self, vars, modules, exogenous_components={}):
+        self._submodel_tracker = {
+            "vars": vars,
+            "modules": modules
+        }
         deps = self.get_dependencies(vars, modules)
         warnings.warn(
             "Selecting submodel, "
@@ -1814,6 +1961,7 @@ class Model(Macro):
         # re-assign the cache_type and initialization order
         self._assign_cache_type()
         self._get_initialize_order()
+        return self
 
     def get_dependencies(self, vars=[], modules=[]):
         """
@@ -1837,8 +1985,8 @@ class Model(Macro):
         dependencies: pysd.py_backend.utils.Dependencies
             Dependencies data object.
 
-        Notes
-        -----
+        Note
+        ----
         modules can be only passed when the model has been split in
         different files during translation.
 
@@ -1872,6 +2020,10 @@ class Model(Macro):
         Stateful objects integrated with the selected variables (total 1):
             _integ_stock1, _integ_stock3, _delay_fixed_delay1
 
+        See also
+        --------
+        :func:`pysd.py_backend.model.Model.get_vars_in_module`
+
         """
         def check_dep(deps_obj, deps, initial=False):
             for dep in deps:
@@ -1928,6 +2080,10 @@ class Model(Macro):
         vars: set
             Set of varible names in the given module.
 
+        See also
+        --------
+        :func:`pysd.py_backend.model.Model.get_dependencies`
+
         """
         if self._modules:
             module_content = self._modules.copy()
@@ -1958,10 +2114,74 @@ class Model(Macro):
 
         return vars
 
+    def copy(self, reload=False):
+        """
+        Create a copy of the current model.
+
+        Parameters
+        ----------
+        reload: bool (optional)
+            If True the model will be copied without applying to it any
+            change, the copy will simply load the model again from the
+            translated file. This would be equivalent to doing
+            :py:func:`pysd.load` with the same arguments. Otherwise, it
+            will apply the same changes that have been applied to the
+            original model and update the states (faithful copy).
+            Default is False.
+
+        Warning
+        -------
+        The copy function will load a new model from the file and apply
+        the same changes to it. If any of these changes have replaced a
+        variable with a function that references other variables in the
+        model, the copy will not work properly since the function will
+        still reference the variables in the original model, in which
+        case the function should be redefined.
+
+        See also
+        --------
+        :func:`pysd.py_backend.model.Model.reload`
+
+        """
+        # initialize the new model?
+        initialize = self.time.stage != 'Load'
+        # create a new model
+        new_model = type(self)(
+            py_model_file=deepcopy(self.py_model_file),
+            data_files=deepcopy(self.data_files),
+            initialize=initialize,
+            missing_values=deepcopy(self.missing_values)
+        )
+        if reload:
+            # return reloaded copy
+            return new_model
+        # copy the values of the stateful objects
+        if initialize:
+            new_model._set_stateful(deepcopy(super().export()))
+        # copy time object values
+        new_model.time._set_time(deepcopy(self.time.export()))
+        # set other components
+        with warnings.catch_warnings():
+            # filter warnings that have been already shown in original model
+            warnings.simplefilter("ignore")
+            # substract submodel
+            if self._submodel_tracker:
+                new_model._select_submodel(**self._submodel_tracker)
+            # copy modified parameters
+            new_model.set_components(self._components_setter_tracker)
+
+        return new_model
+
     def reload(self):
         """
         Reloads the model from the translated model file, so that all the
         parameters are back to their original value.
+
+        See also
+        --------
+        :func:`pysd.py_backend.model.Model.copy`
+        :func:`pysd.py_backend.model.Model.initialize`
+
         """
         self.__init__(self.py_model_file, data_files=self.data_files,
                       initialize=True,
@@ -2044,9 +2264,9 @@ class Model(Macro):
         >>> model.set_initial_condition('exported_pickle.pic')
         >>> model.set_initial_condition((10, {'teacup_temperature': 50}))
 
-        See Also
+        See also
         --------
-        model.set_initial_value()
+        :func:`pysd.py_backend.model.Macro.set_initial_value`
 
         """
         if isinstance(initial_condition, str)\
@@ -2059,6 +2279,7 @@ class Model(Macro):
             self.set_initial_value(*initial_condition)
         elif isinstance(initial_condition, Path):
             self.import_pickle(initial_condition)
+            self.time.set_control_vars(initial_time=self.time())
         elif isinstance(initial_condition, str):
             if initial_condition.lower() in ["original", "o"]:
                 self.time.set_control_vars(
@@ -2118,7 +2339,62 @@ class Model(Macro):
         progressbar.finish()
 
     def _integrate_step(self):
-
         self._euler_step(self.time.time_step())
         self.time.update(self.time()+self.time.time_step())
         self.clean_caches()
+
+    def export(self, file_name):
+        """
+        Export stateful values to pickle file.
+
+        Parameters
+        ----------
+        file_name: str or pathlib.Path
+          Name of the file to export the values.
+
+        See also
+        --------
+        :func:`pysd.py_backend.model.Model.import_pickle`
+
+        """
+        warnings.warn(
+            "\nCompatibility of exported states could be broken between"
+            " different versions of PySD or xarray, current versions:\n"
+            f"\tPySD {__version__}\n\txarray {xr.__version__}\n"
+        )
+        with open(file_name, 'wb') as file:
+            pickle.dump(
+                (self.time.export(),
+                 super().export(),
+                 {'pysd': __version__, 'xarray': xr.__version__}
+                 ), file)
+
+    def import_pickle(self, file_name):
+        """
+        Import stateful values from pickle file.
+
+        Parameters
+        ----------
+        file_name: str or pathlib.Path
+          Name of the file to import the values from.
+
+        See also
+        --------
+        :func:`pysd.py_backend.model.Model.export_pickle`
+
+        """
+        with open(file_name, 'rb') as file:
+            time_dict, stateful_dict, metadata = pickle.load(file)
+
+        if __version__ != metadata['pysd']\
+           or xr.__version__ != metadata['xarray']:  # pragma: no cover
+            warnings.warn(
+                "\nCompatibility of exported states could be broken between"
+                " different versions of PySD or xarray. Current versions:\n"
+                f"\tPySD {__version__}\n\txarray {xr.__version__}\n"
+                "Loaded versions:\n"
+                f"\tPySD {metadata['pysd']}\n\txarray {metadata['xarray']}\n"
+                )
+
+        self.time._set_time(time_dict)
+        self._set_stateful(stateful_dict)
diff --git a/pysd/py_backend/output.py b/pysd/py_backend/output.py
index abfd845..18ed23b 100644
--- a/pysd/py_backend/output.py
+++ b/pysd/py_backend/output.py
@@ -10,6 +10,7 @@ import time as t
 
 from csv import QUOTE_NONE
 from pathlib import Path
+from collections import defaultdict
 
 import regex as re
 
@@ -46,8 +47,7 @@ class ModelOutput():
         self.handler = DataFrameHandler(DatasetHandler(None)).handle(out_file)
 
     def set_capture_elements(self, capture_elements):
-        self.handler.capture_elements_step = capture_elements["step"] + \
-            ["time"]
+        self.handler.capture_elements_step = capture_elements["step"]
         self.handler.capture_elements_run = capture_elements["run"]
 
     def initialize(self, model):
@@ -214,23 +214,9 @@ class DatasetHandler(OutputHandlerInterface):
         super().__init__(next)
         self.out_file = None
         self.ds = None
-        self._step = 0
+        self.__step = 0
         self.nc = __import__("netCDF4")
 
-    @property
-    def step(self):
-        """
-        Used as time index for the output Dataset. Increases by one
-        at each iteration.
-        """
-        return self._step
-
-    def __update_step(self):
-        """
-        Increases the _step attribute by 1 at each model iteration.
-        """
-        self._step = self.step + 1
-
     def process_output(self, out_file):
         """
         If out_file can be handled by this concrete handler, it returns
@@ -267,6 +253,7 @@ class DatasetHandler(OutputHandlerInterface):
         None
 
         """
+        self.__step = 0
         self.ds = self.nc.Dataset(self.out_file, "w")
 
         # defining global attributes
@@ -295,7 +282,7 @@ class DatasetHandler(OutputHandlerInterface):
         # creating the time dimension as unlimited
         self.ds.createDimension("time", None)
         # creating variables
-        self.__create_ds_vars(model, self.capture_elements_step)
+        self.__create_ds_vars(model, self.capture_elements_step + ['time'])
 
     def update(self, model):
         """
@@ -312,19 +299,20 @@ class DatasetHandler(OutputHandlerInterface):
         None
 
         """
+        self.ds['time'][self.__step] = model.time.round()
         for key in self.capture_elements_step:
             comp = model[key]
             if isinstance(comp, xr.DataArray):
-                self.ds[key][self.step, :] = comp.values
+                self.ds[key][self.__step, :] = comp.values
             else:
-                self.ds[key][self.step] = comp
+                self.ds[key][self.__step] = comp
 
-        self.__update_step()
+        self.__step += 1
 
     def __update_run_elements(self, model):
         """
-        Writes values of cache run elements from the cature_elements set
-        in the netCDF4 Dataset.
+        Writes values of cache run elements from the capture_elements
+        set in the netCDF4 Dataset.
         Cache run elements do not have the time dimension.
 
         Parameters
@@ -429,6 +417,7 @@ class DataFrameHandler(OutputHandlerInterface):
         super().__init__(next)
         self.ds = None
         self.out_file = None
+        self.__step = 0
 
     def process_output(self, out_file):
         """
@@ -457,7 +446,7 @@ class DataFrameHandler(OutputHandlerInterface):
 
     def initialize(self, model):
         """
-        Creates a pandas DataFrame and adds model variables as columns.
+        Creates an empty dictionary to save the outputs.
 
         Parameters
         ----------
@@ -469,12 +458,12 @@ class DataFrameHandler(OutputHandlerInterface):
         None
 
         """
-        self.ds = pd.DataFrame(columns=self.capture_elements_step)
+        self.ds = defaultdict(list)
+        self.__step = 0
 
     def update(self, model):
         """
-        Add a row to the results pandas DataFrame with the values of the
-        variables listed in capture_elements.
+        Add new values to the data dictionary.
 
         Parameters
         ----------
@@ -486,13 +475,15 @@ class DataFrameHandler(OutputHandlerInterface):
         None
 
         """
-        self.ds.loc[model.time.round()] = [
-            getattr(model.components, key)()
-            for key in self.capture_elements_step]
+        self.ds['time'].append(model.time.round())
+        for key in self.capture_elements_step:
+            self.ds[key].append(getattr(model.components, key)())
+
+        self.__step += 1
 
     def postprocess(self, **kwargs):
         """
-        Delete time column from the pandas DataFrame and flatten
+        Convert the output dictionary to a pandas DataFrame and flatten
         xarrays if required.
 
         Returns
@@ -501,15 +492,15 @@ class DataFrameHandler(OutputHandlerInterface):
             Simulation results stored as a pandas DataFrame.
 
         """
-        # delete time column as it was created only for avoiding errors
-        # of appending data. See previous TODO.
-        del self.ds["time"]
+        # create the dataframe
+        df = pd.DataFrame.from_dict(self.ds)
+        df.set_index('time', inplace=True)
 
         # enforce flattening if df is to be saved to csv or tab file
         flatten = True if self.out_file else kwargs.get("flatten", None)
 
         df = DataFrameHandler.make_flat_df(
-            self.ds, kwargs["return_addresses"], flatten
+            df, kwargs["return_addresses"], flatten
             )
         if self.out_file:
             NCFile.df_to_text_file(df, self.out_file)
@@ -518,7 +509,7 @@ class DataFrameHandler(OutputHandlerInterface):
 
     def add_run_elements(self, model):
         """
-        Adds constant elements to a dataframe.
+        Adds constant elements to the output data dictionary.
 
         Parameters
         ----------
@@ -530,9 +521,8 @@ class DataFrameHandler(OutputHandlerInterface):
         None
 
         """
-        nx = len(self.ds.index)
-        for element in self.capture_elements_run:
-            self.ds[element] = [getattr(model.components, element)()] * nx
+        for key in self.capture_elements_run:
+            self.ds[key] = [getattr(model.components, key)()] * self.__step
 
     @staticmethod
     def make_flat_df(df, return_addresses, flatten=False):
diff --git a/pysd/py_backend/statefuls.py b/pysd/py_backend/statefuls.py
index b987059..41e490e 100644
--- a/pysd/py_backend/statefuls.py
+++ b/pysd/py_backend/statefuls.py
@@ -384,6 +384,9 @@ class DelayFixed(DynamicStateful):
         # The order can only be set once
         self.order = round(order + SMALL_VENSIM)
 
+        # set the pointer to 0
+        self.pointer = 0
+
         if init_val is None:
             init_state_value = self.init_func()
         else:


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index fe3f39d..0cbe8b5 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -46,3 +46,12 @@ def model(_root, tmp_path, model_path):
         return load(new_path)
     else:
         return ValueError("Invalid model")
+
+
+@pytest.fixture(scope="session")
+def ignore_warns():
+    # warnings to be ignored in the integration tests
+    return [
+        "numpy.ndarray size changed, may indicate binary incompatibility.",
+        "Creating an ndarray from ragged nested sequences.*"
+    ]
diff --git a/tests/pytest_integration/pytest_integration_vensim_pathway.py b/tests/pytest_integration/pytest_integration_vensim_pathway.py
index a4afd46..de24aac 100644
--- a/tests/pytest_integration/pytest_integration_vensim_pathway.py
+++ b/tests/pytest_integration/pytest_integration_vensim_pathway.py
@@ -1,10 +1,9 @@
-import warnings
+import re
 import shutil
 import pytest
 
 from pysd.tools.benchmarking import runner, assert_frames_close
 
-# TODO add warnings catcher per test
 
 vensim_test = {
     "abs": {
@@ -22,7 +21,10 @@ vensim_test = {
     "allocate_available": {
         "folder": "allocate_available",
         "file": "test_allocate_available.mdl",
-        "rtol": 2e-2
+        "rtol": 2e-2,
+        "warns": [
+            "Variable '.*' is defined with different .*types: '.*'",
+        ]
     },
     "allocate_by_priority": {
         "folder": "allocate_by_priority",
@@ -39,7 +41,10 @@ vensim_test = {
     "arguments": {
         "folder": "arguments",
         "file": "test_arguments.mdl",
-        "rtol": 1e-2  # TODO test why it is failing with smaller tolerance
+        "rtol": 1e-2,  # TODO test why it is failing with smaller tolerance
+        "warns": [
+            "_delay.*\nDelay time very small, casting delay order from .*"
+        ]
     },
     "array_with_line_break": {
         "folder": "array_with_line_break",
@@ -76,7 +81,10 @@ vensim_test = {
     },
     "delay_fixed": {
         "folder": "delay_fixed",
-        "file": "test_delay_fixed.mdl"
+        "file": "test_delay_fixed.mdl",
+        "warns": [
+            ".*\nCasting delay order from.*"
+        ]
     },
     "delay_numeric_error": {
         "folder": "delay_numeric_error",
@@ -88,7 +96,11 @@ vensim_test = {
     },
     "delay_pipeline": {
         "folder": "delay_pipeline",
-        "file": "test_pipeline_delays.mdl"
+        "file": "test_pipeline_delays.mdl",
+        "warns": [
+            "_delay.*\nDelay time very small, casting delay order from .*",
+            ".*\nCasting delay order from.*"
+        ]
     },
     "delays": {
         "folder": "delays",
@@ -116,7 +128,10 @@ vensim_test = {
     },
     "except_subranges": {
         "folder": "except_subranges",
-        "file": "test_except_subranges.mdl"
+        "file": "test_except_subranges.mdl",
+        "warns": [
+            "Variable '.*' is defined with different .*types: '.*'",
+        ]
     },
     "exp": {
         "folder": "exp",
@@ -148,7 +163,10 @@ vensim_test = {
     },  marks=pytest.mark.xfail(reason="csv files not implemented")),
     "get_constants_subranges": {
         "folder": "get_constants_subranges",
-        "file": "test_get_constants_subranges.mdl"
+        "file": "test_get_constants_subranges.mdl",
+        "warns": [
+            "Variable '.*' is defined with different .*types: '.*'",
+        ]
     },
     "get_data": pytest.param({
         "folder": "get_data",
@@ -164,7 +182,10 @@ vensim_test = {
     },
     "get_lookups_subscripted_args": {
         "folder": "get_lookups_subscripted_args",
-        "file": "test_get_lookups_subscripted_args.mdl"
+        "file": "test_get_lookups_subscripted_args.mdl",
+        "warns": [
+            "_ext.*\nextrapolating data .* the m.*mum value of the .*",
+        ]
     },
     "get_lookups_subset": {
         "folder": "get_lookups_subset",
@@ -172,7 +193,11 @@ vensim_test = {
     },
     "get_mixed_definitions": {
         "folder": "get_mixed_definitions",
-        "file": "test_get_mixed_definitions.mdl"
+        "file": "test_get_mixed_definitions.mdl",
+        "warns": [
+            "Variable '.*' is defined with different .*types: '.*'",
+            "_ext.*\nextrapolating data .* the m.*mum value of the .*",
+        ]
     },
     "get_subscript_3d_arrays_xls": {
         "folder": "get_subscript_3d_arrays_xls",
@@ -188,7 +213,11 @@ vensim_test = {
     },
     "get_with_missing_values_xlsx": {
         "folder": "get_with_missing_values_xlsx",
-        "file": "test_get_with_missing_values_xlsx.mdl"
+        "file": "test_get_with_missing_values_xlsx.mdl",
+        "warns": [
+            "_ext.*\nD.* value missing or non-valid in:.*",
+            "_ext.*\nextrapolating data .* the m.*mum value of the .*",
+        ]
     },
     "get_xls_cellrange": {
         "folder": "get_xls_cellrange",
@@ -308,7 +337,11 @@ vensim_test = {
     },
     "odd_number_quotes": {
         "folder": "odd_number_quotes",
-        "file": "teacup_3quotes.mdl"
+        "file": "teacup_3quotes.mdl",
+        "warns": [
+            "No encoding specified or detected to translate the model "
+            "file. 'UTF-8' encoding will be used."
+        ]
     },
     "parentheses": {
         "folder": "parentheses",
@@ -316,7 +349,10 @@ vensim_test = {
     },
     "partial_range_definitions": {
         "folder": "partial_range_definitions",
-        "file": "test_partial_range_definitions.mdl"
+        "file": "test_partial_range_definitions.mdl",
+        "warns": [
+            "\nDimension given by subscripts.*"
+        ]
     },
     "power": {
         "folder": "power",
@@ -324,7 +360,11 @@ vensim_test = {
     },
     "reality_checks": {
         "folder": "reality_checks",
-        "file": "test_reality_checks.mdl"
+        "file": "test_reality_checks.mdl",
+        "warns": [
+            "':.*:' detected. The expression content is not parsed and "
+            "will be ignored.",
+        ]
     },
     "reference_capitalization": {
         "folder": "reference_capitalization",
@@ -332,7 +372,10 @@ vensim_test = {
     },
     "repeated_subscript": {
         "folder": "repeated_subscript",
-        "file": "test_repeated_subscript.mdl"
+        "file": "test_repeated_subscript.mdl",
+        "warns": [
+            "\nAdding new subscript range to subscript_dict:.*"
+        ]
     },
     "rounding": {
         "folder": "rounding",
@@ -344,7 +387,12 @@ vensim_test = {
     },
     "smaller_range": {
         "folder": "smaller_range",
-        "file": "test_smaller_range.mdl"
+        "file": "test_smaller_range.mdl",
+        "warns": [
+            "\nDimension given by subscripts.*",
+            "\nAdding new subscript range to subscript_dict:.*",
+            "Variable '.*' is defined with different .*types: '.*'",
+        ]
     },
     "smooth": {
         "folder": "smooth",
@@ -364,7 +412,10 @@ vensim_test = {
     },
     "subrange_merge": {
         "folder": "subrange_merge",
-        "file": "test_subrange_merge.mdl"
+        "file": "test_subrange_merge.mdl",
+        "warns": [
+            "Variable '.*' is defined with different .*types: '.*'",
+        ]
     },
     "subscript_1d_arrays": {
         "folder": "subscript_1d_arrays",
@@ -444,7 +495,10 @@ vensim_test = {
     },
     "subscript_mixed_assembly": {
         "folder": "subscript_mixed_assembly",
-        "file": "test_subscript_mixed_assembly.mdl"
+        "file": "test_subscript_mixed_assembly.mdl",
+        "warns": [
+            "Variable '.*' is defined with different .*types: '.*'",
+        ]
     },
     "subscript_multiples": {
         "folder": "subscript_multiples",
@@ -496,7 +550,10 @@ vensim_test = {
     },
     "subscripted_lookups": {
         "folder": "subscripted_lookups",
-        "file": "test_subscripted_lookups.mdl"
+        "file": "test_subscripted_lookups.mdl",
+        "warns": [
+            "_.*\nextrapolating data .* the m.*mum value of the .*",
+        ]
     },
     "subscripted_ramp_step": {
         "folder": "subscripted_ramp_step",
@@ -504,7 +561,10 @@ vensim_test = {
     },
     "subscripted_round": {
         "folder": "subscripted_round",
-        "file": "test_subscripted_round.mdl"
+        "file": "test_subscripted_round.mdl",
+        "warns": [
+            "Variable '.*' is defined with different .*types: '.*'",
+        ]
     },
     "subscripted_smooth": {
         "folder": "subscripted_smooth",
@@ -560,7 +620,11 @@ vensim_test = {
     },
     "vector_select": {
         "folder": "vector_select",
-        "file": "test_vector_select.mdl"
+        "file": "test_vector_select.mdl",
+        "warns": [
+            "Vensim's help says that numerical_action=5 computes the "
+            "product of selection_array . expression_array..*"
+        ]
     },
     "with_lookup": {
         "folder": "with_lookup",
@@ -623,6 +687,11 @@ class TestIntegrateVensim:
         else:
             return None
 
+    @pytest.fixture
+    def warns(self, test_data, ignore_warns):
+        ign_warns = test_data.get('warns', []) + ignore_warns
+        return [re.compile(w) for w in ign_warns]
+
     @pytest.fixture
     def kwargs(self, test_data):
         """Fixture for atol and rtol"""
@@ -633,8 +702,12 @@ class TestIntegrateVensim:
             kwargs["rtol"] = test_data["rtol"]
         return kwargs
 
-    def test_read_vensim_file(self, model_path, data_path, kwargs):
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore")
-            output, canon = runner(model_path, data_files=data_path)
+    def test_read_vensim_file(self, model_path, data_path, kwargs,
+                              recwarn, warns):
+        output, canon = runner(model_path, data_files=data_path)
+        for warn in recwarn:
+            warn = str(warn.message)
+            assert any([re.match(pwarn, warn) for pwarn in warns]), \
+                f"Couldn't match warning:\n{warn}"
+
         assert_frames_close(output, canon, verbose=True, **kwargs)
diff --git a/tests/pytest_integration/pytest_integration_xmile_pathway.py b/tests/pytest_integration/pytest_integration_xmile_pathway.py
index d452a44..a422bb8 100644
--- a/tests/pytest_integration/pytest_integration_xmile_pathway.py
+++ b/tests/pytest_integration/pytest_integration_xmile_pathway.py
@@ -1,10 +1,9 @@
-import warnings
+import re
 import shutil
 import pytest
 
 from pysd.tools.benchmarking import runner, assert_frames_close
 
-# TODO add warnings catcher per test
 
 xmile_test = {
     "abs": {
@@ -187,10 +186,10 @@ xmile_test = {
         "folder": "smooth_and_stock",
         "file": "test_smooth_and_stock.xmile"
     },  marks=pytest.mark.xfail(reason="failing originally")),
-    "special_characters": pytest.param({
-        "folder": "special_characters",
+    "special_characters": {
+        "folder": "special_characters_xmile",
         "file": "test_special_variable_names.xmile"
-    }, marks=pytest.mark.xfail(reason="failing originally")),
+    },
     "sqrt": {
         "folder": "sqrt",
         "file": "test_sqrt.xmile"
@@ -278,6 +277,11 @@ class TestIntegrateXmile:
         """Return model path"""
         return test_folder.joinpath(test_data["file"])
 
+    @pytest.fixture
+    def warns(self, test_data, ignore_warns):
+        ign_warns = test_data.get('warns', []) + ignore_warns
+        return [re.compile(w) for w in ign_warns]
+
     @pytest.fixture
     def kwargs(self, test_data):
         """Fixture for atol and rtol"""
@@ -288,8 +292,11 @@ class TestIntegrateXmile:
             kwargs["rtol"] = test_data["rtol"]
         return kwargs
 
-    def test_read_vensim_file(self, model_path, kwargs):
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore")
-            output, canon = runner(model_path)
+    def test_read_xmile_file(self, model_path, kwargs, recwarn, warns):
+        output, canon = runner(model_path)
+        for warn in recwarn:
+            warn = str(warn.message)
+            assert any([re.match(pwarn, warn) for pwarn in warns]), \
+                f"Couldn't match warning:\n{warn}"
+
         assert_frames_close(output, canon, **kwargs)
diff --git a/tests/pytest_pysd/pytest_output.py b/tests/pytest_pysd/pytest_output.py
index 341524b..36c4fce 100644
--- a/tests/pytest_pysd/pytest_output.py
+++ b/tests/pytest_pysd/pytest_output.py
@@ -312,15 +312,6 @@ class TestOutput():
         output.set_capture_elements(capture_elements)
         output.initialize(model)
 
-        # Dataset handler step cannot be modified from the outside
-        with pytest.raises(AttributeError):
-            output.handler.step = 5
-
-        with pytest.raises(AttributeError):
-            output.handler.__update_step()
-
-        assert output.handler.step == 0
-
     def test_make_flat_df(self):
 
         df = pd.DataFrame(index=[1], columns=['elem1'])
diff --git a/tests/pytest_pysd/pytest_pysd.py b/tests/pytest_pysd/pytest_pysd.py
index 3dfd434..cb5b85f 100644
--- a/tests/pytest_pysd/pytest_pysd.py
+++ b/tests/pytest_pysd/pytest_pysd.py
@@ -22,6 +22,8 @@ test_model_look = Path(
     + "test_get_lookups_subscripted_args.mdl")
 test_model_data = Path(
     "test-models/tests/get_data_args_3d_xls/test_get_data_args_3d_xls.mdl")
+test_macro_stock = Path("test-models/tests/macro_stock/test_macro_stock.mdl")
+test_model_query = Path("test-models/samples/Query_file/Query_file.mdl")
 
 more_tests = Path("more-tests")
 
@@ -31,6 +33,7 @@ test_model_constant_pipe = more_tests.joinpath(
 test_model_stepper = more_tests.joinpath(
     "stepper_cache/stepper_cache.mdl")
 
+
 class TestPySD():
 
     @pytest.mark.parametrize("model_path", [test_model])
@@ -306,11 +309,14 @@ class TestPySD():
             index=timeseries,
             data=(50 + np.random.rand(len(timeseries)).cumsum())
         )
-        res = model.run(
-            params={"room_temperature": temp_timeseries},
-            return_columns=["room_temperature"],
-            return_timestamps=timeseries,
-        )
+        warn_message = "Replacing a constant value with a time-dependent "\
+                       "value. The value will be interpolated over time."
+        with pytest.warns(UserWarning, match=warn_message):
+            res = model.run(
+                params={"room_temperature": temp_timeseries},
+                return_columns=["room_temperature"],
+                return_timestamps=timeseries,
+            )
         assert (res["room_temperature"] == temp_timeseries).all()
 
     @pytest.mark.parametrize("model_path", [test_model])
@@ -320,7 +326,11 @@ class TestPySD():
             index=timeseries,
             data=(50 + np.random.rand(len(timeseries)).cumsum())
         )
-        model.components.room_temperature = temp_timeseries
+        warn_message = "Replacing a constant value with a time-dependent "\
+                       "value. The value will be interpolated over time."
+        with pytest.warns(UserWarning, match=warn_message):
+            model.components.room_temperature = temp_timeseries
+
         res = model.run(
             return_columns=["room_temperature"],
             return_timestamps=timeseries,
@@ -349,8 +359,10 @@ class TestPySD():
     def test_set_components_with_function(self, model):
         def test_func():
             return 5
-
-        model.set_components({"Room Temperature": test_func})
+        warn_message = "Replacing a constant value with a callable. "\
+                       "The value may not be constant anymore."
+        with pytest.warns(UserWarning, match=warn_message):
+            model.set_components({"Room Temperature": test_func})
         res = model.run(return_columns=["Room Temperature"])
         assert test_func() == res["Room Temperature"].iloc[0]
 
@@ -546,12 +558,15 @@ class TestPySD():
         xr_series = [xr.DataArray(val, coords, dims) for val in val_series]
 
         temp_timeseries = pd.Series(index=timeseries, data=val_series)
-        res = model.run(
-            params={"initial_values": temp_timeseries, "final_time": 10},
-            return_columns=["initial_values"],
-            return_timestamps=timeseries,
-            flatten_output=False
-        )
+        warn_message = "Replacing a constant value with a time-dependent "\
+                       "value. The value will be interpolated over time."
+        with pytest.warns(UserWarning, match=warn_message):
+            res = model.run(
+                params={"initial_values": temp_timeseries, "final_time": 10},
+                return_columns=["initial_values"],
+                return_timestamps=timeseries,
+                flatten_output=False
+            )
 
         assert np.all(
             [
@@ -580,8 +595,11 @@ class TestPySD():
                                                               ).cumsum()]
         temp_timeseries = pd.Series(index=timeseries, data=val_series)
         out_series = [out_b + val for val in val_series]
-        model.set_components({"initial_values": temp_timeseries,
-                              "final_time": 10})
+        warn_message = "Replacing a constant value with a time-dependent "\
+                       "value. The value will be interpolated over time."
+        with pytest.warns(UserWarning, match=warn_message):
+            model.set_components({"initial_values": temp_timeseries,
+                                  "final_time": 10})
         res = model.run(
             return_columns=["initial_values"], flatten_output=False)
         assert np.all(
@@ -607,12 +625,15 @@ class TestPySD():
             data=[init_val + rd for rd in np.random.rand(len(timeseries)
                                                          ).cumsum()],
         )
-        res = model.run(
-            params={"initial_values": temp_timeseries, "final_time": 10},
-            return_columns=["initial_values"],
-            return_timestamps=timeseries,
-            flatten_output=False
-        )
+        warn_message = "Replacing a constant value with a time-dependent "\
+                       "value. The value will be interpolated over time."
+        with pytest.warns(UserWarning, match=warn_message):
+            res = model.run(
+                params={"initial_values": temp_timeseries, "final_time": 10},
+                return_columns=["initial_values"],
+                return_timestamps=timeseries,
+                flatten_output=False
+            )
 
         assert np.all(
             [
@@ -944,10 +965,24 @@ class TestPySD():
         with pytest.raises(TypeError):
             model.set_initial_value(new_time + 2, {'_integ_stock_a': input3})
 
+    @pytest.mark.parametrize("model_path", [test_model])
+    def test_replace_stateful(self, model):
+        warn_message = "Replacing the value of Stateful variable with "\
+                       "an expression. To set initial conditions use "\
+                       "`set_initial_condition` instead..."
+        with pytest.warns(UserWarning, match=warn_message):
+            model.components.teacup_temperature = 3
+
+        stocks = model.run()
+        assert np.all(stocks["Teacup Temperature"] == 3)
+
     @pytest.mark.parametrize("model_path", [test_model])
     def test_replace_element(self, model):
         stocks1 = model.run()
-        model.components.characteristic_time = lambda: 3
+        warn_message = "Replacing a constant value with a callable. "\
+                       "The value may not be constant anymore."
+        with pytest.warns(UserWarning, match=warn_message):
+            model.components.characteristic_time = lambda: 3
         stocks2 = model.run()
         assert stocks1["Teacup Temperature"].loc[10]\
             > stocks2["Teacup Temperature"].loc[10]
@@ -964,12 +999,12 @@ class TestPySD():
         set_temp = model.components.teacup_temperature()
         set_time = model.components.time()
 
-        assert set_temp != initial_temp,\
+        assert set_temp != initial_temp, \
             "Test definition is wrong, please change configuration"
 
         assert set_temp == 500
 
-        assert initial_time != new_time,\
+        assert initial_time != new_time, \
             "Test definition is wrong, please change configuration"
         assert new_time == set_time
 
@@ -992,12 +1027,12 @@ class TestPySD():
         set_temp = model.components.teacup_temperature()
         set_time = model.components.time()
 
-        assert set_temp != initial_temp,\
+        assert set_temp != initial_temp, \
             "Test definition is wrong, please change configuration"
 
         assert set_temp == 500
 
-        assert initial_time != new_time,\
+        assert initial_time != new_time, \
             "Test definition is wrong, please change configuration"
         assert new_time == set_time
 
@@ -1020,12 +1055,12 @@ class TestPySD():
         set_temp = model.components.teacup_temperature()
         set_time = model.components.time()
 
-        assert set_temp != initial_temp,\
+        assert set_temp != initial_temp, \
             "Test definition is wrong, please change configuration"
 
         assert set_temp == 500
 
-        assert initial_time != 10,\
+        assert initial_time != 10, \
             "Test definition is wrong, please change configuration"
 
         assert set_time == 10
@@ -1087,27 +1122,27 @@ class TestPySD():
                     "Initial Values": {
                         "One Dimensional Subscript": ["Entry 1", "Entry 2",
                                                       "Entry 3"],
-                        "Second Dimension Subscript":["Column 1", "Column 2"]
+                        "Second Dimension Subscript": ["Column 1", "Column 2"]
                     },
                     "initial_values": {
                         "One Dimensional Subscript": ["Entry 1", "Entry 2",
                                                       "Entry 3"],
-                        "Second Dimension Subscript":["Column 1", "Column 2"]
+                        "Second Dimension Subscript": ["Column 1", "Column 2"]
                     },
                     "Stock A": {
                         "One Dimensional Subscript": ["Entry 1", "Entry 2",
                                                       "Entry 3"],
-                        "Second Dimension Subscript":["Column 1", "Column 2"]
+                        "Second Dimension Subscript": ["Column 1", "Column 2"]
                     },
                     "stock_a": {
                         "One Dimensional Subscript": ["Entry 1", "Entry 2",
                                                       "Entry 3"],
-                        "Second Dimension Subscript":["Column 1", "Column 2"]
+                        "Second Dimension Subscript": ["Column 1", "Column 2"]
                     },
                     "_integ_stock_a": {
                         "One Dimensional Subscript": ["Entry 1", "Entry 2",
                                                       "Entry 3"],
-                        "Second Dimension Subscript":["Column 1", "Column 2"]
+                        "Second Dimension Subscript": ["Column 1", "Column 2"]
                     }
                 }
             )
@@ -1228,9 +1263,9 @@ class TestPySD():
         model.output = output
         model._integrate()
         res = model.output.handler.ds
-        assert isinstance(res, pd.DataFrame)
+        assert isinstance(res, dict)
         assert 'teacup_temperature' in res
-        assert all(res.index.values == list(range(0, 5, 2)))
+        assert np.array_equal(res['time'], np.arange(0, 5, 2))
 
         model.reload()
         model.time.add_return_timestamps(list(range(0, 5, 2)))
@@ -1290,6 +1325,62 @@ class TestPySD():
         # Check mdl_file
         assert model.mdl_file == str(path)
 
+    @pytest.mark.parametrize("model_path", [test_model_query])
+    def test_performance_dataframe(self, model):
+        model.run()
+
+    @pytest.mark.parametrize("model_path", [test_model])
+    def test_copy(self, model):
+        # check copying the model and running one of them
+        model2 = model.copy()
+        assert model['room_temperature'] == model2['room_temperature']
+        assert model['room_temperature'] != 3
+        model2.components.room_temperature = 3
+        assert model['room_temperature'] != model2['room_temperature']
+        model.run()
+        assert model.time() == 30
+        assert np.isclose(model['teacup_temperature'], 75.374)
+        assert model2.time() == 0
+        assert np.isclose(model2['teacup_temperature'], 180)
+
+        # check that time is copied
+        model3 = model.copy()
+        assert model3.time() == 30
+        assert np.isclose(model3['teacup_temperature'], 75.374)
+
+        # check that changes in variables are done again
+        assert model['room_temperature'] == 70
+        model.components.room_temperature = 200
+        model4 = model.copy()
+        assert model4['room_temperature'] == 200
+
+        # check that changes in control vars are copied
+        model.run(initial_condition=(3, {}), final_time=5)
+        model5 = model.copy()
+        assert model5.time.initial_time() == 3
+        assert model5.time.final_time() == 5
+        assert model5.time.time_step() == model.time.time_step()
+        assert model5.time.saveper() == model.time.saveper()
+        model.run(time_step=3, saveper=6, final_time=16)
+        model6 = model.copy()
+        assert model6.time.initial_time() == model.time.initial_time()
+        assert model6.time.final_time() == 16
+        assert model6.time.time_step() == 3
+        assert model6.time.saveper() == 6
+
+        # check a reloaded model
+        modelr = model.copy(reload=True)
+        assert modelr['room_temperature'] == 70
+        assert modelr['teacup_temperature'] == 180
+        assert modelr.time() == 0
+
+    @pytest.mark.parametrize("model_path", [test_macro_stock])
+    def test_copy_macro(self, model):
+        # check that macro's statefuls are copied recursively
+        out1 = model.run(return_columns=['macro_output']).values[-1]
+        model2 = model.copy()
+        assert model2['macro_output'] == out1
+
 
 class TestModelInteraction():
     """ The tests in this class test pysd's interaction with itself
@@ -1334,7 +1425,10 @@ class TestModelInteraction():
         model_2 = pysd.read_vensim(
             _root.joinpath("test-models/samples/SIR/SIR.mdl"))
 
-        model_1.components.initial_time = lambda: 10
+        warn_message = "Replacing a constant value with a callable. "\
+                       "The value may not be constant anymore."
+        with pytest.warns(UserWarning, match=warn_message):
+            model_1.components.initial_time = lambda: 10
         assert model_2.components.initial_time != 10
 
         # check that the model time is not shared between the two objects
@@ -1473,9 +1567,13 @@ class TestDependencies():
 
         # we should ensure that the constant_cache is removed
         # when passing new param
-        out2 = model.run(params={"constant1": new_var})
+        warn_message = "Replacing a constant value with a "\
+                       "time-dependent value. The value will "\
+                       "be interpolated over time."
+        with pytest.warns(UserWarning, match=warn_message):
+            out2 = model.run(params={"constant1": new_var})
 
-        assert not np.all(out1 - out2 == 0)
+        assert not np.array_equal(out1, out2)
 
         for key in pipe:
             assert model.cache_type[key] == "step"
@@ -1536,10 +1634,15 @@ class TestExportImport():
                 Path('test-models/tests/initial_function/test_initial.mdl'),
                 ([8, 20], [8], [20]),
                 (None, 15, None)
+            ),
+            (
+                test_macro_stock,
+                ([8, 20], [8], [20]),
+                (None, 15, None)
             )
         ],
         ids=["integ", "delays", "delay_fixed", "forecast", "sample_if_true",
-             "smooth", "trend", "initial"]
+             "smooth", "trend", "initial", "macro_stock"]
     )
     @pytest.mark.filterwarnings("ignore")
     def test_run_export_import(self, tmp_path, model, return_ts, final_t):
@@ -1582,6 +1685,16 @@ class TestExportImport():
         if return_ts[2]:
             assert_frames_close(stocks2, stocks.loc[return_ts[2]])
 
+        # test that time component is exported
+        model.export(export_path)
+        model.reload()
+        model.import_pickle(export_path)
+        assert model.time._time == finals[2]
+        assert model.time.final_time() == finals[2]
+        assert model.time.initial_time() == finals[1]
+        assert model.time.stage == "Run"
+
+
 class TestStepper():
 
     @pytest.mark.parametrize("model_path", [test_model_stepper])
@@ -1603,7 +1716,7 @@ class TestStepper():
         steps = int((model["final_time"] - model["initial_time"]) /
                     model["time_step"])
 
-        for i in range(1,steps + 1):
+        for i in range(1, steps+1):
             model.step(1, {"foo": i})
 
         res = output.collect(model)
@@ -1630,5 +1743,6 @@ class TestStepper():
 
         result_df = output.collect(model)
 
-        np.array_equal(result_df["Room Temperature"].values, result_temperatures)
+        np.array_equal(
+            result_df["Room Temperature"].values, result_temperatures)
         assert np.floor(result_df.loc[5, "Teacup Temperature"]) == 144
diff --git a/tests/pytest_pysd/pytest_select_submodel.py b/tests/pytest_pysd/pytest_select_submodel.py
index 470fa89..fba7702 100644
--- a/tests/pytest_pysd/pytest_select_submodel.py
+++ b/tests/pytest_pysd/pytest_select_submodel.py
@@ -177,7 +177,10 @@ class TestSubmodel:
                 assert var in str(record[-1].message)
             assert np.any(np.isnan(model.run()))
             # redefine dependencies
-            assert not np.any(np.isnan(model.run(params=dep_vars)))
+            warn_message = "Replacing a variable by a constant value."
+            with pytest.warns(UserWarning, match=warn_message):
+                out = model.run(params=dep_vars)
+            assert not np.any(np.isnan(out))
 
         # select submodel using contour values
         model.reload()
@@ -187,6 +190,111 @@ class TestSubmodel:
 
         assert not np.any(np.isnan(model.run()))
 
+    def test_select_submodel_copy(self, model, variables, modules,
+                                  n_deps, dep_vars):
+
+        # assert original stateful elements
+        assert len(model._dynamicstateful_elements) == 2
+        assert "_integ_other_stock" in model._stateful_elements
+        assert "_integ_other_stock" in model._dependencies
+        assert "other_stock" in model._dependencies
+        assert "other stock" in model._namespace
+        assert "other stock" in model._doc["Real Name"].to_list()
+        assert "other_stock" in model._doc["Py Name"].to_list()
+        assert "_integ_stock" in model._stateful_elements
+        assert "_integ_stock" in model._dependencies
+        assert "stock" in model._dependencies
+        assert "Stock" in model._namespace
+        assert "Stock" in model._doc["Real Name"].to_list()
+        assert "stock" in model._doc["Py Name"].to_list()
+
+        # select submodel
+        with pytest.warns(UserWarning) as record:
+            model2 = model.select_submodel(vars=variables, modules=modules,
+                                           inplace=False)
+
+        # assert warning
+        assert str(record[0].message) == self.warning
+
+        # assert original stateful elements
+        assert len(model._dynamicstateful_elements) == 2
+        assert "_integ_other_stock" in model._stateful_elements
+        assert "_integ_other_stock" in model._dependencies
+        assert "other_stock" in model._dependencies
+        assert "other stock" in model._namespace
+        assert "other stock" in model._doc["Real Name"].to_list()
+        assert "other_stock" in model._doc["Py Name"].to_list()
+        assert "_integ_stock" in model._stateful_elements
+        assert "_integ_stock" in model._dependencies
+        assert "stock" in model._dependencies
+        assert "Stock" in model._namespace
+        assert "Stock" in model._doc["Real Name"].to_list()
+        assert "stock" in model._doc["Py Name"].to_list()
+
+        # assert stateful elements change
+        assert len(model2._dynamicstateful_elements) == 1
+        assert "_integ_other_stock" not in model2._stateful_elements
+        assert "_integ_other_stock" not in model2._dependencies
+        assert "other_stock" not in model2._dependencies
+        assert "other stock" not in model2._namespace
+        assert "other stock" not in model2._doc["Real Name"].to_list()
+        assert "other_stock" not in model2._doc["Py Name"].to_list()
+        assert "_integ_stock" in model2._stateful_elements
+        assert "_integ_stock" in model2._dependencies
+        assert "stock" in model2._dependencies
+        assert "Stock" in model2._namespace
+        assert "Stock" in model2._doc["Real Name"].to_list()
+        assert "stock" in model2._doc["Py Name"].to_list()
+
+        if not dep_vars:
+            # totally independent submodels can run without producing
+            # nan values
+            assert not np.any(np.isnan(model2.run()))
+        else:
+            # running the model without redefining dependencies will
+            # produce nan values
+            assert "Exogenous components for the following variables are"\
+                + " necessary but not given:" in str(record[-1].message)
+            assert "Please, set them before running the model using "\
+                + "set_components method..." in str(record[-1].message)
+            for var in dep_vars:
+                assert var in str(record[-1].message)
+            assert np.any(np.isnan(model2.run()))
+            # redefine dependencies
+            warn_message = "Replacing a variable by a constant value."
+            with pytest.warns(UserWarning, match=warn_message):
+                out = model2.run(params=dep_vars)
+            assert not np.any(np.isnan(out))
+
+        # select submodel using contour values
+        with pytest.warns(UserWarning) as record:
+            model2 = model.select_submodel(vars=variables, modules=modules,
+                                           exogenous_components=dep_vars,
+                                           inplace=False)
+
+        assert not np.any(np.isnan(model2.run()))
+
+        # copy of submodel
+        model3 = model2.copy()
+        # assert stateful elements change
+        assert len(model3._dynamicstateful_elements) == 1
+        assert "_integ_other_stock" not in model3._stateful_elements
+        assert "_integ_other_stock" not in model3._dependencies
+        assert "other_stock" not in model3._dependencies
+        assert "other stock" not in model3._namespace
+        assert "other stock" not in model3._doc["Real Name"].to_list()
+        assert "other_stock" not in model3._doc["Py Name"].to_list()
+        assert "_integ_stock" in model3._stateful_elements
+        assert "_integ_stock" in model3._dependencies
+        assert "stock" in model3._dependencies
+        assert "Stock" in model3._namespace
+        assert "Stock" in model3._doc["Real Name"].to_list()
+        assert "stock" in model3._doc["Py Name"].to_list()
+
+        # dep_vars should be already set
+        out = model3.run()
+        assert not np.any(np.isnan(out))
+
 
 @pytest.mark.parametrize(
     "model_path,split_views,module,raise_type,error_message",
diff --git a/tests/pytest_types/data/pytest_data_with_model.py b/tests/pytest_types/data/pytest_data_with_model.py
index 7f5e8cc..9f3c02b 100644
--- a/tests/pytest_types/data/pytest_data_with_model.py
+++ b/tests/pytest_types/data/pytest_data_with_model.py
@@ -120,10 +120,11 @@ class TestPySDData:
             expected)
 
     def test_modify_data(self, model, expected):
-        out = model.run(params={
-            "var1": pd.Series(index=[1, 3, 7], data=[10, 20, 30]),
-            "var2": 10
-        })
+        with pytest.warns(UserWarning, match="Replacing .*"):
+            out = model.run(params={
+                "var1": pd.Series(index=[1, 3, 7], data=[10, 20, 30]),
+                "var2": 10
+            })
 
         assert (out["var2"] == 10).all()
         assert (
diff --git a/tests/pytest_types/external/pytest_external.py b/tests/pytest_types/external/pytest_external.py
index 5c4b2ae..8ec5a18 100644
--- a/tests/pytest_types/external/pytest_external.py
+++ b/tests/pytest_types/external/pytest_external.py
@@ -66,9 +66,8 @@ class TestExcels():
         Excels.clean()
         assert list(Excels._Excels_opyxl) == []
 
-    @pytest.mark.skipif(
-        sys.platform.startswith("win"),
-        reason="not working on Windows"
+    @pytest.mark.skip(
+        reason="may fail in parallel running"
     )
     def test_close_file(self, _root):
         """
diff --git a/tests/pytest_types/statefuls/pytest_statefuls.py b/tests/pytest_types/statefuls/pytest_statefuls.py
index bd7e244..c7060f6 100644
--- a/tests/pytest_types/statefuls/pytest_statefuls.py
+++ b/tests/pytest_types/statefuls/pytest_statefuls.py
@@ -1,5 +1,6 @@
 import pytest
 
+import numpy as np
 import xarray as xr
 
 from pysd.py_backend.statefuls import Stateful, Integ, Delay, DelayN,\
@@ -96,6 +97,22 @@ class TestStateful():
         delay_c.initialize(6)
         assert delay_c() == 6
 
+        for p, i in enumerate(np.arange(0, 2.5, 0.5)):
+            delay_c.update(None)
+            assert delay_c() == 6
+            assert delay_c.pointer == p+1
+
+        delay_c.update(None)
+        assert delay_c() == 5
+        assert delay_c.pointer == 0
+
+        # check that the pointer is set to 0 after initialization
+        delay_c.update(None)
+        assert delay_c.pointer != 0
+        delay_c.initialize()
+        assert delay_c.pointer == 0
+
+
     def test_delay_subscript(self):
         """
         Test for subscripted delay
diff --git a/tests/test-models b/tests/test-models
index 87d980f..ac2edbc 160000
--- a/tests/test-models
+++ b/tests/test-models
@@ -1,1 +1,1 @@
-Subproject commit 87d980f61c730efffc829e565cf375795f5af36a
+Subproject commit ac2edbcde87aec2dfd8cc338d0574970c549a2d9