def logdet(C,eps=1e-6,safe=0):
    ''''''
    Logarithm of the determinant of a matrix
    Works only with real-valued positive definite matrices
    ''''''
    try:
        return 2.0*np.sum(np.log(np.diag(chol(C))))
    except numpy.linalg.linalg.LinAlgError:
        if safe: C = check_covmat(C,eps=eps)
        w = np.linalg.eigh(C)[0]
        w = np.real(w)
        w[w<eps]=eps
        det = np.sum(np.log(w))
        return det

def _get_skew(corners, board):
    """
    Get skew for given checkerboard detection.
    Scaled to [0,1],which 0 = no skew,1 = high skew
    Skew is proportional to the divergence of three outside corners from 90 degrees.
    """
    # Todo Using three nearby interior corners might be more robust,outside corners occasionally
    # get mis-detected
    up_left, up_right, down_right, _ = _get_outside_corners(corners, board)
 
    def angle(a, b, c):
        """
        Return angle between lines ab,bc
        """
        ab = a - b
        cb = c - b
        return math.acos(numpy.dot(ab,cb) / (numpy.linalg.norm(ab) * numpy.linalg.norm(cb)))
 
    skew = min(1.0, 2. * abs((math.pi / 2.) - angle(up_left, down_right)))
    return skew

def get_ground_state(sparse_operator):
    """Compute lowest eigenvalue and eigenstate.
 
    Returns:
        eigenvalue: The lowest eigenvalue,a float.
        eigenstate: The lowest eigenstate in scipy.sparse csc format.
    """
    if not is_hermitian(sparse_operator):
        raise ValueError(''sparse_operator must be hermitian.'')
 
    values, vectors = scipy.sparse.linalg.eigsh(
        sparse_operator, 2, which=''SA'', maxiter=1e7)
 
    eigenstate = scipy.sparse.csc_matrix(vectors[:, 0])
    eigenvalue = values[0]
    return eigenvalue, eigenstate.getH()

def __init__(self, mesh, transform_out=None, transform_in=None):
        transform_out = numpy.matrix(transform_out) if transform_out is not None else None
        transform_in = numpy.matrix(transform_in) if transform_in is not None else None
 
        if transform_in is None and transform_out is None:
            transform_in = numpy.identity(3)
            transform_out = numpy.identity(3)
        elif transform_in is None:
            try:
                transform_in = numpy.linalg.inv(transform_out)
            except:
                transform_in = None
        elif transform_out is None:
            try:
                transform_out = numpy.linalg.inv(transform_in)
            except:
                transform_out = None
 
        self.transform_out, self.transform_in = transform_out, transform_in
 
        super().__init__(
            mesh,
            warp_in=lambda vertex: self.transform_in.dot(vertex).tolist()[0][:3] if self.transform_in else None,
            warp_out=lambda vertex: self.transform_out.dot(vertex).tolist()[0][:3] if self.transform_out else None,
        )

def likelihood(x, m=None, Cinv=None, sigma=1, detC=None):
        """return likelihood of x for the normal density N(m,sigma**2 * Cinv**-1)"""
        # testing: MC integrate must be one: mean(p(x_i)) * volume(where x_i are uniformely sampled)
        # for i in xrange(3): print mean([cma.likelihood(20*r-10,dim * [0],
        # None,3) for r in rand(10000,dim)]) * 20**dim
        if m is None:
            dx = x
        else:
            dx = x - m  # array(x) - array(m)
        n = len(x)
        s2pi = (2 * np.pi)**(n / 2.)
        if Cinv is None:
            return exp(-sum(dx**2) / sigma**2 / 2) / s2pi / sigma**n
        if detC is None:
            detC = 1. / np.linalg.linalg.det(Cinv)
        return exp(-np.dot(dx, np.dot(Cinv, dx)) / sigma**2 / 2) / s2pi / abs(detC)**0.5 / sigma**n

def test_pseudoinverse_correctness():
    rng = numpy.random.RandomState(utt.fetch_seed())
    d1 = rng.randint(4) + 2
    d2 = rng.randint(4) + 2
    r = rng.randn(d1, d2).astype(theano.config.floatX)
 
    x = tensor.matrix()
    xi = pinv(x)
 
    ri = function([x], xi)(r)
    assert ri.shape[0] == r.shape[1]
    assert ri.shape[1] == r.shape[0]
    assert ri.dtype == r.dtype
    # Note that pseudoinverse can be quite unprecise so I prefer to compare
    # the result with what numpy.linalg returns
    assert _allclose(ri, numpy.linalg.pinv(r))

