Python scipy.linalg 模块，LinAlgError() 实例源码

def _calculate(self, X, y, categorical, metafeatures, helpers):
        import sklearn.lda
        if len(y.shape) == 1 or y.shape[1] == 1:
            kf = cross_validation.StratifiedKFold(y, n_folds=10)
        else:
            kf = cross_validation.KFold(y.shape[0], n_folds=10)

        accuracy = 0.
        try:
            for train, test in kf:
                lda = sklearn.lda.LDA()

                if len(y.shape) == 1 or y.shape[1] == 1:
                    lda.fit(X[train], y[train])
                else:
                    lda = OneVsRestClassifier(lda)
                    lda.fit(X[train], y[train])

                predictions = lda.predict(X[test])
                accuracy += sklearn.metrics.accuracy_score(predictions, y[test])
            return accuracy / 10
        except LinAlgError as e:
            self.logger.warning("LDA failed: %s Returned 0 instead!" % e)
            return np.NaN
        except ValueError as e:
            self.logger.warning("LDA failed: %s Returned 0 instead!" % e)
            return np.NaN

def _calculate_sparse(self, X, y, categorical, metafeatures, helpers):
        import sklearn.decomposition
        rs = np.random.RandomState(42)
        indices = np.arange(X.shape[0])
        # This is expensive, but necessary with scikit-learn 0.15
        Xt = X.astype(np.float64)
        for i in range(10):
            try:
                rs.shuffle(indices)
                truncated_svd = sklearn.decomposition.TruncatedSVD(
                        n_components=X.shape[1] - 1, random_state=i,
                        algorithm="randomized")
                truncated_svd.fit(Xt[indices])
                return truncated_svd
            except LinAlgError as e:
                pass
        self.logger.warning("Failed to compute a Truncated SVD")
        return None


# Maybe define some more...

def _covariance_factor(Sigma):
  # Assume it is positive-definite and try Cholesky decomposition.
  try:
    return Covariance_Cholesky(Sigma)
  except la.LinAlgError:
    pass

  # XXX In the past, we tried LU decomposition if, owing to
  # floating-point rounding error, the matrix is merely nonsingular,
  # not positive-definite.  However, empirically, that seemed to lead
  # to bad numerical results.  Until we have better numerical analysis
  # of the situation, let's try just falling back to least-squares
  # pseudoinverse approximation.

  # Otherwise, fall back to whatever heuristics scipy can manage.
  return Covariance_Loser(Sigma)

def optimize(self, optimizer=None, start=None, **kwargs):
        self._IN_OPTIMIZATION_ = True
        if self.mpi_comm is None:
            super(DeepGP, self).optimize(optimizer,start,**kwargs)
        elif self.mpi_comm.rank==self.mpi_root:
            super(DeepGP, self).optimize(optimizer,start,**kwargs)
            self.mpi_comm.Bcast(np.int32(-1),root=self.mpi_root)
        elif self.mpi_comm.rank!=self.mpi_root:
            x = self.optimizer_array.copy()
            flag = np.empty(1,dtype=np.int32)
            while True:
                self.mpi_comm.Bcast(flag,root=self.mpi_root)
                if flag==1:
                    try:
                        self.optimizer_array = x
                    except (LinAlgError, ZeroDivisionError, ValueError):
                        pass
                elif flag==-1:
                    break
                else:
                    self._IN_OPTIMIZATION_ = False
                    raise Exception("Unrecognizable flag for synchronization!")
        self._IN_OPTIMIZATION_ = False

def _log_multivariate_normal_density_full(X, means, covars, min_covar=1.e-7):
    """Log probability for full covariance matrices."""
    n_samples, n_dim = X.shape
    nmix = len(means)
    log_prob = np.empty((n_samples, nmix))
    for c, (mu, cv) in enumerate(zip(means, covars)):
        try:
            cv_chol = linalg.cholesky(cv, lower=True)
        except linalg.LinAlgError:
            # The model is most probably stuck in a component with too
            # few observations, we need to reinitialize this components
            try:
                cv_chol = linalg.cholesky(cv + min_covar * np.eye(n_dim),
                                          lower=True)
            except linalg.LinAlgError:
                raise ValueError("'covars' must be symmetric, "
                                 "positive-definite")

        cv_log_det = 2 * np.sum(np.log(np.diagonal(cv_chol)))
        cv_sol = linalg.solve_triangular(cv_chol, (X - mu).T, lower=True).T
        log_prob[:, c] = - .5 * (np.sum(cv_sol ** 2, axis=1) +
                                 n_dim * np.log(2 * np.pi) + cv_log_det)

    return log_prob

