Python scipy.sparse.linalg.lsqr() Examples
The following are 30
code examples of scipy.sparse.linalg.lsqr().
You can vote up the ones you like or vote down the ones you don't like,
and go to the original project or source file by following the links above each example.
You may also want to check out all available functions/classes of the module
scipy.sparse.linalg
, or try the search function
.
.
Example #1
| Source File: test_dwts.py From pylops with GNU Lesser General Public License v3.0 | 6 votes |
|
def test_DWT_1dsignal(par):
"""Dot-test and inversion for DWT operator for 1d signal
"""
DWTop = DWT(dims=[par['nt']], dir=0, wavelet='haar', level=3)
x = np.random.normal(0., 1., par['nt']) + \
par['imag'] * np.random.normal(0., 1., par['nt'])
assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
complexflag=0 if par['imag'] == 0 else 3)
y = DWTop * x
xadj = DWTop.H * y # adjoint is same as inverse for dwt
xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]
assert_array_almost_equal(x, xadj, decimal=8)
assert_array_almost_equal(x, xinv, decimal=8)
Example #2
| Source File: ridge.py From twitter-stock-recommendation with MIT License | 6 votes |
|
def _solve_lsqr(X, y, alpha, max_iter=None, tol=1e-3):
n_samples, n_features = X.shape
coefs = np.empty((y.shape[1], n_features), dtype=X.dtype)
n_iter = np.empty(y.shape[1], dtype=np.int32)
# According to the lsqr documentation, alpha = damp^2.
sqrt_alpha = np.sqrt(alpha)
for i in range(y.shape[1]):
y_column = y[:, i]
info = sp_linalg.lsqr(X, y_column, damp=sqrt_alpha[i],
atol=tol, btol=tol, iter_lim=max_iter)
coefs[i] = info[0]
n_iter[i] = info[2]
return coefs, n_iter
Example #3
| Source File: optimization_internal.py From chumpy with MIT License | 6 votes |
|
def updateGN(self):
tm = timer()
if sp.issparse(self.A):
self.A.eliminate_zeros()
pif('sparse solve...sparsity infill is %.3f%% (hessian %dx%d)' % (100. * self.A.nnz / (self.A.shape[0] * self.A.shape[1]), self.A.shape[0], self.A.shape[1]))
if self.g.size > 1:
self.d_gn = self.solve(self.A, self.g).ravel()
if np.any(np.isnan(self.d_gn)) or np.any(np.isinf(self.d_gn)):
from scipy.sparse.linalg import lsqr
warnings.warn("sparse solve failed, falling back to lsqr")
self.d_gn = lsqr(self.A, self.g)[0].ravel()
else:
self.d_gn = np.atleast_1d(self.g.ravel()[0]/self.A[0,0])
pif('sparse solve...done in %.2fs' % tm())
else:
pif('dense solve...')
try:
self.d_gn = np.linalg.solve(self.A, self.g).ravel()
except Exception:
warnings.warn("dense solve failed, falling back to lsqr")
self.d_gn = np.linalg.lstsq(self.A, self.g)[0].ravel()
pif('dense solve...done in %.2fs' % tm())
Example #4
| Source File: ridge.py From Splunking-Crime with GNU Affero General Public License v3.0 | 6 votes |
|
def _solve_lsqr(X, y, alpha, max_iter=None, tol=1e-3):
n_samples, n_features = X.shape
coefs = np.empty((y.shape[1], n_features), dtype=X.dtype)
n_iter = np.empty(y.shape[1], dtype=np.int32)
# According to the lsqr documentation, alpha = damp^2.
sqrt_alpha = np.sqrt(alpha)
for i in range(y.shape[1]):
y_column = y[:, i]
info = sp_linalg.lsqr(X, y_column, damp=sqrt_alpha[i],
atol=tol, btol=tol, iter_lim=max_iter)
coefs[i] = info[0]
n_iter[i] = info[2]
return coefs, n_iter
Example #5
| Source File: randomwalk.py From pykernels with MIT License | 6 votes |
|
def _compute(self, data_1, data_2):
data_1 = basic.graphs_to_adjacency_lists(data_1)
data_2 = basic.graphs_to_adjacency_lists(data_2)
res = np.zeros((len(data_1), len(data_2)))
N = len(data_1) * len(data_2)
for i, graph1 in enumerate(data_1):
for j, graph2 in enumerate(data_2):
# norm1, norm2 - normalized adjacency matrixes
norm1 = _norm(graph1)
norm2 = _norm(graph2)
# if graph is unweighted, W_prod = kron(a_norm(g1)*a_norm(g2))
w_prod = kron(lil_matrix(norm1), lil_matrix(norm2))
starting_prob = np.ones(w_prod.shape[0]) / (w_prod.shape[0])
stop_prob = starting_prob
# first solve (I - lambda * W_prod) * x = p_prod
A = identity(w_prod.shape[0]) - (w_prod * self._lmb)
x = lsqr(A, starting_prob)
res[i, j] = stop_prob.T.dot(x[0])
# print float(len(data_2)*i + j)/float(N), "%"
return res
Example #6
| Source File: inference.py From ektelo with Apache License 2.0 | 6 votes |
|
def infer(self, Ms, ys, scale_factors=None):
''' Either:
1) Ms is a single M and ys is a single y
(scale_factors ignored) or
2) Ms and ys are lists of M matrices and y vectors
and scale_factors is a list of the same length.
'''
A, y = _apply_scales(Ms, ys, scale_factors)
if self.method == 'standard':
assert self.l2_reg == 0, 'l2 reg not supported with method=standard'
(x_est, _, rank, _) = linalg.lstsq(A.dense_matrix(), y, lapack_driver='gelsy')
elif self.method == 'lsmr':
res = lsmr(A, y, atol=0, btol=0, damp=self.l2_reg)
x_est = res[0]
elif self.method == 'lsqr':
res = lsqr(A, y, atol=0, btol=0, damp=self.l2_reg)
x_est = res[0]
return x_est
Example #7
| Source File: ridge.py From Mastering-Elasticsearch-7.0 with MIT License | 6 votes |
|
def _solve_lsqr(X, y, alpha, max_iter=None, tol=1e-3):
n_samples, n_features = X.shape
coefs = np.empty((y.shape[1], n_features), dtype=X.dtype)
n_iter = np.empty(y.shape[1], dtype=np.int32)
# According to the lsqr documentation, alpha = damp^2.
sqrt_alpha = np.sqrt(alpha)
for i in range(y.shape[1]):
y_column = y[:, i]
info = sp_linalg.lsqr(X, y_column, damp=sqrt_alpha[i],
atol=tol, btol=tol, iter_lim=max_iter)
coefs[i] = info[0]
n_iter[i] = info[2]
return coefs, n_iter
Example #8
| Source File: OLSUngnarStyle.py From scattertext with Apache License 2.0 | 6 votes |
|
def get_scores_and_p_values(self, tdm, category): ''' Parameters ---------- tdm: TermDocMatrix category: str, category name Returns ------- pd.DataFrame(['coef', 'p-val']) ''' X = tdm._X y = self._make_response_variable_1_or_negative_1(category, tdm) pX = X / X.sum(axis=1) ansX = self._anscombe_transform(pX.copy()) B, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var\ = lsqr(A=ansX, b=y, calc_var=True)
Example #9
| Source File: test_dwts.py From pylops with GNU Lesser General Public License v3.0 | 6 votes |
|
def test_DWT2D_3dsignal(par):
"""Dot-test and inversion for DWT operator for 3d signal
"""
for dirs in [(0, 1), (0, 2), (1, 2)]:
DWTop = DWT2D(dims=(par['nt'], par['nx'], par['ny']),
dirs=dirs, wavelet='haar', level=3)
x = np.random.normal(0., 1., (par['nt'], par['nx'], par['ny'])) + \
par['imag'] * np.random.normal(0., 1., (par['nt'], par['nx'], par['ny']))
assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
complexflag=0 if par['imag'] == 0 else 3)
y = DWTop * x.ravel()
xadj = DWTop.H * y # adjoint is same as inverse for dwt
xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]
assert_array_almost_equal(x.ravel(), xadj, decimal=8)
assert_array_almost_equal(x.ravel(), xinv, decimal=8)
Example #10
| Source File: test_dwts.py From pylops with GNU Lesser General Public License v3.0 | 6 votes |
|
def test_DWT2D_2dsignal(par):
"""Dot-test and inversion for DWT2D operator for 2d signal
"""
DWTop = DWT2D(dims=(par['nt'], par['nx']),
dirs=(0, 1), wavelet='haar', level=3)
x = np.random.normal(0., 1., (par['nt'], par['nx'])) + \
par['imag'] * np.random.normal(0., 1., (par['nt'], par['nx']))
assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
complexflag=0 if par['imag'] == 0 else 3)
y = DWTop * x.ravel()
xadj = DWTop.H * y # adjoint is same as inverse for dwt
xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]
assert_array_almost_equal(x.ravel(), xadj, decimal=8)
assert_array_almost_equal(x.ravel(), xinv, decimal=8)
Example #11
| Source File: test_dwts.py From pylops with GNU Lesser General Public License v3.0 | 6 votes |
|
def test_DWT_3dsignal(par):
"""Dot-test and inversion for DWT operator for 3d signal
"""
for dir in [0, 1, 2]:
DWTop = DWT(dims=(par['nt'], par['nx'], par['ny']),
dir=dir, wavelet='haar', level=3)
x = np.random.normal(0., 1., (par['nt'], par['nx'], par['ny'])) + \
par['imag'] * np.random.normal(0., 1., (par['nt'], par['nx'], par['ny']))
assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
complexflag=0 if par['imag'] == 0 else 3)
y = DWTop * x.ravel()
xadj = DWTop.H * y # adjoint is same as inverse for dwt
xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]
assert_array_almost_equal(x.ravel(), xadj, decimal=8)
assert_array_almost_equal(x.ravel(), xinv, decimal=8)
Example #12
| Source File: test_dwts.py From pylops with GNU Lesser General Public License v3.0 | 6 votes |
|
def test_DWT_2dsignal(par):
"""Dot-test and inversion for DWT operator for 2d signal
"""
for dir in [0, 1]:
DWTop = DWT(dims=(par['nt'], par['nx']),
dir=dir, wavelet='haar', level=3)
x = np.random.normal(0., 1., (par['nt'], par['nx'])) + \
par['imag'] * np.random.normal(0., 1., (par['nt'], par['nx']))
assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
complexflag=0 if par['imag'] == 0 else 3)
y = DWTop * x.ravel()
xadj = DWTop.H * y # adjoint is same as inverse for dwt
xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]
assert_array_almost_equal(x.ravel(), xadj, decimal=8)
assert_array_almost_equal(x.ravel(), xinv, decimal=8)
Example #13
| Source File: LaplacianMesh.py From laplacian-meshes with GNU General Public License v3.0 | 6 votes |
|
def solveLaplacianMesh(mesh, anchors, anchorsIdx, cotangent=True):
n = mesh.VPos.shape[0] # N x 3
k = anchorsIdx.shape[0]
operator = (getLaplacianMatrixUmbrella, getLaplacianMatrixCotangent)
L = operator[1](mesh, anchorsIdx) if cotangent else operator[0](mesh, anchorsIdx)
delta = np.array(L.dot(mesh.VPos))
# augment delta solution matrix with weighted anchors
for i in range(k):
delta[n + i, :] = WEIGHT * anchors[i, :]
# update mesh vertices with least-squares solution
for i in range(3):
mesh.VPos[:, i] = lsqr(L, delta[:, i])[0]
return mesh
#Purpose: Given a few RGB colors on a mesh, smoothly interpolate those colors
#by using their values as anchors and
#Inputs: mesh (polygon mesh object), anchors (a K x 3 numpy array of anchor
#coordinates), colorsIdx (a parallel array of the indices of the RGB anchor indices)
Example #14
| Source File: test_diagonal.py From pylops with GNU Lesser General Public License v3.0 | 6 votes |
|
def test_Diagonal_3dsignal(par):
"""Dot-test and inversion for Diagonal operator for 3d signal
"""
for idim, ddim in enumerate((par['ny'], par['nx'], par['nt'])):
d = np.arange(ddim) + 1. +\
par['imag'] * (np.arange(ddim) + 1.)
Dop = Diagonal(d, dims=(par['ny'], par['nx'], par['nt']),
dir=idim, dtype=par['dtype'])
assert dottest(Dop, par['ny']*par['nx']*par['nt'],
par['ny']*par['nx']*par['nt'],
complexflag=0 if par['imag'] == 0 else 3)
x = np.ones((par['ny'], par['nx'], par['nt'])) + \
par['imag']*np.ones((par['ny'], par['nx'], par['nt']))
xlsqr = lsqr(Dop, Dop * x.ravel(), damp=1e-20, iter_lim=300, show=0)[0]
assert_array_almost_equal(x.ravel(), xlsqr.ravel(), decimal=4)
Example #15
| Source File: test_diagonal.py From pylops with GNU Lesser General Public License v3.0 | 6 votes |
|
def test_Diagonal_2dsignal(par):
"""Dot-test and inversion for Diagonal operator for 2d signal
"""
for idim, ddim in enumerate((par['nx'], par['nt'])):
d = np.arange(ddim) + 1. +\
par['imag'] * (np.arange(ddim) + 1.)
Dop = Diagonal(d, dims=(par['nx'], par['nt']),
dir=idim, dtype=par['dtype'])
assert dottest(Dop, par['nx']*par['nt'], par['nx']*par['nt'],
complexflag=0 if par['imag'] == 0 else 3)
x = np.ones((par['nx'], par['nt'])) + \
par['imag']*np.ones((par['nx'], par['nt']))
xlsqr = lsqr(Dop, Dop * x.ravel(), damp=1e-20, iter_lim=300, show=0)[0]
assert_array_almost_equal(x.ravel(), xlsqr.ravel(), decimal=4)
Example #16
| Source File: test_fredholm.py From pylops with GNU Lesser General Public License v3.0 | 6 votes |
|
def test_Fredholm1(par):
"""Dot-test and inversion for Fredholm1 operator
"""
np.random.seed(10)
_F = np.arange(par['nsl'] * par['nx'] * par['ny']).reshape(par['nsl'],
par['nx'],
par['ny'])
F = _F - par['imag'] * _F
x = np.ones((par['nsl'], par['ny'], par['nz'])) + \
par['imag'] * np.ones((par['nsl'], par['ny'], par['nz']))
Fop = Fredholm1(F, nz=par['nz'], saveGt=par['saveGt'],
usematmul=par['usematmul'], dtype=par['dtype'])
assert dottest(Fop, par['nsl']*par['nx']*par['nz'],
par['nsl']*par['ny']*par['nz'],
complexflag=0 if par['imag'] == 0 else 3)
xlsqr = lsqr(Fop, Fop * x.flatten(), damp=1e-20,
iter_lim=30, show=0)[0]
xlsqr = xlsqr.reshape(par['nsl'], par['ny'], par['nz'])
assert_array_almost_equal(x, xlsqr, decimal=3)
Example #17
| Source File: LaplacianMesh.py From laplacian-meshes with GNU General Public License v3.0 | 6 votes |
|
def smoothColors(mesh, anchors, colorsIdx):
colorsIdx = np.array(colorsIdx)
n = mesh.VPos.shape[0]
k = anchors.shape[0]
colors = np.zeros((n, 3))
delta = np.zeros((n + k, 3))
L = getLaplacianMatrixUmbrella(mesh, colorsIdx);
# augment delta solution matrix with weighted anchors
for i in range(k):
delta[n + i, :] = WEIGHT * anchors[i, :]
# update RGB values with least-squares solution
for i in range(3):
colors[:, i] = lsqr(L, delta[:, i])[0]
return colors
#Purpose: Given a mesh, to smooth it by subtracting off the delta coordinates
#from each vertex, normalized by the degree of that vertex
#Inputs: mesh (polygon mesh object)
#Returns: Nothing (should update mesh.VPos)
Example #18
| Source File: test_basicoperators.py From pylops with GNU Lesser General Public License v3.0 | 6 votes |
|
def test_Regression(par):
"""Dot-test, inversion and apply for Regression operator
"""
np.random.seed(10)
order = 4
t = np.arange(par['ny'], dtype=np.float32)
LRop = Regression(t, order=order, dtype=par['dtype'])
assert dottest(LRop, par['ny'], order+1)
x = np.array([1., 2., 0., 3., -1.], dtype=np.float32)
xlsqr = lsqr(LRop, LRop*x, damp=1e-10, iter_lim=300, show=0)[0]
assert_array_almost_equal(x, xlsqr, decimal=3)
y = LRop * x
y1 = LRop.apply(t, x)
assert_array_almost_equal(y, y1, decimal=3)
Example #19
| Source File: test_kronecker.py From pylops with GNU Lesser General Public License v3.0 | 6 votes |
|
def test_Kroneker(par):
"""Dot-test, inversion and comparison with np.kron for Kronecker operator
"""
np.random.seed(10)
G1 = np.random.normal(0, 10, (par['ny'], par['nx'])).astype(par['dtype'])
G2 = np.random.normal(0, 10, (par['ny'], par['nx'])).astype(par['dtype'])
x = np.ones(par['nx']**2) + par['imag']*np.ones(par['nx']**2)
Kop = Kronecker(MatrixMult(G1, dtype=par['dtype']),
MatrixMult(G2, dtype=par['dtype']),
dtype=par['dtype'])
assert dottest(Kop, par['ny']**2, par['nx']**2,
complexflag=0 if par['imag'] == 0 else 3)
xlsqr = lsqr(Kop, Kop * x, damp=1e-20, iter_lim=300, show=0)[0]
assert_array_almost_equal(x, xlsqr, decimal=2)
# Comparison with numpy
assert_array_almost_equal(np.kron(G1, G2), Kop * np.eye(par['nx']**2),
decimal=3)
Example #20
| Source File: test_basicoperators.py From pylops with GNU Lesser General Public License v3.0 | 6 votes |
|
def test_MatrixMult_repeated(par):
"""Dot-test and inversion for test_MatrixMult operator repeated
along another dimension
"""
np.random.seed(10)
G = np.random.normal(0, 10, (par['ny'], par['nx'])).astype('float32') + \
par['imag'] * np.random.normal(0, 10, (par['ny'],
par['nx'])).astype('float32')
Gop = MatrixMult(G, dims=5, dtype=par['dtype'])
assert dottest(Gop, par['ny']*5, par['nx']*5,
complexflag=0 if par['imag'] == 1 else 3)
x = (np.ones((par['nx'], 5)) +
par['imag'] * np.ones((par['nx'], 5))).flatten()
xlsqr = lsqr(Gop, Gop*x, damp=1e-20, iter_lim=300, show=0)[0]
assert_array_almost_equal(x, xlsqr, decimal=4)
Example #21
| Source File: test_lsqr.py From Computable with MIT License | 5 votes |
|
def test_gh_2466():
row = np.array([0, 0])
col = np.array([0, 1])
val = np.array([1, -1])
A = scipy.sparse.coo_matrix((val, (row, col)), shape=(1, 2))
b = np.asarray([4])
lsqr(A, b)
Example #22
| Source File: test_lsqr.py From Computable with MIT License | 5 votes |
|
def test_basic():
svx = np.linalg.solve(G, b)
X = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=maxit)
xo = X[0]
assert_(norm(svx - xo) < 1e-5)
Example #23
| Source File: test_smoothing.py From pylops with GNU Lesser General Public License v3.0 | 5 votes |
|
def test_Smoothing2D(par):
"""Dot-test for smoothing
"""
# 2d kernel on 2d signal
if par['dir'] < 2:
D2op = Smoothing2D(nsmooth=(5, 5), dims=(par['ny'], par['nx']), dtype='float32')
assert dottest(D2op, par['ny']*par['nx'], par['ny']*par['nx'])
x = np.zeros((par['ny'], par['nx']))
x[par['ny']//2, par['nx']//2] = 1.
x = x.flatten()
y = D2op*x
xlsqr = lsqr(D2op, y, damp=1e-10, iter_lim=400, show=0)[0]
assert_array_almost_equal(x, xlsqr, decimal=1)
# 2d kernel on 3d signal
D2op = Smoothing2D(nsmooth=(5, 5), dims=(par['nz'], par['ny'], par['nx']),
nodir=par['dir'], dtype='float32')
assert dottest(D2op, par['nz'] * par['ny'] * par['nx'], par['nz'] * par['ny'] * par['nx'])
x = np.zeros((par['nz'], par['ny'], par['nx']))
x[par['nz']//2, par['ny']//2, par['nx']//2] = 1.
x = x.flatten()
y = D2op*x
xlsqr = lsqr(D2op, y, damp=1e-10, iter_lim=400, show=0)[0]
assert_array_almost_equal(x, xlsqr, decimal=1)
Example #24
| Source File: test_smoothing.py From pylops with GNU Lesser General Public License v3.0 | 5 votes |
|
def test_Smoothing1D(par):
"""Dot-test and inversion for smoothing
"""
# 1d kernel on 1d signal
D1op = Smoothing1D(nsmooth=5, dims=par['nx'], dtype='float32')
assert dottest(D1op, par['nx'], par['nx'])
x = np.random.normal(0, 1, par['nx'])
y = D1op * x
xlsqr = lsqr(D1op, y, damp=1e-10, iter_lim=100, show=0)[0]
assert_array_almost_equal(x, xlsqr, decimal=3)
# 1d kernel on 2d signal
D1op = Smoothing1D(nsmooth=5, dims=(par['ny'], par['nx']), dir=par['dir'], dtype='float32')
assert dottest(D1op, par['ny']*par['nx'], par['ny']*par['nx'])
x = np.random.normal(0, 1, (par['ny'], par['nx'])).flatten()
y = D1op*x
xlsqr = lsqr(D1op, y, damp=1e-10, iter_lim=100, show=0)[0]
assert_array_almost_equal(x, xlsqr, decimal=3)
# 1d kernel on 3d signal
D1op = Smoothing1D(nsmooth=5, dims=(par['nz'], par['ny'], par['nx']),
dir=par['dir'], dtype='float32')
assert dottest(D1op, par['nz'] * par['ny'] * par['nx'], par['nz'] * par['ny'] * par['nx'])
x = np.random.normal(0, 1, (par['nz'], par['ny'], par['nx'])).flatten()
y = D1op*x
xlsqr = lsqr(D1op, y, damp=1e-10, iter_lim=100, show=0)[0]
assert_array_almost_equal(x, xlsqr, decimal=3)
Example #25
| Source File: test_lsqr.py From GraphicDesignPatternByPython with MIT License | 5 votes |
|
def test_basic():
b_copy = b.copy()
X = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=maxit)
assert_(np.all(b_copy == b))
svx = np.linalg.solve(G, b)
xo = X[0]
assert_(norm(svx - xo) < 1e-5)
Example #26
| Source File: test_lsqr.py From GraphicDesignPatternByPython with MIT License | 5 votes |
|
def test_gh_2466():
row = np.array([0, 0])
col = np.array([0, 1])
val = np.array([1, -1])
A = scipy.sparse.coo_matrix((val, (row, col)), shape=(1, 2))
b = np.asarray([4])
lsqr(A, b)
Example #27
| Source File: test_lsqr.py From GraphicDesignPatternByPython with MIT License | 5 votes |
|
def test_well_conditioned_problems():
# Test that sparse the lsqr solver returns the right solution
# on various problems with different random seeds.
# This is a non-regression test for a potential ZeroDivisionError
# raised when computing the `test2` & `test3` convergence conditions.
n = 10
A_sparse = scipy.sparse.eye(n, n)
A_dense = A_sparse.toarray()
with np.errstate(invalid='raise'):
for seed in range(30):
rng = np.random.RandomState(seed + 10)
beta = rng.rand(n)
beta[beta == 0] = 0.00001 # ensure that all the betas are not null
b = A_sparse * beta[:, np.newaxis]
output = lsqr(A_sparse, b, show=show)
# Check that the termination condition corresponds to an approximate
# solution to Ax = b
assert_equal(output[1], 1)
solution = output[0]
# Check that we recover the ground truth solution
assert_array_almost_equal(solution, beta)
# Sanity check: compare to the dense array solver
reference_solution = np.linalg.solve(A_dense, b).ravel()
assert_array_almost_equal(solution, reference_solution)
Example #28
| Source File: test_lsqr.py From GraphicDesignPatternByPython with MIT License | 5 votes |
|
def test_b_shapes():
# Test b being a scalar.
A = np.array([[1.0, 2.0]])
b = 3.0
x = lsqr(A, b)[0]
assert_almost_equal(norm(A.dot(x) - b), 0)
# Test b being a column vector.
A = np.eye(10)
b = np.ones((10, 1))
x = lsqr(A, b)[0]
assert_almost_equal(norm(A.dot(x) - b.ravel()), 0)
Example #29
| Source File: test_lsqr.py From GraphicDesignPatternByPython with MIT License | 5 votes |
|
def test_initialization():
# Test the default setting is the same as zeros
b_copy = b.copy()
x_ref = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=maxit)
x0 = np.zeros(x_ref[0].shape)
x = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=maxit, x0=x0)
assert_(np.all(b_copy == b))
assert_array_almost_equal(x_ref[0], x[0])
# Test warm-start with single iteration
x0 = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=1)[0]
x = lsqr(G, b, show=show, atol=tol, btol=tol, iter_lim=maxit, x0=x0)
assert_array_almost_equal(x_ref[0], x[0])
assert_(np.all(b_copy == b))
Example #30
| Source File: mesh_edit.py From hmd with MIT License | 5 votes |
|
def solveLaplacianMesh(mesh, anchors, anchorsIdx, cotangent=True):
n = mesh.n_vertices()
k = anchorsIdx.shape[0]
operator = (getLaplacianMatrixUmbrella, getLaplacianMatrixCotangent)
L = operator[1](mesh, anchorsIdx) if cotangent else operator[0](mesh, anchorsIdx)
delta = np.array(L.dot(mesh.points()))
# augment delta solution matrix with weighted anchors
for i in range(k):
delta[n + i, :] = WEIGHT * anchors[i, :]
# update mesh vertices with least-squares solution
for i in range(3):
#mesh.points()[:, i] = lsqr(L, delta[:, i])[0]
mesh.points()[:, i] = sparseqr.solve(L, delta[:, i], tolerance = 1e-8)
return mesh
##############################################################
## High Speed Laplacian Mesh Editing ##
##############################################################
# using umbrella weights for higher speed
