rs-station / reciprocalspaceship

Tools for exploring reciprocal space
https://rs-station.github.io/reciprocalspaceship/
MIT License
29 stars 13 forks source link

Support clipping Sigma to avoid negative values in French-Wilson method #237

Closed DHekstra closed 10 months ago

DHekstra commented 11 months ago

I added a flag permitting clipping of Iobs to positive values for the purpose of calculating Sigma during anisotropic FW calculation. The option is off by default.

codecov-commenter commented 11 months ago

Codecov Report

Attention: 1 lines in your changes are missing coverage. Please review.

Comparison is base (0e2f5e2) 92.39% compared to head (a34429a) 92.36%. Report is 7 commits behind head on main.

:exclamation: Current head a34429a differs from pull request most recent head ca49796. Consider uploading reports for the commit ca49796 to get more accurate results

Files Patch % Lines
...alspaceship/algorithms/scale_merged_intensities.py 66.66% 1 Missing :warning:
Additional details and impacted files ```diff @@ Coverage Diff @@ ## main #237 +/- ## ========================================== - Coverage 92.39% 92.36% -0.04% ========================================== Files 37 37 Lines 2434 2436 +2 ========================================== + Hits 2249 2250 +1 - Misses 185 186 +1 ``` | [Flag](https://app.codecov.io/gh/rs-station/reciprocalspaceship/pull/237/flags?src=pr&el=flags&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=rs-station) | Coverage Δ | | |---|---|---| | [unittests](https://app.codecov.io/gh/rs-station/reciprocalspaceship/pull/237/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=rs-station) | `92.36% <66.66%> (-0.04%)` | :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=rs-station#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.

kmdalton commented 11 months ago

Wouldn't it make more sense to clip the returned values rather than the input intensities? This version of the code adds small errors to all miller indices rather than just the ones which would be assigned negative average intensities.

DHekstra commented 11 months ago

are you proposing clipping Sigma or the FW'ed intensities?

kmdalton commented 11 months ago

Sigma

DHekstra commented 11 months ago

cool. makes sense.

DHekstra commented 11 months ago

OK, this seems to do the trick. Clipping sigma to a user-specified small positive value (too small results in numerical warnings).

DHekstra commented 11 months ago
import reciprocalspaceship as rs
import numpy as np

ds=rs.read_mtz('pipeline/data/mtzs_origin/111.mtz')
ds_out=rs.algorithms.scale_merged_intensities(\
        ds, \
        intensity_key="IMEAN", \
        sigma_key="SIGIMEAN", \
        output_columns=["IP","SIGIP","FP","SIGFP"], \
        dropna=True, \
        inplace=False, \
        mean_intensity_method='anisotropic', \
        bw=2.0,
        clip_neg_Sigma=False,
        eps=0.0,
    )

produces

print(np.sum(np.isnan(ds_out["FP"].to_numpy())))

22238

ds_out=rs.algorithms.scale_merged_intensities(\ ds, \ intensity_key="IMEAN", \ sigma_key="SIGIMEAN", \ output_columns=["IP","SIGIP","FP","SIGFP"], \ dropna=True, \ inplace=False, \ mean_intensity_method='anisotropic', \ bw=2.0, clip_neg_Sigma=True, eps=0.1, )

produces

print(np.sum(np.isnan(ds_out["FP"].to_numpy()))) 0

kmdalton commented 11 months ago

I think it's better to just have a single new parameter, minimum_sigma=-np.inf

DHekstra commented 11 months ago

Let me know how this looks. I switched to using a single flag minimum_sigma as @kmdalton suggested, with a default value of -np.inf (no minimum). Because some of the ~440 datasets I tested this on also yielded some negative Sigma with the isotropic method, I moved the clipping into scale_merged_intensities() itself. I am applying the clipping before multiplication by the multiplicity as that strikes me as most appropriate. Let me know how this looks.

kmdalton commented 11 months ago

I am happy with this change. There might be a better parameter name (mea culpa), but otherwise this is the least invasive patch I can think of. Let's ask @JBGreisman for comments before merging, but I'm plenty satisfied.