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

Issue:
pin pandas < 2
this will pandas to be version < 2


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index c5c508ca6..2750def85 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,18 @@ 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.10.0] - 2023-04-21
+
+### Added
+- Added Galveston Capital Shock and CGE models as a submodule[#239](https://github.com/IN-CORE/pyincore/issues/239)
+
+### Fixed
+- CGE output post process util function [#298](https://github.com/IN-CORE/pyincore/issues/298)
+- Population Dislocation utility function arbitrarily assumes there will be dislocated and non-dislocated [#301](https://github.com/IN-CORE/pyincore/issues/301)
+- Seaside & Joplin cge uses a fixed location for temporary files [#312](https://github.com/IN-CORE/pyincore/issues/312)
+- Functional vs non-functional calculation based of failure sample now [#300](https://github.com/IN-CORE/pyincore/issues/300)
+- Exposed timeout and kwargs parameter for incore client methods [#295](https://github.com/IN-CORE/pyincore/issues/295)
+
 ## [1.9.0] - 2023-03-15
 
 ### Added
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 5bb583fdd..2800bd4ad 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -33,9 +33,9 @@
 author = ''
 
 # The short X.Y version
-version = '1.9'
+version = '1.10'
 # The full version, including alpha/beta/rc tags
-release = '1.9.0'
+release = '1.10.0'
 
 # -- General configuration ---------------------------------------------------
 
diff --git a/docs/source/modules.rst b/docs/source/modules.rst
index 2b278a61d..b0992fcb8 100644
--- a/docs/source/modules.rst
+++ b/docs/source/modules.rst
@@ -80,6 +80,29 @@ analyses/example
 ..  autoclass:: example.exampleanalysis.ExampleAnalysis
     :members:
 
+analyses/galvestoncge
+==================
+..  autoclass:: galvestoncge.galvestoncge.GalvestonCGEModel
+    :members:
+..  autoclass:: galvestoncge.equationlib.VarContainer
+    :members:
+..  autoclass:: galvestoncge.equationlib.Variable
+    :members:
+..  autoclass:: galvestoncge.equationlib.ExprItem
+    :members:
+..  autoclass:: galvestoncge.equationlib.Expr
+    :members:
+..  autoclass:: galvestoncge.equationlib.ExprM
+    :members:
+..  autofunction:: galvestoncge.outputfunctions.baseValue
+    :members:
+..  autofunction:: galvestoncge.outputfunctions.newValue
+    :members:
+..  autofunction:: galvestoncge.outputfunctions.getDiff
+    :members:
+..  autofunction:: galvestoncge.galvestonoutput.gams_to_dataframes
+    :members:
+
 analyses/housingrecovery
 ========================
 ..  autoclass:: housingrecovery.housingrecovery.HousingRecovery
@@ -206,6 +229,29 @@ analyses/seasidecge
 ..  autoclass:: seasidecge.solverconstants.equationlib.ExprM
     :members:
 
+analyses/saltlakecge
+==================
+..  autoclass:: saltlakecge.saltlakecge.SaltLakeCGEModel
+    :members:
+..  autoclass:: saltlakecge.equationlib.VarContainer
+    :members:
+..  autoclass:: saltlakecge.equationlib.Variable
+    :members:
+..  autoclass:: saltlakecge.equationlib.ExprItem
+    :members:
+..  autoclass:: saltlakecge.equationlib.Expr
+    :members:
+..  autoclass:: saltlakecge.equationlib.ExprM
+    :members:
+..  autofunction:: saltlakecge.outputfunctions.baseValue
+    :members:
+..  autofunction:: saltlakecge.outputfunctions.newValue
+    :members:
+..  autofunction:: saltlakecge.outputfunctions.getDiff
+    :members:
+..  autofunction:: saltlakecge.saltlakecge.gams_to_dataframes
+    :members:
+
 analyses/socialvulnerability
 ============================
 ..  autoclass:: socialvulnerability.socialvulnerability.SocialVulnerability
diff --git a/environment.yml b/environment.yml
index e28ace1f1..99dc8f8d0 100644
--- a/environment.yml
+++ b/environment.yml
@@ -9,7 +9,7 @@ dependencies:
   - matplotlib>=2.1.0
   - networkx>=2.2
   - numpy>=1.16.6,<2.0a0
-  - pandas>=0.24.1
+  - pandas>=0.24.1,<2
   - pycodestyle>=2.6.0
   - pyomo>=5.6
   - pyproj>=1.9.6
diff --git a/pyincore/analyses/galvestoncge/__init__.py b/pyincore/analyses/galvestoncge/__init__.py
new file mode 100644
index 000000000..6a77dda3d
--- /dev/null
+++ b/pyincore/analyses/galvestoncge/__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.galvestoncge.galvestoncge import GalvestonCGEModel
+from pyincore.analyses.galvestoncge.equationlib import *
diff --git a/pyincore/analyses/galvestoncge/equationlib.py b/pyincore/analyses/galvestoncge/equationlib.py
new file mode 100644
index 000000000..f5d39a0ae
--- /dev/null
+++ b/pyincore/analyses/galvestoncge/equationlib.py
@@ -0,0 +1,846 @@
+# 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="float64").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="float64").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/galvestoncge/galvestoncge.py b/pyincore/analyses/galvestoncge/galvestoncge.py
new file mode 100644
index 000000000..9c10efd15
--- /dev/null
+++ b/pyincore/analyses/galvestoncge/galvestoncge.py
@@ -0,0 +1,2118 @@
+# 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 sys
+from pyincore import globals as pyglobals
+from pyincore import BaseAnalysis
+from pyincore.analyses.galvestoncge.equationlib import *
+from pyincore.analyses.galvestoncge.outputfunctions import *
+from pyincore.analyses.galvestoncge.galvestonoutput import gams_to_dataframes
+
+import os
+import pandas as pd
+from pyomo.environ import *
+from pyomo.opt import SolverFactory
+
+class GalvestonCGEModel(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(GalvestonCGEModel, 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.galveston_cge(iNum, sam, bb, jobcr, misch, employ, outcr, sector_shocks)
+
+    def galveston_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 = ['IAGMIN', 'IUTIL', 'ICONS', 'IMANU', 'IWHOLE', 'IRETAIL', 'ITRANS', 'IPROFSER', 'IREALE', 'IEDU',
+             'IHEALTH', 'IART', 'IACCO',
+             'MAGMIN', 'MUTIL', 'MCONS', 'MMANU', 'MWHOLE', 'MRETAIL', 'MTRANS', 'MPROFSER', 'MREALE', 'MEDU',
+             'MHEALTH', 'MART', 'MACCO',
+             'IL1', 'IL2', 'IL3', 'IL4',
+             'ML1', 'ML2', 'ML3', 'ML4',
+             'KAP', 'INVES',
+             'HS1I', 'HS2I', 'HS3I', 'HS1M', 'HS2M', 'HS3M',
+             'IHH1', 'IHH2', 'IHH3', 'IHH4', 'IHH5',
+             'MHH1', 'MHH2', 'MHH3', 'MHH4', 'MHH5',
+             'IUSSOC1', 'IUSSOC2', 'IUSSOC3', 'IUSSOC4',
+             'MUSSOC1', 'MUSSOC2', 'MUSSOC3', 'MUSSOC4',
+             'PROTAX', 'TXSTX', 'COUNTYSTX', 'USPIT', 'OTX', 'CYGF', 'COUNTY', 'STATE', 'FED',
+             'OUT1', 'OUT2', 'OUT3', 'OUT4', 'ROW']
+
+        # FACTORS
+        F = ['IL1', 'IL2', 'IL3', 'IL4', 'ML1', 'ML2', 'ML3', 'ML4', 'KAP']
+        FW1 = ['IL1']
+        FW2 = ['IL2']
+        FW3 = ['IL3']
+        FW4 = ['IL4']
+
+        # OUT COMMUTERS
+        CM = ['OUT1', 'OUT2', 'OUT3', 'OUT4']
+        CM1 = ['OUT1']
+        CM2 = ['OUT2']
+        CM3 = ['OUT3']
+        CM4 = ['OUT4']
+
+        # LABOR SUPPLY
+        LT = ['IL1', 'IL2', 'IL3', 'IL4', 'ML1', 'ML2', 'ML3', 'ML4', 'OUT1', 'OUT2', 'OUT3', 'OUT4']
+
+        # LABOR WORKING IN GALVESTON
+        L = ['IL1', 'IL2', 'IL3', 'IL4', 'ML1', 'ML2', 'ML3', 'ML4']
+
+        # CAPITAL
+        K = ['KAP']
+
+        # LAND
+        # LA = ['LAND']
+
+        # GOVERNMENTS
+        G = ['IUSSOC1', 'IUSSOC2', 'IUSSOC3', 'IUSSOC4', 'MUSSOC1', 'MUSSOC2', 'MUSSOC3', 'MUSSOC4',
+             'PROTAX', 'TXSTX', 'COUNTYSTX', 'USPIT', 'OTX', 'CYGF', 'COUNTY', 'STATE', 'FED']
+
+        # ENDOGENOUS GOVERNMENTS
+        GN = ['COUNTY', 'STATE', 'FED']
+
+        # LOCAL ENDOGENOUS GOVERNMENTS
+        GNL = ['COUNTY']
+
+        # EXOGENOUS GOVERMENTS
+        GX = ['IUSSOC1', 'IUSSOC2', 'IUSSOC3', 'IUSSOC4', 'MUSSOC1', 'MUSSOC2', 'MUSSOC3', 'MUSSOC4',
+              'PROTAX', 'TXSTX', 'COUNTYSTX', 'USPIT', 'OTX']
+
+        # SALES OR EXCISE TAXES
+        GS = ['TXSTX', 'COUNTYSTX', 'OTX']
+
+        # LOCAL TAX
+        GC = ['PROTAX', 'TXSTX', 'COUNTYSTX', 'OTX']
+
+        # LAND TAXES
+        GL = ['PROTAX']
+
+        # FACTOR TAXES
+        GF = ['IUSSOC1', 'IUSSOC2', 'IUSSOC3', 'IUSSOC4', 'MUSSOC1', 'MUSSOC2', 'MUSSOC3', 'MUSSOC4',
+              'PROTAX']
+
+        # INCOME TAX UNITS
+        GI = ['USPIT']
+
+        # HOUSEHOLD TAX UNITS
+        GH = ['PROTAX']
+
+        # SOCIAL SECURITY PAYMENT
+        GFUS = ['IUSSOC1', 'IUSSOC2', 'IUSSOC3', 'IUSSOC4', 'MUSSOC1', 'MUSSOC2', 'MUSSOC3', 'MUSSOC4']
+
+        # EXOGNOUS TRANSFER PMT
+        GY = ['IUSSOC1', 'IUSSOC2', 'IUSSOC3', 'IUSSOC4', 'MUSSOC1', 'MUSSOC2', 'MUSSOC3', 'MUSSOC4',
+              'PROTAX', 'TXSTX', 'COUNTYSTX', 'USPIT', 'OTX', 'CYGF', 'COUNTY', 'STATE', 'FED']
+
+        # EXOGNOUS TRANSFER PMT
+        GTA = ['IUSSOC1', 'IUSSOC2', 'IUSSOC3', 'IUSSOC4', 'MUSSOC1', 'MUSSOC2', 'MUSSOC3', 'MUSSOC4',
+               'PROTAX', 'TXSTX', 'COUNTYSTX', 'USPIT', 'OTX', 'CYGF', 'COUNTY', 'STATE', 'FED']
+
+        # ENDOGENOUS TRANSFER PMT
+        GT = ['CYGF', 'STATE', 'FED']
+
+        # ALL HOUSEHOLDS IN GALVESTON
+        H = ['IHH1', 'IHH2', 'IHH3', 'IHH4', 'IHH5', 'MHH1', 'MHH2', 'MHH3', 'MHH4', 'MHH5']
+
+        # HOUSEHOLD (INCOME i, i = 1,2,3,4,5)
+        HH1 = ['IHH1', 'MHH1']
+        HH2 = ['IHH2', 'MHH2']
+        HH3 = ['IHH3', 'MHH3']
+        HH4 = ['IHH4', 'MHH4']
+        HH5 = ['IHH5', 'MHH5']
+
+        # I+G SECTORS
+        IG = ['IAGMIN', 'IUTIL', 'ICONS', 'IMANU', 'IWHOLE', 'IRETAIL', 'ITRANS', 'IPROFSER', 'IREALE', 'IEDU',
+              'IHEALTH', 'IART', 'IACCO',
+              'MAGMIN', 'MUTIL', 'MCONS', 'MMANU', 'MWHOLE', 'MRETAIL', 'MTRANS', 'MPROFSER', 'MREALE', 'MEDU',
+              'MHEALTH', 'MART', 'MACCO',
+              'HS1I', 'HS2I', 'HS3I', 'HS1M', 'HS2M', 'HS3M', 'COUNTY', 'STATE', 'FED']
+
+        # INDUSTRY SECTORS
+        I = ['IAGMIN', 'IUTIL', 'ICONS', 'IMANU', 'IWHOLE', 'IRETAIL', 'ITRANS', 'IPROFSER', 'IREALE', 'IEDU',
+             'IHEALTH', 'IART', 'IACCO',
+             'MAGMIN', 'MUTIL', 'MCONS', 'MMANU', 'MWHOLE', 'MRETAIL', 'MTRANS', 'MPROFSER', 'MREALE', 'MEDU',
+             'MHEALTH', 'MART', 'MACCO',
+             'HS1I', 'HS2I', 'HS3I', 'HS1M', 'HS2M', 'HS3M']
+
+        # ENDOGENOUS GOVERNMENTS
+        IG2 = ['COUNTY', 'STATE', 'FED']
+
+        # PRODUCTION SECTORS
+        IP = ['IAGMIN', 'IUTIL', 'ICONS', 'IMANU', 'IWHOLE', 'IRETAIL', 'ITRANS', 'IPROFSER', 'IREALE', 'IEDU',
+              'IHEALTH', 'IART', 'IACCO',
+              'MAGMIN', 'MUTIL', 'MCONS', 'MMANU', 'MWHOLE', 'MRETAIL', 'MTRANS', 'MPROFSER', 'MREALE', 'MEDU',
+              'MHEALTH', 'MART', 'MACCO']
+
+        # PRODUCTION GOV
+        FG = ['COUNTY', 'STATE', 'FED']
+
+        # HOUSING SERVICE DEMAND
+        HSD = ['HS1I', 'HS2I', 'HS3I', 'HS1M', 'HS2M', 'HS3M']
+
+        # HOUSING SERVICES
+        HSSET = ['HS1I', 'HS2I', 'HS3I', 'HS1M', 'HS2M', 'HS3M']
+
+        # HOUSING SERVICES 2 & 3
+        HS23SET = ['HS2I', 'HS3I', 'HS2M', 'HS3M']
+
+        # 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
+
+        # ----------------------------------------------------------------
+        # IMPORT ADDITIONAL DATA FILES
+        # ----------------------------------------------------------------
+
+        # SOCIAL ACCOUNTING MATRIX
+        SAM = sam
+
+        # CAPITAL COMP
+        BB = bb
+
+        # MISC TABLES
+        JOBCR = jobcr
+        MISCH = misch
+        EMPLOY = employ
+        OUTCR = outcr
+
+        # ----------------------------------------------------------------
+        # 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)
+        LAMBDAH = pd.DataFrame(index=HSD, columns=HSD).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)
+        PH0 = pd.Series(index=HSD, 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', 'TXSTX'] = 1
+        # IGTD.loc['STATE', 'UTITX'] = 1
+
+        TPC.loc[H, G] = 0
+        TPC.loc[H, GFUS] = 1
+
+        TAUFF.loc[G, F] = 0
+        TAUFF.loc['IUSSOC1', 'IL1'] = 1
+        TAUFF.loc['IUSSOC2', 'IL2'] = 1
+        TAUFF.loc['IUSSOC3', 'IL3'] = 1
+        TAUFF.loc['IUSSOC4', 'IL4'] = 1
+
+        TAUFF.loc['MUSSOC1', 'ML1'] = 1
+        TAUFF.loc['MUSSOC2', 'ML2'] = 1
+        TAUFF.loc['MUSSOC3', 'ML3'] = 1
+        TAUFF.loc['MUSSOC4', 'ML4'] = 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[HH5, '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[HH5, 'ETAYD'] = 1.8
+
+        MISCH.loc[HH1, 'ETAU'] = -0.8
+        MISCH.loc[HH2, 'ETAU'] = -0.8
+        MISCH.loc[HH3, 'ETAU'] = -0.8
+        MISCH.loc[HH4, 'ETAU'] = -0.6
+        MISCH.loc[HH5, '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']
+
+        LAMBDAH.loc[HSD, HSD] = MISC.loc[HSD, '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
+        PH0.loc[HSD] = 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['IL1', G] = SAM.loc['IL1', G] + SAM.loc['IUSSOC1', G]
+        AGFS.loc['IL2', G] = SAM.loc['IL2', G] + SAM.loc['IUSSOC2', G]
+        AGFS.loc['IL3', G] = SAM.loc['IL3', G] + SAM.loc['IUSSOC3', G]
+        AGFS.loc['IL4', G] = SAM.loc['IL4', G] + SAM.loc['IUSSOC4', G]
+
+        AGFS.loc['ML1', G] = SAM.loc['ML1', G] + SAM.loc['MUSSOC1', G]
+        AGFS.loc['ML2', G] = SAM.loc['ML2', G] + SAM.loc['MUSSOC2', G]
+        AGFS.loc['ML3', G] = SAM.loc['ML3', G] + SAM.loc['MUSSOC3', G]
+        AGFS.loc['ML4', G] = SAM.loc['ML4', G] + SAM.loc['MUSSOC4', 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)
+
+        # COMMUTING OUT
+        '''
+        CMO0.loc[['OUT1']] = 5498
+        CMO0.loc[['OUT2']] = 8635
+        CMO0.loc[['OUT3']] = 5660
+        CMO0.loc[['OUT4']] = 270
+        '''
+        # COMMUTING IN
+        '''
+        CMI0.loc[['L1']] = FD0.loc[['L1'], IG].sum(1) - JOBCOR.loc[H, ['L1']].mul(HW0.loc[H], axis = 'index').sum(0)
+        CMI0.loc[['L2']] = FD0.loc[['L2'], IG].sum(1) - JOBCOR.loc[H, ['L2']].mul(HW0.loc[H], axis = 'index').sum(0)
+        CMI0.loc[['L3']] = FD0.loc[['L3'], IG].sum(1) - JOBCOR.loc[H, ['L3']].mul(HW0.loc[H], axis = 'index').sum(0)
+        CMI0.loc[['L4']] = FD0.loc[['L4'], IG].sum(1) - JOBCOR.loc[H, ['L4']].mul(HW0.loc[H], axis = 'index').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))/8)/ (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))/32 *(SUM(CM, CMO(CM))/(SUM(H1, HW(H1)* SUM(L,JOBCOR(H1,L)))+ SUM(CM, CMO(CM)) +SUM(L,CMI(L)))))** (ETARA(H)*0.1)
+            #                              * ( 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('IL1') ))** 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('IL2') ))** 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('IL3') ))** 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('IL4') ))** 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))/10))** ECOMI(L)));
+            print('LSEQ3')
+            line = ExprM(vars, m=CMI0.loc[L]) * (
+                        (RA.loc(L) / (CPI.loc(H).sum(H) / 10)) ** (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)
+
+            '''
+            #GSEQJ1('CYGF').. S('CYGF')=E= Y('CYGF') -  Y('CYGF');
+            print('GSEQJ1(\'CYGF\')')
+            GSEQJ1 = S.loc(['CYGF']) - Y.loc(['CYGF']) + Y.loc(['CYGF'])
+            GSEQJ1.write(count,filename)
+            print(GSEQJ1.test(vars.initialVals))
+            '''
+
+            #  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 = '/Users/mohanar2/opt/anaconda3/envs/pyincore/bin/ipopt'
+            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'
+            # opt.options['line_search_method'] = 'penalty'
+            # opt.options['watchdog_shortened_iter_trigger'] = 0
+
+            ###
+
+            ### 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. ########  
+        '''
+        # iNum = 1 # dynamic model itterations
+        sims = sector_shocks
+        all_sectors = ['IAGMIN', 'IUTIL', 'ICONS', 'IMANU', 'IWHOLE', 'IRETAIL', 'ITRANS', 'IPROFSER', 'IREALE', 'IEDU',
+             'IHEALTH', 'IART', 'IACCO',
+             'MAGMIN', 'MUTIL', 'MCONS', 'MMANU', 'MWHOLE', 'MRETAIL', 'MTRANS', 'MPROFSER', 'MREALE', 'MEDU',
+             'MHEALTH', 'MART', 'MACCO',
+             'HS1I', 'HS2I', 'HS3I', 'HS1M', 'HS2M', 'HS3M']
+
+
+        missing_sectors = list(set(all_sectors) - set(sims.index.unique()))
+
+        for i in missing_sectors:
+            data = {'shock':1}
+            index = [i]
+            df = pd.DataFrame(data=data, index=index)
+
+            sims = sims.append(df)
+
+
+        # === begin replacing the average shocks on housing services
+        HSmean = sims.loc[['HS1I', 'HS2I', 'HS3I']].mean().mean()
+        sims.loc[['HS1I', 'HS2I', 'HS3I']] = HSmean
+        # === end replacing the average shocks on housing services
+
+        iNum = 1
+        # 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': 'Galveston-cge',
+            'description': 'CGE model for Galveston.',
+            '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 Galveston 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/galvestoncge/galvestonoutput.py b/pyincore/analyses/galvestoncge/galvestonoutput.py
new file mode 100644
index 000000000..83bba5e10
--- /dev/null
+++ b/pyincore/analyses/galvestoncge/galvestonoutput.py
@@ -0,0 +1,225 @@
+"""
+# After running a model, use solution dictionaries to replicate GAMS
+# results with DataFrames.
+#
+# NOTE:
+# Solution variables stored in "vars" and "soln" objects should
+# be the primary source for model evaluation.
+# This was only created to assist those unfamiliar with python objects.
+"""
+
+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'
+    domestic_supply_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'
+    gross_income_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'
+    household_count_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'
+    pre_disaster_demand_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'
+    post_disaster_demand_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/analyses/galvestoncge/outputfunctions.py b/pyincore/analyses/galvestoncge/outputfunctions.py
new file mode 100644
index 000000000..b47a3afc4
--- /dev/null
+++ b/pyincore/analyses/galvestoncge/outputfunctions.py
@@ -0,0 +1,72 @@
+# 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.galvestoncge.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
+
+
diff --git a/pyincore/analyses/joplincge/joplincge.py b/pyincore/analyses/joplincge/joplincge.py
index 442d8efd7..6cd8ffa1b 100644
--- a/pyincore/analyses/joplincge/joplincge.py
+++ b/pyincore/analyses/joplincge/joplincge.py
@@ -7,7 +7,6 @@
 from pyomo.opt import SolverFactory
 
 import os
-import tempfile
 import pandas as pd
 
 logger = pyglobals.LOGGER
@@ -1752,15 +1751,9 @@ def run_solver(cons_filename, temp_file_name):
         # TODO: we need to generate the "solverconstatnt" folder with username since it uses system tmp
         # TODO: there is a situation that multiple users on system can run this together
 
-        cge_tmp_folder = os.path.join(tempfile.gettempdir(), "solverconstants")
-        if not os.path.isdir(cge_tmp_folder):  # create the folder if there is no folder
-            os.mkdir(cge_tmp_folder)
-        logger.debug(cge_tmp_folder)
-
-        filename = os.path.join(cge_tmp_folder, "ipopt_cons.py")
-        tmp = os.path.join(cge_tmp_folder, "tmp.py")
-
-        # print("Calibration: ")
+        filename = "ipopt_cons.py"
+        tmp = "tmp.py"
+        print("Calibration: ")
         run_solver(filename, tmp)
 
         '''
diff --git a/pyincore/analyses/saltlakecge/equationlib.py b/pyincore/analyses/saltlakecge/equationlib.py
index 686981c93..4bef6d0a1 100644
--- a/pyincore/analyses/saltlakecge/equationlib.py
+++ b/pyincore/analyses/saltlakecge/equationlib.py
@@ -561,7 +561,6 @@ def clear_empty(self):
         #     self.itemList.append(ExprItem(1))
 
 
-
 class ExprM:
     """
     Three ways to create a ExprMatrix:
diff --git a/pyincore/analyses/saltlakecge/outputfunctions.py b/pyincore/analyses/saltlakecge/outputfunctions.py
index 7c579ea22..6d4b205b2 100644
--- a/pyincore/analyses/saltlakecge/outputfunctions.py
+++ b/pyincore/analyses/saltlakecge/outputfunctions.py
@@ -7,6 +7,7 @@
 import pandas as pd
 import numpy as np
 
+
 # Return initial values function
 def baseValue(vars, soln, eqName):
     if vars.getInfo(eqName)['size'] == 1:
@@ -66,4 +67,4 @@ def newValue(vars, soln, eqName, ittr):
 # 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
+    return difference
diff --git a/pyincore/analyses/seasidecge/seasidecge.py b/pyincore/analyses/seasidecge/seasidecge.py
index 487bd903b..d9f6e677f 100644
--- a/pyincore/analyses/seasidecge/seasidecge.py
+++ b/pyincore/analyses/seasidecge/seasidecge.py
@@ -5,7 +5,6 @@
 # and is available at https://www.mozilla.org/en-US/MPL/2.0/
 
 import os
-import tempfile
 
 from pyomo.environ import *
 from pyomo.opt import SolverFactory
@@ -1817,15 +1816,9 @@ def run_solver(cons_filename, temp_file_name):
         # TODO: we need to generate the "solverconstatnt" folder with username since it uses system tmp
         # TODO: there is a situation that multiple users on system can run this together
 
-        cge_tmp_folder = os.path.join(tempfile.gettempdir(), "solverconstants")
-        if not os.path.isdir(cge_tmp_folder):  # create the folder if there is no folder
-            os.mkdir(cge_tmp_folder)
-        logger.debug(cge_tmp_folder)
-
-        filename = os.path.join(cge_tmp_folder, "ipopt_cons.py")
-        tmp = os.path.join(cge_tmp_folder, "tmp.py")
-
-        logger.debug("Calibration: ")
+        filename = "ipopt_cons.py"
+        tmp = "tmp.py"
+        print("Calibration: ")
         run_solver(filename, tmp)
 
         '''
diff --git a/pyincore/dataservice.py b/pyincore/dataservice.py
index eca329efd..dbf4c8472 100644
--- a/pyincore/dataservice.py
+++ b/pyincore/dataservice.py
@@ -32,41 +32,46 @@ def __init__(self, client: IncoreClient):
         self.base_earthquake_url = urljoin(client.service_url, 'hazard/api/earthquakes/')
         self.base_tornado_url = urljoin(client.service_url, 'hazard/api/tornadoes/')
 
-    def get_dataset_metadata(self, dataset_id: str):
+    def get_dataset_metadata(self, dataset_id: str, timeout=(30, 600), **kwargs):
         """Retrieve metadata from data service. Dataset API endpoint is called.
 
         Args:
             dataset_id (str): ID of the Dataset.
-
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
         Returns:
             obj: HTTP response containing the metadata.
 
         """
         # construct url with service, dataset api, and id
         url = urljoin(self.base_url, dataset_id)
-        r = self.client.get(url)
+        r = self.client.get(url, timeout=timeout, **kwargs)
         return r.json()
 
-    def get_dataset_files_metadata(self, dataset_id: str):
+    def get_dataset_files_metadata(self, dataset_id: str, timeout=(30, 600), **kwargs):
         """Retrieve metadata of all files associated with the dataset. Files API endpoint is called.
 
         Args:
             dataset_id (str): ID of the Dataset.
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
 
         Returns:
             obj: HTTP response containing the metadata.
 
         """
         url = urljoin(self.base_url, dataset_id + '/files')
-        r = self.client.get(url)
+        r = self.client.get(url, timeout=timeout, **kwargs)
         return r.json()
 
-    def get_dataset_file_metadata(self, dataset_id: str, file_id: str):
+    def get_dataset_file_metadata(self, dataset_id: str, file_id: str, timeout=(30, 600), **kwargs):
         """Retrieve metadata of all files associated with the dataset. Files API endpoint is called.
 
         Args:
             dataset_id (str): ID of the Dataset.
             file_id (str): ID of the File.
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
 
         Returns:
             obj: HTTP response containing the metadata.
@@ -74,15 +79,17 @@ def get_dataset_file_metadata(self, dataset_id: str, file_id: str):
         """
         url = urljoin(self.base_url,
                                    dataset_id + "/files/" + file_id)
-        r = self.client.get(url)
+        r = self.client.get(url, timeout=timeout, **kwargs)
         return r.json()
 
-    def get_dataset_blob(self, dataset_id: str, join=None):
+    def get_dataset_blob(self, dataset_id: str, join=None, timeout=(30, 600), **kwargs):
         """Retrieve a blob of the dataset. Blob API endpoint is called.
 
         Args:
             dataset_id (str): ID of the Dataset.
             join (bool): Add join parameter if True. Default None.
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
 
         Returns:
             str: Folder or file name.
@@ -99,7 +106,7 @@ def get_dataset_blob(self, dataset_id: str, join=None):
             # for consistency check to ensure the repository hash is recorded in service.json
             self.client.create_service_json_entry()
 
-            local_filename = self.download_dataset_blob(cache_data_dir, dataset_id)
+            local_filename = self.download_dataset_blob(cache_data_dir, dataset_id, timeout=timeout, **kwargs )
 
         # if cache_data_dir exist, check if id folder and zip file exist inside
         else:
@@ -108,7 +115,7 @@ def get_dataset_blob(self, dataset_id: str, join=None):
                     local_filename = os.path.join(cache_data_dir, fname)
                     print('Dataset already exists locally. Reading from local cached zip.')
             if not local_filename:
-                local_filename = self.download_dataset_blob(cache_data_dir, dataset_id)
+                local_filename = self.download_dataset_blob(cache_data_dir, dataset_id, timeout=timeout, **kwargs)
 
         folder = self.unzip_dataset(local_filename)
         if folder is not None:
@@ -116,19 +123,19 @@ def get_dataset_blob(self, dataset_id: str, join=None):
         else:
             return local_filename
 
-    def download_dataset_blob(self, cache_data_dir: str, dataset_id: str, join=None):
+    def download_dataset_blob(self, cache_data_dir: str, dataset_id: str, join=None,timeout=(30, 600), **kwargs):
         # construct url for file download
         url = urljoin(self.base_url, dataset_id + '/blob')
-        kwargs = {"stream": True}
+        kwargs["stream"] = True
         if join is None:
-            r = self.client.get(url, **kwargs)
+            r = self.client.get(url, timeout=timeout, **kwargs)
         else:
             payload = {}
             if join is True:
                 payload['join'] = 'true'
             elif join is False:
                 payload['join'] = 'false'
-            r = self.client.get(url, params=payload, **kwargs)
+            r = self.client.get(url, params=payload, timeout=timeout, **kwargs)
 
         # extract filename
         disposition = r.headers['content-disposition']
@@ -145,7 +152,7 @@ def download_dataset_blob(self, cache_data_dir: str, dataset_id: str, join=None)
         return local_filename
 
     def get_datasets(self, datatype: str = None, title: str = None, creator: str = None, skip: int = None,
-                     limit: int = None, space: str = None):
+                     limit: int = None, space: str = None, timeout=(30, 600), **kwargs):
         """Function to get datasets. Blob API endpoint is called.
 
         Args:
@@ -156,6 +163,9 @@ def get_datasets(self, datatype: str = None, title: str = None, creator: str = N
             skip (str): Skip the first n results, default None.
             limit (int): Limit number of results to return, default None.
             space (str): Name of space, default None.
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
+
 
         Returns:
             obj: HTTP response containing the metadata.
@@ -176,33 +186,36 @@ def get_datasets(self, datatype: str = None, title: str = None, creator: str = N
         if space is not None:
             payload['space'] = space
 
-        r = self.client.get(url, params=payload)
+        r = self.client.get(url, params=payload, timeout=timeout, **kwargs)
         # need to handle there is no datasets
         return r.json()
 
-    def create_dataset(self, properties: dict):
+    def create_dataset(self, properties: dict, timeout=(30, 600), **kwargs):
         """Create datasets. Post API endpoint is called.
 
         Args:
             properties (dict): Metadata and files of the dataset to be created.
-
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
         Returns:
             obj: HTTP POST Response. Json of the dataset posted to the server.
 
         """
         payload = {'dataset': json.dumps(properties)}
         url = self.base_url
-        kwargs = {"files": payload}
-        r = self.client.post(url, **kwargs)
+        kwargs["files"] = payload
+        r = self.client.post(url, timeout=timeout, **kwargs)
         return r.json()
 
     def update_dataset(self, dataset_id, property_name: str,
-                       property_value: str):
+                       property_value: str, timeout=(30, 600), **kwargs):
         """Update dataset. Put API endpoint is called.
 
         Args:
             dataset_id (str): ID of the Dataset.
             property_name (str): Property parameters such as name and value.
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
 
         Returns:
             obj: HTTP PUT Response. Json of the dataset updated on the server.
@@ -211,16 +224,18 @@ def update_dataset(self, dataset_id, property_name: str,
         url = urljoin(self.base_url, dataset_id)
         payload = {'update': json.dumps({"property name": property_name,
                                          "property value": property_value})}
-        kwargs = {"files": payload}
-        r = self.client.put(url, **kwargs)
+        kwargs["files"] = payload
+        r = self.client.put(url, timeout=timeout, **kwargs)
         return r.json()
 
-    def add_files_to_dataset(self, dataset_id: str, filepaths: list):
+    def add_files_to_dataset(self, dataset_id: str, filepaths: list, timeout=(30, 600), **kwargs):
         """Add files to the dataset. Post API endpoint is called.
 
         Args:
             dataset_id (str): ID of the Dataset.
             filepaths (list): Path to the files.
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
 
         Returns:
             obj: HTTP POST Response. Json of the files added to the dataset.
@@ -232,8 +247,8 @@ def add_files_to_dataset(self, dataset_id: str, filepaths: list):
             file = open(filepath, 'rb')
             tuple = ('file', file)
             listfiles.append(tuple)
-        kwargs = {"files": listfiles}
-        r = self.client.post(url, **kwargs)
+        kwargs['files'] = listfiles
+        r = self.client.post(url, timeout=timeout, **kwargs)
 
         # close files
         for tuple in listfiles:
@@ -242,7 +257,7 @@ def add_files_to_dataset(self, dataset_id: str, filepaths: list):
         return r.json()
 
     def add_files_to_network_dataset(self, dataset_id: str, filepaths: list,
-                                     nodename: str, linkname: str, graphname: str):
+                                     nodename: str, linkname: str, graphname: str, timeout=(30, 600), **kwargs):
         """Add files to the network dataset. Post API endpoint is called.
 
         Args:
@@ -251,6 +266,8 @@ def add_files_to_network_dataset(self, dataset_id: str, filepaths: list,
             nodename (str): node shapefile name.
             linkname (str): link shapefile name.
             graphname (str): graph file name.
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
 
         Returns:
             obj: HTTP POST Response. Json of the files added to the dataset.
@@ -273,8 +290,8 @@ def add_files_to_network_dataset(self, dataset_id: str, filepaths: list,
                 bodyname = 'graph-file'
             tuple = (bodyname, file)
             listfiles.append(tuple)
-        kwargs = {"files": listfiles}
-        r = self.client.post(url, **kwargs)
+        kwargs["files"] = listfiles
+        r = self.client.post(url, timeout=timeout, **kwargs)
 
         # close files
         for tuple in listfiles:
@@ -282,50 +299,58 @@ def add_files_to_network_dataset(self, dataset_id: str, filepaths: list,
 
         return r.json()
 
-    def delete_dataset(self, dataset_id: str):
+    def delete_dataset(self, dataset_id: str, timeout=(30, 600), **kwargs):
         """Delete dataset. Delete API endpoint is called.
 
         Args:
             dataset_id (str): ID of the Dataset to be deleted.
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
 
         Returns:
             obj: HTTP DELETE Response. Json model of the delete action.
 
         """
         url = urljoin(self.base_url, dataset_id)
-        r = self.client.delete(url)
+        r = self.client.delete(url, timeout=timeout, **kwargs)
         return r.json()
 
-    def get_files(self):
+    def get_files(self, timeout=(30, 600), **kwargs):
         """Get all files. Files API endpoint is called.
-
+        Args:
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
         Returns:
             obj: HTTP response containing the files.
 
         """
         url = self.files_url
-        r = self.client.get(url)
+        r = self.client.get(url, timeout=timeout, **kwargs)
         return r.json()
 
-    def get_file_metadata(self, file_id: str):
+    def get_file_metadata(self, file_id: str, timeout=(30, 600), **kwargs):
         """Function to retrieve metadata of a file defined by id. Files API endpoint is called.
 
         Args:
             file_id (str): ID of the File.
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
 
         Returns:
             obj: HTTP response containing the metadata.
 
         """
         url = urljoin(self.files_url, file_id)
-        r = self.client.get(url)
+        r = self.client.get(url, timeout=timeout, **kwargs)
         return r.json()
 
-    def get_file_blob(self, file_id: str):
+    def get_file_blob(self, file_id: str, timeout=(30, 600), **kwargs):
         """Function to retrieve a blob of the file. Blob API endpoint is called.
 
         Args:
             file_id (str): ID of the Dataset.
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
 
         Returns:
             str: Local file name.
@@ -333,8 +358,8 @@ def get_file_blob(self, file_id: str):
         """
         # construct url for file download
         url = urljoin(self.files_url, file_id + '/blob')
-        kwargs = {"stream": True}
-        r = self.client.get(url, **kwargs)
+        kwargs["stream"] = True
+        r = self.client.get(url, timeout=timeout, **kwargs)
 
         # extract filename
         disposition = r.headers['content-disposition']
@@ -379,12 +404,14 @@ def unzip_dataset(self, local_filename: str):
         zip_ref.close()
         return foldername
 
-    def get_shpfile_from_service(self, fileid, dirname):
+    def get_shpfile_from_service(self, fileid, dirname, timeout=(30, 600), **kwargs):
         """Function to obtain a shape file from Data service.
 
         Args:
             fileid (str):  ID of the Shape file.
             dirname (str): Directory the files are being extracted.
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
 
         Returns:
             str: Filename with shape files.
@@ -394,23 +421,25 @@ def get_shpfile_from_service(self, fileid, dirname):
         request_str_zip = request_str + '/blob'
 
         # obtain file name
-        r = self.client.get(request_str)
+        r = self.client.get(request_str, timeout=timeout, **kwargs)
         first_filename = r.json()['fileDescriptors'][0]['filename']
         filename = os.path.splitext(first_filename)[0]
-        kwargs = {"stream": True}
+        kwargs["stream"] = True
 
-        r = self.client.get(request_str_zip, **kwargs)
+        r = self.client.get(request_str_zip, timeout=timeout, **kwargs)
         z = zipfile.ZipFile(io.BytesIO(r.content))
         z.extractall(dirname)
         # print(r.status_code)
 
         return filename
 
-    def get_tornado_dataset_id_from_service(self, fileid):
+    def get_tornado_dataset_id_from_service(self, fileid, timeout=(30, 600), **kwargs):
         """Function to obtain a tornado dataset Id from Data service.
 
         Args:
             fileid (str): ID of the Tornado file.
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
 
         Returns:
             obj: HTTP response containing the tornado Dataset ID.
@@ -418,17 +447,19 @@ def get_tornado_dataset_id_from_service(self, fileid):
         """
         # dataset
         request_str = self.base_tornado_url + fileid
-        r = self.client.get(request_str)
+        r = self.client.get(request_str, timeout=timeout, **kwargs)
 
         return r.json()['tornadoDatasetId']
 
-    def search_datasets(self, text: str, skip: int = None, limit: int = None):
+    def search_datasets(self, text: str, skip: int = None, limit: int = None, timeout=(30, 600), **kwargs):
         """Function to search datasets.
 
         Args:
             text (str): Text to search by, passed to the parameter "text".
             skip (int):  Skip the first n results, passed to the parameter "skip". Dafault None.
             limit (int):  Limit number of results to return. Passed to the parameter "limit". Dafault None.
+            timeout (tuple[int,int]): Session timeout.
+            **kwargs: A dictionary of external parameters.
 
         Returns:
             obj: HTTP response with search results.
@@ -441,6 +472,6 @@ def search_datasets(self, text: str, skip: int = None, limit: int = None):
         if limit is not None:
             payload['limit'] = limit
 
-        r = self.client.get(url, params=payload)
+        r = self.client.get(url, params=payload, timeout=timeout, **kwargs)
 
         return r.json()
diff --git a/pyincore/globals.py b/pyincore/globals.py
index f17efacbc..32e4799da 100644
--- a/pyincore/globals.py
+++ b/pyincore/globals.py
@@ -10,7 +10,7 @@
 import os
 import shutil
 
-PACKAGE_VERSION = "1.9.0"
+PACKAGE_VERSION = "1.10.0"
 
 INCORE_API_PROD_URL = "https://incore.ncsa.illinois.edu"
 INCORE_API_DEV_URL = "https://incore-dev.ncsa.illinois.edu"
diff --git a/pyincore/utils/cgeoutputprocess.py b/pyincore/utils/cgeoutputprocess.py
index e34018215..4a90dfd19 100644
--- a/pyincore/utils/cgeoutputprocess.py
+++ b/pyincore/utils/cgeoutputprocess.py
@@ -39,19 +39,20 @@ def get_cge_household_count(household_count, household_count_path=None, filename
         else:
             household_group_count = household_count.get_dataframe_from_csv()
 
-        before_values = household_group_count["HH0"]
-        after_values = household_group_count["HHL"]
-
         before_event = {}
         after_event = {}
         pct_change = {}
-        for i in range(len(income_categories)):
-            before_event[income_categories[i]] = before_values[i]
-            after_event[income_categories[i]] = after_values[i]
-            if before_values[i]:
-                pct_change[income_categories[i]] = 100 * ((after_values[i] - before_values[i]) / abs(before_values[i]))
+        for income_category in income_categories:
+            before_event[income_category] = household_group_count[household_group_count["Household Group"] ==
+                                                                  income_category]["HH0"].values[0]
+            after_event[income_category] = household_group_count[household_group_count["Household Group"] ==
+                                                                 income_category]["HHL"].values[0]
+
+            if before_event[income_category]:
+                pct_change[income_category] = 100 * ((after_event[income_category] - before_event[income_category]) /
+                                                     abs(before_event[income_category]))
             else:
-                pct_change[income_categories[i]] = None
+                pct_change[income_category] = None
 
         cge_total_household_count = {"beforeEvent": before_event, "afterEvent": after_event, "%_change": pct_change}
 
@@ -89,19 +90,20 @@ def get_cge_gross_income(gross_income, gross_income_path=None, filename_json=Non
         else:
             household_income = gross_income.get_dataframe_from_csv()
 
-        before_values = household_income["Y0"]
-        after_values = household_income["YL"]
-
         before_event = {}
         after_event = {}
         pct_change = {}
-        for i in range(len(income_categories)):
-            before_event[income_categories[i]] = before_values[i]
-            after_event[income_categories[i]] = after_values[i]
-            if before_values[i]:
-                pct_change[income_categories[i]] = 100 * ((after_values[i] - before_values[i]) / abs(before_values[i]))
+        for income_category in income_categories:
+            before_event[income_category] = household_income[household_income["Household Group"] == income_category][
+                "Y0"].values[0]
+            after_event[income_category] = household_income[household_income["Household Group"] == income_category][
+                "YL"].values[0]
+
+            if before_event[income_category]:
+                pct_change[income_category] = 100 * ((after_event[income_category] - before_event[income_category]) /
+                                                     abs(before_event[income_category]))
             else:
-                pct_change[income_categories[i]] = None
+                pct_change[income_category] = None
 
         cge_total_household_income = {"beforeEvent": before_event, "afterEvent": after_event, "%_change": pct_change}
 
@@ -150,19 +152,18 @@ def get_cge_employment(pre_demand, post_demand, pre_demand_path=None, post_deman
             pre_disaster_demand = pre_demand.get_dataframe_from_csv()
             post_disaster_demand = post_demand.get_dataframe_from_csv()
 
-        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 = {}
         pct_change = {}
-        for i in range(len(demand_categories)):
-            before_event[demand_categories[i]] = before_values[i]
-            after_event[demand_categories[i]] = after_values[i]
-            if before_values[i]:
-                pct_change[demand_categories[i]] = 100 * ((after_values[i] - before_values[i]) / abs(before_values[i]))
+        for demand_category in demand_categories:
+            before_event[demand_category] = pre_disaster_demand[demand_category].sum()
+            after_event[demand_category] = post_disaster_demand[demand_category].sum()
+
+            if before_event[demand_category]:
+                pct_change[demand_category] = 100 * ((after_event[demand_category] - before_event[demand_category]) /
+                                                     abs(before_event[demand_category]))
             else:
-                pct_change[demand_categories[i]] = None
+                pct_change[demand_category] = None
 
         cge_employment = {"beforeEvent": before_event, "afterEvent": after_event, "%_change": pct_change}
 
@@ -203,19 +204,17 @@ def get_cge_domestic_supply(domestic_supply, domestic_supply_path=None, filename
         else:
             sector_supply = domestic_supply.get_dataframe_from_csv()
 
-        before_values = sector_supply["DS0"]
-        after_values = sector_supply["DSL"]
-
         before_event = {}
         after_event = {}
         pct_change = {}
-        for i in range(len(supply_categories)):
-            before_event[supply_categories[i]] = before_values[i]
-            after_event[supply_categories[i]] = after_values[i]
-            if before_values[i]:
-                pct_change[supply_categories[i]] = 100 * ((after_values[i] - before_values[i]) / abs(before_values[i]))
+        for supply_category in supply_categories:
+            before_event[supply_category] = sector_supply[sector_supply["Sectors"] == supply_category]["DS0"].values[0]
+            after_event[supply_category] = sector_supply[sector_supply["Sectors"] == supply_category]["DSL"].values[0]
+            if before_event[supply_category]:
+                pct_change[supply_category] = 100 * ((after_event[supply_category] - before_event[supply_category]) /
+                                                     abs(before_event[supply_category]))
             else:
-                pct_change[supply_categories[i]] = None
+                pct_change[supply_category] = None
 
         cge_domestic_supply = {"beforeEvent": before_event, "afterEvent": after_event, "%_change": pct_change}
 
diff --git a/pyincore/utils/dataprocessutil.py b/pyincore/utils/dataprocessutil.py
index 6724650f8..7208762d2 100644
--- a/pyincore/utils/dataprocessutil.py
+++ b/pyincore/utils/dataprocessutil.py
@@ -7,8 +7,10 @@
 import geopandas as gpd
 import json
 import pandas as pd
+import numpy as np
 
-from pyincore import Dataset, DataService
+from pyincore import Dataset, DataService, IncoreClient
+from functools import reduce
 
 
 class DataProcessUtil:
@@ -16,7 +18,7 @@ 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", arch_col='archetype'
+                                        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.
@@ -43,7 +45,7 @@ def get_mapped_result_from_analysis(client, inventory_id: str, dmg_result_datase
         dmg_result = dmg_result_dataset.get_dataframe_from_csv()
 
         bldg_func_df = bldg_func_dataset.get_dataframe_from_csv()
-        bldg_func_df.rename(columns={'building_guid': 'guid'}, inplace=True)
+        bldg_func_df.rename(columns={"building_guid": "guid", "samples": "failure"}, inplace=True)
 
         arch_mapping = Dataset.from_data_service(archetype_mapping_id, DataService(client)).get_dataframe_from_csv()
 
@@ -57,7 +59,7 @@ def get_mapped_result_from_analysis(client, inventory_id: str, dmg_result_datase
     @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", arch_col='archetype'):
+                                          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.
 
@@ -84,7 +86,7 @@ def get_mapped_result_from_dataset_id(client, inventory_id: str, dmg_result_id:
 
         bldg_func_dataset = Dataset.from_data_service(bldg_func_id, DataService(client))
         bldg_func_df = bldg_func_dataset.get_dataframe_from_csv()
-        bldg_func_df.rename(columns={'building_guid': 'guid'}, inplace=True)
+        bldg_func_df.rename(columns={"building_guid": "guid", "samples": "failure"}, inplace=True)
 
         archtype_mapping_dataset = Dataset.from_data_service(archetype_mapping_id, DataService(client))
         arch_mapping = archtype_mapping_dataset.get_dataframe_from_csv()
@@ -100,7 +102,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, arch_col='archetype'):
+                                    groupby_col_name: str, arch_col="archetype"):
         """
 
         Args:
@@ -121,7 +123,7 @@ def get_mapped_result_from_path(inventory_path: str, dmg_result_path: str,
         inventory = pd.DataFrame(gpd.read_file("zip://" + inventory_path))
         dmg_result = pd.read_csv(dmg_result_path)
         bldg_func_df = pd.read_csv(func_result_path)
-        bldg_func_df.rename(columns={'building_guid': 'guid'}, inplace=True)
+        bldg_func_df.rename(columns={"building_guid": "guid", "samples": "failure"}, inplace=True)
         arch_mapping = pd.read_csv(archetype_mapping_path)
 
         max_state_df = DataProcessUtil.get_max_damage_state(dmg_result)
@@ -134,7 +136,7 @@ def get_mapped_result_from_path(inventory_path: str, dmg_result_path: str,
 
     @staticmethod
     def create_mapped_dmg_result(inventory, dmg_result, arch_mapping, groupby_col_name="max_state",
-                                 arch_col='archetype'):
+                                 arch_col="archetype"):
         """
 
         Args:
@@ -147,26 +149,26 @@ 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')
+        dmg_merged = pd.merge(inventory, dmg_result, on="guid")
         mapped_df = pd.merge(dmg_merged, arch_mapping, on=arch_col)
-        unique_categories = arch_mapping.groupby(by=['cluster', 'category'], sort=False).count().reset_index()
+        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()
-        group_by = group_by.loc[:, ['guid', groupby_col_name, 'cluster', 'category']]
-        group_by.rename(columns={'guid': 'count'}, inplace=True)
+        group_by = mapped_df.groupby(by=[groupby_col_name, "cluster", "category"]).count().reset_index()
+        group_by = group_by.loc[:, ["guid", groupby_col_name, "cluster", "category"]]
+        group_by.rename(columns={"guid": "count"}, inplace=True)
 
-        pivot = group_by.pivot_table(values='count', index=['cluster', 'category'], columns=groupby_col_name,
+        pivot = group_by.pivot_table(values="count", index=["cluster", "category"], columns=groupby_col_name,
                                      fill_value=0)
 
         table = pd.DataFrame()
-        table[['category', 'cluster']] = unique_categories[['category', 'cluster']]
-        result_by_cluster = pd.merge(table, pivot, how='left', on=['cluster', 'category'])
+        table[["category", "cluster"]] = unique_categories[["category", "cluster"]]
+        result_by_cluster = pd.merge(table, pivot, how="left", on=["cluster", "category"])
 
         # Add missing max damage states. Handles case when no inventory fall under some damage states.
         result_by_cluster = result_by_cluster.reindex(result_by_cluster.columns.union(
             dmg_states, sort=False), axis=1, fill_value=0)
 
-        result_by_category = result_by_cluster.groupby(by=['category'], sort=False).sum(min_count=1).reset_index()
+        result_by_category = result_by_cluster.groupby(by=["category"], sort=False).sum(min_count=1).reset_index()
 
         result_by_cluster[dmg_states] = result_by_cluster[dmg_states].fillna(-1).astype(int)
         result_by_category[dmg_states] = result_by_category[dmg_states].fillna(-1).astype(int)
@@ -179,12 +181,12 @@ 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, arch_col='archetype'):
+    def create_mapped_func_result(inventory, bldg_func, arch_mapping, arch_col="archetype"):
         """
 
         Args:
             inventory: dataframe represent inventory
-            bldg_func: building func dataset
+            bldg_func: building func state dataset
             arch_mapping: Path to the archetype mappings
             arch_col: archetype column to use for the clustering
 
@@ -192,49 +194,55 @@ def create_mapped_func_result(inventory, bldg_func, arch_mapping, arch_col='arch
             ret_json: JSON of the results ordered by cluster and category.
 
         """
+        def _sum_average(series):
+            return reduce(lambda x, y: np.mean(x + y).round(0), series)
+
         func_state = ["percent_functional", "percent_non_functional", "num_functional", "num_non_functional"]
-        func_merged = pd.merge(inventory, bldg_func, on='guid')
+
+        # unify mcs and bldg func naming
+        bldg_func.rename(columns={"building_guid": "guid", "samples": "failure"}, inplace=True)
+
+        func_merged = pd.merge(inventory, bldg_func, on="guid")
         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']]
-
-        # group by cluster
-        result_by_cluster = mapped_df.groupby(by=['cluster', 'category'], sort=False, as_index=False).agg(
-            {'guid': 'count',
-             'probability': 'mean'})
-        result_by_cluster.rename(columns={'guid': 'tot_count', 'probability': 'percent_functional'}, inplace=True)
-        result_by_cluster["percent_non_functional"] = 1 - result_by_cluster["percent_functional"]
-        result_by_cluster["num_functional"] = (result_by_cluster["tot_count"] * result_by_cluster[
-            "percent_functional"]).round(0)
-        result_by_cluster["num_non_functional"] = (result_by_cluster["tot_count"] * result_by_cluster[
-            "percent_non_functional"]).round(0)
-        result_by_cluster = result_by_cluster.drop('tot_count', 1)
-        result_by_cluster = pd.merge(unique_cluster, result_by_cluster, how='left', on=['cluster', 'category'])
-        # Add missing max damage states. Handles case when no inventory fall under some damage states.
-        result_by_cluster = result_by_cluster.reindex(result_by_cluster.columns.union(
-            func_state, sort=False), axis=1, fill_value=0)
-        # replace NaN
-        result_by_cluster[func_state] = result_by_cluster[func_state].fillna(-1)
-        result_by_cluster[["num_functional", "num_non_functional"]] = result_by_cluster[["num_functional",
-                                                                                         "num_non_functional"]].astype(
-            int)
-
-        # group by category
-        result_by_category = mapped_df.groupby(by=['category'], sort=False, as_index=False).agg({'guid': 'count',
-                                                                                                 'probability': 'mean'})
-        result_by_category.rename(columns={'guid': 'tot_count', 'probability': 'percent_functional'}, inplace=True)
-        result_by_category["percent_non_functional"] = 1 - result_by_category["percent_functional"]
-        result_by_category["num_functional"] = (
-                result_by_category["tot_count"] * result_by_category["percent_functional"]).round(0)
-        result_by_category["num_non_functional"] = (
-                result_by_category["tot_count"] * result_by_category["percent_non_functional"]).round(0)
-        result_by_category = result_by_category.drop('tot_count', 1)
-        result_by_category = pd.merge(unique_categories, result_by_category, how='left', on=['category'])
-        # replace NaN
-        result_by_category[func_state] = result_by_category[func_state].fillna(-1)
-        result_by_category[["num_functional", "num_non_functional"]] = result_by_category[
-            ["num_functional", "num_non_functional"]].astype(int)
+        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"]]
+
+        mapped_df = mapped_df[["guid", "failure", "category", "cluster"]]
+        mapped_df["failure_array"] = mapped_df["failure"].apply(lambda x: np.array([int(x) for x in x.split(",")]))
+
+        def _group_by(by_column, unique):
+            # group by cluster
+            result = mapped_df.groupby(by=by_column, sort=False, as_index=False).agg(
+                {"guid": "count", "failure_array": [_sum_average]})
+
+            # clean up
+            result.rename(columns={"guid": "tot_count", "failure_array": "num_functional"}, inplace=True)
+
+            # 0 (failed), 1 (not failed). MCS
+            # 0 otherwise (not functional), 1 (functional),  Functionality
+            result["num_non_functional"] = result["tot_count"].squeeze() - result["num_functional"].squeeze()
+            result["percent_functional"] = result["num_functional"].squeeze() / result["tot_count"].squeeze()
+            result["percent_non_functional"] = 1 - result["percent_functional"]
+
+            # remove the tuples in column
+            result.columns = [x[0] if len(x) > 1 else x for x in result.columns]
+
+            # more clean up
+            result = pd.merge(unique, result, how="left", on=by_column)
+
+            # Add missing max damage states. Handles case when no inventory fall under some damage states.
+            result = result.reindex(result.columns.union(func_state, sort=False), axis=1, fill_value=0)
+
+            # replace NaN
+            result[func_state] = result[func_state].fillna(-1)
+            result["tot_count"] = result["tot_count"].fillna(-1)
+            result[["num_functional", "num_non_functional"]] = result[["num_functional", "num_non_functional"]].astype(int)
+
+            return result
+
+        result_by_cluster = _group_by(by_column=["cluster", "category"], unique=unique_cluster)
+        result_by_category = _group_by(by_column=["category"], unique=unique_categories)
 
         cluster_records = result_by_cluster.to_json(orient="records")
         category_records = result_by_category.to_json(orient="records")
@@ -254,20 +262,20 @@ def get_max_damage_state(dmg_result):
             pd.DataFrame: Pandas dataframe that has column GUID and column max_state.
 
         """
-        if all(column in dmg_result.columns for column in ['DS_0', 'DS_1', 'DS_2', 'DS_3']):
-            dmg_states = ['DS_0', 'DS_1', 'DS_2', 'DS_3']
-        elif all(column in dmg_result.columns for column in ['insignific', 'moderate', 'heavy', 'complete']):
-            dmg_states = ['insignific', 'moderate', 'heavy', 'complete']
+        if all(column in dmg_result.columns for column in ["DS_0", "DS_1", "DS_2", "DS_3"]):
+            dmg_states = ["DS_0", "DS_1", "DS_2", "DS_3"]
+        elif all(column in dmg_result.columns for column in ["insignific", "moderate", "heavy", "complete"]):
+            dmg_states = ["insignific", "moderate", "heavy", "complete"]
         elif all(column in dmg_result.columns for column in ["ds-none", "ds-slight", "ds-moderat", "ds-extensi",
                                                              "ds-complet"]):
             dmg_states = ["ds-none", "ds-slight", "ds-moderat", "ds-extensi", "ds-complet"]
         else:
             raise ValueError("Invalid damage state names. Cannot create mapped max damage state.")
 
-        guids = dmg_result[['guid']]
+        guids = dmg_result[["guid"]]
         max_val = dmg_result[dmg_states].max(axis=1)
         max_key = dmg_result[dmg_states].idxmax(axis=1)
         dmg_concat = pd.concat([guids, max_val, max_key], axis=1)
-        dmg_concat.rename(columns={0: 'max_prob', 1: 'max_state'}, inplace=True)
+        dmg_concat.rename(columns={0: "max_prob", 1: "max_state"}, inplace=True)
 
         return dmg_concat
diff --git a/pyincore/utils/popdisloutputprocess.py b/pyincore/utils/popdisloutputprocess.py
index 6866667da..262012d0a 100644
--- a/pyincore/utils/popdisloutputprocess.py
+++ b/pyincore/utils/popdisloutputprocess.py
@@ -514,7 +514,19 @@ def pd_total(self, filename_json=None):
         # Dislocated by race and ethnicity
         hud = self.pop_disl_result
         hud_vals = hud["dislocated"].value_counts()
-        hua_disl = [int(hud_vals[False]), int(hud_vals[True])]
+
+        hud_vals_false = 0
+        hud_vals_true = 0
+
+        # if there is "False" category, assign the value
+        if False in hud_vals.keys().values:
+            hud_vals_false = int(hud_vals[False])
+
+        # if there is "True" category, assign the value
+        if True in hud_vals.keys().values:
+            hud_vals_true = int(hud_vals[True])
+
+        hua_disl = [hud_vals_false, hud_vals_true]
 
         pd_disl = [int(hud["numprec"].where(hud["dislocated"] == 0).sum()),
                    int(hud["numprec"].where(hud["dislocated"] == 1).sum())]
diff --git a/recipes/meta.yaml b/recipes/meta.yaml
index 1c8d673f8..5a0057224 100644
--- a/recipes/meta.yaml
+++ b/recipes/meta.yaml
@@ -45,7 +45,7 @@ requirements:
     - geopandas>=0.6.1
     - matplotlib>=2.1.0
     - networkx>=2.2
-    - pandas>=0.24.1
+    - pandas>=0.24.1,<2
     - pyomo>=5.6
     - pyproj>=1.9.6
     - rasterio>=1.0.18
diff --git a/requirements.min b/requirements.min
index eebcf9676..e53900f05 100644
--- a/requirements.min
+++ b/requirements.min
@@ -4,7 +4,7 @@ ipopt>=3.11
 matplotlib>=2.1.0
 networkx>=2.2
 numpy>=1.16.6,<2.0a0
-pandas>=0.24.1
+pandas>=0.24.1,<2
 pycodestyle>=2.6.0
 pyomo>=5.6
 pyproj>=1.9.6
diff --git a/requirements.txt b/requirements.txt
index 378132a66..853930a8b 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -4,7 +4,7 @@ geopandas>=0.6.1
 matplotlib>=2.1.0
 networkx>=2.2
 numpy>=1.16.6,<2.0a0
-pandas>=0.24.1
+pandas>=0.24.1,<2
 pycodestyle>=2.6.0
 pyomo>=5.6
 pyproj>=1.9.6
diff --git a/setup.py b/setup.py
index 3148d3321..a73dfb521 100644
--- a/setup.py
+++ b/setup.py
@@ -7,7 +7,7 @@
 from setuptools import setup, find_packages
 
 # version number of pyincore
-version = '1.9.0'
+version = '1.10.0'
 
 with open("README.rst", encoding="utf-8") as f:
     readme = f.read()
@@ -60,7 +60,7 @@
         'matplotlib>=2.1.0',
         'networkx>=2.2',
         'numpy>=1.16.6,<2.0a0',
-        'pandas>=0.24.1',
+        'pandas>=0.24.1,<2',
         'pyomo>=5.6',
         'pyproj>=1.9.6',
         'rasterio>=1.0.18',
diff --git a/versions.md b/versions.md
index 40474d587..31e18525a 100644
--- a/versions.md
+++ b/versions.md
@@ -4,47 +4,49 @@ This file is autogenerated and lists the latest versions of all dependencies. Th
 
 | PIP Package | Conda Package | Required Version | Matched Version (pip) |
 |-------------|---------------|------------------|-----------------------|
-| fiona | | \>=1.8.4| 1.8.21 |
-| geopandas | | \>=0.6.1| 0.11.0 |
-| matplotlib | | \>=2.1.0| 3.5.2 |
-| networkx | | \>=2.2| 2.8.5 |
-| numpy | | \>=1.16.6,\<2.0a0| 1.23.1 |
-| pandas | | \>=0.24.1| 1.4.3 |
-| pycodestyle | | \>=2.6.0| 2.8.0 |
-| pyomo | | \>=5.6| 6.4.1 |
-| pyproj | | \>=1.9.6| 3.3.1 |
-| pytest | | \>=3.9.0| 7.1.2 |
+| fiona | | \>=1.8.4| 1.9.2 |
+| geopandas | | \>=0.6.1| 0.12.2 |
+| matplotlib | | \>=2.1.0| 3.7.1 |
+| networkx | | \>=2.2| 3.1 |
+| numpy | | \>=1.16.6,\<2.0a0| 1.24.2 |
+| pandas | | \>=0.24.1,\<2| 1.5.3 |
+| pycodestyle | | \>=2.6.0| 2.10.0 |
+| pyomo | | \>=5.6| 6.5.0 |
+| pyproj | | \>=1.9.6| 3.5.0 |
+| pytest | | \>=3.9.0| 7.2.2 |
 | python-jose | | \>=3.0| 3.3.0 |
-| rasterio | | \>=1.0.18| 1.3.0 |
-| requests | | \>=2.20.0| 2.28.1 |
-| rtree | | \>=0.8.3| 1.0.0 |
-| scipy | | \>=1.2.0| 1.8.1 |
-| shapely | | \>=1.6.4.post1| 1.8.2 |
-| wntr | | \>=0.1.6| 0.4.2 |
-| affine | | | 2.3.1 |
-| attrs | | | 21.4.0 |
-| certifi | | | 2022.6.15 |
-| charset-normalizer | | | 2.1.0 |
+| rasterio | | \>=1.0.18| 1.3.6 |
+| requests | | \>=2.20.0| 2.28.2 |
+| rtree | | \>=0.8.3| 1.0.1 |
+| scipy | | \>=1.2.0| 1.10.1 |
+| shapely | | \>=1.6.4.post1| 2.0.1 |
+| wntr | | \>=0.1.6| 1.0.0 |
+| affine | | | 2.4.0 |
+| attrs | | | 22.2.0 |
+| certifi | | | 2022.12.7 |
+| charset-normalizer | | | 3.1.0 |
 | click | | | 8.1.3 |
 | click-plugins | | | 1.1.1 |
 | cligj | | | 0.7.2 |
+| contourpy | | | 1.0.7 |
 | cycler | | | 0.11.0 |
 | ecdsa | | | 0.18.0 |
-| fonttools | | | 4.34.4 |
-| idna | | | 3.3 |
-| iniconfig | | | 1.1.1 |
+| exceptiongroup | | | 1.1.1 |
+| fonttools | | | 4.39.3 |
+| idna | | | 3.4 |
+| iniconfig | | | 2.0.0 |
 | kiwisolver | | | 1.4.4 |
 | munch | | | 2.5.0 |
-| packaging | | | 21.3 |
-| pillow | | | 9.2.0 |
+| packaging | | | 23.0 |
+| pillow | | | 9.5.0 |
 | pluggy | | | 1.0.0 |
 | ply | | | 3.11 |
 | pyasn1 | | | 0.4.8 |
 | pyparsing | | | 3.0.9 |
 | python-dateutil | | | 2.8.2 |
-| pytz | | | 2022.1 |
-| rsa | | | 4.8 |
+| pytz | | | 2023.3 |
+| rsa | | | 4.9 |
 | six | | | 1.16.0 |
 | snuggs | | | 1.4.7 |
 | tomli | | | 2.0.1 |
-| urllib3 | | | 1.26.10 |
+| urllib3 | | | 1.26.15 |


Test patch:
diff --git a/tests/pyincore/analyses/galvestoncge/test_galvestoncge.py b/tests/pyincore/analyses/galvestoncge/test_galvestoncge.py
new file mode 100644
index 000000000..aecefcf6a
--- /dev/null
+++ b/tests/pyincore/analyses/galvestoncge/test_galvestoncge.py
@@ -0,0 +1,26 @@
+from pyincore import IncoreClient, Dataset
+from pyincore.analyses.galvestoncge import GalvestonCGEModel
+import pyincore.globals as pyglobals
+import time
+
+# This script runs GalvestonCGEModel analysis with input files from
+# IN-CORE development services
+
+def run_base_analysis():
+    client = IncoreClient(pyglobals.INCORE_API_DEV_URL)
+    galveston_cge = GalvestonCGEModel(client)
+
+    galveston_cge.set_parameter("model_iterations", 1)
+    galveston_cge.load_remote_input_dataset("SAM", "6420c377b18d026e7c7dc327")
+    galveston_cge.load_remote_input_dataset("BB", "6420c3d2b18d026e7c7dc328")
+    galveston_cge.load_remote_input_dataset("MISCH", "6420c434b18d026e7c7dc32e")
+    galveston_cge.load_remote_input_dataset("EMPLOY", "6420c474b18d026e7c7dc334")
+    galveston_cge.load_remote_input_dataset("JOBCR", "6420c4d9b18d026e7c7dc33a")
+    galveston_cge.load_remote_input_dataset("OUTCR", "6420c511b18d026e7c7dc340")
+    galveston_cge.load_remote_input_dataset("sector_shocks", "64219be5b18d026e7c7e1534")
+
+    galveston_cge.run_analysis()
+
+
+if __name__ == '__main__':
+    run_base_analysis()
diff --git a/tests/pyincore/utils/test_cgecsvoutputjson.py b/tests/pyincore/utils/test_cgecsvoutputjson.py
index 24d681593..9916b17e6 100644
--- a/tests/pyincore/utils/test_cgecsvoutputjson.py
+++ b/tests/pyincore/utils/test_cgecsvoutputjson.py
@@ -1,7 +1,10 @@
 # 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 pytest
 
+from pyincore.analyses.galvestoncge.galvestoncge import GalvestonCGEModel
+from pyincore.analyses.saltlakecge.saltlakecge import SaltLakeCGEModel
 from pyincore import IncoreClient
 from pyincore.analyses.joplincge import JoplinCGEModel
 from pyincore.utils.cgeoutputprocess import CGEOutputProcess
@@ -12,129 +15,176 @@
 # This script runs JoplinCGEModel analysis with input files from
 # IN-CORE development services. The output csv files are converted to json
 # format suitable for the IN-CORE Playbook tool.
+def run_convert_SLC_cge_json_path(testpath):
+    # test the external file with a path
 
-def run_convert_cge_json_chained():
-    client = IncoreClient(pyglobals.INCORE_API_DEV_URL)
-    joplin_cge = JoplinCGEModel(client)
-
-    # SAM
-    sam = "5cdc7b585648c4048fb53062"
-
-    # CAPITAL COMP
-    bb = "5cdc7d295648c4048fb53089"
-
-    # MISC TABLES
-    iout = "5cdc7e965648c4048fb530ef"
-    misc = "5cdc7f4f5648c4048fb53150"
-    misch = "5cdc7fa05648c4048fb53172"
-    landcap = "5cdc7f0a5648c4048fb5312e"
-    employ = "5cdc7df65648c4048fb530ab"
-    igtd = "5cdc7e405648c4048fb530cd"
-    tauff = "5cdc81da5648c4048fb532b7"
-    jobcr = "5cdc7ed25648c4048fb5310c"
-    outcr = "5cdc7fde5648c4048fb53194"
-    sector_shocks = "5f20653e7887544479c6b94a"
-
-    joplin_cge.set_parameter("model_iterations", 1)
-
-    joplin_cge.load_remote_input_dataset("SAM", sam)
-    joplin_cge.load_remote_input_dataset("BB", bb)
-    joplin_cge.load_remote_input_dataset("IOUT", iout)
-    joplin_cge.load_remote_input_dataset("MISC", misc)
-    joplin_cge.load_remote_input_dataset("MISCH", misch)
-    joplin_cge.load_remote_input_dataset("LANDCAP", landcap)
-    joplin_cge.load_remote_input_dataset("EMPLOY", employ)
-    joplin_cge.load_remote_input_dataset("IGTD", igtd)
-    joplin_cge.load_remote_input_dataset("TAUFF", tauff)
-    joplin_cge.load_remote_input_dataset("JOBCR", jobcr)
-    joplin_cge.load_remote_input_dataset("OUTCR", outcr)
-    joplin_cge.load_remote_input_dataset("sector_shocks", sector_shocks)
+    cge_json = CGEOutputProcess()
+    region = ["R1", "R2"]
 
-    joplin_cge.run_analysis()
+    categories = []
+    for h in ["HH1", "HH2", "HH3", "HH4"]:
+        for r in region:
+            categories.append(h + "_" + r)
 
-    domestic_supply_result = joplin_cge.get_output_dataset("domestic-supply")
-    # household_count_result = joplin_cge.get_output_dataset("household-count")
-    gross_income_result = joplin_cge.get_output_dataset("gross-income")
-    pre_demand_result = joplin_cge.get_output_dataset("pre-disaster-factor-demand")
-    post_demand_result = joplin_cge.get_output_dataset("post-disaster-factor-demand")
+    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)
 
-    cge_json = CGEOutputProcess()
-    cge_json.get_cge_domestic_supply(domestic_supply_result, None, "cge_domestic_supply.json")
-    # cge_json.get_cge_household_count(household_count_result, None, "cge_total_household_count.json")
-    cge_json.get_cge_gross_income(gross_income_result, None, "cge_total_household_income.json")
-    cge_json.get_cge_employment(pre_demand_result, post_demand_result, None, None, "cge_employment.json")
+    categories = []
+    for d in ["AG_MI", "UTIL", "CONS", "MANU", "COMMER", "EDU", "HEALTH", "ART_ACC", "RELIG"]:
+        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", "HS3"]:
+        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
 
 
-def run_convert_cge_json_path(testpath):
+def run_convert_Joplin_cge_json_path(testpath):
     # test the external file with a path
 
     cge_json = CGEOutputProcess()
+    categories = ["HH1", "HH2", "HH3", "HH4", "HH5"]
     cge_json.get_cge_household_count(None,
                                      os.path.join(testpath, "household-count.csv"),
-                                     "cge_total_household_count.json")
+                                     "joplin_cge_total_household_count.json", income_categories=categories)
     cge_json.get_cge_gross_income(None,
                                   os.path.join(testpath, "gross-income.csv"),
-                                  "cge_total_household_income.json")
-    cge_json.get_cge_employment(None, None,
-                                os.path.join(testpath, "pre-disaster-factor-demand.csv"),
+                                  "joplin_cge_total_household_income.json", income_categories=categories)
+
+    categories = ["GOODS", "TRADE", "OTHER"]
+    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"),
-                                "cge_employment.json")
+                                "joplin_cge_employment.json", demand_categories=categories)
+
+    categories = ["Goods", "Trades", "Others", "HS1", "HS2", "HS3"]
     cge_json.get_cge_domestic_supply(None,
-                                     os.path.join(testpath, "domestic-supply.csv"),
-                                     "cge_domestic_supply.json")
+                                     os.path.join(testpath, "domestic-supply.csv"), "joplin_cge_domestic_supply.json",
+                                     supply_categories=categories)
     return True
 
 
-def run_convert_SLC_cge_json_path(testpath):
+def run_convert_Galveston_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)
+    region = ["I", "M"]
 
+    categories = ["IHH1", "IHH2", "IHH3", "IHH4", "IHH5", "MHH1", "MHH2", "MHH3", "MHH4", "MHH5"]
     cge_json.get_cge_household_count(None,
                                      os.path.join(testpath, "household-count.csv"),
-                                     "slc_cge_total_household_count.json", income_categories=categories)
+                                     "galveston_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)
+                                  "galveston_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 d in ["AGMIN", "UTIL", "CONS", "MANU", "WHOLE", "RETAIL", "TRANS", "PROFSER", "REALE", "EDU", "HEALTH", "ART",
+              "ACCO"]:
         for r in region:
-            categories.append(d + "_" + r)
+            categories.append(r + d)
     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)
+                                "galveston_cge_employment.json", demand_categories=categories)
 
     categories = []
-    for d in ["AG_MI", "UTIL", "CONS", "MANU", "COMMER", "EDU", "HEALTH", "ART_ACC", "RELIG",
-              "HS1", "HS2"]:
+    for d in ["AGMIN", "UTIL", "CONS", "MANU", "WHOLE", "RETAIL", "TRANS", "PROFSER", "REALE", "EDU", "HEALTH", "ART",
+              "ACCO"]:
         for r in region:
-            categories.append(d + "_" + r)
+            categories.append(r + d)  # e.g.iAGMIN
+    for d in ["HS1", "HS2", "HS3"]:
+        for r in region:
+            categories.append(d + r)  # e.g.HS1I
     cge_json.get_cge_domestic_supply(None,
-                                     os.path.join(testpath, "domestic-supply.csv"), "slc_cge_domestic_supply.json",
+                                     os.path.join(testpath, "domestic-supply.csv"),
+                                     "galveston_cge_domestic_supply.json",
                                      supply_categories=categories)
     return True
 
 
 if __name__ == '__main__':
-    # test chaining with Joplin CGE analysis
-    run_convert_cge_json_chained()
 
-    # test the external file with a path
-    testpath = ""
-    # testpath = "/Users/<user>/<path_to_pyincore>/pyincore/tests/pyincore/utils"
-    if testpath:
-        run_convert_cge_json_path(testpath)
+    # # run slc cge
+    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")
+    saltlake_cge.load_remote_input_dataset("sector_shocks", "64075ec46121f9438872a802")
+
+    saltlake_cge.run_analysis()
+    run_convert_SLC_cge_json_path(testpath="./")
+    print("Salt lake city post processing done.")
+
+    # run joplin cge
+    joplin_cge = JoplinCGEModel(client)
+
+    # SAM
+    sam = "5cdc7b585648c4048fb53062"
+
+    # CAPITAL COMP
+    bb = "5cdc7d295648c4048fb53089"
+
+    # MISC TABLES
+    iout = "5cdc7e965648c4048fb530ef"
+    misc = "5cdc7f4f5648c4048fb53150"
+    misch = "5cdc7fa05648c4048fb53172"
+    landcap = "5cdc7f0a5648c4048fb5312e"
+    employ = "5cdc7df65648c4048fb530ab"
+    igtd = "5cdc7e405648c4048fb530cd"
+    tauff = "5cdc81da5648c4048fb532b7"
+    jobcr = "5cdc7ed25648c4048fb5310c"
+    outcr = "5cdc7fde5648c4048fb53194"
+    sector_shocks = "5f20653e7887544479c6b94a"
+
+    joplin_cge.set_parameter("model_iterations", 1)
+
+    joplin_cge.load_remote_input_dataset("SAM", sam)
+    joplin_cge.load_remote_input_dataset("BB", bb)
+    joplin_cge.load_remote_input_dataset("IOUT", iout)
+    joplin_cge.load_remote_input_dataset("MISC", misc)
+    joplin_cge.load_remote_input_dataset("MISCH", misch)
+    joplin_cge.load_remote_input_dataset("LANDCAP", landcap)
+    joplin_cge.load_remote_input_dataset("EMPLOY", employ)
+    joplin_cge.load_remote_input_dataset("IGTD", igtd)
+    joplin_cge.load_remote_input_dataset("TAUFF", tauff)
+    joplin_cge.load_remote_input_dataset("JOBCR", jobcr)
+    joplin_cge.load_remote_input_dataset("OUTCR", outcr)
+    joplin_cge.load_remote_input_dataset("sector_shocks", sector_shocks)
+
+    joplin_cge.run_analysis()
+    run_convert_Joplin_cge_json_path("./")
+    print("Joplin post processing done.")
+
+    galveston_cge = GalvestonCGEModel(client)
+
+    galveston_cge.set_parameter("model_iterations", 1)
+    galveston_cge.load_remote_input_dataset("SAM", "6420c377b18d026e7c7dc327")
+    galveston_cge.load_remote_input_dataset("BB", "6420c3d2b18d026e7c7dc328")
+    galveston_cge.load_remote_input_dataset("MISCH", "6420c434b18d026e7c7dc32e")
+    galveston_cge.load_remote_input_dataset("EMPLOY", "6420c474b18d026e7c7dc334")
+    galveston_cge.load_remote_input_dataset("JOBCR", "6420c4d9b18d026e7c7dc33a")
+    galveston_cge.load_remote_input_dataset("OUTCR", "6420c511b18d026e7c7dc340")
+    galveston_cge.load_remote_input_dataset("sector_shocks", "64219be5b18d026e7c7e1534")
 
-    run_convert_SLC_cge_json_path(testpath="../../data/SLC_CGE")
+    galveston_cge.run_analysis()
 
-    print("DONE")
+    run_convert_Galveston_cge_json_path("./")
+    print("Galveston post processing done.")
diff --git a/tests/pyincore/utils/test_dataprocessutil.py b/tests/pyincore/utils/test_dataprocessutil.py
index d2adefe66..db0d20d79 100644
--- a/tests/pyincore/utils/test_dataprocessutil.py
+++ b/tests/pyincore/utils/test_dataprocessutil.py
@@ -1,7 +1,10 @@
+import json
+
+import pandas as pd
 import pytest
 
-from pyincore.utils.dataprocessutil import DataProcessUtil as util
 from pyincore import Dataset, DataService
+from pyincore.utils.dataprocessutil import DataProcessUtil as util
 
 
 @pytest.fixture
@@ -12,11 +15,11 @@ def client():
 def test_get_mapped_result(client):
     bldg_dataset_id = "5f9091df3e86721ed82f701d"
     bldg_dmg_dataset_id = "5f9868c00ace240b22a7f2a5"
-    bldg_func_dataset_id = "60674c89c57ada48e492b97b"
+    bldg_func_state_dataset_id = "642f66e5f27db6680103c4ad"
     archetype_id = "5fca915fb34b193f7a44059b"
 
     dmg_ret_json, func_ret_json, max_state_df = util.get_mapped_result_from_dataset_id(
-        client, bldg_dataset_id, bldg_dmg_dataset_id, bldg_func_dataset_id, archetype_id)
+        client, bldg_dataset_id, bldg_dmg_dataset_id, bldg_func_state_dataset_id, archetype_id)
 
     assert "by_cluster" in dmg_ret_json.keys() and "by_category" in dmg_ret_json.keys()
     assert "by_cluster" in func_ret_json.keys() and "by_category" in func_ret_json.keys()
@@ -32,12 +35,72 @@ def test_get_mapped_result_from_analysis(client):
 
     archetype_id = "5fca915fb34b193f7a44059b"
 
-    bldg_func_dataset_id = "60674c89c57ada48e492b97b"
-    bldg_func_dataset = Dataset.from_data_service(bldg_func_dataset_id, DataService(client))
+    bldg_func_state_dataset_id = "642f66e5f27db6680103c4ad"
+    bldg_func_state_dataset = Dataset.from_data_service(bldg_func_state_dataset_id, DataService(client))
 
     dmg_ret_json, func_ret_json, max_state_df = util.get_mapped_result_from_analysis(
-        client, bldg_dataset_id, dmg_result_dataset, bldg_func_dataset, archetype_id)
+        client, bldg_dataset_id, dmg_result_dataset, bldg_func_state_dataset, archetype_id)
 
     assert "by_cluster" in dmg_ret_json.keys() and "by_category" in dmg_ret_json.keys()
     assert "by_cluster" in func_ret_json.keys() and "by_category" in func_ret_json.keys()
     assert "max_state" in max_state_df._info_axis.values
+
+
+def _functionality_cluster(client, archetype_mapping="5fca915fb34b193f7a44059b",
+                           building_dataset_id="5fa0b132cc6848728b66948d",
+                           bldg_func_state_id="5f0f6fbfb922f96f4e989ed8",
+                           arch_column="archetype",
+                           title="joplin_mcs"):
+    dataservice = DataService(client)
+
+    # Archetype mapping file
+    archetype_mapping_dataset = Dataset.from_data_service(archetype_mapping, dataservice)
+    archetype_mapping_path = archetype_mapping_dataset.get_file_path()
+    arch_mapping = pd.read_csv(archetype_mapping_path)
+
+    # Building dataset id
+    bldg_dataset = Dataset.from_data_service(building_dataset_id, dataservice)
+    buildings = bldg_dataset.get_dataframe_from_shapefile()
+
+    # Cluster the mcs building failure probability - essentially building functionality without electric power being
+    # considered
+    bldg_func_state_dataset = Dataset.from_data_service(bldg_func_state_id, dataservice)
+    bldg_func_state_dataset_path = bldg_func_state_dataset.get_file_path()
+    bldg_func_state = pd.read_csv(bldg_func_state_dataset_path)
+
+    ret_json = util.create_mapped_func_result(buildings, bldg_func_state, arch_mapping, arch_column)
+
+    with open(title + "_cluster.json", "w") as f:
+        json.dump(ret_json, f, indent=2)
+
+
+def test_joplin_mcs_cluster(client):
+    _functionality_cluster(client, archetype_mapping="5fca915fb34b193f7a44059b",
+                           building_dataset_id="5fa0b132cc6848728b66948d",
+                           bldg_func_state_id="5f0f6fbfb922f96f4e989ed8",
+                           arch_column="archetype",
+                           title="joplin_mcs")
+
+
+def test_joplin_bldg_func_cluster(client):
+    _functionality_cluster(client, archetype_mapping="5fca915fb34b193f7a44059b",
+                           building_dataset_id="5fa0b132cc6848728b66948d",
+                           bldg_func_state_id="642f66e5f27db6680103c4ad",
+                           arch_column="archetype",
+                           title="joplin_bldg_func")
+
+
+def test_galveston_mcs_cluster(client):
+    _functionality_cluster(client, archetype_mapping="6420befbb18d026e7c7dbafc",
+                           building_dataset_id="63ff69a96d3b2a308baaca12",
+                           bldg_func_state_id="642f6230f27db66801038964",
+                           arch_column="arch_wind",
+                           title="galveston_mcs")
+
+
+def test_galveston_bldg_func_cluster(client):
+    _functionality_cluster(client, archetype_mapping="6420befbb18d026e7c7dbafc",
+                           building_dataset_id="63ff69a96d3b2a308baaca12",
+                           bldg_func_state_id="642f6495f27db66801039434",
+                           arch_column="arch_wind",
+                           title="galveston_bldg_func")