def test_endian(self):
     msg = "big endian"
     a = np.arange(6, dtype='>i4').reshape((2,3))
     assert_array_equal(umt.inner1d(a,a), np.sum(a*a,axis=-1), err_msg=msg)
     msg = "little endian"
     a = np.arange(6, dtype='<i4').reshape((2,3))
     assert_array_equal(umt.inner1d(a,a), np.sum(a*a,axis=-1), err_msg=msg)

    def diffract_photons(self, photons, intersect, interpos, intercoos):
        '''Vectorized implementation'''
        p = norm_vector(photons['dir'].data[intersect])
        l, d, n = self.e_groove_coos(intercoos[intersect])
        # Minus sign here because we want n, l, d to be a right-handed coordinate system
        d = - d

        wave = energy2wave / photons['energy'].data[intersect]
        # calculate angle between normal and (ray projected in plane perpendicular to groove)
        # -> this is the blaze angle
        p_perp_to_grooves = norm_vector(p - inner1d(p, l)[:, None] * l)
        # Use abs here so that blaze angle is always in 0..pi/2
        # independent of the relative orientation of p and n.
        blazeangle = np.arccos(np.abs(inner1d(p_perp_to_grooves, n)))
        blazeangle += self.blaze_angle_modifier(intercoos[intersect, :])
        m, prob = self.order_selector(photons['energy'].data[intersect],
                                      photons['polarization'].data[intersect],
                                      blazeangle)

        # The idea to calculate the components in the (d,l,n) system separately
        # is taken from MARX
        sign = self.order_sign_convention(p, d)
        p_d = inner1d(p, d) + sign * m * wave / self.d(intercoos[intersect, :])
        p_l = inner1d(p, l)
        # The norm for p_n can be derived, but the direction needs to be chosen.
        p_n = np.sqrt(1. - p_d**2 - p_l**2)
        # Check if the photons have same direction compared to normal before
        direction = np.sign(inner1d(p, n), dtype=np.float)
        if not self.transmission:
            direction *= -1
        dir = p_d[:, None] * d + p_l[:, None] * l + (direction * p_n)[:, None] * n
        return dir, m, prob, blazeangle

 def test_type_cast(self):
     msg = "type cast"
     a = np.arange(6, dtype="short").reshape((2, 3))
     assert_array_equal(umt.inner1d(a, a), np.sum(a * a, axis=-1), err_msg=msg)
     msg = "type cast on one argument"
     a = np.arange(6).reshape((2, 3))
     b = a + 0.1
     assert_array_almost_equal(umt.inner1d(a, a), np.sum(a * a, axis=-1), err_msg=msg)

    def boundedness_selection(self, lagrangecut = 0.9):
        """get a subset of stars that are probably bound.
        
        Arguments:
        lagrangecut -- lagrangian radius outside which stars may be escapers. default: 0.9
        
        Returns:
        Nbody6Subset containing the stars that are probably bound.    
        """
        """# everything here is using the new density center as the origin.
        # if lagrangian radii aren't calculated, get the desired one.
        if self.Lagr_fractions is None:
            self.calc_lagrangian_radii([lagrangecut])
            
        # if the desired cut isn't in the already calculated ones, get it
        if lagrangecut not in self.Lagr_fractions:
            newfracs = np.hstack((self.Lagr_fractions, lagrangecut))
            newfracs.sort()
            self.calc_lagrangian_radii(newfracs)
        """
        self.calc_lagrangian_radii([lagrangecut])
        # get the desired lagrangian radius
        #rcut = self.Lagr_rads[np.where(self.Lagr_fractions == lagrangecut)]
        rcut = self.Lagr_rads[0]
        # get stars inside the cutoff
        selection = (self.Radii_dcm_new <= rcut)
        # and outside
        outside = (self.Radii_dcm_new > rcut)

        # get stars with positive v.r
        vdotr = inner1d(self.Pos - self.dc_pos_new, 
                        self.Vel - self.dc_vel_new)
        streamers = (vdotr > 0)
       
        
        # safe ones: inside or outside and not streaming
        selection = -outside | -(outside & streamers)
        possibles = outside & streamers
        possiblenames = self.Names[possibles]    
                
                
        # get velocities of stars outside the cut
        velrels = self.Vel[possibles] - self.dc_vel_new
        vel2 = np.array(inner1d(velrels, velrels))
        # get squared escape velocity from the masscut at each radius in nbody units
        vesc2 = 2 * self.Masses.sum() * lagrangecut / self.Radii_dcm_new[possibles]
        keepers = (vel2 < vesc2)
        keepnames = possiblenames[keepers] 
        selectionnames = np.vstack((self.Names[selection], keepnames))    
            
        # update radii from density center
        poscm = self.Pos - self.dc_pos
        self.Radii_dcm_new = np.array(np.sqrt(inner1d(poscm, poscm)))  
        if selection.sum() == 0:
            print 'No stars are in the sphere!'
            return None
        else:
            return Nbody6Subset(self, selectionnames) 

 def test_incontiguous_array(self):
     msg = "incontiguous memory layout of array"
     x = np.arange(64).reshape((2, 2, 2, 2, 2, 2))
     a = x[:, 0,:, 0,:, 0]
     b = x[:, 1,:, 1,:, 1]
     a[0, 0, 0] = -1
     msg2 = "make sure it references to the original array"
     assert_equal(x[0, 0, 0, 0, 0, 0], -1, err_msg=msg2)
     assert_array_equal(umt.inner1d(a, b), np.sum(a*b, axis=-1), err_msg=msg)
     x = np.arange(24).reshape(2, 3, 4)
     a = x.T
     b = x.T
     a[0, 0, 0] = -1
     assert_equal(x[0, 0, 0], -1, err_msg=msg2)
     assert_array_equal(umt.inner1d(a, b), np.sum(a*b, axis=-1), err_msg=msg)

