Repository: SDXorg/pysd (language: python)

Issue:
Command-line usage
- Command-line usage added. Now pysd can be called from command line. (python -m pysd).
- Test added (they take some time due to subprocess.run speed).
- Documentation updated.
- Correct bug of initial time in outputs when initial time is modified,
- Some other minor fixes

This PR branch [cli] is a sub-branch of the branch from #269 [bench], as I needed to use totally flat dataframes to save the outputs.
@JamesPHoughton if you want you can accept this PR and close #269 
Fixed issues
This pull request fixes issue #244. It also throws an error if the file to be translated is not a Vensim file (neither mdl nor MDL). I added two tests to check this in the file ‘unit_test_pysd.py’: 
· ‘test_read_not_model_vensim’ checks if PySD throws an exception when the file to translate is not a Vensim model. 
· “test_run_upper” checks if, with an uppercase extension file, Vensim generates the expected output. The new model with uppercase extension is the teacup model.
The model tests used are in the new folder named ‘more-tests’ because I cannot add them to the submodule. 

Fix issue #245. Added RANDOM 0 1 function and documentation for RANDOM UNIFORM function. 

Fix issue #246. I updated model_structure_grammar to define a comment of a Vensim sentence as anything that is not a tilde nor a pipe. I added a test in ‘unit_test_pysd’ file to check the ‘teacup_3quotes.mdl’ model. That model is in ‘more-test’ folder and current version of PySD does not support it.


Discussion / hints:

