diff --git a/.github/workflows/pr-messages.yml b/.github/workflows/pr-messages.yml
index 2e4d38077..47aea7bd4 100644
--- a/.github/workflows/pr-messages.yml
+++ b/.github/workflows/pr-messages.yml
@@ -18,10 +18,12 @@ jobs:
           pullRequestOpened: >
             Thank you for opening this PR. Each PR into dev requires a code review. For the code review, look at the following:
 
-            - [ ] Reviewer (someone other than author) should look for bugs, efficiency, readability, testing, and coverage in examples (if relevant).
+            - [ ] Reviewer should look for bugs, efficiency, readability, testing, and coverage in examples (if relevant).
 
             - [ ] Ensure that each PR adding a new feature should include a test verifying that feature.
 
+            - [ ] All tests must be passing.
+
             - [ ] All errors from static analysis must be resolved.
             
             - [ ] Review the test coverage reports (if there is a change) - will be added as comment on PR if there is a change
diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index cbf3b13ba..2ed76d6be 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -110,42 +110,31 @@ jobs:
     - name: Run copyright check
       run: |
         python scripts/test_copyright.py
-  coverage:
-    timeout-minutes: 30
-    runs-on: ubuntu-latest
-    strategy:
-      matrix:
-        python-version: ['3.10']
-    steps:
-    - uses: actions/checkout@v3
-    - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
-      with:
-        python-version: ${{ matrix.python-version }}
-    - name: Install dependencies
-      run: |
-        python -m pip install --upgrade pip
-        python -m pip install -e .
-        pip install coverage
-        pip install notebook
-        pip install testbook
-        pip install requests
-        pip install importlib_metadata
-    - name: Run coverage
-      run: |
-        coverage run -m tests.test_base_models
-        coverage run -a -m tests.test_battery
-        coverage run -a -m tests.test_centrifugal_pump
-        coverage run -a -m tests.test_composite
-        coverage run -a -m tests.test_datasets
-        coverage run -a -m tests.test_dict_like_matrix_wrapper
-        coverage run -a -m tests.test_direct
-        coverage run -a -m tests.test_ensemble
-        coverage run -a -m tests.test_linear_model
-        coverage run -a -m tests.test_pneumatic_valve
-        coverage run -a -m tests.test_powertrain
-        coverage run -a -m tests.test_serialization
-        coverage run -a -m tests.test_sim_result
-        coverage xml -i
-    - name: "Upload coverage to Codecov"
-      uses: codecov/codecov-action@v3
+  # coverage:
+  #   timeout-minutes: 30
+  #   runs-on: ubuntu-latest
+  #   strategy:
+  #     matrix:
+  #       python-version: ['3.9']
+  #   steps:
+  #   - uses: actions/checkout@v3
+  #   - name: Set up Python ${{ matrix.python-version }}
+  #     uses: actions/setup-python@v4
+  #     with:
+  #       python-version: ${{ matrix.python-version }}
+  #   - name: Install dependencies
+  #     run: |
+  #       python -m pip install --upgrade pip
+  #       python -m pip install -e .
+  #       pip install coverage
+  #       pip install notebook
+  #       pip install testbook
+  #       pip install requests
+  #   - name: Run coverage
+  #     run: |
+  #       coverage run -m tests
+  #       coverage xml
+  #   - name: "Upload coverage to Codecov"
+  #     uses: codecov/codecov-action@v3
+  #     with:
+  #       fail_ci_if_error: true
diff --git a/codecov.yml b/codecov.yml
deleted file mode 100644
index 990f060b4..000000000
--- a/codecov.yml
+++ /dev/null
@@ -1,5 +0,0 @@
-ignore:
-  - "tests"  # Don't include tests
-  - "src/prog_models/models/test_models"  # Dont include test models
-  - "src/prog_models/datasets" # Dont include datasets (too long a test)
-  - "src/prog_models/data_models" # Skip data models (doesn't work with coverage)
diff --git a/examples/param_est.py b/examples/param_est.py
index 8248a8df1..5de47ea7b 100644
--- a/examples/param_est.py
+++ b/examples/param_est.py
@@ -48,11 +48,5 @@ def run_example():
     
     # Sure enough- parameter estimation determined that the thrower's height wasn't 20 m, instead was closer to 1.9m, a much more reasonable height!
 
-    # Note: You can also adjust the metric that is used to estimate parameters.
-    # This is done by setting the "error_method" argument. 
-    # e.g., m.estimate_params([(times, inputs, outputs)], keys, dt=0.01, error_method='MAX_E')
-    # Default is Mean Squared Error (MSE)
-    # See calc_error method for list of options.
-
 if __name__=='__main__':
     run_example()
diff --git a/src/prog_models/__init__.py b/src/prog_models/__init__.py
index da1aff45f..bf9119046 100644
--- a/src/prog_models/__init__.py
+++ b/src/prog_models/__init__.py
@@ -2,11 +2,10 @@
 # National Aeronautics and Space Administration.  All Rights Reserved.
 
 # PrognosticsModel must be first, since the others build on this
-from prog_models.prognostics_model import PrognosticsModel
-from prog_models.ensemble_model import EnsembleModel
-from prog_models.composite_model import CompositeModel
-from prog_models.linear_model import LinearModel
-from prog_models.models.thrown_object import LinearThrownObject
-from prog_models.exceptions import ProgModelException, ProgModelInputException, ProgModelTypeError
+from .prognostics_model import PrognosticsModel
+from .ensemble_model import EnsembleModel
+from .composite_model import CompositeModel
+from .linear_model import LinearModel
+from .exceptions import ProgModelException, ProgModelInputException, ProgModelTypeError
 
 __version__ = '1.5.0.pre'
diff --git a/src/prog_models/composite_model.py b/src/prog_models/composite_model.py
index a772a5eca..c204e27ac 100644
--- a/src/prog_models/composite_model.py
+++ b/src/prog_models/composite_model.py
@@ -3,7 +3,7 @@
 
 from collections.abc import Iterable
 
-from prog_models import PrognosticsModel
+from . import PrognosticsModel
 
 DIVIDER = '.'
 
@@ -32,7 +32,8 @@ class CompositeModel(PrognosticsModel):
         outputs (list[str]):
             Model outputs in format "model_name.output_name". Must be subset of all outputs from models. If not provided, all outputs will be included. 
     """
-    def __init__(self, models, connections = [], **kwargs):
+
+    def __init__(self, models, connections=[], **kwargs):
         # General Input Validation
         if not isinstance(models, Iterable):
             raise ValueError('The models argument must be a list')
@@ -54,8 +55,9 @@ def __init__(self, models, connections = [], **kwargs):
         # Handle models
         for m in models:
             if isinstance(m, Iterable):
-                if len(m) != 2: 
-                    raise ValueError('Each model tuple must be of the form (name: str, model). For example ("Batt1", BatteryElectroChem())')
+                if len(m) != 2:
+                    raise ValueError(
+                        'Each model tuple must be of the form (name: str, model). For example ("Batt1", BatteryElectroChem())')
                 if not isinstance(m[0], str):
                     raise ValueError('The first element of each model tuple must be a string')
                 if not isinstance(m[1], PrognosticsModel):
@@ -73,7 +75,8 @@ def __init__(self, models, connections = [], **kwargs):
                 self.model_names.add(m[0])
                 kwargs['models'].append(m)
             else:
-                raise ValueError(f'Each model must be a PrognosticsModel or tuple (name: str, PrognosticsModel), was {type(m)}')
+                raise ValueError(
+                    f'Each model must be a PrognosticsModel or tuple (name: str, PrognosticsModel), was {type(m)}')
 
         for (name, m) in kwargs['models']:
             self.inputs |= set([name + DIVIDER + u for u in m.inputs])
@@ -81,7 +84,7 @@ def __init__(self, models, connections = [], **kwargs):
             self.outputs |= set([name + DIVIDER + z for z in m.outputs])
             self.events |= set([name + DIVIDER + e for e in m.events])
             self.performance_metric_keys |= set([name + DIVIDER + p for p in m.performance_metric_keys])
-        
+
         # Handle outputs
         if 'outputs' in kwargs:
             if isinstance(kwargs['outputs'], str):
@@ -91,11 +94,11 @@ def __init__(self, models, connections = [], **kwargs):
             if not set(kwargs['outputs']).issubset(self.outputs):
                 raise ValueError('The outputs of the composite model must be a subset of the outputs of the models')
             self.outputs = kwargs['outputs']
-        
+
         # Handle Connections
         kwargs['connections'] = []
-        self.__to_input_connections = {m_name : [] for m_name in self.model_names}
-        self.__to_state_connections = {m_name : [] for m_name in self.model_names}
+        self.__to_input_connections = {m_name: [] for m_name in self.model_names}
+        self.__to_state_connections = {m_name: [] for m_name in self.model_names}
 
         for connection in connections:
             # Input validation
@@ -107,13 +110,14 @@ def __init__(self, models, connections = [], **kwargs):
             in_key, out_key = connection
             # Validation
             if out_key not in self.inputs:
-                raise ValueError(f'The output key, {out_key}, must be an input to one of the composite models. Options include {self.inputs}')
+                raise ValueError(
+                    f'The output key, {out_key}, must be an input to one of the composite models. Options include {self.inputs}')
 
             # Remove the out_key from inputs
             # These no longer are an input to the composite model
             # as they are now satisfied internally
             self.inputs.remove(out_key)
-                
+
             # Split the keys into parts (model, key_part)
             (in_model, in_key_part) = in_key.split('.')
             (out_model, out_key_part) = out_key.split('.')
@@ -125,7 +129,7 @@ def __init__(self, models, connections = [], **kwargs):
                 raise ValueError('The input model must be one of the models in the composite model')
             if out_model not in self.model_names:
                 raise ValueError('The output model must be one of the models in the composite model')
-            
+
             # Add to connections
             if in_key in self.states:
                 self.__to_input_connections[out_model].append((in_key, out_key_part))
@@ -138,11 +142,11 @@ def __init__(self, models, connections = [], **kwargs):
                 self.states.add(in_key)
             else:
                 raise ValueError('The input key must be an output or state of one of the composite models')
-        
+
         # Finish initialization
         super().__init__(**kwargs)
 
-    def initialize(self, u = {}, z = {}):
+    def initialize(self, u={}, z={}):
         if u is None:
             u = {}
         if z is None:
@@ -156,7 +160,7 @@ def initialize(self, u = {}, z = {}):
             x_i = m.initialize(u_i, z_i)
             for key, value in x_i.items():
                 x_0[name + '.' + key] = value
-        
+
             # Process connections
             # This initializes the states that are connected to outputs
             for (in_key_part, in_key) in self.__to_state_connections[name]:
@@ -165,7 +169,7 @@ def initialize(self, u = {}, z = {}):
                 else:  # Missing from z, so estimate using initial state
                     z_ii = m.output(x_i)
                     x_0[in_key] = z_ii.get(in_key_part, None)
-                
+
         return self.StateContainer(x_0)
 
     def next_state(self, x, u, dt):
@@ -177,17 +181,17 @@ def next_state(self, x, u, dt):
             for (in_key, out_key_part) in self.__to_input_connections[name]:
                 u_i[out_key_part] = x[in_key]
             u_i = m.InputContainer(u_i)
-            
+
             # Prepare state
             x_i = m.StateContainer({key: x[name + '.' + key] for key in m.states})
 
-            # Propogate state
+            # Propagate state
             x_next_i = m.next_state(x_i, u_i, dt)
 
             # Save to super state
             for key, value in x_next_i.items():
                 x[name + '.' + key] = value
-            
+
             # Process connections
             # This updates the states that are connected to outputs
             if len(self.__to_state_connections[name]) > 0:
diff --git a/src/prog_models/ensemble_model.py b/src/prog_models/ensemble_model.py
index fd351674e..c554d353f 100644
--- a/src/prog_models/ensemble_model.py
+++ b/src/prog_models/ensemble_model.py
@@ -1,10 +1,9 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the
 # National Aeronautics and Space Administration.  All Rights Reserved.
 
-from collections.abc import Sequence
 import numpy as np
 
-from prog_models import PrognosticsModel
+from . import PrognosticsModel
 
 
 class EnsembleModel(PrognosticsModel):
@@ -20,7 +19,7 @@ class EnsembleModel(PrognosticsModel):
     See example :download:`examples.ensemble <../../../../prog_models/examples/ensemble.py>`
 
     Args:
-        models (list[PrognosticsModel]): List of at least 2 models that form the ensemble
+        models (list): List of models that form the ensemble
 
     Keyword Arguments:
         aggregation_method (function): Function that aggregates the outputs of the models in the ensemble. Default is np.mean
@@ -31,14 +30,6 @@ class EnsembleModel(PrognosticsModel):
     }
 
     def __init__(self, models, **kwargs):
-        if not isinstance(models, Sequence):
-            raise TypeError(f'EnsembleModel must be initialized with a list of models, got {type(models)}')
-        if len(models) < 2:
-            raise ValueError('EnsembleModel requires at least two models')
-        for i, m in enumerate(models):
-            if not isinstance(m, PrognosticsModel):
-                raise TypeError(f'EnsembleModel requires all models to be PrognosticsModel instances. models[{i}] was {type(m)}')
-
         inputs = set()
         states = set()
         outputs = set()
@@ -56,12 +47,8 @@ def __init__(self, models, **kwargs):
         super().__init__(**kwargs)
         self.parameters['models'] = models
 
-    def initialize(self, u=None, z=None):
-        xs = [
-            m.initialize(
-                m.InputContainer(u) if u is not None else None, 
-                m.OutputContainer(z) if z is not None else None
-            ) for m in self.parameters['models']]
+    def initialize(self, u, z=None):
+        xs = [m.initialize(m.InputContainer(u), m.OutputContainer(z) if z is not None else None) for m in self.parameters['models']]
         x0 = {}
         for x in xs:
             for key in x.keys():
diff --git a/src/prog_models/linear_model.py b/src/prog_models/linear_model.py
index 3c1533e95..eae7c6cee 100644
--- a/src/prog_models/linear_model.py
+++ b/src/prog_models/linear_model.py
@@ -3,12 +3,13 @@
 
 from abc import ABC, abstractmethod
 import numpy as np
-from prog_models import PrognosticsModel
 
+from prog_models.prognostics_model import PrognosticsModel
 
 class LinearModel(PrognosticsModel, ABC):
     """
     A linear prognostics :term:`model`. Used when behavior can be described using a simple linear time-series model defined by the following equations:
+
     .. math::
         \dfrac{dx}{dt} = Ax + Bu + E
 
@@ -21,8 +22,8 @@ class LinearModel(PrognosticsModel, ABC):
     Linear Models must inherit from this class and define the following properties:
         * A: 2-d np.array[float], dimensions: n_states x n_states
         * B: 2-d np.array[float], optional (zeros by default), dimensions: n_states x n_inputs
-        * C: 2-d np.array[float], dimensions: n_outputs x n_states
-        * D: 1-d np.array[float], optional (zeros by default), dimensions: n_outputs x 1
+        * C: 2-d np.array[float], optional (zeros by default), dimensions: n_outputs x n_states
+        * D: 1-d np.array[float], dimensions: n_outputs x 1
         * E: 1-d np.array[float], optional (zeros by default), dimensions: n_states x 1
         * F: 2-d np.array[float], dimensions: n_es x n_states
         * G: 1-d np.array[float], optional (zeros by default), dimensions: n_es x 1
@@ -32,80 +33,28 @@ class LinearModel(PrognosticsModel, ABC):
         * events:  list[str] - :term:`event` keys
     """
 
-    # Default Values are set to None
-    default_parameters = {
-        '_B': None,
-        '_D': None,
-        '_E': None,
-        '_G': None
-    }
-
     def __init__(self, **kwargs):
-        params = LinearModel.default_parameters.copy()
-        params.update(self.default_parameters)
-        params.update(kwargs)
-        super().__init__(**params)
-
-        # Set each property to itself
-        # This triggers the default value logic in the setter
-        # for cases where the property has not been overwritten
-        self.B = self.B
-        self.D = self.D
-        self.E = self.E
-        self.G = self.G
-
-        if self.F is None and type(self).event_state == LinearModel.event_state:
-            raise AttributeError(
-                'LinearModel must define F if event_state is not defined. Either override event_state or define F.')
-        
+        super().__init__(**kwargs)
         self.matrixCheck()
 
-# check to see if attributes are different and if functions within the models are overriden as well (i.e threshold_met and event_state)
-    def __eq__(self, other):
-        return isinstance(other, LinearModel) \
-            and np.all(self.A == other.A) \
-            and np.all(self.B == other.B) \
-            and np.all(self.C == other.C) \
-            and np.all(self.D == other.D) \
-            and np.all(self.E == other.E) \
-            and np.all(self.F == other.F) \
-            and np.all(self.G == other.G) \
-            and self.inputs == other.inputs \
-            and self.outputs == other.outputs \
-            and self.events == other.events \
-            and self.states == other.states \
-            and self.performance_metric_keys == other.performance_metric_keys \
-            and self.parameters == other.parameters \
-            and self.state_limits == other.state_limits \
-            and type(self).threshold_met == type(other).threshold_met \
-            and type(self).event_state == type(other).event_state \
-            and type(self).dx == type(other).dx \
-            and type(self).next_state == type(other).next_state \
-            and type(self).output == type(other).output \
-            and type(self).performance_metrics == type(other).performance_metrics
+        if self.F is None and type(self).event_state == LinearModel.event_state:
+            raise AttributeError('LinearModel must define F if event_state is not defined. Either override event_state or define F.')
 
     def matrixCheck(self) -> None:
         """
         Public class method for checking matrices dimensions across all properties of the model.
         """
-        self._propertyCheck(self.n_states, self.n_states,
-                            ["A", "states", "states"])
-        self._propertyCheck(self.n_states, self.n_inputs,
-                            ["B", "states", "inputs"])
-        self._propertyCheck(self.n_outputs, self.n_states,
-                            ["C", "outputs", "states"])
-        self._propertyCheck(self.n_outputs, 1,
-                            ["D", "outputs", "1"])
-        self._propertyCheck(self.n_states, 1, 
-                            ["E", "states", "1"])
-        self._propertyCheck(self.n_events, 1, 
-                            ["G", "events", "1"])
+        self._propertyCheck(self.n_states, self.n_states, ["A","states","states"])
+        self._propertyCheck(self.n_states, self.n_inputs, ["B","states","inputs"])
+        self._propertyCheck(self.n_outputs, self.n_states, ["C","outputs","states"])
+        self._propertyCheck(self.n_outputs, 1, ["D","outputs","1"])
+        self._propertyCheck(self.n_states, 1, ["E","states","1"])
+        self._propertyCheck(self.n_events, 1, ["G","events","1"])
 
         if self.F is not None:
-            self._propertyCheck(self.n_events, self.n_states, [
-                                "F", "events", "states"])
+            self._propertyCheck(self.n_events, self.n_states, ["F","events","states"])
 
-    def _propertyCheck(self, rowsCount: int, colsCount: int, notes: list) -> None:
+    def _propertyCheck(self, rowsCount : int, colsCount : int, notes : list) -> None:
         """
         matrix: Input matrix to check dimensions of (e.g. self.A, self.B, etc)
         rowsCount: Row count to check matrix against
@@ -115,19 +64,17 @@ def _propertyCheck(self, rowsCount: int, colsCount: int, notes: list) -> None:
         target_property = getattr(self, notes[0])
         if isinstance(target_property, list):
             setattr(self, notes[0], np.array(target_property))
-        matrix = getattr(self, notes[0])
+        matrix = getattr(self, notes[0]) 
         if not isinstance(matrix, np.ndarray):
-            raise TypeError(
-                "Matrix type check failed: @property {} dimensions is not of type list or NumPy array.".format(notes[0]))
+            raise TypeError("Matrix type check failed: @property {} dimensions is not of type list or NumPy array.".format(notes[0]))
 
         matrixShape = matrix.shape
-
-        if (matrix.ndim != 2 or  # Checks to see if matrix is two-dimensional
-            matrixShape[0] != rowsCount or  # checks if matrix has correct row count
-            matrixShape[1] != colsCount):  # check all rows are equal to correct column count
-            raise AttributeError(
-                "Matrix size check failed: @property {} dimensions improperly formed along {} x {}.".format(notes[0], notes[1], notes[2]))
-
+        if (matrixShape[0] != rowsCount or # check matrix is 2 dimensional, correspond to rows count
+            len(matrixShape) == 1 or # check .shape returns 2-tuple, meaning all rows are of equal length
+            matrixShape[1] != colsCount or # check all rows are equal to correct column count
+            matrix.ndim != 2): # check matrix is 2 dimensional
+            raise AttributeError("Matrix size check failed: @property {} dimensions improperly formed along {} x {}.".format(notes[0],notes[1],notes[2]))
+   
     @property
     @abstractmethod
     def A(self):
@@ -135,22 +82,7 @@ def A(self):
 
     @property
     def B(self):
-        return self.parameters['_B']
-
-    @B.setter
-    def B(self, value):
-        if (value is None):
-            self.parameters['_B'] = np.zeros((self.n_states, self.n_inputs))
-        else:
-            prev_value = self.parameters['_B']
-            self.parameters['_B'] = value
-            try:
-                self._propertyCheck(self.n_states, self.n_inputs, [
-                                    "B", "states", "inputs"])
-            except (TypeError, AttributeError) as ex:
-                # Unacceptable value, reset and re-raise
-                self.parameters['_B'] = prev_value
-                raise ex
+        return np.zeros((self.n_states, self.n_inputs))
 
     @property
     @abstractmethod
@@ -159,39 +91,11 @@ def C(self):
 
     @property
     def D(self):
-        return self.parameters['_D']
-
-    @D.setter
-    def D(self, value):
-        if (value is None):
-            self.parameters['_D'] = np.zeros((self.n_outputs, 1))
-        else:
-            prev_value = self.parameters['_D']
-            self.parameters['_D'] = value
-            try:
-                self._propertyCheck(self.n_outputs, 1, ["D", "outputs", "1"])
-            except (TypeError, AttributeError) as ex:
-                # Unacceptable value, reset and re-raise
-                self.parameters['_D'] = prev_value
-                raise ex
+        return np.zeros((self.n_outputs, 1))
 
     @property
     def E(self):
-        return self.parameters['_E']
-
-    @E.setter
-    def E(self, value):
-        if (value is None):
-            self.parameters['_E'] = np.zeros((self.n_states, 1))
-        else:
-            prev_value = self.parameters['_E']
-            self.parameters['_E'] = value
-            try:
-                self._propertyCheck(self.n_states, 1, ["E", "states", "1"])
-            except (TypeError, AttributeError) as ex:
-                # Unacceptable value, reset and re-raise
-                self.parameters['_E'] = prev_value
-                raise ex
+        return np.zeros((self.n_states, 1))
 
     @property
     @abstractmethod
@@ -200,28 +104,14 @@ def F(self):
 
     @property
     def G(self):
-        return self.parameters['_G']
-
-    @G.setter
-    def G(self, value):
-        if (value is None):
-            self.parameters['_G'] = np.zeros((self.n_events, 1))
-        else:
-            prev_value = self.parameters['_G']
-            self.parameters['_G'] = value
-            try:
-                self._propertyCheck(self.n_events, 1, ["G", "events", "1"])
-            except (TypeError, AttributeError) as ex:
-                # Unacceptable value, reset and re-raise
-                self.parameters['_G'] = prev_value
-                raise ex
+        return np.zeros((self.n_events, 1))
 
     def dx(self, x, u):
         dx_array = np.matmul(self.A, x.matrix) + self.E
         if len(u.matrix) > 0:
             dx_array += np.matmul(self.B, u.matrix)
         return self.StateContainer(dx_array)
-
+        
     def output(self, x):
         z_array = np.matmul(self.C, x.matrix) + self.D
         return self.OutputContainer(z_array)
@@ -229,7 +119,7 @@ def output(self, x):
     def event_state(self, x):
         es_array = np.matmul(self.F, x.matrix) + self.G
         return {key: value[0] for key, value in zip(self.events, es_array)}
-
+    
     def threshold_met(self, x):
         es_array = np.matmul(self.F, x.matrix) + self.G
         return {key: value[0] <= 0 for key, value in zip(self.events, es_array)}
diff --git a/src/prog_models/models/__init__.py b/src/prog_models/models/__init__.py
index aaf3eded2..b72f90790 100644
--- a/src/prog_models/models/__init__.py
+++ b/src/prog_models/models/__init__.py
@@ -1,15 +1,14 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the
 # National Aeronautics and Space Administration.  All Rights Reserved.
 
-from prog_models.models.battery_circuit import BatteryCircuit
-from prog_models.models.battery_electrochem import BatteryElectroChem, BatteryElectroChemEOD, BatteryElectroChemEOL, BatteryElectroChemEODEOL
-from prog_models.models.centrifugal_pump import CentrifugalPump, CentrifugalPumpBase, CentrifugalPumpWithWear
-from prog_models.models.pneumatic_valve import PneumaticValve, PneumaticValveBase, PneumaticValveWithWear
-from prog_models.models.dcmotor import DCMotor
-from prog_models.models.dcmotor_singlephase import DCMotorSP
-from prog_models.models.esc import ESC
-from prog_models.models.powertrain import Powertrain
-from prog_models.models.propeller_load import PropellerLoad
-from prog_models.models.thrown_object import ThrownObject
-from prog_models.models.experimental.paris_law import ParisLawCrackGrowth
-from prog_models.models.thrown_object import LinearThrownObject
+from .battery_circuit import BatteryCircuit
+from .battery_electrochem import BatteryElectroChem, BatteryElectroChemEOD, BatteryElectroChemEOL, BatteryElectroChemEODEOL
+from .centrifugal_pump import CentrifugalPump, CentrifugalPumpBase, CentrifugalPumpWithWear
+from .pneumatic_valve import PneumaticValve, PneumaticValveBase, PneumaticValveWithWear
+from .dcmotor import DCMotor
+from .dcmotor_singlephase import DCMotorSP
+from .esc import ESC
+from .powertrain import Powertrain
+from .propeller_load import PropellerLoad
+from .thrown_object import ThrownObject
+from .experimental.paris_law import ParisLawCrackGrowth
diff --git a/src/prog_models/models/battery_circuit.py b/src/prog_models/models/battery_circuit.py
index 7b3cc426f..c231d04b9 100644
--- a/src/prog_models/models/battery_circuit.py
+++ b/src/prog_models/models/battery_circuit.py
@@ -4,7 +4,7 @@
 from math import inf
 import numpy as np
 
-from prog_models import PrognosticsModel
+from .. import PrognosticsModel
 
 
 class BatteryCircuit(PrognosticsModel):
diff --git a/src/prog_models/models/battery_electrochem.py b/src/prog_models/models/battery_electrochem.py
index 7922ea9b4..379d8f69c 100644
--- a/src/prog_models/models/battery_electrochem.py
+++ b/src/prog_models/models/battery_electrochem.py
@@ -6,7 +6,7 @@
 from scipy.optimize import fsolve
 import warnings
 
-from prog_models import PrognosticsModel
+from .. import PrognosticsModel
 
 # Constants of nature
 R = 8.3144621  # universal gas constant, J/K/mol
@@ -500,7 +500,7 @@ def event_state(self, x: dict) -> dict:
         charge_EOD = (x['qnS'] + x['qnB'])/self.parameters['qnMax']
         voltage_EOD = (v - self.parameters['VEOD'])/self.parameters['VDropoff'] 
         return {
-            'EOD': np.clip(min(charge_EOD, voltage_EOD), 0, 1)
+            'EOD': min(charge_EOD, voltage_EOD)
         }
 
     def output(self, x: dict):
diff --git a/src/prog_models/models/centrifugal_pump.py b/src/prog_models/models/centrifugal_pump.py
index ac36f22df..b0905a699 100644
--- a/src/prog_models/models/centrifugal_pump.py
+++ b/src/prog_models/models/centrifugal_pump.py
@@ -5,10 +5,10 @@
 import numpy as np
 import warnings
 
-from prog_models import PrognosticsModel
+from .. import prognostics_model
 
 
-class CentrifugalPumpBase(PrognosticsModel):
+class CentrifugalPumpBase(prognostics_model.PrognosticsModel):
     """
     Prognostics :term:`model` for a Centrifugal Pump as described in [0]_.
 
diff --git a/src/prog_models/models/pneumatic_valve.py b/src/prog_models/models/pneumatic_valve.py
index a8c5a2ced..af7456247 100644
--- a/src/prog_models/models/pneumatic_valve.py
+++ b/src/prog_models/models/pneumatic_valve.py
@@ -5,7 +5,7 @@
 import numpy as np
 import warnings
 
-from prog_models import PrognosticsModel
+from .. import prognostics_model
 
 
 def calc_x(x: float, forces: float, Ls: float, new_x: float) -> float:
@@ -26,7 +26,7 @@ def calc_v(x: float, v: float, dv: float, forces: float, Ls: float, new_x: float
     return v + dv
 
 
-class PneumaticValveBase(PrognosticsModel):
+class PneumaticValveBase(prognostics_model.PrognosticsModel):
     """
     Prognostics :term:`model` for a Pneumatic Valve model as described in the following paper:
     `M. Daigle and K. Goebel, "A Model-based Prognostics Approach Applied to Pneumatic Valves," International Journal of Prognostics and Health Management, vol. 2, no. 2, August 2011. https://papers.phmsociety.org/index.php/ijphm/article/view/1359`
diff --git a/src/prog_models/models/propeller_load.py b/src/prog_models/models/propeller_load.py
index 9622a5a7b..09b45240c 100644
--- a/src/prog_models/models/propeller_load.py
+++ b/src/prog_models/models/propeller_load.py
@@ -3,7 +3,7 @@
 
 import numpy as np
 
-from prog_models import PrognosticsModel
+from .. import PrognosticsModel
 
 
 def update_Cq(params):
diff --git a/src/prog_models/models/test_models/linear_models.py b/src/prog_models/models/test_models/linear_models.py
index 14fee8535..41fec63ef 100644
--- a/src/prog_models/models/test_models/linear_models.py
+++ b/src/prog_models/models/test_models/linear_models.py
@@ -24,6 +24,7 @@ class FNoneNoEventStateLM(LinearModel):
         }
     }
 
