Repository: IN-CORE/pyincore (language: python)

Issue:
changed number of cpu to 4
I haven't edited changelog because it was the missing piece that was supposed to be in the previous PR.
Clean up Salt Lake City CGE formatting
Spacing is set to 2 instead of 4.
Building economic loss analysis breaks when no occupancy multiplier is provided
The buildingeconloss analysis breaks when no occupancy multiplier is provided. The [documentation](https://incore.ncsa.illinois.edu/doc/incore/analyses/building_loss.html) says that occupancy_multiplier is not required for this analysis


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7415e3f50..c5c508ca6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,27 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](http://keepachangelog.com/)
 and this project adheres to [Semantic Versioning](http://semver.org/).
 
+## [1.9.0] - 2023-03-15
+
+### Added
+- Method in space service to add a dataset by space name [#273](https://github.com/IN-CORE/pyincore/issues/273)
+- Method in space service to get space id by space name [#272](https://github.com/IN-CORE/pyincore/issues/272)
+- Method in space service to remove dataset from the space [#283](https://github.com/IN-CORE/pyincore/issues/283)
+- Method in space service to remove dataset by space name [#284](https://github.com/IN-CORE/pyincore/issues/284)
+- Combined wind, surge-wave, and flood building loss [#276](https://github.com/IN-CORE/pyincore/issues/276)
+
+### Changed
+- Rewrote clustering utility function to use flexible archetype column [#247](https://github.com/IN-CORE/pyincore/issues/247)
+- Made documentation containter to use requirements instead of environemt [#257](https://github.com/IN-CORE/pyincore/issues/257)
+- Parallelized the HHRS analysis [#268](https://github.com/IN-CORE/pyincore/issues/268)
+- Updated Salt Lake City CGE [#281](https://github.com/IN-CORE/pyincore/issues/281)
+- Tested hurricane windfield test methods [#100](https://github.com/IN-CORE/incore-services/issues/100)
+- Updatd Salt Lake City CGE formatting and handled infeasible case  [#287](https://github.com/IN-CORE/pyincore/issues/287)
+
+### Fixed
+- Duplicate input spec for housing recovery sequential model [#263](https://github.com/IN-CORE/pyincore/issues/263)
+- Updated building economic loss analysis to handle case when no occupancy multiplier is provided [#274](https://github.com/IN-CORE/pyincore/issues/274)
+
 ## [1.8.0] - 2022-11-16
 
 ### Changed
diff --git a/Dockerfile b/Dockerfile
index 9445a2668..284adc242 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -7,11 +7,12 @@ USER root
 
 # install packages
 WORKDIR /src
-COPY environment.yml .
+COPY requirements.txt .
+
 ENV PATH "$MAMBA_ROOT_PREFIX/bin:$PATH"
 RUN micromamba install -y -n base -c conda-forge \
     sphinx sphinx_rtd_theme \
-    -f environment.yml
+    -f requirements.txt
 
 # copy code and generate documentation
 COPY . ./
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 83eff8a82..5bb583fdd 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -33,9 +33,9 @@
 author = ''
 
 # The short X.Y version
-version = '1.8'
+version = '1.9'
 # The full version, including alpha/beta/rc tags
-release = '1.8.0'
+release = '1.9.0'
 
 # -- General configuration ---------------------------------------------------
 
diff --git a/docs/source/modules.rst b/docs/source/modules.rst
index b73eab9f0..2b278a61d 100644
--- a/docs/source/modules.rst
+++ b/docs/source/modules.rst
@@ -42,9 +42,15 @@ analyses/capitalshocks
     :members:
 
 analyses/combinedwindwavesurgebuildingdamage
+============================================
 ..  autoclass:: pyincore.analyses.combinedwindwavesurgebuildingdamage.CombinedWindWaveSurgeBuildingDamage
     :members:
 
+analyses/combinedwindwavesurgebuildingloss
+============================================
+..  autoclass:: pyincore.analyses.combinedwindwavesurgebuildingloss.CombinedWindWaveSurgeBuildingLoss
+    :members:
+
 analyses/cumulativebuildingdamage
 =================================
 ..  autoclass:: cumulativebuildingdamage.cumulativebuildingdamage.CumulativeBuildingDamage
diff --git a/pyincore/__init__.py b/pyincore/__init__.py
index d0e15722f..d525d4669 100644
--- a/pyincore/__init__.py
+++ b/pyincore/__init__.py
@@ -22,6 +22,7 @@
 from pyincore.utils.analysisutil import AnalysisUtil
 from pyincore.utils.popdisloutputprocess import PopDislOutputProcess
 from pyincore.utils.cgeoutputprocess import CGEOutputProcess
+from pyincore.utils.hhrsoutputprocess import HHRSOutputProcess
 from pyincore.dataset import Dataset, InventoryDataset, DamageRatioDataset
 from pyincore.models.fragilitycurveset import FragilityCurveSet
 from pyincore.models.repaircurveset import RepairCurveSet
diff --git a/pyincore/analyses/buildingeconloss/buildingeconloss.py b/pyincore/analyses/buildingeconloss/buildingeconloss.py
index 52338fbc6..75b795420 100755
--- a/pyincore/analyses/buildingeconloss/buildingeconloss.py
+++ b/pyincore/analyses/buildingeconloss/buildingeconloss.py
@@ -40,8 +40,12 @@ def run(self):
         bldg_set = self.get_input_dataset("buildings").get_inventory_reader()
 
         # Occupancy type of the exposure
-        occ_multiplier = self.get_input_dataset("occupancy_multiplier").get_csv_reader()
-        occ_mult_df = pd.DataFrame(occ_multiplier)
+        occ_multiplier = self.get_input_dataset("occupancy_multiplier")
+        if occ_multiplier is not None:
+            occ_multiplier = occ_multiplier.get_csv_reader()
+            occ_mult_df = pd.DataFrame(occ_multiplier)
+        else:
+            occ_mult_df = None
 
         try:
             prop_select = []
@@ -112,7 +116,7 @@ def add_multipliers(self, dmg_set_df, occ_mult_df):
             dmg_set_df = pd.merge(dmg_set_df, occ_mult_df, how="left", left_on="occ_type",
                                   right_on="occ_type", sort=True, copy=True)
         else:
-            dmg_set_df = dmg_set_df["Multiplier"] = 1.0
+            dmg_set_df["Multiplier"] = 1.0
 
         return dmg_set_df
 
diff --git a/pyincore/analyses/buildingfunctionality/buildingfunctionality.py b/pyincore/analyses/buildingfunctionality/buildingfunctionality.py
index 8d0a3c673..2ec438ce1 100644
--- a/pyincore/analyses/buildingfunctionality/buildingfunctionality.py
+++ b/pyincore/analyses/buildingfunctionality/buildingfunctionality.py
@@ -108,11 +108,9 @@ def run(self):
         else:
             poles_df = None
 
-        # All three above needs to be present at the same time or none
-        if not all(epf is None for epf in [interdependency_dict, substations_df, poles_df]) \
-                and not all(epf is not None for epf in [interdependency_dict, substations_df, poles_df]):
-            raise ValueError("To consider electric power availability in the analysis, please provide pole damage "
-                             "samples, substation damage samples, and interdependency.")
+        if (poles_dataset is not None or substations_dataset is not None) and interdependency_dataset is None:
+            raise ValueError("Please provide interdependency table if pole or substation damage is "
+                             "considered in the building functionality calculation.")
 
         functionality_probabilities = []
         functionality_samples = []
@@ -168,30 +166,53 @@ def functionality(self, building_guid, buildings, substations, poles, interdepen
         if interdependency is not None:
 
             if building_guid in interdependency.keys():
-                substations_mc_samples = substations.loc[interdependency[building_guid]["substations_guid"]]
-                poles_mc_samples = poles.loc[interdependency[building_guid]["poles_guid"]]
-
-                substation_list = []
-                try:
-                    substation_list = substations_mc_samples["failure"].split(",")
-                except IndexError:
-                    print("error with substations")
-                    print(interdependency[building_guid]["substations_guid"])
-                    return {building_guid: -1}
-
-                pole_list = []
-                try:
-                    pole_list = poles_mc_samples["failure"].split(",")
-                except IndexError:
-                    print("error with poles")
-                    print(interdependency[building_guid]["poles_guid"])
-                    return {building_guid: -1}
-
-                functionality_samples = [BuildingFunctionality._calc_functionality_samples(building_sample,
-                                                                                           substation_sample,
-                                                                                           pole_sample)
-                                         for building_sample, substation_sample, pole_sample in
-                                         zip(building_list, substation_list, pole_list)]
+                if substations is not None:
+                    substations_mc_samples = substations.loc[interdependency[building_guid]["substations_guid"]]
+                    substation_list = []
+                    try:
+                        substation_list = substations_mc_samples["failure"].split(",")
+                    except IndexError:
+                        print("error with substations")
+                        print(interdependency[building_guid]["substations_guid"])
+                        return {building_guid: -1}
+                else:
+                    substation_list = None
+
+                if poles is not None:
+                    poles_mc_samples = poles.loc[interdependency[building_guid]["poles_guid"]]
+                    pole_list = []
+                    try:
+                        pole_list = poles_mc_samples["failure"].split(",")
+                    except IndexError:
+                        print("error with poles")
+                        print(interdependency[building_guid]["poles_guid"])
+                        return {building_guid: -1}
+                else:
+                    pole_list = None
+
+                if substation_list is not None and pole_list is not None:
+                    functionality_samples = [BuildingFunctionality._calc_functionality_samples(building_sample,
+                                                                                               substation_sample,
+                                                                                               pole_sample)
+                                             for building_sample, substation_sample, pole_sample in
+                                             zip(building_list, substation_list, pole_list)]
+                elif substation_list is not None:
+                    functionality_samples = [BuildingFunctionality._calc_functionality_samples(building_sample,
+                                                                                               substation_sample,
+                                                                                               None)
+                                             for building_sample, substation_sample in
+                                             zip(building_list, substation_list)]
+                elif pole_list is not None:
+                    functionality_samples = [BuildingFunctionality._calc_functionality_samples(building_sample,
+                                                                                               None,
+                                                                                               pole_sample)
+                                             for building_sample, pole_sample in
+                                             zip(building_list, pole_list)]
+                else:
+                    functionality_samples = [BuildingFunctionality._calc_functionality_samples(building_sample,
+                                                                                               None,
+                                                                                               None)
+                                             for building_sample in building_list]
                 probability = BuildingFunctionality._calc_functionality_probability(functionality_samples)
                 return building_guid, ",".join([str(sample) for sample in functionality_samples]), probability
 
@@ -212,8 +233,8 @@ def _calc_functionality_samples(building_sample, substation_sample=None, pole_sa
 
         Args:
             building_sample (str): Monte Carlo samples of building functionality.
-            substation_sample (str): Monte Carlo samples of substation functionality.
-            pole_sample (str): Monte Carlo samples of pole functionality.
+            substation_sample (str|None): Monte Carlo samples of substation functionality.
+            pole_sample (str|None): Monte Carlo samples of pole functionality.
 
         Returns:
             int: 1 if building is functional, 0 otherwise
diff --git a/pyincore/analyses/combinedwindwavesurgebuildingdamage/combinedwindwavesurgebuildingdamage.py b/pyincore/analyses/combinedwindwavesurgebuildingdamage/combinedwindwavesurgebuildingdamage.py
index 0035a37c6..9b85cb78d 100755
--- a/pyincore/analyses/combinedwindwavesurgebuildingdamage/combinedwindwavesurgebuildingdamage.py
+++ b/pyincore/analyses/combinedwindwavesurgebuildingdamage/combinedwindwavesurgebuildingdamage.py
@@ -11,6 +11,7 @@
 
 class CombinedWindWaveSurgeBuildingDamage(BaseAnalysis):
     """ Determines overall building maximum damage state from wind, flood and surge-wave damage
+    and uses the maximum damage probabilities from the 3 damages to determine overall damage
     
     Args:
         incore_client (IncoreClient): Service authentication.
@@ -63,10 +64,70 @@ def run(self):
         combined_output['f_max_ds'] = combined_output['f_max_ds'].replace(old_ds_vals, new_ds_vals)
         combined_output['max_state'] = combined_output['max_state'].replace(old_ds_vals, new_ds_vals)
 
+        # Find combined damage
+        combined_bldg_dmg = self.get_combined_damage(wind_damage, surge_wave_damage, flood_damage)
+
+        # Create the result containing the 3 combined damages into a single damage
+        self.set_result_csv_data("ds_result", combined_bldg_dmg, self.get_parameter("result_name") + "_combined_dmg",
+                                 "dataframe")
+
         # Create the result dataset
         self.set_result_csv_data("result", combined_output, self.get_parameter("result_name"), "dataframe")
 
         return True
+    
+    def get_combined_damage(self, wind_dmg: pd.DataFrame, sw_dmg: pd.DataFrame, flood_dmg: pd.DataFrame):
+        """Calculates overall building damage
+                Determines the overall building damage probabilities from the 3 hazards by taking the maximum.
+
+                Args:
+                    wind_dmg (pd.DataFrame): Table of wind damage for the building inventory
+                    sw_dmg (pd.DataFrame): Table of surge-wave damage for the building inventory
+                    flood_dmg (pd.DataFrame): Table of flood damage for the building inventory
+
+                Returns:
+                    pd.DataFrame: An table of combined damage probabilities for the building inventory
+
+                """
+        flood_dmg.rename(columns={'LS_0': 'f_LS_0', 'LS_1': 'f_LS_1', 'LS_2': 'f_LS_2', 'DS_0': 'f_DS_0', 
+                                  'DS_1': 'f_DS_1', 'DS_2': 'f_DS_2', 'DS_3': 'f_DS_3', 'haz_expose': 'f_haz_expose'}, 
+                         inplace=True)
+
+        sw_dmg.rename(columns={'LS_0': 'sw_LS_0', 'LS_1': 'sw_LS_1', 'LS_2': 'sw_LS_2', 'DS_0': 'sw_DS_0',
+                               'DS_1': 'sw_DS_1', 'DS_2': 'sw_DS_2', 'DS_3': 'sw_DS_3', 'haz_expose': 
+                                   'sw_haz_expose'}, inplace=True)
+
+        wind_dmg.rename(columns={'LS_0': 'w_LS_0', 'LS_1': 'w_LS_1', 'LS_2': 'w_LS_2', 'DS_0': 'w_DS_0',
+                                 'DS_1': 'w_DS_1', 'DS_2': 'w_DS_2', 'DS_3': 'w_DS_3', 'haz_expose': 'w_haz_expose'},
+                        inplace=True)
+
+        combined_df = pd.merge(pd.merge(wind_dmg, sw_dmg, on='guid'), flood_dmg, on='guid')
+
+        def find_match(row, col_name):
+            max_finder = {
+                row['f_DS_3']: 'f_',
+                row['w_DS_3']: 'w_',
+                row['sw_DS_3']: 'sw_'
+            }
+
+            return row[max_finder[max(max_finder.keys())] + col_name]
+
+        combined_df['LS_0'] = combined_df.apply(lambda x: find_match(x, col_name="LS_0"), axis=1)
+        combined_df['LS_1'] = combined_df.apply(lambda x: find_match(x, col_name="LS_1"), axis=1)
+        combined_df['LS_2'] = combined_df.apply(lambda x: find_match(x, col_name="LS_2"), axis=1)
+        combined_df['DS_0'] = combined_df.apply(lambda x: find_match(x, col_name="DS_0"), axis=1)
+        combined_df['DS_1'] = combined_df.apply(lambda x: find_match(x, col_name="DS_1"), axis=1)
+        combined_df['DS_2'] = combined_df.apply(lambda x: find_match(x, col_name="DS_2"), axis=1)
+        combined_df['DS_3'] = combined_df.apply(lambda x: find_match(x, col_name="DS_3"), axis=1)
+        combined_df['haz_expose'] = combined_df.apply(lambda x: find_match(x, col_name="haz_expose"), axis=1)
+
+        # Remove extra columns that are no longer needed
+        combined_df.drop(['w_LS_0', 'w_LS_1', 'w_LS_2', 'sw_LS_0', 'sw_LS_1', 'sw_LS_2', 'f_LS_0', 'f_LS_1',
+                          'f_LS_2', 'w_DS_0', 'w_DS_1', 'w_DS_2', 'w_DS_3', 'sw_DS_0', 'sw_DS_1', 'sw_DS_2', 'sw_DS_3',
+                          'f_DS_0', 'f_DS_1', 'f_DS_2', 'f_DS_3', 'w_haz_expose', 'sw_haz_expose', 'f_haz_expose'],
+                         axis=1, inplace=True)
+
+        return combined_df
 
     def get_spec(self):
         """Get specifications of the combined wind, wave, and surge building damage analysis.
@@ -108,6 +169,12 @@ def get_spec(self):
 
             ],
             'output_datasets': [
+                {
+                    'id': 'ds_result',
+                    'parent_type': 'buildings',
+                    'description': 'CSV file of damage states for building structural damage',
+                    'type': 'ergo:buildingDamageVer6'
+                },
                 {
                     'id': 'result',
                     'parent_type': 'buildings',
diff --git a/pyincore/analyses/combinedwindwavesurgebuildingloss/__init__.py b/pyincore/analyses/combinedwindwavesurgebuildingloss/__init__.py
new file mode 100644
index 000000000..47de6aefe
--- /dev/null
+++ b/pyincore/analyses/combinedwindwavesurgebuildingloss/__init__.py
@@ -0,0 +1,10 @@
+# Copyright (c) 2023 University of Illinois and others. All rights reserved.
+#
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+
+
+from pyincore.analyses.combinedwindwavesurgebuildingloss.combinedwindwavesurgebuildingloss import \
+    CombinedWindWaveSurgeBuildingLoss
+
diff --git a/pyincore/analyses/combinedwindwavesurgebuildingloss/combinedwindwavesurgebuildingloss.py b/pyincore/analyses/combinedwindwavesurgebuildingloss/combinedwindwavesurgebuildingloss.py
new file mode 100755
index 000000000..84919fadb
--- /dev/null
+++ b/pyincore/analyses/combinedwindwavesurgebuildingloss/combinedwindwavesurgebuildingloss.py
@@ -0,0 +1,209 @@
+# Copyright (c) 2023 University of Illinois and others. All rights reserved.
+#
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+
+import pandas as pd
+from pyincore import BaseAnalysis
+import numpy as np
+
+
+class CombinedWindWaveSurgeBuildingLoss(BaseAnalysis):
+    """ 
+    This analysis computes the building structural and content loss from wind, flood and surge-wave damage
+
+    Contributors
+        | Science: Omar Nofal, John W. van de Lindt, Trung Do, Guirong Yan, Sara Hamideh, Daniel Cox, Joel Dietrich
+        | Implementation: Jiate Li, Chris Navarro and NCSA IN-CORE Dev Team
+
+    Related publications
+        Nofal, Omar & Lindt, John & Do, Trung & Yan, Guirong & Hamideh, Sara & Cox, Daniel & Dietrich, Joel. (2021).
+        Methodology for Regional Multi-Hazard Hurricane Damage and Risk Assessment. Journal of Structural Engineering.
+        147. 04021185. 10.1061/(ASCE)ST.1943-541X.0003144.
+        
+    Args:
+        incore_client (IncoreClient): Service authentication.
+    """
+
+    def __init__(self, incore_client):
+        super(CombinedWindWaveSurgeBuildingLoss, self).__init__(incore_client)
+
+    def run(self):
+        """Executes combined wind, wave, surge building loss analysis."""
+
+        # Read buildings
+        buildings = self.get_input_dataset("buildings").get_dataframe_from_shapefile()
+
+        # Read Building wind damage
+        wind_damage = self.get_input_dataset("wind_damage").get_dataframe_from_csv()
+
+        # Read Building surge-wave damage
+        surge_wave_damage = self.get_input_dataset("surge_wave_damage").get_dataframe_from_csv()
+
+        # Read Building flood damage
+        flood_damage = self.get_input_dataset("flood_damage").get_dataframe_from_csv()
+
+        # Read cumulative replacement cost ratio of building content
+        content_cost = self.get_input_dataset("content_cost").get_dataframe_from_csv()
+
+        # Read cumulative replacement cost ratio of structural damage
+        structure_cost = self.get_input_dataset("structural_cost").get_dataframe_from_csv()
+
+        combined_loss = self.get_combined_loss(wind_damage, surge_wave_damage, flood_damage, buildings, content_cost,
+                                               structure_cost)
+
+        # Create the result dataset
+        self.set_result_csv_data("result", combined_loss, self.get_parameter("result_name") + "_loss",
+                                 "dataframe")
+
+        return True
+    
+    def get_combined_loss(self, wind_dmg: pd.DataFrame, sw_dmg: pd.DataFrame, flood_dmg: pd.DataFrame, buildings:
+                          pd.DataFrame, content_cost: pd.DataFrame, structure_cost: pd.DataFrame):
+        """Calculates structural and content loss from wind, surge-wave and flood damage
+
+                Args:
+                    wind_dmg (pd.DataFrame): Table of wind damage for the building inventory
+                    sw_dmg (pd.DataFrame): Table of surge-wave damage for the building inventory
+                    flood_dmg (pd.DataFrame): Table of flood damage for the building inventory
+                    buildings (pd.DataFrame): Table of building attributes
+                    content_cost (pd.DataFrame): Table of content cost ratios for each archetype  
+                    structure_cost (pd.DataFrame): Table of structural cost ratio for each archetype and loss type 
+
+                Returns:
+                    pd.DataFrame: An table of structural and content loss for each building
+
+                """
+        # Rename the columns so there are no overlapping names among the hazard damage states
+        flood_dmg.rename(columns={'LS_0': 'f_LS_0', 'LS_1': 'f_LS_1', 'LS_2': 'f_LS_2', 'DS_0': 'f_DS_0',
+                                  'DS_1': 'f_DS_1', 'DS_2': 'f_DS_2', 'DS_3': 'f_DS_3', 'haz_expose': 'f_haz_expose'},
+                         inplace=True)
+        
+        sw_dmg.rename(columns={'LS_0': 'sw_LS_0', 'LS_1': 'sw_LS_1', 'LS_2': 'sw_LS_2', 'DS_0': 'sw_DS_0',
+                               'DS_1': 'sw_DS_1', 'DS_2': 'sw_DS_2', 'DS_3': 'sw_DS_3', 'haz_expose': 
+                                   'sw_haz_expose'}, inplace=True)
+
+        wind_dmg.rename(columns={'LS_0': 'w_LS_0', 'LS_1': 'w_LS_1', 'LS_2': 'w_LS_2', 'DS_0': 'w_DS_0',
+                                 'DS_1': 'w_DS_1', 'DS_2': 'w_DS_2', 'DS_3': 'w_DS_3', 'haz_expose': 'w_haz_expose'},
+                        inplace=True)
+
+        # Combine all three sets of damages into a single data frame
+        combined_df = pd.merge(pd.merge(wind_dmg, sw_dmg, on='guid'), flood_dmg, on='guid')
+
+        # Add the flood archetype to the combined damage since this is used to find loss multipliers
+        new_combined_df = pd.merge(combined_df, buildings[["guid", "arch_flood"]], on="guid")
+
+        # Create a result data frame for the loss calculations
+        loss_df = new_combined_df[['guid']].copy()
+
+        # Compute content and structural loss
+        loss_df['cont_loss'] = new_combined_df.apply(lambda row: np.dot(row[['f_DS_0', 'f_DS_1', 'f_DS_2']],
+                                                                        content_cost.iloc[int(row['arch_flood']-1),
+                                                                        1:4]), axis=1)
+
+        loss_df['roof_loss'] = new_combined_df.apply(lambda row: np.dot(row[['w_DS_0', 'w_DS_1', 'w_DS_2', 'w_DS_3']],
+                                                                        np.array([0.15, 0.5, 0.75, 1])) *
+                                                     structure_cost.loc[int(row['arch_flood']-1), "Roofing"], axis=1)
+
+        loss_df['ff_loss'] = new_combined_df.apply(lambda row: np.dot(row[['sw_DS_0', 'sw_DS_1', 'sw_DS_2',
+                                                                           'sw_DS_3']], np.array([0.1, 0.5, 0.75, 1])) *
+                                                   structure_cost.loc[int(row['arch_flood'] - 1),
+                                                   "Flooring and Foundation"], axis=1)
+
+        # Computes frame loss from the dominant hazard between wind and surge-wave intensities
+        def compute_frame_loss(row):
+            if row['w_DS_3'] >= row['sw_DS_3']:
+                return np.dot(row[['w_DS_0', 'w_DS_1', 'w_DS_2', 'w_DS_3']],np.array([0.25, 0.5, 0.75, 1])) * \
+                    structure_cost.loc[int(row['arch_flood'] - 1), "Wood Framing"]
+            elif row['sw_DS_3'] > row['w_DS_3']:
+                return np.dot(row[['sw_DS_0', 'sw_DS_1', 'sw_DS_2', 'sw_DS_3']], np.array([0.25, 0.5, 0.75, 1])) * \
+                    structure_cost.loc[int(row['arch_flood'] - 1), "Wood Framing"]
+            else:
+                # If one or both are nan, they can't be compared
+                if pd.isnull(row['sw_DS_3']) and pd.isnull(row['w_DS_3']):
+                    return 0
+                elif not pd.isnull(row['sw_DS_3']):
+                    return np.dot(row[['sw_DS_0', 'sw_DS_1', 'sw_DS_2', 'sw_DS_3']], np.array([0.25, 0.5, 0.75, 1])) * \
+                                  structure_cost.loc[int(row['arch_flood'] - 1), "Wood Framing"]
+                else:
+                    return np.dot(row[['w_DS_0', 'w_DS_1', 'w_DS_2', 'w_DS_3']], np.array([0.25, 0.5, 0.75,
+                                                                                           1])) * structure_cost.loc[
+                        int(row['arch_flood'] - 1), "Wood Framing"]
+
+        loss_df['frame_loss'] = new_combined_df.apply(lambda row: compute_frame_loss(row), axis=1)
+
+        # Fill NA values with 0 otherwise we can't compute the total loss since empty are NaN values and won't add
+        loss_df.fillna(0, inplace=True)
+        loss_df['total_loss'] = loss_df.apply(lambda row: row['cont_loss'] + row['frame_loss'] + row[
+            'roof_loss'] + row['ff_loss'], axis=1)
+
+        return loss_df
+    
+    def get_spec(self):
+        """Get specifications of the combined wind, wave, and surge building loss analysis.
+
+        Returns:
+            obj: A JSON object of specifications of the combined wind, wave, and surge building damage analysis.
+
+        """
+        return {
+            'name': 'combined-wind-wave-surge-building-loss',
+            'description': 'Combined wind wave and surge building loss analysis',
+            'input_parameters': [
+                {
+                    'id': 'result_name',
+                    'required': True,
+                    'description': 'result dataset name',
+                    'type': str
+                },
+            ],
+            'input_datasets': [
+                {
+                    'id': 'buildings',
+                    'required': True,
+                    'description': 'Building Inventory',
+                    'type': ['ergo:buildingInventoryVer4', 'ergo:buildingInventoryVer5',
+                             'ergo:buildingInventoryVer6', 'ergo:buildingInventoryVer7'],
+                },
+                {
+                    'id': 'wind_damage',
+                    'required': True,
+                    'description': 'Wind damage result that has damage intervals in it',
+                    'type': ['ergo:buildingDamageVer6']
+                },
+                {
+                    'id': 'surge_wave_damage',
+                    'required': True,
+                    'description': 'Surge-wave damage result that has damage intervals in it',
+                    'type': ['ergo:buildingDamageVer6']
+                },
+                {
+                    'id': 'flood_damage',
+                    'required': True,
+                    'description': 'Flood damage result that has damage intervals in it',
+                    'type': ['ergo:buildingDamageVer6']
+                },
+                {
+                    'id': 'structural_cost',
+                    'required': True,
+                    'description': 'Structural cost ratio for each archetype and loss type',
+                    'type': ['incore:structuralCostRatio']
+                },
+                {
+                    'id': 'content_cost',
+                    'required': True,
+                    'description': 'Content cost ratio for each archetype and damage state',
+                    'type': ['incore:contentCostRatio']
+                }
+
+            ],
+            'output_datasets': [
+                {
+                    'id': 'result',
+                    'parent_type': 'buildings',
+                    'description': 'CSV file of building structural and content loss',
+                    'type': 'incore:buildingLoss'
+                }
+            ]
+        }
diff --git a/pyincore/analyses/housingrecoverysequential/housingrecoverysequential.py b/pyincore/analyses/housingrecoverysequential/housingrecoverysequential.py
index 60f251c4c..a7b72b991 100644
--- a/pyincore/analyses/housingrecoverysequential/housingrecoverysequential.py
+++ b/pyincore/analyses/housingrecoverysequential/housingrecoverysequential.py
@@ -3,10 +3,12 @@
 # This program and the accompanying materials are made available under the
 # terms of the Mozilla Public License v2.0 which accompanies this distribution,
 # and is available at https://www.mozilla.org/en-US/MPL/2.0/
-from pyincore import BaseAnalysis
+from pyincore import BaseAnalysis, AnalysisUtil
 
 import numpy as np
 import pandas as pd
+import concurrent.futures
+from itertools import repeat
 
 
 class HousingRecoverySequential(BaseAnalysis):
@@ -91,19 +93,15 @@ def run(self):
         t_final = self.get_parameter('t_final')
 
         # Load population block data from IN-CORE
-        pop_dis_selectors = ['guid', 'huid', 'blockid',
-                             'numprec', 'race', 'hispan',
-                             'ownershp', 'randincome', 'dislocated']
+        pop_dis_selectors = ['guid', 'huid', 'blockid', 'race', 'hispan', 'ownershp', 'dislocated']
         households_csv = self.get_input_dataset('population_dislocation_block').get_csv_reader()
         households_df = (pd.DataFrame(households_csv))[pop_dis_selectors]
 
         # Perform  conversions across the dataset from object type into the appropriate type
         households_df['blockid'] = households_df['blockid'].astype('int64')
-        households_df['numprec'] = households_df['numprec'].astype(int)
         households_df['race'] = pd.to_numeric(households_df['race'])
         households_df['hispan'] = pd.to_numeric(households_df['hispan'])
         households_df['ownershp'] = pd.to_numeric(households_df['ownershp'])
-        households_df['randincome'] = pd.to_numeric(households_df['randincome'])
         households_df['dislocated'] = (households_df['dislocated'] == 'True')
 
         # Load the transition probability matrix from IN-CORE
@@ -122,19 +120,68 @@ def run(self):
         initial_prob = pd.DataFrame(list(initial_prob_csv))
         initial_prob['value'] = pd.to_numeric(initial_prob['value'])
 
+        # Obtain the number of CPUs (cores) to execute the analysis with
+        user_defined_cpu = 4
+
+        if not self.get_parameter("num_cpu") is None and self.get_parameter("num_cpu") > 0:
+            user_defined_cpu = self.get_parameter("num_cpu")
+
+        num_workers = AnalysisUtil.determine_parallelism_locally(self, len(households_df), user_defined_cpu)
+
+        # Chop dataset into `num_size` chunks
+        max_chunk_size = int(np.ceil(len(households_df)/num_workers))
+
+        households_df_list = list(
+            map(lambda x: households_df[x * max_chunk_size:x * max_chunk_size + max_chunk_size],
+                list(range(num_workers)))
+        )
+
         # Run the analysis
-        self.housing_serial_recovery_model(t_delta, t_final, tpm, initial_prob, households_df)
+        result = self.hhrs_concurrent_future(self.housing_serial_recovery_model,
+                                                 num_workers,
+                                                 households_df_list,
+                                                 repeat(t_delta),
+                                                 repeat(t_final),
+                                                 repeat(tpm),
+                                                 repeat(initial_prob))
+
+        result_name = self.get_parameter("result_name")
+        self.set_result_csv_data("ds_result", result, name=result_name, source="dataframe")
+
+        return True
+
+    def hhrs_concurrent_future(self, function_name, parallelism, *args):
+        """Utilizes concurrent.future module.
+
+        Args:
+            function_name (function): The function to be parallelized.
+            parallelism (int): Number of workers in parallelization.
+            *args: All the arguments in order to pass into parameter function_name.
+
+        Returns:
+            pd.DataFrame: outcome DataFrame containing the results from the concurrent function.
+
+        """
+        output_ds = pd.DataFrame()
+
+        with concurrent.futures.ProcessPoolExecutor(max_workers=parallelism) as executor:
+            for ret in executor.map(function_name, *args):
+                output_ds = pd.concat([output_ds, ret], ignore_index=True)
 
-    def housing_serial_recovery_model(self, t_delta, t_final, tpm, initial_prob, households_df):
+        return output_ds
+
+    def housing_serial_recovery_model(self, households_df, t_delta, t_final, tpm, initial_prob):
         """Performs the computation of the model as indicated in Sutley and Hamide (2020).
 
         Args:
+            households_df (pd.DataFrame): Households with population dislocation data.
             t_delta (float): Time step size.
             t_final (float): Final time.
             tpm (np.Array): Transition probability matrix.
             initial_prob (pd.DataFrame): Initial probability Markov vector.
-            households_df (pd.DataFrame): Households with population dislocation data.
 
+        Returns:
+            pd.DataFrame: outcome of the HHRS model for a given household dataset.
         """
         seed = self.get_parameter('seed')
         rng = np.random.RandomState(seed)
@@ -270,8 +317,7 @@ def housing_serial_recovery_model(self, t_delta, t_final, tpm, initial_prob, hou
         for i, c in enumerate(column_names):
             result[c] = markov_stages[i]
 
-        result_name = self.get_parameter("result_name")
-        self.set_result_csv_data("ds_result", result, name=result_name, source="dataframe")
+        return result
 
     @staticmethod
     def compute_social_vulnerability_zones(sv_result, households_df):
@@ -352,6 +398,7 @@ def compute_regressions(markov_stages, household, lower, upper):
         regressions = 0
 
         # Check the previous timesteps between [lower, upper] prior to the current timestep.
+
         for t in range(lower, upper):
             if markov_stages[t][household] < markov_stages[t - 1][household]:
                 regressions += 1
@@ -393,6 +440,12 @@ def get_spec(self):
                     'description': 'Seed to ensure replication of the Markov Chain path'
                                    'in connection with Population Dislocation.',
                     'type': int
+                },
+                {
+                    'id': 'num_cpu',
+                    'required': False,
+                    'description': 'If using parallel execution, the number of cpus to request',
+                    'type': int
                 }
             ],
             'input_datasets': [
@@ -416,12 +469,6 @@ def get_spec(self):
                     'description': 'initial mass probability function for stage 0 of the Markov Chain',
                     'type': 'incore:houseRecInitialStageProbability'
                 },
-                {
-                    'id': 'initial_stage_probabilities',
-                    'required': True,
-                    'description': 'initial mass probability function for stage 0 of the Markov Chain',
-                    'type': 'incore:houseRecInitialStageProbability'
-                },
                 {
                     'id': 'sv_result',
                     'required': True,
diff --git a/pyincore/analyses/joplincge/joplincge.py b/pyincore/analyses/joplincge/joplincge.py
index 9a8d813cb..442d8efd7 100644
--- a/pyincore/analyses/joplincge/joplincge.py
+++ b/pyincore/analyses/joplincge/joplincge.py
@@ -59,7 +59,7 @@ def get_spec(self):
                     'description': 'Social accounting matrix (SAM) contains data for firms, '
                                    'households and government which are organized in a way to '
                                    'represent the interactions of all three entities in a typical economy.',
-                    'type': ['incore:JoplinCGEsam']
+                    'type': ['incore:CGEsam']
                 },
                 {
                     'id': 'BB',
@@ -67,14 +67,14 @@ def get_spec(self):
                     'description': 'BB is a matrix which describes how investment in physical infrastructure is'
                                    ' transformed into functioning capital such as commercial and residential buildings.'
                                    ' These data are collected from the Bureau of Economic Analysis (BEA).',
-                    'type': ['incore:JoplinCGEbb']
+                    'type': ['incore:CGEbb']
                 },
                 {
                     'id': 'IOUT',
                     'required': True,
                     'description': 'IOUT is a matrix that describes the transfer of tax revenue collected by the local'
                                    ' government to help finance local government expenditures.',
-                    'type': ['incore:JoplinCGEiout']
+                    'type': ['incore:CGEiout']
                 },
                 {
                     'id': 'MISC',
@@ -82,27 +82,27 @@ def get_spec(self):
                     'description': 'MISC is the name of a file that contains data for commercial sector employment'
                                    ' and physical capital. It also contains data for the number of households and'
                                    ' working households in the economy.',
-                    'type': ['incore:JoplinCGEmisc']
+                    'type': ['incore:CGEmisc']
                 },
                 {
                     'id': 'MISCH',
                     'required': True,
                     'description': 'MISCH is a file that contains elasticities for the supply of labor with'
                                    ' respect to paying income taxes.',
-                    'type': ['incore:JoplinCGEmisch']
+                    'type': ['incore:CGEmisch']
                 },
                 {
                     'id': 'LANDCAP',
                     'required': True,
                     'description': 'LANDCAP contains information regarding elasticity values for the response of '
                                    'changes in the price of physical capital with respect to the supply of investment.',
-                    'type': ['incore:JoplinCGElandcap']
+                    'type': ['incore:CGElandcap']
                 },
                 {
                     'id': 'EMPLOY',
                     'required': True,
                     'description': 'EMPLOY is a table name containing data for commercial sector employment.',
-                    'type': ['incore:JoplinCGEemploy']
+                    'type': ['incore:CGEemploy']
                 },
                 {
                     'id': 'IGTD',
@@ -110,27 +110,27 @@ def get_spec(self):
                     'description': 'IGTD variable represents a matrix describing the transfer of taxes collected'
                                    ' to a variable which permits governments to spend the tax revenue on workers and'
                                    ' intermediate inputs.',
-                    'type': ['incore:JoplinCGEigtd']
+                    'type': ['incore:CGEigtd']
                 },
                 {
                     'id': 'TAUFF',
                     'required': True,
                     'description': 'TAUFF represents social security tax rates',
-                    'type': ['incore:JoplinCGEtauff']
+                    'type': ['incore:CGEtauff']
                 },
                 {
                     'id': 'JOBCR',
                     'required': True,
                     'description': 'JOBCR is a matrix describing the supply of workers'
                                    ' coming from each household group in the economy.',
-                    'type': ['incore:JoplinCGEjobcr']
+                    'type': ['incore:CGEjobcr']
                 },
                 {
                     'id': 'OUTCR',
                     'required': True,
                     'description': 'OUTCR is a matrix describing the number of workers who'
                                    ' live in Joplin but commute outside of town to work.',
-                    'type': ['incore:JoplinCGEoutcr']
+                    'type': ['incore:CGEoutcr']
                 },
                 {
                     'id': 'sector_shocks',
diff --git a/pyincore/analyses/saltlakecge/__init__.py b/pyincore/analyses/saltlakecge/__init__.py
new file mode 100644
index 000000000..361c79902
--- /dev/null
+++ b/pyincore/analyses/saltlakecge/__init__.py
@@ -0,0 +1,7 @@
+# Copyright (c) 2021 University of Illinois and others. All rights reserved.
+#
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+from pyincore.analyses.saltlakecge.saltlakecge import SaltLakeCGEModel
+from pyincore.analyses.saltlakecge.equationlib import *
\ No newline at end of file
diff --git a/pyincore/analyses/saltlakecge/equationlib.py b/pyincore/analyses/saltlakecge/equationlib.py
new file mode 100644
index 000000000..686981c93
--- /dev/null
+++ b/pyincore/analyses/saltlakecge/equationlib.py
@@ -0,0 +1,851 @@
+# Copyright (c) 2021 University of Illinois and others. All rights reserved.
+#
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+from copy import deepcopy
+import pandas as pd
+import operator as op
+import math
+
+
+class VarContainer:
+    """
+     All matrix variable(tables) in the GAMS model is flatten to a array to make a better
+     interface to the solver.
+
+     AllVarList stores all initial values of varibles used in the GAMS model in an array.
+     It also has a indexing system for looking up.
+
+     Attributes:
+         namelist: A dictionary with all stored GAMS variables and its information.
+         nvars: The length of the array, i.e. the size of all matrix variables summed up.
+         initialVals: Stored initial values of all variables
+
+    """
+
+    def __init__(self):
+        """
+          Initialize to an empty list
+        """
+        self.namelist = {}
+        self.initialVals = []
+        self.nvars = 0
+        self.LO = []
+        self.UP = []
+
+    def add(self, name, rows=None, cols=None):
+        """
+        :param name:
+        :param rows:
+        :param cols:
+        :return:
+        """
+        if rows is not None and cols is not None:
+            size = len(rows) * len(cols)
+            self.namelist[name] = {
+                'nrows': len(rows),
+                'ncols': len(cols),
+                'rows': rows,
+                'cols': cols,
+                'start': self.nvars,
+                'size': size
+            }
+        elif rows is not None and cols is None:
+            size = len(rows)
+            self.namelist[name] = {
+                'nrows': len(rows),
+                'rows': rows,
+                'start': self.nvars,
+                'size': size
+            }
+        else:
+            size = 1
+            self.namelist[name] = {'start': self.nvars, 'size': 1}
+
+        self.initialVals.extend([None] * size)
+        self.LO.extend([None] * size)
+        self.UP.extend([None] * size)
+        self.nvars += size
+
+        return ExprM(self, name=name, rows=rows, cols=cols)
+
+    def set_value(self, name, values, target):
+        """
+        An internal method for setting the initial values or UPs and LOs for variables
+
+        :param name: Name of the variable in GAMS
+        :param value: a pandas DataFrame, pandas Series, int or float with initial values
+        :param target: target array to be set
+
+        :return: None
+        """
+        if type(values) == int or type(values) == float:
+            info = self.namelist[name]
+            if 'nrows' in info and 'ncols' in info:
+                values = pd.DataFrame(index=info['rows'], columns=info['cols']).fillna(values)
+            elif 'nrows' in info and 'ncols' not in info:
+                values = pd.Series(index=info['rows'], dtype=pd.Float64Dtype).fillna(values)
+
+        if type(values) == pd.DataFrame:
+            rows = values.index.tolist()
+            cols = values.columns.tolist()
+            for i in rows:
+                for j in cols:
+                    target[self.getIndex(name, row=i, col=j)] = float(values.loc[i][j])
+        elif type(values) == pd.Series:
+            rows = values.index.tolist()
+            for i in rows:
+                target[self.getIndex(name, row=i)] = float(values.loc[i])
+        else:
+            target[self.getIndex(name)] = values
+
+    def init(self, name, initialValue):
+        """
+        Flatten the table variable and add to the list.
+        Also set the initial variable values array.
+
+        :param name: Name of the variable in GAMS
+        :param initialValue: a pandas DataFrame or pandas Series with initial values
+
+        :return: None.
+        """
+        self.set_value(name, initialValue, self.initialVals)
+
+    def inList(self, name):
+        """
+        Check if a GAMS varible is added to the container
+
+        :param name(str): name of GAMS variable you want to look up
+        :return: Boolean, whether the variable is added.
+        """
+        return name in self.namelist
+
+    def getInfo(self, name):
+        """
+        Get the information about a GAMS variable
+
+        :param name(str): name of GAMS variable you want to look up
+        :return: a dictionary with all information
+        """
+        return self.namelist[name]
+
+    def getIndex(self, name, row=None, col=None):
+        """
+        Look up the index by providing the variable name and label information
+
+        :param name: name of GAMS variable you want to look up
+        :param row: row label of the position you want to look up index for(if it has row labels)
+        :param col: column label of the position you want to look up index for(if it has column labels)
+        :return: the index of the position in the array
+        """
+        info = self.namelist[name]
+        result = info['start']
+        if row is not None and col is not None:
+            result += info['rows'].index(row) * info['ncols'] + info['cols'].index(col)
+        elif row is not None:
+            result += info['rows'].index(row)
+        return result
+
+    def getLabel(self, index):
+        """
+        Look up variable name and label information by providing the index
+
+        :param index: the index in the array
+        :return: its information including the variable name, row label and column label if applicable
+        """
+        result = []
+        for i in self.namelist.keys():
+            if index >= self.namelist[i]['start'] and index < self.namelist[i]['start'] + self.namelist[i]['size']:
+                result.append(i)
+                if self.namelist[i]['size'] > 1:
+                    diff = index - self.namelist[i]['start']
+                    if 'ncols' in self.namelist[i]:
+                        result.append(self.namelist[i]['rows'][int(diff / self.namelist[i]['ncols'])])
+                        result.append(self.namelist[i]['cols'][diff % self.namelist[i]['ncols']])
+                    else:
+                        result.append(self.namelist[i]['rows'][diff])
+                return result
+
+    def get_all_variable_printed(self, output=None):
+        if output is None:
+            output = self.initialVals
+        for i in range(len(output)):
+            label = self.getLabel(i)
+            if len(label) == 1:
+                print(label[0] + '=' + '{0:.7f}'.format(output[i]) + ';')
+            elif len(label) == 2:
+                print(label[0] + '(\'' + label[1] + '\')' + '=' + '{0:.7f}'.format(output[i]) + ';')
+            elif len(label) == 3:
+                print(
+                    label[0] + '(\'' + label[1] + '\',\'' + label[2] + '\')' + '=' + '{0:.7f}'.format(output[i]) + ';')
+
+    def get(self, name, x=None):
+        """
+          Returns a Dataframe, Series, or a variable based on the given name
+          and the result array returned from the solver
+
+          :param name: GAMS variable name
+          :return: if x is not given, it returns the initial values
+            if x is set to the result, returns the result variable value
+        """
+        if x is None:
+            x = self.initialVals
+
+        info = self.namelist[name]
+        if 'nrows' in info and 'ncols' in info:
+            ret = pd.DataFrame(index=info['rows'], columns=info['cols']).fillna(0.0)
+            for i in info['rows']:
+                for j in info['cols']:
+                    ret.at[i, j] = x[self.getIndex(name, row=i, col=j)]
+        elif 'nrows' in info and 'ncols' not in info:
+            ret = pd.Series(index=info['rows'], dtype=pd.Float64Dtype).fillna(0.0)
+            for i in info['rows']:
+                ret.at[i] = x[self.getIndex(name, row=i)]
+        elif 'nrows' not in info and 'ncols' not in info:
+            ret = x[self.getIndex(name)]
+
+        return ret
+
+    def lo(self, name, value):
+        """
+          Set the LOs of a GAMS variable providing the LOs with a Dataframe, Series, int or float
+
+          :param name: GAMS variable name
+          :param value: The lower bound to be set
+          :return: None
+        """
+        self.set_value(name, value, self.LO)
+
+    def up(self, name, value):
+        """
+          Set the UPs of a GAMS variable providing the LOs with a Dataframe, Series, int or float
+
+          :param name: GAMS variable name
+          :param value: The upper bound to be set
+          :return: None
+        """
+        self.set_value(name, value, self.UP)
+
+    def write(self, filename):
+        """
+          Write(append) the variables to a file, in the format of setting ipopt model variables
+
+          :param filename: the output filename
+          :return: None
+        """
+        with open(filename, 'a') as f:
+            for i in range(self.nvars):
+                lower = -1e20 if self.LO[i] == None else self.LO[i]
+                upper = 1e20 if self.UP[i] == None else self.UP[i]
+                value = 0 if math.isnan(self.initialVals[i]) else self.initialVals[i]
+
+                # tmp_list = ['CG','CH','CN','CPI','CPIN','CPIH','CX','D','DD','DS','FD','GCP','GCP1',
+                # 'HH','HN','HW','IGT','KS']
+                '''tmp_list = ['CG','CH','CN','CPI','CPIN','CPIH','CX','D','DD','DS','FD','GCP','GCP1',
+                'HH','HN','HW','IGT','KS','LAS','M','N','NKI','LNFOR','KPFOR','GVFOR','P','PD','PVA']                        
+                if self.getLabel(i)[0] in tmp_list:
+                    f.write('model.x' + str(i) + ' = Var(bounds=(' + str(lower) + ',' + str(upper) + '),initialize=' + 
+                    str(self.initialVals[i]) + ')' + '\n')
+                else:
+                    f.write('model.x' + str(i) + ' = Var(bounds=(' + str(self.initialVals[i]-1e-3) + ',' + 
+                    str(self.initialVals[i]+1e-3) + '),initialize=' + str(self.initialVals[i]) + ')' + '\n')
+                '''
+                # if i == 1507:
+                #  f.write('model.x' + str(i) + ' = Var(bounds=(' + str(value) + ',' + str(value*10) + '),initialize=
+                #  ' + str(value) + ')' + '\n')
+                # else:
+                #  f.write('model.x' + str(i) + ' = Var(bounds=(' + str(lower) + ',' + str(upper) + '),initialize=' +
+                #  str(value) + ')' + '\n')
+                f.write('model.x' + str(i) + ' = Var(bounds=(' + str(lower) + ',' + str(upper) + '),initialize=' +
+                        str(value) + ')' + '\n')
+
+
+class variable:
+    """
+      A GMAS variable, initialized by given the GAMS variable value and its label
+    """
+
+    def __init__(self, vars, name, row=None, col=None):
+        """
+          Initialize it with a variable container, the GAMS name, the labels
+
+          :param vars: the variable container that already added the GAMS variable
+          :param name: GAMS variable name
+          :param row: GAMS row label if there is
+          :param col: GAMS col label if there is
+          :return: None
+        """
+        try:
+            self.index = vars.getIndex(name, row, col)
+        except Exception as e:
+            print(e)
+            # print("invalid name for a variable")
+
+    def __str__(self):
+        """
+          returns the variable in the format of "model.x#" if gets printed,
+          with # being the index in the array in the container
+
+          :return: String
+        """
+        return 'model.x' + str(self.index) + ''
+        # return 'x[' + str(self.index) + ']'
+
+    def debug_test_str(self):
+        return 'x[' + str(self.index) + ']'
+
+
+class ExprItem:
+    '''
+      You can construct it with a variable, a constant or a deepcopy of another ExprItem
+    '''
+
+    def __init__(self, v, const=1):
+        self.varList = []
+        self.const = const
+        if type(v) == variable:
+            self.varList.append(v)
+        elif type(v) == int or type(v) == float:
+            self.const = v
+        elif type(v) == ExprItem:
+            self.const = deepcopy(v.const)
+            self.varList = deepcopy(v.varList)
+        else:
+            print("invalid parameter to create a item")
+
+    '''
+      You could multiply it with a number, a variable, a ExprItem or Expression
+    '''
+
+    def __mul__(self, rhs):
+        copy = ExprItem(self)
+        if type(rhs) == int or type(rhs) == float:
+            copy.const = copy.const * rhs
+        elif type(rhs) == variable:
+            copy.varList.append(rhs)
+        elif type(rhs) == ExprItem:
+            copy.const *= rhs.const
+            copy.varList.extend(deepcopy(rhs.varList))
+        elif type(rhs) == Expr:
+            if rhs.isComposite:
+                if rhs.operator == '/':
+                    return rhs * copy
+                else:
+                    copy.varList.append(Expr(rhs))
+            else:
+                copyrhs = Expr(rhs)
+                copyrhs.itemList = [i * copy for i in copyrhs.itemList]
+                return copyrhs
+        return copy
+
+    def __str__(self):
+        # if abs(self.const) > 0.1:
+        #  result = '' + str(round(self.const,3))
+        # else:
+        #  result = '' + str(self.const)
+        # result = '' + str(round(self.const,7))
+        result = '' + str(self.const)
+        for i in range(len(self.varList)):
+            result += "*" + self.varList[i].__str__()
+        return result
+
+    def debug_test_str(self):
+        # if abs(self.const) > 0.1:
+        #  result = '' + str(round(self.const,3))
+        # else:
+        #  result = '' + str(self.const)
+        result = '' + str(self.const)
+        for i in range(len(self.varList)):
+            result += "*" + self.varList[i].debug_test_str()
+        return result
+
+    def is_empty(self):
+        if abs(self.const) < 0.00000001:
+            return True
+        else:
+            return False
+
+    def is_one(self):
+        if abs(self.const - 1) < 0.0000001 and len(self.varList) == 0:
+            return True
+        else:
+            return False
+
+
+class Expr:
+
+    def __init__(self, item):
+        self.itemList = []
+        self.isComposite = False
+        if type(item) == ExprItem or type(item) == variable or type(item) == int or type(item) == float:
+            self.itemList.append(ExprItem(item))
+        elif type(item) == Expr:
+            self.itemList = [
+                ExprItem(item.itemList[i]) if type(item.itemList[i]) == ExprItem else Expr(item.itemList[i]) for i in
+                range(len(item.itemList))]
+            try:
+                self.isComposite = deepcopy(item.isComposite)
+            except Exception as e:
+                pass
+            try:
+                self.operator = deepcopy(item.operator)
+            except Exception as e:
+                pass
+            try:
+                self.first = Expr(item.first)
+            except Exception as e:
+                pass
+            try:
+                self.second = Expr(item.second)
+            except Exception as e:
+                pass
+        else:
+            print("invalid parameter for creating a Expr")
+
+        self.clear_empty()
+
+    def __add__(self, rhs):
+        copy = Expr(self)
+        if self.isComposite:
+            tmp = Expr(self)
+            copy.itemList = [tmp]
+            copy.isComposite = False
+        if type(rhs) == Expr:
+            if rhs.isComposite:
+                rhscopy = Expr(rhs)
+                copy.itemList.append(rhscopy)
+            else:
+                rhscopy = Expr(rhs)
+                copy.itemList = copy.itemList + rhscopy.itemList
+        elif type(rhs) == ExprItem or type(rhs) == int or type(rhs) == float or type(rhs) == variable:
+            copy.itemList.append(ExprItem(rhs))
+        return copy
+
+    def __sub__(self, rhs):
+        copy = Expr(self)
+        return copy + rhs * -1
+
+    def __mul__(self, rhs):
+        copy = Expr(self)
+        if type(rhs) == int or type(rhs) == float or type(rhs) == ExprItem:
+            if copy.isComposite:
+                if copy.operator == '/':
+                    copy.first = copy.first * rhs
+                elif copy.opeartor == '**':
+                    return ExprItem(rhs) * copy
+            else:
+                result = []
+                for i in copy.itemList:
+                    result.append(i * rhs)
+                copy.itemList = result
+
+
+        elif type(rhs) == Expr:
+            if copy.isComposite and not rhs.isComposite:
+                if copy.operator == '/':
+                    copy.first = copy.first * rhs
+                else:
+                    copy.itemList = [i * copy for i in rhs.itemList]
+                    copy.isComposite = False
+            elif not copy.isComposite and not rhs.isComposite:
+                # print('copy',copy)
+                # print('rhs',rhs)
+
+                if len(copy.itemList) * len(rhs.itemList) > 10:
+
+                    tmpItem = ExprItem(1)
+                    tmpItem.varList.append(Expr(copy))
+                    tmpItem.varList.append(Expr(rhs))
+                    copy.itemList = [tmpItem]
+
+                else:
+
+                    # The expanding way is super slow for long equations
+
+                    result = []
+                    for i in copy.itemList:
+                        for j in rhs.itemList:
+                            if type(i) == Expr:
+                                tmp = Expr(i)
+                            elif type(i) == ExprItem:
+                                tmp = ExprItem(i)
+                            result.append(tmp * j)
+                    copy.itemList = result
+
+            elif not self.isComposite and rhs.isComposite:
+                return rhs * copy
+            else:
+                # both are composite
+                if copy.operator == '/' and rhs.operator == '/':
+                    copy.first = copy.first * rhs.first
+                    copy.second = copy.second * rhs.second
+                elif copy.operator == '**' and rhs.operator == '**':
+                    copy.first = copy.first * rhs.first
+                    copy.second = copy.second + rhs.second
+                elif copy.operator == '/' and rhs.operator == '**':
+                    copy.first = copy.first * rhs
+                elif copy.operator == '**' and rhs.operator == '/':
+                    copyrhs = Expr(rhs)
+                    copyrhs.first = copyrhs.first * copy
+        return copy
+
+    def __truediv__(self, rhs):
+        copy = Expr(self)
+        tmp = Expr(copy)
+        copy.isComposite = True
+        copy.operator = '/'
+        copy.first = tmp
+        copy.second = Expr(rhs)
+        return copy
+
+    def __pow__(self, rhs):
+        copy = Expr(self)
+        tmp = Expr(copy)
+        copy.isComposite = True
+        copy.operator = '**'
+        copy.first = tmp
+        copy.second = Expr(rhs)
+        return copy
+
+    def __str__(self):
+        if self.isComposite:
+            result = "(" + self.first.__str__() + self.operator + self.second.__str__() + ")"
+        else:
+            if len(self.itemList) == 0:
+                return '0'
+            result = ''
+            for i in self.itemList[0:-1]:
+                result += "(" + i.__str__() + ")" + "+"
+            result += "(" + self.itemList[-1].__str__() + ")"
+            if len(self.itemList) > 1:
+                result = "(" + result + ")"
+        return result
+
+    def debug_test_str(self):
+        if self.isComposite:
+            result = "(" + self.first.debug_test_str() + self.operator + self.second.debug_test_str() + ")"
+        else:
+            if len(self.itemList) == 0:
+                return '0'
+            result = ''
+            for i in self.itemList[0:-1]:
+                result += "(" + i.debug_test_str() + ")" + "+"
+            result += "(" + self.itemList[-1].debug_test_str() + ")"
+            if len(self.itemList) > 1:
+                result = "(" + result + ")"
+        return result
+
+    def is_empty(self):
+        if not self.isComposite and len(self.itemList) == 0:
+            return True
+        else:
+            return False
+
+    def clear_empty(self):
+        result = []
+        if self.isComposite == False:
+            for i in range(len(self.itemList)):
+                if type(self.itemList[i]) == ExprItem:
+                    if not self.itemList[i].is_empty():
+                        result.append(self.itemList[i])
+                if type(self.itemList[i]) == Expr:
+                    self.itemList[i].clear_empty()
+                    if not self.itemList[i].is_empty():
+                        result.append(self.itemList[i])
+            self.itemList = result
+        # else:
+        #   if self.operator == '**' and is_empty(self.second):
+        #     self.isComposite = False
+        #     self.itemList = []
+        #     self.itemList.append(ExprItem(1))
+
+
+
+class ExprM:
+    """
+    Three ways to create a ExprMatrix:
+    1. Give it the variable name, selected rows and cols(could be empty),
+      The constructor will create a Expression matrix from the variable matrix
+    2. Give it a pandas Series or DataFrame, it will create the Expression matrix
+      with the content in the Series or DataFrame as constants
+    3. Give it a ExprMatrix, will return a deep copy of it
+    """
+    def __init__(self, vars, name=None, rows=None, cols=None, m=None, em=None):
+        self.vars = vars
+        self.hasCondition = False
+        if em is None:
+            # if these are the variables, we need to create an Expression by the variable name
+            if name is not None:
+                if self.vars.inList(name):
+                    self.info = deepcopy(self.vars.getInfo(name))
+                else:
+                    print("Can't find this variable in the all variable list")
+
+                self.info['height'] = 1
+                self.info['width'] = 1
+                self.info['rows'] = deepcopy(rows)
+                self.info['cols'] = deepcopy(cols)
+                if cols is not None:
+                    self.info['width'] = len(cols)
+                if rows is not None:
+                    self.info['height'] = len(rows)
+
+                if cols is not None:  # if it is a DataFrame
+                    self.m = [[Expr(variable(self.vars, name, i, j)) for j in cols] for i in rows]
+                elif rows is not None:  # if it is a Series
+                    self.m = [[Expr(variable(self.vars, name, i))] for i in rows]
+                else:  # if it is a variable
+                    self.m = [[Expr(variable(self.vars, name))]]
+
+            # otherwise these are just constants
+            else:
+                self.info = {}
+                self.info['height'] = 1
+                self.info['width'] = 1
+                self.info['rows'] = None
+                self.info['cols'] = None
+                if type(m) == pd.DataFrame:
+                    self.info['rows'] = m.index.tolist()
+                    self.info['height'] = len(self.info['rows'])
+                    self.info['cols'] = m.columns.tolist()
+                    self.info['width'] = len(self.info['cols'])
+                    self.m = [[Expr(float(m.loc[i, j])) for j in self.info['cols']] for i in self.info['rows']]
+                elif type(m) == pd.Series:
+                    self.info['rows'] = m.index.tolist()
+                    self.info['height'] = len(self.info['rows'])
+                    self.m = [[Expr(float(m.loc[i]))] for i in self.info['rows']]
+                else:
+                    self.m = [[Expr(float(m))]]
+        else:
+            # make deep copy
+            self.info = deepcopy(em.info)
+            self.m = [[Expr(em.m[i][j]) for j in range(em.info['width'])] for i in range(em.info['height'])]
+
+    def operation(self, rhs, oper):
+        copy = ExprM(self.vars, em=self)
+
+        if type(rhs) == int or type(rhs) == float:
+            for i in range(self.info['height']):
+                for j in range(self.info['width']):
+                    copy.m[i][j] = oper(copy.m[i][j], Expr(rhs))
+
+        elif type(rhs) == pd.DataFrame or type(rhs) == pd.Series:
+            return self.operation(ExprM(self.vars, m=rhs), oper)
+
+        elif type(rhs) == ExprM:
+            if copy.info['rows'] == rhs.info['rows'] and copy.info['cols'] == rhs.info['cols']:
+                # same size, apply operation element-wise
+                for i in range(self.info['height']):
+                    for j in range(self.info['width']):
+                        copy.m[i][j] = oper(copy.m[i][j], rhs.m[i][j])
+            elif copy.info['rows'] == rhs.info['rows'] and copy.info['width'] == 1:
+                # for each column in rhs
+                result = [[Expr(0) for j in range(rhs.info['width'])] for i in range(self.info['height'])]
+                for i in range(self.info['height']):
+                    for j in range(rhs.info['width']):
+                        result[i][j] = oper(copy.m[i][0], rhs.m[i][j])
+                copy.m = result
+                copy.info['width'] = rhs.info['width']
+                copy.info['cols'] = deepcopy(rhs.info['cols'])
+            elif copy.info['rows'] == rhs.info['rows'] and rhs.info['width'] == 1:
+                # for each column in lhs
+                for i in range(self.info['height']):
+                    for j in range(self.info['width']):
+                        copy.m[i][j] = oper(copy.m[i][j], rhs.m[i][0])
+            elif copy.info['height'] == 1 and copy.info['cols'] == rhs.info['cols']:
+                # for each row in rhs
+                result = [[Expr(0) for j in range(rhs.info['width'])] for i in range(rhs.info['height'])]
+                for i in range(rhs.info['height']):
+                    for j in range(rhs.info['width']):
+                        result[i][j] = oper(copy.m[0][j], rhs.m[i][j])
+                copy.m = result
+                copy.info['height'] = rhs.info['height']
+                copy.info['rows'] = deepcopy(rhs.info['rows'])
+            elif rhs.info['height'] == 1 and copy.info['cols'] == rhs.info['cols']:
+                # for each row in lhs
+                for i in range(self.info['height']):
+                    for j in range(self.info['width']):
+                        copy.m[i][j] = oper(copy.m[i][j], rhs.m[0][j])
+            elif (self.info['width'] == rhs.info['height'] and self.info['cols'] == rhs.info['rows']) or \
+                    (self.info['height'] == rhs.info['width'] and self.info['rows'] == rhs.info['cols']):
+                # flip the matrix
+                return oper(self, ~rhs)
+
+            else:
+                print(copy.info['rows'], copy.info['cols'], rhs.info['rows'], rhs.info['cols'])
+                print(copy.info['height'], copy.info['width'], rhs.info['height'], rhs.info['width'])
+                print("Invalid size for ", str(oper))
+        return copy
+
+    def __add__(self, rhs):
+        return self.operation(rhs, op.add)
+
+    def __radd__(self, lhs):
+        return self + lhs
+
+    def __sub__(self, rhs):
+        return self.operation(rhs, op.sub)
+
+    def __rsub__(self, lhs):
+        return self * -1 + lhs
+
+    def __mul__(self, rhs):
+        return self.operation(rhs, op.mul)
+
+    def __truediv__(self, rhs):
+        return self.operation(rhs, op.truediv)
+
+    def __pow__(self, rhs):
+        return self.operation(rhs, op.pow)
+
+    def __xor__(self, rhs):
+        """
+          create 2d list out of 2 single lists
+        """
+        # has to be 2 single lists
+        if self.info['width'] != 1 or rhs.info['width'] != 1:
+            print("Invalid size for creating a 2-D matrix")
+        else:
+            copy = ExprM(self.vars, em=self)
+            copy.m = [[copy.m[i][0] * rhs.m[j][0] for j in range(rhs.info['height'])] for i in
+                      range(copy.info['height'])]
+            copy.info['cols'] = deepcopy(rhs.info['rows'])
+            copy.info['width'] = deepcopy(rhs.info['height'])
+            return copy
+
+    def __invert__(self):
+        """
+          Return the transpose of a Expression matrix
+        """
+        copy = ExprM(self.vars, em=self)
+        result = [[copy.m[i][j] for i in range(copy.info['height'])] for j in range(copy.info['width'])]
+        copy.info['height'], copy.info['width'] = copy.info['width'], copy.info['height']
+        copy.info['rows'], copy.info['cols'] = copy.info['cols'], copy.info['rows']
+        copy.m = result
+        return copy
+
+    def __str__(self):
+        result = ''
+        for i in self.m:
+            for j in i:
+                result += j.__str__() + '\n'
+            result += '///////////////////\n'
+        return result
+
+    def loc(self, rows=None, cols=None):
+        """
+          get a subset of the matrix by labels
+        """
+        copy = ExprM(self.vars, em=self)
+        if cols is not None:
+            result = [[Expr(copy.m[self.info['rows'].index(i)][self.info['cols'].index(j)]) for j in cols] for i in
+                      rows]
+            copy.m = result
+            copy.info['rows'] = deepcopy(rows)
+            copy.info['height'] = len(copy.info['rows'])
+            copy.info['cols'] = deepcopy(cols)
+            copy.info['width'] = len(copy.info['cols'])
+            return copy
+        elif rows is not None:
+            result = [[Expr(copy.m[self.info['rows'].index(i)][0])] for i in rows]
+            copy.m = result
+            copy.info['rows'] = deepcopy(rows)
+            copy.info['height'] = len(copy.info['rows'])
+            return copy
+        else:
+            return copy
+
+    def sum(self, label=None):
+        copy = ExprM(self.vars, em=self)
+        if label == None:
+            result = Expr(0)
+            for i in range(self.info['height']):
+                for j in range(self.info['width']):
+                    result = result + copy.m[i][j]
+            copy.m = [[result]]
+            copy.info['width'] = 1
+            copy.info['height'] = 1
+            copy.info['rows'] = None
+            copy.info['cols'] = None
+
+        elif label == self.info['rows'] or label == 0:
+            result = [[Expr(0) for j in range(self.info['width'])]]
+            for i in range(self.info['height']):
+                for j in range(self.info['width']):
+                    if not self.hasCondition or self.hasCondition and self.mark[i][j]:
+                        result[0][j] = result[0][j] + copy.m[i][j]
+            copy.m = result
+            copy.info['height'] = 1
+            copy.info['rows'] = None
+
+        elif label == self.info['cols'] or label == 1:
+            result = [[Expr(0)] for i in range(self.info['height'])]
+            for i in range(self.info['height']):
+                for j in range(self.info['width']):
+                    if not self.hasCondition or self.hasCondition and self.mark[i][j]:
+                        result[i][0] = result[i][0] + copy.m[i][j]
+            copy.m = result
+            copy.info['width'] = 1
+            copy.info['cols'] = None
+        return copy
+
+    def prod(self, label):
+        copy = ExprM(self.vars, em=self)
+        if label == self.info['rows'] or label == 0:
+            result = [[Expr(1) for j in range(self.info['width'])]]
+            for i in range(self.info['height']):
+                for j in range(self.info['width']):
+                    if not self.hasCondition or self.hasCondition and self.mark[i][j]:
+                        result[0][j] = result[0][j] * copy.m[i][j]
+            copy.m = result
+            copy.info['height'] = 1
+            copy.info['rows'] = None
+        elif label == self.info['cols'] or label == 1:
+            result = [[Expr(1)] for i in range(self.info['height'])]
+            for i in range(self.info['height']):
+                for j in range(self.info['width']):
+                    if not self.hasCondition or self.hasCondition and self.mark[i][j]:
+                        result[i][0] = result[i][0] * copy.m[i][j]
+            copy.m = result
+            copy.info['width'] = 1
+            copy.info['cols'] = None
+        return copy
+
+    def setCondition(self, matrix, operator=None, value=None):
+        self.hasCondition = True
+        mappings = {'LT': op.lt, 'LE': op.le, 'EQ': op.eq, 'INEQ': op.ne}  # not complete list
+        if type(matrix) == pd.DataFrame:
+            self.mark = [[False for j in range(self.info['width'])] for i in range(self.info['height'])]
+            for i in range(self.info['height']):
+                for j in range(self.info['width']):
+                    if operator == None:
+                        self.mark[i][j] = op.ne(matrix.loc[self.info['rows'][i]][self.info['cols'][j]], 0)
+                    else:
+                        self.mark[i][j] = mappings[operator](matrix.loc[self.info['rows'][i]][self.info['cols'][j]],
+                                                             value)
+
+        if type(matrix) == pd.Series:
+            self.mark = [[False] for i in range(self.info['height'])]
+            for i in range(self.info['height']):
+                if operator == None:
+                    self.mark[i][0] = op.ne(matrix.loc[self.info['rows'][i]], 0)
+                else:
+                    self.mark[i][0] = mappings[operator](matrix.loc[self.info['rows'][i]], value)
+
+    def write(self, count, filename):
+        f = open(filename, 'a')
+        for i in range(self.info['height']):
+            for j in range(self.info['width']):
+                if not self.hasCondition or self.hasCondition and self.mark[i][j]:
+                    f.write(
+                        'model.equality' + str(count[0]) + ' = Constraint(expr=' + self.m[i][j].__str__() + ' == 0)\n')
+                    count[0] += 1
+        f.close()
+
+    def test(self, x):
+        for i in range(self.info['height']):
+            for j in range(self.info['width']):
+                if not self.hasCondition or self.hasCondition and self.mark[i][j]:
+                    fun = lambda x: eval(self.m[i][j].debug_test_str())
+                    print(i, j, fun(x))
diff --git a/pyincore/analyses/saltlakecge/outputfunctions.py b/pyincore/analyses/saltlakecge/outputfunctions.py
new file mode 100644
index 000000000..7c579ea22
--- /dev/null
+++ b/pyincore/analyses/saltlakecge/outputfunctions.py
@@ -0,0 +1,69 @@
+# Copyright (c) 2021 University of Illinois and others. All rights reserved.
+#
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+from pyincore.analyses.saltlakecge.equationlib import *
+import pandas as pd
+import numpy as np
+
+# Return initial values function
+def baseValue(vars, soln, eqName):
+    if vars.getInfo(eqName)['size'] == 1:
+        basevalues = vars.initialVals[vars.getIndex(eqName)]
+
+    else:
+        rows = vars.getInfo(eqName)['rows']
+        rws = []
+        if vars.getInfo(eqName)['size'] == vars.getInfo(eqName)['nrows']:
+            for rr in rows:
+                number = vars.initialVals[vars.getIndex(eqName, row=rr)]
+                rws.append(number)
+            basevalues = pd.DataFrame(rws, index=rows)
+
+        else:
+            cols = vars.getInfo(eqName)['cols']
+            for rr in rows:
+                clms = []
+                for cc in cols:
+                    number = vars.initialVals[vars.getIndex(eqName, row=rr, col=cc)]
+                    clms.append(number)
+                rws.append(clms)
+            basevalues = pd.DataFrame(rws, index=rows, columns=cols)
+
+    return basevalues
+
+
+# Get the new values function
+def newValue(vars, soln, eqName, ittr):
+    #    for scalar
+    if vars.getInfo(eqName)['size'] == 1:
+        newvalue = soln[ittr][vars.getIndex(eqName)]
+
+    else:
+        rows = vars.getInfo(eqName)['rows']
+        rws = []
+        #       vectors
+        if vars.getInfo(eqName)['size'] == vars.getInfo(eqName)['nrows']:
+            for rr in rows:
+                number = soln[ittr][vars.getIndex(eqName, row=rr)]
+                rws.append(number)
+            newvalue = pd.DataFrame(rws, index=rows)
+        #       data frames
+        else:
+            for rr in rows:
+                cols = vars.getInfo(eqName)['cols']
+                clms = []
+                for cc in cols:
+                    number = soln[ittr][vars.getIndex(eqName, row=rr, col=cc)]
+                    clms.append(number)
+                rws.append(clms)
+            newvalue = pd.DataFrame(rws, index=rows, columns=cols)
+
+    return newvalue
+
+
+# Get differences function
+def getDiff(vars, soln, eqName, ittr):
+    difference = newValue(vars, soln, eqName, ittr) - baseValue(vars, soln, eqName)
+    return difference
\ No newline at end of file
diff --git a/pyincore/analyses/saltlakecge/saltlakecge.py b/pyincore/analyses/saltlakecge/saltlakecge.py
new file mode 100644
index 000000000..2bd83a0c4
--- /dev/null
+++ b/pyincore/analyses/saltlakecge/saltlakecge.py
@@ -0,0 +1,2088 @@
+# Copyright (c) 2021 University of Illinois and others. All rights reserved.
+#
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+
+# increase recursion limit so that LSEQ1 can be solved
+import os
+from pyincore import globals as pyglobals
+from pyincore import BaseAnalysis
+from pyincore.analyses.saltlakecge.equationlib import *
+from pyincore.analyses.saltlakecge.outputfunctions import *
+from pyincore.analyses.saltlakecge.saltlakeoutput import gams_to_dataframes
+
+from pyomo.environ import *
+from pyomo.opt import SolverFactory, SolverStatus, TerminationCondition
+
+import pandas as pd
+import sys
+
+
+class SaltLakeCGEModel(BaseAnalysis):
+    """A computable general equilibrium (CGE) model is based on fundamental economic principles.
+    A CGE model uses multiple data sources to reflect the interactions of households,
+    firms and relevant government entities as they contribute to economic activity.
+    The model is based on (1) utility-maximizing households that supply labor and capital,
+    using the proceeds to pay for goods and services (both locally produced and imported)
+    and taxes; (2) the production sector, with perfectly competitive, profit-maximizing firms
+    using intermediate inputs, capital, land and labor to produce goods and services for both
+    domestic consumption and export; (3) the government sector that collects taxes and uses
+    tax revenues in order to finance the provision of public services; and (4) the rest of the world.
+
+
+
+    Args:
+        incore_client (IncoreClient): Service authentication.
+
+    """
+
+    def __init__(self, incore_client):
+        super(SaltLakeCGEModel, self).__init__(incore_client)
+        sys.setrecursionlimit(10000)
+
+    def run(self):
+        """
+
+        Returns:
+
+        """
+        iNum = self.get_parameter("model_iterations")
+
+        # TODO: Update SPEC data types to be generic
+        SAM = pd.read_csv(self.get_input_dataset("SAM").get_file_path('csv'), index_col=0)
+        BB = pd.read_csv(self.get_input_dataset("BB").get_file_path('csv'), index_col=0)
+        JOBCR = pd.read_csv(self.get_input_dataset("JOBCR").get_file_path('csv'), index_col=0)
+        MISCH = pd.read_csv(self.get_input_dataset("MISCH").get_file_path('csv'), index_col=0)
+        EMPLOY = pd.read_csv(self.get_input_dataset("EMPLOY").get_file_path('csv'), index_col=0)
+        OUTCR = pd.read_csv(self.get_input_dataset("OUTCR").get_file_path('csv'), index_col=0)
+        sector_shocks = pd.read_csv(self.get_input_dataset("sector_shocks").get_file_path('csv'), index_col=0)
+
+        self.salt_lake_city_cge(iNum, SAM, BB, JOBCR, MISCH, EMPLOY, OUTCR, sector_shocks)
+
+    def salt_lake_city_cge(self, iNum, SAM, BB, JOBCR, MISCH, EMPLOY, OUTCR, sector_shocks):
+        """
+
+        Args:
+            iNum (int):
+            SAM (pd.DataFrame):
+            BB (str):
+            JOBCR:
+            MISCH:
+            EMPLOY:
+            OUTCR:
+            sector_shocks:
+
+        Returns:
+
+        """
+
+        def _(x):
+            return ExprM(vars, m=x)
+
+        # ----------------------------------------------------------------
+        # define sets
+        # ----------------------------------------------------------------
+
+        # ALL ACCOUNTS IN SOCIAL ACCOUNTING MATRIX
+        Z = [
+            'AG_MI_R1',
+            'UTIL_R1',
+            'CONS_R1',
+            'MANU_R1',
+            'COMMER_R1',
+            'EDU_R1',
+            'HEALTH_R1',
+            'ART_ACC_R1',
+            'RELIG_R1',
+            'AG_MI_R2',
+            'UTIL_R2',
+            'CONS_R2',
+            'MANU_R2',
+            'COMMER_R2',
+            'EDU_R2',
+            'HEALTH_R2',
+            'ART_ACC_R2',
+            'RELIG_R2',
+            'L1W_R1',
+            'L2W_R1',
+            'L3W_R1',
+            'L4W_R1',
+            'L1W_R2',
+            'L2W_R2',
+            'L3W_R2',
+            'L4W_R2',
+            'KAP',
+            'INVES',
+            'HS1_R1',
+            'HS2_R1',
+            'HS3_R1',
+            'HS1_R2',
+            'HS2_R2',
+            'HS3_R2',
+            'HH1_R1',
+            'HH2_R1',
+            'HH3_R1',
+            'HH4_R1',
+            'HH1_R2',
+            'HH2_R2',
+            'HH3_R2',
+            'HH4_R2',
+            'USSOC1_R1',
+            'USSOC2_R1',
+            'USSOC3_R1',
+            'USSOC4_R1',
+            'USSOC1_R2',
+            'USSOC2_R2',
+            'USSOC3_R2',
+            'USSOC4_R2',
+            'PROTAX',
+            'UTSTX',
+            'COUNTYSTX',
+            'UTITX',
+            'USPIT',
+            'OTX',
+            'CYGF',
+            'COUNTY',
+            'STATE',
+            'FED',
+            'OUT1',
+            'OUT2',
+            'OUT3',
+            'OUT4',
+            'ROW']
+
+        # FACTORS
+        F = ['L1W_R1', 'L2W_R1', 'L3W_R1', 'L4W_R1', 'L1W_R2', 'L2W_R2', 'L3W_R2', 'L4W_R2', 'KAP']
+        FW1 = ['L1W_R1']
+        FW2 = ['L2W_R1']
+        FW3 = ['L3W_R1']
+        FW4 = ['L4W_R1']
+
+        # OUT COMMUTERS
+        CM = ['OUT1', 'OUT2', 'OUT3', 'OUT4']
+        CM1 = ['OUT1']
+        CM2 = ['OUT2']
+        CM3 = ['OUT3']
+        CM4 = ['OUT4']
+
+        # LABOR SUPPLY
+        LT = ['L1W_R1', 'L2W_R1', 'L3W_R1', 'L4W_R1', 'L1W_R2', 'L2W_R2', 'L3W_R2', 'L4W_R2',
+              'OUT1', 'OUT2', 'OUT3', 'OUT4']
+
+        # LABOR WORKING IN GALVESTON
+        L = ['L1W_R1', 'L2W_R1', 'L3W_R1', 'L4W_R1', 'L1W_R2', 'L2W_R2', 'L3W_R2', 'L4W_R2']
+
+        # CAPITAL
+        K = ['KAP']
+
+        # LAND
+        # LA = ['LAND']
+
+        # GOVERNMENTS
+        G = ['USSOC1_R1', 'USSOC2_R1', 'USSOC3_R1', 'USSOC4_R1', 'USSOC1_R2', 'USSOC2_R2', 'USSOC3_R2', 'USSOC4_R2',
+             'PROTAX', 'UTSTX', 'COUNTYSTX', 'UTITX', 'USPIT', 'OTX', 'CYGF', 'COUNTY', 'STATE', 'FED']
+
+        # ENDOGENOUS GOVERNMENTS
+        GN = ['COUNTY', 'STATE', 'FED']
+
+        # LOCAL ENDOGENOUS GOVERNMENTS
+        GNL = ['COUNTY']
+
+        # EXOGENOUS GOVERMENTS
+        GX = ['USSOC1_R1', 'USSOC2_R1', 'USSOC3_R1', 'USSOC4_R1', 'USSOC1_R2', 'USSOC2_R2', 'USSOC3_R2', 'USSOC4_R2',
+              'PROTAX', 'UTSTX', 'COUNTYSTX', 'UTITX', 'USPIT', 'OTX']
+
+        # SALES OR EXCISE TAXES
+        GS = ['UTSTX', 'COUNTYSTX', 'OTX']
+
+        # LOCAL TAX
+        GC = ['UTSTX', 'COUNTYSTX', 'OTX']
+
+        # LAND TAXES
+        GL = ['PROTAX']
+
+        # FACTOR TAXES
+        GF = ['USSOC1_R1', 'USSOC2_R1', 'USSOC3_R1', 'USSOC4_R1', 'USSOC1_R2', 'USSOC2_R2', 'USSOC3_R2',
+              'USSOC4_R2', 'PROTAX']
+
+        # INCOME TAX UNITS
+        GI = ['UTITX', 'USPIT']
+
+        # HOUSEHOLD TAX UNITS
+        GH = ['PROTAX']
+
+        # SOCIAL SECURITY PAYMENT
+        GFUS = ['USSOC1_R1', 'USSOC2_R1', 'USSOC3_R1', 'USSOC4_R1', 'USSOC1_R2', 'USSOC2_R2', 'USSOC3_R2', 'USSOC4_R2']
+
+        # EXOGNOUS TRANSFER PMT
+        GY = ['USSOC1_R1', 'USSOC2_R1', 'USSOC3_R1', 'USSOC4_R1', 'USSOC1_R2', 'USSOC2_R2', 'USSOC3_R2', 'USSOC4_R2',
+              'PROTAX', 'UTSTX', 'COUNTYSTX', 'UTITX', 'USPIT', 'OTX', 'CYGF', 'COUNTY', 'STATE', 'FED']
+
+        # EXOGNOUS TRANSFER PMT
+        GTA = ['USSOC1_R1', 'USSOC2_R1', 'USSOC3_R1', 'USSOC4_R1', 'USSOC1_R2', 'USSOC2_R2', 'USSOC3_R2', 'USSOC4_R2',
+               'PROTAX', 'UTSTX', 'COUNTYSTX', 'UTITX', 'USPIT', 'OTX', 'CYGF', 'COUNTY', 'STATE', 'FED']
+
+        # ENDOGENOUS TRANSFER PMT
+        GT = ['CYGF', 'STATE', 'FED']
+
+        # ALL HOUSEHOLDS IN GALVESTON
+        H = ['HH1_R1', 'HH2_R1', 'HH3_R1', 'HH4_R1', 'HH1_R2', 'HH2_R2', 'HH3_R2', 'HH4_R2']
+
+        # HOUSEHOLD (INCOME i, i = 1,2,3,4)
+        HH1 = ['HH1_R1', 'HH1_R2']
+        HH2 = ['HH2_R1', 'HH2_R2']
+        HH3 = ['HH3_R1', 'HH3_R2']
+        HH4 = ['HH4_R1', 'HH4_R2']
+
+        # I+G SECTORS
+        IG = ['AG_MI_R1', 'UTIL_R1', 'CONS_R1', 'MANU_R1', 'COMMER_R1', 'EDU_R1', 'HEALTH_R1', 'ART_ACC_R1', 'RELIG_R1',
+              'AG_MI_R2', 'UTIL_R2', 'CONS_R2', 'MANU_R2', 'COMMER_R2', 'EDU_R2', 'HEALTH_R2', 'ART_ACC_R2', 'RELIG_R2',
+              'HS1_R1', 'HS2_R1', 'HS3_R1', 'HS1_R2', 'HS2_R2', 'HS3_R2',
+              'COUNTY', 'STATE', 'FED']
+
+        # INDUSTRY SECTORS
+        I = ['AG_MI_R1', 'UTIL_R1', 'CONS_R1', 'MANU_R1', 'COMMER_R1', 'EDU_R1', 'HEALTH_R1', 'ART_ACC_R1', 'RELIG_R1',
+             'AG_MI_R2', 'UTIL_R2', 'CONS_R2', 'MANU_R2', 'COMMER_R2', 'EDU_R2', 'HEALTH_R2', 'ART_ACC_R2', 'RELIG_R2',
+             'HS1_R1', 'HS2_R1', 'HS3_R1', 'HS1_R2', 'HS2_R2', 'HS3_R2']
+
+        # ENDOGENOUS GOVERNMENTS
+        IG2 = ['COUNTY', 'STATE', 'FED']
+
+        # PRODUCTION SECTORS
+        IP = ['AG_MI_R1', 'UTIL_R1', 'CONS_R1', 'MANU_R1', 'COMMER_R1', 'EDU_R1', 'HEALTH_R1', 'ART_ACC_R1', 'RELIG_R1',
+              'AG_MI_R2', 'UTIL_R2', 'CONS_R2', 'MANU_R2', 'COMMER_R2', 'EDU_R2', 'HEALTH_R2', 'ART_ACC_R2', 'RELIG_R2']
+
+        # PRODUCTION GOV
+        FG = ['COUNTY', 'STATE', 'FED']
+
+        # HOUSING SERVICE DEMAND
+        HSD = ['HS1_R1', 'HS2_R1', 'HS3_R1', 'HS1_R2', 'HS2_R2', 'HS3_R2', ]
+
+        # ETA ELASTICITIES
+        ETA = ['ETAL1', 'ETAI1', 'ETALB1', 'ETAPIT', 'ETAPT', 'ETARA', 'NRPG', 'ETAYD', 'ETAU', 'ETAM', 'ETAE', 'ETAY',
+               'ETAOP']
+
+        # ETALANDCAP(ETA) LANDCAP TABLE ELASTICITIES
+        ETALANDCAP = ['ETAL1', 'ETAI1', 'ETALB1']
+
+        # ETAMISCH(ETA) MISCH TABLE ELASTICITIES
+        ETAMISCH = ['ETAPIT', 'ETAPT', 'ETARA', 'NRPG', 'ETAYD', 'ETAU']
+
+        # ETAMISC(ETA) MISC TABLE ELASTICITIES
+        ETAMISC = ['ETAM', 'ETAE', 'ETAY', 'ETAOP']
+
+        # SIMMLOOP
+        SM = ['BASE', 'TODAY', 'simm']
+
+        # EPORT 1 FOR SCALARS
+        R1H = ['GFREV', 'SFREV', 'PIT',
+               'DGF', 'DSF', 'DDRE', 'PDRE', 'SPI', 'COMM', 'COMMO',
+               'GN', 'NKI', 'HH', 'W', 'W1', 'W2', 'W3',
+               'R', 'RL', 'L', 'K', 'HN', 'HW', 'GFSAV', 'LD',
+               'CMO', 'CMI', 'HC', 'SSC', 'LAND', 'LAS']
+
+        # REPORT 2FOR STATUS
+        R2H = ['M-STAT', 'S-STAT']
+
+        # LABELS FOR MODEL STATUS
+        MS = ['OPTIMAL', 'LOCALOP', 'UNBOUND',
+              'INFSBLE', 'INFSLOC', 'INFSINT',
+              'NOOPTML', 'MIPSOLN', 'NOINTGR',
+              'INFSMIP', 'UNUSED', 'UNKNOWN',
+              'NOSOLUT']
+
+        # LABELS FOR SOLVER STATUS
+        SS = ['OK', 'ITERATE', 'RESRCE',
+              'SOLVER', 'EVALUATE', 'NOTKNWN',
+              'NOTUSED', 'PRE-PROC', 'SETUP',
+              'SLVFAIL', 'SLVINTER', 'POST-PROC',
+              'METSYS']
+
+        # ----------------------------------------------------------------
+        # SET ALIASES
+        # ----------------------------------------------------------------
+
+        # ALIAS
+        J = I
+        I1 = I
+        Z1 = Z
+        F1 = F
+        G1 = G
+        G2 = G
+        GI1 = GI
+        GS1 = GS
+        GX1 = GX
+        GN1 = GN
+        GH1 = GH
+        GF1 = GF
+        H1 = H
+        HSD1 = HSD
+        JP = IP
+        JG = IG
+        GT1 = GT
+        GNL1 = GNL
+        L1 = L
+
+        # ----------------------------------------------------------------
+        # PARAMETER DECLARATION
+        # ----------------------------------------------------------------
+
+        # These are data frames with zeros to be filled during calibration.
+        A = pd.DataFrame(index=Z, columns=Z).fillna(0.0)
+        AD = pd.DataFrame(index=Z, columns=Z).fillna(0.0)
+        AG = pd.DataFrame(index=Z, columns=G).fillna(0.0)
+        ALPHA = pd.DataFrame(index=Z, columns=I).fillna(0.0)
+        B = pd.DataFrame(index=I, columns=IG).fillna(0.0)
+        CMOWAGE = pd.Series(index=CM, dtype="float64").fillna(0.0)
+        CMIWAGE = pd.Series(index=L, dtype="float64").fillna(0.0)
+        FCONST = pd.DataFrame(index=F, columns=I).fillna(0.0)
+        DELTA = pd.Series(index=I, dtype="float64").fillna(0.0)
+        GAMMA = pd.Series(index=I, dtype="float64").fillna(0.0)
+        PIT = pd.DataFrame(index=G, columns=H).fillna(0.0)
+        PIT0 = pd.DataFrame(index=G, columns=H).fillna(0.0)
+        PRIVRET = pd.Series(index=H, dtype="float64").fillna(0.0)
+        KFOR = pd.Series(index=K, dtype="float64").fillna(0.0)
+        # LFOR = pd.Series(index=LA, dtype="float64").fillna(0.0)
+        GFOR = pd.Series(index=G, dtype="float64").fillna(0.0)
+        NRPG = pd.Series(index=H, dtype="float64").fillna(0.0)
+        RHO = pd.Series(index=I, dtype="float64").fillna(0.0)
+        SIGMA = pd.Series(index=I, dtype=float).fillna(0.0)
+
+        # TAXES
+        TAUFH = pd.DataFrame(index=G, columns=Z).fillna(0.0)
+        TAUFL = pd.DataFrame(index=G, columns=L).fillna(0.0)
+        TAUFK = pd.DataFrame(index=G, columns=K).fillna(0.0)
+        # TAUFLA = pd.DataFrame(index=G, columns=LA).fillna(0.0)
+        TAUH = pd.DataFrame(index=G, columns=H).fillna(0.0)
+        TAUH0 = pd.DataFrame(index=G, columns=H).fillna(0.0)
+        TAUM = pd.DataFrame(index=G, columns=IG).fillna(0.0)
+        TAUQ = pd.DataFrame(index=G, columns=IG).fillna(0.0)
+        TAUC = pd.DataFrame(index=G, columns=I).fillna(0.0)
+        TAUCH = pd.DataFrame(index=G, columns=HSD).fillna(0.0)
+        TAUV = pd.DataFrame(index=G, columns=I).fillna(0.0)
+        TAUN = pd.DataFrame(index=G, columns=IG).fillna(0.0)
+        TAUX = pd.DataFrame(index=G, columns=IG).fillna(0.0)
+        TAUG = pd.DataFrame(index=G, columns=I).fillna(0.0)
+        TAXS = pd.DataFrame(index=G, columns=GX).fillna(0.0)
+        TAXS1 = pd.Series(index=GNL, dtype="float64").fillna(0.0)
+
+        # ELASTICITIES AND TAX DATA IMPOSED
+        LAMBDA = pd.DataFrame(index=I, columns=I).fillna(0.0)
+        BETA = pd.DataFrame(index=I, columns=H).fillna(0.0)
+        BETAH = pd.DataFrame(index=HSD, columns=H).fillna(0.0)
+        ETAD = pd.Series(index=I, dtype="float64").fillna(0.0)
+        ETAE = pd.Series(index=I, dtype="float64").fillna(0.0)
+        ETAI = pd.Series(index=IG, dtype="float64").fillna(0.0)
+        ETAIX = pd.DataFrame(index=K, columns=IG).fillna(0.0)
+        # ETAL = pd.DataFrame(index=LA, columns=IG).fillna(0.0)
+        ETAL1 = pd.Series(index=IG, dtype="float64").fillna(0.0)
+        ETALB1 = pd.Series(index=IG, dtype="float64").fillna(0.0)
+        ETALB = pd.DataFrame(index=L, columns=IG).fillna(0.0)
+        ETAM = pd.Series(index=I, dtype="float64").fillna(0.0)
+        ETARA = pd.Series(index=H, dtype="float64").fillna(0.0)
+        ETAYD = pd.Series(index=H, dtype="float64").fillna(0.0)
+        ETAU = pd.Series(index=H, dtype="float64").fillna(0.0)
+        ETAYDI = pd.Series(index=H, dtype="float64").fillna(0.0)
+        ETAYDO = pd.Series(index=H, dtype="float64").fillna(0.0)
+        ETAUO = pd.Series(index=H, dtype="float64").fillna(0.0)
+        ETAUI = pd.Series(index=H, dtype="float64").fillna(0.0)
+        ETAPT = pd.Series(index=H, dtype="float64").fillna(0.0)
+        ETAPIT = pd.Series(index=H, dtype="float64").fillna(0.0)
+        ECOMI = pd.Series(index=L, dtype="float64").fillna(0.0)
+        ECOMO = pd.Series(index=CM, dtype="float64").fillna(0.0)
+
+        # VARIABLES FOR INTERMEDIATE CALCULATIONS
+        AGFS = pd.DataFrame(index=Z, columns=G).fillna(0.0)
+        ALPHA1 = pd.DataFrame(index=F, columns=I).fillna(0.0)
+        B1 = pd.DataFrame(index=I, columns=I).fillna(0.0)
+        PRIVRET1 = pd.Series(index=H, dtype="float64").fillna(0.0)
+        OUT = pd.DataFrame(index=G, columns=G).fillna(0.0)
+        KPFOR01 = pd.Series(index=K, dtype="float64").fillna(0.0)
+        TT = pd.Series(index=I, dtype="float64").fillna(0.0)
+
+        # ARRAYS BUILT TO EXPORT RESULTS TO SEPARATE FILE
+        R1 = pd.DataFrame(index=R1H, columns=SM).fillna(0.0)
+        R2 = pd.DataFrame(index=R2H, columns=SM).fillna(0.0)
+
+        # INITIAL VALUES OF ENDOGENOUS VARIABLES
+        CG0 = pd.DataFrame(index=I, columns=G).fillna(0.0)
+        CG0T = pd.DataFrame(index=I, columns=G).fillna(0.0)
+        CH0 = pd.DataFrame(index=I, columns=H).fillna(0.0)
+        CH0T = pd.DataFrame(index=I, columns=H).fillna(0.0)
+        CMI0 = pd.Series(index=L, dtype="float64").fillna(0.0)
+        CMO0 = pd.Series(index=CM, dtype="float64").fillna(0.0)
+        CN0 = pd.Series(index=I, dtype="float64").fillna(0.0)
+        CN0T = pd.Series(index=I, dtype="float64").fillna(0.0)
+        CPI0 = pd.Series(index=H, dtype="float64").fillna(0.0)
+        CPIN0 = pd.Series(index=H, dtype="float64").fillna(0.0)
+        CPIH0 = pd.Series(index=H, dtype="float64").fillna(0.0)
+        CX0 = pd.Series(index=I, dtype="float64").fillna(0.0)
+        D0 = pd.Series(index=I, dtype="float64").fillna(0.0)
+        DD0 = pd.Series(index=Z, dtype="float64").fillna(0.0)
+        DS0 = pd.Series(index=Z, dtype="float64").fillna(0.0)
+        DQ0 = pd.Series(index=Z, dtype="float64").fillna(0.0)
+
+        FD0 = pd.DataFrame(index=F, columns=Z).fillna(0.0)
+        IGT0 = pd.DataFrame(index=G, columns=GX).fillna(0.0)
+        KS0 = pd.DataFrame(index=K, columns=IG).fillna(0.0)
+        # LAS0 = pd.DataFrame(index=LA, columns=IG).fillna(0.0)
+        HH0 = pd.Series(index=H, dtype="float64").fillna(0.0)
+        HN0 = pd.Series(index=H, dtype="float64").fillna(0.0)
+        HW0 = pd.Series(index=H, dtype="float64").fillna(0.0)
+        EXWGEO = pd.Series(index=CM, dtype="float64").fillna(0.0)
+
+        HOUSECOR = pd.DataFrame(index=H, columns=HSD).fillna(0.0)
+        JOBCOR = pd.DataFrame(index=H, columns=L).fillna(0.0)
+        OUTCOR = pd.DataFrame(index=H, columns=CM).fillna(0.0)
+
+        M0 = pd.Series(index=I, dtype="float64").fillna(0.0)
+        M01 = pd.Series(index=Z, dtype="float64").fillna(0.0)
+        MI0 = pd.Series(index=H, dtype="float64").fillna(0.0)
+        MO0 = pd.Series(index=H, dtype="float64").fillna(0.0)
+        N0 = pd.DataFrame(index=K, columns=IG).fillna(0.0)
+        # NKI0
+        KPFOR0 = pd.Series(index=K, dtype="float64").fillna(0.0)
+        # LNFOR0 = pd.Series(index=LA, dtype="float64").fillna(0.0)
+        # LNFOR01 = pd.Series(index=LA, dtype="float64").fillna(0.0)
+
+        GVFOR0 = pd.Series(index=G, dtype="float64").fillna(0.0)
+        P0 = pd.Series(index=IG, dtype="float64").fillna(0.0)
+        PD0 = pd.Series(index=I, dtype="float64").fillna(0.0)
+        PVA0 = pd.Series(index=I, dtype="float64").fillna(0.0)
+        PW0 = pd.Series(index=I, dtype="float64").fillna(0.0)
+        PWM0 = pd.Series(index=I, dtype="float64").fillna(0.0)
+        Q0 = pd.Series(index=Z, dtype="float64").fillna(0.0)
+        Q10 = pd.Series(index=Z, dtype="float64").fillna(0.0)
+        R0 = pd.DataFrame(index=F, columns=Z).fillna(1.0)
+        RA0 = pd.Series(index=Z, dtype="float64").fillna(0.0)
+
+        S0 = pd.Series(index=Z, dtype="float64").fillna(0.0)
+        # SPI0
+        V0 = pd.Series(index=I, dtype="float64").fillna(0.0)
+        V0T = pd.Series(index=I, dtype="float64").fillna(0.0)
+        TP = pd.DataFrame(index=H, columns=G).fillna(0.0)
+        # TAUF0 = Table(G,F,Z)
+        YD0 = pd.Series(index=H, dtype="float64").fillna(0.0)
+        Y0 = pd.Series(index=Z, dtype="float64").fillna(0.0)
+        Y01 = pd.Series(index=H, dtype="float64").fillna(0.0)
+        YT0 = pd.Series(index=G, dtype="float64").fillna(0.0)
+        GCP10 = pd.Series(index=I, dtype="float64").fillna(0.0)
+        # GCP0
+        SD3 = pd.Series(index=GX, dtype="float64").fillna(0.0)
+        DDCX = pd.Series(index=I, dtype="float64").fillna(0.0)
+
+        TESTA1 = pd.DataFrame(index=F, columns=I).fillna(0.0)
+        TESTA2 = pd.DataFrame(index=F, columns=I).fillna(0.0)
+        TESTA3 = pd.DataFrame(index=F, columns=I).fillna(0.0)
+
+        # Lstore(Z,L,t)
+        # DYSTORE = pd.DataFrame(index=Z, columns=T).fillna(0.0)
+        # DSSTORE = pd.DataFrame(index=I, columns=T).fillna(0.0)
+        # HHSTORE = pd.DataFrame(index=H, columns=T).fillna(0.0)
+        # FSTORE(Z,F,T)
+        # CPISTORE = pd.DataFrame(index=H, columns=T).fillna(0.0)
+
+        TPC = pd.DataFrame(index=H, columns=G).fillna(0.0)
+        IGTD = pd.DataFrame(index=G, columns=G1).fillna(0.0)
+        TAUFF = pd.DataFrame(index=G, columns=F).fillna(0.0)
+        IOUT = pd.DataFrame(index=G1, columns=G1).fillna(0.0)
+        LANDCAP = pd.DataFrame(index=IG, columns=ETALANDCAP).fillna(0.0)
+        MISC = pd.DataFrame(index=IG, columns=ETAMISC).fillna(0.0)
+        # MISCH= pd.DataFrame(index=H, columns=ETAMISCH).fillna(0.0)
+
+        # -------------------------------------------------------------------------------------------------------------
+        # SIMPLIFYING TABLES AND DOING AWAY WITH MISC FILES
+        # -------------------------------------------------------------------------------------------------------------
+
+        IGTD.loc[G, G1] = 0
+        IGTD.loc['FED', GFUS] = 1
+        IGTD.loc['FED', 'USPIT'] = 1
+
+        IGTD.loc['CYGF', 'PROTAX'] = 1
+        IGTD.loc['CYGF', 'COUNTYSTX'] = 1
+
+        IGTD.loc['CYGF', 'OTX'] = 1
+
+        IGTD.loc['STATE', 'UTSTX'] = 1
+        IGTD.loc['STATE', 'UTITX'] = 1
+
+        TPC.loc[H, G] = 0
+        TPC.loc[H, GFUS] = 1
+
+        TAUFF.loc[G, F] = 0
+        TAUFF.loc['USSOC1_R1', 'L1W_R1'] = 1
+        TAUFF.loc['USSOC2_R1', 'L2W_R1'] = 1
+        TAUFF.loc['USSOC3_R1', 'L3W_R1'] = 1
+        TAUFF.loc['USSOC4_R1', 'L4W_R1'] = 1
+
+        TAUFF.loc['USSOC1_R2', 'L1W_R2'] = 1
+        TAUFF.loc['USSOC2_R2', 'L2W_R2'] = 1
+        TAUFF.loc['USSOC3_R2', 'L3W_R2'] = 1
+        TAUFF.loc['USSOC4_R2', 'L4W_R2'] = 1
+
+        TAUFF.loc['PROTAX', 'KAP'] = 1
+
+        for label in G1:
+            IOUT.loc[label, label] = 0
+        # IOUT.loc['OTHER_LOCAL', 'CYGF'] = 1
+        IOUT.loc['COUNTY', 'CYGF'] = 1
+
+        LANDCAP.loc[IG, ETALANDCAP] = 1
+
+        MISCH.loc[H, ETAMISCH] = 0
+
+        MISCH.loc[HH1, 'ETAPIT'] = -0.5
+        MISCH.loc[HH2, 'ETAPIT'] = -0.15
+        MISCH.loc[HH3, 'ETAPIT'] = -0.15
+        MISCH.loc[HH4, 'ETAPIT'] = -0.15
+
+        MISCH.loc[HH1, 'ETAPT'] = -0.5
+        MISCH.loc[H, 'NRPG'] = 1
+        MISCH.loc[H, 'ETARA'] = 1
+
+        MISCH.loc[HH1, 'ETAYD'] = 1.3
+        MISCH.loc[HH2, 'ETAYD'] = 1.3
+        MISCH.loc[HH3, 'ETAYD'] = 1.6
+        MISCH.loc[HH4, 'ETAYD'] = 1.6
+
+        MISCH.loc[HH1, 'ETAU'] = -0.8
+        MISCH.loc[HH2, 'ETAU'] = -0.8
+        MISCH.loc[HH3, 'ETAU'] = -0.8
+        MISCH.loc[HH4, 'ETAU'] = -0.6
+
+        MISC.loc[IG, ETAMISC] = 0
+        MISC.loc[IP, 'ETAM'] = 1
+        MISC.loc[IP, 'ETAE'] = -1.65
+        MISC.loc[I, 'ETAY'] = 1
+        MISC.loc[I, 'ETAOP'] = -1
+
+        # -------------------------------------------------------------------------------------------------------------
+        # CALCULATIONS OF PARAMETERS AND INITIAL VALUES
+        # -------------------------------------------------------------------------------------------------------------
+
+        # CALCULATE COLUMN AND ROW TOTALS OF SAM TO COMPARE FOR BALANCE
+        Q10.loc[Z] = SAM.loc[Z].sum(0)  # Column Totals of SAM table
+
+        Q0.loc[Z] = SAM.loc[Z].sum(1)  # Row Totals of SAM table
+
+        DQ0.loc[Z] = Q10.loc[Z] - Q0.loc[Z]  # difference of SAM row and coloumn totals
+
+        B1.loc[I, I] = SAM.loc[I, I]
+
+        # READ IN ELASTICITY PARAMETERS FROM MISC.PRN
+        OUT.loc[G, G] = IOUT.loc[G, G]
+
+        BETA.loc[I, H] = MISC.loc[I, 'ETAY']
+
+        BETAH.loc[HSD, H] = MISC.loc[HSD, 'ETAY']
+
+        LAMBDA.loc[I, I] = MISC.loc[I, 'ETAOP']
+
+        ETAE.loc[I] = MISC.loc[I, 'ETAE']
+
+        ETAM.loc[I] = MISC.loc[I, 'ETAM']
+
+        ETARA.loc[H] = MISCH.loc[H, 'ETARA']
+
+        ETAPIT.loc[H] = MISCH.loc[H, 'ETAPIT']
+
+        ETAPT.loc[H] = MISCH.loc[H, 'ETAPT']
+
+        ETAYD.loc[H] = MISCH.loc[H, 'ETAYD']
+
+        NRPG.loc[H] = MISCH.loc[H, 'NRPG']
+
+        ETAU.loc[H] = MISCH.loc[H, 'ETAU']
+
+        ECOMI.loc[L] = 1.0
+        ECOMO.loc[CM] = 1.0
+        EXWGEO.loc[CM] = 1.0
+
+        ETAI.loc[IG] = LANDCAP.loc[IG, 'ETAI1']
+        ETAL1.loc[IG] = LANDCAP.loc[IG, 'ETAL1']
+
+        ETAIX.loc['KAP', IG] = ETAI.loc[IG]
+        # ETAL.loc['LAND',IG] = ETAL1.loc[IG]
+
+        # CALCULATE TAX RATES FROM SAM INFORMATION
+        TAUQ_1 = SAM.loc[GS, I]
+        TAUQ_2 = SAM.loc[I, I].sum(1)
+        TAUQ_3 = SAM.loc[I, H].sum(1)
+        TAUQ_4 = SAM.loc[I, ['INVES']].sum(1)
+        TAUQ_5 = SAM.loc[I, G].sum(1)
+        TAUQ_6 = SAM.loc[I, ['ROW']].sum(1)
+        TAUQ_7 = SAM.loc[GS1, I].sum(0)
+
+        TAUQ.loc[GS, I] = TAUQ_1 / (TAUQ_2 + TAUQ_3 + TAUQ_4 + TAUQ_5 + TAUQ_6 - TAUQ_7)
+
+        TAUC.loc[GS, I] = TAUQ.loc[GS, I]
+        TAUV.loc[GS, I] = TAUQ.loc[GS, I]
+        TAUN.loc[GS, I] = TAUQ.loc[GS, I]
+        TAUG.loc[GS, I] = TAUQ.loc[GS, I]
+        TAUX.loc[GS, I] = TAUQ.loc[GS, I]
+
+        # FACTOR TAX RATES
+        TAUFarray = [[[0 for k in range(len(Z))] for j in range(len(F))] for i in range(len(G))]
+        for i in range(len(GF)):
+            for j in range(len(F)):
+                for k in range(len(I)):
+                    if SAM.loc[F[j], I[k]] != 0 and TAUFF.loc[GF[i], F[j]] != 0:
+                        TAUFarray[G.index(GF[i])][j][Z.index(I[k])] = SAM.loc[GF[i], I[k]] / SAM.loc[F[j], I[k]]
+
+        for i in range(len(GF)):
+            for j in range(len(F)):
+                for k in range(len(G)):
+                    if SAM.loc[F[j], G[k]] != 0 and TAUFF.loc[GF[i], F[j]] != 0:
+                        TAUFarray[G.index(GF[i])][j][Z.index(G[k])] = SAM.loc[GF[i], G[k]] / SAM.loc[F[j], G[k]]
+
+        TAUFX_SUM_array = [[0 for j in range(len(Z))] for i in range(len(F))]
+        for i in range(len(F)):
+            for j in range(len(Z)):
+                tmp = 0
+                for k in range(len(G)):
+                    tmp += TAUFarray[k][i][j]
+                TAUFX_SUM_array[i][j] = tmp
+
+        # TAUFX summed over GX
+        TAUFX_SUM = pd.DataFrame(TAUFX_SUM_array, index=F, columns=Z).fillna(0.0)
+
+        TAUFX_GF = pd.DataFrame(TAUFX_SUM_array, index=F, columns=Z).fillna(0.0)
+
+        TAUFXgx = {}
+        for i in range(len(GX)):
+            TAUFXgx[GX[i]] = pd.DataFrame(TAUFarray[i], index=F, columns=Z).fillna(0.0)
+
+        # SS TAX RATES
+        # TAUFH(GF,F)$(TAUFF(GF,F)) =SAM(GF,F) / SUM(Z, SAM(Z,F));
+        for i in GF:
+            for j in F:
+                if TAUFF.loc[i, j] != 0:
+                    #            TAUFH.set_value(i, j, SAM.loc[i, j] / SAM.loc[Z, F].sum(0).loc[j])
+                    TAUFH.at[i, j] = SAM.at[i, j] / SAM.loc[Z, F].sum(0).at[j]
+
+        # TAUFH(GFUS,L)  =SAM(GFUS,L) / SUM(IG, SAM(L,IG));
+        for i in GFUS:
+            for j in L:
+                if TAUFF.loc[i, j] != 0:
+                    #            TAUFH.set_value(i, j, SAM.loc[i, j] / SAM.loc[L, IG].sum(1).loc[j])
+                    TAUFH.at[i, j] = SAM.at[i, j] / SAM.loc[L, IG].sum(1).at[j]
+
+        # EMPLOYEE PORTION OF FACTOR TAXES
+        TAUFL.loc[GF, L] = SAM.loc[GF, L] / SAM.loc[Z, L].sum(0)
+        TAUFK.loc[GF, K] = SAM.loc[GF, K] / SAM.loc[Z, K].sum(0)
+
+        # SHARES OF ENDOGENOUS GOVERNMENTS TRANFERS TO REVENUE
+        TAXS.loc[G, GX] = SAM.loc[G, GX] / SAM.loc[G, GX].sum(0)
+        TAXS1.loc[GNL] = SAM.loc[GNL, ['CYGF']].sum(1) / SAM.loc[GNL, ['CYGF']].sum(1).sum(0)
+
+        # SET INITIAL INTER GOVERNMENTAL TRANSFERS
+        IGT0.loc[G, GX] = SAM.loc[G, GX]
+
+        PW0.loc[I] = 1.0
+        PWM0.loc[I] = 1.0
+        P0.loc[I] = 1.0
+        PD0.loc[I] = 1.0
+        CPI0.loc[H] = 1.0
+        CPIN0.loc[H] = 1.0
+        CPIH0.loc[H] = 1.0
+
+        # HOUSEHOLD TRANSFER PAYMENTS AND PERSONAL INCOME TAXES
+
+        # TOTAL HH
+        HH0.loc[H] = MISCH.loc[H, 'HH0']
+
+        # TOTAL WORKING HH
+        HW0.loc[H] = MISCH.loc[H, 'HW0']
+
+        # NON WORKING HH
+        HN0.loc[H] = HH0.loc[H] - HW0.loc[H]
+
+        # NOMINAL GOVERNMENT SS PAYMENTS
+        TP.loc[H, G] = SAM.loc[H, G].div(HH0.loc[H], axis='index').fillna(0.0)
+
+        # FACTOR RENTALS
+        JOBCOR.loc[H, L] = JOBCR.loc[H, L]
+
+        # COMMUTERS
+        CMO0.loc[CM] = (OUTCR.loc[H, CM].mul(HW0.loc[H], axis='index')).sum(0)
+
+        # RENTAL RATE FOR FACTORS
+        R0.loc[F, Z] = 1.0
+
+        R0.loc[F, IG] = (SAM.loc[F, IG] / EMPLOY.loc[IG, F].T).fillna(1.0)
+
+        # REAL FACTOR DEMAND
+        FD0.loc[F, IG] = EMPLOY.loc[IG, F].T
+        KS0.loc[K, IG] = FD0.loc[K, IG]
+
+        # SHARES FOUND IN THE SOCIAL ACCOUNTING MATRIX DATA
+        A.loc[Z, Z] = SAM.loc[Z, Z].div(Q0.loc[Z], axis='columns')
+
+        # AGFS: LABOR PAYMENTS BY G SECTOR + USSOC PAYMENTS BY LABOR (GROSS LABOR PAYMENTS)
+        AGFS.loc['L1W_R1', G] = SAM.loc['L1W_R1', G] + SAM.loc['USSOC1_R1', G]
+        AGFS.loc['L2W_R1', G] = SAM.loc['L2W_R1', G] + SAM.loc['USSOC2_R1', G]
+        AGFS.loc['L3W_R1', G] = SAM.loc['L3W_R1', G] + SAM.loc['USSOC3_R1', G]
+        AGFS.loc['L4W_R1', G] = SAM.loc['L4W_R1', G] + SAM.loc['USSOC4_R1', G]
+
+        AGFS.loc['L1W_R2', G] = SAM.loc['L1W_R2', G] + SAM.loc['USSOC1_R2', G]
+        AGFS.loc['L2W_R2', G] = SAM.loc['L2W_R2', G] + SAM.loc['USSOC2_R2', G]
+        AGFS.loc['L3W_R2', G] = SAM.loc['L3W_R2', G] + SAM.loc['USSOC3_R2', G]
+        AGFS.loc['L4W_R2', G] = SAM.loc['L4W_R2', G] + SAM.loc['USSOC4_R2', G]
+
+        # AG - GOVERNMENT SPENDING SHARES OF NET INCOME
+        AG_1 = SAM.loc[I, G]
+        AG_2 = SAM.loc[I, G].sum(0) + SAM.loc[F, G].sum(0) + SAM.loc[GF, G].sum(0)
+        AG.loc[I, G] = AG_1 / AG_2
+
+        AG_1 = SAM.loc[F, G]
+        AG.loc[F, G] = AG_1 / AG_2
+
+        AG_1 = AGFS.loc[L, G]
+        AG.loc[L, G] = AG_1 / AG_2
+        AG = AG.fillna(0.0)
+
+        # TRADE INTERMEDIATES CONSUMPTION INVESTMENT INITIAL LEVELS
+
+        # REAL EXPORT CONSUMPTION
+        CX0.loc[I] = SAM.loc[I, ["ROW"]].div(P0.loc[I], axis='index').div(1.0 + TAUQ.loc[GS, I].sum(0),
+                                                                          axis='index').sum(1)
+
+        # REAL IMPORTS
+        # M01.loc[I] = SAM.loc[["ROW"], I].sum(0) / PWM0[I].T
+        M0.loc[IP] = SAM.loc[IP, Z].sum(1) - (B1.loc[I, IP].sum(0) + SAM.loc[F, IP].sum(0) + SAM.loc[G, IP].sum(0))
+        M0.loc[I] = (M0[I] / PWM0[I])
+        M0 = M0.fillna(0.0)
+
+        # REAL INTERMEDIATE DEMAND
+        V0.loc[I] = SAM.loc[I, I].sum(1) / P0.loc[I] / (1.0 + TAUQ.loc[GS, I].sum(0))
+        V0T.loc[I] = SAM.loc[I, I].sum(1) / P0.loc[I]
+
+        # REAL PRIVATE CONSUMPTION
+        CH0.loc[I, H] = SAM.loc[I, H].div(P0.loc[I], axis='index').div(1.0 + TAUQ.loc[GS, I].sum(0), axis='index')
+
+        CH0T.loc[I, H] = SAM.loc[I, H].div(P0[I], axis='index')
+
+        CG0.loc[I, GN] = SAM.loc[I, GN].div(P0.loc[I], axis='index').div((1.0 + TAUQ.loc[GS, I].sum(0)), axis='index')
+
+        CG0T.loc[I, GN] = SAM.loc[I, GN].div(P0.loc[I], axis='index')
+
+        DEPR = float((SAM.loc[IG, ["INVES"]].sum(0)) / (KS0.loc[K, IG].sum(1).sum(0)))
+
+        N0.loc[K, IG] = KS0.loc[K, IG] * DEPR
+
+        # INVESTMENT BY SECTOR OF SOURCE
+        CN0.loc[I] = 0.0
+
+        B.loc[I, IG] = BB.loc[I, IG].fillna(0.0)
+
+        CN0.loc[I] = (B.loc[I, IG].mul(N0.loc[K, IG].sum(0), axis='columns')).sum(1).div(P0.loc[I], axis='index').div(
+            1.0 + TAUN.loc[GS, I].sum(0), axis='index').transpose()
+
+        CN0T.loc[I] = (B.loc[I, IG].mul(N0.loc[K, IG].sum(0), axis='columns')).sum(1).div(P0.loc[I], axis='index')
+
+        DD0.loc[I] = CH0.loc[I, H].sum(1) + CG0.loc[I, G].sum(1) + CN0.loc[I] + V0.loc[I]
+
+        D0.loc[I] = 1.0 - (M0.loc[I] / DD0.loc[I])
+
+        # CORRECT IMPORT ELASTICITY TO DOMESTIC SHARE ELASTICITY
+        ETAD.loc[I] = -1.0 * ETAM.loc[I] * M0.loc[I] / (DD0.loc[I] * D0.loc[I])
+
+        # PRODUCTION DATA
+        DS0.loc[I] = DD0.loc[I] + CX0.loc[I] - M0.loc[I];
+
+        AD.loc[I, I] = SAM.loc[I, I].div(P0.loc[I], axis='index').div(1.0 + TAUQ.loc[GS, I].sum(0), axis='index') / \
+                       DS0.loc[I]
+
+        PVA0.loc[I] = PD0.loc[I] - (
+            (AD.loc[I, I].mul(P0.loc[I], axis='index').mul(1.0 + TAUQ.loc[GS, I].sum(0).T, axis='index')).sum(0).T)
+
+        # AVERAGE RENTAL RATES FOR FACTORS (NORMALIZED)
+        RA0.loc[F] = 1.0
+
+        '''
+        # CALIBRATION OF PRODUCTION EXPONENTS FOR COBB DOUGLAS
+        a = pd.Series(index=I, dtype="float64").fillna(0.0)
+        a = SAM.loc[USSOCL, I].append(a, ignore_index=True).append(SAM.loc[GL, I]) # labor, land, capital
+        a.index = F
+        
+        ALPHA.loc[F, I] = (SAM.loc[F, I] + a.loc[F, I]) / (SAM.loc[F, I].sum(0) + SAM.loc[GF, I].sum(0))
+        ALPHA.loc[F, I] = ALPHA.loc[F, I] / ALPHA.loc[F, I].sum(0)
+        '''
+
+        # DELTA(I) = DS0(I)/ (PROD(F$ALPHA(F,I),FD0(F,I)**ALPHA(F,I)));
+        # replace takes care of multiplying by zeros, by changing zeros to ones.
+        # DELTA.loc[I] = DS0.loc[I] / (FD0.loc[F, I] ** ALPHA.loc[F, I]).replace({0: 1}).product(0)
+
+        SIGMA.loc[I] = 0.67
+        RHO.loc[I] = (SIGMA.loc[I] - 1) / SIGMA.loc[I]
+
+        TESTA1.loc[F, I] = 0.0
+
+        # TESTA1(F,I)      = R0(F,I) * RA0(F)*(1 + SUM(GF,TAUFX(GF,F,I)))*((FD0(F,I))**(1-RHO(I)));
+        TESTA1.loc[F, I] = R0.loc[F, I].mul(RA0.loc[F], axis='index') * (1 + TAUFX_SUM.loc[F, I]) * (
+                    FD0.loc[F, I] ** (1 - RHO.loc[I]))
+
+        # TESTA2(F,I)      = SUM(F1, R0(F1,I) * RA0(F1) * (1 + SUM(GF,TAUFX(GF,F1,I)))*((FD0(F1,I))**(1-RHO(I))));
+        TESTA2_pre = (R0.loc[F, I].mul(RA0.loc[F], axis='index') * (1 + TAUFX_SUM.loc[F, I]) * (
+                    FD0.loc[F, I] ** (1 - RHO.loc[I]))).sum(0)
+        TESTA2.loc[F, I] = 1
+        TESTA2.loc[F, I] = TESTA2.loc[F, I] * TESTA2_pre
+
+        TESTA3.loc[F, I] = TESTA1.loc[F, I] / TESTA2.loc[F, I]
+
+        ALPHA.loc[F, I] = TESTA3.loc[F, I]
+
+        # add THETA to make FDEQ work
+        THETA = ALPHA.where(ALPHA == 0, 1)
+
+        ATEST1 = pd.Series(index=I, dtype=float).fillna(0.0)
+        ATEST1.loc[I] = ALPHA.loc[F, I].sum(0)
+
+        # GAMMA(I)    = DS0(I) / ( SUM(F, ALPHA(F,I) *  FD0(F,I)) ** ( RHO(I) ) ) ** ( 1 / RHO(I) );
+        GAMMA.loc[I] = DS0.loc[I] / (((ALPHA.loc[F, I] * FD0.loc[F, I] ** (RHO.loc[I])).sum(0)) ** (1 / RHO.loc[I]))
+
+        # OTHER DATA
+        # HH INVESTMENT INCOME FROM ROW
+        PRIVRET.loc[H] = SAM.loc[Z, H].sum(0) - (SAM.loc[H, F].sum(1) + SAM.loc[H, CM].sum(1) + SAM.loc[H, GX].sum(1))
+
+        PRIVRET.loc[H] = PRIVRET.loc[H] / HH0.loc[H]
+
+        # TOTAL OUTPUT
+        Y0.loc[F] = SAM.loc[F, IG].sum(1)
+
+        KPFOR01.loc[K] = SAM.loc[K, ["ROW"]].sum(1)
+
+        # NOMINAL CAPITAL OUTFLOW
+        KPFOR0[K] = SAM.loc[Z, K].sum(0).T - SAM.loc[K, IG].sum(1)
+
+        # NOMINAL GOVERNMENT OUTFLOWS
+        GVFOR0.loc[G] = SAM.loc[G, ["ROW"]].sum(1)
+
+        '''
+        GVFOR0.loc[GT] = SAM.loc[Z, GT].sum(0) - (
+                SAM.loc[GT, I].sum(1)
+                + SAM.loc[GT, F].sum(1)
+                + SAM.loc[GT, H].sum(1)
+                + SAM.loc[GT, G].sum(1)
+                )
+        '''
+        # original:
+        # A.loc[H, L] = SAM.loc[H, L].div(HW0.loc[H], axis='index') / (
+        #        Y0.loc[L] + SAM.loc[Z, L].sum(0) - SAM.loc[L, IG].sum(1)) * (1.0 - TAUFL.loc[G, L].sum(0))
+        A.loc[H, L] = SAM.loc[H, L].div(HW0.loc[H], axis='index') / (Y0.loc[L] + SAM.loc[L, ["ROW"]].sum(1)) * (
+                    1.0 - TAUFL.loc[G, L].sum(0))
+
+        A.loc[H, K] = SAM.loc[H, K].div(HW0.loc[H], axis='index') / (
+                    Y0.loc[K] + SAM.loc[Z, K].sum(0) - SAM.loc[K, IG].sum(1))
+
+        # HH TAXES OTHER THAN PIT
+        TAUH.loc[GH, H] = SAM.loc[GH, H].div(HH0.loc[H], axis='columns')
+
+        TAUH0.loc[GH, H] = SAM.loc[GH, H] / HH0.loc[H]
+
+        S0.loc[H] = SAM.loc[["INVES"], H].T.sum(1)
+
+        YD0.loc[H] = SAM.loc[I, H].sum(0).T + S0.loc[H]
+
+        Y0.loc[G] = SAM.loc[G, Z].sum(1) - SAM.loc[G, ["ROW"]].sum(1)
+
+        S0.loc[G] = SAM.loc[["INVES"], G].sum(0)
+
+        CMI0.loc[L] = EMPLOY.loc[IG, L].sum(0) - (JOBCOR.loc[H, L].mul(HW0.loc[H], axis='index')).sum(0)
+
+        # AVERAGE WAGE FLOWING INTO
+        CMOWAGE.loc[CM] = SAM.loc[CM, ["ROW"]].sum(1).div(CMO0.loc[CM], axis='index').fillna(0.0)
+
+        # AVERAGE WAGES FLOWING OUT OF Robeson
+        CMIWAGE.loc[L] = (-1) * (SAM.loc[L, ["ROW"]].sum(1).div(CMI0.loc[L], axis='index').fillna(0.0))
+
+        # PROPORTION OF CAPITAL INCOME OUTFLOW
+        KFOR.loc[K] = KPFOR0.loc[K] / SAM.loc[["KAP"], IG].sum(1)
+
+        # PROPORTION OF GOVERNMENT INCOME OUTFLOW
+        GFOR.loc[G] = GVFOR0.loc[G] / Y0.loc[G]
+
+        A.loc[H, CM] = SAM.loc[H, CM] / SAM.loc[Z, CM].sum(0)
+
+        # NOMINAL NET CAPITAL INFLOW
+        NKI0 = (M0.loc[I] * PWM0.loc[I]).sum(0) - (CX0.loc[I] * PD0.loc[I]).sum(0) - \
+               (PRIVRET.loc[H] * HH0.loc[H]).sum(0) - KPFOR0.loc[K].sum(0) - GVFOR0.loc[G].sum(0) - \
+               (CMOWAGE.loc[CM] * CMO0.loc[CM]).sum(0) + \
+               (CMIWAGE.loc[L] * CMI0.loc[L]).sum(0)
+
+        # REAL HH NET INCOME
+        Y0.loc[H] = (A.loc[H, L].mul(HW0[H], axis='index').div(A.loc[H, L].mul(HW0[H], axis='index').sum(0),
+                                                               axis='columns') \
+                     .mul(Y0.loc[L] - (CMIWAGE.loc[L] * CMI0.loc[L]), axis='columns').mul(1.0 - TAUFL.loc[G, L].sum(0),
+                                                                                          axis='columns')).sum(1) \
+                    + (A.loc[H, CM].mul((CMOWAGE.loc[CM] * CMO0.loc[CM]), axis='columns')).sum(1) \
+                    + (A.loc[H, K].mul(HW0[H], axis='index') / A.loc[H, K].mul(HW0[H], axis='index').sum(0) \
+                       * (Y0[K] * (1.0 - TAUFK.loc[G, K].sum(0)) + KPFOR0.loc[K])).sum(1)
+
+        # PERSONAL INCOME OBJECTIVE FUNCTION
+        SPI0 = (Y0.loc[H].sum(0) +
+                TP.loc[H, G].mul(HH0.loc[H], axis='index').sum(1).sum(0) +
+                (PRIVRET[H] * HH0[H]).sum(0))
+
+        # PERSONAL INCOME TAX
+        PIT.loc[GI, H] = SAM.loc[GI, H].div(Y0[H], axis='columns')
+        PIT0.loc[GI, H] = SAM.loc[GI, H].div(Y0[H], axis='columns')
+
+        MI0.loc[H] = HH0.loc[H] * 0.04
+
+        MO0.loc[H] = HH0.loc[H] * 0.04
+
+        GCP0 = CH0.loc[I, H].sum(1).sum(0) + CN0.loc[I].sum(0) + CG0.loc[I, GN].sum(1).sum(0) + CX0.loc[I].sum(0) - \
+               M0.loc[I].sum(0)
+
+        GCP10.loc[I] = CH0.loc[I, H].sum(1) + CN0.loc[I] + CG0.loc[I, GN].sum(1) + CX0.loc[I] - M0.loc[I]
+
+        ###########################################
+        # VARIABLE DECLARATION
+        ###########################################
+
+        vars = VarContainer()
+
+        # PUBLIC CONSUMPTION
+        CG = vars.add('CG', rows=I, cols=G)
+
+        # PRIVATE CONSUMPTION
+        CH = vars.add('CH', rows=I, cols=H)
+
+        # GROSS INVESTMENT BY SECTOR OF SOURCE
+        CN = vars.add('CN', rows=I)
+
+        # CONSUMER PRICE INDEX
+        CPI = vars.add('CPI', rows=H)
+
+        # NONHOUSING CONSUMER PRICE INDEX
+        CPIN = vars.add('CPIN', rows=H)
+
+        # NONHOUSING CONSUMER PRICE INDEX
+        CPIH = vars.add('CPIH', rows=H)
+
+        # EXPORT DEMAND
+        CX = vars.add('CX', rows=I)
+
+        # DOMESTIC SHARE OF DOMESTIC DEMAND
+        D = vars.add('D', rows=I)
+
+        # DOMESTIC DEMAND
+        DD = vars.add('DD', rows=I)
+
+        # DOMESTIC SUPPLY
+        DS = vars.add('DS', rows=I)
+
+        # SECTORAL FACTOR DEMAND
+        FD = vars.add('FD', rows=F, cols=Z)
+
+        # GROSS AGGREGATE CITY PRODUCT
+        GCP = vars.add('GCP')
+
+        # GROSS CITY PRODUCT BY SECTOR
+        GCP1 = vars.add('GCP1', rows=I)
+
+        # NUMBER OF HOUSEHOLDS
+        HH = vars.add('HH', rows=H)
+
+        # NUMBER OF NONWORKING HOUSEHOLDS
+        HN = vars.add('HN', rows=H)
+
+        # NUMBER OF WORKING HOUSEHOLDS
+        HW = vars.add('HW', rows=H)
+
+        # INTER GOVERNMENTAL TRANSFERS
+        IGT = vars.add('IGT', rows=G, cols=GX)
+
+        # CAPITAL FLOW
+        KS = vars.add('KS', rows=K, cols=IG)
+
+        # LAND FLOW
+        # LAS = vars.add('LAS', rows=LA, cols=IG)
+
+        # IMPORTS
+        M = vars.add('M', rows=I)
+
+        # GROSS INVESTMENT BY SECTOR OF DESTINATION
+        N = vars.add('N', rows=K, cols=IG)
+
+        # NET CAPITAL INFLOW
+        NKI = vars.add('NKI')
+
+        # LAND OUTFLOW
+        # LNFOR = vars.add('LNFOR', rows=LA)
+
+        # CAPITAL OUTFLOW
+        KPFOR = vars.add('KPFOR', rows=K)
+
+        # GOVT OUTFLOW
+        GVFOR = vars.add('GVFOR', rows=G)
+
+        # AGGREGATE DOMESTIC PRICE PAID BY PURCHASERS
+        P = vars.add('P', rows=I)
+
+        # DOMESTIC PRICE RECEIVED BY SUPPLIERS
+        PD = vars.add('PD', rows=I)
+
+        # VALUE ADDED PRICE
+        PVA = vars.add('PVA', rows=I)
+
+        # ECONOMY WIDE SCALAR RENTAL RATES OF FACTORS
+        RA = vars.add('RA', rows=F)
+
+        # SECTORAL RENTAL RATES
+        R = vars.add('R', rows=F, cols=Z)
+
+        # SAVINGS
+        S = vars.add('S', rows=Z)
+
+        # PERSONAL INCOME (OBJECTIVE FUNCTION)
+        SPI = vars.add('SPI')
+
+        # INTERMEDIATE GOODS
+        V = vars.add('V', rows=I)
+
+        # GROSS INCOMES
+        Y = vars.add('Y', rows=Z)
+
+        # AFTER TAX TOTAL HOUSEHOLD INCOMES
+        YD = vars.add('YD', rows=H)
+
+        # OUT COMMUTERS
+        CMO = vars.add('CMO', rows=CM)
+
+        # IN COMMUTERS
+        CMI = vars.add('CMI', rows=L)
+
+        ###########################################
+        # INITIALIZE VARIABLES FOR SOLVER
+        ###########################################
+
+        vars.init('CG', CG0.loc[I, G])
+        vars.init('CH', CH0.loc[I, H])
+        vars.init('CN', CN0.loc[I])
+        vars.init('CPI', CPI0.loc[H])
+        vars.init('CPIN', CPIN0.loc[H])
+        vars.init('CPIH', CPIH0.loc[H])
+        vars.init('CX', CX0.loc[I])
+        vars.init('D', D0.loc[I])
+        vars.init('DD', DD0.loc[I])
+        vars.init('DS', DS0.loc[I])
+        vars.init('FD', FD0.loc[F, Z])
+        vars.init('GCP', GCP0)
+        vars.init('GCP1', GCP10.loc[I])
+        vars.init('HH', HH0.loc[H])
+        vars.init('HN', HN0.loc[H])
+        vars.init('HW', HW0.loc[H])
+        vars.init('IGT', IGT0.loc[G, GX])
+        vars.init('KS', KS0.loc[K, IG])
+        # vars.init('LAS', LAS0.loc[LA, IG])
+        vars.init('M', M0.loc[I])
+        vars.init('N', N0.loc[K, IG])
+        vars.init('NKI', NKI0)
+        # vars.init('LNFOR', LNFOR0.loc[LA])
+        vars.init('KPFOR', KPFOR0.loc[K])
+        vars.init('GVFOR', GVFOR0.loc[G])
+        vars.init('P', P0.loc[I])
+        vars.init('PD', PD0.loc[I])
+        vars.init('PVA', PVA0.loc[I])
+        vars.init('RA', RA0.loc[F])
+        vars.init('R', R0.loc[F, Z])
+        vars.init('S', S0.loc[Z])
+        vars.init('SPI', SPI0)
+        vars.init('V', V0.loc[I])
+        vars.init('Y', Y0.loc[Z])
+        vars.init('YD', YD0.loc[H])
+        vars.init('CMO', CMO0.loc[CM])
+        vars.init('CMI', CMI0.loc[L])
+
+        # -------------------------------------------------------------------------------------------------------------
+        # DEFINE BOUNDS FOR VARIABLES
+        # -------------------------------------------------------------------------------------------------------------
+
+        vars.lo('P', vars.get('P') / 1000)
+        vars.up('P', vars.get('P') * 1000)
+        vars.lo('PD', vars.get('PD') / 1000)
+        vars.up('PD', vars.get('PD') * 1000)
+        vars.lo('PVA', vars.get('PVA') / 1000)
+        vars.up('PVA', vars.get('PVA') * 1000)
+        vars.lo('RA', vars.get('RA') / 1000)
+        vars.up('RA', vars.get('RA') * 1000)
+        vars.lo('CPI', vars.get('CPI') / 1000)
+        vars.up('CPI', vars.get('CPI') * 1000)
+        vars.lo('DS', vars.get('DS') / 1000)
+        vars.up('DS', vars.get('DS') * 1000)
+        vars.lo('DD', vars.get('DD') / 1000)
+        vars.up('DD', vars.get('DD') * 1000)
+        vars.lo('D', vars.get('D') / 1000)
+        vars.up('D', vars.get('D') * 1000)
+        vars.lo('V', vars.get('V') / 1000)
+        vars.up('V', vars.get('V') * 1000)
+        vars.lo('FD', vars.get('FD') / 1000)
+        vars.up('FD', vars.get('FD') * 1000)
+        vars.lo('HH', vars.get('HH') / 1000)
+        vars.up('HH', vars.get('HH') * 1000)
+        vars.lo('HW', vars.get('HW') / 1000)
+        vars.up('HW', vars.get('HW') * 1000)
+        vars.lo('HN', vars.get('HN') / 1000)
+        vars.up('HN', vars.get('HN') * 1000)
+        vars.lo('KS', vars.get('KS') / 1000)
+        vars.up('KS', vars.get('KS') * 1000)
+        # vars.lo('LAS',vars.get('LAS')/1000); vars.up('LAS',vars.get('LAS')*1000)
+        vars.lo('M', vars.get('M') / 1000)
+        vars.up('M', vars.get('M') * 1000)
+        vars.lo('Y', vars.get('Y') / 1000)
+        vars.up('Y', vars.get('Y') * 1000)
+        vars.lo('YD', vars.get('YD') / 1000)
+        vars.up('YD', vars.get('YD') * 1000)
+        vars.lo('CH', vars.get('CH') / 1000)
+        vars.up('CH', vars.get('CH') * 1000)
+        vars.lo('CG', vars.get('CG') / 1000)
+        vars.up('CG', vars.get('CG') * 1000)
+        vars.lo('CX', vars.get('CX') / 1000)
+        vars.up('CX', vars.get('CX') * 1000)
+        vars.lo('R', vars.get('R') / 1000)
+        vars.up('R', vars.get('R') * 1000)
+        vars.lo('CMI', vars.get('CMI') / 1000)
+        vars.up('CMI', vars.get('CMI') * 1000)
+        vars.lo('CMO', vars.get('CMO') / 1000)
+        vars.up('CMO', vars.get('CMO') * 1000)
+        vars.lo('N', 0)
+        vars.lo('CN', 0)
+
+        def set_variable(filname):
+
+            # clear the file before write the variables
+            with open(filename, 'w') as f:
+                f.write("")
+
+            vars.write(filename)
+
+        # -------------------------------------------------------------------------------------------------------------
+        # DEFINE EQUATIONS AND SET CONDITIONS
+        # -------------------------------------------------------------------------------------------------------------
+
+        def set_equation(filename):
+
+            count = [0]
+
+            #   CPIEQ(H)..
+            #   CPI(H)=E=
+            #   SUM(I, P(I) * ( 1 + SUM(GS, TAUC(GS,I) ) ) * CH(I,H) )
+            #       / SUM(I, P0(I) * ( 1 + SUM(GS, TAUQ(GS,I) ) ) * CH(I,H) );
+            print('CPIEQ(H)')
+            line1 = (P.loc(I) * ExprM(vars, m=1 + TAUC.loc[GS, I].sum(0)) * CH.loc(I, H)).sum(I)
+            line2 = (ExprM(vars, m=P0.loc[I] * (1 + TAUQ.loc[GS, I].sum(0))) * CH.loc(I, H)).sum(I)
+
+            CPIEQ = (line1 / line2 - ~CPI.loc(H))
+            print(CPIEQ.test(vars.initialVals))
+            CPIEQ.write(count, filename)
+
+            # YEQ(H).. Y(H) =E= SUM(L,  A(H,L) * HW(H) / SUM(H1, A(H1,L) * HW(H1) ) * (Y(L)-(CMIWAGE(L)*CMI(L))) * ( 1 - SUM(G, TAUFL(G,L))))
+            #                + SUM(CM, A(H,CM)*(CMOWAGE(CM)*CMO(CM)))
+            #                + SUM(K,  A(H,K) * HW(H) / SUM(H1, A(H1,K) * HW(H1)) * (Y(K) + KPFOR(K)) * ( 1 - SUM(G, TAUFK(G,K) ) ) ) ;
+            print('YEQ(H)')
+            line1 = (ExprM(vars, m=A.loc[H, L]) * HW.loc(H) / (ExprM(vars, m=A.loc[H1, L]) * HW.loc(H1)).sum(H1) * (
+                        Y.loc(L) - ExprM(vars, m=CMIWAGE.loc[L]) * CMI.loc(L)) * ExprM(vars, m=1 - TAUFL.loc[G, L].sum(
+                0))).sum(L)
+            line2 = (ExprM(vars, m=A.loc[H, CM] * CMOWAGE.loc[CM]) * CMO.loc(CM)).sum(CM)
+            # line3 = (ExprM(vars, m= A.loc[H,LA]) * HW.loc(H) / (ExprM(vars, m= A.loc[H1,LA]) * HW.loc(H1)).sum(H1) * ExprM(vars, m= 1 - TAUFLA.loc[G,LA].sum(0))).sum(LA)
+            line4 = (ExprM(vars, m=A.loc[H, K]) * HW.loc(H) / (ExprM(vars, m=A.loc[H1, K]) * HW.loc(H1)).sum(H1) * (
+                        Y.loc(K) + KPFOR.loc(K)) * ExprM(vars, m=1 - TAUFK.loc[G, K].sum(0))).sum(K)
+
+            YEQ = ((line1 + line2 + line4) - Y.loc(H))
+            print(YEQ.test(vars.initialVals))
+            YEQ.write(count, filename)
+            # print(YEQ)
+
+            #  YDEQ(H).. YD(H)          =E=   Y(H) + (PRIVRET(H) * HH(H))
+            #                                         + SUM(G, TP(H,G) * HH(H))
+            #                                         - SUM(GI, PIT0(GI,H)  * Y(H))
+            #                                         - SUM(G, TAUH(G,H)  * HH(H));
+            print('YDEQ(H)')
+            line1 = Y.loc(H) + ExprM(vars, m=PRIVRET.loc[H]) * HH.loc(H)
+            line2 = (ExprM(vars, m=TP.loc[H, G]) * HH.loc(H)).sum(G)
+            line3 = ~(ExprM(vars, m=PIT0.loc[GI, H]) * Y.loc(H)).sum(GI)
+            line4 = ~(ExprM(vars, m=TAUH.loc[G, H]) * HH.loc(H)).sum(G)
+
+            YDEQ = ((line1 + line2 - line3 - line4) - YD.loc(H))
+            YDEQ.write(count, filename)
+            print(YDEQ.test(vars.initialVals))
+            # print(YDEQ)
+
+            #  CHEQ(I,H).. CH(I,H)      =E= CH0(I,H)* ((YD(H) / YD0(H)) / ( CPI(H) / CPI0(H)))**(BETA(I,H))
+            #                                 * PROD(J, ((P(J)*( 1 + SUM(GS, TAUC(GS,J))))/ (P0(J)*(1 + SUM(GS, TAUQ(GS,J)))))** (LAMBDA(J,I)));
+            print('CHEQ(I,H)')
+            line1 = ExprM(vars, m=CH0.loc[I, H]) * (
+                        (YD.loc(H) / ExprM(vars, m=YD0.loc[H])) / (CPI.loc(H) / ExprM(vars, m=CPI0.loc[H]))) ** ExprM(
+                vars, m=BETA.loc[I, H])
+            line2 = (((P.loc(J) * ExprM(vars, m=1 + TAUC.loc[GS, J].sum(0))) / (
+                        ExprM(vars, m=P0.loc[J]) * ExprM(vars, m=1 + TAUQ.loc[GS, J].sum(0)))) ** ExprM(vars,
+                                                                                                        m=LAMBDA.loc[
+                                                                                                            J, I])).prod(
+                0)
+
+            CHEQ = ((line1 * line2) - CH.loc(I, H))
+            CHEQ.write(count, filename)
+            print(CHEQ.test(vars.initialVals))
+            # print(CHEQ)
+
+            #  SHEQ(H).. S(H)           =E= YD(H) - SUM(I, P(I) * CH(I,H) * ( 1 + SUM(GS, TAUC(GS,I))));
+            print('SHEQ(H)')
+            line = YD.loc(H) - ~((P.loc(I) * CH.loc(I, H) * ExprM(vars, m=1 + TAUC.loc[GS, I].sum(0))).sum(I))
+
+            SHEQ = (line - S.loc(H))
+            SHEQ.write(count, filename)
+            print(SHEQ.test(vars.initialVals))
+            # print(SHEQ)
+
+            #  PVAEQ(I).. PVA(I)        =E= PD(I) - SUM(J, AD(J,I) * P(J) * (1 + SUM(GS, TAUQ(GS, J))));
+            print('PVAEQ(I)')
+            line = PD.loc(I) - ~(
+                (ExprM(vars, m=AD.loc[J, I]) * P.loc(J) * ExprM(vars, m=1 + TAUQ.loc[GS, J].sum(0))).sum(0))
+
+            PVAEQ = (line - PVA.loc(I))
+            PVAEQ.write(count, filename)
+            print(PVAEQ.test(vars.initialVals))
+            # print(PVAEQ)
+
+            # Cobb-Douglas production function
+
+            #  PFEQ(I)..DS(I)           =E= DELTA(I)*PROD(F, (FD(F,I))**ALPHA(F,I));
+            #  print('PFEQ(I)')
+            #  line =  ExprM(vars, m =DELTA.loc[I]) *((FD.loc(F, I)) ** ExprM(vars, m = ALPHA.loc[F, I])).prod(F)
+
+            #  PFEQ = (line - DS.loc(I))
+            #  PFEQ.write(count,filename)
+            # print(PFEQ)
+
+            # CES production function
+            # PFEQ(I)..DS(I) =E= GAMMA(I)*(SUM(F, ALPHA(F,I)*(FD(F,I)**(RHO(I)))))**(1/RHO(I));
+            print('PFEQ(I)')
+            line = ExprM(vars, m=GAMMA.loc[I]) * (
+                (ExprM(vars, m=ALPHA.loc[F, I]) * (FD.loc(F, I) ** ExprM(vars, m=RHO.loc[I]))).sum(F)) ** ExprM(vars,
+                                                                                                                m=1 /
+                                                                                                                  RHO.loc[
+                                                                                                                      I])
+
+            PFEQ = (line - DS.loc(I))
+            PFEQ.write(count, filename)
+            print(PFEQ.test(vars.initialVals))
+            # print(PFEQ)
+
+            # Cobb-Douglas factor demand function
+            # FDEQ(F,I).. R(F,I) * RA(F) * (1 + SUM(GF,TAUFX(GF,F,I) ) )* (FD(F,I))
+            #                         =E= PVA(I) * DS(I) * ALPHA(F,I);
+            #  print('FDEQ(F,I)')
+            #  left = R.loc(F,I) * RA.loc(F) * ExprM(vars, m = 1 + TAUFX_SUM.loc[F,I])* (FD.loc(F,I))
+            #  right = ~(PVA.loc(I) * DS.loc(I) * ExprM(vars, m = ALPHA.loc[F, I]))
+
+            #  FDEQ = (right - left)
+
+            # FDEQ.test(vars.initialVals)
+            #  FDEQ.write(count,filename)
+            # print(FDEQ)
+
+            # CES factor demand function
+            # FDEQ(F,I).. (R(F,I) * RA(F)*(1 + SUM(GF,TAUFX(GF,F,I))))* (FD(F,I)**(1-RHO(I)))
+            #               =E= PVA(I)* ALPHA(F,I)*(GAMMA(I)**RHO(I))*(DS(I)**(1-RHO(I)));
+            print('FDEQ(F,I)')
+            left = (R.loc(F, I) * RA.loc(F) * ExprM(vars, m=1 + TAUFX_SUM.loc[F, I])) * (
+                        FD.loc(F, I) ** ExprM(vars, m=1 - RHO.loc[I]))
+            right = ~(PVA.loc(I) * (ExprM(vars, m=ALPHA.loc[F, I] * (GAMMA.loc[I] ** RHO.loc[I])) * (
+                        DS.loc(I) ** ExprM(vars, m=1 - RHO.loc[I]))))
+            FDEQ = (right - left)
+            FDEQ.setCondition(FD0.loc[F, I], 'INEQ', 0)
+
+            # FDEQ.test(vars.initialVals)
+            FDEQ.write(count, filename)
+            print(FDEQ.test(vars.initialVals))
+            # print(FDEQ)
+
+            #   VEQ(I).. V(I) =E= SUM(J, AD(I,J) * DS(J) );
+            print('VEQ(I)')
+            line = (ExprM(vars, m=AD.loc[I, J]) * ~DS.loc(J)).sum(1)
+
+            VEQ = (line - V.loc(I))
+            VEQ.write(count, filename)
+            print(VEQ.test(vars.initialVals))
+            # print(VEQ)
+
+            #   YFEQL(F).. Y(F) =E= SUM(IG, R(F,IG)* RA(F)*FD(F,IG));
+            print('YFEQL(F)')
+            line = (R.loc(F, IG) * RA.loc(F) * FD.loc(F, IG)).sum(IG)
+
+            YFEQL = (line - Y.loc(F))
+            YFEQL.write(count, filename)
+            print(YFEQL.test(vars.initialVals))
+            # print(YFEQL)
+
+            # YFEQK(K).. Y('KAP') =E= SUM(IG, R('KAP',IG) * RA('KAP') * FD('KAP',IG));
+            print('YFEQK(K)')
+            line = (R.loc(['KAP'], IG) * RA.loc(['KAP']) * FD.loc(['KAP'], IG)).sum(IG)
+
+            YFEQK = (line - Y.loc(['KAP']))
+            YFEQK.write(count, filename)
+            print(YFEQK.test(vars.initialVals))
+            # print(YFEQK)
+
+            #  YFEQLA(LA).. Y('LAND')   =E= SUM(IG, R('LAND',IG) * RA('LAND') * FD('LAND',IG));
+            #  print('YFEQLA(LA)')
+            #  line = (R.loc(['LAND'],IG) * RA.loc(['LAND']) * FD.loc(['LAND'],IG)).sum(IG)
+
+            #  YFEQLA = (line - Y.loc(['LAND']))
+            #  YFEQLA.write(count,filename)
+            #  print(YFEQLA.test(vars.initialVals))
+            # print(YFEQLA)
+
+            #  LANFOR(LA).. LNFOR(LA)   =E= LFOR(LA) * Y(LA);
+            #  print('LANFOR(LA)')
+            #  line = ExprM(vars, m = LFOR.loc[LA])*Y.loc(LA)
+
+            #  LANFOR = (line - LNFOR.loc(LA))
+            #  LANFOR.write(count,filename)
+            #  print(LANFOR.test(vars.initialVals))
+            # print(LANFOR)
+
+            #  KAPFOR(K).. KPFOR(K)     =E= KFOR(K) * Y(K);
+            print('KAPFOR(K)')
+            line = ExprM(vars, m=KFOR.loc[K]) * Y.loc(K)
+
+            KAPFOR = (line - KPFOR.loc(K))
+            KAPFOR.write(count, filename)
+            print(KAPFOR.test(vars.initialVals))
+            # print(KAPFOR)
+
+            # XEQ(I).. CX(I) =E= CX0(I)*( (PD(I)*(1+SUM(GS,TAUX(GS,I))))
+            #                  /(PW0(I)*(1+SUM(GS,TAUQ(GS,I))))) **(ETAE(I));
+            print('XEQ(I)')
+            line = ExprM(vars, m=CX0.loc[I]) * ((PD.loc(I) * ExprM(vars, m=1 + TAUX.loc[GS, I].sum(0))) / ExprM(vars, m=
+            PW0.loc[I] * (1 + TAUQ.loc[GS, I].sum(0)))) ** ExprM(vars, m=ETAE.loc[I])
+
+            XEQ = (line - CX.loc(I))
+            XEQ.write(count, filename)
+            print(XEQ.test(vars.initialVals))
+            # print(XEQ)
+
+            #  DEQ(I)$PWM0(I).. D(I)    =E= D0(I) *(PD(I)/PWM0(I))**(ETAD(I)*0.1);
+            print('DEQ(I)$PWM0(I)')
+            line = ExprM(vars, m=D0.loc[I]) * (PD.loc(I) / ExprM(vars, m=PWM0.loc[I])) ** ExprM(vars,
+                                                                                                m=ETAD.loc[I] * 0.1)
+
+            DEQ = (line - D.loc(I))
+            #  DEQ.setCondition(PWM0.loc[I])
+            DEQ.write(count, filename)
+            print(DEQ.test(vars.initialVals))
+            # print(DEQ)
+
+            #  PEQ(I)..  P(I)           =E= D(I) * PD(I) + ( 1 - D(I) ) * PWM0(I);
+            print('PEQ(I)')
+            line = (D.loc(I) * PD.loc(I) + (1 - D.loc(I)) * ExprM(vars, m=PWM0.loc[I]))
+
+            PEQ = (line - P.loc(I))
+            PEQ.write(count, filename)
+            print(PEQ.test(vars.initialVals))
+            # print(PEQ)
+
+            #  MEQ(I).. M(I)            =E= ( 1 - D(I) ) * DD(I);
+            print('MEQ(I)')
+            line = (1 - D.loc(I)) * DD.loc(I)
+
+            MEQ = (line - M.loc(I))
+            MEQ.write(count, filename)
+            print(MEQ.test(vars.initialVals))
+            # print(MEQ)
+
+            #  NKIEQ.. NKI              =E= SUM(I, M(I) * PWM0(I) )
+            #                                 - SUM(I, CX(I) * PD(I) )
+            #                                 - SUM(H, PRIVRET(H)*HH(H))
+            #                                 - SUM(K, KPFOR(K))
+            #                                 - SUM(G, GVFOR(G))
+            #                                 - SUM(CM,CMOWAGE(CM)*CMO(CM))
+            #                                 + SUM(L,CMIWAGE(L)*CMI(L));
+
+            print('NKIEQ')
+            line1 = (M.loc(I) * ExprM(vars, m=PWM0.loc[I])).sum(I)
+            line2 = (CX.loc(I) * PD.loc(I)).sum(I)
+            line3 = (ExprM(vars, m=PRIVRET.loc[H]) * HH.loc(H)).sum(H)
+            line5 = KPFOR.loc(K).sum(K)
+            line6 = GVFOR.loc(G).sum(G)
+            line7 = (ExprM(vars, m=CMOWAGE.loc[CM]) * CMO.loc(CM)).sum(CM)
+            line8 = (ExprM(vars, m=CMIWAGE.loc[L]) * CMI.loc(L)).sum(L)
+
+            NKIEQ = ((line1 - line2 - line3 - line5 - line6 - line7 + line8) - NKI)
+            NKIEQ.write(count, filename)
+            print(NKIEQ.test(vars.initialVals))
+            # print(NKIEQ)
+
+            #  NEQ(K,I).. N(K,I)        =E= N0(K,I)*(R(K,I)/R0(K,I))**(ETAIX(K,I)*0.1);
+            print('NEQ(K,I)')
+            line = ExprM(vars, m=N0.loc[K, I]) * (R.loc(K, I) / ExprM(vars, m=R0.loc[K, I])) ** ExprM(vars, m=ETAIX.loc[
+                                                                                                                  K, I] * 0.1)
+
+            NEQ = (line - N.loc(K, I))
+            NEQ.write(count, filename)
+            print(NEQ.test(vars.initialVals))
+            # print(NEQ)
+
+            #  CNEQ(I).. P(I)*(1 + SUM(GS, TAUN(GS,I)))*CN(I)
+            #                         =E= SUM(IG, B(I,IG)*(SUM(K, N(K,IG))));
+            print('CNEQ(I)')
+            left = P.loc(I) * ExprM(vars, m=1 + TAUN.loc[GS, I].sum(0)) * CN.loc(I)
+            right = (ExprM(vars, m=B.loc[I, IG]) * N.loc(K, IG).sum(K)).sum(IG)
+
+            CNEQ = (right - left)
+            CNEQ.write(count, filename)
+            print(CNEQ.test(vars.initialVals))
+            # print(CNEQ)
+
+            #  KSEQ(K,IG).. KS(K,IG)    =E= KS0(K,IG) * ( 1 - DEPR) + N(K,IG) ;
+            print('KSEQ(K,IG)')
+            line = ExprM(vars, m=KS0.loc[K, IG] * (1 - DEPR)) + N.loc(K, IG)
+
+            KSEQ = (line - KS.loc(K, IG))
+            KSEQ.write(count, filename)
+            print(KSEQ.test(vars.initialVals))
+            # print(KSEQ)
+
+            # LSEQ1(H).. HW(H)/HH(H)   =E= (HW0(H)/HH0(H))
+            #                              *((SUM(L, RA(L) / RA0(L))/24)/ (CPI(H) / CPI0(H))*(SUM((Z,L), FD(L,Z))/(SUM(H1, HW(H1)* SUM(L, JOBCOR(H1,L)))+ SUM(CM, CMO(CM)) + SUM(L,CMI(L))))+ SUM((CM,L), EXWGEO(CM)/RA(L))/96 *(SUM(CM, CMO(CM))/(SUM(H1, HW(H1)* SUM(L,JOBCOR(H1,L)))+ SUM(CM, CMO(CM)) +SUM(L,CMI(L)))))** (ETARA(H))
+            #                              * ( SUM(G, TP(H,G) / CPI(H) )/ SUM(G, TP(H,G) / CPI0(H) )) ** ETAPT(H)
+            #                              *  ((SUM(GI, PIT0(GI,H)* HH0(H))+ SUM(G, TAUH(G,H)*HH0(H)))/(SUM(GI, PIT(GI,H)* HH(H))+ SUM(G, TAUH(G,H)*HH(H))))**(ETAPIT(H));
+            print('LSEQ1(H)')
+            line1 = ExprM(vars, m=HW0.loc[H] / HH0.loc[H])
+            LSEQ1line2pre = FD.loc(L, Z).sum(1)
+            line2 = (((RA.loc(L) / ExprM(vars, m=RA0.loc[L])).sum(L) / 8) / (CPI.loc(H) / ExprM(vars, m=CPI0[H])) \
+                     * (LSEQ1line2pre.sum(0) / (
+                                (HW.loc(H1) * ExprM(vars, m=JOBCOR.loc[H1, L].sum(1))).sum(H1) + CMO.loc(CM).sum(
+                            CM) + CMI.loc(L).sum(L))) \
+                     + ((ExprM(vars, m=EXWGEO.loc[CM].sum(0)) / RA.loc(L)).sum(0) / 32) \
+                     * (CMO.loc(CM).sum(CM) / (
+                                (HW.loc(H1) * ExprM(vars, m=JOBCOR.loc[H1, L].sum(1))).sum(H1) + CMO.loc(CM).sum(
+                            CM) + CMI.loc(L).sum(L)))) \
+                    ** (ExprM(vars, m=ETARA.loc[H]) * 0.1)
+            line3 = ((ExprM(vars, m=TP.loc[H, G]) / CPI.loc(H)).sum(G) / (
+                        ExprM(vars, m=TP.loc[H, G]) / ExprM(vars, m=CPI0.loc[H])).sum(G)) ** ExprM(vars, m=ETAPT.loc[H])
+            line4 = (((ExprM(vars, m=PIT0.loc[GI, H]) * ExprM(vars, m=HH0.loc[H])).sum(GI) + (
+                        ExprM(vars, m=TAUH0.loc[G, H]) * ExprM(vars, m=HH0.loc[H])).sum(G)) \
+                     / ((ExprM(vars, m=PIT.loc[GI, H]) * HH.loc(H)).sum(GI) + (
+                                ExprM(vars, m=TAUH.loc[G, H]) * HH.loc(H)).sum(G))) ** ExprM(vars, m=ETAPIT.loc[H])
+
+            LSEQ1 = ((line1 * line2 * line3 * line4) - HW.loc(H) / HH.loc(H))
+            LSEQ1.write(count, filename)
+            print(LSEQ1.test(vars.initialVals))
+            # print(LSEQ1)
+
+            # LSEQ2A('OUT1').. CMO('OUT1')=E= CMO0('OUT1')* (((EXWGEO('OUT1') /RA('L1WA') ))** ECOMO('OUT1'));
+            print('LSEQ2A')
+            line = ExprM(vars, m=CMO0.loc[CM1]) * (
+                        (ExprM(vars, m=EXWGEO.loc[CM1].sum(0)) / RA.loc(FW1)) ** (ExprM(vars, m=ECOMO.loc[CM1])))
+
+            LSEQ2A = (line - CMO.loc(CM1))
+            LSEQ2A.write(count, filename)
+            print(LSEQ2A.test(vars.initialVals))
+            # print(LSEQ2A)
+
+            # LSEQ2B('OUT2').. CMO('OUT2')=E= CMO0('OUT2')* (((EXWGEO('OUT2') /RA('L2WA') ))** ECOMO('OUT2'));
+            print('LSEQ2B')
+            line = ExprM(vars, m=CMO0.loc[CM2]) * (
+                        (ExprM(vars, m=EXWGEO.loc[CM2].sum(0)) / RA.loc(FW2)) ** (ExprM(vars, m=ECOMO.loc[CM2])))
+
+            LSEQ2B = (line - CMO.loc(CM2))
+            LSEQ2B.write(count, filename)
+            print(LSEQ2B.test(vars.initialVals))
+            # print(LSEQ2B)
+
+            # LSEQ2C('OUT3').. CMO('OUT3')=E= CMO0('OUT3')* (((EXWGEO('OUT3') /RA('L3WA') ))** ECOMO('OUT3'));
+            print('LSEQ2C')
+            line = ExprM(vars, m=CMO0.loc[CM3]) * (
+                        (ExprM(vars, m=EXWGEO.loc[CM3].sum(0)) / RA.loc(FW3)) ** (ExprM(vars, m=ECOMO.loc[CM3])))
+
+            LSEQ2C = (line - CMO.loc(CM3))
+            LSEQ2C.write(count, filename)
+            print(LSEQ2C.test(vars.initialVals))
+            # print(LSEQ2C)
+
+            # LSEQ2D('OUT4').. CMO('OUT4')=E= CMO0('OUT4')* (((EXWGEO('OUT4') /RA('L4WA') ))** ECOMO('OUT4'));
+            print('LSEQ2D')
+            line = ExprM(vars, m=CMO0.loc[CM4]) * (
+                        (ExprM(vars, m=EXWGEO.loc[CM4].sum(0)) / RA.loc(FW4)) ** (ExprM(vars, m=ECOMO.loc[CM4])))
+
+            LSEQ2D = (line - CMO.loc(CM4))
+            LSEQ2D.write(count, filename)
+            print(LSEQ2D.test(vars.initialVals))
+            # print(LSEQ2D)
+
+            # LSEQ3(L).. CMI(L)  =E= CMI0(L)* ((( RA(L)/(SUM( H, CPI(H))/30))** ECOMI(L)));
+            print('LSEQ3')
+            line = ExprM(vars, m=CMI0.loc[L]) * ((RA.loc(L) / (CPI.loc(H).sum(H) / 8)) ** (ExprM(vars, m=ECOMI.loc[L])))
+
+            LSEQ3 = (line - CMI.loc(L))
+            LSEQ3.write(count, filename)
+            print(LSEQ3.test(vars.initialVals))
+            # print(LSEQ3)
+
+            #  LASEQ1(LA,I).. LAS(LA,I) =E= LAS0(LA,I)*(R(LA, I)/R0(LA, I))**(ETAL(LA,I));
+            #  print('LASEQ1(LA,I)')
+            #  line = ExprM(vars, m = LAS0.loc[LA, I]) * (R.loc(LA,I) / ExprM(vars, m = R0.loc[LA, I]))**ExprM(vars, m = ETAL.loc[LA,I])
+
+            #  LASEQ1 = (line - LAS.loc(LA, I))
+            #  LASEQ1.write(count,filename)
+            #  print(LASEQ1.test(vars.initialVals))
+            # print(LASEQ1)
+
+            #  POPEQ(H).. HH(H)         =E= HH0(H) * NRPG(H)
+            #                                + MI0(H) * ((YD(H)/HH(H))/(YD0(H)/HH0(H))/(CPI(H)/CPI0(H))) ** (ETAYD(H))
+            #                                         *((HN(H)/HH(H))/(HN0(H)/HH0(H))) ** (ETAU(H))
+            #                                - MO0(H) *((YD0(H)/HH0(H))/(YD(H)/HH(H))/(CPI0(H)/CPI(H))) ** (ETAYD(H))
+            #                                         *((HN0(H)/HH0(H))/(HN(H)/HH(H))) ** (ETAU(H))  ;
+
+            print('POPEQ(H)')
+            line1 = ExprM(vars, m=HH0.loc[H] * NRPG.loc[H])
+            line2 = ExprM(vars, m=MI0.loc[H]) * ((YD.loc(H) / HH.loc(H)) / ExprM(vars, m=YD0.loc[H] / HH0.loc[H]) / (
+                        CPI.loc(H) / ExprM(vars, m=CPI0.loc[H]))) ** ExprM(vars, m=ETAYD.loc[H])
+            line3 = ((HN.loc(H) / HH.loc(H)) / ExprM(vars, m=HN0.loc[H] / HH0.loc[H])) ** ExprM(vars, m=ETAU.loc[H])
+            line4 = ExprM(vars, m=MO0.loc[H]) * (ExprM(vars, m=YD0.loc[H] / HH0.loc[H]) / (YD.loc(H) / HH.loc(H)) / (
+                        ExprM(vars, m=CPI0.loc[H]) / CPI.loc(H))) ** ExprM(vars, m=ETAYD.loc[H])
+            line5 = (ExprM(vars, m=HN0.loc[H] / HH0.loc[H]) / (HN.loc(H) / HH.loc(H))) ** ExprM(vars, m=ETAU.loc[H])
+
+            POPEQ = (line1 + line2 * line3 - line4 * line5 - HH.loc(H))
+            POPEQ.write(count, filename)
+            print(POPEQ.test(vars.initialVals))
+            # print(POPEQ)
+
+            #  ANEQ(H).. HN(H)          =E= HH(H) - HW(H);
+            print('ANEQ(H)')
+            line = HH.loc(H) - HW.loc(H)
+
+            ANEQ = (line - HN.loc(H))
+            ANEQ.write(count, filename)
+            print(ANEQ.test(vars.initialVals))
+            # print(ANEQ)
+
+            #  YGEQ(GX).. Y(GX)         =E=   SUM(I, TAUV(GX,I) * V(I) * P(I) )
+            #                                 + SUM(I, TAUX(GX,I)* CX(I) *PD(I))
+            #                                 + SUM((H,I), TAUC(GX,I) * CH(I,H) * P(I) )
+            #                                 + SUM(I, TAUN(GX,I) * CN(I) * P(I) )
+            #                                 + SUM((GN,I), TAUG(GX,I) * CG(I,GN) * P(I) )
+            #                                 + SUM((F,I), TAUFX(GX,F,I) * RA(F) * R(F,I) * FD(F,I) )
+            #                                 + SUM((F,GN), TAUFX(GX,F,GN) * RA(F) * R(F,GN) * FD(F,GN) )
+            #                                 + SUM(F, TAUFH(GX,F) * (Y(F)))
+            #                                 + SUM(H, PIT0(GX,H) * Y(H) )
+            #                                 + SUM(H, TAUH(GX,H) * HH(H) )
+            #                                 + SUM(GX1, IGT(GX,GX1));
+            print('YGEQ')
+            line1 = (ExprM(vars, m=TAUV.loc[GX, I]) * V.loc(I) * P.loc(I)).sum(I)
+            line2 = (ExprM(vars, m=TAUX.loc[GX, I]) * CX.loc(I) * PD.loc(I)).sum(I)
+
+            YGEQline3pre = ExprM(vars,
+                                 m=pd.DataFrame(index=GX, columns=H).fillna(0.0))  # first just set it to correct size
+            for label in GX:
+                for hlabel in H:
+                    YGEQline3pre.m[GX.index(label)][H.index(hlabel)] = Expr(
+                        (ExprM(vars, m=TAUC.loc[[label], I]) * CH.loc(I, [hlabel]) * P.loc(I)).sum().m[0][0])
+
+            line3 = YGEQline3pre.sum(1)
+            line4 = (ExprM(vars, m=TAUN.loc[GX, I]) * CN.loc(I) * P.loc(I)).sum(I)
+
+            YGEQline5pre = ExprM(vars, m=TAUG.loc[GX, I]).sum(I)  # first just set it to correct size
+            for label in GX:
+                YGEQline5pre.m[GX.index(label)][0] = Expr(
+                    (ExprM(vars, m=TAUG.loc[[label], I]) * CG.loc(I, GN) * P.loc(I)).sum().m[0][0])
+            line5 = YGEQline5pre
+
+            line6 = ExprM(vars, m=TAUG.loc[GX, I]).sum(I)
+            for label in GX:
+                line6.m[GX.index(label)][0] = Expr(
+                    (ExprM(vars, m=TAUFXgx[label].loc[F, I]) * RA.loc(F) * R.loc(F, I) * FD.loc(F, I)).sum().m[0][0])
+
+            line7 = ExprM(vars, m=TAUG.loc[GX, I]).sum(I)
+            for label in GX:
+                line7.m[GX.index(label)][0] = Expr(
+                    (ExprM(vars, m=TAUFXgx[label].loc[F, GN]) * RA.loc(F) * R.loc(F, GN) * FD.loc(F, GN)).sum().m[0][0])
+
+            line8 = (ExprM(vars, m=TAUFH.loc[GX, F]) * Y.loc(F)).sum(F)
+            line9 = (ExprM(vars, m=PIT0.loc[GX, H]) * Y.loc(H)).sum(H)
+            line10 = (ExprM(vars, m=TAUH.loc[GX, H]) * HH.loc(H)).sum(H)
+            line11 = IGT.loc(GX, GX1).sum(1)
+
+            YGEQ = ((line1 + line2 + line3 + line4 + line5 + line6 + line7 + line8 + line9 + line10 + line11) - Y.loc(
+                GX))
+            YGEQ.write(count, filename)
+            print(YGEQ.test(vars.initialVals))
+            # print(YGEQ)
+
+            #    YGEQ2(GT).. Y(GT)        =E= SUM(GX, IGT(GT,GX));
+            print('YGEQ2')
+            line = IGT.loc(GT, GX).sum(GX)
+            YGEQ2 = (line - Y.loc(GT))
+            YGEQ2.write(count, filename)
+            print(YGEQ2.test(vars.initialVals))
+            # print(YGEQ2)
+
+            #  YGEQL(GNL).. Y(GNL)      =E= TAXS1(GNL)*Y('CYGF');
+            print('YGEQL(GNL)')
+            line = ExprM(vars, m=TAXS1.loc[GNL]) * Y.loc(['CYGF'])
+            YGEQ1 = (line - Y.loc(GNL))
+            YGEQ1.write(count, filename)
+            print(YGEQ1.test(vars.initialVals))
+            # print(YGEQ1)
+
+            #  GOVFOR(G).. GVFOR(G)     =E= GFOR(G)*Y(G);
+            print('GOVFOR(G)')
+            line = ExprM(vars, m=GFOR.loc[G]) * Y.loc(G)
+
+            GOVFOR = (line - GVFOR.loc(G))
+            GOVFOR.write(count, filename)
+            print(GOVFOR.test(vars.initialVals))
+            # print(GOVFOR)
+
+            #  CGEQ(I,GN).. P(I)*(1 + SUM(GS, TAUG(GS,I))) * CG(I,GN)
+            #                         =E= AG(I,GN) * (Y(GN)+ GFOR(GN)*Y(GN));
+            print('CGEQ(I,GN)')
+            left = P.loc(I) * ExprM(vars, m=1 + TAUG.loc[GS, I].sum(0)) * CG.loc(I, GN)
+            right = ExprM(vars, m=AG.loc[I, GN]) * (Y.loc(GN) + ExprM(vars, m=GFOR.loc[GN]) * Y.loc(GN))
+
+            CGEQ = (right - left)
+            CGEQ.write(count, filename)
+            print(CGEQ.test(vars.initialVals))
+            # print(CGEQ)
+
+            #  GFEQ(F,GN)..  FD(F,GN) * R(F,GN) * RA(F)*( 1 + SUM(GF, TAUFX(GF,F,GN)))
+            #                         =E= AG(F,GN) * (Y(GN)+ GFOR(GN)*Y(GN));
+            print('GFEQ(F,GN)')
+            left = FD.loc(F, GN) * R.loc(F, GN) * RA.loc(F) * (1 + ExprM(vars, m=TAUFX_SUM.loc[F, GN]))
+            right = ExprM(vars, m=AG.loc[F, GN]) * (Y.loc(GN) + ExprM(vars, m=GFOR.loc[GN]) * Y.loc(GN))
+
+            GFEQ = left - right
+            GFEQ.write(count, filename)
+            print(GFEQ.test(vars.initialVals))
+            # print(GFEQ)
+
+            #  GSEQL(GN).. S(GN)        =E= (Y(GN)+ GVFOR(GN))
+            #                                 - SUM(I, CG(I,GN)*P(I)*(1 + SUM(GS, TAUG(GS,I))))
+            #                                 - SUM(F, FD(F,GN)*R(F,GN)*RA(F)*(1 + SUM(GF, TAUFX(GF,F,GN))));
+            print('GSEQL(GN)')
+            line1 = Y.loc(GN) + GVFOR.loc(GN)
+            line2 = (CG.loc(I, GN) * P.loc(I) * (1 + ExprM(vars, m=TAUG.loc[GS, I]).sum(GS))).sum(I)
+            line3 = (FD.loc(F, GN) * R.loc(F, GN) * RA.loc(F) * (1 + ExprM(vars, m=TAUFX_SUM.loc[F, GN]))).sum(F)
+
+            GSEQL = ((line1 - ~line2 - ~line3) - S.loc(GN))
+            GSEQL.write(count, filename)
+            print(GSEQL.test(vars.initialVals))
+            # print(GSEQL)
+
+            #  GSEQ(GX).. S(GX)         =E= (Y(GX) + GFOR(GX)*Y(GX)) - SUM(H, (TP(H,GX)*HH(H))) - SUM(G,IGT(G,GX));
+            print('GSEQ(GX)')
+            line1 = (Y.loc(GX) + ExprM(vars, m=GFOR.loc[GX]) * Y.loc(GX))
+            line2 = (ExprM(vars, m=TP.loc[H, GX]) * HH.loc(H)).sum(H)
+            line3 = IGT.loc(G, GX).sum(G)
+
+            GSEQ = ((line1 - ~line2 - ~line3) - S.loc(GX))
+            GSEQ.write(count, filename)
+            print(GSEQ.test(vars.initialVals))
+            # print(GSEQ)
+
+            #  TDEQ(G,GX)$(IGTD(G,GX) EQ 1).. IGT(G,GX)
+            #                         =E= TAXS(G,GX)*(Y(GX) + GVFOR(GX)- SUM(H, (TP(H,GX)*HH(H))));
+            print('TDEQ(G,GX)$(IGTD(G,GX) EQ 1)')
+            line = ExprM(vars, m=TAXS.loc[G, GX]) * (
+                        Y.loc(GX) + GVFOR.loc(GX) - ~(ExprM(vars, m=TP.loc[H, GX]) * HH.loc(H)).sum(H))
+
+            TDEQ = line - IGT.loc(G, GX)
+            TDEQ.setCondition(IGTD.loc[G, GX], 'EQ', 1)
+            TDEQ.write(count, filename)
+            print(TDEQ.test(vars.initialVals))
+            # print(TDEQ)
+
+            #  SPIEQ.. SPI =E= SUM(H, Y(H)) + SUM((H,G), TP(H,G)*HH(H)) + SUM(H, PRIVRET(H)*HH(H));
+            print('SPIEQ')
+            line = Y.loc(H).sum(H) + (ExprM(vars, m=TP.loc[H, G]) * HH.loc(H)).sum() + (
+                        ExprM(vars, m=PRIVRET.loc[H]) * HH.loc(H)).sum(H)
+
+            SPIEQ = (line - SPI)
+            SPIEQ.write(count, filename)
+            print(SPIEQ.test(vars.initialVals))
+            # print(SPIEQ)
+
+            #  LMEQ1(L).. SUM(H, HW(H)* JOBCOR(H,L)) + CMI(L) =E= SUM(Z, FD(L,Z)));
+            print('LMEQ1')
+            left = (ExprM(vars, m=JOBCOR.loc[H, L]) * HW.loc(H)).sum(H) + CMI.loc(L)
+            right = FD.loc(L, Z).sum(Z)
+            # right = FD.loc(['L1'], Z).sum(Z) + CMO.loc(CM1).sum(CM1)
+
+            LMEQ1 = (right - left)
+            LMEQ1.write(count, filename)
+            print(LMEQ1.test(vars.initialVals))
+            # print(LMEQ1)
+
+            #  KMEQ(K,IG).. KS(K,IG)    =E= FD(K,IG);
+            print('KMEQ(K,IG)')
+            KMEQ = (FD.loc(K, IG) - KS.loc(K, IG))
+            KMEQ.write(count, filename)
+            print(KMEQ.test(vars.initialVals))
+            # print(KMEQ)
+
+            #  LAMEQ(LA,IG).. LAS(LA,IG)=E= FD(LA,IG);
+            #  print('LAMEQ(LA,IG)')
+            #  LAMEQ = (FD.loc(LA, IG) - LAS.loc(LA, IG))
+            #  LAMEQ.write(count,filename)
+            #  print(LAMEQ.test(vars.initialVals))
+            # print(LAMEQ)
+
+            #  GMEQ(I).. DS(I)          =E= DD(I) + CX(I) - M(I);
+            print('GMEQ(I)')
+            GMEQ = (DD.loc(I) + CX.loc(I) - M.loc(I) - DS.loc(I))
+            GMEQ.write(count, filename)
+            print(GMEQ.test(vars.initialVals))
+            # print(GMEQ)
+
+            #  DDEQ(I).. DD(I)          =E= V(I) + SUM(H, CH(I,H) ) + SUM(G, CG(I,G) ) + CN(I);
+            print('DDEQ(I)')
+            DDEQ = (V.loc(I) + CH.loc(I, H).sum(H) + CG.loc(I, G).sum(G) + CN.loc(I) - DD.loc(I))
+            DDEQ.write(count, filename)
+            print(DDEQ.test(vars.initialVals))
+            # print(DDEQ)
+
+            # IGT.FX(G,GX)$(NOT IGT0(G,GX))=0;
+            print('IGT.FX(G,GX)$(NOT IGT0(G,GX))=0')
+            FX1 = IGT.loc(G, GX)
+            FX1.setCondition(IGT0.loc[G, GX], 'EQ', 0)
+            FX1.write(count, filename)
+            # print(FX1)
+
+            # IGT.FX(G,GX)$(IGTD(G,GX) EQ 2)=IGT0(G,GX);
+            print('IGT.FX(G,GX)$(IGTD(G,GX) EQ 2)=IGT0(G,GX)')
+            FX2 = IGT.loc(G, GX) - ExprM(vars, m=IGT0.loc[G, GX])
+            FX2.setCondition(IGTD.loc[G, GX], 'EQ', 2)
+            FX2.write(count, filename)
+            # print(FX2)
+
+            # R.FX(L,Z)=R0(L,Z);
+            print('R.FX(L,Z)=R0(L,Z)')
+            FX3 = R.loc(L, Z) - ExprM(vars, m=R0.loc[L, Z])
+            FX3.write(count, filename)
+            # print(FX3)
+
+            '''
+            #RA.FX(LA)=RA0(LA);
+            print('RA.FX(LA)=RA0(LA)')
+            FX4 = RA.loc(LA) - ExprM(vars, m = RA0.loc[LA])
+            FX4.write(count,filename)
+            # print(FX4)
+            '''
+
+            # RA.FX(K)=RA0(K);
+            print('RA.FX(K)=RA0(K)')
+            FX5 = RA.loc(K) - ExprM(vars, m=RA0.loc[K])
+            FX5.write(count, filename)
+            # print(FX5)
+
+            print("Objective")
+            obj = vars.getIndex('SPI')
+
+            with open(filename, 'a') as f:
+                f.write('model.obj = Objective(expr=-1*model.x' + str(obj) + ')')
+
+        def run_solver(cons_filename, temp_file_name="tmp.py"):
+            solver = 'ipopt'
+            solver_io = 'nl'
+            stream_solver = True  # True prints solver output to screen
+            keepfiles = False  # True prints intermediate file names (.nl,.sol,...)
+
+            executable_path = self.get_parameter("solver_path") \
+                if self.get_parameter("solver_path") is not None else pyglobals.IPOPT_PATH
+            if not os.path.exists(executable_path):
+                print("Invalid executable path, please make sure you have Pyomo installed.")
+
+            opt = SolverFactory(solver, solver_io=solver_io, executable=executable_path)
+
+            if opt is None:
+                print("")
+                print("ERROR: Unable to create solver plugin for %s "
+                      "using the %s interface" % (solver, solver_io))
+                print("")
+                exit(1)
+
+            ### Create the model
+            model = ConcreteModel()
+            set_variable(cons_filename)
+            set_equation(cons_filename)
+            ###
+
+            ### read the model
+            exec(open(cons_filename).read())
+            ###
+
+            ### Declare all suffixes
+
+            # Ipopt bound multipliers (obtained from solution)
+            model.ipopt_zL_out = Suffix(direction=Suffix.IMPORT)
+            model.ipopt_zU_out = Suffix(direction=Suffix.IMPORT)
+
+            # Ipopt bound multipliers (sent to solver)
+            model.ipopt_zL_in = Suffix(direction=Suffix.EXPORT)
+            model.ipopt_zU_in = Suffix(direction=Suffix.EXPORT)
+
+            # Obtain dual solutions from first solve and send to warm start
+            model.dual = Suffix(direction=Suffix.IMPORT_EXPORT)
+
+            # opt.options['halt_on_ampl_error'] = 'yes'
+            # opt.options['acceptable_tol'] = '1e-3'
+            # Send the model to ipopt and collect the solution
+            # results = opt.solve(model,keepfiles=keepfiles,tee=stream_solver)
+
+            # IMPORT / IMPORT_EXPORT Suffix components)
+            # model.solutions.load_from(results)
+
+            # Set Ipopt options for warm-start
+            # The current values on the ipopt_zU_out and
+            # ipopt_zL_out suffixes will be used as initial
+            # conditions for the bound multipliers to solve
+            # the new problem
+            model.ipopt_zL_in.update(model.ipopt_zL_out)
+            model.ipopt_zU_in.update(model.ipopt_zU_out)
+            opt.options['warm_start_init_point'] = 'yes'
+            opt.options['warm_start_bound_push'] = 1e-6
+            opt.options['warm_start_mult_bound_push'] = 1e-6
+            opt.options['mu_init'] = 1e-6
+            # opt.options['acceptable_tol'] = 10
+            # opt.options['max_iter'] = 15
+            # opt.options['compl_inf_tol'] = 1e-3
+            # opt.options['bound_relax_factor'] = 0
+            # opt.options['start_with_resto'] = 'yes'
+            # opt.options['acceptable_iter'] = 15
+            # opt.options['halt_on_ampl_error'] = 'yes'
+            # opt.options['fixed_variable_treatment'] = 'relax_bounds'
+            # opt.options['print_options_documentation'] = 'yes'
+
+            ###
+
+            # Send the model and suffix information to ipopt and collect the solution
+            # The solver plugin will scan the model for all active suffixes
+            # valid for importing, which it will store into the results object
+
+            results = opt.solve(model, keepfiles=keepfiles, tee=stream_solver)
+
+            if results.solver.status == SolverStatus.ok and results.solver.termination_condition == TerminationCondition.optimal:
+
+                x = [None for i in range(vars.nvars)]
+
+                with open(temp_file_name, 'w') as f:
+                    for i in range(vars.nvars):
+                        f.write('x[' + str(i) + ']=value(model.x' + str(i) + ')\n')
+
+                exec(open(temp_file_name).read())
+
+                soln.append(x[:])
+
+                return True
+            elif results.solver.termination_condition == TerminationCondition.infeasible:
+                print("Solution is not feasible.")
+                return False
+            else:
+                print("Solver Status: ", results.solver.status)
+                return False
+
+        '''
+        Calibrate the model 
+        '''
+
+        soln = []
+        # TODO: I am not sure this is needed here. We might want to leave the python version
+        # of the models out for now
+        filename = "ipopt_cons.py"
+        tmp = "tmp.py"
+        print("Calibration: ")
+        run_solver(filename, tmp)
+
+        # sys.exit()
+
+        '''
+        Simulation code below:
+        In each simulation:
+        
+        1. Apply simulation code (for instance PI(I) = 1.02).
+        2. Rewrite all equations
+        3. Solve the new model with the result from last run as initial guess.
+        
+        '''
+
+        '''
+        ######## The following is for individual shocks. ######## 
+        
+        iNum = 1 # dynamic model itterations
+        
+        for ittr in range(iNum):
+          print("Simulation: ", ittr+1)
+          if ittr == 0: # if it is the first simulation, apply the shock
+        
+              #DELTA.loc[I] = 1.02 * DELTA.loc[I]
+        
+              KS0.loc[K, I] = KS0.loc[K, I]*0.7
+        
+              #KS0.loc[K, ['HS1']] = KS0.loc[K, ['HS1']] * 0.675
+              #KS0.loc[K, ['HS2']] = KS0.loc[K, ['HS2']] * 0.739
+              #KS0.loc[K, ['HS3']] = KS0.loc[K, ['HS3']] * 0.958
+              #KS0.loc[K, ['GOODS']] = KS0.loc[K, ['GOODS']] * 0.658
+              #KS0.loc[K, ['TRADE']] = KS0.loc[K, ['TRADE']] * 0.961
+              #KS0.loc[K, ['OTHER']] = KS0.loc[K, ['OTHER']] * 0.673
+        
+        
+          else: # other simulations
+        
+              #KS0 = KSNEW*(1-DEPR)+vars.get('N', x=soln[-1])
+              KS0 = vars.get('KS', x=soln[-1])
+        
+          run_solver(filename, tmp)
+        
+        '''
+
+        # The following is for random shocks. #
+        sims = sector_shocks
+        iNum = 1  # dynamic model itterations
+        # iNum = len(sims.columns)
+        KS00 = KS0.copy()
+
+        solver_status = False
+        for num in range(iNum):
+            KS0.loc[K, I] = KS00.loc[K, I].mul(sims.iloc[:, num])
+            KS0 = KS0.fillna(0.0)
+            solver_status = run_solver(filename, tmp)
+
+        if solver_status:
+            domestic_supply, gross_income, household_count, pre_disaster_demand, post_disaster_demand = \
+                gams_to_dataframes(iNum, vars, H, L, soln)
+
+            self.set_result_csv_data("domestic-supply", domestic_supply, name="domestic-supply", source="dataframe",
+                                     index=True)
+            self.set_result_csv_data("pre-disaster-factor-demand", pre_disaster_demand,
+                                     name="pre-disaster-factor-demand", source="dataframe", index=True)
+            self.set_result_csv_data("post-disaster-factor-demand", post_disaster_demand,
+                                     name="post-disaster-factor-demand", source="dataframe", index=True)
+            self.set_result_csv_data("gross-income", gross_income, name="gross-income", source="dataframe", index=True)
+            self.set_result_csv_data("household-count", household_count, name="household-count", source="dataframe",
+                                     index=True)
+        else:
+            raise ValueError("Solution infeasible - no output to save")
+
+    def get_spec(self):
+        return {
+            'name': 'Salt-Lake-cge',
+            'description': 'CGE model for Salt Lake City.',
+            'input_parameters': [
+                {
+                    'id': 'model_iterations',
+                    'required': True,
+                    'description': 'Number of dynamic model iterations.',
+                    'type': int
+                },
+                {
+                    'id': 'solver_path',
+                    'required': False,
+                    'description': 'Path to ipopt package. If none is provided, it will default to your environment\'ts'
+                                   'path to the package.',
+                    'type': str
+                }
+
+            ],
+            'input_datasets': [
+                {
+                    'id': 'SAM',
+                    'required': True,
+                    'description': 'Social accounting matrix (SAM) contains data for firms, '
+                                   'households and government which are organized in a way to '
+                                   'represent the interactions of all three entities in a typical economy.',
+                    'type': ['incore:CGEsam']
+                },
+                {
+                    'id': 'BB',
+                    'required': True,
+                    'description': 'BB is a matrix which describes how investment in physical infrastructure is'
+                                   ' transformed into functioning capital such as commercial and residential buildings.'
+                                   ' These data are collected from the Bureau of Economic Analysis (BEA).',
+                    'type': ['incore:CGEbb']
+                },
+                {
+                    'id': 'IOUT',
+                    'required': False,
+                    'description': 'IOUT is a matrix that describes the transfer of tax revenue collected by the local'
+                                   ' government to help finance local government expenditures.',
+                    'type': ['incore:CGEiout']
+                },
+                {
+                    'id': 'MISC',
+                    'required': False,
+                    'description': 'MISC is the name of a file that contains data for commercial sector employment'
+                                   ' and physical capital. It also contains data for the number of households and'
+                                   ' working households in the economy.',
+                    'type': ['incore:CGEmisc']
+                },
+                {
+                    'id': 'MISCH',
+                    'required': True,
+                    'description': 'MISCH is a file that contains elasticities for the supply of labor with'
+                                   ' respect to paying income taxes.',
+                    'type': ['incore:CGEmisch']
+                },
+                {
+                    'id': 'EMPLOY',
+                    'required': True,
+                    'description': 'EMPLOY is a table name containing data for commercial sector employment.',
+                    'type': ['incore:CGEemploy']
+                },
+                {
+                    'id': 'JOBCR',
+                    'required': True,
+                    'description': 'JOBCR is a matrix describing the supply of workers'
+                                   ' coming from each household group in the economy.',
+                    'type': ['incore:CGEjobcr']
+                },
+                {
+                    'id': 'OUTCR',
+                    'required': True,
+                    'description': 'OUTCR is a matrix describing the number of workers who'
+                                   ' live in Salt Lake City but commute outside of town to work.',
+                    'type': ['incore:CGEoutcr']
+                },
+                {
+                    'id': 'sector_shocks',
+                    'required': True,
+                    'description': 'Aggregation of building functionality states to capital shocks per sector',
+                    'type': ['incore:capitalShocks']
+                }
+            ],
+            'output_datasets': [
+                {
+                    'id': 'domestic-supply',
+                    'parent_type': '',
+                    'description': 'CSV file of resulting domestic supply',
+                    'type': 'incore:Employment'
+                },
+                {
+                    'id': 'gross-income',
+                    'parent_type': '',
+                    'description': 'CSV file of resulting gross income',
+                    'type': 'incore:Employment'
+                },
+                {
+                    'id': 'pre-disaster-factor-demand',
+                    'parent_type': '',
+                    'description': 'CSV file of factor demand before disaster',
+                    'type': 'incore:FactorDemand'
+                },
+                {
+                    'id': 'post-disaster-factor-demand',
+                    'parent_type': '',
+                    'description': 'CSV file of resulting factor-demand',
+                    'type': 'incore:FactorDemand'
+                },
+                {
+                    'id': 'household-count',
+                    'parent_type': '',
+                    'description': 'CSV file of household count',
+                    'type': 'incore:HouseholdCount'
+                }
+            ]
+        }
diff --git a/pyincore/analyses/saltlakecge/saltlakeoutput.py b/pyincore/analyses/saltlakecge/saltlakeoutput.py
new file mode 100644
index 000000000..b01c1203c
--- /dev/null
+++ b/pyincore/analyses/saltlakecge/saltlakeoutput.py
@@ -0,0 +1,220 @@
+# Copyright (c) 2021 University of Illinois and others. All rights reserved.
+#
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+import pandas as pd
+
+
+def gams_to_dataframes(iNum, vars, H, L, soln):
+    # CG0
+    CG0 = vars.get('CG', x=soln[0])
+
+    # CH0
+    CH0 = vars.get('CH', x=soln[0])
+
+    # CMI0
+    CMI0 = vars.get('CMI', x=soln[0])
+
+    # CMO0
+    CMO0 = vars.get('CMO', x=soln[0])
+
+    # CN0
+    CN0 = vars.get('CN', x=soln[0])
+
+    # CPI0
+    CPI0 = vars.get('CPI', x=soln[0])
+
+    # CX0
+    CX0 = vars.get('CX', x=soln[0])
+
+    # D0
+    D0 = vars.get('D', x=soln[0])
+
+    # DD0
+    DD0 = vars.get('DD', x=soln[0])
+
+    # DS0
+    DS0 = vars.get('DS', x=soln[0])
+
+    # FD
+    FD0 = vars.get('FD', x=soln[0])
+
+    # IGT
+    IGT0 = vars.get('IGT', x=soln[0])
+
+    # KS
+    KS0 = vars.get('KS', x=soln[0])
+
+    # LAS
+    # LAS0 = vars.get('LAS', x=soln[0])
+
+    # HH
+    HH0 = vars.get('HH', x=soln[0])
+
+    # HN
+    HN0 = vars.get('HN', x=soln[0])
+
+    # HW
+    HW0 = vars.get('HW', x=soln[0])
+
+    # M
+    M0 = vars.get('M', x=soln[0])
+
+    # N
+    N0 = vars.get('N', x=soln[0])
+
+    # NKI
+    NKI0 = vars.get('NKI', x=soln[0])
+
+    # LNFOR
+    # LNFOR0 = vars.get('LNFOR', x=soln[0])
+
+    # KPFOR
+    KPFOR0 = vars.get('KPFOR', x=soln[0])
+
+    # GVFOR
+    GVFOR0 = vars.get('GVFOR', x=soln[0])
+
+    # P
+    P0 = vars.get('P', x=soln[0])
+
+    # PD
+    PD0 = vars.get('PD', x=soln[0])
+
+    # PVA
+    PVA0 = vars.get('PVA', x=soln[0])
+
+    # RA
+    RA0 = vars.get('RA', x=soln[0])
+
+    # R
+    R0 = vars.get('R', x=soln[0])
+
+    # S
+    S0 = vars.get('S', x=soln[0])
+
+    # SPI
+    SPI0 = vars.get('SPI', x=soln[0])
+
+    # V
+    V0 = vars.get('V', x=soln[0])
+
+    # Y
+    Y0 = vars.get('Y', x=soln[0])
+
+    # Yd
+    YD0 = vars.get('YD', x=soln[0])
+
+    for i in range(iNum):
+        CGL = vars.get('CG', x=soln[i + 1])
+        CHL = vars.get('CH', x=soln[i + 1])
+        CMIL = vars.get('CMI', x=soln[i + 1])
+        CMOL = vars.get('CMO', x=soln[i + 1])
+        CNL = vars.get('CN', x=soln[i + 1])
+        CPIL = vars.get('CPI', x=soln[i + 1])
+        CXL = vars.get('CX', x=soln[i + 1])
+        DL = vars.get('D', x=soln[i + 1])
+        DDL = vars.get('DD', x=soln[i + 1])
+        DSL = vars.get('DS', x=soln[i + 1])
+        FDL = vars.get('FD', x=soln[i + 1])
+        IGTL = vars.get('IGT', x=soln[i + 1])
+        KSL = vars.get('KS', x=soln[i + 1])
+        # LASL = vars.get('LAS', x=soln[i+1])
+        HHL = vars.get('HH', x=soln[i + 1])
+        HNL = vars.get('HN', x=soln[i + 1])
+        HWL = vars.get('HW', x=soln[i + 1])
+        ML = vars.get('M', x=soln[i + 1])
+        NL = vars.get('N', x=soln[i + 1])
+        NKIL = vars.get('NKI', x=soln[i + 1])
+        # LNFORL = vars.get('LNFOR', x=soln[i+1])
+        KPFORL = vars.get('KPFOR', x=soln[i + 1])
+        GVFORL = vars.get('GVFOR', x=soln[i + 1])
+        PL = vars.get('P', x=soln[i + 1])
+        PDL = vars.get('PD', x=soln[i + 1])
+        PVAL = vars.get('PVA', x=soln[i + 1])
+        RAL = vars.get('RA', x=soln[i + 1])
+        RL = vars.get('R', x=soln[i + 1])
+        SL = vars.get('S', x=soln[i + 1])
+        SPIL = vars.get('SPI', x=soln[i + 1])
+        VL = vars.get('V', x=soln[i + 1])
+        YL = vars.get('Y', x=soln[i + 1])
+        YDL = vars.get('YD', x=soln[i + 1])
+
+        # DFCG.L(I,G)      = CG.L(I,G)-CG0(I,G);
+        DFCG = CGL - CG0
+
+        # DFFD.L(F,Z)      = FD.L(F,Z)-FD0(F,Z);
+        DFFD = FDL - FD0
+
+        # DK.L(K,Z)        = FD.L(K,Z)-FD0(K,Z);
+        DK = KSL - KS0
+
+        # DY.L(Z)          = Y.L(Z)-Y0(Z);
+        # DY = YL - Y0
+        DY = (YL / CPIL) - Y0
+
+        # DDS.L(I)         = DS.L(I)-DS0(I);
+        DDS = DSL - DS0
+
+        # DDD.L(I)         = DD.L(I) - DD0(I);
+        DDD = DDL - DD0
+
+        # DCX.L(I)         = CX.L(I) -CX0(I);
+        DCX = CXL - CX0
+
+        # DCH.L(I,H)       = CH.L(I,H) - CH0(I,H);
+        DCH = CHL - CH0
+
+        # DRR.L(F,Z)       = R.L(F,Z)-R0(F,Z);
+        DR = RL - R0
+
+        # DCMI.L(L)        = CMI.L(L) - CMI0(L);
+        DCMI = CMIL - CMI0
+
+        # DCMO.L(CM)       = CMO.L(CM) - CMO0(CM);
+        DCMO = CMOL - CMO0
+
+        DM = ML - M0
+
+        DV = VL - V0
+
+        DN = NL - N0
+
+        s_name = 'Simulation ' + str(i + 1)
+
+    '''
+    export domestic supply, household income (gross income), number of household, and factor demand
+    '''
+    # domestic supply
+    cols = {'DS0': DS0,
+            'DSL': DSL}
+    domestic_supply_df = pd.DataFrame.from_dict(cols)
+    domestic_supply_df.index.name = 'Sectors'
+    # df.to_csv('domestic-supply.csv')
+
+    # gross income
+    cols = {'Y0': Y0.loc[H],
+            'YL': YL.loc[H]}
+    gross_income_df = pd.DataFrame.from_dict(cols)
+    gross_income_df.index.name = 'Household Group'
+    # df.to_csv('gross-income.csv')
+
+    # household count
+    cols = {'HH0': HH0,
+            'HHL': HHL}
+    household_count_df = pd.DataFrame.from_dict(cols)
+    household_count_df.index.name = 'Household Group'
+    # df.to_csv('household-count.csv')
+
+    # pre-disaster-factor-demand
+    pre_disaster_demand_df = FD0.loc[L]
+    pre_disaster_demand_df.index.name = 'Labor Group'
+    # df.to_csv('pre-disaster-factor-demand.csv')
+
+    # post-disaster-factor-demand
+    post_disaster_demand_df = FDL.loc[L]
+    post_disaster_demand_df.index.name = 'Labor Group'
+    # df.to_csv('post-disaster-factor-demand.csv')
+
+    return domestic_supply_df, gross_income_df, household_count_df, pre_disaster_demand_df, post_disaster_demand_df
diff --git a/pyincore/baseanalysis.py b/pyincore/baseanalysis.py
index 34784c32c..b63d1adec 100644
--- a/pyincore/baseanalysis.py
+++ b/pyincore/baseanalysis.py
@@ -241,7 +241,7 @@ def validate_output_dataset(dataset_spec, dataset):
 
     """ convenience function(s) for setting result data as a csv """
 
-    def set_result_csv_data(self, result_id, result_data, name, source='file'):
+    def set_result_csv_data(self, result_id, result_data, name, source='file', index=False):
         if name is None:
             name = self.spec["name"] + "-result"
 
@@ -254,7 +254,7 @@ def set_result_csv_data(self, result_id, result_data, name, source='file'):
         if source == 'file':
             dataset = Dataset.from_csv_data(result_data, name, dataset_type)
         elif source == 'dataframe':
-            dataset = Dataset.from_dataframe(result_data, name, dataset_type)
+            dataset = Dataset.from_dataframe(result_data, name, dataset_type, index)
 
         self.set_output_dataset(result_id, dataset)
 
diff --git a/pyincore/dataset.py b/pyincore/dataset.py
index da51932ab..ac587800e 100644
--- a/pyincore/dataset.py
+++ b/pyincore/dataset.py
@@ -108,19 +108,20 @@ def from_file(cls, file_path, data_type):
         return instance
 
     @classmethod
-    def from_dataframe(cls, dataframe, name, data_type):
+    def from_dataframe(cls, dataframe, name, data_type, index=False):
         """Get Dataset from Panda's DataFrame.
 
         Args:
             dataframe (obj): Panda's DataFrame.
             name (str): filename.
             data_type (str): Incore data type, e.g. incore:xxxx or ergo:xxxx
+            index (bool): Store the index column
 
         Returns:
             obj: Dataset from file.
 
         """
-        dataframe.to_csv(name, index=False)
+        dataframe.to_csv(name, index=index)
         return Dataset.from_file(name, data_type)
 
     @classmethod
diff --git a/pyincore/globals.py b/pyincore/globals.py
index b98aa8e56..f17efacbc 100644
--- a/pyincore/globals.py
+++ b/pyincore/globals.py
@@ -10,7 +10,7 @@
 import os
 import shutil
 
-PACKAGE_VERSION = "1.8.0"
+PACKAGE_VERSION = "1.9.0"
 
 INCORE_API_PROD_URL = "https://incore.ncsa.illinois.edu"
 INCORE_API_DEV_URL = "https://incore-dev.ncsa.illinois.edu"
diff --git a/pyincore/spaceservice.py b/pyincore/spaceservice.py
index 2b59785e5..0199ff770 100644
--- a/pyincore/spaceservice.py
+++ b/pyincore/spaceservice.py
@@ -73,6 +73,24 @@ def get_space_by_id(self, space_id: str):
 
         return response
 
+    def get_space_by_name(self, space_name: str):
+        """Get space information by querying the name of space.
+
+        Args:
+            space_name (str): A space representation. Name of the Space.
+
+        Returns:
+            obj: HTTP response with the returned space information.
+
+        """
+        r = self.client.get(self.base_space_url, params={"name": space_name})
+        if r.status_code == 200:
+            return r.json()
+        else:
+            # throw an error instead of returning the empty result
+            # return []
+            raise r.exceptions.HTTPError("There is no matching name or you don't have a privilege to view the space.")
+
     def update_space(self, space_id: str, space_json):
         """Updates a Space.
 
@@ -91,6 +109,57 @@ def update_space(self, space_id: str, space_json):
         response = r.json()
         return response
 
+    def add_to_space_by_name(self, space_name: str, dataset_id: str):
+        """Add dataset to a space by using space name and dataset id.
+
+        Args:
+            space_name (str): Name of the space that dataset id will go to.
+            dataset_id (json): ID of the dataset that will be copied in to the space.
+
+        Returns:
+            obj: HTTP response with the updated space.
+
+        """
+        space_id = self.get_space_by_name(space_name)[0]
+
+        response = self.add_dataset_to_space(space_id["id"], dataset_id)
+
+        return response
+
+    def remove_from_space_by_name(self, space_name: str, dataset_id: str):
+        """Remove dataset from a space by using space name and dataset id.
+
+        Args:
+            space_name (str): Name of the space that dataset id will go to.
+            dataset_id (json): ID of the dataset that will be copied in to the space.
+
+        Returns:
+            obj: HTTP response with the updated space.
+
+        """
+        space_id = self.get_space_by_name(space_name)[0]
+
+        response = self.remove_dataset_from_space(space_id["id"], dataset_id)
+
+        return response
+
+    def remove_dataset_from_space(self, space_id: str, dataset_id: str):
+        """Remove dataset from the space using dataset id and space id
+
+        Args:
+            space_id (str): ID of the space to update.
+            dataset_id (str): ID of the dataset to be removed from the space.
+
+        Returns:
+            obj: HTTP response with the updated space.
+
+        """
+        url = urljoin(self.base_space_url, space_id + "/members/" + dataset_id)
+
+        r = self.client.delete(url)
+        response = r.json()
+        return response
+
     def add_dataset_to_space(self, space_id: str, dataset_id: str):
         """Add member to a Space.
 
diff --git a/pyincore/utils/cgeoutputprocess.py b/pyincore/utils/cgeoutputprocess.py
index 35dafe89c..e34018215 100644
--- a/pyincore/utils/cgeoutputprocess.py
+++ b/pyincore/utils/cgeoutputprocess.py
@@ -12,7 +12,8 @@ class CGEOutputProcess:
     """This class converts csv results outputs of Joplin CGE analysis to json format."""
 
     @staticmethod
-    def get_cge_household_count(household_count, household_count_path=None, filename_json=None):
+    def get_cge_household_count(household_count, household_count_path=None, filename_json=None,
+                                income_categories=("HH1", "HH2", "HH3", "HH4", "HH5")):
         """Calculate income results from the output files of the Joplin CGE analysis and convert the results
         to json format.
         {
@@ -26,12 +27,12 @@ def get_cge_household_count(household_count, household_count_path=None, filename
             household_count_path (obj): A fallback for the case that household count object of CGE is not provided.
                  For example a user wants to directly pass in csv files, a path to CGE household count result.
             filename_json (str): Path and name to save json output file in. E.g "cge_total_household_count.json"
+            income_categories (list): A list of income categories to partition the data
 
         Returns:
             obj: CGE total household count. A JSON of the total household count results ordered by category.
 
         """
-        income_categories = ["HH1", "HH2", "HH3", "HH4", "HH5"]
 
         if household_count_path:
             household_group_count = pd.read_csv(household_count_path)
@@ -56,12 +57,13 @@ def get_cge_household_count(household_count, household_count_path=None, filename
 
         if filename_json:
             with open(filename_json, "w") as outfile:
-                json.dump(cge_total_household_count, outfile)
+                json.dump(cge_total_household_count, outfile, indent=2)
         # Serializing json
         return json.dumps(cge_total_household_count)
 
     @staticmethod
-    def get_cge_gross_income(gross_income, gross_income_path=None, filename_json=None):
+    def get_cge_gross_income(gross_income, gross_income_path=None, filename_json=None,
+                             income_categories=("HH1", "HH2", "HH3", "HH4", "HH5")):
         """Calculate household gross income results from the output files of the Joplin CGE analysis
         and convert the results to json format.
         {
@@ -75,12 +77,12 @@ def get_cge_gross_income(gross_income, gross_income_path=None, filename_json=Non
             gross_income_path (obj): A fallback for the case that gross_income object of CGE is not provided.
                  For example a user wants to directly pass in csv files, a path to CGE gross income result.
             filename_json (str): Path and name to save json output file in. E.g "cge_total_house_income.json"
+            income_categories (list): A list of income categories to partition the data
 
         Returns:
             obj: CGE total house income. A JSON of the total household income results ordered by category.
 
         """
-        income_categories = ["HH1", "HH2", "HH3", "HH4", "HH5"]
 
         if gross_income_path:
             household_income = pd.read_csv(gross_income_path)
@@ -105,14 +107,13 @@ def get_cge_gross_income(gross_income, gross_income_path=None, filename_json=Non
 
         if filename_json:
             with open(filename_json, "w") as outfile:
-                json.dump(cge_total_household_income, outfile)
+                json.dump(cge_total_household_income, outfile, indent=2)
         # Serializing json
         return json.dumps(cge_total_household_income)
 
     @staticmethod
-    def get_cge_employment(pre_demand, post_demand,
-                           pre_demand_path=None, post_demand_path=None,
-                           filename_json=None):
+    def get_cge_employment(pre_demand, post_demand, pre_demand_path=None, post_demand_path=None,
+                           filename_json=None, demand_categories=("GOODS", "TRADE", "OTHER")):
         """Calculate employment results from the output files of the Joplin CGE analysis and convert the results
         to json format. The value is a sum of L1, L2 and L3 Labor groups numbers.
         {
@@ -135,12 +136,12 @@ def get_cge_employment(pre_demand, post_demand,
                 of CGE is not provided. For example a user wants to directly pass in csv files, a path to CGE
                 household count result.
             filename_json (str): Path and name to save json output file in. E.g "cge_employment.json"
+            demand_categories (list): demand categories to partition data with.
 
         Returns:
             obj: CGE total employment. A JSON of the employment results ordered by category.
 
         """
-        demand_categories = ["Goods", "Trade", "Other"]
 
         if pre_demand_path and post_demand_path:
             pre_disaster_demand = pd.read_csv(pre_demand_path)
@@ -149,12 +150,8 @@ def get_cge_employment(pre_demand, post_demand,
             pre_disaster_demand = pre_demand.get_dataframe_from_csv()
             post_disaster_demand = post_demand.get_dataframe_from_csv()
 
-        before_values = [pre_disaster_demand["GOODS"].sum(),
-                         pre_disaster_demand["TRADE"].sum(),
-                         pre_disaster_demand["OTHER"].sum()]
-        after_values = [post_disaster_demand["GOODS"].sum(),
-                        post_disaster_demand["TRADE"].sum(),
-                        post_disaster_demand["OTHER"].sum()]
+        before_values = [pre_disaster_demand[demand_category].sum() for demand_category in demand_categories]
+        after_values = [post_disaster_demand[demand_category].sum() for demand_category in demand_categories]
 
         before_event = {}
         after_event = {}
@@ -171,12 +168,13 @@ def get_cge_employment(pre_demand, post_demand,
 
         if filename_json:
             with open(filename_json, "w") as outfile:
-                json.dump(cge_employment, outfile)
+                json.dump(cge_employment, outfile, indent=2)
         # Serializing json
         return json.dumps(cge_employment)
 
     @staticmethod
-    def get_cge_domestic_supply(domestic_supply, domestic_supply_path=None, filename_json=None):
+    def get_cge_domestic_supply(domestic_supply, domestic_supply_path=None, filename_json=None,
+                                supply_categories=("Goods", "Trade", "Other", "HS1", "HS2", "HS3")):
         """Calculate domestic supply results from the output files of the Joplin CGE analysis and convert the results
         to json format.
         {
@@ -190,6 +188,7 @@ def get_cge_domestic_supply(domestic_supply, domestic_supply_path=None, filename
 
         Args:
             domestic_supply (obj): IN-CORE dataset for CGE domestic supply result.
+            supply_categories (list): demand categories to partition data with.
             domestic_supply_path (obj): A fallback for the case that domestic supply object of CGE is not provided.
                  For example a user wants to directly pass in csv files, a path to CGE household count result.
             filename_json (str): Path and name to save json output file in. E.g "cge_domestic_supply"
@@ -198,7 +197,7 @@ def get_cge_domestic_supply(domestic_supply, domestic_supply_path=None, filename
             obj: CGE total domestic supply. A JSON of the total domestic supply results ordered by category.
 
         """
-        supply_categories = ["Goods", "Trade", "Other", "HS1", "HS2", "HS3"]
+
         if domestic_supply_path:
             sector_supply = pd.read_csv(domestic_supply_path)
         else:
@@ -222,6 +221,6 @@ def get_cge_domestic_supply(domestic_supply, domestic_supply_path=None, filename
 
         if filename_json:
             with open(filename_json, "w") as outfile:
-                json.dump(cge_domestic_supply, outfile)
+                json.dump(cge_domestic_supply, outfile, indent=2)
         # Serializing json
         return json.dumps(cge_domestic_supply)
diff --git a/pyincore/utils/dataprocessutil.py b/pyincore/utils/dataprocessutil.py
index f6c679a3f..6724650f8 100644
--- a/pyincore/utils/dataprocessutil.py
+++ b/pyincore/utils/dataprocessutil.py
@@ -16,7 +16,8 @@ class DataProcessUtil:
     @staticmethod
     def get_mapped_result_from_analysis(client, inventory_id: str, dmg_result_dataset,
                                         bldg_func_dataset, archetype_mapping_id: str,
-                                        groupby_col_name: str = "max_state"):
+                                        groupby_col_name: str = "max_state", arch_col='archetype'
+                                        ):
         """Use this if you want to load results directly from the output files of the analysis, than storing the results
         to data service and loading from there using ids.
         It takes the static inputs: inventory & archetypes as dataset ids. The result inputs are taken as
@@ -33,6 +34,7 @@ def get_mapped_result_from_analysis(client, inventory_id: str, dmg_result_datase
             dmg_ret_json: JSON of the damage state results ordered by cluster and category.
             func_ret_json: JSON of the building functionality results ordered by cluster and category.
             mapped_df: Dataframe of max damage state
+            arch_col: column name for the archetype to perform the merge
 
         """
         bldg_inv = Dataset.from_data_service(inventory_id, DataService(client))
@@ -47,15 +49,15 @@ def get_mapped_result_from_analysis(client, inventory_id: str, dmg_result_datase
 
         max_state_df = DataProcessUtil.get_max_damage_state(dmg_result)
         dmg_ret_json = DataProcessUtil.create_mapped_dmg_result(inventory, max_state_df, arch_mapping,
-                                                                groupby_col_name)
-        func_ret_json = DataProcessUtil.create_mapped_func_result(inventory, bldg_func_df, arch_mapping)
+                                                                groupby_col_name, arch_col)
+        func_ret_json = DataProcessUtil.create_mapped_func_result(inventory, bldg_func_df, arch_mapping, arch_col)
 
         return dmg_ret_json, func_ret_json, max_state_df
 
     @staticmethod
     def get_mapped_result_from_dataset_id(client, inventory_id: str, dmg_result_id: str, bldg_func_id,
                                           archetype_mapping_id: str,
-                                          groupby_col_name: str = "max_state"):
+                                          groupby_col_name: str = "max_state", arch_col='archetype'):
         """Use this if your damage results are already stored in the data service and you have their dataset ids.
         All the inputs (except groupby_col_name) are dataset ids.
 
@@ -66,6 +68,7 @@ def get_mapped_result_from_dataset_id(client, inventory_id: str, dmg_result_id:
             bldg_func_id: Incore dataset for building func id
             archetype_mapping_id: Mapping id dataset for archetype
             groupby_col_name: column name to group by, default to max_state
+            arch_col: column name for the archetype to perform the merge
 
         Returns:
             dmg_ret_json: JSON of the damage state results ordered by cluster and category.
@@ -88,8 +91,8 @@ def get_mapped_result_from_dataset_id(client, inventory_id: str, dmg_result_id:
 
         max_state_df = DataProcessUtil.get_max_damage_state(dmg_result)
         dmg_ret_json = DataProcessUtil.create_mapped_dmg_result(inventory, max_state_df, arch_mapping,
-                                                                groupby_col_name)
-        func_ret_json = DataProcessUtil.create_mapped_func_result(inventory, bldg_func_df, arch_mapping)
+                                                                groupby_col_name, arch_col)
+        func_ret_json = DataProcessUtil.create_mapped_func_result(inventory, bldg_func_df, arch_mapping, arch_col)
 
         return dmg_ret_json, func_ret_json, max_state_df
 
@@ -97,7 +100,7 @@ def get_mapped_result_from_dataset_id(client, inventory_id: str, dmg_result_id:
     def get_mapped_result_from_path(inventory_path: str, dmg_result_path: str,
                                     func_result_path: str,
                                     archetype_mapping_path: str,
-                                    groupby_col_name: str):
+                                    groupby_col_name: str, arch_col='archetype'):
         """
 
         Args:
@@ -107,6 +110,7 @@ def get_mapped_result_from_path(inventory_path: str, dmg_result_path: str,
             func_result_path: Path to the bldg functionality result output file
             archetype_mapping_path: Path to the arechetype mappings
             groupby_col_name: column name to group by, default to max_state
+            arch_col: column name for the archetype to perform the merge
 
         Returns:
             dmg_ret_json: JSON of the damage state results ordered by cluster and category.
@@ -122,20 +126,21 @@ def get_mapped_result_from_path(inventory_path: str, dmg_result_path: str,
 
         max_state_df = DataProcessUtil.get_max_damage_state(dmg_result)
         dmg_ret_json = DataProcessUtil.create_mapped_dmg_result(inventory, max_state_df, arch_mapping,
-                                                                groupby_col_name)
+                                                                groupby_col_name, arch_col)
 
-        func_ret_json = DataProcessUtil.create_mapped_func_result(inventory, bldg_func_df, arch_mapping)
+        func_ret_json = DataProcessUtil.create_mapped_func_result(inventory, bldg_func_df, arch_mapping, arch_col)
 
         return dmg_ret_json, func_ret_json, max_state_df
 
     @staticmethod
-    def create_mapped_dmg_result(inventory, dmg_result, arch_mapping, groupby_col_name="max_state"):
+    def create_mapped_dmg_result(inventory, dmg_result, arch_mapping, groupby_col_name="max_state",
+                                 arch_col='archetype'):
         """
 
         Args:
             inventory: dataframe represent inventory
             dmg_result: dmg_result that need to merge with inventory and be grouped
-            arch_mapping: Path to the arechetype mappings
+            arch_mapping: Path to the archetype mappings
 
         Returns:
             ret_json: JSON of the results ordered by cluster and category.
@@ -143,7 +148,7 @@ def create_mapped_dmg_result(inventory, dmg_result, arch_mapping, groupby_col_na
         """
         dmg_states = dmg_result[groupby_col_name].unique().tolist()  # get unique damage states
         dmg_merged = pd.merge(inventory, dmg_result, on='guid')
-        mapped_df = pd.merge(dmg_merged, arch_mapping, on='archetype')
+        mapped_df = pd.merge(dmg_merged, arch_mapping, on=arch_col)
         unique_categories = arch_mapping.groupby(by=['cluster', 'category'], sort=False).count().reset_index()
 
         group_by = mapped_df.groupby(by=[groupby_col_name, 'cluster', 'category']).count().reset_index()
@@ -174,13 +179,14 @@ def create_mapped_dmg_result(inventory, dmg_result, arch_mapping, groupby_col_na
         return {"by_cluster": json_by_cluster, "by_category": json_by_category}
 
     @staticmethod
-    def create_mapped_func_result(inventory, bldg_func, arch_mapping):
+    def create_mapped_func_result(inventory, bldg_func, arch_mapping, arch_col='archetype'):
         """
 
         Args:
             inventory: dataframe represent inventory
             bldg_func: building func dataset
-            arch_mapping: Path to the arechetype mappings
+            arch_mapping: Path to the archetype mappings
+            arch_col: archetype column to use for the clustering
 
         Returns:
             ret_json: JSON of the results ordered by cluster and category.
@@ -188,7 +194,7 @@ def create_mapped_func_result(inventory, bldg_func, arch_mapping):
         """
         func_state = ["percent_functional", "percent_non_functional", "num_functional", "num_non_functional"]
         func_merged = pd.merge(inventory, bldg_func, on='guid')
-        mapped_df = pd.merge(func_merged, arch_mapping, on='archetype')
+        mapped_df = pd.merge(func_merged, arch_mapping, on=arch_col)
         unique_categories = arch_mapping.groupby(by=['category'], sort=False, as_index=False).count()['category']
         unique_cluster = arch_mapping.groupby(by=['cluster', 'category'], sort=False, as_index=False).count()[[
             'cluster', 'category']]
diff --git a/pyincore/utils/hhrsoutputprocess.py b/pyincore/utils/hhrsoutputprocess.py
new file mode 100644
index 000000000..593e1959b
--- /dev/null
+++ b/pyincore/utils/hhrsoutputprocess.py
@@ -0,0 +1,36 @@
+# Copyright (c) 2021 University of Illinois and others. All rights reserved.
+#
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+import json
+
+
+class HHRSOutputProcess:
+    """This class converts results outputs of housing household recovery sequential analysis to certain format."""
+
+    @staticmethod
+    def get_hhrs_stage_count(timesteps, hhrs_df, filename_json="hhrs_stage_count.json"):
+        """
+
+        Args:
+            timesteps: timesteps in the unit of month ["0", "6", "12", "24", "48", etc]
+            hhrs_df: pandas dataframe of the output of housingrecoverysequential
+            filename_json: the name of the json file to store the output
+
+        Returns:
+
+        """
+
+        hhrs_stage_count = {}
+        for t in timesteps:
+            stage = hhrs_df[t]
+            hhrs_stage_count[t] = [int((stage == 1.0).sum()), int((stage == 2.0).sum()), int((stage == 3.0).sum()),
+                                   int((stage == 4.0).sum()), int((stage == 5.0).sum())]
+
+        if filename_json:
+            with open(filename_json, "w") as outfile:
+                json.dump(hhrs_stage_count, outfile, indent=2)
+
+        return hhrs_stage_count
+
diff --git a/recipes/meta.yaml b/recipes/meta.yaml
index c9d912b54..1c8d673f8 100644
--- a/recipes/meta.yaml
+++ b/recipes/meta.yaml
@@ -1,5 +1,5 @@
 {% set name = "pyincore" %}
-{% set version = "1.8.0" %}
+{% set version = "1.9.0" %}
 
 package:
   name: {{ name|lower }}
diff --git a/scripts/release-packages.yml b/scripts/release-packages.yml
index 3b727a881..5126ea33a 100644
--- a/scripts/release-packages.yml
+++ b/scripts/release-packages.yml
@@ -1,7 +1,7 @@
 # List of pyincore analyses to either include or exclude for release
 # If provided, 'excludes' takes preference over 'includeOnly'
 
-version: '1.8.0'
+version: '1.9.0'
 
 includeOnly:
   - buildingdamage
diff --git a/setup.py b/setup.py
index a0e3dda48..3148d3321 100644
--- a/setup.py
+++ b/setup.py
@@ -7,7 +7,7 @@
 from setuptools import setup, find_packages
 
 # version number of pyincore
-version = '1.8.0'
+version = '1.9.0'
 
 with open("README.rst", encoding="utf-8") as f:
     readme = f.read()


Test patch:
diff --git a/tests/pyincore/analyses/combinedwindwavesurgebuildingloss/test_combinedwindwavesurgebuildingloss.py b/tests/pyincore/analyses/combinedwindwavesurgebuildingloss/test_combinedwindwavesurgebuildingloss.py
new file mode 100644
index 000000000..7a6c6c38e
--- /dev/null
+++ b/tests/pyincore/analyses/combinedwindwavesurgebuildingloss/test_combinedwindwavesurgebuildingloss.py
@@ -0,0 +1,87 @@
+from pyincore import IncoreClient, FragilityService, MappingSet
+from pyincore.analyses.combinedwindwavesurgebuildingloss import CombinedWindWaveSurgeBuildingLoss
+from pyincore.analyses.buildingdamage import BuildingDamage
+import pyincore.globals as pyglobals
+from timeit import default_timer as timer
+
+
+def run_with_base_class():
+
+    client = IncoreClient(pyglobals.INCORE_API_DEV_URL)
+
+    # Galveston building inventory
+    bldg_dataset_id = "63053e74c8f8b7614f6cff9f"
+
+    # Hazard information
+    hazard_id = "5fa472033c1f0c73fe81461a"
+    hazard_type = "hurricane"
+
+    # Fragility Service
+    fragility_service = FragilityService(client)
+
+    # Surge-wave mapping
+    mapping_id = "6303e51bd76c6d0e1f6be080"
+    mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))
+
+    # surge-wave building damage
+    sw_bldg_dmg = BuildingDamage(client)
+    sw_bldg_dmg.load_remote_input_dataset("buildings", bldg_dataset_id)
+    sw_bldg_dmg.set_input_dataset('dfr3_mapping_set', mapping_set)
+    sw_bldg_dmg.set_parameter("result_name", "Galveston-sw-dmg")
+    sw_bldg_dmg.set_parameter("hazard_type", hazard_type)
+    sw_bldg_dmg.set_parameter("hazard_id", hazard_id)
+    sw_bldg_dmg.set_parameter("num_cpu", 4)
+    sw_bldg_dmg.run_analysis()
+
+    # wind mapping
+    mapping_id = "62fef3a6cef2881193f2261d"
+    mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))
+
+    # wind building damage
+    w_bldg_dmg = BuildingDamage(client)
+    w_bldg_dmg.load_remote_input_dataset("buildings", bldg_dataset_id)
+    w_bldg_dmg.set_input_dataset('dfr3_mapping_set', mapping_set)
+    w_bldg_dmg.set_parameter("result_name", "Galveston-wind-dmg")
+    w_bldg_dmg.set_parameter("hazard_type", hazard_type)
+    w_bldg_dmg.set_parameter("hazard_id", hazard_id)
+    w_bldg_dmg.set_parameter("num_cpu", 4)
+    w_bldg_dmg.run_analysis()
+
+    # flood mapping
+    mapping_id = "62fefd688a30d30dac57bbd7"
+    mapping_set = MappingSet(fragility_service.get_mapping(mapping_id))
+
+    # flood building damage
+    f_bldg_dmg = BuildingDamage(client)
+    f_bldg_dmg.load_remote_input_dataset("buildings", bldg_dataset_id)
+    f_bldg_dmg.set_input_dataset('dfr3_mapping_set', mapping_set)
+    f_bldg_dmg.set_parameter("result_name", "Galveston-flood-dmg")
+    f_bldg_dmg.set_parameter("hazard_type", hazard_type)
+    f_bldg_dmg.set_parameter("hazard_id", hazard_id)
+    f_bldg_dmg.set_parameter("num_cpu", 4)
+    f_bldg_dmg.run_analysis()
+
+    # Get damage outputs from different hazards
+    surge_wave_damage = sw_bldg_dmg.get_output_dataset("ds_result")
+    wind_damage = w_bldg_dmg.get_output_dataset("ds_result")
+    flood_damage = f_bldg_dmg.get_output_dataset("ds_result")
+
+    start = timer()
+
+    # Combined building loss from each hazard
+    combined_bldg_loss = CombinedWindWaveSurgeBuildingLoss(client)
+    combined_bldg_loss.load_remote_input_dataset("buildings", bldg_dataset_id)
+    combined_bldg_loss.set_input_dataset("surge_wave_damage", surge_wave_damage)
+    combined_bldg_loss.set_input_dataset("wind_damage", wind_damage)
+    combined_bldg_loss.set_input_dataset("flood_damage", flood_damage)
+    combined_bldg_loss.load_remote_input_dataset("structural_cost", "63fd15716d3b2a308ba914c8")
+    combined_bldg_loss.load_remote_input_dataset("content_cost", "63fd16956d3b2a308ba9269a")
+    combined_bldg_loss.set_parameter("result_name", "Galveston")
+    combined_bldg_loss.run_analysis()
+
+    end = timer()
+    print(f'Elapsed time: {end - start:.3f} seconds')
+
+
+if __name__ == '__main__':
+    run_with_base_class()
diff --git a/tests/pyincore/analyses/housingrecoveryserial/test_housingrecoverysequential.py b/tests/pyincore/analyses/housingrecoveryserial/test_housingrecoverysequential.py
index d4f058f05..348439b50 100644
--- a/tests/pyincore/analyses/housingrecoveryserial/test_housingrecoverysequential.py
+++ b/tests/pyincore/analyses/housingrecoveryserial/test_housingrecoverysequential.py
@@ -2,14 +2,17 @@
 # terms of the Mozilla Public License v2.0 which accompanies this distribution,
 # and is available at https://www.mozilla.org/en-US/MPL/2.0/
 
-from pyincore import IncoreClient, FragilityService, MappingSet
+from pyincore import IncoreClient, HHRSOutputProcess
 from pyincore.analyses.housingrecoverysequential import HousingRecoverySequential
+from timeit import default_timer as timer
 import pyincore.globals as pyglobals
 
 
 def run_with_base_class():
     client = IncoreClient(pyglobals.INCORE_API_DEV_URL)
 
+    start = timer()
+
     population_dislocation = "60f098f502897f12fcda19ec"  # dev Galveston testbed
     transition_probability_matrix = "60ef513802897f12fcd9765c"
     initial_probability_vector = "60ef532e02897f12fcd9ac63"
@@ -23,9 +26,11 @@ def run_with_base_class():
     housing_recovery = HousingRecoverySequential(client)
 
     # Parameter setup
+    housing_recovery.set_parameter('num_cpu', 4)
     housing_recovery.set_parameter('seed', seed)
     housing_recovery.set_parameter('t_delta', t_delta)
     housing_recovery.set_parameter('t_final', t_final)
+    housing_recovery.set_parameter('result_name', 'results_hhrs_galveston.csv')
 
     # Dataset inputs
     housing_recovery.load_remote_input_dataset("population_dislocation_block", population_dislocation)
@@ -35,6 +40,18 @@ def run_with_base_class():
 
     housing_recovery.run()
 
+    # test the utilility
+    hhr_result = housing_recovery.get_output_dataset("ds_result")
+    hhrs_df = hhr_result.get_dataframe_from_csv()
+
+    end = timer()
+
+    print(f'Elapsed time: {end - start:.3f} seconds')
+
+    timesteps = ["1", "7", "13", "25", "49"]  # t0, t6, t12, t24, t48
+
+    print(HHRSOutputProcess.get_hhrs_stage_count(timesteps, hhrs_df))
+
 
 if __name__ == '__main__':
     run_with_base_class()
diff --git a/tests/pyincore/analyses/saltlakecge/test_saltlakecge.py b/tests/pyincore/analyses/saltlakecge/test_saltlakecge.py
new file mode 100644
index 000000000..64c14c183
--- /dev/null
+++ b/tests/pyincore/analyses/saltlakecge/test_saltlakecge.py
@@ -0,0 +1,32 @@
+from pyincore import IncoreClient, Dataset
+from pyincore.analyses.saltlakecge import SaltLakeCGEModel
+import pyincore.globals as pyglobals
+
+# This script runs SaltLakeCGEModel analysis with input files from
+# IN-CORE development services.
+
+
+def run_base_analysis():
+    client = IncoreClient(pyglobals.INCORE_API_DEV_URL)
+    saltlake_cge = SaltLakeCGEModel(client)
+
+    saltlake_cge.set_parameter("model_iterations", 1)
+
+    saltlake_cge.load_remote_input_dataset("SAM", "640758d66121f943887299a2")
+    saltlake_cge.load_remote_input_dataset("BB", "640759b56121f94388729d01")
+    saltlake_cge.load_remote_input_dataset("MISCH", "64075a5d6121f94388729f8a")
+    saltlake_cge.load_remote_input_dataset("EMPLOY", "64075b116121f94388729f91")
+    saltlake_cge.load_remote_input_dataset("JOBCR", "64075d1c6121f9438872a648")
+    saltlake_cge.load_remote_input_dataset("OUTCR", "64075e306121f9438872a7fb")
+
+    # This sector shock should work and print results
+    saltlake_cge.load_remote_input_dataset("sector_shocks", "64075ec46121f9438872a802")
+
+    # Bad shock data - should produce infeasible result and not write output files
+    # saltlake_cge.load_remote_input_dataset("sector_shocks", "640b94348bec9f280fdf9164")
+
+    saltlake_cge.run_analysis()
+
+
+if __name__ == '__main__':
+    run_base_analysis()
diff --git a/tests/pyincore/test_hazardservice.py b/tests/pyincore/test_hazardservice.py
index 9b0ce0c98..9ec375e4c 100644
--- a/tests/pyincore/test_hazardservice.py
+++ b/tests/pyincore/test_hazardservice.py
@@ -427,12 +427,15 @@ def test_search_floods(hazardsvc):
 
 
 @pytest.mark.skip(reason="performance issues")
-def test_create_hurricane_windfield(hazardsvc):
+def test_create_and_delete_hurricane_windfield(hazardsvc):
     with open(os.path.join(pyglobals.TEST_DATA_DIR, "hurricanewf.json"), 'r') as file:
         hurr_wf_inputs = file.read()
 
-    response = hazardsvc.create_hurricane_windfield(hurr_wf_inputs)
-    assert response["id"] is not None
+    post_response = hazardsvc.create_hurricane_windfield(hurr_wf_inputs)
+    assert post_response["id"] is not None
+
+    del_response = hazardsvc.delete_hurricanewf(post_response["id"])
+    assert del_response["id"] is not None
 
 
 @pytest.mark.skip(reason="performance issues")
diff --git a/tests/pyincore/test_spaceservice.py b/tests/pyincore/test_spaceservice.py
index 5469fd285..52b4bd88d 100644
--- a/tests/pyincore/test_spaceservice.py
+++ b/tests/pyincore/test_spaceservice.py
@@ -55,3 +55,20 @@ def test_get_space(spacesvc):
     space_id = "5c813be55648c42a9168d5c1"
     space = spacesvc.get_space_by_id(space_id)
     assert "id" in space and space["id"] == space_id
+
+
+def test_get_space_by_name(spacesvc):
+    space_name = "incore"
+    space = spacesvc.get_space_by_name(space_name)
+    assert "id" in space[0]
+
+
+def test_add_to_space_by_name(spacesvc):
+    space_name = "incore"
+    dataset_id = "640128d5b40ec85799d46039"
+    space = spacesvc.add_to_space_by_name(space_name, dataset_id)
+    assert "id" in space and dataset_id in space["members"]
+
+    # delete the dataset id from the space after testing add
+    space = spacesvc.remove_from_space_by_name(space_name, dataset_id)
+    assert "id" in space and dataset_id not in space["members"]
diff --git a/tests/pyincore/utils/test_cgecsvoutputjson.py b/tests/pyincore/utils/test_cgecsvoutputjson.py
index 637f4af5a..24d681593 100644
--- a/tests/pyincore/utils/test_cgecsvoutputjson.py
+++ b/tests/pyincore/utils/test_cgecsvoutputjson.py
@@ -87,6 +87,44 @@ def run_convert_cge_json_path(testpath):
     return True
 
 
+def run_convert_SLC_cge_json_path(testpath):
+    # test the external file with a path
+
+    cge_json = CGEOutputProcess()
+    region = ["R1", "R2", "R3", "R4", "R5", "R6", "R7"]
+
+    categories = []
+    for h in ["HH1", "HH2", "HH3", "HH4"]:
+        for r in region:
+            categories.append(h + "_" + r)
+
+    cge_json.get_cge_household_count(None,
+                                     os.path.join(testpath, "household-count.csv"),
+                                     "slc_cge_total_household_count.json", income_categories=categories)
+    cge_json.get_cge_gross_income(None,
+                                  os.path.join(testpath, "gross-income.csv"),
+                                  "slc_cge_total_household_income.json", income_categories=categories)
+
+    categories = []
+    for d in ["AG_MI", "UTIL", "CONS", "MANU", "COMMER", "EDU", "HEALTH", "ART_ACC", "RELIG", "L1W", "L2W", "L3W",
+              "L4W", "HS1", "HS2", "HH1", "HH2", "HH3", "HH4", "USSOC1", "USSOC2", "USSOC3", "USSOC4"]:
+        for r in region:
+            categories.append(d + "_" + r)
+    cge_json.get_cge_employment(None, None, os.path.join(testpath, "pre-disaster-factor-demand.csv"),
+                                os.path.join(testpath, "post-disaster-factor-demand.csv"),
+                                "slc_cge_employment.json", demand_categories=categories)
+
+    categories = []
+    for d in ["AG_MI", "UTIL", "CONS", "MANU", "COMMER", "EDU", "HEALTH", "ART_ACC", "RELIG",
+              "HS1", "HS2"]:
+        for r in region:
+            categories.append(d + "_" + r)
+    cge_json.get_cge_domestic_supply(None,
+                                     os.path.join(testpath, "domestic-supply.csv"), "slc_cge_domestic_supply.json",
+                                     supply_categories=categories)
+    return True
+
+
 if __name__ == '__main__':
     # test chaining with Joplin CGE analysis
     run_convert_cge_json_chained()
@@ -97,4 +135,6 @@ def run_convert_cge_json_path(testpath):
     if testpath:
         run_convert_cge_json_path(testpath)
 
+    run_convert_SLC_cge_json_path(testpath="../../data/SLC_CGE")
+
     print("DONE")
diff --git a/tests/pyincore/utils/test_csvoutputjson.py b/tests/pyincore/utils/test_csvoutputjson.py
index 790c3bedb..6d5632762 100644
--- a/tests/pyincore/utils/test_csvoutputjson.py
+++ b/tests/pyincore/utils/test_csvoutputjson.py
@@ -5,7 +5,7 @@
 from pyincore.utils.cgeoutputprocess import CGEOutputProcess
 import os
 
-PYINCOREPATH = "path-to-pyincore"
+PYINCOREPATH = ""
 TESTAPATH = "pyincore/tests/pyincore/analyses/"
 
 def run_convert_cge_json():
@@ -14,18 +14,24 @@ def run_convert_cge_json():
     filepath = os.path.join(PYINCOREPATH, TESTAPATH, "joplincge")
 
     cge_json.get_cge_household_count(None,
-                                     os.path.join(filepath, "joplin-pop-disl-results.csv"),
-                                     "cge_total_household_count.json")
+                                     os.path.join(filepath, "household-count.csv"),
+                                     "cge_total_household_count.json",
+                                     ["HH1", "HH2", "HH3", "HH4", "HH5"]
+                                     )
     cge_json.get_cge_gross_income(None,
                                   os.path.join(filepath, "gross-income.csv"),
-                                  "cge_total_household_income.json")
+                                  "cge_total_household_income.json",
+                                  ["HH1", "HH2", "HH3", "HH4", "HH5"])
     cge_json.get_cge_employment(None, None,
                                 os.path.join(filepath, "pre-disaster-factor-demand.csv"),
                                 os.path.join(filepath, "post-disaster-factor-demand.csv"),
-                                "cge_employment.json")
+                                "cge_employment.json",
+                                ["GOODS", "TRADE", "OTHER"]
+                                )
     cge_json.get_cge_domestic_supply(None,
                                      os.path.join(filepath, "domestic-supply.csv"),
-                                     "cge_domestic_supply.json")
+                                     "cge_domestic_supply.json",
+                                     ["Goods", "trade", "other", "HS1", "HS2", "HS3"])
     return True