diff --git a/pysd/building/__init__.py b/pysd/building/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/pysd/building/python/__init__.py b/pysd/building/python/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/pysd/building/python/imports.py b/pysd/building/python/imports.py
new file mode 100644
index 00000000..bd69864d
--- /dev/null
+++ b/pysd/building/python/imports.py
@@ -0,0 +1,78 @@
+
+class ImportsManager():
+    """
+    Class to save the imported modules information for intelligent import
+    """
+    _external_libs = {"numpy": "np", "xarray": "xr"}
+    _external_submodules = ["scipy"]
+    _internal_libs = [
+        "functions", "statefuls", "external", "data", "lookups", "utils"
+    ]
+
+    def __init__(self):
+        self._numpy, self._xarray, self._subs = False, False, False
+        self._functions, self._statefuls, self._external, self._data,\
+            self._lookups, self._utils, self._scipy =\
+            set(), set(), set(), set(), set(), set(), set()
+
+    def add(self, module, function=None):
+        """
+        Add a function from module.
+
+        Parameters
+        ----------
+        module: str
+          module name.
+
+        function: str or None
+          function name. If None module will be set to true.
+
+        """
+        if function:
+            getattr(self, f"_{module}").add(function)
+        else:
+            setattr(self, f"_{module}", True)
+
+    def get_header(self, outfile):
+        """
+        Returns the importing information to print in the model file
+
+        Parameters
+        ----------
+        outfile: str
+            Name of the outfile to print in the header.
+
+        Returns
+        -------
+        text: str
+            Header of the translated model file.
+
+        """
+        text =\
+            f'"""\nPython model \'{outfile}\'\nTranslated using PySD\n"""\n\n'
+
+        text += "from pathlib import Path\n"
+
+        for module, shortname in self._external_libs.items():
+            if getattr(self, f"_{module}"):
+                text += f"import {module} as {shortname}\n"
+
+        for module in self._external_submodules:
+            if getattr(self, f"_{module}"):
+                text += "%(module)s import %(submodules)s\n" % {
+                    "module": module,
+                    "submodules": ", ".join(getattr(self, f"_{module}"))}
+
+        text += "\n"
+
+        for module in self._internal_libs:
+            if getattr(self, f"_{module}"):
+                text += "from pysd.py_backend.%(module)s import %(methods)s\n"\
+                        % {
+                            "module": module,
+                            "methods": ", ".join(getattr(self, f"_{module}"))}
+
+        if self._subs:
+            text += "from pysd import subs\n"
+
+        return text
diff --git a/pysd/building/python/namespace.py b/pysd/building/python/namespace.py
new file mode 100644
index 00000000..71250ec6
--- /dev/null
+++ b/pysd/building/python/namespace.py
@@ -0,0 +1,145 @@
+import re
+
+from unicodedata import normalize
+
+# used to create python safe names with the variable reserved_words
+from keyword import kwlist
+from builtins import __dir__ as bidir
+from pysd.py_backend.components import __dir__ as cdir
+from pysd.py_backend.data import __dir__ as ddir
+from pysd.py_backend.decorators import __dir__ as dedir
+from pysd.py_backend.external import __dir__ as edir
+from pysd.py_backend.functions import __dir__ as fdir
+from pysd.py_backend.statefuls import __dir__ as sdir
+from pysd.py_backend.utils import __dir__ as udir
+
+
+class NamespaceManager:
+    reserved_words = set(
+        dir() + bidir() + cdir() + ddir() + dedir() + edir() + fdir()
+        + sdir() + udir()).union(kwlist)
+
+    def __init__(self, parameters=[]):
+        self.used_words = self.reserved_words.copy()
+        self.namespace = {"Time": "time"}
+        self.cleanspace = {"time": "time"}
+        for parameter in parameters:
+            self.add_to_namespace(parameter)
+
+    def add_to_namespace(self, string):
+        self.make_python_identifier(string, add_to_namespace=True)
+
+    def make_python_identifier(self, string, prefix=None, add_to_namespace=False):
+        """
+        Takes an arbitrary string and creates a valid Python identifier.
+
+        If the input string is in the namespace, return its value.
+
+        If the python identifier created is already in the namespace,
+        but the input string is not (ie, two similar strings resolve to
+        the same python identifier)
+
+        or if the identifier is a reserved word in the reserved_words
+        list, or is a python default reserved word,
+        adds _1, or if _1 is in the namespace, _2, etc.
+
+        Parameters
+        ----------
+        string: str
+            The text to be converted into a valid python identifier.
+
+        namespace: dict
+            Map of existing translations into python safe identifiers.
+            This is to ensure that two strings are not translated into
+            the same python identifier. If string is already in the namespace
+            its value will be returned. Otherwise, namespace will be mutated
+            adding string as a new key and its value.
+
+        Returns
+        -------
+        identifier: str
+            A vaild python identifier based on the input string.
+
+        Examples
+        --------
+        >>> make_python_identifier('Capital')
+        'capital'
+
+        >>> make_python_identifier('multiple words')
+        'multiple_words'
+
+        >>> make_python_identifier('multiple     spaces')
+        'multiple_spaces'
+
+        When the name is a python keyword, add '_1' to differentiate it
+        >>> make_python_identifier('for')
+        'for_1'
+
+        Remove leading and trailing whitespace
+        >>> make_python_identifier('  whitespace  ')
+        'whitespace'
+
+        Remove most special characters outright:
+        >>> make_python_identifier('H@t tr!ck')
+        'ht_trck'
+
+        add valid string to leading digits
+        >>> make_python_identifier('123abc')
+        'nvs_123abc'
+
+        already in namespace
+        >>> make_python_identifier('Var$', namespace={'Var$': 'var'})
+        'var'
+
+        namespace conflicts
+        >>> make_python_identifier('Var@', namespace={'Var$': 'var'})
+        'var_1'
+
+        >>> make_python_identifier('Var$', namespace={'Var@': 'var',
+        ...                                           'Var%':'var_1'})
+        'var_2'
+
+        References
+        ----------
+        Identifiers must follow the convention outlined here:
+            https://docs.python.org/2/reference/lexical_analysis.html#identifiers
+
+        """
+        s = string.lower()
+        clean_s = s.replace(" ", "_")
+
+        if prefix is None and clean_s in self.cleanspace:
+            return self.cleanspace[clean_s]
+
+        # Make spaces into underscores
+        s = re.sub(r"[\s\t\n_]+", "_", s)
+
+        # remove accents, diaeresis and others ó -> o
+        s = normalize("NFD", s).encode("ascii", "ignore").decode("utf-8")
+
+        # Remove invalid characters
+        s = re.sub(r"[^0-9a-zA-Z_]", "", s)
+
+        # If leading character is not a letter add nvs_.
+        # Only letters can be leading characters.
+        if prefix is not None:
+            s = prefix + "_" + s
+        elif re.findall(r"^[0-9]", s):
+            s = "nvs_" + s
+        elif re.findall(r"^_", s):
+            s = "nvs" + s
+
+        # Check that the string is not a python identifier
+        identifier = s
+        i = 1
+        while identifier in self.used_words:
+            identifier = s + '_' + str(i)
+            i += 1
+
+        self.used_words.add(identifier)
+
+        if add_to_namespace:
+            self.namespace[string] = identifier
+            self.cleanspace[clean_s] = identifier
+
+        return identifier
diff --git a/pysd/building/python/python_builder.py b/pysd/building/python/python_builder.py
new file mode 100644
index 00000000..10b70aa5
--- /dev/null
+++ b/pysd/building/python/python_builder.py
@@ -0,0 +1,550 @@
+import textwrap
+import black
+import json
+
+from pysd.translation.structures.abstract_model import\
+    AbstractComponent, AbstractElement, AbstractModel, AbstractSection
+
+from . import visitors as vs
+from .namespace import NamespaceManager
+from .subscripts import SubscriptManager
+from .imports import ImportsManager
+from pysd._version import __version__
+
+
+class ModelBuilder:
+
+    def __init__(self, abstract_model: AbstractModel):
+        self.__dict__ = abstract_model.__dict__.copy()
+        self.sections = [
+            SectionBuilder(section)
+            for section in abstract_model.sections
+        ]
+        self.macrospace = {
+            section.name: section for section in self.sections[1:]}
+
+    def build_model(self):
+        # TODO: add special building for main
+        for section in self.sections:
+            section.macrospace = self.macrospace
+            section.build_section()
+
+        return self.sections[0].path
+
+
+class SectionBuilder:
+
+    def __init__(self, abstract_section: AbstractSection):
+        self.__dict__ = abstract_section.__dict__.copy()
+        self.root = self.path.parent
+        self.model_name = self.path.with_suffix("").name
+        self.subscripts = SubscriptManager(
+            abstract_section.subscripts, self.root)
+        self.elements = [
+            ElementBuilder(element, self)
+            for element in abstract_section.elements
+        ]
+        self.namespace = NamespaceManager(self.params)
+        self.imports = ImportsManager()
+        self.macrospace = {}
+        self.dependencies = {}
+
+    def __str__(self):
+        return "SectionBuilder " + self.path.name
+
+    def build_section(self):
+        # Create namespace
+        for element in self.elements:
+            self.namespace.add_to_namespace(element.name)
+            identifier = self.namespace.namespace[element.name]
+            element.identifier = identifier
+            self.subscripts.elements[identifier] = element.subscripts
+
+        # TODO
+        # 1. split control variables, main element, elements from other modules
+        # 2. build elements (only build 1 time!)
+        # 3. write model
+
+        for element in self.elements:
+            element.build_element()
+            self.dependencies[element.identifier] = element.dependencies
+            for subelement in element.objects.values():
+                if "calls" in subelement:
+                    self.dependencies[subelement["name"]] = subelement["calls"]
+
+        if self.split:
+            self._build_modular(self.views_dict)
+        else:
+            self._build()
+
+    def process_views_tree(self, view_name, view_content, wdir):
+        """
+        Creates a directory tree based on the elements_per_view dictionary.
+        If it's the final view, it creates a file, if not, it creates a folder.
+        """
+        if isinstance(view_content, set):
+            # will become a module
+
+            # convert subview elements names to python names
+            view_content = {
+                self.namespace.cleanspace[var] for var in view_content
+            }
+
+            # get subview elements
+            subview_elems = [
+                element for element in self.elements_remaining
+                if element.identifier in view_content
+            ]
+
+            # remove elements from remaining ones
+            [
+                self.elements_remaining.remove(element)
+                for element in subview_elems
+            ]
+
+            self._build_separate_module(subview_elems, view_name, wdir)
+
+            return sorted(view_content)
+
+        else:
+            # the current view has subviews
+            wdir = wdir.joinpath(view_name)
+            wdir.mkdir(exist_ok=True)
+            return {
+                subview_name:
+                self.process_views_tree(subview_name, subview_content, wdir)
+                for subview_name, subview_content in view_content.items()
+            }
+
+    def _build_modular(self, elements_per_view):
+        self.elements_remaining = self.elements.copy()
+        elements_per_view = self.process_views_tree(
+            "modules_" + self.model_name, elements_per_view, self.root)
+        # building main file using the build function
+        self._build_main_module(self.elements_remaining)
+
+        for file, values in {
+          "modules_%s/_modules": elements_per_view,
+          "_namespace_%s": self.namespace.namespace,
+          "_subscripts_%s": self.subscripts.subscripts,
+          "_dependencies_%s": self.dependencies}.items():
+
+            with self.root.joinpath(
+                file % self.model_name).with_suffix(
+                    ".json").open("w") as outfile:
+                json.dump(values, outfile, indent=4, sort_keys=True)
+
+    def _build_separate_module(self, elements, module_name, module_dir):
+        """
+        Constructs and writes the python representation of a specific model
+        module, when the split_views=True in the read_vensim function.
+
+        Parameters
+        ----------
+        elements: list
+            Elements belonging to the module module_name.
+
+        module_name: str
+            Name of the module
+
+        module_dir: str
+            Path of the directory where module files will be stored.
+
+        Returns
+        -------
+        None
+
+        """
+        text = textwrap.dedent('''
+        """
+        Module %(module_name)s
+        Translated using PySD version %(version)s
+        """
+        ''' % {
+            "module_name": module_name,
+            "version": __version__,
+        })
+        funcs = self._generate_functions(elements)
+        text += funcs
+        text = black.format_file_contents(
+            text, fast=True, mode=black.FileMode())
+
+        outfile_name = module_dir.joinpath(module_name + ".py")
+
+        with outfile_name.open("w", encoding="UTF-8") as out:
+            out.write(text)
+
+    def _build_main_module(self, elements):
+        """
+        Constructs and writes the python representation of the main model
+        module, when the split_views=True in the read_vensim function.
+
+        Parameters
+        ----------
+        elements: list
+            Elements belonging to the main module. Ideally, there should only be
+            the initial_time, final_time, saveper and time_step, functions, though
+            there might be others in some situations. Each element is a
+            dictionary, with the various components needed to assemble a model
+            component in python syntax. This will contain multiple entries for
+            elements that have multiple definitions in the original file, and
+            which need to be combined.
+
+        Returns
+        -------
+        None or text: None or str
+            If file_name="return" it will return the content of the output file
+            instead of saving it. It is used for testing.
+
+        """
+        # separating between control variables and rest of variables
+        control_vars, funcs = self._build_variables(elements)
+
+        self.imports.add("utils", "load_model_data")
+        self.imports.add("utils", "load_modules")
+
+        # import of needed functions and packages
+        text = self.imports.get_header(self.path.name)
+
+        # import namespace from json file
+        text += textwrap.dedent("""
+        __pysd_version__ = '%(version)s'
+
+        __data = {
+            'scope': None,
+            'time': lambda: 0
+        }
+
+        _root = Path(__file__).parent
+
+        _namespace, _subscript_dict, _dependencies, _modules = load_model_data(
+            _root, "%(model_name)s")
+        """ % {
+            "model_name": self.model_name,
+            "version": __version__
+        })
+
+        text += self._get_control_vars(control_vars)
+
+        text += textwrap.dedent("""
+            # load modules from modules_%(model_name)s directory
+            exec(load_modules("modules_%(model_name)s", _modules, _root, []))
+
+            """ % {
+                "model_name": self.model_name,
+            })
+
+        text += funcs
+        text = black.format_file_contents(text, fast=True, mode=black.FileMode())
+
+        with self.path.open("w", encoding="UTF-8") as out:
+            out.write(text)
+
+    def _build(self):
+        control_vars, funcs = self._build_variables(self.elements)
+
+        text = self.imports.get_header(self.path.name)
+        text += textwrap.dedent("""
+        __pysd_version__ = '%(version)s'
+
+        __data = {
+            'scope': None,
+            'time': lambda: 0
+        }
+
+        _root = Path(__file__).parent
+
+        _subscript_dict = %(subscript_dict)s
+
+        _namespace = %(namespace)s
+
+        _dependencies = %(dependencies)s
+        """ % {
+            "subscript_dict": repr(self.subscripts.subscripts),
+            "namespace": repr(self.namespace.namespace),
+            "dependencies": repr(self.dependencies),
+            "version": __version__,
+        })
+
+        text += self._get_control_vars(control_vars) + funcs
+
+        text = black.format_file_contents(
+            text, fast=True, mode=black.FileMode())
+
+        # this is used for testing
+        if not self.path:
+            return text
+
+        with self.path.open("w", encoding="UTF-8") as out:
+            out.write(text)
+
+    def _build_variables(self, elements):
+        """
+        Build model variables (functions) and separate then in control variables
+        and regular variables.
+
+        Returns
+        -------
+        control_vars, regular_vars: tuple, str
+            control_vars is a tuple of length 2. First element is the dictionary
+            of original control vars. Second is the string to add the control
+            variables' functions. regular_vars is the string to add the regular
+            variables' functions.
+
+        """
+        # returns of the control variables
+        control_vars_dict = {
+            "initial_time": "__data['time'].initial_time()",
+            "final_time": "__data['time'].final_time()",
+            "time_step": "__data['time'].time_step()",
+            "saveper": "__data['time'].saveper()"
+        }
+        regular_vars = []
+        control_vars = []
+
+        for element in elements:
+            if element.identifier in control_vars_dict:
+                # change the return expression in the element and update the dict
+                # with the original expression
+                control_vars_dict[element.identifier], element.expression =\
+                    element.expression, control_vars_dict[element.identifier]
+                control_vars.append(element)
+            else:
+                regular_vars.append(element)
+
+        if len(control_vars) == 0:
+            # macro objects, no control variables
+            control_vars_dict = ""
+        else:
+            control_vars_dict = """
+        _control_vars = {
+            "initial_time": lambda: %(initial_time)s,
+            "final_time": lambda: %(final_time)s,
+            "time_step": lambda: %(time_step)s,
+            "saveper": lambda: %(saveper)s
+        }
+        """ % control_vars_dict
+
+        return (control_vars_dict,
+                self._generate_functions(control_vars)),\
+            self._generate_functions(regular_vars)
+
+    def _generate_functions(self, elements):
+        """
+        Builds all model elements as functions in string format.
+        NOTE: this function calls the build_element function, which updates the
+        import_modules.
+        Therefore, it needs to be executed before the_generate_automatic_imports
+        function.
+
+        Parameters
+        ----------
+        elements: dict
+            Each element is a dictionary, with the various components needed to
+            assemble a model component in python syntax. This will contain
+            multiple entries for elements that have multiple definitions in the
+            original file, and which need to be combined.
+
+        Returns
+        -------
+        funcs: str
+            String containing all formated model functions
+
+        """
+        return "\n".join([element.build_element_out() for element in elements])
+
+    def _get_control_vars(self, control_vars):
+        """
+        Create the section of control variables
+
+        Parameters
+        ----------
+        control_vars: str
+            Functions to define control variables.
+
+        Returns
+        -------
+        text: str
+            Control variables section and header of model variables section.
+
+        """
+        text = textwrap.dedent("""
+        #######################################################################
+        #                          CONTROL VARIABLES                          #
+        #######################################################################
+        %(control_vars_dict)s
+        def _init_outer_references(data):
+            for key in data:
+                __data[key] = data[key]
+
+
+        def time():
+            return __data['time']()
+
+        """ % {"control_vars_dict": control_vars[0]})
+
+        text += control_vars[1]
+
+        text += textwrap.dedent("""
+        #######################################################################
+        #                           MODEL VARIABLES                           #
+        #######################################################################
+        """)
+
+        return text
+
+
+class SubSectionBuilder(SectionBuilder):
+    def __init__(self, abstract_section: AbstractSection):
+        pass
+        # TODO Use an intermediate class to split model, this calls could be inexistent and point to Section
+        # Namespace, subscripts and imports should point to parent section, others should remain in subsection
+
+
+class ElementBuilder:
+
+    def __init__(self, abstract_element: AbstractElement, section: SectionBuilder):
+        self.__dict__ = abstract_element.__dict__.copy()
+        self.type = None
+        self.subtype = None
+        self.arguments = getattr(self.components[0], "arguments", "")
+        self.components = [
+            ComponentBuilder(component, self, section)
+            for component in abstract_element.components
+        ]
+        self.section = section
+        self.subscripts = section.subscripts.make_merge_list(
+            [component.subscripts[0] for component in self.components])
+        self.subs_dict = section.subscripts.make_coord_dict(self.subscripts)
+        self.dependencies = {}
+        self.objects = {}
+
+    def build_element(self):
+        # TODO think better how to build the components at once to build
+        # in one declaration the external objects
+        # TODO include some kind of magic vectorization to identify patterns
+        # that can be easily vecorized (GET, expressions, Stocks...)
+        expressions = []
+        for component in self.components:
+            expr, subs = component.build_component()
+            if expr is None:
+                continue
+            else:
+                subs = {
+                    esubs: subs[csubs]
+                    for csubs, esubs in zip(subs, self.subscripts)
+                }
+                expressions.append({"expr": expr, "subs": subs})
+
+        if len(expressions) > 1:
+            # NUMPY: xrmerge would be sustitute by a multiple line definition
+            # e.g.:
+            # value = np.empty((len(dim1), len(dim2)))
+            # value[:, 0] = expression1
+            # value[:, 1] = expression2
+            # return value
+            # This allows reference to the same variable
+            # from: VAR[A] = 5; VAR[B] = 2*VAR[A]
+            # to: value[0] = 5; value[1] = 2*value[0]
+            self.section.imports.add("numpy")
+            self.pre_expression =\
+                "value = xr.DataArray(np.nan, {%s}, %s)\n" % (
+                    ", ".join("'%(dim)s': _subscript_dict['%(dim)s']" %
+                              {"dim": subs} for subs in self.subscripts),
+                    self.subscripts)
+            for expression in expressions:
+                if expression["expr"].subscripts:
+                    # get the values
+                    # NUMPY not necessary
+                    expression["expr"].lower_order(0, force_0=True)
+                    expression["expr"].expression += ".values"
+                self.pre_expression += "value.loc[%(subs)s] = %(expr)s\n" % (
+                    expression)
+            self.expression = "value"
+        else:
+            self.pre_expression = ""
+            self.expression = expressions[0]["expr"]
+
+        self.type = ", ".join(
+            set(component.type for component in self.components)
+        )
+        self.subtype = ", ".join(
+            set(component.subtype for component in self.components)
+        )
+
+    def build_element_out(self):
+        """
+        Returns a string that has processed a single element dictionary.
+
+        Returns
+        -------
+        func: str
+            The function to write in the model file.
+
+        """
+        # TODO: merge with the previous build to do all at once
+        contents = self.pre_expression + "return %s" % self.expression
+
+        self.subs_dec = ""
+        self.subs_doc = "None"
+
+        if self.subscripts:
+            # We add the list of the subs to the __doc__ of the function
+            # this will give more information to the user and make possible
+            # to rewrite subscripted values with model.run(params=X) or
+            # model.run(initial_condition=(n,x))
+            self.subs_doc = "%s" % self.subscripts
+            self.subs_dec =\
+                "@subs(%s, _subscript_dict)" % self.subscripts
+            self.section.imports.add("subs")
+
+        objects = "\n\n".join([
+            value["expression"] for value in self.objects.values()
+            if value["expression"] is not None
+        ])
+
+        indent = 12
+
+        self.contents = contents.replace("\n", "\n" + " " * (indent+4))
+        self.objects = objects.replace("\n", "\n" + " " * indent)
+
+        # convert newline indicator and add expected level of indentation
+        # TODO check if this is neccessary
+        self.documentation = self.documentation.replace(
+            "\\", "\n").replace("\n", "\n" + "" * indent)
+
+        return textwrap.dedent('''
+            %(subs_dec)s
+            def %(identifier)s(%(arguments)s):
+                """
+                Real Name: %(name)s
+                Original Eqn:
+                Units: %(units)s
+                Limits: %(range)s
+                Type: %(type)s
+                Subtype: %(subtype)s
+                Subs: %(subscripts)s
+
+                %(documentation)s
+                """
+                %(contents)s
+
+
+            %(objects)s
+            ''' % self.__dict__)
+
+
+class ComponentBuilder:
+
+    def __init__(self, abstract_component: AbstractComponent,
+                 element: ElementBuilder, section: SectionBuilder):
+        self.__dict__ = abstract_component.__dict__.copy()
+        self.element = element
+        self.section = section
+        if not hasattr(self, "keyword"):
+            self.keyword = None
+
+    def build_component(self):
+        self.subscripts_dict = self.section.subscripts.make_coord_dict(
+            self.subscripts[0])
+        return (vs.ASTVisitor(self).visit(), self.subscripts_dict)
diff --git a/pysd/building/python/python_functions.py b/pysd/building/python/python_functions.py
new file mode 100644
index 00000000..bdde7130
--- /dev/null
+++ b/pysd/building/python/python_functions.py
@@ -0,0 +1,88 @@
+
+# functions that can be diretcly applied over an array
+functionspace = {
+    # directly build functions without dependencies
+    "elmcount": ("len(_subscript_dict['%(0)s'])", None),
+
+    # directly build numpy based functions
+    "abs": ("np.abs(%(0)s)", ("numpy",)),
+    "min": ("np.minimum(%(0)s, %(1)s)", ("numpy",)),
+    "max": ("np.maximum(%(0)s, %(1)s)", ("numpy",)),
+    "exp": ("np.exp(%(0)s)", ("numpy",)),
+    "sin": ("np.sin(%(0)s)", ("numpy",)),
+    "cos": ("np.cos(%(0)s)", ("numpy",)),
+    "tan": ("np.tan(%(0)s)", ("numpy",)),
+    "arcsin": ("np.arcsin(%(0)s)", ("numpy",)),
+    "arccos": ("np.arccos(%(0)s)", ("numpy",)),
+    "arctan": ("np.arctan(%(0)s)", ("numpy",)),
+    "sinh": ("np.sinh(%(0)s)", ("numpy",)),
+    "cosh": ("np.cosh(%(0)s)", ("numpy",)),
+    "tanh": ("np.tanh(%(0)s)", ("numpy",)),
+    "sqrt": ("np.sqrt(%(0)s)", ("numpy",)),
+    "ln": ("np.log(%(0)s)", ("numpy",)),
+    "log": ("(np.log(%(0)s)/np.log(%(1)s))", ("numpy",)),
+    # NUMPY: "invert_matrix": ("np.linalg.inv(%(0)s)", ("numpy",)),
+
+    # vector functions with axis to apply over
+    # NUMPY:
+    # "prod": "np.prod(%(0)s, axis=%(axis)s)", ("numpy",)),
+    # "sum": "np.sum(%(0)s, axis=%(axis)s)", ("numpy",)),
+    # "vmax": "np.max(%(0)s, axis=%(axis)s)", ("numpy", )),
+    # "vmin": "np.min(%(0)s, axis=%(axis)s)", ("numpy",))
+    "prod": ("prod(%(0)s, dim=%(axis)s)", ("functions", "prod")),
+    "sum": ("sum(%(0)s, dim=%(axis)s)", ("functions", "sum")),
+    "vmax": ("vmax(%(0)s, dim=%(axis)s)", ("functions", "vmax")),
+    "vmin": ("vmin(%(0)s, dim=%(axis)s)", ("functions", "vmin")),
+
+    # functions defined in pysd.py_bakcend.functions
+    "active_initial": (  # TODO replace time by stage when doing a non compatible version
+        "active_initial(__data['time'], lambda: %(0)s, %(1)s)",
+        ("functions", "active_initial")),
+    "if_then_else": (
+        "if_then_else(%(0)s, lambda: %(1)s, lambda: %(2)s)",
+        ("functions", "if_then_else")),
+    "integer": (
+        "integer(%(0)s)",
+        ("functions", "integer")),
+    "invert_matrix": (  # NUMPY: remove
+        "invert_matrix(%(0)s)",
+        ("functions", "invert_matrix")),  # NUMPY: remove
+    "modulo": (
+        "modulo(%(0)s, %(1)s)",
+        ("functions", "modulo")),
+    "pulse": (
+        "pulse(__data['time'], %(0)s, %(1)s)",
+        ("functions", "pulse")),
+    "pulse_train": (
+        "pulse_train(__data['time'], %(0)s, %(1)s, %(2)s, %(3)s)",
+        ("functions", "pulse_train")),
+    "quantum": (
+        "quantum(%(0)s, %(1)s)",
+        ("functions", "quantum")),
+    "ramp": (
+        "ramp(__data['time'], %(0)s, %(1)s, %(2)s)",
+        ("functions", "ramp")),
+    "step": (
+        "step(__data['time'], %(0)s, %(1)s)",
+        ("functions", "step")),
+    "xidz": (
+        "xidz(%(0)s, %(1)s, %(2)s)",
+        ("functions", "xidz")),
+    "zidz": (
+        "zidz(%(0)s, %(1)s)",
+        ("functions", "zidz")),
+
+    # random functions must have the shape of the component subscripts
+    # most of them are shifted, scaled and truncated
+    # TODO: it is difficult to find same parametrization in python,
+    # maybe build a new model
+    "random_0_1": (
+        "np.random.uniform(0, 1, size=%(size)s)",
+        ("numpy",)),
+    "random_uniform": (
+        "np.random.uniform(%(0)s, %(1)s, size=%(size)s)",
+        ("numpy",)),
+    "random_normal": (
+        "stats.truncnorm.rvs(%(0)s, %(1)s, loc=%(2)s, scale=%(3)s, size=%(size)s))",
+        ("scipy", "stats")),
+}
diff --git a/pysd/building/python/python_utils.py b/pysd/building/python/python_utils.py
new file mode 100644
index 00000000..1836bdac
--- /dev/null
+++ b/pysd/building/python/python_utils.py
@@ -0,0 +1,61 @@
+import re
+import warnings
+import numpy as np
+
+# used to create python safe names with the variable reserved_words
+from keyword import kwlist
+from builtins import __dir__ as bidir
+from pysd.py_backend.components import __dir__ as cdir
+from pysd.py_backend.data import __dir__ as ddir
+from pysd.py_backend.decorators import __dir__ as dedir
+from pysd.py_backend.external import __dir__ as edir
+from pysd.py_backend.functions import __dir__ as fdir
+from pysd.py_backend.statefuls import __dir__ as sdir
+from pysd.py_backend.utils import __dir__ as udir
+
+
+reserved_words = set(
+    dir() + bidir() + cdir() + ddir() + dedir() + edir() + fdir()
+    + sdir() + udir()).union(kwlist)
+
+
+def simplify_subscript_input(coords, subscript_dict, return_full, merge_subs):
+    """
+    Parameters
+    ----------
+    coords: dict
+        Coordinates to write in the model file.
+
+    subscript_dict: dict
+        The subscript dictionary of the model file.
+
+    return_full: bool
+        If True the when coords == subscript_dict, '_subscript_dict'
+        will be returned
+
+    merge_subs: list of strings
+        List of the final subscript range of the python array after
+        merging with other objects
+
+    Returns
+    -------
+    coords: str
+        The equations to generate the coord dicttionary in the model file.
+
+    """
+
+    if coords == subscript_dict and return_full:
+        # variable defined with all the subscripts
+        return "_subscript_dict"
+
+    coordsp = []
+    for ndim, (dim, coord) in zip(merge_subs, coords.items()):
+        # find dimensions can be retrieved from _subscript_dict
+        if coord == subscript_dict[dim]:
+            # use _subscript_dict
+            coordsp.append(f"'{ndim}': _subscript_dict['{dim}']")
+        else:
+            # write whole dict
+            coordsp.append(f"'{ndim}': {coord}")
+
+    return "{" + ", ".join(coordsp) + "}"
diff --git a/pysd/building/python/subscripts.py b/pysd/building/python/subscripts.py
new file mode 100644
index 00000000..d004f0a9
--- /dev/null
+++ b/pysd/building/python/subscripts.py
@@ -0,0 +1,350 @@
+from multiprocessing.sharedctypes import Value
+import warnings
+from pathlib import Path
+import numpy as np
+from pysd.translation.structures.abstract_model import AbstractSubscriptRange
+from pysd.py_backend.external import ExtSubscript
+from typing import List
+
+
+class SubscriptManager:
+    def __init__(self, abstract_subscripts: List[AbstractSubscriptRange],
+                 _root: Path):
+        self._root = _root
+        self._copied = []
+        self.mapping = {}
+        self.subscripts = abstract_subscripts
+        self.elements = {}
+        self.subranges = self._get_main_subscripts()
+        self.subscript2num = self._get_subscript2num()
+        # TODO: manage subscript mapping
+
+    @property
+    def subscripts(self):
+        return self._subscripts
+
+    @subscripts.setter
+    def subscripts(self, abstract_subscripts):
+        self._subscripts = {}
+        missing = []
+        for sub in abstract_subscripts:
+            self.mapping[sub.name] = sub.mapping
+            if isinstance(sub.subscripts, list):
+                # regular definition of subscripts
+                self._subscripts[sub.name] = sub.subscripts
+            elif isinstance(sub.subscripts, str):
+                # copied subscripts, this will be always a subrange,
+                # then we need to prevent them of being saved as a main range
+                self._copied.append(sub.name)
+                self.mapping[sub.name].append(sub.subscripts)
+                if sub.subscripts in self._subscripts:
+                    self._subscripts[sub.name] =\
+                        self._subscripts[sub.subscripts]
+                else:
+                    missing.append(sub)
+            elif isinstance(sub.subscripts, dict):
+                # subscript from file
+                self._subscripts[sub.name] = ExtSubscript(
+                    file_name=sub.subscripts["file"],
+                    sheet=sub.subscripts["tab"],
+                    firstcell=sub.subscripts["firstcell"],
+                    lastcell=sub.subscripts["lastcell"],
+                    prefix=sub.subscripts["prefix"],
+                    root=self._root).subscript
+            else:
+                raise ValueError(
+                    f"Invalid definition of subscript {sub.name}:\n\t"
+                    + str(sub.subscripts))
+
+        while missing:
+            # second loop for copied subscripts
+            sub = missing.pop()
+            self._subscripts[sub.name] =\
+                self._subscripts[sub.subscripts]
+
+    def _get_main_subscripts(self):
+        """
+        Reutrns a dictionary with the main ranges as keys and their
+        subranges as values.
+        """
+        subscript_sets = {
+            name: set(subs) for name, subs in self.subscripts.items()}
+
+        subranges = {}
+        for range, subs in subscript_sets.items():
+            # current subscript range
+            subranges[range] = []
+            for subrange, subs2 in subscript_sets.items():
+                if range == subrange:
+                    # pass current range
+                    continue
+                elif subs == subs2:
+                    # range is equal to the subrange, as Vensim does
+                    # the main range will be the first one alphabetically
+                    # make it case insensitive
+                    range_l = range.replace(" ", "_").lower()
+                    subrange_l = subrange.replace(" ", "_").lower()
+                    if range_l < subrange_l and range not in self._copied:
+                        subranges[range].append(subrange)
+                    else:
+                        # copied subscripts ranges or subscripts ranges
+                        # that come later alphabetically
+                        del subranges[range]
+                        break
+                elif subs2.issubset(subs):
+                    # subrange is a subset of range, append it to the list
+                    subranges[range].append(subrange)
+                elif subs2.issuperset(subs):
+                    # it exist a range that contents the elements of the range
+                    del subranges[range]
+                    break
+
+        return subranges
+
+    def _get_subscript2num(self):
+        """
+        Build a dictionary to return the numeric value or values of a
+        subscript or subscript range.
+        """
+        s2n = {}
+        for range, subranges in self.subranges.items():
+            # a main range is direct to return
+            s2n[range.replace(" ", "_").lower()] = (
+                f"np.arange(1, len(_subscript_dict['{range}'])+1)",
+                {range: self.subscripts[range]}
+            )
+            for i, sub in enumerate(self.subscripts[range], start=1):
+                # a subscript must return its numeric position
+                # in the main range
+                s2n[sub.replace(" ", "_").lower()] = (str(i), {})
+            for subrange in subranges:
+                # subranges may return the position of each subscript
+                # in the main range
+                sub_index = [
+                    self.subscripts[range].index(sub)+1
+                    for sub in self.subscripts[subrange]]
+
+                if np.all(
+                  sub_index
+                  == np.arange(sub_index[0], sub_index[0]+len(sub_index))):
+                    # subrange definition can be simplified with a range
+                    subsarray = f"np.arange({sub_index[0]}, "\
+                        f"len(_subscript_dict['{subrange}'])+{sub_index[0]})"
+                else:
+                    # subrange definition cannot be simplified
+                    subsarray = f"np.array({sub_index})"
+
+                s2n[subrange.replace(" ", "_").lower()] = (
+                    subsarray,
+                    {subrange: self.subscripts[subrange]}
+                )
+
+        return s2n
+
+    def find_subscript_name(self, element, avoid=[]):
+        """
+        Given a subscript dictionary, and a member of a subscript family,
+        return the first key of which the member is within the value list.
+        If element is already a subscript name, return that.
+
+        Parameters
+        ----------
+        element: str
+            Subscript or subscriptrange name to find.
+        avoid: list (optional)
+            List of subscripts to avoid. Default is an empty list.
+
+        Returns
+        -------
+
+        Examples
+        --------
+        >>> find_subscript_name('D')
+        'Dim2'
+        >>> find_subscript_name('B')
+        'Dim1'
+        >>> find_subscript_name('B', avoid=['Dim1'])
+        'Dim2'
+
+        """
+        if element in self.subscripts.keys():
+            return element
+
+        for name, elements in self.subscripts.items():
+            if element in elements and name not in avoid:
+                return name
+
+    def make_coord_dict(self, subs):
+        """
+        This is for assisting with the lookup of a particular element.
+
+        Parameters
+        ----------
+        subs: list of strings
+            Coordinates, either as names of dimensions, or positions within
+            a dimension.
+
+        Returns
+        -------
+        coordinates: dict
+            Coordinates needed to access the xarray quantities we are
+            interested in.
+
+        Examples
+        --------
+        >>> make_coord_dict(['Dim1', 'D'])
+        {'Dim1': ['A', 'B', 'C'], 'Dim2': ['D']}
+
+        """
+        sub_elems_list = [y for x in self.subscripts.values() for y in x]
+        coordinates = {}
+        for sub in subs:
+            if sub in sub_elems_list:
+                name = self.find_subscript_name(
+                    sub, avoid=subs + list(coordinates))
+                coordinates[name] = [sub]
+            else:
+                if sub.endswith("!"):
+                    coordinates[sub] = self.subscripts[sub[:-1]]
+                else:
+                    coordinates[sub] = self.subscripts[sub]
+        return coordinates
+
+    def make_merge_list(self, subs_list, element=""):
+        """
+        This is for assisting when building xrmerge. From a list of subscript
+        lists returns the final subscript list after mergin. Necessary when
+        merging variables with subscripts comming from different definitions.
+
+        Parameters
+        ----------
+        subs_list: list of lists of strings
+            Coordinates, either as names of dimensions, or positions within
+            a dimension.
+        element: str (optional)
+            Element name, if given it will be printed with any error or
+            warning message. Default is "".
+
+        Returns
+        -------
+        dims: list
+            Final subscripts after merging.
+
+        Examples
+        --------
+        >>> sm = SubscriptManager()
+        >>> sm.subscripts = {"upper": ["A", "B"], "all": ["A", "B", "C"]}
+        >>> sm.make_merge_list([['upper'], ['C']])
+        ['all']
+
+        """
+        coords_set = [set() for i in range(len(subs_list[0]))]
+        coords_list = [
+            self.make_coord_dict(subs)
+            for subs in subs_list
+        ]
+
+        # update coords set
+        [[coords_set[i].update(coords[dim]) for i, dim in enumerate(coords)]
+         for coords in coords_list]
+
+        dims = [None] * len(coords_set)
+        # create an array with the name of the subranges for all
+        # merging elements
+        dims_list = np.array([
+            list(coords) for coords in coords_list]).transpose()
+        indexes = np.arange(len(dims))
+
+        for i, coord2 in enumerate(coords_set):
+            dims1 = [
+                dim for dim in dims_list[i]
+                if dim is not None and set(self.subscripts[dim]) == coord2
+            ]
+            if dims1:
+                # if the given coordinate already matches return it
+                dims[i] = dims1[0]
+            else:
+                # find a suitable coordinate
+                other_dims = dims_list[indexes != i]
+                for name, elements in self.subscripts.items():
+                    if coord2 == set(elements) and name not in other_dims:
+                        dims[i] = name
+                        break
+
+                if not dims[i]:
+                    # the dimension is incomplete use the smaller
+                    # dimension that completes it
+                    for name, elements in self.subscripts.items():
+                        if coord2.issubset(set(elements))\
+                           and name not in other_dims:
+                            dims[i] = name
+                            warnings.warn(
+                                element
+                                + "\nDimension given by subscripts:"
+                                + "\n\t{}\nis incomplete ".format(coord2)
+                                + "using {} instead.".format(name)
+                                + "\nSubscript_dict:"
+                                + "\n\t{}".format(self.subscripts)
+                            )
+                            break
+
+                if not dims[i]:
+                    for name, elements in self.subscripts.items():
+                        if coord2 == set(elements):
+                            j = 1
+                            while name + str(j) in self.subscripts.keys():
+                                j += 1
+                            self.subscripts[name + str(j)] = elements
+                            dims[i] = name + str(j)
+                            warnings.warn(
+                                element
+                                + "\nAdding new subscript range to"
+                                + " subscript_dict:\n"
+                                + name + str(j) + ": " + ', '.join(elements))
+                            break
+
+                if not dims[i]:
+                    # not able to find the correct dimension
+                    raise ValueError(
+                        element
+                        + "\nImpossible to find the dimension that contains:"
+                        + "\n\t{}\nFor subscript_dict:".format(coord2)
+                        + "\n\t{}".format(self.subscripts)
+                    )
+
+        return dims
+
+    def simplify_subscript_input(self, coords, merge_subs):
+        """
+        Parameters
+        ----------
+        coords: dict
+            Coordinates to write in the model file.
+
+        merge_subs: list of strings
+            List of the final subscript range of the python array after
+            merging with other objects
+
+        Returns
+        -------
+        final_subs, coords: dict, str
+            Final subscripts and the equations to generate the coord
+            dicttionary in the model file.
+
+        """
+        coordsp = []
+        final_subs = {}
+        for ndim, (dim, coord) in zip(merge_subs, coords.items()):
+            # find dimensions can be retrieved from _subscript_dict
+            final_subs[ndim] = coord
+            if dim.endswith("!") and coord == self.subscripts[dim[:-1]]:
+                # use _subscript_dict
+                coordsp.append(f"'{ndim}': _subscript_dict['{dim[:-1]}']")
+            elif not dim.endswith("!") and coord == self.subscripts[dim]:
+                # use _subscript_dict
+                coordsp.append(f"'{ndim}': _subscript_dict['{dim}']")
+            else:
+                # write whole dict
+                coordsp.append(f"'{ndim}': {coord}")
+
+        return final_subs, "{" + ", ".join(coordsp) + "}"
\ No newline at end of file
diff --git a/pysd/building/python/visitors.py b/pysd/building/python/visitors.py
new file mode 100644
index 00000000..3034ad18
--- /dev/null
+++ b/pysd/building/python/visitors.py
@@ -0,0 +1,1257 @@
+from re import X
+import warnings
+from dataclasses import dataclass
+
+import numpy as np
+from pysd.py_backend.utils import compute_shape
+
+from pysd.translation.structures import components as ct
+from .python_functions import functionspace
+
+
+@dataclass
+class BuildAST:
+    expression: str
+    calls: dict
+    subscripts: dict
+    order: int
+
+    def __str__(self):
+        # makes easier building
+        return self.expression
+
+    def reshape(self, subscripts, final_subscripts):
+        subscripts_out = subscripts.simplify_subscript_input(
+                final_subscripts, list(final_subscripts))[1]
+        if not final_subscripts or (
+          self.subscripts == final_subscripts
+          and list(self.subscripts) == list(final_subscripts)):
+            # same dictionary in the same orde, do nothing
+            pass
+        elif not self.subscripts:
+            # original expression is not an array
+            # NUMPY: object.expression = np.full(%s, %(shape)s)
+            self.expression = "xr.DataArray(%s, %s, %s)" % (
+                self.expression, subscripts_out, list(final_subscripts)
+            )
+            self.order = 0
+        else:
+            # original expression is an array
+            # NUMPY: reorder dims if neccessary with np.moveaxis or similar
+            # NUMPY: add new axis with [:, None, :] or np.tile,
+            #        depending on an input argument
+            # NUMPY: if order is not 0 need to lower the order to 0
+            # using force!
+            self.expression = "(xr.DataArray(0, %s, %s) + %s)" % (
+                subscripts_out, list(final_subscripts), self.expression
+                )
+            self.order = 0
+        self.subscripts = final_subscripts
+
+    def lower_order(self, new_order, force_0=False):
+        if self.order >= new_order and self.order != 0\
+          and (new_order != 0 or force_0):
+            # if current operator order is 0 do not need to do anything
+            # if the order of operations conflicts add parenthesis
+            # if new order is 0 do not need to do anything, as it may be
+            # an argument to a function, unless force_0 is True which
+            # will force the parenthesis (necessary to reshape some
+            # numpy arrays)
+            self.expression = "(%s)" % self.expression
+            self.order = 0
+
+
+class StructureBuilder:
+    def __init__(self, value, component):
+        self.value = value
+        self.arguments = {}
+        self.component = component
+        self.element = component.element
+        self.section = component.section
+        self.def_subs = component.subscripts_dict
+
+    def build(self, arguments):
+        return BuildAST(
+            expression=repr(self.value),
+            calls={},
+            subscripts={},
+            order=0)
+
+    def join_calls(self, arguments):
+        if len(arguments) == 0:
+            return {}
+        elif len(arguments) == 1:
+            return arguments["0"].calls
+        else:
+            return merge_dependencies(
+                *[val.calls for val in arguments.values()])
+
+    def reorder(self, arguments, def_subs=None, force=None):
+
+        if force == "component":
+            final_subscripts = def_subs or {}
+        else:
+            final_subscripts = self.get_final_subscripts(
+                arguments, def_subs)
+
+        [arguments[key].reshape(self.section.subscripts, final_subscripts)
+         for key in arguments
+         if arguments[key].subscripts or force == "equal"]
+
+        return final_subscripts
+
+    def get_final_subscripts(self, arguments, def_subs):
+        if len(arguments) == 0:
+            return {}
+        elif len(arguments) == 1:
+            return arguments["0"].subscripts
+        else:
+            return self._compute_final_subscripts(
+                [arg.subscripts for arg in arguments.values()],
+                def_subs)
+
+    def _compute_final_subscripts(self, subscripts_list, def_subs):
+        expression = {}
+        [expression.update(subscript)
+         for subscript in subscripts_list if subscript]
+        # TODO reorder final_subscripts taking into account def_subs
+        return expression
+
+
+class OperationBuilder(StructureBuilder):
+    operators_build = {
+        "^": ("%(left)s**%(right)s", None, 1),
+        "*": ("%(left)s*%(right)s", None, 2),
+        "/": ("%(left)s/%(right)s", None, 2),
+        "+": ("%(left)s + %(right)s", None, 3),
+        "-": ("%(left)s - %(right)s", None, 3),
+        "=": ("%(left)s == %(right)s", None, 4),
+        "<>": ("%(left)s != %(right)s", None, 4),
+        ">=": ("%(left)s >= %(right)s", None, 4),
+        ">": ("%(left)s > %(right)s", None, 4),
+        "<=": ("%(left)s <= %(right)s", None, 4),
+        "<": ("%(left)s < %(right)s", None, 4),
+        ":NOT:": ("np.logical_not(%s)", ("numpy",), 0),
+        ":AND:": ("np.logical_and(%(left)s, %(right)s)", ("numpy",), 0),
+        ":OR:": ("np.logical_or(%(left)s, %(right)s)", ("numpy",), 0),
+        "negative": ("-%s", None, 3),
+    }
+
+    def __init__(self, operation, component):
+        super().__init__(None, component)
+        self.operators = operation.operators.copy()
+        self.arguments = {
+            str(i): arg for i, arg in enumerate(operation.arguments)}
+
+    def build(self, arguments):
+        operands = {}
+        calls = self.join_calls(arguments)
+        final_subscripts = self.reorder(arguments, def_subs=self.def_subs)
+        arguments = [arguments[str(i)] for i in range(len(arguments))]
+        dependencies, order = self.operators_build[self.operators[-1]][1:]
+
+        if dependencies:
+            self.section.imports.add(*dependencies)
+
+        if self.operators[-1] == "^":
+            # right side of the exponential can be from higher order
+            arguments[-1].lower_order(2)
+        else:
+            arguments[-1].lower_order(order)
+
+        if len(arguments) == 1:
+            # not and negative operations (only 1 element)
+            if self.operators[0] == "negative":
+                order = 1
+            expression = self.operators_build[self.operators[0]][0]
+            return BuildAST(
+                expression=expression % arguments[0],
+                calls=calls,
+                subscripts=final_subscripts,
+                order=order)
+
+        operands["right"] = arguments.pop()
+        while arguments or self.operators:
+            expression = self.operators_build[self.operators.pop()][0]
+            operands["left"] = arguments.pop()
+            operands["left"].lower_order(order)
+            operands["right"] = expression % operands
+
+        return BuildAST(
+            expression=operands["right"],
+            calls=calls,
+            subscripts=final_subscripts,
+            order=order)
+
+
+class GameBuilder(StructureBuilder):
+    def __init__(self, game_str, component):
+        super().__init__(None, component)
+        self.arguments = {"expr": game_str.expression}
+
+    def build(self, arguments):
+        return arguments["expr"]
+
+
+class CallBuilder(StructureBuilder):
+    def __init__(self, call_str, component):
+        super().__init__(None, component)
+        function_name = call_str.function.reference
+        self.arguments = {
+            str(i): arg for i, arg in enumerate(call_str.arguments)}
+
+        # move this to a setter
+        if function_name in self.section.macrospace:
+            # build macro
+            self.macro_name = function_name
+            self.build = self.build_macro_call
+        elif function_name in self.section.namespace.cleanspace:
+            # build lookupcall
+            self.arguments["function"] = call_str.function
+            self.build = self.build_lookups_call
+        elif function_name in functionspace:
+            # build direct function
+            self.function = function_name
+            self.build = self.build_function_call
+        elif function_name == "a_function_of":
+            self.build = self.build_incomplete_call
+        else:
+            # error
+            raise ValueError("Undefined function %s" % function_name)
+
+    def build_macro_call(self, arguments):
+        self.section.imports.add("statefuls", "Macro")
+        macro = self.section.macrospace[self.macro_name]
+
+        calls = self.join_calls(arguments)
+        final_subscripts = self.reorder(arguments, def_subs=self.def_subs)
+
+        arguments["name"] = self.section.namespace.make_python_identifier(
+            self.macro_name + "_" + self.element.identifier, prefix="_macro")
+        arguments["file"] = macro.path.name
+        arguments["macro_name"] = macro.name
+        arguments["args"] = "{%s}" % ", ".join([
+            "'%s': lambda: %s" % (key, val)
+            for key, val in zip(macro.params, arguments.values())
+        ])
+
+        self.element.objects[arguments["name"]] = {
+            "name": arguments["name"],
+            "expression": "%(name)s = Macro(_root.joinpath('%(file)s'), "
+                          "%(args)s, '%(macro_name)s', "
+                          "time_initialization=lambda: __data['time'], "
+                          "py_name='%(name)s')" % arguments,
+            "calls": {
+                "initial": calls,
+                "step": calls
+            }
+        }
+        return BuildAST(
+            expression="%s()" % arguments["name"],
+            calls={arguments["name"]: 1},
+            subscripts=final_subscripts,
+            order=0)
+
+    def build_incomplete_call(self, arguments):
+        warnings.warn(
+            "%s has no equation specified" % self.element.name,
+            SyntaxWarning, stacklevel=2
+        )
+        self.section.imports.add("functions", "incomplete")
+        return BuildAST(
+            expression="incomplete(%s)" % ", ".join(
+                arg.expression for arg in arguments.values()),
+            calls=self.join_calls(arguments),
+            subscripts=self.def_subs,
+            order=0)
+
+    def build_lookups_call(self, arguments):
+        expression = arguments["function"].expression.replace("()", "(%(0)s)")
+        final_subscripts = self.get_final_subscripts(arguments, self.def_subs)
+        # NUMPY: we need to manage inside lookup with subscript and later
+        #        return the values in a correct ndarray
+        return BuildAST(
+            expression=expression % arguments,
+            calls=self.join_calls(arguments),
+            subscripts=final_subscripts,
+            order=0)
+
+    def build_function_call(self, arguments):
+        expression, modules = functionspace[self.function]
+        if modules:
+            self.section.imports.add(*modules)
+
+        calls = self.join_calls(arguments)
+
+        if "__data['time']" in expression:
+            merge_dependencies(calls, {"time": 1}, inplace=True)
+
+        # TODO modify dimensions of BuildAST
+        if "%(axis)s" in expression:
+            final_subscripts, arguments["axis"] = self.compute_axis(arguments)
+
+        elif "%(size)s" in expression:
+            final_subscripts = self.reorder(
+                arguments,
+                def_subs=self.def_subs,
+                force="component"
+            )
+            arguments["size"] = compute_shape(final_subscripts)
+
+        elif self.function == "active_initial":
+            # we need to ensure that active initial outputs are always the
+            # same and update dependencies as stateful object
+            # TODO: update calls as statefull object
+            name = self.section.namespace.make_python_identifier(
+                self.element.identifier, prefix="_active_initial")
+            final_subscripts = self.reorder(
+                arguments,
+                def_subs=self.def_subs,
+                force="equal"
+            )
+            self.element.objects[name] = {
+                "name": name,
+                "expression": None,
+                "calls": {
+                    "initial": arguments["1"].calls,
+                    "step": arguments["0"].calls
+                }
+
+            }
+            calls = {name: 1}
+        else:
+            final_subscripts = self.reorder(
+                arguments,
+                def_subs=self.def_subs
+            )
+            if self.function == "xidz" and final_subscripts:
+                if not arguments["1"].subscripts:
+                    new_args = {"0": arguments["0"], "2": arguments["2"]}
+                    self.reorder(
+                        new_args,
+                        def_subs=self.def_subs,
+                        force="equal"
+                    )
+                    arguments.update(new_args)
+            if self.function == "if_then_else" and final_subscripts:
+                if not arguments["0"].subscripts:
+                    # NUMPY: we need to ensure that if_then_else always returs
+                    # the same shape object
+                    new_args = {"1": arguments["1"], "2": arguments["2"]}
+                    self.reorder(
+                        new_args,
+                        def_subs=self.def_subs,
+                        force="equal"
+                    )
+                    arguments.update(new_args)
+                else:
+                    self.reorder(
+                        arguments,
+                        def_subs=self.def_subs,
+                        force="equal"
+                    )
+
+        return BuildAST(
+            expression=expression % arguments,
+            calls=calls,
+            subscripts=final_subscripts,
+            order=0)
+
+    def compute_axis(self, arguments):
+        subscripts = arguments["0"].subscripts
+        axis = []
+        coords = {}
+        for subs in subscripts:
+            if subs.endswith("!"):
+                # dimensions to apply along
+                axis.append(subs)
+            else:
+                # dimensions remaining
+                coords[subs] = subscripts[subs]
+        return coords, axis
+
+
+class ExtLookupBuilder(StructureBuilder):
+    def __init__(self, getlookup_str,  component):
+        super().__init__(None, component)
+        self.file = getlookup_str.file
+        self.tab = getlookup_str.tab
+        self.x_row_or_col = getlookup_str.x_row_or_col
+        self.cell = getlookup_str.cell
+        self.arguments = {}
+
+    def build(self, arguments):
+        self.component.type = "Lookup"
+        self.component.subtype = "External"
+        arguments["params"] = "'%s', '%s', '%s', '%s'" % (
+            self.file, self.tab, self.x_row_or_col, self.cell
+        )
+        final_subs, arguments["subscripts"] =\
+            self.section.subscripts.simplify_subscript_input(
+                self.def_subs, self.element.subscripts)
+
+        if "ext_lookups" in self.element.objects:
+            # object already exists
+            self.element.objects["ext_lookups"]["expression"] += "\n\n"\
+                + self.element.objects["ext_lookups"]["name"]\
+                + ".add(%(params)s, %(subscripts)s)" % arguments
+
+            return None
+        else:
+            # create a new object
+            self.section.imports.add("external", "ExtLookup")
+
+            arguments["name"] = self.section.namespace.make_python_identifier(
+                self.element.identifier, prefix="_ext_lookup")
+
+            self.element.objects["ext_lookups"] = {
+                "name": arguments["name"],
+                "expression": "%(name)s = ExtLookup(%(params)s, "
+                              "%(subscripts)s, "
+                              "_root, '%(name)s')" % arguments
+            }
+
+            return BuildAST(
+                expression=arguments["name"] + "(x)",
+                calls={"__external__": None, "__lookup__": None},
+                subscripts=final_subs,
+                order=0)
+
+class ExtDataBuilder(StructureBuilder):
+    def __init__(self, getdata_str,  component):
+        super().__init__(None, component)
+        self.file = getdata_str.file
+        self.tab = getdata_str.tab
+        self.time_row_or_col = getdata_str.time_row_or_col
+        self.cell = getdata_str.cell
+        self.keyword = component.keyword
+        self.arguments = {}
+
+    def build(self, arguments):
+        self.component.type = "Data"
+        self.component.subtype = "External"
+        arguments["params"] = "'%s', '%s', '%s', '%s'" % (
+            self.file, self.tab, self.time_row_or_col, self.cell
+        )
+        final_subs, arguments["subscripts"] =\
+            self.section.subscripts.simplify_subscript_input(
+                self.def_subs, self.element.subscripts)
+        arguments["method"] = "'%s'" % self.keyword if self.keyword else None
+
+        if "ext_data" in self.element.objects:
+            # object already exists
+            self.element.objects["ext_data"]["expression"] += "\n\n"\
+                + self.element.objects["ext_data"]["name"]\
+                + ".add(%(params)s, %(method)s, %(subscripts)s)" % arguments
+
+            return None
+        else:
+            # create a new object
+            self.section.imports.add("external", "ExtData")
+
+            arguments["name"] = self.section.namespace.make_python_identifier(
+                self.element.identifier, prefix="_ext_data")
+
+            self.element.objects["ext_data"] = {
+                "name": arguments["name"],
+                "expression": "%(name)s = ExtData(%(params)s, "
+                              " %(method)s, %(subscripts)s, "
+                              "_root, '%(name)s')" % arguments
+            }
+
+            return BuildAST(
+                expression=arguments["name"] + "(time())",
+                calls={"__external__": None, "time": 1},
+                subscripts=final_subs,
+                order=0)
+
+
+class ExtConstantBuilder(StructureBuilder):
+    def __init__(self, getlookup_str,  component):
+        super().__init__(None, component)
+        self.file = getlookup_str.file
+        self.tab = getlookup_str.tab
+        self.cell = getlookup_str.cell
+        self.arguments = {}
+
+    def build(self, arguments):
+        self.component.type = "Constant"
+        self.component.subtype = "External"
+        arguments["params"] = "'%s', '%s', '%s'" % (
+            self.file, self.tab, self.cell
+        )
+        final_subs, arguments["subscripts"] =\
+            self.section.subscripts.simplify_subscript_input(
+                self.def_subs, self.element.subscripts)
+
+        if "constants" in self.element.objects:
+            # object already exists
+            self.element.objects["constants"]["expression"] += "\n\n"\
+                + self.element.objects["constants"]["name"]\
+                + ".add(%(params)s, %(subscripts)s)" % arguments
+
+            return None
+        else:
+            # create a new object
+            self.section.imports.add("external", "ExtConstant")
+
+            arguments["name"] = self.section.namespace.make_python_identifier(
+                self.element.identifier, prefix="_ext_constant")
+
+            self.element.objects["constants"] = {
+                "name": arguments["name"],
+                "expression": "%(name)s = ExtConstant(%(params)s, "
+                              "%(subscripts)s, _root, '%(name)s')" % arguments
+            }
+
+            return BuildAST(
+                expression=arguments["name"] + "()",
+                calls={"__external__": None},
+                subscripts=final_subs,
+                order=0)
+
+
+class TabDataBuilder(StructureBuilder):
+    def __init__(self, data_str,  component):
+        super().__init__(None, component)
+        self.keyword = component.keyword
+        self.arguments = {}
+
+    def build(self, arguments):
+        self.section.imports.add("data", "TabData")
+
+        final_subs, arguments["subscripts"] =\
+            self.section.subscripts.simplify_subscript_input(
+                self.def_subs, self.element.subscripts)
+
+        arguments["real_name"] = self.element.name
+        arguments["py_name"] =\
+            self.section.namespace.namespace[self.element.name]
+        arguments["subscripts"] = self.def_subs
+        arguments["method"] = "'%s'" % self.keyword if self.keyword else None
+
+        arguments["name"] = self.section.namespace.make_python_identifier(
+            self.element.identifier, prefix="_data")
+
+        self.element.objects["tab_data"] = {
+            "name": arguments["name"],
+            "expression": "%(name)s = TabData('%(real_name)s', '%(py_name)s', "
+                          "%(subscripts)s,  %(method)s)" % arguments
+        }
+
+        return BuildAST(
+            expression=arguments["name"] + "(time())",
+            calls={"time": 1, "__data__": None},
+            subscripts=final_subs,
+            order=0)
+
+
+class InitialBuilder(StructureBuilder):
+    def __init__(self, initial_str, component):
+        super().__init__(None, component)
+        self.arguments = {
+            "initial": initial_str.initial
+        }
+
+    def build(self, arguments):
+        self.component.type = "Stateful"
+        self.component.subtype = "Initial"
+        self.section.imports.add("statefuls", "Initial")
+        arguments["initial"].reshape(self.section.subscripts, self.def_subs)
+        arguments["name"] = self.section.namespace.make_python_identifier(
+            self.element.identifier, prefix="_initial")
+
+        self.element.objects[arguments["name"]] = {
+            "name": arguments["name"],
+            "expression": "%(name)s = Initial(lambda: %(initial)s, "
+                          "'%(name)s')" % arguments,
+            "calls": {
+                "initial": arguments["initial"].calls,
+                "step": {}
+            }
+
+        }
+        return BuildAST(
+            expression=arguments["name"] + "()",
+            calls={arguments["name"]: 1},
+            subscripts=self.def_subs,
+            order=0)
+
+
+class IntegBuilder(StructureBuilder):
+    def __init__(self, integ_str, component):
+        super().__init__(None, component)
+        self.arguments = {
+            "flow": integ_str.flow,
+            "initial": integ_str.initial
+        }
+
+    def build(self, arguments):
+        self.component.type = "Stateful"
+        self.component.subtype = "Integ"
+        self.section.imports.add("statefuls", "Integ")
+        arguments["initial"].reshape(self.section.subscripts, self.def_subs)
+        arguments["flow"].reshape(self.section.subscripts, self.def_subs)
+        arguments["name"] = self.section.namespace.make_python_identifier(
+            self.element.identifier, prefix="_integ")
+
+        self.element.objects[arguments["name"]] = {
+            "name": arguments["name"],
+            "expression": "%(name)s = Integ(lambda: %(flow)s, "
+                          "lambda: %(initial)s, '%(name)s')" % arguments,
+            "calls": {
+                "initial": arguments["initial"].calls,
+                "step": arguments["flow"].calls
+            }
+
+        }
+        return BuildAST(
+            expression=arguments["name"] + "()",
+            calls={arguments["name"]: 1},
+            subscripts=self.def_subs,
+            order=0)
+
+
+class DelayBuilder(StructureBuilder):
+    def __init__(self, dtype, delay_str,  component):
+        super().__init__(None, component)
+        self.arguments = {
+            "input": delay_str.input,
+            "delay_time": delay_str.delay_time,
+            "initial": delay_str.initial,
+            "order": delay_str.order
+        }
+        self.dtype = dtype
+
+    def build(self, arguments):
+        self.component.type = "Stateful"
+        self.component.subtype = "Delay"
+        self.section.imports.add("statefuls", self.dtype)
+        arguments["input"].reshape(self.section.subscripts, self.def_subs)
+        arguments["delay_time"].reshape(self.section.subscripts, self.def_subs)
+        arguments["initial"].reshape(self.section.subscripts, self.def_subs)
+        arguments["name"] = self.section.namespace.make_python_identifier(
+            self.element.identifier, prefix=f"_{self.dtype.lower()}")
+        arguments["dtype"] = self.dtype
+
+        self.element.objects[arguments["name"]] = {
+            "name": arguments["name"],
+            "expression": "%(name)s = %(dtype)s(lambda: %(input)s, "
+                          "lambda: %(delay_time)s, lambda: %(initial)s, "
+                          "lambda: %(order)s, "
+                          "time_step, '%(name)s')" % arguments,
+            "calls": {
+                "initial": merge_dependencies(
+                    arguments["initial"].calls,
+                    arguments["delay_time"].calls,
+                    arguments["order"].calls),
+                "step": merge_dependencies(
+                    arguments["input"].calls,
+                    arguments["delay_time"].calls)
+
+            }
+        }
+        return BuildAST(
+            expression=arguments["name"] + "()",
+            calls={arguments["name"]: 1},
+            subscripts=self.def_subs,
+            order=0)
+
+
+class DelayFixedBuilder(StructureBuilder):
+    def __init__(self, delay_str,  component):
+        super().__init__(None, component)
+        self.arguments = {
+            "input": delay_str.input,
+            "delay_time": delay_str.delay_time,
+            "initial": delay_str.initial,
+        }
+
+    def build(self, arguments):
+        self.component.type = "Stateful"
+        self.component.subtype = "DelayFixed"
+        self.section.imports.add("statefuls", "DelayFixed")
+        arguments["input"].reshape(self.section.subscripts, self.def_subs)
+        arguments["initial"].reshape(self.section.subscripts, self.def_subs)
+        arguments["name"] = self.section.namespace.make_python_identifier(
+            self.element.identifier, prefix="_delayfixed")
+
+        self.element.objects[arguments["name"]] = {
+            "name": arguments["name"],
+            "expression": "%(name)s = DelayFixed(lambda: %(input)s, "
+                          "lambda: %(delay_time)s, lambda: %(initial)s, "
+                          "time_step, '%(name)s')" % arguments,
+            "calls": {
+                "initial": merge_dependencies(
+                    arguments["initial"].calls,
+                    arguments["delay_time"].calls),
+                "step": arguments["input"].calls
+            }
+        }
+        return BuildAST(
+            expression=arguments["name"] + "()",
+            calls={arguments["name"]: 1},
+            subscripts=self.def_subs,
+            order=0)
+
+
+class SmoothBuilder(StructureBuilder):
+    def __init__(self, smooth_str,  component):
+        super().__init__(None, component)
+        self.arguments = {
+            "input": smooth_str.input,
+            "smooth_time": smooth_str.smooth_time,
+            "initial": smooth_str.initial,
+            "order": smooth_str.order
+        }
+
+    def build(self, arguments):
+        self.component.type = "Stateful"
+        self.component.subtype = "Smooth"
+        self.section.imports.add("statefuls", "Smooth")
+        arguments["input"].reshape(self.section.subscripts, self.def_subs)
+        arguments["smooth_time"].reshape(self.section.subscripts, self.def_subs)
+        arguments["initial"].reshape(self.section.subscripts, self.def_subs)
+        arguments["name"] = self.section.namespace.make_python_identifier(
+            self.element.identifier, prefix="_smooth")
+
+        # TODO in the future we need to ad timestep to show warnings about
+        # the smooth time as its done with delays (see vensim help for smooth)
+        # TODO in the future we may want to have 2 py_backend classes for
+        # smooth as the behaviour is different for SMOOTH and SMOOTH N when
+        # using RingeKutta scheme
+        self.element.objects[arguments["name"]] = {
+            "name": arguments["name"],
+            "expression": "%(name)s = Smooth(lambda: %(input)s, "
+                          "lambda: %(smooth_time)s, lambda: %(initial)s, "
+                          "lambda: %(order)s, '%(name)s')" % arguments,
+            "calls": {
+                "initial": merge_dependencies(
+                    arguments["initial"].calls,
+                    arguments["smooth_time"].calls,
+                    arguments["order"].calls),
+                "step": merge_dependencies(
+                    arguments["input"].calls,
+                    arguments["smooth_time"].calls)
+            }
+
+        }
+        return BuildAST(
+            expression=arguments["name"] + "()",
+            calls={arguments["name"]: 1},
+            subscripts=self.def_subs,
+            order=0)
+
+
+class TrendBuilder(StructureBuilder):
+    def __init__(self, trend_str,  component):
+        super().__init__(None, component)
+        self.arguments = {
+            "input": trend_str.input,
+            "average_time": trend_str.average_time,
+            "initial": trend_str.initial,
+        }
+
+    def build(self, arguments):
+        self.component.type = "Stateful"
+        self.component.subtype = "Trend"
+        self.section.imports.add("statefuls", "Trend")
+        arguments["input"].reshape(self.section.subscripts, self.def_subs)
+        arguments["average_time"].reshape(self.section.subscripts, self.def_subs)
+        arguments["initial"].reshape(self.section.subscripts, self.def_subs)
+        arguments["name"] = self.section.namespace.make_python_identifier(
+            self.element.identifier, prefix="_trend")
+
+        self.element.objects[arguments["name"]] = {
+            "name": arguments["name"],
+            "expression": "%(name)s = Trend(lambda: %(input)s, "
+                          "lambda: %(average_time)s, lambda: %(initial)s, "
+                          "'%(name)s')" % arguments,
+            "calls": {
+                "initial": merge_dependencies(
+                    arguments["initial"].calls,
+                    arguments["input"].calls,
+                    arguments["average_time"].calls),
+                "step": merge_dependencies(
+                    arguments["input"].calls,
+                    arguments["average_time"].calls)
+            }
+
+        }
+        return BuildAST(
+            expression=arguments["name"] + "()",
+            calls={arguments["name"]: 1},
+            subscripts=self.def_subs,
+            order=0)
+
+
+class ForecastBuilder(StructureBuilder):
+    def __init__(self, forecast_str,  component):
+        super().__init__(None, component)
+        self.arguments = {
+            "input": forecast_str.input,
+            "average_time": forecast_str.average_time,
+            "horizon": forecast_str.horizon,
+        }
+
+    def build(self, arguments):
+        self.component.type = "Stateful"
+        self.component.subtype = "Forecast"
+        self.section.imports.add("statefuls", "Forecast")
+        arguments["input"].reshape(self.section.subscripts, self.def_subs)
+        arguments["average_time"].reshape(self.section.subscripts, self.def_subs)
+        arguments["horizon"].reshape(self.section.subscripts, self.def_subs)
+        arguments["name"] = self.section.namespace.make_python_identifier(
+            self.element.identifier, prefix="_forecast")
+
+        self.element.objects[arguments["name"]] = {
+            "name": arguments["name"],
+            "expression": "%(name)s = Forecast(lambda: %(input)s, "
+                          "lambda: %(average_time)s, lambda: %(horizon)s, "
+                          "'%(name)s')" % arguments,
+            "calls": {
+                "initial":
+                    arguments["input"].calls,
+                "step": merge_dependencies(
+                    arguments["input"].calls,
+                    arguments["average_time"].calls,
+                    arguments["horizon"].calls)
+            }
+
+        }
+        return BuildAST(
+            expression=arguments["name"] + "()",
+            calls={arguments["name"]: 1},
+            subscripts=self.def_subs,
+            order=0)
+
+
+class SampleIfTrueBuilder(StructureBuilder):
+    def __init__(self, sampleiftrue_str,  component):
+        super().__init__(None, component)
+        self.arguments = {
+            "condition": sampleiftrue_str.condition,
+            "input": sampleiftrue_str.input,
+            "initial": sampleiftrue_str.initial,
+        }
+
+    def build(self, arguments):
+        self.component.type = "Stateful"
+        self.component.subtype = "SampleIfTrue"
+        self.section.imports.add("statefuls", "SampleIfTrue")
+        arguments["condition"].reshape(self.section.subscripts, self.def_subs)
+        arguments["input"].reshape(self.section.subscripts, self.def_subs)
+        arguments["initial"].reshape(self.section.subscripts, self.def_subs)
+        arguments["name"] = self.section.namespace.make_python_identifier(
+            self.element.identifier, prefix="_sampleiftrue")
+
+        self.element.objects[arguments["name"]] = {
+            "name": arguments["name"],
+            "expression": "%(name)s = SampleIfTrue(lambda: %(condition)s, "
+                          "lambda: %(input)s, lambda: %(initial)s, "
+                          "'%(name)s')" % arguments,
+            "calls": {
+                "initial":
+                    arguments["initial"].calls,
+                "step": merge_dependencies(
+                    arguments["condition"].calls,
+                    arguments["input"].calls)
+            }
+
+        }
+        return BuildAST(
+            expression=arguments["name"] + "()",
+            calls={arguments["name"]: 1},
+            subscripts=self.def_subs,
+            order=0)
+
+
+class LookupsBuilder(StructureBuilder):
+    def __init__(self, lookups_str,  component):
+        super().__init__(None, component)
+        self.arguments = {}
+        self.x = lookups_str.x
+        self.y = lookups_str.y
+
+    def build(self, arguments):
+        self.component.type = "Lookup"
+        self.component.subtype = "Normal"
+        arguments["x"] = np.array2string(
+            np.array(self.x),
+            separator=",",
+            threshold=len(self.x)
+        )
+        arguments["y"] = np.array2string(
+            np.array(self.y),
+            separator=",",
+            threshold=len(self.y)
+        )
+        arguments["subscripts"] = self.def_subs
+
+        if "hardcoded_lookups" in self.element.objects:
+            # object already exists
+            self.element.objects["hardcoded_lookups"]["expression"] += "\n\n"\
+                + self.element.objects["hardcoded_lookups"]["name"]\
+                + ".add(%(x)s, %(y)s, %(subscripts)s)" % arguments
+
+            return None
+        else:
+            # create a new object
+            self.section.imports.add("lookups", "HardcodedLookups")
+
+            arguments["name"] = self.section.namespace.make_python_identifier(
+                self.element.identifier, prefix="_hardcodedlookup")
+
+            self.element.objects["hardcoded_lookups"] = {
+                "name": arguments["name"],
+                "expression": "%(name)s = HardcodedLookups(%(x)s, %(y)s, "
+                              "%(subscripts)s, '%(name)s')" % arguments
+            }
+
+            return BuildAST(
+                expression=arguments["name"] + "(x)",
+                calls={"__lookup__": None},
+                subscripts=self.def_subs,
+                order=0)
+
+
+class InlineLookupsBuilder(StructureBuilder):
+    def __init__(self, inlinelookups_str,  component):
+        super().__init__(None, component)
+        self.arguments = {
+            "value": inlinelookups_str.argument
+        }
+        self.lookups = inlinelookups_str.lookups
+
+    def build(self, arguments):
+        self.component.type = "Auxiliary"
+        self.component.subtype = "with Lookup"
+        self.section.imports.add("numpy")
+        arguments["x"] = np.array2string(
+            np.array(self.lookups.x),
+            separator=",",
+            threshold=len(self.lookups.x)
+        )
+        arguments["y"] = np.array2string(
+            np.array(self.lookups.y),
+            separator=",",
+            threshold=len(self.lookups.y)
+        )
+        return BuildAST(
+            expression="np.interp(%(value)s, %(x)s, %(y)s)" % arguments,
+            calls=arguments["value"].calls,
+            subscripts=arguments["value"].subscripts,
+            order=0)
+
+
+class ReferenceBuilder(StructureBuilder):
+    def __init__(self, reference_str,  component):
+        super().__init__(None, component)
+        self.mapping_subscripts = {}
+        self.reference = reference_str.reference
+        self.subscripts = reference_str.subscripts
+        self.arguments = {}
+        self.section.imports.add("xarray")
+
+    @property
+    def subscripts(self):
+        return self._subscripts
+
+    @subscripts.setter
+    def subscripts(self, subscripts):
+        """Get subscript dictionary from reference"""
+        self._subscripts = self.section.subscripts.make_coord_dict(
+            getattr(subscripts, "subscripts", {}))
+
+        # get the subscripts after applying the mapping if necessary
+        for dim, coordinates in self._subscripts.items():
+            if len(coordinates) > 1:
+                # we create the mapping only with those subscripts that are
+                # ranges as we need to ignore singular subscripts because
+                # that dimension is removed from final element
+                if dim not in self.def_subs and not dim.endswith("!"):
+                    # the reference has a subscripts which is it not
+                    # applied (!) and does not appear in the definition
+                    # of the variable
+                    for mapped in self.section.subscripts.mapping[dim]:
+                        # check the mapped subscripts
+                        # TODO update this and the parser to make it
+                        # compatible with more complex mappings
+                        if mapped in self.def_subs\
+                          and mapped not in self._subscripts:
+                            # the mapped subscript appears in the definition
+                            # and it is not already in the variable
+                            self.mapping_subscripts[mapped] =\
+                                self.section.subscripts.subscripts[mapped]
+                            break
+                else:
+                    # the subscript is in the variable definition,
+                    # do not change it
+                    self.mapping_subscripts[dim] = coordinates
+
+    def build(self, arguments):
+        if self.reference not in self.section.namespace.cleanspace:
+            # Manage references to subscripts (subscripts used as variables)
+            expression, subscripts =\
+                self.section.subscripts.subscript2num[self.reference]
+            subscripts_out = self.section.subscripts.simplify_subscript_input(
+                subscripts, list(subscripts))[1]
+            if subscripts:
+                self.section.imports.add("numpy")
+                # NUMPY: not need this if
+                expression = "xr.DataArray(%s, %s, %s)" % (
+                    expression, subscripts_out, list(subscripts))
+            return BuildAST(
+                expression=expression,
+                calls={},
+                subscripts=subscripts,
+                order=0)
+
+        reference = self.section.namespace.cleanspace[self.reference]
+
+        # TODO lookups are passed as a reference first, in that case we will
+        # need to replace () in the lookup call
+        expression = reference + "()"
+
+        if not self.subscripts:
+            return BuildAST(
+                expression=expression,
+                calls={reference: 1},
+                subscripts={},
+                order=0)
+
+        original_subs = self.section.subscripts.make_coord_dict(
+                    self.section.subscripts.elements[reference])
+
+        expression, final_subs = self.visit_subscripts(expression, original_subs)
+
+        return BuildAST(
+            expression=expression,
+            calls={reference: 1},
+            subscripts=final_subs,
+            order=0)
+
+    def visit_subscripts(self, expression, original_subs):
+        final_subs, rename, loc, reset_coords, float = {}, {}, [], False, True
+        for (dim, coord), (orig_dim, orig_coord)\
+          in zip(self.subscripts.items(), original_subs.items()):
+            if len(coord) == 1:
+                # subset a 1 dimension value
+                # NUMPY: subset value [:, N, :, :]
+                loc.append(repr(coord[0]))
+                reset_coords = True
+            elif len(coord) < len(orig_coord):
+                # subset a subrange
+                # NUMPY: subset value [:, :, np.array([1, 0]), :]
+                # NUMPY: as order may change we need to check if dim != orig_dim
+                # NUMPY: use also ranges [:, :, 2:5, :] when possible
+                loc.append("_subscript_dict['%s']" % dim)
+                final_subs[dim] = coord
+                float = False
+            else:
+                # do nothing
+                # NUMPY: same, we can remove float = False
+                loc.append(":")
+                final_subs[dim] = coord
+                float = False
+
+            if dim != orig_dim and len(coord) != 1:
+                # NUMPY: check order of dimensions, make all subranges work with the same dimensions?
+                # NUMPY: this could be solved in the previous if/then/else
+                rename[orig_dim] = dim
+
+        if any(dim != ":" for dim in loc):
+            # NUMPY: expression += "[%s]" % ", ".join(loc)
+            expression += ".loc[%s]" % ", ".join(loc)
+        if reset_coords and float:
+            # NUMPY: Not neccessary
+            expression = "float(" + expression + ")"
+        elif reset_coords:
+            # NUMPY: Not neccessary
+            expression += ".reset_coords(drop=True)"
+        if rename:
+            # NUMPY: Not neccessary
+            expression += ".rename(%s)" % rename
+
+        # NUMPY: This will not be necessary, we only need to return
+        #        self.mapping_subscripts
+        if self.mapping_subscripts != final_subs:
+            subscripts_out = self.section.subscripts.simplify_subscript_input(
+                self.mapping_subscripts, list(self.mapping_subscripts))[1]
+            expression = "xr.DataArray(%s.values, %s, %s)" % (
+                expression, subscripts_out, list(self.mapping_subscripts)
+            )
+
+        return expression, self.mapping_subscripts
+
+
+class NumericBuilder(StructureBuilder):
+    # Standard class, inherit all from StructureBuilder
+    pass
+
+
+class ArrayBuilder(StructureBuilder):
+    # Standard class, inherit all from StructureBuilder
+    def build(self, arguments):
+        self.value = np.array2string(
+            self.value.reshape(compute_shape(self.def_subs)),
+            separator=",",
+            threshold=np.prod(self.value.shape)
+        )
+        self.component.type = "Constant"
+        self.component.subtype = "Normal"
+
+        final_subs, subscripts_out =\
+            self.section.subscripts.simplify_subscript_input(
+                self.def_subs, self.element.subscripts)
+
+        return BuildAST(
+            expression="xr.DataArray(%s, %s, %s)" % (
+                self.value, subscripts_out, list(final_subs)),
+            calls={},
+            subscripts=final_subs,
+            order=0)
+
+
+def merge_dependencies(*dependencies, inplace=False):
+    # TODO improve dependencies in the next major release, include info
+    # about external objects and simplify the stateful objects, think about
+    # how to include data/lookups objects
+    current = dependencies[0]
+    if inplace:
+        current = dependencies[0]
+    else:
+        current = dependencies[0].copy()
+    for new in dependencies[1:]:
+        if not current:
+            current.update(new)
+        elif new:
+            # regular element
+            _merge_dependencies(current, new)
+
+    return current
+
+
+def _merge_dependencies(current, new):
+    """
+    Merge two dependencies dicts of an element.
+
+    Parameters
+    ----------
+    current: dict
+        Current dependencies of the element. It will be mutated.
+
+    new: dict
+        New dependencies to add.
+
+    Returns
+    -------
+    None
+
+    """
+    current_set, new_set = set(current), set(new)
+    for dep in current_set.intersection(new_set):
+        if dep.startswith("__"):
+            # if it is special (__lookup__, __external__) continue
+            continue
+        # if dependency is in both sum the number of calls
+        if dep in ["initial", "step"]:
+            _merge_dependencies(current[dep], new[dep])
+        else:
+            current[dep] += new[dep]
+    for dep in new_set.difference(current_set):
+        # if dependency is only in new copy it
+        current[dep] = new[dep]
+
+
+class ASTVisitor:
+    builders = {
+        ct.InitialStructure: InitialBuilder,
+        ct.IntegStructure: IntegBuilder,
+        ct.DelayStructure: lambda x, y: DelayBuilder("Delay", x, y),
+        ct.DelayNStructure: lambda x, y: DelayBuilder("DelayN", x, y),
+        ct.DelayFixedStructure: DelayFixedBuilder,
+        ct.SmoothStructure: SmoothBuilder,
+        ct.SmoothNStructure: SmoothBuilder,
+        ct.TrendStructure: TrendBuilder,
+        ct.ForecastStructure: ForecastBuilder,
+        ct.SampleIfTrueStructure: SampleIfTrueBuilder,
+        ct.GetConstantsStructure: ExtConstantBuilder,
+        ct.GetDataStructure: ExtDataBuilder,
+        ct.GetLookupsStructure: ExtLookupBuilder,
+        ct.LookupsStructure: LookupsBuilder,
+        ct.InlineLookupsStructure: InlineLookupsBuilder,
+        ct.DataStructure: TabDataBuilder,
+        ct.ReferenceStructure: ReferenceBuilder,
+        ct.CallStructure: CallBuilder,
+        ct.GameStructure: GameBuilder,
+        ct.LogicStructure: OperationBuilder,
+        ct.ArithmeticStructure: OperationBuilder,
+        int: NumericBuilder,
+        float: NumericBuilder,
+        np.ndarray: ArrayBuilder
+    }
+
+    def __init__(self, component):
+        self.ast = component.ast
+        self.subscripts = component.subscripts_dict
+        self.component = component
+        # TODO add a attribute for "new structures"
+
+    def visit(self):
+        # TODO: if final_subscripts == self.subscripts OK, else -> redimension
+        visit_out = self._visit(self.ast)
+
+        if not visit_out:
+            # external objects that are declared with other expression
+            return None
+
+        if not visit_out.calls and self.component.type == "Auxiliary":
+            self.component.type = "Constant"
+            self.component.subtype = "Normal"
+
+        # include dependencies of the current component in the element
+        merge_dependencies(
+            self.component.element.dependencies,
+            visit_out.calls,
+            inplace=True)
+
+        if not visit_out.subscripts\
+          and self.subscripts != self.component.element.subs_dict:
+            return visit_out
+
+        # NUMPY not needed
+        # get subscript in elements as name of the ranges may change
+        subscripts_in_element = {
+            dim: coords
+            for dim, coords
+            in zip(self.component.element.subscripts, self.subscripts.values())
+        }
+
+        reshape = (
+            (visit_out.subscripts != self.subscripts
+             or list(visit_out.subscripts) != list(self.subscripts))
+            and
+            (visit_out.subscripts != subscripts_in_element
+             or list(visit_out.subscripts) != list(subscripts_in_element))
+        )
+
+        if reshape:
+            # We are only comparing the dictionaries (set of dimensions)
+            # and not the list (order).
+            # With xarray we don't need to compare the order because the
+            # decorator @subs will reorder the objects
+            # NUMPY: in this case we need to tile along dims if neccessary
+            #        or reorder the dimensions
+            # NUMPY: if the output is a float or int and they are several
+            #        definitions we can return float or int as we can
+            #        safely do "var[:, 1, :] = 3"
+            visit_out.reshape(
+                self.component.section.subscripts, self.subscripts)
+
+        return visit_out
+
+    def _visit(self, ast_object):
+        builder = self.builders[type(ast_object)](ast_object, self.component)
+        arguments = {name: self._visit(value) for name, value in builder.arguments.items()}
+        return builder.build(arguments)
diff --git a/pysd/py_backend/components.py b/pysd/py_backend/components.py
index 5eccf370..6d4c41f0 100644
--- a/pysd/py_backend/components.py
+++ b/pysd/py_backend/components.py
@@ -4,6 +4,7 @@
 
 import os
 import random
