Changed default diagonalization_method to shur_forcedtriu in DORT

smrt/rtsolver/dort.py

@@ -21,17 +21,17 @@
     The DORT solver is very robust in passive mode but may raise exception in active mode due to a matrix
     diagonalisation problem. The exception provides detailed information on how to address this issue. Two new
     diagonalisation approaches were added in January 2024. They are activated by setting the `diagonalization_method`
-    optional argument (see :py:mod:`smrt.core.make_model`). The first method (``diagonalization_method='shur'``)
-    replaces the scipy.linalg.eig function by a shur decomposition followed by a diagonalisation of the shur matrix.
-    While scipy.linalg.eig performs such a shur decomposition internally in any case, it seems that explicitly calling
-    the shur decomposition beforehand improves the stability. Nevertheless to really solve the problem, the second
-    method (``diagonalization_method='shur_forcedtriu'``) consists in removing the 2x2 and 3x3 blocks from the shur
-    matrix, i.e. forcing the shur matrix to be upper triangular (triu in numpy jargon = zeroing the lower part of this
+    optional argument (see :py:mod:`smrt.core.make_model`). The first method (``diagonalization_method='schur'``)
+    replaces the scipy.linalg.eig function by a schur decomposition followed by a diagonalisation of the schur matrix.
+    While scipy.linalg.eig performs such a schur decomposition internally in any case, it seems that explicitly calling
+    the schur decomposition beforehand improves the stability. Nevertheless to really solve the problem, the second
+    method (``diagonalization_method='schur_forcedtriu'``) consists in removing the 2x2 and 3x3 blocks from the schur
+    matrix, i.e. forcing the schur matrix to be upper triangular (triu in numpy jargon = zeroing the lower part of this
     matrix). This problem is due to the structure of the matrix to be diagonalized and the formulation of the DORT
     method in the polarimetric configuration. The eigenvalues come by triplets and can be very close to each other for
     the three H, V, U Stokes components when scattering is becoming small (or equiv. the azimuth mode 'm' is large). As
     a consequence of the Gershgorin theorem, this results in slightly complex eigenvalues (i.e. eigenvalues with very
-    small imaginary part) that comes from 2x2 or 3x3 blocks in the shur decomposition. This would not be a problem if
+    small imaginary part) that comes from 2x2 or 3x3 blocks in the schur decomposition. This would not be a problem if
     the eigenvectors were correctly estimated, but this is not the case. It is indeed difficult to find the correct
     orientation of eigenvectors associated to very close eigenvalues. To overcome the problem, the solution is to remove
     the 2x2 and 3x3 blocks. In principle, it would be safer to check that these blocks are nearly diagonal but this is