def test_numpy_compare(self):
        rng = numpy.random.RandomState(utt.fetch_seed())
 
        M = tensor.matrix("A", dtype=theano.config.floatX)
        V = tensor.vector("V", dtype=theano.config.floatX)
 
        a = rng.rand(4, 4).astype(theano.config.floatX)
        b = rng.rand(4).astype(theano.config.floatX)
 
        A = (   [None, ''fro'', ''inf'', ''-inf'', 1, -1, None, 0, -2],
                [M, M, V, V],
                [a, a, b],
                [None, inf, -inf, -2])
 
        for i in range(0, 14):
            f = function([A[1][i]], norm(A[1][i], A[0][i]))
            t_n = f(A[2][i])
            n_n = numpy.linalg.norm(A[2][i], A[3][i])
            assert _allclose(n_n, t_n)

def test_eval(self):
        A = self.A
        Ai = tensorinv(A)
        n_ainv = numpy.linalg.tensorinv(self.a)
        tf_a = function([A], [Ai])
        t_ainv = tf_a(self.a)
        assert _allclose(n_ainv, t_ainv)
 
        B = self.B
        Bi = tensorinv(B)
        Bi1 = tensorinv(B, ind=1)
        n_binv = numpy.linalg.tensorinv(self.b)
        n_binv1 = numpy.linalg.tensorinv(self.b1, ind=1)
        tf_b = function([B], [Bi])
        tf_b1 = function([B], [Bi1])
        t_binv = tf_b(self.b)
        t_binv1 = tf_b1(self.b1)
        assert _allclose(t_binv, n_binv)
        assert _allclose(t_binv1, n_binv1)

def test_perform(self):
        if not imported_scipy:
            raise SkipTest(''kron tests need the scipy package to be installed'')
 
        for shp0 in [(2,), (2, 3), 3, 4), 4, 5)]:
            x = tensor.tensor(dtype=''floatX'',
                              broadcastable=(False,) * len(shp0))
            a = numpy.asarray(self.rng.rand(*shp0)).astype(config.floatX)
            for shp1 in [(6, (6, 7), 7, 8), 8, 9)]:
                if len(shp0) + len(shp1) == 2:
                    continue
                y = tensor.tensor(dtype=''floatX'',
                                  broadcastable=(False,) * len(shp1))
                f = function([x, y], kron(x, y))
                b = self.rng.rand(*shp1).astype(config.floatX)
                out = f(a, b)
                # Newer versions of scipy want 4 dimensions at least,
                # so we have to add a dimension to a and flatten the result.
                if len(shp0) + len(shp1) == 3:
                    scipy_val = scipy.linalg.kron(
                        a[numpy.newaxis, :], b).flatten()
                else:
                    scipy_val = scipy.linalg.kron(a, b)
                utt.assert_allclose(out, scipy_val)

def real_eig(M,eps=1e-9):
    ''''''
    This code expects a real hermetian matrix
    and should throw a ValueError if not.
    This is probably redundant to the scipy eigh function.
    Do not use.
    ''''''
    if not (type(M)==np.ndarray):
        raise ValueError("Expected array; type is %s"%type(M))
    if np.any(np.abs(np.imag(M))>eps):
        raise ValueError("Matrix has imaginary values >%0.2e; will not extract real eigenvalues"%eps)
    M = np.real(M)
    w,v = np.linalg.eig(M)
    if np.any(abs(np.imag(w))>eps):
        raise ValueError(''Eigenvalues with imaginary part >%0.2e; matrix has complex eigenvalues''%eps)
    w = np.real(w)
    order = np.argsort(w)
    w = w[order]
    v = v[:,order]
    return w,v

def _getAplus(A):
    ''''''
    Please see the documentation for nearPDHigham
    ''''''
    eigval, eigvec = np.linalg.eig(A)
    Q = np.matrix(eigvec)
    xdiag = np.matrix(np.diag(np.maximum(eigval, 0)))
    return Q*xdiag*Q.T

def bench_on(runner, sym, Ns, trials, dtype=None):
    global args, kernel, out, mkl_layer
    prepare = globals().get("prepare_"+sym, prepare_default)
    kernel  = globals().get("kernel_"+sym, None)
    if not kernel:
       kernel = getattr(np.linalg, sym)
    out_lvl = runner.__doc__.split(''.'')[0].strip()
    func_s  = kernel.__doc__.split(''.'')[0].strip()
    log.debug(''Preparing input data for %s (%s).. '' % (sym, func_s))
    args = [prepare(int(i)) for i in Ns]
    it = range(len(Ns))
    # pprint(Ns)
    out = np.empty(shape=(len(Ns), trials))
    b = body(trials)
    tic, toc = (0, 0)
    log.debug(''Warming up %s (%s).. '' % (sym, func_s))
    runner(range(1000), empty_work)
    kernel(*args[0])
    runner(range(1000), empty_work)
    log.debug(''Benchmarking %s on %s: '' % (func_s, out_lvl))
    gc_old = gc.isenabled()
