def main(self, *args, **kwds):
        if self.inventory.dumpregistry:
            self._dumpRegsitry()
            return
        
        import mcni
        instrument = self._createInstrument()
        
        geometer = self.geometer

        context = self._makeSimContext()
        
        if self.ncount < self.buffer_size:
            self.buffer_size = int(self.ncount)
        n = int(self.ncount / self.buffer_size)
        
        from mcni import journal
        logger = journal.logger(
            'info', 'instrument', header='', footer='', format='-> %s')
        for i in range(n):
            logger("mpi node %s at loop %s" % (self.mpi.rank, i))
            neutrons = mcni.neutron_buffer( self.buffer_size )
            context.iteration_no = i
            mcni.simulate( instrument, geometer, neutrons, context=context)
            continue
        
        remain = int(self.ncount % self.buffer_size)
        if remain:
            logger("mpi node %s at last loop" % (self.mpi.rank,))
            neutrons = mcni.neutron_buffer(remain)
            context.iteration_no = n
            mcni.simulate( instrument, geometer, neutrons, context=context)
            
        print os.times()
        return

 def test4(self):
     "neutron buffer: swap"
     from mcni import neutron_buffer
     b1 = neutron_buffer(1)
     b2 = neutron_buffer(2)
     b1.swap(b2)
     self.assertEqual(len(b1), 2)
     self.assertEqual(len(b2), 1)
     return

 def test5(self):
     "neutron buffer: appendNeutrons"
     from mcni import neutron_buffer
     b1 = neutron_buffer(1)
     b2 = neutron_buffer(2)
     
     b1.appendNeutrons(b2, 0, len(b2))
     self.assertEqual(len(b1), 3)
     
     b1.appendNeutrons(b2)
     self.assertEqual(len(b1), 5)
     return

 def test2a(self):
     'normalize: large buffer'
     import mcni
     # create a dummy input
     neutrons = mcni.neutron_buffer( int(3e6) )
     narr = neutrons.to_npyarr()
     narr[:, -1] = 1
     neutrons.from_npyarr(narr)
     out = 'tmp-nst-test2a.ns'
     mns.dump(neutrons, out)
     del neutrons
     
     # try normalize out-of-place
     out2 = 'tmp-nst-test2a-normalized.ns'
     if os.path.exists(out2): os.remove(out2)
     mns.normalize(out, 10., out2)
     neutrons2 = mns.load(out2)
     # and see if it is done correctly
     narr = neutrons2.to_npyarr()
     self.assertTrue((narr[:, -1] == .1).all())
     del neutrons2, narr
     
     # try normalize in-place
     mns.normalize(out, 10.)
     neutrons2 = mns.load(out)
     # and see if it is done correctly
     narr = neutrons2.to_npyarr()
     self.assertTrue((narr[:, -1] == .1).all())
     return

    def test2(self):
        'normalize'
        import mcni
        # create a dummy input
        neutrons = mcni.neutron_buffer( 10 )
        for n in neutrons:
            n.probability = 1
            continue
        out = 'tmp-nst-test2.ns'
        mns.dump(neutrons, out)

        # try normalize out-of-place
        out2 = 'tmp-nst-test2-normalized.ns'
        if os.path.exists(out2): os.remove(out2)
        mns.normalize(out, 10., out2)
        neutrons2 = mns.load(out2)
        # and see if it is done correctly
        for n in neutrons2:
            self.assertAlmostEqual(n.probability, .1)
            continue

        # try normalize in-place
        mns.normalize(out, 10.)
        neutrons2 = mns.load(out)
        # and see if it is done correctly
        for n in neutrons2:
            self.assertAlmostEqual(n.probability, .1)
            continue
        return

    def test1(self):
        'mccomponents.sample.samplecomponent: IsotropicKernel'
        import mcni
        neutron = mcni.neutron( r = (0,0,0), v = (0,0,3000), time = 0, prob = 1 )
        from mcni.components.MonochromaticSource import MonochromaticSource
        component1 = MonochromaticSource('source', neutron)
        from mccomponents.sample import samplecomponent
        component2 = samplecomponent( 'Al', 'sampleassemblies/Al-isotropickernel/sampleassembly.xml' )
        instrument = mcni.instrument( [component1, component2] )
        
        geometer = mcni.geometer()
        geometer.register( component1, (0,0,0), (0,0,0) )
        geometer.register( component2, (0,0,1), (0,0,0) )

        N0 = 1000
        neutrons = mcni.neutron_buffer(N0)

        mcni.simulate( instrument, geometer, neutrons )

        N = len(neutrons)

        for i in range(10):
            neutron = neutrons[i]
            print neutron
            continue

        # N should be about 2/3 of N0. this is determined by
        # the mc weights in Al-isotropic-kernel-plate-scatterer.xml
        self.assert_( N < 0.72*N0 and N > 0.6*N0)

        return

    def test(self):
        'mcni.pyre_components.NeutronsOnCone_FixedQE'
        from mcni.pyre_components.NeutronsOnCone_FixedQE import NeutronsOnCone_FixedQE as factory
        component = factory( 'source' )
        component.inventory.Q = Q
        component.inventory.E = E
        component.inventory.Ei = Ei
        component.inventory.L1 = L1
        component._configure()
        component._init()

        import mcni
        neutrons = mcni.neutron_buffer( 10 )
        component.process( neutrons )
        from numpy.linalg import norm
        from mcni.utils import v2e, e2v, v2k
        vi = e2v( Ei )
        t = L1/vi
        for n in neutrons:
            vv = n.state.velocity
            vQv = -vv[0], -vv[1], vi-vv[2]
            vQ = norm(vQv)
            Q1 = v2k(vQ)
            v = norm( vv )
            self.assertAlmostEqual( v2e(v), Ei-E, 4)
            self.assertAlmostEqual( n.time, t, 4 )
            self.assertAlmostEqual( Q, Q1, 4 )
            continue

        component._fini()
        return

    def test1(self):
        'NDMonitor'
        
        from mcni.components.NDMonitor import NDMonitor, Axis
        xaxis = Axis(
            name = 'x', expression='x',
            bins = 100, range=(0, 1000.),
            unit = 'meter',
            )
        m = NDMonitor('abc', [xaxis])
        
        N = 100
        from mcni import neutron_buffer, neutron
        b = neutron_buffer(N)
        for i in range(N):
            b[i] = neutron()
            continue

        m.process(b)

        self.assertEqual(m.histogram.I[0], N)
        self.assertEqual(m.histogram.E2[0], N)
        self.assertEqual(m.histogram.I[1], 0.)
        self.assertEqual(m.histogram.E2[1], 0.)

        if interactive:
            from histogram.plotter import defaultPlotter as plotter
            plotter.plot(m.histogram)
        return

