MESMER-group / mesmer

spatially-resolved ESM-specific multi-scenario initial-condition ensemble emulator
https://mesmer-emulator.readthedocs.io/en/latest/
GNU General Public License v3.0
22 stars 15 forks source link

align equality check in yeo johnson transform #436

Closed veni-vidi-vici-dormivi closed 1 month ago

veni-vidi-vici-dormivi commented 2 months ago

Changing the comparison to eps in _yeo_johnson_transform to be consistent with sklearn's yet-johnsons power transform.

codecov[bot] commented 2 months ago

Codecov Report

Attention: Patch coverage is 0% with 2 lines in your changes are missing coverage. Please review.

Project coverage is 81.94%. Comparing base (9b0b76b) to head (0bfe836). Report is 37 commits behind head on main.

Files Patch % Lines
mesmer/mesmer_m/power_transformer.py 0.00% 2 Missing :warning:
Additional details and impacted files ```diff @@ Coverage Diff @@ ## main #436 +/- ## ========================================== - Coverage 87.90% 81.94% -5.96% ========================================== Files 40 44 +4 Lines 1745 1939 +194 ========================================== + Hits 1534 1589 +55 - Misses 211 350 +139 ``` | [Flag](https://app.codecov.io/gh/MESMER-group/mesmer/pull/436/flags?src=pr&el=flags&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=MESMER-group) | Coverage Δ | | |---|---|---| | [unittests](https://app.codecov.io/gh/MESMER-group/mesmer/pull/436/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=MESMER-group) | `81.94% <0.00%> (-5.96%)` | :arrow_down: | Flags with carried forward coverage won't be shown. [Click here](https://docs.codecov.io/docs/carryforward-flags?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=MESMER-group#carryforward-flags-in-the-pull-request-comment) to find out more.

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

mathause commented 2 months ago

Can you define

eps = np.spacing(np.float64(1))

or actually, if we don't do it for an array

eps = np.finfo(np.float64).eps

? That would be clearer to me (not sure if here or in another PR, though).

veni-vidi-vici-dormivi commented 2 months ago

I think another PR would make sense.

veni-vidi-vici-dormivi commented 2 months ago

Then let's merge this then?

mathause commented 2 months ago

Can you merge main - I need to think about this.

veni-vidi-vici-dormivi commented 2 months ago

The question is, if when lambda is exactly eps, if we consider that to be zero, or not. I just wanted to make it consistent with the sklearn function. They however are not consistent themselves I think because they write if abs(lmbda) < np.spacing(1.0) for $\lambda = 0$ but then if not abs(lmbda - 2) > np.spacing(1.0) for $\lambda = 2$...

# when x >= 0
if abs(lmbda) < np.spacing(1.0):
    out[pos] = np.log1p(x[pos])
else:  # lmbda != 0
    out[pos] = (np.power(x[pos] + 1, lmbda) - 1) / lmbda

# when x < 0
if abs(lmbda - 2) > np.spacing(1.0):
    out[~pos] = -(np.power(-x[~pos] + 1, 2 - lmbda) - 1) / (2 - lmbda)
else:  # lmbda == 2
    out[~pos] = -np.log1p(-x[~pos])

Maybe that's the reason why Shruti wrote it herself? Because for the inverse transform she actually just uses the once from sklearn...

https://github.com/MESMER-group/mesmer/blob/03ab48cacc324782f7e12b1a0227abc87b7b95b4/mesmer/mesmer_m/power_transformer.py#L272

mathause commented 2 months ago

Shruti wrote this herself so she could have variable (or dependent) $\lambda$ values.

I think the >= / > mess comes because the original author also confused that $\lambda$ can be any real value and does not have to be in 0..2 (also https://github.com/MESMER-group/mesmer/pull/430#discussion_r1590938808). The value problem was fixed but this particular inconsistency remained - see scikit-learn/scikit-learn#12522 (Assuming 0 <= lambda <= 2 the operators make sense.)

mathause commented 2 months ago

scipy does the same but I think it's written by the same author: https://github.com/scipy/scipy/blob/7dcd8c59933524986923cde8e9126f5fc2e6b30b/scipy/stats/_morestats.py#L1572

veni-vidi-vici-dormivi commented 2 months ago

Shruti wrote this herself so she could have variable (or dependent) values.

Hm, I mean we could also pass every (value, lambda) pair to the sklearn power transform no? Like we do for the inverse transform.

$\lambda$ can be any real value and does not have to be in 0..2

In our case it is because lambda is derived from a logistic function.

veni-vidi-vici-dormivi commented 2 months ago

I don't get how it makes sense that in one case $\lambda == eps$ is not contained in the case and in the other it is?

mathause commented 2 months ago

Shruti wrote this herself so she could have variable (or dependent) values.

Hm, I mean we could also pass every (value, lambda) pair to the sklearn power transform no? Like we do for the inverse transform.

Yes, but then we have to check if this is vectorized - otherwise it will be too slow.

λ can be any real value and does not have to be in 0..2

In our case it is because lambda is derived from a logistic function.

Ah ok, sorry - it's too many open PRs and comments. But then this is a function of the covariate function and it's not optimal if this is in

https://github.com/MESMER-group/mesmer/blob/10c8c432ba21b98d75dd05fbfe70b7d7d30404d5/mesmer/mesmer_m/power_transformer.py#L219

(Technically the user will not be able to easily replace the lambda_function but it's misleading when the bounds are associated with the Yeo-Johnson transform and not with the covariate function. And almost impossible to find out if the lambda_function should ever be changed...)

veni-vidi-vici-dormivi commented 2 months ago

I agree, it is pretty hard to see through it all. The get_yeo_johnson_lambdas is pretty horrible. I think we could get rid of it when rewriting the whole thing for xarrays. Maybe we should think about passing the lambda function as argument... But in this PR I actually just wanted to fix the operator, to be the same as in sklearn. Am I missing something?

mathause commented 2 months ago

But in this PR I actually just wanted to fix the operator, to be the same as in sklearn. Am I missing something?

No - I was trying to understand why it's inconsistent in scikit-learn (and I think I do now).