Repository: SDXorg/pysd (language: python)

Issue:
Add support for subscripted Smoth and Trend
Something similar to subscripted delays can be done. Test models already added:
```
integration_test_vensim_pathway.py: test_subscripted_smooth
test-models/tests/subscripted_smooth/test_subscripted_smooth.mdl
integration_test_vensim_pathway.py: test_subscripted_trend
test-models/tests/subscripted_trend/test_subscripted_trend.mdl

```
Lookup Pandas Series positional arguments
@enekomartinmartinez 
Sorry for my amateuer questions.
Its the same model like i ask you a few days ago. i changed a few of my model and the new reduktion() looks like:
```
@cache.step
def zunahme():
    """
    Real Name: Zunahme
    Original Eqn: max(FaktorK*Boote - if then else (Time - INITIAL TIME< 50, 0,Reduktion(Time-INITIAL TIME)),-Boote)
    Units:
    Limits: (None, None)
    Type: component
    Subs: None

    2 Zunahme beim betriebsw. Erfolg
    """
    return np.maximum(
        faktork() * boote()
        - if_then_else(
            time() - initial_time() < 50,
            lambda: 0,
            lambda: reduktion(time() - initial_time()),
        ),
        -boote(),
    )


def reduktion(x):
    """
    Real Name: Reduktion
    Original Eqn: ( [(0,0)-(100,10)],(0,0),(50,0),(100,0))
    Units:
    Limits: (None, None)
    Type: lookup
    Subs: None

    7
    """
    return lookup(x, [0, 50, 100], [0, 0, 0])
```

how can i add a new element in python?
with this call (e.g. 120,4 or  something else)
```
temp = pd.Series(index=[0, 50,100, 120], data=[0,0,0, 4])
values = model.run(params={'Reduktion':temp}, return_columns=['TIME','Umsatzrendite', 'Fang', 'Boote', 'Fische'])
```

i got this error
TypeError: <lambda>() takes 0 positional arguments but 1 was given

`68             time() - initial_time() < 50,
     69             lambda: 0,
---> 70             lambda: reduktion(time() - initial_time()),
     71         ),
     72         -boote(),

TypeError: <lambda>() takes 0 positional arguments but 1 was given`

I guess it had to do something with the x in def reduktion(x):

Best regards


Discussion / hints:

I see the problem. That's a bug we may need to solve.

