numpy.r_

Here are the examples of the python api numpy.r_ taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.

def create_rope_demo(env, rope_poss):
    rope_sim_obj = create_rope(rope_poss)
    env.sim.add_objects([rope_sim_obj])
    env.sim.settle()
    scene_state = env.observe_scene()
    env.sim.remove_objects([rope_sim_obj])
    
    pick_pos = rope_poss[0] + .1 * (rope_poss[1] - rope_poss[0])
    drop_pos = rope_poss[3] + .1 * (rope_poss[2] - rope_poss[3]) + np.r_[0, .2, 0]
    pick_R = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]])
    drop_R = np.array([[0, 1, 0], [0, 0, -1], [-1, 0, 0]])
    move_height = .2
    aug_traj = create_augmented_traj(env.sim.robot, pick_pos, drop_pos, pick_R, drop_R, move_height)
    
    demo = Demonstration("rope_demo", scene_state, aug_traj)
    return demo

def draw_finger_pts_traj(sim_env, flr2finger_pts_traj, color):
    handles = []
    for finger_lr, pts_traj in flr2finger_pts_traj.items():
        for pts in pts_traj:
            handles.append(sim_env.env.drawlinestrip(np.r_[pts, pts[0][None,:]], 1, color))
    return handles

def retime2(positions, vel_limits_j):
    good_inds = np.r_[True,(abs(positions[1:] - positions[:-1]) >= 1e-6).any(axis=1)]
    positions = positions[good_inds]
    
    move_nj = positions[1:] - positions[:-1]
    dur_n = (np.abs(move_nj) / vel_limits_j[None,:]).max(axis=1) # time according to velocity limit
    times = np.cumsum(np.r_[0,dur_n])
    return times, good_inds    

def make_traj_with_limits(positions, vel_limits_j, acc_limits_j):
    times, inds = retime(positions, vel_limits_j, acc_limits_j)
    positions = positions[inds]

    deg = min(3, len(positions) - 1)
    s = len(positions)*.001**2
    (tck, _) = si.splprep(positions.T, s=s, u=times, k=deg)
    smooth_positions = np.r_[si.splev(times,tck,der=0)].T
    smooth_velocities = np.r_[si.splev(times,tck,der=1)].T    
    return smooth_positions, smooth_velocities, times    

    @staticmethod
    def get_objective2(x_nd, y_md, f, g, corr_nm, rad):
        r"""Returns the following 10 objectives:
        
            - :math:`\frac{1}{n} \sum_{i=1}^n \sum_{j=1}^m m_{ij} ||y_j - f(x_i)||_2^2`
            - :math:`\lambda Tr(A_f^\top K A_f)`
            - :math:`Tr((B_f - I) R (B_f - I))`
            - :math:`\frac{2T}{n} \sum_{i=1}^n \sum_{j=1}^m m_{ij} \log m_{ij}`
            - :math:`-\frac{2T}{n} \sum_{i=1}^n \sum_{j=1}^m m_{ij}`
            - :math:`\frac{1}{m} \sum_{j=1}^m \sum_{i=1}^n m_{ij} ||x_i - g(y_j)||_2^2`
            - :math:`\lambda Tr(A_g^\top K A_g)`
            - :math:`Tr((B_g - I) R (B_g - I))`
            - :math:`\frac{2T}{m} \sum_{j=1}^m \sum_{i=1}^n m_{ij} \log m_{ij}`
            - :math:`-\frac{2T}{m} \sum_{j=1}^m \sum_{i=1}^n m_{ij}`
        """
        cost = np.r_[TpsRpmRegistration.get_objective2(x_nd, y_md, f, corr_nm, rad), 
                     TpsRpmRegistration.get_objective2(y_md, x_nd, g, corr_nm.T, rad)]
        return cost

    def cost(self, demo, test_scene_state):
        """Gets the costs of the forward and backward thin plate spline 
        objective of the resulting registration
        
        Args:
            demo: Demonstration which has the demonstration scene
            test_scene_state: SceneState of the test scene
        
        Returns:
            A 1-dimensional numpy.array containing the residual, bending and 
            rotation cost of the forward and backward spline, each already 
            premultiplied by the respective coefficients.
        """
        reg = self.register(demo, test_scene_state, callback=None)
        cost = np.r_[reg.f.get_objective(), reg.g.get_objective()]
        return cost

