Python numpy.linalg 模块，solve() 实例源码

def test_0_size_k(self):
        # test zero multiple equation (K=0) case.
        class ArraySubclass(np.ndarray):
            pass
        a = np.arange(4).reshape(1, 2, 2)
        b = np.arange(6).reshape(3, 2, 1).view(ArraySubclass)

        expected = linalg.solve(a, b)[:, :, 0:0]
        result = linalg.solve(a, b[:, :, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        # test both zero.
        expected = linalg.solve(a, b)[:, 0:0, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

def TPS_test(self,N,sizex,sizey,sizez,control_num,show_times,stop_time,typeofT,colors):
        for i in range(show_times):
            # source control points
            x, y, z = np.meshgrid(np.linspace(-1, 1, control_num),np.linspace(-1, 1, control_num), np.linspace(-1, 1, control_num))
            xs = x.flatten();ys = y.flatten();zs = z.flatten()
            cps = np.vstack([xs, ys, zs]).T
            # target control points
            xt = xs + np.random.uniform(-0.3, 0.3, size=xs.size);yt = ys + np.random.uniform(-0.3, 0.3, size=ys.size); zt = zs + np.random.uniform(-0.3, 0.3, size=zs.size)
            # construct T
            T = self.makeT(cps)
            # solve cx, cy (coefficients for x and y)
            xtAug = np.concatenate([xt, np.zeros(4)]);ytAug = np.concatenate([yt, np.zeros(4)]);ztAug = np.concatenate([zt, np.zeros(4)])
            cx = nl.solve(T, xtAug); cy = nl.solve(T, ytAug); cz = nl.solve(T, ztAug)
            # dense grid
            x = np.linspace(-sizex, sizex, 2*sizex); y = np.linspace(-sizey, sizey, 2*sizey); z = np.linspace(-sizez, sizez, 2*sizez)
            x, y, z = np.meshgrid(x, y, z)
            xgs, ygs, zgs = x.flatten(), y.flatten(), z.flatten()
            gps = np.vstack([xgs, ygs, zgs]).T
            # transform
            pgLift = self.liftPts(gps, cps) # [N x (K+3)]
            xgt = np.dot(pgLift, cx.T);ygt = np.dot(pgLift, cy.T); zgt = np.dot(pgLift,cz.T)
            # display
            showIm = ShowImage()
            showIm.Show_transform(xgs,ygs,zgs,xgt,ygt,zgt,sizex,sizey,sizez,stop_time,typeofT,N,colors)

def test_0_size_k(self):
        # test zero multiple equation (K=0) case.
        class ArraySubclass(np.ndarray):
            pass
        a = np.arange(4).reshape(1, 2, 2)
        b = np.arange(6).reshape(3, 2, 1).view(ArraySubclass)

        expected = linalg.solve(a, b)[:, :, 0:0]
        result = linalg.solve(a, b[:, :, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        # test both zero.
        expected = linalg.solve(a, b)[:, 0:0, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

def _single_sample_update(self, X, i, w):
        n_features = X.shape[1]
        if X.indptr[i + 1] - X.indptr[i] != 0:
            X_subset = X.data[X.indptr[i]:X.indptr[i + 1]]
            subset = X.indices[X.indptr[i]:X.indptr[i + 1]]
            len_subset = subset.shape[0]
            reduction = n_features / len_subset
            self.feature_n_iter_[subset] += 1
            components_subset = self.components_[:, subset]
            Dx = components_subset.dot(X_subset)
            G = components_subset.dot(components_subset.T)
            G.flat[::self.n_components + 1] += self.alpha / reduction
            self.code_[i] = linalg.solve(G, Dx)
            code = self.code_[i]
            w_B = np.minimum(1,
                             w * self.n_iter_ / self.feature_n_iter_[subset])
            self.B_[:, subset] *= 1 - w_B
            self.B_[:, subset] += np.outer(code, X_subset * w_B)

def enet_regression_multi_gram_(G, Dx, X, code, l1_ratio, alpha,
                                positive):
    batch_size = code.shape[0]
    if l1_ratio == 0:
        n_components = G.shape[1]
        for i in range(batch_size):
            G.flat[::n_components + 1] += alpha
            code[i] = linalg.solve(G[i], Dx[i])
            G.flat[::n_components + 1] -= alpha
    else:
        # Unused but unfortunate API
        random_state = check_random_state(0)
        for i in range(batch_size):
            cd_fast.enet_coordinate_descent_gram(
                code[i],
                alpha * l1_ratio,
                alpha * (
                    1 - l1_ratio),
                G[i], Dx[i], X[i], 100, 1e-2,
                random_state,
                False, positive)
    return code

def enet_regression_single_gram_(G, Dx, X, code, code_l1_ratio, code_alpha,
                                 code_pos):
    batch_size = code.shape[0]
    if code_l1_ratio == 0:
        n_components = G.shape[0]
        G = G.copy()
        G.flat[::n_components + 1] += code_alpha
        code[:] = linalg.solve(G, Dx.T).T
        G.flat[::n_components + 1] -= code_alpha
    else:
        # Unused but unfortunate API
        random_state = check_random_state(0)
        for i in range(batch_size):
            cd_fast.enet_coordinate_descent_gram(
                code[i],
                code_alpha * code_l1_ratio,
                code_alpha * (
                    1 - code_l1_ratio),
                G, Dx[i], X[i], 100, 1e-2,
                random_state,
                False, code_pos)
    return code

def _uwham_obj_grad(self, f_i):
        """Return the log-likelihood (scalar objective function) and its
        gradient (wrt the free energies) for a given value of the free
        energies.  Note that this expects one less free energy than there are
        states, since we always solve under the constraint that f_1 = 0.
        """
        _f_i = hstack((zeros(1), asarray(f_i)))
        # For numerical stability, use log quantities.
        logQ_nj = _f_i - self._u_nj_m
        logNorm_n = logsumexp(logQ_nj, 1, self.PIsdiag[newaxis, :])
        W_nj = exp(logQ_nj - logNorm_n[:, newaxis])
        # Compute matrix products and sums (equivalent to multiplying by
        # appropriate vector of ones). Note that using dot() with ndarrays is
        # _much_ faster than multiplying matrix objects.
        PIW = self.PIs.dot(W_nj.T)
        WPI = W_nj.dot(self.PIs)
        g = PIW.mean(axis=1) # used in gradient and Hessian computation
        kappa = logNorm_n.mean() - (self.PIsdiag.dot(_f_i)).sum()
        grad = (g - self.PIsdiag)[1:]
        self._hess = (diagflat(g) - PIW.dot(WPI) / self.total_samples)[1:, 1:]
        return kappa, grad

def solve(self):
        # known and unknown values
        self.VU = [node[key] for node in self.U.values() for key in ("ux",)]
        self.VF = [node[key] for node in self.F.values() for key in ("fx",)]
        knw = [pos for pos,value in enumerate(self.VU) if not value is np.nan]
        unknw = [pos for pos,value in enumerate(self.VU) if value is np.nan]
        # Matrices to solve
        self.K2S = np.delete(np.delete(self.KG,knw,0),knw,1)
        self.F2S = np.delete(self.VF,knw,0)
        # For displacements
        self.solved_u = la.solve(self.K2S,self.F2S)
        # Updating U (displacements vector)
        for k,ic in enumerate(unknw):
            nd, var = self.index2key(ic)
            self.U[nd][var] = self.solved_u[k]
            self.nodes[ic].ux = self.solved_u[k]
        # For nodal forces/reactions
        self.NF = self.F.copy()
        self.VU = [node[key] for node in self.U.values() for key in ("ux",)]
        nf_calc = np.dot(self.KG, self.VU)
        for k,ic in enumerate(range(self.get_number_of_nodes())):
            nd, var = self.index2key(ic, ("fx",))
            self.NF[nd][var] = nf_calc[k]
            self.nodes[ic].fx = nf_calc[k]

def getBarycentricCoords(A, B, C, X, checkValidity = True):
    T = np.array( [ [A.x - C.x, B.x - C.x ], [A.y - C.y, B.y - C.y] ] )
    y = np.array( [ [X.x - C.x], [X.y - C.y] ] )
    lambdas = linalg.solve(T, y)
    lambdas = lambdas.flatten()
    lambdas = np.append(lambdas, 1 - (lambdas[0] + lambdas[1]))
    if checkValidity:
        if (lambdas[0] < 0 or lambdas[1] < 0 or lambdas[2] < 0):
            print "ERROR: Not a convex combination; lambda = %s"%lambdas
            print "pointInsideConvexPolygon2D = %s"%pointInsideConvexPolygon2D([A, B, C], X, 0)
            plt.plot([A.x, B.x, C.x, A.x], [A.y, B.y, C.y, A.y], 'r')
            plt.hold(True)
            plt.plot([X.x], [X.y], 'b.')
            plt.show()
        assert (lambdas[0] >= 0 and lambdas[1] >= 0 and lambdas[2] >= 0)
    else:
        lambdas[0] = max(lambdas[0], 0)
        lambdas[1] = max(lambdas[1], 0)
        lambdas[2] = max(lambdas[2], 0)
    return lambdas

def get_mvdr_vector(atf_vector, noise_psd_matrix):
    """
    Returns the MVDR beamforming vector.

    :param atf_vector: Acoustic transfer function vector
        with shape (..., bins, sensors)
    :param noise_psd_matrix: Noise PSD matrix
        with shape (bins, sensors, sensors)
    :return: Set of beamforming vectors with shape (..., bins, sensors)
    """

    while atf_vector.ndim > noise_psd_matrix.ndim - 1:
        noise_psd_matrix = np.expand_dims(noise_psd_matrix, axis=0)

    # Make sure matrix is hermitian
    noise_psd_matrix = 0.5 * (
        noise_psd_matrix + np.conj(noise_psd_matrix.swapaxes(-1, -2)))

    numerator = solve(noise_psd_matrix, atf_vector)
    denominator = np.einsum('...d,...d->...', atf_vector.conj(), numerator)
    beamforming_vector = numerator / np.expand_dims(denominator, axis=-1)

    return beamforming_vector

def get_mvdr_vector(atf_vector, noise_psd_matrix):
    """
    Returns the MVDR beamforming vector: h = (Npsd^-1)*A / (A^H*(Npsd^-1)A)
    :param atf_vector: Acoustic transfer function vector
        with shape (..., bins, sensors)
    :param noise_psd_matrix: Noise PSD matrix
        with shape (bins, sensors, sensors)
    :return: Set of beamforming vectors with shape (..., bins, sensors)
    """
    while atf_vector.ndim > noise_psd_matrix.ndim - 1:
        noise_psd_matrix = np.expand_dims(noise_psd_matrix, axis=0)

    # Make sure matrix is hermitian
    noise_psd_matrix = 0.5 * (
        noise_psd_matrix + np.conj(noise_psd_matrix.swapaxes(-1, -2)))

    numerator = solve(noise_psd_matrix, atf_vector)
    denominator = np.einsum('...d,...d->...', atf_vector.conj(), numerator)
    beamforming_vector = numerator / np.expand_dims(denominator, axis=-1)

    return beamforming_vector

def apply_sdw_mwf(mix, target_psd_matrix, noise_psd_matrix, mu=1, corr=None):
    """
    Apply speech distortion weighted MWF: h = Tpsd * e1 / (Tpsd + mu*Npsd) 
    :param mix: the signal complex FFT
    :param target_psd_matrix (bins, sensors, sensors) 
    :param noise_psd_matrix
    :param mu: the lagrange factor
    :return 
    """
    bins, sensors, frames = mix.shape
    ref_vector = np.zeros((sensors,1), dtype=np.float)    
    if corr is None:
        ref_ch = 0
    else: # choose the channel with highest correlation with the others
        corr=corr.tolist()        
        while len(corr) > sensors:
            corr.remove(np.min(corr))
        ref_ch=np.argmax(corr)
    ref_vector[ref_ch,0]=1 

    mwf_vector = solve(target_psd_matrix + mu*noise_psd_matrix, target_psd_matrix[:,:,ref_ch])
    return np.einsum('...a,...at->...t', mwf_vector.conj(), mix)

def test_0_size_k(self):
        # test zero multiple equation (K=0) case.
        class ArraySubclass(np.ndarray):
            pass
        a = np.arange(4).reshape(1, 2, 2)
        b = np.arange(6).reshape(3, 2, 1).view(ArraySubclass)

        expected = linalg.solve(a, b)[:, :, 0:0]
        result = linalg.solve(a, b[:, :, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        # test both zero.
        expected = linalg.solve(a, b)[:, 0:0, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

def test_0_size_k(self):
        # test zero multiple equation (K=0) case.
        class ArraySubclass(np.ndarray):
            pass
        a = np.arange(4).reshape(1, 2, 2)
        b = np.arange(6).reshape(3, 2, 1).view(ArraySubclass)

        expected = linalg.solve(a, b)[:, :, 0:0]
        result = linalg.solve(a, b[:, :, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        # test both zero.
        expected = linalg.solve(a, b)[:, 0:0, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

def get_mvdr_vector(atf_vector, noise_psd_matrix):
    """
    Returns the MVDR beamforming vector.

    :param atf_vector: Acoustic transfer function vector
        with shape (..., bins, sensors)
    :param noise_psd_matrix: Noise PSD matrix
        with shape (bins, sensors, sensors)
    :return: Set of beamforming vectors with shape (..., bins, sensors)
    """

    while atf_vector.ndim > noise_psd_matrix.ndim - 1:
        noise_psd_matrix = np.expand_dims(noise_psd_matrix, axis=0)

    # Make sure matrix is hermitian
    noise_psd_matrix = 0.5 * (
        noise_psd_matrix + np.conj(noise_psd_matrix.swapaxes(-1, -2)))

    numerator = solve(noise_psd_matrix, atf_vector)
    denominator = np.einsum('...d,...d->...', atf_vector.conj(), numerator)
    beamforming_vector = numerator / np.expand_dims(denominator, axis=-1)

    return beamforming_vector

def test_0_size_k(self):
        # test zero multiple equation (K=0) case.
        class ArraySubclass(np.ndarray):
            pass
        a = np.arange(4).reshape(1, 2, 2)
        b = np.arange(6).reshape(3, 2, 1).view(ArraySubclass)

        expected = linalg.solve(a, b)[:, :, 0:0]
        result = linalg.solve(a, b[:, :, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        # test both zero.
        expected = linalg.solve(a, b)[:, 0:0, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

def getBarycentricCoords(A, B, C, X, checkValidity = True):
    T = np.array( [ [A.x - C.x, B.x - C.x ], [A.y - C.y, B.y - C.y] ] )
    y = np.array( [ [X.x - C.x], [X.y - C.y] ] )
    lambdas = linalg.solve(T, y)
    lambdas = lambdas.flatten()
    lambdas = np.append(lambdas, 1 - (lambdas[0] + lambdas[1]))
    if checkValidity:
        if (lambdas[0] < 0 or lambdas[1] < 0 or lambdas[2] < 0):
            print "ERROR: Not a convex combination; lambda = %s"%lambdas
            print "pointInsideConvexPolygon2D = %s"%pointInsideConvexPolygon2D([A, B, C], X, 0)
            plt.plot([A.x, B.x, C.x, A.x], [A.y, B.y, C.y, A.y], 'r')
            plt.hold(True)
            plt.plot([X.x], [X.y], 'b.')
            plt.show()
        assert (lambdas[0] >= 0 and lambdas[1] >= 0 and lambdas[2] >= 0)
    else:
        lambdas[0] = max(lambdas[0], 0)
        lambdas[1] = max(lambdas[1], 0)
        lambdas[2] = max(lambdas[2], 0)
    return lambdas

def update_V(m_opts, m_vars):
    P, N = E_x_omega_col(m_opts, m_vars) # expectation of xi_{nl} for n = i, expecation of omega_{nl} for n = i
    Kappa = PG_col(m_opts, m_vars) # polyagamma kappa_{nl} for n = i
    PN = P*N
    PK = P*Kappa

    for i in range(m_vars['n_labels']):
        PN_i = PN[i][:,np.newaxis]
        PK_i = PK[i]

        sigma = m_vars['U_batch'].T.dot(PN_i*m_vars['U_batch'])# + m_opts['lam_v']*np.eye(m_opts['n_components'])
        x = m_vars['U_batch'].T.dot(PK_i)

        m_vars['sigma_v'][i] = (1-m_vars['gamma'])*m_vars['sigma_v'][i] + m_vars['gamma']*sigma
        m_vars['x_v'][i] = (1-m_vars['gamma'])*m_vars['x_v'][i] + m_vars['gamma']*x
        m_vars['V'][i] = linalg.solve(m_vars['sigma_v'][i], m_vars['x_v'][i])

def test_0_size_k(self):
        # test zero multiple equation (K=0) case.
        class ArraySubclass(np.ndarray):
            pass
        a = np.arange(4).reshape(1, 2, 2)
        b = np.arange(6).reshape(3, 2, 1).view(ArraySubclass)

        expected = linalg.solve(a, b)[:, :, 0:0]
        result = linalg.solve(a, b[:, :, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        # test both zero.
        expected = linalg.solve(a, b)[:, 0:0, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

def solve(self, A, b):
        epsilon = self.epsilon
        force = self.force
        while True:
            try:
                rtn = super(InsufficientRankSolver, self).solve(A, b)
                break;
            except LinAlgError as e:
                if epsilon is None and force is None:
                    raise RuntimeError("Did not provide a way to resolve sigular matrix")

            if epsilon is not None:
                diagval = np.sum(np.diagonal(A))
                E = np.ones(A.shape, A.dtype)*epsilon
                if diagval==0:
                    E-=np.diag(np.diag(E))
                A+=E
                epsilon = None
            else:
                A[-1,:] = np.ones(A.shape[0], A.dtype)
                b[-1] = 0
                force = None

        return rtn;

def test_0_size_k(self):
        # test zero multiple equation (K=0) case.
        class ArraySubclass(np.ndarray):
            pass
        a = np.arange(4).reshape(1, 2, 2)
        b = np.arange(6).reshape(3, 2, 1).view(ArraySubclass)

        expected = linalg.solve(a, b)[:, :, 0:0]
        result = linalg.solve(a, b[:, :, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        # test both zero.
        expected = linalg.solve(a, b)[:, 0:0, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

def gaussian2D(window):
    """Real 2D Gaussian interpolation for sub pixel shift"""
    #ref on paper
    w = np.ones((3, 3))*(1./9)
    rhs = np.zeros(6)
    M = np.zeros((6,6))
    for i in [-1, 0, 1]:
        for j in [-1, 0, 1]:
            rhs = rhs     +     np.array([i*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                          j*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                        i*j*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                        i*i*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                        j*j*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                            w[i+1, j+1]*np.log(np.abs(window[i+1, j+1]))], dtype='float')

            M = M + w[i+1, j+1]*np.array([[  i*i,   i*j,   i*i*j,   i*i*i,   i*j*j,   i],
                                          [  i*j,   j*j,   i*j*j,   i*i*j,   j*j*j,   j],
                                          [i*i*j, i*j*j, i*i*j*j, i*i*i*j, i*j*j*j, i*j],
                                          [i*i*i, i*i*j, i*i*i*j, i*i*i*i, i*i*j*j, i*i],
                                          [i*j*j, j*j*j, i*j*j*j, i*i*j*j, j*j*j*j, j*j],
                                          [    i,     j,     i*j,     i*i,     j*j,   1]], dtype='float')
    solution = nl.solve(M, rhs)

    dx = (    solution[2]*solution[1] - 2.0*solution[0]*solution[4])/ \
         (4.0*solution[3]*solution[4] -     solution[2]*solution[2])

    dy = (    solution[2]*solution[0] - 2.0*solution[1]*solution[3])/ \
         (4.0*solution[3]*solution[4] -     solution[2]*solution[2])

    return dx, dy

def mittelpunkt(A,B,C,D):
    ''' schnitt der Diagonalen AC und BD '''
    M=np.array([A-C,D-B])
    R=D-C
    x=npl.solve(M.T,R)
    (np.dot(M,x)-R)
    S=x[0]*A+(1-x[0])*C
    return S

def mittelpunkt(A,B,C,D):
    ''' schnitt der Diagonalen AC und BD '''
    M=np.array([A-C,D-B])
    R=D-C
    x=npl.solve(M.T,R)
    (np.dot(M,x)-R)
    S=x[0]*A+(1-x[0])*C
    return S

def do(self, a, b):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(imply(isinstance(b, matrix), isinstance(x, matrix)))

def test_types(self):
        def check(dtype):
            x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
            assert_equal(linalg.solve(x, x).dtype, dtype)
        for dtype in [single, double, csingle, cdouble]:
            yield check, dtype

def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

        expected = linalg.solve(a, b)[:, 0:0, :]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0, :])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2)  # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b)

def line_segments_intersect(p1, p2, q1, q2, magn_delta_p):
    # solve equation line1 = line2
    # p1 + lambda * (p2-p1) = q1 + my * (q2-q1)
    # lambda * (p2-p1) - my * (q2-q1) = q1 - p1
    # normalize deltas to 1
    # lambda * (p2-p1)/|p1-p2| - my * (q2-q1)/|q1-q2| = q1 - p1
    # magn_delta_p = |p1-p2|  (magnitude)
    # magn_delta_q = euclidean_distance(q1[0], q1[1], q2[0], q2[1])
    if max(p1[0], p2[0]) < min(q1[0], q2[0]) or max(q1[0], q2[0]) < min(p1[0], p2[0]) or \
            max(p1[1], p2[1]) < min(q1[1], q2[1]) or max(q1[1], q2[1]) < min(p1[1], p2[1]):
        return False

    dif_p_x = p2[0] - p1[0]
    dif_p_y = p2[1] - p1[1]
    dif_q_x = q2[0] - q1[0]
    dif_q_y = q2[1] - q1[1]
    a = array([[dif_p_x, -dif_q_x], [dif_p_y, -dif_q_y]])
    # [[dif_p_x / magn_delta_p, -dif_q_x / magn_delta_q], [dif_p_y / magn_delta_p, -dif_q_y / magn_delta_q]])
    b = array([q1[0] - p1[0], q1[1] - p1[1]])
    try:
        x = linalg.solve(a, b)
    except linalg.linalg.LinAlgError:
        # Singular matrix (lines parallel, there is not intersection)
        return False

    if x[0] < 0 or x[0] > 1 or x[1] < 0 or x[1] > 1:
        # intersection of the two lines appears before or after actual line segments
        # in this use case it is important to include the points themselves when checking for intersections
        # that way a new edge that barely touches the old polygon is also legal
        return False

    return True

def AffineFromPoints(src1,src2,dst1,dst2,new_size,filter=BILINEAR):
    ''' 
    An affine transform that will rotate, translate, and scale to map one 
    set of points to the other. For example, to align eye coordinates in face images.

    Find a transform (a,b,tx,ty) such that it maps the source points to the 
    destination points::

        a*x1-b*y1+tx = x2
        b*x1+a*y1+ty = y2

    The mapping between the two points creates a set of  four linear equations 
    with four unknowns. This set of equations is solved to find the transform.

    @param src1: the first link.Point in the source image.
    @param src2: the second link.Point in the source image.
    @param dst1: the first link.Point in the destination image.
    @param dst2: the second link.Point in the destination image.
    @param new_size: new size for the image.
    @param filter: PIL filter to use.
    '''

    # Compute the transformation parameters
    A = [[src1.X(),-src1.Y(),1,0],
         [src1.Y(),src1.X(),0,1],
         [src2.X(),-src2.Y(),1,0],
         [src2.Y(),src2.X(),0,1]]
    b = [dst1.X(),dst1.Y(),dst2.X(),dst2.Y()]
    A = array(A)
    b = array(b)
    result = solve(A,b)

    a,b,tx,ty = result    
    # Create the transform matrix
    matrix = array([[a,-b,tx],[b,a,ty],[0,0,1]],'d')

    return AffineTransform(matrix,new_size,filter)

def _solve_theta(lo, hi, topic, X, BTBpR, c0, c1, f):
    theta_batch = np.empty((hi - lo, f), dtype=topic.dtype)
    for ib, u in enumerate(range(lo, hi)):
        x_u, idx_u = get_row(X, u)
        B_u = topic[idx_u]
        a = x_u.dot(c1 * B_u)
        '''
        non-zero elements handled in this loop
        '''
        B = BTBpR + B_u.T.dot((c1 - c0) * B_u)#B_u only contains vectors corresponding to non-zero doc-word occurence
        theta_batch[ib] = LA.solve(B, a)
    theta_batch = theta_batch.clip(0)
    return theta_batch

def _solve_topic(lo, hi, theta, alpha, gamma, XT, BTBpR, c0, c1, f, lam_d, lam_t):
    topic_batch = np.empty((hi - lo, f), dtype=theta.dtype)
    for ib, u in enumerate(range(lo, hi)):
        x_u, idx_u = get_row(XT, u)
        B_u = theta[idx_u]
        cpAT = gamma[u].dot(alpha.T)
        a = lam_d * x_u.dot(c1 * B_u) + lam_t * cpAT
        '''
        non-zero elements handled in this loop
        '''
        B = BTBpR + B_u.T.dot((c1 - c0) * B_u)#B_u only contains vectors corresponding to non-zero doc-word occurence
        topic_batch[ib] = LA.solve(B, a)
    topic_batch = topic_batch.clip(0)
    return topic_batch

def _solve_gamma(lo, hi, alpha, beta, topic, M, B, c0, c1, f, lam_w, lam_t):
    gamma_batch = np.empty((hi - lo, f), dtype=alpha.dtype)
    for ib, i in enumerate(range(lo, hi)):
        t_i = topic[i,:]
        m_i, idx_m_i = get_row(M, i)
        B_i = beta[idx_m_i]
        '''
        the reason why they put G_i in the loop instead of calculate GTG = gamma.T * gamma is that in the objective function,
        we currently only consider the non-zero elements in matrix W.
        '''
        a = lam_t * np.dot(t_i, alpha) + lam_w * np.dot(m_i, B_i)
        gamma_batch[ib] = LA.solve(B, a)
    return gamma_batch

def _solve_beta(lo, hi, gamma, MT, B, f):
    beta_batch = np.empty((hi - lo, f), dtype=gamma.dtype)
    for ib, j in enumerate(range(lo, hi)):
        m_j, idx_m_j = get_row(MT, j)
        C_j = gamma[idx_m_j]
        a = np.dot(m_j, C_j)
        beta_batch[ib] = LA.solve(B, a)
    return beta_batch

def _solve_alpha(lo, hi, gamma, topicT, B, f):
    alpha_batch = np.empty((hi - lo, f), dtype=gamma.dtype)
    for ib, j in enumerate(range(lo, hi)):
        t_j = topicT[j,:]        
        a = np.dot(t_j, gamma)
        alpha_batch[ib] = LA.solve(B, a)
    return alpha_batch

def _estimate_gradient(self, policy):
        env = self.environment
        var = self.var
        # store current policy parameter
        parameter = policy.parameters
        par_dim = policy.parameter_space.dimension

        # using forward differences
        trace = env.rollout(policy)
        j_ref = sum([x[2] for x in trace]) / len(trace)

        dj = np.zeros((2 * par_dim))
        dv = np.append(np.eye(par_dim), -np.eye(par_dim), axis=0)
        dv *= var

        for n in range(par_dim):
            variation = dv[n]

            policy.parameters = parameter + variation
            trace_n = env.rollout(policy)

            jn = sum([x[2] for x in trace]) / len(trace_n)

            dj[n] = j_ref - jn

        grad = solve(dv.T.dot(dv), dv.T.dot(dj))

        # reset current policy parameter
        policy.parameters = parameter

        return grad

def _estimate_gradient(self, policy):
        env = self.environment

        parameter = policy.parameters
        par_dim = policy.parameter_space.dimension

        dj = np.zeros((par_dim,))
        dv = np.eye(par_dim) * self.var / 2

        for n in range(par_dim):
            variation = dv[n]

            policy.parameters = parameter + variation
            trace_n = env.rollout(policy)

            policy.parameters = parameter - variation
            trace_n_ref = env.rollout(policy)

            jn = sum([x[2] for x in trace_n]) / len(trace_n)
            jn_ref = sum([x[2] for x in trace_n_ref]) / len(trace_n_ref)

            dj[n] = jn - jn_ref

        grad = solve(dv.T.dot(dv), dv.T.dot(dj))
        policy.parameters = parameter

        return grad

def _omp(x, D, Gram, alpha, n_nonzero_coefs=None, tol=None):
    _, n_atoms = D.shape
    # the dict indexes of the atoms this datapoint uses
    Dx = np.array([]).astype(int)
    z = np.zeros(n_atoms)
    # the residual
    r = np.copy(x)
    i = 0
    if n_nonzero_coefs is not None:
        tol = 1e-10
        def cont_criterion():
            not_reached_sparsity = i < n_nonzero_coefs
            return (not_reached_sparsity and norm(r) > tol)
    else:
        cont_criterion = lambda: norm(r) >= tol

    while (cont_criterion()):

        # find the atom that correlates the
        # most with the residual
        k = np.argmax(np.abs(alpha))
        if k in Dx:
            break
        Dx = np.append(Dx, k)
        # solve the Least Squares problem
        # to find the coefs z
        DI = D[:, Dx]
        G = Gram[Dx, :][:, Dx]
        G = np.atleast_2d(G)
        try:
            G_inv = inv(G)
        except LinAlgError:
            print gram_singular_msg
            break

        z[Dx] = np.dot(G_inv, np.dot(D.T, x)[Dx])
        r = x - np.dot(D[:, Dx], z[Dx])
        alpha = np.dot(D.T, r)
        i += 1

    return z

def llc(X, D, knn=5):
    # the sparse coder introduced in
    # "Locality-constrained Linear Coding for Image Classification"

    n_samples = X.shape[1]
    n_atoms = D.shape[1]
    # has the distance of
    # each sample to each atom
    dist = np.zeros((n_samples, n_atoms))
    # calculate the distances
    for i in range(n_samples):
        for j in range(n_atoms):
            dist[i, j] = norm(X[:, i] - D[:, j])

    # has the indices of the atoms
    # that are nearest neighbour to each sample
    knn_idx = np.zeros((n_samples, knn)).astype(int)
    for i in xrange(n_samples):
        knn_idx[i, :] = np.argsort(dist[i, :])[:knn]
    # the sparse coding matrix
    Z = np.zeros((n_atoms, n_samples))
    II = np.eye(knn)
    beta = 1e-4
    b = np.ones(knn)
    for i in xrange(n_samples):
        idx = knn_idx[i, :]
        z = D.T[idx, :] - repmat(X.T[i, :], knn, 1)
        C = np.dot(z, z.T)
        C = C + II * beta * np.trace(C)
        # solve the linear system C*c=b
        c = solve(C, b)
        # enforce the constraint on the sparse codes
        # such that sum(c)=1
        c = c / float(np.sum(c))
        Z[idx, i] = c

    return Z

def do(self, a, b):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(imply(isinstance(b, matrix), isinstance(x, matrix)))

def test_types(self):
        def check(dtype):
            x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
            assert_equal(linalg.solve(x, x).dtype, dtype)
        for dtype in [single, double, csingle, cdouble]:
            yield check, dtype

def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

        expected = linalg.solve(a, b)[:, 0:0, :]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0, :])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2)  # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b)

def _refit(self, X):
        n_samples, n_features = X.shape
        for i in range(n_samples):
            X_subset = X.data[X.indptr[i]:X.indptr[i + 1]]
            subset = X.indices[X.indptr[i]:X.indptr[i + 1]]
            len_subset = subset.shape[0]
            reduction = n_features / len_subset
            components_subset = self.components_[:, subset]
            Dx = components_subset.dot(X_subset)
            G = components_subset.dot(components_subset.T)
            G.flat[::self.n_components + 1] += self.alpha / reduction
            self.code_[i] = linalg.solve(G, Dx)

def fvar(self):
        """Estimate the variances of the reduced free energies."""
        # Shorthand indices - notation similiar to Ref. 1
        n, m, mpk = self.total_samples, self.nstates_sampled, self.nstates
        k = mpk - m
        mask0, maskn0 = self.mask_zero, self.mask_nonzero
        _W_nj = self._W_nj
        O = _W_nj.T.dot(_W_nj) / n
        Os = O[:, :m]
        B1 = (hstack((Os.dot(self.PIs), zeros((mpk, k))))
              - identity(mpk))[1:, 1:]
        A1 = (O - Os.dot(self.PIs).dot(Os.T))[1:, 1:]
        V = solve(B1, A1).dot(inv(B1.T)) / n
#        if any(diagonal(V) < 0.):
#            D = _W_nj.T.dot(self._R_nj) / n
#            A1 = (O - D.dot(self.PIs).dot(D.T))[1:, 1:]
#            V = solve(B1, A1).dot(inv(B1.T)) / n
        # Unshuffle the state indices. Note that the variance of state 0 is
        # zero by construction and thus is omitted from V - since the user
        # selected state 0 may be different from the lowest indexed sample
        # state, re-adjust so that the actual state 0 has zero variance.
        #
        V_full = zeros((mpk, mpk))
        V_full[1:, 1:] = V
        var = zeros(mpk)
        var[maskn0] += diagonal(V_full)[:m]
        var[mask0] += diagonal(V_full)[m:]
        var += var[0] 
        var[0] = 0.0
        return var