[![Coverage Status](https://coveralls.io/builds/40781542/badge)](https://coveralls.io/builds/40781542)

Coverage decreased (-4.5%) to 91.004% when pulling **4efee97249ca5b8730644252b9ef5009c6b00d43 on enekomartinmartinez:cli** into **85606265aa92878c35a41dd81ce9588d23350e19 on JamesPHoughton:master**.


[![Coverage Status](https://coveralls.io/builds/38646341/badge)](https://coveralls.io/builds/38646341)

Coverage decreased (-0.1%) to 94.033% when pulling **7368e3bc08a6d1f9cad3dd2bc204c73fa1225647 on marrobl:master** into **03a92578337c8834d4ef2b183263c530682c1a1c on JamesPHoughton:master**.

Hey @marrobl, 

Thanks for the PR, and welcome to the community! =)

To update the test suite, would you mind adding the model files for testing an odd number of quotes, and the case sensitivity of the model filename as examples in the [test models](https://github.com/SDXorg/test-models) repo? We're collecting tests there for behaviors that are general to vensim model tools, so that they can be used for any project that implements a model translator. Then references for those can go in [the vensim integration test file](https://github.com/JamesPHoughton/pysd/blob/master/tests/integration_test_vensim_pathway.py).

Also, can you ignore changes in your `.gitmodules` fileso that we still point to the original test models repo?

Thanks!
Hi, @JamesPHoughton! 
Finally, I changed what you asked me and created a pull request in test-models repository. 
I hope everything is correct and if not let me know and i will change it.   
> Hey @marrobl , I'm working out why the tests aren't running properly, but in the meantime, a few questions for you. =)

Sorry @JamesPHoughton, it was my fault. In my PySD repository, the submodule was not pointing to master branch. The submodule pointed to a branch that did not have the new test models.

Fix patch:
diff --git a/docs/advanced_usage.rst b/docs/advanced_usage.rst
index 750501da..e5588717 100644
--- a/docs/advanced_usage.rst
+++ b/docs/advanced_usage.rst
@@ -41,6 +41,32 @@ We can substitute this function directly for the heat_loss_to_room model compone
 
 If you want to replace a subscripted variable, you need to ensure that the output from the new function is the same as the previous one. You can check the current coordinates and dimensions of a component by using :py:data:`.get_coords(variable_name)` as it is explained in :doc:`basic usage <../basic_usage>`.
 
-Supplying additional arguments to the integrator
-------------------------------------------------
-the :py:func:`run` function's argument `intg_kwargs` is a pass-through for keyword arguments to scipy's odeint function, and as such can take on any of the keywords that odeint recognizes.
+Starting simulations from an end-state of another simulation
+------------------------------------------------------------
+The current state of a model can be saved in a pickle file using the :py:data:`.export()`method::
+
+   import pysd
+   model1 = pysd.read_vensim("my_model.mdl")
+   model1.run(final_time=50)
+   model1.export("final_state.pic")
+
+Then the exported data can be used in another session::
+
+   import pysd
+   model2 = pysd.load("my_model.py")
+   model2 = run(initial_condition="final_state.pic", return_timestamps=[55, 60])
+
+the new simulation will have initial time equal to 50 with the saved values from the previous one.
+
+.. note::
+   You can set the exact final time of the simulation using the *final_time* argument.
+   If you want to avoid returning the dataframe of the stocks you can use *return_timestamps=[]*::
+
+     model1.run(final_time=50, return_timestamps=[])
+
+.. note::
+   The changes done with *params* arguments are not ported to the new model (*model2*) object that you initialize with *final_state.pic*. If you want to keep them, you need to call run with the same *params* values as in the original model (*model1*).
+
+.. warning::
+  Exported data is saved and loaded using `pickle <https://docs.python.org/3/library/pickle.html>`_, this data can be not compatible with future versions of
+  *PySD* or *xarray*. In order to prevent data losses save always the source code.
diff --git a/pysd/_version.py b/pysd/_version.py
index f49459c7..14d9d2f5 100644
--- a/pysd/_version.py
+++ b/pysd/_version.py
@@ -1,1 +1,1 @@
-__version__ = "1.6.1"
+__version__ = "1.7.0"
diff --git a/pysd/cli/help.txt b/pysd/cli/help.txt
deleted file mode 100644
index e69de29b..00000000
diff --git a/pysd/cli/main.py b/pysd/cli/main.py
index 9b0304fe..aff09e9d 100644
--- a/pysd/cli/main.py
+++ b/pysd/cli/main.py
@@ -31,10 +31,13 @@ def main(args):
         print("\nFinished!")
         sys.exit()
 
-    initialize(model, options)
+    model.initialize()
 
     output = model.run(**create_configuration(model, options))
 
+    if options.export_file:
+        model.export(options.export_file)
+
     save(output, options)
     print("\nFinished!")
     sys.exit()
@@ -67,28 +70,6 @@ def load(model_file, missing_values):
                          missing_values=missing_values)
 
 
-def initialize(model, options):
-    """
-    Update model control variables and initialize the model.
-
-    Parameters
-    ----------
-    model: pysd.model object
-
-    options: argparse.Namespace
-
-    Returns
-    -------
-    None
-
-    """
-    for func_name in ['initial_time', 'time_step', 'final_time', 'saveper']:
-        if getattr(options, func_name):
-            model.set_components({func_name: getattr(options, func_name)})
-
-    model.initialize()
-
-
 def create_configuration(model, options):
     """
     create configuration dict to pass to the run method.
@@ -107,12 +88,19 @@ def create_configuration(model, options):
     conf_dict = {
         "progress": options.progress,
         "params": options.new_values['param'],
-        "initial_condition": (model.time(), options.new_values['initial']),
+        "initial_condition": (options.initial_time or model.time(),
+                              options.new_values['initial']),
         "return_columns": options.return_columns,
+        "final_time": options.final_time,
+        "time_step": options.time_step,
+        "saveper": options.saveper,
         "flatten_output": True,  # need to return totally flat DF
         "return_timestamps": options.return_timestamps  # given or None
     }
 
+    if options.import_file:
+        conf_dict["initial_condition"] = options.import_file
+
     return conf_dict
 
 
diff --git a/pysd/cli/parser.py b/pysd/cli/parser.py
index fc4e27c0..b460cfec 100644
--- a/pysd/cli/parser.py
+++ b/pysd/cli/parser.py
@@ -175,6 +175,17 @@ def __call__(self, parser, namespace, values, option_string=None):
     help='provide the return columns separated by commas or a .txt file'
          ' where each row is a variable')
 
+parser.add_argument(
+    '-e', '--export', dest='export_file',
+    type=str, metavar='FILE',
+    help='export to a pickle stateful objects states at the end of the '
+         'simulation')
+
+parser.add_argument(
+    '-i', '--import-initial', dest='import_file',
+    type=str, metavar='FILE',
+    help='import stateful objects states from a pickle file,'
+         'if given initial conditions from var:value will be ignored')
 
 ###################
 # Model arguments #
@@ -209,7 +220,7 @@ def __call__(self, parser, namespace, values, option_string=None):
     action='store', type=split_timestamps,
     metavar='\'value1, value2, .., valueN\'',
     help='provide the return time stamps separated by commas, if given '
-         '--final-time and --saveper will be ignored')
+         '--saveper will be ignored')
 
 
 #######################
diff --git a/pysd/py_backend/builder.py b/pysd/py_backend/builder.py
index 506d695e..34113077 100644
--- a/pysd/py_backend/builder.py
+++ b/pysd/py_backend/builder.py
@@ -134,7 +134,7 @@ def build(elements, subscript_dict, namespace, outfile_name):
     text, root = Imports.get_header(os.path.basename(outfile_name))
 
     text += textwrap.dedent('''
-    __pysd_version__ = "%(version)s"
+    __pysd_version__ = '%(version)s'
 
     _subscript_dict = %(subscript_dict)s
 
@@ -753,7 +753,8 @@ def add_n_delay(identifier, delay_input, delay_time, initial_value, order,
     return "%s()" % py_name, new_structure
 
 
-def add_sample_if_true(identifier, condition, actual_value, initial_value):
+def add_sample_if_true(identifier, condition, actual_value, initial_value,
+                       subs):
     """
     Creates code to instantiate a stateful 'SampleIfTrue' object,
     and provides reference to that object's output.
@@ -774,6 +775,10 @@ def add_sample_if_true(identifier, condition, actual_value, initial_value):
     initial_value: <string>
         This is used to initialize the state of the sample if true function.
 
+    subs: list of strings
+        List of strings of subscript indices that correspond to the
+        list of expressions, and collectively define the shape of the output
+
     Returns
     -------
     reference: basestring
@@ -786,25 +791,49 @@ def add_sample_if_true(identifier, condition, actual_value, initial_value):
     """
     Imports.add('functions', 'SampleIfTrue')
 
+    new_structure = []
     py_name = '_sample_if_true_%s' % identifier
 
+    if len(subs) == 0:
+        stateful_py_expr = "SampleIfTrue(lambda: %s, lambda: %s,"\
+                           "lambda: %s, '%s')" % (
+                               condition, actual_value, initial_value, py_name)
+
+    else:
+        stateful_py_expr = "SampleIfTrue(lambda: %s, lambda: %s,"\
+                           "_sample_if_true_init_%s, '%s')" % (
+                               condition, actual_value, identifier, py_name)
+
+        # following elements not specified in the model file, but must exist
+        # create the delay initialization element
+        new_structure.append({
+            'py_name': '_sample_if_true_init_%s' % identifier,
+            'real_name': 'Implicit',
+            'kind': 'setup',  # not specified in the model file, but must exist
+            'py_expr': initial_value,
+            'subs': subs,
+            'doc': 'Provides initial conditions for %s function' % identifier,
+            'unit': 'See docs for %s' % identifier,
+            'lims': 'None',
+            'eqn': 'None',
+            'arguments': ''
+        })
     # describe the stateful object
-    stateful = {
+    new_structure.append({
         'py_name': py_name,
         'real_name': 'Sample if true of %s' % identifier,
         'doc': 'Initial value: %s \n  Input: %s \n Condition: %s' % (
             initial_value, actual_value, condition),
-        'py_expr': "SampleIfTrue(lambda: %s, lambda: %s, lambda: %s, '%s')" % (
-                    condition, actual_value, initial_value, py_name),
+        'py_expr': stateful_py_expr,
         'unit': 'None',
         'lims': 'None',
         'eqn': 'None',
         'subs': '',
         'kind': 'stateful',
         'arguments': ''
-    }
+    })
 
-    return "%s()" % stateful['py_name'], [stateful]
+    return "%s()" % py_name, new_structure
 
 
 def add_n_smooth(identifier, smooth_input, smooth_time, initial_value, order,
diff --git a/pysd/py_backend/functions.py b/pysd/py_backend/functions.py
index 78a89615..27fe8a77 100644
--- a/pysd/py_backend/functions.py
+++ b/pysd/py_backend/functions.py
@@ -8,6 +8,7 @@
 import inspect
 import os
 import re
+import pickle
 import random
 import warnings
 from importlib.machinery import SourceFileLoader
@@ -89,12 +90,20 @@ def __init__(self, ddt, initial_value, py_name="Integ object"):
         self.shape_info = None
         self.py_name = py_name
 
-    def initialize(self):
-        self.state = self.init_func()
+    def initialize(self, init_val=None):
+        if init_val is None:
+            self.state = self.init_func()
+        else:
+            self.state = init_val
         if isinstance(self.state, xr.DataArray):
             self.shape_info = {'dims': self.state.dims,
                                'coords': self.state.coords}
 
+    def export(self):
+        return {self.py_name: {
+            'state': self.state,
+            'shape_info': self.shape_info}}
+
 
 class Delay(DynamicStateful):
     """
@@ -130,7 +139,7 @@ def __init__(self, delay_input, delay_time, initial_value, order,
         self.shape_info = None
         self.py_name = py_name
 
-    def initialize(self):
+    def initialize(self, init_val=None):
         order = self.order_func()
 
         if order != int(order):
@@ -146,7 +155,10 @@ def initialize(self):
                           'Delay time very small, casting delay order '
                           + f'from {int(order)} to {self.order}')
 
-        init_state_value = self.init_func() * self.delay_time_func()
+        if init_val is None:
+            init_state_value = self.init_func() * self.delay_time_func()
+        else:
+            init_state_value = init_val * self.delay_time_func()
 
         if isinstance(init_state_value, xr.DataArray):
             # broadcast self.state
@@ -173,6 +185,11 @@ def ddt(self):
             inflows[0] = self.input_func()
         return (inflows - outflows) * self.order
 
+    def export(self):
+        return {self.py_name: {
+            'state': self.state,
+            'shape_info': self.shape_info}}
+
 
 class DelayN(DynamicStateful):
     """
@@ -209,7 +226,7 @@ def __init__(self, delay_input, delay_time, initial_value, order,
         self.shape_info = None
         self.py_name = py_name
 
-    def initialize(self):
+    def initialize(self, init_val=None):
         order = self.order_func()
 
         if order != int(order):
@@ -225,7 +242,10 @@ def initialize(self):
                           'Delay time very small, casting delay order '
                           + f'from {int(order)} to {self.order}')
 
-        init_state_value = self.init_func() * self.delay_time_func()
+        if init_val is None:
+            init_state_value = self.init_func() * self.delay_time_func()
+        else:
+            init_state_value = init_val * self.delay_time_func()
 
         if isinstance(init_state_value, xr.DataArray):
             # broadcast self.state
@@ -263,6 +283,12 @@ def ddt(self):
         inflows[0] = self.input_func()
         return (inflows - outflows)*self.order
 
+    def export(self):
+        return {self.py_name: {
+            'state': self.state,
+            'times': self.times,
+            'shape_info': self.shape_info}}
+
 
 class DelayFixed(DynamicStateful):
     """
@@ -291,7 +317,7 @@ def __init__(self, delay_input, delay_time, initial_value, tstep,
         self.pointer = 0
         self.py_name = py_name
 
-    def initialize(self):
+    def initialize(self, init_val=None):
         order = max(self.delay_time_func()/self.tstep(), 1)
 
         if order != int(order):
@@ -303,7 +329,10 @@ def initialize(self):
         # need to add a small decimal to ensure that 0.5 is rounded to 1
         self.order = round(order + small_vensim)  # The order can only be set once
 
-        init_state_value = self.init_func()
+        if init_val is None:
+            init_state_value = self.init_func()
+        else:
+            init_state_value = init_val
 
         self.state = init_state_value
         self.pipe = [init_state_value] * self.order
@@ -319,6 +348,12 @@ def update(self, state):
         self.pointer = (self.pointer + 1) % self.order
         self.state = self.pipe[self.pointer]
 
+    def export(self):
+        return {self.py_name: {
+            'state': self.state,
+            'pointer': self.pointer,
+            'pipe': self.pipe}}
+
 
 class Smooth(DynamicStateful):
     """
@@ -346,9 +381,14 @@ def __init__(self, smooth_input, smooth_time, initial_value, order,
         self.shape_info = None
         self.py_name = py_name
 
-    def initialize(self):
+    def initialize(self, init_val=None):
         self.order = self.order_func()  # The order can only be set once
-        init_state_value = self.init_func()
+
+        if init_val is None:
+            init_state_value = self.init_func()
+        else:
+            init_state_value = init_val
+
         if isinstance(init_state_value, xr.DataArray):
             # broadcast self.state
             self.state = init_state_value.expand_dims({
@@ -356,7 +396,7 @@ def initialize(self):
             self.shape_info = {'dims': self.state.dims,
                                'coords': self.state.coords}
         else:
-            self.state = np.array([self.init_func()] * self.order)
+            self.state = np.array([init_state_value] * self.order)
 
     def __call__(self):
         if self.shape_info:
@@ -372,6 +412,11 @@ def ddt(self):
             targets[0] = self.input_func()
         return (targets - self.state) * self.order / self.smooth_time_func()
 
+    def export(self):
+        return {self.py_name: {
+            'state': self.state,
+            'shape_info': self.shape_info}}
+
 
 class Trend(DynamicStateful):
     """
@@ -396,9 +441,13 @@ def __init__(self, trend_input, average_time, initial_trend,
         self.input_func = trend_input
         self.py_name = py_name
 
-    def initialize(self):
-        self.state = self.input_func()\
-            / (1 + self.init_func() * self.average_time_function())
+    def initialize(self, init_val=None):
+        if init_val is None:
+            self.state = self.input_func()\
+                / (1 + self.init_func()*self.average_time_function())
+        else:
+            self.state = self.input_func()\
+                / (1 + init_val*self.average_time_function())
 
         if isinstance(self.state, xr.DataArray):
             self.shape_info = {'dims': self.state.dims,
@@ -411,8 +460,13 @@ def __call__(self):
     def ddt(self):
         return (self.input_func() - self.state) / self.average_time_function()
 
+    def export(self):
+        return {self.py_name: {
+            'state': self.state,
+            'shape_info': self.shape_info}}
 
-class SampleIfTrue(Stateful):
+
+class SampleIfTrue(DynamicStateful):
     def __init__(self, condition, actual_value, initial_value,
                  py_name="SampleIfTrue object"):
         """
@@ -428,38 +482,63 @@ def __init__(self, condition, actual_value, initial_value,
         super().__init__()
         self.condition = condition
         self.actual_value = actual_value
-        self.initial_value = initial_value
+        self.init_func = initial_value
         self.py_name = py_name
 
-    def initialize(self):
-        self.state = self.initial_value()
+    def initialize(self, init_val=None):
+        if init_val is None:
+            self.state = self.init_func()
+        else:
+            self.state = init_val
+        if isinstance(self.state, xr.DataArray):
+            self.shape_info = {'dims': self.state.dims,
+                               'coords': self.state.coords}
 
     def __call__(self):
-        self.state = if_then_else(self.condition(),
-                                  self.actual_value,
-                                  lambda: self.state)
-        return self.state
+        return if_then_else(self.condition(),
+                            self.actual_value,
+                            lambda: self.state)
+
+    def ddt(self):
+        return np.nan
+
+    def update(self, state):
+        self.state = self.state*0 + if_then_else(self.condition(),
+                                                 self.actual_value,
+                                                 lambda: self.state)
+
+    def export(self):
+        return {self.py_name: {
+            'state': self.state,
+            'shape_info': self.shape_info}}
 
 
 class Initial(Stateful):
     """
     Implements INITIAL function
     """
-    def __init__(self, func, py_name="Initial object"):
+    def __init__(self, initial_value, py_name="Initial object"):
         """
 
         Parameters
         ----------
-        func: function
+        initial_value: function
         py_name: str
           Python name to identify the object
         """
         super().__init__()
-        self.func = func
+        self.init_func = initial_value
         self.py_name = py_name
 
-    def initialize(self):
-        self.state = self.func()
+    def initialize(self, init_val=None):
+        if init_val is None:
+            self.state = self.init_func()
+        else:
+            self.state = init_val
+
+    def export(self):
+        return {self.py_name: {
+            'state': self.state}}
 
 
 class Macro(DynamicStateful):
@@ -493,6 +572,7 @@ def __init__(self, py_model_file, params=None, return_func=None,
         self.time = time
         self.time_initialization = time_initialization
         self.py_name = py_name
+        self.initialize_order = None
 
         # need a unique identifier for the imported module.
         module_name = os.path.splitext(py_model_file)[0]\
@@ -572,10 +652,7 @@ def initialize(self, initialization_order=None):
 
         # Initialize time
         if self.time is None:
-            if self.time_initialization is None:
-                self.time = Time()
-            else:
-                self.time = self.time_initialization()
+            self.time = self.time_initialization()
 
         self.components.cache.clean()
         self.components.cache.time = self.time()
@@ -591,6 +668,33 @@ def initialize(self, initialization_order=None):
 
         Excels.clean()
 
+        remaining = set(self._stateful_elements)
+        if len(set([element.py_name for element in self._stateful_elements]))\
+           == len(set(self._stateful_elements)) and self.initialize_order:
+            # use elements names to initialize them, this is available
+            # after the model is initialized one time
+            # solves issue #247 until we have a dependency dictionary
+            try:
+                for element_name in self.initialize_order:
+                    for element in remaining:
+                        if element.py_name == element_name:
+                            element.initialize()
+                            break
+                    remaining.remove(element)
+                assert len(remaining) == 0
+                return
+            except Exception as err:
+                # if user includes new stateful objects or some other
+                # dependencies the previous initialization order may
+                # not be keept
+                warnings.warn(
+                    err.args[0] +
+                    "\n\nNot able to initialize statefull elements "
+                    "with the same order as before..."
+                    "Trying to find a new order.")
+                # initialize as always
+
+        self.initialize_order = []
         # Initialize stateful elements
         remaining = set(self._stateful_elements)
         while remaining:
@@ -599,6 +703,7 @@ def initialize(self, initialization_order=None):
                 try:
                     element.initialize()
                     progress.add(element)
+                    self.initialize_order.append(element.py_name)
                 except (KeyError, TypeError, AttributeError):
                     pass
 
@@ -625,6 +730,59 @@ 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
+          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 = {}
+        [stateful_elements.update(component.export()) for component
+         in self._stateful_elements]
+
+        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
+          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']:
+            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.update(time)
+        self.components.cache.reset(time)
+
     def get_args(self, param):
         """
         Returns the arguments of a model element.
@@ -648,7 +806,12 @@ def get_args(self, param):
                 param,
                 self.components._namespace) or param
 
-            func = getattr(self.components, func_name)
+            if hasattr(self.components, func_name):
+                func = getattr(self.components, func_name)
+            else:
+                NameError(
+                    "\n'%s' is not recognized as a model component."
+                    % param)
         else:
             func = param
 
@@ -657,7 +820,10 @@ def get_args(self, param):
             return func.args
         else:
             # regular functions
-            return inspect.getfullargspec(func)[0]
+            args = inspect.getfullargspec(func)[0]
+            if 'self' in args:
+                args.remove('self')
+            return args
 
     def get_coords(self, param):
         """
@@ -683,7 +849,12 @@ def get_coords(self, param):
                 param,
                 self.components._namespace) or param
 
-            func = getattr(self.components, func_name)
+            if hasattr(self.components, func_name):
+                func = getattr(self.components, func_name)
+            else:
+                NameError(
+                    "\n'%s' is not recognized as a model component."
+                    % param)
         else:
             func = param
 
@@ -717,10 +888,9 @@ def set_components(self, params):
 
 
         """
-        # It might make sense to allow the params argument to take a pandas series, where
-        # the indices of the series are variable names. This would make it easier to
-        # do a Pandas apply on a DataFrame of parameter values. However, this may conflict
-        # with a pandas series being passed in as a dictionary element.
+        # 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.
 
         for key, value in params.items():
             func_name = utils.get_value_by_insensitive_key_or_value(
@@ -728,15 +898,16 @@ def set_components(self, params):
                 self.components._namespace)
 
             if isinstance(value, np.ndarray) or isinstance(value, list):
-                raise ValueError('When setting ' + key +'\n'
-                                 + 'Setting subscripted must be done'
-                                 + 'using a xarray.DataArray with the '
-                                 + 'correct dimensions or a constant value '
-                                 + '(https://pysd.readthedocs.io/en/master/basic_usage.html)')
+                raise TypeError(
+                    'When setting ' + key + '\n'
+                    'Setting subscripted must be done using a xarray.DataArray'
+                    ' with the correct dimensions or a constant value '
+                    '(https://pysd.readthedocs.io/en/master/basic_usage.html)')
 
             if func_name is None:
-                raise NameError('%s is not recognized as a model component'
-                                % key)
+                raise NameError(
+                    "\n'%s' is not recognized as a model component."
+                    % key)
 
             try:
                 func = getattr(self.components, func_name)
@@ -837,73 +1008,104 @@ def _constant_component(self, value, dims, args=[]):
         else:
             return lambda: value
 
-    def set_state(self, t, state):
-        """ Set the system state.
+    def set_state(self, t, initial_value):
+        """ Old set_state method use set_initial_value"""
+        warnings.warn(
+            "\nset_state will be deprecated, use set_initial_value instead.",
+            FutureWarning)
+        self.set_initial_value(t, initial_value)
+
+    def set_initial_value(self, t, initial_value):
+        """ Set the system initial value.
 
         Parameters
         ----------
         t : numeric
             The system time
 
-        state : dict
-            A (possibly partial) dictionary of the system state.
+        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
+
         """
         self.time.update(t)
         self.components.cache.reset(t)
+        stateful_name = "_NONE"
+        # TODO make this more solid, link with builder or next TODO?
+        stateful_init = [
+            "_integ_", "_delay_", "_delayfixed_", "_delayn_",
+            "_sample_if_true_", "_smooth_", "_trend_", "_initial_"]
 
-        for key, value in state.items():
-            # TODO Implement map with reference between component and stateful element?
+        for key, value in initial_value.items():
             component_name = utils.get_value_by_insensitive_key_or_value(
                 key, self.components._namespace)
             if component_name is not None:
-                stateful_name = '_integ_%s' % component_name
+                for element in self._stateful_elements:
+                    # TODO make this more solid, add link between stateful
+                    # objects and model vars
+                    for init in stateful_init:
+                        if init + component_name == element.py_name:
+                            stateful_name = element.py_name
             else:
                 component_name = key
                 stateful_name = key
 
             try:
-                if component_name.startswith('_integ_'):
-                    # we need to check the original expression to retrieve
-                    # the dimensions
-                    _, dims = self.get_coords(component_name[7:])
-                else:
-                    _, dims = self.get_coords(component_name)
-            except (AttributeError, TypeError):
+                _, dims = self.get_coords(component_name)
+            except TypeError:
                 dims = None
 
+            if isinstance(value, xr.DataArray)\
+               and not set(value.dims).issubset(set(dims)):
+                raise ValueError(
+                    f"\nInvalid dimensions for {component_name}."
+                    f"It should be a subset of {dims}, "
+                    f"but passed value has {list(value.dims)}")
+
             if isinstance(value, np.ndarray) or isinstance(value, list):
-                raise ValueError('When setting ' + key + '\n'
-                                 + 'Setting subscripted must be done'
-                                 + 'using a xarray.DataArray with the '
-                                 + 'correct dimensions or a constant value '
-                                 + '(https://pysd.readthedocs.io/en/master/basic_usage.html)')
+                raise TypeError(
+                    'When setting ' + key + '\n'
+                    'Setting subscripted must be done using a xarray.DataArray'
+                    ' with the correct dimensions or a constant value '
+                    '(https://pysd.readthedocs.io/en/master/basic_usage.html)')
 
             # Try to update stateful component
             if hasattr(self.components, stateful_name):
-                try:
-                    element = getattr(self.components, stateful_name)
-                    if dims:
-                        value = utils.rearrange(
-                            value, dims,
-                            self.components._subscript_dict)
-                    element.update(value)
-                    self.components.cache.clean()
-                except AttributeError:
-                    print("'%s' has no state elements, assignment failed")
-                    raise
+                element = getattr(self.components, stateful_name)
+                if dims:
+                    value = utils.rearrange(
+                        value, dims,
+                        self.components._subscript_dict)
+                element.initialize(value)
+                self.components.cache.clean()
             else:
                 # Try to override component
-                try:
-                    setattr(self.components, component_name,
-                            self._constant_component(
-                                value, dims,
-                                self.get_args(component_name)))
-                    self.components.cache.clean()
-                except AttributeError:
-                    print("'%s' has no component, assignment failed")
-                    raise
+                warnings.warn(
+                    f"\nSetting {component_name} to a constant value with "
+                    "initial_conditions will be deprecated. Use params={"
+                    f"'{component_name}': {value}"+"} instead.",
+                    FutureWarning)
+
+                setattr(self.components, component_name,
+                        self._constant_component(
+                            value, dims,
+                            self.get_args(component_name)))
+                self.components.cache.clean()
+
+    def set_stateful(self, stateful_dict):
+        """
+        Set stateful values.
+
+        Parameters
+        ----------
+        stateful_dict: dict
+          Dictionary of the stateful elements and the attributes to change.
+
+        """
+        for element, attrs in stateful_dict.items():
+            for attr, value in attrs.items():
+                setattr(getattr(self.components, element), attr, value)
 
     def doc(self):
         """
@@ -1004,66 +1206,92 @@ def initialize(self):
         External.missing = self.missing_values
         super().initialize()
 
-    def _build_euler_timeseries(self, return_timestamps=None):
+    def _build_euler_timeseries(self, return_timestamps=None, final_time=None):
         """
         - The integration steps need to include the return values.
         - There is no point running the model past the last return value.
         - The last timestep will be the last in that requested for return
-        - Spacing should be at maximum what is specified by the integration time step.
-        - The initial time should be the one specified by the model file, OR
-          it should be the initial condition.
-        - This function needs to be called AFTER the model is set in its initial state
+        - Spacing should be at maximum what is specified by the integration
+          time step.
+        - The initial time should be the one specified by the model file,
+          OR it should be the initial condition.
+        - This function needs to be called AFTER the model is set in its
+          initial state
         Parameters
         ----------
         return_timestamps: numpy array
-          Must be specified by user or built from model file before this function is called.
+            Must be specified by user or built from model file before this
+            function is called.
+
+        final_time: float or None
+            Final time of the simulation. If float, the given final time
+            will be used. If None, the last return_timestamps will be used.
+            Default is None.
 
         Returns
         -------
         ts: numpy array
             The times that the integrator will use to compute time history
+
         """
         t_0 = self.time()
-        t_f = return_timestamps[-1]
-        dt = self.components.time_step()
-        ts = np.arange(t_0, t_f, dt, dtype=np.float64)
+        try:
+            t_f = return_timestamps[-1]
+        except IndexError:
+            # return_timestamps is an empty list
+            # model default final time or passed argument value
+            t_f = self.final_time
+
+        if final_time is not None:
+            t_f = max(final_time, t_f)
+
+        ts = np.arange(
+            t_0, t_f+self.time_step/2, self.time_step, dtype=np.float64)
 
-        # Add the returned time series into the integration array. Best we can do for now.
-        # This does change the integration ever so slightly, but for well-specified
-        # models there shouldn't be sensitivity to a finer integration time step.
-        ts = np.sort(np.unique(np.append(ts, return_timestamps)))
-        return ts
+        # Add the returned time series into the integration array.
+        # Best we can do for now. This does change the integration ever
+        # so slightly, but for well-specified models there shouldn't be
+        # sensitivity to a finer integration time step.
+        return np.sort(np.unique(np.append(ts, return_timestamps)))
 
     def _format_return_timestamps(self, return_timestamps=None):
         """
         Format the passed in return timestamps value as a numpy array.
         If no value is passed, build up array of timestamps based upon
         model start and end times, and the 'saveper' value.
+
+        Parameters
+        ----------
+        return_timestamps: float, iterable of floats or None (optional)
+          Iterable of timestamps to return or None. Default is None.
+
+        Returns
+        -------
+        ndarray (float)
+
         """
         if return_timestamps is None:
             # Build based upon model file Start, Stop times and Saveper
-            # Vensim's standard is to expect that the data set includes the `final time`,
-            # so we have to add an extra period to make sure we get that value in what
-            # numpy's `arange` gives us.
-            return_timestamps_array = np.arange(
+            # Vensim's standard is to expect that the data set includes
+            # the `final time`, so we have to add an extra period to
+            # make sure we get that value in what numpy's `arange` gives us.
+
+            return np.arange(
                 self.time(),
-                self.components.final_time() + self.components.saveper()/2,
-                self.components.saveper(), dtype=float
+                self.final_time + self.saveper/2,
+                self.saveper, dtype=float
             )
-        elif inspect.isclass(range) and isinstance(return_timestamps, range):
-            return_timestamps_array = np.array(return_timestamps, ndmin=1)
-        elif isinstance(return_timestamps, (list, int, float, np.ndarray)):
-            return_timestamps_array = np.array(return_timestamps, ndmin=1)
-        elif isinstance(return_timestamps, pd.Series):
-            return_timestamps_array = return_timestamps.as_matrix()
-        else:
-            raise TypeError('`return_timestamps` expects a list, array,'
-                            'pandas Series, or numeric value')
-        return return_timestamps_array
+
+        try:
+            return np.array(return_timestamps, ndmin=1, dtype=float)
+        except Exception:
+            raise TypeError(
+                '`return_timestamps` expects an iterable of numeric values'
+                ' or a single numeric value')
 
     def run(self, params=None, return_columns=None, return_timestamps=None,
-            initial_condition='original', reload=False, progress=False,
-            flatten_output=False):
+            initial_condition='original', final_time=None, time_step=None,
+            saveper=None, reload=False, progress=False, flatten_output=False):
         """
         Simulate the model's behavior over time.
         Return a pandas dataframe with timestamps as rows,
@@ -1071,30 +1299,48 @@ def run(self, params=None, return_columns=None, return_timestamps=None,
 
         Parameters
         ----------
-        params : dictionary (optional)
+        params: dict (optional)
             Keys are strings of model component names.
             Values are numeric or pandas Series.
             Numeric values represent constants over the model integration.
             Timeseries will be interpolated to give time-varying input.
 
-        return_timestamps : list, numeric, ndarray (1D) (optional)
+        return_timestamps: list, numeric, ndarray (1D) (optional)
             Timestamps in model execution at which to return state information.
             Defaults to model-file specified timesteps.
 
-        return_columns : list, 'step' or None (optional)
+        return_columns: list, 'step' or None (optional)
             List of string model component names, returned dataframe
             will have corresponding columns. If 'step' only variables with
             cache step will be returned. If None, variables with cache step
             and run will be returned. Default is None.
 
-        initial_condition : 'original'/'o', 'current'/'c' or (t, {state}) (optional)
-            The starting time, and the state of the system (the values of all the stocks)
-            at that starting time.
-
-            * 'original' (default) uses model-file specified initial condition
-            * 'current' uses the state of the model after the previous execution
-            * (t, {state}) lets the user specify a starting time and (possibly partial)
-              list of stock values.
+        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
+            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'.
+
+        final_time: float or None
+            Final time of the simulation. If float, the given value will be
+            used to compute the return_timestamps (if not given) and as a
+            final time. If None the last value of return_timestamps will be
+            used as a final time. Default is None.
+
+        time_step: float or None
+            Time step of the simulation. If float, the given value will be
+            used to compute the return_timestamps (if not given) and
+            euler time series. If None the default value from components
+            will be used. Default is None.
+
+        saveper: float or None
+            Saving step of the simulation. If float, the given value will be
+            used to compute the return_timestamps (if not given). If None
+            the default value from components will be used. Default is None.
 
         reload : bool (optional)
             If True, reloads the model from the translated model file
@@ -1131,14 +1377,21 @@ def run(self, params=None, return_columns=None, return_timestamps=None,
 
         if params:
             self.set_components(params)
+
         self.set_initial_condition(initial_condition)
 
-        # save initial time for the output
+        # save control variables
         self.initial_time = self.time()
+        self.final_time = final_time or self.components.final_time()
+        self.time_step = time_step or self.components.time_step()
+        self.saveper = saveper or max(self.time_step,
+                                      self.components.saveper())
+        # need to take bigger saveper if time_step is > saveper
 
         return_timestamps = self._format_return_timestamps(return_timestamps)
 
-        t_series = self._build_euler_timeseries(return_timestamps)
+        t_series = self._build_euler_timeseries(return_timestamps, final_time)
+        self.final_time = t_series[-1]
 
         if return_columns is None or isinstance(return_columns, str):
             return_columns = self._default_return_columns(return_columns)
@@ -1240,37 +1493,37 @@ def set_initial_condition(self, initial_condition):
 
         Parameters
         ----------
-        initial_condition : <string> or <tuple>
-            Takes on one of the following sets of values:
-
-            * 'original'/'o' : Reset to the model-file specified initial condition.
-            * 'current'/'c' : Use the current state of the system to start
-              the next simulation. This includes the simulation time, so this
-              initial condition must be paired with new return timestamps
-            * (t, {state}) : Lets the user specify a starting time and list of stock values.
+        initial_condition : str or (float, dict)
+            The starting time, and the state of the system (the values of
+            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.
 
         >>> model.set_initial_condition('original')
         >>> model.set_initial_condition('current')
+        >>> model.set_initial_condition('exported_pickle.pic')
         >>> model.set_initial_condition((10, {'teacup_temperature': 50}))
 
         See Also
         --------
-        PySD.set_state()
+        PySD.set_initial_value()
 
         """
 
         if isinstance(initial_condition, tuple):
-            # Todo: check the values more than just seeing if they are a tuple.
-            self.set_state(*initial_condition)
+            self.initialize()
+            self.set_initial_value(*initial_condition)
         elif isinstance(initial_condition, str):
             if initial_condition.lower() in ['original', 'o']:
                 self.initialize()
             elif initial_condition.lower() in ['current', 'c']:
                 pass
             else:
-                raise ValueError('Valid initial condition strings include:  \n' +
-                                 '    "original"/"o",                       \n' +
-                                 '    "current"/"c"')
+                self.import_pickle(initial_condition)
         else:
             raise TypeError('Check documentation for valid entries')
 
@@ -1356,10 +1609,15 @@ def _add_run_elements(self, df, capture_elements):
             df[element] = [getattr(self.components, element)()] * nt
 
         # update initial time values in df (necessary if initial_conditions)
-        for it in ['INITIAL TIME', 'INITIAL_TIME',
-                   'initial time', 'initial_time']:
+        for it in ['initial_time', 'final_time', 'saveper', 'time_step']:
             if it in df:
-                df[it] = self.initial_time
+                df[it] = getattr(self, it)
+            elif it.upper() in df:
+                df[it.upper()] = getattr(self, it)
+            elif it.replace('_', ' ') in df:
+                df[it.replace('_', ' ')] = getattr(self, it)
+            elif it.replace('_', ' ').upper() in df:
+                df[it.replace('_', ' ').upper()] = getattr(self, it)
 
 
 def ramp(time, slope, start, finish=0):
@@ -1664,6 +1922,7 @@ def active_initial(time, expr, init_val):
     else:
         return expr()
 
+
 def bounded_normal(minimum, maximum, mean, std, seed):
     """
     Implements vensim's BOUNDED NORMAL function
@@ -1699,10 +1958,14 @@ def random_uniform(m, x, s):
 
     Returns
     -------
-    A random number from the uniform distribution between m and x (exclusive of the endpoints).
+    A random number from the uniform distribution between m and x
+    (exclusive of the endpoints).
+
     """
-    if(s!=0):
-        warnings.warn("Random uniform with a nonzero seed value, may not give the same result as vensim", RuntimeWarning)
+    if s != 0:
+        warnings.warn(
+            "Random uniform with a nonzero seed value, may not give the "
+            "same result as vensim", RuntimeWarning)
 
     return np.random.uniform(m, x)
 
@@ -1835,4 +2098,3 @@ def vmax(x, dim=None):
         return float(x.max())
 
     return x.max(dim=dim)
-
diff --git a/pysd/py_backend/vensim/vensim2py.py b/pysd/py_backend/vensim/vensim2py.py
index 7fc8e96b..95766a8b 100644
--- a/pysd/py_backend/vensim/vensim2py.py
+++ b/pysd/py_backend/vensim/vensim2py.py
@@ -701,7 +701,8 @@ def parse_units(units_str):
         identifier=element['py_name'],
         condition=args[0],
         actual_value=args[1],
-        initial_value=args[2]
+        initial_value=args[2],
+        subs=element['subs']
     ),
 
     "smooth": lambda element, subscript_dict, args: builder.add_n_smooth(


Test patch:
diff --git a/tests/more-tests/initialization_order/test_initialization_order.py b/tests/more-tests/initialization_order/test_initialization_order.py
new file mode 100644
index 00000000..fb65a776
--- /dev/null
+++ b/tests/more-tests/initialization_order/test_initialization_order.py
@@ -0,0 +1,146 @@
+"""
+Python model 'trying_to_reproduce_bug.py'
+Translated using PySD
+"""
+
+
+from pysd.py_backend.functions import Integ
+from pysd import cache
+
+__pysd_version__ = "1.6.0"
+
+_subscript_dict = {}
+
+_namespace = {
+    "TIME": "time",
+    "Time": "time",
+    "Stock B": "stock_b",
+    "Stock A": "stock_a",
+    "Initial Parameter": "initial_parameter",
+    "FINAL TIME": "final_time",
+    "INITIAL TIME": "initial_time",
+    "SAVEPER": "saveper",
+    "TIME STEP": "time_step",
+}
+
+__data = {"scope": None, "time": lambda: 0}
+
+
+def _init_outer_references(data):
+    for key in data:
+        __data[key] = data[key]
+
+
+def time():
+    return __data["time"]()
+
+
+@cache.step
+def stock_b():
+    """
+    Real Name: Stock B
+    Original Eqn: INTEG(1, Stock A)
+    Units:
+    Limits: (None, None)
+    Type: component
+    Subs: None
+
+
+    """
+    return _integ_stock_b()
+
+
+@cache.step
+def stock_a():
+    """
+    Real Name: Stock A
+    Original Eqn: INTEG (1, Initial Parameter)
+    Units:
+    Limits: (None, None)
+    Type: component
+    Subs: None
+
+
+    """
+    return _integ_stock_a()
+
+
+@cache.run
+def initial_parameter():
+    """
+    Real Name: Initial Parameter
+    Original Eqn: 42
+    Units:
+    Limits: (None, None)
+    Type: constant
+    Subs: None
+
+
+    """
+    return 42
+
+
+@cache.run
+def final_time():
+    """
+    Real Name: FINAL TIME
+    Original Eqn: 1
+    Units: Month
+    Limits: (None, None)
+    Type: constant
+    Subs: None
+
+    The final time for the simulation.
+    """
+    return 1
+
+
+@cache.run
+def initial_time():
+    """
+    Real Name: INITIAL TIME
+    Original Eqn: 0
+    Units: Month
+    Limits: (None, None)
+    Type: constant
+    Subs: None
+
+    The initial time for the simulation.
+    """
+    return 0
+
+
+@cache.step
+def saveper():
+    """
+    Real Name: SAVEPER
+    Original Eqn: TIME STEP
+    Units: Month
+    Limits: (0.0, None)
+    Type: component
+    Subs: None
+
+    The frequency with which output is stored.
+    """
+    return time_step()
+
+
+@cache.run
+def time_step():
+    """
+    Real Name: TIME STEP
+    Original Eqn: 1
+    Units: Month
+    Limits: (0.0, None)
+    Type: constant
+    Subs: None
+
+    The time step for the simulation.
+    """
+    return 1
+
+
+_integ_stock_b = Integ(lambda: 1, lambda: stock_a(), "_integ_stock_b")
+
+
+_integ_stock_a = Integ(lambda: 1, lambda: initial_parameter(), "_integ_stock_a")
diff --git a/tests/unit_test_cli.py b/tests/unit_test_cli.py
index 627167c7..6af5dbcf 100644
--- a/tests/unit_test_cli.py
+++ b/tests/unit_test_cli.py
@@ -289,6 +289,8 @@ def test_model_arguments(self):
         self.assertEqual(out.returncode, 0)
         stocks = load_outputs(out_tab_file)
         self.assertTrue((np.diff(stocks.index.values) == time_step).all())
+        self.assertTrue((stocks['SAVEPER'] == time_step).all().all())
+        self.assertTrue((stocks['TIME STEP'] == time_step).all().all())
         os.remove(out_tab_file)
 
         # check saveper
@@ -300,6 +302,8 @@ def test_model_arguments(self):
         self.assertEqual(out.returncode, 0)
         stocks = load_outputs(out_tab_file)
         self.assertTrue((np.diff(stocks.index.values) == saveper).all())
+        self.assertTrue((stocks['SAVEPER'] == saveper).all().all())
+        self.assertTrue((stocks['TIME STEP'] == time_step).all().all())
         os.remove(out_tab_file)
 
         # check all
@@ -385,3 +389,43 @@ def test_set_timeseries_parameter_lookup(self):
         self.assertTrue(
             (stocks['lookup_2d_time[Row2]'] == temp_timeseries).all())
         os.remove(out_tab_file)
+
+    def test_export_import(self):
+        import pysd
+        from pysd.tools.benchmarking import assert_frames_close
+
+        exp_file = 'teacup15.pic'
+        model = pysd.read_vensim(test_model)
+        stocks = model.run(return_timestamps=[0, 10, 20, 30])
+        self.assertTrue((stocks['INITIAL TIME'] == 0).all().all())
+
+        command = f'{call} -o {out_tab_file} -e {exp_file} -F 15 -R 0,10'\
+                  f' {test_model}'
+
+        command2 = f'{call} -o {out_tab_file} -i {exp_file} -R 20,30'\
+                   f' {test_model}'
+
+        out = subprocess.run(split_bash(command), capture_output=True)
+        self.assertEqual(out.returncode, 0)
+        stocks1 = load_outputs(out_tab_file)
+
+        out = subprocess.run(split_bash(command2), capture_output=True)
+        self.assertEqual(out.returncode, 0)
+        stocks2 = load_outputs(out_tab_file)
+
+        os.remove(exp_file)
+        os.remove(out_tab_file)
+
+        print(stocks1['FINAL TIME'])
+        self.assertTrue((stocks1['INITIAL TIME'] == 0).all().all())
+        self.assertTrue((stocks1['FINAL TIME'] == 15).all().all())
+        self.assertTrue((stocks2['INITIAL TIME'] == 15).all().all())
+        stocks.drop('INITIAL TIME', axis=1, inplace=True)
+        stocks1.drop('INITIAL TIME', axis=1, inplace=True)
+        stocks2.drop('INITIAL TIME', axis=1, inplace=True)
+        stocks.drop('FINAL TIME', axis=1, inplace=True)
+        stocks1.drop('FINAL TIME', axis=1, inplace=True)
+        stocks2.drop('FINAL TIME', axis=1, inplace=True)
+
+        assert_frames_close(stocks1, stocks.loc[[0, 10]])
+        assert_frames_close(stocks2, stocks.loc[[20, 30]])
diff --git a/tests/unit_test_functions.py b/tests/unit_test_functions.py
index 085b2c65..c5cf2e97 100644
--- a/tests/unit_test_functions.py
+++ b/tests/unit_test_functions.py
@@ -337,6 +337,9 @@ def init_func():
         stock.initialize()
         self.assertEqual(stock(), 11)
 
+        stock.initialize(30)
+        self.assertEqual(stock(), 30)
+
     def test_stateful_identification(self):
         import pysd
 
@@ -354,7 +357,7 @@ def test_delay(self):
                                        initial_value=lambda: 4.234,
                                        order=lambda: 3,
                                        tstep=lambda: 0.5,
-                                       py_name='delay')
+                                       py_name='delay_a')
 
         delay_a.initialize()
 
@@ -362,6 +365,40 @@ def test_delay(self):
 
         self.assertEqual(float(delay_a.ddt()[0]), (5-4.234)*3)
 
+        delay_a.initialize(6)
+        self.assertEqual(delay_a(), 6)
+        self.assertEqual(float(delay_a.ddt()[0]), (5-6)*3)
+
+        delay_b = pysd.functions.DelayN(delay_input=lambda: 5,
+                                        delay_time=lambda: 3,
+                                        initial_value=lambda: 4.234,
+                                        order=lambda: 3,
+                                        tstep=lambda: 0.5,
+                                        py_name='delay_b')
+
+        delay_b.initialize()
+
+        self.assertEqual(delay_b(), 4.234)
+
+        self.assertEqual(float(delay_b.ddt()[0]), (5-4.234)*3)
+
+        delay_b.initialize(6)
+        self.assertEqual(delay_b(), 6)
+        self.assertEqual(float(delay_b.ddt()[0]), (5-6)*3)
+
+        delay_c = pysd.functions.DelayFixed(delay_input=lambda: 5,
+                                            delay_time=lambda: 3,
+                                            initial_value=lambda: 4.234,
+                                            tstep=lambda: 0.5,
+                                            py_name='delay_c')
+
+        delay_c.initialize()
+
+        self.assertEqual(delay_c(), 4.234)
+
+        delay_c.initialize(6)
+        self.assertEqual(delay_c(), 6)
+
     def test_delay_subscript(self):
         """
         Test for subscripted delay
@@ -507,6 +544,62 @@ def func2():
         self.assertNotEqual(f2_1, f2_i1)
         self.assertEqual(f2_i1, f2_0)
 
+        # test change initial condition
+        initial_f20.initialize(123)
+        self.assertEqual(initial_f20(), 123)
+
+    def test_smooth(self):
+        import pysd
+        smooth = pysd.functions.Smooth(smooth_input=lambda: 5,
+                                       smooth_time=lambda: 3,
+                                       initial_value=lambda: 4.234,
+                                       order=lambda: 3,
+                                       py_name='smooth')
+
+        smooth.initialize()
+        self.assertEqual(smooth(), 4.234)
+
+        smooth.initialize(6)
+        self.assertEqual(smooth(), 6)
+
+    def test_trend(self):
+        import pysd
+        trend = pysd.functions.Trend(trend_input=lambda: 5,
+                                     average_time=lambda: 3,
+                                     initial_trend=lambda: 4.234,
+                                     py_name='trend')
+
+        trend.initialize()
+        self.assertEqual(trend(), 4.234)
+
+        trend.initialize(6)
+        self.assertEqual(round(trend(), 8), 6)
+
+    def test_sampleiftrue(self):
+        import pysd
+
+        bool_v = False
+
+        def condition():
+            return bool_v
+
+        sif = pysd.functions.SampleIfTrue(condition=condition,
+                                          actual_value=lambda: 3,
+                                          initial_value=lambda: 4.234,
+                                          py_name='sampleiftrue')
+
+        sif.initialize()
+        self.assertEqual(sif(), 4.234)
+
+        sif.initialize(6)
+        self.assertEqual(sif(), 6)
+
+        bool_v = True
+        self.assertEqual(sif(), 3)
+
+
+class TestFunctions(unittest.TestCase):
+
     def test_sum(self):
         """
         Test for sum function
diff --git a/tests/unit_test_pysd.py b/tests/unit_test_pysd.py
index 59b1b58e..ee658b1c 100644
--- a/tests/unit_test_pysd.py
+++ b/tests/unit_test_pysd.py
@@ -1,4 +1,5 @@
 import unittest
+import os
 import warnings
 import pandas as pd
 import numpy as np
@@ -118,6 +119,10 @@ def test_run_return_timestamps(self):
         stocks = model.run(return_timestamps=5)
         self.assertEqual(stocks.index[0], 5)
 
+        timestamps = ['A', 'B']
+        with self.assertRaises(TypeError):
+            model.run(return_timestamps=timestamps)
+
     def test_run_return_timestamps_past_final_time(self):
         """ If the user enters a timestamp that is longer than the euler
         timeseries that is defined by the normal model file, should
@@ -155,7 +160,6 @@ def test_run_return_columns_step(self):
         import pysd
         model = pysd.read_vensim(test_model)
         result = model.run(return_columns='step')
-        print(result.columns)
         self.assertEqual(
             set(result.columns),
             {'Teacup Temperature', 'SAVEPER', 'Heat Loss to Room'})
@@ -181,11 +185,165 @@ def test_run_return_columns_pysafe_names(self):
         result = model.run(return_columns=return_columns)
         self.assertEqual(set(result.columns), set(return_columns))
 
+    def test_run_export_import(self):
+        import pysd
+        from pysd.tools.benchmarking import assert_frames_close
+
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            model = pysd.read_vensim(test_model)
+            stocks = model.run(return_timestamps=[0, 10, 20, 30])
+            self.assertTrue((stocks['INITIAL TIME'] == 0).all().all())
+            self.assertTrue((stocks['FINAL TIME'] == 30).all().all())
+
+            model.initialize()
+            stocks1 = model.run(return_timestamps=[0, 10], final_time=12)
+            self.assertTrue((stocks1['INITIAL TIME'] == 0).all().all())
+            self.assertTrue((stocks1['FINAL TIME'] == 12).all().all())
+            model.export('teacup12.pic')
+            model.initialize()
+            stocks2 = model.run(initial_condition='teacup12.pic',
+                                return_timestamps=[20, 30])
+            self.assertTrue((stocks2['INITIAL TIME'] == 12).all().all())
+            self.assertTrue((stocks2['FINAL TIME'] == 30).all().all())
+            stocks.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks1.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks2.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks.drop('FINAL TIME', axis=1, inplace=True)
+            stocks1.drop('FINAL TIME', axis=1, inplace=True)
+            stocks2.drop('FINAL TIME', axis=1, inplace=True)
+            os.remove('teacup12.pic')
+
+            assert_frames_close(stocks1, stocks.loc[[0, 10]])
+            assert_frames_close(stocks2, stocks.loc[[20, 30]])
+
+            # delays
+            test_delays = 'test-models/tests/delays/test_delays.mdl'
+            model = pysd.read_vensim(test_delays)
+            stocks = model.run(return_timestamps=20)
+            model.initialize()
+            model.run(return_timestamps=[], final_time=7)
+            model.export('delays7.pic')
+            stocks2 = model.run(initial_condition='delays7.pic',
+                                return_timestamps=20)
+            self.assertTrue((stocks['INITIAL TIME'] == 0).all().all())
+            self.assertTrue((stocks2['INITIAL TIME'] == 7).all().all())
+            stocks.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks2.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks.drop('FINAL TIME', axis=1, inplace=True)
+            stocks2.drop('FINAL TIME', axis=1, inplace=True)
+            os.remove('delays7.pic')
+
+            assert_frames_close(stocks2, stocks)
+
+            # delay fixed
+            test_delayf = 'test-models/tests/delay_fixed/test_delay_fixed.mdl'
+            model = pysd.read_vensim(test_delayf)
+            stocks = model.run(return_timestamps=20)
+            model.initialize()
+            model.run(return_timestamps=7)
+            model.export('delayf7.pic')
+            stocks2 = model.run(initial_condition='delayf7.pic',
+                                return_timestamps=20)
+            self.assertTrue((stocks['INITIAL TIME'] == 0).all().all())
+            self.assertTrue((stocks2['INITIAL TIME'] == 7).all().all())
+            stocks.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks2.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks.drop('FINAL TIME', axis=1, inplace=True)
+            stocks2.drop('FINAL TIME', axis=1, inplace=True)
+            os.remove('delayf7.pic')
+
+            assert_frames_close(stocks2, stocks)
+
+            # smooth
+            test_smooth = 'test-models/tests/subscripted_smooth/'\
+                          'test_subscripted_smooth.mdl'
+            model = pysd.read_vensim(test_smooth)
+            stocks = model.run(return_timestamps=20, flatten_output=True)
+            model.initialize()
+            model.run(return_timestamps=7)
+            model.export('smooth7.pic')
+            stocks2 = model.run(initial_condition='smooth7.pic',
+                                return_timestamps=20,
+                                flatten_output=True)
+            self.assertTrue((stocks['INITIAL TIME'] == 0).all().all())
+            self.assertTrue((stocks2['INITIAL TIME'] == 7).all().all())
+            stocks.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks2.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks.drop('FINAL TIME', axis=1, inplace=True)
+            stocks2.drop('FINAL TIME', axis=1, inplace=True)
+            os.remove('smooth7.pic')
+
+            assert_frames_close(stocks2, stocks)
+
+            # trend
+            test_trend = 'test-models/tests/subscripted_trend/'\
+                         'test_subscripted_trend.mdl'
+            model = pysd.read_vensim(test_trend)
+            stocks = model.run(return_timestamps=20, flatten_output=True)
+            model.initialize()
+            model.run(return_timestamps=7)
+            model.export('trend7.pic')
+            stocks2 = model.run(initial_condition='trend7.pic',
+                                return_timestamps=20,
+                                flatten_output=True)
+            self.assertTrue((stocks['INITIAL TIME'] == 0).all().all())
+            self.assertTrue((stocks2['INITIAL TIME'] == 7).all().all())
+            stocks.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks2.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks.drop('FINAL TIME', axis=1, inplace=True)
+            stocks2.drop('FINAL TIME', axis=1, inplace=True)
+            os.remove('trend7.pic')
+
+            assert_frames_close(stocks2, stocks)
+
+            # initial
+            test_initial = 'test-models/tests/initial_function/'\
+                           'test_initial.mdl'
+            model = pysd.read_vensim(test_initial)
+            stocks = model.run(return_timestamps=20)
+            model.initialize()
+            model.run(return_timestamps=7)
+            model.export('initial7.pic')
+            stocks2 = model.run(initial_condition='initial7.pic',
+                                return_timestamps=20)
+            self.assertTrue((stocks['INITIAL TIME'] == 0).all().all())
+            self.assertTrue((stocks2['INITIAL TIME'] == 7).all().all())
+            stocks.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks2.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks.drop('FINAL TIME', axis=1, inplace=True)
+            stocks2.drop('FINAL TIME', axis=1, inplace=True)
+            os.remove('initial7.pic')
+
+            assert_frames_close(stocks2, stocks)
+
+            # sample if true
+            test_sample_if_true = 'test-models/tests/sample_if_true/'\
+                                  'test_sample_if_true.mdl'
+            model = pysd.read_vensim(test_sample_if_true)
+            stocks = model.run(return_timestamps=20, flatten_output=True)
+            model.initialize()
+            model.run(return_timestamps=7)
+            model.export('sample_if_true7.pic')
+            stocks2 = model.run(initial_condition='sample_if_true7.pic',
+                                return_timestamps=20,
+                                flatten_output=True)
+            self.assertTrue((stocks['INITIAL TIME'] == 0).all().all())
+            self.assertTrue((stocks2['INITIAL TIME'] == 7).all().all())
+            stocks.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks2.drop('INITIAL TIME', axis=1, inplace=True)
+            stocks.drop('FINAL TIME', axis=1, inplace=True)
+            stocks2.drop('FINAL TIME', axis=1, inplace=True)
+            os.remove('sample_if_true7.pic')
+
+            assert_frames_close(stocks2, stocks)
+
     def test_initial_conditions_tuple_pysafe_names(self):
         import pysd
         model = pysd.read_vensim(test_model)
-        stocks = model.run(initial_condition=(3000, {'teacup_temperature': 33}),
-                           return_timestamps=list(range(3000, 3010)))
+        stocks = model.run(
+            initial_condition=(3000, {'teacup_temperature': 33}),
+            return_timestamps=list(range(3000, 3010)))
         self.assertEqual(stocks.index[0], 3000)
         self.assertEqual(stocks['Teacup Temperature'].iloc[0], 33)
 
@@ -210,7 +368,7 @@ def test_initial_conditions_current(self):
     def test_initial_condition_bad_value(self):
         import pysd
         model = pysd.read_vensim(test_model)
-        with self.assertRaises(ValueError):
+        with self.assertRaises(FileNotFoundError):
             model.run(initial_condition='bad value')
 
     def test_initial_conditions_subscripted_value_with_constant(self):
@@ -223,9 +381,18 @@ def test_initial_conditions_subscripted_value_with_constant(self):
         output = xr.DataArray([[5, 5], [5, 5], [5, 5]], coords, dims)
 
         model = pysd.read_vensim(test_model_subs)
-        res = model.run(initial_condition=(5, {'initial_values': 5}),
-                        return_columns=['Initial Values'],
-                        return_timestamps=list(range(5, 10)))
+
+        with warnings.catch_warnings(record=True) as ws:
+            res = model.run(initial_condition=(5, {'initial_values': 5}),
+                            return_columns=['Initial Values'],
+                            return_timestamps=list(range(5, 10)))
+            # use only future warnings
+            wf = [w for w in ws if issubclass(w.category, FutureWarning)]
+            self.assertEqual(len(wf), 1)
+            self.assertIn(
+                "a constant value with initial_conditions will be deprecated",
+                str(wf[0].message))
+
         self.assertTrue(output.equals(res['Initial Values'].iloc[0]))
         self.assertEqual(res.index[0], 5)
 
@@ -237,14 +404,23 @@ def test_initial_conditions_subscripted_value_with_partial_xarray(self):
                   'Second Dimension Subscript': ['Column 1', 'Column 2']}
         dims = ['One Dimensional Subscript', 'Second Dimension Subscript']
         output = xr.DataArray([[5, 3], [5, 3], [5, 3]], coords, dims)
-        input_val = xr.DataArray([5, 3],
+        input_val = xr.DataArray(
+            [5, 3],
             {'Second Dimension Subscript': ['Column 1', 'Column 2']},
             ['Second Dimension Subscript'])
 
         model = pysd.read_vensim(test_model_subs)
-        res = model.run(initial_condition=(5, {'Initial Values': input_val}),
-                        return_columns=['Initial Values'],
-                        return_timestamps=list(range(5, 10)))
+        with warnings.catch_warnings(record=True) as ws:
+            res = model.run(initial_condition=(5, {'Initial Values': input_val}),
+                            return_columns=['Initial Values'],
+                            return_timestamps=list(range(5, 10)))
+            # use only future warnings
+            wf = [w for w in ws if issubclass(w.category, FutureWarning)]
+            self.assertEqual(len(wf), 1)
+            self.assertIn(
+                "a constant value with initial_conditions will be deprecated",
+                str(wf[0].message))
+
         self.assertTrue(output.equals(res['Initial Values'].iloc[0]))
         self.assertEqual(res.index[0], 5)
 
@@ -258,23 +434,30 @@ def test_initial_conditions_subscripted_value_with_xarray(self):
         output = xr.DataArray([[5, 3], [4, 8], [9, 3]], coords, dims)
 
         model = pysd.read_vensim(test_model_subs)
-        res = model.run(initial_condition=(5, {'initial_values': output}),
-                        return_columns=['Initial Values'],
-                        return_timestamps=list(range(5, 10)))
+        with warnings.catch_warnings(record=True) as ws:
+            res = model.run(initial_condition=(5, {'initial_values': output}),
+                            return_columns=['Initial Values'],
+                            return_timestamps=list(range(5, 10)))
+            # use only future warnings
+            wf = [w for w in ws if issubclass(w.category, FutureWarning)]
+            self.assertEqual(len(wf), 1)
+            self.assertIn(
+                "a constant value with initial_conditions will be deprecated",
+                str(wf[0].message))
+
         self.assertTrue(output.equals(res['Initial Values'].iloc[0]))
         self.assertEqual(res.index[0], 5)
 
     def test_initial_conditions_subscripted_value_with_numpy_error(self):
-        import xarray as xr
         import pysd
 
         input_ = np.array([[5, 3], [4, 8], [9, 3]])
 
         model = pysd.read_vensim(test_model_subs)
-        with self.assertRaises(ValueError):
-            res = model.run(initial_condition=(5, {'initial_values': input_}),
-                            return_columns=['Initial Values'],
-                            return_timestamps=list(range(5, 10)))
+        with self.assertRaises(TypeError):
+            model.run(initial_condition=(5, {'initial_values': input_}),
+                      return_columns=['Initial Values'],
+                      return_timestamps=list(range(5, 10)))
 
     def test_set_constant_parameter(self):
         """ In response to: re: https://github.com/JamesPHoughton/pysd/issues/5"""
@@ -286,6 +469,9 @@ def test_set_constant_parameter(self):
         model.run(params={'room_temperature': 70})
         self.assertEqual(model.components.room_temperature(), 70)
 
+        with self.assertRaises(NameError):
+            model.set_components({'not_a_var': 20})
+
     def test_set_timeseries_parameter(self):
         import pysd
         model = pysd.read_vensim(test_model)
@@ -480,7 +666,7 @@ def test_set_subscripted_value_with_numpy_error(self):
         input_ = np.array([[5, 3], [4, 8], [9, 3]])
 
         model = pysd.read_vensim(test_model_subs)
-        with self.assertRaises(ValueError):
+        with self.assertRaises(TypeError):
             model.set_components({'initial_values': input_, 'final_time': 10})
 
     def test_set_subscripted_timeseries_parameter_with_constant(self):
@@ -723,29 +909,95 @@ def test_initialize(self):
         self.assertNotEqual(initial_temp, final_temp)
         self.assertEqual(initial_temp, reset_temp)
 
+    def test_initialize_order(self):
+        import pysd
+        model = pysd.load('more-tests/initialization_order/'
+                          'test_initialization_order.py')
+
+        if model._stateful_elements[0].py_name.endswith('stock_a'):
+            # we want to have stock b first always
+            model._stateful_elements.reverse()
+
+        self.assertEqual(model.components.stock_b(), 42)
+        self.assertEqual(model.components.stock_a(), 42)
+        model.components.initial_parameter = lambda: 1
+        model.initialize()
+        self.assertEqual(model.components.stock_b(), 1)
+        self.assertEqual(model.components.stock_a(), 1)
+
     def test_set_state(self):
         import pysd
         model = pysd.read_vensim(test_model)
 
         initial_temp = model.components.teacup_temperature()
-        initial_time = model.components.time()
+
+        new_time = np.random.rand()
+
+        with warnings.catch_warnings(record=True) as ws:
+            # Test that we can set with real names
+            model.set_state(new_time, {'Teacup Temperature': 500})
+            self.assertNotEqual(initial_temp, 500)
+            self.assertEqual(model.components.teacup_temperature(), 500)
+            self.assertEqual(model.components.time(), new_time)
+            # use only future warnings
+            wf = [w for w in ws if issubclass(w.category, FutureWarning)]
+            self.assertEqual(len(wf), 1)
+            self.assertIn(
+                "set_state will be deprecated, use set_initial_value instead.",
+                str(wf[0].message))
+
+        with warnings.catch_warnings(record=True) as ws:
+            # Test setting with pysafe names
+            model.set_state(new_time + 1, {'teacup_temperature': 202})
+            self.assertEqual(model.components.teacup_temperature(), 202)
+            self.assertEqual(model.components.time(), new_time + 1)
+            # use only future warnings
+            wf = [w for w in ws if issubclass(w.category, FutureWarning)]
+            self.assertEqual(len(wf), 1)
+            self.assertIn(
+                "set_state will be deprecated, use set_initial_value instead.",
+                str(wf[0].message))
+
+        with warnings.catch_warnings(record=True) as ws:
+            # Test setting with stateful object name
+            model.set_state(new_time + 2, {'_integ_teacup_temperature': 302})
+            self.assertEqual(model.components.teacup_temperature(), 302)
+            self.assertEqual(model.components.time(), new_time + 2)
+            # use only future warnings
+            wf = [w for w in ws if issubclass(w.category, FutureWarning)]
+            self.assertEqual(len(wf), 1)
+            self.assertIn(
+                "set_state will be deprecated, use set_initial_value instead.",
+                str(wf[0].message))
+
+    def test_set_initial_value(self):
+        import pysd
+        model = pysd.read_vensim(test_model)
+
+        initial_temp = model.components.teacup_temperature()
 
         new_time = np.random.rand()
 
         # Test that we can set with real names
-        model.set_state(new_time, {'Teacup Temperature': 500})
+        model.set_initial_value(new_time, {'Teacup Temperature': 500})
         self.assertNotEqual(initial_temp, 500)
         self.assertEqual(model.components.teacup_temperature(), 500)
+        self.assertEqual(model.components.time(), new_time)
 
         # Test setting with pysafe names
-        model.set_state(new_time + 1, {'teacup_temperature': 202})
+        model.set_initial_value(new_time + 1, {'teacup_temperature': 202})
         self.assertEqual(model.components.teacup_temperature(), 202)
+        self.assertEqual(model.components.time(), new_time + 1)
 
         # Test setting with stateful object name
-        model.set_state(new_time + 2, {'_integ_teacup_temperature': 302})
+        model.set_initial_value(new_time + 2, {'_integ_teacup_temperature': 302})
         self.assertEqual(model.components.teacup_temperature(), 302)
+        self.assertEqual(model.components.time(), new_time + 2)
 
-    def test_set_state_lookup(self):
+        with self.assertRaises(NameError):
+            model.set_initial_value(new_time, {'not_a_var': 500})
+
+    def test_set_initial_value_lookup(self):
         import xarray as xr
         import pysd
         model = pysd.read_vensim(test_model_look)
@@ -753,10 +1005,26 @@ def test_set_state_lookup(self):
         new_time = np.random.rand()
 
         # Test that we can set with real names
-        model.set_state(new_time, {'lookup 1d': 500})
+        with warnings.catch_warnings(record=True) as ws:
+            model.set_initial_value(new_time, {'lookup 1d': 500})
+            # use only future warnings
+            wf = [w for w in ws if issubclass(w.category, FutureWarning)]
+            self.assertEqual(len(wf), 1)
+            self.assertIn(
+                "a constant value with initial_conditions will be deprecated",
+                str(wf[0].message))
+
         self.assertEqual(model.components.lookup_1d(0), 500)
         self.assertEqual(model.components.lookup_1d(100), 500)
-        model.set_state(new_time, {'lookup 2d': 520})
+
+        with warnings.catch_warnings(record=True) as ws:
+            model.set_initial_value(new_time, {'lookup 2d': 520})
+            # use only future warnings
+            wf = [w for w in ws if issubclass(w.category, FutureWarning)]
+            self.assertEqual(len(wf), 1)
+            self.assertIn(
+                "a constant value with initial_conditions will be deprecated",
+                str(wf[0].message))
 
         expected = xr.DataArray(520, {'Rows': ['Row1', 'Row2']}, ['Rows'])
         self.assertTrue(model.components.lookup_2d(0).equals(expected) )
@@ -767,7 +1035,7 @@ def test_set_state_lookup(self):
             warnings.simplefilter("ignore")
             model.run()
 
-    def test_set_state_subscripted_value_with_constant(self):
+    def test_set_initial_value_subscripted_value_with_constant(self):
         import xarray as xr
         import pysd
 
@@ -782,19 +1050,26 @@ def test_set_state_subscripted_value_with_constant(self):
         initial_stock = model.components.stock_a()
 
         # Test that we can set with real names
-        model.set_state(new_time, {'Stock A': 500})
+        model.set_initial_value(new_time, {'Stock A': 500})
         self.assertFalse(initial_stock.equals(output_b + 500))
         self.assertTrue(model.components.stock_a().equals(output_b + 500))
 
         # Test setting with pysafe names
-        model.set_state(new_time + 1, {'stock_a': 202})
+        model.set_initial_value(new_time + 1, {'stock_a': 202})
         self.assertTrue(model.components.stock_a().equals(output_b + 202))
 
         # Test setting with stateful object name
-        model.set_state(new_time + 2, {'_integ_stock_a': 302})
+        model.set_initial_value(new_time + 2, {'_integ_stock_a': 302})
         self.assertTrue(model.components.stock_a().equals(output_b + 302))
 
-    def test_set_state_subscripted_value_with_partial_xarray(self):
+        # Test error when coords are not a subset
+        with self.assertRaises(ValueError):
+            model.set_initial_value(
+                new_time + 2,
+                {'_integ_stock_a': xr.DataArray(302, {'D': ['A', 'B']}, ['D'])}
+            )
+
+    def test_set_initial_value_subscripted_value_with_partial_xarray(self):
         import xarray as xr
         import pysd
 
@@ -822,19 +1097,19 @@ def test_set_state_subscripted_value_with_partial_xarray(self):
         initial_stock = model.components.stock_a()
 
         # Test that we can set with real names
-        model.set_state(new_time, {'Stock A': input_val1})
+        model.set_initial_value(new_time, {'Stock A': input_val1})
         self.assertFalse(initial_stock.equals(output1))
         self.assertTrue(model.components.stock_a().equals(output1))
 
         # Test setting with pysafe names
-        model.set_state(new_time + 1, {'stock_a': input_val2})
+        model.set_initial_value(new_time + 1, {'stock_a': input_val2})
         self.assertTrue(model.components.stock_a().equals(output2))
 
         # Test setting with stateful object name
-        model.set_state(new_time + 2, {'_integ_stock_a': input_val3})
+        model.set_initial_value(new_time + 2, {'_integ_stock_a': input_val3})
         self.assertTrue(model.components.stock_a().equals(output3))
 
-    def test_set_state_subscripted_value_with_xarray(self):
+    def test_set_initial_value_subscripted_value_with_xarray(self):
         import xarray as xr
         import pysd
 
@@ -851,19 +1126,19 @@ def test_set_state_subscripted_value_with_xarray(self):
         initial_stock = model.components.stock_a()
 
         # Test that we can set with real names
-        model.set_state(new_time, {'Stock A': output1})
+        model.set_initial_value(new_time, {'Stock A': output1})
         self.assertFalse(initial_stock.equals(output1))
         self.assertTrue(model.components.stock_a().equals(output1))
 
         # Test setting with pysafe names
-        model.set_state(new_time + 1, {'stock_a': output2})
+        model.set_initial_value(new_time + 1, {'stock_a': output2})
         self.assertTrue(model.components.stock_a().equals(output2))
 
         # Test setting with stateful object name
-        model.set_state(new_time + 2, {'_integ_stock_a': output3})
+        model.set_initial_value(new_time + 2, {'_integ_stock_a': output3})
         self.assertTrue(model.components.stock_a().equals(output3))
 
-    def test_set_state_subscripted_value_with_numpy_error(self):
+    def test_set_initial_value_subscripted_value_with_numpy_error(self):
         import pysd
         input1 = np.array([[5, 3], [4, 8], [9, 3]])
         input2 = np.array([[53, 43], [84, 80], [29, 63]])
@@ -874,16 +1149,16 @@ def test_set_state_subscripted_value_with_numpy_error(self):
         model = pysd.read_vensim(test_model_subs)
 
         # Test that we can set with real names
-        with self.assertRaises(ValueError):
-            model.set_state(new_time, {'Stock A': input1})
+        with self.assertRaises(TypeError):
+            model.set_initial_value(new_time, {'Stock A': input1})
 
         # Test setting with pysafe names
-        with self.assertRaises(ValueError):
-            model.set_state(new_time + 1, {'stock_a': input2})
+        with self.assertRaises(TypeError):
+            model.set_initial_value(new_time + 1, {'stock_a': input2})
 
         # Test setting with stateful object name
-        with self.assertRaises(ValueError):
-            model.set_state(new_time + 2, {'_integ_stock_a': input3})
+        with self.assertRaises(TypeError):
+            model.set_initial_value(new_time + 2, {'_integ_stock_a': input3})
 
     def test_replace_element(self):
         import pysd
@@ -974,7 +1249,15 @@ def test_set_initial_condition_for_constant_component(self):
         new_state = {'Room Temperature': 100}
         new_time = 10
 
-        model.set_initial_condition((new_time, new_state))
+        with warnings.catch_warnings(record=True) as ws:
+            model.set_initial_condition((new_time, new_state))
+            # use only future warnings
+            wf = [w for w in ws if issubclass(w.category, FutureWarning)]
+            self.assertEqual(len(wf), 1)
+            self.assertIn(
+                "a constant value with initial_conditions will be deprecated",
+                str(wf[0].message))
+
         set_temp = model.components.room_temperature()
         set_time = model.components.time()
 
@@ -984,6 +1267,25 @@ def test_set_initial_condition_for_constant_component(self):
         self.assertNotEqual(initial_time, 10, "Test definition is wrong, please change configuration")
         self.assertEqual(set_time, 10)
 
+    def test_get_args(self):
+        import pysd
+
+        model = pysd.read_vensim(test_model)
+        model2 = pysd.read_vensim(test_model_look)
+
+        self.assertEqual(model.get_args('Room Temperature'), [])
+        self.assertEqual(model.get_args('room_temperature'), [])
+        self.assertEqual(model.get_args('teacup_temperature'), [])
+        self.assertEqual(model.get_args('_integ_teacup_temperature'), [])
+
+        self.assertEqual(model2.get_args('lookup 1d'), ['x'])
+        self.assertEqual(model2.get_args('lookup_1d'), ['x'])
+        self.assertEqual(model2.get_args('lookup 2d'), ['x'])
+        self.assertEqual(model2.get_args('lookup_2d'), ['x'])
+
+        with self.assertRaises(NameError):
+            model.get_args('not_a_var')
+
     def test_get_coords(self):
         import pysd
 
@@ -1007,18 +1309,26 @@ def test_get_coords(self):
         self.assertEqual(model2.get_coords('stock_a'), coords_dims)
         self.assertEqual(model2.get_coords('_integ_stock_a'), coords_dims)
 
+        with self.assertRaises(NameError):
+            model.get_coords('not_a_var')
+
     def test__build_euler_timeseries(self):
         import pysd
         model = pysd.read_vensim(test_model)
         model.components.initial_time = lambda: 3
-        model.components.final_time = lambda: 10
-        model.components.time_step = lambda: 1
+        model.final_time = 50
+        model.time_step = 1
         model.initialize()
 
         actual = list(model._build_euler_timeseries(return_timestamps=[10]))
         expected = range(3, 11, 1)
         self.assertSequenceEqual(actual, expected)
 
+        actual = list(model._build_euler_timeseries(return_timestamps=[10],
+                                                    final_time=50))
+        expected = range(3, 51, 1)
+        self.assertSequenceEqual(actual, expected)
+
     def test__integrate(self):
         import pysd
         # Todo: think through a stronger test here...