Repository : ssh://g18-sc-serv-04.diamond.ac.uk/dials_scratch
On branch : master
commit ca62fda314b554c4528d2491cb9ca019fb644e51
Author: James Parkhurst <[log in to unmask]>
Date: Thu Jul 12 13:31:31 2018 +0100
Added models for angular and simple mosaicity models
ca62fda314b554c4528d2491cb9ca019fb644e51
jmp/potato/model.py | 404 ++++++++++++++++++++++++++++++++++++++++++++--------
1 file changed, 342 insertions(+), 62 deletions(-)
diff --git a/jmp/potato/model.py b/jmp/potato/model.py
index 8f5d397..f528d1f 100644
--- a/jmp/potato/model.py
+++ b/jmp/potato/model.py
@@ -1,26 +1,164 @@
from __future__ import division
import scitbx.linalg
from dials_scratch.jmp.potato.parameterisation import SimpleMosaicityParameterisation
+from dials_scratch.jmp.potato.parameterisation import AngularMosaicityParameterisation
+from scitbx.linalg import eigensystem, l_l_transpose_cholesky_decomposition_in_place
+from dials_scratch.jmp.potato import PredictorAngular
+from dials_scratch.jmp.potato import PredictorSimple
+from dials_scratch.jmp.potato import BBoxCalculatorAngular
+from dials_scratch.jmp.potato import BBoxCalculatorSimple
+from dials_scratch.jmp.potato import MaskCalculatorAngular
+from dials_scratch.jmp.potato import MaskCalculatorSimple
from dials.array_family import flex
from scitbx import matrix
+from math import exp
-class SimpleMosaicityModel(object):
+class Simple6ProfileModel(object):
'''
- A class to hold a simple 6 parameter mosaicity model
+ Class to store profile model
'''
- def __init__(self, sigma):
+ def __init__(self, params):
'''
- Initialise the with mosaicity matrix
+ Initialise the class
'''
- self._sigma = sigma
+ self.params = params
def parameterisation(self):
'''
- Return the parameterisation
+ Get the parameterisation
+
+ '''
+ return SimpleMosaicityParameterisation(self.params)
+
+ def sigma(self):
+ '''
+ Get the sigma
+
+ '''
+ return self.parameterisation().sigma()
+
+ def update_model_state_parameters(self, state):
+ '''
+ Update the model state with the parameters
+
+ '''
+ state.set_M_params(self.params)
+
+ def update_model(self, state):
+ '''
+ Update the model
+
+ '''
+
+ # Compute the eigen decomposition of the covariance matrix and check
+ # largest eigen value
+ eigen_decomposition = eigensystem.real_symmetric(
+ state.get_M().as_flex_double_matrix())
+ L = eigen_decomposition.values()
+ if L[0] > 1e-5:
+ raise RuntimeError("Mosaicity matrix is unphysically large")
+
+ self.params = state.get_M_params()
+
+ def predict_reflections(self,
+ experiments,
+ miller_indices,
+ probability=0.997):
+ '''
+ Predict the reflections
+
+ '''
+ predictor = PredictorSimple(
+ experiments[0],
+ self.sigma(),
+ probability)
+ return predictor.predict(miller_indices)
+
+ def compute_bbox(self, experiments, reflections):
+ '''
+ Compute the bounding box
+
+ '''
+ calculator = BBoxCalculatorSimple(
+ experiments[0],
+ self.sigma(),
+ 0.999999998,
+ 4)
+ calculator.compute(reflections)
+
+ def compute_mask(self, experiments, reflections):
+ '''
+ Compute the mask
+
+ '''
+ calculator = MaskCalculatorSimple(
+ experiments[0],
+ self.sigma(),
+ 0.999999998)
+ calculator.compute(reflections)
+
+ def sigma_for_reflection(self, s0, r):
+ '''
+ Get sigma for a reflections
+
+ '''
+ return self.sigma()
+
+ def compute_partiality(self, experiments, reflections):
+ '''
+ Compute the partiality
+
+ '''
+ s0 = matrix.col(experiments[0].beam.get_s0())
+ partiality = flex.double(len(reflections))
+ for k in range(len(reflections)):
+ s1 = matrix.col(reflections[k]['s1'])
+ s2 = matrix.col(reflections[k]['s2'])
+ sbox = reflections[k]['shoebox']
+
+ r = s2 - s0
+ sigma = experiments[0].crystal.mosaicity.sigma()
+ R = compute_change_of_basis_operation(s0, s2)
+ S = R*(sigma)*R.transpose()
+ mu = R*s2
+ assert(abs(1-mu.normalize().dot(matrix.col((0,0,1)))) < 1e-7)
+
+ S11 = matrix.sqr((
+ S[0], S[1],
+ S[3], S[4]))
+ S12 = matrix.col((S[2], S[5]))
+ S21 = matrix.col((S[6], S[7])).transpose()
+ S22 = S[8]
+
+ mu1 = matrix.col((mu[0], mu[1]))
+ mu2 = mu[2]
+ partiality[k] = exp(-0.5*(s0.length()-mu2) * (1/S22) * (s0.length()-mu2))
+ reflections['partiality'] = partiality
+
+ @classmethod
+ def from_sigma_d(Class, sigma_d):
+ '''
+ Create the profile model from sigma_d estimate
+
+ '''
+ return Class.from_params(flex.double((sigma_d, 0, sigma_d, 0, 0, sigma_d)))
+
+ @classmethod
+ def from_params(Class, params):
+ '''
+ Create the class from some parameters
+
+ '''
+ return Class(params)
+
+ @classmethod
+ def from_sigma(Class, sigma):
+ '''
+ Construct the profile model from the sigma
'''
@@ -28,33 +166,208 @@ class SimpleMosaicityModel(object):
LL = flex.double()
for j in range(3):
for i in range(j+1):
- LL.append(self._sigma[j*3+i])
+ LL.append(sigma[j*3+i])
# Do the cholesky decomposition
- ll = scitbx.linalg.l_l_transpose_cholesky_decomposition_in_place(LL)
+ ll = l_l_transpose_cholesky_decomposition_in_place(LL)
# Setup the parameters
- params = (
+ return Class.from_params(flex.double((
LL[0],
LL[1], LL[2],
- LL[3], LL[4], LL[5])
+ LL[3], LL[4], LL[5])))
+
- # Return the parameterisation
- return SimpleMosaicityParameterisation(params)
- def compose(self, parameterisation):
+class Angular2ProfileModel(object):
+ '''
+ Class to store profile model
+
+ '''
+
+ def __init__(self, params):
'''
- Compose the model from the parameterisation
+ Initialise the class
'''
- self._sigma = parameterisation.sigma()
+ self.params = params
+
+ def parameterisation(self):
+ '''
+ Get the parameterisation
+
+ '''
+ return AngularMosaicityParameterisation(self.params)
def sigma(self):
'''
- Get the covariance matrix
+ Get the sigma
+
+ '''
+ return self.parameterisation().sigma()
+
+ def update_model_state_parameters(self, state):
+ '''
+ Update the model state with the parameters
+
+ '''
+ state.set_M_params(self.params)
+
+ def update_model(self, state):
+ '''
+ Update the model
+
+ '''
+
+ # Compute the eigen decomposition of the covariance matrix and check
+ # largest eigen value
+ eigen_decomposition = eigensystem.real_symmetric(
+ state.get_M().as_flex_double_matrix())
+ L = eigen_decomposition.values()
+ if L[0] > 1e-5:
+ raise RuntimeError("Mosaicity matrix is unphysically large")
+
+ self.params = state.get_M_params()
+
+ def sigma_for_reflection(self, s0, r):
+ '''
+ Get sigma for a reflections
+
+ '''
+ Q = compute_change_of_basis_operation(s0, r)
+ return Q.transpose()*self.sigma()*Q
+
+ def predict_reflections(self,
+ experiments,
+ miller_indices,
+ probability=0.997):
+ '''
+ Predict the reflections
+
+ '''
+ predictor = PredictorAngular(
+ experiments[0],
+ self.sigma(),
+ probability)
+ return predictor.predict(miller_indices)
+
+ def compute_bbox(self, experiments, reflections):
+ '''
+ Compute the bounding box
+
+ '''
+ calculator = BBoxCalculatorAngular(
+ experiments[0],
+ self.sigma(),
+ 0.999999998,
+ 4)
+ calculator.compute(reflections)
+
+ def compute_mask(self, experiments, reflections):
+ '''
+ Compute the mask
+
+ '''
+ calculator = MaskCalculatorAngular(
+ experiments[0],
+ self.sigma(),
+ 0.999999998)
+ calculator.compute(reflections)
+
+ def compute_partiality(self, experiments, reflections):
+ '''
+ Compute the partiality
+
+ '''
+ s0 = matrix.col(experiments[0].beam.get_s0())
+ partiality = flex.double(len(reflections))
+ for k in range(len(reflections)):
+ s1 = matrix.col(reflections[k]['s1'])
+ s2 = matrix.col(reflections[k]['s2'])
+ sbox = reflections[k]['shoebox']
+
+ r = s2 - s0
+ sigma = experiments[0].crystal.mosaicity.sigma()
+ R = compute_change_of_basis_operation(s0, s2)
+ Q = compute_change_of_basis_operation(s0, r)
+ S = R*(Q.transpose()*sigma*Q)*R.transpose()
+ mu = R*s2
+ assert(abs(1-mu.normalize().dot(matrix.col((0,0,1)))) < 1e-7)
+
+ S11 = matrix.sqr((
+ S[0], S[1],
+ S[3], S[4]))
+ S12 = matrix.col((S[2], S[5]))
+ S21 = matrix.col((S[6], S[7])).transpose()
+ S22 = S[8]
+
+ mu1 = matrix.col((mu[0], mu[1]))
+ mu2 = mu[2]
+ partiality[k] = exp(-0.5*(s0.length()-mu2) * (1/S22) * (s0.length()-mu2))
+ reflections['partiality'] = partiality
+
+
+ @classmethod
+ def from_sigma_d(Class, sigma_d):
+ '''
+ Create the profile model from sigma_d estimate
+
+ '''
+ return Class.from_params(flex.double((sigma_d, sigma_d)))
+
+ @classmethod
+ def from_params(Class, params):
+ '''
+ Create the class from some parameters
+
+ '''
+ return Class(params)
+
+ @classmethod
+ def from_sigma(Class, sigma):
+ '''
+ Construct the profile model from the sigma
+
+ '''
+
+ # Construct triangular matrix
+ LL = flex.double()
+ for j in range(3):
+ for i in range(j+1):
+ LL.append(sigma[j*3+i])
+
+ # Do the cholesky decomposition
+ ll = l_l_transpose_cholesky_decomposition_in_place(LL)
+
+ # Check the sigma is as we expect
+ TINY = 1e-10
+ assert abs(LL[1] - 0) < TINY
+ assert abs(LL[2] - LL[0]) < TINY
+ assert abs(LL[3] - 0) < TINY
+ assert abs(LL[4] - 0) < TINY
+
+ # Setup the parameters
+ return Class.from_params(flex.double((LL[0], LL[5])))
+
+
+class ProfileModelFactory(object):
+ '''
+ Class to create profile models
+
+ '''
+ @classmethod
+ def from_sigma_d(Class, params, sigma_d):
+ '''
+ Construct a profile model from an initial sigma estimate
'''
- return self._sigma
+ if params.profile.rlp_mosaicity.model == "simple6":
+ return Simple6ProfileModel.from_sigma_d(sigma_d)
+ elif params.profile.rlp_mosaicity.model == "angular2":
+ return Angular2ProfileModel.from_sigma_d(sigma_d)
+
+ raise RuntimeError("Unknown profile model: %s" %
+ params.profile.rlp_mosaicity.model)
def compute_change_of_basis_operation(s0, s2):
@@ -71,51 +384,18 @@ def compute_change_of_basis_operation(s0, s2):
e3.elems)
return R
+def compute_change_of_basis_operation2(s0, r):
+ '''
+ Compute the change of basis operation that puts the r vector along the z axis
-# class ReflectionModel(object):
-
-# def __init__(self, model, s2):
-
-# # Save some stuff
-# self.model = model
-# self.s0 = s0
-# self.s2 = s2
-
-# # Compute the change of basis
-# self.R = compute_change_of_basis_operation(s0, s2)
-
-# # Rotate the covariance matrix
-# self.S = self.R*self.model.sigma()*self.R.transpose()
-
-# # Construct the marginal distribution
-# self._marginal = MarginalDistribution(self.S)
-
-# # Construct the conditional distribution
-# self._conditional = ConditionalDistribution(self.S, s0.length(), s2.length())
-
-# def marginal(self):
-# '''
-# Return the marginal
-
-# '''
-# return self._marginal
-
-# def conditional(self):
-# '''
-# Return the conditional
-
-# '''
-# return self._conditional
-
-# def s1(self):
-# '''
-# Compute the diffracted beam vector
+ '''
+ e1 = r.cross(s0).normalize()
+ e2 = r.cross(e1).normalize()
+ e3 = r.normalize()
+ R = matrix.sqr(
+ e1.elems +
+ e2.elems +
+ e3.elems)
+ return R
-# '''
-# mubar = self._conditional.mean()
-# v = matrix.col((
-# mubar[0],
-# mubar[1],
-# self.s0.length())).normalize() * self.s0.length()
-# return self.R.transpose() * v