def check_pd(A, lower=True):
    """
    Checks if A is PD.
    If so returns True and Cholesky decomposition,
    otherwise returns False and None
    """
    try:
        return True, np.tril(cho_factor(A, lower=lower)[0])
    except LinAlgError as err:
        if 'not positive definite' in str(err):
            return False, None

def _calculate(self, X, y, categorical, metafeatures, helpers):
        import sklearn.decomposition
        pca = sklearn.decomposition.PCA(copy=True)
        rs = np.random.RandomState(42)
        indices = np.arange(X.shape[0])
        for i in range(10):
            try:
                rs.shuffle(indices)
                pca.fit(X[indices])
                return pca
            except LinAlgError as e:
                pass
        self.logger.warning("Failed to compute a Principle Component Analysis")
        return None

def make_from_matrix(self, cov_mat, check_singular=False):
        """
        Sets the covariance matrix manually.

        Parameters
        ----------

        **cov_mat** : numpy.matrix
            A *square*, *symmetric* (and usually *regular*) matrix.

        Keyword Arguments
        -----------------

        check_singular : boolean, optional
            Whether to force singularity check. Defaults to ``False``.
        """

        # Check matrix suitability
        _mat = np.asmatrix(cov_mat)  # cast to matrix
        _shp = _mat.shape
        if not _shp[0] == _shp[1]:   # check square shape
            raise ValueError("Failed to make ErrorSource: matrix must be "
                               "square, got shape %r" % (_shp,))
        if (_mat == _mat.T).all():  # check if symmetric
            if check_singular:
                try:
                    _mat.I  # try to invert matrix
                except LinAlgError:
                    raise ValueError("Failed to make ErrorSource: singular "
                                       "matrix!")
        else:
            raise ValueError("Failed to make ErrorSource: covariance "
                               "matrix not symmetric")

        self.error_type = 'matrix'
        self.error_value = _mat
        self.size = _shp[0]
        self.has_correlations = (np.diag(np.diag(_mat)) == _mat).all()

def solve(A,b,rtol=10**-8):
    '''
    Solve the matrix equation A*x=b by QR decomposition.

    Parameters
    ----------
    A : 2d ndarray
        The coefficient matrix.
    b : 1d ndarray
        The ordinate values.
    rtol : np.float64
        The relative tolerance of the solution.

    Returns
    -------
    1d ndarray
        The solution.

    Raises
    ------
    LinAlgError
        When no solution exists.
    '''
    assert A.ndim==2
    nrow,ncol=A.shape
    if nrow>=ncol:
        result=np.zeros(ncol,dtype=np.find_common_type([],[A.dtype,b.dtype]))
        q,r=sl.qr(A,mode='economic',check_finite=False)
        temp=q.T.dot(b)
        for i,ri in enumerate(r[::-1]):
            result[-1-i]=(temp[-1-i]-ri[ncol-i:].dot(result[ncol-i:]))/ri[-1-i]
    else:
        temp=np.zeros(nrow,dtype=np.find_common_type([],[A.dtype,b.dtype]))
        q,r=sl.qr(dagger(A),mode='economic',check_finite=False)
        for i,ri in enumerate(dagger(r)):
            temp[i]=(b[i]-ri[:i].dot(temp[:i]))/ri[i]
        result=q.dot(temp)
    if not np.allclose(A.dot(result),b,rtol=rtol):
        raise sl.LinAlgError('solve error: no solution.')
    return result

