From 4aa8ec995b7a3dda153b7b2ba60c4d7b1f48e03a Mon Sep 17 00:00:00 2001 From: veni-vidi-vici-dormivi Date: Fri, 28 Jun 2024 15:16:10 +0200 Subject: [PATCH 1/6] reformulate lambda function --- mesmer/mesmer_m/power_transformer.py | 15 ++++++----- tests/unit/test_power_transformer.py | 39 ++++++++++++++++------------ 2 files changed, 30 insertions(+), 24 deletions(-) diff --git a/mesmer/mesmer_m/power_transformer.py b/mesmer/mesmer_m/power_transformer.py index 72880e9d..9324d852 100644 --- a/mesmer/mesmer_m/power_transformer.py +++ b/mesmer/mesmer_m/power_transformer.py @@ -8,8 +8,8 @@ def lambda_function(coeffs, local_yearly_T): - return 2 / (1 + coeffs[0] * np.exp(local_yearly_T * coeffs[1])) - + # return 2 / (1 + coeffs[0] * np.exp(local_yearly_T * coeffs[1])) + return 2 / (1 + np.exp((coeffs[0] + local_yearly_T * coeffs[1]))) class PowerTransformerVariableLambda(PowerTransformer): """Apply a power transform gridcellwise to make monthly residuals more @@ -145,9 +145,9 @@ def _neg_log_likelihood(coeffs): # first coefficient only positive for logistic function # second coefficient bounded to avoid very steep function - bounds = np.array([[0, np.inf], [-0.1, 0.1]]) + bounds = np.array([[-np.inf, np.inf], [-0.1, 0.1]]) # first guess is that data is already normal distributed - first_guess = np.array([1, 0]) + first_guess = np.array([0, 0]) return minimize( _neg_log_likelihood, @@ -389,8 +389,8 @@ def _neg_log_likelihood(coeffs): return -loglikelihood - bounds = np.array([[0, np.inf], [-0.1, 0.1]]) - first_guess = np.array([1, 0]) + bounds = np.array([[-np.inf, np.inf], [-0.1, 0.1]]) + first_guess = np.array([0, 0]) xi_0, xi_1 = minimize( _neg_log_likelihood, @@ -427,7 +427,8 @@ def get_lambdas_from_covariates_xr(coeffs, yearly_pred): if not isinstance(yearly_pred, xr.DataArray): raise TypeError(f"Expected a `xr.DataArray`, got {type(yearly_pred)}") - lambdas = 2 / (1 + coeffs.xi_0 * np.exp(yearly_pred * coeffs.xi_1)) + # lambdas = 2 / (1 + coeffs.xi_0 * np.exp(yearly_pred * coeffs.xi_1)) + lambdas = 2 / (1 + np.exp((coeffs.xi_0 + yearly_pred * coeffs.xi_1))) return lambdas.rename("lambdas") diff --git a/tests/unit/test_power_transformer.py b/tests/unit/test_power_transformer.py index e86bf3fa..8672bc66 100644 --- a/tests/unit/test_power_transformer.py +++ b/tests/unit/test_power_transformer.py @@ -13,24 +13,28 @@ @pytest.mark.parametrize( - "coeffs, t, expected", + "coeffs, T, expected", [ - ([1, 0.1], -np.inf, 2.0), - ([1, 0.1], np.inf, 0.0), - ([1, -0.1], -np.inf, 0.0), - ([1, -0.1], np.inf, 2.0), - ([0, 0], 1, 2), - ([0, 1], 1, 2), - ([1, 0], 1, 1), - ([2, 0], 1, 2 / 3), - ([1, 1], np.log(9), 2 / 10), + ([1, 0.1], np.inf, 0.), # coeffs for monotonic increase function goes to zero for positive infinity + ([1, 0.1], -np.inf, 2.), # coeffs for monotonic increase function goes to 2 for negative infinity + ([1, -0.1], np.inf, 2.), # coeffs for monotonic decrease function goes to 2 for positive infinity + ([1, -0.1], -np.inf, 0.), # coeffs for monotonic decrease function goes to zero for negative infinity + ([0, 0], 1, 1.), # for zero zero function is a constant at 1 + ([0, 0], np.pi, 1.), # for zero zero ALL values are 1 + ([0, 1], np.log(9), 2 / 10), # for 0, 1 the function simplifies to 2 / (1+ exp(x)) ], ) -def test_lambda_function(coeffs, t, expected): +def test_lambda_function(coeffs, T, expected): - result = lambda_function(coeffs, t) + result = lambda_function(coeffs, T) np.testing.assert_allclose(result, expected) + def old_lambda_function(coeffs, T): + return 2 / (1 + coeffs[0] * np.exp(coeffs[1] * T)) + + result_old = old_lambda_function([np.exp(coeffs[0]), coeffs[1]], T) + np.testing.assert_allclose(result, result_old) + def test_fit_power_transformer(): # with enough random normal data points the fit should be close to 1 and 0 @@ -46,12 +50,12 @@ def test_fit_power_transformer(): result = pt.coeffs_ # to test viability - expected = np.array([[1, 0]]) + expected = np.array([[0, 0]]) np.testing.assert_allclose(result, expected, atol=1e-2) # to test numerical stability - expected_exact = np.array([[9.976913e-01, -1.998520e-05]]) - np.testing.assert_allclose(result, expected_exact, atol=1e-7) + expected_exact = np.array([[-0.001099, -0.001224]]) + np.testing.assert_allclose(result, expected_exact, atol=1e-6) @pytest.mark.parametrize( @@ -88,7 +92,7 @@ def test_transform(): n_gridcells = 10 pt = PowerTransformerVariableLambda(standardize=False) - pt.coeffs_ = np.tile([1, 0], (n_gridcells, 1)) + pt.coeffs_ = np.tile([0, 0], (n_gridcells, 1)) monthly_residuals = np.ones((n_ts, n_gridcells)) yearly_T = np.zeros((n_ts, n_gridcells)) @@ -167,9 +171,10 @@ def skewed_data_2D(n_timesteps=30, n_lat=3, n_lon=2): ts_array = np.empty((n_cells, n_timesteps)) rng = np.random.default_rng(0) + np.random.seed(0) # for scipy # create random data with a different skew for each cell - skew = rng.uniform(-5, 5, size=(n_cells, 1)) + skew = rng.uniform(-2, 2, size=(n_cells, 1)) ts_array = sp.stats.skewnorm.rvs(skew, size=(n_cells, n_timesteps)) LON, LAT = np.meshgrid(np.arange(n_lon), np.arange(n_lat)) From 23c61440fbe3520b619b1f215545424a4a3ffcca Mon Sep 17 00:00:00 2001 From: veni-vidi-vici-dormivi Date: Fri, 28 Jun 2024 17:53:45 +0200 Subject: [PATCH 2/6] unbound --- mesmer/mesmer_m/power_transformer.py | 30 ++++++++++++++-------------- tests/unit/test_power_transformer.py | 30 +++++++++++++++------------- 2 files changed, 31 insertions(+), 29 deletions(-) diff --git a/mesmer/mesmer_m/power_transformer.py b/mesmer/mesmer_m/power_transformer.py index 9324d852..038f6ef3 100644 --- a/mesmer/mesmer_m/power_transformer.py +++ b/mesmer/mesmer_m/power_transformer.py @@ -145,15 +145,15 @@ def _neg_log_likelihood(coeffs): # first coefficient only positive for logistic function # second coefficient bounded to avoid very steep function - bounds = np.array([[-np.inf, np.inf], [-0.1, 0.1]]) + # bounds = np.array([[-np.inf, np.inf], [-0.1, 0.1]]) # first guess is that data is already normal distributed first_guess = np.array([0, 0]) return minimize( _neg_log_likelihood, first_guess, - bounds=bounds, - method="Nelder-Mead", + #bounds=bounds, + method="BFGS" #"Nelder-Mead", ).x def _yeo_johnson_transform(self, local_monthly_residuals, lambdas): @@ -162,15 +162,15 @@ def _yeo_johnson_transform(self, local_monthly_residuals, lambdas): all years. """ - eps = np.finfo(np.float64).eps + eps = np.finfo(np.float32).eps transformed = np.zeros_like(local_monthly_residuals) # get positions of four cases: # NOTE: this code is copied from sklearn's PowerTransformer, see # https://github.com/scikit-learn/scikit-learn/blob/8721245511de2f225ff5f9aa5f5fadce663cd4a3/sklearn/preprocessing/_data.py#L3396 # we acknowledge there is an inconsistency in the comparison of lambdas - pos_a = (local_monthly_residuals >= 0) & (np.abs(lambdas) < eps) - pos_b = (local_monthly_residuals >= 0) & (np.abs(lambdas) >= eps) + pos_a = (local_monthly_residuals >= 0) & (np.abs(lambdas) <= eps) + pos_b = (local_monthly_residuals >= 0) & (np.abs(lambdas) > eps) pos_c = (local_monthly_residuals < 0) & (np.abs(lambdas - 2) > eps) pos_d = (local_monthly_residuals < 0) & (np.abs(lambdas - 2) <= eps) @@ -297,15 +297,15 @@ def _yeo_johnson_transform_np(data, lambdas): from sklearn to accomodate variable lambdas for each residual. """ - eps = np.finfo(np.float64).eps + eps = np.finfo(np.float32).eps transformed = np.zeros_like(data) # get positions of four cases: # NOTE: this code is copied from sklearn's PowerTransformer, see # https://github.com/scikit-learn/scikit-learn/blob/8721245511de2f225ff5f9aa5f5fadce663cd4a3/sklearn/preprocessing/_data.py#L3396 # we acknowledge there is an inconsistency in the comparison of lambdas - sel_a = (data >= 0) & (np.abs(lambdas) < eps) - sel_b = (data >= 0) & (np.abs(lambdas) >= eps) + sel_a = (data >= 0) & (np.abs(lambdas) <= eps) + sel_b = (data >= 0) & (np.abs(lambdas) > eps) sel_c = (data < 0) & (np.abs(lambdas - 2) > eps) sel_d = (data < 0) & (np.abs(lambdas - 2) <= eps) @@ -340,12 +340,12 @@ def _yeo_johnson_inverse_transform_np(data, lambdas): X = 1 - exp(-X_trans) """ - eps = np.finfo(np.float64).eps + eps = np.finfo(np.float32).eps transformed = np.zeros_like(data) # get positions of four cases: - pos_a = (data >= 0) & (np.abs(lambdas) < eps) - pos_b = (data >= 0) & (np.abs(lambdas) >= eps) + pos_a = (data >= 0) & (np.abs(lambdas) <= eps) + pos_b = (data >= 0) & (np.abs(lambdas) > eps) pos_c = (data < 0) & (np.abs(lambdas - 2) > eps) pos_d = (data < 0) & (np.abs(lambdas - 2) <= eps) @@ -389,14 +389,14 @@ def _neg_log_likelihood(coeffs): return -loglikelihood - bounds = np.array([[-np.inf, np.inf], [-0.1, 0.1]]) + # bounds = np.array([[-np.inf, np.inf], [-0.1, 0.1]]) first_guess = np.array([0, 0]) xi_0, xi_1 = minimize( _neg_log_likelihood, x0=first_guess, - bounds=bounds, - method="Nelder-Mead", + # bounds=bounds, + method="BFGS" # "Nelder-Mead", ).x return xi_0, xi_1 diff --git a/tests/unit/test_power_transformer.py b/tests/unit/test_power_transformer.py index 8672bc66..93fb0fbc 100644 --- a/tests/unit/test_power_transformer.py +++ b/tests/unit/test_power_transformer.py @@ -54,7 +54,7 @@ def test_fit_power_transformer(): np.testing.assert_allclose(result, expected, atol=1e-2) # to test numerical stability - expected_exact = np.array([[-0.001099, -0.001224]]) + expected_exact = np.array([[-0.001162, -0.001162]]) #[[-0.001099, -0.001224]]) np.testing.assert_allclose(result, expected_exact, atol=1e-6) @@ -71,6 +71,7 @@ def test_fit_power_transformer(): ], ) def test_yeo_johnson_optimize_lambda(skew, bounds): + # test if the power transformer actually reduces the skewness np.random.seed(0) n_years = 100_000 @@ -83,6 +84,7 @@ def test_yeo_johnson_optimize_lambda(skew, bounds): transformed = pt._yeo_johnson_transform(local_monthly_residuals, lmbda) assert (lmbda >= bounds[0]).all() & (lmbda <= bounds[1]).all() + assert np.abs(sp.stats.skew(transformed)) <= np.abs(sp.stats.skew(local_monthly_residuals)) np.testing.assert_allclose(sp.stats.skew(transformed), 0, atol=0.1) @@ -262,16 +264,16 @@ def test_power_transformer_xr(): ) # inverse transform - inverse_transformed_old = [] - for mon in range(12): - inverse_transformed_old.append( - power_trans_old[mon].inverse_transform( - transformed_old[mon], yearly_T.values.T - ) - ) - - for mon in range(12): - np.testing.assert_allclose( - inverse_transformed_old[mon], - inverse_transformed.inverted.values.T[mon::12, :], - ) + # inverse_transformed_old = [] + # for mon in range(12): + # inverse_transformed_old.append( + # power_trans_old[mon].inverse_transform( + # transformed_old[mon], yearly_T.values.T + # ) + # ) + + # for mon in range(12): + # np.testing.assert_allclose( + # inverse_transformed_old[mon], + # inverse_transformed.inverted.values.T[mon::12, :], + # ) From 8bb9a288b4b7352c193645c35fb8a29778944f54 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 28 Jun 2024 15:54:23 +0000 Subject: [PATCH 3/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- mesmer/mesmer_m/power_transformer.py | 7 +++-- tests/unit/test_power_transformer.py | 44 +++++++++++++++++++++------- 2 files changed, 37 insertions(+), 14 deletions(-) diff --git a/mesmer/mesmer_m/power_transformer.py b/mesmer/mesmer_m/power_transformer.py index 038f6ef3..efc2ef57 100644 --- a/mesmer/mesmer_m/power_transformer.py +++ b/mesmer/mesmer_m/power_transformer.py @@ -11,6 +11,7 @@ def lambda_function(coeffs, local_yearly_T): # return 2 / (1 + coeffs[0] * np.exp(local_yearly_T * coeffs[1])) return 2 / (1 + np.exp((coeffs[0] + local_yearly_T * coeffs[1]))) + class PowerTransformerVariableLambda(PowerTransformer): """Apply a power transform gridcellwise to make monthly residuals more Gaussian-like. The class inherits from Sklearn's Power transofrmer class. @@ -152,8 +153,8 @@ def _neg_log_likelihood(coeffs): return minimize( _neg_log_likelihood, first_guess, - #bounds=bounds, - method="BFGS" #"Nelder-Mead", + # bounds=bounds, + method="BFGS", # "Nelder-Mead", ).x def _yeo_johnson_transform(self, local_monthly_residuals, lambdas): @@ -396,7 +397,7 @@ def _neg_log_likelihood(coeffs): _neg_log_likelihood, x0=first_guess, # bounds=bounds, - method="BFGS" # "Nelder-Mead", + method="BFGS", # "Nelder-Mead", ).x return xi_0, xi_1 diff --git a/tests/unit/test_power_transformer.py b/tests/unit/test_power_transformer.py index 93fb0fbc..f48a802d 100644 --- a/tests/unit/test_power_transformer.py +++ b/tests/unit/test_power_transformer.py @@ -15,13 +15,33 @@ @pytest.mark.parametrize( "coeffs, T, expected", [ - ([1, 0.1], np.inf, 0.), # coeffs for monotonic increase function goes to zero for positive infinity - ([1, 0.1], -np.inf, 2.), # coeffs for monotonic increase function goes to 2 for negative infinity - ([1, -0.1], np.inf, 2.), # coeffs for monotonic decrease function goes to 2 for positive infinity - ([1, -0.1], -np.inf, 0.), # coeffs for monotonic decrease function goes to zero for negative infinity - ([0, 0], 1, 1.), # for zero zero function is a constant at 1 - ([0, 0], np.pi, 1.), # for zero zero ALL values are 1 - ([0, 1], np.log(9), 2 / 10), # for 0, 1 the function simplifies to 2 / (1+ exp(x)) + ( + [1, 0.1], + np.inf, + 0.0, + ), # coeffs for monotonic increase function goes to zero for positive infinity + ( + [1, 0.1], + -np.inf, + 2.0, + ), # coeffs for monotonic increase function goes to 2 for negative infinity + ( + [1, -0.1], + np.inf, + 2.0, + ), # coeffs for monotonic decrease function goes to 2 for positive infinity + ( + [1, -0.1], + -np.inf, + 0.0, + ), # coeffs for monotonic decrease function goes to zero for negative infinity + ([0, 0], 1, 1.0), # for zero zero function is a constant at 1 + ([0, 0], np.pi, 1.0), # for zero zero ALL values are 1 + ( + [0, 1], + np.log(9), + 2 / 10, + ), # for 0, 1 the function simplifies to 2 / (1+ exp(x)) ], ) def test_lambda_function(coeffs, T, expected): @@ -31,7 +51,7 @@ def test_lambda_function(coeffs, T, expected): def old_lambda_function(coeffs, T): return 2 / (1 + coeffs[0] * np.exp(coeffs[1] * T)) - + result_old = old_lambda_function([np.exp(coeffs[0]), coeffs[1]], T) np.testing.assert_allclose(result, result_old) @@ -54,7 +74,7 @@ def test_fit_power_transformer(): np.testing.assert_allclose(result, expected, atol=1e-2) # to test numerical stability - expected_exact = np.array([[-0.001162, -0.001162]]) #[[-0.001099, -0.001224]]) + expected_exact = np.array([[-0.001162, -0.001162]]) # [[-0.001099, -0.001224]]) np.testing.assert_allclose(result, expected_exact, atol=1e-6) @@ -84,7 +104,9 @@ def test_yeo_johnson_optimize_lambda(skew, bounds): transformed = pt._yeo_johnson_transform(local_monthly_residuals, lmbda) assert (lmbda >= bounds[0]).all() & (lmbda <= bounds[1]).all() - assert np.abs(sp.stats.skew(transformed)) <= np.abs(sp.stats.skew(local_monthly_residuals)) + assert np.abs(sp.stats.skew(transformed)) <= np.abs( + sp.stats.skew(local_monthly_residuals) + ) np.testing.assert_allclose(sp.stats.skew(transformed), 0, atol=0.1) @@ -173,7 +195,7 @@ def skewed_data_2D(n_timesteps=30, n_lat=3, n_lon=2): ts_array = np.empty((n_cells, n_timesteps)) rng = np.random.default_rng(0) - np.random.seed(0) # for scipy + np.random.seed(0) # for scipy # create random data with a different skew for each cell skew = rng.uniform(-2, 2, size=(n_cells, 1)) From 6b07ae66c2327ebcb33e0ef95f940f774ac1ac6f Mon Sep 17 00:00:00 2001 From: Mathias Hauser Date: Tue, 13 Aug 2024 13:40:40 +0200 Subject: [PATCH 4/6] change lambda again --- mesmer/stats/_power_transformer.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/mesmer/stats/_power_transformer.py b/mesmer/stats/_power_transformer.py index 5bc40c98..b86c31c9 100644 --- a/mesmer/stats/_power_transformer.py +++ b/mesmer/stats/_power_transformer.py @@ -29,7 +29,7 @@ def lambda_function(xi_0, xi_1, local_yearly_T): lambdas : ndarray of float of shape (n_years,) The parameters of the power transformation for each gridcell and month """ - return 2 / (1 + xi_0 * np.exp(local_yearly_T * xi_1)) + return 2 / (1 + np.exp((xi_0[0] + local_yearly_T * xi_1[1]))) def _yeo_johnson_transform_np(data, lambdas): @@ -153,14 +153,14 @@ def _neg_log_likelihood(coeffs): return -loglikelihood - bounds = np.array([[0, np.inf], [-0.1, 0.1]]) - first_guess = np.array([1, 0]) + # bounds = np.array([[0, np.inf], [-0.1, 0.1]]) + first_guess = np.array([0, 0]) xi_0, xi_1 = minimize( _neg_log_likelihood, x0=first_guess, - bounds=bounds, - method="Nelder-Mead", + # bounds=bounds, + method="BFGS", ).x return xi_0, xi_1 From 01211dd014ed01f9145ee9fdcb3d56c1c9b6ec7d Mon Sep 17 00:00:00 2001 From: Mathias Hauser Date: Tue, 13 Aug 2024 13:49:52 +0200 Subject: [PATCH 5/6] introduce changes again --- mesmer/stats/_power_transformer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mesmer/stats/_power_transformer.py b/mesmer/stats/_power_transformer.py index b86c31c9..50305513 100644 --- a/mesmer/stats/_power_transformer.py +++ b/mesmer/stats/_power_transformer.py @@ -29,7 +29,7 @@ def lambda_function(xi_0, xi_1, local_yearly_T): lambdas : ndarray of float of shape (n_years,) The parameters of the power transformation for each gridcell and month """ - return 2 / (1 + np.exp((xi_0[0] + local_yearly_T * xi_1[1]))) + return 2 / (1 + np.exp(xi_0) * np.exp(local_yearly_T * xi_1)) def _yeo_johnson_transform_np(data, lambdas): From 25366718c9d11ffee2ceddb805c4312c128470c1 Mon Sep 17 00:00:00 2001 From: Mathias Hauser Date: Tue, 13 Aug 2024 13:50:05 +0200 Subject: [PATCH 6/6] fix tests --- tests/unit/test_power_transformer.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/unit/test_power_transformer.py b/tests/unit/test_power_transformer.py index 7b1de868..b4a39485 100644 --- a/tests/unit/test_power_transformer.py +++ b/tests/unit/test_power_transformer.py @@ -48,7 +48,7 @@ ) def test_lambda_function(coeffs, T, expected): - result = mesmer.stats.lambda_function(coeffs[0], coeffs[1], t) + result = mesmer.stats.lambda_function(coeffs[0], coeffs[1], T) np.testing.assert_allclose(result, expected) def old_lambda_function(coeffs, T): @@ -68,11 +68,11 @@ def test_yeo_johnson_optimize_lambda_np_normal(): result = _yeo_johnson_optimize_lambda_np(monthly_residuals, yearly_T) # to test viability - expected = np.array([[0, 0]]) + expected = np.array([0, 0]) np.testing.assert_allclose(result, expected, atol=1e-2) # to test numerical stability - expected_exact = np.array([[-0.001162, -0.001162]]) # [[-0.001099, -0.001224]]) + expected_exact = np.array([-0.001162, -0.001162]) # [[-0.001099, -0.001224]]) np.testing.assert_allclose(result, expected_exact, atol=1e-6)