Here are the examples of the python api numpy.linalg.solve taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.
93 Examples
4
Example 1
def solve(self, wt_n, y_nd, bend_coef, f_res):
if y_nd.shape[0] != self.n or y_nd.shape[1] != self.d:
raise RuntimeError("The dimensions of y_nd doesn't match the dimensions of x_nd")
WQN = wt_n[:, None] * self.QN
lhs = self.QN.T.dot(WQN) + bend_coef * self.NKN + self.NRN
rhs = self.NR + WQN.T.dot(y_nd)
z = np.linalg.solve(lhs, rhs)
theta = self.N.dot(z)
f_res.update(self.x_nd, y_nd, bend_coef, self.rot_coef, wt_n, theta, N=self.N, z=z)
3
Example 2
Project: pydy Source File: ode_function_generators.py
@linear_sys_solver.setter
def linear_sys_solver(self, v):
if isinstance(v, type(lambda x: x)):
self._solve_linear_system = v
elif v == 'numpy':
self._solve_linear_system = numpy.linalg.solve
elif v == 'scipy':
self._solve_linear_system = scipy.linalg.solve
else:
msg = '{} is not a valid solver.'
raise ValueError(msg.format(self.linear_sys_solver))
3
Example 3
Project: hedge Source File: local.py
@memoize_method
def stiffness_t_matrices(self):
return [numpy.asarray(
la.solve(
self.ldis.vandermonde().T,
numpy.dot(
diff_vdm.T,
numpy.diag(self.volume_weights))),
order="C")
for diff_vdm in self.diff_vandermonde_matrices()]
3
Example 4
Project: AWS-Lambda-ML-Microservice-Skeleton Source File: test_linalg.py
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))
3
Example 5
Project: FreeCAD_drawing_dimensioning Source File: svgLib_dd.py
def t_of_position( self, pos ):
pos_element = numpy.linalg.solve( self.T_element, pos - self.c_element )
A = numpy.linalg.solve(self.T_ellipse, pos_element)
B = self.center_circular
C = B + array([cos(self.theta_start), sin(self.theta_start)])
AB = (A - B) / norm(A-B)
BC = C - B
theta = arccos2(dot(AB,BC))
D = array([ A[0] - B[0], A[1] - B[1], 0.0 ])
E = array([ A[0] - C[0], A[1] - C[1], 0.0 ])
F = cross(D,E)
if F[2] < 0:
theta = -theta
#print(theta, self.dtheta, theta / self.dtheta )
return theta / self.dtheta
3
Example 6
def beta(self):
'''\
The linear estimation of the parameter vector :math:`\beta` given by
.. math::
\beta = (X^T X)^-1 X^T y
'''
t = self.X.transpose()
XX = dot(t,self.X)
XY = dot(t,self.y)
return linalg.solve(XX,XY)
3
Example 7
Project: pymatgen Source File: analyzer.py
def get_facet_chempots(self, facet):
"""
Calculates the chemical potentials for each element within a facet.
Args:
facet: Facet of the phase diagram.
Returns:
{ element: chempot } for all elements in the phase diagram.
"""
complist = [self._pd.qhull_entries[i].composition for i in facet]
energylist = [self._pd.qhull_entries[i].energy_per_atom for i in facet]
m = self._make_comp_matrix(complist)
chempots = np.linalg.solve(m, energylist)
return dict(zip(self._pd.elements, chempots))
3
Example 8
Project: hedge Source File: local.py
@memoize_method
def face_mass_matrix(self):
return numpy.asarray(
la.solve(
self.ldis.face_vandermonde().T,
numpy.dot(
self.face_vandermonde().T,
numpy.diag(self.face_weights))),
order="C")
3
Example 9
def do_problem(elements):
Kg = build_global_stiffness() # global stiffness (Kg)
Kr = apply_boundary_conditions() # reduced stiffness (Kr)
#F = K*x
F = get_forces()
x = solve(Kr, F) # deflections
(stress, strain) = recover_stress_Strain(elements, x)
3
Example 10
def test_solve(self, mat, b, x, f):
"""Solve a linear system where the solution is equal to the right-hand
side and check the result."""
mat.assemble()
x = np.linalg.solve(mat.values, b.data)
eps = 1.e-8
assert_allclose(x, f.data, eps)
3
Example 11
def E_dist_mahalanobis(self, X ):
'''Calculate Mahalanobis distance to x
dist(x) = dX'*E[Lambda]*dX
If X is a matrix, we compute a vector of distances to each row vector
dist(X)[n] = dX[n]'*E[Lambda]*dX[n]
'''
xWx = np.linalg.solve(self.cholinvW(), X.T)
xWx *= xWx
return self.v * np.sum(xWx, axis=0)
3
Example 12
Project: GPflow Source File: test_variational.py
def referenceMultivariatePriorKL(meanA, covA, meanB, covB):
# KL[ qA | qB ] = E_{qA} \log [qA / qB] where qA and aB are K dimensional multivariate normal distributions.
# Analytically tractable and equal to...
# 0.5 * (Tr(covB^{-1} covA) + (meanB - meanA)^T covB^{-1} (meanB - meanA) - K + log(det(covB)) - log (det(covA)))
K = covA.shape[0]
traceTerm = 0.5 * np.trace(np.linalg.solve(covB, covA))
delta = meanB - meanA
mahalanobisTerm = 0.5 * np.dot(delta.T, np.linalg.solve(covB, delta))
constantTerm = -0.5 * K
priorLogDeterminantTerm = 0.5*np.linalg.slogdet(covB)[1]
variationalLogDeterminantTerm = -0.5 * np.linalg.slogdet(covA)[1]
return traceTerm + mahalanobisTerm + constantTerm + priorLogDeterminantTerm + variationalLogDeterminantTerm
3
Example 13
def test_exp(self):
man = self.man
x = man.rand()
u = man.randvec(x)
e = sp.linalg.expm(la.solve(x, u))
np_testing.assert_allclose(multiprod(x, e), man.exp(x, u))
u = u * 1e-6
np_testing.assert_allclose(man.exp(x, u), x + u)
3
Example 14
Project: ignition Source File: test_iterative_methods.py
def cg_alg_driver(alg, n=100):
A, b = get_fd_poisson(n)
iters, x = alg(A, b, n)
x_numpy = solve(A, b)
diff = x - x_numpy
l_inf = max(diff)
l_2 = sqrt(dot(diff, diff))
# assert(l_inf < 1e-10)
# assert(l_2 < 1e-10)
return iters, l_2, l_inf
3
Example 15
Project: pymanopt Source File: grassmann.py
def log(self, X, Y):
ytx = multiprod(multitransp(Y), X)
At = multitransp(Y) - multiprod(ytx, multitransp(X))
Bt = np.linalg.solve(ytx, At)
u, s, vt = svd(multitransp(Bt), full_matrices=False)
arctan_s = np.expand_dims(np.arctan(s), -2)
U = multiprod(u * arctan_s, vt)
return U
3
Example 16
Project: hedge Source File: gmsh.py
def __init__(self, nodes, ldis):
self.nodes = nodes
self.ldis = ldis
node_src_indices = np.array(
ldis.get_lexicographic_gmsh_node_indices(),
dtype=np.intp)
nodes = np.array(nodes, dtype=np.float64)
reordered_nodes = nodes[node_src_indices, :]
self.modal_coeff = la.solve(
ldis.equidistant_vandermonde(), reordered_nodes)
# axis 0: node number, axis 1: xyz axis
if False:
for i, c in zip(ldis.generate_mode_identifiers(), self.modal_coeff):
print i, c
3
Example 17
Project: pymatgen Source File: analyzer.py
def get_decomposition(self, comp):
"""
Provides the decomposition at a particular composition.
Args:
comp: A composition
Returns:
Decomposition as a dict of {Entry: amount}
"""
facet = self._get_facet(comp)
comp_list = [self._pd.qhull_entries[i].composition for i in facet]
m = self._make_comp_matrix(comp_list)
compm = self._make_comp_matrix([comp])
decomp_amts = np.linalg.solve(m.T, compm.T)
return {self._pd.qhull_entries[f]: amt[0]
for f, amt in zip(facet, decomp_amts)
if abs(amt[0]) > PDAnalyzer.numerical_tol}
3
Example 18
Project: pyhsmm Source File: general.py
def solve_psd(A,b,chol=None,overwrite_b=False,overwrite_A=False):
if A.shape[0] < 5000 and chol is None:
return np.linalg.solve(A,b)
else:
if chol is None:
chol = np.linalg.cholesky(A)
return scipy.linalg.solve_triangular(
chol.T,
scipy.linalg.solve_triangular(chol,b,lower=True,overwrite_b=overwrite_b),
lower=False,overwrite_b=True)
3
Example 19
@memoize_method
def mass_matrix(self):
return numpy.asarray(
la.solve(
self.ldis.vandermonde().T,
numpy.dot(
self.vandermonde().T,
numpy.diag(self.volume_weights))),
order="C")
3
Example 20
Project: elephant Source File: icsd.py
def get_csd(self, ):
'''
Perform the CSD estimate from the LFP and forward matrix F, i.e as
CSD=F**-1*LFP
Arguments
---------
Returns
-------
csd : np.ndarray * quantity.Quantity
Array with the csd estimate
'''
csd = np.linalg.solve(self.f_matrix, self.lfp)
return csd * (self.f_matrix.units**-1*self.lfp.units).simplified
3
Example 21
def solve(a, b):
"""Returns the solution of A X = B."""
try:
return linalg.solve(a, b)
except linalg.LinAlgError:
return np.dot(linalg.pinv(a), b)
3
Example 22
def dist_mahalanobis(self, X):
'''Calculate Mahalanobis distance to x
dist(x) = (x-m)'*W*(x-m)
If X is a matrix, we compute a vector of distances to each row vector
dist(X)[n] = (x_n-m)'*W*(x_n-m)
'''
dX = (X-self.m).T
Q = np.linalg.solve( self.cholinvW(), dX )
Q *= Q
return np.sum( Q, axis=0)
3
Example 23
def estimate(xlist, tlist):
"""訓練データからパラメータを推定"""
# M次多項式のときはM+1個のパラメータがある
A = []
for i in range(M + 1):
for j in range(M + 1):
temp = (xlist ** (i + j)).sum()
A.append(temp)
A = array(A).reshape(M + 1, M + 1)
T = []
for i in range(M + 1):
T.append(((xlist ** i) * tlist).sum())
T = array(T)
# パラメータwは線形連立方程式の解
wlist = np.linalg.solve(A, T)
return wlist
3
Example 24
def dist_mahalanobis(self, X):
''' Given NxD matrix X, compute Nx1 vector Dist
Dist[n] = ( X[n]-m )' L (X[n]-m)
'''
if self.doSigma:
Q = np.linalg.solve(self.cholSigma(), X.T)
else:
Q = dotABT(self.cholL(), X)
Q *= Q
return np.sum(Q, axis=0)
3
Example 25
def horner(self, s):
'''Evaluate the systems's transfer function for a complex variable
Returns a matrix of values evaluated at complex variable s.
'''
resp = self.C * solve(s * eye(self.states) - self.A,
self.B) + self.D
return array(resp)
3
Example 26
@memoized
def risky_ss(self):
oo = self.order
if oo == 1:
return self['ys']
elif oo in (2,3):
#TODO: RSS for order 3 should be computed differently
import numpy.linalg
A = self['g_a']
I = np.eye(A.shape[0])
D = self['g_ss']/2
dx = numpy.linalg.solve( I - A, D)
return self['ys'] + dx
3
Example 27
def estimate(xlist, tlist):
A = []
for i in range(M + 1):
for j in range(M + 1):
temp = (xlist ** (i + j)).sum()
A.append(temp)
A = array(A).reshape(M + 1, M + 1)
T = []
for i in range(M + 1):
T.append(((xlist ** i) * tlist).sum())
T = array(T)
# パラメータwは線形連立方程式の解
wlist = np.linalg.solve(A, T)
return wlist
3
Example 28
Project: pybasicbayes Source File: regression.py
def resample(self,data=[],stats=None):
if len(data) > 0 or stats is not None:
stats = self._get_statistics(data) if stats is None else stats
for itr in range(self.niter):
self.A, self.sigma = \
sample_mniw(*self._natural_to_standard(
self.natural_hypparam + stats))
mat = self.M_0 - self.A
self.resample_K(1./2*np.einsum(
'ij,ij->j',mat,np.linalg.solve(self.sigma,mat)))
else:
self.resample_K()
super(ARDRegression,self).resample()
3
Example 29
Project: FreeCAD_assembly2 Source File: degreesOfFreedom.py
def determine_R_about_axis(self, R_effective, checkAnswer=True, tol=10**-12): #not used anymore
'determine R_about_axis so that R_effective = R_about_axis * R_to_align_axis'
A = self.R_to_align_axis.transpose()
X = numpy.array([
numpy.linalg.solve(A, R_effective[row,:]) for row in range(3)
])
#prettyPrintArray(X)
if checkAnswer:
print(' determine_R_about_axis: diff between R_effective and R_about_axis * R_to_align_axis (should be all close to zero):')
error = R_effective - dotProduct(X, self.R_to_align_axis)
assert norm(error) <= tol
return X
3
Example 30
def test_exp(self):
# Test against manopt implementation, test that for small vectors
# exp(x, u) = x + u.
man = self.man
x = man.rand()
u = man.randvec(x)
e = np.zeros((self.k, self.n, self.n))
for i in range(self.k):
e[i] = sp.linalg.expm(la.solve(x[i], u[i]))
np_testing.assert_allclose(multiprod(x, e), man.exp(x, u))
u = u * 1e-6
np_testing.assert_allclose(man.exp(x, u), x + u)
3
Example 31
Project: pgmult Source File: distributions.py
def conditional_psi(self, x):
"""
Compute the conditional distribution over psi given observation x and omega
:param x:
:return:
"""
assert x.ndim == 2
Omega = np.diag(self.omega)
Sigma_cond = inv(Omega + inv(self.Sigma))
# kappa is the mean dot precision, i.e. the sufficient statistic of a Gaussian
# therefore we can sum over datapoints
kappa = kappa_vec(x).sum(0)
mu_cond = Sigma_cond.dot(kappa +
solve(self.Sigma, self.mu))
return mu_cond, Sigma_cond
3
Example 32
def test_dist(self):
k = self.k
# n = self.n
man = self.man
x = man.rand()
y = man.rand()
# Compare to the implementation in manopt
d = la.solve(x, y)
for i in range(k):
d[i] = sp.linalg.logm(d[i])
np_testing.assert_almost_equal(man.dist(x, y), la.norm(d))
3
Example 33
def _precompute(self):
# call super method
super(InverseCompositional, self)._precompute()
# compute appearance model mean gradient
nabla_t = self.interface.gradient(self.template)
# compute masked inverse Jacobian
self.J_m = self.interface.steepest_descent_images(-nabla_t, self.dW_dp)
# compute masked inverse Hessian
self.JJ_m = self.J_m.T.dot(self.J_m)
# compute masked Jacobian pseudo-inverse
self.pinv_J_m = np.linalg.solve(self.JJ_m, self.J_m.T)
3
Example 34
Project: pyBAST Source File: classes.py
def llik(self,P=np.array( [0., 0., 0., 0., 0., 1., 1.] ) ):
""" Compute log-likelihood of parameter set P within
likelihood distribution bgmap object.
"""
delta = self.mu - P
sigma = self.sigma.copy()
# Interpret infs in covariance matrix as contributing
# zero to the chi^2 (set delta[i] = 0).
I = np.nonzero(np.isinf(np.diag(sigma)))[0]
delta[I] = 0
sigma[I,:] = 0
sigma[:,I] = 0
sigma[I,I] = 1
return -0.5 * delta.dot( solve( sigma, delta ) )
3
Example 35
def exp(self, x, u):
# TODO: Check which method is faster depending on n, k.
x_inv_u = np.linalg.solve(x, u)
if self._k > 1:
e = np.zeros(np.shape(x))
for i in range(self._k):
e[i] = sp.linalg.expm(x_inv_u[i])
else:
e = sp.linalg.expm(x_inv_u)
return multiprod(x, e)
3
Example 36
@memoized
def risky_ss(self):
oo = self.order
if oo == 1:
return self['ys']
elif oo in (2,3):
#TODO: RSS for order 3 should be computed differently
import numpy.linalg
A = self['g_a']
I = np.eye(A.shape[0])
D = self['g_ss']/2
dx = numpy.linalg.solve( I - A, D)
return self['ys'] + dx
2
Example 37
Project: tfdeploy Source File: tfdeploy.py
@Operation.factory(types=("MatrixSolve", "BatchMatrixSolve"))
def MatrixSolve(a, b):
"""
Matrix solve op.
"""
return np.linalg.solve(a, b),
0
Example 38
def dcgain(self):
"""Return the zero-frequency gain
The zero-frequency gain of a state-space system is given by:
.. math: G(0) = - C A^{-1} B + D
Returns
-------
gain : ndarray
The zero-frequency gain, or np.nan if the system has a pole
at the origin
"""
try:
gain = np.asarray(self.D -
self.C.dot(np.linalg.solve(self.A, self.B)))
except LinAlgError:
# zero eigenvalue: singular matrix
return np.nan
return np.squeeze(gain)
0
Example 39
Project: python-control Source File: statesp.py
def _rss_generate(states, inputs, outputs, type):
"""Generate a random state space.
This does the actual random state space generation expected from rss and
drss. type is 'c' for continuous systems and 'd' for discrete systems.
"""
# Probability of repeating a previous root.
pRepeat = 0.05
# Probability of choosing a real root. Note that when choosing a complex
# root, the conjugate gets chosen as well. So the expected proportion of
# real roots is pReal / (pReal + 2 * (1 - pReal)).
pReal = 0.6
# Probability that an element in B or C will not be masked out.
pBCmask = 0.8
# Probability that an element in D will not be masked out.
pDmask = 0.3
# Probability that D = 0.
pDzero = 0.5
# Check for valid input arguments.
if states < 1 or states % 1:
raise ValueError("states must be a positive integer. states = %g." %
states)
if inputs < 1 or inputs % 1:
raise ValueError("inputs must be a positive integer. inputs = %g." %
inputs)
if outputs < 1 or outputs % 1:
raise ValueError("outputs must be a positive integer. outputs = %g." %
outputs)
# Make some poles for A. Preallocate a complex array.
poles = zeros(states) + zeros(states) * 0.j
i = 0
while i < states:
if rand() < pRepeat and i != 0 and i != states - 1:
# Small chance of copying poles, if we're not at the first or last
# element.
if poles[i-1].imag == 0:
# Copy previous real pole.
poles[i] = poles[i-1]
i += 1
else:
# Copy previous complex conjugate pair of poles.
poles[i:i+2] = poles[i-2:i]
i += 2
elif rand() < pReal or i == states - 1:
# No-oscillation pole.
if type == 'c':
poles[i] = -exp(randn()) + 0.j
elif type == 'd':
poles[i] = 2. * rand() - 1.
i += 1
else:
# Complex conjugate pair of oscillating poles.
if type == 'c':
poles[i] = complex(-exp(randn()), 3. * exp(randn()))
elif type == 'd':
mag = rand()
phase = 2. * pi * rand()
poles[i] = complex(mag * cos(phase),
mag * sin(phase))
poles[i+1] = complex(poles[i].real, -poles[i].imag)
i += 2
# Now put the poles in A as real blocks on the diagonal.
A = zeros((states, states))
i = 0
while i < states:
if poles[i].imag == 0:
A[i, i] = poles[i].real
i += 1
else:
A[i, i] = A[i+1, i+1] = poles[i].real
A[i, i+1] = poles[i].imag
A[i+1, i] = -poles[i].imag
i += 2
# Finally, apply a transformation so that A is not block-diagonal.
while True:
T = randn(states, states)
try:
A = dot(solve(T, A), T) # A = T \ A * T
break
except LinAlgError:
# In the unlikely event that T is rank-deficient, iterate again.
pass
# Make the remaining matrices.
B = randn(states, inputs)
C = randn(outputs, states)
D = randn(outputs, inputs)
# Make masks to zero out some of the elements.
while True:
Bmask = rand(states, inputs) < pBCmask
if any(Bmask): # Retry if we get all zeros.
break
while True:
Cmask = rand(outputs, states) < pBCmask
if any(Cmask): # Retry if we get all zeros.
break
if rand() < pDzero:
Dmask = zeros((outputs, inputs))
else:
Dmask = rand(outputs, inputs) < pDmask
# Apply masks.
B = B * Bmask
C = C * Cmask
D = D * Dmask
return StateSpace(A, B, C, D)
0
Example 40
Project: QuantEcon.py Source File: matrix_eqn.py
def solve_discrete_riccati(A, B, Q, R, N=None, tolerance=1e-10, max_iter=500):
"""
Solves the discrete-time algebraic Riccati equation
X = A'XA - (N + B'XA)'(B'XB + R)^{-1}(N + B'XA) + Q
via a modified structured doubling algorithm. An explanation of the
algorithm can be found in the reference below.
Note that SciPy also has a discrete riccati equation solver. However it
cannot handle the case where R is not invertible, or when N is nonzero.
Both of these cases can be handled in the algorithm implemented below.
Parameters
----------
A : array_like(float, ndim=2)
k x k array.
B : array_like(float, ndim=2)
k x n array
Q : array_like(float, ndim=2)
k x k, should be symmetric and non-negative definite
R : array_like(float, ndim=2)
n x n, should be symmetric and positive definite
N : array_like(float, ndim=2)
n x k array
tolerance : scalar(float), optional(default=1e-10)
The tolerance level for convergence
max_iter : scalar(int), optional(default=500)
The maximum number of iterations allowed
Returns
-------
X : array_like(float, ndim=2)
The fixed point of the Riccati equation; a k x k array
representing the approximate solution
References
----------
Chiang, Chun-Yueh, Hung-Yuan Fan, and Wen-Wei Lin. "STRUCTURED DOUBLING
ALGORITHM FOR DISCRETE-TIME ALGEBRAIC RICCATI EQUATIONS WITH SINGULAR
CONTROL WEIGHTING MATRICES." Taiwanese Journal of Mathematics 14, no. 3A
(2010): pp-935.
"""
# == Set up == #
error = tolerance + 1
fail_msg = "Convergence failed after {} iterations."
# == Make sure that all array_likes are np arrays, two-dimensional == #
A, B, Q, R = np.atleast_2d(A, B, Q, R)
n, k = R.shape[0], Q.shape[0]
I = np.identity(k)
if N is None:
N = np.zeros((n, k))
else:
N = np.atleast_2d(N)
# == Choose optimal value of gamma in R_hat = R + gamma B'B == #
current_min = np.inf
candidates = (0.0, 0.01, 0.1, 0.25, 0.5, 1.0, 2.0, 10.0, 100.0, 10e5)
BB = dot(B.T, B)
BTA = dot(B.T, A)
for gamma in candidates:
Z = R + gamma * BB
cn = np.linalg.cond(Z)
if np.isfinite(cn):
Q_tilde = - Q + dot(N.T, solve(Z, N + gamma * BTA)) + gamma * I
G0 = dot(B, solve(Z, B.T))
A0 = dot(I - gamma * G0, A) - dot(B, solve(Z, N))
H0 = gamma * dot(A.T, A0) - Q_tilde
f1 = np.linalg.cond(Z, np.inf)
f2 = gamma * f1
f3 = np.linalg.cond(I + dot(G0, H0))
f_gamma = max(f1, f2, f3)
if f_gamma < current_min:
best_gamma = gamma
current_min = f_gamma
# == If no candidate successful then fail == #
if current_min == np.inf:
msg = "Unable to initialize routine due to ill conditioned arguments"
raise ValueError(msg)
gamma = best_gamma
R_hat = R + gamma * BB
# == Initial conditions == #
Q_tilde = - Q + dot(N.T, solve(R_hat, N + gamma * BTA)) + gamma * I
G0 = dot(B, solve(R_hat, B.T))
A0 = dot(I - gamma * G0, A) - dot(B, solve(R_hat, N))
H0 = gamma * dot(A.T, A0) - Q_tilde
i = 1
# == Main loop == #
while error > tolerance:
if i > max_iter:
raise ValueError(fail_msg.format(i))
else:
A1 = dot(A0, solve(I + dot(G0, H0), A0))
G1 = G0 + dot(dot(A0, G0), solve(I + dot(H0, G0), A0.T))
H1 = H0 + dot(A0.T, solve(I + dot(H0, G0), dot(H0, A0)))
error = np.max(np.abs(H1 - H0))
A0 = A1
G0 = G1
H0 = H1
i += 1
return H1 + gamma * I # Return X
0
Example 41
Project: ray Source File: linalg.py
@ray.remote
def solve(a, b):
return np.linalg.solve(a, b)
0
Example 42
Project: pymanopt Source File: psd.py
def inner(self, x, u, v):
return np.tensordot(la.solve(x, u), la.solve(x, v), axes=x.ndim)
0
Example 43
def energy(self, x):
L1x = solve(self.L, x)
return .5 * dot(L1x.T, L1x)
0
Example 44
def solve(self):
assert(self.system_force_vector is not None), "There are no forces on the structure"
self.__assemble_system_matrix()
self.__process_conditions()
# solution of the reduced system (reduced due to support conditions)
reduced_displacement_vector = np.linalg.solve(self.reduced_system_matrix, self.reduced_force_vector)
# add the solution of the reduced system in the complete system displacement vector
self.system_displacement_vector = np.zeros(self.shape_system_matrix)
count = 0
for i in self.remainder_indexes:
self.system_displacement_vector[i] = reduced_displacement_vector[count]
count += 1
# determine the displacement vector of the elements
for el in self.elements:
index_node_1 = (el.node_1.ID - 1) * 3
index_node_2 = (el.node_2.ID - 1) * 3
for i in range(3): # node 1 ux, uz, phi
el.element_displacement_vector[i] = self.system_displacement_vector[index_node_1 + i]
for i in range(3): # node 2 ux, uz, phi
el.element_displacement_vector[3 + i] = self.system_displacement_vector[index_node_2 + i]
el.determine_force_vector()
# determining the node results in post processing class
self.post_processor.node_results()
self.post_processor.reaction_forces()
self.post_processor.element_results()
# check the values in the displacement vector for extreme values, indicating a flawed calculation
for value in np.nditer(self.system_displacement_vector):
assert(value < 1e6), "The displacements of the structure exceed 1e6. Check your support conditions," \
"or your elements Young's modulus"
return self.system_displacement_vector
0
Example 45
def random(self):
n = normal(size=self.L.shape[0])
return solve(self.L.T, n)
0
Example 46
def norm(self, x, u):
return la.norm(la.solve(x, u))
0
Example 47
def feedback(self, other=1, sign=-1):
"""Feedback interconnection between two FRD objects."""
other = _convertToFRD(other, omega=self.omega)
if (self.outputs != other.inputs or
self.inputs != other.outputs):
raise ValueError(
"FRD.feedback, inputs/outputs mismatch")
fresp = empty((self.outputs, self.inputs, len(other.omega)),
dtype=complex)
# TODO: vectorize this
# TODO: handle omega re-mapping
for k, w in enumerate(other.omega):
fresp[:, :, k] = self.fresp[:, :, k].view(type=matrix)* \
linalg.solve(
eye(self.inputs) +
other.fresp[:, :, k].view(type=matrix) *
self.fresp[:, :, k].view(type=matrix),
eye(self.inputs))
return FRD(fresp, other.omega, smooth=(self.ifunc is not None))
0
Example 48
Project: python-whiteboard Source File: wiimote.py
def calibrate(self, p_screen, p_wii):
l = []
for i in range(0,4):
l.append( [p_wii[i][0], p_wii[i][1], 1, 0, 0, 0,
(-p_screen[i][0] * p_wii[i][0]),
(-p_screen[i][0] * p_wii[i][1])] )
l.append( [0, 0, 0, p_wii[i][0], p_wii[i][1], 1,
(-p_screen[i][1] * p_wii[i][0]),
(-p_screen[i][1] * p_wii[i][1])] )
A = matrix(l)
x = matrix( [
[p_screen[0][0]],
[p_screen[0][1]],
[p_screen[1][0]],
[p_screen[1][1]],
[p_screen[2][0]],
[p_screen[2][1]],
[p_screen[3][0]],
[p_screen[3][1]],
])
self.hCoefs = linalg.solve(A, x)
self.h11 = self.hCoefs[0]
self.h12 = self.hCoefs[1]
self.h13 = self.hCoefs[2]
self.h21 = self.hCoefs[3]
self.h22 = self.hCoefs[4]
self.h23 = self.hCoefs[5]
self.h31 = self.hCoefs[6]
self.h32 = self.hCoefs[7]
self.calibrationPoints = list(p_wii)
self.screenPoints = list(p_screen)
self.state = Wiimote.CALIBRATED
area_inside = calculateArea(self.calibrationPoints)
total_area = Wiimote.MAX_X * Wiimote.MAX_Y
self.utilization = float(area_inside)/float(total_area)
0
Example 49
Project: pydy Source File: test_ode_function_generator.py
def test_set_linear_system_solver(self):
g = ODEFunctionGenerator(self.sys.eom_method.forcing_full,
self.sys.coordinates,
self.sys.speeds,
self.sys.constants_symbols,
mass_matrix=self.sys.eom_method.mass_matrix_full)
assert g._solve_linear_system == np.linalg.solve
g = ODEFunctionGenerator(self.sys.eom_method.forcing_full,
self.sys.coordinates,
self.sys.speeds,
self.sys.constants_symbols,
mass_matrix=self.sys.eom_method.mass_matrix_full,
linear_sys_solver='numpy')
assert g._solve_linear_system == np.linalg.solve
g = ODEFunctionGenerator(self.sys.eom_method.forcing_full,
self.sys.coordinates,
self.sys.speeds,
self.sys.constants_symbols,
mass_matrix=self.sys.eom_method.mass_matrix_full,
linear_sys_solver='scipy')
assert g._solve_linear_system == sp.linalg.solve
solver = lambda A, b: np.dot(np.inv(A), b)
g = ODEFunctionGenerator(self.sys.eom_method.forcing_full,
self.sys.coordinates,
self.sys.speeds,
self.sys.constants_symbols,
mass_matrix=self.sys.eom_method.mass_matrix_full,
linear_sys_solver=solver)
assert g._solve_linear_system == solver
0
Example 50
Project: python-control Source File: statesp.py
def feedback(self, other=1, sign=-1):
"""Feedback interconnection between two LTI systems."""
other = _convertToStateSpace(other)
# Check to make sure the dimensions are OK
if ((self.inputs != other.outputs) or (self.outputs != other.inputs)):
raise ValueError("State space systems don't have compatible \
inputs/outputs for feedback.")
# Figure out the sampling time to use
if (self.dt == None and other.dt != None):
dt = other.dt # use dt from second argument
elif (other.dt == None and self.dt != None) or \
timebaseEqual(self, other):
dt = self.dt # use dt from first argument
else:
raise ValueError("Systems have different sampling times")
A1 = self.A
B1 = self.B
C1 = self.C
D1 = self.D
A2 = other.A
B2 = other.B
C2 = other.C
D2 = other.D
F = eye(self.inputs) - sign * D2 * D1
if matrix_rank(F) != self.inputs:
raise ValueError("I - sign * D2 * D1 is singular.")
# Precompute F\D2 and F\C2 (E = inv(F))
E_D2 = solve(F, D2)
E_C2 = solve(F, C2)
T1 = eye(self.outputs) + sign * D1 * E_D2
T2 = eye(self.inputs) + sign * E_D2 * D1
A = concatenate(
(concatenate(
(A1 + sign * B1 * E_D2 * C1, sign * B1 * E_C2), axis=1),
concatenate(
(B2 * T1 * C1, A2 + sign * B2 * D1 * E_C2), axis=1)),
axis=0)
B = concatenate((B1 * T2, B2 * D1 * T2), axis=0)
C = concatenate((T1 * C1, sign * D1 * E_C2), axis=1)
D = D1 * T2
return StateSpace(A, B, C, D, dt)