Use only cc_half by default in resolution analysis

[rjgildea]

CC1/2 is now accepted as the most appropriate for estimating the resolution
cutoff (other than paired refinement). I/sigI-based cutoffs can be unreliable
due to dependence of the sigmas on error modelling in integration/scaling.
This is now consistent with xia2 default behaviour.

See also xia2/xia2#374

[codecov]

Codecov Report

Merging #1492 (bd2377a) into master (a583e91) will not change coverage.
The diff coverage is 100.00%.

@@           Coverage Diff           @@
##           master    #1492   +/-   ##
=======================================
  Coverage   65.57%   65.57%           
=======================================
  Files         614      614           
  Lines       68956    68956           
  Branches     9519     9519           
=======================================
  Hits        45219    45219           
  Misses      21898    21898           
  Partials     1839     1839           

[graeme-winter]

Makes sense, though I wonder what will happen if $USER passes in a data set with multiplicity ~ 1? Overall though good that this enforces current best practice for consistency.

[graeme-winter]

possibly "change default to only use cc_half in resolution analysis"?

[rjgildea]

Makes sense, though I wonder what will happen if $USER passes in a data set with multiplicity ~ 1? Overall though good that this enforces current best practice for consistency.

It would only use those reflections with multiplicity > 1, obviously how many of these there are will depend on the symmetry etc. The cc_half_significance will also come into play here, as it will only use bins where the calculated CC1/2 is significant given the number of reflections in the bin. This could be a use case for cc_half_method=sigma_tau.

First scale a single narrow wedge of data:

$ dials.scale $(dials.data get -q multi_crystal_proteinase_k)/{experiments_1.json,reflections_1.pickle}

Run dials.estimate_resolution with defaults:

dials.estimate_resolution scaled.{expt,refl}

Using cc_half_method=sigma_tau:

dials.estimate_resolution scaled.{expt,refl} cc_half_method=sigma_tau

Using misigma=1 (I don't know how reliable the sigmas will be after scaling with such low multiplicity):

dials.estimate_resolution scaled.{expt,refl} misigma=1

[dagewa]

It makes sense for this to be consistent with xia2.

The issue of what cut off value to use has come up again recently in CCP4 land. 0.143 has been suggested as being better than 0.3, because a CC1/2 of 0.143 corresponds to CC* of 0.5. But having a CC* of 0.5 is also an arbitrary value. The right value of the cut off is out of scope here, but I thought it might be worth noting the discussion.

[jbeilstenedmands]

Makes sense to me too.

To test a really low-multiplicity dataset, you can do something like this:
$ dials.scale $(dials.data get -q l_cysteine_dials_output)/{11_integrated.expt,11_integrated.refl,23_integrated.refl,23_integrated.expt}
$ dials.estimate_resolution scaled.{expt,refl}

N.B. This fails both on master and this branch with a stacktrace, but this shouldn't be a blocker for this work:

    merging_stats, metric, model, limit, sel=sel
  File "/dials/modules/dials/util/resolution_analysis.py", line 143, in resolution_fit_from_merging_stats
    return resolution_fit(d_star_sq, y_obs, model, limit, sel=sel)
  File "/Users/whi10850/dials/modules/dials/util/resolution_analysis.py", line 178, in resolution_fit
    raise RuntimeError("No reflections left for fitting")
RuntimeError: Please report this error to dials-support@lists.sourceforge.net: No reflections left for fitting```

[rjgildea]

The issue here is that even with just 4 bins, there is only 1 reflection in each of 2 bins with a multiplicity > 1:

  Statistics by resolution bin:
 d_max  d_min   #obs  #uniq   mult.  %comp       <I>  <I/sI>    r_mrg   r_meas    r_pim   r_anom   cc1/2   cc_ano
  1.94   1.01     11     10    1.10   3.36    3856.6    39.3    0.004    0.005    0.004    0.007   0.000   0.000
  1.01   0.77      9      9    1.00   2.47    1965.0    24.0    0.000    0.000    0.000     0.000   0.000
  0.77   0.64     11     10    1.10   2.06    1170.4    19.1    0.130    0.184    0.130    0.251   0.000   0.000
  0.64   0.62      9      9    1.00   6.12     502.5    10.1    0.000    0.000    0.000     0.000   0.000
  1.94   0.62     40     38    1.05   2.82    1907.3    23.5    0.032    0.045    0.032    0.063   0.995   0.000

Even the sigma-tau CC1/2 method requires > 1 unique reflection with multiplicity > 1 in order to calculate a non-zero CC1/2 value:
https://github.com/cctbx/cctbx_project/blob/5ae04d2955fa897b34ed4ec5bd8aea427aeeeafc/cctbx/miller/__init__.py#L5127-L5158

Granted, we should probably gracefully handle the resulting error, but as you say, that is outside the context of the PR.