Calculate R_d for XDS unmerged reflections

[DHekstra]

We received a request to implement R_d for unmerged integrated intensities before/after correction using XSCALE as described in
Diederichs, Acta D, 2005.

[kmdalton]

This looks relatively straight forward. I'd like to put the actual calculation in rs.algorithms and include a test.

[DHekstra]

from Tom Caradoc Davies:

Apologies for the delay in getting back to you. We had our first users on the beamline this week and it was a bit hectic.

I have uploaded two data files to google drive here:
https://drive.google.com/drive/folders/15H03HNwcSfxk74cFS-XQ72tmkyO6761_?usp=drive_link

There are two datasets from the same insulin crystal. Both processed the same way with XDS, one with Friedel’s law True and one False.

I can also send the unmerged reflection files (INTEGRATE.HKL) files if that is more useful?

The dataset is 720 degrees and processing the first 180 vrs the last there is ~50% drop in signal in the highest resolution bin. Hopefully this is enough radiation damage to show up on an Rd plot.

If not let me know and I can roast some crystals for you!

[DHekstra]

(I requested access to that GDrive folder and will add files when I've got them)

[DHekstra]

import argparse
import pandas as pd
import numpy as np
import reciprocalspaceship as rs
from rsbooster.stats.parser import BaseParser
from careless.io import xds
from matplotlib import pyplot as plt
from tqdm import tqdm

plt.rcParams.update({'font.size': 16})


class ArgumentParser(BaseParser):
    def __init__(self):
        super().__init__(
            description=__doc__
        )
        
        # Required arguments
        self.add_argument(
            "hkl",
            help="HKL containing unmerged, scaled reflections from XDS",
        )
        self.add_argument(
            "--ranges",
            nargs="+",
            default=None,
            help="Resolution ranges specified as 'low,high' (e.g. '100,3'). Multiple ranges allowed."
        )
        self.add_argument(
            "--n-samples",
            type=int,
            default=10000,
            help="Number of samples to process (default: 10000)"
        )
        self.add_argument(
            "--n-dist-bins",
            type=int,
            default=100,
            help="Number of distribution bins (default: 100)"
        )
        self.add_argument(
            "--anomalous",
            action="store_true",
            default=False,
            help="If passed, keep Friedel mates separate"
        )


def parse_resolution_ranges(range_strings):
    """
    Converts list of strings like ['10,20', '30,40'] into list of tuples [(10, 20), (30, 40)].
    """
    ranges = []
    for r in range_strings:
        try:
            low, high = map(float, r.split(','))
            ranges.append((low, high))
        except ValueError:
            raise argparse.ArgumentTypeError(f"Invalid range format: '{r}'. Expected format 'low,high'.")
    return ranges


def calculate_Rd_per_HKL(group, key):
    # Implement Diederichs, Acta Cryst. (2006). D62, 96–101, p. 98 top.
    result=[]
    for i in range(len(group.index)):
        for j in range(len(group.index)):
           if i > j:
                d_ij = np.abs(group.iloc[i]["BATCH"] - group.iloc[j]["BATCH"])
                DI_ij= np.abs(group.iloc[i]["IOBS"]  - group.iloc[j]["IOBS"])
                MI_ij= np.abs(group.iloc[i]["IOBS"]  + group.iloc[j]["IOBS"])/2.0
                result.append([key, d_ij, DI_ij,MI_ij])
           else:
                pass
    return rs.DataSet(result, columns=["key","d_ij","DI_ij","MI_ij"])  # Return as a Series or DataFrame

def moving_average(data, window_size):
    return np.convolve(data, np.ones(window_size)/window_size, mode='valid')

def run_analysis(args):
    print("Input file:             ", args.hkl)
    print("Resolution ranges:      ", args.ranges)
    print("Number of samples:      ", args.n_samples)
    print("Number of bins:         ", args.n_dist_bins)
    print("Keep Friedel separate:  ", args.anomalous)    
    parsed_ranges = parse_resolution_ranges(args.ranges) if args.ranges else None

    if parsed_ranges is None:
        parsed_ranges = [(1000,0.01)]

    results = []
    labels = None
    for res_range in parsed_ranges:
        print(f"Processing range {res_range} A:")
        ds=xds.read_hkl(args.hkl)
        ds.set_index(['H', 'K', 'L'], inplace=True)
        ds.hkl_to_asu(inplace=True, anomalous=args.anomalous)
        ds.compute_dHKL(inplace=True)

        ds=ds.loc[(ds.dHKL>=res_range[1]) & (ds.dHKL<=res_range[0]),:]
        ds=ds.sample(args.n_samples)

        results = []
        for key, group in tqdm(ds.loc[:,["BATCH","IOBS"]].groupby(ds.index), desc="Processing Miller indices"):
            result = calculate_Rd_per_HKL(group, key)
            results.append(result)
            
        result_df=pd.concat(results)
        result_df["d_ij_binned"]=pd.qcut(result_df["d_ij"].astype(float), args.n_dist_bins)
        
        result_df_mean=result_df.groupby("d_ij_binned").mean()
        result_df_mean.reset_index(inplace=True)
        
        x=result_df_mean.d_ij
        y=result_df_mean.DI_ij/result_df_mean.MI_ij
        plt.plot(x,y,'o',alpha=0.3,label=res_range)
        
        z = np.polyfit(x,y,3)  # Degree 1 for linear trendline
        p = np.poly1d(z)
        plt.plot(x, p(x), label='Trendline')

        # Smooth the data
        window_size = np.ceil(args.n_dist_bins/20).astype(int)  # Adjust as needed for more or less smoothing
        y_smooth = moving_average(y, window_size)
        x_smooth = x[(window_size - 1) // 2 : -(window_size // 2)]
        plt.plot(x_smooth, y_smooth, label=f'Moving Average (window={window_size})')
        
    plt.xlabel("Frame (discretized)")
    plt.ylabel(r"$R_d$")
    plt.ylim([0,1.1*plt.ylim()[1]])
    plt.grid()
    plt.show()
    
    if args.image is not None:
        plt.savefig(args.image)

    if args.show:
        plt.show()


def parse_arguments():
    return ArgumentParser()

def main():
    run_analysis(parse_arguments().parse_args())

if __name__ == "__main__":
    main()
'''

[kmdalton]

@DHekstra -- Is this completed? There's no PR referenced. Should we work on putting your script in the stats submodule?

[DHekstra]

Good question. This is specific to XDS output. I'm happy to start a PR to put this in the stats module rather than have this be a one-off script.

[DHekstra]

It's a little slow (~1 min) and could benefit from multiprocessing but otherwise fully functional.

[kmdalton]

I see several place with for loops that can probably be vectorized. There's most likely no need to parallelize.