ConsIndShock: cFunc extends below mNrmMin
A student at Goethe University flagged this issue for us.

The consumption function in ConsIndShockModel exists below mNrmMin, but it should not.

Explanation: cFunc is a LowerEnvelope instance containing cFuncCnst and cFuncUnc.  cFuncCnst hits zero and ceases to exist at mNrmMin, but cFuncUnc continues further down to BoroCnstNat.  LowerEnvelope uses nanmax rather than max, so it throws out the nan from cFuncCnst and thinks it should use cFuncUnc below mNrmMin.

Solution: Add optional boolean input for LowerEnvelope (and related classes) to indicate whether to use max/min or nanmax/nanmin.  Default should be nanmax.
diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py
index 4ffd9a6d..76923542 100644
--- a/HARK/ConsumptionSaving/ConsIndShockModel.py
+++ b/HARK/ConsumptionSaving/ConsIndShockModel.py
@@ -893,7 +893,8 @@ class ConsIndShockSolverBasic(ConsIndShockSetup):
         cFuncNowUnc = interpolator(mNrm,cNrm)
 
         # Combine the constrained and unconstrained functions into the true consumption function
-        cFuncNow = LowerEnvelope(cFuncNowUnc,self.cFuncNowCnst)
+        cFuncNow = LowerEnvelope(cFuncNowUnc, self.cFuncNowCnst,
+                                 nan_bool = False)
 
         # Make the marginal value function and the marginal marginal value function
         vPfuncNow = MargValueFunc(cFuncNow,self.CRRA)