def angle(A, B, C, degrees=True):
    """
    Returns the angle at *B* between *AB* and *BC*.

    Parameters
    ----------
    A, B, C : (*x*, *y*, *z*) triples or Nx3 arrays of triples
        Points on sphere.

    degrees : bool, optional
        If `True` (default) the result is returned in decimal degrees,
        otherwise radians.

    Returns
    -------
    angle : float or array of floats
        The angle at *B* between *AB* and *BC*, in range 0 to 2π.

    References
    ----------

    .. [1] Miller, Robert D.  Computing the area of a spherical
       polygon.  Graphics Gems IV.  1994.  Academic Press.
    """
    if HAS_C_UFUNCS:
        angle = math_util.angle(A, B, C)
    else:
        A = np.asanyarray(A)
        B = np.asanyarray(B)
        C = np.asanyarray(C)

        A, B, C = np.broadcast_arrays(A, B, C)

        ABX = _fast_cross(A, B)
        ABX = _cross_and_normalize(B, ABX)
        BCX = _fast_cross(C, B)
        BCX = _cross_and_normalize(B, BCX)
        X = _cross_and_normalize(ABX, BCX)
        diff = inner1d(B, X)
        inner = inner1d(ABX, BCX)
        with np.errstate(invalid='ignore'):
            angle = np.arccos(inner)
        angle = np.where(diff < 0.0, (2.0 * np.pi) - angle, angle)

    if degrees:
        angle = np.rad2deg(angle)

    return angle

 def calc_new_density_center(self, axes = [0, 1, 2]): 
     """Calculate the density center of a set of stars
     
     Uses the method of Casertano & Hut ApJ 1985, 298, 80. 
     
     Arguments:
     axes -- use the full 3d information (default mode) or pass in two axes to 
             use the projection along the missing axis.
     """
     positions = self.Pos[:, axes]
     #  get a nearest neighbor tree  
     kdtree = spatial.cKDTree(positions)    
     # the first result is the point itself, so sixth neighbor is the seventh result
     near6 = kdtree.query(positions, 7)[0][:,6]
     vols = np.pi * near6**2
     densities = 5.0 / vols
 
     # density center is density weighted radius of the stars
     self.dc_pos_new = (densities[:,np.newaxis] * positions).sum(0) / densities.sum()    
     self.dc_vel_new = (densities[:,np.newaxis] * 
         self.Vel[:, axes]).sum(0) / densities.sum()
         
     # update radii from density center
     poscm = self.Pos - self.dc_pos_new
     self.Radii_dcm_new = np.array(np.sqrt(inner1d(poscm, poscm)))  

 def STRONG(self, lam_max, lam, resid, X):
     """
     Screen variables using the STRONG rule.
     """
     resid_prod = np.fabs( inner1d(X.T,resid) )
     idx = resid_prod >= 2*lam_max - lam
     return np.where(idx)[0]

 def SAFE(self, lam_max, lam, y, X):
     """
     Screen variables using the SAFE rule.
     """
     resid_prod = np.fabs( inner1d(X.T,resid) )
     idx = resid_prod >= lam - la.norm(X[:,i])*la.norm(y)*((lam_max-lam)/lam_max)
     return np.where(idx)[0]