+import numpy as np
 from importlib.machinery import SourceFileLoader
 
 from pysd._version import __version__
@@ -84,6 +85,8 @@ def _set_component(self, name, value):
 
 
 class Time(object):
+    rprec = 1e-10  # relative precission for final time and saving time
+
     def __init__(self):
         self._time = None
         self.stage = None
@@ -135,7 +138,7 @@ def in_bounds(self):
             True if time is smaller than final time. Otherwise, returns Fase.
 
         """
-        return self._time < self.final_time()
+        return self._time + self.time_step()*self.rprec < self.final_time()
 
     def in_return(self):
         """ Check if current time should be returned """
@@ -144,9 +147,15 @@ def in_return(self):
 
         time_delay = self._time - self._initial_time
         save_per = self.saveper()
-        prec = self.time_step() * 1e-10
+        prec = self.time_step() * self.rprec
         return time_delay % save_per < prec or -time_delay % save_per < prec
 
+    def round(self):
+        """ Return rounded time to outputs to avoid float precission error"""
+        return np.round(
+            self._time,
+            -int(np.log10(self.time_step()*self.rprec)))
+
     def add_return_timestamps(self, return_timestamps):
         """ Add return timestamps """
         if return_timestamps is None or hasattr(return_timestamps, '__len__'):
diff --git a/pysd/py_backend/data.py b/pysd/py_backend/data.py
index 4a69d6fa..a46af99f 100644
--- a/pysd/py_backend/data.py
+++ b/pysd/py_backend/data.py
@@ -1,5 +1,6 @@
 import warnings
 import re
+import random
 from pathlib import Path
 
 import numpy as np
@@ -50,6 +51,8 @@ def read_file(cls, file_name, encoding=None):
             indicate if the output file is transposed.
 
         """