def plot_polygon(ax, polygon):
    from sfepy.postprocess.plot_dofs import _get_axes

    dim = polygon.shape[1]
    ax = _get_axes(ax, dim)

    pp = nm.r_[polygon, polygon[:1]]
    px, py = pp[:, 0], pp[:, 1]
    if dim == 2:
        ax.plot(px, py)

    else:
        pz = pp[:, 2]
        ax.plot(px, py, pz)

    return ax

    def deriv(self, m, v=None):
        nC = self.nP/2
        shp = (nC, nC*2)

        def fwd(v):
            return v[:nC] + v[nC:]*1j

        def adj(v):
            return np.r_[v.real, v.imag]
        if v is not None:
            return LinearOperator(shp, matvec=fwd, rmatvec=adj) * v
        return LinearOperator(shp, matvec=fwd, rmatvec=adj)

    @property
    def vectorNy(self):
        """Nodal grid vector (1D) in the y direction."""
        if self.isSymmetric:
            # There aren't really any nodes, but all the grids need
            # somewhere to live, why not zero?!
            return np.r_[0]
        return np.r_[0, self.hy[:-1].cumsum()] + self.hy[0]*0.5

    def test_timeProblem_setTimeSteps(self):
        self.prob.timeSteps = [(1e-6, 3), 1e-5, (1e-4, 2)]
        trueTS = np.r_[1e-6, 1e-6, 1e-6, 1e-5, 1e-4, 1e-4]
        self.assertTrue(np.all(trueTS == self.prob.timeSteps))

        self.prob.timeSteps = trueTS
        self.assertTrue(np.all(trueTS == self.prob.timeSteps))

        self.assertTrue(self.prob.nT == 6)

        self.assertTrue(np.all(self.prob.times == np.r_[0, trueTS].cumsum()))

    def test_slices(self):
        expMap = Maps.ExpMap(Mesh.TensorMesh((3,)))
        PM = MyPropMap({'maps':[('sigma', expMap)], 'slices':{'sigma':[2,1,0]}})
        assert PM.sigmaIndex == [2,1,0]
        m = PM(np.r_[1.,2,3])
        assert np.all(m.sigmaModel == np.r_[3,2,1])
        assert np.all(m.sigma == np.exp(np.r_[3,2,1]))

    def test_Projections(self):
        m = Mesh.TensorMesh((3,))
        iMap = Maps.IdentityMap(m)
        PM = MyReciprocalPropMap([('sigma', iMap)])
        v = np.r_[1,2.,3]
        pm = PM(v)

        assert pm.sigmaProj is not None
        assert pm.rhoProj   is None
        assert pm.muProj    is None
        assert pm.muiProj   is None

        assert np.all(pm.sigmaProj * v == pm.sigmaModel)

    def test_sourceIndex(self):
        survey = self.D.survey
        srcs = survey.srcList
        assert survey.getSourceIndex([srcs[1],srcs[0]]) == [1,0]
        assert survey.getSourceIndex([srcs[1],srcs[2],srcs[2]]) == [1,2,2]
        SrcNotThere = Survey.BaseSrc(srcs[0].rxList, loc=np.r_[0,0,0])
        self.assertRaises(KeyError, survey.getSourceIndex, [SrcNotThere])
        self.assertRaises(KeyError, survey.getSourceIndex, [srcs[1],srcs[2],SrcNotThere])

    def test_edgeCurl(self):

        hx, hy, hz = np.r_[1.,2,3,4], np.r_[5.,6,7,8], np.r_[9.,10,11,12]
        M = Mesh.TreeMesh([hx, hy, hz], levels=2)
        M.refine(lambda xc:2)
        # M.plotGrid(showIt=True)
        Mr = Mesh.TensorMesh([hx, hy, hz])

        # plt.subplot(211).spy(Mr.faceDiv)
        # plt.subplot(212).spy(M.permuteCC.T*M.faceDiv*M.permuteF)
        # plt.show()

        assert (Mr.edgeCurl - M.permuteF*M.edgeCurl*M.permuteE.T).nnz == 0

    def test_faceInnerProduct(self):

        hx, hy, hz = np.r_[1.,2,3,4], np.r_[5.,6,7,8], np.r_[9.,10,11,12]
        # hx, hy, hz = [[(1,4)], [(1,4)], [(1,4)]]

        M = Mesh.TreeMesh([hx, hy, hz], levels=2)
        M.refine(lambda xc:2)
        # M.plotGrid(showIt=True)
        Mr = Mesh.TensorMesh([hx, hy, hz])

        # plt.subplot(211).spy(Mr.getFaceInnerProduct())
        # plt.subplot(212).spy(M.getFaceInnerProduct())
        # plt.show()

        # print(M.nC, M.nF, M.getFaceInnerProduct().shape, M.permuteF.shape)

        assert np.allclose(Mr.getFaceInnerProduct().todense(), (M.permuteF * M.getFaceInnerProduct() * M.permuteF.T).todense())
        assert np.allclose(Mr.getEdgeInnerProduct().todense(), (M.permuteE * M.getEdgeInnerProduct() * M.permuteE.T).todense())

    def test_VectorIdenties(self):
        hx, hy, hz = [[(1,4)], [(1,4)], [(1,4)]]

        M = Mesh.TreeMesh([hx, hy, hz], levels=2)
        Mr = Mesh.TensorMesh([hx, hy, hz])

        assert (M.faceDiv * M.edgeCurl).nnz == 0
        assert (Mr.faceDiv * Mr.edgeCurl).nnz == 0

        hx, hy, hz = np.r_[1.,2,3,4], np.r_[5.,6,7,8], np.r_[9.,10,11,12]

        M = Mesh.TreeMesh([hx, hy, hz], levels=2)
        Mr = Mesh.TensorMesh([hx, hy, hz])

        assert np.max(np.abs((M.faceDiv * M.edgeCurl).todense().flatten())) < TOL
        assert np.max(np.abs((Mr.faceDiv * Mr.edgeCurl).todense().flatten())) < TOL

    def _raveled_index(self):
        """
        Flattened array of ints, specifying the index of this object.
        This has to account for shaped parameters!
        """
        return np.r_[:self.size]

    def test_remove(self):
        removed = self.param_index.remove(three, np.r_[3:13])
        self.assertListEqual(removed.tolist(), [4,7,10])
        self.assertListEqual(self.param_index[three].tolist(), [2])
        removed = self.param_index.remove(one, [1])
        self.assertListEqual(removed.tolist(), [])
        self.assertListEqual(self.param_index[one].tolist(), [3,9])
        self.assertListEqual(self.param_index.remove('not in there', []).tolist(), [])
        removed = self.param_index.remove(one, [9])
        self.assertListEqual(removed.tolist(), [9])
        self.assertListEqual(self.param_index[one].tolist(), [3])
        self.assertListEqual(self.param_index.remove('not in there', [2,3,4]).tolist(), [])
        self.assertListEqual(self.view.remove('not in there', [2,3,4]).tolist(), [])

