Repository: SDXorg/pysd (language: python)

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.

Improve memory allocation for saving results
Maybe the dataframe is not the best object to save the outputs.
We should look for the fastest way to do so the integration doesn't have to wait.

We could preallocate memory for output of size `(final_time-initial_time)/saveper` however, we may need to update the allocated space as `saveper` and `final_time `values may change during integration, we could also do something like `list.append` as[ space allocation is done in chunks](https://www.opensourceforu.com/2021/05/memory-management-in-lists-and-tuples/).

We could also save the data for the last `saveper` and launch a subprocess that saves the data in a file while the model continues integrating.

We should also consider saving the data in the same format the user had asked for (for example flattening it during integration if this doesn't take too much extra time). This would avoid iterating over the dataframe's columns after the integration.
DelayFixed returning a list assignment index out of range due to wrong pointer value
I tried to run .mdl file using pysd
```jupyter-notebook
#context: when I changed a few parameters in the model, it showed list assignment index out of range, but these values were working #fine in vensim model

!pip -q install pySD
import numpy as np
import pandas as pd
import pysd
import matplotlib.pyplot as plt

vmodel = pysd.read_vensim('Query_file.mdl')
SD_result = vmodel.run()
#the values of LT1, LT2 parameters in the original model are 30 and 3 respectively
#in the next step I tried to change the parameter (LT1, LT2) values to 60 and 6 respectively 

SD_result = vmodel.run(params={ 'LT1':60, 'LT2':6})

#I again tried to change the parameter values to 30 and 3, but it showed me an error: **list assignment index out of range**

SD_result = vmodel.run(params={ 'LT1':30, 'LT2':3})
``` 
[Query_file.zip](https://github.com/quaquel/EMAworkbench/files/13623164/Query_file.zip)
please let me know why is it happening with my code
I have attached vensim model in the zip file for reference





Discussion / hints:
And by 'hypothetical user' I really mean: I just did this, and couldn't figure out what was happening.

Hi @JamesPHoughton 
I am thinking about including some small improvements in 3.8.0 to close some old issues. I could probably include different warnings in params let me know if they suit you (I will probably include also information about the var names and an option to don't show the warning messages):

- When setting a Data o Lookup object to a constant "Warning: Replacing interpolation data with constant values."
- When setting a cache constant value with a type different than a constant (i.e. a constant that is computed with other constant values): "Warning: Replacing a constant value that is not primitive by another constant value."
- When replacing a constant with a pandas.Series "Warning: replacing a constant value with a time-dependent value. The value will be interpolated over time."
- When replacing a time-dependent variable with a constant or a pandas.Series "Warning: replacing the value of a variable by a constant/time dependent value"
- When replacing a stateful object with a constant or a pandas.Series "Warning: replacing the value of Stateful variable by a constant/time-dependent value. To set initial conditions use xxxxxx"
Move to numpy backend first #328
HI @sanketme 

Thank you for opening the issue. I will try to reproduce this bug and see if there is something that we need to correct.
Bug found!

`self.pointer` wasn't being set to 0 in the initialization. I will try to update it in the master soon.

Fix patch:
diff --git a/docs/advanced_usage.rst b/docs/advanced_usage.rst
index 0c8f7756..c583b194 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 b7ec3172..b1646125 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 55f918d3..3da4238f 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 f0506cec..4f524802 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 d1a7f1e0..62ee17b8 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 ce432835..019a69c5 100644
--- a/pysd/builders/python/python_model_builder.py
+++ b/pysd/builders/python/python_model_builder.py
@@ -663,7 +663,7 @@ def build_element(self) -> None:
         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 @@ def build_element(self) -> None:
         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 140ac9b2..8eb79d9b 100644
--- a/pysd/py_backend/components.py
+++ b/pysd/py_backend/components.py
@@ -7,6 +7,7 @@
 import random
 import inspect
 import importlib.util
+from copy import deepcopy
 
 import numpy as np
 
@@ -129,13 +130,45 @@ def __init__(self):
         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 @@ def set_control_vars(self, **kwargs):
             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 @@ def _convert_value(value):
                 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 @@ def in_return(self):
 
         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 @@ def add_return_timestamps(self, return_timestamps):
            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 @@ def update(self, value):
     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 4f52b330..3f69a189 100644
--- a/pysd/py_backend/model.py
+++ b/pysd/py_backend/model.py
@@ -10,6 +10,7 @@
 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 @@ def __init__(self, py_model_file, params=None, return_func=None,
         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 @@ def __init__(self, py_model_file, params=None, return_func=None,
             # 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 @@ def modules(self) -> Union[dict, None]:
 
     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 @@ def state(self, new_value):
         [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 @@ def initialize_external_data(self, externals=None):
         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 @@ def initialize_external_data(self, externals=None):
         -------
         None
 
+        See also
+        --------
+        :func:`pysd.py_backend.model.Macro.serialize_externals`
+
         """
 
         if not externals:
@@ -612,6 +576,9 @@ def serialize_externals(self, export_path="externals.nc",
         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 @@ def serialize_externals(self, export_path="externals.nc",
         -------
         None
 
+        See also
+        --------
+        :func:`pysd.py_backend.model.Macro.initialize_external_data`
+
         """
         data = {}
         metadata = {}
@@ -726,7 +697,9 @@ def __include_for_serialization(self, ext, py_name_clean, data, metadata,
         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 @@ def get_args(self, param):
         >>> 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 @@ def get_coords(self, param):
         >>> 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 @@ def __getitem__(self, param):
         ----
         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 @@ def get_series_data(self, param):
                 "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 @@ def set_components(self, params, new=False):
                     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 @@ def set_components(self, params, new=False):
 
                 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 @@ def set_components(self, params, new=False):
                 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 @@ def _timeseries_component(self, series, dims):
 
         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 @@ def _constant_component(self, value, dims):
         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 @@ def _get_elements_to_initialize(self, modified_statefuls):
 
         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 @@ def set_stateful(self, stateful_dict):
 
         """
         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 @@ def run(self, params=None, return_columns=None, return_timestamps=None,
             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 @@ def run(self, params=None, return_columns=None, return_timestamps=None,
             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 @@ def run(self, params=None, return_columns=None, return_timestamps=None,
 
         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 @@ def run(self, params=None, return_columns=None, return_timestamps=None,
         >>> 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 @@ def set_stepper(self, output_obj, params=None, step_vars=[],
                     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 @@ def set_stepper(self, output_obj, params=None, step_vars=[],
            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 @@ def step(self, num_steps=1, step_vars={}):
         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 @@ def step(self, num_steps=1, step_vars={}):
         -------
         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 @@ def _config_simulation(self, params, return_columns, return_timestamps,
         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 @@ def _config_simulation(self, params, return_columns, return_timestamps,
         # 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 @@ def _set_control_vars(self, return_timestamps, final_time, time_step,
             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 @@ def select_submodel(self, vars=[], modules=[], exogenous_components={}):
             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 @@ def select_submodel(self, vars=[], modules=[], exogenous_components={}):
         ...         "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 @@ def select_submodel(self, vars=[], modules=[], exogenous_components={}):
         # 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 @@ def get_dependencies(self, vars=[], modules=[]):
         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 @@ def get_dependencies(self, vars=[], modules=[]):
         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 @@ def get_vars_in_module(self, module):
         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 @@ def get_vars_in_module(self, module):
 
         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 @@ def set_initial_condition(self, initial_condition):
         >>> 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 @@ def set_initial_condition(self, initial_condition):
             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 @@ def _integrate(self):
         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 abfd8451..18ed23b1 100644
--- a/pysd/py_backend/output.py
+++ b/pysd/py_backend/output.py
@@ -10,6 +10,7 @@
 
 from csv import QUOTE_NONE
 from pathlib import Path
+from collections import defaultdict
 
 import regex as re
 
@@ -46,8 +47,7 @@ def __init__(self, out_file=None):
         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 @@ def __init__(self, next):
         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 @@ def initialize(self, model):
         None
 
         """
+        self.__step = 0
         self.ds = self.nc.Dataset(self.out_file, "w")
 
         # defining global attributes
@@ -295,7 +282,7 @@ def initialize(self, model):
         # 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 @@ def update(self, model):
         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 @@ def __init__(self, next):
         super().__init__(next)
         self.ds = None
         self.out_file = None
+        self.__step = 0
 
     def process_output(self, out_file):
         """
@@ -457,7 +446,7 @@ def process_output(self, out_file):
 
     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 @@ def initialize(self, model):
         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 @@ def update(self, model):
         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 @@ def postprocess(self, **kwargs):
             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 @@ def postprocess(self, **kwargs):
 
     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 @@ def add_run_elements(self, model):
         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 b987059f..41e490e9 100644
--- a/pysd/py_backend/statefuls.py
+++ b/pysd/py_backend/statefuls.py
@@ -384,6 +384,9 @@ def initialize(self, init_val=None):
         # 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 fe3f39d2..0cbe8b5b 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 a4afd463..de24aac3 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 @@
     "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 @@
     "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 @@
     },
     "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 @@
     },
     "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 @@
     },
     "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 @@
     },  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 @@
     },
     "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 @@
     },
     "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 @@
     },
     "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 @@
     },
     "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 @@
     },
     "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 @@
     },
     "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 @@
     },
     "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 @@
     },
     "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 @@
     },
     "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 @@
     },
     "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 @@
     },
     "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 @@
     },
     "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 @@
     },
     "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 @@ def data_path(self, test_folder, test_data):
         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 @@ def kwargs(self, test_data):
             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 d452a441..a422bb89 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 @@
         "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 @@ def model_path(self, test_folder, test_data):
         """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 @@ def kwargs(self, test_data):
             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 341524bf..36c4fce8 100644
--- a/tests/pytest_pysd/pytest_output.py
+++ b/tests/pytest_pysd/pytest_output.py
@@ -312,15 +312,6 @@ def test_dataset_handler_step_setter(self, tmp_path, model):
         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 3dfd434c..cb5b85fa 100644
--- a/tests/pytest_pysd/pytest_pysd.py
+++ b/tests/pytest_pysd/pytest_pysd.py
@@ -22,6 +22,8 @@
     + "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_stepper = more_tests.joinpath(
     "stepper_cache/stepper_cache.mdl")
 
+
 class TestPySD():
 
     @pytest.mark.parametrize("model_path", [test_model])
@@ -306,11 +309,14 @@ def test_set_timeseries_parameter(self, model):
             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 @@ def test_set_timeseries_parameter_inline(self, model):
             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 @@ def test_set_components_warnings(self, model):
     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 @@ def test_set_subscripted_timeseries_parameter_with_constant(self, model):
         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 @@ def test_set_subscripted_timeseries_parameter_with_partial_xarray(self,
                                                               ).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 @@ def test_set_subscripted_timeseries_parameter_with_xarray(self, model):
             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 @@ def test_set_initial_value_subscripted_value_with_numpy_error(self, model):
         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 @@ def test_set_initial_condition_origin_full(self, model):
         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 @@ def test_set_initial_condition_origin_short(self, model):
         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 @@ def test_set_initial_condition_for_stock_component(self, model):
         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 @@ def test_get_args(self, model, args):
                     "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 @@ def test__integrate(self, tmp_path, model):
         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 @@ def test_files(self, model, model_path, tmp_path):
         # 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 @@ def test_no_crosstalk(self, _root):
         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 @@ def test_change_constant_pipe(self, model):
 
         # 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 @@ def test_run_export_import(self, tmp_path, model, return_ts, final_t):
         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 @@ def test_stepper_cache_assignment(self, model):
         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 @@ def test_stepper_workflow(self, model):
 
         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 470fa89d..fba77029 100644
--- a/tests/pytest_pysd/pytest_select_submodel.py
+++ b/tests/pytest_pysd/pytest_select_submodel.py
@@ -177,7 +177,10 @@ def test_select_submodel(self, model, variables, modules,
                 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 @@ def test_select_submodel(self, model, variables, modules,
 
         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 7f5e8cc6..9f3c02be 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 @@ def test_get_data_and_run(self, model, expected):
             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 5c4b2ae2..8ec5a184 100644
--- a/tests/pytest_types/external/pytest_external.py
+++ b/tests/pytest_types/external/pytest_external.py
@@ -66,9 +66,8 @@ def test_read_clean_opyxl(self, _root):
         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 bd7e244f..c7060f64 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 @@ def test_delay(self):
         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 87d980f6..ac2edbcd 160000
--- a/tests/test-models
+++ b/tests/test-models
@@ -1,1 +1,1 @@
-Subproject commit 87d980f61c730efffc829e565cf375795f5af36a
+Subproject commit ac2edbcde87aec2dfd8cc338d0574970c549a2d9