+        # in the most cases variables will be split per columns, then
+        # read the first row to have all the column names
         out = cls.read_line(file_name, encoding)
         if out is None:
             raise ValueError(
@@ -59,10 +62,16 @@ def read_file(cls, file_name, encoding=None):
         transpose = False
 
         try:
-            [float(col) for col in out]
-            out = cls.read_row(file_name, encoding)
+            # if we fail converting columns to float then they are
+            # not numeric values, so current direction is okay
+            [float(col) for col in random.sample(out, 3)]
+            # we did not fail, read the first column to see if variables
+            # are split per rows
+            out = cls.read_col(file_name, encoding)
             transpose = True
-            [float(col) for col in out]
+            # if we still are able to transform values to float the
+            # file is not valid
+            [float(col) for col in random.sample(out, 3)]
         except ValueError:
             return out, transpose
         else:
@@ -91,7 +100,7 @@ def read_line(cls, file_name, encoding=None):
             return None
 
     @classmethod
-    def read_row(cls, file_name, encoding=None):
+    def read_col(cls, file_name, encoding=None):
         """
         Read the firts column and return a set of it.
         """
@@ -190,9 +199,9 @@ def __call__(self, time):
                 outdata = self.data[0]
             elif self.interp == "interpolate":
                 outdata = self.data.interp(time=time)
-            elif self.interp == 'look forward':
+            elif self.interp == 'look_forward':
                 outdata = self.data.sel(time=time, method="backfill")
-            elif self.interp == 'hold backward':
+            elif self.interp == 'hold_backward':
                 outdata = self.data.sel(time=time, method="pad")
 
             if self.is_float:
@@ -214,16 +223,23 @@ def __call__(self, time):
 
 class TabData(Data):
     """
-    Data from tabular file tab/cls, it could be from Vensim output.
+    Data from tabular file tab/csv, it could be from Vensim output.
     """
     def __init__(self, real_name, py_name, coords, interp="interpolate"):
         self.real_name = real_name
         self.py_name = py_name
         self.coords = coords
-        self.interp = interp
+        self.interp = interp.replace(" ", "_") if interp else None
         self.is_float = not bool(coords)
         self.data = None
 
+        if self.interp not in ["interpolate", "raw",
+                               "look_forward", "hold_backward"]:
+            raise ValueError(self.py_name + "\n"
+                             + " The interpolation method (interp) must be "
+                             + "'raw', 'interpolate', "
+                             + "'look_forward' or 'hold_backward")
+
     def load_data(self, file_names):
         """
         Load data values from files.
diff --git a/pysd/py_backend/decorators.py b/pysd/py_backend/decorators.py
index f796dee1..195a5c77 100644
--- a/pysd/py_backend/decorators.py
+++ b/pysd/py_backend/decorators.py
@@ -14,6 +14,7 @@ def subs(dims, subcoords):
     """
     def decorator(function):
         function.dims = dims
+        function.args = inspect.getfullargspec(function)[0]
 
         @wraps(function)
         def wrapper(*args):
diff --git a/pysd/py_backend/external.py b/pysd/py_backend/external.py
index c54a815a..9da02cb9 100644
--- a/pysd/py_backend/external.py
+++ b/pysd/py_backend/external.py
@@ -13,6 +13,7 @@
 from openpyxl import load_workbook
 from . import utils
 from .data import Data
+from .lookups import Lookups
 
 
 class Excels():
@@ -551,9 +552,9 @@ def _interpolate_missing(self, x, xr, yr):
                 y[i] = yr[-1]
             elif value <= xr[0]:
                 y[i] = yr[0]
-            elif self.interp == 'look forward':
+            elif self.interp == 'look_forward':
                 y[i] = yr[xr >= value][0]
-            elif self.interp == 'hold backward':
+            elif self.interp == 'hold_backward':
                 y[i] = yr[xr <= value][-1]
             else:
                 y[i] = np.interp(value, xr, yr)
@@ -705,7 +706,8 @@ def __init__(self, file_name, sheet, time_row_or_col, cell,
         self.cells = [cell]
         self.coordss = [coords]
         self.root = root
-        self.interp = interp
+        # TODO remove in 3.0.0 (self.interp = interp)
+        self.interp = interp.replace(" ", "_") if interp else None
         self.is_float = not bool(coords)
 
         # check if the interpolation method is valid
@@ -713,11 +715,11 @@ def __init__(self, file_name, sheet, time_row_or_col, cell,
             self.interp = "interpolate"
 
         if self.interp not in ["interpolate", "raw",
-                               "look forward", "hold backward"]:
+                               "look_forward", "hold_backward"]:
             raise ValueError(self.py_name + "\n"
                              + " The interpolation method (interp) must be "
                              + "'raw', 'interpolate', "
-                             + "'look forward' or 'hold backward")
+                             + "'look_forward' or 'hold_backward")
 
     def add(self, file_name, sheet, time_row_or_col, cell,
             interp, coords):
@@ -732,7 +734,7 @@ def add(self, file_name, sheet, time_row_or_col, cell,
 
         if not interp:
             interp = "interpolate"
-        if interp != self.interp:
+        if interp.replace(" ", "_") != self.interp:
             raise ValueError(self.py_name + "\n"
                              + "Error matching interpolation method with "
                              + "previously defined one")
@@ -753,7 +755,7 @@ def initialize(self):
                    self.cells, self.coordss)])
 
 