@@ -1214,7 +1215,7 @@ def solveConsIndShock(solution_next,IncomeDstn,LivPrb,DiscFac,CRRA,Rfree,PermGro
         included in the reported solution.
     CubicBool: boolean
         Indicator for whether the solver should use cubic or linear interpolation.
-
+    
     Returns
     -------
     solution_now : ConsumerSolution
diff --git a/HARK/interpolation.py b/HARK/interpolation.py
index ff138a81..45b69f6f 100644
--- a/HARK/interpolation.py
+++ b/HARK/interpolation.py
@@ -1642,7 +1642,7 @@ class LowerEnvelope(HARKinterpolator1D):
     '''
     distance_criteria = ['functions']
 
-    def __init__(self,*functions):
+    def __init__(self, *functions, nan_bool = True):
         '''
         Constructor to make a new lower envelope iterpolation.
 
@@ -1650,11 +1650,21 @@ class LowerEnvelope(HARKinterpolator1D):
         ----------
         *functions : function
             Any number of real functions; often instances of HARKinterpolator1D
-
+        nan_bool : boolean
+            An indicator for whether the solver should exclude NA's when forming
+            the lower envelope.
         Returns
         -------
         new instance of LowerEnvelope
         '''
+        
+        if nan_bool:
+            self.compare = np.nanmin
+            self.argcompare = np.nanargmin
+        else:
+            self.compare = np.min
+            self.argcompare = np.argmin
+        
         self.functions = []
         for function in functions:
             self.functions.append(function)
@@ -1665,14 +1675,16 @@ class LowerEnvelope(HARKinterpolator1D):
         Returns the level of the function at each value in x as the minimum among
         all of the functions.  Only called internally by HARKinterpolator1D.__call__.
         '''
+        
         if _isscalar(x):
-            y = np.nanmin([f(x) for f in self.functions])
+            y = self.compare([f(x) for f in self.functions])
         else:
             m = len(x)
             fx = np.zeros((m,self.funcCount))
             for j in range(self.funcCount):
                 fx[:,j] = self.functions[j](x)
-            y = np.nanmin(fx,axis=1)
+            y = self.compare(fx,axis=1)
+
         return y
 
     def _der(self,x):
@@ -1680,7 +1692,7 @@ class LowerEnvelope(HARKinterpolator1D):
         Returns the first derivative of the function at each value in x.  Only
         called internally by HARKinterpolator1D.derivative.
         '''
-        y,dydx = self.eval_with_derivative(x)
+        y,dydx = self._evalAndDer(x)
         return dydx  # Sadly, this is the fastest / most convenient way...
 
     def _evalAndDer(self,x):
@@ -1692,8 +1704,7 @@ class LowerEnvelope(HARKinterpolator1D):
         fx = np.zeros((m,self.funcCount))
         for j in range(self.funcCount):
             fx[:,j] = self.functions[j](x)
-        fx[np.isnan(fx)] = np.inf
-        i = np.argmin(fx,axis=1)
+        i = self.argcompare(fx,axis=1)
         y = fx[np.arange(m),i]
         dydx = np.zeros_like(y)
         for j in range(self.funcCount):
@@ -1701,7 +1712,6 @@ class LowerEnvelope(HARKinterpolator1D):
             dydx[c] = self.functions[j].derivative(x[c])
         return y,dydx
 
-
 class UpperEnvelope(HARKinterpolator1D):
     '''
     The upper envelope of a finite set of 1D functions, each of which can be of
@@ -1710,7 +1720,7 @@ class UpperEnvelope(HARKinterpolator1D):
     '''
     distance_criteria = ['functions']
 
-    def __init__(self,*functions):
+    def __init__(self,*functions, nan_bool=True):
         '''
         Constructor to make a new upper envelope iterpolation.
 
@@ -1718,11 +1728,21 @@ class UpperEnvelope(HARKinterpolator1D):
         ----------
         *functions : function
             Any number of real functions; often instances of HARKinterpolator1D
+        nan_bool : boolean
+            An indicator for whether the solver should exclude NA's when forming
+            the lower envelope.
 
         Returns
         -------
         new instance of UpperEnvelope
         '''
+        if nan_bool:
+            self.compare = np.nanmax
+            self.argcompare = np.nanargmax
+        else:
+            self.compare = np.max
+            self.argcompare = np.argmax
+        
         self.functions = []
         for function in functions:
             self.functions.append(function)
@@ -1734,13 +1754,14 @@ class UpperEnvelope(HARKinterpolator1D):
         all of the functions.  Only called internally by HARKinterpolator1D.__call__.
         '''
         if _isscalar(x):
-            y = np.nanmax([f(x) for f in self.functions])
+            y = self.compare([f(x) for f in self.functions])
         else:
             m = len(x)
             fx = np.zeros((m,self.funcCount))
             for j in range(self.funcCount):
                 fx[:,j] = self.functions[j](x)
-            y = np.nanmax(fx,axis=1)
+            y = self.compare(fx,axis=1)
+        
         return y
 
     def _der(self,x):
@@ -1748,7 +1769,7 @@ class UpperEnvelope(HARKinterpolator1D):
         Returns the first derivative of the function at each value in x.  Only
         called internally by HARKinterpolator1D.derivative.
         '''
-        y,dydx = self.eval_with_derivative(x)
+        y,dydx = self._evalAndDer(x)
         return dydx  # Sadly, this is the fastest / most convenient way...
 
     def _evalAndDer(self,x):
@@ -1760,8 +1781,7 @@ class UpperEnvelope(HARKinterpolator1D):
         fx = np.zeros((m,self.funcCount))
         for j in range(self.funcCount):
             fx[:,j] = self.functions[j](x)
-        fx[np.isnan(fx)] = np.inf
-        i = np.argmax(fx,axis=1)
+        i = self.argcompare(fx,axis=1)
         y = fx[np.arange(m),i]
         dydx = np.zeros_like(y)
         for j in range(self.funcCount):
@@ -1778,7 +1798,7 @@ class LowerEnvelope2D(HARKinterpolator2D):
     '''
     distance_criteria = ['functions']
 
-    def __init__(self,*functions):
+    def __init__(self,*functions, nan_bool = True):
         '''
         Constructor to make a new lower envelope iterpolation.
 
@@ -1786,11 +1806,21 @@ class LowerEnvelope2D(HARKinterpolator2D):
         ----------
         *functions : function
             Any number of real functions; often instances of HARKinterpolator2D
+         nan_bool : boolean
+            An indicator for whether the solver should exclude NA's when forming
+            the lower envelope.
 
         Returns
         -------
         new instance of LowerEnvelope2D
         '''
+        
+        if nan_bool:
+            self.compare = np.nanmin
+            self.argcompare = np.nanargmin
+        else:
+            self.compare = np.min
+            self.argcompare = np.argmin
         self.functions = []
         for function in functions:
             self.functions.append(function)
@@ -1802,14 +1832,16 @@ class LowerEnvelope2D(HARKinterpolator2D):
         among all of the functions.  Only called internally by
         HARKinterpolator2D.__call__.
         '''
+        
         if _isscalar(x):
-            f = np.nanmin([f(x,y) for f in self.functions])
+            f = self.compare([f(x,y) for f in self.functions])
         else:
             m = len(x)
             temp = np.zeros((m,self.funcCount))
             for j in range(self.funcCount):
                 temp[:,j] = self.functions[j](x,y)
-            f = np.nanmin(temp,axis=1)
+            f = self.compare(temp,axis=1)
+        
         return f
 
     def _derX(self,x,y):
@@ -1821,8 +1853,7 @@ class LowerEnvelope2D(HARKinterpolator2D):
         temp = np.zeros((m,self.funcCount))
         for j in range(self.funcCount):
             temp[:,j] = self.functions[j](x,y)
-        temp[np.isnan(temp)] = np.inf
-        i = np.argmin(temp,axis=1)
+        i = self.argcompare(temp,axis=1)
         dfdx = np.zeros_like(x)
         for j in range(self.funcCount):
             c = i == j
@@ -1838,8 +1869,7 @@ class LowerEnvelope2D(HARKinterpolator2D):
         temp = np.zeros((m,self.funcCount))
         for j in range(self.funcCount):
             temp[:,j] = self.functions[j](x,y)
-        temp[np.isnan(temp)] = np.inf
-        i = np.argmin(temp,axis=1)
+        i = self.argcompare(temp,axis=1)
         y = temp[np.arange(m),i]
         dfdy = np.zeros_like(x)
         for j in range(self.funcCount):
@@ -1856,7 +1886,7 @@ class LowerEnvelope3D(HARKinterpolator3D):
     '''
     distance_criteria = ['functions']
 
-    def __init__(self,*functions):
+    def __init__(self,*functions, nan_bool = True):
         '''
         Constructor to make a new lower envelope iterpolation.
 
@@ -1864,11 +1894,20 @@ class LowerEnvelope3D(HARKinterpolator3D):
         ----------
         *functions : function
             Any number of real functions; often instances of HARKinterpolator3D
+        nan_bool : boolean
+            An indicator for whether the solver should exclude NA's when forming
+            the lower envelope.
 
         Returns
         -------
         None
         '''
+        if nan_bool:
+            self.compare = np.nanmin
+            self.argcompare = np.nanargmin
+        else:
+            self.compare = np.min
+            self.argcompare = np.argmin
         self.functions = []
         for function in functions:
             self.functions.append(function)
@@ -1880,14 +1919,16 @@ class LowerEnvelope3D(HARKinterpolator3D):
         among all of the functions.  Only called internally by
         HARKinterpolator3D.__call__.
         '''
+        
         if _isscalar(x):
-            f = np.nanmin([f(x,y,z) for f in self.functions])
+            f = self.compare([f(x,y,z) for f in self.functions])
         else:
             m = len(x)
             temp = np.zeros((m,self.funcCount))
             for j in range(self.funcCount):
                 temp[:,j] = self.functions[j](x,y,z)
-            f = np.nanmin(temp,axis=1)
+            f = self.compare(temp,axis=1)
+            
         return f
 
     def _derX(self,x,y,z):
@@ -1899,8 +1940,7 @@ class LowerEnvelope3D(HARKinterpolator3D):
         temp = np.zeros((m,self.funcCount))
         for j in range(self.funcCount):
             temp[:,j] = self.functions[j](x,y,z)
-        temp[np.isnan(temp)] = np.inf
-        i = np.argmin(temp,axis=1)
+        i = self.argcompare(temp,axis=1)
         dfdx = np.zeros_like(x)
         for j in range(self.funcCount):
             c = i == j
@@ -1916,8 +1956,7 @@ class LowerEnvelope3D(HARKinterpolator3D):
         temp = np.zeros((m,self.funcCount))
         for j in range(self.funcCount):
             temp[:,j] = self.functions[j](x,y,z)
-        temp[np.isnan(temp)] = np.inf
-        i = np.argmin(temp,axis=1)
+        i = self.argcompare(temp,axis=1)
         y = temp[np.arange(m),i]
         dfdy = np.zeros_like(x)
         for j in range(self.funcCount):
@@ -1934,8 +1973,7 @@ class LowerEnvelope3D(HARKinterpolator3D):
         temp = np.zeros((m,self.funcCount))
         for j in range(self.funcCount):
             temp[:,j] = self.functions[j](x,y,z)
-        temp[np.isnan(temp)] = np.inf
-        i = np.argmin(temp,axis=1)
+        i = self.argcompare(temp,axis=1)
         y = temp[np.arange(m),i]
         dfdz = np.zeros_like(x)
         for j in range(self.funcCount):
['HARK/ConsumptionSaving/ConsIndShockModel.py', 'HARK/interpolation.py']