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))
python类solve()的实例源码
transformation_tests_func.py 文件源码
项目：3D_Dense_Transformer_Networks
作者: JohnYC1995
项目源码
文件源码
阅读 27
收藏 0
点赞 0
评论 0
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)
test_linalg.py 文件源码
项目：PyDataLondon29-EmbarrassinglyParallelDAWithAWSLambda
作者: SignalMedia
项目源码
文件源码
阅读 28
收藏 0
点赞 0
评论 0
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