-class ExtLookup(External):
+class ExtLookup(External, Lookups):
     """
     Class for Vensim GET XLS LOOKUPS/GET DIRECT LOOKUPS
     """
@@ -768,6 +770,7 @@ def __init__(self, file_name, sheet, x_row_or_col, cell,
         self.root = root
         self.coordss = [coords]
         self.interp = "interpolate"
+        self.is_float = not bool(coords)
 
     def add(self, file_name, sheet, x_row_or_col, cell, coords):
         """
@@ -794,58 +797,6 @@ def initialize(self):
             in zip(self.files, self.sheets, self.x_row_or_cols,
                    self.cells, self.coordss)])
 
-    def __call__(self, x):
-        return self._call(self.data, x)
-
-    def _call(self, data, x):
-        if isinstance(x, xr.DataArray):
-            if not x.dims:
-                # shape 0 xarrays
-                return self._call(data, float(x))
-            if np.all(x > data['lookup_dim'].values[-1]):
-                outdata, _ = xr.broadcast(data[-1], x)
-                warnings.warn(
-                  self.py_name + "\n"
-                  + "extrapolating data above the maximum value of the series")
-            elif np.all(x < data['lookup_dim'].values[0]):
-                outdata, _ = xr.broadcast(data[0], x)
-                warnings.warn(
-                  self.py_name + "\n"
-                  + "extrapolating data below the minimum value of the series")
-            else:
-                data, _ = xr.broadcast(data, x)
-                outdata = data[0].copy()
-                for a in utils.xrsplit(x):
-                    outdata.loc[a.coords] = self._call(
-                        data.loc[a.coords],
-                        float(a))
-            # the output will be always an xarray
-            return outdata.reset_coords('lookup_dim', drop=True)
-
-        else:
-            if x in data['lookup_dim'].values:
-                outdata = data.sel(lookup_dim=x)
-            elif x > data['lookup_dim'].values[-1]:
-                outdata = data[-1]
-                warnings.warn(
-                  self.py_name + "\n"
-                  + "extrapolating data above the maximum value of the series")
-            elif x < data['lookup_dim'].values[0]:
-                outdata = data[0]
-                warnings.warn(
-                  self.py_name + "\n"
-                  + "extrapolating data below the minimum value of the series")
-            else:
-                outdata = data.interp(lookup_dim=x)
-
-            # the output could be a float or an xarray
-            if self.coordss[0]:
-                # Remove lookup dimension coord from the DataArray
-                return outdata.reset_coords('lookup_dim', drop=True)
-            else:
-                # if lookup has no-coords return a float
-                return float(outdata)
-
 
 class ExtConstant(External):
     """
diff --git a/pysd/py_backend/functions.py b/pysd/py_backend/functions.py
index b609e48c..ea751e69 100644
--- a/pysd/py_backend/functions.py
+++ b/pysd/py_backend/functions.py
@@ -240,7 +240,19 @@ def if_then_else(condition, val_if_true, val_if_false):
     The value depending on the condition.
 
     """
+    # NUMPY: replace xr by np
     if isinstance(condition, xr.DataArray):
+        # NUMPY: neccessarry for keep the same shape always
+        # if condition.all():
+        #    value = val_if_true()
+        # elif not condition.any():
+        #    value = val_if_false()
+        # else:
+        #    return np.where(condition, val_if_true(), val_if_false())
+        #
+        # if isinstance(value, np.ndarray):
+        #    return value
+        # return np.full_like(condition, value)
         if condition.all():
             return val_if_true()
         elif not condition.any():
@@ -293,7 +305,7 @@ def logical_or(*args):
     return current
 
 
-def xidz(numerator, denominator, value_if_denom_is_zero):
+def xidz(numerator, denominator, x):
     """
     Implements Vensim's XIDZ function.
     https://www.vensim.com/documentation/fn_xidz.htm
@@ -304,26 +316,34 @@ def xidz(numerator, denominator, value_if_denom_is_zero):
     Parameters
     ----------
     numerator: float or xarray.DataArray
+        Numerator of the operation.
     denominator: float or xarray.DataArray
-        Components of the division operation
-    value_if_denom_is_zero: float or xarray.DataArray
-        The value to return if the denominator is zero
+        Denominator of the operation.
+    x: float or xarray.DataArray
+        The value to return if the denominator is zero.
 
     Returns
     -------
-    numerator / denominator if denominator > 1e-6
+    numerator/denominator if denominator > small_vensim
     otherwise, returns value_if_denom_is_zero
 
     """
+    # NUMPY: replace DataArray by np.ndarray, xr.where -> np.where
     if isinstance(denominator, xr.DataArray):
         return xr.where(np.abs(denominator) < small_vensim,
-                        value_if_denom_is_zero,
-                        numerator * 1.0 / denominator)
+                        x,
+                        numerator/denominator)
 
     if abs(denominator) < small_vensim:
-        return value_if_denom_is_zero
+        # NUMPY: neccessarry for keep the same shape always
+        # if isinstance(numerator, np.ndarray):
+        #    return np.full_like(numerator, x)
+        return x
     else:
-        return numerator * 1.0 / denominator
+        # NUMPY: neccessarry for keep the same shape always
+        # if isinstance(x, np.ndarray):
+        #    return np.full_like(x, numerator/denominator)
+        return numerator/denominator
 
 
 def zidz(numerator, denominator):
@@ -345,15 +365,21 @@ def zidz(numerator, denominator):
     otherwise zero.
 
     """
+    # NUMPY: replace DataArray by np.ndarray, xr.where -> np.where
     if isinstance(denominator, xr.DataArray):
         return xr.where(np.abs(denominator) < small_vensim,
                         0,
-                        numerator * 1.0 / denominator)
+                        numerator/denominator)
 
     if abs(denominator) < small_vensim:
+        # NUMPY: neccessarry for keep the same shape always
+        # if isinstance(denominator, np.ndarray):
+        #    return np.zeros_like(denominator)
+        if isinstance(numerator, xr.DataArray):
+            return xr.DataArray(0, numerator.coords, numerator.dims)
         return 0
     else:
-        return numerator * 1.0 / denominator
+        return numerator/denominator
 
 
 def active_initial(time, expr, init_val):
@@ -361,15 +387,19 @@ def active_initial(time, expr, init_val):
     Implements vensim's ACTIVE INITIAL function
     Parameters
     ----------
-    time: function
-        The current time function
-    expr
-    init_val
+    stage: str
+        The stage of the model.
+    expr: function
+        Running stage value
+    init_val: float or xarray.DataArray
+        Initialization stage value.
 
     Returns
     -------
 
     """
+    # TODO replace time by stage when doing a non compatible version
+    # NUMPY: both must have same dimensions in inputs, remove time.stage
     if time.stage == 'Initialization':
         return init_val
     else:
@@ -414,6 +444,7 @@ def log(x, base):
     float
       The log of 'x' in base 'base'.
     """
+    # TODO remove with PySD 3.0.0, log could be directly created in the file
     return np.log(x) / np.log(base)
 
 
@@ -431,6 +462,7 @@ def integer(x):
     Returns integer part of x.
 
     """
+    # NUMPY: replace xr by np
     if isinstance(x, xr.DataArray):
         return x.astype(int)
     else:
@@ -454,6 +486,7 @@ def quantum(a, b):
         If b > 0 returns b * integer(a/b). Otherwise, returns a.
 
     """
+    # NUMPY: replace xr by np
     if isinstance(b, xr.DataArray):
         return xr.where(b < small_vensim, a, b*integer(a/b))
     if b < small_vensim:
@@ -500,6 +533,7 @@ def sum(x, dim=None):
       The result of the sum operation in the given dimensions.
 
     """
+    # NUMPY: replace by np.sum(x, axis=axis) put directly in the file
     # float returned if the function is applied over all the dimensions
     if dim is None or set(x.dims) == set(dim):
         return float(x.sum())
@@ -525,6 +559,7 @@ def prod(x, dim=None):
       The result of the product operation in the given dimensions.
 
     """
+    # NUMPY: replace by np.prod(x, axis=axis) put directly in the file
     # float returned if the function is applied over all the dimensions
     if dim is None or set(x.dims) == set(dim):
         return float(x.prod())
@@ -550,6 +585,7 @@ def vmin(x, dim=None):
       The result of the minimum value over the given dimensions.
 
     """
+    # NUMPY: replace by np.min(x, axis=axis) put directly in the file
     # float returned if the function is applied over all the dimensions
     if dim is None or set(x.dims) == set(dim):
         return float(x.min())
@@ -575,6 +611,7 @@ def vmax(x, dim=None):
       The result of the maximum value over the dimensions.
 
     """
+    # NUMPY: replace by np.max(x, axis=axis) put directly in the file
     # float returned if the function is applied over all the dimensions
     if dim is None or set(x.dims) == set(dim):
         return float(x.max())
@@ -599,4 +636,6 @@ def invert_matrix(mat):
         Inverted matrix.
 
     """
+    # NUMPY: avoid converting to xarray, put directly the expression
+    # in the model
     return xr.DataArray(np.linalg.inv(mat.values), mat.coords, mat.dims)
diff --git a/pysd/py_backend/lookups.py b/pysd/py_backend/lookups.py
new file mode 100644
index 00000000..98ef2132
--- /dev/null
+++ b/pysd/py_backend/lookups.py
@@ -0,0 +1,115 @@
+import warnings
+
+import numpy as np
+import xarray as xr
+
+from . import utils
+
+
+class Lookups(object):
+    # TODO add __init__ and use this class for used input pandas.Series
+    # as Lookups
+    # def __init__(self, data, coords, interp="interpolate"):
+
+    def __call__(self, x):
+        return self._call(self.data, x)
+
+    def _call(self, data, x):
+        if isinstance(x, xr.DataArray):
+            if not x.dims:
+                # shape 0 xarrays
+                return self._call(data, float(x))
+            if np.all(x > data['lookup_dim'].values[-1]):
+                outdata, _ = xr.broadcast(data[-1], x)
+                warnings.warn(
+                  self.py_name + "\n"
+                  + "extrapolating data above the maximum value of the series")
+            elif np.all(x < data['lookup_dim'].values[0]):
+                outdata, _ = xr.broadcast(data[0], x)
+                warnings.warn(
+                  self.py_name + "\n"
+                  + "extrapolating data below the minimum value of the series")
+            else:
+                data, _ = xr.broadcast(data, x)
+                outdata = data[0].copy()
+                for a in utils.xrsplit(x):
+                    outdata.loc[a.coords] = self._call(
+                        data.loc[a.coords],
+                        float(a))
+            # the output will be always an xarray
+            return outdata.reset_coords('lookup_dim', drop=True)
+
+        else:
+            if x in data['lookup_dim'].values:
+                outdata = data.sel(lookup_dim=x)
+            elif x > data['lookup_dim'].values[-1]:
+                outdata = data[-1]
+                warnings.warn(
+                  self.py_name + "\n"
+                  + "extrapolating data above the maximum value of the series")
+            elif x < data['lookup_dim'].values[0]:
+                outdata = data[0]
+                warnings.warn(
+                  self.py_name + "\n"
+                  + "extrapolating data below the minimum value of the series")
+            else:
+                outdata = data.interp(lookup_dim=x)
+
+            # the output could be a float or an xarray
+            if self.is_float:
+                # if lookup has no-coords return a float
+                return float(outdata)
+            else:
+                # Remove lookup dimension coord from the DataArray
+                return outdata.reset_coords('lookup_dim', drop=True)
+
+
+class HardcodedLookups(Lookups):
+    """Class for lookups defined in the file"""
+
+    def __init__(self, x, y, coords, py_name):
+        # TODO: avoid add and merge all declarations in one definition
+        self.is_float = not bool(coords)
+        self.py_name = py_name
+        self.data = xr.DataArray(
+            np.array(y).reshape(tuple([len(x)] + utils.compute_shape(coords))),
+            {"lookup_dim": x, **coords},
+            ["lookup_dim"] + list(coords)
+        )
+        self.x = set(x)
+
+    def add(self, x, y, coords):
+        self.data = self.data.combine_first(
+            xr.DataArray(
+                np.array(y).reshape(tuple([len(x)] + utils.compute_shape(coords))),
+                {"lookup_dim": x, **coords},
+                ["lookup_dim"] + list(coords)
+            ))
+
+        if np.any(np.isnan(self.data)):
+            # fill missing values of different input lookup_dim values
+            values = self.data.values
+            self._fill_missing(self.data.lookup_dim.values, values)
+            self.data = xr.DataArray(values, self.data.coords, self.data.dims)
+
+    def _fill_missing(self, series,  data):
+        """
+        Fills missing values in lookups to have a common series.
+        Mutates the values in data.
+
+        Returns
+        -------
+        None
+
+        """
+        if len(data.shape) > 1:
+            # break the data array until arrive to a vector
+            for i in range(data.shape[1]):
+                if np.any(np.isnan(data[:, i])):
+                    self._fill_missing(series, data[:, i])
+        elif not np.all(np.isnan(data)):
+            # interpolate missing values
+            data[np.isnan(data)] = np.interp(
+                series[np.isnan(data)],
+                series[~np.isnan(data)],
+                data[~np.isnan(data)])
diff --git a/pysd/py_backend/statefuls.py b/pysd/py_backend/statefuls.py
index 6567224e..8d36b743 100644
--- a/pysd/py_backend/statefuls.py
+++ b/pysd/py_backend/statefuls.py
@@ -889,7 +889,6 @@ def _isdynamic(self, dependencies):
             return True
         for dep in dependencies:
             if dep.startswith("_") and not dep.startswith("_initial_")\
-               and not dep.startswith("_active_initial_")\
                and not dep.startswith("__"):
                 return True
         return False
@@ -1148,6 +1147,7 @@ def get_series_data(self, param):
         func_name = utils.get_key_and_value_by_insensitive_key_or_value(
             param,
             self.components._namespace)[1] or param
+        print(func_name, self.get_args(getattr(self.components, func_name)))
 
         try:
             if func_name.startswith("_ext_"):
@@ -1427,6 +1427,26 @@ def doc(self):
         variable names, and understand how they are translated into
         python safe names.
 
+        Returns
+        -------
+        docs_df: pandas dataframe
+            Dataframe with columns for the model components:
+                - Real names
+                - Python safe identifiers (as used in model.components)
+                - Units string
+                - Documentation strings from the original model file
+        """
+        warnings.warn(
+            "doc method will become an attribute in version 3.0.0...",
+            FutureWarning)
+        return self._doc
+
+    def _build_doc(self):
+        """
+        Formats a table of documentation strings to help users remember
+        variable names, and understand how they are translated into
+        python safe names.
+
         Returns
         -------
         docs_df: pandas dataframe
@@ -1454,25 +1474,43 @@ def doc(self):
                     eqn = '; '.join(
                         [line.strip() for line in lines[3:unit_line]])
 
-                collector.append(
-                    {'Real Name': name,
-                     'Py Name': varname,
-                     'Eqn': eqn,
-                     'Unit': lines[unit_line].replace("Units:", "").strip(),
-                     'Lims': lines[unit_line+1].replace("Limits:", "").strip(),
-                     'Type': lines[unit_line+2].replace("Type:", "").strip(),
-                     'Subs': lines[unit_line+3].replace("Subs:", "").strip(),
-                     'Comment': '\n'.join(lines[(unit_line+4):]).strip()})
+                vardoc = {
+                    'Real Name': name,
+                    'Py Name': varname,
+                    'Eqn': eqn,
+                    'Unit': lines[unit_line].replace("Units:", "").strip(),
+                    'Lims': lines[unit_line+1].replace("Limits:", "").strip(),
+                    'Type': lines[unit_line+2].replace("Type:", "").strip()
+                }
+
+                if "Subtype:" in lines[unit_line+3]:
+                    vardoc["Subtype"] =\
+                        lines[unit_line+3].replace("Subtype:", "").strip()
+                    vardoc["Subs"] =\
+                        lines[unit_line+4].replace("Subs:", "").strip()
+                    vardoc["Comment"] =\
+                        '\n'.join(lines[(unit_line+5):]).strip()
+                else:
+                    vardoc["Subtype"] = None
+                    vardoc["Subs"] =\
+                        lines[unit_line+3].replace("Subs:", "").strip()
+                    vardoc["Comment"] =\
+                        '\n'.join(lines[(unit_line+4):]).strip()
+
+                collector.append(vardoc)
             except Exception:
                 pass
 
-        docs_df = pd.DataFrame(collector)
-        docs_df.fillna('None', inplace=True)
-
-        order = ['Real Name', 'Py Name', 'Unit', 'Lims',
-                 'Type', 'Subs', 'Eqn', 'Comment']
-        return docs_df[order].sort_values(
-            by='Real Name').reset_index(drop=True)
+        if collector:
+            docs_df = pd.DataFrame(collector)
+            docs_df.fillna("None", inplace=True)
+            order = ["Real Name", "Py Name", "Unit", "Lims",
+                     "Type", "Subtype", "Subs", "Eqn", "Comment"]
+            return docs_df[order].sort_values(
+                by="Real Name").reset_index(drop=True)
+        else:
+            # manage models with no documentation (mainly test models)
+            return None
 
     def __str__(self):
         """ Return model source files """
@@ -1496,6 +1534,7 @@ def __init__(self, py_model_file, data_files, initialize, missing_values):
         self.time.set_control_vars(**self.components._control_vars)
         self.data_files = data_files
         self.missing_values = missing_values
+        self._doc = self._build_doc()
         if initialize:
             self.initialize()
 
@@ -2125,8 +2164,9 @@ def _integrate(self, capture_elements):
 
         while self.time.in_bounds():
             if self.time.in_return():
-                outputs.at[self.time()] = [getattr(self.components, key)()
-                                           for key in capture_elements]
+                outputs.at[self.time.round()] = [
+                    getattr(self.components, key)()
+                    for key in capture_elements]
             self._euler_step(self.time.time_step())
             self.time.update(self.time()+self.time.time_step())
             self.clean_caches()
@@ -2135,8 +2175,8 @@ def _integrate(self, capture_elements):
         # need to add one more time step, because we run only the state
         # updates in the previous loop and thus may be one short.
         if self.time.in_return():
-            outputs.at[self.time()] = [getattr(self.components, key)()
-                                       for key in capture_elements]
+            outputs.at[self.time.round()] = [getattr(self.components, key)()
+                                             for key in capture_elements]
 
         progressbar.finish()
 
diff --git a/pysd/pysd.py b/pysd/pysd.py
index 69eaccf9..1de73930 100644
--- a/pysd/pysd.py
+++ b/pysd/pysd.py
@@ -8,6 +8,7 @@
 import sys
 from .py_backend.statefuls import Model
 
+
 if sys.version_info[:2] < (3, 7):  # pragma: no cover
     raise RuntimeError(
         "\n\n"
@@ -69,7 +70,7 @@ def read_xmile(xmile_file, data_files=None, initialize=True,
 
 def read_vensim(mdl_file, data_files=None, initialize=True,
                 missing_values="warning", split_views=False,
-                encoding=None, **kwargs):
+                encoding=None, old=False, **kwargs):
     """
     Construct a model from Vensim `.mdl` file.
 