def gen_gendat_ar0(ar):
    def gendat_ar0(msg = False):
        ars = AR_simulator()
        ars.ngroups = 200
        ars.params = np.r_[0, -1, 1, 0, 0.5]
        ars.error_sd = 2
        ars.dep_params = [ar,]
        ars.simulate()
        return ars, Autoregressive()
    return gendat_ar0

def gen_gendat_ar1(ar):
    def gendat_ar1():
        ars = AR_simulator()
        ars.ngroups = 200
        ars.params = np.r_[0, -0.8, 1.2, 0, 0.5]
        ars.error_sd = 2
        ars.dep_params = [ar,]
        ars.simulate()
        return ars, Autoregressive()
    return gendat_ar1

def gendat_nested1():
    ns = Nested_simulator()
    ns.error_sd = 2.
    ns.params = np.r_[0, 1, 1.3, -0.8, -1.2]
    ns.ngroups = 50
    ns.nest_sizes = [10, 5]
    ns.dep_params = [1., 3.]
    ns.simulate()
    return ns, Nested(ns.id_matrix)

def gendat_nested1():
    ns = Nested_simulator()
    ns.error_sd = 2.
    ns.params = np.r_[0, 1, 1.3, -0.8, -1.2]
    ns.ngroups = 50
    ns.nest_sizes = [10, 5]
    ns.dparams = [1., 3.]
    ns.simulate()
    return ns, Nested(ns.id_matrix)

