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
python类LinAlgError()的实例源码
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