@@ -124,9 +125,21 @@ def read_vensim(mdl_file, data_files=None, initialize=True,
     >>> model = read_vensim('../tests/test-models/samples/teacup/teacup.mdl')
 
     """
-    from .translation.vensim.vensim2py import translate_vensim
+    if old:
+        from .translation.vensim.vensim2py import translate_vensim
+        py_model_file = translate_vensim(mdl_file, split_views, encoding, **kwargs)
+    else:
+        from pysd.translation.vensim.vensin_file import VensimFile
+        from pysd.building.python.python_builder import ModelBuilder
+        ven_file = VensimFile(mdl_file)
+        ven_file.parse()
+        if split_views:
+            subview_sep = kwargs.get("subview_sep", "")
+            ven_file.parse_sketch(subview_sep)
+
+        abs_model = ven_file.get_abstract_model()
+        py_model_file = ModelBuilder(abs_model).build_model()
 
-    py_model_file = translate_vensim(mdl_file, split_views, encoding, **kwargs)
     model = load(py_model_file, data_files, initialize, missing_values)
     model.mdl_file = str(mdl_file)
     return model
diff --git a/pysd/tools/benchmarking.py b/pysd/tools/benchmarking.py
index 8012e35a..16754ac0 100644
--- a/pysd/tools/benchmarking.py
+++ b/pysd/tools/benchmarking.py
@@ -2,18 +2,18 @@
 Benchmarking tools for testing and comparing outputs between different files.
 Some of these functions are also used for testing.
 """
-
-import os.path
 import warnings
+from pathlib import Path
 
 import numpy as np
 import pandas as pd
 
-from pysd import read_vensim, read_xmile
+from pysd import read_vensim, read_xmile, load
 from ..py_backend.utils import load_outputs, detect_encoding
 
 
-def runner(model_file, canonical_file=None, transpose=False, data_files=None):
+def runner(model_file, canonical_file=None, transpose=False, data_files=None,
+           old=False):
     """
     Translates and runs a model and returns its output and the
     canonical output.
@@ -34,34 +34,42 @@ def runner(model_file, canonical_file=None, transpose=False, data_files=None):
     data_files: list (optional)
         List of the data files needed to run the model.
 
+    old: bool(optional)
+        If True use old translation method, used for testing backward compatibility.
+
     Returns
     -------
     output, canon: (pandas.DataFrame, pandas.DataFrame)
         pandas.DataFrame of the model output and the canonical output.
 
     """
-    directory = os.path.dirname(model_file)
+    if isinstance(model_file, str):
+        model_file = Path(model_file)
+
+    directory = model_file.parent
 
     # load canonical output
     if not canonical_file:
-        if os.path.isfile(os.path.join(directory, 'output.csv')):
-            canonical_file = os.path.join(directory, 'output.csv')
-        elif os.path.isfile(os.path.join(directory, 'output.tab')):
-            canonical_file = os.path.join(directory, 'output.tab')
+        if directory.joinpath('output.csv').is_file():
+            canonical_file = directory.joinpath('output.csv')
+        elif directory.joinpath('output.tab').is_file():
+            canonical_file = directory.joinpath('output.tab')
         else:
-            raise FileNotFoundError('\nCanonical output file not found.')
+            raise FileNotFoundError("\nCanonical output file not found.")
 
     canon = load_outputs(canonical_file,
                          transpose=transpose,
                          encoding=detect_encoding(canonical_file))
 
     # load model
-    if model_file.lower().endswith('.mdl'):
-        model = read_vensim(model_file, data_files)
-    elif model_file.lower().endswith(".xmile"):
+    if model_file.suffix.lower() == ".mdl":
+        model = read_vensim(model_file, data_files, old=old)
+    elif model_file.suffix.lower() == ".xmile":
         model = read_xmile(model_file, data_files)
+    elif model_file.suffix.lower() == ".py":
+        model = load(model_file, data_files)
     else:
-        raise ValueError('\nModelfile should be *.mdl or *.xmile')
+        raise ValueError("\nModelfile should be *.mdl, *.xmile, or *.py")
 
     # run model and return the result
 
@@ -87,8 +95,8 @@ def assert_frames_close(actual, expected, assertion="raise",
 
     assertion: str (optional)
         "raise" if an error should be raised when not able to assert
-        that two frames are close. Otherwise, it will show a warning
-        message. Default is "raise".
+        that two frames are close. If "warning", it will show a warning
+        message. If "return" it will return information. Default is "raise".
 
     verbose: bool (optional)
         If True, if any column is not close the actual and expected values
@@ -166,15 +174,17 @@ def assert_frames_close(actual, expected, assertion="raise",
         message = ""
 
         if actual_cols.difference(expected_cols):
-            columns = ["'" + col + "'" for col
-                       in actual_cols.difference(expected_cols)]
+            columns = sorted([
+                "'" + col + "'" for col
+                in actual_cols.difference(expected_cols)])
             columns = ", ".join(columns)
             message += '\nColumns ' + columns\
                        + ' from actual values not found in expected values.'
 
         if expected_cols.difference(actual_cols):
-            columns = ["'" + col + "'" for col
-                       in expected_cols.difference(actual_cols)]
+            columns = sorted([
+                "'" + col + "'" for col
+                in expected_cols.difference(actual_cols)])
             columns = ", ".join(columns)
             message += '\nColumns ' + columns\
                        + ' from expected values not found in actual values.'
@@ -190,8 +200,8 @@ def assert_frames_close(actual, expected, assertion="raise",
 
     # TODO let compare dataframes with different timestamps if "warn"
     assert np.all(np.equal(expected.index.values, actual.index.values)), \
-        'test set and actual set must share a common index' \
-        'instead found' + expected.index.values + 'vs' + actual.index.values
+        "test set and actual set must share a common index, "\
+        "instead found %s vs %s" % (expected.index.values, actual.index.values)
 
     # if for Vensim outputs where constant values are only in the first row
     _remove_constant_nan(expected)
@@ -201,13 +211,25 @@ def assert_frames_close(actual, expected, assertion="raise",
                         actual[columns],
                         **kwargs)
 
-    if c.all():
-        return
+    if c.all().all():
+        return (set(), np.nan, set()) if assertion == "return" else None
 
-    columns = np.array(columns, dtype=str)[~c.values]
+    # Get the columns that have the first different value, useful for
+    # debugging
+    false_index = c.apply(
+        lambda x: np.where(~x)[0][0] if not x.all() else np.nan)
+    index_first_false = int(np.nanmin(false_index))
+    time_first_false = c.index[index_first_false]
+    variable_first_false = sorted(
+        false_index.index[false_index == index_first_false])
+
+    columns = sorted(np.array(columns, dtype=str)[~c.all().values])
 
     assertion_details = "\nFollowing columns are not close:\n\t"\
-                        + ", ".join(columns)
+                        + ", ".join(columns) + "\n\n"\
+                        + f"First false values ({time_first_false}):\n\t"\
+                        + ", ".join(variable_first_false)
+
     if verbose:
         for col in columns:
             assertion_details += '\n\n'\
@@ -229,13 +251,15 @@ def assert_frames_close(actual, expected, assertion="raise",
 
     if assertion == "raise":
         raise AssertionError(assertion_details)
+    elif assertion == "return":
+        return (set(columns), time_first_false, set(variable_first_false))
     else:
         warnings.warn(assertion_details)
 
 
 def assert_allclose(x, y, rtol=1.e-5, atol=1.e-5):
     """
-    Asserts if all numeric values from two arrays are close.
+    Asserts if numeric values from two arrays are close.
 
     Parameters
     ----------
@@ -253,7 +277,7 @@ def assert_allclose(x, y, rtol=1.e-5, atol=1.e-5):
     None
 
     """
-    return ((abs(x - y) <= atol + rtol * abs(y)) + x.isna()*y.isna()).all()
+    return ((abs(x - y) <= atol + rtol * abs(y)) + x.isna()*y.isna())
 
 
 def _remove_constant_nan(df):
diff --git a/pysd/translation/structures/__init__.py b/pysd/translation/structures/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/pysd/translation/structures/abstract_model.py b/pysd/translation/structures/abstract_model.py
new file mode 100644
index 00000000..b295eba9
--- /dev/null
+++ b/pysd/translation/structures/abstract_model.py
@@ -0,0 +1,168 @@
+from dataclasses import dataclass
+from typing import Tuple, List, Union
+from pathlib import Path
+
+
+@dataclass
+class AbstractComponent:
+    subscripts: Tuple[str]
+    ast: object
+    type: str = "Auxiliary"
+    subtype: str = "Normal"
+
+    def __str__(self) -> str:
+        return "AbstractComponent %s\n" % (
+            "%s" % repr(list(self.subscripts)) if self.subscripts else "")
+
+    def dump(self, depth=None, indent="") -> str:
+        if depth == 0:
+            return self.__str__()
+
+        return self.__str__() + "\n" + self._str_child(depth, indent)
+
+    def _str_child(self, depth, indent) -> str:
+        return str(self.ast).replace("\t", indent).replace("\n", "\n" + indent)
+
+
+@dataclass
+class AbstractUnchangeableConstant(AbstractComponent):
+    subscripts: Tuple[str]
+    ast: object
+    type: str = "Constant"
+    subtype: str = "Unchangeable"
+
+    def __str__(self) -> str:
+        return "AbstractLookup %s\n" % (
+            "%s" % repr(list(self.subscripts)) if self.subscripts else "")
+
+
+@dataclass
+class AbstractLookup(AbstractComponent):
+    subscripts: Tuple[str]
+    ast: object
+    arguments: str = "x"
+    type: str = "Lookup"
+    subtype: str = "Hardcoded"
+
+    def __str__(self) -> str:
+        return "AbstractLookup %s\n" % (
+            "%s" % repr(list(self.subscripts)) if self.subscripts else "")
+
+
+@dataclass
+class AbstractData(AbstractComponent):
+    subscripts: Tuple[str]
+    ast: object
+    keyword: Union[str, None] = None
+    type: str = "Data"
+    subtype: str = "Normal"
+
+    def __str__(self) -> str:
+        return "AbstractData (%s)  %s\n" % (
+            self.keyword,
+            "%s" % repr(list(self.subscripts)) if self.subscripts else "")
+
+    def dump(self, depth=None, indent="") -> str:
+        if depth == 0:
+            return self.__str__()
+
+        return self.__str__() + "\n" + self._str_child(depth, indent)
+
+    def _str_child(self, depth, indent) -> str:
+        return str(self.ast).replace("\n", "\n" + indent)
+
+
+@dataclass
+class AbstractElement:
+    name: str
+    components: List[AbstractComponent]
+    units: str = ""
+    range: tuple = (None, None)
+    documentation: str = ""
+
+    def __str__(self) -> str:
+        return "AbstractElement:\t%s (%s, %s)\n%s\n" % (
+            self.name, self.units, self.range, self.documentation)
+
+    def dump(self, depth=None, indent="") -> str:
+        if depth == 0:
+            return self.__str__()
+        elif depth is not None:
+            depth -= 1
+
+        return self.__str__() + "\n" + self._str_child(depth, indent)
+
+    def _str_child(self, depth, indent) -> str:
+        return "\n".join([
+            component.dump(depth, indent) for component in self.components
+            ]).replace("\n", "\n" + indent)
+
+
+@dataclass
+class AbstractSubscriptRange:
+    name: str
+    subscripts: Tuple[str]
+    mapping: Tuple[str]
+
+    def __str__(self) -> str:
+        return "AbstractSubscriptRange:\t%s\n\t%s\n" % (
+            self.name,
+            "%s <- %s" % (self.subscripts, self.mapping)
+            if self.mapping else self.subscripts)
+
+    def dump(self, depth=None, indent="") -> str:
+        return self.__str__()
+
+
+@dataclass
+class AbstractSection:
+    name: str
+    path: Path
+    type: str  # main, macro or module
+    params: List[str]
+    returns: List[str]
+    subscripts: Tuple[AbstractSubscriptRange]
+    elements: Tuple[AbstractElement]
+    split: bool
+    views_dict: Union[dict, None]
+
+    def __str__(self) -> str:
+        return "AbstractSection (%s):\t%s (%s)\n" % (
+            self.type, self.name, self.path)
+
+    def dump(self, depth=None, indent="") -> str:
+        if depth == 0:
+            return self.__str__()
+        elif depth is not None:
+            depth -= 1
+
+        return self.__str__() + "\n" + self._str_child(depth, indent)
+
+    def _str_child(self, depth, indent) -> str:
+        return "\n".join([
+            element.dump(depth, indent) for element in self.subscripts
+            ] + [
+            element.dump(depth, indent) for element in self.elements
+            ]).replace("\n", "\n" + indent)
+
+
+@dataclass
+class AbstractModel:
+    original_path: Path
+    sections: Tuple[AbstractSection]
+
+    def __str__(self) -> str:
+        return "AbstractModel:\t%s\n" % self.original_path
+
+    def dump(self, depth=None, indent="") -> str:
+        if depth == 0:
+            return self.__str__()
+        elif depth is not None:
+            depth -= 1
+
+        return self.__str__() + "\n" + self._str_child(depth, indent)
+
+    def _str_child(self, depth, indent) -> str:
+        return "\n".join([
+            section.dump(depth, indent) for section in self.sections
+            ]).replace("\n", "\n" + indent)
diff --git a/pysd/translation/structures/components.py b/pysd/translation/structures/components.py
new file mode 100644
index 00000000..7996269a
--- /dev/null
+++ b/pysd/translation/structures/components.py
@@ -0,0 +1,272 @@
+from dataclasses import dataclass
+from typing import Union
+
+
+@dataclass
+class ArithmeticStructure:
+    operators: str
+    arguments: tuple
+
+    def __str__(self) -> str:
+        return "ArithmeticStructure:\n\t %s %s" % (
+            self.operators, self.arguments)
+
+
+@dataclass
+class LogicStructure:
+    operators: str
+    arguments: tuple
+
+    def __str__(self) -> str:
+        return "LogicStructure:\n\t %s %s" % (
+            self.operators, self.arguments)
+
+
+@dataclass
+class SubscriptsReferenceStructure:
+    subscripts: tuple
+
+    def __str__(self) -> str:
+        return "SubscriptReferenceStructure:\n\t %s" % self.subscripts
+
+
+@dataclass
+class ReferenceStructure:
+    reference: str
+    subscripts: Union[SubscriptsReferenceStructure, None] = None
+
+    def __str__(self) -> str:
+        return "ReferenceStructure:\n\t %s%s" % (
+            self.reference,
+            "\n\t" + str(self.subscripts or "").replace("\n", "\n\t"))
+
+
+@dataclass
+class CallStructure:
+    function: Union[str, object]
+    arguments: tuple
+
+    def __str__(self) -> str:
+        return "CallStructure:\n\t%s(%s)" % (
+            self.function,
+            "\n\t\t,".join([
+                "\n\t\t" + str(arg).replace("\n", "\n\t\t")
+                for arg in self.arguments
+            ]))
+
+
+@dataclass
+class GameStructure:
+    expression: object
+
+    def __str__(self) -> str:
+        return "GameStructure:\n\t%s" % self.expression
+
+
+@dataclass
+class InitialStructure:
+    initial: object
+
+    def __str__(self) -> str:
+        return "InitialStructure:\n\t%s" % (
+            self.initial)
+
+
+@dataclass
+class IntegStructure:
+    flow: object
+    initial: object
+
+    def __str__(self) -> str:
+        return "IntegStructure:\n\t%s,\n\t%s" % (
+            self.flow,
+            self.initial)
+
+
+@dataclass
+class DelayStructure:
+    input: object
+    delay_time: object
+    initial: object
+    order: float
+
+    def __str__(self) -> str:
+        return "DelayStructure (order %s):\n\t%s,\n\t%s,\n\t%s" % (
+            self.order,
+            self.input,
+            self.delay_time,
+            self.initial)
+
+
+@dataclass
+class DelayNStructure:
+    input: object
+    delay_time: object
+    initial: object
+    order: object
+
+    # DELAY N may behave different than other delays when the delay time
+    # changes during integration
+
+    def __str__(self) -> str:
+        return "DelayNStructure (order %s):\n\t%s,\n\t%s,\n\t%s" % (
+            self.order,
+            self.input,
+            self.delay_time,
+            self.initial)
+
+
+@dataclass
+class DelayFixedStructure:
+    input: object
+    delay_time: object
+    initial: object
+
+    def __str__(self) -> str:
+        return "DelayFixedStructure:\n\t%s,\n\t%s,\n\t%s" % (
+            self.input,
+            self.delay_time,
+            self.initial)
+
+
+@dataclass
+class SmoothStructure:
+    input: object
+    smooth_time: object
+    initial: object
+    order: float
+
+    def __str__(self) -> str:
+        return "SmoothStructure (order %s):\n\t%s,\n\t%s,\n\t%s" % (
+            self.order,
+            self.input,
+            self.smooth_time,
+            self.initial)
+
+
+@dataclass
+class SmoothNStructure:
+    input: object
+    smooth_time: object
+    initial: object
+    order: object
+
+    # SMOOTH N may behave different than other smooths with RungeKutta
+    # integration
+
+    def __str__(self) -> str:
+        return "SmoothNStructure (order %s):\n\t%s,\n\t%s,\n\t%s" % (
+            self.order,
+            self.input,
+            self.smooth_time,
+            self.initial)
+
+
+@dataclass
+class TrendStructure:
+    input: object
+    average_time: object
+    initial: object
+
+    def __str__(self) -> str:
+        return "TrendStructure:\n\t%s,\n\t%s,\n\t%s" % (
+            self.input,
+            self.average_time,
+            self.initial)
+
+
+@dataclass
+class ForecastStructure:
+    input: object
+    average_time: object
+    horizon: object
+
+    def __str__(self) -> str:
+        return "ForecastStructure:\n\t%s,\n\t%s,\n\t%s" % (
+            self.input,
+            self.average_time,
+            self.horizon)
+
+
+@dataclass
+class SampleIfTrueStructure:
+    condition: object
+    input: object
+    initial: object
+
+    def __str__(self) -> str:
+        return "SampleIfTrueStructure:\n\t%s,\n\t%s,\n\t%s" % (
+            self.condition,
+            self.input,
+            self.initial)
+
+
+@dataclass
+class LookupsStructure:
+    x: tuple
+    y: tuple
+    x_range: tuple
+    y_range: tuple
+
+    def __str__(self) -> str:
+        return "LookupStructure:\n\tx %s = %s\n\ty %s = %s\n" % (
+            self.x_range, self.x, self.y_range, self.y
+        )
+
+
+@dataclass
+class InlineLookupsStructure:
+    argument: None
+    lookups: LookupsStructure
+
+    def __str__(self) -> str:
+        return "InlineLookupsStructure:\n\t%s\n\t%s" % (
+            str(self.argument).replace("\n", "\n\t"),
+            str(self.lookups).replace("\n", "\n\t")
+        )
+
+
+@dataclass
+class DataStructure:
+    pass
+
+    def __str__(self) -> str:
+        return "DataStructure"
+
+
+@dataclass
+class GetLookupsStructure:
+    file: str
+    tab: str
+    x_row_or_col: str
+    cell: str
+
+    def __str__(self) -> str:
+        return "GetLookupStructure:\n\t'%s', '%s', '%s', '%s'\n" % (
+            self.file, self.tab, self.x_row_or_col, self.cell
+        )
+
+
+@dataclass
+class GetDataStructure:
+    file: str
+    tab: str
+    time_row_or_col: str
+    cell: str
+
+    def __str__(self) -> str:
+        return "GetDataStructure:\n\t'%s', '%s', '%s', '%s'\n" % (
+            self.file, self.tab, self.time_row_or_col, self.cell
+        )
+
+
+@dataclass
+class GetConstantsStructure:
+    file: str
+    tab: str
+    cell: str
+
+    def __str__(self) -> str:
+        return "GetConstantsStructure:\n\t'%s', '%s', '%s'\n" % (
+            self.file, self.tab, self.cell
+        )
diff --git a/pysd/translation/utils.py b/pysd/translation/utils.py
index 601926f6..d967c14d 100644
--- a/pysd/translation/utils.py
+++ b/pysd/translation/utils.py
@@ -487,9 +487,9 @@ def merge_nested_dicts(original_dict, dict_to_merge):
 
     Returns
     -------
-        None
-    """
+    None
 
+    """
     for k, v in dict_to_merge.items():
         if (k in original_dict and isinstance(original_dict[k], dict)
                 and isinstance(dict_to_merge[k], Mapping)):
diff --git a/pysd/translation/vensim/parsing_expr/common_grammar.peg b/pysd/translation/vensim/parsing_expr/common_grammar.peg
new file mode 100644
index 00000000..329bba92
--- /dev/null
+++ b/pysd/translation/vensim/parsing_expr/common_grammar.peg
@@ -0,0 +1,19 @@
+# Parsing Expression Grammar: common_grammar
+
+name = basic_id / escape_group
+
+# This takes care of models with Unicode variable names
+basic_id = id_start id_continue*
+
+id_start = ~r"[\w]"IU
+id_continue = id_start / ~r"[0-9\'\$\s\_]"
+
+# between quotes, either escaped quote or character that is not a quote
+escape_group = "\"" ( "\\\"" / ~r"[^\"]" )* "\""
+
+number = raw_number
+raw_number = ("+"/"-")? ~r"\d+\.?\d*([eE][+-]?\d+)?"
+string = "\'" (~r"[^\']"IU)* "\'"
+range = _ "[" ~r"[^\]]*" "]" _ ","
+
+_ = ~r"[\s\\]*"
\ No newline at end of file
diff --git a/pysd/translation/vensim/parsing_expr/components.peg b/pysd/translation/vensim/parsing_expr/components.peg
new file mode 100644
index 00000000..09afc0f0
--- /dev/null
+++ b/pysd/translation/vensim/parsing_expr/components.peg
@@ -0,0 +1,36 @@
+# Parsing Expression Grammar: components
+
+expr_type = array / final_expr / empty
+
+final_expr = logic_expr _ (logic_oper _ logic_expr)*  # logic operators (:and:, :or:)
+logic_expr = not_oper? _ comp_expr  # :not: operator
+comp_expr = add_expr _ (comp_oper _ add_expr)?  # comparison (e.g. '<', '=>')
+add_expr = prod_expr _ (add_oper _ prod_expr)*  # addition and substraction
+prod_expr = exp_expr _ (prod_oper _ exp_expr)*  # product and division
+exp_expr = neg_expr _ (exp_oper _ neg_expr)*  # exponential
+neg_expr = pre_oper? _ expr  # pre operators (-, +)
+expr = lookup_with_def / call / parens / number / reference / nan
+
+lookup_with_def = ~r"(WITH\ LOOKUP)"I _ "(" _ final_expr _ "," _ "(" _  range?  ( "(" _ raw_number _ "," _ raw_number _ ")" _ ","? _ )+ _ ")" _ ")"
+
+nan = ":NA:"
+
+arguments = ((string / final_expr) _ ","? _)*
+parens   = "(" _ final_expr _ ")"
+
+call = reference _ "(" _ arguments _ ")"
+
+reference = (name _ subscript_list) / name  # check first for subscript
+subscript_list = "[" _ (name _ "!"? _ ","? _)+ _ "]"
+
+array = (raw_number _ ("," / ";")? _)+ !~r"."  # negative lookahead for
+
+logic_oper = ~r"(%(logic_ops)s)"IU
+not_oper = ~r"(%(not_ops)s)"IU
+comp_oper = ~r"(%(comp_ops)s)"IU
+add_oper = ~r"(%(add_ops)s)"IU
+prod_oper = ~r"(%(prod_ops)s)"IU
+exp_oper = ~r"(%(exp_ops)s)"IU
+pre_oper = ~r"(%(pre_ops)s)"IU
+
+empty = "" # empty string
\ No newline at end of file
diff --git a/pysd/translation/vensim/parsing_expr/element_object.peg b/pysd/translation/vensim/parsing_expr/element_object.peg
new file mode 100644
index 00000000..58fb8310
--- /dev/null
+++ b/pysd/translation/vensim/parsing_expr/element_object.peg
@@ -0,0 +1,47 @@
+# Parsing Expression Grammar: element_object
+
+entry = unchangeable_constant / component / data_definition / subscript_definition / lookup_definition / subscript_copy
+
+# Regular component definition "="
+component = name _ subscript_component? _ "=" _ expression
+
+# Unchangeable constant definition "=="
+unchangeable_constant = name _ subscript_component? _ "==" _ expression
+
+# Lookup definition "()", uses lookahead assertion to capture whole group
+lookup_definition = name _ subscript_component? &"(" _ expression
+
+
+# Data type definition ":=" or empty with keyword
+data_definition = component_data_definition / empty_data_definition
+component_data_definition = name _ subscript_component? _ keyword? _ ":=" _ expression
+empty_data_definition = name _ subscript_component? _ keyword
+
+# Subscript ranges
+# Subcript range  regular definition ":"
+subscript_definition = name _ ":" _ (imported_subscript / literal_subscript) _ subscript_mapping_list?
+imported_subscript = basic_id _ "(" _ (string _ ","? _)* ")"
+literal_subscript = (subscript_range / subscript) _ ("," _ (subscript_range / subscript) _)*
+subscript_range = "(" _ basic_id _ "-" _ basic_id _ ")"
+
+# Subcript range  definition by copy "<->"
+subscript_copy = name _ "<->" _ name_mapping
+
+# Subscript mapping
+subscript_mapping_list = "->" _ subscript_mapping _ ("," _ subscript_mapping _)*
+subscript_mapping = (_ name_mapping _) / (_ "(" _ name_mapping _ ":" _ index_list _")" )
+name_mapping = basic_id / escape_group
+
+# Subscript except match
+subscript_list_except = ":EXCEPT:" _ '[' _ subscript_except _ ("," _ subscript_except _)* _ ']'
+subscript_except = basic_id / escape_group
+
+# Subscript match
+subscript_list = "[" _ index_list _ "]"
+index_list = subscript _ ("," _ subscript _)*
+subscript = basic_id / escape_group
+
+# Other definitions
+subscript_component = subscript_list _ subscript_list_except?
+expression = ~r".*"  # expression could be anything, at this point.
+keyword = ":" _ basic_id _ ":"
\ No newline at end of file
diff --git a/pysd/translation/vensim/parsing_expr/file_sections.peg b/pysd/translation/vensim/parsing_expr/file_sections.peg
new file mode 100644
index 00000000..2dd9c463
--- /dev/null
+++ b/pysd/translation/vensim/parsing_expr/file_sections.peg
@@ -0,0 +1,15 @@
+# Parsing Expression Grammar: file_sections
+
+# full file
+file = encoding? _ ((macro / main) _)+
+
+# macro definition
+macro = ":MACRO:" _ name _ "(" _ (name _ ","? _)+ _ ":"? _ (name _ ","? _)* _ ")" ~r".+?(?=:END OF MACRO:)" ":END OF MACRO:"
+
+# regular expressions
+main = main_part / main_end
+main_part = !":MACRO:" ~r".+(?=:MACRO:)"
+main_end = !":MACRO:" ~r".+"
+
+# encoding
+encoding = ~r"\{[^\}]*\}"
diff --git a/pysd/translation/vensim/parsing_expr/lookups.peg b/pysd/translation/vensim/parsing_expr/lookups.peg
new file mode 100644
index 00000000..ef088c9e
--- /dev/null
+++ b/pysd/translation/vensim/parsing_expr/lookups.peg
@@ -0,0 +1,7 @@
+# Parsing Expression Grammar: lookups
+
+lookup = _ "(" _ (regularLookup / excelLookup) _ ")"
+regularLookup = range? _ ( "(" _ number _ "," _ number _ ")" _ ","? _ )+
+excelLookup = ~"GET( |_)(XLS|DIRECT)( |_)LOOKUPS"I _ "(" (args _ ","? _)+ ")"
+
+args = ~r"[^,()]*"
diff --git a/pysd/translation/vensim/parsing_expr/section_elements.peg b/pysd/translation/vensim/parsing_expr/section_elements.peg
new file mode 100644
index 00000000..40250fa5
--- /dev/null
+++ b/pysd/translation/vensim/parsing_expr/section_elements.peg
@@ -0,0 +1,12 @@
+# Parsing Expression Grammar: section_elements
+
+model = (entry / section)+ sketch?
+entry = element "~" element "~" doc ("~" element)? "|"
+section = element "~" element "|"
+sketch = ~r".*"  #anything
+
+# Either an escape group, or a character that is not tilde or pipe
+element = ( escape_group / ~r"[^~|]")*
+
+# Anything other that is not a tilde or pipe
+doc = (~r"[^~|]")*
diff --git a/pysd/translation/vensim/parsing_expr/sketch.peg b/pysd/translation/vensim/parsing_expr/sketch.peg
new file mode 100644
index 00000000..b4dd0546
--- /dev/null
+++ b/pysd/translation/vensim/parsing_expr/sketch.peg
@@ -0,0 +1,58 @@
+# Parsing Expression Grammar: sketch
+
+line = var_definition / view_intro / view_title / view_definition / arrow / flow / other_objects / anything
+view_intro = ~r"\s*Sketch.*?names$" / ~r"^V300.*?ignored$"
+view_title = "*" view_name
+view_name = ~r"(?<=\*)[^\n]+$"
+view_definition = "$" color "," digit "," font_properties "|" ( ( color / ones_and_dashes ) "|")* view_code
+var_definition = var_code "," var_number "," var_name "," position "," var_box_type "," arrows_in_allowed "," hide_level "," var_face "," var_word_position "," var_thickness "," var_rest_conf ","? ( ( ones_and_dashes / color) ",")* font_properties? ","? extra_bytes?
+
+# elements used in a line defining the properties of a variable or stock
+var_name = element
+var_name = ~r"(?<=,)[^,]+(?=,)"
+var_number = digit
+var_box_type = ~r"(?<=,)\d+,\d+,\d+(?=,)" # improve this regex
+arrows_in_allowed = ~r"(?<=,)\d+(?=,)" # if this is an even number it's a shadow variable
+hide_level = digit
+var_face = digit
+var_word_position = ~r"(?<=,)\-*\d+(?=,)"
+var_thickness = digit
+var_rest_conf = digit "," ~r"\d+"
+extra_bytes = ~r"\d+,\d+,\d+,\d+,\d+,\d+" # required since Vensim 8.2.1
+arrow = arrow_code "," digit "," origin_var "," destination_var "," (digit ",")+ (ones_and_dashes ",")?  ((color ",") / ("," ~r"\d+") / (font_properties "," ~r"\d+"))* "|(" position ")|"
+
+# arrow origin and destination (this may be useful if further parsing is required)
+origin_var = digit
+destination_var = digit
+
+# flow arrows
+flow = source_or_sink_or_plot / flow_arrow
+
+# if you want to extend the parsing, these three would be a good starting point (they are followed by "anything")
+source_or_sink_or_plot = multipurpose_code "," anything
+flow_arrow =  flow_arrow_code "," anything
+other_objects = other_objects_code "," anything
+
+# fonts
+font_properties = font_name? "|" font_size "|" font_style? "|" color
+font_style =  "B" / "I" / "U" / "S" / "V"  # italics, bold, underline, etc
+font_size =  ~r"\d+"  # this needs to be made a regex to match any font
+font_name = ~r"(?<=,)[^\|\d]+(?=\|)"
+
+# x and y within the view layout. This may be useful if further parsing is required
+position = ~r"-*\d+,-*\d+"
+
+# rgb color (e.g. 255-255-255)
+color = ~r"((?<!\d|\.)([0-9]?[0-9]|1[0-9][0-9]|2[0-4][0-9]|25[0-5])(?!\d|\.) *[-] *){2}(?<!\d|\.)([0-9]?[0-9]|1[0-9][0-9]|2[0-4][0-9]|25[0-5])(?!\d|\.)"
+
+# lines that start with specific numbers (1:arrows, 11:flow_arrow, etc.)
+arrow_code = ~r"^1(?=,)"
+flow_arrow_code = ~r"^11(?=,)"
+var_code = ~r"^10(?=,)"
+multipurpose_code = ~r"^12(?=,)" # source, sink, plot, comment
+other_objects_code = ~r"^(30|31)(?=,)"
+view_code = ~r"\d+" "," digit "," digit "," ~r"\d+" # code at
+
+digit = ~r"(?<=,)\d+(?=,)" # comma separated value/s
+ones_and_dashes = ~r"\-1\-\-1\-\-1"
+anything = ~r".*"
diff --git a/pysd/translation/vensim/vensim2py.py b/pysd/translation/vensim/vensim2py.py
index 71005fbf..8bcec591 100644
--- a/pysd/translation/vensim/vensim2py.py
+++ b/pysd/translation/vensim/vensim2py.py
@@ -68,7 +68,6 @@ def get_file_sections(file_str):
     encoding = ~r"\{[^\}]*\}"
     """
     )
-
     parser = parsimonious.Grammar(file_structure_grammar)
     tree = parser.parse(file_str)
 
diff --git a/pysd/translation/vensim/vensim_element.py b/pysd/translation/vensim/vensim_element.py
new file mode 100644
index 00000000..1bdcf86e
--- /dev/null
+++ b/pysd/translation/vensim/vensim_element.py
@@ -0,0 +1,480 @@
+import re
+import warnings
+import parsimonious
+import numpy as np
+
+from ..structures.abstract_model import AbstractData, AbstractLookup, AbstractComponent, AbstractUnchangeableConstant
+
+from . import vensim_utils as vu
+from .vensim_structures import structures, operators, parsing_ops
+
+
+class Element():  # File section dataclass
+
+    def __init__(self, equation, units, documentation):
+        self.equation = equation
+        self.units, self.limits = self._parse_units(units)
+        self.documentation = documentation
+
+    def __str__(self):
+        return "Model element:\n\t%s\nunits: %s\ndocs: %s\n" % (
+            self.equation, self.units, self.documentation)
+
+    @property
+    def _verbose(self):
+        return self.__str__()
+
+    @property
+    def verbose(self):
+        print(self._verbose)
+
+    def _parse_units(self, units_str):
+        # TODO improve units parsing: move to _parse_section_elements
+        if not units_str:
+            return "", (None, None)
+
+        if units_str.endswith("]"):
+            units, lims = units_str.rsplit("[")  # types: str, str
+        else:
+            units = units_str
+            lims = "?, ?]"
+
+        lims = tuple(
+            [
+                float(x) if x.strip() != "?" else None
+                for x in lims.strip("]").split(",")
+            ]
+        )
+        return units, lims
+
+    def _parse(self):
+        tree = vu.Grammar.get("element_object").parse(self.equation)
+        self.component = ElementsComponentParser(tree).component
+        self.component.units = self.units
+        self.component.limits = self.limits
+        self.component.documentation = self.documentation
+        return self.component
+
+
+class ElementsComponentParser(parsimonious.NodeVisitor):
+    def __init__(self, ast):
+        self.mapping = []
+        self.subscripts = []
+        self.subscripts_except = []
+        self.name = None
+        self.expression = None
+        self.keyword = None
+        self.visit(ast)
+
+    def visit_subscript_definition(self, n, vc):
+        self.component = SubscriptRange(
+            self.name, self.subscripts, self.mapping)
+
+    def visit_lookup_definition(self, n, vc):
+        self.component = Lookup(
+            self.name,
+            (self.subscripts, self.subscripts_except),
+            self.expression
+        )
+
+    def visit_unchangeable_constant(self, n, vc):
+        self.component = UnchangeableConstant(
+            self.name,
+            (self.subscripts, self.subscripts_except),
+            self.expression
+        )
+
+    def visit_component(self, n, vc):
+        self.component = Component(
+            self.name,
+            (self.subscripts, self.subscripts_except),
+            self.expression
+        )
+
+    def visit_data_definition(self, n, vc):
+        self.component = Data(
+            self.name,
+            (self.subscripts, self.subscripts_except),
+            self.keyword,
+            self.expression
+        )
+
+    def visit_keyword(self, n, vc):
+        self.keyword = n.text.strip()[1:-1].lower().replace(" ", "_")
+
+    def visit_imported_subscript(self, n, vc):
+        self.subscripts = {
+            arg_name: argument.strip().strip("'")
+            for arg_name, argument
+            in zip(
+                ("file", "tab", "firstcell", "lastcell", "prefix"),
+                vc[4].split(",")
+            )
+        }
+
+    def visit_subscript_copy(self, n, vc):
+        self.component = SubscriptRange(self.name, vc[4].strip())
+
+    def visit_subscript_mapping(self, n, vc):
+
+        warnings.warn(
+            "\n Subscript mapping detected."
+            + "This feature works only in some simple cases."
+        )
+
+        if ":" in str(vc):
+            # TODO: add test for this condition
+            # Obtain subscript name and split by : and (
+            self.mapping.append(str(vc).split(":")[0].split("(")[1].strip())
+        else:
+            self.mapping.append(vc[0].strip())
+
+    def visit_subscript_range(self, n, vc):
+        subs_start = re.findall(r"\d+|\D+", vc[2].strip())
+        subs_end = re.findall(r"\d+|\D+", vc[6].strip())
+        prefix_start, num_start = "".join(subs_start[:-1]), int(subs_start[-1])
+        prefix_end, num_end = "".join(subs_end[:-1]), int(subs_end[-1])
+
+        if not prefix_start or not prefix_end:
+            raise ValueError(
+                "\nA numeric range must contain at least one letter.")
+        elif num_start >= num_end:
+            raise ValueError(
+                "\nThe number of the first subscript value must be "
+                "lower than the second subscript value in a "
+                "subscript numeric range.")
+        elif prefix_start != prefix_end:
+            raise ValueError(
+                "\nOnly matching names ending in numbers are valid.")
+
+        self.subscripts += [
+            prefix_start + str(i) for i in range(num_start, num_end + 1)
+            ]
+
+    def visit_name(self, n, vc):
+        self.name = vc[0].strip()
+
+    def visit_subscript(self, n, vc):
+        self.subscripts.append(n.text.strip())
+
+    def visit_subscript_except(self, n, vc):
+        self.subscripts_except.append(n.text.strip())
+
+    def visit_expression(self, n, vc):
+        self.expression = n.text.strip()
+
+    def generic_visit(self, n, vc):
+        return "".join(filter(None, vc)) or n.text
+
+    def visit__(self, n, vc):
+        # TODO check if necessary when finished
+        return " "
+
+
+class SubscriptRange():
+    """Subscript range definition, defined by ":" or "<->" in Vensim."""
+
+    def __init__(self, name, definition, mapping=[]):
+        self.name = name
+        self.definition = definition
+        self.mapping = mapping
+
+    def __str__(self):
+        return "\nSubscript range definition:  %s\n\t%s\n" % (
+            self.name,
+            "%s <- %s" % (self.definition, self.mapping)
+            if self.mapping else self.definition)
+
+    @property
+    def _verbose(self):
+        return self.__str__()
+
+    @property
+    def verbose(self):
+        print(self._verbose)
+
+
+class Component():
+    """Model component defined by "name = expr" in Vensim."""
+    kind = "Model component"
+
+    def __init__(self, name, subscripts, expression):
+        self.name = name
+        self.subscripts = subscripts
+        self.expression = expression
+
+    def __str__(self):
+        text = "\n%s definition: %s" % (self.kind, self.name)
+        text += "\nSubscrips: %s" % repr(self.subscripts[0])\
+            if self.subscripts[0] else ""
+        text += "  EXCEPT  %s" % repr(self.subscripts[1])\
+            if self.subscripts[1] else ""
+        text += "\n\t%s" % self._expression
+        return text
+
+    @property
+    def _expression(self):
+        if hasattr(self, "ast"):
+            return str(self.ast).replace("\n", "\n\t")
+
+        else:
+            return self.expression.replace("\n", "\n\t")
+
+    @property
+    def _verbose(self):
+        return self.__str__()
+
+    @property
+    def verbose(self):
+        print(self._verbose)
+
+    def _parse(self):
+        tree = vu.Grammar.get("components", parsing_ops).parse(self.expression)
+        self.ast = ComponentsParser(tree).translation
+        if isinstance(self.ast, structures["get_xls_lookups"]):
+            self.lookup = True
+        else:
+            self.lookup = False
+
+    def get_abstract_component(self):
+        if self.lookup:
+            return AbstractLookup(subscripts=self.subscripts, ast=self.ast)
+        else:
+            return AbstractComponent(subscripts=self.subscripts, ast=self.ast)
+
+
+class UnchangeableConstant(Component):
+    """Unchangeable constant defined by "name == expr" in Vensim."""
+    kind = "Unchangeable constant component"
+
+    def __init__(self, name, subscripts, expression):
+        super().__init__(name, subscripts, expression)
+
+    def get_abstract_component(self):
+        return AbstractUnchangeableConstant(
+            subscripts=self.subscripts, ast=self.ast)
+
+
+class Lookup(Component):
+    """Lookup variable, defined by "name(expr)" in Vensim."""
+    kind = "Lookup component"
+
+    def __init__(self, name, subscripts, expression):
+        super().__init__(name, subscripts, expression)
+
+    def _parse(self):
+        tree = vu.Grammar.get("lookups").parse(self.expression)
+        self.ast = LookupsParser(tree).translation
+
+    def get_abstract_component(self):
+        return AbstractLookup(subscripts=self.subscripts, ast=self.ast)
+
+
+class Data(Component):
+    """Data variable, defined by "name := expr" in Vensim."""
+    kind = "Data component"
+
+    def __init__(self, name, subscripts, keyword, expression):
+        super().__init__(name, subscripts, expression)
+        self.keyword = keyword
+
+    def __str__(self):
+        text = "\n%s definition: %s" % (self.kind, self.name)
+        text += "\nSubscrips: %s" % repr(self.subscripts[0])\
+            if self.subscripts[0] else ""
+        text += "  EXCEPT  %s" % repr(self.subscripts[1])\
+            if self.subscripts[1] else ""
+        text += "\nKeyword: %s" % self.keyword if self.keyword else ""
+        text += "\n\t%s" % self._expression
+        return text
+
+    def _parse(self):
+        if not self.expression:
+            # empty data vars, read from vdf file
+            self.ast = structures["data"]()
+        else:
+            super()._parse()
+
+    def get_abstract_component(self):
+        return AbstractData(
+            subscripts=self.subscripts, ast=self.ast, keyword=self.keyword)
+
+
+class LookupsParser(parsimonious.NodeVisitor):
+    def __init__(self, ast):
+        self.translation = None
+        self.visit(ast)
+
+    def visit_range(self, n, vc):
+        return n.text.strip()[:-1].replace(")-(", "),(")
+
+    def visit_regularLookup(self, n, vc):
+        if vc[0]:
+            xy_range = np.array(eval(vc[0]))
+        else:
+            xy_range = np.full((2, 2), np.nan)
+
+        values = np.array((eval(vc[2])))
+        values = values[np.argsort(values[:, 0])]
+
+        self.translation = structures["lookup"](
+            x=tuple(values[:, 0]),
+            y=tuple(values[:, 1]),
+            x_range=tuple(xy_range[:, 0]),
+            y_range=tuple(xy_range[:, 1])
+        )
+
+    def visit_excelLookup(self, n, vc):
+        arglist = vc[3].split(",")
+
+        self.translation = structures["get_xls_lookups"](
+            file=eval(arglist[0]),
+            tab=eval(arglist[1]),
+            x_row_or_col=eval(arglist[2]),
+            cell=eval(arglist[3])
+        )
+
+    def generic_visit(self, n, vc):
+        return "".join(filter(None, vc)) or n.text
+
+
+class ComponentsParser(parsimonious.NodeVisitor):
+    def __init__(self, ast):
+        self.translation = None
+        self.elements = {}
+        self.subs = None  # the subscripts if given
+        self.negatives = set()
+        self.visit(ast)
+
+    def visit_expr_type(self, n, vc):
+        self.translation = self.elements[vc[0]]
+
+    def visit_final_expr(self, n, vc):
+        return vu.split_arithmetic(
+            structures["logic"], parsing_ops["logic_ops"],
+            "".join(vc).strip(), self.elements)
+
+    def visit_logic_expr(self, n, vc):
+        id = vc[2]
+        if vc[0].lower() == ":not:":
+            id = self.add_element(structures["logic"](
+                [":NOT:"],
+                (self.elements[id],)
+                ))
+        return id
+
+    def visit_comp_expr(self, n, vc):
+        return vu.split_arithmetic(
+            structures["logic"], parsing_ops["comp_ops"],
+            "".join(vc).strip(), self.elements)
+
+    def visit_add_expr(self, n, vc):
+        return vu.split_arithmetic(
+            structures["arithmetic"], parsing_ops["add_ops"],
+            "".join(vc).strip(), self.elements)
+
+    def visit_prod_expr(self, n, vc):
+        return vu.split_arithmetic(
+            structures["arithmetic"], parsing_ops["prod_ops"],
+            "".join(vc).strip(), self.elements)
+
+    def visit_exp_expr(self, n, vc):
+        return vu.split_arithmetic(
+            structures["arithmetic"], parsing_ops["exp_ops"],
+            "".join(vc).strip(), self.elements, self.negatives)
+
+    def visit_neg_expr(self, n, vc):
+        id = vc[2]
+        if vc[0] == "-":
+            if isinstance(self.elements[id], (float, int)):
+                self.elements[id] = -self.elements[id]
+            else:
+                self.negatives.add(id)
+        return id
+
+    def visit_call(self, n, vc):
+        func = self.elements[vc[0]]
+        args = self.elements[vc[4]]
+        if func.reference in structures:
+            return self.add_element(structures[func.reference](*args))
+        else:
+            return self.add_element(structures["call"](func, args))
+
+    def visit_reference(self, n, vc):
+        id = self.add_element(structures["reference"](
+            vc[0].lower().replace(" ", "_"), self.subs))
+        self.subs = None
+        return id
+
+    def visit_range(self, n, vc):
+        return self.add_element(n.text.strip()[:-1].replace(")-(", "),("))
+
+    def visit_lookup_with_def(self, n, vc):
+        if vc[10]:
+            xy_range = np.array(eval(self.elements[vc[10]]))
+        else:
+            xy_range = np.full((2, 2), np.nan)
+
+        values = np.array((eval(vc[11])))
+        values = values[np.argsort(values[:, 0])]
+
+        lookup = structures["lookup"](
+            x=tuple(values[:, 0]),
+            y=tuple(values[:, 1]),
+            x_range=tuple(xy_range[:, 0]),
+            y_range=tuple(xy_range[:, 1])
+        )
+
+        return self.add_element(structures["with_lookup"](
+            self.elements[vc[4]], lookup))
+
+    def visit_array(self, n, vc):
+        if ";" in n.text or "," in n.text:
+            return self.add_element(np.squeeze(np.array(
+                [row.split(",") for row in n.text.strip(";").split(";")],
+                dtype=float)))
+        else:
+            return self.add_element(eval(n.text))
+
+    def visit_subscript_list(self, n, vc):
+        subs = [x.strip() for x in vc[2].split(",")]
+        self.subs = structures["subscripts_ref"](subs)
+        return ""
+
+    def visit_name(self, n, vc):
+        return n.text.strip()
+
+    def visit_expr(self, n, vc):
+        if vc[0] not in self.elements:
+            return self.add_element(eval(vc[0]))
+        else:
+            return vc[0]
+
+    def visit_string(self, n, vc):
+        return self.add_element(eval(n.text))
+
+    def visit_arguments(self, n, vc):
+        arglist = tuple(x.strip(",") for x in vc)
+        return self.add_element(tuple(
+            self.elements[arg] if arg in self.elements
+            else eval(arg) for arg in arglist))
+
+    def visit_parens(self, n, vc):
+        return vc[2]
+
+    def visit__(self, n, vc):
+        """Handles whitespace characters"""
+        return ""
+
+    def visit_nan(self, n, vc):
+        return "np.nan"
+
+    def visit_empty(self, n, vc):
+        #warnings.warn(f"Empty expression for '{element['real_name']}''.")
+        return self.add_element(None)
+
+    def generic_visit(self, n, vc):
+        return "".join(filter(None, vc)) or n.text
+
+    def add_element(self, element):
+        return vu.add_element(self.elements, element)
diff --git a/pysd/translation/vensim/vensim_section.py b/pysd/translation/vensim/vensim_section.py
new file mode 100644
index 00000000..30293662
--- /dev/null
+++ b/pysd/translation/vensim/vensim_section.py
@@ -0,0 +1,124 @@
+from typing import List, Union
+from pathlib import Path
+import parsimonious
+
+from ..structures.abstract_model import\
+    AbstractElement, AbstractSubscriptRange,  AbstractSection
+
+from . import vensim_utils as vu
+from .vensim_element import Element, SubscriptRange, Component
+
+
+class FileSection():  # File section dataclass
+
+    def __init__(self, name: str, path: Path, type: str,
+                 params: List[str], returns: List[str],
+                 content: str, split: bool, views_dict: Union[dict, None]
+                 ) -> object:
+        self.name = name
+        self.path = path
+        self.type = type
+        self.params = params
+        self.returns = returns
+        self.content = content
+        self.split = split
+        self.views_dict = views_dict
+        self.elements = None
+
+    def __str__(self):
+        return "\nFile section: %s\n" % self.name
+
+    @property
+    def _verbose(self):
+        text = self.__str__()
+        if self.elements:
+            for element in self.elements:
+                text += element._verbose
+        else:
+            text += self.content
+
+        return text
+
+    @property
+    def verbose(self):
+        print(self._verbose)
+
+    def _parse(self):
+        tree = vu.Grammar.get("section_elements").parse(self.content)
+        self.elements = SectionElementsParser(tree).entries
+        self.elements = [element._parse() for element in self.elements]
+        # split subscript from other components
+        self.subscripts = [
+            element for element in self.elements
+            if isinstance(element, SubscriptRange)
+        ]
+        self.components = [
+            element for element in self.elements
+            if isinstance(element, Component)
+        ]
+        # reorder element list for better printing
+        self.elements = self.subscripts + self.components
+
+        [component._parse() for component in self.components]
+
+    def get_abstract_section(self):
+        return AbstractSection(
+            name=self.name,
+            path=self.path,
+            type=self.type,
+            params=self.params,
+            returns=self.returns,
+            subscripts=self.solve_subscripts(),
+            elements=self.merge_components(),
+            split=self.split,
+            views_dict=self.views_dict
+        )
+
+    def solve_subscripts(self):
+        return [AbstractSubscriptRange(
+            name=subs_range.name,
+            subscripts=subs_range.definition,
+            mapping=subs_range.mapping
+        ) for subs_range in self.subscripts]
+
+    def merge_components(self):
+        merged = {}
+        for component in self.components:
+            name = component.name.lower().replace(" ", "_")
+            if name not in merged:
+                merged[name] = AbstractElement(
+                    name=component.name,
+                    components=[])
+
+            if component.units:
+                merged[name].units = component.units
+            if component.limits[0] is not None\
+              or component.limits[1] is not None:
+                merged[name].range = component.limits
+            if component.documentation:
+                merged[name].documentation = component.documentation
+
+            merged[name].components.append(component.get_abstract_component())
+
+
+
+        return list(merged.values())
+
+
+class SectionElementsParser(parsimonious.NodeVisitor):
+    # TODO include units parsing
+    def __init__(self, ast):
+        self.entries = []
+        self.visit(ast)
+
+    def visit_entry(self, n, vc):
+        self.entries.append(
+            Element(
+                equation=vc[0].strip(),
+                units=vc[2].strip(),
+                documentation=vc[4].strip(),
+            )
+        )
+
+    def generic_visit(self, n, vc):
+        return "".join(filter(None, vc)) or n.text or ""
diff --git a/pysd/translation/vensim/vensim_structures.py b/pysd/translation/vensim/vensim_structures.py
new file mode 100644
index 00000000..5341358c
--- /dev/null
+++ b/pysd/translation/vensim/vensim_structures.py
@@ -0,0 +1,54 @@
+import re
+from ..structures import components as cs
+
+
+structures = {
+    "reference":  cs.ReferenceStructure,
+    "subscripts_ref": cs.SubscriptsReferenceStructure,
+    "arithmetic": cs.ArithmeticStructure,
+    "logic": cs.LogicStructure,
+    "with_lookup": cs.InlineLookupsStructure,
+    "call": cs.CallStructure,
+    "game": cs.GameStructure,
+    "get_xls_lookups": cs.GetLookupsStructure,
+    "get_direct_lookups": cs.GetLookupsStructure,
+    "get_xls_data": cs.GetDataStructure,
+    "get_direct_data": cs.GetDataStructure,
+    "get_xls_constants": cs.GetConstantsStructure,
+    "get_direct_constants": cs.GetConstantsStructure,
+    "initial": cs.InitialStructure,
+    "integ": cs.IntegStructure,
+    "delay1": lambda x, y: cs.DelayStructure(x, y, x, 1),
+    "delay1i": lambda x, y, z: cs.DelayStructure(x, y, z, 1),
+    "delay3": lambda x, y: cs.DelayStructure(x, y, x, 3),
+    "delay3i": lambda x, y, z: cs.DelayStructure(x, y, z, 3),
+    "delay_n": cs.DelayNStructure,
+    "delay_fixed": cs.DelayFixedStructure,
+    "smooth": lambda x, y: cs.SmoothStructure(x, y, x, 1),
+    "smoothi": lambda x, y, z: cs.SmoothStructure(x, y, z, 1),
+    "smooth3": lambda x, y: cs.SmoothStructure(x, y, x, 3),
+    "smooth3i": lambda x, y, z: cs.SmoothStructure(x, y, z, 3),
+    "smooth_n": cs.SmoothNStructure,
+    "trend": cs.TrendStructure,
+    "forecast": cs.ForecastStructure,
+    "sample_if_true": cs.SampleIfTrueStructure,
+    "lookup": cs.LookupsStructure,
+    "data": cs.DataStructure
+}
+
+
+operators = {
+    "logic_ops": [":AND:", ":OR:"],
+    "not_ops": [":NOT:"],
+    "comp_ops": ["=", "<>", "<=", "<", ">=", ">"],
+    "add_ops": ["+", "-"],
+    "prod_ops": ["*", "/"],
+    "exp_ops": ["^"],
+    "pre_ops": ["+", "-"]
+}
+
+
+parsing_ops = {
+    key: "|".join(re.escape(x) for x in values)
+    for key, values in operators.items()
+}
diff --git a/pysd/translation/vensim/vensim_utils.py b/pysd/translation/vensim/vensim_utils.py
new file mode 100644
index 00000000..ca369555
--- /dev/null
+++ b/pysd/translation/vensim/vensim_utils.py
@@ -0,0 +1,115 @@
+import re
+import warnings
+import uuid
+
+import parsimonious
+from typing import Dict
+from pathlib import Path
+from chardet import detect
+
+
+class Grammar():
+    _common_grammar = None
+    _grammar_path: Path = Path(__file__).parent.joinpath("parsing_expr")
+    _grammar: Dict = {}
+
+    @classmethod
+    def get(cls, grammar: str, subs: dict = {}) -> parsimonious.Grammar:
+        """Get parsimonious grammar for parsing"""
+        if grammar not in cls._grammar:
+            # include grammar in the class singleton
+            cls._grammar[grammar] = parsimonious.Grammar(
+                cls._read_grammar(grammar) % subs
+            )
+
+        return cls._grammar[grammar]
+
+    @classmethod
+    def _read_grammar(cls, grammar: str) -> str:
+        """Read grammar from a file and include common grammar"""
+        with cls._gpath(grammar).open(encoding="ascii") as gfile:
+            source_grammar: str = gfile.read()
+
+        return cls._include_common_grammar(source_grammar)
+
+    @classmethod
+    def _include_common_grammar(cls, source_grammar: str) -> str:
+        """Include common grammar"""
+        if not cls._common_grammar:
+            with cls._gpath("common_grammar").open(encoding="ascii") as gfile:
+                cls._common_grammar: str = gfile.read()
+
+        return r"{source_grammar}{common_grammar}".format(
+            source_grammar=source_grammar, common_grammar=cls._common_grammar
+        )
+
+    @classmethod
+    def _gpath(cls, grammar: str) -> Path:
+        """Get the grammar file path"""
+        return cls._grammar_path.joinpath(grammar).with_suffix(".peg")
+
+    @classmethod
+    def clean(cls) -> None:
+        """Clean the saved grammars (used for debugging)"""
+        cls._common_grammar = None
+        cls._grammar: Dict = {}
+
+
+def _detect_encoding_from_file(mdl_file: Path) -> str:
+    """Detect and return the encoding from a Vensim file"""
+    try:
+        with mdl_file.open("rb") as in_file:
+            f_line: bytes = in_file.readline()
+        f_line: str = f_line.decode(detect(f_line)['encoding'])
+        return re.search(r"(?<={)(.*)(?=})", f_line).group()
+    except (AttributeError, UnicodeDecodeError):
+        warnings.warn(
+            "No encoding specified or detected to translate the model "
+            "file. 'UTF-8' encoding will be used.")
+        return "UTF-8"
+
+
+def split_arithmetic(structure: object, parsing_ops: dict,
+                     expression: str, elements: dict,
+                     negatives: set = set()) -> object:
+    pattern = re.compile(parsing_ops)
+    parts = pattern.split(expression)
+    ops = pattern.findall(expression)
+    if not ops:
+        if parts[0] in negatives:
+            negatives.remove(parts[0])
+            return add_element(
+                elements,
+                structure(["negative"], (elements[parts[0]],)))
+        else:
+            return expression
+    else:
+        if not negatives:
+            return add_element(
+                elements,
+                structure(
+                    ops,
+                    tuple([elements[id] for id in parts])))
+        else:
+            # manage negative expressions
+            current_id = parts.pop()
+            current = elements[current_id]
+            if current_id in negatives:
+                negatives.remove(current_id)
+                current = structure(["negative"], (current,))
+            while ops:
+                current_id = parts.pop()
+                current = structure(
+                    [ops.pop()],
+                    (elements[current_id], current))
+                if current_id in negatives:
+                    negatives.remove(current_id)
+                    current = structure(["negative"], (current,))
+
+            return add_element(elements, current)
+
+
+def add_element(elements: dict, element: object) -> str:
+    id = uuid.uuid4().hex
+    elements[id] = element
+    return id
diff --git a/pysd/translation/vensim/vensin_file.py b/pysd/translation/vensim/vensin_file.py
new file mode 100644
index 00000000..fbe63f31
--- /dev/null
+++ b/pysd/translation/vensim/vensin_file.py
@@ -0,0 +1,276 @@
+import re
+from pathlib import Path
+import warnings
+import parsimonious
+from collections.abc import Mapping
+
+from ..structures.abstract_model import AbstractModel
+
+from . import vensim_utils as vu
+from .vensim_section import FileSection
+
+
+class VensimFile():
+    """
+    Create a VensimFile object which allows parsing a mdl file.
+
+    Parameters
+    ----------
+    mdl_path: str or pathlib.Path
+        Path to the Vensim model.
+
+    encoding: str or None (optional)
+        Encoding of the source model file. If None, the encoding will be
+        read from the model, if the encoding is not defined in the model
+        file it will be set to 'UTF-8'. Default is None.
+
+    """
+    def __init__(self, mdl_path, encoding=None):
+        self.mdl_path = Path(mdl_path)
+        self.root_path = self.mdl_path.parent
+        self.model_text = self._read(encoding)
+        self.sketch = ""
+        self.view_elements = None
+        self._split_sketch()
+
+    def __str__(self):
+        return "\nVensim model file, loaded from:\n\t%s\n" % self.mdl_path
+
+    @property
+    def _verbose(self):
+        text = self.__str__()
+        for section in self.sections:
+            text += section._verbose
+
+        return text
+
+    @property
+    def verbose(self):
+        print(self._verbose)
+
+    def _read(self, encoding):
+        """Read a Vensim file and assign its content to self.model_text"""
+        # check for model extension
+        if self.mdl_path.suffix.lower() != ".mdl":
+            raise ValueError(
+                "The file to translate, '%s' " % self.mdl_path
+                + "is not a vensim model. It must end with mdl extension."
+            )
+
+        if encoding is None:
+            encoding = vu._detect_encoding_from_file(self.mdl_path)
+
+        with self.mdl_path.open("r", encoding=encoding,
+                                errors="ignore") as in_file:
+            model_text = in_file.read()
+
+        return model_text
+
+    def _split_sketch(self):
+        """Split model from the sketch"""
+        try:
+            split_model = self.model_text.split("\\\\\\---///", 1)
+            self.model_text = self._clean(split_model[0])
+            # remove plots section, if it exists
+            self.sketch = split_model[1].split("///---\\\\\\")[0]
+        except LookupError:
+            pass
+
+    def _clean(self, text):
+        return re.sub(r"[\n\t\s]+", " ", re.sub(r"\\\n\t", " ", text))
+
+    def parse(self):
+        tree = vu.Grammar.get("file_sections").parse(self.model_text)
+        self.sections = FileSectionsParser(tree).entries
+        self.sections[0].path = self.mdl_path.with_suffix(".py")
+        for section in self.sections[1:]:
+            section.path = self.mdl_path.parent.joinpath(
+                self.clean_file_names(section.name)[0]
+                ).with_suffix(".py")
+        # TODO modify names and paths of macros
+        for section in self.sections:
+            section._parse()
+
+    def parse_sketch(self, subview_sep):
+        if self.sketch:
+            sketch = list(map(
+                lambda x: x.strip(),
+                self.sketch.split("\\\\\\---/// ")
+            ))
+        else:
+            warnings.warn(
+                "No sketch detected. The model will be built in a "
+                "single file.")
+            return None
+
+        grammar = vu.Grammar.get("sketch")
+        view_elements = {}
+        for module in sketch:
+            for sketch_line in module.split("\n"):
+                # parsed line could have information about new view name
+                # or of a variable inside a view
+                parsed = SketchParser(grammar.parse(sketch_line))
+
+                if parsed.view_name:
+                    view_name = parsed.view_name
+                    view_elements[view_name] = set()
+
+                elif parsed.variable_name:
+                    view_elements[view_name].add(parsed.variable_name)
+
+        # removes views that do not include any variable in them
+        non_empty_views = {
+            key: value for key, value in view_elements.items() if value
+        }
+
+        # split into subviews, if subview_sep is provided
+        views_dict = {}
+
+        if len(non_empty_views) == 1:
+            warnings.warn(
+                "Only a single view with no subviews was detected. The model"
+                " will be built in a single file.")
+            return
+        elif subview_sep and any(
+          sep in view for sep in subview_sep for view in non_empty_views):
+            escaped_separators = list(map(lambda x: re.escape(x), subview_sep))
+            for full_name, values in non_empty_views.items():
+                # split the full view name using the separator and make the
+                # individual parts safe file or directory names
+                clean_view_parts = self.clean_file_names(
+                    *re.split("|".join(escaped_separators), full_name))
+                # creating a nested dict for each view.subview
+                # (e.g. {view_name: {subview_name: [values]}})
+                nested_dict = values
+
+                for item in reversed(clean_view_parts):
+                    nested_dict = {item: nested_dict}
+                # merging the new nested_dict into the views_dict, preserving
+                # repeated keys
+                self.merge_nested_dicts(views_dict, nested_dict)
+        else:
+            # view names do not have separators or separator characters
+            # not provided
+
+            if subview_sep and not any(
+              sep in view for sep in subview_sep for view in non_empty_views):
+                warnings.warn(
+                    "The given subview separators were not matched in "
+                    "any view name.")
+
+            for view_name, elements in non_empty_views.items():
+                views_dict[self.clean_file_names(view_name)[0]] = elements
+
+        self.sections[0].split = True
+        self.sections[0].views_dict = views_dict
+
+    def get_abstract_model(self):
+        return AbstractModel(
+            original_path=self.mdl_path,
+            sections=tuple(section.get_abstract_section()
+                           for section in self.sections))
+
+    @staticmethod
+    def clean_file_names(*args):
+        """
+        Removes special characters and makes clean file names.
+
+        Parameters
+        ----------
+        *args: tuple
+            Any number of strings to to clean.
+
+        Returns
+        -------
+        clean: list
+            List containing the clean strings.
+
+        """
+        return [
+            re.sub(
+                r"[\W]+", "",
+                name.replace(" ", "_")
+            ).lstrip("0123456789")
+            for name in args]
+
+    def merge_nested_dicts(self, original_dict, dict_to_merge):
+        """
+        Merge dictionaries recursively, preserving common keys.
+
+        Parameters
+        ----------
+        original_dict: dict
+            Dictionary onto which the merge is executed.
+
+        dict_to_merge: dict
+            Dictionary to be merged to the original_dict.
+
+        Returns
+        -------
+        None
+
+        """
+        for key, value in dict_to_merge.items():
+            if (key in original_dict and isinstance(original_dict[key], dict)
+                    and isinstance(value, Mapping)):
+                self.merge_nested_dicts(original_dict[key], value)
+            else:
+                original_dict[key] = value
+
+
+class FileSectionsParser(parsimonious.NodeVisitor):
+    """Parse file sections"""
+    def __init__(self, ast):
+        self.entries = [None]
+        self.visit(ast)
+
+    def visit_main(self, n, vc):
+        # main will be always stored as the first entry
+        if self.entries[0] is None:
+            self.entries[0] = FileSection(
+                name="__main__",
+                path=Path("."),
+                type="main",
+                params=[],
+                returns=[],
+                content=n.text.strip(),
+                split=False,
+                views_dict=None
+            )
+        else:
+            # this is needed when macro parts are in the middle of the file
+            self.entries[0].content += n.text.strip()
+
+    def visit_macro(self, n, vc):
+        self.entries.append(
+            FileSection(
+                name=vc[2].strip().lower().replace(" ", "_"),
+                path=Path("."),
+                type="macro",
+                params=[x.strip() for x in vc[6].split(",")] if vc[6] else [],
+                returns=[x.strip() for x in vc[10].split(",")] if vc[10] else [],
+                content=vc[13].strip(),
+                split=False,
+                views_dict=None
+            )
+        )
+
+    def generic_visit(self, n, vc):
+        return "".join(filter(None, vc)) or n.text or ""
+
+
+class SketchParser(parsimonious.NodeVisitor):
+    def __init__(self, ast):
+        self.variable_name = None
+        self.view_name = None
+        self.visit(ast)
+
+    def visit_view_name(self, n, vc):
+        self.view_name = n.text.lower()
+
+    def visit_var_definition(self, n, vc):
+        if int(vc[10]) % 2 != 0:  # not a shadow variable
+            self.variable_name = vc[4].replace(" ", "_").lower()
+
+    def generic_visit(self, n, vc):
+        return "".join(filter(None, vc)) or n.text.strip() or ""
diff --git a/pysd/translation/xmile/xmile2py.py b/pysd/translation/xmile/xmile2py.py
index e714d4ed..e2b973fd 100644
--- a/pysd/translation/xmile/xmile2py.py
+++ b/pysd/translation/xmile/xmile2py.py
@@ -21,6 +21,8 @@ def translate_xmile(xmile_file):
     Functionality is currently limited.
 
     """
+    if not isinstance(xmile_file, str):
+        xmile_file = str(xmile_file)
     # process xml file
     xml_parser = etree.XMLParser(encoding="utf-8", recover=True)
     root = etree.parse(xmile_file, parser=xml_parser).getroot()
diff --git a/tests/conftest.py b/tests/conftest.py
index a6be8df1..b4fd902b 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -8,6 +8,12 @@ def _root():
     return Path(__file__).parent.resolve()
 
 
+@pytest.fixture(scope="session")
+def _test_models(_root):
+    # test-models directory
+    return _root.joinpath("test-models/tests")
+
+
 @pytest.fixture(scope="class")
 def shared_tmpdir(tmpdir_factory):
     # shared temporary directory for each class
diff --git a/tests/integration_test_factory.py b/tests/integration_test_factory.py
index ea85aeab..7ad08458 100644
--- a/tests/integration_test_factory.py
+++ b/tests/integration_test_factory.py
@@ -1,53 +1,32 @@
-from __future__ import print_function
+import os.path
+import glob
 
-run = False
-
-if run:
-
-
-    import os.path
-    import textwrap
-    import glob
-    from pysd import utils
-
-    test_dir = 'test-models/'
-    vensim_test_files = glob.glob(test_dir+'tests/*/*.mdl')
+if False:
+    vensim_test_files = glob.glob("test-models/tests/*/*.mdl")
+    vensim_test_files.sort()
 
     tests = []
     for file_path in vensim_test_files:
-        (path, file_name) = os.path.split(file_path)
-        (name, ext) = os.path.splitext(file_name)
-
-        test_name = utils.make_python_identifier(path.split('/')[-1])[0]
+        path, file_name = os.path.split(file_path)
+        folder = path.split("/")[-1]
 
         test_func_string = """
-            def test_%(test_name)s(self):
-                output, canon = runner('%(file_path)s')
-                assert_frames_close(output, canon, rtol=rtol)
-            """ % {
-            'file_path': file_path,
-            'test_name': test_name,
+    "%(test_name)s": {
+        "folder": "%(folder)s",
+        "file": "%(file_name)s"
+    },""" % {
+            "folder": folder,
+            "test_name": folder,
+            "file_name": file_name,
         }
         tests.append(test_func_string)
 
-    file_string = textwrap.dedent('''
-        """
-        Note that this file is autogenerated by `integration_test_factory.py`
-        and changes are likely to be overwritten.
-        """
-
-        import unittest
-        from pysd.tools.benchmarking import runner, assert_frames_close
-
-        rtol = .05
-
-
-        class TestIntegrationExamples(unittest.TestCase):
-        %(tests)s
-
-        ''' % {'tests': ''.join(tests)})
+    file_string = """
+    vensim_test = {%(tests)s
+    }
+    """ % {"tests": "".join(tests)}
 
-    with open('integration_test_pysd.py', 'w', encoding='UTF-8') as ofile:
+    with open("test_factory_result.py", "w", encoding="UTF-8") as ofile:
         ofile.write(file_string)
 
-    print('generated %i integration tests' % len(tests))
+    print("Generated %i integration tests" % len(tests))
diff --git a/tests/integration_test_vensim_pathway.py b/tests/integration_test_vensim_pathway.py
index e2efaff2..fb564e1c 100644
--- a/tests/integration_test_vensim_pathway.py
+++ b/tests/integration_test_vensim_pathway.py
@@ -17,119 +17,119 @@
 class TestIntegrationExamples(unittest.TestCase):
 
     def test_abs(self):
-        output, canon = runner(test_models + '/abs/test_abs.mdl')
+        output, canon = runner(test_models + '/abs/test_abs.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_active_initial(self):
-        output, canon = runner(test_models + '/active_initial/test_active_initial.mdl')
+        output, canon = runner(test_models + '/active_initial/test_active_initial.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_active_initial_circular(self):
-        output, canon = runner(test_models + '/active_initial_circular/test_active_initial_circular.mdl')
+        output, canon = runner(test_models + '/active_initial_circular/test_active_initial_circular.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_arguments(self):
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
-            output, canon = runner(test_models + '/arguments/test_arguments.mdl')
+            output, canon = runner(test_models + '/arguments/test_arguments.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_array_with_line_break(self):
-        output, canon = runner(test_models + '/array_with_line_break/test_array_with_line_break.mdl')
+        output, canon = runner(test_models + '/array_with_line_break/test_array_with_line_break.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_builtin_max(self):
-        output, canon = runner(test_models + '/builtin_max/builtin_max.mdl')
+        output, canon = runner(test_models + '/builtin_max/builtin_max.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_builtin_min(self):
-        output, canon = runner(test_models + '/builtin_min/builtin_min.mdl')
+        output, canon = runner(test_models + '/builtin_min/builtin_min.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_chained_initialization(self):
-        output, canon = runner(test_models + '/chained_initialization/test_chained_initialization.mdl')
+        output, canon = runner(test_models + '/chained_initialization/test_chained_initialization.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     @unittest.skip("Working on it")
     def test_conditional_subscripts(self):
-        output, canon = runner(test_models + '/conditional_subscripts/test_conditional_subscripts.mdl')
+        output, canon = runner(test_models + '/conditional_subscripts/test_conditional_subscripts.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_control_vars(self):
-        output, canon = runner(test_models + '/control_vars/test_control_vars.mdl')
+        output, canon = runner(test_models + '/control_vars/test_control_vars.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_constant_expressions(self):
-        output, canon = runner(test_models + '/constant_expressions/test_constant_expressions.mdl')
+        output, canon = runner(test_models + '/constant_expressions/test_constant_expressions.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_data_from_other_model(self):
         output, canon = runner(
             test_models + '/data_from_other_model/test_data_from_other_model.mdl',
-            data_files=test_models + '/data_from_other_model/data.tab')
+            data_files=test_models + '/data_from_other_model/data.tab', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_delay_fixed(self):
         # issue https://github.com/JamesPHoughton/pysd/issues/147
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
-            output, canon = runner(test_models + '/delay_fixed/test_delay_fixed.mdl')
+            output, canon = runner(test_models + '/delay_fixed/test_delay_fixed.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_delay_numeric_error(self):
         # issue https://github.com/JamesPHoughton/pysd/issues/225
-        output, canon = runner(test_models + '/delay_numeric_error/test_delay_numeric_error.mdl')
+        output, canon = runner(test_models + '/delay_numeric_error/test_delay_numeric_error.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_delay_parentheses(self):
-        output, canon = runner(test_models + '/delay_parentheses/test_delay_parentheses.mdl')
+        output, canon = runner(test_models + '/delay_parentheses/test_delay_parentheses.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_delay_pipeline(self):
         # issue https://github.com/JamesPHoughton/pysd/issues/147
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
-            output, canon = runner(test_models + '/delay_pipeline/test_pipeline_delays.mdl')
+            output, canon = runner(test_models + '/delay_pipeline/test_pipeline_delays.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_delays(self):
         # issue https://github.com/JamesPHoughton/pysd/issues/147
-        output, canon = runner(test_models + '/delays/test_delays.mdl')
+        output, canon = runner(test_models + '/delays/test_delays.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_dynamic_final_time(self):
         # issue https://github.com/JamesPHoughton/pysd/issues/278
-        output, canon = runner(test_models + '/dynamic_final_time/test_dynamic_final_time.mdl')
+        output, canon = runner(test_models + '/dynamic_final_time/test_dynamic_final_time.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_euler_step_vs_saveper(self):
-        output, canon = runner(test_models + '/euler_step_vs_saveper/test_euler_step_vs_saveper.mdl')
+        output, canon = runner(test_models + '/euler_step_vs_saveper/test_euler_step_vs_saveper.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_exp(self):
-        output, canon = runner(test_models + '/exp/test_exp.mdl')
+        output, canon = runner(test_models + '/exp/test_exp.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_exponentiation(self):
-        output, canon = runner(test_models + '/exponentiation/exponentiation.mdl')
+        output, canon = runner(test_models + '/exponentiation/exponentiation.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_forecast(self):
-        output, canon = runner(test_models + '/forecast/test_forecast.mdl')
+        output, canon = runner(test_models + '/forecast/test_forecast.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_function_capitalization(self):
-        output, canon = runner(test_models + '/function_capitalization/test_function_capitalization.mdl')
+        output, canon = runner(test_models + '/function_capitalization/test_function_capitalization.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_game(self):
-        output, canon = runner(test_models + '/game/test_game.mdl')
+        output, canon = runner(test_models + '/game/test_game.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_get_constants_subrange(self):
         output, canon = runner(
           test_models + '/get_constants_subranges/'
-          + 'test_get_constants_subranges.mdl'
+          + 'test_get_constants_subranges.mdl', old=True
         )
         assert_frames_close(output, canon, rtol=rtol)
 
@@ -142,7 +142,7 @@ def test_get_data_args_3d_xls(self):
         """
         output, canon = runner(
           test_models + '/get_data_args_3d_xls/'
-          + 'test_get_data_args_3d_xls.mdl'
+          + 'test_get_data_args_3d_xls.mdl', old=True
         )
         assert_frames_close(output, canon, rtol=rtol)
 
@@ -155,7 +155,7 @@ def test_get_lookups_data_3d_xls(self):
         """
         output, canon = runner(
           test_models + '/get_lookups_data_3d_xls/'
-          + 'test_get_lookups_data_3d_xls.mdl'
+          + 'test_get_lookups_data_3d_xls.mdl', old=True
         )
         assert_frames_close(output, canon, rtol=rtol)
 
@@ -164,14 +164,14 @@ def test_get_lookups_subscripted_args(self):
             warnings.simplefilter("ignore")
             output, canon = runner(
               test_models + '/get_lookups_subscripted_args/'
-              + 'test_get_lookups_subscripted_args.mdl'
+              + 'test_get_lookups_subscripted_args.mdl', old=True
             )
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_get_lookups_subset(self):
         output, canon = runner(
           test_models + '/get_lookups_subset/'
-          + 'test_get_lookups_subset.mdl'
+          + 'test_get_lookups_subset.mdl', old=True
         )
         assert_frames_close(output, canon, rtol=rtol)
 
@@ -180,7 +180,7 @@ def test_get_with_missing_values_xlsx(self):
             warnings.simplefilter("ignore")
             output, canon = runner(
                 test_models + '/get_with_missing_values_xlsx/'
-                + 'test_get_with_missing_values_xlsx.mdl'
+                + 'test_get_with_missing_values_xlsx.mdl', old=True
             )
 
         assert_frames_close(output, canon, rtol=rtol)
@@ -188,7 +188,7 @@ def test_get_with_missing_values_xlsx(self):
     def test_get_mixed_definitions(self):
         output, canon = runner(
           test_models + '/get_mixed_definitions/'
-          + 'test_get_mixed_definitions.mdl'
+          + 'test_get_mixed_definitions.mdl', old=True
         )
         assert_frames_close(output, canon, rtol=rtol)
 
@@ -201,343 +201,343 @@ def test_get_subscript_3d_arrays_xls(self):
         """
         output, canon = runner(
           test_models + '/get_subscript_3d_arrays_xls/'
-          + 'test_get_subscript_3d_arrays_xls.mdl'
+          + 'test_get_subscript_3d_arrays_xls.mdl', old=True
         )
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_get_xls_cellrange(self):
         output, canon = runner(
           test_models + '/get_xls_cellrange/'
-          + 'test_get_xls_cellrange.mdl'
+          + 'test_get_xls_cellrange.mdl', old=True
         )
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_if_stmt(self):
-        output, canon = runner(test_models + '/if_stmt/if_stmt.mdl')
+        output, canon = runner(test_models + '/if_stmt/if_stmt.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_initial_function(self):
-        output, canon = runner(test_models + '/initial_function/test_initial.mdl')
+        output, canon = runner(test_models + '/initial_function/test_initial.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_input_functions(self):
-        output, canon = runner(test_models + '/input_functions/test_inputs.mdl')
+        output, canon = runner(test_models + '/input_functions/test_inputs.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscripted_round(self):
-        output, canon = runner(test_models + '/subscripted_round/test_subscripted_round.mdl')
+        output, canon = runner(test_models + '/subscripted_round/test_subscripted_round.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_invert_matrix(self):
-        output, canon = runner(test_models + '/invert_matrix/test_invert_matrix.mdl')
+        output, canon = runner(test_models + '/invert_matrix/test_invert_matrix.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_limits(self):
-        output, canon = runner(test_models + '/limits/test_limits.mdl')
+        output, canon = runner(test_models + '/limits/test_limits.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_line_breaks(self):
-        output, canon = runner(test_models + '/line_breaks/test_line_breaks.mdl')
+        output, canon = runner(test_models + '/line_breaks/test_line_breaks.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_line_continuation(self):
-        output, canon = runner(test_models + '/line_continuation/test_line_continuation.mdl')
+        output, canon = runner(test_models + '/line_continuation/test_line_continuation.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_ln(self):
-        output, canon = runner(test_models + '/ln/test_ln.mdl')
+        output, canon = runner(test_models + '/ln/test_ln.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_log(self):
-        output, canon = runner(test_models + '/log/test_log.mdl')
+        output, canon = runner(test_models + '/log/test_log.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_logicals(self):
-        output, canon = runner(test_models + '/logicals/test_logicals.mdl')
+        output, canon = runner(test_models + '/logicals/test_logicals.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_lookups(self):
-        output, canon = runner(test_models + '/lookups/test_lookups.mdl')
+        output, canon = runner(test_models + '/lookups/test_lookups.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_lookups_without_range(self):
-        output, canon = runner(test_models + '/lookups_without_range/test_lookups_without_range.mdl')
+        output, canon = runner(test_models + '/lookups_without_range/test_lookups_without_range.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_lookups_funcnames(self):
-        output, canon = runner(test_models + '/lookups_funcnames/test_lookups_funcnames.mdl')
+        output, canon = runner(test_models + '/lookups_funcnames/test_lookups_funcnames.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_lookups_inline(self):
-        output, canon = runner(test_models + '/lookups_inline/test_lookups_inline.mdl')
+        output, canon = runner(test_models + '/lookups_inline/test_lookups_inline.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_lookups_inline_bounded(self):
-        output, canon = runner(test_models + '/lookups_inline_bounded/test_lookups_inline_bounded.mdl')
+        output, canon = runner(test_models + '/lookups_inline_bounded/test_lookups_inline_bounded.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_lookups_with_expr(self):
-        output, canon = runner(test_models + '/lookups_with_expr/test_lookups_with_expr.mdl')
+        output, canon = runner(test_models + '/lookups_with_expr/test_lookups_with_expr.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_macro_cross_reference(self):
-        output, canon = runner(test_models + '/macro_cross_reference/test_macro_cross_reference.mdl')
+        output, canon = runner(test_models + '/macro_cross_reference/test_macro_cross_reference.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_macro_expression(self):
-        output, canon = runner(test_models + '/macro_expression/test_macro_expression.mdl')
+        output, canon = runner(test_models + '/macro_expression/test_macro_expression.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_macro_multi_expression(self):
-        output, canon = runner(test_models + '/macro_multi_expression/test_macro_multi_expression.mdl')
+        output, canon = runner(test_models + '/macro_multi_expression/test_macro_multi_expression.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_macro_multi_macros(self):
-        output, canon = runner(test_models + '/macro_multi_macros/test_macro_multi_macros.mdl')
+        output, canon = runner(test_models + '/macro_multi_macros/test_macro_multi_macros.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     @unittest.skip('working')
     def test_macro_output(self):
-        output, canon = runner(test_models + '/macro_output/test_macro_output.mdl')
+        output, canon = runner(test_models + '/macro_output/test_macro_output.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_macro_stock(self):
-        output, canon = runner(test_models + '/macro_stock/test_macro_stock.mdl')
+        output, canon = runner(test_models + '/macro_stock/test_macro_stock.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
-    @unittest.skip('do we need this?')
+    @unittest.skip("Working on it")
     def test_macro_trailing_definition(self):
-        output, canon = runner(test_models + '/macro_trailing_definition/test_macro_trailing_definition.mdl')
+        output, canon = runner(test_models + '/macro_trailing_definition/test_macro_trailing_definition.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_model_doc(self):
-        output, canon = runner(test_models + '/model_doc/model_doc.mdl')
+        output, canon = runner(test_models + '/model_doc/model_doc.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_nested_functions(self):
-        output, canon = runner(test_models + '/nested_functions/test_nested_functions.mdl')
+        output, canon = runner(test_models + '/nested_functions/test_nested_functions.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_number_handling(self):
-        output, canon = runner(test_models + '/number_handling/test_number_handling.mdl')
+        output, canon = runner(test_models + '/number_handling/test_number_handling.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_parentheses(self):
-        output, canon = runner(test_models + '/parentheses/test_parens.mdl')
+        output, canon = runner(test_models + '/parentheses/test_parens.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     @unittest.skip('low priority')
     def test_reference_capitalization(self):
         """A properly formatted Vensim model should never create this failure"""
-        output, canon = runner(test_models + '/reference_capitalization/test_reference_capitalization.mdl')
+        output, canon = runner(test_models + '/reference_capitalization/test_reference_capitalization.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_repeated_subscript(self):
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
-            output, canon = runner(test_models + '/repeated_subscript/test_repeated_subscript.mdl')
+            output, canon = runner(test_models + '/repeated_subscript/test_repeated_subscript.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_rounding(self):
-        output, canon = runner(test_models + '/rounding/test_rounding.mdl')
+        output, canon = runner(test_models + '/rounding/test_rounding.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_sample_if_true(self):
-        output, canon = runner(test_models + '/sample_if_true/test_sample_if_true.mdl')
+        output, canon = runner(test_models + '/sample_if_true/test_sample_if_true.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_smooth(self):
-        output, canon = runner(test_models + '/smooth/test_smooth.mdl')
+        output, canon = runner(test_models + '/smooth/test_smooth.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_smooth_and_stock(self):
-        output, canon = runner(test_models + '/smooth_and_stock/test_smooth_and_stock.mdl')
+        output, canon = runner(test_models + '/smooth_and_stock/test_smooth_and_stock.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_special_characters(self):
-        output, canon = runner(test_models + '/special_characters/test_special_variable_names.mdl')
+        output, canon = runner(test_models + '/special_characters/test_special_variable_names.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_sqrt(self):
-        output, canon = runner(test_models + '/sqrt/test_sqrt.mdl')
+        output, canon = runner(test_models + '/sqrt/test_sqrt.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subrange_merge(self):
-        output, canon = runner(test_models + '/subrange_merge/test_subrange_merge.mdl')
+        output, canon = runner(test_models + '/subrange_merge/test_subrange_merge.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_logicals(self):
-        output, canon = runner(test_models + '/subscript_logicals/test_subscript_logicals.mdl')
+        output, canon = runner(test_models + '/subscript_logicals/test_subscript_logicals.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_multiples(self):
-        output, canon = runner(test_models + '/subscript_multiples/test_multiple_subscripts.mdl')
+        output, canon = runner(test_models + '/subscript_multiples/test_multiple_subscripts.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_1d_arrays(self):
-        output, canon = runner(test_models + '/subscript_1d_arrays/test_subscript_1d_arrays.mdl')
+        output, canon = runner(test_models + '/subscript_1d_arrays/test_subscript_1d_arrays.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_2d_arrays(self):
-        output, canon = runner(test_models + '/subscript_2d_arrays/test_subscript_2d_arrays.mdl')
+        output, canon = runner(test_models + '/subscript_2d_arrays/test_subscript_2d_arrays.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_3d_arrays(self):
-        output, canon = runner(test_models + '/subscript_3d_arrays/test_subscript_3d_arrays.mdl')
+        output, canon = runner(test_models + '/subscript_3d_arrays/test_subscript_3d_arrays.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_3d_arrays_lengthwise(self):
-        output, canon = runner(test_models + '/subscript_3d_arrays_lengthwise/test_subscript_3d_arrays_lengthwise.mdl')
+        output, canon = runner(test_models + '/subscript_3d_arrays_lengthwise/test_subscript_3d_arrays_lengthwise.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_3d_arrays_widthwise(self):
-        output, canon = runner(test_models + '/subscript_3d_arrays_widthwise/test_subscript_3d_arrays_widthwise.mdl')
+        output, canon = runner(test_models + '/subscript_3d_arrays_widthwise/test_subscript_3d_arrays_widthwise.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_aggregation(self):
-        output, canon = runner(test_models + '/subscript_aggregation/test_subscript_aggregation.mdl')
+        output, canon = runner(test_models + '/subscript_aggregation/test_subscript_aggregation.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_constant_call(self):
-        output, canon = runner(test_models + '/subscript_constant_call/test_subscript_constant_call.mdl')
+        output, canon = runner(test_models + '/subscript_constant_call/test_subscript_constant_call.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_copy(self):
-        output, canon = runner(test_models + '/subscript_copy/test_subscript_copy.mdl')
+        output, canon = runner(test_models + '/subscript_copy/test_subscript_copy.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_docs(self):
-        output, canon = runner(test_models + '/subscript_docs/subscript_docs.mdl')
+        output, canon = runner(test_models + '/subscript_docs/subscript_docs.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_element_name(self):
         # issue https://github.com/JamesPHoughton/pysd/issues/216
-        output, canon = runner(test_models + '/subscript_element_name/test_subscript_element_name.mdl')
+        output, canon = runner(test_models + '/subscript_element_name/test_subscript_element_name.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_individually_defined_1_of_2d_arrays(self):
-        output, canon = runner(test_models + '/subscript_individually_defined_1_of_2d_arrays/subscript_individually_defined_1_of_2d_arrays.mdl')
+        output, canon = runner(test_models + '/subscript_individually_defined_1_of_2d_arrays/subscript_individually_defined_1_of_2d_arrays.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_individually_defined_1_of_2d_arrays_from_floats(self):
-        output, canon = runner(test_models + '/subscript_individually_defined_1_of_2d_arrays_from_floats/subscript_individually_defined_1_of_2d_arrays_from_floats.mdl')
+        output, canon = runner(test_models + '/subscript_individually_defined_1_of_2d_arrays_from_floats/subscript_individually_defined_1_of_2d_arrays_from_floats.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_individually_defined_1d_arrays(self):
-        output, canon = runner(test_models + '/subscript_individually_defined_1d_arrays/subscript_individually_defined_1d_arrays.mdl')
+        output, canon = runner(test_models + '/subscript_individually_defined_1d_arrays/subscript_individually_defined_1d_arrays.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_individually_defined_stocks(self):
-        output, canon = runner(test_models + '/subscript_individually_defined_stocks/test_subscript_individually_defined_stocks.mdl')
+        output, canon = runner(test_models + '/subscript_individually_defined_stocks/test_subscript_individually_defined_stocks.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_mapping_simple(self):
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
-            output, canon = runner(test_models + '/subscript_mapping_simple/test_subscript_mapping_simple.mdl')
+            output, canon = runner(test_models + '/subscript_mapping_simple/test_subscript_mapping_simple.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_mapping_vensim(self):
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
-            output, canon = runner(test_models + '/subscript_mapping_vensim/test_subscript_mapping_vensim.mdl')
+            output, canon = runner(test_models + '/subscript_mapping_vensim/test_subscript_mapping_vensim.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_mixed_assembly(self):
-        output, canon = runner(test_models + '/subscript_mixed_assembly/test_subscript_mixed_assembly.mdl')
+        output, canon = runner(test_models + '/subscript_mixed_assembly/test_subscript_mixed_assembly.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_selection(self):
-        output, canon = runner(test_models + '/subscript_selection/subscript_selection.mdl')
+        output, canon = runner(test_models + '/subscript_selection/subscript_selection.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_numeric_range(self):
-        output, canon = runner(test_models + '/subscript_numeric_range/test_subscript_numeric_range.mdl')
+        output, canon = runner(test_models + '/subscript_numeric_range/test_subscript_numeric_range.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_subranges(self):
-        output, canon = runner(test_models + '/subscript_subranges/test_subscript_subrange.mdl')
+        output, canon = runner(test_models + '/subscript_subranges/test_subscript_subrange.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_subranges_equal(self):
-        output, canon = runner(test_models + '/subscript_subranges_equal/test_subscript_subrange_equal.mdl')
+        output, canon = runner(test_models + '/subscript_subranges_equal/test_subscript_subrange_equal.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_switching(self):
-        output, canon = runner(test_models + '/subscript_switching/subscript_switching.mdl')
+        output, canon = runner(test_models + '/subscript_switching/subscript_switching.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_transposition(self):
-        output, canon = runner(test_models + '/subscript_transposition/test_subscript_transposition.mdl')
+        output, canon = runner(test_models + '/subscript_transposition/test_subscript_transposition.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscript_updimensioning(self):
-        output, canon = runner(test_models + '/subscript_updimensioning/test_subscript_updimensioning.mdl')
+        output, canon = runner(test_models + '/subscript_updimensioning/test_subscript_updimensioning.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscripted_delays(self):
-        output, canon = runner(test_models + '/subscripted_delays/test_subscripted_delays.mdl')
+        output, canon = runner(test_models + '/subscripted_delays/test_subscripted_delays.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscripted_flows(self):
-        output, canon = runner(test_models + '/subscripted_flows/test_subscripted_flows.mdl')
+        output, canon = runner(test_models + '/subscripted_flows/test_subscripted_flows.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscripted_if_then_else(self):
-        output, canon = runner(test_models + '/subscripted_if_then_else/test_subscripted_if_then_else.mdl')
+        output, canon = runner(test_models + '/subscripted_if_then_else/test_subscripted_if_then_else.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscripted_logicals(self):
-        output, canon = runner(test_models + '/subscripted_logicals/test_subscripted_logicals.mdl')
+        output, canon = runner(test_models + '/subscripted_logicals/test_subscripted_logicals.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscripted_smooth(self):
         # issue https://github.com/JamesPHoughton/pysd/issues/226
-        output, canon = runner(test_models + '/subscripted_smooth/test_subscripted_smooth.mdl')
+        output, canon = runner(test_models + '/subscripted_smooth/test_subscripted_smooth.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscripted_trend(self):
         # issue https://github.com/JamesPHoughton/pysd/issues/226
-        output, canon = runner(test_models + '/subscripted_trend/test_subscripted_trend.mdl')
+        output, canon = runner(test_models + '/subscripted_trend/test_subscripted_trend.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subscripted_xidz(self):
-        output, canon = runner(test_models + '/subscripted_xidz/test_subscripted_xidz.mdl')
+        output, canon = runner(test_models + '/subscripted_xidz/test_subscripted_xidz.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_subset_duplicated_coord(self):
         output, canon = runner(test_models + '/subset_duplicated_coord/'
-                               + 'test_subset_duplicated_coord.mdl')
+                               + 'test_subset_duplicated_coord.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_time(self):
-        output, canon = runner(test_models + '/time/test_time.mdl')
+        output, canon = runner(test_models + '/time/test_time.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_trend(self):
-        output, canon = runner(test_models + '/trend/test_trend.mdl')
+        output, canon = runner(test_models + '/trend/test_trend.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_trig(self):
-        output, canon = runner(test_models + '/trig/test_trig.mdl')
+        output, canon = runner(test_models + '/trig/test_trig.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_variable_ranges(self):
-        output, canon = runner(test_models + '/variable_ranges/test_variable_ranges.mdl')
+        output, canon = runner(test_models + '/variable_ranges/test_variable_ranges.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_unicode_characters(self):
-        output, canon = runner(test_models + '/unicode_characters/unicode_test_model.mdl')
+        output, canon = runner(test_models + '/unicode_characters/unicode_test_model.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_xidz_zidz(self):
-        output, canon = runner(test_models + '/xidz_zidz/xidz_zidz.mdl')
+        output, canon = runner(test_models + '/xidz_zidz/xidz_zidz.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_run_uppercase(self):
-        output, canon = runner(test_models + '/case_sensitive_extension/teacup-upper.MDL')
+        output, canon = runner(test_models + '/case_sensitive_extension/teacup-upper.MDL', old=True)
         assert_frames_close(output, canon, rtol=rtol)
 
     def test_odd_number_quotes(self):
-        output, canon = runner(test_models + '/odd_number_quotes/teacup_3quotes.mdl')
+        output, canon = runner(test_models + '/odd_number_quotes/teacup_3quotes.mdl', old=True)
         assert_frames_close(output, canon, rtol=rtol)
diff --git a/tests/more-tests/subscript_individually_defined_stocks2/test_subscript_individually_defined_stocks2.mdl b/tests/more-tests/subscript_individually_defined_stocks2/test_subscript_individually_defined_stocks2.mdl
index aa30749e..e9470cf1 100644
--- a/tests/more-tests/subscript_individually_defined_stocks2/test_subscript_individually_defined_stocks2.mdl
+++ b/tests/more-tests/subscript_individually_defined_stocks2/test_subscript_individually_defined_stocks2.mdl
@@ -23,9 +23,9 @@ Third Dimension Subscript:
 	~		|
 	
 Initial Values[One Dimensional Subscript,Second Dimension Subscript,Depth 1]=
-	Initial Values A ~~|
+	Initial Values A[One Dimensional Subscript,Second Dimension Subscript] ~~|
 Initial Values[One Dimensional Subscript,Second Dimension Subscript,Depth 2]=
-	Initial Values B
+	Initial Values B[One Dimensional Subscript,Second Dimension Subscript]
 	~	
 	~		|
 	
diff --git a/tests/pytest_integration/vensim_pathway/pytest_integration_vensim_pathway.py b/tests/pytest_integration/vensim_pathway/pytest_integration_vensim_pathway.py
new file mode 100644
index 00000000..22876b26
--- /dev/null
+++ b/tests/pytest_integration/vensim_pathway/pytest_integration_vensim_pathway.py
@@ -0,0 +1,523 @@
+import pytest
+from pysd.tools.benchmarking import runner, assert_frames_close
+
+# TODO add warnings catcher per test
+
+vensim_test = {
+    "abs": {
+        "folder": "abs",
+        "file": "test_abs.mdl"
+    },
+    "active_initial": {
+        "folder": "active_initial",
+        "file": "test_active_initial.mdl"
+    },
+    "active_initial_circular": {
+        "folder": "active_initial_circular",
+        "file": "test_active_initial_circular.mdl"
+    },
+    "arithmetics": {
+        "folder": "arithmetics",
+        "file": "test_arithmetics.mdl"
+    },
+    "arguments": {
+        "folder": "arguments",
+        "file": "test_arguments.mdl",
+        "rtol": 1e-2 # TODO test why it is failing with smaller tolerance
+    },
+    "array_with_line_break": {
+        "folder": "array_with_line_break",
+        "file": "test_array_with_line_break.mdl"
+    },
+    "builtin_max": {
+        "folder": "builtin_max",
+        "file": "builtin_max.mdl"
+    },
+    "builtin_min": {
+        "folder": "builtin_min",
+        "file": "builtin_min.mdl"
+    },
+    "chained_initialization": {
+        "folder": "chained_initialization",
+        "file": "test_chained_initialization.mdl"
+    },
+    "conditional_subscripts": {
+        "folder": "conditional_subscripts",
+        "file": "test_conditional_subscripts.mdl"
+    },
+    "constant_expressions": {
+        "folder": "constant_expressions",
+        "file": "test_constant_expressions.mdl"
+    },
+    "control_vars": {
+        "folder": "control_vars",
+        "file": "test_control_vars.mdl"
+    },
+    "data_from_other_model": {
+        "folder": "data_from_other_model",
+        "file": "test_data_from_other_model.mdl",
+        "data_files": "data.tab"
+    },
+    "delay_fixed": {
+        "folder": "delay_fixed",
+        "file": "test_delay_fixed.mdl"
+    },
+    "delay_numeric_error": {
+        "folder": "delay_numeric_error",
+        "file": "test_delay_numeric_error.mdl"
+    },
+    "delay_parentheses": {
+        "folder": "delay_parentheses",
+        "file": "test_delay_parentheses.mdl"
+    },
+    "delay_pipeline": {
+        "folder": "delay_pipeline",
+        "file": "test_pipeline_delays.mdl"
+    },
+    "delays": {
+        "folder": "delays",
+        "file": "test_delays.mdl"
+    },
+    "dynamic_final_time": {
+        "folder": "dynamic_final_time",
+        "file": "test_dynamic_final_time.mdl"
+    },
+    "euler_step_vs_saveper": {
+        "folder": "euler_step_vs_saveper",
+        "file": "test_euler_step_vs_saveper.mdl"
+    },
+    "except": {
+        "folder": "except",
+        "file": "test_except.mdl"
+    },
+    "exp": {
+        "folder": "exp",
+        "file": "test_exp.mdl"
+    },
+    "exponentiation": {
+        "folder": "exponentiation",
+        "file": "exponentiation.mdl"
+    },
+    "forecast": {
+        "folder": "forecast",
+        "file": "test_forecast.mdl"
+    },
+    "function_capitalization": {
+        "folder": "function_capitalization",
+        "file": "test_function_capitalization.mdl"
+    },
+    "game": {
+        "folder": "game",
+        "file": "test_game.mdl"
+    },
+    "get_constants": pytest.param({
+        "folder": "get_constants",
+        "file": "test_get_constants.mdl"
+    },  marks=pytest.mark.xfail(reason="csv files not implemented")),
+    "get_constants_subranges": {
+        "folder": "get_constants_subranges",
+        "file": "test_get_constants_subranges.mdl"
+    },
+    "get_data": pytest.param({
+        "folder": "get_data",
+        "file": "test_get_data.mdl"
+    },  marks=pytest.mark.xfail(reason="csv files not implemented")),
+    "get_data_args_3d_xls": {
+        "folder": "get_data_args_3d_xls",
+        "file": "test_get_data_args_3d_xls.mdl"
+    },
+    "get_lookups_data_3d_xls": {
+        "folder": "get_lookups_data_3d_xls",
+        "file": "test_get_lookups_data_3d_xls.mdl"
+    },
+    "get_lookups_subscripted_args": {
+        "folder": "get_lookups_subscripted_args",
+        "file": "test_get_lookups_subscripted_args.mdl"
+    },
+    "get_lookups_subset": {
+        "folder": "get_lookups_subset",
+        "file": "test_get_lookups_subset.mdl"
+    },
+    "get_mixed_definitions": {
+        "folder": "get_mixed_definitions",
+        "file": "test_get_mixed_definitions.mdl"
+    },
+    "get_subscript_3d_arrays_xls": {
+        "folder": "get_subscript_3d_arrays_xls",
+        "file": "test_get_subscript_3d_arrays_xls.mdl"
+    },
+    "get_with_missing_values_xlsx": {
+        "folder": "get_with_missing_values_xlsx",
+        "file": "test_get_with_missing_values_xlsx.mdl"
+    },
+    "get_xls_cellrange": {
+        "folder": "get_xls_cellrange",
+        "file": "test_get_xls_cellrange.mdl"
+    },
+    "if_stmt": {
+        "folder": "if_stmt",
+        "file": "if_stmt.mdl"
+    },
+    "initial_function": {
+        "folder": "initial_function",
+        "file": "test_initial.mdl"
+    },
+    "input_functions": {
+        "folder": "input_functions",
+        "file": "test_inputs.mdl"
+    },
+    "invert_matrix": {
+        "folder": "invert_matrix",
+        "file": "test_invert_matrix.mdl"
+    },
+    "limits": {
+        "folder": "limits",
+        "file": "test_limits.mdl"
+    },
+    "line_breaks": {
+        "folder": "line_breaks",
+        "file": "test_line_breaks.mdl"
+    },
+    "line_continuation": {
+        "folder": "line_continuation",
+        "file": "test_line_continuation.mdl"
+    },
+    "ln": {
+        "folder": "ln",
+        "file": "test_ln.mdl"
+    },
+    "log": {
+        "folder": "log",
+        "file": "test_log.mdl"
+    },
+    "logicals": {
+        "folder": "logicals",
+        "file": "test_logicals.mdl"
+    },
+    "lookups": {
+        "folder": "lookups",
+        "file": "test_lookups.mdl"
+    },
+    "lookups_funcnames": {
+        "folder": "lookups_funcnames",
+        "file": "test_lookups_funcnames.mdl"
+    },
+    "lookups_inline": {
+        "folder": "lookups_inline",
+        "file": "test_lookups_inline.mdl"
+    },
+    "lookups_inline_bounded": {
+        "folder": "lookups_inline_bounded",
+        "file": "test_lookups_inline_bounded.mdl"
+    },
+    "lookups_with_expr": {
+        "folder": "lookups_with_expr",
+        "file": "test_lookups_with_expr.mdl"
+    },
+    "lookups_without_range": {
+        "folder": "lookups_without_range",
+        "file": "test_lookups_without_range.mdl"
+    },
+    "macro_cross_reference": {
+        "folder": "macro_cross_reference",
+        "file": "test_macro_cross_reference.mdl"
+    },
+    "macro_expression": {
+        "folder": "macro_expression",
+        "file": "test_macro_expression.mdl"
+    },
+    "macro_multi_expression": {
+        "folder": "macro_multi_expression",
+        "file": "test_macro_multi_expression.mdl"
+    },
+    "macro_multi_macros": {
+        "folder": "macro_multi_macros",
+        "file": "test_macro_multi_macros.mdl"
+    },
+    "macro_stock": {
+        "folder": "macro_stock",
+        "file": "test_macro_stock.mdl"
+    },
+    "macro_trailing_definition": {
+        "folder": "macro_trailing_definition",
+        "file": "test_macro_trailing_definition.mdl"
+    },
+    "model_doc": {
+        "folder": "model_doc",
+        "file": "model_doc.mdl"
+    },
+    "multiple_lines_def": {
+        "folder": "multiple_lines_def",
+        "file": "test_multiple_lines_def.mdl"
+    },
+    "nested_functions": {
+        "folder": "nested_functions",
+        "file": "test_nested_functions.mdl"
+    },
+    "number_handling": {
+        "folder": "number_handling",
+        "file": "test_number_handling.mdl"
+    },
+    "odd_number_quotes": {
+        "folder": "odd_number_quotes",
+        "file": "teacup_3quotes.mdl"
+    },
+    "parentheses": {
+        "folder": "parentheses",
+        "file": "test_parens.mdl"
+    },
+    "reference_capitalization": {
+        "folder": "reference_capitalization",
+        "file": "test_reference_capitalization.mdl"
+    },
+    "repeated_subscript": {
+        "folder": "repeated_subscript",
+        "file": "test_repeated_subscript.mdl"
+    },
+    "rounding": {
+        "folder": "rounding",
+        "file": "test_rounding.mdl"
+    },
+    "sample_if_true": {
+        "folder": "sample_if_true",
+        "file": "test_sample_if_true.mdl"
+    },
+    "smooth": {
+        "folder": "smooth",
+        "file": "test_smooth.mdl"
+    },
+    "smooth_and_stock": {
+        "folder": "smooth_and_stock",
+        "file": "test_smooth_and_stock.mdl"
+    },
+    "special_characters": {
+        "folder": "special_characters",
+        "file": "test_special_variable_names.mdl"
+    },
+    "sqrt": {
+        "folder": "sqrt",
+        "file": "test_sqrt.mdl"
+    },
+    "subrange_merge": {
+        "folder": "subrange_merge",
+        "file": "test_subrange_merge.mdl"
+    },
+    "subscript_1d_arrays": {
+        "folder": "subscript_1d_arrays",
+        "file": "test_subscript_1d_arrays.mdl"
+    },
+    "subscript_2d_arrays": {
+        "folder": "subscript_2d_arrays",
+        "file": "test_subscript_2d_arrays.mdl"
+    },
+    "subscript_3d_arrays": {
+        "folder": "subscript_3d_arrays",
+        "file": "test_subscript_3d_arrays.mdl"
+    },
+    "subscript_3d_arrays_lengthwise": {
+        "folder": "subscript_3d_arrays_lengthwise",
+        "file": "test_subscript_3d_arrays_lengthwise.mdl"
+    },
+    "subscript_3d_arrays_widthwise": {
+        "folder": "subscript_3d_arrays_widthwise",
+        "file": "test_subscript_3d_arrays_widthwise.mdl"
+    },
+    "subscript_aggregation": {
+        "folder": "subscript_aggregation",
+        "file": "test_subscript_aggregation.mdl"
+    },
+    "subscript_constant_call": {
+        "folder": "subscript_constant_call",
+        "file": "test_subscript_constant_call.mdl"
+    },
+    "subscript_copy": {
+        "folder": "subscript_copy",
+        "file": "test_subscript_copy.mdl"
+    },
+    "subscript_docs": {
+        "folder": "subscript_docs",
+        "file": "subscript_docs.mdl"
+    },
+    "subscript_element_name": {
+        "folder": "subscript_element_name",
+        "file": "test_subscript_element_name.mdl"
+    },
+    "subscript_individually_defined_1_of_2d_arrays": {
+        "folder": "subscript_individually_defined_1_of_2d_arrays",
+        "file": "subscript_individually_defined_1_of_2d_arrays.mdl"
+    },
+    "subscript_individually_defined_1_of_2d_arrays_from_floats": {
+        "folder": "subscript_individually_defined_1_of_2d_arrays_from_floats",
+        "file": "subscript_individually_defined_1_of_2d_arrays_from_floats.mdl"
+    },
+    "subscript_individually_defined_1d_arrays": {
+        "folder": "subscript_individually_defined_1d_arrays",
+        "file": "subscript_individually_defined_1d_arrays.mdl"
+    },
+    "subscript_individually_defined_stocks": {
+        "folder": "subscript_individually_defined_stocks",
+        "file": "test_subscript_individually_defined_stocks.mdl"
+    },
+    "subscript_logicals": {
+        "folder": "subscript_logicals",
+        "file": "test_subscript_logicals.mdl"
+    },
+    "subscript_mapping_simple": {
+        "folder": "subscript_mapping_simple",
+        "file": "test_subscript_mapping_simple.mdl"
+    },
+    "subscript_mapping_vensim": {
+        "folder": "subscript_mapping_vensim",
+        "file": "test_subscript_mapping_vensim.mdl"
+    },
+    "subscript_mixed_assembly": {
+        "folder": "subscript_mixed_assembly",
+        "file": "test_subscript_mixed_assembly.mdl"
+    },
+    "subscript_multiples": {
+        "folder": "subscript_multiples",
+        "file": "test_multiple_subscripts.mdl"
+    },
+    "subscript_numeric_range": {
+        "folder": "subscript_numeric_range",
+        "file": "test_subscript_numeric_range.mdl"
+    },
+    "subscript_selection": {
+        "folder": "subscript_selection",
+        "file": "subscript_selection.mdl"
+    },
+    "subscript_subranges": {
+        "folder": "subscript_subranges",
+        "file": "test_subscript_subrange.mdl"
+    },
+    "subscript_subranges_equal": {
+        "folder": "subscript_subranges_equal",
+        "file": "test_subscript_subrange_equal.mdl"
+    },
+    "subscript_switching": {
+        "folder": "subscript_switching",
+        "file": "subscript_switching.mdl"
+    },
+    "subscript_transposition": {
+        "folder": "subscript_transposition",
+        "file": "test_subscript_transposition.mdl"
+    },
+    "subscript_updimensioning": {
+        "folder": "subscript_updimensioning",
+        "file": "test_subscript_updimensioning.mdl"
+    },
+    "subscripted_delays": {
+        "folder": "subscripted_delays",
+        "file": "test_subscripted_delays.mdl"
+    },
+    "subscripted_flows": {
+        "folder": "subscripted_flows",
+        "file": "test_subscripted_flows.mdl"
+    },
+    "subscripted_if_then_else": {
+        "folder": "subscripted_if_then_else",
+        "file": "test_subscripted_if_then_else.mdl"
+    },
+    "subscripted_logicals": {
+        "folder": "subscripted_logicals",
+        "file": "test_subscripted_logicals.mdl"
+    },
+    "subscripted_lookups": {
+        "folder": "subscripted_lookups",
+        "file": "test_subscripted_lookups.mdl"
+    },
+    "subscripted_round": {
+        "folder": "subscripted_round",
+        "file": "test_subscripted_round.mdl"
+    },
+    "subscripted_smooth": {
+        "folder": "subscripted_smooth",
+        "file": "test_subscripted_smooth.mdl"
+    },
+    "subscripted_trend": {
+        "folder": "subscripted_trend",
+        "file": "test_subscripted_trend.mdl"
+    },
+    "subscripted_xidz": {
+        "folder": "subscripted_xidz",
+        "file": "test_subscripted_xidz.mdl"
+    },
+    "subset_duplicated_coord": {
+        "folder": "subset_duplicated_coord",
+        "file": "test_subset_duplicated_coord.mdl"
+    },
+    "time": {
+        "folder": "time",
+        "file": "test_time.mdl"
+    },
+    "trend": {
+        "folder": "trend",
+        "file": "test_trend.mdl"
+    },
+    "trig": {
+        "folder": "trig",
+        "file": "test_trig.mdl"
+    },
+    "unicode_characters": {
+        "folder": "unicode_characters",
+        "file": "unicode_test_model.mdl"
+    },
+    "variable_ranges": {
+        "folder": "variable_ranges",
+        "file": "test_variable_ranges.mdl"
+    },
+    "xidz_zidz": {
+        "folder": "xidz_zidz",
+        "file": "xidz_zidz.mdl"
+    }
+}
+
+
+@pytest.mark.parametrize(
+    "test_data",
+    [item for item in vensim_test.values()],
+    ids=list(vensim_test)
+)
+class TestIntegrateVensim:
+    """
+    Test for splitting Vensim views in modules and submodules
+    """
+
+    @pytest.fixture
+    def model_path(self, _test_models, test_data):
+        return _test_models.joinpath(
+            test_data["folder"]).joinpath(test_data["file"])
+
+    @pytest.fixture
+    def data_path(self, _test_models, test_data):
+        """Fixture for models with data_path"""
+        if "data_files" in test_data:
+            if isinstance(test_data["data_files"], str):
+                return _test_models.joinpath(
+                    test_data["folder"]).joinpath(test_data["data_files"])
+            elif isinstance(test_data["data_files"], list):
+                return [
+                    _test_models.joinpath(test_data["folder"]).joinpath(file)
+                    for file in test_data["data_files"]
+                ]
+            else:
+                return {
+                    _test_models.joinpath(test_data["folder"]).joinpath(file):
+                    values for file, values in test_data["data_files"].items()
+                }
+        else:
+            return None
+
+    @pytest.fixture
+    def kwargs(self, test_data):
+        """Fixture for atol and rtol"""
+        kwargs = {}
+        if "atol" in test_data:
+            kwargs["atol"] = test_data["atol"]
+        if "rtol" in test_data:
+            kwargs["rtol"] = test_data["rtol"]
+        return kwargs
+
+    def test_read_vensim_file(self, model_path, data_path, kwargs):
+        output, canon = runner(model_path, data_files=data_path)
+        assert_frames_close(output, canon, **kwargs)
diff --git a/tests/pytest_pysd/user_interaction/pytest_select_submodel.py b/tests/pytest_pysd/user_interaction/pytest_select_submodel.py
index 0e101fe0..e868e06d 100644
--- a/tests/pytest_pysd/user_interaction/pytest_select_submodel.py
+++ b/tests/pytest_pysd/user_interaction/pytest_select_submodel.py
@@ -146,18 +146,16 @@ def test_select_submodel(self, model, variables, modules,
         if not dep_vars:
             # totally independent submodels can run without producing
             # nan values
-            assert len(record) == 1
             assert not np.any(np.isnan(model.run()))
         else:
             # running the model without redefining dependencies will
             # produce nan values
-            assert len(record) == 2
             assert "Exogenous components for the following variables are"\
-                + " necessary but not given:" in str(record[1].message)
+                + " necessary but not given:" in str(record[-1].message)
             assert "Please, set them before running the model using "\
-                + "set_components method..." in str(record[1].message)
+                + "set_components method..." in str(record[-1].message)
             for var in dep_vars:
-                assert var in str(record[1].message)
+                assert var in str(record[-1].message)
             assert np.any(np.isnan(model.run()))
             # redefine dependencies
             assert not np.any(np.isnan(model.run(params=dep_vars)))
@@ -168,7 +166,6 @@ def test_select_submodel(self, model, variables, modules,
             model.select_submodel(vars=variables, modules=modules,
                                   exogenous_components=dep_vars)
 
-        assert len(record) == 1
         assert not np.any(np.isnan(model.run()))
 
 
diff --git a/tests/pytest_translation/vensim_parser/pytest_vensim_file.py b/tests/pytest_translation/vensim_parser/pytest_vensim_file.py
new file mode 100644
index 00000000..135dd7d8
--- /dev/null
+++ b/tests/pytest_translation/vensim_parser/pytest_vensim_file.py
@@ -0,0 +1,54 @@
+
+import pytest
+from pathlib import Path
+
+from pysd.translation.vensim.vensin_file import VensimFile
+
+
+@pytest.mark.parametrize(
+    "path",
+    [
+        (  # teacup
+            "test-models/samples/teacup/teacup.mdl"
+        ),
+        (  # macros
+            "test-models/tests/macro_multi_expression/test_macro_multi_expression.mdl"
+        ),
+        (  # mapping
+            "test-models/tests/subscript_mapping_vensim/test_subscript_mapping_vensim.mdl"
+        ),
+        (  # data
+            "test-models/tests/data_from_other_model/test_data_from_other_model.mdl"
+        ),
+        (  # except
+            "test-models/tests/except/test_except.mdl"
+        )
+    ],
+    ids=["teacup", "macros", "mapping", "data", "except"]
+)
+class TestVensimFile:
+    """
+    Test for splitting Vensim views in modules and submodules
+    """
+    @pytest.fixture
+    def model_path(self, _root, path):
+        return _root.joinpath(path)
+
+    @pytest.mark.dependency(name="read_vensim_file")
+    def test_read_vensim_file(self, request, path, model_path):
+        # assert that the files don't exist in the temporary directory
+        ven_file = VensimFile(model_path)
+
+        assert hasattr(ven_file, "mdl_path")
+        assert hasattr(ven_file, "root_path")
+        assert hasattr(ven_file, "model_text")
+
+        assert isinstance(getattr(ven_file, "mdl_path"), Path)
+        assert isinstance(getattr(ven_file, "root_path"), Path)
+        assert isinstance(getattr(ven_file, "model_text"), str)
+
+    @pytest.mark.dependency(depends=["read_vensim_file"])
+    def test_file_split_file_sections(self, request, path, model_path):
+        ven_file = VensimFile(model_path)
+        ven_file.parse()
+        print(ven_file.verbose)
\ No newline at end of file
diff --git a/tests/test-models b/tests/test-models
index de294a8a..75ea19ba 160000
--- a/tests/test-models
+++ b/tests/test-models
@@ -1 +1 @@
-Subproject commit de294a8ad0f2c1a2bf41c351cfc4ab637bc39825
+Subproject commit 75ea19badf2ed6e94aba0f707ef5d6c97d80195b
diff --git a/tests/unit_test_benchmarking.py b/tests/unit_test_benchmarking.py
index 696c6c46..2bc74fc2 100644
--- a/tests/unit_test_benchmarking.py
+++ b/tests/unit_test_benchmarking.py
@@ -28,7 +28,7 @@ def test_non_valid_model(self):
                 "more-tests/not_vensim/test_not_vensim.txt"))
 
         self.assertIn(
-            'Modelfile should be *.mdl or *.xmile',
+            'Modelfile should be *.mdl, *.xmile, or *.py',
             str(err.exception))
 
     def test_different_frames_error(self):
diff --git a/tests/unit_test_external.py b/tests/unit_test_external.py
index 78f0aae3..0b2cadbe 100644
--- a/tests/unit_test_external.py
+++ b/tests/unit_test_external.py
@@ -207,12 +207,12 @@ def test_fill_missing(self):
         interp = np.array([1., 1., 1., 3., 3.5, 4.,
                            5., 6., 7., 8., 8., 8.])
 
-        ext.interp = "hold backward"
+        ext.interp = "hold_backward"
         datac = data.copy()
         ext._fill_missing(series, datac)
         self.assertTrue(np.all(hold_back == datac))
 
-        ext.interp = "look forward"
+        ext.interp = "look_forward"
         datac = data.copy()
         ext._fill_missing(series, datac)
         self.assertTrue(np.all(look_for == datac))
@@ -456,7 +456,7 @@ def test_data_interp_vn1d(self):
 
     def test_data_forward_h1d(self):
         """
-        ExtData test for 1d horizontal series look forward
+        ExtData test for 1d horizontal series look_forward
         """
         import pysd
 
@@ -465,7 +465,7 @@ def test_data_forward_h1d(self):
         time_row_or_col = "4"
         cell = "C5"
         coords = {}
-        interp = "look forward"
+        interp = "look_forward"
         py_name = "test_data_forward_h1d"
 
         data = pysd.external.ExtData(file_name=file_name,
@@ -486,7 +486,7 @@ def test_data_forward_h1d(self):
 
     def test_data_forward_v1d(self):
         """
-        ExtData test for 1d vertical series look forward
+        ExtData test for 1d vertical series look_forward
         """
         import pysd
 
@@ -495,7 +495,7 @@ def test_data_forward_v1d(self):
         time_row_or_col = "B"
         cell = "C5"
         coords = {}
-        interp = "look forward"
+        interp = "look_forward"
         py_name = "test_data_forward_v1d"
 
         data = pysd.external.ExtData(file_name=file_name,
@@ -516,7 +516,7 @@ def test_data_forward_v1d(self):
 
     def test_data_forward_hn1d(self):
         """
-        ExtData test for 1d horizontal series look forward by cell range names
+        ExtData test for 1d horizontal series look_forward by cell range names
         """
         import pysd
 
@@ -525,7 +525,7 @@ def test_data_forward_hn1d(self):
         time_row_or_col = "time"
         cell = "data_1d"
         coords = {}
-        interp = "look forward"
+        interp = "look_forward"
         py_name = "test_data_forward_hn1d"
 
         data = pysd.external.ExtData(file_name=file_name,
@@ -546,7 +546,7 @@ def test_data_forward_hn1d(self):
 
     def test_data_forward_vn1d(self):
         """
-        ExtData test for 1d vertical series look forward by cell range names
+        ExtData test for 1d vertical series look_forward by cell range names
         """
         import pysd
 
@@ -555,7 +555,7 @@ def test_data_forward_vn1d(self):
         time_row_or_col = "time"
         cell = "data_1d"
         coords = {}
-        interp = "look forward"
+        interp = "look_forward"
         py_name = "test_data_forward_vn1d"
 
         data = pysd.external.ExtData(file_name=file_name,
@@ -576,7 +576,7 @@ def test_data_forward_vn1d(self):
 
     def test_data_backward_h1d(self):
         """
-        ExtData test for 1d horizontal series hold backward
+        ExtData test for 1d horizontal series hold_backward
         """
         import pysd
 
@@ -585,7 +585,7 @@ def test_data_backward_h1d(self):
         time_row_or_col = "4"
         cell = "C5"
         coords = {}
-        interp = "hold backward"
+        interp = "hold_backward"
         py_name = "test_data_backward_h1d"
 
         data = pysd.external.ExtData(file_name=file_name,
@@ -606,7 +606,7 @@ def test_data_backward_h1d(self):
 
     def test_data_backward_v1d(self):
         """
-        ExtData test for 1d vertical series hold backward by cell range names
+        ExtData test for 1d vertical series hold_backward by cell range names
         """
         import pysd
 
@@ -615,7 +615,7 @@ def test_data_backward_v1d(self):
         time_row_or_col = "B"
         cell = "C5"
         coords = {}
-        interp = "hold backward"
+        interp = "hold_backward"
         py_name = "test_data_backward_v1d"
 
         data = pysd.external.ExtData(file_name=file_name,
@@ -636,7 +636,7 @@ def test_data_backward_v1d(self):
 
     def test_data_backward_hn1d(self):
         """
-        ExtData test for 1d horizontal series hold backward by cell range names
+        ExtData test for 1d horizontal series hold_backward by cell range names
         """
         import pysd
 
@@ -645,7 +645,7 @@ def test_data_backward_hn1d(self):
         time_row_or_col = "time"
         cell = "data_1d"
         coords = {}
-        interp = "hold backward"
+        interp = "hold_backward"
         py_name = "test_data_backward_hn1d"
 
         data = pysd.external.ExtData(file_name=file_name,
@@ -666,7 +666,7 @@ def test_data_backward_hn1d(self):
 
     def test_data_backward_vn1d(self):
         """
-        ExtData test for 1d vertical series hold backward by cell range names
+        ExtData test for 1d vertical series hold_backward by cell range names
         """
         import pysd
 
@@ -675,7 +675,7 @@ def test_data_backward_vn1d(self):
         time_row_or_col = "time"
         cell = "data_1d"
         coords = {}
-        interp = "hold backward"
+        interp = "hold_backward"
         py_name = "test_data_backward_vn1d"
 
         data = pysd.external.ExtData(file_name=file_name,
@@ -726,7 +726,7 @@ def test_data_interp_vn2d(self):
 
     def test_data_forward_hn2d(self):
         """
-        ExtData test for 2d vertical series look forward by cell range names
+        ExtData test for 2d vertical series look_forward by cell range names
         """
         import pysd
 
@@ -735,7 +735,7 @@ def test_data_forward_hn2d(self):
         time_row_or_col = "time"
         cell = "data_2d"
         coords = {'ABC': ['A', 'B', 'C']}
-        interp = "look forward"
+        interp = "look_forward"
         py_name = "test_data_forward_hn2d"
 
         data = pysd.external.ExtData(file_name=file_name,
@@ -756,7 +756,7 @@ def test_data_forward_hn2d(self):
 
     def test_data_backward_v2d(self):
         """
-        ExtData test for 2d vertical series hold backward
+        ExtData test for 2d vertical series hold_backward
         """
         import pysd
 
@@ -765,7 +765,7 @@ def test_data_backward_v2d(self):
         time_row_or_col = "B"
         cell = "C5"
         coords = {'ABC': ['A', 'B', 'C']}
-        interp = "hold backward"
+        interp = "hold_backward"
         py_name = "test_data_backward_v2d"
 
         data = pysd.external.ExtData(file_name=file_name,
@@ -827,7 +827,7 @@ def test_data_interp_h3d(self):
 
     def test_data_forward_v3d(self):
         """
-        ExtData test for 3d vertical series look forward
+        ExtData test for 3d vertical series look_forward
         """
         import pysd
 
@@ -838,7 +838,7 @@ def test_data_forward_v3d(self):
         cell_2 = "F5"
         coords_1 = {'XY': ['X'], 'ABC': ['A', 'B', 'C']}
         coords_2 = {'XY': ['Y'], 'ABC': ['A', 'B', 'C']}
-        interp = "look forward"
+        interp = "look_forward"
         py_name = "test_data_forward_v3d"
 
         data = pysd.external.ExtData(file_name=file_name,
@@ -867,7 +867,7 @@ def test_data_forward_v3d(self):
 
     def test_data_backward_hn3d(self):
         """
-        ExtData test for 3d horizontal series hold backward by cellrange names
+        ExtData test for 3d horizontal series hold_backward by cellrange names
         """
         import pysd
 
@@ -878,7 +878,7 @@ def test_data_backward_hn3d(self):
         cell_2 = "data_2db"
         coords_1 = {'XY': ['X'], 'ABC': ['A', 'B', 'C']}
         coords_2 = {'XY': ['Y'], 'ABC': ['A', 'B', 'C']}
-        interp = "hold backward"
+        interp = "hold_backward"
         py_name = "test_data_backward_hn3d"
 
         data = pysd.external.ExtData(file_name=file_name,
@@ -2890,7 +2890,7 @@ def test_data_h3d_interp(self):
         coords_1 = {'ABC': ['A', 'B', 'C'], 'XY': ['X']}
         coords_2 = {'ABC': ['A', 'B', 'C'], 'XY': ['Y']}
         interp = None
-        interp2 = "look forward"
+        interp2 = "look_forward"
         py_name = "test_data_h3d_interp"
 
         data = pysd.external.ExtData(file_name=file_name,
diff --git a/tests/unit_test_pysd.py b/tests/unit_test_pysd.py
index f3cafa88..1502a4a3 100644
--- a/tests/unit_test_pysd.py
+++ b/tests/unit_test_pysd.py
@@ -581,6 +581,55 @@ def test_docs(self):
         model = pysd.read_vensim(test_model)
         self.assertIsInstance(str(model), str)  # tests string conversion of
         # model
+        print(model.doc().columns)
+
+        doc = model._doc
+        self.assertIsInstance(doc, pd.DataFrame)
+        self.assertSetEqual(
+            {
+                "Characteristic Time",
+                "Teacup Temperature",
+                "FINAL TIME",
+                "Heat Loss to Room",
+                "INITIAL TIME",
+                "Room Temperature",
+                "SAVEPER",
+                "TIME STEP",
+            },
+            set(doc["Real Name"].values),
+        )
+
+        self.assertEqual(
+            doc[doc["Real Name"] == "Heat Loss to Room"]["Unit"].values[0],
+            "Degrees Fahrenheit/Minute",
+        )
+        self.assertEqual(
+            doc[doc["Real Name"] == "Teacup Temperature"]["Py Name"].values[0],
+            "teacup_temperature",
+        )
+        self.assertEqual(
+            doc[doc["Real Name"] == "INITIAL TIME"]["Comment"].values[0],
+            "The initial time for the simulation.",
+        )
+        self.assertEqual(
+            doc[doc["Real Name"] == "Characteristic Time"]["Type"].values[0],
+            "Constant"
+        )
+        self.assertEqual(
+            doc[doc["Real Name"] == "Characteristic Time"]["Subtype"].values[0],
+            "Normal"
+        )
+        self.assertEqual(
+            doc[doc["Real Name"] == "Teacup Temperature"]["Lims"].values[0],
+            "(32.0, 212.0)",
+        )
+
+    def test_docs_old(self):
+        """ Test that the model prints some documentation """
+
+        model = pysd.read_vensim(test_model, old=True)
+        self.assertIsInstance(str(model), str)  # tests string conversion of
+        # model
 
         doc = model.doc()
         self.assertIsInstance(doc, pd.DataFrame)
@@ -636,8 +685,9 @@ def test_docs_multiline_eqn(self):
         self.assertEqual(
             doc[doc["Real Name"] == "price"]["Subs"].values[0], "['fruits']"
         )
-        self.assertEqual(doc[doc["Real Name"] == "price"]["Eqn"].values[0],
-                         "1.2; .; .; .; 1.4")
+        # TODO: keep eqn?
+        #self.assertEqual(doc[doc["Real Name"] == "price"]["Eqn"].values[0],
+        #                 "1.2; .; .; .; 1.4")
 
     def test_stepwise_cache(self):
         from pysd.py_backend.decorators import Cache
@@ -1382,7 +1432,7 @@ def test_multiple_deps(self):
                 + "test_subscript_individually_defined_stocks2.mdl"))
 
         expected_dep = {
-            "stock_a": {"_integ_stock_a": 2},
+            "stock_a": {"_integ_stock_a": 1, "_integ_stock_a_1": 1},
             "inflow_a": {"rate_a": 1},
             "inflow_b": {"rate_a": 1},
             "initial_values": {"initial_values_a": 1, "initial_values_b": 1},
@@ -1394,9 +1444,13 @@ def test_multiple_deps(self):
             "saveper": {"time_step": 1},
             "time_step": {},
             "_integ_stock_a": {
-                "initial": {"initial_values": 2},
-                "step": {"inflow_a": 1, "inflow_b": 1}
+                "initial": {"initial_values": 1},
+                "step": {"inflow_a": 1}
             },
+            '_integ_stock_a_1': {
+                'initial': {'initial_values': 1},
+                'step': {'inflow_b': 1}
+            }
         }
         self.assertEqual(model.components._dependencies, expected_dep)