@dec.skipif(not have_matplotlib)
def test_plot_acf():
    # Just test that it runs.
    fig = plt.figure()
    ax = fig.add_subplot(111)

    ar = np.r_[1., -0.9]
    ma = np.r_[1., 0.9]
    armaprocess = tsp.ArmaProcess(ar, ma)
    rs = np.random.RandomState(1234)
    acf = armaprocess.generate_sample(100, distrvs=rs.standard_normal)
    plot_acf(acf, ax=ax, lags=10)
    plot_acf(acf, ax=ax)
    plot_acf(acf, ax=ax, alpha=None)

    plt.close(fig)

@dec.skipif(not have_matplotlib)
def test_plot_acf_irregular():
    # Just test that it runs.
    fig = plt.figure()
    ax = fig.add_subplot(111)

    ar = np.r_[1., -0.9]
    ma = np.r_[1., 0.9]
    armaprocess = tsp.ArmaProcess(ar, ma)
    rs = np.random.RandomState(1234)
    acf = armaprocess.generate_sample(100, distrvs=rs.standard_normal)
    plot_acf(acf, ax=ax, lags=np.arange(1, 11))
    plot_acf(acf, ax=ax, lags=10, zero=False)
    plot_acf(acf, ax=ax, alpha=None, zero=False)

    plt.close(fig)

@dec.skipif(not have_matplotlib)
def test_plot_pacf():
    # Just test that it runs.
    fig = plt.figure()
    ax = fig.add_subplot(111)

    ar = np.r_[1., -0.9]
    ma = np.r_[1., 0.9]
    armaprocess = tsp.ArmaProcess(ar, ma)
    rs = np.random.RandomState(1234)
    pacf = armaprocess.generate_sample(100, distrvs=rs.standard_normal)
    plot_pacf(pacf, ax=ax)
    plot_pacf(pacf, ax=ax, alpha=None)

    plt.close(fig)

@dec.skipif(not have_matplotlib)
def test_plot_pacf_irregular():
    # Just test that it runs.
    fig = plt.figure()
    ax = fig.add_subplot(111)

    ar = np.r_[1., -0.9]
    ma = np.r_[1., 0.9]
    armaprocess = tsp.ArmaProcess(ar, ma)
    rs = np.random.RandomState(1234)
    pacf = armaprocess.generate_sample(100, distrvs=rs.standard_normal)
    plot_pacf(pacf, ax=ax, lags=np.arange(1, 11))
    plot_pacf(pacf, ax=ax, lags=10, zero=False)
    plot_pacf(pacf, ax=ax, alpha=None, zero=False)

    plt.close(fig)

