Skip to content

Commit d4148bf

Browse files
fengdu78fengdu78
committed
修改稀疏矩阵部分
1 parent b24ff61 commit d4148bf

File tree

1 file changed

+104
-152
lines changed

1 file changed

+104
-152
lines changed

‎4.scipy/8.scipy-sparse.ipynb

+104-152
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
"cells": [
33
{
44
"cell_type": "code",
5-
"execution_count": 3,
5+
"execution_count": 1,
66
"metadata": {},
77
"outputs": [],
88
"source": [
@@ -29,7 +29,7 @@
2929
},
3030
{
3131
"cell_type": "code",
32-
"execution_count": 10,
32+
"execution_count": 2,
3333
"metadata": {},
3434
"outputs": [
3535
{
@@ -54,7 +54,7 @@
5454
},
5555
{
5656
"cell_type": "code",
57-
"execution_count": 11,
57+
"execution_count": 3,
5858
"metadata": {},
5959
"outputs": [
6060
{
@@ -79,7 +79,7 @@
7979
},
8080
{
8181
"cell_type": "code",
82-
"execution_count": 12,
82+
"execution_count": 4,
8383
"metadata": {},
8484
"outputs": [
8585
{
@@ -108,229 +108,181 @@
108108
"cell_type": "markdown",
109109
"metadata": {},
110110
"source": [
111-
"### 最短路径"
111+
"### 矩阵向量相乘"
112112
]
113113
},
114114
{
115115
"cell_type": "code",
116-
"execution_count": 14,
117-
"metadata": {},
118-
"outputs": [
119-
{
120-
"ename": "SyntaxError",
121-
"evalue": "invalid syntax (<ipython-input-14-adc4f84eb492>, line 3)",
122-
"output_type": "error",
123-
"traceback": [
124-
"\u001b[1;36m File \u001b[1;32m\"<ipython-input-14-adc4f84eb492>\"\u001b[1;36m, line \u001b[1;32m3\u001b[0m\n\u001b[1;33m digraph graph1 {\u001b[0m\n\u001b[1;37m ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
125-
]
126-
}
127-
],
128-
"source": [
129-
"\n",
130-
"#%figonly=使用稀疏矩阵可以表示有向图\n",
131-
"digraph graph1 {\n",
132-
" rankdir=LR;\n",
133-
" size=\"8,5\"\n",
134-
" node [shape = circle];\n",
135-
" A -> B [ label = \"10\" ];\n",
136-
" B -> C [ label = \"5\" ];\n",
137-
" A -> C [ label = \"3\" ];\n",
138-
" C -> D [ label = \"7\" ];\n",
139-
" D -> A [ label = \"4\" ];\n",
140-
" D -> C [ label = \"6\" ];\n",
141-
"}"
142-
]
143-
},
144-
{
145-
"cell_type": "code",
146-
"execution_count": 15,
116+
"execution_count": 5,
147117
"metadata": {},
148118
"outputs": [
149119
{
150120
"data": {
151121
"text/plain": [
152-
"matrix([[ 0., 10., 3., 0.],\n",
153-
" [ 0., 0., 5., 0.],\n",
154-
" [ 0., 0., 0., 7.],\n",
155-
" [ 4., 0., 6., 0.]])"
122+
"array([ 1, -3, -1], dtype=int32)"
156123
]
157124
},
158-
"execution_count": 15,
125+
"execution_count": 5,
159126
"metadata": {},
160127
"output_type": "execute_result"
161128
}
162129
],
163130
"source": [
164-
"w = sparse.dok_matrix((4, 4))\n",
165-
"\n",
166-
"edges = [(0, 1, 10), (1, 2, 5), (0, 2, 3),\n",
167-
" (2, 3, 7), (3, 0, 4), (3, 2, 6)]\n",
168-
"\n",
169-
"for i, j, v in edges:\n",
170-
" w[i, j] = v\n",
171-
"\n",
172-
"w.todense()"
131+
" import numpy as np\n",
132+
"from scipy.sparse import csr_matrix\n",
133+
"A = csr_matrix([[1, 2, 0], [0, 0, 3], [4, 0, 5]])\n",
134+
"v = np.array([1, 0, -1])\n",
135+
"A.dot(v)"
136+
]
137+
},
138+
{
139+
"cell_type": "markdown",
140+
"metadata": {},
141+
"source": [
142+
"#### 示例1"
143+
]
144+
},
145+
{
146+
"cell_type": "markdown",
147+
"metadata": {},
148+
"source": [
149+
"构造一个1000x1000 lil_matrix并添加值："
173150
]
174151
},
175152
{
176153
"cell_type": "code",
177-
"execution_count": 17,
154+
"execution_count": 6,
178155
"metadata": {},
179-
"outputs": [
180-
{
181-
"ename": "IndexError",
182-
"evalue": "index 0 is out of bounds for axis 0 with size 0",
183-
"output_type": "error",
184-
"traceback": [
185-
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
186-
"\u001b[1;31mIndexError\u001b[0m Traceback (most recent call last)",
187-
"\u001b[1;32m<ipython-input-17-c7c2ac13b02d>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m 5\u001b[0m \u001b[0mH\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mW\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mbimg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 6\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 7\u001b[1;33m \u001b[0mx0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwhere\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbimg\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mH\u001b[0m\u001b[1;33m//\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m==\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;31m#❷\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 8\u001b[0m \u001b[0mbimg\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mH\u001b[0m\u001b[1;33m//\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m:\u001b[0m\u001b[0mx0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 9\u001b[0m \u001b[0mbimg\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mH\u001b[0m\u001b[1;33m//\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
188-
"\u001b[1;31mIndexError\u001b[0m: index 0 is out of bounds for axis 0 with size 0"
189-
]
190-
}
191-
],
156+
"outputs": [],
192157
"source": [
193-
"img = pl.imread(\"maze.jpg\")\n",
194-
"sx, sy = (400, 979)\n",
195-
"ex, ey = (398, 25)\n",
196-
"bimg = np.all(img > 0.81, axis=2) #❶\n",
197-
"H, W = bimg.shape\n",
198-
"\n",
199-
"x0, x1 = np.where(bimg[H//2, :]==0)[0][[0, -1]] #❷\n",
200-
"bimg[H//2, :x0] = 0\n",
201-
"bimg[H//2, x1:] = 0"
158+
"from scipy.sparse import lil_matrix\n",
159+
"from scipy.sparse.linalg import spsolve\n",
160+
"from numpy.linalg import solve, norm\n",
161+
"from numpy.random import rand"
202162
]
203163
},
204164
{
205165
"cell_type": "code",
206-
"execution_count": 18,
166+
"execution_count": 7,
207167
"metadata": {},
208168
"outputs": [],
209169
"source": [
210-
"#上下相邻白色像素\n",
211-
"mask = (bimg[1:, :] & bimg[:-1, :]) \n",
212-
"idx = np.where(mask.ravel())[0]\n",
213-
"vedge = np.c_[idx, idx + W]\n",
214-
"pl.imsave(\"tmp.png\", mask, cmap=\"gray\")\n",
215-
"\n",
216-
"#左右相邻白色像素\n",
217-
"mask = (bimg[:, 1:] & bimg[:, :-1])\n",
218-
"y, x = np.where(mask)\n",
219-
"idx = y * W + x\n",
220-
"hedge = np.c_[idx, idx + 1]\n",
221-
"\n",
222-
"edges = np.vstack([vedge, hedge]) #❶\n",
223-
"\n",
224-
"values = np.ones(edges.shape[0])\n",
225-
"w = sparse.coo_matrix((values, (edges[:, 0], edges[:, 1])), #❷\n",
226-
" shape=(bimg.size, bimg.size))"
170+
"A = lil_matrix((1000, 1000))\n",
171+
"A[0, :100] = rand(100)\n",
172+
"A[1, 100:200] = A[0, :100]\n",
173+
"A.setdiag(rand(1000))"
174+
]
175+
},
176+
{
177+
"cell_type": "markdown",
178+
"metadata": {},
179+
"source": [
180+
"现在将其转换为CSR格式，并求解$A x = b$的$x$："
227181
]
228182
},
229183
{
230184
"cell_type": "code",
231-
"execution_count": 19,
185+
"execution_count": 8,
186+
"metadata": {},
187+
"outputs": [],
188+
"source": [
189+
"A = A.tocsr()\n",
190+
"b = rand(1000)\n",
191+
"x = spsolve(A, b)"
192+
]
193+
},
194+
{
195+
"cell_type": "markdown",
196+
"metadata": {},
197+
"source": [
198+
"将其转换为密集矩阵并求解，并检查结果是否相同："
199+
]
200+
},
201+
{
202+
"cell_type": "code",
203+
"execution_count": 9,
204+
"metadata": {},
205+
"outputs": [],
206+
"source": [
207+
"x_ = solve(A.toarray(), b)"
208+
]
209+
},
210+
{
211+
"cell_type": "markdown",
232212
"metadata": {},
233-
"outputs": [
234-
{
235-
"name": "stdout",
236-
"output_type": "stream",
237-
"text": [
238-
"(1, 801600)\n",
239-
"(1, 801600)\n"
240-
]
241-
}
242-
],
243213
"source": [
244-
"from scipy.sparse import csgraph\n",
245-
"startid = sy * W + sx\n",
246-
"endid = ey * W + ex\n",
247-
"d, p = csgraph.dijkstra(w, indices=[startid], return_predecessors=True, directed=False)\n",
248-
"print(d.shape)\n",
249-
"print(p.shape)"
214+
"现在我们可以使用以下公式计算错误的范数："
250215
]
251216
},
252217
{
253218
"cell_type": "code",
254-
"execution_count": 20,
219+
"execution_count": 10,
255220
"metadata": {},
256221
"outputs": [
257222
{
258223
"data": {
259224
"text/plain": [
260-
"110"
225+
"True"
261226
]
262227
},
263-
"execution_count": 20,
228+
"execution_count": 10,
264229
"metadata": {},
265230
"output_type": "execute_result"
266231
}
267232
],
268233
"source": [
269-
"np.isinf(d[0]).sum()"
234+
"err = norm(x-x_)\n",
235+
"err < 1e-10"
270236
]
271237
},
272238
{
273-
"cell_type": "code",
274-
"execution_count": 21,
239+
"cell_type": "markdown",
275240
"metadata": {},
276-
"outputs": [],
277241
"source": [
278-
"path = []\n",
279-
"node_id = endid\n",
280-
"while True:\n",
281-
" path.append(node_id)\n",
282-
" if node_id == startid or node_id < 0:\n",
283-
" break\n",
284-
" node_id = p[0, node_id]\n",
285-
"path = np.array(path)"
242+
"#### 示例2\n",
243+
"构造COO格式的矩阵："
286244
]
287245
},
288246
{
289247
"cell_type": "code",
290-
"execution_count": 22,
248+
"execution_count": 11,
291249
"metadata": {},
292-
"outputs": [
293-
{
294-
"ename": "ValueError",
295-
"evalue": "assignment destination is read-only",
296-
"output_type": "error",
297-
"traceback": [
298-
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
299-
"\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)",
300-
"\u001b[1;32m<ipython-input-22-8e148d29628c>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[1;31m#%figonly=用dijkstra计算最短路径\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 2\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0my\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpath\u001b[0m \u001b[1;33m%\u001b[0m \u001b[0mW\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpath\u001b[0m \u001b[1;33m//\u001b[0m \u001b[0mW\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0mimg\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 4\u001b[0m \u001b[0mfig\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxes\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpl\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfigsize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m16\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m12\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 5\u001b[0m \u001b[0maxes\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mimshow\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mimg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
301-
"\u001b[1;31mValueError\u001b[0m: assignment destination is read-only"
302-
]
303-
}
304-
],
250+
"outputs": [],
305251
"source": [
306-
"#%figonly=用dijkstra计算最短路径\n",
307-
"x, y = path % W, path // W\n",
308-
"img[y, x, :] = 0\n",
309-
"fig, axes = pl.subplots(1, 2, figsize=(16, 12))\n",
310-
"axes[0].imshow(img)\n",
311-
"axes[1].imshow(bimg, cmap=\"gray\")\n",
312-
"for ax in axes:\n",
313-
" ax.axis(\"off\")\n",
314-
"fig.subplots_adjust(0, 0, 1, 1, 0, 0)"
252+
"from scipy import sparse\n",
253+
"from numpy import array\n",
254+
"I = array([0,3,1,0])\n",
255+
"J = array([0,3,1,2])\n",
256+
"V = array([4,5,7,9])\n",
257+
"A = sparse.coo_matrix((V,(I,J)),shape=(4,4))"
315258
]
316259
},
317260
{
318261
"cell_type": "markdown",
319262
"metadata": {},
320263
"source": [
321-
"> **SOURCE**\n",
264+
"注意，索引不需要排序。\n",
322265
"\n",
323-
"> `scpy2.scipy.hrd_solver`使用`csgraph`计算华容道游戏【横刀立马】布局步数最少的解法。"
266+
"转换为CSR或CSC时，将对重复的（i，j）条目进行求和。"
324267
]
325268
},
326269
{
327270
"cell_type": "code",
328-
"execution_count": 1,
271+
"execution_count": 12,
329272
"metadata": {},
330273
"outputs": [],
331274
"source": [
332-
"#%hide\n",
333-
"%exec_python -m scpy2.scipy.hrd_solver"
275+
"I = array([0,0,1,3,1,0,0])\n",
276+
"J = array([0,2,1,3,1,0,0])\n",
277+
"V = array([1,1,1,1,1,1,1])\n",
278+
"B = sparse.coo_matrix((V,(I,J)),shape=(4,4)).tocsr()"
279+
]
280+
},
281+
{
282+
"cell_type": "markdown",
283+
"metadata": {},
284+
"source": [
285+
"这对于构造有限元刚度矩阵和质量矩阵是有用的。"
334286
]
335287
}
336288
],
@@ -350,7 +302,7 @@
350302
"name": "python",
351303
"nbconvert_exporter": "python",
352304
"pygments_lexer": "ipython3",
353-
"version": "3.6.4"
305+
"version": "3.7.3"
354306
}
355307
},
356308
"nbformat": 4,

0 commit comments

Comments
 (0)