violini in workflow

.github/workflows/CI.yml

@@ -38,6 +38,8 @@ jobs:
         run: |
           unzip tests/data/PbTe-CNCS-HH4-Ei_12-E_5-q_0/dhkls.npy.zip -d tests/data/PbTe-CNCS-HH4-Ei_12-E_5-q_0
           ls -tl tests/data/beam/ARCS/100meV/out/neutrons
+          mkdir tests/data/SEQUOIA_data/Mn3Si2Te6
+          tar xfz tests/data/SEQUOIA_data/Mn3Si2Te6-mcvine-res-sim.tgz --directory tests/data/SEQUOIA_data/Mn3Si2Te6/
 
       - name: ~/.mantid
         run: git clone https://github.com/yxqd/dotmantid ~/.mantid

dgsres/instruments.py

@@ -24,6 +24,7 @@ class sequoia:
         name = 'SEQ',
         detsys_radius = "5.5*meter",
         L_m2s = "20.05*meter",
+        L_m2fc = "18*meter",
         offset_sample2beam = "-0.15*meter" # offset from sample to saved beam. don't change this unless you are sure what you are doing
     )
 
@@ -37,6 +38,13 @@ class sequoia:
     def scattering_angle_constraints(cls, theta, phi):
         return ((theta<60.) * (theta>-30)) * (phi<18) * (phi>-18)
 