+
 class OneInputNoOutputNoEventLM(LinearModel):
     """
     Simple model that increases state by u1 every step. 
@@ -55,7 +56,7 @@ class OneInputOneOutputNoEventLM(LinearModel):
     A = np.array([[0]])
     B = np.array([[1]])
     C = np.array([[1]])
-    F = np.empty((0, 1))
+    F = np.empty((0,1))
 
     default_parameters = {
         'process_noise': 0,
@@ -65,43 +66,6 @@ class OneInputOneOutputNoEventLM(LinearModel):
     }
 
 
-class OneInputOneOutputOneEventLM(OneInputOneOutputNoEventLM):
-    events = ['x1 == 10']
-    performance_metric_keys = ['pm1']
-
-    F = np.array([[-0.1]])
-    G = np.array([[1]])
-
-    def performance_metrics(self, x) -> dict:
-        return {'pm1': x['x1'] + 1}
-
-class OneInputOneOutputOneEventAltLM(LinearModel):
-    """
-    Simple model that increases state by u1 every step. Event occurs when state == 10
-    """
-    inputs = ['u2']
-    states = ['x2']
-    outputs = ['z2']
-    performance_metric_keys = ['pm2']
-    events = ['x2 == 5']
-
-    A = np.array([[0]])
-    B = np.array([[1]])
-    C = np.array([[1]])
-    F = np.array([[-0.2]])
-    G = np.array([[1]])
-
-    default_parameters = {
-        'process_noise': 0,
-        'x0': {
-            'x2': 0
-        }
-    }
-
-    def performance_metrics(self, x) -> dict:
-        return {'pm2': x['x2'] + 1}
-
-
 class OneInputOneOutputNoEventLMPM(OneInputOneOutputNoEventLM):
     """
     Same as OneInputOneOutputNoEventLM, but with performance metrics defined as a function. Has a single performance metric that is always the state, plus 1
@@ -133,7 +97,6 @@ class OneInputNoOutputTwoEventLM(LinearModel):
     A = np.array([[0]])
     B = np.array([[1, 0.5]])
     C = np.empty((0,1))
-    D = np.empty((0,1))
     F = np.array([[-0.1], [-0.2]])
     G = np.array([[1], [1]])
 
@@ -156,7 +119,6 @@ class TwoInputNoOutputOneEventLM(LinearModel):
     A = np.array([[0]])
     B = np.array([[1, 0.5]])
     C = np.empty((0,1))
-    D = np.empty((0,1))
     F = np.array([[-0.1]])
     G = np.array([[1]])
 
@@ -179,7 +141,6 @@ class TwoInputNoOutputTwoEventLM(LinearModel):
     A = np.array([[0]])
     B = np.array([[1, 0.5]])
     C = np.empty((0,1))
-    D = np.empty((0,1))
     F = np.array([[-0.1], [-0.2]])
     G = np.array([[1], [1]])
 
diff --git a/src/prog_models/models/test_models/linear_thrown_object.py b/src/prog_models/models/test_models/linear_thrown_object.py
index e804819c3..87dd066cd 100644
--- a/src/prog_models/models/test_models/linear_thrown_object.py
+++ b/src/prog_models/models/test_models/linear_thrown_object.py
@@ -13,8 +13,8 @@ class LinearThrownObject(LinearModel):
     events = ['impact']
 
     A = np.array([[0, 1], [0, 0]])
-    C = np.array([[1, 0]])
     E = np.array([[0], [-9.81]])
+    C = np.array([[1, 0]])
     F = None # Will override method
 
     default_parameters = {
@@ -41,91 +41,3 @@ def event_state(self, x):
             'falling': np.maximum(x['v']/self.parameters['throwing_speed'],0),  # Throwing speed is max speed
             'impact': np.maximum(x['x']/x_max,0) if x['v'] < 0 else 1  # 1 until falling begins, then it's fraction of height
         }
-
-class LinearThrownObjectNoE(LinearThrownObject):
-    E = np.array([[0], [-9.81]])
-
-class LinearThrownDiffThrowingSpeed(LinearThrownObject):
-    inputs = [] 
-    states = ['x', 'v']
-    outputs = ['x']
-    events = ['impact']
-
-    A = np.array([[0, 1], [0, 0]])
-    C = np.array([[1, 0]])
-    D = np.array([[1]])
-    E = np.array([[0], [-9.81]])
-    F = None # Will override method
-
-    default_parameters = {
-        'thrower_height': 1.83,  # m
-        'throwing_speed': 20,  # m/s
-        'g': -9.81  # Acceleration due to gravity in m/s^2
-    }
-
-class LinearThrownObjectWrongB(LinearThrownObject):
-    inputs = [] 
-    states = ['x', 'v']
-    outputs = ['x']
-    events = ['impact']
-
-    A = np.array([[0, 1], [0, 0]])
-    B = np.array([[1, 0], [0, 1]])
-    C = np.array([[1, 0]])
-    D = np.array([[1]])
-    E = np.array([[0], [-9.81]])
-    F = None # Will override method
-
-
-# Wrong x statecontainer parameter. Has Throwing_speed when it should be thrower_height
-class LinearThrownObjectUpdatedInitalizedMethod(LinearThrownObject):
-    inputs = [] 
-    states = ['x', 'v']
-    outputs = ['x']
-    events = ['impact']
-
-    A = np.array([[0, 1], [0, 0]])
-    C = np.array([[1, 0]])
-    D = np.array([[1]])
-    E = np.array([[0], [-9.81]])
-    F = None # Will override method
-
-    default_parameters = {
-        'thrower_height': 1.83,  # m
-        'throwing_speed': 40,  # m/s
-        'g': -9.81  # Acceleration due to gravity in m/s^2
-    }
-
-    def initialize(self, u=None, z=None):
-        return self.StateContainer({
-            'x': self.parameters['throwing_speed'],  # Thrown, so initial altitude is height of thrower
-            })
-
-class LinearThrownObjectDiffDefaultParams(LinearThrownObject):
-    inputs = [] 
-    states = ['x', 'v']
-    outputs = ['x']
-    events = ['impact']
-
-    A = np.array([[0, 1], [0, 0]])
-    C = np.array([[1, 0]])
-    D = np.array([[1]])
-    E = np.array([[0], [-9.81]])
-    F = None # Will override method
-
-    default_parameters = {
-        'thrower_height': 1.83,  # m
-        'throwing_speed': 40,  # m/s
-        'g': -9.81,  # Acceleration due to gravity in m/s^2
-        'x': 1111
-    }
-
-class LinearThrownObjectFourStates(LinearThrownObject):
-    inputs = []
-    states = ['x', 'v', 'y' ,'z']
-    outputs = ['x']
-    events = ['impact']
-
-    A = np.array([[0, 1, 2, 3], [0, 1, 2, 3], [0, 1, 2, 3], [0, 1, 2, 3]])
-    C = np.array([[0, 1, 2, 3]])
-    E = np.array([[0], [1], [2], [3]])
diff --git a/src/prog_models/models/thrown_object.py b/src/prog_models/models/thrown_object.py
index e551768d8..d383838ae 100644
--- a/src/prog_models/models/thrown_object.py
+++ b/src/prog_models/models/thrown_object.py
@@ -1,8 +1,10 @@
-# Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
+# Copyright © 2021 United States Government as represented by the Administrator of the
+# National Aeronautics and Space Administration.  All Rights Reserved.
 
 import numpy as np
 
-from prog_models import PrognosticsModel, LinearModel
+from .. import PrognosticsModel
+
 
 def calc_lumped_param(params):
     return {
@@ -123,41 +125,3 @@ def event_state(self, x: dict) -> dict:
             'falling': np.maximum(x['v']/self.parameters['throwing_speed'], 0),  # Throwing speed is max speed
             'impact': np.maximum(x['x']/x_max, 0)  # then it's fraction of height
         }
-
-class LinearThrownObject(LinearModel):
-    inputs = [] 
-    states = ['x', 'v']
-    outputs = ['x']
-    events = ['impact']
-
-    A = np.array([[0, 1], [0, 0]])
-    D = np.array([[1]])
-    E = np.array([[0], [-9.81]])
-    C = np.array([[1, 0]])
-    F = None # Will override method
-    
-
-    default_parameters = {
-        'thrower_height': 1.83,  # m
-        'throwing_speed': 40,  # m/s
-        'g': -9.81  # Acceleration due to gravity in m/s^2
-    }
-
-    def initialize(self, u=None, z=None):
-        return self.StateContainer({
-            'x': self.parameters['thrower_height'],  # Thrown, so initial altitude is height of thrower
-            'v': self.parameters['throwing_speed']  # Velocity at which the ball is thrown - this guy is a professional baseball pitcher
-            })
-    
-    def threshold_met(self, x):
-        return {
-            'falling': x['v'] < 0,
-            'impact': x['x'] <= 0
-        }
-
-    def event_state(self, x): 
-        x_max = x['x'] + np.square(x['v'])/(-self.parameters['g']*2) # Use speed and position to estimate maximum height
-        return {
-            'falling': np.maximum(x['v']/self.parameters['throwing_speed'],0),  # Throwing speed is max speed
-            'impact': np.maximum(x['x']/x_max,0) if x['v'] < 0 else 1  # 1 until falling begins, then it's fraction of height
-        }
diff --git a/src/prog_models/prognostics_model.py b/src/prog_models/prognostics_model.py
index 0a6d4fa9b..7522e3164 100644
--- a/src/prog_models/prognostics_model.py
+++ b/src/prog_models/prognostics_model.py
@@ -14,11 +14,9 @@
 from prog_models.exceptions import ProgModelInputException, ProgModelTypeError, ProgModelException, ProgModelStateLimitWarning
 from prog_models.sim_result import SimResult, LazySimResult
 from prog_models.utils import ProgressBar
-from prog_models.utils import calc_error
 from prog_models.utils.containers import DictLikeMatrixWrapper
-from prog_models.utils.next_state import euler_next_state, rk4_next_state, euler_next_state_wrapper, rk4_next_state_wrapper
 from prog_models.utils.parameters import PrognosticsModelParameters
-from prog_models.utils.serialization import CustomEncoder, custom_decoder
+from prog_models.utils.serialization import *
 from prog_models.utils.size import getsizeof
 
 
@@ -67,11 +65,11 @@ class PrognosticsModel(ABC):
             Limits on the state variables format {'state_name': (lower_limit, upper_limit)}
         param_callbacks : dict[str, list[function]], optional
             Callbacks for derived parameters
-        inputs: list[str], optional
+        inputs: list[str]
             Identifiers for each :term:`input`
         states: list[str]
             Identifiers for each :term:`state`
-        outputs: list[str], optional
+        outputs: list[str]
             Identifiers for each :term:`output`
         performance_metric_keys: list[str], optional
             Identifiers for each performance metric
@@ -104,9 +102,7 @@ class PrognosticsModel(ABC):
     # events = []       # Identifiers for each event
     param_callbacks = {}  # Callbacks for derived parameters
 
-    SimulationResults = namedtuple(
-        'SimulationResults', 
-        ['times', 'inputs', 'states', 'outputs', 'event_states'])
+    SimulationResults = namedtuple('SimulationResults', ['times', 'inputs', 'states', 'outputs', 'event_states'])
 
     def __init__(self, **kwargs):
         # Default params for any model
@@ -123,7 +119,7 @@ def __init__(self, **kwargs):
 
         PrognosticsModel.__setstate__(self, params)
 
-    def __eq__(self, other: "PrognosticsModel") -> bool:
+    def __eq__(self, other : "PrognosticsModel") -> bool:
         """
         Check if two models are equal
         """
@@ -135,9 +131,9 @@ def __str__(self) -> str:
     def __getstate__(self) -> dict:
         return self.parameters.data
 
-    def __setstate__(self, params: dict) -> None:
+    def __setstate__(self, params : dict) -> None:
         # This method is called when depickling and in construction. It builds the model from the parameters
-        
+
         if not hasattr(self, 'inputs'):
             self.inputs = []
         self.n_inputs = len(self.inputs)
@@ -167,9 +163,7 @@ def __setstate__(self, params: dict) -> None:
         self.n_performance = len(self.performance_metric_keys)
 
         # Setup Containers
-        # These containers should be used instead of dictionaries for models 
-        # that use the internal matrix state
-
+        # These containers should be used instead of dictionaries for models that use the internal matrix state
         states = self.states
 
         class StateContainer(DictLikeMatrixWrapper):
@@ -193,7 +187,7 @@ def __init__(self, data):
 
         self.parameters = PrognosticsModelParameters(self, params, self.param_callbacks)
 
-    def initialize(self, u=None, z=None):
+    def initialize(self, u = None, z = None):
         """
         Calculate initial state given inputs and outputs. If not defined for a model, it will return parameters['x0']
 
@@ -251,7 +245,7 @@ def apply_measurement_noise(self, z):
         z.matrix += np.random.normal(0, self.parameters['measurement_noise'].matrix, size=z.matrix.shape)
         return z
 
-    def apply_process_noise(self, x, dt: float = 1):
+    def apply_process_noise(self, x, dt : int = 1):
         """
         Apply process noise to the state
 
@@ -322,7 +316,7 @@ def dx(self, x, u):
         """
         raise ProgModelException('dx not defined - please use next_state()')
 
-    def next_state(self, x, u, dt: float):
+    def next_state(self, x, u, dt : float):
         """
         State transition equation: Calculate next state
 
@@ -414,6 +408,51 @@ def apply_limits(self, x):
                 x[key] = np.minimum(x[key], limit[1])
         return x
 
+    def __next_state(self, x, u, dt : float):
+        """
+        State transition equation: Calls next_state(), calculating the next state, and then adds noise
+
+        Parameters
+        ----------
+        x : StateContainer
+            state, with keys defined by model.states \n
+            e.g., x = m.StateContainer({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
+        u : InputContainer
+            Inputs, with keys defined by model.inputs \n
+            e.g., u = m.InputContainer({'i':3.2}) given inputs = ['i']
+        dt : float
+            Timestep size in seconds (≥ 0) \n
+            e.g., dt = 0.1
+
+        Returns
+        -------
+        x : StateContainer
+            Next state, with keys defined by model.states
+            e.g., x = m.StateContainer({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
+
+        Example
+        -------
+        | m = PrognosticsModel() # Replace with specific model being simulated
+        | u = m.InputContainer({'u1': 3.2})
+        | z = m.OutputContainer({'z1': 2.2})
+        | x = m.initialize(u, z) # Initialize first state
+        | x = m.__next_state(x, u, 0.1) # Returns state, with noise, at 3.1 seconds given input u
+
+        See Also
+        --------
+        next_state
+
+        Note
+        ----
+        A model should not overwrite '__next_state'
+        A model should overwrite either `next_state` or `dx`. Override `dx` for continuous models, and `next_state` for discrete, where the behavior cannot be described by the first derivative.
+        """
+        # Calculate next state and add process noise
+        next_state = self.apply_process_noise(self.next_state(x, u, dt), dt)
+
+        # Apply Limits
+        return self.apply_limits(next_state)
+
     def performance_metrics(self, x) -> dict:
         """
         Calculate performance metrics where
@@ -677,7 +716,7 @@ def time_of_event(self, x, future_loading_eqn = lambda t,x=None: {}, **kwargs) -
             t = result.times[-1]
         return time_of_event
 
-    def simulate_to(self, time : float, future_loading_eqn: Callable = lambda t,x=None: {}, first_output=None, **kwargs) -> namedtuple:
+    def simulate_to(self, time : float, future_loading_eqn : Callable = lambda t,x=None: {}, first_output = None, **kwargs) -> namedtuple:
         """
         Simulate prognostics model for a given number of seconds
 