def lnlike_spec(spec_mu, obs=None, spec_noise=None, **vectors):
        """Calculate the likelihood of the spectroscopic data given the
        spectroscopic model.  Allows for the use of a gaussian process
        covariance matrix for multiplicative residuals.

        :param spec_mu:
            The mean model spectrum, in linear or logarithmic units, including
            e.g. calibration and sky emission.

        :param obs: (optional)
            A dictionary of the observational data, including the keys
            *``spectrum`` a numpy array of the observed spectrum, in linear or
             logarithmic units (same as ``spec_mu``).
            *``unc`` the uncertainty of same length as ``spectrum``
            *``mask`` optional boolean array of same length as ``spectrum``
            *``wavelength`` if using a GP, the metric that is used in the
             kernel generation, of same length as ``spectrum`` and typically
             giving the wavelength array.

        :param spec_noise: (optional)
            A NoiseModel object with the methods `compute` and `lnlikelihood`.
            If ``spec_noise`` is supplied, the `wavelength` entry in the obs
            dictionary must exist.

        :param vectors: (optional)
            A dictionary of vectors of same length as ``wavelength`` giving
            possible weghting functions for the kernels

        :returns lnlikelhood:
            The natural logarithm of the likelihood of the data given the mean
            model spectrum.
        """
        if obs['spectrum'] is None:
            return 0.0

        mask = obs.get('mask', slice(None))
        vectors['mask'] = mask
        vectors['wavelength'] = obs['wavelength']

        delta = (obs['spectrum'] - spec_mu)[mask]

        if spec_noise is not None:
            try:
                spec_noise.compute(**vectors)
                return spec_noise.lnlikelihood(delta)
            except(LinAlgError):
                return np.nan_to_num(-np.inf)
        else:
            # simple noise model
            var = (obs['unc'][mask])**2
            lnp = -0.5*( (delta**2/var).sum() + np.log(2*np.pi*var).sum() )
            return lnp

def lnlike_phot(phot_mu, obs=None, phot_noise=None, **vectors):
    """Calculate the likelihood of the photometric data given the spectroscopic
    model.  Allows for the use of a gaussian process covariance matrix.

    :param phot_mu:
        The mean model sed, in linear flux units (i.e. maggies).

    :param obs: (optional)
        A dictionary of the observational data, including the keys
          *``maggies`` a numpy array of the observed SED, in linear flux
           units
          *``maggies_unc`` the uncertainty of same length as ``maggies``
          *``phot_mask`` optional boolean array of same length as
           ``maggies``
          *``filters`` optional list of sedpy.observate.Filter objects,
           necessary if using fixed filter groups with different gp
           amplitudes for each group.
       If not supplied then the obs dictionary given at initialization will
       be used.

    :param phot_noise: (optional)
        A ``prospect.likelihood.NoiseModel`` object with the methods
        ``compute()`` and ``lnlikelihood()``.  If not supplied a simple chi^2
        likelihood will be evaluated.

    :param vectors:
        A dictionary of possibly relevant vectors of same length as maggies
        that will be passed to the NoiseModel object for constructing weighted
        covariance matrices.

    :returns lnlikelhood:
        The natural logarithm of the likelihood of the data given the mean
        model spectrum.
    """
    if obs['maggies'] is None:
        return 0.0

    mask = obs.get('phot_mask', slice(None))
    delta = (obs['maggies'] - phot_mu)[mask]

    if phot_noise is not None:
        filternames = [f.name for f in obs['filters']]
        vectors['mask'] = mask
        vectors['filternames'] = np.array(filternames)
        try:
            phot_noise.compute(**vectors)
            return phot_noise.lnlikelihood(delta)
        except(LinAlgError):
            return np.nan_to_num(-np.inf)
    else:
        # simple noise model
        var = np.ma.masked_where(not mask, obs['maggies_unc'])**2
        #var = (obs['maggies_unc'][mask])**2
        lnp = -0.5*( (delta**2/var).sum() + np.log(2*np.pi*var).sum() )
        return lnp