def neutrons_from_npyarr(arr, neutrons=None):
    """copy data from a numpy array to a boost python instance of
    Neutron::Events.

    arr: the numpy array
    neutrons: the Neutron::Events instance where data will be copied.
      if None, a new instance will be created.
    """
    shape = arr.shape
    assert shape[1] == ndblsperneutron
    n = len(arr)

    if neutrons is None:
        import mcni

        neutrons = mcni.neutron_buffer(n)
        pass

    n = min(n, len(neutrons))

    cevents = binding.cevents_from_npyarr(arr)

    neutrons.fromCevents(cevents, n)

    return neutrons

 def test2(self):
     'shape positioning: cylinder with axis along beam'
     # source
     from mcni.components.MonochromaticSource import MonochromaticSource
     import mcni
     neutron = mcni.neutron(r=(0,0,-1), v=(0,0,1000), prob=1)
     source = MonochromaticSource('s', neutron, dx=0.09, dy=0.09, dE=0)
     # sample
     from mccomponents.sample import samplecomponent
     scatterer = samplecomponent('sa', 'cyl-along-beam/sampleassembly.xml' )
     # neutrons
     N = 1000
     neutrons = mcni.neutron_buffer(N)
     neutrons = source.process(neutrons)
     # find neutrons out of target
     arr = neutrons.to_npyarr()
     x,y,z = arr[:, :3].T
     missing = x*x+y*y > 0.04**2
     # print neutrons
     scatterer.process(neutrons)
     # print neutrons
     arr = neutrons.to_npyarr()
     x,y,z = arr[:, :3].T
     assert (z[missing] < -.9).all()
     hit = arr[np.logical_not(missing), :3]
     x,y,z = hit.T
     assert (z > -.1).all()
     assert (np.isclose((x*x + y*y)**.5, 0.04) | np.isclose(np.abs(z), 0.005)).all()
     return

    def test(self):
        "Source_simple --> E_monitor"
        from mcstas2 import componentfactory
        ssimplefac = componentfactory( 'sources', 'Source_simple' )
        ssimple = ssimplefac(
            'ssimple',
            radius=0.1, dist=2, xw=0.1, yh=0.1, E0=55, dE=2)
        
        from mcstas2.wrappers import wrap
        wrap( 'E_monitor.comp', 'monitors' ) 
        emonfac = componentfactory( 'monitors', 'E_monitor' )
        emon = emonfac(
            'emon',
            nchan=20, filename="e.dat",
            xmin=-0.2, xmax=0.2,
            ymin=-0.2, ymax=0.2,
            Emin=50, Emax=60)

        import mcni
        instrument = mcni.instrument( [ssimple, emon] )

        geometer = mcni.geometer()
        geometer.register( ssimple, (0,0,0), (0,0,0) )
        geometer.register( emon, (0,0,1), (0,0,0) )

        neutrons = mcni.neutron_buffer( 100 )

        mcni.simulate( instrument, geometer, neutrons )

        return

    def test3(self):
        'neutron_storage.Storage: wrap-reading (nread>ntotal)'

        path = 'test-storage-3'
        if os.path.exists(path):
            os.remove( path )
        
        from mcni.neutron_storage.Storage import Storage

        #open storage for writing
        s = Storage( path, 'w' )

        #create neutrons
        import mcni
        neutrons = mcni.neutron_buffer( 7 )
        for i in range(7):
            neutrons[i] = mcni.neutron( v = (i,0,0) )

        #write 
        s.write( neutrons )

        # flush
        del s
        
        #open the storage for reading
        sr = Storage( path, 'r', packetsize=10)

        neutrons = sr.read()
        self.assertEqual( len(neutrons), 10 )
        self.assertAlmostEqual( neutrons[5].state.velocity[0] , 5 )
        self.assertAlmostEqual( neutrons[9].state.velocity[0] , 2 )

        neutrons = sr.read()
        self.assertAlmostEqual( neutrons[0].state.velocity[0] , 3 )
        return

    def test(self):
        'neutron_storage.Storage: write and then read'

        path = 'test-storage'
        if os.path.exists(path):
            os.remove( path )
        
        from mcni.neutron_storage.Storage import Storage

        #open storage for writing
        s = Storage( path, 'w' )

        #create neutrons
        import mcni
        neutrons = mcni.neutron_buffer( 7 )
        neutrons[5] = mcni.neutron( v = (8,9,10) )

        #write neutrons
        s.write( neutrons )

        #delete the storage to make sure it flushes all neutrons
        del s

        #open the storage for reading
        sr = Storage( path, 'r')
        neutrons = sr.read()
        self.assertEqual( len(neutrons), 7 )

        self.assertAlmostEqual( neutrons[5].state.velocity[0] , 8 )
        return

    def test5(self):
        'neutron_storage.Storage: wrap-reading (nread>>ntotal)'

        path = 'test-storage-4'
        if os.path.exists(path):
            os.remove( path )
        
        from mcni.neutron_storage.Storage import Storage

        #open storage for writing
        s = Storage( path, 'w' )

        #create neutrons
        import mcni
        neutrons = mcni.neutron_buffer( 5 )
        for i in range(5):
            neutrons[i] = mcni.neutron( v = (i,0,0) )

        #write 
        s.write( neutrons )

        # flush
        del s
        
        #open the storage for reading
        sr = Storage( path, 'r')

        neutrons = sr.read(100)
        self.assertEqual( len(neutrons), 100 )
        self.assertAlmostEqual( neutrons[3].state.velocity[0] , 3 )
        self.assertAlmostEqual( neutrons[4].state.velocity[0] , 4 )
        self.assertAlmostEqual( neutrons[6].state.velocity[0] , 1 )
        self.assertAlmostEqual( neutrons[7].state.velocity[0] , 2 )

        return

 def test1(self):
     'shape positioning: plate perp to beam'
     # source
     from mcni.components.MonochromaticSource import MonochromaticSource
     import mcni
     neutron = mcni.neutron(r=(0,0,-1), v=(0,0,1000), prob=1)
     source = MonochromaticSource('s', neutron, dx=0.07, dy=0.09, dE=0)
     # sample
     from mccomponents.sample import samplecomponent
     scatterer = samplecomponent('sa', 'plate/sampleassembly.xml' )
     # neutrons
     N = 1000
     neutrons = mcni.neutron_buffer(N)
     neutrons = source.process(neutrons)
     # find neutrons out of target
     arr = neutrons.to_npyarr()
     x,y,z = arr[:, :3].T
     missing = (x>0.03) | (x<-0.03) | (y>0.04) | (y<-0.04)
     # print neutrons
     scatterer.process(neutrons)
     # print neutrons
     arr = neutrons.to_npyarr()
     x,y,z = arr[:, :3].T
     assert (z[missing] < -.9).all()
     hit = arr[np.logical_not(missing), :3]
     x,y,z = hit.T
     assert (z > -.1).all()
     assert (np.isclose(np.abs(x), 0.03) | np.isclose(np.abs(y), 0.04) | np.isclose(np.abs(z), 0.005)).all()
     return