Two possible fast solutions for you @foxy-munich :
1. Make the needed change in the model file
2. Or do as explained in [docs > advanced usage](https://pysd.readthedocs.io/en/master/advanced_usage.html#replacing-model-components-with-more-complex-objects):
```

model = pysd.load(...)

def new_reduktion(x):
    return np.interp(x, [0, 50, 100, 120], [0, 0, 0, 4])

model.components.reduktion = new_reduktion
```
I hope that can help you.

Bug to correct (for me and @JamesPHoughton):
When setting a pandas.Series it is made for a constant, as the time is passed as first argument to interpolate,
[see the code](https://github.com/JamesPHoughton/pysd/blob/3597d9ef17113d8d96d1f5a083285af4cbe1b291/pysd/py_backend/functions.py#L481).
```
return lambda: np.interp(self.time(), series.index, series.values)
```

We need to check if the overwritten function is a constant or a lookup and in the second case return it with an argument instead:
```
return lambda x: np.interp(x, series.index, series.values)
```

I will busy the following weeks but I am interested in correcting this bug and I may be able to do it by the end of the month.

Ok thank you very much.

in vensim the model runs from 1950 to 1996 (timesetp 46) without a point in the lookup table of reduction. then you add a point to the lookup table of reduction for the year 1998 (timestep 48). for example 4. then for the year 2000 (timestep 50) and so on. every two timesteps a value 0 to 9 is added to the lookup table.

in python i want to do this in the form of a repeat call. Maybe there is another way to do this without changing the model.py file

UPDATE: i found a soultion. i did it like you said in your second point "Or do as explained in docs > advanced usage:"


Fix patch:
diff --git a/README.md b/README.md
index 93e165f6..41ec2f23 100644
--- a/README.md
+++ b/README.md
@@ -3,12 +3,14 @@ PySD
 [![Build Status](https://travis-ci.com/JamesPHoughton/pysd.svg?branch=master)](https://travis-ci.com/JamesPHoughton/pysd)
 [![Coverage Status](https://coveralls.io/repos/github/JamesPHoughton/pysd/badge.svg?branch=master)](https://coveralls.io/github/JamesPHoughton/pysd?branch=master)
 [![PyPI version](https://badge.fury.io/py/pysd.svg)](https://badge.fury.io/py/pysd)
+[![PyPI status](https://img.shields.io/pypi/status/pysd.svg)](https://pypi.python.org/pypi/pysd/)
+[![Py version](https://img.shields.io/pypi/pyversions/pysd.svg)](https://pypi.python.org/pypi/pysd/)
 [![Docs](https://readthedocs.org/projects/pysd/badge/?version=latest)](https://pysd.readthedocs.io/en/latest/?badge=latest)
 
 
 Simulating System Dynamics Models in Python
 
-This project is a simple library for running [System Dynamics](http://en.wikipedia.org/wiki/System_dynamics) models in python, with the purpose of improving integration of *Big Data* and *Machine Learning* into the SD workflow. 
+This project is a simple library for running [System Dynamics](http://en.wikipedia.org/wiki/System_dynamics) models in python, with the purpose of improving integration of *Big Data* and *Machine Learning* into the SD workflow.
 
 **The current version needs to run at least Python 3.7. If you need support for Python 2, please use the release here: https://github.com/JamesPHoughton/pysd/releases/tag/LastPy2**
 
@@ -32,9 +34,9 @@ If you use PySD in any published work, consider citing the [PySD Introductory Pa
 
 ### Why create a new SD simulation engine?
 
-There are a number of great SD programs out there ([Vensim](http://vensim.com/), [iThink](http://www.iseesystems.com/Softwares/Business/ithinkSoftware.aspx), [AnyLogic](http://www.anylogic.com/system-dynamics), [Insight Maker](http://insightmaker.com/), and [others](http://en.wikipedia.org/wiki/List_of_system_dynamics_software)). In order not to waste our effort, or fall victim to the [Not-Invented-Here](http://en.wikipedia.org/wiki/Not_invented_here) fallacy, we should have a very good reason for starting a new project. 
+There are a number of great SD programs out there ([Vensim](http://vensim.com/), [iThink](http://www.iseesystems.com/Softwares/Business/ithinkSoftware.aspx), [AnyLogic](http://www.anylogic.com/system-dynamics), [Insight Maker](http://insightmaker.com/), and [others](http://en.wikipedia.org/wiki/List_of_system_dynamics_software)). In order not to waste our effort, or fall victim to the [Not-Invented-Here](http://en.wikipedia.org/wiki/Not_invented_here) fallacy, we should have a very good reason for starting a new project.
 
-That reason is this: There is a whole world of computational tools being developed in the larger data science community. **System dynamicists should directly use the tools that other people are building, instead of replicating their functionality in SD specific software.** The best way to do this is to bring specific SD functionality to the domain where those other tools are being developed. 
+That reason is this: There is a whole world of computational tools being developed in the larger data science community. **System dynamicists should directly use the tools that other people are building, instead of replicating their functionality in SD specific software.** The best way to do this is to bring specific SD functionality to the domain where those other tools are being developed.
 
 This approach allows SD modelers to take advantage of the most recent developments in data science, and focus our efforts on improving the part of the stack that is unique to System Dynamics modeling.
 
@@ -54,7 +56,7 @@ You can use PySD in [R](https://www.r-project.org/) via the [PySD2R](https://git
 ### Contributors
 
 Many people have contributed to developing this project - by
-[submitting code](https://github.com/JamesPHoughton/pysd/graphs/contributors), bug reports, and advice. 
+[submitting code](https://github.com/JamesPHoughton/pysd/graphs/contributors), bug reports, and advice.
 
 Special thanks to the [sdCloud.io](http://sdcloud.io) development team, who have
 made great contributions to XMILE support, and for integrating PySD into their cloud-based model simulation environment.
diff --git a/pysd/py_backend/builder.py b/pysd/py_backend/builder.py
index 7e430b86..f85379f7 100644
--- a/pysd/py_backend/builder.py
+++ b/pysd/py_backend/builder.py
@@ -347,7 +347,7 @@ def merge_partial_elements(element_list):
 
 
 def add_stock(identifier, expression, initial_condition,
-              subs, subscript_dict):
+              subs):
     """
     Creates new model element dictionaries for the model elements associated
     with a stock.
@@ -367,9 +367,6 @@ def add_stock(identifier, expression, initial_condition,
         List of strings of subscript indices that correspond to the
         list of expressions, and collectively define the shape of the output
 
-    subscript_dict: dictionary
-        Dictionary describing the possible dimensions of the stock's subscripts
-
     Returns
     -------
     reference: string
@@ -385,18 +382,19 @@ def add_stock(identifier, expression, initial_condition,
     import_modules['functions'].add("Integ")
 
     new_structure = []
+    py_name = '_integ_%s' % identifier
 
     if len(subs) == 0:
-        stateful_py_expr = 'Integ(lambda: %s, lambda: %s)' % (
-            expression, initial_condition)
+        stateful_py_expr = "Integ(lambda: %s, lambda: %s, '%s')" % (
+            expression, initial_condition, py_name)
     else:
-        stateful_py_expr = 'Integ(lambda: _d%s_dt(), lambda: '\
-                           '_init_%s())' % (identifier, identifier)
+        stateful_py_expr = "Integ(_integ_input_%s, _integ_init_%s, '%s')" % (
+            identifier, identifier, py_name)
 
         # following elements not specified in the model file, but must exist
         # create the stock initialization element
         new_structure.append({
-            'py_name': '_init_%s' % identifier,
+            'py_name': '_integ_init_%s' % identifier,
             'real_name': 'Implicit',
             'kind': 'setup',
             'py_expr': initial_condition,
@@ -409,7 +407,7 @@ def add_stock(identifier, expression, initial_condition,
         })
 
         new_structure.append({
-            'py_name': '_d%s_dt' % identifier,
+            'py_name': '_integ_input_%s' % identifier,
             'real_name': 'Implicit',
             'kind': 'component',
             'doc': 'Provides derivative for %s function' % identifier,
@@ -422,8 +420,8 @@ def add_stock(identifier, expression, initial_condition,
         })
 
     # describe the stateful object
-    stateful = {
-        'py_name': '_integ_%s' % identifier,
+    new_structure.append({
+        'py_name': py_name,
         'real_name': 'Representation of  %s' % identifier,
         'doc': 'Integrates Expression %s' % expression,
         'py_expr': stateful_py_expr,
@@ -433,14 +431,13 @@ def add_stock(identifier, expression, initial_condition,
         'subs': '',
         'kind': 'stateful',
         'arguments': ''
-    }
+    })
 
-    new_structure.append(stateful)
-    return "%s()" % stateful['py_name'], new_structure
+    return "%s()" % py_name, new_structure
 
 
-def add_n_delay(identifier, delay_input, delay_time, initial_value, order,
-                subs, subscript_dict):
+def add_delay(identifier, delay_input, delay_time, initial_value, order,
+              subs):
     """
     Creates code to instantiate a stateful 'Delay' object,
     and provides reference to that object's output.
@@ -477,9 +474,6 @@ def add_n_delay(identifier, delay_input, delay_time, initial_value, order,
         List of strings of subscript indices that correspond to the
         list of expressions, and collectively define the shape of the output
 
-    subscript_dict: dictionary
-        Dictionary describing the possible dimensions of the stock's subscripts
-
     Returns
     -------
     reference: basestring
@@ -493,22 +487,24 @@ def add_n_delay(identifier, delay_input, delay_time, initial_value, order,
     import_modules['functions'].add("Delay")
 
     new_structure = []
+    py_name = '_delay_%s' % identifier
 
     if len(subs) == 0:
-        stateful_py_expr = 'Delay(lambda: %s, lambda: %s,'\
-                           'lambda: %s, lambda: %s)' % (
-                               delay_input, delay_time, initial_value, order)
+        stateful_py_expr = "Delay(lambda: %s, lambda: %s,"\
+                           "lambda: %s, lambda: %s, time_step, '%s')" % (
+                               delay_input, delay_time,
+                               initial_value, order, py_name)
 
     else:
-        stateful_py_expr = 'Delay(lambda: _delinput_%s(),'\
-                           'lambda: _deltime_%s(), lambda: _init_%s(),'\
-                           'lambda: %s)' % (
-                               identifier, identifier, identifier, order)
+        stateful_py_expr = "Delay(_delay_input_%s, lambda: %s, _delay_init_%s,"\
+                           "lambda: %s, time_step, '%s')" % (
+                               identifier, delay_time, identifier,
+                               order, py_name)
 
         # following elements not specified in the model file, but must exist
         # create the delay initialization element
         new_structure.append({
-            'py_name': '_init_%s' % identifier,
+            'py_name': '_delay_init_%s' % identifier,
             'real_name': 'Implicit',
             'kind': 'setup',  # not specified in the model file, but must exist
             'py_expr': initial_value,
@@ -521,20 +517,7 @@ def add_n_delay(identifier, delay_input, delay_time, initial_value, order,
         })
 
         new_structure.append({
-            'py_name': '_deltime_%s' % identifier,
-            'real_name': 'Implicit',
-            'kind': 'component',
-            'doc': 'Provides delay time for %s function' % identifier,
-            'subs': subs,
-            'unit': 'See docs for %s' % identifier,
-            'lims': 'None',
-            'eqn': 'None',
-            'py_expr': delay_time,
-            'arguments': ''
-        })
-
-        new_structure.append({
-            'py_name': '_delinput_%s' % identifier,
+            'py_name': '_delay_input_%s' % identifier,
             'real_name': 'Implicit',
             'kind': 'component',
             'doc': 'Provides input for %s function' % identifier,
@@ -547,8 +530,8 @@ def add_n_delay(identifier, delay_input, delay_time, initial_value, order,
         })
 
     # describe the stateful object
-    stateful = {
-        'py_name': '_delay_%s' % identifier,
+    new_structure.append({
+        'py_name': py_name,
         'real_name': 'Delay of %s' % delay_input,
         'doc': 'Delay time: %s \n Delay initial value %s \n Delay order %s' % (
             delay_time, initial_value, order),
@@ -559,93 +542,158 @@ def add_n_delay(identifier, delay_input, delay_time, initial_value, order,
         'subs': '',
         'kind': 'stateful',
         'arguments': ''
+    })
+
+    return "%s()" % py_name, new_structure
+
+
+def add_delay_f(identifier, delay_input, delay_time, initial_value):
+    """
+    Creates code to instantiate a stateful 'DelayFixed' object,
+    and provides reference to that object's output.
+
+    The name of the stateful object is based upon the passed in parameters,
+    so if there are multiple places where identical delay functions are
+    referenced, the translated python file will only maintain one stateful
+    object, and reference it multiple times.
+
+    Parameters
+    ----------
+    identifier: basestring
+        the python-safe name of the stock
+
+    delay_input: <string>
+        Reference to the model component that is the input to the delay
+
+    delay_time: <string>
+        Can be a number (in string format) or a reference to another model
+        element which will calculate the delay. This is calculated throughout
+        the simulation at runtime.
+
+    initial_value: <string>
+        This is used to initialize the stocks that are present in the delay.
+        We initialize the stocks with equal values so that the outflow in
+        the first timestep is equal to this value.
+
+    Returns
+    -------
+    reference: basestring
+        reference to the delay object `__call__` method, which will return
+        the output of the delay process
+
+    new_structure: list
+        list of element construction dictionaries for the builder to assemble
+
+    """
+    import_modules['functions'].add("DelayFixed")
+
+    py_name = '_delayfixed_%s' % identifier
+
+    stateful_py_expr = "DelayFixed(lambda: %s, lambda: %s,"\
+                       "lambda: %s, time_step, '%s')" % (
+                           delay_input, delay_time,
+                           initial_value, py_name)
+
+    # describe the stateful object
+    stateful = {
+        'py_name': py_name,
+        'real_name': 'Delay fixed  of %s' % delay_input,
+        'doc': 'DelayFixed time: %s \n Delay initial value %s' % (
+            delay_time, initial_value),
+        'py_expr': stateful_py_expr,
+        'unit': 'None',
+        'lims': 'None',
+        'eqn': 'None',
+        'subs': '',
+        'kind': 'stateful',
+        'arguments': ''
     }
-    new_structure.append(stateful)
 
-    return "%s()" % stateful['py_name'], new_structure
+    return "%s()" % py_name, [stateful]
 
 
-def add_sample_if_true(identifier, condition, actual_value, initial_value,
-                subs, subscript_dict):
+def add_n_delay(identifier, delay_input, delay_time, initial_value, order,
+                subs):
     """
-    Creates code to instantiate a stateful 'SampleIfTrue' object,
+    Creates code to instantiate a stateful 'DelayN' object,
     and provides reference to that object's output.
 
+    The name of the stateful object is based upon the passed in parameters,
+    so if there are multiple places where identical delay functions are
+    referenced, the translated python file will only maintain one stateful
+    object, and reference it multiple times.
+
     Parameters
     ----------
     identifier: basestring
         the python-safe name of the stock
 
-    condition: <string>
-        Reference to another model element that is the condition to the
-        'sample if true' function
+    delay_input: <string>
+        Reference to the model component that is the input to the delay
 
-    actual_value: <string>
+    delay_time: <string>
         Can be a number (in string format) or a reference to another model
-        element which is calculated throughout the simulation at runtime.
+        element which will calculate the delay. This is calculated throughout
+        the simulation at runtime.
 
     initial_value: <string>
-        This is used to initialize the state of the sample if true function.
+        This is used to initialize the stocks that are present in the delay.
+        We initialize the stocks with equal values so that the outflow in
+        the first timestep is equal to this value.
+
+    order: string
+        The number of stocks in the delay pipeline. As we construct the
+        delays at build time, this must be an integer and cannot be calculated
+        from other model components. Anything else will yield a ValueError.
 
     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
 
-    subscript_dict: dictionary
-        Dictionary describing the possible dimensions of the stock's subscripts
-
     Returns
     -------
     reference: basestring
-        reference to the sample if true object `__call__` method,
-        which will return the output of the sample if true process
+        reference to the delay object `__call__` method, which will return
+        the output of the delay process
 
     new_structure: list
         list of element construction dictionaries for the builder to assemble
 
     """
-    import_modules['functions'].add("SampleIfTrue")
+    import_modules['functions'].add("DelayN")
 
     new_structure = []
+    py_name = '_delayn_%s' % identifier
 
     if len(subs) == 0:
-        stateful_py_expr = 'SampleIfTrue(lambda: %s, lambda: %s,'\
-                           'lambda: %s)' % (condition, actual_value, initial_value)
+        stateful_py_expr = "DelayN(lambda: %s, lambda: %s,"\
+                           "lambda: %s, lambda: %s, time_step, '%s')" % (
+                               delay_input, delay_time,
+                               initial_value, order, py_name)
 
     else:
-        stateful_py_expr = 'SampleIfTrue(lambda: _condition_%s(),'\
-                           'lambda: _input_%s(), lambda: _init_%s(),)' % (
-                               identifier, identifier, identifier)
+        stateful_py_expr = "DelayN(_delayn_input_%s, lambda: %s,"\
+                           " _delayn_init_%s, lambda: %s, time_step, '%s')" % (
+                               identifier, delay_time, identifier,
+                               order, py_name)
+
         # following elements not specified in the model file, but must exist
-        # create the sample if true initialization element
+        # create the delay initialization element
         new_structure.append({
-            'py_name': '_init_%s' % identifier,
+            'py_name': '_delayn_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 value for %s function' % identifier,
-            'unit': 'See docs for %s' % identifier,
-            'lims': 'None',
-            'eqn': 'None',
-            'arguments': ''
-        })
-
-        new_structure.append({
-            'py_name': '_condition_%s' % identifier,
-            'real_name': 'Implicit',
-            'kind': 'component',
-            'doc': 'Provides condition for %s function' % identifier,
-            'subs': subs,
+            'doc': 'Provides initial conditions for %s function' % identifier,
             'unit': 'See docs for %s' % identifier,
             'lims': 'None',
             'eqn': 'None',
-            'py_expr': condition,
             'arguments': ''
         })
 
         new_structure.append({
-            'py_name': '_input_%s' % identifier,
+            'py_name': '_delayn_input_%s' % identifier,
             'real_name': 'Implicit',
             'kind': 'component',
             'doc': 'Provides input for %s function' % identifier,
@@ -653,17 +701,71 @@ def add_sample_if_true(identifier, condition, actual_value, initial_value,
             'unit': 'See docs for %s' % identifier,
             'lims': 'None',
             'eqn': 'None',
-            'py_expr': actual_value,
+            'py_expr': delay_input,
             'arguments': ''
         })
 
+    # describe the stateful object
+    new_structure.append({
+        'py_name': py_name,
+        'real_name': 'DelayN of %s' % delay_input,
+        'doc': 'DelayN time: %s \n DelayN initial value %s \n DelayN order %s' % (
+            delay_time, initial_value, order),
+        'py_expr': stateful_py_expr,
+        'unit': 'None',
+        'lims': 'None',
+        'eqn': 'None',
+        'subs': '',
+        'kind': 'stateful',
+        'arguments': ''
+    })
+
+    return "%s()" % py_name, new_structure
+
+
+def add_sample_if_true(identifier, condition, actual_value, initial_value):
+    """
+    Creates code to instantiate a stateful 'SampleIfTrue' object,
+    and provides reference to that object's output.
+
+    Parameters
+    ----------
+    identifier: basestring
+        the python-safe name of the stock
+
+    condition: <string>
+        Reference to another model element that is the condition to the
+        'sample if true' function
+
+    actual_value: <string>
+        Can be a number (in string format) or a reference to another model
+        element which is calculated throughout the simulation at runtime.
+
+    initial_value: <string>
+        This is used to initialize the state of the sample if true function.
+
+    Returns
+    -------
+    reference: basestring
+        reference to the sample if true object `__call__` method,
+        which will return the output of the sample if true process
+
+    new_structure: list
+        list of element construction dictionaries for the builder to assemble
+
+    """
+    import_modules['functions'].add("SampleIfTrue")
+
+    py_name = '_sample_if_true_%s' % identifier
+
     # describe the stateful object
     stateful = {
-        'py_name': '_sample_if_true_%s' % identifier,
+        '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': stateful_py_expr,
+        'py_expr': "SampleIfTrue(lambda: %s, lambda: %s, lambda: %s, '%s')" % (
+                    condition, actual_value, initial_value, py_name),
         'unit': 'None',
         'lims': 'None',
         'eqn': 'None',
@@ -671,13 +773,12 @@ def add_sample_if_true(identifier, condition, actual_value, initial_value,
         'kind': 'stateful',
         'arguments': ''
     }
-    new_structure.append(stateful)
 
-    return "%s()" % stateful['py_name'], new_structure
+    return "%s()" % stateful['py_name'], [stateful]
 
 
 def add_n_smooth(identifier, smooth_input, smooth_time, initial_value, order,
-                 subs, subscript_dict):
+                 subs):
     """
     Constructs stock and flow chains that implement the calculation of
     a smoothing function.
@@ -710,9 +811,6 @@ def add_n_smooth(identifier, smooth_input, smooth_time, initial_value, order,
         List of strings of subscript indices that correspond to the
         list of expressions, and collectively define the shape of the output
 
-    subscript_dict: dictionary
-        Dictionary describing the possible dimensions of the stock's subscripts
-
     Returns
     -------
     reference: basestring
@@ -725,26 +823,69 @@ def add_n_smooth(identifier, smooth_input, smooth_time, initial_value, order,
     """
     import_modules['functions'].add("Smooth")
 
-    stateful = {
-        'py_name': '_smooth_%s' % identifier,
+    new_structure = []
+    py_name = '_smooth_%s' % identifier
+
+    if len(subs) == 0:
+        stateful_py_expr = "Smooth(lambda: %s, lambda: %s,"\
+                           "lambda: %s, lambda: %s, '%s')" % (
+                               smooth_input, smooth_time, initial_value,
+                               order, py_name)
+
+    else:
+        # only need to re-dimension init and input as xarray will take care of other
+        stateful_py_expr = "Smooth(_smooth_input_%s, lambda: %s,"\
+                           " _smooth_init_%s, lambda: %s, '%s')" % (
+                               identifier, smooth_time, identifier,
+                               order, py_name)
+
+        # following elements not specified in the model file, but must exist
+        # create the delay initialization element
+        new_structure.append({
+            'py_name': '_smooth_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': ''
+        })
+
+        new_structure.append({
+            'py_name': '_smooth_input_%s' % identifier,
+            'real_name': 'Implicit',
+            'kind': 'component',
+            'doc': 'Provides input for %s function' % identifier,
+            'subs': subs,
+            'unit': 'See docs for %s' % identifier,
+            'lims': 'None',
+            'eqn': 'None',
+            'py_expr': smooth_input,
+            'arguments': ''
+        })
+
+    new_structure.append({
+        'py_name': py_name,
         'real_name': 'Smooth of %s' % smooth_input,
         'doc': 'Smooth time: %s \n Smooth initial value %s \n Smooth order %s' % (
             smooth_time, initial_value, order),
-        'py_expr': 'Smooth(lambda: %s, lambda: %s, lambda: %s, lambda: %s)' % (
-            smooth_input, smooth_time, initial_value, order),
+        '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_trend(identifier, trend_input, average_time, initial_trend,
-                subs, subscript_dict):
+                subs):
     """
     Trend.
 
@@ -764,9 +905,6 @@ def add_n_trend(identifier, trend_input, average_time, initial_trend,
         List of strings of subscript indices that correspond to the
         list of expressions, and collectively define the shape of the output
 
-    subscript_dict: dictionary
-        Dictionary describing the possible dimensions of the stock's subscripts
-
     Returns
     -------
     reference: basestring
@@ -780,22 +918,51 @@ def add_n_trend(identifier, trend_input, average_time, initial_trend,
 
     import_modules['functions'].add("Trend")
 
-    stateful = {
-        'py_name': '_trend_%s' % identifier,
+    new_structure = []
+    py_name = '_trend_%s' % identifier
+
+    if len(subs) == 0:
+        stateful_py_expr = "Trend(lambda: %s, lambda: %s,"\
+                           " lambda: %s, '%s')" % (
+                               trend_input, average_time, initial_trend,
+                               py_name)
+
+    else:
+        # only need to re-dimension init as xarray will take care of other
+        stateful_py_expr = "Trend(lambda: %s, lambda: %s,"\
+                           " _trend_init_%s, '%s')" % (
+                               trend_input, average_time, identifier, py_name)
+
+        # following elements not specified in the model file, but must exist
+        # create the delay initialization element
+        new_structure.append({
+            'py_name': '_trend_init_%s' % identifier,
+            'real_name': 'Implicit',
+            'kind': 'setup',  # not specified in the model file, but must exist
+            'py_expr': initial_trend,
+            'subs': subs,
+            'doc': 'Provides initial conditions for %s function' % identifier,
+            'unit': 'See docs for %s' % identifier,
+            'lims': 'None',
+            'eqn': 'None',
+            'arguments': ''
+        })
+
+    new_structure.append({
+        'py_name': py_name,
         'real_name': 'trend of %s' % trend_input,
         'doc': 'Trend average time: %s \n Trend initial value %s' % (
             average_time, initial_trend),
-        'py_expr': 'Trend(lambda: %s, lambda: %s, lambda: %s)' % (
-            trend_input, average_time, initial_trend),
+        '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_initial(initial_input):
@@ -820,14 +987,15 @@ def add_initial(initial_input):
     """
 
     import_modules['functions'].add("Initial")
+    py_name = utils.make_python_identifier('_initial_%s'
+                                           % initial_input)[0]
 
     stateful = {
-        'py_name': utils.make_python_identifier('_initial_%s'
-                                                % initial_input)[0],
-        'real_name': 'Smooth of %s' % initial_input,
+        'py_name': py_name,
+        'real_name': 'Initial %s' % initial_input,
         'doc': 'Returns the value taken on during the initialization phase',
-        'py_expr': 'Initial(lambda: %s)' % (
-            initial_input),
+        'py_expr': "Initial(lambda: %s, '%s')" % (
+            initial_input, py_name),
         'unit': 'None',
         'lims': 'None',
         'eqn': 'None',
@@ -1080,18 +1248,20 @@ def add_macro(macro_name, filename, arg_names, arg_vals):
     """
     import_modules['functions'].add("Macro")
 
+    py_name = '_macro_' + macro_name + '_' + '_'.join(
+        [utils.make_python_identifier(f)[0] for f in arg_vals])
     func_args = '{ %s }' % ', '.join(["'%s': lambda: %s" % (key, val)
                                       for key, val in
                                       zip(arg_names, arg_vals)])
 
     stateful = {
-        'py_name': '_macro_' + macro_name + '_' + '_'.join(
-            [utils.make_python_identifier(f)[0] for f in arg_vals]),
+        'py_name': py_name,
         'real_name': 'Macro Instantiation of ' + macro_name,
         'doc': 'Instantiates the Macro',
         'py_expr': "Macro('%s', %s, '%s',"
-                   "time_initialization=lambda: __data['time'])" % (
-                   filename, func_args, macro_name),
+                   " time_initialization=lambda: __data['time'],"
+                   " py_name='%s')" % (
+                   filename, func_args, macro_name, py_name),
         'unit': 'None',
         'lims': 'None',
         'eqn': 'None',
diff --git a/pysd/py_backend/functions.py b/pysd/py_backend/functions.py
index f2106451..9dbb2f4c 100644
--- a/pysd/py_backend/functions.py
+++ b/pysd/py_backend/functions.py
@@ -7,17 +7,15 @@
 
 import inspect
 import os
-import sys
 import re
 import random
 import warnings
-import pathlib
 from importlib.machinery import SourceFileLoader
-from ast import literal_eval
 
 import numpy as np
 import pandas as pd
 import xarray as xr
+import scipy.stats as stats
 from funcsigs import signature
 
 from . import utils
@@ -25,21 +23,6 @@
 
 from pysd._version import __version__
 
-try:
-    import scipy.stats as stats
-
-    def bounded_normal(minimum, maximum, mean, std, seed):
-        """ Implements vensim's BOUNDED NORMAL function """
-        # np.random.seed(seed)  # we could bring this back later, but for now, ignore
-        return stats.truncnorm.rvs(minimum, maximum, loc=mean, scale=std)
-
-except ImportError:
-    warnings.warn("Scipy required for functions:"
-                  "- Bounded Normal (falling back to unbounded normal)")
-
-    def bounded_normal(minimum, maximum, mean, std, seed):
-        """ Warning: using unbounded normal due to no scipy """
-        return np.random.normal(mean, std)
 
 small_vensim = 1e-6  # What is considered zero according to Vensim Help
 
@@ -51,41 +34,61 @@ class Stateful(object):
     # algebraic operations
     def __init__(self):
         self._state = None
+        self.shape_info = None
 
     def __call__(self, *args, **kwargs):
         return self.state
 
-    def initialize(self):
-        raise NotImplementedError
-
     @property
     def state(self):
         if self._state is None:
-            raise AttributeError('Attempt to call stateful element before it is initialized.')
+            raise AttributeError('Attempt to call stateful element'
+                                 + ' before it is initialized.')
         return self._state
 
     @state.setter
     def state(self, new_value):
-        self._state = new_value
+        if self.shape_info:
+            self._state = xr.DataArray(data=new_value, **self.shape_info)
+        else:
+            self._state = new_value
+
+
+class DynamicStateful(Stateful):
+
+    def __init__(self):
+        super().__init__()
 
     def update(self, state):
-        self.state = state
+        try:
+            self.state = state
+        except Exception as err:
+            raise ValueError(err.args[0] + "\n\n"
+                             + "Could not update the value of "
+                             + self.py_name)
 
 
-class Integ(Stateful):
-    def __init__(self, ddt, initial_value):
+class Integ(DynamicStateful):
+    """
+    Implements INTEG function
+    """
+    def __init__(self, ddt, initial_value, py_name="Integ object"):
         """
 
         Parameters
         ----------
         ddt: function
-            This will become an attribute of the object
-        initial_value
+          This will become an attribute of the object
+        initial_value: function
+          Initial value
+        py_name: str
+          Python name to identify the object
         """
         super().__init__()
         self.init_func = initial_value
         self.ddt = ddt
         self.shape_info = None
+        self.py_name = py_name
 
     def initialize(self):
         self.state = self.init_func()
@@ -93,19 +96,89 @@ def initialize(self):
             self.shape_info = {'dims': self.state.dims,
                                'coords': self.state.coords}
 
-    @property
-    def state(self):
-        return self._state
 
-    @state.setter
-    def state(self, new_value):
-        if self.shape_info is not None:
-            self._state = xr.DataArray(data=new_value, **self.shape_info)
+class Delay(DynamicStateful):
+    """
+    Implements DELAY function
+    """
+    # note that we could have put the `delay_input` argument as a parameter to
+    # the `__call__` function, and more closely mirrored the vensim syntax.
+    # However, people may get confused this way in thinking that they need
+    # only one delay object and can call it with various arguments to delay
+    # whatever is convenient. This method forces them to acknowledge that
+    # additional structure is being created in the delay object.
+
+    def __init__(self, delay_input, delay_time, initial_value, order,
+                 tstep=lambda: 0, py_name="Delay object"):
+        """
+
+        Parameters
+        ----------
+        delay_input: function
+        delay_time: function
+        initial_value: function
+        order: function
+        py_name: str
+          Python name to identify the object
+        """
+        super().__init__()
+        self.init_func = initial_value
+        self.delay_time_func = delay_time
+        self.input_func = delay_input
+        self.order_func = order
+        self.order = None
+        self.tstep = tstep
+        self.shape_info = None
+        self.py_name = py_name
+
+    def initialize(self):
+        order = self.order_func()
+
+        if order != int(order):
+            warnings.warn(self.py_name + '\n' +
+                          'Casting delay order '
+                          + f'from {order} to {int(order)}')
+
+        self.order = int(order)  # The order can only be set once
+        if self.order*self.tstep() > np.min(self.delay_time_func()):
+            while self.order*self.tstep() > np.min(self.delay_time_func()):
+                self.order -= 1
+            warnings.warn(self.py_name + '\n' +
+                          '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 isinstance(init_state_value, xr.DataArray):
+            # broadcast self.state
+            self.state = init_state_value.expand_dims({
+                '_delay': np.arange(self.order)}, axis=0)
         else:
-            self._state = new_value
+            self.state = xr.DataArray(
+                init_state_value,
+                {'_delay': np.arange(self.order)},
+                ['_delay'])
+        self.shape_info = {'dims': self.state.dims,
+                           'coords': self.state.coords}
+
+    def __call__(self):
+        if self.shape_info['dims'] == ('_delay',):
+            return float(self.state[-1] / self.delay_time_func())
+        else:
+            return self.state[-1].reset_coords('_delay', drop=True)\
+                / self.delay_time_func()
+
+    def ddt(self):
+        outflows = self.state / self.delay_time_func()
+        inflows = outflows.roll({'_delay': 1}, False)
+        inflows[0] = self.input_func()
+        return (inflows - outflows)*self.order
 
 
-class Delay(Stateful):
+class DelayN(DynamicStateful):
+    """
+    Implements DELAY N function
+    """
     # note that we could have put the `delay_input` argument as a parameter to
     # the `__call__` function, and more closely mirrored the vensim syntax.
     # However, people may get confused this way in thinking that they need
@@ -113,7 +186,8 @@ class Delay(Stateful):
     # whatever is convenient. This method forces them to acknowledge that
     # additional structure is being created in the delay object.
 
-    def __init__(self, delay_input, delay_time, initial_value, order):
+    def __init__(self, delay_input, delay_time, initial_value, order,
+                 tstep, py_name):
         """
 
         Parameters
@@ -122,8 +196,8 @@ def __init__(self, delay_input, delay_time, initial_value, order):
         delay_time: function
         initial_value: function
         order: function
-        coords: dictionary (optional)
-        dims: list (optional)
+        py_name: str
+          Python name to identify the object
         """
         super().__init__()
         self.init_func = initial_value
@@ -131,135 +205,263 @@ def __init__(self, delay_input, delay_time, initial_value, order):
         self.input_func = delay_input
         self.order_func = order
         self.order = None
+        self.times = None
+        self.tstep = tstep
         self.shape_info = None
+        self.py_name = py_name
 
     def initialize(self):
         order = self.order_func()
 
         if order != int(order):
-            warnings.warn('Casting delay order from %f to %i' % (
-                order, int(order)))
+            warnings.warn(self.py_name + '\n' +
+                          'Casting delay order '
+                          + f'from {order} to {int(order)}')
 
         self.order = int(order)  # The order can only be set once
+        if self.order*self.tstep() > np.min(self.delay_time_func()):
+            while self.order*self.tstep() > np.min(self.delay_time_func()):
+                self.order -= 1
+            warnings.warn(self.py_name + '\n' +
+                          'Delay time very small, casting delay order '
+                          + f'from {int(order)} to {self.order}')
 
-        init_state_value = self.init_func() * self.delay_time_func()\
-                           / self.order
+        init_state_value = self.init_func() * self.delay_time_func()
 
         if isinstance(init_state_value, xr.DataArray):
             # broadcast self.state
             self.state = init_state_value.expand_dims({
-                'delay': np.arange(self.order)}, axis=0)
-            self.shape_info = {'dims': self.state.dims,
-                               'coords': self.state.coords}
+                '_delay': np.arange(self.order)}, axis=0)
+            self.times = self.delay_time_func().expand_dims({
+                '_delay': np.arange(self.order)}, axis=0)
         else:
-            self.state = np.array([init_state_value] * self.order)
+            self.state = xr.DataArray(
+                init_state_value,
+                {'_delay': np.arange(self.order)},
+                ['_delay'])
+            self.times = xr.DataArray(
+                self.delay_time_func(),
+                {'_delay': np.arange(self.order)},
+                ['_delay'])
+        self.shape_info = {'dims': self.state.dims,
+                           'coords': self.state.coords}
 
     def __call__(self):
-        if self.shape_info:
-            return self.state[-1].reset_coords('delay', drop=True)\
-                   / (self.delay_time_func() / self.order)
+        if self.shape_info['dims'] == ('_delay',):
+            return float((self.state[-1]/self.times[0].values))
         else:
-            return self.state[-1] / (self.delay_time_func() / self.order)
+            return self.state[-1].reset_coords('_delay', drop=True)\
+                / self.times[0].reset_coords('_delay', drop=True)
 
     def ddt(self):
-        outflows = self.state / (self.delay_time_func() / self.order)
-        inflows = np.roll(outflows, 1, axis=0)
-        if self.shape_info:
-            inflows[0] = self.input_func().values
-        else:
-            inflows[0] = self.input_func()
-        return inflows - outflows
+        self.times = self.times.roll({'_delay': 1}, False)
+        self.times[0] = self.delay_time_func()
+        outflows = self.state / self.times
+        inflows = outflows.roll({'_delay': 1}, False)
+        inflows[0] = self.input_func()
+        return (inflows - outflows)*self.order
 
-class SampleIfTrue(Stateful):
-    def __init__(self, condition, actual_value, initial_value):
+
+class DelayFixed(DynamicStateful):
+    """
+    Implements DELAY FIXED function
+    """
+
+    def __init__(self, delay_input, delay_time, initial_value, tstep,
+                 py_name):
         """
 
         Parameters
         ----------
-        condition: function
-        actual_value: function
+        delay_input: function
+        delay_time: function
         initial_value: function
+        order: function
+        py_name: str
+          Python name to identify the object
         """
         super().__init__()
-        self.condition = condition
-        self.actual_value = actual_value
-        self.initial_value = initial_value
+        self.init_func = initial_value
+        self.delay_time_func = delay_time
+        self.input_func = delay_input
+        self.tstep = tstep
+        self.order = None
+        self.pointer = 0
+        self.py_name = py_name
 
     def initialize(self):
-        self.state = self.initial_value()
-        if isinstance(self.state, xr.DataArray):
-            self.shape_info = {'dims': self.state.dims,
-                               'coords': self.state.coords}
+        order = max(self.delay_time_func()/self.tstep(), 1)
+
+        if order != int(order):
+            warnings.warn(
+                self.py_name + '\n'
+                + 'Casting delay order from %f to %i' % (
+                    order, round(order + small_vensim)))
+
+        # 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()
+
+        self.state = init_state_value
+        self.pipe = [init_state_value] * self.order
 
     def __call__(self):
-        self.state = if_then_else(self.condition(), self.actual_value, lambda: self.state)
         return self.state
 
     def ddt(self):
-        return 0
+        return np.nan
 
     def update(self, state):
-        # this doesn't change once it's set up.
-        pass
+        self.pipe[self.pointer] = self.input_func()
+        self.pointer = (self.pointer + 1) % self.order
+        self.state = self.pipe[self.pointer]
 
-class Smooth(Stateful):
-    def __init__(self, smooth_input, smooth_time, initial_value, order):
+
+class Smooth(DynamicStateful):
+    """
+    Implements SMOOTH function
+    """
+    def __init__(self, smooth_input, smooth_time, initial_value, order,
+                 py_name="Smooth object"):
+        """
+
+        Parameters
+        ----------
+        smooth_input: function
+        smooth_time: function
+        initial_value: function
+        order: function
+        py_name: str
+          Python name to identify the object
+        """
         super().__init__()
         self.init_func = initial_value
         self.smooth_time_func = smooth_time
         self.input_func = smooth_input
         self.order_func = order
         self.order = None
+        self.shape_info = None
+        self.py_name = py_name
 
     def initialize(self):
         self.order = self.order_func()  # The order can only be set once
-        self.state = np.array([self.init_func()] * self.order)
+        init_state_value = self.init_func()
+        if isinstance(init_state_value, xr.DataArray):
+            # broadcast self.state
+            self.state = init_state_value.expand_dims({
+                '_smooth': np.arange(self.order)}, axis=0)
+        else:
+            self.state = xr.DataArray(
+                init_state_value,
+                {'_smooth': np.arange(self.order)},
+                ['_smooth'])
+
+        self.shape_info = {'dims': self.state.dims,
+                           'coords': self.state.coords}
 
     def __call__(self):
-        return self.state[-1]
+        if self.shape_info['dims'] == ('_smooth',):
+            return float(self.state[-1])
+        else:
+            return self.state[-1].reset_coords('_smooth', drop=True)
 
     def ddt(self):
-        targets = np.roll(self.state, 1)
+        targets = self.state.roll({'_smooth': 1}, False)
         targets[0] = self.input_func()
         return (targets - self.state) * self.order / self.smooth_time_func()
 
 
-class Trend(Stateful):
-    def __init__(self, trend_input, average_time, initial_trend):
+class Trend(DynamicStateful):
+    """
+    Implements TREND function
+    """
+    def __init__(self, trend_input, average_time, initial_trend,
+                 py_name="Trend object"):
+        """
+
+        Parameters
+        ----------
+        trend_input: function
+        average_time: function
+        initial_trend: function
+        py_name: str
+          Python name to identify the object
+        """
+
         super().__init__()
         self.init_func = initial_trend
         self.average_time_function = average_time
         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())
 
+        if isinstance(self.state, xr.DataArray):
+            self.shape_info = {'dims': self.state.dims,
+                               'coords': self.state.coords}
+
     def __call__(self):
         return zidz(self.input_func() - self.state,
-                    self.average_time_function() * abs(self.state))
+                    self.average_time_function() * np.abs(self.state))
 
     def ddt(self):
         return (self.input_func() - self.state) / self.average_time_function()
 
 
+class SampleIfTrue(Stateful):
+    def __init__(self, condition, actual_value, initial_value,
+                 py_name="SampleIfTrue object"):
+        """
+
+        Parameters
+        ----------
+        condition: function
+        actual_value: function
+        initial_value: function
+        py_name: str
+          Python name to identify the object
+        """
+        super().__init__()
+        self.condition = condition
+        self.actual_value = actual_value
+        self.initial_value = initial_value
+        self.py_name = py_name
+
+    def initialize(self):
+        self.state = self.initial_value()
+
+    def __call__(self):
+        self.state = if_then_else(self.condition(),
+                                  self.actual_value,
+                                  lambda: self.state)
+        return self.state
+
+
 class Initial(Stateful):
-    def __init__(self, func):
+    """
+    Implements INITIAL function
+    """
+    def __init__(self, func, py_name="Initial object"):
+        """
+
+        Parameters
+        ----------
+        func: function
+        py_name: str
+          Python name to identify the object
+        """
         super().__init__()
         self.func = func
+        self.py_name = py_name
 
     def initialize(self):
         self.state = self.func()
 
-    def ddt(self):
-        return 0
-
-    def update(self, state):
-        # this doesn't change once it's set up.
-        pass
-
 
-class Macro(Stateful):
+class Macro(DynamicStateful):
     """
     The Model class implements a stateful representation of the system,
     and contains the majority of methods for accessing and modifying model
@@ -271,7 +473,7 @@ class Macro(Stateful):
     """
 
     def __init__(self, py_model_file, params=None, return_func=None,
-                 time=None, time_initialization=None):
+                 time=None, time_initialization=None, py_name=None):
         """
         The model object will be created with components drawn from a
         translated python model file.
@@ -289,10 +491,11 @@ def __init__(self, py_model_file, params=None, return_func=None,
         super().__init__()
         self.time = time
         self.time_initialization = time_initialization
+        self.py_name = py_name
 
         # need a unique identifier for the imported module.
         module_name = os.path.splitext(py_model_file)[0]\
-                      + str(random.randint(0, 1000000))
+            + str(random.randint(0, 1000000))
         try:
             self.components = SourceFileLoader(module_name,
                                                py_model_file).load_module()
@@ -327,6 +530,10 @@ def __init__(self, py_model_file, params=None, return_func=None,
             getattr(self.components, name) for name in dir(self.components)
             if isinstance(getattr(self.components, name), Stateful)
         ]
+        self._dynamicstateful_elements = [
+            getattr(self.components, name) for name in dir(self.components)
+            if isinstance(getattr(self.components, name), DynamicStateful)
+        ]
         self._external_elements = [
             getattr(self.components, name) for name in dir(self.components)
             if isinstance(getattr(self.components, name), External)
@@ -397,19 +604,25 @@ def initialize(self, initialization_order=None):
             if progress:
                 remaining.difference_update(progress)
             else:
-                raise KeyError('Unresolvable Reference: Probable circular initialization' +
-                               '\n'.join([repr(e) for e in remaining]))
+                raise ValueError('Unresolvable Reference: '
+                                 + 'Probable circular initialization...\n'
+                                 + 'Not able to initialize the '
+                                 + 'following objects:\n\t'
+                                 + '\n\t'.join([e.py_name for e in remaining]))
 
     def ddt(self):
-        return np.array([component.ddt() for component in self._stateful_elements], dtype=object)
+        return np.array([component.ddt() for component
+                         in self._dynamicstateful_elements], dtype=object)
 
     @property
     def state(self):
-        return np.array([component.state for component in self._stateful_elements], dtype=object)
+        return np.array([component.state for component
+                         in self._dynamicstateful_elements], dtype=object)
 
     @state.setter
     def state(self, new_value):
-        [component.update(val) for component, val in zip(self._stateful_elements, new_value)]
+        [component.update(val) for component, val
+         in zip(self._dynamicstateful_elements, new_value)]
 
     def get_coords(self, param):
         """
@@ -418,18 +631,21 @@ def get_coords(self, param):
         >>> model.set_components('birth_rate')
         >>> model.set_components('Birth Rate')
         """
-        func_name = utils.get_value_by_insensitive_key_or_value(param,
+        func_name = utils.get_value_by_insensitive_key_or_value(
+            param,
             self.components._namespace) or param
 
-        try:
-            # TODO: make this less fragile, may crash if the component
-            # takes arguments and is a xarray
+        if not inspect.getfullargspec(getattr(self.components, func_name))[0]:
             value = getattr(self.components, func_name)()
+        else:
+            value = getattr(self.components, func_name)(0)
+
+        if isinstance(value, xr.DataArray):
             dims = list(value.dims)
             coords = {coord: list(value.coords[coord].values)
                       for coord in value.coords}
             return coords, dims
-        except Exception:
+        else:
             return None
 
     def set_components(self, params):
@@ -458,15 +674,6 @@ def set_components(self, params):
 
             func_name = utils.get_value_by_insensitive_key_or_value(key,
                 self.components._namespace)
-            try:
-                # TODO: make this less fragile, may crash if the component
-                # takes arguments and is an xarray
-                # can use the __doc__ but will not be backward compatible:
-                # dims = literal_eval(re.findall("Subs: .+",
-                # getattr(self.components, func_name).__doc__)[0][6:])
-                dims = getattr(self.components, func_name)().dims
-            except Exception:
-                dims = None
 
             if isinstance(value, np.ndarray) or isinstance(value, list):
                 raise ValueError('When setting ' + key +'\n'
@@ -479,12 +686,23 @@ def set_components(self, params):
                 raise NameError('%s is not recognized as a model component'
                                 % key)
 
+            try:
+                _, dims = self.get_coords(key)
+            except (AttributeError, TypeError):
+                dims = None
+
+            try:
+                args = inspect.getfullargspec(
+                    getattr(self.components, func_name))[0]
+            except (AttributeError, TypeError):
+                args = None
+
             if isinstance(value, pd.Series):
-                new_function = self._timeseries_component(value, dims)
+                new_function = self._timeseries_component(value, dims, args)
             elif callable(value):
                 new_function = value
             else:
-                new_function = self._constant_component(value, dims)
+                new_function = self._constant_component(value, dims, args)
 
             # this won't handle other statefuls...
             if '_integ_' + func_name in dir(self.components):
@@ -493,30 +711,71 @@ def set_components(self, params):
                               stacklevel=2)
 
             setattr(self.components, func_name, new_function)
+            self.components.cache.clean()
 
-    def _timeseries_component(self, series, dims):
+    def _timeseries_component(self, series, dims, args=[]):
         """ Internal function for creating a timeseries model element """
         # this is only called if the set_component function recognizes a
         # pandas series
         # TODO: raise a warning if extrapolating from the end of the series.
-        if isinstance(series.values[0], xr.DataArray):
-            return lambda: utils.rearrange(xr.concat(series.values,
+        if isinstance(series.values[0], xr.DataArray) and args:
+            # the argument is already given in the model when the model
+            # is called
+            return lambda x: utils.rearrange(xr.concat(
+                series.values,
+                series.index).interp(concat_dim=x).reset_coords(
+                'concat_dim', drop=True),
+                dims, self.components._subscript_dict)
+
+        elif isinstance(series.values[0], xr.DataArray):
+            # the interpolation will be time dependent
+            return lambda: utils.rearrange(xr.concat(
+                series.values,
                 series.index).interp(concat_dim=self.time()).reset_coords(
-                'concat_dim', drop=True),dims, self.components._subscript_dict)
+                'concat_dim', drop=True),
+                dims, self.components._subscript_dict)
+
+        elif args and dims:
+            # the argument is already given in the model when the model
+            # is called
+            return lambda x: utils.rearrange(
+                np.interp(x, series.index, series.values),
+                dims, self.components._subscript_dict)
+
+        elif args:
+            # the argument is already given in the model when the model
+            # is called
+            return lambda x: np.interp(x, series.index, series.values)
+
+        elif dims:
+            # the interpolation will be time dependent
+            return lambda: utils.rearrange(
+                np.interp(self.time(), series.index, series.values),
+                dims, self.components._subscript_dict)
 
         else:
-            if dims:
-                return lambda: utils.rearrange(np.interp(self.time(),
-                    series.index, series.values),
-                    dims, self.components._subscript_dict)
+            # the interpolation will be time dependent
             return lambda: np.interp(self.time(), series.index, series.values)
 
-    def _constant_component(self, value, dims):
+    def _constant_component(self, value, dims, args=[]):
         """ Internal function for creating a constant model element """
-        if dims:
-            return lambda: utils.rearrange(value, dims,
-                                           self.components._subscript_dict)
-        return lambda: value
+        if args and dims:
+            # need to pass an argument to keep consistency with the calls
+            # to the function
+            return lambda x: utils.rearrange(
+                value, dims, self.components._subscript_dict)
+
+        elif args:
+            # need to pass an argument to keep consistency with the calls
+            # to the function
+            return lambda x: value
+
+        elif dims:
+            return lambda: utils.rearrange(
+                value, dims, self.components._subscript_dict)
+
+        else:
+            return lambda: value
 
     def set_state(self, t, state):
         """ Set the system state.
@@ -536,8 +795,8 @@ def set_state(self, t, state):
 
         for key, value in state.items():
             # TODO Implement map with reference between component and stateful element?
-            component_name = utils.get_value_by_insensitive_key_or_value(key,
-                self.components._namespace)
+            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
             else:
@@ -545,35 +804,37 @@ def set_state(self, t, state):
                 stateful_name = key
 
             try:
-                # TODO make this less fragile (avoid using the __doc__)
                 if component_name[:7] == '_integ_':
                     # we need to check the original expression to retrieve
                     # the dimensions
-                    dims = literal_eval(re.findall("Subs: .+",
-                        getattr(self.components,
-                                component_name[7:]).__doc__)[0][6:])
+                    _, dims = self.get_coords(component_name[7:])
                 else:
-                    dims = literal_eval(re.findall("Subs: .+",
-                        getattr(self.components,
-                                component_name).__doc__)[0][6:])
-            except Exception:
+                    _, dims = self.get_coords(component_name)
+            except (AttributeError, TypeError):
                 dims = None
 
             if isinstance(value, np.ndarray) or isinstance(value, list):
-                raise ValueError('When setting ' + key +'\n'
+                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)')
 
+            # get the arguments of the function to know if it is a regular
+            # variable or lookup
+            args = inspect.getfullargspec(
+                getattr(self.components, component_name))[0]
+
             # 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,
+                        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
@@ -581,7 +842,8 @@ def set_state(self, t, state):
                 # Try to override component
                 try:
                     setattr(self.components, component_name,
-                            self._constant_component(value, dims))
+                            self._constant_component(value, dims, args))
+                    self.components.cache.clean()
                 except AttributeError:
                     print("'%s' has no component, assignment failed")
                     raise
@@ -1097,6 +1359,7 @@ def lookup_discrete(x, xs, ys):
 def if_then_else(condition, val_if_true, val_if_false):
     """
     Implements Vensim's IF THEN ELSE function.
+    https://www.vensim.com/documentation/20475.htm
 
     Parameters
     ----------
@@ -1125,6 +1388,8 @@ def if_then_else(condition, val_if_true, val_if_false):
 def xidz(numerator, denominator, value_if_denom_is_zero):
     """
     Implements Vensim's XIDZ function.
+    https://www.vensim.com/documentation/fn_xidz.htm
+
     This function executes a division, robust to denominator being zero.
     In the case of zero denominator, the final argument is returned.
 
@@ -1157,6 +1422,7 @@ def zidz(numerator, denominator):
     """
     This function bypasses divide-by-zero errors,
     implementing Vensim's ZIDZ function
+    https://www.vensim.com/documentation/fn_zidz.htm
 
     Parameters
     ----------
@@ -1201,6 +1467,15 @@ def active_initial(time, expr, init_val):
     else:
         return expr()
 
+def bounded_normal(minimum, maximum, mean, std, seed):
+    """
+    Implements vensim's BOUNDED NORMAL function
+    """
+    # np.random.seed(seed)
+    # we could bring this back later, but for now, ignore
+    return stats.truncnorm.rvs(minimum, maximum, loc=mean, scale=std)
+
+
 def random_0_1():
     """
     Implements Vensim's RANDOM 0 1 function.
@@ -1211,6 +1486,7 @@ def random_0_1():
     """
     return np.random.uniform(0,1)
 
+
 def random_uniform(m, x, s):
     """
     Implements Vensim's RANDOM UNIFORM function.
diff --git a/pysd/py_backend/vensim/vensim2py.py b/pysd/py_backend/vensim/vensim2py.py
index 74fc5999..228898b0 100644
--- a/pysd/py_backend/vensim/vensim2py.py
+++ b/pysd/py_backend/vensim/vensim2py.py
@@ -623,62 +623,50 @@ def parse_units(units_str):
         identifier=element['py_name'],
         expression=args[0],
         initial_condition=args[1],
-        subs=element['subs'],
-        subscript_dict=subscript_dict
+        subs=element['subs']
     ),
 
-    "delay1": lambda element, subscript_dict, args: builder.add_n_delay(
+    "delay1": lambda element, subscript_dict, args: builder.add_delay(
         identifier=element['py_name'],
         delay_input=args[0],
         delay_time=args[1],
         initial_value=args[0],
         order='1',
-        subs=element['subs'],
-        subscript_dict=subscript_dict
+        subs=element['subs']
     ),
 
-    "delay1i": lambda element, subscript_dict, args: builder.add_n_delay(
+    "delay1i": lambda element, subscript_dict, args: builder.add_delay(
         identifier=element['py_name'],
         delay_input=args[0],
         delay_time=args[1],
         initial_value=args[2],
         order='1',
-        subs=element['subs'],
-        subscript_dict=subscript_dict
+        subs=element['subs']
     ),
 
-    "delay3": lambda element, subscript_dict, args: builder.add_n_delay(
+    "delay3": lambda element, subscript_dict, args: builder.add_delay(
         identifier=element['py_name'],
         delay_input=args[0],
         delay_time=args[1],
         initial_value=args[0],
         order='3',
-        subs=element['subs'],
-        subscript_dict=subscript_dict
+        subs=element['subs']
     ),
 
-    "delay3i": lambda element, subscript_dict, args: builder.add_n_delay(
+    "delay3i": lambda element, subscript_dict, args: builder.add_delay(
         identifier=element['py_name'],
         delay_input=args[0],
         delay_time=args[1],
         initial_value=args[2],
         order='3',
-        subs=element['subs'],
-        subscript_dict=subscript_dict
+        subs=element['subs']
     ),
 
-    "delay fixed": lambda element, subscript_dict, args: builder.add_n_delay(
+    "delay fixed": lambda element, subscript_dict, args: builder.add_delay_f(
         identifier=element['py_name'],
         delay_input=args[0],
-        delay_time='time_step()' if args[1] == 'time_step()'
-                   else builder.build_function_call(
-                       functions_utils['round'],
-                       [args[1] + ' / time_step()']) + '* time_step()',
-        initial_value=args[2],
-        order='1.' if args[1] == 'time_step()'
-              else args[1] + ' / time_step()',
-        subs=element['subs'],
-        subscript_dict=subscript_dict
+        delay_time=args[1],
+        initial_value=args[2]
     ),
 
     "delay n": lambda element, subscript_dict, args: builder.add_n_delay(
@@ -687,17 +675,14 @@ def parse_units(units_str):
         delay_time=args[1],
         initial_value=args[2],
         order=args[3],
-        subs=element['subs'],
-        subscript_dict=subscript_dict
+        subs=element['subs']
     ),
 
     "sample if true": lambda element, subscript_dict, args: builder.add_sample_if_true(
         identifier=element['py_name'],
         condition=args[0],
         actual_value=args[1],
-        initial_value=args[2],
-        subs=element['subs'],
-        subscript_dict=subscript_dict
+        initial_value=args[2]
     ),
 
     "smooth": lambda element, subscript_dict, args: builder.add_n_smooth(
@@ -706,8 +691,7 @@ def parse_units(units_str):
         smooth_time=args[1],
         initial_value=args[0],
         order='1',
-        subs=element['subs'],
-        subscript_dict=subscript_dict
+        subs=element['subs']
     ),
 
     "smoothi": lambda element, subscript_dict, args: builder.add_n_smooth(
@@ -716,8 +700,7 @@ def parse_units(units_str):
         smooth_time=args[1],
         initial_value=args[2],
         order='1',
-        subs=element['subs'],
-        subscript_dict=subscript_dict
+        subs=element['subs']
     ),
 
     "smooth3": lambda element, subscript_dict, args: builder.add_n_smooth(
@@ -726,8 +709,7 @@ def parse_units(units_str):
         smooth_time=args[1],
         initial_value=args[0],
         order='3',
-        subs=element['subs'],
-        subscript_dict=subscript_dict
+        subs=element['subs']
     ),
 
     "smooth3i": lambda element, subscript_dict, args: builder.add_n_smooth(
@@ -736,8 +718,7 @@ def parse_units(units_str):
         smooth_time=args[1],
         initial_value=args[2],
         order='3',
-        subs=element['subs'],
-        subscript_dict=subscript_dict
+        subs=element['subs']
     ),
 
     "smooth n": lambda element, subscript_dict, args: builder.add_n_smooth(
@@ -746,8 +727,7 @@ def parse_units(units_str):
         smooth_time=args[1],
         initial_value=args[2],
         order=args[3],
-        subs=element['subs'],
-        subscript_dict=subscript_dict
+        subs=element['subs']
     ),
 
     "trend": lambda element, subscript_dict, args: builder.add_n_trend(
@@ -755,8 +735,8 @@ def parse_units(units_str):
         trend_input=args[0],
         average_time=args[1],
         initial_trend=args[2],
-        subs=element['subs'],
-        subscript_dict=subscript_dict),
+        subs=element['subs']
+    ),
 
     "get xls data": lambda element, subscript_dict, args: builder.add_ext_data(
         identifier=element['py_name'],
@@ -810,9 +790,9 @@ def parse_units(units_str):
 )
 
 
-def parse_general_expression(element, namespace=None, subscript_dict=None,
-                             macro_list=None, elements_subs_dict=None,
-                             subs_compatibility=None):
+def parse_general_expression(element, namespace={}, subscript_dict={},
+                             macro_list=None, elements_subs_dict={},
+                             subs_compatibility={}):
     """
     Parses a normal expression
     # its annoying that we have to construct and compile the grammar every time...
@@ -866,12 +846,6 @@ def parse_general_expression(element, namespace=None, subscript_dict=None,
        'real_name': None}])
 
     """
-    if namespace is None:
-        namespace = {}
-    if subscript_dict is None:
-        subscript_dict = {}
-    if subs_compatibility is None:
-        subs_compatibility = {}
 
     # spaces important for word-based operators
     in_ops = {
@@ -898,7 +872,9 @@ def parse_general_expression(element, namespace=None, subscript_dict=None,
 
     expression_grammar = _include_common_grammar(r"""
     expr_type = array / expr / empty
-    expr = _ pre_oper? _ (lookup_with_def / build_call / macro_call / call / lookup_call / parens / number / string / reference) _ (in_oper _ expr)?
+    expr = _ pre_oper? _ (lookup_with_def / build_call / macro_call /
+                          call / lookup_call / parens / number / string / reference)
+           _ (in_oper _ expr)?
 
     lookup_with_def = ~r"(WITH\ LOOKUP)"I _ "(" _ expr _ "," _ "(" _  ("[" ~r"[^\]]*" "]" _ ",")?  ( "(" _ expr _ "," _ expr _ ")" _ ","? _ )+ _ ")" _ ")"
 
@@ -1158,12 +1134,15 @@ def visit_subscript_list(self, n, vc):
             return ""
 
         def visit_build_call(self, n, vc):
-            call = vc[0]
-            arglist = vc[4]
+            # use only the dict with the final subscripts
+            # needed for the good working of externals
+
+            subs = {k: subscript_dict[k]
+                    for k in elements_subs_dict[element['py_name']]}
             self.kind = 'component'
-            builder_name = call.strip().lower()
-            name, structure = builders[builder_name](element, subscript_dict,
-                                                     arglist)
+            builder_name = vc[0].strip().lower()
+            name, structure = builders[builder_name](element, subs,
+                                                     vc[4])
 
             self.new_structure += structure
 
@@ -1176,9 +1155,6 @@ def visit_build_call(self, n, vc):
             elif 'data' in builder_name:
                 # External data
                 self.kind = 'component_ext_data'
-            elif builder_name == 'delay fixed':
-                warnings.warn("Delay fixed only approximates solution,"
-                              " may not give the same result as vensim")
 
             return name
 
@@ -1331,8 +1307,6 @@ def translate_section(section, macro_list, root_path):
                     subscript_dict[compatible] =\
                         subscript_dict[e['subs_compatibility'][compatible][0]]
 
-
-
     elements_subs_dict = {}
     # add model elements
     for element in model_elements:
diff --git a/pysd/py_backend/xmile/SMILE2Py.py b/pysd/py_backend/xmile/SMILE2Py.py
index 91177294..cbb81ee6 100644
--- a/pysd/py_backend/xmile/SMILE2Py.py
+++ b/pysd/py_backend/xmile/SMILE2Py.py
@@ -152,8 +152,7 @@
             delay_time = args[1],
             initial_value = args[2] if len(args) > 2 else args[0],
             order = "1",
-            subs = element['subs'],
-            subscript_dict = subscript_dict
+            subs = element['subs']
         ),
 
     "delay3": lambda element, subscript_dict, args: builder.add_n_delay(
@@ -161,8 +160,7 @@
             delay_time = args[1],
             initial_value = args[2] if len(args) > 2 else args[0],
             order = "3",
-            subs = element['subs'],
-            subscript_dict = subscript_dict
+            subs = element['subs']
         ),
 
     "delayn": lambda element, subscript_dict, args: builder.add_n_delay(
@@ -170,8 +168,7 @@
             delay_time = args[1],
             initial_value = args[2] if len(args) > 3 else args[0],
             order = args[2],
-            subs = element['subs'],
-            subscript_dict = subscript_dict
+            subs = element['subs']
         ),
 
     "smth1": lambda element, subscript_dict, args: builder.add_n_smooth(
@@ -179,8 +176,7 @@
             smooth_time = args[1],
             initial_value = args[2] if len(args) > 2 else args[0],
             order = "1",
-            subs = element['subs'],
-            subscript_dict = subscript_dict
+            subs = element['subs']
         ),
 
     "smth3": lambda element, subscript_dict, args: builder.add_n_smooth(
@@ -188,8 +184,7 @@
             smooth_time = args[1],
             initial_value = args[2] if len(args) > 2 else args[0],
             order = "3",
-            subs = element['subs'],
-            subscript_dict = subscript_dict
+            subs = element['subs']
         ),
 
     "smthn": lambda element, subscript_dict, args: builder.add_n_smooth(
@@ -197,8 +192,7 @@
             smooth_time = args[1],
             initial_value = args[2] if len(args) > 3 else args[0],
             order = args[2],
-            subs = element['subs'],
-            subscript_dict = subscript_dict
+            subs = element['subs']
         ),
 
     # "forcst" !TODO!
@@ -207,8 +201,7 @@
             trend_input = args[0],
             average_time = args[1],
             initial_trend = args[2] if len(args) > 2 else 0,
-            subs = element['subs'],
-            subscript_dict = subscript_dict
+            subs = element['subs']
         ),
 
     "init": lambda element, subscript_dict, args: builder.add_initial(args[0]),
diff --git a/pysd/py_backend/xmile/xmile2py.py b/pysd/py_backend/xmile/xmile2py.py
index c60baf85..b7614d69 100644
--- a/pysd/py_backend/xmile/xmile2py.py
+++ b/pysd/py_backend/xmile/xmile2py.py
@@ -249,8 +249,7 @@ def parse_lookup_xml_node(node):
         py_expr, new_structure = builder.add_stock(identifier=py_name,
                                                    subs=[],  # Todo later
                                                    expression=py_ddt,
-                                                   initial_condition=py_initial_value,
-                                                   subscript_dict={},  # Todo later
+                                                   initial_condition=py_initial_value
                                                    )
         element['py_expr'] = py_expr
         model_elements.append(element)
diff --git a/setup.py b/setup.py
index a232670b..ae19325a 100755
--- a/setup.py
+++ b/setup.py
@@ -3,8 +3,6 @@
 exec(open('pysd/_version.py').read())
 print(__version__)
 
-long_description = """Project Documentation: http://pysd.readthedocs.org/"""
-
 setup(
     name='pysd',
     version=__version__,
@@ -15,7 +13,8 @@
     url='https://github.com/JamesPHoughton/pysd',
     license='LICENSE.txt',
     description='System Dynamics Modeling in Python',
-    long_description=long_description,
+    long_description=open('README.md').read(),
+    long_description_content_type='text/markdown',
     keywords=['System Dynamics', 'Vensim', 'XMILE'],
     classifiers=[
         'Development Status :: 3 - Alpha',
@@ -23,6 +22,10 @@
         'Topic :: Scientific/Engineering :: Mathematics',
         'Topic :: Scientific/Engineering :: Information Analysis',
         'Intended Audience :: Science/Research',
+
+        'Programming Language :: Python :: 3.7',
+        'Programming Language :: Python :: 3.8',
+        'Programming Language :: Python :: 3.9',
     ],
     install_requires=[
         'pandas',
@@ -31,7 +34,7 @@
         'black',
         'xarray',
         'lxml',
-        'xlwt==1.2.0',        
+        'xlwt==1.2.0',
         'funcsigs',
         'pydoe',
         'xlrd',


Test patch:
diff --git a/tests/integration_test_vensim_pathway.py b/tests/integration_test_vensim_pathway.py
index aa5bb848..15e585b3 100644
--- a/tests/integration_test_vensim_pathway.py
+++ b/tests/integration_test_vensim_pathway.py
@@ -4,6 +4,7 @@
 and changes are likely to be overwritten.
 """
 
+import warnings
 import unittest
 from .test_utils import runner, assert_frames_close
 
@@ -21,7 +22,9 @@ def test_active_initial(self):
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_arguments(self):
-        output, canon = runner('test-models/tests/arguments/test_arguments.mdl')
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            output, canon = runner('test-models/tests/arguments/test_arguments.mdl')
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_builtin_max(self):
@@ -40,10 +43,11 @@ def test_constant_expressions(self):
         output, canon = runner('test-models/tests/constant_expressions/test_constant_expressions.mdl')
         assert_frames_close(output, canon, rtol=rtol)
 
-    @unittest.skip('working on it #147')
     def test_delay_fixed(self):
         # issue https://github.com/JamesPHoughton/pysd/issues/147
-        output, canon = runner('test-models/tests/delay_fixed/test_delay_fixed.mdl')
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            output, canon = runner('test-models/tests/delay_fixed/test_delay_fixed.mdl')
         assert_frames_close(output, canon, rtol=rtol)
 
     @unittest.skip('to be fixed #225')
@@ -56,13 +60,13 @@ def test_delay_parentheses(self):
         output, canon = runner('test-models/tests/delay_parentheses/test_delay_parentheses.mdl')
         assert_frames_close(output, canon, rtol=rtol)
 
-    @unittest.skip('to be fixed #147')
     def test_delay_pipeline(self):
         # issue https://github.com/JamesPHoughton/pysd/issues/147
-        output, canon = runner('test-models/tests/delay_pipeline/test_pipeline_delays.mdl')
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            output, canon = runner('test-models/tests/delay_pipeline/test_pipeline_delays.mdl')
         assert_frames_close(output, canon, rtol=rtol)
 
-    @unittest.skip('to be fixed #147')
     def test_delays(self):
         # issue https://github.com/JamesPHoughton/pysd/issues/147
         output, canon = runner('test-models/tests/delays/test_delays.mdl')
@@ -115,8 +119,6 @@ def test_get_lookups_data_3d_xls(self):
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_get_lookups_subscripted_args(self):
-        import warnings
-
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
             output, canon = runner(
@@ -133,8 +135,6 @@ def test_get_lookups_subset(self):
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_get_with_missing_values_xlsx(self):
-        import warnings
-
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
             output, canon = runner(
@@ -426,13 +426,11 @@ def test_subscripted_if_then_else(self):
         output, canon = runner('test-models/tests/subscripted_if_then_else/test_subscripted_if_then_else.mdl')
         assert_frames_close(output, canon, rtol=rtol)
 
-    @unittest.skip('working on it #226')
     def test_subscripted_smooth(self):
         # issue https://github.com/JamesPHoughton/pysd/issues/226
         output, canon = runner('test-models/tests/subscripted_smooth/test_subscripted_smooth.mdl')
         assert_frames_close(output, canon, rtol=rtol)
 
-    @unittest.skip('working on it #226')
     def test_subscripted_trend(self):
         # issue https://github.com/JamesPHoughton/pysd/issues/226
         output, canon = runner('test-models/tests/subscripted_trend/test_subscripted_trend.mdl')
diff --git a/tests/unit_test_functions.py b/tests/unit_test_functions.py
index d1543e5e..085b2c65 100644
--- a/tests/unit_test_functions.py
+++ b/tests/unit_test_functions.py
@@ -1,4 +1,5 @@
 import unittest
+import warnings
 
 import numpy as np
 import xarray as xr
@@ -351,13 +352,15 @@ def test_delay(self):
         delay_a = pysd.functions.Delay(delay_input=lambda: 5,
                                        delay_time=lambda: 3,
                                        initial_value=lambda: 4.234,
-                                       order=lambda: 3)
+                                       order=lambda: 3,
+                                       tstep=lambda: 0.5,
+                                       py_name='delay')
 
         delay_a.initialize()
 
         self.assertEqual(delay_a(), 4.234)
 
-        self.assertEqual(delay_a.ddt()[0], 5-4.234)
+        self.assertEqual(float(delay_a.ddt()[0]), (5-4.234)*3)
 
     def test_delay_subscript(self):
         """
@@ -372,16 +375,97 @@ def test_delay_subscript(self):
         xr_delay_time = xr.DataArray([[3, 2], [3, 2]], coords, dims)
 
         delay = pysd.functions.Delay(delay_input=lambda: xr_input,
-                                       delay_time=lambda: xr_delay_time,
-                                       initial_value=lambda: xr_initial,
-                                       order=lambda: 2)
+                                     delay_time=lambda: xr_delay_time,
+                                     initial_value=lambda: xr_initial,
+                                     order=lambda: 2,
+                                     tstep=lambda: 0.5,
+                                     py_name='delay')
 
         delay.initialize()
 
         self.assertTrue(delay().equals(xr_initial))
-        delay_ddt = delay.ddt()[0].reset_coords('delay', drop=True)
+        delay_ddt = delay.ddt()[0].reset_coords('_delay', drop=True)
+        print(delay_ddt)
 
-        self.assertTrue(delay_ddt.equals(xr_input-xr_initial))
+        self.assertTrue(delay_ddt.equals((xr_input-xr_initial)*2))
+
+    def test_delay_order(self):
+        import pysd
+
+        # order 3 to 2
+        delay1 = pysd.functions.Delay(delay_input=lambda: 10,
+                                      delay_time=lambda: 1,
+                                      initial_value=lambda: 0,
+                                      order=lambda: 3,
+                                      tstep=lambda: 0.4,
+                                      py_name='delay1')
+
+        # order 3 to 2
+        delay2 = pysd.functions.DelayN(delay_input=lambda: 10,
+                                       delay_time=lambda: 1,
+                                       initial_value=lambda: 0,
+                                       order=lambda: 3,
+                                       tstep=lambda: 0.4,
+                                       py_name='delay2')
+
+        # 1.5 to 1
+        delay3 = pysd.functions.Delay(delay_input=lambda: 10,
+                                      delay_time=lambda: 1,
+                                      initial_value=lambda: 0,
+                                      order=lambda: 1.5,
+                                      tstep=lambda: 0.4,
+                                      py_name='delay3')
+
+        # 1.5  to 1
+        delay4 = pysd.functions.DelayN(delay_input=lambda: 10,
+                                       delay_time=lambda: 1,
+                                       initial_value=lambda: 0,
+                                       order=lambda: 1.5,
+                                       tstep=lambda: 0.4,
+                                       py_name='delay4')
+        # 1.5 rounded to 2
+        delay5 = pysd.functions.DelayFixed(delay_input=lambda: 10,
+                                           delay_time=lambda: 0.75,
+                                           initial_value=lambda: 0,
+                                           tstep=lambda: 0.5,
+                                           py_name='delay5')
+
+        with warnings.catch_warnings(record=True) as ws:
+            delay1.initialize()
+            wu = [w for w in ws if issubclass(w.category, UserWarning)]
+            self.assertEqual(len(wu), 1)
+            self.assertIn('Delay time very small, casting delay order '
+                          + 'from 3 to 2',
+                          str(wu[0].message))
+
+        with warnings.catch_warnings(record=True) as ws:
+            delay2.initialize()
+            wu = [w for w in ws if issubclass(w.category, UserWarning)]
+            self.assertEqual(len(wu), 1)
+            self.assertIn('Delay time very small, casting delay order '
+                          + 'from 3 to 2',
+                          str(wu[0].message))
+
+        with warnings.catch_warnings(record=True) as ws:
+            delay3.initialize()
+            wu = [w for w in ws if issubclass(w.category, UserWarning)]
+            self.assertEqual(len(wu), 1)
+            self.assertIn("Casting delay order from 1.5 to 1",
+                          str(wu[0].message))
+
+        with warnings.catch_warnings(record=True) as ws:
+            delay4.initialize()
+            wu = [w for w in ws if issubclass(w.category, UserWarning)]
+            self.assertEqual(len(wu), 1)
+            self.assertIn("Casting delay order from 1.5 to 1",
+                          str(wu[0].message))
+
+        with warnings.catch_warnings(record=True) as ws:
+            delay5.initialize()
+            wu = [w for w in ws if issubclass(w.category, UserWarning)]
+            self.assertEqual(len(wu), 1)
+            self.assertIn("Casting delay order from 1.500000 to 2",
+                          str(wu[0].message))
 
     def test_initial(self):
         import pysd
diff --git a/tests/unit_test_pysd.py b/tests/unit_test_pysd.py
index 15d9c7c7..5678ec02 100644
--- a/tests/unit_test_pysd.py
+++ b/tests/unit_test_pysd.py
@@ -1,10 +1,15 @@
 import unittest
+import warnings
 import pandas as pd
 import numpy as np
+from xarray.coding.cftime_offsets import YearEnd
 
 test_model = 'test-models/samples/teacup/teacup.mdl'
-test_model_subs = 'test-models/tests/subscript_2d_arrays/test_subscript_2d_arrays.mdl'
-test_not_vensim_model = 'more-tests/Not-Vensim.txt'
+test_model_subs = 'test-models/tests/subscript_2d_arrays/'\
+                  + 'test_subscript_2d_arrays.mdl'
+test_model_look = 'test-models/tests/get_lookups_subscripted_args/'\
+                  + 'test_get_lookups_subscripted_args.mdl'
+
 
 class TestPySD(unittest.TestCase):
 
@@ -66,12 +71,12 @@ def test_load_type_error(self):
             pysd.load("type_error.py")
 
         os.remove("type_error.py")
-    
+
     def test_read_not_model_vensim(self):
         import pysd
         with self.assertRaises(ValueError):
-            model = pysd.read_vensim(test_not_vensim_model)
-        
+            pysd.read_vensim('more-tests/Not-Vensim.txt')
+
     def test_run(self):
         import pysd
         model = pysd.read_vensim(test_model)
@@ -400,10 +405,109 @@ def test_set_subscripted_value_with_xarray(self):
         res = model.run(return_columns=['Initial Values'])
         self.assertTrue(output.equals(res['Initial Values'].iloc[0]))
 
-    def test_set_subscripted_value_with_numpy_error(self):
-        import warnings
+    # TODO
+    def test_set_constant_parameter_lookup(self):
         import xarray as xr
         import pysd
+        model = pysd.read_vensim(test_model_look)
+
+        with warnings.catch_warnings():
+            # avoid warnings related to extrapolation
+            warnings.simplefilter("ignore")
+            model.set_components({'lookup_1d': 20})
+            for i in range(100):
+                self.assertEqual(model.components.lookup_1d(i), 20)
+
+            model.run(params={'lookup_1d': 70})
+            for i in range(100):
+                self.assertEqual(model.components.lookup_1d(i), 70)
+
+            model.set_components({'lookup_2d': 20})
+            for i in range(100):
+                self.assertTrue(
+                    model.components.lookup_2d(i).equals(
+                        xr.DataArray(
+                            20,
+                            {'Rows': ['Row1', 'Row2']},
+                            ['Rows'])))
+
+            model.run(params={'lookup_2d': 70})
+            for i in range(100):
+                self.assertTrue(
+                    model.components.lookup_2d(i).equals(
+                        xr.DataArray(
+                            70,
+                            {'Rows': ['Row1', 'Row2']},
+                            ['Rows'])))
+
+            xr1 = xr.DataArray(
+                [-10, 50],
+                {'Rows': ['Row1', 'Row2']},
+                ['Rows'])
+            model.set_components({'lookup_2d': xr1})
+            for i in range(100):
+                self.assertTrue(
+                    model.components.lookup_2d(i).equals(xr1))
+
+            xr2 = xr.DataArray(
+                [-100, 500],
+                {'Rows': ['Row1', 'Row2']},
+                ['Rows'])
+            model.run(params={'lookup_2d': xr2})
+            for i in range(100):
+                self.assertTrue(
+                    model.components.lookup_2d(i).equals(xr2))
+
+    def test_set_timeseries_parameter_lookup(self):
+        import xarray as xr
+        import pysd
+        model = pysd.read_vensim(test_model_look)
+        timeseries = list(range(30))
+
+        with warnings.catch_warnings():
+            # avoid warnings related to extrapolation
+            warnings.simplefilter("ignore")
+            temp_timeseries = pd.Series(
+                index=timeseries,
+                data=(50 + np.random.rand(len(timeseries)).cumsum()))
+
+            res = model.run(params={'lookup_1d': temp_timeseries},
+                            return_columns=['lookup_1d_time'],
+                            return_timestamps=timeseries)
+
+            self.assertTrue((res['lookup_1d_time'] == temp_timeseries).all())
+
+            res = model.run(params={'lookup_2d': temp_timeseries},
+                            return_columns=['lookup_2d_time'],
+                            return_timestamps=timeseries)
+
+            self.assertTrue(all([
+                a.equals(xr.DataArray(b, {'Rows': ['Row1', 'Row2']}, ['Rows']))
+                for a, b in zip(res['lookup_2d_time'].values, temp_timeseries)
+                ]))
+
+            temp_timeseries2 = pd.Series(
+                index=timeseries,
+                data=[
+                    xr.DataArray(
+                        [50 + x, 20 - y],
+                        {'Rows': ['Row1', 'Row2']}, ['Rows'])
+                        for x, y in
+                        zip(np.random.rand(len(timeseries)).cumsum(),
+                            np.random.rand(len(timeseries)).cumsum())
+                ])
+
+            res = model.run(params={'lookup_2d': temp_timeseries2},
+                            return_columns=['lookup_2d_time'],
+                            return_timestamps=timeseries)
+
+            self.assertTrue(all([
+                a.equals(b)
+                for a, b in zip(res['lookup_2d_time'].values, temp_timeseries2)
+                ]))
+
+    def test_set_subscripted_value_with_numpy_error(self):
+        import pysd
 
         input_ = np.array([[5, 3], [4, 8], [9, 3]])
 
@@ -423,7 +527,7 @@ def test_set_subscripted_timeseries_parameter_with_constant(self):
         timeseries = list(range(10))
         val_series = [50 + rd for rd in
                       np.random.rand(len(timeseries)).cumsum()]
-        xr_series = [xr.DataArray([[val,val],[val,val],[val,val]], coords, dims)
+        xr_series = [xr.DataArray(val, coords, dims)
                      for val in val_series]
 
         temp_timeseries = pd.Series(index=timeseries,
@@ -432,7 +536,8 @@ def test_set_subscripted_timeseries_parameter_with_constant(self):
                         return_columns=['initial_values'],
                         return_timestamps=timeseries)
 
-        self.assertTrue(np.all([r.equals(t) for r,t in
+        self.assertTrue(np.all([
+            r.equals(t) for r, t in
             zip(res['initial_values'].values, xr_series)]))
 
     def test_set_subscripted_timeseries_parameter_with_partial_xarray(self):
@@ -454,13 +559,13 @@ def test_set_subscripted_timeseries_parameter_with_partial_xarray(self):
         temp_timeseries = pd.Series(index=timeseries,
                                     data=val_series)
         out_series = [out_b + val
-                     for val in val_series]
+                      for val in val_series]
         model.set_components({'initial_values': temp_timeseries, 'final_time': 10})
         res = model.run(return_columns=['initial_values'])
-        self.assertTrue(np.all([r.equals(t) for r,t in
+        self.assertTrue(np.all([
+            r.equals(t) for r, t in
             zip(res['initial_values'].values, out_series)]))
 
-
     def test_set_subscripted_timeseries_parameter_with_xarray(self):
         import xarray as xr
         import pysd
@@ -480,7 +585,8 @@ def test_set_subscripted_timeseries_parameter_with_xarray(self):
                         return_columns=['initial_values'],
                         return_timestamps=timeseries)
 
-        self.assertTrue(np.all([r.equals(t) for r,t in
+        self.assertTrue(np.all([
+            r.equals(t) for r, t in
             zip(res['initial_values'].values, temp_timeseries.values)]))
 
     def test_docs(self):
@@ -671,6 +777,28 @@ def test_set_state(self):
         model.set_state(new_time + 2, {'_integ_teacup_temperature': 302})
         self.assertEqual(model.components.teacup_temperature(), 302)
 
+    def test_set_state_lookup(self):
+        import xarray as xr
+        import pysd
+        model = pysd.read_vensim(test_model_look)
+
+        new_time = np.random.rand()
+
+        # Test that we can set with real names
+        model.set_state(new_time, {'lookup 1d': 500})
+        self.assertEqual(model.components.lookup_1d(0), 500)
+        self.assertEqual(model.components.lookup_1d(100), 500)
+        model.set_state(new_time, {'lookup 2d': 520})
+
+        expected = xr.DataArray(520, {'Rows': ['Row1', 'Row2']}, ['Rows'])
+        self.assertTrue(model.components.lookup_2d(0).equals(expected) )
+        self.assertTrue(model.components.lookup_2d(100).equals(expected))
+
+        with warnings.catch_warnings():
+            # avoid warnings related to extrapolation
+            warnings.simplefilter("ignore")
+            model.run()
+
     def test_set_state_subscripted_value_with_constant(self):
         import xarray as xr
         import pysd
@@ -684,7 +812,6 @@ def test_set_state_subscripted_value_with_constant(self):
 
         model = pysd.read_vensim(test_model_subs)
         initial_stock = model.components.stock_a()
-        initial_time = model.components.time()
 
         # Test that we can set with real names
         model.set_state(new_time, {'Stock A': 500})
@@ -699,7 +826,6 @@ def test_set_state_subscripted_value_with_constant(self):
         model.set_state(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):
         import xarray as xr
         import pysd
@@ -770,8 +896,6 @@ def test_set_state_subscripted_value_with_xarray(self):
         self.assertTrue(model.components.stock_a().equals(output3))
 
     def test_set_state_subscripted_value_with_numpy_error(self):
-        import warnings
-        import xarray as xr
         import pysd
         input1 = np.array([[5, 3], [4, 8], [9, 3]])
         input2 = np.array([[53, 43], [84, 80], [29, 63]])
@@ -892,7 +1016,6 @@ 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)
 
-    #TODO
     def test_get_coords(self):
         import pysd
 
@@ -916,7 +1039,6 @@ def test_get_coords(self):
         self.assertEqual(model2.get_coords('stock_a'), coords_dims)
         self.assertEqual(model2.get_coords('_integ_stock_a'), coords_dims)
 
-
     def test__build_euler_timeseries(self):
         import pysd
         model = pysd.read_vensim(test_model)
@@ -1082,6 +1204,9 @@ def _init_outer_references(data):
         def time():
             return __data["time"]()
 
+        def time_step():
+            return 0.5
+
         def initial_time():
             return 0
 
@@ -1091,9 +1216,10 @@ def integ():
         def delay():
             return _delay_delay()
 
-        _integ_integ = Integ(lambda: 2, lambda: delay())
+        _integ_integ = Integ(lambda: 2, lambda: delay(), '_integ_integ')
 
-        _delay_delay = Delay(lambda: 2, lambda: 1, lambda: integ(), 1)
+        _delay_delay = Delay(lambda: 2, lambda: 1,
+                             lambda: integ(), 1, time_step, '_delay_delay')
         """
 
         model_main = model_main.replace("\n        ", "\n")
@@ -1101,11 +1227,37 @@ def delay():
         with open("circular.py", "w") as f:
             f.write(model_main)
 
-        with self.assertRaises(KeyError):
+        with self.assertRaises(ValueError) as err:
             pysd.load("circular.py")
 
         os.remove("circular.py")
 
+        self.assertIn('_integ_integ', str(err.exception))
+        self.assertIn('_delay_delay', str(err.exception))
+        self.assertIn('Unresolvable Reference: '
+                      + 'Probable circular initialization...\n'
+                      + 'Not able to initialize the '
+                      + 'following objects:',
+                      str(err.exception))
+
+    def test_not_able_to_update_stateful_object(self):
+        import xarray as xr
+        import pysd
+
+        integ = pysd.functions.Integ(
+            lambda: xr.DataArray([1, 2], {'Dim': ['A', 'B']}, ['Dim']),
+            lambda: xr.DataArray(0, {'Dim': ['A', 'B']}, ['Dim']),
+            'my_integ_object')
+
+        integ.initialize()
+
+        with self.assertRaises(ValueError) as err:
+            integ.update(np.array([[1, 2], [3, 4]]))
+
+        self.assertIn(
+            'Could not update the value of my_integ_object',
+            str(err.exception))
+
 
 class TestMultiRun(unittest.TestCase):
     def test_delay_reinitializes(self):
diff --git a/tests/unit_test_vensim2py.py b/tests/unit_test_vensim2py.py
index e1e0d09e..c13d831b 100644
--- a/tests/unit_test_vensim2py.py
+++ b/tests/unit_test_vensim2py.py
@@ -6,6 +6,7 @@ class TestGetFileSections(unittest.TestCase):
     def test_normal_load(self):
         """normal model file with no macros"""
         from pysd.py_backend.vensim.vensim2py import get_file_sections
+
         actual = get_file_sections(r'a~b~c| d~e~f| g~h~i|')
         expected = [{'returns': [], 'params': [], 'name': '_main_', 'string': 'a~b~c| d~e~f| g~h~i|'}]
         self.assertEqual(actual, expected)
@@ -13,6 +14,7 @@ def test_normal_load(self):
     def test_macro_only(self):
         """ Macro Only """
         from pysd.py_backend.vensim.vensim2py import get_file_sections
+
         actual = get_file_sections(':MACRO: MAC(z) a~b~c| :END OF MACRO:')
         expected = [{'returns': [], 'params': ['z'], 'name': 'MAC', 'string': 'a~b~c|'}]
         self.assertEqual(actual, expected)
@@ -20,6 +22,7 @@ def test_macro_only(self):
     def test_macro_and_model(self):
         """ basic macro and model """
         from pysd.py_backend.vensim.vensim2py import get_file_sections
+
         actual = get_file_sections(':MACRO: MAC(z) a~b~c| :END OF MACRO: d~e~f| g~h~i|')
         expected = [{'returns': [], 'params': ['z'], 'name': 'MAC', 'string': 'a~b~c|'},
                     {'returns': [], 'params': [], 'name': '_main_', 'string': 'd~e~f| g~h~i|'}]
@@ -28,6 +31,7 @@ def test_macro_and_model(self):
     def test_macro_multiple_inputs(self):
         """ macro with multiple input parameters """
         from pysd.py_backend.vensim.vensim2py import get_file_sections
+
         actual = get_file_sections(':MACRO: MAC(z, y) a~b~c| :END OF MACRO: d~e~f| g~h~i|')
         expected = [{'returns': [], 'params': ['z', 'y'], 'name': 'MAC', 'string': 'a~b~c|'},
                     {'returns': [], 'params': [], 'name': '_main_', 'string': 'd~e~f| g~h~i|'}]
@@ -36,6 +40,7 @@ def test_macro_multiple_inputs(self):
     def test_macro_with_returns(self):
         """ macro with return values """
         from pysd.py_backend.vensim.vensim2py import get_file_sections
+
         actual = get_file_sections(':MACRO: MAC(z, y :x, w) a~b~c| :END OF MACRO: d~e~f| g~h~i|')
         expected = [{'returns': ['x', 'w'],
                      'params': ['z', 'y'],
@@ -50,6 +55,7 @@ def test_macro_with_returns(self):
     def test_handle_encoding(self):
         """ Handle encoding """
         from pysd.py_backend.vensim.vensim2py import get_file_sections
+
         actual = get_file_sections(r'{UTF-8} a~b~c| d~e~f| g~h~i|')
         expected = [{'returns': [], 'params': [], 'name': '_main_', 'string': 'a~b~c| d~e~f| g~h~i|'}]
         self.assertEqual(actual, expected)
@@ -57,6 +63,7 @@ def test_handle_encoding(self):
     def test_handle_encoding_like_strings(self):
         """ Handle encoding-like strings in other places in the file """
         from pysd.py_backend.vensim.vensim2py import get_file_sections
+
         actual = get_file_sections(r'a~b~c| d~e~f{special}| g~h~i|')
         expected = [{'returns': [],
                      'params': [],
@@ -67,8 +74,10 @@ def test_handle_encoding_like_strings(self):
 
 class TestEquationStringParsing(unittest.TestCase):
     """ Tests the 'get_equation_components function """
+
     def test_basics(self):
         from pysd.py_backend.vensim.vensim2py import get_equation_components
+
         self.assertEqual(
             get_equation_components(r'constant = 25'),
             {'expr': '25', 'kind': 'component', 'subs': [], 
@@ -78,6 +87,7 @@ def test_basics(self):
     def test_equals_handling(self):
         """ Parse cases with equal signs within the expression """
         from pysd.py_backend.vensim.vensim2py import get_equation_components
+
         self.assertEqual(
             get_equation_components(r'Boolean = IF THEN ELSE(1 = 1, 1, 0)'),
             {'expr': 'IF THEN ELSE(1 = 1, 1, 0)', 'kind': 'component', 'subs': [],
@@ -105,6 +115,7 @@ def test_whitespace_handling(self):
 
     def test_subscript_definition_parsing(self):
         from pysd.py_backend.vensim.vensim2py import get_equation_components
+
         self.assertEqual(
             get_equation_components(r'''Sub1: Entry 1, Entry 2, Entry 3 '''),
             {'expr': None, 'kind': 'subdef', 'subs': ['Entry 1', 'Entry 2', 'Entry 3'],
@@ -113,6 +124,7 @@ def test_subscript_definition_parsing(self):
 
     def test_subscript_references(self):
         from pysd.py_backend.vensim.vensim2py import get_equation_components
+
         self.assertEqual(
             get_equation_components(r'constant [Sub1, Sub2] = 10, 12; 14, 16;'),
             {'expr': '10, 12; 14, 16;', 'kind': 'component', 'subs': ['Sub1', 'Sub2'],
@@ -139,6 +151,7 @@ def test_subscript_references(self):
 
     def test_lookup_definitions(self):
         from pysd.py_backend.vensim.vensim2py import get_equation_components
+
         self.assertEqual(
             get_equation_components(r'table([(0,-1)-(45,1)],(0,0),(5,0))'),
             {'expr': '([(0,-1)-(45,1)],(0,0),(5,0))', 'kind': 'lookup', 'subs': [],
@@ -155,10 +168,10 @@ def test_get_lookup(self):
         from pysd.py_backend.vensim.vensim2py import parse_lookup_expression
 
         res = parse_lookup_expression({
-	'expr': r"(GET DIRECT LOOKUPS('path2excel.xlsx', "
-                + r"'SheetName', 'index'\ , 'values'))",
-                'py_name': 'get_lookup',
-                'subs': []}, {})[1][0]
+            'expr': r"(GET DIRECT LOOKUPS('path2excel.xlsx', "
+                    + r"'SheetName', 'index'\ , 'values'))",
+            'py_name': 'get_lookup',
+            'subs': []}, {})[1][0]
 
         self.assertEqual(
             res['py_expr'],
@@ -168,6 +181,7 @@ def test_get_lookup(self):
 
     def test_pathological_names(self):
         from pysd.py_backend.vensim.vensim2py import get_equation_components
+
         self.assertEqual(
             get_equation_components(r'"silly-string" = 25'),
             {'expr': '25', 'kind': 'component', 'subs': [], 
@@ -198,11 +212,13 @@ class TestParse_general_expression(unittest.TestCase):
 
     def test_arithmetic(self):
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
+
         res = parse_general_expression({'expr': '-10^3+4'})
         self.assertEqual(res[0]['py_expr'], '-10**3+4')
 
     def test_caps_handling(self):
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
+
         res = parse_general_expression({'expr': 'Abs(-3)'})
         self.assertEqual(res[0]['py_expr'], 'abs(-3)')
 
@@ -214,6 +230,7 @@ def test_caps_handling(self):
 
     def test_function_calls(self):
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
+
         res = parse_general_expression({'expr': 'ABS(StockA)'}, {'StockA': 'stocka'})
         self.assertEqual(res[0]['py_expr'], 'abs(stocka())')
 
@@ -226,11 +243,13 @@ def test_function_calls(self):
 
     def test_id_parsing(self):
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
+
         res = parse_general_expression({'expr': 'StockA'}, {'StockA': 'stocka'})
         self.assertEqual(res[0]['py_expr'], 'stocka()')
 
     def test_number_parsing(self):
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
+
         res = parse_general_expression({'expr': '20'})
         self.assertEqual(res[0]['py_expr'], '20')
 
@@ -247,11 +266,13 @@ def test_stock_construction_function_no_subscripts(self):
         """ stock construction should create a stateful variable and reference it """
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
         from pysd.py_backend.functions import Integ
-        from pysd import functions
+
         res = parse_general_expression({'expr': 'INTEG (FlowA, -10)',
-                                      'py_name': 'test_stock',
-                                      'subs': []},
-                                     {'FlowA': 'flowa'})
+                                        'py_name': 'test_stock',
+                                        'subs': []},
+                                       {'FlowA': 'flowa'},
+                                       elements_subs_dict={'test_stock': []}
+                                       )
 
         self.assertEqual(res[1][0]['kind'], 'stateful')
         a = eval(res[1][0]['py_expr'])
@@ -263,14 +284,19 @@ def test_stock_construction_function_no_subscripts(self):
     def test_delay_construction_function_no_subscripts(self):
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
         from pysd.py_backend.functions import Delay
-        from pysd import functions
+
         res = parse_general_expression({'expr': 'DELAY1(Variable, DelayTime)',
                                         'py_name': 'test_delay',
                                         'subs': []},
                                        {'Variable': 'variable',
-                                        'DelayTime': 'delaytime'},
+                                        'DelayTime': 'delaytime',
+                                        'TIME STEP': 'time_step'},
+                                       elements_subs_dict={'test_delay': {}}
                                        )
 
+        def time_step():
+            return 0.5
+
         self.assertEqual(res[1][0]['kind'], 'stateful')
         a = eval(res[1][0]['py_expr'])
         self.assertIsInstance(a, Delay)
@@ -281,17 +307,18 @@ def test_delay_construction_function_no_subscripts(self):
     def test_smooth_construction_function_no_subscripts(self):
         """ Tests translation of 'smooth'
 
-        This translation should create a new stateful object to hold the delay elements,
-        and then pass back a reference to that value
+        This translation should create a new stateful object to hold the delay
+        elements, and then pass back a reference to that value
         """
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
         from pysd.py_backend.functions import Smooth
-        import pysd.py_backend.functions as functions  # for eval statement
+
         res = parse_general_expression({'expr': 'SMOOTH(Variable, DelayTime)',
                                         'py_name': 'test_smooth',
                                         'subs': []},
                                        {'Variable': 'variable',
                                         'DelayTime': 'delaytime'},
+                                       elements_subs_dict={'test_smooth': []}
                                        )
 
         # check stateful object creation
@@ -304,6 +331,7 @@ def test_smooth_construction_function_no_subscripts(self):
 
     def test_subscript_float_initialization(self):
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
+
         _subscript_dict = {'Dim1': ['A', 'B', 'C'],
                            'Dim2': ['D', 'E']}
         element = parse_general_expression({'expr': '3.32',
@@ -316,7 +344,7 @@ def test_subscript_float_initialization(self):
         self.assertEqual(string, 'xr.DataArray(3.32,{dim: '
                                  + '_subscript_dict[dim] for dim in '
                                  + "['Dim1']},['Dim1'])")
-        a = xr.DataArray(3.32, {dim: _subscript_dict[dim] 
+        a = xr.DataArray(3.32, {dim: _subscript_dict[dim]
                                for dim in ['Dim1']},['Dim1'])
         self.assertDictEqual({key: list(val.values) for key, val in a.coords.items()},
                               {'Dim1': ['A', 'B', 'C']})
@@ -324,6 +352,7 @@ def test_subscript_float_initialization(self):
 
     def test_subscript_1d_constant(self):
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
+
         _subscript_dict = {'Dim1': ['A', 'B', 'C'],
                            'Dim2': ['D', 'E']}
         element = parse_general_expression({'expr': '1, 2, 3',
@@ -344,6 +373,7 @@ def test_subscript_1d_constant(self):
 
     def test_subscript_2d_constant(self):
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
+
         _subscript_dict = {'Dim1': ['A', 'B', 'C'],
                            'Dim2': ['D', 'E']}
         element = parse_general_expression({'expr': '1, 2; 3, 4; 5, 6;',
@@ -359,6 +389,7 @@ def test_subscript_2d_constant(self):
 
     def test_subscript_3d_depth(self):
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
+
         _subscript_dict = {'Dim1': ['A', 'B', 'C'],
                            'Dim2': ['D', 'E']}
         element = parse_general_expression({'expr': '1, 2; 3, 4; 5, 6;',
@@ -374,6 +405,7 @@ def test_subscript_3d_depth(self):
 
     def test_subscript_reference(self):
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
+
         res = parse_general_expression({'expr': 'Var A[Dim1, Dim2]'},
                                      {'Var A': 'var_a'},
                                      {'Dim1': ['A', 'B'],
@@ -410,6 +442,7 @@ def test_subscript_reference(self):
 
     def test_subscript_ranges(self):
         from pysd.py_backend.vensim.vensim2py import parse_general_expression
+
         res = parse_general_expression({'expr': 'Var D[Range1]'},
                                      {'Var D': 'var_c'},
                                      {'Dim1': ['A', 'B', 'C', 'D', 'E', 'F'],
@@ -426,11 +459,14 @@ def test_incomplete_expression(self):
         with catch_warnings(record=True) as w:
             res = parse_general_expression({
                 'expr': 'A FUNCTION OF(Unspecified Eqn,Var A,Var B)',
-                'real_name': 'Incomplete Func'
+                'real_name': 'Incomplete Func',
+                'py_name': 'incomplete_func'
                 },
                 {'Unspecified Eqn': 'unspecified_eqn',
                  'Var A': 'var_a',
-                 'Var B': 'var_b'})
+                 'Var B': 'var_b'},
+                elements_subs_dict={'incomplete_func': []}
+                )
             self.assertEqual(len(w), 1)
             self.assertTrue('Incomplete Func has no equation specified'
                             in str(w[-1].message))