@@ -718,10 +757,10 @@ def simulate_to(self, time : float, future_loading_eqn: Callable = lambda t,x=No
         Example
         -------
         >>> def future_load_eqn(t):
-        >>>    if t< 5.0: # Load is 3.0 for first 5 seconds
-        >>>        return 3.0
-        >>>    else:
-        >>>        return 5.0
+        >>>     if t< 5.0: # Load is 3.0 for first 5 seconds
+        >>>         return 3.0
+        >>>     else:
+        >>>         return 5.0
         >>> first_output = m.OutputContainer({'o1': 3.2, 'o2': 1.2})
         >>> m = PrognosticsModel() # Replace with specific model being simulated
         >>> (times, inputs, states, outputs, event_states) = m.simulate_to(200, future_load_eqn, first_output)
@@ -739,7 +778,7 @@ def simulate_to(self, time : float, future_loading_eqn: Callable = lambda t,x=No
 
         return self.simulate_to_threshold(future_loading_eqn, first_output, **kwargs)
  
-    def simulate_to_threshold(self, future_loading_eqn: Callable = None, first_output = None, threshold_keys: list = None, **kwargs) -> namedtuple:
+    def simulate_to_threshold(self, future_loading_eqn : Callable = None, first_output = None, threshold_keys : list = None, **kwargs) -> namedtuple:
         """
         Simulate prognostics model until any or specified threshold(s) have been met
 
@@ -883,6 +922,7 @@ def simulate_to_threshold(self, future_loading_eqn: Callable = None, first_outpu
             x = self.StateContainer(x)
         
         # Optimization
+        next_state = self.__next_state
         output = self.__output
         threshold_met_eqn = self.threshold_met
         event_state = self.event_state
@@ -934,6 +974,7 @@ def check_thresholds(thresholds_met):
         if config['print']:
             def update_all():
                 saved_times.append(t)
+                print('u, saved_input: ', u)
                 saved_inputs.append(u)
                 saved_states.append(deepcopy(x))  # Avoid optimization where x is not copied
                 saved_outputs.append(output(x))
@@ -986,6 +1027,19 @@ def load_eqn(t, x):
                 u = future_loading_eqn(t, x)
                 return self.InputContainer(u)
 
+        if not isinstance(self.next_state(x.copy(), u, dt0), DictLikeMatrixWrapper):
+            # Wrapper around next_state
+            def next_state(x, u, dt):
+                # Calculate next state, and convert
+                x_new = self.next_state(x, u, dt)
+                x_new = self.StateContainer(x_new)
+
+                # Calculate next state and add process noise
+                next_state = self.apply_process_noise(x_new, dt)
+
+                # Apply Limits
+                return self.apply_limits(next_state)
+
         if not isinstance(self.output(x), DictLikeMatrixWrapper):
             # Wrapper around the output equation
             def output(x):
@@ -1014,16 +1068,21 @@ def output(x):
             apply_limits = self.apply_limits
             apply_process_noise = self.apply_process_noise
             StateContainer = self.StateContainer
-            if not isinstance(self.dx(x.copy(), u), DictLikeMatrixWrapper):
-                next_state = rk4_next_state
-            else:
-                next_state = rk4_next_state_wrapper
-        elif config['integration_method'].lower() == 'euler':
-            if not isinstance(self.next_state(x.copy(), u, dt0), DictLikeMatrixWrapper):
-                next_state = euler_next_state_wrapper
-            else:
-                next_state = euler_next_state
-        else:
+            def next_state(x, u, dt):
+                dx1 = StateContainer(dx(x, u))
+                
+                x2 = StateContainer({key: x[key] + dt*dx_i/2 for key, dx_i in dx1.items()})
+                dx2 = dx(x2, u)
+
+                x3 = StateContainer({key: x[key] + dt*dx_i/2 for key, dx_i in dx2.items()})
+                dx3 = dx(x3, u)
+                
+                x4 = StateContainer({key: x[key] + dt*dx_i for key, dx_i in dx3.items()})   
+                dx4 = dx(x4, u)
+
+                x = StateContainer({key: x[key]+ dt/3*(dx1[key]/2 + dx2[key] + dx3[key] + dx4[key]/2) for key in dx1.keys()})
+                return apply_limits(apply_process_noise(x))
+        elif config['integration_method'].lower() != 'euler':
             raise ProgModelInputException(f"'integration_method' mode {config['integration_method']} not supported. Must be 'euler' or 'rk4'")
        
         while t < horizon:
@@ -1033,7 +1092,7 @@ def output(x):
             # This is sometimes referred to as 'leapfrog integration'
             u = load_eqn(t, x)
             t = t + dt/2
-            x = next_state(self, x, u, dt)
+            x = next_state(x, u, dt)
 
             # Save if at appropriate time
             if (t >= next_save):
@@ -1085,7 +1144,7 @@ def output(x):
     def __sizeof__(self):
         return getsizeof(self)
 
-    def calc_error(self, times: List[float], inputs: List[dict], outputs: List[dict], **kwargs) -> float:
+    def calc_error(self, times : List[float], inputs : List[dict], outputs : List[dict], **kwargs) -> float:
         """Calculate Mean Squared Error (MSE) between simulated and observed
 
         Args:
@@ -1094,47 +1153,51 @@ def calc_error(self, times: List[float], inputs: List[dict], outputs: List[dict]
             outputs (list[dict]): array of output dictionaries where output[x] corresponds to time[x]
         
         Keyword Args:
-            method (str, optional): Error method to use. Supported methods include:
-                * MSE (Mean Squared Error)
-                * RMSE (Root Mean Squared Error)
-                * MAX_E (Maximum Error)
-                * MAE (Mean Absolute Error)
-                * MAPE (Mean Absolute Percentage Error)
             x0 (dict, optional): Initial state
-            dt (float, optional): Minimum time step in simulation. Defaults to 1e99.
-            stability_tol (double, optional): Configurable parameter.
-                Configurable cutoff value, between 0 and 1, that determines the fraction of the data points for which the model must be stable.
-                In some cases, a prognostics model will become unstable under certain conditions, after which point the model can no longer represent behavior. 
-                stability_tol represents the fraction of the provided argument `times` that are required to be met in simulation, 
-                before the model goes unstable in order to produce a valid estimate of mean squared error. 
-
-                If the model goes unstable before stability_tol is met, NaN is returned. 
-                Else, model goes unstable after stability_tol is met, the mean squared error calculated from data up to the instability is returned.
+            dt (double, optional): time step
 
         Returns:
-            float: error
-
-        See Also:
-            :func:`calc_error.MSE`
+            double: Total error
         """
-        method = kwargs.get('method', 'MSE')
-
-        # Call appropriate error calculation method
-        if method.lower() == 'mse':
-            return calc_error.MSE(self, times, inputs, outputs, **kwargs)
-        if method.lower() == 'max_e':
-            return calc_error.MAX_E(self, times, inputs, outputs, **kwargs)
-        if method.lower() == 'rmse':
-            return calc_error.RMSE(self, times, inputs, outputs, **kwargs)
-        if method.lower() == 'mae':
-            return calc_error.MAE(self, times, inputs, outputs, **kwargs)
-        if method.lower() == 'mape':
-            return calc_error.MAPE(self, times, inputs, outputs, **kwargs)
-
-        # If we get here, method is not supported
-        raise ProgModelInputException(f"Error method '{method}' not supported")
+        if isinstance(times[0], Iterable):
+            # Calculate error for each
+            error = [self.calc_error(t, i, z, **kwargs) for (t, i, z) in zip(times, inputs, outputs)]
+            return sum(error)/len(error)
+
+        x = kwargs.get('x0', self.initialize(inputs[0], outputs[0]))
+        dt = kwargs.get('dt', 1e99)
+
+        if not isinstance(x, self.StateContainer):
+            x = [self.StateContainer(x_i) for x_i in x]
+
+        if not isinstance(inputs[0], self.InputContainer):
+            inputs = [self.InputContainer(u_i) for u_i in inputs]
+        
+        if not isinstance(outputs[0], self.OutputContainer):
+            outputs = [self.OutputContainer(z_i) for z_i in outputs]
+
+        counter = 0  # Needed to account for skipped (i.e., none) values
+        t_last = times[0]
+        err_total = 0
+        z_obs = self.output(x)  # Initialize
+        for t, u, z in zip(times, inputs, outputs):
+            while t_last < t:
+                t_new = min(t_last + dt, t)
+                x = self.next_state(x, u, t_new-t_last)
+                t_last = t_new
+                if t >= t_last:
+                    # Only recalculate if required
+                    z_obs = self.output(x)
+            if not (None in z_obs.matrix or None in z.matrix):
+                if any(np.isnan(z_obs.matrix)):
+                    warn("Model unstable- NaN reached in simulation (t={})".format(t))
+                    break
+                err_total += np.sum(np.square(z.matrix - z_obs.matrix), where= ~np.isnan(z.matrix))
+                counter += 1
+
+        return err_total/counter
     
-    def estimate_params(self, runs: List[tuple] = None, keys: List[str] = None, times = None, inputs = None, outputs = None, method = 'nelder-mead', **kwargs) -> None:
+    def estimate_params(self, runs : List[tuple] = None, keys : List[str] = None, times = None, inputs = None, outputs = None, **kwargs) -> None:
         """Estimate the model parameters given data. Overrides model parameters
 
         Keyword Args:
@@ -1148,8 +1211,6 @@ def estimate_params(self, runs: List[tuple] = None, keys: List[str] = None, time
                 Array of output containers where output[x] corresponds to time[x]
             method (str, optional): 
                 Optimization method- see scipy.optimize.minimize for options
-            error_method (str, optional):
-                Method to use in calculating error. See calc_error for options
             bounds (tuple or dict): 
                 Bounds for optimization in format ((lower1, upper1), (lower2, upper2), ...) or {key1: (lower1, upper1), key2: (lower2, upper2), ...}
             options (dict): 
@@ -1166,9 +1227,9 @@ def estimate_params(self, runs: List[tuple] = None, keys: List[str] = None, time
             keys = [key for key in self.parameters.keys() if isinstance(self.parameters[key], Number)]
 
         config = {
+            'method': 'nelder-mead',
             'bounds': tuple((-np.inf, np.inf) for _ in keys),
             'options': {'xatol': 1e-8},
-            'error_method': 'MSE'
         }
         config.update(kwargs)
 
@@ -1220,7 +1281,7 @@ def optimization_fcn(params):
             err = 0
             for run in runs:
                 try:
-                    err += self.calc_error(run[0], run[1], run[2], method = config['error_method'], **kwargs)
+                    err += self.calc_error(run[0], run[1], run[2], **kwargs)
                 except Exception:
                     return 1e99 
                     # If it doesn't work (i.e., throws an error), don't use it
@@ -1228,7 +1289,7 @@ def optimization_fcn(params):
         
         params = np.array([self.parameters[key] for key in keys])
 
-        res = minimize(optimization_fcn, params, method=method, bounds = config['bounds'], options=config['options'])
+        res = minimize(optimization_fcn, params, method=config['method'], bounds = config['bounds'], options=config['options'])
         for x, key in zip(res.x, keys):
             self.parameters[key] = x
 
diff --git a/src/prog_models/sim_result.py b/src/prog_models/sim_result.py
index 5b4297da3..cfdbaa40e 100644
--- a/src/prog_models/sim_result.py
+++ b/src/prog_models/sim_result.py
@@ -4,10 +4,11 @@
 from copy import deepcopy
 from matplotlib.pyplot import figure
 import numpy as np
-from typing import Callable, Dict, List
+import pandas as pd
+from typing import Callable, Dict, List, Union
 
-from prog_models.utils.containers import DictLikeMatrixWrapper
-from prog_models.visualize import plot_timeseries
+from .utils.containers import DictLikeMatrixWrapper
+from .visualize import plot_timeseries
 
 
 class SimResult(UserList):
@@ -17,22 +18,34 @@ class SimResult(UserList):
     Args:
             times (array[float]): Times for each data point where times[n] corresponds to data[n]
             data (array[Dict[str, float]]): Data points where data[n] corresponds to times[n]
+            frame (pd.DataFrame): all times and data sorted in a pandas DataFrame
     """
 
-    __slots__ = ['times', 'data']  # Optimization 
-    
-    def __init__(self, times : list = None, data : list = None, _copy = True):
+    __slots__ = ['times', 'data', 'frame']  # Optimization
+
+    def __init__(self, times: list[float] = None, data: list[Union[DictLikeMatrixWrapper, dict]] = None, _copy=True):
+        # empty lists are passed
         if times is None or data is None:
-            self.times = [] 
+            self.times = []
             self.data = []
+            self.frame = pd.DataFrame()
         else:
             self.times = times.copy()
             if _copy:
                 self.data = deepcopy(data)
             else:
                 self.data = data
+            # creating multi row pd.DataFrame from data list of dict
+            if len(self.data) > 0:
+                self.frame = pd.concat([
+                    pd.DataFrame(dict(dframe), index=[0]) for dframe in self.data
+                ], ignore_index=True, axis=0)
+            else:
+                self.frame = pd.DataFrame(self.data)
+            self.frame.insert(0, "time", self.times)
+            self.frame.reindex()
 
-    def __eq__(self, other : "SimResult") -> bool:
+    def __eq__(self, other: "SimResult") -> bool:
         """Compare 2 SimResults
 
         Args:
@@ -41,11 +54,14 @@ def __eq__(self, other : "SimResult") -> bool:
         Returns:
             bool: If the two SimResults are equal
         """
-        return self.times == other.times and self.data == other.data
+        time_check = self.times == other.times
+        data_check = (self.data == other.data)
+        frame_check = self.frame.equals(other.frame)
+        return time_check and data_check and frame_check
 
-    def index(self, other : dict, *args, **kwargs) -> int:
+    def index(self, other: dict, *args, **kwargs) -> int:
         """
-        Get the index of the first sample where other occurs
+        Get the index of the first location where other occurs
 
         Args:
             other (dict)
@@ -55,7 +71,7 @@ def index(self, other : dict, *args, **kwargs) -> int:
         """
         return self.data.index(other, *args, **kwargs)
 
-    def extend(self, other : "SimResult") -> None:
+    def extend(self, other: "SimResult") -> None:
         """
         Extend the SimResult with another SimResult or LazySimResult object
 
@@ -66,10 +82,11 @@ def extend(self, other : "SimResult") -> None:
         if isinstance(other, SimResult):
             self.times.extend(other.times)
             self.data.extend(other.data)
+            self.frame = pd.concat([self.frame, other.frame], ignore_index=True, axis=0)
         else:
             raise ValueError(f"ValueError: Argument must be of type {self.__class__}")
 
-    def pop(self, index : int = -1) -> dict:
+    def pop(self, index: int = -1) -> dict:
         """Remove and return an element
 
         Args:
@@ -78,10 +95,18 @@ def pop(self, index : int = -1) -> dict:
         Returns:
             dict: Element Removed
         """
+        if index is -1:
+            index_df = len(self.frame.index) - 1
+        else:
+            index_df = index
         self.times.pop(index)
+        temp_df = self.frame.T
+        temp_df.pop(index_df)
+        self.frame = temp_df.T
+        self.frame = self.frame.reset_index(drop=True)
         return self.data.pop(index)
 
-    def remove(self, d : dict = None, t : float = None) -> None:
+    def remove(self, d: Union[float, dict] = None, t: float = None) -> None:
         """Remove an element
 
         Args:
@@ -91,19 +116,25 @@ def remove(self, d : dict = None, t : float = None) -> None:
         if sum([i is None for i in (d, t)]) != 1:
             raise ValueError("ValueError: Only one named argument (d, t) can be specified.")
 
-        if (t is not None):
-            self.data.pop(self.times.index(t))
-            self.times.remove(t)
+        if t is not None:
+            # removes the index of the timestamp meant for removal
+            index = self.times.index(t)
+            self.pop(index)
         else:
-            self.times.pop(self.data.index(d))
-            self.data.remove(d)
+            # finds index of dictionary meant to be removed
+            if isinstance(d, dict) and len(d) > 0:
+                index = self.data.index(d)
+                self.pop(index)
+            else:
+                raise ValueError
 
     def clear(self) -> None:
         """Clear the SimResult"""
         self.times = []
         self.data = []
+        self.frame = pd.DataFrame()
 
-    def time(self, index : int) -> float:
+    def get_time(self, index: int) -> float:
         """Get time for data point at index `index`
 
         Args:
@@ -112,27 +143,23 @@ def time(self, index : int) -> float:
         Returns:
             float: Time for which the data point at index `index` corresponds
         """
+        # return self.frame['time'].iat[index]
         return self.times[index]
 
-    def to_numpy(self, keys = None) -> np.ndarray:
+    def to_numpy(self, keys: list = None) -> np.ndarray:
         """
         Convert from simresult to numpy array
-
         Args:
             keys: Subset of keys to return as part of numpy array (by default, all)
-
         Returns:
             np.ndarray: numpy array representing simresult
         """
-        if len(self.data) == 0:
+        if len(self.data) is 0:
             return np.array([[]], dtype=np.float64)
-        if len(self.data[0]) == 0:
-            return np.array([[] for _ in self.data], dtype=np.float64)
-        if isinstance(self.data[0], DictLikeMatrixWrapper) and keys is None:
-            return np.array([u_i.matrix[:, 0] for u_i in self.data], dtype=np.float64)
-        if keys is None:
-            keys = self.data[0].keys()
-        return np.array([[u_i[key] for key in keys] for u_i in self.data], dtype=np.float64)
+        if keys is not None:
+            with_keys_numpy = self.frame.drop(['time'], axis=1)[keys].to_numpy(dtype=np.float64)
+            return with_keys_numpy
+        return self.frame.drop(['time'], axis=1).to_numpy(dtype=np.float64)
 
     def plot(self, **kwargs) -> figure:
         """
@@ -154,7 +181,7 @@ def plot(self, **kwargs) -> figure:
         Returns:
             Figure
         """
-        return plot_timeseries(self.times, self.data, legend = {'display': True}, options=kwargs)
+        return plot_timeseries(self.times, self.data, legend={'display': True}, options=kwargs)
 
     def monotonicity(self) -> Dict[str, float]:
         """
@@ -170,7 +197,6 @@ def monotonicity(self) -> Dict[str, float]:
         https://www.sciencedirect.com/science/article/pii/S0004370222000078
 
         Args:
-            None
 
         Returns:
             float: Value between [0, 1] indicating monotonicity of a given event for the Prediction.
@@ -185,11 +211,26 @@ def monotonicity(self) -> Dict[str, float]:
         result = {}
         for key, l in by_event.items():
             mono_sum = 0
-            for i in range(len(l)-1): 
-                mono_sum += np.sign(l[i+1] - l[i])
-            result[key] = abs(mono_sum / (len(l)-1))
+            for i in range(len(l) - 1):
+                mono_sum += np.sign(l[i + 1] - l[i])
+            result[key] = abs(mono_sum / (len(l) - 1))
         return result
 
+    def monotonicity_df(self) -> pd.DataFrame:
+        """
+        Returns:
+            pd.DataFrame: values from monotonicity in a DataFrame
+        """
+        # creating multi row pd.DataFrame from data list of dict
+        mono_dict = self.monotonicity()
+        if len(mono_dict) > 0:
+            mono_df = pd.DataFrame(mono_dict, index=[0])
+        else:
+            mono_df = pd.DataFrame()
+        return mono_df
+
+
+
     def __not_implemented(self):  # lgtm [py/inheritance/signature-mismatch]
         raise NotImplementedError("Not Implemented")
 
@@ -207,7 +248,8 @@ class LazySimResult(SimResult):  # lgtm [py/missing-equals]
     """
     Used to store the result of a simulation, which is only calculated on first request
     """
-    def __init__(self, fcn : Callable, times : list = None, states : list = None, _copy = True) -> None:
+
+    def __init__(self, fcn: Callable, times: list = None, states: list = None, _copy: bool = True) -> None:
         """
         Args:
             fcn (callable): function (x) -> z where x is the state and z is the data
@@ -219,15 +261,27 @@ def __init__(self, fcn : Callable, times : list = None, states : list = None, _c
         if times is None or states is None:
             self.times = []
             self.states = []
+            self.frame = pd.DataFrame()
         else:
             self.times = times.copy()
             if _copy:
                 self.states = deepcopy(states)
+
             else:
                 self.states = states
+            if len(self.states) > 0:  # BOOKMARK
+                # state DataFrame
+                self.frame = pd.concat([
+                    pd.DataFrame(dict(dframe), index=[0]) for dframe in self.states
+                ], ignore_index=True, axis=0)
+                # inserting time column
+                self.frame.insert(0, "time", self.times)
+                self.frame.reindex()
+            else:
+                self.frame = pd.DataFrame()
 
     def __reduce__(self):
-        return (self.__class__.__base__, (self.times, self.data))
+        return self.__class__.__base__, (self.times, self.data)
 
     def is_cached(self) -> bool:
         """
@@ -243,33 +297,38 @@ def clear(self) -> None:
         self.times = []
         self.__data = None
         self.states = []
+        self.frame = pd.DataFrame()
 
-    def extend(self, other : "LazySimResult", _copy=True) -> None:
+    def extend(self, other: "LazySimResult", _copy: bool = True) -> None:
         """
         Extend the LazySimResult with another LazySimResult object
         Raise ValueError if SimResult is passed
         Function fcn of other LazySimResult MUST match function fcn of LazySimResult object to be extended
 
         Args:
-            other (LazySimResult)
+            _copy: bool
+            other: (LazySimResult)
 
         """
-        if (isinstance(other, self.__class__)):
+        if isinstance(other, self.__class__):
             self.times.extend(other.times)
             if _copy:
                 self.states.extend(deepcopy(other.states))
+                self.frame = pd.concat([self.frame, deepcopy(other.frame)], ignore_index=True, axis=0)
             else:
                 self.states.extend(other.states)
+                self.frame = pd.concat([self.frame, other.frame], ignore_index=True, axis=0)
             if self.__data is None or not other.is_cached():
                 self.__data = None
             else:
                 self.__data.extend(other.data)
-        elif (isinstance(other, SimResult)):
-            raise ValueError(f"ValueError: {self.__class__} cannot be extended by SimResult. First convert to SimResult using to_simresult() method.")
+        elif isinstance(other, SimResult):
+            raise ValueError(
+                f"ValueError: {self.__class__} cannot be extended by SimResult. First convert to SimResult using to_simresult() method.")
         else:
             raise ValueError(f"ValueError: Argument must be of type {self.__class__}.")
 
-    def pop(self, index : int = -1) -> dict:
+    def pop(self, index: int = -1) -> dict:
         """Remove an element. If data hasn't been cached, remove the state - so it wont be calculated
 
         Args:
@@ -279,12 +338,21 @@ def pop(self, index : int = -1) -> dict:
             dict: Element Removed
         """
         self.times.pop(index)
+        # to pop from self.frame
+        if index is -1:
+            index_df = len(self.frame.index) - 1
+        else:
+            index_df = index
+        temp_df = self.frame.T
+        temp_df.pop(index_df)
+        self.frame = temp_df.T
+        self.frame = self.frame.reset_index(drop=True)
         x = self.states.pop(index)
         if self.__data is not None:
             return self.__data.pop(index)
         return self.fcn(x)
 
-    def remove(self, d : float = None, t : float = None, s = None) -> None:
+    def remove(self, d: float = None, t: float = None, s=None) -> None:
         """Remove an element
 
         Args:
@@ -294,24 +362,14 @@ def remove(self, d : float = None, t : float = None, s = None) -> None:
         """
         if sum([i is None for i in (d, t, s)]) != 2:
             raise ValueError("ValueError: Only one named argument (d, t, s) can be specified.")
-
-        if (t is not None):
-            target_index = self.times.index(t)
-            self.times.pop(target_index)
-            self.states.pop(target_index)
-            if self.__data is not None:
-                self.__data.pop(target_index)
-        elif (s is not None):
-            target_index = self.states.index(s)
-            self.times.pop(target_index)
-            self.states.pop(target_index)
-            if self.__data is not None:
-                self.__data.pop(target_index)
+        # get index value
+        if t is not None:
+            index = self.times.index(t)
+        elif s is not None:
+            index = self.states.index(s)
         else:
-            target_index = self.data.index(d)
-            self.times.pop(target_index)
-            self.states.pop(target_index)
-            self.__data.pop(target_index)
+            index = self.data.index(d)
+        self.pop(index)
 
     def to_simresult(self) -> SimResult:
         return SimResult(self.times, self.data)
@@ -327,3 +385,23 @@ def data(self) -> List[dict]:
         if self.__data is None:
             self.__data = [self.fcn(x) for x in self.states]
         return self.__data
+
+    def get_frame_data(self) -> pd.DataFrame:
+        """
+        place fcn data (list[dict]) into a pd.DataFrame format.
+
+        Returns:
+            pd.DataFrame: frame
+        """
+        if len(self.data) == 0:
+            return pd.DataFrame()
+        else:
+            # creating fcn(x) DataFrame
+            # fcn data DataFrame
+            frame = pd.concat([
+                pd.DataFrame(dict(dframe), index=[0]) for dframe in self.data
+            ], ignore_index=True, axis=0)
+            # inserting time column
+            frame.insert(0, "time", self.times)
+            frame.reindex()
+            return frame
diff --git a/src/prog_models/utils/calc_error.py b/src/prog_models/utils/calc_error.py
deleted file mode 100644
index 86c413093..000000000
--- a/src/prog_models/utils/calc_error.py
+++ /dev/null
@@ -1,363 +0,0 @@
-# Copyright © 2021 United States Government as represented by the Administrator of the
-# National Aeronautics and Space Administration.  All Rights Reserved.
-
-"""
-This file contains functions for calculating error given a model and some data (times, inputs, outputs). This is used by the PrognosticsModel.calc_error() method.
-"""
-
-from collections.abc import Iterable
-from warnings import warn
-import math
-import numpy as np
-
-
-def MAX_E(m, times, inputs, outputs, **kwargs):
-    """
-    Calculate the Maximum Error between model behavior and some collected data.
-
-    Args:
-        m (PrognosticsModel): Model to use for comparison
-        times (list[float]): array of times for each sample
-        inputs (list[dict]): array of input dictionaries where input[x] corresponds to time[x]
-        outputs (list[dict]): array of output dictionaries where output[x] corresponds to time[x]
-
-    Keyword Args:
-        x0 (StateContainer): Current State of the model
-        dt (float, optional): Minimum time step in simulation. Defaults to 1e99.
-        stability_tol (double, optional): Configurable parameter.
-            Configurable cutoff value, between 0 and 1, that determines the fraction of the data points for which the model must be stable.
-            In some cases, a prognostics model will become unstable under certain conditions, after which point the model can no longer represent behavior. 
-            stability_tol represents the fraction of the provided argument `times` that are required to be met in simulation, 
-            before the model goes unstable in order to produce a valid estimate of mean squared error. 
-
-            If the model goes unstable before stability_tol is met, NaN is returned. 
-            Else, model goes unstable after stability_tol is met, the mean squared error calculated from data up to the instability is returned.
-
-    Returns:
-        float: Maximum error between model and data
-    """
-    if isinstance(times[0], Iterable):
-        # Calculate error for each
-        error = [MAX_E(t, i, z, **kwargs) for (t, i, z) in zip(times, inputs, outputs)]
-        return max(error)
-
-    x = kwargs.get('x0', m.initialize(inputs[0], outputs[0]))
-    dt = kwargs.get('dt', 1e99)
-    stability_tol = kwargs.get('stability_tol', 0.95)
-
-    if not isinstance(x, m.StateContainer):
-        x = m.StateContainer(x)
-
-    if not isinstance(inputs[0], m.InputContainer):
-        inputs = [m.InputContainer(u_i) for u_i in inputs]
-    
-    if not isinstance(outputs[0], m.OutputContainer):
-        outputs = [m.OutputContainer(z_i) for z_i in outputs]
-
-    # Checks stability_tol is within bounds
-    # Throwing a default after the warning.
-    if stability_tol >= 1 or stability_tol < 0:
-        warn(f"configurable cutoff must be some float value in the domain (0, 1]. Received {stability_tol}. Resetting value to 0.95")
-        stability_tol = 0.95
-
-    counter = 0
-    t_last = times[0]
-    err_max = 0
-    z_obs = m.output(x)  # Initialize
-    cutoffThreshold = math.floor(stability_tol * len(times))
-
-    for t, u, z in zip(times, inputs, outputs):
-        while t_last < t:
-            t_new = min(t_last + dt, t)
-            x = m.next_state(x, u, t_new-t_last)
-            t_last = t_new
-            if t >= t_last:
-                # Only recalculate if required
-                z_obs = m.output(x)
-        if not (None in z_obs.matrix or None in z.matrix):
-            # The none check above is used to cover the case where the model
-            # is not able to produce an output for a given input yet
-            # For example, in LSTM models, the first few inputs will not 
-            # produce an output until the model has received enough data
-            # This is true for any window-based model
-            if any(np.isnan(z_obs.matrix)):
-                if counter < cutoffThreshold:
-                    raise ValueError(f"Model unstable- NAN reached in simulation (t={t}) before cutoff threshold. Cutoff threshold is {cutoffThreshold}, or roughly {stability_tol * 100}% of the data")                
-                else: 
-                    warn(f"Model unstable- NaN reached in simulation (t={t})")
-                    break
-            err_max = max(err_max, np.max(
-                np.abs(z.matrix - z_obs.matrix)
-            ))
-            counter += 1
-
-    if counter == 0:
-        return np.nan
-
-    return err_max
-
-
-def RMSE(m, times, inputs, outputs, **kwargs):
-    """
-    Calculate the Root Mean Squared Error between model behavior and some collected data.
-
-    Args:
-        m (PrognosticsModel): Model to use for comparison
-        times (list[float]): array of times for each sample
-        inputs (list[dict]): array of input dictionaries where input[x] corresponds to time[x]
-        outputs (list[dict]): array of output dictionaries where output[x] corresponds to time[x]
-
-    Keyword Args:
-        x0 (StateContainer): Current State of the model
-        dt (float, optional): Minimum time step in simulation. Defaults to 1e99.
-        stability_tol (double, optional): Configurable parameter.
-            Configurable cutoff value, between 0 and 1, that determines the fraction of the data points for which the model must be stable.
-            In some cases, a prognostics model will become unstable under certain conditions, after which point the model can no longer represent behavior. 
-            stability_tol represents the fraction of the provided argument `times` that are required to be met in simulation, 
-            before the model goes unstable in order to produce a valid estimate of mean squared error. 
-
-            If the model goes unstable before stability_tol is met, NaN is returned. 
-            Else, model goes unstable after stability_tol is met, the mean squared error calculated from data up to the instability is returned.
-
-    Returns:
-        float: RMSE between model and data
-    """
-    return np.sqrt(MSE(m, times, inputs, outputs, **kwargs))
-
-
-def MSE(self, times, inputs, outputs, **kwargs) -> float:
-    """Calculate Mean Squared Error (MSE) between simulated and observed
-
-    Args:
-        times (list[float]): Array of times for each sample.
-        inputs (list[dict]): Array of input dictionaries where input[x] corresponds to time[x].
-        outputs (list[dict]): Array of output dictionaries where output[x] corresponds to time[x].
-
-    Keyword Args:
-        x0 (dict, optional): Initial state.
-        dt (double, optional): Maximum time step.
-        stability_tol (double, optional): Configurable parameter.
-            Configurable cutoff value, between 0 and 1, that determines the fraction of the data points for which the model must be stable.
-            In some cases, a prognostics model will become unstable under certain conditions, after which point the model can no longer represent behavior. 
-            stability_tol represents the fraction of the provided argument `times` that are required to be met in simulation, 
-            before the model goes unstable in order to produce a valid estimate of mean squared error. 
-
-            If the model goes unstable before stability_tol is met, NaN is returned. 
-            Else, model goes unstable after stability_tol is met, the mean squared error calculated from data up to the instability is returned.
-
-    Returns:
-        double: Total error
-    """
-    if isinstance(times[0], Iterable):
-        # Calculate error for each
-        error = [self.calc_error(t, i, z, **kwargs) for (t, i, z) in zip(times, inputs, outputs)]
-        return sum(error)/len(error)
-
-    x = kwargs.get('x0', self.initialize(inputs[0], outputs[0]))
-    dt = kwargs.get('dt', 1e99)
-    stability_tol = kwargs.get('stability_tol', 0.95)
-
-    if not isinstance(x, self.StateContainer):
-        x = self.StateContainer(x)
-
-    if not isinstance(inputs[0], self.InputContainer):
-        inputs = [self.InputContainer(u_i) for u_i in inputs]
-    
-    if not isinstance(outputs[0], self.OutputContainer):
-        outputs = [self.OutputContainer(z_i) for z_i in outputs]
-
-    # Checks stability_tol is within bounds
-    # Throwing a default after the warning.
-    if stability_tol >= 1 or stability_tol < 0:
-        warn(f"configurable cutoff must be some float value in the domain (0, 1]. Received {stability_tol}. Resetting value to 0.95")
-        stability_tol = 0.95
-
-    counter = 0  # Needed to account for skipped (i.e., none) values
-    t_last = times[0]
-    err_total = 0
-    z_obs = self.output(x)
-    cutoffThreshold = math.floor(stability_tol * len(times))
-
-    for t, u, z in zip(times, inputs, outputs):
-        while t_last < t:
-            t_new = min(t_last + dt, t)
-            x = self.next_state(x, u, t_new-t_last)
-            t_last = t_new
-            if t >= t_last:
-                # Only recalculate if required
-                z_obs = self.output(x)
-        if not (None in z_obs.matrix or None in z.matrix):
-            # The none check above is used to cover the case where the model
-            # is not able to produce an output for a given input yet
-            # For example, in LSTM models, the first few inputs will not 
-            # produce an output until the model has received enough data
-            # This is true for any window-based model
-            if any(np.isnan(z_obs.matrix)):
-                if counter < cutoffThreshold:
-                    raise ValueError(f"Model unstable- NAN reached in simulation (t={t}) before cutoff threshold. Cutoff threshold is {cutoffThreshold}, or roughly {stability_tol * 100}% of the data")     
-                else:
-                    warn(f"Model unstable- NaN reached in simulation (t={t})")
-                    break
-            err_total += np.sum(np.square(z.matrix - z_obs.matrix), where=~np.isnan(z.matrix))
-            counter += 1
-
-    return err_total/counter
-
-def MAE(m, times, inputs, outputs, **kwargs):
-    """
-    Calculate the Mean Absolute Error between model behavior and some collected data.
-
-    Args:
-        m (PrognosticsModel): Model to use for comparison
-        times (list[float]): array of times for each sample
-        inputs (list[dict]): array of input dictionaries where input[x] corresponds to time[x]
-        outputs (list[dict]): array of output dictionaries where output[x] corresponds to time[x]
-
-    Keyword Args:
-        x0 (StateContainer): Current State of the model
-        dt (float, optional): Minimum time step in simulation. Defaults to 1e99.
-        stability_tol (double, optional): Configurable parameter.
-            Configurable cutoff value, between 0 and 1, that determines the fraction of the data points for which the model must be stable.
-            In some cases, a prognostics model will become unstable under certain conditions, after which point the model can no longer represent behavior. 
-            stability_tol represents the fraction of the provided argument `times` that are required to be met in simulation, 
-            before the model goes unstable in order to produce a valid estimate of mean squared error. 
-
-            If the model goes unstable before stability_tol is met, NaN is returned. 
-            Else, model goes unstable after stability_tol is met, the mean squared error calculated from data up to the instability is returned.
-
-    Returns:
-        float: MAE between model and data
-    """
-    if isinstance(times[0], Iterable):
-        # Calculate error for each
-        error = [MAE(t, i, z, **kwargs) for (t, i, z) in zip(times, inputs, outputs)]
-        return sum(error)/len(error)
-
-    x = kwargs.get('x0', m.initialize(inputs[0], outputs[0]))
-    dt = kwargs.get('dt', 1e99)
-    stability_tol = kwargs.get('stability_tol', 0.95)
-
-    if not isinstance(x, m.StateContainer):
-        x = m.StateContainer(x)
-
-    if not isinstance(inputs[0], m.InputContainer):
-        inputs = [m.InputContainer(u_i) for u_i in inputs]
-    
-    if not isinstance(outputs[0], m.OutputContainer):
-        outputs = [m.OutputContainer(z_i) for z_i in outputs]
-
-    # Checks stability_tol is within bounds
-    # Throwing a default after the warning.
-    if stability_tol >= 1 or stability_tol < 0:
-        warn(f"configurable cutoff must be some float value in the domain (0, 1]. Received {stability_tol}. Resetting value to 0.95")
-        stability_tol = 0.95
-
-    counter = 0  # Needed to account for skipped (i.e., none) values
-    t_last = times[0]
-    err_total = 0
-    z_obs = m.output(x)  # Initialize
-    cutoffThreshold = math.floor(stability_tol * len(times))
-
-    for t, u, z in zip(times, inputs, outputs):
-        while t_last < t:
-            t_new = min(t_last + dt, t)
-            x = m.next_state(x, u, t_new-t_last)
-            t_last = t_new
-            if t >= t_last:
-                # Only recalculate if required
-                z_obs = m.output(x)
-        if not (None in z_obs.matrix or None in z.matrix):
-            # The none check above is used to cover the case where the model
-            # is not able to produce an output for a given input yet
-            # For example, in LSTM models, the first few inputs will not 
-            # produce an output until the model has received enough data
-            # This is true for any window-based model
-            if any(np.isnan(z_obs.matrix)):
-                if counter < cutoffThreshold:
-                    raise ValueError(f"Model unstable- NAN reached in simulation (t={t}) before cutoff threshold. Cutoff threshold is {cutoffThreshold}, or roughly {stability_tol * 100}% of the data")     
-                else:
-                    warn(f"Model unstable- NaN reached in simulation (t={t})")
-                    break
-            err_total += np.sum(
-                np.abs(z.matrix - z_obs.matrix))
-            counter += 1
-    return err_total/counter
-
-def MAPE(m, times, inputs, outputs, **kwargs):
-    """
-    Calculate the Mean Absolute Percentage Error between model behavior and some collected data.
-
-    Args:
-        m (PrognosticsModel): Model to use for comparison
-        times (list[float]): array of times for each sample
-        inputs (list[dict]): array of input dictionaries where input[x] corresponds to time[x]
-        outputs (list[dict]): array of output dictionaries where output[x] corresponds to time[x]
-
-    Keyword Args:
-        x0 (StateContainer): Current State of the model
-        dt (float, optional): Minimum time step in simulation. Defaults to 1e99.
-        stability_tol (double, optional): Configurable parameter.
-            Configurable cutoff value, between 0 and 1, that determines the fraction of the data points for which the model must be stable.
-            In some cases, a prognostics model will become unstable under certain conditions, after which point the model can no longer represent behavior. 
-            stability_tol represents the fraction of the provided argument `times` that are required to be met in simulation, 
-            before the model goes unstable in order to produce a valid estimate of mean squared error. 
-
-            If the model goes unstable before stability_tol is met, NaN is returned. 
-            Else, model goes unstable after stability_tol is met, the mean squared error calculated from data up to the instability is returned.
-
-    Returns:
-        float: MAPE between model and data
-    """
-    if isinstance(times[0], Iterable):
-        # Calculate error for each
-        error = [MAPE(t, i, z, **kwargs) for (t, i, z) in zip(times, inputs, outputs)]
-        return sum(error)/len(error)
-
-    x = kwargs.get('x0', m.initialize(inputs[0], outputs[0]))
-    dt = kwargs.get('dt', 1e99)
-    stability_tol = kwargs.get('stability_tol', 0.95)
-
-    if not isinstance(x, m.StateContainer):
-        x = m.StateContainer(x)
-
-    if not isinstance(inputs[0], m.InputContainer):
-        inputs = [m.InputContainer(u_i) for u_i in inputs]
-    
-    if not isinstance(outputs[0], m.OutputContainer):
-        outputs = [m.OutputContainer(z_i) for z_i in outputs]
-
-    # Checks stability_tol is within bounds
-    # Throwing a default after the warning.
-    if stability_tol >= 1 or stability_tol < 0:
-        warn(f"configurable cutoff must be some float value in the domain (0, 1]. Received {stability_tol}. Resetting value to 0.95")
-        stability_tol = 0.95
-
-    counter = 0  # Needed to account for skipped (i.e., none) values
-    t_last = times[0]
-    err_total = 0
-    z_obs = m.output(x)  # Initialize
-    cutoffThreshold = math.floor(stability_tol * len(times))
-
-    for t, u, z in zip(times, inputs, outputs):
-        while t_last < t:
-            t_new = min(t_last + dt, t)
-            x = m.next_state(x, u, t_new-t_last)
-            t_last = t_new
-            if t >= t_last:
-                # Only recalculate if required
-                z_obs = m.output(x)
-        if not (None in z_obs.matrix or None in z.matrix):
-            # The none check above is used to cover the case where the model
-            # is not able to produce an output for a given input yet
-            # For example, in LSTM models, the first few inputs will not 
-            # produce an output until the model has received enough data
-            # This is true for any window-based model
-            if any(np.isnan(z_obs.matrix)):
-                if counter < cutoffThreshold:
-                    raise ValueError(f"Model unstable- NAN reached in simulation (t={t}) before cutoff threshold. Cutoff threshold is {cutoffThreshold}, or roughly {stability_tol * 100}% of the data")     
-                else:
-                    warn(f"Model unstable- NaN reached in simulation (t={t})")
-                    break
-            err_total += np.sum(np.abs(z.matrix - z_obs.matrix)/z.matrix)
-            counter += 1
-    return err_total/counter
diff --git a/src/prog_models/utils/containers.py b/src/prog_models/utils/containers.py
index 6b8012043..75d78470d 100644
--- a/src/prog_models/utils/containers.py
+++ b/src/prog_models/utils/containers.py
@@ -1,39 +1,47 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the
 # National Aeronautics and Space Administration.  All Rights Reserved.
 
-import numpy as np
-from typing import Union
-
+from numpy import float64, matrix, ndarray, array, newaxis, nan, delete, atleast_1d, array_equal
+import pandas as pd
 from prog_models.exceptions import ProgModelTypeError
+from typing import Union
 
 
 class DictLikeMatrixWrapper():
     """
-    A container that behaves like a dictionary, but is backed by a numpy array, which is itself directly accessable. This is used for model states, inputs, and outputs- and enables efficient matrix operations.
-    
+    A container that uses pandas dictionary like data structure, but is backed by a numpy array, which is itself directly accessible. This is used for model states, inputs, and outputs- and enables efficient matrix operations.
+
     Arguments:
         keys -- list: The keys of the dictionary. e.g., model.states or model.inputs
         data -- dict or numpy array: The contained data (e.g., :term:`input`, :term:`state`, :term:`output`). If numpy array should be column vector in same order as keys
     """
 
-    def __init__(self, keys: list, data: Union[dict, np.array]):
-        """ Initializes the container
+    def __init__(self, keys: list, data: Union[dict, array]):
+        """
+        Initializes the container
         """
         if not isinstance(keys, list):
-            keys = list(keys)   # creates list with keys
+            keys = list(keys)  # creates list with keys
         self._keys = keys.copy()
-        if isinstance(data, np.matrix):
-            self.matrix = np.array(data, dtype=np.float64)
-        elif isinstance(data, np.ndarray):
+
+        if isinstance(data, matrix):
+            self.data = pd.DataFrame(array(data, dtype=float64), self._keys, dtype=float64)
+            self.matrix = self.data.to_numpy(dtype=float64)
+        elif isinstance(data, ndarray):
             if data.ndim == 1:
-                data = data[np.newaxis].T
+                data = data[newaxis].T
+                self.data = pd.DataFrame(data, self._keys)
+            self.data = pd.DataFrame(data, self._keys).T
             self.matrix = data
         elif isinstance(data, (dict, DictLikeMatrixWrapper)):
-          # ravel is used to prevent vectorized case, where data[key] returns multiple values,  from resulting in a 3D matrix
-          self.matrix = np.array(
-             [
-                  np.ravel([data[key]]) if key in data else [None] for key in keys
-              ], dtype=np.float64)
+            if data and not isinstance(list(data.values())[0], ndarray):  # len(self.matrix[0]) == 1:
+                if isinstance(data, DictLikeMatrixWrapper):
+                    data = dict(data.copy())
+                self.data = pd.DataFrame(data, columns=self._keys, index=[0], dtype=float64).replace(
+                    nan, None)
+            else:
+                self.data = pd.DataFrame(data, columns=self._keys)
+            self.matrix = self.data.to_numpy(dtype=float64).T if len(data) > 0 else array([])
         else:
             raise ProgModelTypeError(f"Data must be a dictionary or numpy array, not {type(data)}")
 
@@ -41,67 +49,78 @@ def __reduce__(self):
         """
         reduce is overridden for pickles
         """
-        return (DictLikeMatrixWrapper, (self._keys, self.matrix))
+        return DictLikeMatrixWrapper, (self._keys, self.matrix)
 
     def __getitem__(self, key: str) -> int:
         """
         get all values associated with a key, ex: all values of 'i'
         """
-        row = self.matrix[self._keys.index(key)]    #creates list from a row of matrix
-        if len(row) == 1:   #list contains 1 value, returns that value (non-vectorized)
+        row = self.data.loc[:, key].to_list()  # creates list from a column of pandas DF
+        if len(row) == 1:  # list contains 1 value, returns that value (non-vectorized)
             return row[0]
-        return row #returns entire row/list (vectorized case)
+        else:
+            return row  # returns entire row/list (vectorized case)
 
     def __setitem__(self, key: str, value: int) -> None:
         """
         sets a row at the key given
         """
-        index = self._keys.index(key)   #the int value index for the key given
-        self.matrix[index] = np.atleast_1d(value)
+        index = self._keys.index(key)  # the int value index for the key given
+        self.matrix[index] = atleast_1d(value)
 
-    
     def __delitem__(self, key: str) -> None:
         """
         removes row associated with key
         """
-        self.matrix = np.delete(self.matrix, self._keys.index(key), axis=0)
+        # self.matrix = delete(self.matrix, self._keys.index(key), axis=0)
         self._keys.remove(key)
+        self.data = self.data.drop(columns=[key], axis=1)
+        self.matrix = self.data.T.to_numpy()
 
     def __add__(self, other: "DictLikeMatrixWrapper") -> "DictLikeMatrixWrapper":
         """
         add another matrix to the existing matrix
         """
-        return DictLikeMatrixWrapper(self._keys, self.matrix + other.matrix)
+        rowadded = self.data.add(other.data).T.to_numpy()
+        return DictLikeMatrixWrapper(self._keys, rowadded)
 
     def __iter__(self):
         """
         creates iterator object for the list of keys
         """
-        return iter(self._keys)
+        return iter(self.data.keys())
 
     def __len__(self) -> int:
         """
         returns the length of key list
         """
-        return len(self._keys)
+        return len(self.data.keys())
 
     def __eq__(self, other: "DictLikeMatrixWrapper") -> bool:
         """
         Compares two DictLikeMatrixWrappers (i.e. *Containers) or a DictLikeMatrixWrapper and a dictionary
         """
-        if isinstance(other, dict): #checks that the list of keys for each matrix match
-            list_key_check = (list(self.keys()) == list(other.keys()))  #checks that the list of keys for each matrix are equal
-            matrix_check = (self.matrix == np.array([[other[key]] for key in self._keys])).all()    #checks to see that each row matches
-            return list_key_check and matrix_check
+        if isinstance(other, dict):  # checks that the list of keys for each matrix match
+            list_key_check = (list(self.keys()) == list(
+                other.keys()))  # checks that the list of keys for each matrix are equal
+            matrix_check = (self.matrix == array(
+                [[other[key]] for key in self._keys])).all()  # checks to see that each row matches
+            # check if DF is the same or if both are empty
+            df_check = self.data.equals(other.data) or (self.data.empty and other.data.empty)
+            return list_key_check and matrix_check and df_check
         list_key_check = self.keys() == other.keys()
         matrix_check = (self.matrix == other.matrix).all()
-        return list_key_check and matrix_check
+        # check if DF is the same or if both are empty
+        df_check = self.data.equals(other.data) or (self.data.empty and other.data.empty)
+        return list_key_check and matrix_check and df_check
 
     def __hash__(self):
         """
         returns hash value sum for keys and matrix
         """
-        return hash(self.keys) + hash(self.matrix)
+        sum_hash = 0
+        sum_hash = (sum_hash + x for x in pd.util.hash_pandas_object(self.data))
+        return sum_hash
 
     def __str__(self) -> str:
         """
@@ -109,54 +128,60 @@ def __str__(self) -> str:
         """
         return self.__repr__()
 
-    def get(self, key, default=None):
+    def get(self, key: str, default=None):
         """
         gets the list of values associated with the key given
         """
         if key in self._keys:
-            return self[key]
+            return self.data.loc[0, key]
         return default
 
     def copy(self) -> "DictLikeMatrixWrapper":
         """
         creates copy of object
         """
-        return DictLikeMatrixWrapper(self._keys, self.matrix.copy())
+        matrix_df = self.data.T.to_numpy().copy()
+        return DictLikeMatrixWrapper(self._keys, matrix_df)
 
     def keys(self) -> list:
         """
         returns list of keys for container
         """
-        return self._keys
+        return self.data.keys().to_list()
 
-    def values(self) -> np.array:
+    def values(self) -> array:
         """
         returns array of matrix values
         """
-        if len(self.matrix) > 0 and len(self.matrix[0]) == 1:   #if the first row of the matrix has one value (i.e., non-vectorized)
-            return np.array([value[0] for value in self.matrix])    #the value from the first row
-        return self.matrix    # the matrix (vectorized case)
+        matrix_df = self.data.T.to_numpy()
+        if len(matrix_df) > 0 and len(
+                matrix_df[0]) == 1:  # if the first row of the matrix has one value (i.e., non-vectorized)
+            return array([value[0] for value in matrix_df])  # the value from the first row
+        return matrix_df  # the matrix (vectorized case)
 
     def items(self) -> zip:
         """
         returns keys and values as a list of tuples (for iterating)
         """
-        if len(self.matrix) > 0 and len(self.matrix[0]) == 1:   #first row of the matrix has one value (non-vectorized case)
-            return zip(self._keys, np.array([value[0] for value in self.matrix]))
-        return zip(self._keys, self.matrix)
+        matrix_df = self.data.T.to_numpy()
+        if len(matrix_df) > 0 and len(matrix_df[0]) == 1:  # first row of the matrix has one value (non-vectorized case)
+            return zip(self.data.keys(), array([value[0] for value in matrix_df]))
+        return zip(self.data.keys(), matrix_df)
 
     def update(self, other: "DictLikeMatrixWrapper") -> None:
         """
         merges other DictLikeMatrixWrapper, updating values
         """
-        for key in other.keys():   
-            if key in self._keys:   #checks to see if every key in 'other' is in 'self'
+        for key in other.data.index.to_list():
+            if key in self.data.index.to_list():  # checks to see if the key exists
                 # Existing key
-                self[key] = other[key]
-            else:   #else it isn't it is appended to self._keys list
-                # A new key!
-                self._keys.append(key)
-                self.matrix = np.vstack((self.matrix, np.array([other[key]])))
+                self.data.loc[key] = other.data.loc[key]
+            else:  # the key doesn't exist within
+                # the key
+                temp_df = DictLikeMatrixWrapper([key], {key: other.data.loc[key, 0]})
+                self.data = pd.concat([self.data, temp_df.data])
+        self._keys = self.data.index.to_list()
+        self.matrix = self.data.to_numpy()
 
     def __contains__(self, key: str) -> bool:
         """
@@ -166,9 +191,11 @@ def __contains__(self, key: str) -> bool:
         -------
         >>> from prog_models.utils.containers import DictLikeMatrixWrapper
         >>> dlmw = DictLikeMatrixWrapper(['a', 'b', 'c'], {'a': 1, 'b': 2, 'c': 3})
-        >>> 'a' in dlmw  # True
+        >>> 'a' in dlmw
+        True
         """
-        return key in self._keys
+        key_list = self.data.keys()
+        return key in key_list
 
     def __repr__(self) -> str:
         """
@@ -176,6 +203,6 @@ def __repr__(self) -> str:
 
         returns: a string of dictionaries containing all the keys and associated matrix values
         """
-        if len(self.matrix) > 0 and len(self.matrix[0]) == 1:   #the matrix has rows and the first row/list has one value in it
-            return str({key: value[0] for key, value in zip(self._keys, self.matrix)})
-        return str(dict(zip(self._keys, self.matrix)))  
+        if len(self.data.columns) > 0:
+            return str(self.data.to_dict('records')[0])
+        return str(self.data.to_dict())
diff --git a/src/prog_models/utils/next_state.py b/src/prog_models/utils/next_state.py
deleted file mode 100644
index baba77801..000000000
--- a/src/prog_models/utils/next_state.py
+++ /dev/null
@@ -1,146 +0,0 @@
-# Copyright © 2021 United States Government as represented by the Administrator of the
-# National Aeronautics and Space Administration.  All Rights Reserved.
-
-def euler_next_state(model, x, u, dt: float):
-    """
-    State transition equation using simple euler integration: Calls next_state(), calculating the next state, and then adds noise and applies limits
-
-    Parameters
-    ----------
-    x : StateContainer
-        state, with keys defined by model.states \n
-        e.g., x = m.StateContainer({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
-    u : InputContainer
-        Inputs, with keys defined by model.inputs \n
-        e.g., u = m.InputContainer({'i':3.2}) given inputs = ['i']
-    dt : float
-        Timestep size in seconds (≥ 0) \n
-        e.g., dt = 0.1
-
-    Returns
-    -------
-    x : StateContainer
-        Next state, with keys defined by model.states
-        e.g., x = m.StateContainer({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
-
-    See Also
-    --------
-    PrognosticsModel.next_state
-    """
-    # Calculate next state and add process noise
-    next_state = model.apply_process_noise(model.next_state(x, u, dt), dt)
-
-    # Apply Limits
-    return model.apply_limits(next_state)
-
-def euler_next_state_wrapper(model, x, u, dt: float):
-    """
-    State transition equation using simple euler integration: Calls next_state(), calculating the next state, and then adds noise and applies limits
-
-    Parameters
-    ----------
-    x : StateContainer
-        state, with keys defined by model.states \n
-        e.g., x = m.StateContainer({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
-    u : InputContainer
-        Inputs, with keys defined by model.inputs \n
-        e.g., u = m.InputContainer({'i':3.2}) given inputs = ['i']
-    dt : float
-        Timestep size in seconds (≥ 0) \n
-        e.g., dt = 0.1
-
-    Returns
-    -------
-    x : StateContainer
-        Next state, with keys defined by model.states
-        e.g., x = m.StateContainer({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
-
-    See Also
-    --------
-    PrognosticsModel.next_state
-    """
-    # Calculate next state and add process noise
-    next_state = model.StateContainer(model.next_state(x, u, dt))
-    next_state = model.apply_process_noise(next_state, dt)
-
-    # Apply Limits
-    return model.apply_limits(next_state)
-
-def rk4_next_state(model, x, u, dt: float):
-    """
-    State transition equation using rungekutta4 integration: Calls next_state(), calculating the next state, and then adds noise and applies limits
-
-    Parameters
-    ----------
-    x : StateContainer
-        state, with keys defined by model.states \n
-        e.g., x = m.StateContainer({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
-    u : InputContainer
-        Inputs, with keys defined by model.inputs \n
-        e.g., u = m.InputContainer({'i':3.2}) given inputs = ['i']
-    dt : float
-        Timestep size in seconds (≥ 0) \n
-        e.g., dt = 0.1
-
-    Returns
-    -------
-    x : StateContainer
-        Next state, with keys defined by model.states
-        e.g., x = m.StateContainer({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
-
-    See Also
-    --------
-    PrognosticsModel.next_state
-    """
-    dx1 = model.StateContainer(model.dx(x, u))
-    x2 = x.matrix + dx1.matrix*dt/2
-    dx2 = model.dx(x2, u)
-
-    x3 = model.StateContainer({key: x[key] + dt*dx_i/2 for key, dx_i in dx2.items()})
-    dx3 = model.dx(x3, u)
-    
-    x4 = model.StateContainer({key: x[key] + dt*dx_i for key, dx_i in dx3.items()})
-    dx4 = model.dx(x4, u)
-
-    x = model.StateContainer({key: x[key] + dt/3*(dx1[key]/2 + dx2[key] + dx3[key] + dx4[key]/2) for key in dx1.keys()})
-    return model.apply_limits(model.apply_process_noise(x))
-
-def rk4_next_state_wrapper(model, x, u, dt: float):
-    """
-    State transition equation using rungekutta4 integration: Calls next_state(), calculating the next state, and then adds noise and applies limits
-
-    Parameters
-    ----------
-    x : StateContainer
-        state, with keys defined by model.states \n
-        e.g., x = m.StateContainer({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
-    u : InputContainer
-        Inputs, with keys defined by model.inputs \n
-        e.g., u = m.InputContainer({'i':3.2}) given inputs = ['i']
-    dt : float
-        Timestep size in seconds (≥ 0) \n
-        e.g., dt = 0.1
-
-    Returns
-    -------
-    x : StateContainer
-        Next state, with keys defined by model.states
-        e.g., x = m.StateContainer({'abc': 332.1, 'def': 221.003}) given states = ['abc', 'def']
-
-    See Also
-    --------
-    PrognosticsModel.next_state
-    """
-    dx1 = model.StateContainer(model.dx(x, u))
-                
-    x2 = model.StateContainer({key: x[key] + dt*dx_i/2 for key, dx_i in dx1.items()})
-    dx2 = model.dx(x2, u)
-
-    x3 = model.StateContainer({key: x[key] + dt*dx_i/2 for key, dx_i in dx2.items()})
-    dx3 = model.dx(x3, u)
-    
-    x4 = model.StateContainer({key: x[key] + dt*dx_i for key, dx_i in dx3.items()})
-    dx4 = model.dx(x4, u)
-
-    x = model.StateContainer({key: x[key] + dt/3*(dx1[key]/2 + dx2[key] + dx3[key] + dx4[key]/2) for key in dx1.keys()})
-    return model.apply_limits(model.apply_process_noise(x))
diff --git a/src/prog_models/utils/parameters.py b/src/prog_models/utils/parameters.py
index 262543398..2a703d9ac 100644
--- a/src/prog_models/utils/parameters.py
+++ b/src/prog_models/utils/parameters.py
@@ -5,7 +5,6 @@
 from copy import deepcopy
 import json
 from numbers import Number
-import numpy as np
 import types
 from typing import Callable
 
@@ -29,7 +28,7 @@ class PrognosticsModelParameters(UserDict):
         dict_in: Initial parameters
         callbacks: Any callbacks for derived parameters f(parameters) : updates (dict)
     """
-    def __init__(self, model: "PrognosticsModel", dict_in: dict = {}, callbacks: dict = {}, _copy: bool = True):
+    def __init__(self, model : "PrognosticsModel", dict_in : dict = {}, callbacks : dict = {}, _copy: bool = True):
         super().__init__()
         self._m = model
         self.callbacks = {}
@@ -49,29 +48,16 @@ def __init__(self, model: "PrognosticsModel", dict_in: dict = {}, callbacks: dic
     def __sizeof__(self):
         return getsizeof(self)
 
-    def __eq__(self, other):
-        if set(self.data.keys()) != set(other.data.keys()):
-            return False
-        for key, value in self.data.items():
-            if not np.all(value == other[key]): 
-                  # Note: np.all is used to handle numpy array elements
-                  # Otherwise value == other[key] would return a numpy array of bools for each element
-                  return False
-        return True
-        
-
     def copy(self):
         return self.__class__(self._m, self.data, self.callbacks, _copy=False)
 
     def __copy__(self):
         return self.__class__(self._m, self.data, self.callbacks, _copy=False)
     
-    def __deepcopy__(self, memo):
-        result = self.__class__(self._m, self.data, self.callbacks, _copy=True)
-        memo[id(self)] = result
-        return result
+    def __deepcopy__(self):
+        return self.__class__(self._m, self.data, self.callbacks, _copy=True)
 
-    def __setitem__(self, key: str, value: float, _copy: bool = False) -> None:
+    def __setitem__(self, key : str, value : float, _copy : bool = False) -> None:
         """Set model configuration, overrides dict.__setitem__()
 
         Args:
diff --git a/src/prog_models/utils/serialization.py b/src/prog_models/utils/serialization.py
index 19289ba46..5ba6e63dd 100644
--- a/src/prog_models/utils/serialization.py
+++ b/src/prog_models/utils/serialization.py
@@ -4,7 +4,7 @@
 import json
 import numpy as np
 
-from prog_models.utils.containers import DictLikeMatrixWrapper
+from .containers import DictLikeMatrixWrapper
 
 __all__ = ['CustomEncoder', 'custom_decoder']
 
diff --git a/tests/__init__.py b/tests/__init__.py
index f50d9acbc..b2b2929e1 100644
--- a/tests/__init__.py
+++ b/tests/__init__.py
@@ -1,4 +1,4 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the
 # National Aeronautics and Space Administration.  All Rights Reserved.
 
-__all__ = ['test_base_models', 'test_battery', 'test_examples', 'test_pneumatic_valve', 'test_dict_like_matrix_wrapper', 'test_tutorials', 'test_datasets', 'test_manual', 'test_linear_model', 'test_composite']
+__all__ = ['test_base_models', 'test_battery', 'test_examples', 'test_pneumatic_valve', 'test_dict_like_matrix_wrapper', 'test_tutorials', 'test_datasets', 'test_manual']
diff --git a/tests/__main__.py b/tests/__main__.py
index 66db4e2f7..d45e7543f 100644
--- a/tests/__main__.py
+++ b/tests/__main__.py
@@ -1,27 +1,51 @@
 # Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
 
-from tests.test_base_models import main as base_models_main
-from tests.test_sim_result import main as sim_result_main
-from tests.test_dict_like_matrix_wrapper import main as dict_like_matrix_wrapper_main
-from tests.test_examples import main as examples_main
-from tests.test_centrifugal_pump import main as centrifugal_pump_main
-from tests.test_pneumatic_valve import main as pneumatic_valve_main
-from tests.test_battery import main as battery_main
-from tests.test_tutorials import main as tutorials_main
-from tests.test_datasets import main as datasets_main
-from tests.test_powertrain import main as powertrain_main
-from tests.test_surrogates import main as surrogates_main
-from tests.test_data_model import main as lstm_main
-from tests.test_direct import main as direct_main
-from tests.test_linear_model import main as linear_main
-from tests.test_composite import main as composite_main
-from tests.test_serialization import main as serialization_main
-from tests.test_ensemble import main as ensemble_main
+from .test_base_models import main as base_models_main
+from .test_sim_result import main as sim_result_main
+from .test_dict_like_matrix_wrapper import main as dict_like_matrix_wrapper_main
+from .test_examples import main as examples_main
+from .test_centrifugal_pump import main as centrifugal_pump_main
+from .test_pneumatic_valve import main as pneumatic_valve_main
+from .test_battery import main as battery_main
+from .test_tutorials import main as tutorials_main
+from .test_datasets import main as datasets_main
+from .test_powertrain import main as powertrain_main
+from .test_surrogates import main as surrogates_main
+from .test_data_model import main as lstm_main
+from .test_direct import main as direct_main
+from .test_linear_model import main as linear_main
+
+from io import StringIO
+import matplotlib.pyplot as plt
+import sys
+from timeit import timeit
+from unittest.mock import patch
+
+from examples import sim as sim_example
+
+def _test_ex():
+    # Run example
+    sim_example.run_example()
 
 if __name__ == '__main__':
     was_successful = True
 
-    print("\n\nTesting individual execution of test files")
+    try:
+        # set stdout (so it wont print)
+        sys.stdout = StringIO()
+        
+        with patch('matplotlib.pyplot.show'):
+            runtime = timeit(_test_ex, number=10)
+            plt.close('all')
+
+        # Reset stdout 
+        sys.stdout = sys.__stdout__
+        print(f"\nExample Runtime: {runtime}")
+    except Exception as e:
+        print("\Benchmarking Failed: ", e)
+        was_successful = False
+
+    print("\n\nTesting individual exectution of test files")
 
     # Run tests individually to test them and make sure they can be executed individually
     try:
@@ -40,7 +64,6 @@
         was_successful = False
 
     try:
-    
         examples_main()
     except Exception:
         was_successful = False
@@ -94,21 +117,6 @@
         linear_main()
     except Exception:
         was_successful = False
-    
-    try:
-        composite_main()
-    except Exception:
-        was_successful = False
-
-    try:
-        serialization_main()
-    except Exception:
-        was_successful = False
-
-    try:
-        ensemble_main()
-    except Exception:
-        was_successful = False
 
     if not was_successful:
         raise Exception("Failed test")
diff --git a/tests/benchmarking.py b/tests/benchmarking.py
index c408cc3d4..f13a95f80 100644
--- a/tests/benchmarking.py
+++ b/tests/benchmarking.py
@@ -20,7 +20,7 @@
     print(f'{t} |')
 
     print(FORMAT_STR.format('model initialization'), end='')
-    t = timeit.timeit('ThrownObject()', 'from prog_models.models import ThrownObject', number = 1000, timer = process_time)
+    t = timeit.timeit('ThrownObject()', 'from prog_models.models import ThrownObject', timer = process_time)
     print(f'{t} |')
 
     m = ThrownObject()
diff --git a/tests/test_base_models.py b/tests/test_base_models.py
index 90c339cde..2b25f0951 100644
--- a/tests/test_base_models.py
+++ b/tests/test_base_models.py
@@ -7,14 +7,14 @@
 import pickle
 import sys
 import unittest
+import pandas as pd
 
 # This ensures that the directory containing ProgModelTemplate is in the python search directory
 sys.path.append(join(dirname(__file__), ".."))
 
-from prog_models import ProgModelTypeError, ProgModelInputException, ProgModelException, PrognosticsModel, CompositeModel
-from prog_models.models import ThrownObject, BatteryElectroChemEOD
+from prog_models import *
+from prog_models.models import *
 from prog_models.models.test_models.linear_models import (OneInputNoOutputNoEventLM, OneInputOneOutputNoEventLM, OneInputNoOutputOneEventLM, OneInputOneOutputNoEventLMPM)
-from prog_models.models.test_models.linear_thrown_object import (LinearThrownObject, LinearThrownDiffThrowingSpeed, LinearThrownObjectUpdatedInitalizedMethod, LinearThrownObjectDiffDefaultParams)
 
 
 class MockModel():
@@ -41,7 +41,7 @@ def output(self, x):
         return self.OutputContainer({'o1': x['a'] + x['b'] + x['c']})
 
 
-class MockProgModel(MockModel, PrognosticsModel):
+class MockProgModel(MockModel, prognostics_model.PrognosticsModel):
     events = ['e1', 'e2']
 
     def event_state(self, x):
@@ -69,6 +69,42 @@ def derived_callback3(config):
         'p4': -2 * config['p2'], 
     }
 
+class LinearThrownObject(LinearModel):
+    inputs = [] 
+    states = ['x', 'v']
+    outputs = ['x']
+    events = ['impact']
+
+    A = np.array([[0, 1], [0, 0]])
+    E = np.array([[0], [-9.81]])
+    C = np.array([[1, 0]])
+    F = None # Will override method
+
+    default_parameters = {
+        'thrower_height': 1.83,  # m
+        'throwing_speed': 40,  # m/s
+        'g': -9.81  # Acceleration due to gravity in m/s^2
+    }
+
+    def initialize(self, u=None, z=None):
+        return self.StateContainer({
+            'x': self.parameters['thrower_height'],  # Thrown, so initial altitude is height of thrower
+            'v': self.parameters['throwing_speed']  # Velocity at which the ball is thrown - this guy is a professional baseball pitcher
+            })
+    
+    def threshold_met(self, x):
+        return {
+            'falling': x['v'] < 0,
+            'impact': x['x'] <= 0
+        }
+
+    def event_state(self, x): 
+        x_max = x['x'] + np.square(x['v'])/(-self.parameters['g']*2) # Use speed and position to estimate maximum height
+        return {
+            'falling': np.maximum(x['v']/self.parameters['throwing_speed'],0),  # Throwing speed is max speed
+            'impact': np.maximum(x['x']/x_max,0) if x['v'] < 0 else 1  # 1 until falling begins, then it's fraction of height
+        }
+
 class MockModelWithDerived(MockProgModel):
     param_callbacks = {
             'p1': [derived_callback],
@@ -77,34 +113,35 @@ class MockModelWithDerived(MockProgModel):
 
 
 class TestModels(unittest.TestCase):
+    """
     def setUp(self):
         # set stdout (so it wont print)
         sys.stdout = io.StringIO()
 
     def tearDown(self):
-        sys.stdout = sys.__stdout__
+        sys.stdout = sys.__stdout__"""
 
     def test_non_container(self):
         class MockProgModelStateDict(MockProgModel):
             def next_state(self, x, u, dt):
-                return {
+                ns_ret = {
                     'a': x['a'] + u['i1']*dt,
                     'b': x['b'],
                     'c': x['c'] - u['i2'],
                     't': x['t'] + dt
                 }
+                return ns_ret
         
-        m = MockProgModelStateDict(
-            process_noise_dist='none',
-            measurement_noise_dist='none')
-        
+        m = MockProgModelStateDict(process_noise_dist='none', measurement_noise_dist='none')
         def load(t, x=None):
             return {'i1': 1, 'i2': 2.1}
 
         # Any event, default
         (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(load, {'o1': 0.8}, **{'dt': 0.5, 'save_freq': 1.0})
         self.assertAlmostEqual(times[-1], 5.0, 5)
+        self.assertAlmostEqual(outputs.frame['time'].iloc[-1], 5.0, 5)
         self.assertAlmostEqual(outputs[-1]['o1'], -13.2)
+        self.assertAlmostEqual(outputs.get_frame_data()['o1'].iloc[-1], -13.2)
         self.assertIsInstance(outputs[-1], m.OutputContainer)
 
         class MockProgModelStateNdarray(MockProgModel):
@@ -116,17 +153,16 @@ def next_state(self, x, u, dt):
                     [x['t'] + dt]]
                 )
 
-        m = MockProgModelStateNdarray(
-            process_noise_dist='none', 
-            measurement_noise_dist='none')
+        m = MockProgModelStateNdarray(process_noise_dist='none', measurement_noise_dist='none')
 
         # Any event, default
-        config = {'dt': 0.5, 'save_freq': 1.0}
-        (times, _, _, outputs, _) = m.simulate_to_threshold(
-            load,
-            {'o1': 0.8},
-            **config)
+        (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(load, {'o1': 0.8}, **{'dt': 0.5, 'save_freq': 1.0})
+
         self.assertAlmostEqual(times[-1], 5.0, 5)
+        self.assertAlmostEqual(inputs.frame['time'].iloc[-1], 5.0, 5)
+        self.assertAlmostEqual(states.frame['time'].iloc[-1], 5.0, 5)
+        self.assertAlmostEqual(outputs.frame['time'].iloc[-1], 5.0, 5)
+        self.assertAlmostEqual(event_states.frame['time'].iloc[-1], 5.0, 5)
 
     def test_size(self):
         m = MockProgModel()
@@ -220,10 +256,10 @@ def test_derived(self):
     def test_broken_models(self):
 
 
-        class missing_states(PrognosticsModel):
+        class missing_states(prognostics_model.PrognosticsModel):
             inputs = ['i1', 'i2']
             outputs = ['o1']
-            parameters = {'process_noise': 0.1}
+            parameters = {'process_noise':0.1}
             def initialize(self, u, z):
                 pass
             def next_state(self, x, u, dt):
@@ -232,11 +268,11 @@ def output(self, x):
                 pass
         
 
-        class empty_states(PrognosticsModel):
+        class empty_states(prognostics_model.PrognosticsModel):
             states = []
             inputs = ['i1', 'i2']
             outputs = ['o1']
-            parameters = {'process_noise': 0.1}
+            parameters = {'process_noise':0.1}
             def initialize(self, u, z):
                 pass
             def next_state(self, x, u, dt):
@@ -245,10 +281,10 @@ def output(self, x):
                 pass
         
 
-        class missing_inputs(PrognosticsModel):
+        class missing_inputs(prognostics_model.PrognosticsModel):
             states = ['x1', 'x2']
             outputs = ['o1']
-            parameters = {'process_noise': 0.1}
+            parameters = {'process_noise':0.1}
             def initialize(self, u, z):
                 pass
             def next_state(self, x, u, dt):
@@ -257,10 +293,10 @@ def output(self, x):
                 pass
         
 
-        class missing_outputs(PrognosticsModel):
+        class missing_outputs(prognostics_model.PrognosticsModel):
             states = ['x1', 'x2']
             inputs = ['i1']
-            parameters = {'process_noise': 0.1}
+            parameters = {'process_noise':0.1}
             def initialize(self, u, z):
                 pass
             def next_state(self, x, u, dt):
@@ -269,29 +305,32 @@ def output(self, x):
                 pass
         
 
-        class missing_initiialize(PrognosticsModel):
+        class missing_initiialize(prognostics_model.PrognosticsModel):
             inputs = ['i1']
             states = ['x1', 'x2']
             outputs = ['o1']
-            parameters = {'process_noise': 0.1}
+            parameters = {'process_noise':0.1}
             def next_state(self, x, u, dt):
                 pass
             def output(self, x):
                 pass
         
 
-        class missing_output(PrognosticsModel):
+        class missing_output(prognostics_model.PrognosticsModel):
             inputs = ['i1']
             states = ['x1', 'x2']
             outputs = ['o1']
-            parameters = {'process_noise': 0.1}
+            parameters = {'process_noise':0.1}
             def initialize(self, u, z):
                 pass
             def next_state(self, x, u, dt):
                 pass
 
-        with self.assertRaises(ProgModelTypeError):
+        try: 
             m = missing_states()
+            self.fail("Should not have worked, missing 'states'")
+        except ProgModelTypeError:
+            pass
 
         m = empty_states()
         self.assertEqual(len(m.states), 0)
@@ -330,9 +369,12 @@ def add_one(self, x):
         x = getattr(m, "apply_{}".format(noise_key))({key: 1 for key in keys})
         self.assertEqual(x[keys[0]], 2)
 
-        with self.assertRaises(Exception):
+        try:
             noise = []
-            m = MockProgModel(**{noise_key: noise})         
+            m = MockProgModel(**{noise_key: noise})
+            self.fail("Should have raised exception - improper format")
+        except Exception:
+            pass            
 
         # Test that it ignores process_noise_dist in case where process_noise is a function
         m = MockProgModel(**{noise_key: add_one, dist_key: 'invalid one'})
@@ -340,13 +382,19 @@ def add_one(self, x):
         self.assertEqual(x[keys[0]], 2)
 
         # Invalid dist
-        with self.assertRaises(ProgModelTypeError):
-            noise = {key: 0.0 for key in keys}
+        try:
+            noise = {key : 0.0 for key in keys}
             m = MockProgModel(**{noise_key: noise, dist_key: 'invalid one'})
+            self.fail("Invalid noise distribution")
+        except ProgModelTypeError:
+            pass
 
         # Invalid dist
-        with self.assertRaises(ProgModelTypeError):
+        try:
             m = MockProgModel(**{noise_key: 0, dist_key: 'invalid one'})
+            self.fail("Invalid noise distribution")
+        except ProgModelTypeError:
+            pass
 
         # Valid distributions
         m = MockProgModel(**{noise_key: 0, dist_key: 'uniform'})
@@ -415,7 +463,7 @@ def test_prog_model(self):
 
     def test_default_es_and_tm(self):
         # Test 1: TM only
-        class NoES(MockModel, PrognosticsModel):
+        class NoES(MockModel, prognostics_model.PrognosticsModel):
             events = ['e1', 'e2']
 
             def threshold_met(self, _):
@@ -427,7 +475,7 @@ def threshold_met(self, _):
         self.assertDictEqual(m.event_state({}), {'e1': 1.0, 'e2': 0.0})
 
         # Test 2: ES only
-        class NoTM(MockModel, PrognosticsModel):
+        class NoTM(MockModel, prognostics_model.PrognosticsModel):
             events = ['e1', 'e2']
 
             def event_state(self, _):
@@ -439,7 +487,7 @@ def event_state(self, _):
         self.assertDictEqual(m.event_state({}), {'e1': 0.0, 'e2': 1.0})
 
         # Test 3: Neither ES or TM 
-        class NoESTM(MockModel, PrognosticsModel):
+        class NoESTM(MockModel, prognostics_model.PrognosticsModel):
             events = []
 
         m = NoESTM()
@@ -513,8 +561,11 @@ def thresh_met(thresholds):
         (times, inputs, states, outputs, event_states) = m_noevents.simulate_to_threshold(linear_load, {'o1': 0.8}, **{'dt': 0.5, 'save_freq': 1.0, 'thresholds_met_eqn': thresh_met})
         self.assertListEqual(times, [0, 0.5]) # Only one step
 
-        with self.assertRaises(ProgModelInputException):
+        try:
             (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(load, {'o1': 0.8}, threshold_keys=['e1', 'e2', 'e3'], **{'dt': 0.5, 'save_freq': 1.0})
+            self.fail("Should fail- extra threshold key")
+        except ProgModelInputException:
+            pass
 
     def test_sim_past_thresh(self):
         m = MockProgModel(process_noise = 0.0)
@@ -530,11 +581,11 @@ def load(t, x=None):
             return {'i1': 1, 'i2': 2.1}
         (times, inputs, states, outputs, event_states) = m.simulate_to(6, load, {'o1': 0.8}, **{'dt': 0.5, 'save_freq': 1.0})
         named_results = m.simulate_to(6, load, {'o1': 0.8}, **{'dt': 0.5, 'save_freq': 1.0})
-        self.assertEqual(times, named_results.times)
-        self.assertEqual(inputs, named_results.inputs)
-        self.assertEqual(states, named_results.states)
-        self.assertEqual(outputs, named_results.outputs)
-        self.assertEqual(event_states, named_results.event_states)
+        self.assertEquals(times, named_results.times)
+        self.assertEquals(inputs, named_results.inputs)
+        self.assertEquals(states, named_results.states)
+        self.assertEquals(outputs, named_results.outputs)
+        self.assertEquals(event_states, named_results.event_states)
         
     def test_next_time_fcn(self):
         m = MockProgModel(process_noise = 0.0)
@@ -714,17 +765,29 @@ def load(t, x=None):
         (times, inputs, states, outputs, event_states) = m.simulate_to(0, load, {'o1': 0.8})
         self.assertEqual(len(times), 1)
 
-        with self.assertRaises(ProgModelInputException):
+        try:
             m.simulate_to(-30, load, {'o1': 0.8})
+            self.fail("Should have failed- time must be greater than 0")
+        except ProgModelInputException:
+            pass
 
-        with self.assertRaises(ProgModelInputException):
+        try:
             m.simulate_to([12], load, {'o1': 0.8})
+            self.fail("Should have failed- time must be a number")
+        except ProgModelInputException:
+            pass
 
-        with self.assertRaises(ProgModelInputException):
+        try:
             m.simulate_to(12, load, {'o2': 0.9})
+            self.fail("Should have failed- output must contain each field (e.g., o1)")
+        except ProgModelInputException:
+            pass
 
-        with self.assertRaises(ProgModelInputException):
+        try:
             m.simulate_to(12, 132, {'o1': 0.8})
+            self.fail("Should have failed- future_load should be callable")
+        except ProgModelInputException:
+            pass
 
         ## Simulate
         (times, inputs, states, outputs, event_states) = m.simulate_to(3.5, load, {'o1': 0.8}, **{'dt': 0.5, 'save_freq': 1.0})
@@ -814,37 +877,61 @@ def load(t, x=None):
         
         ## Check inputs
         config = {'dt': [1, 2]}
-        with self.assertRaises(ProgModelInputException):
+        try:
             (times, inputs, states, outputs, event_states) = m.simulate_to(0, load, {'o1': 0.8}, **config)
+            self.fail("should have failed - dt must be number")
+        except ProgModelInputException:
+            pass
 
         config = {'dt': -1}
-        with self.assertRaises(ProgModelInputException):
+        try:
             (times, inputs, states, outputs, event_states) = m.simulate_to(0, load, {'o1': 0.8}, **config)
+            self.fail("Should have failed- dt must be positive")
+        except ProgModelInputException:
+            pass
 
         config = {'save_freq': [1, 2]}
-        with self.assertRaises(ProgModelInputException):
+        try:
             (times, inputs, states, outputs, event_states) = m.simulate_to(0, load, {'o1': 0.8}, **config)
+            self.fail("Should have failed- save_freq must be number")
+        except ProgModelInputException:
+            pass
 
         config = {'save_freq': -1}
-        with self.assertRaises(ProgModelInputException):
+        try:
             (times, inputs, states, outputs, event_states) = m.simulate_to(0, load, {'o1': 0.8}, **config)
+            self.fail("Should have failed- save_freq must be positive")
+        except ProgModelInputException:
+            pass
 
         config = {'horizon': [1, 2]}
-        with self.assertRaises(ProgModelInputException):
+        try:
             (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(load, {'o1': 0.8}, **config)
+            self.fail("Should have failed Horizon should be number")
+        except ProgModelInputException:
+            pass
 
         config = {'horizon': -1}
-        with self.assertRaises(ProgModelInputException):
+        try:
             (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(load, {'o1': 0.8}, **config)
+            self.fail("Should have failed- horizon must be positive")
+        except ProgModelInputException:
+            pass
         
         config = {'thresholds_met_eqn': -1}
-        with self.assertRaises(ProgModelInputException):
+        try:
             (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(load, {'o1': 0.8}, **config)
+            self.fail("Should have failed- thresholds_met_eqn must be callable")
+        except ProgModelInputException:
+            pass
 
         # incorrect number of arguments
         config = {'thresholds_met_eqn': lambda a, b: print(a, b)}
-        with self.assertRaises(ProgModelInputException):
+        try:
             (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(load, {'o1': 0.8}, **config)
+            self.fail()
+        except ProgModelInputException:
+            pass
 
     def test_sim_modes(self):
         m = ThrownObject(process_noise = 0, measurement_noise = 0)
@@ -882,10 +969,8 @@ def load(t, x=None):
         
         with self.assertRaises(ProgModelException):
             m.simulate_to_threshold(load, integration_method='rk4')
-
         # With linear model
         m = LinearThrownObject(process_noise = 0, measurement_noise = 0)
-
         result = m.simulate_to_threshold(load, dt = 0.1, integration_method='rk4')
         self.assertAlmostEqual(result.times[-1], 8.3)
 
@@ -955,30 +1040,45 @@ def load(t, x=None):
         self.assertListEqual(list(x['t']), [50, -100, 100])
 
         # when state doesn't exist
-        with self.assertRaises(Exception):
+        try:
             x0['n'] = 0
             (times, inputs, states, outputs, event_states) = m.simulate_to(0.001, load, {'o1': 0.8}, x = x0)
+            self.fail()
+        except Exception:
+            pass
 
         # when state entered incorrectly
-        with self.assertRaises(Exception):
+        try:
             x0['t'] = 'f'
             (times, inputs, states, outputs, event_states) = m.simulate_to(0.001, load, {'o1': 0.8}, x = x0)
+            self.fail()
+        except Exception:
+            pass
 
         # when boundary entered incorrectly
-        with self.assertRaises(Exception):
+        try:
             m.state_limits = { 't': ('f', 100) }
             x0['t'] = 0
             (times, inputs, states, outputs, event_states) = m.simulate_to(0.001, load, {'o1': 0.8}, x = x0)
+            self.fail()
+        except Exception:
+            pass
 
-        with self.assertRaises(Exception):
+        try:
             m.state_limits = { 't': (-100, 'f') }
             x0['t'] = 0
             (times, inputs, states, outputs, event_states) = m.simulate_to(0.001, load, {'o1': 0.8}, x = x0)
+            self.fail()
+        except Exception:
+            pass
 
-        with self.assertRaises(Exception):
+        try:
             m.state_limits = { 't': (100) }
             x0['t'] = 0
             (times, inputs, states, outputs, event_states) = m.simulate_to(0.001, load, {'o1': 0.8}, x = x0)
+            self.fail()
+        except Exception:
+            pass
 
     def test_progress_bar(self):
         m = MockProgModel(process_noise = 0.0)
@@ -989,6 +1089,7 @@ def load(t, x=None):
         capturedOutput = io.StringIO()
         sys.stdout = capturedOutput
 
+
         # Test progress bar matching
         simulate_results = m.simulate_to_threshold(load, {'o1': 0.8}, **{'dt': 0.5, 'save_freq': 1.0}, print=False, progress=True)
         capture_split =  [l+"%" for l in capturedOutput.getvalue().split("%") if l][:11]
@@ -1001,17 +1102,22 @@ def test_containers(self):
         m = ThrownObject()
         c1 = m.StateContainer({'x': 1.7, 'v': 40})
         c2 = m.StateContainer(np.array([[1.7], [40]]))
+        self.assertTrue(c1.data.equals(c2.data))
+        self.assertTrue(c2.data.equals(c1.data))
         self.assertEqual(c1, c2)
         self.assertListEqual(list(c1.keys()), m.states)
 
         input_c1 = m.InputContainer({})
         input_c2 = m.InputContainer(np.array([]))
         self.assertEqual(input_c1, input_c2)
+        # Issue with index in pd.DataFrame, check with empty
+        self.assertTrue(input_c1.data.empty and input_c2.data.empty)
         self.assertListEqual(list(input_c1.keys()), m.inputs)
 
         output_c1 = m.OutputContainer({'x': 1.7})
         output_c2 = m.OutputContainer(np.array([[1.7]]))
         self.assertEqual(output_c1, output_c2)
+        self.assertTrue(output_c1.data.equals(output_c2.data))
         self.assertListEqual(list(output_c1.keys()), m.outputs)
 
     def test_thrown_object_drag(self):
@@ -1034,19 +1140,25 @@ def future_load(t, x=None):
         # Test no drag simulated results different from default
         self.assertNotEqual(simulated_results_nd.times, simulated_results_df.times)
         self.assertNotEqual(simulated_results_nd.states, simulated_results_df.states)
+        self.assertFalse(simulated_results_nd.states.frame.equals(simulated_results_df.states.frame))
         self.assertGreater(simulated_results_nd.times[-1], simulated_results_df.times[-1])
+        self.assertGreater(simulated_results_nd.states.frame['time'].iloc[-1], simulated_results_df.states.frame['time'].iloc[-1])
 
         # Test high drag simulated results different from default
         self.assertNotEqual(simulated_results_hi.times, simulated_results_df.times)
         self.assertNotEqual(simulated_results_hi.states, simulated_results_df.states)
+        self.assertFalse(simulated_results_hi.states.frame.equals(simulated_results_df.states.frame))
         self.assertLess(simulated_results_hi.times[-1], simulated_results_df.times[-1])
+        self.assertLess(simulated_results_hi.states.frame['time'].iloc[-1], simulated_results_df.states.frame['time'].iloc[-1])
 
         # Test high drag simulated results different from no drag
         self.assertNotEqual(simulated_results_hi.times, simulated_results_nd.times)
         self.assertNotEqual(simulated_results_hi.states, simulated_results_nd.states)
+        self.assertFalse(simulated_results_hi.states.frame.equals(simulated_results_nd.states.frame))
 
     def test_composite_broken(self):
         m1 = OneInputOneOutputNoEventLM()
+
         # Insufficient number of models
         with self.assertRaises(ValueError):
             CompositeModel([])
@@ -1105,7 +1217,7 @@ def test_composite_broken(self):
             # extra
             CompositeModel([m1, m1], outputs=['OneInputOneOutputNoEventLM.z1', 'OneInputOneOutputNoEventLM_2.z1', 'z1'])
 
-    def test_composite(self):        
+    def test_composite(self):
         m1 = OneInputOneOutputNoEventLM()
         m2 = OneInputNoOutputOneEventLM()
         m1_withpm = OneInputOneOutputNoEventLMPM()
@@ -1162,22 +1274,37 @@ def test_composite(self):
         self.assertEqual(x0['OneInputOneOutputNoEventLM_2.x1'], 0)
         self.assertEqual(x0['OneInputOneOutputNoEventLM.x1'], 0)
         self.assertEqual(x0['OneInputOneOutputNoEventLM.z1'], 0)
+        # DataFrame check, x0
+        self.assertEqual(x0.data['OneInputOneOutputNoEventLM_2.x1'][0], 0)
+        self.assertEqual(x0.data['OneInputOneOutputNoEventLM.x1'][0], 0)
+        self.assertEqual(x0.data['OneInputOneOutputNoEventLM.z1'][0], 0)
         # Only provide non-zero input for first model
         u = m_composite.InputContainer({'OneInputOneOutputNoEventLM.u1': 1})
         x = m_composite.next_state(x0, u, 1)
         self.assertSetEqual(set(x.keys()), {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1', 'OneInputOneOutputNoEventLM.z1'})
         self.assertEqual(x['OneInputOneOutputNoEventLM_2.x1'], 1) # Propogates through, because of the order. If the connection were the other way it wouldn't
         self.assertEqual(x['OneInputOneOutputNoEventLM.x1'], 1)
+        # DataFrame check, x
+        self.assertSetEqual(set(x.data.columns.to_list()), {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1', 'OneInputOneOutputNoEventLM.z1'})
+        self.assertEqual(x.data['OneInputOneOutputNoEventLM_2.x1'][0], 1)
+        self.assertEqual(x.data['OneInputOneOutputNoEventLM.x1'][0], 1)
         z = m_composite.output(x)
         self.assertSetEqual(set(z.keys()), {'OneInputOneOutputNoEventLM_2.z1', 'OneInputOneOutputNoEventLM.z1'})
         self.assertEqual(z['OneInputOneOutputNoEventLM_2.z1'], 1)
         self.assertEqual(z['OneInputOneOutputNoEventLM.z1'], 1)
+        # DataFrame Check, z
+        self.assertSetEqual(set(z.data.columns.to_list()), {'OneInputOneOutputNoEventLM_2.z1', 'OneInputOneOutputNoEventLM.z1'})
+        self.assertEqual(z.data['OneInputOneOutputNoEventLM_2.z1'][0], 1)
+        self.assertEqual(z.data['OneInputOneOutputNoEventLM.z1'][0], 1)
 
         # Propogate again
         x = m_composite.next_state(x, u, 1)
         self.assertSetEqual(set(x.keys()), {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1', 'OneInputOneOutputNoEventLM.z1'})
         self.assertEqual(x['OneInputOneOutputNoEventLM_2.x1'], 3) # 1 + 2
         self.assertEqual(x['OneInputOneOutputNoEventLM.x1'], 2)
+        # DataFrame Check, x
+        self.assertEqual(x.data['OneInputOneOutputNoEventLM_2.x1'][0], 3)
+        self.assertEqual(x.data['OneInputOneOutputNoEventLM.x1'][0], 2)
 
         # Test with connections - state, no event
         m_composite = CompositeModel([m1, m1], connections=[('OneInputOneOutputNoEventLM.x1', 'OneInputOneOutputNoEventLM_2.u1')])
@@ -1187,7 +1314,6 @@ def test_composite(self):
         self.assertSetEqual(m_composite.inputs, {'OneInputOneOutputNoEventLM.u1',})
         self.assertSetEqual(m_composite.outputs, {'OneInputOneOutputNoEventLM.z1', 'OneInputOneOutputNoEventLM_2.z1'})
         self.assertSetEqual(m_composite.events, set())
-        
         x0 = m_composite.initialize()
         self.assertSetEqual(set(x0.keys()), {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1'})
         self.assertEqual(x0['OneInputOneOutputNoEventLM_2.x1'], 0)
@@ -1197,15 +1323,25 @@ def test_composite(self):
         x = m_composite.next_state(x0, u, 1)
         self.assertEqual(x['OneInputOneOutputNoEventLM_2.x1'], 1) # Propogates through, because of the order. If the connection were the other way it wouldn't
         self.assertEqual(x['OneInputOneOutputNoEventLM.x1'], 1)
+        # DataFrame Check, x
+        self.assertEqual(x.data['OneInputOneOutputNoEventLM_2.x1'][0], 1)
+        self.assertEqual(x.data['OneInputOneOutputNoEventLM.x1'][0], 1)
         z = m_composite.output(x)
         self.assertEqual(z['OneInputOneOutputNoEventLM_2.z1'], 1)
         self.assertEqual(z['OneInputOneOutputNoEventLM.z1'], 1)
+        # DataFrame Check, z
+        self.assertEqual(z.data['OneInputOneOutputNoEventLM_2.z1'][0], 1)
+        self.assertEqual(z.data['OneInputOneOutputNoEventLM.z1'][0], 1)
 
         # Propogate again
         x = m_composite.next_state(x, u, 1)
         self.assertSetEqual(set(x.keys()), {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1'})
         self.assertEqual(x['OneInputOneOutputNoEventLM_2.x1'], 3) # 1 + 2
         self.assertEqual(x['OneInputOneOutputNoEventLM.x1'], 2)
+        # DataFrame checks, x
+        self.assertSetEqual(set(x.data.columns.to_list()), {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1'})
+        self.assertEqual(x.data['OneInputOneOutputNoEventLM_2.x1'][0], 3)
+        self.assertEqual(x.data['OneInputOneOutputNoEventLM.x1'][0], 2)
 
         # Test with connections - two events
         m_composite = CompositeModel([m2, m2], connections=[('OneInputNoOutputOneEventLM.x1', 'OneInputNoOutputOneEventLM_2.u1')])
@@ -1234,7 +1370,6 @@ def test_composite(self):
         self.assertSetEqual(set(tm.keys()), {'OneInputNoOutputOneEventLM.x1 == 10', 'OneInputNoOutputOneEventLM_2.x1 == 10'})
         self.assertFalse(tm['OneInputNoOutputOneEventLM.x1 == 10'])
         self.assertTrue(tm['OneInputNoOutputOneEventLM_2.x1 == 10'])
-
         # Test with outputs specified
         m_composite = CompositeModel([m1, m1], connections=[('OneInputOneOutputNoEventLM.x1', 'OneInputOneOutputNoEventLM_2.u1')], outputs=['OneInputOneOutputNoEventLM_2.z1'])
         self.assertSetEqual(m_composite.states, {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1'})
@@ -1252,32 +1387,6 @@ def test_composite(self):
         self.assertSetEqual(m_composite.inputs, {'m1.u1',})
         self.assertSetEqual(m_composite.outputs, {'m1.z1', })
         self.assertSetEqual(m_composite.events, {'m2.x1 == 10', })
-    
-    # Fill parameters with different types of objects instead
-    def test_parameter_equality(self):
-
-        m1 = LinearThrownObject()
-        m2 = LinearThrownObject()
-
-        self.assertTrue(m1.parameters == m2.parameters) #Checking to see if the parameters are equal
-        self.assertTrue(m2.parameters == m1.parameters) #Parameters should be equal
-
-        m3 = LinearThrownDiffThrowingSpeed() # A model with a different throwing speed
-        self.assertFalse(m1.parameters == m3.parameters)
-        self.assertFalse(m3.parameters == m1.parameters) # Checking both directions 
-
-        m4 = LinearThrownObjectDiffDefaultParams() # Model with an extra default parameter.
-
-        self.assertFalse(m1.parameters == m4.parameters)
-        self.assertFalse(m4.parameters == m1.parameters) # checking both directions
-
-        m5 = LinearThrownObjectUpdatedInitalizedMethod() # Model with incorrectly initalized throwing height, but same parameters
-
-        self.assertFalse(m1.parameters == m5.parameters) 
-        self.assertFalse(m5.parameters == m1.parameters) 
-
-        self.assertTrue(m1.parameters == m2.parameters) # Checking to see previous equal statements stay the same
-        self.assertTrue(m2.parameters == m1.parameters)
 
 # This allows the module to be executed directly
 def run_tests():
diff --git a/tests/test_composite.py b/tests/test_composite.py
deleted file mode 100644
index c10d20d93..000000000
--- a/tests/test_composite.py
+++ /dev/null
@@ -1,243 +0,0 @@
-# Copyright © 2021 United States Government as represented by the Administrator of the National Aeronautics and Space Administration.  All Rights Reserved.
-
-from copy import deepcopy
-import io
-import numpy as np
-from os.path import dirname, join
-import pickle
-import sys
-import unittest
-
-# This ensures that the directory containing ProgModelTemplate is in the python search directory
-sys.path.append(join(dirname(__file__), ".."))
-
-from prog_models import *
-from prog_models.models import *
-from prog_models.models.test_models.linear_models import (
-    OneInputNoOutputNoEventLM, OneInputOneOutputNoEventLM, OneInputNoOutputOneEventLM, OneInputOneOutputNoEventLMPM)
-from prog_models.models.thrown_object import LinearThrownObject
-
-class TestCompositeModel(unittest.TestCase):
-    def test_composite_broken(self):
-        m1 = OneInputOneOutputNoEventLM()
-
-        # Insufficient number of models
-        with self.assertRaises(ValueError):
-            CompositeModel([])
-        with self.assertRaises(ValueError):
-            CompositeModel([m1])
-        
-        # Wrong type
-        with self.assertRaises(ValueError):
-            CompositeModel([m1, m1, 'abc'])
-
-        # Incorrect named format
-        with self.assertRaises(ValueError):
-            # Too many elements
-            CompositeModel([('a', m1, 'Something else'), ('b', m1)])
-        with self.assertRaises(ValueError):
-            # Not a string
-            CompositeModel([(m1, m1)])
-        with self.assertRaises(ValueError):
-            # Not a model
-            CompositeModel([('a', 'b')])
-        with self.assertRaises(ValueError):
-            # Too few elements
-            CompositeModel([(m1, )])
-
-        # Incorrect connections
-        with self.assertRaises(ValueError):
-            # without model name
-            CompositeModel([m1, m1], connections=[('z1', 'u1')])
-        with self.assertRaises(ValueError):
-            # broken in
-            CompositeModel([m1, m1], connections=[('z1', 'OneInputOneOutputNoEventLM.u1')])
-        with self.assertRaises(ValueError):
-            # broken out
-            CompositeModel([m1, m1], connections=[('OneInputOneOutputNoEventLM.z1', 'u1')])
-        with self.assertRaises(ValueError):
-            # Switched
-            CompositeModel([m1, m1], connections=[('OneInputOneOutputNoEventLM.u1', 'OneInputOneOutputNoEventLM_2.z1')])
-        with self.assertRaises(ValueError):
-            # Improper format - too long
-            CompositeModel([m1, m1], connections=[('OneInputOneOutputNoEventLM.z1', 'OneInputOneOutputNoEventLM.u1', 'Something else')])
-        with self.assertRaises(ValueError):
-            # Improper format - not a string
-            CompositeModel([m1, m1], connections=[(m1, m1)])
-        with self.assertRaises(ValueError):
-            # Improper format - too short
-            CompositeModel([m1, m1], connections=[('OneInputOneOutputNoEventLM.z1', )])
-        with self.assertRaises(ValueError):
-            # Improper format - not a tuple
-            CompositeModel([m1, m1], connections=['m1'])
-
-        # Incorrect outputs
-        with self.assertRaises(ValueError):
-            # without model name
-            CompositeModel([m1, m1], outputs=['z1'])
-        with self.assertRaises(ValueError):
-            # extra
-            CompositeModel([m1, m1], outputs=['OneInputOneOutputNoEventLM.z1', 'OneInputOneOutputNoEventLM_2.z1', 'z1'])
-
-    def test_composite(self):
-        m1 = OneInputOneOutputNoEventLM()
-        m2 = OneInputNoOutputOneEventLM()
-        m1_withpm = OneInputOneOutputNoEventLMPM()
-
-        # Test with no connections
-        m_composite = CompositeModel([m1, m1])
-        self.assertSetEqual(m_composite.states, {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1'})
-        self.assertSetEqual(m_composite.inputs, {'OneInputOneOutputNoEventLM.u1', 'OneInputOneOutputNoEventLM_2.u1'})
-        self.assertSetEqual(m_composite.outputs, {'OneInputOneOutputNoEventLM.z1', 'OneInputOneOutputNoEventLM_2.z1'})
-        self.assertSetEqual(m_composite.events, set())
-        self.assertSetEqual(m_composite.performance_metric_keys, set(), "Shouldn't have any performance metrics")
-
-        x0 = m_composite.initialize()
-        self.assertSetEqual(set(x0.keys()), {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1'})
-        self.assertEqual(x0['OneInputOneOutputNoEventLM_2.x1'], 0)
-        self.assertEqual(x0['OneInputOneOutputNoEventLM.x1'], 0)
-        # Only provide non-zero input for the first model
-        u = m_composite.InputContainer({'OneInputOneOutputNoEventLM.u1': 1, 'OneInputOneOutputNoEventLM_2.u1': 0})
-        x = m_composite.next_state(x0, u, 1)
-        self.assertSetEqual(set(x.keys()), {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1'})
-        self.assertEqual(x['OneInputOneOutputNoEventLM_2.x1'], 0)
-        self.assertEqual(x['OneInputOneOutputNoEventLM.x1'], 1)
-        z = m_composite.output(x)
-        self.assertSetEqual(set(z.keys()), {'OneInputOneOutputNoEventLM_2.z1', 'OneInputOneOutputNoEventLM.z1'})
-        self.assertEqual(z['OneInputOneOutputNoEventLM_2.z1'], 0)
-        self.assertEqual(z['OneInputOneOutputNoEventLM.z1'], 1)
-        pm = m_composite.performance_metrics(x)
-        self.assertSetEqual(set(pm.keys()), set())
-
-        # With Performance Metrics
-        # Everything else should behave the same, so we're only testing the performance metrics
-        m_composite = CompositeModel([m1_withpm, m1_withpm])
-        self.assertSetEqual(m_composite.performance_metric_keys, {'OneInputOneOutputNoEventLMPM_2.x1+1', 'OneInputOneOutputNoEventLMPM.x1+1'})
-
-        x0 = m_composite.initialize()
-        u = m_composite.InputContainer({'OneInputOneOutputNoEventLMPM.u1': 1, 'OneInputOneOutputNoEventLMPM_2.u1': 0})
-        x = m_composite.next_state(x0, u, 1)
-        pm = m_composite.performance_metrics(x)
-        self.assertSetEqual(set(pm.keys()), {'OneInputOneOutputNoEventLMPM_2.x1+1', 'OneInputOneOutputNoEventLMPM.x1+1'})
-        self.assertEqual(pm['OneInputOneOutputNoEventLMPM_2.x1+1'], 1)
-        self.assertEqual(pm['OneInputOneOutputNoEventLMPM.x1+1'], 2)
-
-        # Test with connections - output, no event
-        m_composite = CompositeModel([m1, m1], connections=[('OneInputOneOutputNoEventLM.z1', 'OneInputOneOutputNoEventLM_2.u1')])
-        # Additional state to store output
-        self.assertSetEqual(m_composite.states, {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1', 'OneInputOneOutputNoEventLM.z1'})
-        # One less input - since it's internally connected
-        self.assertSetEqual(m_composite.inputs, {'OneInputOneOutputNoEventLM.u1',})
-        self.assertSetEqual(m_composite.outputs, {'OneInputOneOutputNoEventLM.z1', 'OneInputOneOutputNoEventLM_2.z1'})
-        self.assertSetEqual(m_composite.events, set())
-
-        x0 = m_composite.initialize()
-        self.assertSetEqual(set(x0.keys()), {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1', 'OneInputOneOutputNoEventLM.z1'})
-        self.assertEqual(x0['OneInputOneOutputNoEventLM_2.x1'], 0)
-        self.assertEqual(x0['OneInputOneOutputNoEventLM.x1'], 0)
-        self.assertEqual(x0['OneInputOneOutputNoEventLM.z1'], 0)
-        # Only provide non-zero input for first model
-        u = m_composite.InputContainer({'OneInputOneOutputNoEventLM.u1': 1})
-        x = m_composite.next_state(x0, u, 1)
-        self.assertSetEqual(set(x.keys()), {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1', 'OneInputOneOutputNoEventLM.z1'})
-        self.assertEqual(x['OneInputOneOutputNoEventLM_2.x1'], 1) # Propogates through, because of the order. If the connection were the other way it wouldn't
-        self.assertEqual(x['OneInputOneOutputNoEventLM.x1'], 1)
-        z = m_composite.output(x)
-        self.assertSetEqual(set(z.keys()), {'OneInputOneOutputNoEventLM_2.z1', 'OneInputOneOutputNoEventLM.z1'})
-        self.assertEqual(z['OneInputOneOutputNoEventLM_2.z1'], 1)
-        self.assertEqual(z['OneInputOneOutputNoEventLM.z1'], 1)
-
-        # Propogate again
-        x = m_composite.next_state(x, u, 1)
-        self.assertSetEqual(set(x.keys()), {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1', 'OneInputOneOutputNoEventLM.z1'})
-        self.assertEqual(x['OneInputOneOutputNoEventLM_2.x1'], 3) # 1 + 2
-        self.assertEqual(x['OneInputOneOutputNoEventLM.x1'], 2)
-
-        # Test with connections - state, no event
-        m_composite = CompositeModel([m1, m1], connections=[('OneInputOneOutputNoEventLM.x1', 'OneInputOneOutputNoEventLM_2.u1')])
-        # No additional state to store output, since state is used for the connection
-        self.assertSetEqual(m_composite.states, {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1'})
-        # One less input - since it's internally connected
-        self.assertSetEqual(m_composite.inputs, {'OneInputOneOutputNoEventLM.u1',})
-        self.assertSetEqual(m_composite.outputs, {'OneInputOneOutputNoEventLM.z1', 'OneInputOneOutputNoEventLM_2.z1'})
-        self.assertSetEqual(m_composite.events, set())
-        
-        x0 = m_composite.initialize()
-        self.assertSetEqual(set(x0.keys()), {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1'})
-        self.assertEqual(x0['OneInputOneOutputNoEventLM_2.x1'], 0)
-        self.assertEqual(x0['OneInputOneOutputNoEventLM.x1'], 0)
-        # Only provide non-zero input for model 1
-        u = m_composite.InputContainer({'OneInputOneOutputNoEventLM.u1': 1})
-        x = m_composite.next_state(x0, u, 1)
-        self.assertEqual(x['OneInputOneOutputNoEventLM_2.x1'], 1) # Propogates through, because of the order. If the connection were the other way it wouldn't
-        self.assertEqual(x['OneInputOneOutputNoEventLM.x1'], 1)
-        z = m_composite.output(x)
-        self.assertEqual(z['OneInputOneOutputNoEventLM_2.z1'], 1)
-        self.assertEqual(z['OneInputOneOutputNoEventLM.z1'], 1)
-
-        # Propogate again
-        x = m_composite.next_state(x, u, 1)
-        self.assertSetEqual(set(x.keys()), {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1'})
-        self.assertEqual(x['OneInputOneOutputNoEventLM_2.x1'], 3) # 1 + 2
-        self.assertEqual(x['OneInputOneOutputNoEventLM.x1'], 2)
-
-        # Test with connections - two events
-        m_composite = CompositeModel([m2, m2], connections=[('OneInputNoOutputOneEventLM.x1', 'OneInputNoOutputOneEventLM_2.u1')])
-        self.assertSetEqual(m_composite.states, {'OneInputNoOutputOneEventLM_2.x1', 'OneInputNoOutputOneEventLM.x1'})
-        # One less input - since it's internally connected
-        self.assertSetEqual(m_composite.inputs, {'OneInputNoOutputOneEventLM.u1',})
-        self.assertSetEqual(m_composite.outputs, set())
-        self.assertSetEqual(m_composite.events, {'OneInputNoOutputOneEventLM.x1 == 10', 'OneInputNoOutputOneEventLM_2.x1 == 10'})
-
-        x0 = m_composite.initialize()
-        u = m_composite.InputContainer({'OneInputNoOutputOneEventLM.u1': 1})
-        x = m_composite.next_state(x0, u, 1) # 1, 1
-        x = m_composite.next_state(x, u, 1) # 2, 3
-        x = m_composite.next_state(x, u, 1) # 3, 6
-        tm = m_composite.threshold_met(x)
-        self.assertSetEqual(set(tm.keys()), {'OneInputNoOutputOneEventLM.x1 == 10', 'OneInputNoOutputOneEventLM_2.x1 == 10'})
-        self.assertFalse(tm['OneInputNoOutputOneEventLM.x1 == 10'])
-        self.assertFalse(tm['OneInputNoOutputOneEventLM_2.x1 == 10'])
-
-        x = m_composite.next_state(x, u, 1) # 4, 10
-        es = m_composite.event_state(x)
-        self.assertSetEqual(set(es.keys()), {'OneInputNoOutputOneEventLM.x1 == 10', 'OneInputNoOutputOneEventLM_2.x1 == 10'})
-        self.assertEqual(es['OneInputNoOutputOneEventLM.x1 == 10'], 0.6)
-        self.assertEqual(es['OneInputNoOutputOneEventLM_2.x1 == 10'], 0.0)
-        tm = m_composite.threshold_met(x)
-        self.assertSetEqual(set(tm.keys()), {'OneInputNoOutputOneEventLM.x1 == 10', 'OneInputNoOutputOneEventLM_2.x1 == 10'})
-        self.assertFalse(tm['OneInputNoOutputOneEventLM.x1 == 10'])
-        self.assertTrue(tm['OneInputNoOutputOneEventLM_2.x1 == 10'])
-
-        # Test with outputs specified
-        m_composite = CompositeModel([m1, m1], connections=[('OneInputOneOutputNoEventLM.x1', 'OneInputOneOutputNoEventLM_2.u1')], outputs=['OneInputOneOutputNoEventLM_2.z1'])
-        self.assertSetEqual(m_composite.states, {'OneInputOneOutputNoEventLM_2.x1', 'OneInputOneOutputNoEventLM.x1'})
-        self.assertSetEqual(m_composite.inputs, {'OneInputOneOutputNoEventLM.u1',})
-        # One less output
-        self.assertSetEqual(set(m_composite.outputs), {'OneInputOneOutputNoEventLM_2.z1', })
-        self.assertSetEqual(m_composite.events, set())
-        x0 = m_composite.initialize()
-        z = m_composite.output(x0)
-        self.assertSetEqual(set(z.keys()), {'OneInputOneOutputNoEventLM_2.z1', })
-
-        # With Names
-        m_composite = CompositeModel([('m1', m1), ('m2', m2)], connections=[('m1.x1', 'm2.u1')])
-        self.assertSetEqual(m_composite.states, {'m2.x1', 'm1.x1'})
-        self.assertSetEqual(m_composite.inputs, {'m1.u1',})
-        self.assertSetEqual(m_composite.outputs, {'m1.z1', })
-        self.assertSetEqual(m_composite.events, {'m2.x1 == 10', })
-
-def run_tests():
-    unittest.main()
-    
-def main():
-    l = unittest.TestLoader()
-    runner = unittest.TextTestRunner()
-    print("\n\nTesting Composite Models")
-    result = runner.run(l.loadTestsFromTestCase(TestCompositeModel)).wasSuccessful()
-
-    if not result:
-        raise Exception("Failed test")
-
-if __name__ == '__main__':
-    main()
diff --git a/tests/test_data_model.py b/tests/test_data_model.py
index 8b51e070a..4a72e2b2c 100644
--- a/tests/test_data_model.py
+++ b/tests/test_data_model.py
@@ -140,23 +140,22 @@ def test_lstm_simple(self):
                 self.assertIsInstance(m2, LSTMStateTransitionModel)
                 self.assertIsInstance(m2, DataModel)
                 self.assertListEqual(m2.outputs, ['x'])
-
-            # Deepcopy test
-            m3 = deepcopy(m2)
         except:
             warnings.warn("Pickling not supported for LSTMStateTransitionModel on this system")
             pass
 
+        # Deepcopy test
+        m3 = deepcopy(m2)
         # More tests in examples.lstm_model
 
     def test_dmd_simple(self):
         self._test_simple_case(DMDModel, max_error=25)
 
         # Inferring dt
-        self._test_simple_case(DMDModel, max_error=8, WITH_DT=False)
+        self._test_simple_case(DMDModel, max_error=8, WITH_DT = False)
 
         # Without velocity, DMD doesn't perform well
-        m = self._test_simple_case(DMDModel, WITH_STATES=False, max_error=100)
+        m = self._test_simple_case(DMDModel, WITH_STATES = False, max_error=100)
 
         # Test pickling model m
         pickled_m = pickle.dumps(m)
@@ -179,19 +178,18 @@ def future_loading(t, x=None):
         m3 = LSTMStateTransitionModel.from_model(
             m,
             [future_loading for _ in range(5)],
-            dt=[TIMESTEP, TIMESTEP/2, TIMESTEP/4, TIMESTEP*2, TIMESTEP*4],
-            window=2,
-            epochs=20)
+            dt = [TIMESTEP, TIMESTEP/2, TIMESTEP/4, TIMESTEP*2, TIMESTEP*4],
+            window=2, 
+            epochs=20)  
 
         # Should get keys from original model
         self.assertSetEqual(set(m3.inputs), set(['dt', 'x_t-1']))
         self.assertSetEqual(set(m3.outputs), set(m.outputs))
 
-        # Step 3: Use model to simulate_to time of threshold
+         # Step 3: Use model to simulate_to time of threshold
         t_counter = 0
         x_counter = m.initialize()
-
-        def future_loading2(t, x=None):
+        def future_loading2(t, x = None):
             # Future Loading is a bit complicated here 
             # Loading for the resulting model includes the data inputs, 
             # and the output from the last timestep
diff --git a/tests/test_ensemble.py b/tests/test_ensemble.py
deleted file mode 100644
index ce4e216a9..000000000
--- a/tests/test_ensemble.py
+++ /dev/null
@@ -1,213 +0,0 @@
-# Copyright © 2021 United States Government as represented by the Administrator of the
-# National Aeronautics and Space Administration.  All Rights Reserved.
-
-from io import StringIO
-import numpy as np
-import sys
-import unittest
-
-from prog_models import EnsembleModel
-from prog_models.models.test_models.linear_models import OneInputOneOutputOneEventLM, OneInputOneOutputOneEventAltLM
-
-
-class TestEnsemble(unittest.TestCase):
-    def setUp(self):
-        # set stdout (so it wont print)
-        sys.stdout = StringIO()
-
-    def tearDown(self):
-        sys.stdout = sys.__stdout__
-
-    def test_no_model(self):
-        with self.assertRaises(ValueError):
-            EnsembleModel([])
-
-    def test_single_model(self):
-        # An ensemble model with one model should raise an exception
-
-        m = OneInputOneOutputOneEventLM()
-        
-        with self.assertRaises(ValueError):
-            EnsembleModel([m])
-
-    def test_wrong_type(self):
-        # An ensemble model with a non-model should raise an exception
-
-        m = OneInputOneOutputOneEventLM()
-        with self.assertRaises(TypeError):
-            EnsembleModel(m)
-        with self.assertRaises(TypeError):
-            EnsembleModel([m, 1])
-        with self.assertRaises(TypeError):
-            EnsembleModel(77)
-        with self.assertRaises(TypeError):
-            EnsembleModel([m, m, m, 77])
-
-    def test_two_models_identical(self):
-        """
-        This tests that the ensemble model works with two identical model-types with identical states, inputs, etc., one with slightly altered parameters.
-
-        The result is that each state, output, etc. is combined using the specified aggregation method.
-        """
-        m = OneInputOneOutputOneEventLM()
-        m2 = OneInputOneOutputOneEventLM(x0={'x1': 2})
-
-        # Make sure they're not the same - matricies are 3x their original values
-        # The result is a model where state changes 3x as fast.
-        # Event state degrades 9x as fast, since B and F compound
-        m2.B = np.array([[3]])
-        m2.C = np.array([[3]])
-        m2.F = np.array([[-0.3]])
-
-        # An ensemble model with two models should work
-        em = EnsembleModel([m, m2])
-
-        # Since they're identical, the inputs, state, etc. should be the same
-        self.assertEqual(em.inputs, m.inputs)
-        self.assertEqual(em.states, m.states)
-        self.assertEqual(em.outputs, m.outputs)
-        self.assertEqual(em.events, m.events)
-
-        # The resulting initial state should be exactly between the two:
-        x_t0 = em.initialize()
-        self.assertEqual(x_t0['x1'], 1)
-
-        # Same with state transition
-        u = em.InputContainer({'u1': 1})
-        x_t1 = em.next_state(x_t0, u, 1)
-        # m would give 2, m2 would give 4
-        self.assertEqual(x_t1['x1'], 3)
-
-        # Same with output
-        z = em.output(x_t1)
-        # m would give 3, m2 would give 9
-        self.assertEqual(z['z1'], 6)
-
-        # Same with event state
-        es = em.event_state(x_t1)
-        # m would give 0.7, m2 would give 0.1
-        self.assertEqual(es['x1 == 10'], 0.4)
-
-        # performance metrics
-        pm = em.performance_metrics(x_t1)
-        self.assertEqual(pm['pm1'], 4)
-
-        # Time of event 
-        toe = em.time_of_event(x_t0, lambda t, x=None: u, dt=1e-3)
-        self.assertAlmostEqual(toe['x1 == 10'], 4.8895)
-
-        # threshold met should be false
-        self.assertFalse(em.threshold_met(x_t1)['x1 == 10'])
-
-        # Transition again
-        x_t2 = em.next_state(x_t1, u, 1)
-        # threshold met should be true (because one of 2 models says it is)
-        self.assertTrue(em.threshold_met(x_t2)['x1 == 10'])
-
-    def test_two_models_different(self):
-        """
-        This tests that the ensemble model works with two different model-types with different states, inputs, etc. Tests how the ensemble model handles the different values.
-
-        The result is that the different states, inputs, etc. are combined into a single set without being aggregated (unlike test_two_models_identical).
-        """
-        m = OneInputOneOutputOneEventLM()
-        m2 = OneInputOneOutputOneEventAltLM()
-        em = EnsembleModel([m, m2])
-
-        # inputs, states, outputs, events should be a combination of the two models
-        self.assertSetEqual(set(em.inputs), {'u1', 'u2'})
-        self.assertSetEqual(set(em.states), {'x1', 'x2'})
-        self.assertSetEqual(set(em.outputs), {'z1', 'z2'})
-        self.assertSetEqual(set(em.events), {'x1 == 10', 'x2 == 5'})
-
-        # Initialize - should be combination of the two
-        x_t0 = em.initialize()
-        self.assertEqual(x_t0['x1'], 0)
-        self.assertEqual(x_t0['x2'], 0)
-
-        # State transition - should be combination of the two
-        u = em.InputContainer({'u1': 1, 'u2': 2})
-        x_t1 = em.next_state(x_t0, u, 1)
-        self.assertEqual(x_t1['x1'], 1)
-        self.assertEqual(x_t1['x2'], 2)
-
-        # Output - should be combination of the two
-        z = em.output(x_t1)
-        self.assertEqual(z['z1'], 1)
-        self.assertEqual(z['z2'], 2)
-
-        # Event state - should be combination of the two
-        es = em.event_state(x_t1)
-        self.assertEqual(es['x1 == 10'], 0.9)
-        self.assertEqual(es['x2 == 5'], 0.6)
-
-        # Threshold met - should be combination of the two
-        self.assertFalse(em.threshold_met(x_t1)['x1 == 10'])
-        self.assertFalse(em.threshold_met(x_t1)['x2 == 5'])
-
-        # Transition again
-        x_t2 = em.next_state(x_t1, u, 2)
-
-        # Threshold met - should be combination of the two
-        # x1 == 3, x2 == 6
-        self.assertFalse(em.threshold_met(x_t2)['x1 == 10'])
-        self.assertTrue(em.threshold_met(x_t2)['x2 == 5'])
-
-    def test_two_models_alt_aggrigation(self):
-        """
-        This test repeats test_two_models_identical with different aggrigation method.
-        """
-        m = OneInputOneOutputOneEventLM()
-        m2 = OneInputOneOutputOneEventLM(x0={'x1': 2})
-
-        # Make sure they're not the same - 3x the impact
-        m2.B = np.array([[3]])
-        m2.C = np.array([[3]])
-        m2.F = np.array([[-0.3]])
-
-        # An ensemble model with two models should work
-        em = EnsembleModel([m, m2], aggregation_method=np.max)
-
-        # The resulting initial state should be max of the two:
-        x_t0 = em.initialize()
-        self.assertEqual(x_t0['x1'], 2)
-
-        # Same with state transition
-        u = em.InputContainer({'u1': 1})
-        x_t1 = em.next_state(x_t0, u, 1)
-        # m would give 3, m2 would give 5
-        self.assertEqual(x_t1['x1'], 5)
-
-        # Same with output
-        z = em.output(x_t1)
-        # m would give 5, m2 would give 15
-        self.assertEqual(z['z1'], 15)
-
-        # Same with event state
-        es = em.event_state(x_t1)
-        # m would give 0.5, m2 would give -0.5
-        self.assertEqual(es['x1 == 10'], 0.5)
-
-        # threshold met should be false
-        self.assertFalse(em.threshold_met(x_t1)['x1 == 10'])
-
-        # Next state
-        x2 = em.next_state(x_t1, u, 2)
-        # threshold met should be true (because both of the models agree)
-        self.assertTrue(em.threshold_met(x2)['x1 == 10'])
-
-# This allows the module to be executed directly
-def run_tests():
-    unittest.main()
-    
-def main():
-    l = unittest.TestLoader()
-    runner = unittest.TextTestRunner()
-    print("\n\nTesting Ensemble models")
-    result = runner.run(l.loadTestsFromTestCase(TestEnsemble)).wasSuccessful()
-
-    if not result:
-        raise Exception("Failed test")
-
-if __name__ == '__main__':
-    main()
diff --git a/tests/test_linear_model.py b/tests/test_linear_model.py
index f0eb2e85a..f544048c7 100644
--- a/tests/test_linear_model.py
+++ b/tests/test_linear_model.py
@@ -2,27 +2,19 @@
 
 import numpy as np
 import unittest
-import copy
-import pickle
-from prog_models.models.test_models.linear_thrown_object import (LinearThrownObject, LinearThrownObjectNoE, LinearThrownObjectWrongB,
-                                                                 LinearThrownDiffThrowingSpeed, LinearThrownObjectUpdatedInitalizedMethod,
-                                                                 LinearThrownObjectFourStates)
+
 from prog_models.models.test_models.linear_models import FNoneNoEventStateLM
+from prog_models.models.test_models.linear_thrown_object import LinearThrownObject
 
 
 class TestLinearModel(unittest.TestCase):
     def test_linear_model(self):
         m = LinearThrownObject()
-
-        #Checking to see if initalization would error when passing in incorrect parameter forms
-        with self.assertRaises(AttributeError):
-            b = LinearThrownObjectWrongB()
-
+        m.simulate_to_threshold(lambda t, x = None: m.InputContainer({}))
         # len() = events states inputs outputs
         #         1      2      0      1
         # Matrix overwrite type checking (Can't set attributes for B, D, G; not overwritten)
         # when matrix is not of type NumPy ndarray or standard list
-
         # @A
         with self.assertRaises(TypeError):
             m.A = "[[0, 1], [0, 0]]" # string
@@ -48,78 +40,13 @@ def test_linear_model(self):
         with self.assertRaises(TypeError):
             m.A = True # boolean
             m.matrixCheck()
-        # Matrix Dimension Checking
-        # when matrix is not proper dimensional (1-D array = C, D, G; 2-D array = A,B,E; None = F;)
-        with self.assertRaises(AttributeError):
-            m.A = np.array([0, 1]) # 1-D array
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.A = np.array([[[[0, 1], [0, 0], [1, 0]]]]) # 3-D array
-            m.matrixCheck()
-        # When Matrix is improperly formed
-        with self.assertRaises(AttributeError):
-            m.A = np.array([[0, 1, 2, 3], [0, 0, 1, 2]]) # extra column values per row
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.A = np.array([[0], [0]]) # less column values per row
-            m.matrixCheck()
-        with self.assertRaises(AttributeError): 
-            m.A = np.array([[0, 1], [0, 0], [2, 2]]) # extra row
-            m.matrixCheck()
-        with self.assertRaises(AttributeError): 
-            m.A = np.array([[0, 1]]) # less row
-            m.matrixCheck()
-        # Reset Process
-        m.A = np.array([[0, 1], [0, 2]])
-        m.matrixCheck()
-
-        # @B
-        with self.assertRaises(TypeError):
-            m.B = "[[0, 1], [0, 0]]"
-            m.matrixCheck()
-        with self.assertRaises(TypeError):
-            m.B = 0 # int
-            m.matrixCheck()
-        with self.assertRaises(TypeError):
-            m.B = 3.14 # float
-            m.matrixCheck()
-        with self.assertRaises(TypeError):
-            m.B = {} # dict
-            m.matrixCheck() 
-        with self.assertRaises(TypeError):
-            m.B = () # tuple
-            m.matrixCheck() 
-        with self.assertRaises(TypeError):
-            m.B = set() # set
-            m.matrixCheck() 
-        with self.assertRaises(TypeError):
-            m.B = True # boolean
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.B = np.array(2) # 0-D array
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.B = np.array([[0, 0], [1, 1]]) # 2-D array
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.B = np.array([[1, 0, 2]]) # extra column values per row
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.B = np.array([[1]]) # less column values per row
-            m.matrixCheck()
-        with self.assertRaises(AttributeError): 
-            m.B = np.array([[0, 0], [1, 1], [2, 2]]) # extra row
-            m.matrixCheck()
-        with self.assertRaises(AttributeError): 
-            m.B = np.array([[]]) # less row 
-            m.matrixCheck()
-        m.B = None #sets parameter B to default value
-        m.matrixCheck()
-
         # @C
         with self.assertRaises(TypeError):
             m.C = "[[0, 1], [0, 0]]" # string
             m.matrixCheck() 
+        with self.assertRaises(TypeError):
+            m.C = None # None
+            m.matrixCheck() 
         with self.assertRaises(TypeError):
             m.C = 0 # int
             m.matrixCheck()
@@ -138,72 +65,13 @@ def test_linear_model(self):
         with self.assertRaises(TypeError):
             m.C = True # boolean
             m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.C = np.array(2) # 0-D array
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.C = np.array([[0, 0], [1, 1]]) # 2-D array
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.C = np.array([[1, 0, 2]]) # extra column values per row
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.C = np.array([[1]]) # less column values per row
-            m.matrixCheck()
-        with self.assertRaises(AttributeError): 
-            m.C = np.array([[0, 0], [1, 1], [2, 2]]) # extra row
-            m.matrixCheck()
-        with self.assertRaises(AttributeError): 
-            m.C = np.array([[]]) # less row 
-            m.matrixCheck()
-        m.C = np.array([[1, 0]])
-        m.matrixCheck()
-
-
-# Included some tests that are checking if exceptions are being thrown without matrxiCheck being invoked
-        with self.assertRaises(TypeError):
-            m.D = "[[0, 1], [0, 0]]" # string
-            m.matrixCheck() 
-        with self.assertRaises(TypeError):
-            m.D = 0 # int
-            m.matrixCheck()
-        with self.assertRaises(TypeError):
-            m.D = 3.14 # float
-            m.matrixCheck()
-        with self.assertRaises(TypeError):
-            m.D = {} # dict
-            m.matrixCheck() 
-        with self.assertRaises(TypeError):
-            m.D = () # tuple
-        with self.assertRaises(TypeError):
-            m.D = set() # set
-        with self.assertRaises(TypeError):
-            m.D = True # boolean
-            m.matrixCheck()
-        # @D 1x
-        with self.assertRaises(AttributeError):
-            m.D = np.array(1) # 0-D array
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.D = np.array([[2], [1]]) # 2-D array with incorrect values passed in
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.D = np.array([1, 2]) # extra column values per row
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.D = np.array([[]]) # less column values per row
-        with self.assertRaises(AttributeError): 
-            m.D = np.array([[0], [1]]) # extra row
-        with self.assertRaises(AttributeError): 
-            m.D = np.array([[]]) # less row
-        m.D = np.array([[1]]) # sets to Default Value
-        m.D = None
-        m.matrixCheck()
-
         # @E
         with self.assertRaises(TypeError):
             m.E = "[[0, 1], [0, 0]]" # string
             m.matrixCheck() 
+        with self.assertRaises(TypeError):
+            m.E = None # None
+            m.matrixCheck() 
         with self.assertRaises(TypeError):
             m.E = 0 # int
             m.matrixCheck()
@@ -222,28 +90,6 @@ def test_linear_model(self):
         with self.assertRaises(TypeError):
             m.E = True # boolean
             m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.E = np.array([[0]]) # 2-D array
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.E = np.array([[0], [1], [2]]) # 3-D array
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.E = np.array([[0,0], [-9.81, -1]]) # extra column values per row
-            m.matrixCheck()
-        with self.assertRaises(AttributeError):
-            m.E = np.array([[], []]) # less column values per row
-            m.matrixCheck()
-        with self.assertRaises(AttributeError): 
-            m.E = np.array([[0, 1], [0, 0], [2, 2]]) # extra row
-            m.matrixCheck()
-        with self.assertRaises(AttributeError): 
-            m.E = np.array([[0, 1]]) # less row
-            m.matrixCheck()
-        m.E = np.array([[0], [-9.81]])
-        m.matrixCheck()
-        
-
         # @F
         with self.assertRaises(TypeError):
             m.F = "[[0, 1], [0, 0]]" # string
@@ -267,68 +113,114 @@ def test_linear_model(self):
             m.F = True # boolean
             m.matrixCheck()
         with self.assertRaises(AttributeError):
-            m.F = np.array([[0]]) # 2-D array
+            # if F is none, we need to override event_state
+            m_noes = FNoneNoEventStateLM()
+        
+        # Matrix Dimension Checking
+        # when matrix is not proper dimensional (1-D array = C, D, G; 2-D array = A,B,E; None = F;)
+        # @A 2x2
+        with self.assertRaises(AttributeError):
+            m.A = np.array([[0, 1]]) # 1-D array
+            m.matrixCheck()
+        with self.assertRaises(AttributeError):
+            m.A = np.array([[0, 1], [0, 0], [1, 0]]) # 3-D array
+            m.matrixCheck()
+        # @B 2x0
+        with self.assertRaises(AttributeError):
+            m.B = np.array([[]]) # 1-D array
+            m.matrixCheck()
+        with self.assertRaises(AttributeError):
+            m.B = np.array([[], [], []]) # 3-D array
             m.matrixCheck()
+        # @C 1x2
         with self.assertRaises(AttributeError):
-            m.F = np.array([[0], [1], [2]]) # 3-D array
+            m.C = np.array([[]]) # 0-D array
             m.matrixCheck()
         with self.assertRaises(AttributeError):
-            m.F = np.array([[0,0], [-9.81, -1]]) # extra column values per row
+            m.C = np.array([[0, 0], [1, 1]]) # 2-D array
             m.matrixCheck()
+        # @D 1x1
         with self.assertRaises(AttributeError):
-            m.F = np.array([[], []]) # less column values per row
+            m.D = np.array([]) # 0-D array
+            m.matrixCheck()
+        with self.assertRaises(AttributeError):
+            m.D = np.array([[0], [1]]) # 2-D array
+            m.matrixCheck()
+        # E 2x1
+        with self.assertRaises(AttributeError):
+            m.E = np.array([[0]]) # 1-D array
+            m.matrixCheck()
+        with self.assertRaises(AttributeError):
+            m.E = np.array([[0], [1], [2]]) # 3-D array
+            m.matrixCheck()
+        
+        # when matrix is improperly shaped
+        # @A 2x2
+        with self.assertRaises(AttributeError):
+            m.A = np.array([[0, 1, 2, 3], [0, 0, 1, 2]]) # extra column values per row
+            m.matrixCheck()
+        with self.assertRaises(AttributeError):
+            m.A = np.array([[0], [0]]) # less column values per row
             m.matrixCheck()
         with self.assertRaises(AttributeError): 
-            m.F = np.array([[0, 1], [0, 0], [2, 2]]) # extra row
+            m.A = np.array([[0, 1], [0, 0], [2, 2]]) # extra row
             m.matrixCheck()
+        with self.assertRaises(AttributeError): 
+            m.A = np.array([[0, 1]]) # less row
+            m.matrixCheck()
+        # @B 2x0
         with self.assertRaises(AttributeError):
-            # if F is none, we need to override event_state
-            m_noes = FNoneNoEventStateLM()
-        m.F = np.array([[0, 1]]) # less row
-        m.matrixCheck()
-        m.F = None
-        m.matrixCheck()
-
-        #G
-        with self.assertRaises(TypeError):
-            m.G = "[[0, 1], [0, 0]]" # string
-            m.matrixCheck() 
-        with self.assertRaises(TypeError):
-            m.G = 0 # int
+            m.B = np.array([[0, 1 ,2]]) # extra column values per row
             m.matrixCheck()
-        with self.assertRaises(TypeError):
-            m.G = 3.14 # float
+        with self.assertRaises(AttributeError):
+            m.B = np.array([[0]]) # less column values per row
             m.matrixCheck()
-        with self.assertRaises(TypeError):
-            m.G = {} # dict
-            m.matrixCheck() 
-        with self.assertRaises(TypeError):
-            m.G = () # tuple
-            m.matrixCheck() 
-        with self.assertRaises(TypeError):
-            m.G = set() # set
-            m.matrixCheck() 
-        with self.assertRaises(TypeError):
-            m.G = True # boolean
+        with self.assertRaises(AttributeError): 
+            m.B = np.array([[0, 1], [1, 1], [2, 2]]) # extra row
+            m.matrixCheck()
+        with self.assertRaises(AttributeError): 
+            m.B = np.array([[0, 1]]) # less row
             m.matrixCheck()
+        # @C 1x2
         with self.assertRaises(AttributeError):
-            m.G = np.array([0]) # 1-D Array
+            m.C = np.array([[1, 0, 2]]) # extra column values per row
             m.matrixCheck()
         with self.assertRaises(AttributeError):
-            m.G = np.array([[[0], [1], [2]]]) # 3-D array
+            m.C = np.array([[1]]) # less column values per row
             m.matrixCheck()
+        with self.assertRaises(AttributeError): 
+            m.C = np.array([[0, 0], [1, 1], [2, 2]]) # extra row
+            m.matrixCheck()
+        with self.assertRaises(AttributeError): 
+            m.C = np.array([[]]) # less row
+            m.matrixCheck()
+        # @D 1x1
         with self.assertRaises(AttributeError):
-            m.G = np.array([[0,0], [-9.81, -1]]) # extra column values per row
+            m.D = np.array([[1, 2]]) # extra column values per row
             m.matrixCheck()
         with self.assertRaises(AttributeError):
-            m.G = np.array([[], []]) # less column values per row
+            m.D = np.array([[]]) # less column values per row
+            m.matrixCheck()
+        with self.assertRaises(AttributeError): 
+            m.D = np.array([[0], [1]]) # extra row
+            m.matrixCheck()
+        with self.assertRaises(AttributeError): 
+            m.D = np.array([[]]) # less row
+            m.matrixCheck()
+        # @E 2x1
+        with self.assertRaises(TypeError):
+            m.E = np.array([0,0], [-9.81, -1]) # extra column values per row
+            m.matrixCheck()
+        with self.assertRaises(AttributeError):
+            m.E = np.array([[], []]) # less column values per row
             m.matrixCheck()
         with self.assertRaises(AttributeError): 
-            m.G = np.array([[0, 1], [0, 0], [2, 2]]) # extra row
+            m.E = np.array([[0, 1], [0, 0], [2, 2]]) # extra row
             m.matrixCheck()
         with self.assertRaises(AttributeError): 
-            m.G = np.array([[0, 1]]) # less row
+            m.E = np.array([[0, 1]]) # less row
             m.matrixCheck()
+        # @G 1x1
         with self.assertRaises(AttributeError):
             m.G = np.array([0, 1]) # extra column values per row
             m.matrixCheck()
@@ -341,93 +233,6 @@ def test_linear_model(self):
         with self.assertRaises(AttributeError): 
             m.G = np.array([[]]) # less row
             m.matrixCheck()
-        m.G = np.array([[0]]) # 1-D Array
-        m.matrixCheck()
-        m.G = None # sets to Default Value
-        m.matrixCheck()
-
-        # Error Demonstration
-        mTest = LinearThrownObjectFourStates()
-        with self.assertRaises(AttributeError):
-            mTest.B = np.array([[0], [1], [2], [3]])
-
-    @unittest.skip
-    def test_copy_linear(self):
-        m1 = LinearThrownObject()
-        copym1 = copy.copy(m1)  
-
-        self.assertTrue(m1 == copym1) # Testing Copy for Linear Model
-        deepcopym1 = copy.deepcopy(m1)
-        # Checking to see if all the copies are equal to one another before making changes to one or the other
-        self.assertTrue(m1 == deepcopym1)
-        self.assertTrue(copym1 == deepcopym1)
-
-
-        m2 = LinearThrownObject()
-        copym2 = copy.copy(m2)
-        self.assertTrue(m1 == m2)
-        self.assertTrue(m1 == copym2)
-        self.assertTrue(m2 == copym1)
-        self.assertTrue(copym2 == copym1)
-
-        m3 = LinearThrownObjectNoE()
-        m4 = LinearThrownDiffThrowingSpeed()
-
-        copym3 = copy.copy(m3)
-        copym4 = copy.copy(m4)
-
-        deepcopym3 = copy.deepcopy(m3)
-        deepcopym4 = copy.deepcopy(m4)
-
-        self.assertTrue(deepcopym3 == copym3)
-        self.assertFalse(copym4 == copym3)
-        
-        self.assertFalse(deepcopym4 == deepcopym3)
-
-
-        m5 = LinearThrownObjectUpdatedInitalizedMethod()
-        copym5 = copy.copy(m5)
-        deepcopym5 = copy.deepcopy(m5)
-
-        self.assertTrue(m5 == copym5 == deepcopym5)
-        
-        m5.states.append('C')
-        copym5.states.append('D')
-        deepcopym5.states.append('E')
-        # This test should be failing, but it is passing
-        self.assertFalse(copym5 == deepcopym5)
-
-        
-    def test_linear_pickle(self):
-        # future tests can compress, transfer to a file, and see if it still works
-        m1 = LinearThrownObject()
-        m2 = LinearThrownDiffThrowingSpeed()
-
-        # Note: dumps = serializing;
-        #       loads = deserializing
-
-        bytes_m1 = pickle.dumps(m1) #serializing object 
-        loaded_m1 = pickle.loads(bytes_m1) #deserializing the object
-        self.assertTrue(m1 == loaded_m1) # see if serializing and deserializing changes original form
-
-        bytes_m2 = pickle.dumps(m2)
-        loaded_m2 = pickle.loads(bytes_m2)
-        self.assertTrue(m2 == loaded_m2)
-
-        m3 = LinearThrownObject()
-        bytes_m3 = pickle.dumps(m3)
-        loaded_m3 = pickle.loads(bytes_m3)
-        self.assertTrue(m3 == loaded_m3)
-
-        self.assertTrue(bytes_m1 == bytes_m3)
-        self.assertTrue(loaded_m3 == loaded_m1)
-        self.assertTrue(LinearThrownObject, type(loaded_m3))
-
-        l = LinearThrownObjectUpdatedInitalizedMethod()
-        bytes_l = pickle.dumps(l)
-        loaded_l = pickle.loads(bytes_l)
-        self.assertTrue(l == loaded_l)
-        self.assertFalse(bytes_l == bytes_m1)
 
     def test_F_property_not_none(self):
         class ThrownObject(LinearThrownObject):
@@ -475,13 +280,14 @@ def threshold_met(self, x):
                 }
         # Needs more development; test coverage needs testing of event_state not overridden
 
+# This allows the module to be executed directly
 def run_tests():
     unittest.main()
     
 def main():
     l = unittest.TestLoader()
     runner = unittest.TextTestRunner()
-    print("\n\nTesting Linear Models")
+    print("\n\nTesting Base Models")
     result = runner.run(l.loadTestsFromTestCase(TestLinearModel)).wasSuccessful()
 
     if not result:
diff --git a/tests/test_serialization.py b/tests/test_serialization.py
index e470fd040..5fc8e9c94 100644
--- a/tests/test_serialization.py
+++ b/tests/test_serialization.py
@@ -65,7 +65,11 @@ def future_loading_2(t, x=None):
         new_model = DMDModel.from_json(save_json_dict)
 
         # Check serialization
-        self.assertEqual(surrogate_orig.parameters, new_model.parameters)
+        for key in surrogate_orig.parameters.keys():
+            if key != 'dmd_matrix':
+                self.assertEqual(surrogate_orig.parameters[key], new_model.parameters[key]) 
+            else: 
+                self.assertEqual((surrogate_orig.parameters['dmd_matrix']==new_model.parameters['dmd_matrix']).all(), True)
         
         # Check deserialization
         options_sim = {
diff --git a/tests/test_sim_result.py b/tests/test_sim_result.py
index 31a5e3631..0d77fb6f3 100644
--- a/tests/test_sim_result.py
+++ b/tests/test_sim_result.py
@@ -5,6 +5,7 @@
 import pickle
 import sys
 import unittest
+import pandas as pd
 
 from prog_models.models import BatteryElectroChemEOD
 from prog_models.sim_result import SimResult, LazySimResult
@@ -12,161 +13,212 @@
 
 
 class TestSimResult(unittest.TestCase):
-    def setUp(self):
-        # set stdout (so it wont print)
+    """def setUp(self):
+        set stdout (so it won't print)
         sys.stdout = StringIO()
 
     def tearDown(self):
-        sys.stdout = sys.__stdout__
-    
+        sys.stdout = sys.__stdout__"""
+
     def test_sim_result(self):
-        NUM_ELEMENTS = 5
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+        # Variables
+        time = list(range(5))
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(5)]
+        frame = pd.DataFrame(state)
+        frame.insert(0, "time", time)
         result = SimResult(time, state)
+        # Checks values from SimResult object and static variables
         self.assertListEqual(list(result), state)
         self.assertListEqual(result.times, time)
+        self.assertTrue(frame.equals(result.frame))
         for i in range(5):
-            self.assertEqual(result.time(i), time[i])
+            self.assertEqual(result.get_time(i), time[i])
             self.assertEqual(result[i], state[i])
-
+        self.assertTrue(frame.equals(result.frame))
         try:
-            tmp = result[NUM_ELEMENTS]
+            tmp = result[5]
             self.fail("Should be out of range error")
         except IndexError:
             pass
-
         try:
-            tmp = result.time(NUM_ELEMENTS)
+            tmp = result.times[5]
             self.fail("Should be out of range error")
         except IndexError:
             pass
-    
+
     def test_pickle(self):
-        NUM_ELEMENTS = 5
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+        # Variables
+        time = list(range(5))  # list of int, 0 to 4
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(5)]
         result = SimResult(time, state)
         pickle.dump(result, open('model_test.pkl', 'wb'))
         result2 = pickle.load(open('model_test.pkl', 'rb'))
         self.assertEqual(result, result2)
+        self.assertTrue(result.frame.equals(result2.frame))
 
     def test_extend(self):
-        NUM_ELEMENTS = 5  # Creating two result objects
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 2.5} for i in range(NUM_ELEMENTS)]
+        # Variables
+        time = list(range(5))  # list of int from 0 to 4
+        state = [{'a': i * 2.5, 'b': i * 2.5} for i in range(5)]
         result = SimResult(time, state)
-        NUM_ELEMENTS = 10
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
-        result2 = SimResult(time, state)
-        self.assertEqual(result.times, [0, 1, 2, 3, 4])
-        self.assertEqual(result2.times, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0.0}, {'a': 2.5, 'b': 2.5}, {'a': 5.0, 'b': 5.0}, {'a': 7.5, 'b': 7.5}, {'a': 10.0, 'b': 10.0}]) # Assert data is correct before extending
-        self.assertEqual(result2.data, [{'a': 0, 'b': 0}, {'a': 5, 'b': 5}, {'a': 10, 'b': 10}, {'a': 15, 'b': 15}, {'a': 20, 'b': 20}, {'a': 25, 'b': 25}, {'a': 30, 'b': 30}, {'a': 35, 'b': 35}, {'a': 40, 'b': 40}, {'a': 45, 'b': 45}])
-        
+        time2 = list(range(10))  # list of int from 0 to 9
+        state2 = [{'a': i * 5, 'b': i * 5} for i in range(10)]
+        result2 = SimResult(time2, state2)
+        time_extended = time + time2
+        state_extended = state + state2
+        result_df = pd.DataFrame(state)
+        result_df.insert(0, "time", time)
+        result2_df = pd.DataFrame(state2)
+        result2_df.insert(0, "time", time2)
+        result_extend_df = pd.concat([result_df, result2_df], axis=0).reset_index(drop=True)
+
+        self.assertEqual(result.times, time)
+        self.assertEqual(result2.times, time2)
+        self.assertEqual(result.data, state)  # Assert data is correct before extending
+        self.assertEqual(result2.data, state2)
+
         result.extend(result2)  # Extend result with result2
-        self.assertEqual(result.times, [0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0.0}, {'a': 2.5, 'b': 2.5}, {'a': 5.0, 'b': 5.0}, {'a': 7.5, 'b': 7.5}, {'a': 10.0, 'b': 10.0}, {'a': 0, 'b': 0}, {'a': 5, 'b': 5}, {'a': 10, 'b': 10}, {'a': 15, 'b': 15}, {'a': 20, 'b': 20}, {'a': 25, 'b': 25}, {'a': 30, 'b': 30}, {'a': 35, 'b': 35}, {'a': 40, 'b': 40}, {'a': 45, 'b': 45}])
+        self.assertEqual(result.times, time_extended)
+        self.assertEqual(result.data, state_extended)
+        self.assertTrue(result.frame.equals(result_extend_df))
 
         self.assertRaises(ValueError, result.extend, 0)  # Passing non-LazySimResult types to extend method
-        self.assertRaises(ValueError, result.extend, [0,1])
+        self.assertRaises(ValueError, result.extend, [0, 1])
         self.assertRaises(ValueError, result.extend, {})
         self.assertRaises(ValueError, result.extend, set())
         self.assertRaises(ValueError, result.extend, 1.5)
 
     def test_extended_by_lazy(self):
-        NUM_ELEMENTS = 5 
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 2.5} for i in range(NUM_ELEMENTS)]
-        result = SimResult(time, state) # Creating one SimResult object
-        def f(x):
-            return {k:v * 2 for k,v in x.items()}
-        NUM_ELEMENTS = 10
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
-        result2 = LazySimResult(f, time, state) # Creating one LazySimResult object
+        # Variables
+        time = list(range(5))  # list of int, 0 to 4
+        state = [{'a': i * 2.5, 'b': i * 2.5} for i in range(5)]
+        time2 = list(range(10))  # list of int, 0 to 9
+        state2 = [{'a': i * 5, 'b': i * 5} for i in range(10)]
+        data2 = [{'a': i * 10, 'b': i * 10} for i in range(10)]
+        result_df = pd.DataFrame(state)
+        result_df.insert(0, "time", time)
+        result2_df = pd.DataFrame(state2)
+        result2_df.insert(0, "time", time2)
+        result_extend_df = pd.concat([result_df, result2_df], axis=0).reset_index(drop=True)
+        result = SimResult(time, state)  # Creating one SimResult object
 
-        self.assertEqual(result.times, [0, 1, 2, 3, 4])
-        self.assertEqual(result2.times, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0.0}, {'a': 2.5, 'b': 2.5}, {'a': 5.0, 'b': 5.0}, {'a': 7.5, 'b': 7.5}, {'a': 10.0, 'b': 10.0}]) # Assert data is correct before extending
-        self.assertEqual(result2.data, [{'a': 0, 'b': 0}, {'a': 10, 'b': 10}, {'a': 20, 'b': 20}, {'a': 30, 'b': 30}, {'a': 40, 'b': 40}, {'a': 50, 'b': 50}, {'a': 60, 'b': 60}, {'a': 70, 'b': 70}, {'a': 80, 'b': 80}, {'a': 90, 'b': 90}])
-        result.extend(result2) # Extend result with result2
-        self.assertEqual(result.times, [0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0.0}, {'a': 2.5, 'b': 2.5}, {'a': 5.0, 'b': 5.0}, {'a': 7.5, 'b': 7.5}, {'a': 10.0, 'b': 10.0}, {'a': 0, 'b': 0}, {'a': 10, 'b': 10}, {'a': 20, 'b': 20}, {'a': 30, 'b': 30}, {'a': 40, 'b': 40}, {'a': 50, 'b': 50}, {'a': 60, 'b': 60}, {'a': 70, 'b': 70}, {'a': 80, 'b': 80}, {'a': 90, 'b': 90}])
+        def f(x):
+            return {k: v * 2 for k, v in x.items()}
+
+        result2 = LazySimResult(f, time2, state2)  # Creating one LazySimResult object
+        # confirming the data in result and result2 are correct
+        self.assertEqual(result.times, time)
+        self.assertEqual(result2.times, time2)
+        self.assertEqual(result.data, state)  # Assert data is correct before extending
+        self.assertEqual(result2.data, data2)
+        result.extend(result2)  # Extend result with result2
+        # check data when result is extended with result2
+        self.assertEqual(result.times, time + time2)
+        self.assertEqual(result.data, state + data2)
+        self.assertTrue(result.frame.equals(result_extend_df))
 
     def test_pickle_lazy(self):
         def f(x):
-            return {k:v * 2 for k,v in x.items()}
-        NUM_ELEMENTS = 5
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 2.5} for i in range(NUM_ELEMENTS)]
+            return {k: v * 2 for k, v in x.items()}
+
+        # Variables
+        time = list(range(5))  # list of int, 0 to 4
+        state = [{'a': i * 2.5, 'b': i * 2.5} for i in range(5)]
         lazy_result = LazySimResult(f, time, state)  # Ordinary LazySimResult with f, time, state
         sim_result = SimResult(time, state)  # Ordinary SimResult with time,state
 
         converted_lazy_result = SimResult(lazy_result.times, lazy_result.data)
+        self.assertFalse(sim_result.frame.equals(converted_lazy_result))
         self.assertNotEqual(sim_result, converted_lazy_result)  # converted is not the same as the original SimResult
 
         pickle.dump(lazy_result, open('model_test.pkl', 'wb'))
         pickle_converted_result = pickle.load(open('model_test.pkl', 'rb'))
         self.assertEqual(converted_lazy_result, pickle_converted_result)
-    
+        self.assertTrue(converted_lazy_result.frame.equals(pickle_converted_result.frame))
+
     def test_index(self):
-        NUM_ELEMENTS = 5 # Creating two result objects
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+        # Variables
+        time = list(range(5))  # list of int, 0 to 4
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(5)]
         result = SimResult(time, state)
 
         self.assertEqual(result.index({'a': 10, 'b': 20}), 4)
         self.assertEqual(result.index({'a': 2.5, 'b': 5}), 1)
         self.assertEqual(result.index({'a': 0, 'b': 0}), 0)
-        self.assertRaises(ValueError, result.index, 6.0) # Other argument doesn't exist
-        self.assertRaises(ValueError, result.index, -1) # Non-existent data value
-        self.assertRaises(ValueError, result.index, "7.5") # Data specified incorrectly as string
-        self.assertRaises(ValueError, result.index, None) # Not type errors because its simply looking for an object in list
+        self.assertRaises(ValueError, result.index, 6.0)  # Other argument doesn't exist
+        self.assertRaises(ValueError, result.index, -1)  # Non-existent data value
+        self.assertRaises(ValueError, result.index, "7.5")  # Data specified incorrectly as string
+        self.assertRaises(ValueError, result.index,
+                          None)  # Not type errors because its simply looking for an object in list
         self.assertRaises(ValueError, result.index, [1, 2])
         self.assertRaises(ValueError, result.index, {})
         self.assertRaises(ValueError, result.index, set())
 
     def test_pop(self):
-        NUM_ELEMENTS = 5
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+        # Variables
+        time = list(map(float, range(5)))
+        state = [{'a': i * 2.5, 'b': i * 5.0} for i in range(5)]
         result = SimResult(time, state)
-
-        result.pop(2) # Test specified index
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0}, {'a': 2.5, 'b': 5}, {'a': 7.5, 'b': 15}, {'a': 10.0, 'b': 20}])
-        result.pop() # Test default index -1 (last element)
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0}, {'a': 2.5, 'b': 5}, {'a': 7.5, 'b': 15}])
-        result.pop(-1) # Test argument of index -1 (last element)
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0}, {'a': 2.5, 'b': 5}])
-        result.pop(0) # Test argument of 0
-        self.assertEqual(result.data, [{'a': 2.5, 'b': 5}])
-        self.assertRaises(IndexError, result.pop, 5) # Test specifying an invalid index value
+        result_df = pd.DataFrame(state)
+        result_df.insert(0, "time", time)
+
+        result.pop(2)  # Test specified index
+        state.remove({'a': 5.0, 'b': 10})  # update state by removing value
+        self.assertEqual(result.data, state)
+        # removing row from DataFrame
+        temp_df = result_df.T
+        temp_df.pop(2)
+        result_df = temp_df.T.reset_index(drop=True)
+        self.assertTrue(result.frame.equals(result_df))
+        result.pop()  # Test default index -1 (last element)
+        state.pop()  # pop state, removes last item
+        # removing row from DataFrame
+        temp_df = result_df.T
+        temp_df.pop(3)
+        result_df = temp_df.T.reset_index(drop=True)
+        self.assertEqual(result.data, state)
+        self.assertTrue(result.frame.equals(result_df))
+        result.pop(-1)  # Test argument of index -1 (last element)
+        state.pop()  # pop state, removes last item
+        # removing row from DataFrame
+        temp_df = result_df.T
+        temp_df.pop(2)
+        result_df = temp_df.T.reset_index(drop=True)
+        self.assertEqual(result.data, state)
+        self.assertTrue(result.frame.equals(result_df))
+        result.pop(0)  # Test argument of 0
+        state.pop(0)  # pop state, removes first item
+        # removing row from DataFrame
+        temp_df = result_df.T
+        temp_df.pop(0)
+        result_df = temp_df.T.reset_index(drop=True)
+        self.assertTrue(result.frame.equals(result_df))
+        self.assertEqual(result.data, state)
+        self.assertRaises(IndexError, result.pop, 5)  # Test specifying an invalid index value
         self.assertRaises(IndexError, result.pop, 3)
-        self.assertRaises(TypeError, result.pop, "5") # Test specifying an invalid index type
-        self.assertRaises(TypeError, result.pop, [0,1])
+        self.assertRaises(TypeError, result.pop, "5")  # Test specifying an invalid index type
+        self.assertRaises(TypeError, result.pop, [0, 1])
         self.assertRaises(TypeError, result.pop, {})
         self.assertRaises(TypeError, result.pop, set())
         self.assertRaises(TypeError, result.pop, 1.5)
 
     def test_to_numpy(self):
-        NUM_ELEMENTS = 10
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+        # Variables
+        time = list(range(10))  # list of int, 0 to 9
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(10)]
         result = SimResult(time, state)
         np_result = result.to_numpy()
         self.assertIsInstance(np_result, np.ndarray)
-        self.assertEqual(np_result.shape, (NUM_ELEMENTS, 2))
+        self.assertEqual(np_result.shape, (10, 2))
         self.assertEqual(np_result.dtype, np.dtype('float64'))
-        self.assertTrue(np.all(np_result==np.array([[i * 2.5, i * 5] for i in range(NUM_ELEMENTS)])))
+        self.assertTrue(np.all(np_result == np.array([[i * 2.5, i * 5] for i in range(10)])))
 
         # Subset of keys
         result = result.to_numpy(['b'])
         self.assertIsInstance(result, np.ndarray)
-        self.assertEqual(result.shape, (NUM_ELEMENTS, 1))
+        self.assertEqual(result.shape, (10, 1))
         self.assertEqual(result.dtype, np.dtype('float64'))
-        self.assertTrue(np.all(result==np.array([[i * 5] for i in range(NUM_ELEMENTS)])))
+        self.assertTrue(np.all(result == np.array([[i * 5] for i in range(10)])))
 
         # Now test when empty
         result = SimResult([], [])
@@ -180,127 +232,182 @@ def test_to_numpy(self):
         result = SimResult(time, state)
         result = result.to_numpy()
         self.assertIsInstance(result, np.ndarray)
-        self.assertEqual(result.shape, (NUM_ELEMENTS, 2))
+        self.assertEqual(result.shape, (10, 2))
         self.assertEqual(result.dtype, np.dtype('float64'))
-        self.assertTrue(np.all(result==np.array([[i * 2.5, i * 5] for i in range(NUM_ELEMENTS)])))
+        self.assertTrue(np.all(result == np.array([[i * 2.5, i * 5] for i in range(10)])))
 
     def test_remove(self):
-        NUM_ELEMENTS = 5 # Creating two result objects
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+        # Variables
+        time = list(range(5))  # list of int, 0 to 4
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(5)]
         result = SimResult(time, state)
-
-        result.remove( {'a': 5.0, 'b': 10}) # Positional defaults to removing data
-        self.assertEqual(result.times, [0, 1, 3, 4])
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0}, {'a': 2.5, 'b': 5}, {'a': 7.5, 'b': 15}, {'a': 10.0, 'b': 20}])
-        result.remove(d = {'a': 0.0, 'b': 0}) # Testing named removal of data
-        self.assertEqual(result.times, [1, 3, 4])
-        self.assertEqual(result.data, [{'a': 2.5, 'b': 5}, {'a': 7.5, 'b': 15}, {'a': 10.0, 'b': 20}])
-        result.remove(t = 3) # Testing named removal of time
-        self.assertEqual(result.times, [1, 4])
-        self.assertEqual(result.data, [{'a': 2.5, 'b': 5}, {'a': 10.0, 'b': 20}])
-        result.remove(t = 1) 
-        self.assertEqual(result.times, [4])
-        self.assertEqual(result.data, [{'a': 10.0, 'b': 20}])
-
-        self.assertRaises(ValueError, result.remove, ) # If nothing specified, raise ValueError
-        self.assertRaises(ValueError, result.remove, None, None) # Passing both as None
-        self.assertRaises(ValueError, result.remove, 0.0, 1) # Passing arguments to both
-        self.assertRaises(ValueError, result.remove, 7.5) # Test nonexistent data value
-        self.assertRaises(ValueError, result.remove, -1) # Type checking negated as index searches for element in list
-        self.assertRaises(ValueError, result.remove, "5") # Thus all value types allowed to be searched
-        self.assertRaises(ValueError, result.remove, [0,1])
+        result_df = pd.DataFrame(state)
+        result_df.insert(0, "time", time)
+
+        result.remove({'a': 5.0, 'b': 10})  # Positional defaults to removing data
+        # Update Variables
+        time.remove(2)
+        state.remove({'a': 5.0, 'b': 10})
+        temp_df = result_df.T
+        temp_df.pop(2)
+        result_df = temp_df.T.reset_index(drop=True)
+        self.assertTrue(result_df.equals(result.frame))
+        self.assertEqual(result.times, time)
+        self.assertEqual(result.data, state)
+        result.remove(d={'a': 0.0, 'b': 0})  # Testing named removal of data
+        # Update Variables
+        time.remove(0)
+        state.remove({'a': 0.0, 'b': 0})
+        temp_df = result_df.T
+        temp_df.pop(0)
+        result_df = temp_df.T.reset_index(drop=True)
+        self.assertTrue(result_df.equals(result.frame))
+        self.assertEqual(result.times, time)
+        self.assertEqual(result.data, state)
+        result.remove(t=3)  # Testing named removal of time
+        # Update Variables
+        time.remove(3)
+        state.remove({'a': 7.5, 'b': 15})
+        temp_df = result_df.T
+        temp_df.pop(1)
+        result_df = temp_df.T.reset_index(drop=True)
+        self.assertTrue(result_df.equals(result.frame))
+        self.assertEqual(result.times, time)
+        self.assertEqual(result.data, state)
+        result.remove(t=1)
+        # Update Variables
+        time.remove(1)
+        state.remove({'a': 2.5, 'b': 5})
+        temp_df = result_df.T
+        temp_df.pop(0)
+        result_df = temp_df.T.reset_index(drop=True)
+        self.assertTrue(result_df.equals(result.frame))
+        self.assertEqual(result.times, time)
+        self.assertEqual(result.data, state)
+
+        self.assertRaises(ValueError, result.remove, )  # If nothing specified, raise ValueError
+        self.assertRaises(ValueError, result.remove, None, None)  # Passing both as None
+        self.assertRaises(ValueError, result.remove, 0.0, 1)  # Passing arguments to both
+        self.assertRaises(ValueError, result.remove, 7.5)  # Test nonexistent data value
+        self.assertRaises(ValueError, result.remove, -1)  # Type checking negated as index searches for element in list
+        self.assertRaises(ValueError, result.remove, "5")  # Thus all value types allowed to be searched
+        self.assertRaises(ValueError, result.remove, [0, 1])
         self.assertRaises(ValueError, result.remove, {})
         self.assertRaises(ValueError, result.remove, set())
 
     def test_clear(self):
-        NUM_ELEMENTS = 5 # Creating two result objects
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+        # Variables
+        time = list(range(5))  # list of int, 0 to 4
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(5)]
         result = SimResult(time, state)
-        self.assertEqual(result.times, [0, 1, 2, 3, 4])
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0}, {'a': 2.5, 'b': 5}, {'a': 5, 'b': 10}, {'a': 7.5, 'b': 15}, {'a': 10.0, 'b': 20}])
+        result_df = pd.DataFrame(state)
+        result_df.insert(0, "time", time)
+        self.assertEqual(result.times, time)
+        self.assertEqual(result.data, state)
+        self.assertTrue(result.frame.equals(result_df))
         self.assertRaises(TypeError, result.clear, True)
 
         result.clear()
+        self.assertTrue(result.frame.equals(pd.DataFrame()))
         self.assertEqual(result.times, [])
         self.assertEqual(result.data, [])
 
-    def test_time(self):
-        NUM_ELEMENTS = 5 # Creating two result objects
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+    def test_get_time(self):
+        # Variables
+        # Creating two result objects
+        time = list(range(5))  # list of int, 0 to 4
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(5)]
         result = SimResult(time, state)
-        self.assertEqual(result.time(0), result.times[0])
-        self.assertEqual(result.time(1), result.times[1])
-        self.assertEqual(result.time(2), result.times[2])
-        self.assertEqual(result.time(3), result.times[3])
-        self.assertEqual(result.time(4), result.times[4])
-
-        self.assertRaises(TypeError, result.time, ) # Test no input given
-        self.assertRaises(TypeError, result.time, "0") # Tests specifying an invalid index type 
-        self.assertRaises(TypeError, result.time, [0,1])
-        self.assertRaises(TypeError, result.time, {})
-        self.assertRaises(TypeError, result.time, set())
-        self.assertRaises(TypeError, result.time, 1.5)
+        self.assertEqual(result.get_time(0), result.times[0])
+        self.assertEqual(result.get_time(1), result.times[1])
+        self.assertEqual(result.get_time(2), result.times[2])
+        self.assertEqual(result.get_time(3), result.times[3])
+        self.assertEqual(result.get_time(4), result.times[4])
+
+        self.assertRaises(TypeError, result.get_time, )  # Test no input given
+        self.assertRaises(TypeError, result.get_time, "0")  # Tests specifying an invalid index type
+        self.assertRaises(TypeError, result.get_time, [0, 1])
+        self.assertRaises(TypeError, result.get_time, {})
+        self.assertRaises(TypeError, result.get_time, set())
+        self.assertRaises(TypeError, result.get_time, 1.5)
 
     def test_plot(self):
         # Testing model taken from events.py
         YELLOW_THRESH, RED_THRESH, THRESHOLD = 0.15, 0.1, 0.05
+
         class MyBatt(BatteryElectroChemEOD):
             events = BatteryElectroChemEOD.events + ['EOD_warn_yellow', 'EOD_warn_red', 'EOD_requirement_threshold']
+
             def event_state(self, state):
                 event_state = super().event_state(state)
-                event_state['EOD_warn_yellow'] = (event_state['EOD']-YELLOW_THRESH)/(1-YELLOW_THRESH) 
-                event_state['EOD_warn_red'] = (event_state['EOD']-RED_THRESH)/(1-RED_THRESH)
-                event_state['EOD_requirement_threshold'] = (event_state['EOD']-THRESHOLD)/(1-THRESHOLD)
+                event_state['EOD_warn_yellow'] = (event_state['EOD'] - YELLOW_THRESH) / (1 - YELLOW_THRESH)
+                event_state['EOD_warn_red'] = (event_state['EOD'] - RED_THRESH) / (1 - RED_THRESH)
+                event_state['EOD_requirement_threshold'] = (event_state['EOD'] - THRESHOLD) / (1 - THRESHOLD)
                 return event_state
+
             def threshold_met(self, x):
-                t_met =  super().threshold_met(x)
+                t_met = super().threshold_met(x)
                 event_state = self.event_state(x)
                 t_met['EOD_warn_yellow'] = event_state['EOD_warn_yellow'] <= 0
                 t_met['EOD_warn_red'] = event_state['EOD_warn_red'] <= 0
                 t_met['EOD_requirement_threshold'] = event_state['EOD_requirement_threshold'] <= 0
                 return t_met
+
         def future_loading(t, x=None):
-            if (t < 600): i = 2
-            elif (t < 900): i = 1
-            elif (t < 1800): i = 4
-            elif (t < 3000): i = 2     
-            else: i = 3
-            return {'i': i} 
+            if (t < 600):
+                i = 2
+            elif (t < 900):
+                i = 1
+            elif (t < 1800):
+                i = 4
+            elif (t < 3000):
+                i = 2
+            else:
+                i = 3
+            return {'i': i}
+
         m = MyBatt()
-        (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(future_loading, threshold_keys=['EOD'], print = False)
-        plot_test = event_states.plot() # Plot doesn't raise error
+        (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(future_loading, threshold_keys=['EOD'],
+                                                                                 print=False)
+        plot_test = event_states.plot()  # Plot doesn't raise error
 
     def test_namedtuple_access(self):
         # Testing model taken from events.py
         YELLOW_THRESH, RED_THRESH, THRESHOLD = 0.15, 0.1, 0.05
-        
+
         class MyBatt(BatteryElectroChemEOD):
             events = BatteryElectroChemEOD.events + ['EOD_warn_yellow', 'EOD_warn_red', 'EOD_requirement_threshold']
+
             def event_state(self, state):
                 event_state = super().event_state(state)
-                event_state['EOD_warn_yellow'] = (event_state['EOD']-YELLOW_THRESH)/(1-YELLOW_THRESH) 
-                event_state['EOD_warn_red'] = (event_state['EOD']-RED_THRESH)/(1-RED_THRESH)
-                event_state['EOD_requirement_threshold'] = (event_state['EOD']-THRESHOLD)/(1-THRESHOLD)
+                event_state['EOD_warn_yellow'] = (event_state['EOD'] - YELLOW_THRESH) / (1 - YELLOW_THRESH)
+                event_state['EOD_warn_red'] = (event_state['EOD'] - RED_THRESH) / (1 - RED_THRESH)
+                event_state['EOD_requirement_threshold'] = (event_state['EOD'] - THRESHOLD) / (1 - THRESHOLD)
                 return event_state
+
             def threshold_met(self, x):
-                t_met =  super().threshold_met(x)
+                t_met = super().threshold_met(x)
                 event_state = self.event_state(x)
                 t_met['EOD_warn_yellow'] = event_state['EOD_warn_yellow'] <= 0
                 t_met['EOD_warn_red'] = event_state['EOD_warn_red'] <= 0
                 t_met['EOD_requirement_threshold'] = event_state['EOD_requirement_threshold'] <= 0
                 return t_met
+
         def future_loading(t, x=None):
-            if (t < 600): i = 2
-            elif (t < 900): i = 1
-            elif (t < 1800): i = 4
-            elif (t < 3000): i = 2     
-            else: i = 3
-            return {'i': i} 
+            if (t < 600):
+                i = 2
+            elif (t < 900):
+                i = 1
+            elif (t < 1800):
+                i = 4
+            elif (t < 3000):
+                i = 2
+            else:
+                i = 3
+            return {'i': i}
+
         m = MyBatt()
-        named_results = m.simulate_to_threshold(future_loading, threshold_keys=['EOD'], print = False)
+        named_results = m.simulate_to_threshold(future_loading, threshold_keys=['EOD'], print=False)
         times = named_results.times
         inputs = named_results.inputs
         states = named_results.states
@@ -321,154 +428,261 @@ def test_not_implemented(self):
     # Tests for LazySimResult
     def test_lazy_data_fcn(self):
         def f(x):
-            return {k:v * 2 for k,v in x.items()}
-        NUM_ELEMENTS = 5
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+            return {k: v * 2 for k, v in x.items()}
+
+        # Variables
+        time = list(range(5))
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(5)]
+        state2 = [{'a': i * 5.0, 'b': i * 10} for i in range(5)]
         result = LazySimResult(f, time, state)
+        result_df = pd.DataFrame(state)
+        result_df.insert(0, "time", time)
+        result_df2 = pd.DataFrame(state2)
+        result_df2.insert(0, "time", time)
+
         self.assertFalse(result.is_cached())
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0}, {'a': 5.0, 'b': 10}, {'a': 10.0, 'b': 20}, {'a': 15.0, 'b': 30}, {'a': 20.0, 'b': 40}])
+        self.assertEqual(result.data, state2)
+        self.assertTrue(result_df.equals(result.frame))
+        self.assertTrue(result_df2.equals(result.get_frame_data()))
         self.assertTrue(result.is_cached())
 
     def test_lazy_clear(self):
         def f(x):
-            return {k:v * 2 for k,v in x.items()}
-        NUM_ELEMENTS = 5
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+            return {k: v * 2 for k, v in x.items()}
+
+        # Variables
+        time = list(range(5))  # list of int, 0 to 4
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(5)]
+        state2 = [{'a': i * 5.0, 'b': i * 10} for i in range(5)]
+        result_df = pd.DataFrame(state)
+        result_df.insert(0, "time", time)
+        result_df2 = pd.DataFrame(state2)
+        result_df2.insert(0, "time", time)
         result = LazySimResult(f, time, state)
-        self.assertEqual(result.times, [0, 1, 2, 3, 4])
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0}, {'a': 5.0, 'b': 10}, {'a': 10.0, 'b': 20}, {'a': 15.0, 'b': 30}, {'a': 20.0, 'b': 40}])
-        self.assertEqual(result.states, [{'a': 0.0, 'b': 0}, {'a': 2.5, 'b': 5}, {'a': 5.0, 'b': 10}, {'a': 7.5, 'b': 15}, {'a': 10.0, 'b': 20}])
+        self.assertEqual(result.times, time)
+        self.assertEqual(result.data, state2)
+        self.assertEqual(result.states, state)
+        self.assertTrue(result_df.equals(result.frame))
+        self.assertTrue(result_df2.equals(result.get_frame_data()))
         self.assertRaises(TypeError, result.clear, True)
 
         result.clear()
         self.assertEqual(result.times, [])
         self.assertEqual(result.data, [])
         self.assertEqual(result.states, [])
+        self.assertTrue(result.frame.equals(pd.DataFrame()))
+        self.assertTrue(result.get_frame_data().equals(pd.DataFrame()))
 
     def test_lazy_extend(self):
         def f(x):
-            return {k:v * 2 for k,v in x.items()}
-        NUM_ELEMENTS = 5
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+            return {k: v * 2 for k, v in x.items()}
+
+        # Variables
+        time = list(range(5))  # list of int, 0 to 4
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(5)]
         result = LazySimResult(f, time, state)
+        time2 = list(range(10))  # list of int, 0 to 9
+        state2 = [{'a': i * 5, 'b': i * 10} for i in range(10)]
+        data2 = [{'a': i * 25, 'b': i * 50} for i in range(10)]
+        data = [{'a': i * 5.0, 'b': i * 10} for i in range(5)]
+        result_df = pd.DataFrame(state)
+        result_df.insert(0, "time", time)
+        result_df2 = pd.DataFrame(state2)
+        result_df2.insert(0, "time", time2)
+        result_data_df = pd.DataFrame(data)
+        result_data_df.insert(0, "time", time)
+        result_data_df2 = pd.DataFrame(data2)
+        result_data_df2.insert(0, "time", time2)
+        result_df_extended = pd.DataFrame(state+state2)
+        result_df_extended.insert(0, "time", time+time2)
+        result_data_df_extended = pd.DataFrame(data+data2)
+        result_data_df_extended.insert(0, "time", time+time2)
 
         def f2(x):
-            return {k:v * 5 for k,v in x.items()}
-        NUM_ELEMENTS = 10
-        time2 = list(range(NUM_ELEMENTS))
-        state2 = [{'a': i * 5, 'b': i * 10} for i in range(NUM_ELEMENTS)]
+            return {k: v * 5 for k, v in x.items()}
+
         result2 = LazySimResult(f2, time2, state2)
-        self.assertEqual(result.times, [0, 1, 2, 3, 4]) # Assert data is correct before extending
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0}, {'a': 5.0, 'b': 10}, {'a': 10.0, 'b': 20}, {'a': 15.0, 'b': 30}, {'a': 20.0, 'b': 40}])
-        self.assertEqual(result.states, [{'a': 0.0, 'b': 0}, {'a': 2.5, 'b': 5}, {'a': 5.0, 'b': 10}, {'a': 7.5, 'b': 15}, {'a': 10.0, 'b': 20}])
-        self.assertEqual(result2.times, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
-        self.assertEqual(result2.data, [{'a': 0, 'b': 0}, {'a': 25, 'b': 50}, {'a': 50, 'b': 100}, {'a': 75, 'b': 150}, {'a': 100, 'b': 200}, {'a': 125, 'b': 250}, {'a': 150, 'b': 300}, {'a': 175, 'b': 350}, {'a': 200, 'b': 400}, {'a': 225, 'b': 450}])
-        self.assertEqual(result2.states, [{'a': 0, 'b': 0}, {'a': 5, 'b': 10}, {'a': 10, 'b': 20}, {'a': 15, 'b': 30}, {'a': 20, 'b': 40}, {'a': 25, 'b': 50}, {'a': 30, 'b': 60}, {'a': 35, 'b': 70}, {'a': 40, 'b': 80}, {'a': 45, 'b': 90}])
+        self.assertEqual(result.times, time)  # Assert data is correct before extending
+        self.assertEqual(result.data, data)
+        self.assertEqual(result.states, state)
+        self.assertEqual(result2.times, time2)
+        self.assertEqual(result2.data, data2)
+        self.assertEqual(result2.states, state2)
+        self.assertTrue(result.frame.equals(result_df))
+        self.assertTrue(result.get_frame_data().equals(result_data_df))
+        self.assertTrue(result2.frame.equals(result_df2))
+        self.assertTrue(result2.get_frame_data().equals(result_data_df2))
 
         result.extend(result2)
-        self.assertEqual(result.times, [0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) # Assert data is correct after extending
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0}, {'a': 5.0, 'b': 10}, {'a': 10.0, 'b': 20}, {'a': 15.0, 'b': 30}, {'a': 20.0, 'b': 40}, {'a': 0, 'b': 0}, {'a': 25, 'b': 50}, {'a': 50, 'b': 100}, {'a': 75, 'b': 150}, {'a': 100, 'b': 200}, {'a': 125, 'b': 250}, {'a': 150, 'b': 300}, {'a': 175, 'b': 350}, {'a': 200, 'b': 400}, {'a': 225, 'b': 450}])
-        self.assertEqual(result.states, [{'a': 0.0, 'b': 0}, {'a': 2.5, 'b': 5}, {'a': 5.0, 'b': 10}, {'a': 7.5, 'b': 15}, {'a': 10.0, 'b': 20}, {'a': 0, 'b': 0}, {'a': 5, 'b': 10}, {'a': 10, 'b': 20}, {'a': 15, 'b': 30}, {'a': 20, 'b': 40}, {'a': 25, 'b': 50}, {'a': 30, 'b': 60}, {'a': 35, 'b': 70}, {'a': 40, 'b': 80}, {'a': 45, 'b': 90}])
+        self.assertEqual(result.times, time + time2)  # Assert data is correct after extending
+        self.assertEqual(result.data, data + data2)
+        self.assertEqual(result.states, state + state2)
+        self.assertTrue(result.frame.equals(result_df_extended))
+        self.assertTrue(result.get_frame_data().equals(result_data_df_extended))
 
     def test_lazy_extend_cache(self):
         def f(x):
-            return {k:v * 2 for k,v in x.items()}
-        NUM_ELEMENTS = 5
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+            return {k: v * 2 for k, v in x.items()}
+
+        # Variables
+        time = list(range(5))
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(5)]
+        data = [{'a': i * 5.0, 'b': i * 10} for i in range(5)]
         result1 = LazySimResult(f, time, state)
         result2 = LazySimResult(f, time, state)
+        result_df = pd.DataFrame(state)
+        result_df.insert(0, "time", time)
+        result_data_df = pd.DataFrame(data)
+        result_data_df.insert(0, "time", time)
+        result_df_ext = pd.DataFrame(state+state)
+        result_df_ext.insert(0, "time", time+time)
+        result_data_df_ext = pd.DataFrame(data+data)
+        result_data_df_ext.insert(0, "time", time+time)
 
         # Case 1
         result1.extend(result2)
-        self.assertFalse(result1.is_cached()) # False
-        
+        self.assertFalse(result1.is_cached())  # False
+
         # Case 2
-        result1 = LazySimResult(f, time, state) # Reset result1
-        store_test_data = result1.data # Access result1 data
-        result1.extend(result2) 
-        self.assertFalse(result1.is_cached()) # False
+        result1 = LazySimResult(f, time, state)  # Reset result1
+        self.assertTrue(result1.frame.equals(result_df))
+        store_test_data = result1.data  # Access result1 data
+        result1.extend(result2)
+        self.assertTrue(result1.frame.equals(result_df_ext))
+        self.assertFalse(result1.is_cached())  # False
 
         # Case 3
-        result1 = LazySimResult(f, time, state) # Reset result1
-        store_test_data = result2.data # Access result2 data
-        result1.extend(result2) 
-        self.assertFalse(result1.is_cached()) # False
+        result1 = LazySimResult(f, time, state)  # Reset result1
+        store_test_data = result2.data  # Access result2 data
+        result1.extend(result2)
+        self.assertFalse(result1.is_cached())  # False
 
         # Case 4
-        result1 = LazySimResult(f, time, state) # Reset result1
-        result2 = LazySimResult(f, time, state) # Reset result2
-        store_test_data1 = result1.data # Access result1 data
-        store_test_data2 = result2.data # Access result2 data
-        result1.extend(result2) 
-        self.assertTrue(result1.is_cached()) # True
+        result1 = LazySimResult(f, time, state)  # Reset result1
+        result2 = LazySimResult(f, time, state)  # Reset result2
+        store_test_data1 = result1.data  # Access result1 data
+        store_test_data2 = result2.data  # Access result2 data
+        result1.extend(result2)
+        self.assertTrue(result1.is_cached())  # True
+        self.assertTrue(result1.get_frame_data().equals(result_data_df_ext))
 
     def test_lazy_extend_error(self):
         def f(x):
-            return {k:v * 2 for k,v in x.items()}
-        NUM_ELEMENTS = 5
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+            return {k: v * 2 for k, v in x.items()}
+
+        # Variables
+        time = list(range(5))  # list of int, - to 4
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(5)]
         result = LazySimResult(f, time, state)
         sim_result = SimResult(time, state)
 
-        self.assertRaises(ValueError, result.extend, sim_result) # Passing a SimResult to LazySimResult's extend
-        self.assertRaises(ValueError, result.extend, 0) # Passing non-LazySimResult types to extend method
-        self.assertRaises(ValueError, result.extend, [0,1])
+        self.assertRaises(ValueError, result.extend, sim_result)  # Passing a SimResult to LazySimResult's extend
+        self.assertRaises(ValueError, result.extend, 0)  # Passing non-LazySimResult types to extend method
+        self.assertRaises(ValueError, result.extend, [0, 1])
         self.assertRaises(ValueError, result.extend, {})
         self.assertRaises(ValueError, result.extend, set())
         self.assertRaises(ValueError, result.extend, 1.5)
 
     def test_lazy_pop(self):
         def f(x):
-            return {k:v * 2 for k,v in x.items()}
-        NUM_ELEMENTS = 5
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
-        result = LazySimResult(f, time, state)
+            return {k: v * 2 for k, v in x.items()}
 
-        result.pop(1) # Test specified index
-        self.assertEqual(result.times, [0, 2, 3, 4])
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0}, {'a': 10.0, 'b': 20}, {'a': 15.0, 'b': 30}, {'a': 20.0, 'b': 40}])
-        self.assertEqual(result.states, [{'a': 0.0, 'b': 0}, {'a': 5.0, 'b': 10}, {'a': 7.5, 'b': 15}, {'a': 10.0, 'b': 20}])
-
-        result.pop() # Test default index -1 (last element)
-        self.assertEqual(result.times, [0, 2, 3])
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0}, {'a': 10.0, 'b': 20}, {'a': 15.0, 'b': 30}])
-        self.assertEqual(result.states, [{'a': 0.0, 'b': 0}, {'a': 5.0, 'b': 10}, {'a': 7.5, 'b': 15}])
-
-        result.pop(-1) # Test argument of index -1 (last element)
-        self.assertEqual(result.times, [0, 2])
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0}, {'a': 10.0, 'b': 20}])
-        self.assertEqual(result.states, [{'a': 0.0, 'b': 0}, {'a': 5.0, 'b': 10}])
-        result.pop(0) # Test argument of 0
-        self.assertEqual(result.times, [2])
-        self.assertEqual(result.data, [{'a': 10.0, 'b': 20}])
-        self.assertEqual(result.states, [{'a': 5.0, 'b': 10}])
+        # Variables
+        time = list(map(float, range(5)))  # list of int, 0 to 4
+        state = [{'a': i * 2.5, 'b': i * 5.0} for i in range(5)]
+        data = [{'a': i * 5.0, 'b': i * 10.0} for i in range(5)]
+        result = LazySimResult(f, time, state)
+        result_df = pd.DataFrame(state)
+        result_df.insert(0, "time", time)
+        result_data_df = pd.DataFrame(data)
+        result_data_df.insert(0, "time", time)
+
+        result.pop(1)  # Test specified index
+        time.remove(1)  # remove value '1' to check time values after pop
+        # removing row from DataFrame
+        temp_df = result_df.T
+        temp_df.pop(1)
+        result_df = temp_df.T.reset_index(drop=True)
+        self.assertTrue(result.frame.equals(result_df))
+        self.assertEqual(result.times, time)
+        # removing row from fcn DataFrame
+        temp_df = result_data_df.T
+        temp_df.pop(1)
+        result_data_df = temp_df.T.reset_index(drop=True)
+        self.assertTrue(result.get_frame_data().equals(result_data_df))
+        data.remove({'a': 5.0, 'b': 10})  # removes index 1 value from data list
+        self.assertEqual(result.data, data)
+        state.remove({'a': 2.5, 'b': 5})  # removes index 1 value from state list
+        self.assertEqual(result.states, state)
+        # removing row from fcn DataFrame
+        self.assertTrue(result.get_frame_data().equals(result_data_df))
+
+        result.pop()  # Test default index -1 (last element)
+        time.pop()
+        data.pop()
+        state.pop()
+        # removing row from fcn DataFrame
+        temp_df = result_data_df.T
+        temp_df.pop(3)
+        result_data_df = temp_df.T.reset_index(drop=True)
+        self.assertEqual(result.times, time)
+        self.assertEqual(result.data, data)
+        self.assertEqual(result.states, state)
+        self.assertTrue(result.get_frame_data().equals(result_data_df))
+
+        result.pop(-1)  # Test argument of index -1 (last element)
+        time.pop(-1)
+        data.pop(-1)
+        state.pop(-1)
+        # removing row from fcn DataFrame
+        temp_df = result_data_df.T
+        temp_df.pop(2)
+        result_data_df = temp_df.T.reset_index(drop=True)
+        self.assertEqual(result.times, time)
+        self.assertEqual(result.data, data)
+        self.assertEqual(result.states, state)
+        self.assertTrue(result.get_frame_data().equals(result_data_df))
+        result.pop(0)  # Test argument of 0
+        time.pop(0)
+        data.pop(0)
+        state.pop(0)
+        # removing row from fcn DataFrame
+        temp_df = result_data_df.T
+        temp_df.pop(0)
+        result_data_df = temp_df.T.reset_index(drop=True)
+        self.assertTrue(result.get_frame_data().equals(result_data_df))
+        self.assertEqual(result.times, time)
+        self.assertEqual(result.data, data)
+        self.assertEqual(result.states, state)
         # Test erroneous input
-        self.assertRaises(IndexError, result.pop, 5) # Test specifying an invalid index value
+        self.assertRaises(IndexError, result.pop, 5)  # Test specifying an invalid index value
         self.assertRaises(IndexError, result.pop, 3)
-        self.assertRaises(TypeError, result.pop, "5") # Test specifying an invalid index type
-        self.assertRaises(TypeError, result.pop, [0,1])
+        self.assertRaises(TypeError, result.pop, "5")  # Test specifying an invalid index type
+        self.assertRaises(TypeError, result.pop, [0, 1])
         self.assertRaises(TypeError, result.pop, {})
         self.assertRaises(TypeError, result.pop, set())
         self.assertRaises(TypeError, result.pop, 1.5)
 
     def test_cached_sim_result(self):
         def f(x):
-            return {k:v * 2 for k,v in x.items()}
+            return {k: v * 2 for k, v in x.items()}
+
         NUM_ELEMENTS = 5
         time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+        state = [{'a': i * 2.5, 'b': i * 5.0} for i in range(NUM_ELEMENTS)]
+        data = [{'a': i * 5.0, 'b': i * 10.0} for i in range(NUM_ELEMENTS)]
+        result_df = pd.DataFrame(state)
+        result_df.insert(0, "time", time)
+        result_data_df = pd.DataFrame(data)
+        result_data_df.insert(0, "time", time)
         result = LazySimResult(f, time, state)
+        result_df_ext = pd.DataFrame(state+state)
+        result_df_ext.insert(0, "time", time+time)
         self.assertFalse(result.is_cached())
         self.assertListEqual(result.times, time)
         for i in range(5):
-            self.assertEqual(result.time(i), time[i])
-            self.assertEqual(result[i], {k:v*2 for k,v in state[i].items()})
+            self.assertEqual(result.get_time(i), time[i])
+            self.assertEqual(result[i], {k: v * 2 for k, v in state[i].items()})
         self.assertTrue(result.is_cached())
 
         try:
@@ -478,7 +692,7 @@ def f(x):
             pass
 
         try:
-            tmp = result.time(NUM_ELEMENTS)
+            tmp = result.get_time(NUM_ELEMENTS)
             self.fail("Should be out of range error")
         except IndexError:
             pass
@@ -488,60 +702,105 @@ def f(x):
         result = LazySimResult(f, time, state)
         result2 = LazySimResult(f, time, state)
         self.assertTrue(result == result2)
+        self.assertTrue(result.frame.equals(result2.frame))
         self.assertEqual(len(result), len(result2))
         result.extend(LazySimResult(f, time, state))
         self.assertFalse(result == result2)
         self.assertNotEqual(len(result), len(result2))
+        self.assertTrue(result.frame.equals(result_df_ext))
 
     def test_lazy_remove(self):
         def f(x):
-            return {k:v * 2 for k,v in x.items()}
-        NUM_ELEMENTS = 10
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
-        result = LazySimResult(f, time, state)
+            return {k: v * 2 for k, v in x.items()}
 
-        result.remove({'a': 5.0, 'b': 10}) # Unnamed default positional argument removal of data value
-        self.assertEqual(result.times, [0, 2, 3, 4, 5, 6, 7, 8, 9])
-        self.assertEqual(result.data, [{'a': 0.0, 'b': 0}, {'a': 10.0, 'b': 20}, {'a': 15.0, 'b': 30}, {'a': 20.0, 'b': 40}, {'a': 25.0, 'b': 50}, {'a': 30.0, 'b': 60}, {'a': 35.0, 'b': 70}, {'a': 40.0, 'b': 80}, {'a': 45.0, 'b': 90}])
-        self.assertEqual(result.states, [{'a': 0.0, 'b': 0}, {'a': 5.0, 'b': 10}, {'a': 7.5, 'b': 15}, {'a': 10.0, 'b': 20}, {'a': 12.5, 'b': 25}, {'a': 15.0, 'b': 30}, {'a': 17.5, 'b': 35}, {'a': 20.0, 'b': 40}, {'a': 22.5, 'b': 45}])
-        result.remove(d = {'a': 0.0, 'b': 0}) # Named argument removal of data value
-        self.assertEqual(result.times, [2, 3, 4, 5, 6, 7, 8, 9])
-        self.assertEqual(result.data, [{'a': 10.0, 'b': 20}, {'a': 15.0, 'b': 30}, {'a': 20.0, 'b': 40}, {'a': 25.0, 'b': 50}, {'a': 30.0, 'b': 60}, {'a': 35.0, 'b': 70}, {'a': 40.0, 'b': 80}, {'a': 45.0, 'b': 90}])
-        self.assertEqual(result.states, [{'a': 5.0, 'b': 10}, {'a': 7.5, 'b': 15}, {'a': 10.0, 'b': 20}, {'a': 12.5, 'b': 25}, {'a': 15.0, 'b': 30}, {'a': 17.5, 'b': 35}, {'a': 20.0, 'b': 40}, {'a': 22.5, 'b': 45}])
-        result.remove(t = 7) # Named argument removal of times value
-        self.assertEqual(result.times, [2, 3, 4, 5, 6, 8, 9])
-        self.assertEqual(result.data, [{'a': 10.0, 'b': 20}, {'a': 15.0, 'b': 30}, {'a': 20.0, 'b': 40}, {'a': 25.0, 'b': 50}, {'a': 30.0, 'b': 60}, {'a': 40.0, 'b': 80}, {'a': 45.0, 'b': 90}])
-        self.assertEqual(result.states, [{'a': 5.0, 'b': 10}, {'a': 7.5, 'b': 15}, {'a': 10.0, 'b': 20}, {'a': 12.5, 'b': 25}, {'a': 15.0, 'b': 30}, {'a': 20.0, 'b': 40}, {'a': 22.5, 'b': 45}])
-        result.remove(s = {'a': 12.5, 'b': 25}) # Named argument removal of states value
-        self.assertEqual(result.times, [2, 3, 4, 6, 8, 9])
-        self.assertEqual(result.data, [{'a': 10.0, 'b': 20}, {'a': 15.0, 'b': 30}, {'a': 20.0, 'b': 40}, {'a': 30.0, 'b': 60}, {'a': 40.0, 'b': 80}, {'a': 45.0, 'b': 90}])
-        self.assertEqual(result.states, [{'a': 5.0, 'b': 10}, {'a': 7.5, 'b': 15}, {'a': 10.0, 'b': 20}, {'a': 15.0, 'b': 30}, {'a': 20.0, 'b': 40}, {'a': 22.5, 'b': 45}])
-
-        self.assertRaises(ValueError, result.remove, ) # Test no values specified
-        self.assertRaises(ValueError, result.remove, 90.0, 2) # Test two values specified positionally
-        self.assertRaises(ValueError, result.remove, 90.0, 2, 15.0) # Test three values specified positionally
-        self.assertRaises(ValueError, result.remove, d=90.0, t=2) # Test d,t values specified by name
-        self.assertRaises(ValueError, result.remove, t=2, s=15.0) # Test s,t values specified by name
-        self.assertRaises(ValueError, result.remove, d=90.0, s=15.0) # Test d,s values specified by name
-        self.assertRaises(ValueError, result.remove, d=90.0, t=2, s=15.0) # Test three values specified by name
-        self.assertRaises(ValueError, result.remove, 90.0) # Test nonexistent data value
-        self.assertRaises(ValueError, result.remove, d=90.0) # Test nonexistent data value
-        self.assertRaises(ValueError, result.remove, t=90.0) # Test nonexistent times value
-        self.assertRaises(ValueError, result.remove, s=90.0) # Test nonexistent states value
-        self.assertRaises(ValueError, result.remove, -1) # Type checking negated as index searches for element in list
-        self.assertRaises(ValueError, result.remove, "5") # Thus all value types allowed to be searched
-        self.assertRaises(ValueError, result.remove, [0,1])
+        # Variables
+        time = list(range(10))  # list of int, 0 to 9
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(10)]
+        result = LazySimResult(f, time, state)
+        data = [{'a': i * 5.0, 'b': i * 10} for i in range(10)]
+        result_df = pd.DataFrame(state)
+        result_df.insert(0, "time", time)
+        result_data_df = pd.DataFrame(state)
+        result_data_df.insert(0, "time", time)
+
+        result.remove({'a': 5.0, 'b': 10})  # Unnamed default positional argument removal of data value
+        # Update Variables
+        state.remove({'a': 2.5, 'b': 5})
+        time.remove(1)
+        data.remove({'a': 5.0, 'b': 10})
+        # removing row from DataFrame
+        temp_df = result_df.T
+        temp_df.pop(1)
+        result_df = temp_df.T.reset_index(drop=True)
+        self.assertTrue(result.frame.equals(result_df))
+        self.assertEqual(result.times, time)
+        self.assertEqual(result.data, data)
+        self.assertEqual(result.states, state)
+        result.remove(d={'a': 0.0, 'b': 0})  # Named argument removal of data value
+        # Update Variables
+        state.remove({'a': 0.0, 'b': 0})
+        time.remove(0)
+        data.remove({'a': 0.0, 'b': 0})
+        # removing row from DataFrame
+        temp_df = result_df.T
+        temp_df.pop(0)
+        result_df = temp_df.T.reset_index(drop=True)
+        self.assertTrue(result.frame.equals(result_df))
+        self.assertEqual(result.times, time)
+        self.assertEqual(result.data, data)
+        self.assertEqual(result.states, state)
+        result.remove(t=7)  # Named argument removal of times value
+        # Update Variables
+        state.remove({'a': 17.5, 'b': 35})
+        time.remove(7)
+        data.remove({'a': 35.0, 'b': 70})
+        # removing row from DataFrame
+        temp_df = result_df.T
+        temp_df.pop(5)
+        result_df = temp_df.T.reset_index(drop=True)
+        self.assertTrue(result.frame.equals(result_df))
+        self.assertEqual(result.times, time)
+        self.assertEqual(result.data, data)
+        self.assertEqual(result.states, state)
+        result.remove(s={'a': 12.5, 'b': 25})  # Named argument removal of states value
+        # Update Variables
+        state.remove({'a': 12.5, 'b': 25})
+        time.remove(5)
+        data.remove({'a': 25, 'b': 50})
+        # removing row from DataFrame
+        temp_df = result_df.T
+        temp_df.pop(3)
+        result_df = temp_df.T.reset_index(drop=True)
+        self.assertTrue(result.frame.equals(result_df))
+        self.assertEqual(result.times, time)
+        self.assertEqual(result.data, data)
+        self.assertEqual(result.states, state)
+
+        self.assertRaises(ValueError, result.remove, )  # Test no values specified
+        self.assertRaises(ValueError, result.remove, 90.0, 2)  # Test two values specified positionally
+        self.assertRaises(ValueError, result.remove, 90.0, 2, 15.0)  # Test three values specified positionally
+        self.assertRaises(ValueError, result.remove, d=90.0, t=2)  # Test d,t values specified by name
+        self.assertRaises(ValueError, result.remove, t=2, s=15.0)  # Test s,t values specified by name
+        self.assertRaises(ValueError, result.remove, d=90.0, s=15.0)  # Test d,s values specified by name
+        self.assertRaises(ValueError, result.remove, d=90.0, t=2, s=15.0)  # Test three values specified by name
+        self.assertRaises(ValueError, result.remove, 90.0)  # Test nonexistent data value
+        self.assertRaises(ValueError, result.remove, d=90.0)  # Test nonexistent data value
+        self.assertRaises(ValueError, result.remove, t=90.0)  # Test nonexistent times value
+        self.assertRaises(ValueError, result.remove, s=90.0)  # Test nonexistent states value
+        self.assertRaises(ValueError, result.remove, -1)  # Type checking negated as index searches for element in list
+        self.assertRaises(ValueError, result.remove, "5")  # Thus all value types allowed to be searched
+        self.assertRaises(ValueError, result.remove, [0, 1])
         self.assertRaises(ValueError, result.remove, {})
         self.assertRaises(ValueError, result.remove, set())
 
     def test_lazy_not_implemented(self):
         # Not implemented functions, should raise errors
         def f(x):
-            return {k:v * 2 for k,v in x.items()}
-        NUM_ELEMENTS = 5
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+            return {k: v * 2 for k, v in x.items()}
+
+        # Variables
+        time = list(range(5))  # list of int, 0 to 4
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(5)]
         result = LazySimResult(f, time, state)
         self.assertRaises(NotImplementedError, result.append)
         self.assertRaises(NotImplementedError, result.count)
@@ -550,42 +809,54 @@ def f(x):
 
     def test_lazy_to_simresult(self):
         def f(x):
-            return {k:v * 2 for k,v in x.items()}
-        NUM_ELEMENTS = 5
-        time = list(range(NUM_ELEMENTS))
-        state = [{'a': i * 2.5, 'b': i * 5} for i in range(NUM_ELEMENTS)]
+            return {k: v * 2 for k, v in x.items()}
+
+        # Variables
+        time = list(map(float, range(5)))  # list of int, 0 to 4
+        state = [{'a': i * 2.5, 'b': i * 5} for i in range(5)]
+        data = [{'a': i * 5.0, 'b': i * 10} for i in range(5)]
+        result_df = pd.DataFrame(state)
+        result_df.insert(0, "time", time)
+        result_data_df = pd.DataFrame(data)
+        result_data_df.insert(0, "time", time)
         result = LazySimResult(f, time, state)
 
         converted_result = result.to_simresult()
-        self.assertTrue(isinstance(converted_result, SimResult)) # Ensure type is SimResult
-        self.assertEqual(converted_result.times, result.times) # Compare to original LazySimResult
+        self.assertTrue(isinstance(converted_result, SimResult))  # Ensure type is SimResult
+        self.assertEqual(converted_result.times, result.times)  # Compare to original LazySimResult
         self.assertEqual(converted_result.data, result.data)
-        self.assertEqual(converted_result.times, [0, 1, 2, 3, 4]) # Compare to expected values
-        self.assertEqual(converted_result.data, [{'a': 0.0, 'b': 0}, {'a': 5.0, 'b': 10}, {'a': 10.0, 'b': 20}, {'a': 15.0, 'b': 30}, {'a': 20.0, 'b': 40}])
+        self.assertEqual(converted_result.times, time)  # Compare to expected values
+        self.assertEqual(converted_result.data, data)
+        self.assertTrue(converted_result.frame.equals(result.get_frame_data()))
 
     def test_monotonicity(self):
-        NUM_ELEMENTS = 5
-        time = list(range(NUM_ELEMENTS))
+        # Variables
+        time = list(range(5))
 
         # Test monotonically increasing, decreasing
-        states = [{'a': 1+i/10, 'b': 2-i/5} for i in range(NUM_ELEMENTS)]
+        states = [{'a': 1 + i / 10, 'b': 2 - i / 5} for i in range(5)]
         result = SimResult(time, states)
         self.assertDictEqual(result.monotonicity(), {'a': 1.0, 'b': 1.0})
+        self.assertTrue(result.monotonicity_df().equals(pd.DataFrame({'a': 1.0, 'b': 1.0}, index=[0])))
 
         # Test monotonicity between range [0,1]
-        states = [{'a': i*(i%3-1), 'b': i*(i%3-1)} for i in range(NUM_ELEMENTS)]
+        states = [{'a': i * (i % 3 - 1), 'b': i * (i % 3 - 1)} for i in range(5)]
         result = SimResult(time, states)
         self.assertDictEqual(result.monotonicity(), {'a': 0.25, 'b': 0.25})
+        self.assertTrue(result.monotonicity_df().equals(pd.DataFrame({'a': 0.25, 'b': 0.25}, index=[0])))
 
         # # Test no monotonicity
-        states = [{'a': i*(i%2), 'b': i*(i%2)} for i in range(NUM_ELEMENTS)]
+        states = [{'a': i * (i % 2), 'b': i * (i % 2)} for i in range(5)]
         result = SimResult(time, states)
         self.assertDictEqual(result.monotonicity(), {'a': 0.0, 'b': 0.0})
-        
+        self.assertTrue(result.monotonicity_df().equals(pd.DataFrame(pd.DataFrame({'a': 0.0, 'b': 0.0}, index=[0]))))
+
+
 # This allows the module to be executed directly
 def run_tests():
     unittest.main()
-    
+
+
 def main():
     l = unittest.TestLoader()
     runner = unittest.TextTestRunner()
@@ -593,7 +864,8 @@ def main():
     result = runner.run(l.loadTestsFromTestCase(TestSimResult)).wasSuccessful()
 
     if not result:
-         raise Exception("Failed test")
+        raise Exception("Failed test")
+
 
 if __name__ == '__main__':
     main()