我们从Python开源项目中，提取了以下41个代码示例，用于说明如何使用scipy.sparse.linalg.spsolve()。
def ac_solve(net): """ :param net: :return: """ net.conductance_matrix() net.dynamic_matrix() net.rhs_matrix() # frequency f = float(net.analysis[-1]) # linear system definition net.x = spsolve(net.G + 1j * 2 * np.pi * f* net.C, net.rhs)
def solve_linear(model:Model.fem_model): K_bar,F_bar,index=model.K_,model.F_,model.index Dvec=model.D Logger.info('Solving linear model with %d DOFs...'%model.DOF) n_nodes=model.node_count try: #sparse matrix solution delta_bar = sl.spsolve(sp.csr_matrix(K_bar),F_bar,sym_pos=True) delta = delta_bar #fill original displacement vector prev = 0 for idx in index: gap=idx-prev if gap>0: delta=np.insert(delta,prev,[0]*gap) prev = idx + 1 if idx==index[-1] and idx!=n_nodes-1: delta = np.insert(delta,prev, [0]*(n_nodes*6-prev)) delta += Dvec except Exception as e: print(e) return None model.is_solved=True return delta
def solve_linear2(model:Model.fem_model): K_bar,F_bar,index=model.K_,model.F_,model.index Dvec=model.D Logger.info('Solving linear model with %d DOFs...'%model.DOF) n_nodes=model.node_count #sparse matrix solution delta_bar = sl.spsolve(sp.csc_matrix(K_bar),F_bar) #delta_bar=linalg.solve(K_bar,F_bar,sym_pos=True) delta = delta_bar #fill original displacement vector prev = 0 for idx in index: gap=idx-prev if gap>0: delta=np.insert(delta,prev,[0]*gap) prev = idx + 1 if idx==index[-1] and idx!=n_nodes-1: delta = np.insert(delta,prev, [0]*(n_nodes*6-prev)) delta += Dvec model.is_solved=True return delta
def raise_order(self, amount): """ Raise the polynomial order of the curve. :param int amount: Number of times to raise the order :return: self """ if amount < 0: raise ValueError('Raise order requires a non-negative parameter') elif amount == 0: return # create the new basis newBasis = self.bases[0].raise_order(amount) # set up an interpolation problem. This is in projective space, so no problems for rational cases interpolation_pts_t = newBasis.greville() # parametric interpolation points (t) N_old = self.bases[0].evaluate(interpolation_pts_t) N_new = newBasis.evaluate(interpolation_pts_t, sparse=True) interpolation_pts_x = N_old * self.controlpoints # projective interpolation points (x,y,z,w) # solve the interpolation problem self.controlpoints = np.array(splinalg.spsolve(N_new, interpolation_pts_x)) self.bases = [newBasis] return self
def rebuild(self, p, n): """ Creates an approximation to this curve by resampling it using a uniform knot vector of order *p* with *n* control points. :param int p: Polynomial discretization order :param int n: Number of control points :return: A new approximate curve :rtype: Curve """ # establish uniform open knot vector knot = [0] * p + list(range(1, n - p + 1)) + [n - p + 1] * p basis = BSplineBasis(p, knot) # set parametric range of the new basis to be the same as the old one basis.normalize() t0 = self.bases[0].start() t1 = self.bases[0].end() basis *= (t1 - t0) basis += t0 # fetch evaluation points and solve interpolation problem t = basis.greville() N = basis.evaluate(t, sparse=True) controlpoints = splinalg.spsolve(N, self.evaluate(t)) # return new resampled curve return Curve(basis, controlpoints)
def interpolate(x, basis, t=None): """ Perform general spline interpolation on a provided basis. :param matrix-like x: Matrix *X[i,j]* of interpolation points *xi* with components *j* :param BSplineBasis basis: Basis on which to interpolate :param array-like t: parametric values at interpolation points; defaults to Greville points if not provided :return: Interpolated curve :rtype: Curve """ # wrap x into a numpy matrix x = np.matrix(x) # evaluate all basis functions in the interpolation points if t is None: t = basis.greville() N = basis.evaluate(t, sparse=True) # solve interpolation problem controlpoints = splinalg.spsolve(N, x) return Curve(basis, controlpoints)
def do_newton_fmg_cycle(self, prob, rhs, level, nu0, nu1, nu2, max_inner=1,max_outer=20): self.fh[level] = rhs current_ndofs = int(math.sqrt(rhs.shape[0])) mgrid = MyMultigrid(current_ndofs,int(np.log2(current_ndofs+1))) mgrid.attach_transfer(LinearTransfer2D) M = specificJacobi(current_ndofs,prob.gamma, np.ones(rhs.shape[0])) mgrid.attach_smoother(WeightedJacobi,M,omega=2.0/3.0) if (level < self.nlevels - 1): self.fh[level + 1] = mgrid.trans[0].restrict(self.fh[level]) self.vh[level + 1] = self.do_newton_fmg_cycle(prob,self.fh[level + 1], level + 1, nu0, nu1, nu2) else: self.vh[-1] = self.fh[-1]/M[0,0]#self.do_newton_cycle(prob,self.vh[-1], self.fh[-1], nu1, nu2, level, n_v_cycles=1, max_outer=1) #sLA.spsolve(M, self.fh[-1]) return self.vh[-1] for i in range(nu0): self.vh[level] = self.do_newton_cycle(prob,self.vh[level], self.fh[level], nu1, nu2, level, max_inner, max_outer)[0] return self.vh[level]
def do_fmg_cycle_recursive(self, rhs, h, nu0, nu1, nu2): # set intial conditions (note: resetting vectors here is important!) self.fh[0] = rhs # downward cycle if (h < self.nlevels - 1): self.fh[h + 1] = self.trans[h].restrict(self.fh[h]) # plt.plot(h, self.fh[h]) self.vh[h + 1] = self.do_fmg_cycle_recursive(self.fh[h + 1], h + 1, nu0, nu1, nu2) else: self.vh[-1] = sLA.spsolve(self.Acoarse, self.fh[-1]) return self.vh[-1] # correct self.vh[h] = self.trans[h].prolong(self.vh[h + 1]) for i in range(nu0): self.vh[h] = self.do_v_cycle(self.vh[h], self.fh[h], nu1, nu2, h) return self.vh[h]
def exact_is_fixpoint(smoo_class): ndofs = 127 prob = Poisson1D(ndofs) smoo = smoo_class(prob.A, 2.0/3.0) k = 6 xvalues = np.array([(i + 1) * prob.dx for i in range(prob.ndofs)]) rhs = np.sin(np.pi * k * xvalues) uex = spLA.spsolve(smoo.A, rhs) u = uex for i in range(10): u = smoo.smooth(rhs, u) err = np.linalg.norm(u - uex, np.inf) assert err <= 1E-14, 'Exact solution is not a fixpoint of the iteration!'
def computeX(template, L, H, D, landmarkAll, landmarkIndex): X = np.array(template.v) ptsNum = len(X) left = L.dot(L) right = L.dot(H) landmark3D = X[landmarkIndex] rowForP = np.array(map(lambda x:x, xrange(2*ptsNum))).repeat(3) colForP = (np.tile(np.array(map(lambda x:x, xrange(3*ptsNum))).reshape((ptsNum,3)),2)).flatten() for landmark2D in landmarkAll: P = calP(landmark2D[19:], landmark3D) valForP = np.tile(P[0:2,0:3].flatten(),ptsNum) Pplus = csr_matrix((valForP, (rowForP, colForP)), shape=(2*ptsNum, 3*ptsNum)) W = np.zeros((2*vCount,1)) for count, index in enumerate(landmarkIndex): W[2*index] = landmark2D[count, 0] - P[0, 3] W[2*index+1] = landmark2D[count, 1] - P[1, 3] tempL = Pplus.dot(D) left += lam*tempL.T.dot(tempL) right += lam*Pplus.T.dot(W) #return lsqr(left, np.array(right))[0] return spsolve(left, right)
def itera(template): #L = computeL(template) vertex = np.array(template.v) ptsNum = len(vertex) X = vertex.reshape((3*vCount,1)) landmark3D = vertex[landmarkIndex] sumL = L.dot(L) sumR = sumL.dot(X) rowForP = np.array(map(lambda x:x, xrange(2*ptsNum))).repeat(3) colForP = (np.tile(np.array(map(lambda x:x, xrange(3*ptsNum))).reshape((ptsNum,3)),2)).flatten() for landmark2D in landmarkAll: P = calP(landmark2D, landmark3D) valForP = np.tile(P[0:2,0:3].flatten(),ptsNum) Pplus = csr_matrix((valForP, (rowForP, colForP)), shape=(2*ptsNum, 3*ptsNum)) W = np.zeros((2*vCount,1)) for count, index in enumerate(landmarkIndex): W[2*index] = landmark2D[count, 0] - P[0, 3] W[2*index+1] = landmark2D[count, 1] - P[1, 3] tempL = Pplus.dot(D) sumL = sumL + lamda * tempL.T.dot(tempL) sumR = sumR + lamda * Pplus.T.dot(W) newV = spsolve(sumL, sumR) template.v = newV.reshape((len(newV)/3, 3)) return template
def iteration(self, user, fixed_vecs): num_solve = self.num_users if user else self.num_items num_fixed = fixed_vecs.shape[0] YTY = fixed_vecs.T.dot(fixed_vecs) eye = sparse.eye(num_fixed) lambda_eye = self.reg_param * sparse.eye(self.num_factors) solve_vecs = np.zeros((num_solve, self.num_factors)) for i in xrange(num_solve): if user: matrix_i = self.matrix[i].toarray() else: matrix_i = self.matrix[:, i].T.toarray() CuI = sparse.diags(matrix_i, [0]) pu = matrix_i.copy() pu[np.where(pu != 0)] = 1.0 YTCuIY = fixed_vecs.T.dot(CuI).dot(fixed_vecs) YTCupu = fixed_vecs.T.dot(CuI + eye).dot(sparse.csr_matrix(pu).T) xu = spsolve(YTY + YTCuIY + lambda_eye, YTCupu) solve_vecs[i] = xu return solve_vecs
def dc_solve(net): """ "dc_solve" solves DC network :return: * x, solution """ net.conductance_matrix() net.rhs_matrix() # linear system definition net.x = spsolve(net.G, net.rhs)
def linear_alg_solver(A, B): return linalg.spsolve(A, B, use_umfpack=True)
def solve(A, b): """ A generic linear solver for sparse and dense matrices Args: A : array_like, possibly sparse b : array_like Returns: x such that Ax = b """ if sparse.issparse(A): return spsolve(A, b) else: return np.linalg.solve(A, b)
def test_solve_identity_ones(dtype): b = np.ones(solver.N, dtype=dtype) mat = sp.eye(solver.N, format='csr', dtype=dtype) obs = solver.solve(mat, b) exp = spla.spsolve(mat, b) assert np.allclose(exp, obs)
def test_solve_identity_range(dtype): b = np.arange(solver.N, dtype=dtype) mat = sp.eye(solver.N, format='csr', dtype=dtype) obs = solver.solve(mat, b) exp = spla.spsolve(mat, b) assert np.allclose(exp, obs)
def test_solve_ones_ones(dtype): b = np.ones(solver.N, dtype=dtype) mat = solver.ones(dtype=dtype) + 9*sp.eye(solver.N, format='csr', dtype=dtype) obs = solver.solve(mat, b) exp = spla.spsolve(mat, b) assert np.allclose(exp, obs)
def test_solve_ones_range(dtype): b = np.arange(solver.N, dtype=dtype) mat = solver.ones(dtype=dtype) + 9*sp.eye(solver.N, format='csr', dtype=dtype) obs = solver.solve(mat, b) exp = spla.spsolve(mat, b) assert np.allclose(exp, obs)
def test_solve_range_range(dtype): b = np.arange(solver.N, dtype=dtype) mat = solver.ones(dtype=dtype) + sp.diags([b], offsets=[0], shape=(solver.N, solver.N), format='csr', dtype=dtype) obs = solver.solve(mat, b) exp = spla.spsolve(mat, b) assert np.allclose(exp, obs)
def propagate_information_filter_SLOW(x_analysis, P_analysis, P_analysis_inverse, M_matrix, Q_matrix): """Information filter state propagation using the INVERSER state covariance matrix and a linear state transition model. This function returns `None` for the forecast covariance matrix (as this takes forever). This method is based on the approximation to the inverse of the KF covariance matrix. Parameters ----------- x_analysis : array The analysis state vector. This comes either from the assimilation or directly from a previoulsy propagated state. P_analysis : 2D sparse array The analysis covariance matrix (typically will be a sparse matrix). As this is an information filter update, you will typically pass `None` to it, as it is unused. P_analysis_inverse : 2D sparse array The INVERSE analysis covariance matrix (typically a sparse matrix). M_matrix : 2D array The linear state propagation model. Q_matrix: 2D array (sparse) The state uncertainty inflation matrix that is added to the covariance matrix. Returns ------- x_forecast (forecast state vector), `None` and P_forecast_inverse (forecast inverse covariance matrix)""" logging.info("Starting the propagation...") x_forecast = M_matrix.dot(x_analysis) n, n = P_analysis_inverse.shape S= P_analysis_inverse.dot(Q_matrix) A = (sp.eye(n) + S).tocsc() P_forecast_inverse = spl.spsolve(A, P_analysis_inverse) logging.info("DOne with propagation") return x_forecast, None, P_forecast_inverse
def test_exact_is_fixpoint_of_vcycle(self): k = 6 xvalues = np.array([(i+1) * self.prob.dx for i in range(self.prob.ndofs)]) self.prob.rhs = (np.pi*k)**2 * np.sin(np.pi*k*xvalues) uex = spLA.spsolve(self.prob.A, self.prob.rhs) u = uex for i in range(10): u = self.mymg.do_v_cycle(u, self.prob.rhs, nu1=self.nu1, nu2=self.nu2, lstart=0) err = np.linalg.norm(u - uex, np.inf) assert err <= 1E-14, 'Exact solution is not a fixpoint of the V-cycle iteration!' + str(err)
def newton_raphson_sparse(f, guess, dfdx, x_tol=1e-10, lim_iter=100): """Solve f(x) = 0 with initial guess for x and dfdx(x). dfdx(x) should return a sparse Jacobian. Terminate if error on norm of f(x) is < x_tol or there were more than lim_iter iterations. """ converged = False n_iter = 0 F = f(guess) diff = norm(F,np.Inf) logger.debug("Error at iteration %d: %f", n_iter, diff) while diff > x_tol and n_iter < lim_iter: n_iter +=1 guess = guess - spsolve(dfdx(guess),F) F = f(guess) diff = norm(F,np.Inf) logger.debug("Error at iteration %d: %f", n_iter, diff) if diff > x_tol: logger.warn("Warning, we didn't reach the required tolerance within %d iterations, error is at %f. See the section \"Troubleshooting\" in the documentation for tips to fix this. ", n_iter, diff) elif not np.isnan(diff): converged = True return guess, n_iter, diff, converged
def calculate_PTDF(sub_network,skip_pre=False): """ Calculate the Power Transfer Distribution Factor (PTDF) for sub_network. Sets sub_network.PTDF as a (dense) numpy array. Parameters ---------- sub_network : pypsa.SubNetwork skip_pre: bool, default False Skip the preliminary steps of computing topology, calculating dependent values, finding bus controls and computing B and H. """ if not skip_pre: calculate_B_H(sub_network) #calculate inverse of B with slack removed n_pvpq = len(sub_network.pvpqs) index = np.r_[:n_pvpq] I = csc_matrix((np.ones(n_pvpq), (index, index))) B_inverse = spsolve(sub_network.B[1:, 1:],I) #exception for two-node networks, where B_inverse is a 1d array if issparse(B_inverse): B_inverse = B_inverse.toarray() elif B_inverse.shape == (1,): B_inverse = B_inverse.reshape((1,1)) #add back in zeroes for slack B_inverse = np.hstack((np.zeros((n_pvpq,1)),B_inverse)) B_inverse = np.vstack((np.zeros(n_pvpq+1),B_inverse)) sub_network.PTDF = sub_network.H*B_inverse
def predictor(wh, dt): """ Returns the Galerkin predictor, given the WENO reconstruction at tn """ nx, ny, nz, = wh.shape[:3] wh = wh.reshape([nx, ny, nz, (N+1)**ndim, n]) qh = zeros([nx, ny, nz, NT, n]) for i, j, k in product(range(nx), range(ny), range(nz)): w = wh[i, j, k] Ww = dot(W, w) if hidalgo: q = hidalgo_initial_guess(w, dt*gaps) else: q = standard_initial_guess(w) if stiff: func = lambda X: X - spsolve(U, rhs(X, Ww, dt)) qh[i, j, k] = newton_krylov(func, q, f_tol=TOL, method='bicgstab') else: for count in range(MAX_ITER): qNew = spsolve(U, rhs(q, Ww, dt)) if isnan(qNew).any(): print("DG root finding failed") break elif (absolute(q-qNew) > TOL * (1 + absolute(q))).any(): # Check convergence q = qNew continue else: qh[i, j, k] = qNew break else: print("Maximum iterations reached") return qh
def test_basic_spsolve_vector(): ps.remove_stored_factorization() ps.free_memory() A, b = create_test_A_b_rand() xpp = spsolve(A,b) xscipy = scipyspsolve(A,b) np.testing.assert_array_almost_equal(xpp, xscipy)
def test_basic_spsolve_matrix(): ps.remove_stored_factorization() ps.free_memory() A, b = create_test_A_b_rand(matrix=True) xpp = spsolve(A,b) xscipy = scipyspsolve(A,b) np.testing.assert_array_almost_equal(xpp, xscipy)
def linearMethod(self): """ Determines ``pi`` by solving a system of linear equations using :func:`spsolve`. The method has no parameters since it is an exact method. The result is stored in the class attribute ``pi``. Example ------- >>> P = np.array([[0.5,0.5],[0.6,0.4]]) >>> mc = markovChain(P) >>> mc.linearMethod() >>> print(mc.pi) [ 0.54545455 0.45454545] Remarks ------- For large state spaces, the linear algebra solver may not work well due to memory overflow. Code due to http://stackoverflow.com/questions/21308848/ """ P = self.getIrreducibleTransitionMatrix() #if P consists of one element, then set self.pi = 1.0 if P.shape == (1, 1): self.pi = np.array([1.0]) return size = P.shape[0] dP = P - eye(size) #Replace the first equation by the normalizing condition. A = vstack([np.ones(size), dP.T[1:,:]]).tocsr() rhs = np.zeros((size,)) rhs[0] = 1 self.pi = spsolve(A, rhs)
def __init__(self, training_points, training_values, num_leaves=2, n=5, comp=2): """ Initialize all attributes. Parameters ---------- training_points : ndarray training_values : ndarray num_leaves : int n : int comp : int """ super(RBFInterpolator, self).__init__(training_points, training_values, num_leaves) if self._ntpts < n: raise ValueError('RBFInterpolator only given {0} training points, but requested n={1}.' .format(self._ntpts, n)) # Comp is an arbitrary value that picks a function to use self.comp = comp # For weights, first find the training points radial neighbors tdist, tloc = self._KData.query(self._tp, n) Tt = tdist[:, :-1] / tdist[:, -1:] # Next determine weight matrix Rt = self._find_R(self._ntpts, Tt, tloc) weights = (spsolve(csc_matrix(Rt), self._tv))[..., np.newaxis] self.N = n self.weights = weights
def solve(a, b, silent=False, **kwargs): """Wrapper for spsolve removing null columns The null columns of matrix ``a`` is removed such and the linear system of equations is solved. The corresponding values of the solution ``x`` where the columns are null will also be null values. Parameters ---------- a : ndarray or sparse matrix A square matrix that will be converted to CSR form in the solution. b : scipy sparse matrix The matrix or vector representing the right hand side of the equation. silent : bool, optional A boolean to tell whether the msg messages should be printed. kwargs : keyword arguments, optional Other arguments directly passed to :func:`spsolve`. Returns ------- x : ndarray or sparse matrix The solution of the sparse linear equation. If ``b`` is a vector, then ``x`` is a vector of size ``a.shape[1]``. If ``b`` is a sparse matrix, then ``x`` is a matrix of size ``(a.shape[1], b.shape[1])``. """ a, used_cols = remove_null_cols(a, silent=silent) px = spsolve(a, b[used_cols], **kwargs) x = np.zeros(b.shape[0], dtype=b.dtype) x[used_cols] = px return x
def __call__(self, x, gradient, hessian, support=None): """ :param np.ndarray x: current estimate of the parameter :param np.ndarray gradient: parameter gradient. :param np.ndarray hessian: parameter Hessian. :param tuple(float) support: the bounds of the variable being updated, used to truncate the updated value :return: the new estimate after moving in the direction of the Newton step. :rtype: np.ndarray """ if hessian is None: raise ValueError('Hessian required for second order methods') else: if np.isscalar(hessian): step_vec = -gradient / hessian elif isinstance(hessian, np.ndarray): # dense matrix if hessian.size == gradient.size: # assume Hessian diagonal is stored step_vec = -gradient / np.asarray(hessian) else: step_vec = -np.linalg.solve(hessian, gradient.ravel(order=self.ravel_order)) else: # sparse matrix if hessian.shape[0] == 1: # sp.linalg.spsolve cannot handle 1D matrices step_vec = -gradient / hessian.toarray() else: step_vec = -spsolve(hessian, gradient.ravel(order=self.ravel_order)) self.step = step_vec.reshape(x.shape, order=self.ravel_order) value = x + self.step_size * self.step if np.any(~np.isfinite(value)): raise RuntimeError("Newly computed values are not all finite!") if support is not None: np.clip(value, support[0], support[1], out=value) return value
def calculate(model): boundaries = model['boundaries'] nodes = model['nodes'] unfixed = unfixed_coos(keys(nodes), values(boundaries)) coo_indexes = index_dict(unfixed) section_keys = 'Ax', 'Iz', 'Iy', 'Ay', 'Az', 'theta', 'J' sections = calculated_sections(items(model['sections']), section_keys) material_keys = 'E', 'G' materials = calculated_materials(items(model['materials']), material_keys) a = dok_matrix((len(coo_indexes),) * 2) for ln in values(model['lines']): v = line_vector(*line_nodes(ln, nodes)) p = dict(sections[ln['section']], **materials[ln['material']]) for k, (n1, n2) in zip(line.stiffness_global(*v, **p), stiffness_node_ids(ln)): for i, row in get_indexes(n1, coo_indexes): for j, col in get_indexes(n2, coo_indexes): a[row, col] += k[i][j] b = zeros(len(coo_indexes)) for ld in values(model['nodeloads']): for i, row in get_indexes(ld['node'], coo_indexes): if ld[coos[i]]: b[row] += ld[coos[i]] dis_array = spsolve(a, b) dis = {} for (node_id, coo), i in coo_indexes.items(): if node_id in dis: dis[node_id][coo] = dis_array[i] else: dis[node_id] = {coo: dis_array[i]} return { 'displacements': dis }
def solve(self, b, u=None, axis=0): """Solve matrix system Au = b where A is the current matrix (self) args: b (input/output) Vector of right hand side on entry. Solution on exit unless u is provided. u (output) Optional output vector Vectors may be one- or multidimensional. """ assert self.shape[0] == self.shape[1] assert self.shape[0] == b.shape[axis] if u is None: u = b else: assert u.shape == b.shape # Roll relevant axis to first if axis > 0: u = np.moveaxis(u, axis, 0) if not u is b: b = np.moveaxis(b, axis, 0) if b.ndim == 1: u[:] = spsolve(self.diags(), b) else: N = b.shape[0] P = np.prod(b.shape[1:]) u[:] = spsolve(self.diags(), b.reshape((N, P))).reshape(u.shape) if axis > 0: u = np.moveaxis(u, 0, axis) if not u is b: b = np.moveaxis(b, 0, axis) return u
def main(): """Main demo.""" # Load survey data llh, data = get_flightlines() # llh in n*3 lon-lat-hei format # data in n*3 k th u format # Crop to ROI ROI = (120.4, 120.5, -27.4, -27.3) keep = ((llh[:, 0] > ROI[0]) * (llh[:, 0] < ROI[1]) * (llh[:, 1] > ROI[2]) * (llh[:, 1] < ROI[3])) llh = llh[keep] data = data[keep] # Change reference frame to local North/East/Up in metres frame = LocalFrame(ROI) sensor_xyz = frame.to_xyz(llh) sensor_range = 1000 # metres sensor_scale = 1e2 # sensor parameter. Unknown? sensor_tol = 1e-12 # Sensor noise level pre-scaling noise_var = sensor_scale**2 * 1e-9 # known or learn? res = 20 # Output Grid resolution ground = make_ground(sensor_xyz, pad=2*sensor_range, res=res, frame=frame) n_grid = np.prod(ground.shape[:2]) n_sensor = sensor_xyz.shape[0] sensor_gain = sensor_scale * res ** 2 S = sensor_gain * sensitivity_matrix(ground, sensor_xyz, sensor_range, sensor_tol) # Investigate the sensor itself - we assume a white noise spatial prior # for now... G = sparse.eye(n_grid) # Ground sparsity? K = S.dot(G.dot(S.T)) + noise_var * sparse.eye(n_sensor) y = data[:, 2] # uranium column mu = S.T.dot(spsolve(K, y)).reshape(ground.shape[:2]) # Display the prediction pl.figure() # its not quite 2d but we can approximately show as a flat image gx = np.mean(ground[:, :, 0], axis=1) gy = np.mean(ground[:, :, 1], axis=0) pl.imshow(mu.T, interpolation='none', extent=(gx[0], gx[-1], gy[0], gy[-1])) # Display the points used in the calculation fig = pl.figure() ax = fig.add_subplot(111, projection='3d') ax.plot(ground[::2, ::2, 0].ravel(), ground[::2, ::2, 1].ravel(), 'b.', zs=ground[::2, ::2, 2].ravel()) ax.plot(sensor_xyz[:, 0], sensor_xyz[:, 1], 'k.', zs=sensor_xyz[:, 2]) ax.set_zlim((-200, 250)) ax.set_xlabel('Eastings (m)') ax.set_ylabel('Northings (m)') pl.axis('equal') pl.show() exit()
def do_v_cycle(self, v0, rhs, nu1, nu2, lstart): """Straightforward implementation of a V-cycle This can also be used inside an FMG-cycle! Args: v0 (numpy.array): initial values on finest level rhs (numpy.array): right-hand side on finest level nu1 (int): number of downward smoothing steps nu2 (int): number of upward smoothing steps lstart (int): starting level Returns: numpy.array: solution vector on finest level """ assert self.nlevels >= lstart >= 0 assert v0.size == self.vh[lstart].size # set intial conditions (note: resetting vectors here is important!) self.reset_vectors(lstart) self.vh[lstart] = v0 self.fh[lstart] = rhs # downward cycle for l in range(lstart, self.nlevels-1): # print('V-down: %i -> %i' %(l,l+1)) # pre-smoothing for i in range(nu1): self.vh[l] = self.smoo[l].smooth(self.fh[l], self.vh[l]) # restrict self.fh[l+1] = self.trans[l].restrict(self.fh[l] - self.smoo[l].A.dot(self.vh[l])) # solve on coarsest level self.vh[-1] = sLA.spsolve(self.Acoarse, self.fh[-1]) # upward cycle for l in reversed(range(lstart, self.nlevels-1)): # print('V-up: %i -> %i' %(l+1,l)) # correct self.vh[l] += self.trans[l].prolong(self.vh[l+1]) # post-smoothing for i in range(nu2): self.vh[l] = self.smoo[l].smooth(self.fh[l], self.vh[l]) return self.vh[lstart]
def do_v_cycle_recursive(self, v0, rhs, nu1, nu2, level): """Recursive implementation of a V-cycle This can also be used inside an FMG-cycle! Args: v0 (numpy.array): initial values on finest level rhs (numpy.array): right-hand side on finest level nu1 (int): number of downward smoothing steps nu2 (int): number of upward smoothing steps level (int): current level Returns: numpy.array: solution vector on current level """ assert self.nlevels > level >= 0 assert v0.size == self.vh[level].size # set intial conditions self.vh[level] = v0 self.fh[level] = rhs # downward cycle if level < self.nlevels-1: # pre-smoothing for i in range(nu1): self.vh[level] = self.smoo[level].smooth(self.fh[level], self.vh[level]) # restrict self.fh[level+1] = self.trans[level].restrict(self.fh[level] - self.smoo[level].A.dot(self.vh[level])) # recursive call to v-cycle self.vh[level+1] = self.do_v_cycle_recursive(np.zeros(self.vh[level+1].size), self.fh[level+1], nu1, nu2, level+1) # on coarsest level else: # solve on coarsest level self.vh[level] = sLA.spsolve(self.Acoarse, self.fh[level]) return self.vh[level] # correct self.vh[level] += self.trans[level].prolong(self.vh[level+1]) # post-smoothing for i in range(nu2): self.vh[level] = self.smoo[level].smooth(self.fh[level], self.vh[level]) return self.vh[level]
def solve(A, b, method, tol=1e-3): """ General sparse solver interface. method can be one of - spsolve_umfpack_mmd_ata - spsolve_umfpack_colamd - spsolve_superlu_mmd_ata - spsolve_superlu_colamd - bicg - bicgstab - cg - cgs - gmres - lgmres - minres - qmr - lsqr - lsmr """ if method == 'spsolve_umfpack_mmd_ata': return spla.spsolve(A,b,use_umfpack=True, permc_spec='MMD_ATA') elif method == 'spsolve_umfpack_colamd': return spla.spsolve(A,b,use_umfpack=True, permc_spec='COLAMD') elif method == 'spsolve_superlu_mmd_ata': return spla.spsolve(A,b,use_umfpack=False, permc_spec='MMD_ATA') elif method == 'spsolve_superlu_colamd': return spla.spsolve(A,b,use_umfpack=False, permc_spec='COLAMD') elif method == 'bicg': res = spla.bicg(A,b,tol=tol) return res[0] elif method == 'bicgstab': res = spla.bicgstab(A,b,tol=tol) return res[0] elif method == 'cg': res = spla.cg(A,b,tol=tol) return res[0] elif method == 'cgs': res = spla.cgs(A,b,tol=tol) return res[0] elif method == 'gmres': res = spla.gmres(A,b,tol=tol) return res[0] elif method == 'lgmres': res = spla.lgmres(A,b,tol=tol) return res[0] elif method == 'minres': res = spla.minres(A,b,tol=tol) return res[0] elif method == 'qmr': res = spla.qmr(A,b,tol=tol) return res[0] elif method == 'lsqr': res = spla.lsqr(A,b,atol=tol,btol=tol) return res[0] elif method == 'lsmr': res = spla.lsmr(A,b,atol=tol,btol=tol) return res[0] else: raise Exception('UnknownSolverType')
def itera(template): #L = calL(template) vertex = np.array(template.v) X = vertex.reshape((3*vCount,1))#3p vector #X = X0 landmark3D = vertex[landmarkIndex] imgList = os.listdir(imgSetDir) Pset = [] Wset = [] pMatrix = [] for im in imgList: imgPath = os.path.join(imgSetDir, im) landmark2D = landmarkFromFacepp(imgPath)[19:] #landmark2D = landmarkFromFacepp(imgPath) P = calP(landmark2D, landmark3D) pMatrix.append(P) W = np.zeros((2*vCount,1)) Pplus = np.zeros((2*vCount, 3*vCount)) count = 0 for index in landmarkIndex: Pplus[2*index:2*index+2, 3*index:3*index+3] = P[0:2,0:3] W[2*index] = landmark2D[count, 0] - P[0, 3] W[2*index+1] = landmark2D[count, 1] - P[1, 3] #W[2*index] = landmark2D[count, 0] #W[2*index+1] = landmark2D[count, 1] count = count + 1 Pplus = csr_matrix(Pplus) W = csr_matrix(W) Pset.append(Pplus) Wset.append(W) sumL = L.dot(L) sumR = sumL.dot(X) costVal2 = 0 for i in range(len(Pset)): #sumL = sumL + D.dot(Pset[i].T).dot(Pset[i]).dot(D) tempL = Pset[i].dot(D) costVal2 = costVal2 + np.linalg.norm(tempL.dot(X) - Wset[i]) sumL = sumL + lamda * tempL.T.dot(tempL) sumR = sumR + lamda * (Pset[i].T).dot(Wset[i]) costV2.append(costVal2) newV = spsolve(sumL, sumR) template.v = newV.reshape((len(newV)/3, 3)) return template, pMatrix, Wset, Pset #def main():
def solve_helmholtz(degree, resolution, analytic=False, return_error=False): """Solve a model Helmholtz problem on a unit square mesh with ``resolution`` elements in each direction, using equispaced Lagrange elements of degree ``degree``.""" # Set up the mesh, finite element and function space required. mesh = UnitSquareMesh(resolution, resolution) fe = LagrangeElement(mesh.cell, degree) fs = FunctionSpace(mesh, fe) # Create a function to hold the analytic solution for comparison purposes. analytic_answer = Function(fs) analytic_answer.interpolate(lambda x: cos(4*pi*x[0])*x[1]**2*(1.-x[1])**2) # If the analytic answer has been requested then bail out now. if analytic: return analytic_answer, 0.0 # Create the right hand side function and populate it with the # correct values. f = Function(fs) f.interpolate(lambda x: ((16*pi**2 + 1)*(x[1] - 1)**2*x[1]**2 - 12*x[1]**2 + 12*x[1] - 2) * cos(4*pi*x[0])) # Assemble the finite element system. A, l = assemble(fs, f) # Create the function to hold the solution. u = Function(fs) # Cast the matrix to a sparse format and use a sparse solver for # the linear system. This is vastly faster than the dense # alternative. A = sp.csr_matrix(A) u.values[:] = splinalg.spsolve(A, l) # Compute the L^2 error in the solution for testing purposes. error = errornorm(analytic_answer, u) if return_error: u.values -= analytic_answer.values # Return the solution and the error in the solution. return u, error
def solve_poisson(degree, resolution, analytic=False, return_error=False): """Solve a model Poisson problem on a unit square mesh with ``resolution`` elements in each direction, using equispaced Lagrange elements of degree ``degree``.""" # Set up the mesh, finite element and function space required. mesh = UnitSquareMesh(resolution, resolution) fe = LagrangeElement(mesh.cell, degree) fs = FunctionSpace(mesh, fe) # Create a function to hold the analytic solution for comparison purposes. analytic_answer = Function(fs) analytic_answer.interpolate(lambda x: sin(4*pi*x[0])*x[1]**2*(1.-x[1])**2) # If the analytic answer has been requested then bail out now. if analytic: return analytic_answer, 0.0 # Create the right hand side function and populate it with the # correct values. f = Function(fs) f.interpolate(lambda x: (16*pi**2*(x[1] - 1)**2*x[1]**2 - 2*(x[1] - 1)**2 - 8*(x[1] - 1)*x[1] - 2*x[1]**2) * sin(4*pi*x[0])) # Assemble the finite element system. A, l = assemble(fs, f) # Create the function to hold the solution. u = Function(fs) # Cast the matrix to a sparse format and use a sparse solver for # the linear system. This is vastly faster than the dense # alternative. A = sp.csr_matrix(A) u.values[:] = splinalg.spsolve(A, l) # Compute the L^2 error in the solution for testing purposes. error = errornorm(analytic_answer, u) if return_error: u.values -= analytic_answer.values # Return the solution and the error in the solution. return u, error
def CRAM16(A, n0, dt): """ Chebyshev Rational Approximation Method, order 16 Algorithm is the 16th order Chebyshev Rational Approximation Method, implemented in the more stable incomplete partial fraction (IPF) form [cram16]_. .. [cram16] Pusa, Maria. "Higher-Order Chebyshev Rational Approximation Method and Application to Burnup Equations." Nuclear Science and Engineering 182.3 (2016). Parameters ---------- A : scipy.linalg.csr_matrix Matrix to take exponent of. n0 : numpy.array Vector to operate a matrix exponent on. dt : float Time to integrate to. Returns ------- numpy.array Results of the matrix exponent. """ alpha = np.array([+2.124853710495224e-16, +5.464930576870210e+3 - 3.797983575308356e+4j, +9.045112476907548e+1 - 1.115537522430261e+3j, +2.344818070467641e+2 - 4.228020157070496e+2j, +9.453304067358312e+1 - 2.951294291446048e+2j, +7.283792954673409e+2 - 1.205646080220011e+5j, +3.648229059594851e+1 - 1.155509621409682e+2j, +2.547321630156819e+1 - 2.639500283021502e+1j, +2.394538338734709e+1 - 5.650522971778156e+0j], dtype=np.complex128) theta = np.array([+0.0, +3.509103608414918 + 8.436198985884374j, +5.948152268951177 + 3.587457362018322j, -5.264971343442647 + 16.22022147316793j, +1.419375897185666 + 10.92536348449672j, +6.416177699099435 + 1.194122393370139j, +4.993174737717997 + 5.996881713603942j, -1.413928462488886 + 13.49772569889275j, -10.84391707869699 + 19.27744616718165j], dtype=np.complex128) n = A.shape[0] alpha0 = 2.124853710495224e-16 k = 8 y = np.array(n0, dtype=np.float64) for l in range(1, k+1): y = 2.0*np.real(alpha[l]*sla.spsolve(A*dt - theta[l]*sp.eye(n), y)) + y y *= alpha0 return y
