Calculate RMSDs from CrystFEL 0.11 Millepede bin files

rmsd_from_millepede.py

import numpy as np
import sys

# Usage: cat indexing5/mille-run*/mille-data-*.bin | dials.python ~/scripts/rmsd_from_millepede.py | tee indexing5/RMSD
# TODO: Use online standard deviation calculation
# TODO: Make this more general

NPANELS = 9
residuals = [[] for x in range(NPANELS * 2)]

rawfloat = np.frombuffer(sys.stdin.buffer.read(), dtype=np.float32)
rawint = rawfloat.view(dtype=np.int32)

pos = 0

while pos != len(rawint):
    blocksize = rawint[pos]
    itemcount = (blocksize - 2) // 2

    pos += 2 # size and 0
    grads = rawfloat[pos:(pos + itemcount)]
    pos += itemcount + 1 # float array and 0
    ids = rawint[pos:(pos + itemcount)]
    pos += itemcount # int array

    zeros = np.nonzero(ids == 0)[0]
    assert len(zeros) % 6 == 0
#    print(blocksize, pos, zeros)

    for idx in range(len(zeros) // 6):
        residual_x = grads[zeros[idx * 6]] 
        residual_y = grads[zeros[idx * 6 + 2]]
        panel = ids[zeros[idx * 6 + 1] + 1] // 100
        residuals[panel * 2].append(residual_x)
        residuals[panel * 2 + 1].append(residual_y)
        residuals[0].append(residual_x) # global
        residuals[1].append(residual_y)

for i, errors in enumerate(residuals):
     panelname = "Panel %d" % (i // 2)
     if (i // 2) == 0:
          panelname = "Global"

     direction = "fast scan"
     if i % 2 == 1:
         direction = "slow scan"

     rmsd = np.std(errors)

     print("%-7s %s: %.3f px %-8d spots" % (panelname, direction, rmsd, len(errors)))