def pdf_between_data(train_data, input_data, std):
    """
    Compute PDF between two samples.

    Parameters
    ----------
    train_data : array
        Training dataset.

    input_data : array
        Input dataset

    std : float
        Standard deviation for Probability Density
        Function (PDF).

    Returns
    -------
    array-like
    """
    n_train_samples = train_data.shape[0]
    n_samples = input_data.shape[0]

    results = np.zeros((n_train_samples, n_samples))
    variance = std ** 2
    const = std * math.sqrt(2 * math.pi)

    for i, input_row in enumerate(input_data):
        inputs = np.tile(input_row, (n_train_samples, 1))
        class_difference = (train_data - inputs)
        total_distance = inner1d(class_difference, class_difference)
        results[:, i] = np.exp(-total_distance / variance) / const

    return results

def pdf_between_data(train_data, input_data, std):
    """ Compute PDF between two samples.

    Parameters
    ----------
    train_data : array
        train sample
    input_data : array
        input sample
    std : float
        standard deviation for PDF
    """
    # Note: This implementation works faster than 3D arrays
    # and use less memory.
    results = zeros((train_data.shape[0], input_data.shape[0]))
    variance = std ** 2
    function_const = std * sqrt(2 * pi)
    train_data_size = train_data.shape[0]

    for i, input_row in enumerate(input_data):
        inputs = tile(input_row, (train_data_size, 1))
        class_difference = (train_data - inputs)
        total_distance = inner1d(class_difference, class_difference)
        results[:, i] = exp(-total_distance / variance) / function_const

    return results

    def _boost_real(self, iboost, X, y, sample_weight, X_argsorted=None):
        """Implement a single boost using the SAMME.R real algorithm."""
        estimator = self._make_estimator()

        if X_argsorted is not None:
            estimator.fit(X, y, sample_weight=sample_weight,
                          X_argsorted=X_argsorted)
        else:
            estimator.fit(X, y, sample_weight=sample_weight)

        y_predict_proba = estimator.predict_proba(X)

        if iboost == 0:
            self.classes_ = getattr(estimator, 'classes_', None)
            self.n_classes_ = len(self.classes_)

        y_predict = self.classes_.take(np.argmax(y_predict_proba, axis=1),
                                       axis=0)

        # Instances incorrectly classified
        incorrect = y_predict != y

        # Error fraction
        estimator_error = np.mean(
            np.average(incorrect, weights=sample_weight, axis=0))

        # Stop if classification is perfect
        if estimator_error <= 0:
            return sample_weight, 1., 0.

        # Construct y coding as described in Zhu et al [2]:
        #
        #    y_k = 1 if c == k else -1 / (K - 1)
        #
        # where K == n_classes_ and c, k in [0, K) are indices along the second
        # axis of the y coding with c being the index corresponding to the true
        # class label.
        n_classes = self.n_classes_
        classes = self.classes_
        y_codes = np.array([-1. / (n_classes - 1), 1.])
        y_coding = y_codes.take(classes == y[:, np.newaxis])

        # Displace zero probabilities so the log is defined.
        # Also fix negative elements which may occur with
        # negative sample weights.
        y_predict_proba[y_predict_proba <= 0] = 1e-5

        # Boost weight using multi-class AdaBoost SAMME.R alg
        estimator_weight = (-1. * self.learning_rate
                                * (((n_classes - 1.) / n_classes) *
                                   inner1d(y_coding, np.log(y_predict_proba))))

        # Only boost the weights if it will fit again
        if not iboost == self.n_estimators - 1:
            # Only boost positive weights
            sample_weight *= np.exp(estimator_weight *
                                    ((sample_weight > 0) |
                                     (estimator_weight < 0)))

        return sample_weight, 1., estimator_error

 def _boost_real(self, iboost, X, y, sample_weight):
     estimator = self._make_estimator()
     try:
         estimator.set_params(random_state=self.random_state)
     except ValueError:
         pass
     estimator.fit(X, y, sample_weight=sample_weight)
     y_predict_proba = estimator.predict_proba(X)
     if iboost == 0:
         self.classes_ = getattr(estimator, 'classes_', None)
         self.n_classes_ = len(self.classes_)
     y_predict = self.classes_.take(np.argmax(y_predict_proba, axis=1), axis=0)
     # Instances incorrectly classified
     incorrect = y_predict != y
     # Error fraction
     estimator_error = np.mean(
         np.average(incorrect, weights=sample_weight, axis=0))
     # Stop if classification is perfect
     if estimator_error <= 0:
         return sample_weight, 1., 0.
     n_classes = self.n_classes_
     classes = self.classes_
     y_codes = np.array([-1. / (n_classes - 1), 1.])
     y_coding = y_codes.take(classes == y[:, np.newaxis])
     y_predict_proba[y_predict_proba <= 0] = 1e-5
     # Boost weight using multi-class AdaBoost SAMME.R alg
     estimator_weight = (-1. * self.learning_rate
                             * (((n_classes - 1.) / n_classes) *
                                inner1d(y_coding, np.log(y_predict_proba))))
     if not iboost == self.n_estimators - 1:
         # Only boost positive weights
         sample_weight *= np.exp(estimator_weight *
                                 ((sample_weight > 0) |
                                  (estimator_weight < 0)))
     return sample_weight, 1., estimator_error