def _many_sample(self, xloc, yloc, xa, ya, vra):
        "Solve the Kriging System with more than one sample"
        na = len(vra)
        # number of equations, for simple kriging there're na,
        # for ordinary there're na + 1
        neq = na + self.ktype

        # Establish left hand side covariance matrix:
        left = np.full((neq, neq), np.nan)
        for i, j in product(xrange(na), xrange(na)):
            if np.isnan(left[j, i]):
                left[i, j] = self._cova2(xa[i], ya[i], xa[j], ya[j])
            else:
                left[i, j] = left[j, i]

        # Establish the Right Hand Side Covariance:
        right = list()

        for j in xrange(na):
            xx = xa[j] - xloc
            yy = ya[j] - yloc
            if not self.block_kriging:
                cb = self._cova2(xx, yy, self.xdb[0], self.ydb[0])
            else:
                cb = 0.0
                for i in xrange(self.ndb):
                    cb += self._cova2(xx, yy, self.xdb[i], self.ydb[i])
                    dx = xx - self.xdb[i]
                    dy = yy - self.ydb[i]
                    if dx*dx + dy*dy < np.finfo('float').eps:
                        cb -= self.c0
                cb /= self.ndb
            right.append(cb)

        if self.ktype == 1:  # for ordinary kriging
            # Set the unbiasedness constraint
            left[neq-1, :-1] = self.unbias
            left[:-1, neq-1] = self.unbias
            left[-1, -1] = 0
            right.append(self.unbias)

        # Solve the kriging system
        s = None
        try:
            s = linalg.solve(left, right)
        except linalg.LinAlgError as inst:
            print("Warning kb2d: singular matrix for block " + \
                  "{},{}".format((xloc-self.xmn)/self.xsiz,
                                 (yloc-self.ymn)/self.ysiz))
            return np.nan, np.nan

        estv = self.block_covariance
        if self.ktype == 1:  # ordinary kriging
            estv -= s[-1]*self.unbias  # s[-1] is mu
        est = np.sum(s[:na]*vra[:na])
        estv -= np.sum(s[:na]*right[:na])
        if self.ktype == 0:  # simple kriging
            est += (1 - np.sum(s[:na])) * self.skmean
        return est, estv

def _solve_cholesky_kernel(K, y, alpha, sample_weight=None, copy=False):
    # dual_coef = inv(X X^t + alpha*Id) y
    n_samples = K.shape[0]
    n_targets = y.shape[1]

    if copy:
        K = K.copy()

    alpha = np.atleast_1d(alpha)
    one_alpha = (alpha == alpha[0]).all()
    has_sw = isinstance(sample_weight, np.ndarray) \
        or sample_weight not in [1.0, None]

    if has_sw:
        # Unlike other solvers, we need to support sample_weight directly
        # because K might be a pre-computed kernel.
        sw = np.sqrt(np.atleast_1d(sample_weight))
        y = y * sw[:, np.newaxis]
        K *= np.outer(sw, sw)

    if one_alpha:
        # Only one penalty, we can solve multi-target problems in one time.
        K.flat[::n_samples + 1] += alpha[0]

        try:
            # Note: we must use overwrite_a=False in order to be able to
            #       use the fall-back solution below in case a LinAlgError
            #       is raised
            dual_coef = linalg.solve(K, y, sym_pos=True,
                                     overwrite_a=False)
        except np.linalg.LinAlgError:
            warnings.warn("Singular matrix in solving dual problem. Using "
                          "least-squares solution instead.")
            dual_coef = linalg.lstsq(K, y)[0]

        # K is expensive to compute and store in memory so change it back in
        # case it was user-given.
        K.flat[::n_samples + 1] -= alpha[0]

        if has_sw:
            dual_coef *= sw[:, np.newaxis]

        return dual_coef
    else:
        # One penalty per target. We need to solve each target separately.
        dual_coefs = np.empty([n_targets, n_samples])

        for dual_coef, target, current_alpha in zip(dual_coefs, y.T, alpha):
            K.flat[::n_samples + 1] += current_alpha

            dual_coef[:] = linalg.solve(K, target, sym_pos=True,
                                        overwrite_a=False).ravel()

            K.flat[::n_samples + 1] -= current_alpha

        if has_sw:
            dual_coefs *= sw[np.newaxis, :]

        return dual_coefs.T