Mccomposite

.gitignore

@@ -2,6 +2,7 @@
 *.dat
 *.mcv
 out.*
+log.*
 compiled_composite*
 composite_kernel_*
 
@@ -164,4 +165,4 @@ dmypy.json
 cython_debug/
 
 # PyCharm
-.idea/
+.idea/

acc/_numba.py

@@ -1,3 +1,4 @@
+import os, tempfile
 import numpy as np
 from numba.core import config
 if config.ENABLE_CUDASIM:
@@ -6,3 +7,38 @@ def xoroshiro128p_uniform_float32(rng_states, threadindex):
     from numba.cuda.random import xoroshiro128p_type
 else:
     from numba.cuda.random import xoroshiro128p_uniform_float32, xoroshiro128p_type
+
+class Coder:
+
+    def __init__(self, workdir=None):
+        self.workdir = workdir = workdir or os.path.abspath(".mcvine.acc.coder")
+        if not os.path.exists(workdir):
+            os.makedirs(workdir)
+        self.modules = dict()
+
+    def getModule(self, type, N):
+        container = coder.createDir(type)
+        modulepath = os.path.join(container, f'compiled_{N}.py')
+        key = type, N
+        if os.path.exists(modulepath) and key not in self.modules:
+            self.modules[key] = modulepath
+        return modulepath
+
+    def createUniqueDir(self, prefix):
+        return tempfile.mkdtemp(prefix=prefix, dir=self.workdir)
+
+    def createDir(self, name):
+        wd = os.path.join(self.workdir, name)
+        if not os.path.exists(wd):
+            os.makedirs(wd)
+        return wd
+
+    @classmethod
+    def unrollLoop(cls, N, indent=4*' ', before_loop=None, in_loop=None, after_loop=None):
+        lines=[indent+line for line in before_loop or []]
+        for i in range(N):
+            lines+=[indent+line.format(i=i) for line in in_loop or []]
+        lines+=[indent+line for line in after_loop or []]
+        return lines
+
+coder = Coder()

acc/components/samples/composite.py

@@ -0,0 +1,48 @@
+import os
+import numpy as np, numba
+from numba import cuda, void, int64
+from numba.cuda.random import xoroshiro128p_uniform_float32, xoroshiro128p_type
+from ...config import get_numba_floattype
+NB_FLOAT = get_numba_floattype()
+
+from .SampleBase import SampleBase
+
+category = 'samples'
+
+def sampleassembly_from_xml(name, samplexml):
+    """factory method to create sample assembly component.
+    can only be used in nonacc instrument script
+    """
+    klass = SampleAssemblyFromXml(samplexml)
+    return klass(name)
+
+def SampleAssemblyFromXml(samplexml):
+    from . import loadScattererComposite
+    composite = loadScattererComposite(samplexml)
+    return factory(composite)
+
+def factory(composite):
+    from ...scatterers import scatter_func_factory
+    methods = scatter_func_factory.render(composite)
+    scatter = methods['scatter']
+    class Composite(SampleBase):
+
+        def __init__(self, name):
+            self.name = name
+            self.propagate_params = ()
+
+            # Aim neutrons toward the sample to cause JIT compilation.
+            import mcni
+            neutrons = mcni.neutron_buffer(1)
+            neutrons[0] = mcni.neutron(r=(0, 0, -1), v=(0, 0, 1), prob=1, time=0)
+            self.process(neutrons)
+
+        @cuda.jit(
+            void(int64, xoroshiro128p_type[:], NB_FLOAT[:]),
+            device=True, inline=True,
+        )
+        def propagate(threadindex, rng_states, neutron):
+            return
+    Composite.propagate = scatter
+    Composite.register_propagate_method(scatter)
+    return Composite

acc/geometry/arrow_intersect.py

@@ -26,6 +26,18 @@ def render(self, shape):
         return shape.identify(self)
 
     def onUnion(self, u):
+        nelements = len(u.shapes)
+        if nelements == 1:
+            return u.shapes[0].identify(self)
+        elif nelements == 2:
+            return self.onUnion2(u)
+        elif nelements == 3:
+            from .composite_3 import createRayTracingMethods_NonOverlappingShapes as createMethods
+            return createMethods(u.shapes)['intersect_all']
+        else:
+            raise NotImplementedError(f"locate for union of {nelements} elements")
+
+    def onUnion2(self, u):
         locate1 = self.locate_func_factory.onUnion(u)
         s1, s2 = u.shapes
         f1 = s1.identify(self)
@@ -85,6 +97,22 @@ def intersectBlock(x, y, z, vx, vy, vz, ts):
 
         return intersectBlock
 
+    def onTranslation(self, t):
+        x0,y0,z0 = [_/units.length.meter for _ in t.vector]
+        body = t.body
+        intersect_body = self.render(body)
+        @cuda.jit(device=True, inline=True)
+        def intersectTranslation(x,y,z, vx,vy,vz, ts):
+            return intersect_body(x-x0,y-y0,z-z0, vx,vy,vz, ts)
+        return intersectTranslation
+
+    def onRotation(self, r):
+        euler_angles = [_/units.angle.rad for _ in r.euler_angles]
+        # hack
+        if not np.allclose(euler_angles, [0,0,0]):
+            raise NotImplementedError
+        return self.render(r.body)
+
     def onDifference(self, s):
         locate1 = self.locate_func_factory.onDifference(s)
         f1 = s.op1.identify(self)
@@ -150,8 +178,13 @@ def remove_item(idx, l, N):
 
 @cuda.jit(device=True)
 def insert_into_sorted_list(d, l, N):
-    'insert data "d" into existing sorted array (low to high) of length N'
-    if N>=len(l): return N
+    '''insert data "d" into existing sorted array (low to high) of length N
+    and keep the lower number when there are too many elements
+    '''
+    if N>=len(l):
+        if l[N-1] <= d:
+            return N
+        N = N - 1
     if N==0:
         l[0] = d
         return 1