diff --git a/smrt/rtsolver/dort.py b/smrt/rtsolver/dort.py index a2848fa..d1e3e8d 100644 --- a/smrt/rtsolver/dort.py +++ b/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). diff --git a/smrt/rtsolver/test_dort.py b/smrt/rtsolver/test_dort.py index 85d1679..27aad6a 100644 --- a/smrt/rtsolver/test_dort.py +++ b/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, diff --git a/smrt/test/test_dmrtdort.py b/smrt/test/test_dmrtdort.py index ba960c1..d84bf1a 100644 --- a/smrt/test/test_dmrtdort.py +++ b/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) diff --git a/smrt/test/test_integration_iba.py b/smrt/test/test_integration_iba.py index eae8f9e..c32ff10 100644 --- a/smrt/test/test_integration_iba.py +++ b/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