+    class violini:
+
+        tau_P = 10
+        tau_M = 8
+        sigma_thetai = 0.01
+        sigma_phii = 0.01
+
 class cncs:
     instrument = sx.instrument(
         name = 'CNCS',

dgsres/singlextal/plot.py

@@ -1,6 +1,12 @@
 import os, numpy as np
 from . import use_covmat
+from .pointcloud import PointCloud
+from .simdata import McvineResolutionData
+from .violini import VioliniModel
+# backward compatible
+AnalyticalModel = VioliniModel
 
+# obsolete
 def createARCSAnalyticalModel(
     tau_P, tau_M, pix_r, pix_h,
     sample_thickness, 
@@ -24,147 +30,22 @@ def createARCSAnalyticalModel(
         height = "meter*%s" % pix_h,
         pressure = "10*atm",
         )
-    return AnalyticalModel(
+    return VioliniModel(
         instrument, pixel, tofwidths, beamdivs, 
         sample_yml, samplethickness,
         Ei, psi_scan)
 
 
-class PointCloud:
-
-    """resolution point cloud. can be regarded as "events" with dh,dk,dl,dE,weight data.
-    These data could be generated from mcvine simulation of the resolution function,
-    or generated from a simple normal distribution using a 4D cov matrix (see AnalyticalModel).
-    """
-
-    def __init__(self, dhs, dks, dls, dEs, weights):
-        self.dhs = dhs
-        self.dks = dks
-        self.dls = dls
-        self.dEs = dEs
-        self.weights = weights
-        return
-
-    def getThinSlice(
-        self, 
-        axis1=('h', -0.1,0.1, 0.002), axis2=('E', -3, 3., 0.04),
-        axis3=('k', -0.006, 0.006), axis4=('l', -0.006, 0.006)
-    ):
-        axis1_name = axis1[0]; axis1_ticks = np.arange(*axis1[1:])
-        axis2_name = axis2[0]; axis2_ticks = np.arange(*axis2[1:])
-        condition = True
-        for ax in [axis3, axis4]:
-            name = ax[0]
-            min, max = ax[1:]
-            arr = self._evts(name)
-            condition *= (arr<max) * (arr>min)
-            continue
-        Ixy, xedges, yedges = np.histogram2d(
-            self._evts(axis1_name)[condition], self._evts(axis2_name)[condition],
-            bins=(axis1_ticks, axis2_ticks), weights=self.weights[condition]) 
-        xbc = (xedges[:-1] + xedges[1:])/2
-        ybc = (yedges[:-1] + yedges[1:])/2
-        xg, yg = np.meshgrid(xbc, ybc)
-        return xg,yg,Ixy
-
-    def _evts(self, name):
-        return getattr(self, 'd%ss' % name)
-
-
-class McvineResolutionData:
-
-    """resolution data simulated by mcvine. 
-    The data was simulated earlier and saved to disk. 
-    This class handles locating the data files and reading the data, 
-    and convert the data to different forms: point cloud, cov matrix, etc.
-    """
-
-    def __init__(self, parent_dir, dirname_template='E%s_hkl%s'):
-        self.parent_dir = parent_dir
-        self.dirname_template = dirname_template
-        return
-
-    def path(self, hkl, E):
-        return os.path.join(self.parent_dir, self.dirname_template % (E, '%s,%s,%s' % tuple(hkl)))
-
-    def loadData(self, hkl, E):
-        p = self.path(hkl, E)
-        return self._loadData(p)
-
-    def loadPointCloud(self, hkl, E):
-        dhs, dks, dls, dEs, probs = self.loadData(hkl, E)
-        return PointCloud(dhs, dks, dls, dEs, probs)
-
-    def _loadData(self, outdir1):
-        dhkls = np.load('%s/dhkls.npy' % outdir1)
-        dEs = np.load('%s/dEs.npy' % outdir1)
-        probs = np.load('%s/probs.npy' % outdir1)
-        dhs,dks,dls = dhkls.T
-        # there might be unreasonable data points with unreasonable weights
-        mask = (dhs> -2.)*(dhs<2.) \
-            * (dks> -2.)*(dks<2.) \
-            * (dls> -2.)*(dls<2.) 
-        return np.array([dhs[mask], dks[mask], dls[mask], dEs[mask], probs[mask]])
-
-    def computeCovMat(self, hkl, E):
-        data = self.loadData(hkl, E)
-        Data = data[:4]; probs = data[-1]
-        return np.cov(Data, aweights=probs)
-
-class AnalyticalModel:
-
-    """analytical resolution model based on paper by Violini et al.
-    The main calculation is done in module .use_covmat.
-    """
-
-    def __init__(
-        self, instrument, pixel, tofwidths, beamdivs, 
-        sample_yml, samplethickness,
-        Ei, psi_scan):
-        self.instrument = instrument
-        self.pixel = pixel
-        self.tofwidths = tofwidths
-        self.beamdivs = beamdivs
-        self.sample_yml = sample_yml
-        self.samplethickness = samplethickness
-        self.Ei = Ei
-        self.psi_scan = psi_scan
-        return
-
-    def computePointCloud(self, hkl, E, N=int(1e6)):
-        covmat = self.computeCovMat(hkl, E)
-        events = np.random.multivariate_normal(np.zeros(4), covmat, size=N)
-        dhs, dks, dls, dEs = events.T
-        ws = np.ones(dhs.shape)
-        return PointCloud(dhs, dks, dls, dEs, ws)
-
-    def computeCovMat(self, hkl, E):
-        class dynamics:
-            hkl_dir = np.array([1.,0.,0.])
-            dq = 0
-        dynamics.hkl0 = hkl
-        dynamics.E = E
-        cm_res = use_covmat.compute(
-            self.sample_yml, self.Ei, 
-            dynamics, 
-            self.psi_scan,
-            self.instrument, self.pixel,
-            self.tofwidths, self.beamdivs, self.samplethickness,
-            plot=False)
-        # ellipsoid_trace = cm_res['u']
-        InvCov4D = cm_res['hklE_inv_cov']
-        return np.linalg.inv(InvCov4D)/2.355
-
-def plotEllipsoid(covmat, q, symbol='.'):
+def plot_qE_ellipse(covmat, q, symbol='.', **kwds):
     """plot 2d ellipsoid along a particular hkl direction, given the cov matrix.
     """
     invcm = np.linalg.inv(covmat)
-    _, u = computeEllipsoid(invcm, q)
+    _, u = compute_qE_ellipse(invcm, q)
     from matplotlib import pyplot as plt
-    plt.plot(u[:,0], u[:,1], symbol)
+    plt.plot(u[:,0], u[:,1], symbol, **kwds)
     return
 
-def computeEllipsoid(InvCov4D, q):
+def compute_qE_ellipse(InvCov4D, q):
     "compute ellipsoid along a q direction"
     qE2qE = np.array(
         [np.hstack([q, [0]]),
@@ -181,8 +62,36 @@ def computeEllipsoid(InvCov4D, q):
     u = np.dot(up, mR.T)
     return inv_cov_qE, u
 
-def computeEllipsoids(InvCov4D, directions=None):
+def compute_qE_ellipses(InvCov4D, directions=None):
     if directions is None:
         directions = np.eye(3, dtype=float)
-    return [(q, computeEllipsoid(InvCov4D, q)) for q in directions]
+    return [(q, compute_qE_ellipse(InvCov4D, q)) for q in directions]
+
+def plot_qq_ellipse(covmat, q1, q2, symbol='.', **kwds):
+    """plot 2d ellipsoid along two hkl directions, given the cov matrix.
+    """
+    invcm = np.linalg.inv(covmat)
+    _, u = compute_qq_ellipse(invcm, q1, q2)
+    from matplotlib import pyplot as plt
+    plt.plot(u[:,0], u[:,1], symbol, **kwds)
+    return
+
+def compute_qq_ellipse(InvCov4D, q1, q2):
+    "compute ellipsoid along a q direction"
+    q1q2_to_qE = np.array([
+        np.hstack([q1, [0]]),
+        np.hstack([q2, [0]]),
+    ])
+    inv_cov = np.dot(q1q2_to_qE, np.dot(InvCov4D, q1q2_to_qE.T))
+    # print inv_cov_hE
+    r = np.linalg.eig(inv_cov)
+    mR = r[1]; lambdas = r[0]
+    RR = 2*np.log(2)
+    theta = np.arange(0, 360, 1.)*np.pi/180
+    u1p = np.sqrt(RR/lambdas[0])*np.cos(theta)
+    u2p = np.sqrt(RR/lambdas[1])*np.sin(theta)
+    up = np.array([u1p, u2p]).T
+    u = np.dot(up, mR.T)
+    return inv_cov, u
+
 

dgsres/singlextal/pointcloud.py

@@ -0,0 +1,53 @@
+import os, numpy as np
+
+class PointCloud:
+
+    """resolution point cloud. can be regarded as "events" with dh,dk,dl,dE,weight data.
+    These data could be generated from mcvine simulation of the resolution function,
+    or generated from a simple normal distribution using a 4D cov matrix (see AnalyticalModel).
+    """
+
+    def __init__(self, dhs, dks, dls, dEs, weights):
+        self.dhs = dhs
+        self.dks = dks
+        self.dls = dls
+        self.dEs = dEs
+        self.weights = weights
+        return
+
+    def getThinSlice(
+        self, 
+        axis1=('h', -0.1,0.1, 0.002), axis2=('E', -3, 3., 0.04),
+        axis3=('k', -0.006, 0.006), axis4=('l', -0.006, 0.006)
+    ):
+        axis1_name = axis1[0]; axis1_ticks = np.arange(*axis1[1:])
+        axis2_name = axis2[0]; axis2_ticks = np.arange(*axis2[1:])
+        condition = True
+        for ax in [axis3, axis4]:
+            name = ax[0]
+            min, max = ax[1:]
+            arr = self._evts(name)
+            condition *= (arr<max) * (arr>min)
+            continue
+        Ixy, xedges, yedges = np.histogram2d(
+            self._evts(axis1_name)[condition], self._evts(axis2_name)[condition],
+            bins=(axis1_ticks, axis2_ticks), weights=self.weights[condition]) 
+        xbc = (xedges[:-1] + xedges[1:])/2
+        ybc = (yedges[:-1] + yedges[1:])/2
+        xg, yg = np.meshgrid(xbc, ybc)
+        return xg,yg,Ixy
+
+    def getIE(self, Emin, Emax, dE):
+        bins = np.arange(Emin, Emax, dE)
+        I, edges = np.histogram(self.dEs, bins=bins, weights=self.weights)
+        return edges, I
+
+    def getIq(self, direction, qmin, qmax, dq):
+        "getIq('h', -1., 1., 0.02)"
+        bins = np.arange(qmin, qmax, dq)
+        I, edges = np.histogram(self._evts(direction), bins=bins, weights=self.weights)
+        return edges, I
+
+    def _evts(self, name):
+        return getattr(self, 'd%ss' % name)
+

dgsres/singlextal/simdata.py

@@ -0,0 +1,46 @@
+import os, numpy as np
+from .pointcloud import PointCloud
+
+class McvineResolutionData:
+
+    """resolution data simulated by mcvine. 
+    The data was simulated earlier and saved to disk. 
+    This class handles locating the data files and reading the data, 
+    and convert the data to different forms: point cloud, cov matrix, etc.
+    """
+
+    def __init__(self, parent_dir, dirname_template='E%s_hkl%s'):
+        self.parent_dir = parent_dir
+        self.dirname_template = dirname_template
+        return
+
+    def path(self, hkl, E):
+        return os.path.join(self.parent_dir, self.dirname_template % (E, '%s,%s,%s' % tuple(hkl)))
+
+    def loadData(self, hkl, E, dhkl_ranges=None):
+        p = self.path(hkl, E)
+        return self._loadData(p, dhkl_ranges=dhkl_ranges)
+
+    def loadPointCloud(self, hkl, E):
+        dhs, dks, dls, dEs, probs = self.loadData(hkl, E)
+        return PointCloud(dhs, dks, dls, dEs, probs)
+
+    def _loadData(self, outdir1, dhkl_ranges=None):
+        if dhkl_ranges is None:
+            dhkl_ranges = [ (-2., 2.) ] * 3
+        dhkls = np.load('%s/dhkls.npy' % outdir1)
+        dEs = np.load('%s/dEs.npy' % outdir1)
+        probs = np.load('%s/probs.npy' % outdir1)
+        dhs,dks,dls = dhkls.T
+        # there might be unreasonable data points with unreasonable weights
+        (dh_min, dh_max), (dk_min, dk_max), (dl_min, dl_max) = dhkl_ranges
+        mask = (dhs>dh_min)*(dhs<dh_max) \
+            * (dks>dk_min)*(dks<dk_max) \
+            * (dls>dl_min)*(dls<dl_max) 
+        return np.array([dhs[mask], dks[mask], dls[mask], dEs[mask], probs[mask]])
+
+    def computeCovMat(self, hkl, E, dhkl_ranges=None):
+        data = self.loadData(hkl, E, dhkl_ranges=dhkl_ranges)
+        Data = data[:4]; probs = data[-1]
+        return np.cov(Data, aweights=probs)
+

dgsres/singlextal/spinw.py

@@ -41,5 +41,5 @@ def wrap_spinw_disp_func(f):
     import matlab.engine, matlab
     def _(start, end, Nq_disp):
         sp = f(matlab.double(list(start)), matlab.double(list(end)), Nq_disp)
-        return sp['omega'], sp['swInt']
+        return np.array(sp['omega']), np.array(sp['swInt'])
     return _

dgsres/singlextal/use_covmat.py

@@ -32,8 +32,9 @@ def __init__(self, theta, phi):
         return
 
 def compute(
-        sample_yml, Ei, dynamics, psi_scan, instrument, pixel, tofwidths, beamdivs, samplethickness,
-        plot=False):
+        sample_yml, Ei, dynamics, psi_scan, instrument, pixel,
+        tofwidths, beamdivs, samplethickness,
+        plot=False, verbose=False):
     # XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
     # should P be the T0 chopper?
     L_PM=mcvine.units.parse(instrument.L_m2fc)/mcvine.units.meter # P chopper to M chopper distance
@@ -70,8 +71,9 @@ def compute(
     from mcvine.workflow.singlextal.coords_transform import hkl2Q
     for r in results:
         xtalori.psi = r*np.pi/180
-        print("psi=%s, Q=%s" % (r, hkl2Q(hkl, xtalori)))
-        print("hkl2Q=%r\n(Q = hkl dot hkl2Q)" % (xtalori.hkl2cartesian_mat(),))
+        if verbose:
+            print("psi=%s, Q=%s" % (r, hkl2Q(hkl, xtalori)))
+            print("hkl2Q=%r\n(Q = hkl dot hkl2Q)" % (xtalori.hkl2cartesian_mat(),))
     # these are the psi angles that the particular point of interest will be measured
     # print results
     assert len(results)