@@ -119,8 +119,8 @@ class DORT(RTSolverBase, CoherentLayerMixin, DiscreteOrdinatesMixin, PlanckMixin
             contribution of the lowest layers is neglegible. The optical depth is a good criteria to determine this
             limit. A value of about 6 is recommended. Use with care, especially values lower than 6.
         diagonalization_method: This value set the method for the diagonalization in the eigenvalue solver. The defaut
-            is "eig", it uses the scipy.linalg.eig function. The "shur" replaces the scipy.linalg.eig function by a shur
-            decomposition followed by a diagonalisation of the shur matrix. The "shur_forcedtriu" forces the shur matrix
+            is "eig", it uses the scipy.linalg.eig function. The "schur" replaces the scipy.linalg.eig function by a schur
+            decomposition followed by a diagonalisation of the schur matrix. The "schur_forcedtriu" forces the schur matrix
             to be upper triangular. The "half_rank_eig" is the fastest method but requires symmetry and energy
             conservation which may fail with some EMModels and for some parameters. The "stamnes88" is another half rank
             fast method.
@@ -155,7 +155,7 @@ def __init__(
         error_handling="exception",
         process_coherent_layers=False,
         prune_deep_snowpack=None,
-        diagonalization_method="eig",
+        diagonalization_method="schur_forcedtriu",
         diagonalization_cache=False,
         rayleigh_jeans_approximation=False,
     ):
@@ -778,8 +778,8 @@ def __init__(self, ke, ks, ft_even_phase_function, mu, weight, m_max, normalizat
             normalization (str or bool): Phase-function normalization option (e.g. ``'auto'``, ``'forced'``,
                 or ``False``).
             symmetrization (bool): If True, enforce phase-function symmetry.
-            method (str): Diagonalization method. Supported values: ``'eig'``, ``'shur'``,
-                ``'shur_forcedtriu'``, ``'half_rank_eig'``, ``'stamnes88'``.
+            method (str): Diagonalization method. Supported values: ``'eig'``, ``'schur'``,
+                ``'schur_forcedtriu'``, ``'half_rank_eig'``, ``'stamnes88'``.
             cache (bool): If True, cache diagonalization results to avoid recomputation.
 
         """
@@ -808,10 +808,10 @@ def __init__(self, ke, ks, ft_even_phase_function, mu, weight, m_max, normalizat
         match self.method:
             case "eig":
                 self.diagonalize_function = self.diagonalize_eig
-            case "shur":
-                self.diagonalize_function = partial(self.diagonalize_shur, force_triu=False)
-            case "shur_forcedtriu":
-                self.diagonalize_function = partial(self.diagonalize_shur, force_triu=True)
+            case "schur":
+                self.diagonalize_function = partial(self.diagonalize_schur, force_triu=False)
+            case "schur_forcedtriu":
+                self.diagonalize_function = partial(self.diagonalize_schur, force_triu=True)
             case "half_rank_eig":
                 self.compute_half_rank_phase = True
                 self.diagonalize_function = self.diagonalize_half_rank_eig
@@ -977,19 +977,19 @@ def diagonalize_eig(self, m, A):
 
         return self.validate_eigen(beta, Eu, Ed, m)
 
-    def diagonalize_shur(self, m, A, force_triu=False):
-        # diagonalise the matrix. Eq (13) using Shur decomposition. This avoids some instabilities with the direct eig
+    def diagonalize_schur(self, m, A, force_triu=False):
+        # diagonalise the matrix. Eq (13) using schur decomposition. This avoids some instabilities with the direct eig
         # function
         # in addition it is possible to remove the 2x2 or 3x3 blocks that occurs when eigenvalues are close
-        # forcing the lower triangular part of the shur matrix to zero solves this problem but is radical
+        # forcing the lower triangular part of the schur matrix to zero solves this problem but is radical
         # a better algorithm would first check that the 2x2 nd 3x3 blocks are nearly diagional (values are very small)
-
         try:
             T, Z = scipy.linalg.schur(A)
         except scipy.linalg.LinAlgError:
             raise SMRTError("Schur decomposition failed.\n" + self.diagonalization_error_message())
         try:
             if force_triu:
+                # Zero out strictly lower triangular values
                 T[np.tril_indices(T.shape[0], k=-1)] = 0
 
             beta, E = scipy.linalg.eig(T, overwrite_a=True)
@@ -1256,9 +1256,9 @@ def diagonalization_error_message(self):
 
 - almost diagonal matrix. Such a matrix often arises in active mode when m_max is quite high. However it can
 also arises in passive mode or with low m_max. To solve this issue  you can try to activate the
-diagonalization_method="shur" option and if it does not work the more radical diagonalization_method="shur_forcedtriu".
+diagonalization_method="schur" option and if it does not work the more radical diagonalization_method="schur_forcedtriu".
 These options are experimental, please report your results (both success and failure).
-diagonalization_method="shur_forcedtriu" should become the default if success are reported. Alternatively you could
+diagonalization_method="schur_forcedtriu" should become the default if success are reported. Alternatively you could
 reduce the m_max option progressively but high values of m_max give more accurate results in active mode (but tends to
 produce almost diagonal matrix).
 

smrt/rtsolver/test_dort.py

@@ -127,9 +127,9 @@ def test_shallow_snowpack():
 
 @pytest.mark.parametrize(
     "microstructure_model,m_max,emmodel,diagonalization_method",
-    [("independent_sphere", 6, Rayleigh, "shur"), ("exponential", 16, IBA, "shur_forcedtriu")],
+    [("independent_sphere", 6, Rayleigh, "schur"), ("exponential", 16, IBA, "schur_forcedtriu")],
 )
-def test_shur_based_diagonalisation(microstructure_model, m_max, emmodel, diagonalization_method):
+def test_schur_based_diagonalisation(microstructure_model, m_max, emmodel, diagonalization_method):
     sp = make_snowpack(
         thickness=[1000],
         microstructure_model=microstructure_model,
@@ -142,7 +142,7 @@ def test_shur_based_diagonalisation(microstructure_model, m_max, emmodel, diagon
     nstreams = 32
 
     # this setting fails when DORT  use scipy.linalg.eig
-    # but this works with the shur diagonalization. Let check this:
+    # but this works with the schur diagonalization. Let check this:
 
     m = Model(
         emmodel,

smrt/test/test_dmrtdort.py

@@ -86,8 +86,8 @@ def test_less_refringent_bottom_layer():
         stickiness=0.2,
         substrate=make_soil("transparent", 1, 270),
     )
-    # this test fails with some version of scipy if not using the shur method
-    m = make_model("dmrt_qcacp_shortrange", "dort", rtsolver_options=dict(diagonalization_method="shur_forcedtriu"))
+    # this test fails with some version of scipy if not using the schur method
+    m = make_model("dmrt_qcacp_shortrange", "dort", rtsolver_options=dict(diagonalization_method="schur_forcedtriu"))
     scat = active(10e9, 45)
     warnings.simplefilter("ignore", category=SMRTWarning)
     res = m.run(scat, snowpack)
@@ -99,7 +99,7 @@ def test_less_refringent_bottom_layer():
 def test_less_refringent_bottom_layer_VH():
     # Regression test 19-03-2018: value may change if other bugs found
     snowpack = make_snowpack([0.2, 0.3], "sticky_hard_spheres", density=[290.0, 250.0], radius=1e-4, stickiness=0.2)
-    m = make_model("dmrt_qcacp_shortrange", "dort", rtsolver_options=dict(diagonalization_method="shur_forcedtriu"))
+    m = make_model("dmrt_qcacp_shortrange", "dort", rtsolver_options=dict(diagonalization_method="schur_forcedtriu"))
     scat = active(10e9, 45)
     warnings.simplefilter("ignore", category=SMRTWarning)
     res = m.run(scat, snowpack)

smrt/test/test_integration_iba.py

@@ -31,7 +31,7 @@ def setup_snowpack_2():
 
 
 @pytest.mark.parametrize(
-    "method,atol", [("eig", 0.0001), ("shur", 0.0001), ("half_rank_eig", 0.0001), ("stamnes88", 0.4)]
+    "method,atol", [("eig", 0.0001), ("schur", 0.0001), ("half_rank_eig", 0.0001), ("stamnes88", 0.4)]
 )
 def test_iba_dort_oneconfig_passive(setup_snowpack_2, method, atol):
     # prepare inputs
@@ -50,7 +50,7 @@ def test_iba_dort_oneconfig_passive(setup_snowpack_2, method, atol):
 
 
 @pytest.mark.parametrize(
-    "method,atol", [("eig", 0.001), ("shur", 0.0001), ("half_rank_eig", 0.001), ("stamnes88", 1.2)]
+    "method,atol", [("eig", 0.001), ("schur", 0.0001), ("half_rank_eig", 0.001), ("stamnes88", 1.2)]
 )
 def test_iba_dort_oneconfig_active(setup_snowpack_2, method, atol):
     # prepare inputs