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

[DHekstra]

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]

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.

❗ 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 ⚠️

Additional details and impacted files

@@            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	Coverage Δ
unittests	92.36% <66.66%> (-0.04%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[kmdalton]

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]

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

[kmdalton]

Sigma

[DHekstra]

cool. makes sense.

[DHekstra]

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

[DHekstra]

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]

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

[DHekstra]

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]

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.

[JBGreisman]

This looks fine to me -- I agree that it makes sense to have the default behavior as unchanged relative to before, and the proposed call to np.clip seems fine for getting this functionality.