我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.linalg.solve()。
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 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 _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)
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 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 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 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
def _solve_topic(lo, hi, theta, alpha, gamma, XT, BTBpR, c0, c1, f, lam_d, lam_t): topic_batch = np.empty((hi - lo, f), dtype=theta.dtype) for ib, u in enumerate(range(lo, hi)): x_u, idx_u = get_row(XT, u) B_u = theta[idx_u] cpAT = gamma[u].dot(alpha.T) a = lam_d * x_u.dot(c1 * B_u) + lam_t * cpAT ''' 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 topic_batch[ib] = LA.solve(B, a) topic_batch = topic_batch.clip(0) return topic_batch
def _solve_gamma(lo, hi, alpha, beta, topic, M, B, c0, c1, f, lam_w, lam_t): gamma_batch = np.empty((hi - lo, f), dtype=alpha.dtype) for ib, i in enumerate(range(lo, hi)): t_i = topic[i,:] m_i, idx_m_i = get_row(M, i) B_i = beta[idx_m_i] ''' the reason why they put G_i in the loop instead of calculate GTG = gamma.T * gamma is that in the objective function, we currently only consider the non-zero elements in matrix W. ''' a = lam_t * np.dot(t_i, alpha) + lam_w * np.dot(m_i, B_i) gamma_batch[ib] = LA.solve(B, a) return gamma_batch
def _solve_beta(lo, hi, gamma, MT, B, f): beta_batch = np.empty((hi - lo, f), dtype=gamma.dtype) for ib, j in enumerate(range(lo, hi)): m_j, idx_m_j = get_row(MT, j) C_j = gamma[idx_m_j] a = np.dot(m_j, C_j) beta_batch[ib] = LA.solve(B, a) return beta_batch
def _solve_alpha(lo, hi, gamma, topicT, B, f): alpha_batch = np.empty((hi - lo, f), dtype=gamma.dtype) for ib, j in enumerate(range(lo, hi)): t_j = topicT[j,:] a = np.dot(t_j, gamma) alpha_batch[ib] = LA.solve(B, a) return alpha_batch
def _estimate_gradient(self, policy): env = self.environment var = self.var # store current policy parameter parameter = policy.parameters par_dim = policy.parameter_space.dimension # using forward differences trace = env.rollout(policy) j_ref = sum([x[2] for x in trace]) / len(trace) dj = np.zeros((2 * par_dim)) dv = np.append(np.eye(par_dim), -np.eye(par_dim), axis=0) dv *= var for n in range(par_dim): variation = dv[n] policy.parameters = parameter + variation trace_n = env.rollout(policy) jn = sum([x[2] for x in trace]) / len(trace_n) dj[n] = j_ref - jn grad = solve(dv.T.dot(dv), dv.T.dot(dj)) # reset current policy parameter policy.parameters = parameter return grad
def _estimate_gradient(self, policy): env = self.environment parameter = policy.parameters par_dim = policy.parameter_space.dimension dj = np.zeros((par_dim,)) dv = np.eye(par_dim) * self.var / 2 for n in range(par_dim): variation = dv[n] policy.parameters = parameter + variation trace_n = env.rollout(policy) policy.parameters = parameter - variation trace_n_ref = env.rollout(policy) jn = sum([x[2] for x in trace_n]) / len(trace_n) jn_ref = sum([x[2] for x in trace_n_ref]) / len(trace_n_ref) dj[n] = jn - jn_ref grad = solve(dv.T.dot(dv), dv.T.dot(dj)) policy.parameters = parameter return grad
def _omp(x, D, Gram, alpha, n_nonzero_coefs=None, tol=None): _, n_atoms = D.shape # the dict indexes of the atoms this datapoint uses Dx = np.array([]).astype(int) z = np.zeros(n_atoms) # the residual r = np.copy(x) i = 0 if n_nonzero_coefs is not None: tol = 1e-10 def cont_criterion(): not_reached_sparsity = i < n_nonzero_coefs return (not_reached_sparsity and norm(r) > tol) else: cont_criterion = lambda: norm(r) >= tol while (cont_criterion()): # find the atom that correlates the # most with the residual k = np.argmax(np.abs(alpha)) if k in Dx: break Dx = np.append(Dx, k) # solve the Least Squares problem # to find the coefs z DI = D[:, Dx] G = Gram[Dx, :][:, Dx] G = np.atleast_2d(G) try: G_inv = inv(G) except LinAlgError: print gram_singular_msg break z[Dx] = np.dot(G_inv, np.dot(D.T, x)[Dx]) r = x - np.dot(D[:, Dx], z[Dx]) alpha = np.dot(D.T, r) i += 1 return z
def llc(X, D, knn=5): # the sparse coder introduced in # "Locality-constrained Linear Coding for Image Classification" n_samples = X.shape[1] n_atoms = D.shape[1] # has the distance of # each sample to each atom dist = np.zeros((n_samples, n_atoms)) # calculate the distances for i in range(n_samples): for j in range(n_atoms): dist[i, j] = norm(X[:, i] - D[:, j]) # has the indices of the atoms # that are nearest neighbour to each sample knn_idx = np.zeros((n_samples, knn)).astype(int) for i in xrange(n_samples): knn_idx[i, :] = np.argsort(dist[i, :])[:knn] # the sparse coding matrix Z = np.zeros((n_atoms, n_samples)) II = np.eye(knn) beta = 1e-4 b = np.ones(knn) for i in xrange(n_samples): idx = knn_idx[i, :] z = D.T[idx, :] - repmat(X.T[i, :], knn, 1) C = np.dot(z, z.T) C = C + II * beta * np.trace(C) # solve the linear system C*c=b c = solve(C, b) # enforce the constraint on the sparse codes # such that sum(c)=1 c = c / float(np.sum(c)) Z[idx, i] = c return Z
def _refit(self, X): n_samples, n_features = X.shape for i in range(n_samples): 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 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)
def fvar(self): """Estimate the variances of the reduced free energies.""" # Shorthand indices - notation similiar to Ref. 1 n, m, mpk = self.total_samples, self.nstates_sampled, self.nstates k = mpk - m mask0, maskn0 = self.mask_zero, self.mask_nonzero _W_nj = self._W_nj O = _W_nj.T.dot(_W_nj) / n Os = O[:, :m] B1 = (hstack((Os.dot(self.PIs), zeros((mpk, k)))) - identity(mpk))[1:, 1:] A1 = (O - Os.dot(self.PIs).dot(Os.T))[1:, 1:] V = solve(B1, A1).dot(inv(B1.T)) / n # if any(diagonal(V) < 0.): # D = _W_nj.T.dot(self._R_nj) / n # A1 = (O - D.dot(self.PIs).dot(D.T))[1:, 1:] # V = solve(B1, A1).dot(inv(B1.T)) / n # Unshuffle the state indices. Note that the variance of state 0 is # zero by construction and thus is omitted from V - since the user # selected state 0 may be different from the lowest indexed sample # state, re-adjust so that the actual state 0 has zero variance. # V_full = zeros((mpk, mpk)) V_full[1:, 1:] = V var = zeros(mpk) var[maskn0] += diagonal(V_full)[:m] var[mask0] += diagonal(V_full)[m:] var += var[0] var[0] = 0.0 return var
def solve(self): # Solve LS self.VU = [node[key] for node in self.U.values() for key in ("ux","uy")] self.VF = [node[key] for node in self.F.values() for key in ("fx","fy")] 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] 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) for k,ic in enumerate(unknw): nd, var = self.index2key(ic) self.U[nd][var] = self.solved_u[k] # Updating nodes displacements for nd in self.nodes.values(): if np.isnan(nd.ux): nd.ux = self.U[nd.label]["ux"] if np.isnan(nd.uy): nd.uy = self.U[nd.label]["uy"] # For nodal forces/reactions self.NF = self.F.copy() self.VU = [node[key] for node in self.U.values() for key in ("ux","uy")] nf_calc = np.dot(self.KG, self.VU) for k in range(2*self.get_number_of_nodes()): nd, var = self.index2key(k, ("fx","fy")) self.NF[nd][var] = nf_calc[k] cnlab = np.floor(k/float(self.dof)) if var=="fx": self.nodes[cnlab].fx = nf_calc[k] elif var=="fy": self.nodes[cnlab].fy = nf_calc[k]
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] if len(unknw)==1: _k = unknw[0] _rowtmp = self.KG[_k,:] _ftmp = self.VF[_k] _fk = _ftmp - np.dot(np.delete(_rowtmp,_k), np.delete(self.VU,_k)) _uk = _fk / self.KG[_k, _k] # Then self.solved_u = np.array([_uk]) else: # "Normal" case self.K2S = np.delete(np.delete(self.KG,knw,0),knw,1) self.F2S = np.delete(self.VF,knw,0) self.solved_u = la.solve(self.K2S,self.F2S) # For displacements # 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 solve(self): # Solve LS self.VU = [node[key] for node in self.U.values() for key in ("uy","ur")] self.VF = [node[key] for node in self.F.values() for key in ("fy","m")] 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] 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) for k,ic in enumerate(unknw): nd, var = self.index2key(ic) self.U[nd][var] = self.solved_u[k] # Updating nodes displacements for nd in self.nodes.values(): if np.isnan(nd.uy): nd.uy = self.U[nd.label]["uy"] if np.isnan(nd.ur): nd.ur = self.U[nd.label]["ur"] # For nodal forces/reactions self.NF = self.F.copy() self.VU = [node[key] for node in self.U.values() for key in ("uy","ur")] nf_calc = np.dot(self.KG, self.VU) for k in range(2*self.get_number_of_nodes()): nd, var = self.index2key(k, ("fy","m")) self.NF[nd][var] = nf_calc[k] cnlab = np.floor(k/float(self.dof)) if var=="fy": self.nodes[cnlab].fy = nf_calc[k] elif var=="m": self.nodes[cnlab].m = nf_calc[k]
def solve(self): self._check_nodes() # Solve LS self.VU = [node[key] for node in self.U.values() for key in ("ux","uy")] self.VF = [node[key] for node in self.F.values() for key in ("fx","fy")] 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] self.K2S = np.delete(np.delete(self.KG,knw,0),knw,1) self.F2S = np.delete(self.VF,knw,0) # For displacements try: self.solved_u = la.solve(self.K2S,self.F2S) except: print("Solved using LSTSQ") self.solved_u = la.lstsq(self.K2S, self.F2S)[0] for k,ic in enumerate(unknw): nd, var = self.index2key(ic) self.U[nd][var] = self.solved_u[k] # Updating nodes displacements for nd in self.nodes.values(): if np.isnan(nd.ux): nd.ux = self.U[nd.label]["ux"] if np.isnan(nd.uy): nd.uy = self.U[nd.label]["uy"] # For nodal forces/reactions self.NF = self.F.copy() self.VU = [node[key] for node in self.U.values() for key in ("ux","uy")] nf_calc = np.dot(self.KG, self.VU) for k in range(2*self.get_number_of_nodes()): nd, var = self.index2key(k, ("fx","fy")) self.NF[nd][var] = nf_calc[k] cnlab = np.floor(k/float(self.dof)) if var=="fx": self.nodes[cnlab].fx = nf_calc[k] elif var=="fy": self.nodes[cnlab].fy = nf_calc[k]
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)