def _normalize(input, N, output):
    """normalize the neutrons in the given file by the factor N
    and save it in the output
    """
    import os
    input = os.path.abspath(input)
    output = os.path.abspath(output)
    # guards
    assert input != output
    if os.path.exists(output):
        raise IOError('%s already exists' % output)
    # total # of neutrons
    remain = count(input)
    # output storage
    out = storage(output, mode='w')
    # input storage
    input = storage(input)
    #
    chunk = int(1e6)
    # 
    from mcni import neutron_buffer
    nb = neutron_buffer(0)
    while remain:
        n = min(remain, chunk)
        # read
        neutrons = input.read(n, asnpyarr=True)
        # normalize
        neutrons[:, -1] /= N
        # write
        nb.from_npyarr(neutrons)
        out.write(nb)
        #
        remain -= n
        continue
    return

 def test2(self):
     'shape positioning: rotated sphere'
     # source
     from mcni.components.MonochromaticSource import MonochromaticSource
     import mcni
     neutron = mcni.neutron(r=(0,0,-1), v=(0,0,1000), prob=1)
     source = MonochromaticSource('s', neutron, dx=0.10, dy=0.09, dE=0)
     # sample
     from mccomponents.sample import samplecomponent
     scatterer = samplecomponent('sa', 'sphere-rotated-arbitrarily/sampleassembly.xml' )
     # neutrons
     N = 1000
     neutrons = mcni.neutron_buffer(N)
     neutrons = source.process(neutrons)
     # find neutrons out of target
     arr = neutrons.to_npyarr()
     x,y,z = arr[:, :3].T
     missing = x*x+y*y>0.05*0.05
     # print "missed:", missing.sum()
     # print neutrons
     scatterer.process(neutrons)
     # print neutrons
     arr = neutrons.to_npyarr()
     x,y,z = arr[:, :3].T
     assert (z[missing] < -.9).all()
     hit = arr[np.logical_not(missing), :3]
     x,y,z = hit.T
     assert (z > -.1).all()
     assert np.isclose(np.abs(x*x+y*y+z*z), 0.05*0.05).all()
     return

    def test1(self):
        'mccomponents.sample.samplecomponent isotropic kernel, multiple-scattering'
        import mcni
        neutron = mcni.neutron( r = (0,0,0), v = (0,0,3000), time = 0, prob = 1 )
        from mcni.components.MonochromaticSource import MonochromaticSource
        component1 = MonochromaticSource('source', neutron)
        from mccomponents.sample import samplecomponent
        component2 = samplecomponent( 'Al', 'sampleassemblies/Al-isotropickernel/sampleassembly.xml' )
        instrument = mcni.instrument( [component1, component2] )
        
        geometer = mcni.geometer()
        geometer.register( component1, (0,0,0), (0,0,0) )
        geometer.register( component2, (0,0,1), (0,0,0) )

        N0 = 1
        neutrons = mcni.neutron_buffer(N0)

        mcni.simulate( instrument, geometer, neutrons, multiple_scattering=True)

        N = len(neutrons)

        for i in range(N):
            neutron = neutrons[i]
            print neutron
            continue

        return

 def make_neutrons(N):
     neutrons = mcni.neutron_buffer(N)
     #
     # randomly select E, the energy transfer
     E = core.Emin + np.random.random(N) * (core.Emax-core.Emin)
     # the final energy
     Ef = core.Ei - E
     # the final velocity
     vf = conversion.e2v(Ef)
     # choose cos(theta) between -1 and 1
     cos_t = np.random.random(N) * 2 - 1
     # theta
     theta = np.arccos(cos_t)
     # sin(theta)
     sin_t = np.sin(theta)
     # phi: 0 - 2pi
     phi = np.random.random(N) * 2 * np.pi
     # compute final velocity vector
     vx,vy,vz = vf*sin_t*np.cos(phi), vf*sin_t*np.sin(phi), vf*cos_t
     # neutron position, spin, tof are set to zero
     x = y = z = sx = sy = t = np.zeros(N, dtype="float64")
     # probability
     prob = np.ones(N, dtype="float64") * (vf/vi)
     # XXX: this assumes a specific data layout of neutron struct
     n_arr = np.array([x,y,z,vx,vy,vz, sx,sy, t, prob]).T.copy()
     neutrons.from_npyarr(n_arr)
     return neutrons

    def test1(self):
        'mccomponents.sample.samplecomponent'
        import mcni
        neutron = mcni.neutron( r = (0,0,0), v = (0,0,4149.48), time = 0, prob = 1 )
        from mcni.components.MonochromaticSource import MonochromaticSource
        component1 = MonochromaticSource('source', neutron)
        from mccomponents.sample import samplecomponent
        component2 = samplecomponent( 'V-constantE', 'sampleassemblies/V-constantE/sampleassembly.xml' )
        instrument = mcni.instrument( [component1, component2] )
        
        geometer = mcni.geometer()
        geometer.register( component1, (0,0,0), (0,0,0) )
        geometer.register( component2, (0,0,1), (0,0,0) )
        
        neutrons = mcni.neutron_buffer( 1 )
        
        from mcni.pyre_support.ConsoleNeutronTracer import ConsoleNeutronTracer
        tracer = ConsoleNeutronTracer()
        mcni.simulate( 
            instrument, geometer, neutrons, 
            multiple_scattering=True,
            tracer = tracer
            )

        return