def test_local_coordsys(geom):
    '''Ensure local coordiante systems are orthonormal'''
    rot = transforms3d.euler.euler2mat(*(np.pi * 2 * np.random.rand(3)))
    g = geom({'rotation': rot,
              'position': np.random.rand(3)})

    x, y, z = g.get_local_euklid_bases(np.random.rand(5, 2))

    assert np.allclose(inner1d(x, y), 0)
    assert np.allclose(inner1d(x, z), 0)

    for vec in [x, y, z]:
        assert np.allclose(np.linalg.norm(vec, axis=1), 1.)

    # Check it's a right-handed coordinage system
    assert np.allclose(np.cross(x[:, :3], y[:, :3]), z[:, :3])

def cal_rms( matA, matB ):
    L = matA.shape[0]
    if (L<1): return 0
    dx = matA - matB
    dx2 = inner1d(dx,dx)

    return np.sqrt( sum(dx2)/L )

 def sphere_selection(self, origin = None, radius = 1.0):
     """get a spherical subset of stars
     
     Arguments:
     origin -- origin of the spherical region. default: [0, 0, 0]
     radius -- radius in code units of the spherical region
     
     Returns:
     Nbody6Subset containing the stars in the requested mass range
     """
 
     if origin is not None and np.allclose(origin, self.dc_pos_new):
         rads = self.Radii_dcm_new
     if origin is not None and np.allclose(origin, self.dc_pos):
         rads = self.Radii_dcm
     if origin is None:
         rads = self.Radii
     else:
         posorigin = self.Pos - origin
         rads = np.array(np.sqrt(inner1d(posorigin, posorigin)))  
     selection = (rads < radius)
     if selection.sum() == 0:
         print 'No stars are in the sphere!'
         return None
     else:
         return Nbody6Subset(self, self.Names[selection])            

    def approx_predictive_ll(self, Arow, Acol, Wrow, Wcol, M=100):
        """
        Approximate the (marginal) predictive probability by averaging over M
        samples of the predictive parameters
        """
        from scipy.stats import multivariate_normal
        from numpy.core.umath_tests import inner1d
        import scipy.linalg

        N, B = self.N, self.B
        assert Arow.shape == Acol.shape == (N + 1,)

        # Get the predictive parameters
        lps = np.zeros(M)
        for m in xrange(M):
            Murow, Mucol, Lrow, Lcol = self.sample_predictive_parameters()

            # for n in xrange(N+1):
            #     if Arow[n]:
            #         lps[m] += multivariate_normal(Murow[n], Sigrow[n]).pdf(Wrow[n])
            #
            #     if n < N and Acol[n]:
            #         lps[m] += multivariate_normal(Mucol[n], Sigcol[n]).pdf(Wcol[n])

            # Sigrow_chol = np.array([np.linalg.cholesky(S) for S in Sigrow])
            # Sigcol_chol = np.array([np.linalg.cholesky(S) for S in Sigcol])

            for n in xrange(N + 1):
                if Arow[n]:
                    L = Lrow[n]
                    x = Wrow[n] - Murow[n]
                    xs = scipy.linalg.solve_triangular(L, x.T, lower=True)
                    lps[m] += (
                        -1.0 / 2.0 * inner1d(xs.T, xs.T) - B / 2.0 * np.log(2 * np.pi) - np.log(L.diagonal()).sum()
                    )

                if n < N and Acol[n]:
                    L = Lcol[n]
                    x = Wcol[n] - Mucol[n]
                    xs = scipy.linalg.solve_triangular(L, x.T, lower=True)
                    lps[m] += (
                        -1.0 / 2.0 * inner1d(xs.T, xs.T) - B / 2.0 * np.log(2 * np.pi) - np.log(L.diagonal()).sum()
                    )

        # Compute average log probability
        lp = -np.log(M) + logsumexp(lps)
        return lp