def revrt(X,m=None):
    """
    Inverse of forrt. Equivalent to Munro (1976) REVRT routine.
    """
    if m is None:
        m = len(X)
    i = int(m // 2+1)
    y = X[:i] + np.r_[0,X[i:],0]*1j
    return np.fft.irfft(y)*m

    def unmapped_indices(self, region_mapping):
        """
        Compute vector of indices of regions in connectivity which are not in the given
        region mapping.

        """

        return numpy.setdiff1d(numpy.r_[:self.number_of_regions], region_mapping)

    def _lri(self, nnz_row_el_idx):
        "Flat array of indices afferent, non-zero-weight connections."
        if not hasattr(self, '_cached_lri'):
            rows = numpy.r_[-1, nnz_row_el_idx]
            self._cached_lri, = numpy.argwhere(numpy.diff(rows)).T
            self._cached_nzr = numpy.unique(nnz_row_el_idx)
            LOG.debug('lri.size %d nzr.size %d', self._cached_lri.size, self._cached_nzr.size)
        return self._cached_lri, self._cached_nzr

    def config_for_sim(self, simulator):
        """Configure monitor for given simulator.

        Grab the Simulator's integration step size. Set the monitor's variables
        of interest based on the Monitor's 'variables_of_interest' attribute, if
        it was specified, otherwise use the 'variables_of_interest' specified 
        for the Model. Calculate the number of integration steps (isteps)
        between returns by the record method. This method is called from within
        the the Simulator's configure() method.

        """
        self.dt = simulator.integrator.dt
        self.istep = iround(self.period / self.dt)
        self.voi = self.variables_of_interest
        if self.voi is None or self.voi.size == 0:
            self.voi = numpy.r_[:len(simulator.model.variables_of_interest)]

    def _prepare_stimulus(self):
        if self.stimulus is None:
            stimulus = 0.0
        else:
            time = numpy.r_[0.0 : self.simulation_length : self.integrator.dt]
            self.stimulus.configure_time(time.reshape((1, -1)))
            stimulus = numpy.zeros((self.model.nvar, self.number_of_nodes, 1))
            LOG.debug("stimulus shape is: %s", stimulus.shape)
        return stimulus

    def test_clamp(self):
        vode = integrators.VODE(
            clamped_state_variable_indices=numpy.r_[0, 3],
            clamped_state_variable_values=numpy.array([[42.0, 24.0]])
        )
        x = numpy.ones((5, 4, 2))
        for i in range(10):
            x = vode.scheme(x, self._dummy_dfun, 0.0, 0.0, 0.0)
        for idx, val in zip(vode.clamped_state_variable_indices, vode.clamped_state_variable_values):
            self.assertTrue(numpy.allclose(x[idx], val))

def smooth_hanning(x, size=11):
    """smooth a 1D array using a hanning window with requested size."""

    if x.ndim != 1:
        raise ValueError, "smooth_hanning only accepts 1-D arrays."
    if x.size < size:
        raise ValueError, "Input vector needs to be bigger than window size."
    if size < 3:
        return x

    s = np.r_[x[size - 1:0:-1], x, x[-1:-size:-1]]
    w = np.hanning(size)
    y = np.convolve(w / w.sum(), s, mode='valid')
    return y

def Concat(y,x=None):
	if x is None:
		x = Tensor.context
	#print x.value.mean()
	#print y.value.sum()
	dat = np.r_[x.value,y.value];
	#print dat.sum();
	x = Variable(dat, volatile=True)
	t = ChainerTensor(x	)
	t.use()
	return t

    def setUp(self):
        spread = N.r_[50:101:10]
        self.pos = N.zeros((2*len(spread), 3))
        self.pos[:len(spread), 0] = spread
        self.pos[len(spread):, 1] = spread
        self.pos[:,2] = 4.5
        
        self.field = HeliostatField(self.pos, 8., 8., 0., 90.)
        s2 = N.sqrt(2)/2
        self.sunvec = N.r_[-s2, 0, s2] # Noon, winterish.
        
        ray_pos = (self.pos + self.sunvec).T
        ray_dir = N.tile(-self.sunvec, (self.pos.shape[0], 1)).T
        self.rays = RayBundle(ray_pos, ray_dir, energy=N.ones(self.pos.shape[0]))

    def test_vertical(self):
        """Finite cylinder parallel to the Z axis"""
        correct_prm = N.r_[0.5, 0, 0.5, 0, 0.50062461, 1.50187383, 0.30594117, 0.]
        prm = self.gm.find_intersections(N.eye(4), self.bund)
        N.testing.assert_array_equal(prm[[1,3, 7]], N.ones(3)*N.inf)
        prm[[1,3, 7]] = 0.
        N.testing.assert_array_almost_equal(prm, correct_prm)
        
        self.gm.select_rays(N.r_[0, 2])
        
        correct_norm = N.c_[[0.,-1., 0.], [0.,1., 0.]]
        norm = self.gm.get_normals()
        N.testing.assert_array_almost_equal(norm, correct_norm)
        
        correct_pts = N.c_[[0., 0.5, 0.], [0., 0.5, 0.]]
        pts = self.gm.get_intersection_points_global()
        N.testing.assert_array_almost_equal(pts, correct_pts)

    def test_select_rays_normals(self):
        """A tilted flat geometry manager returns normals only for selected rays"""
        s2 = math.sqrt(2)
        self.gm.select_rays(N.r_[1,3])
        n = self.gm.get_normals()
        N.testing.assert_array_almost_equal(n, N.tile(N.c_[[0,1/s2,1/s2]], (1,2)))

    def setUp(self):
        self.eighth_circle_trans = generate_transform(N.r_[1., 0, 0], N.pi/4, 
            N.c_[[0., 1, 0]])
        
        self.surf = Surface(flat_surface.FlatGeometryManager(), \
            optics_callables.perfect_mirror)
        self.obj = AssembledObject(surfs=[self.surf])
        self.sub_assembly = Assembly()
        self.sub_assembly.add_object(self.obj, self.eighth_circle_trans)
        self.assembly = Assembly()
        self.assembly.add_assembly(self.sub_assembly, self.eighth_circle_trans)

    def test_multiple_ref_idxs(self):
        n1 = N.r_[1., 1.]
        n2 = N.r_[1., 1.5]
        dir = N.c_[[0, -1, 1], [0, -1, -1]]/math.sqrt(2)
        norm = N.c_[[0, 1, 0]]
        refr, dirs = optics.refractions(n1, n2, dir, norm)
        
        self.failUnless(refr.all(), "Some rays did not refract")
        correct_refrs = N.c_[dir[:,0], -N.sqrt([0, 7./9, 2./9])]
        N.testing.assert_array_almost_equal(dirs, correct_refrs)

    def test_TIR(self):
        """With rays exiting a glass, some are TIR'ed, some not"""
        n1 = 1.5
        n2 = 1.
        dir = N.c_[[0, -1, 0], N.r_[0, -1, -1]/math.sqrt(2)]
        norm = N.c_[[0, -1, 0]]
        refr, dirs = optics.refractions(n1, n2, dir, norm)
        
        N.testing.assert_array_equal(refr, N.r_[True, False])
        N.testing.assert_array_equal(dirs, N.c_[[0,-1,0]])

    def test_transformed(self):
        """Translated and rotated dish rejects missing rays"""
        trans = generate_transform(N.r_[1., 0., 0.], N.pi/4., N.c_[[0., 0., 1.]])
        
        self.surf.transform_frame(trans)
        misses = N.isinf(self.surf.register_incoming(self.bund))
        N.testing.assert_array_equal(misses, N.r_[False, False, True, True])

    def setUp(self):
        surface = Surface(HemisphereGM(1.), opt.perfect_mirror,
            rotation=general_axis_rotation(N.r_[1,0,0], N.pi))
        self._bund = RayBundle(energy=N.ones(3))
        self._bund.set_directions(N.c_[[0,1,0],[0,1,0],[0,-1,0]])
        self._bund.set_vertices(N.c_[[0,-2.,0.001],[0,0,0.001],[0,2,0.001]])

        assembly = Assembly()
        object = AssembledObject()
        object.add_surface(surface)
        assembly.add_object(object)
        
        self.engine = TracerEngine(assembly)

    def setUp(self):
        self.assembly = Assembly()

        surface1 = Surface(FlatGeometryManager(), opt.perfect_mirror)
        self.object1 = AssembledObject()
        self.object1.add_surface(surface1)  
        
        boundary = BoundarySphere(location=N.r_[0,0.,3], radius=3.)
        surface3 = Surface(CutSphereGM(2., boundary), opt.perfect_mirror)
        self.object2 = AssembledObject()
        self.object2.add_surface(surface3)

        self.transform1 = generate_transform(N.r_[1.,0,0], N.pi/4, N.c_[[0,0,-1.]])
        self.transform2 = translate(0., 0., 2.)
        self.assembly.add_object(self.object1, self.transform1)
        self.assembly.add_object(self.object2, self.transform2)

    def test_assembly3(self):
        """Tests the assembly after three iterations"""
        self.engine.ray_tracer(self._bund,3,.05)[0]
        params = self.engine.tree.ordered_parents()
        correct_params = [N.r_[1,2],N.r_[0,0,1,1],N.r_[1,2,1,2]]
        N.testing.assert_equal(params, correct_params)

    def test_final_order(self):
        """Rays come out of a homogenizer with the right parents order"""
        self.engine.ray_tracer(self.bund, 300, .05)
        parents = self.engine.tree.ordered_parents()
        correct_parents = [N.r_[0, 1, 2, 3], N.r_[1, 0, 3, 2]]
        N.testing.assert_equal(parents, correct_parents)

    def setUp(self):
        absorptive = Surface(RectPlateGM(1., 1.), opt.Reflective(1.),
            location=N.r_[ 0.5, 0., 1.])
        reflective = Surface(RectPlateGM(1., 1.), opt.Reflective(0.),
            location=N.r_[-0.5, 0., 1.])
        self.assembly = Assembly(
            objects=[AssembledObject(surfs=[absorptive, reflective])])
        
        # 4 rays: two toward absorptive, two toward reflective.
        pos = N.zeros((3,4))
        pos[0] = N.r_[0.5, 0.25, -0.25, -0.5]
        direct = N.zeros((3,4))
        direct[2] = 1.
        self.bund = RayBundle(pos, direct, energy=N.ones(4))

    def test_absorbed_to_back(self):
        """Absorbed rays moved to back of recorded bundle"""
        engine = TracerEngine(self.assembly)
        engine.ray_tracer(self.bund, 300, .05)
        
        parents = engine.tree.ordered_parents()
        N.testing.assert_equal(parents, [N.r_[2,3, 0, 1]])

    def test_single_vec(self):
        """A rotation into one vector is correct"""
        vec = np.r_[1., 1., 1.]/np.sqrt(3)
        rot = sg.rotation_to_z(vec)
        np.testing.assert_array_almost_equal(rot, np.c_[
            np.r_[1., -1, 0.]/np.sqrt(2), np.r_[1., 1., -2.]/np.sqrt(6), vec])

    def test_two_vecs(self):
        """Vectorization of rotation into a vector"""
        vecs = np.vstack((np.r_[1., 1., 1.]/np.sqrt(3), np.r_[1., 0., 0.]))
        rots = sg.rotation_to_z(vecs)
        
        np.testing.assert_array_almost_equal(rots[0], np.c_[
            np.r_[1., -1, 0.]/np.sqrt(2), np.r_[1., 1., -2.]/np.sqrt(6), vecs[0]])
        np.testing.assert_array_almost_equal(rots[1], np.c_[
            [0., -1., 0.], [0., 0., -1.], [1., 0., 0.]])

    def test_complex_matlab_blobs(self):
        blobs = Blob().fetch.order_by('id')['blob']
        assert_equal(blobs[0][0], 'character string')
        assert_true(np.array_equal(blobs[1][0], np.r_[1:180:15]))
        assert_list_equal([r[0] for r in blobs[2]], ['string1', 'string2'])
        assert_list_equal([r[0, 0] for r in blobs[3]['a'][0]], [1, 2])
        assert_tuple_equal(blobs[3]['b'][0, 0]['c'][0, 0].shape, (3, 3))
        assert_true(np.array_equal(blobs[4], np.r_[1:25].reshape((2, 3, 4), order='F')))
        assert_true(blobs[4].dtype == 'float64')
        assert_true(np.array_equal(blobs[5], np.r_[1:25].reshape((2, 3, 4), order='F')))
        assert_true(blobs[5].dtype == 'uint8')
        assert_tuple_equal(blobs[6].shape, (2, 3, 4))
        assert_true(blobs[6].dtype == 'complex128')