Repository: gbrammer/eazy-py

Issue:
Asada+ CGM implementation will always be used to generate the tempfilt file
In photoz.py, when the TemplateGrid is [initialized on line 359](https://github.com/gbrammer/eazy-py/blob/5f96fb6176729032838f456ec57af1215c34dbe0/eazy/photoz.py#L359):

```
            self.tempfilt = TemplateGrid(self.zgrid, self.templates, 
                                        RES=self.param['FILTERS_RES'], 
                                        f_numbers=self.f_numbers, 
                                        add_igm=self.param['IGM_SCALE_TAU'], 
                                    galactic_ebv=self.MW_EBV, 
                                    Eb=self.param['SCALE_2175_BUMP'], 
                                    n_proc=n_proc, cosmology=self.cosmology, 
                                    array_dtype=self.ARRAY_DTYPE, 
                                    tempfilt_data=tempfilt_data)
```

The `add_cgm` flag is not set, and since [it's set to true by default](https://github.com/gbrammer/eazy-py/blob/5f96fb6176729032838f456ec57af1215c34dbe0/eazy/photoz.py#L5655), the Asada+ CGM prescription will be run regardless of whether the user sets ADD_CGM to be true or false. I have updated this by adding an extra line:

```
            self.tempfilt = TemplateGrid(self.zgrid, self.templates, 
                                        RES=self.param['FILTERS_RES'], 
                                        f_numbers=self.f_numbers, 
                                        add_igm=self.param['IGM_SCALE_TAU'],
                                        add_cgm=self.param['ADD_CGM'],
                                    galactic_ebv=self.MW_EBV, 
                                    Eb=self.param['SCALE_2175_BUMP'], 
                                    n_proc=n_proc, cosmology=self.cosmology, 
                                    array_dtype=self.ARRAY_DTYPE, 
                                    tempfilt_data=tempfilt_data)

```

In line 359+, but [also here as well](https://github.com/gbrammer/eazy-py/blob/5f96fb6176729032838f456ec57af1215c34dbe0/eazy/photoz.py#L2531C1-L2540C1).

Hints from issue comments:

Solution patch:
diff --git a/eazy/igm.py b/eazy/igm.py
index 9c0d190..d4b3c11 100644
--- a/eazy/igm.py
+++ b/eazy/igm.py
@@ -7,7 +7,7 @@ from . import __file__ as filepath
 __all__ = ["Asada24","Inoue14"]
 
 class Asada24(object):
-    def __init__(self, sigmoid_params=(3.48347968, 1.25809685, 18.24922789), scale_tau=1., add_cgm=True):
+    def __init__(self, sigmoid_params=(3.5918, 1.8414, 18.001), scale_tau=1., add_cgm=True):
         """
         Compute IGM+CGM transmission from Asada et al. 2024, in prep.
         The IGM model is from Inoue+ (2014).
diff --git a/eazy/photoz.py b/eazy/photoz.py
index ab97c33..37dd69e 100644
--- a/eazy/photoz.py
+++ b/eazy/photoz.py
@@ -359,8 +359,10 @@ class PhotoZ(object):
             self.tempfilt = TemplateGrid(self.zgrid, self.templates, 
                                         RES=self.param['FILTERS_RES'], 
                                         f_numbers=self.f_numbers, 
-                                        add_igm=self.param['IGM_SCALE_TAU'], 
-                                    galactic_ebv=self.MW_EBV, 
+                                        add_igm=self.param['IGM_SCALE_TAU'],
+                                    add_cgm=IGM_OBJECT.add_cgm,
+                                    sigmoid_params=IGM_OBJECT.sigmoid_params,
+                                    galactic_ebv=self.MW_EBV,
                                     Eb=self.param['SCALE_2175_BUMP'], 
                                     n_proc=n_proc, cosmology=self.cosmology, 
                                     array_dtype=self.ARRAY_DTYPE, 
@@ -2532,8 +2534,10 @@ class PhotoZ(object):
             self.tempfilt = TemplateGrid(self.zgrid, self.templates, 
                                          RES=self.RES, 
                                          f_numbers=self.f_numbers, 
-                                         add_igm=self.param['IGM_SCALE_TAU'], 
-                                         galactic_ebv=self.MW_EBV, 
+                                         add_igm=self.param['IGM_SCALE_TAU'],
+                                         add_cgm=IGM_OBJECT.add_cgm,
+                                         sigmoid_params=IGM_OBJECT.sigmoid_params,
+                                         galactic_ebv=self.MW_EBV,
                                          Eb=self.param['SCALE_2175_BUMP'], 
                                          n_proc=0, cosmology=self.cosmology, 
                                          array_dtype=self.ARRAY_DTYPE)
@@ -3387,7 +3391,9 @@ class PhotoZ(object):
                                    RES=self.RES, 
                                    f_numbers=np.array(f_numbers), 
                                    add_igm=self.param['IGM_SCALE_TAU'],
-                                   galactic_ebv=self.MW_EBV, 
+                                   add_cgm=IGM_OBJECT.add_cgm,
+                                   sigmoid_params=IGM_OBJECT.sigmoid_params,
+                                   galactic_ebv=self.MW_EBV,
                                    Eb=self.param['SCALE_2175_BUMP'], 
                                    n_proc=n_proc, verbose=verbose, 
                                    cosmology=self.cosmology,
@@ -5094,8 +5100,10 @@ class PhotoZ(object):
             tempfilt = TemplateGrid(self.zgrid, template_list, 
                                     RES=self.RES, 
                                     f_numbers=self.f_numbers, 
-                                    add_igm=self.param['IGM_SCALE_TAU'], 
-                                    galactic_ebv=self.MW_EBV, 
+                                    add_igm=self.param['IGM_SCALE_TAU'],
+                                    add_cgm=IGM_OBJECT.add_cgm,
+                                    sigmoid_params=IGM_OBJECT.sigmoid_params,
+                                    galactic_ebv=self.MW_EBV,
                                     Eb=self.param['SCALE_2175_BUMP'], 
                                     cosmology=self.cosmology, 
                                     array_dtype=self.ARRAY_DTYPE)
@@ -5652,7 +5660,7 @@ def _obj_nnls(coeffs, A, fnu_i, efnu_i):
 
 
 class TemplateGrid(object):
-    def __init__(self, zgrid, templates, RES='FILTERS.RES.latest', f_numbers=[156], add_igm=True, add_cgm=True, sigmoid_params=(3.48347968, 1.25809685, 18.24922789), galactic_ebv=0, Eb=0, n_proc=4, interpolator=None, filters=None, verbose=2, cosmology=None, array_dtype=np.float32, tempfilt_data=None):
+    def __init__(self, zgrid, templates, RES='FILTERS.RES.latest', f_numbers=[156], add_igm=True, add_cgm=True, sigmoid_params=(3.5918, 1.8414, 18.001), galactic_ebv=0, Eb=0, n_proc=4, interpolator=None, filters=None, verbose=2, cosmology=None, array_dtype=np.float32, tempfilt_data=None):
         """
         Integrate filters through filters on a redshift grid
         
@@ -6542,9 +6550,14 @@ def template_lsq(fnu_i, efnu_i, Ain, TEFz, zp, ndraws, fitter, renorm_t, hess_th
         #try:
         if ok_temp.sum() > 0:
             covar = utils.safe_invert(mat)
-            draws = np.random.multivariate_normal((coeffs_i * An)[ok_temp], 
-                                                  covar, 
-                                                  size=ndraws)
+            if np.isfinite(covar).sum() != covar.size:
+                draws = np.full((ndraws, ok_temp.sum()), np.nan)
+            else:
+                draws = np.random.multivariate_normal(
+                    (coeffs_i * An)[ok_temp],
+                    covar,
+                    size=ndraws
+                )
             coeffs_draw[:, ok_temp] = draws / An[ok_temp]
         else:
             #print('Error getting coeffs draws')
diff --git a/eazy/templates.py b/eazy/templates.py
index dc69e27..6b6dd23 100644
--- a/eazy/templates.py
+++ b/eazy/templates.py
@@ -756,7 +756,7 @@ class Template():
                             redshifts=self.redshifts)
 
 
-    def to_observed_frame(self, z=0, scalar=1., extra_sigma=0, lsf_func=None, to_air=False, wavelengths=None, smoothspec_kwargs={'fftsmooth':False}, include_igm=True, sigmoid_params=(3.48347968, 1.25809685, 18.24922789), scale_tau=1., add_cgm=True, clip_wavelengths=[4500,9400], as_template=True):
+    def to_observed_frame(self, z=0, scalar=1., extra_sigma=0, lsf_func=None, to_air=False, wavelengths=None, smoothspec_kwargs={'fftsmooth':False}, include_igm=True, sigmoid_params=(3.5918, 1.8414, 18.001), scale_tau=1., add_cgm=True, clip_wavelengths=[4500,9400], as_template=True):
         """
         Smooth and resample to observed-frame wavelengths, including an
         optional Line Spread Function (LSF)
@@ -1040,7 +1040,7 @@ class Template():
         return iz
 
 
-    def zscale(self, z, scalar=1, include_igm=True, sigmoid_params=(3.48347968, 1.25809685, 18.24922789), scale_tau=1., add_cgm=True, **kwargs):
+    def zscale(self, z, scalar=1, include_igm=True, sigmoid_params=(3.5918, 1.8414, 18.001), scale_tau=1., add_cgm=True, **kwargs):
         """Redshift the template and multiply by a scalar.
 
         Parameters
@@ -1090,7 +1090,7 @@ class Template():
                         name=f'{self.name} z={z}')
 
 
-    def integrate_filter(self, filt, flam=False, scale=1., z=0, include_igm=False, sigmoid_params=(3.48347968, 1.25809685, 18.24922789), scale_tau=1., add_cgm=True, redshift_type='nearest', iz=None):
+    def integrate_filter(self, filt, flam=False, scale=1., z=0, include_igm=False, sigmoid_params=(3.5918, 1.8414, 18.001), scale_tau=1., add_cgm=True, redshift_type='nearest', iz=None):
         """
         Integrate the template through a `FilterDefinition` filter object.
         
@@ -1192,7 +1192,7 @@ class Template():
             return np.array(fluxes)
 
 
-    def igm_absorption(self, z, sigmoid_params=(3.48347968, 1.25809685, 18.24922789), scale_tau=1., add_cgm=True):
+    def igm_absorption(self, z, sigmoid_params=(3.5918, 1.8414, 18.001), scale_tau=1., add_cgm=True):
         """
         Compute IGM absorption with `eazy.igm.Asada24`.
         
diff --git a/eazy/utils.py b/eazy/utils.py
index 2e1d36e..86802b0 100644
--- a/eazy/utils.py
+++ b/eazy/utils.py
@@ -369,14 +369,37 @@ def fill_between_steps(x, y, z, ax=None, *args, **kwargs):
     ax.fill_between(xfull[so], yfull[so], zfull[so], *args, **kwargs)
 
 
-def safe_invert(arr):
+def safe_invert(arr, tol=None):
     """
     Version-safe matrix inversion using `numpy.linalg.inv` or `numpy.matrix.I`
+
+    Parameters
+    ----------
+    arr : array-like
+        Square matrix to invert
+
+    tol : float
+        Tolerance parameter for calculating if arr is singular using the rank from
+        ``np.linalg.matrix_rank(arr, tol=tol)``.  If None, use the numpy default
+        derived from the matrix size and machine epsilon.  If ``arr`` is found to be
+        singular (rank < max(arr.shape)), don't try to invert it and just return an
+        array of NaN values.
+
+    Returns
+    -------
+    inv : array-like
+        Inverse of ``arr``
+
     """
+    # Check that matrix is not singular
+    rank = np.linalg.matrix_rank(arr, tol=tol)
+    if rank < np.max(arr.shape):
+        return arr * np.nan
+
     try:
         from numpy.linalg import inv
         _inv = inv(arr)
-    except:
+    except ImportError:
         _inv = np.matrix(arr).I.A
     
     return _inv


Test patch:
diff --git a/eazy/tests/test_utils.py b/eazy/tests/test_utils.py
index a14504a..879b340 100644
--- a/eazy/tests/test_utils.py
+++ b/eazy/tests/test_utils.py
@@ -92,6 +92,11 @@ def test_invert():
     ainv = utils.safe_invert(a)
     assert(np.allclose(np.dot(a, ainv), np.eye(2)))
 
+    # Singular matrix returns np.nan inverse without raising exception
+    a = np.ones((2,2))
+    ainv = utils.safe_invert(a)
+    assert np.all(~np.isfinite(utils.safe_invert(a)))
+
 
 def test_query_string():
     """