def multivariate_t_loglik(y,nu,mu,lmbda):
    # returns the log value
    d = len(mu)
    yc = np.array(y-mu,ndmin=2)
    ys, LT = general.solve_chofactor_system(lmbda,yc.T,overwrite_b=True)
    return scipy.special.gammaln((nu+d)/2.) - scipy.special.gammaln(nu/2.) \
            - (d/2.)*np.log(nu*np.pi) - np.log(LT.diagonal()).sum() \
            - (nu+d)/2.*np.log1p(1./nu*inner1d(ys.T,ys.T))

def get_lambda_max(X,y):
    """ 
    Find the value of lambda at which all coefficients are set to zero
    by finding the minimum value such that 0 is in the subdifferential
    and the coefficients are all zero.
    """
    subgrads = np.fabs( inner1d(X.T, y))
    return np.max( subgrads )

def multivariate_t_loglik(y,nu,mu,lmbda):
    # returns the log value
    d = len(mu)
    yc = np.array(y-mu,ndmin=2)
    L = np.linalg.cholesky(lmbda)
    ys = scipy.linalg.solve_triangular(L,yc.T,overwrite_b=True,lower=True)
    return scipy.special.gammaln((nu+d)/2.) - scipy.special.gammaln(nu/2.) \
            - (d/2.)*np.log(nu*np.pi) - np.log(L.diagonal()).sum() \
            - (nu+d)/2.*np.log1p(1./nu*inner1d(ys.T,ys.T))