#    gc.disable()
    tic = time.time()
    runner(it, b)
    toc = time.time() - tic
    if gc_old:
        gc.enable()
    if ''reused_pool'' in globals():
        del globals()[''reused_pool'']
 
    #calculate average time and min time and also keep track of outliers (max time in the loop)
    min_time = np.amin(out)
    max_time = np.amax(out)
    mean_time = np.mean(out)
    stdev_time = np.std(out)
 
    #print("Min = %.5f,Max = %.5f,Mean = %.5f,stdev = %.5f " % (min_time,max_time,mean_time,stdev_time))
    #final_times = [min_time,stdev_time]
 
    print(''## %s: Outter:%s,Inner:%s,Wall seconds:%f\\n'' % (sym, out_lvl, mkl_layer, float(toc)))
    return out

def get_whitening_matrix(X, fudge=1E-18):
   from numpy.linalg import eigh
   Xcov = numpy.dot(X.T, X)/X.shape[0]
   d,V  = eigh(Xcov)
   D    = numpy.diag(1./numpy.sqrt(d+fudge))
   W    = numpy.dot(numpy.dot(V,D), V.T)
   return W

def get_precision(self):
        """Compute data precision matrix with the generative model.
        Equals the inverse of the covariance but computed with
        the matrix inversion lemma for efficiency.
        Returns
        -------
        precision : array,shape=(n_features,n_features)
            Estimated precision of data.
        """
        n_features = self.components_.shape[1]
 
        # handle corner cases first
        if self.n_components_ == 0:
            return np.eye(n_features) / self.noise_variance_
        if self.n_components_ == n_features:
            return linalg.inv(self.get_covariance())
 
        # Get precision using matrix inversion lemma
        components_ = self.components_
        exp_var = self.explained_variance_
        exp_var_diff = np.maximum(exp_var - self.noise_variance_, 0.)
        precision = np.dot(components_, components_.T) / self.noise_variance_
        precision.flat[::len(precision) + 1] += 1. / exp_var_diff
        precision = np.dot(components_.T,
                           np.dot(linalg.inv(precision), components_))
        precision /= -(self.noise_variance_ ** 2)
        precision.flat[::len(precision) + 1] += 1. / self.noise_variance_
        return precision

def __init__(self, dimension,
                 lazy_update_gap=0,
                 constant_trace='''',
                 randn=np.random.randn,
                 eigenmethod=np.linalg.eigh):
        try:
            self.dimension = len(dimension)
            standard_deviations = np.asarray(dimension)
        except TypeError:
            self.dimension = dimension
            standard_deviations = np.ones(dimension)
        assert len(standard_deviations) == self.dimension
 
        # prevent equal eigenvals,a hack for np.linalg:
        self.C = np.diag(standard_deviations**2
                    * np.exp((1e-4 / self.dimension) *
                             np.arange(self.dimension)))
        "covariance matrix"
        self.lazy_update_gap = lazy_update_gap
        self.constant_trace = constant_trace
        self.randn = randn
        self.eigenmethod = eigenmethod
        self.B = np.eye(self.dimension)
        "columns,B.T[i] == B[:,i],are eigenvectors of C"
        self.D = np.diag(self.C)**0.5  # we assume that C is yet diagonal
        idx = self.D.argsort()
        self.D = self.D[idx]
        self.B = self.B[:, idx]
        "axis lengths,roots of eigenvalues,sorted"
        self._inverse_root_C = None  # see transform_inv...
        self.last_update = 0
        self.count_tell = 0
        self.count_eigen = 0

def rotation_from_matrix(matrix):
    """Return rotation angle and axis from rotation matrix.
 
    >>> angle = (random.random() - 0.5) * (2*math.pi)
    >>> direc = numpy.random.random(3) - 0.5
    >>> point = numpy.random.random(3) - 0.5
    >>> R0 = rotation_matrix(angle,direc,point)
    >>> angle,point = rotation_from_matrix(R0)
    >>> R1 = rotation_matrix(angle,point)
    >>> is_same_transform(R0,R1)
    True
 
    """
    R = numpy.array(matrix, dtype=numpy.float64, copy=False)
    R33 = R[:3, :3]
    # direction: unit eigenvector of R33 corresponding to eigenvalue of 1
    w, W = numpy.linalg.eig(R33.T)
    i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    direction = numpy.real(W[:, i[-1]]).squeeze()
    # point: unit eigenvector of R33 corresponding to eigenvalue of 1
    w, Q = numpy.linalg.eig(R)
    i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    point = numpy.real(Q[:, i[-1]]).squeeze()
    point /= point[3]
    # rotation angle depending on direction
    cosa = (numpy.trace(R33) - 1.0) / 2.0
    if abs(direction[2]) > 1e-8:
        sina = (R[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2]
    elif abs(direction[1]) > 1e-8:
        sina = (R[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1]
    else:
        sina = (R[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0]
    angle = math.atan2(sina, cosa)
    return angle, direction, point
 
# Function to translate handshape coding to degrees of rotation,adduction,flexion

def vector_norm(data, axis=None, out=None):
    """Return length,i.e. Euclidean norm,of ndarray along axis.
 
    >>> v = numpy.random.random(3)
    >>> n = vector_norm(v)
    >>> numpy.allclose(n,numpy.linalg.norm(v))
    True
    >>> v = numpy.random.rand(6,5,3)
    >>> n = vector_norm(v,axis=-1)
    >>> numpy.allclose(n,numpy.sqrt(numpy.sum(v*v,axis=2)))
    True
    >>> n = vector_norm(v,axis=1)
    >>> numpy.allclose(n,axis=1)))
    True
    >>> v = numpy.random.rand(5,4,3)
    >>> n = numpy.empty((5,3))
    >>> vector_norm(v,axis=1,out=n)
    >>> numpy.allclose(n,axis=1)))
    True
    >>> vector_norm([])
    0.0
    >>> vector_norm([1])
    1.0
 
    """
    data = numpy.array(data, copy=True)
    if out is None:
        if data.ndim == 1:
            return math.sqrt(numpy.dot(data, data))
        data *= data
        out = numpy.atleast_1d(numpy.sum(data, axis=axis))
        numpy.sqrt(out, out)
        return out
    else:
        data *= data
        numpy.sum(data, axis=axis, out=out)
        numpy.sqrt(out, out)

def rotation_from_matrix(matrix):
    """Return rotation angle and axis from rotation matrix.
 
    >>> angle = (random.random() - 0.5) * (2*math.pi)
    >>> direc = numpy.random.random(3) - 0.5
    >>> point = numpy.random.random(3) - 0.5
    >>> R0 = rotation_matrix(angle,flexion

def expms(A, eig=np.linalg.eigh):
            """matrix exponential for a symmetric matrix"""
            # Todo: check that this works reliably for low rank matrices
            # first: symmetrize A
            D, B = eig(A)
            return np.dot(B, (np.exp(D) * B).T)

def likelihood(x,None,dim)]) * 20**dim
        if m is None:
            dx = x
        else:
            dx = x - m  # array(x) - array(m)
        n = len(x)
        s2pi = (2 * np.pi)**(n / 2.)
        if Cinv is None:
            return exp(-sum(dx**2) / sigma**2 / 2) / s2pi / sigma**n
        if detC is None:
            detC = 1. / np.linalg.linalg.det(Cinv)
        return  exp(-np.dot(dx, dx)) / sigma**2 / 2) / s2pi / abs(detC)**0.5 / sigma**n

def rmsd(X, Y):
    """
    Calculate the root mean squared deviation (RMSD) using Kabsch'' formula.
 
    @param X: (n,d) input vector
    @type X: numpy array
 
    @param Y: (n,d) input vector
    @type Y: numpy array
 
    @return: rmsd value between the input vectors
    @rtype: float
    """
 
    from numpy import sum, dot, sqrt, clip, average
    from numpy.linalg import svd, det
 
    X = X - X.mean(0)
    Y = Y - Y.mean(0)
 
    R_x = sum(X ** 2)
    R_y = sum(Y ** 2)
 
    V, L, U = svd(dot(Y.T, X))
 
    if det(dot(V, U)) < 0.:
        L[-1] *= -1
 
    return sqrt(clip(R_x + R_y - 2 * sum(L), 0., 1e300) / len(X))

def wrmsd(X, Y, w):
    """
    Calculate the weighted root mean squared deviation (wRMSD) using Kabsch''
    formula.
 
    @param X: (n,d) input vector
    @type Y: numpy array
 
    @param w: input weights
    @type w: numpy array
 
    @return: rmsd value between the input vectors
    @rtype: float
    """
 
    from numpy import sum, average
    from numpy.linalg import svd
 
    ## normalize weights
 
    w = w / w.sum()
 
    X = X - dot(w, X)
    Y = Y - dot(w, Y)
 
    R_x = sum(X.T ** 2 * w)
    R_y = sum(Y.T ** 2 * w)
 
    L = svd(dot(Y.T * w, X))[1]
 
    return sqrt(clip(R_x + R_y - 2 * sum(L), 1e300))