Skip to content

Commit 0a0e151

Browse files
mathematicalmichaelmathematicalmichael
authored
use SVD to craft QoI (#51)
* new notebook * minor changes * run notebook * clean notbeook
1 parent 1580c30 commit 0a0e151

File tree

4 files changed

+396
-3
lines changed

4 files changed

+396
-3
lines changed

‎.gitignore‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ __pycache__/*
1212
.cache/*
1313
.*.swp
1414
*/.ipynb_checkpoints/*
15+
.ipynb_checkpoints/*
1516
*/*/.ipynb_checkpoints/*
1617
.DS_Store
1718

‎notebooks/MUD_SVD.ipynb‎

Lines changed: 388 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,388 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"id": "df13b096-f876-44d9-b80b-10d303e02793",
7+
"metadata": {},
8+
"outputs": [],
9+
"source": [
10+
"import numpy as np\n",
11+
"import scipy.stats.distributions as ds\n",
12+
"import matplotlib.pyplot as plt\n",
13+
"\n",
14+
"from mud.base import DensityProblem\n",
15+
"\n",
16+
"import mud_examples.poisson as ps"
17+
]
18+
},
19+
{
20+
"cell_type": "code",
21+
"execution_count": null,
22+
"id": "f1f73fca-73b5-429d-87ae-273204dc22df",
23+
"metadata": {},
24+
"outputs": [],
25+
"source": [
26+
"sample_dist = \"u\"\n",
27+
"input_dim = 2\n",
28+
"prefix = 1000\n",
29+
"fdir = f\"pde_{input_dim}D\"\n",
30+
"fname = f\"data/{fdir}/ref_{prefix}_{input_dim}{sample_dist}.pkl\"\n",
31+
"\n",
32+
"P = ps.pdeProblem(fname)\n",
33+
"P.load()\n",
34+
"if sample_dist == \"n\":\n",
35+
" P.dist = ds.norm # required for generator functions\n",
36+
" loc = -2.0\n",
37+
" scale = 0.2\n",
38+
"else:\n",
39+
" P.dist = ds.uniform\n",
40+
" loc = -4.0\n",
41+
" scale = 4.0"
42+
]
43+
},
44+
{
45+
"cell_type": "code",
46+
"execution_count": null,
47+
"id": "c25cda4f-3a26-4c8b-a9d6-4d4b7ae48923",
48+
"metadata": {},
49+
"outputs": [],
50+
"source": [
51+
"print(fname)"
52+
]
53+
},
54+
{
55+
"cell_type": "code",
56+
"execution_count": null,
57+
"id": "041017c4-ac99-4a65-8c4d-0bfc5746c4d2",
58+
"metadata": {},
59+
"outputs": [],
60+
"source": [
61+
"num_samples, num_max_qoi = P.qoi.shape\n",
62+
"print(f\"There are {num_max_qoi} sensors available and {num_samples} samples of {input_dim}-D parameter space.\")"
63+
]
64+
},
65+
{
66+
"cell_type": "code",
67+
"execution_count": null,
68+
"id": "36fdea18-087b-46a3-af14-db51a4a0921a",
69+
"metadata": {},
70+
"outputs": [],
71+
"source": [
72+
"# modify this method for new plot based on singular vectors\n",
73+
"num_qoi = 2\n",
74+
"ps.plot_without_fenics(\n",
75+
" fname,\n",
76+
" num_sensors=100,\n",
77+
" mode=\"hor\",\n",
78+
" num_qoi=num_qoi,\n",
79+
" example=\"mud\", # TODO: rename this variable\n",
80+
")"
81+
]
82+
},
83+
{
84+
"cell_type": "markdown",
85+
"id": "1af0b9c2-1586-4ca5-91ce-e02a178279e9",
86+
"metadata": {},
87+
"source": [
88+
"# Baseline Solution"
89+
]
90+
},
91+
{
92+
"cell_type": "code",
93+
"execution_count": null,
94+
"id": "7d4c5999-c869-48f2-b21a-3cd09c5e6217",
95+
"metadata": {},
96+
"outputs": [],
97+
"source": [
98+
"# generator function which takes (num_obs, sd) as arguments and returns a mud_problem\n",
99+
"mud_generator = P.mud_vector_horizontal(num_qoi, loc=loc, scale=scale)"
100+
]
101+
},
102+
{
103+
"cell_type": "code",
104+
"execution_count": null,
105+
"id": "2e1fa06a-9ef0-4e24-8172-85e0d6f3e003",
106+
"metadata": {},
107+
"outputs": [],
108+
"source": [
109+
"mud_problem = mud_generator(500, 0.01) # simulates noisy measurements"
110+
]
111+
},
112+
{
113+
"cell_type": "code",
114+
"execution_count": null,
115+
"id": "b2df0883-0c3a-4514-940b-c533ee40acf4",
116+
"metadata": {},
117+
"outputs": [],
118+
"source": [
119+
"mud_problem"
120+
]
121+
},
122+
{
123+
"cell_type": "code",
124+
"execution_count": null,
125+
"id": "8d5c8654-c7bc-4fb5-849d-700f75cb82ed",
126+
"metadata": {},
127+
"outputs": [],
128+
"source": [
129+
"mud_problem.mud_point()"
130+
]
131+
},
132+
{
133+
"cell_type": "markdown",
134+
"id": "5fc70123-6aec-4bba-9326-5eed435abc74",
135+
"metadata": {},
136+
"source": [
137+
"# Reference Solution"
138+
]
139+
},
140+
{
141+
"cell_type": "code",
142+
"execution_count": null,
143+
"id": "579f2d79-4f16-4d89-ace8-da0e02e3cc8a",
144+
"metadata": {},
145+
"outputs": [],
146+
"source": [
147+
"# TODO : turn this snippet into a method\n",
148+
"closest_fit_index_out = np.argmin(\n",
149+
" np.linalg.norm(P.qoi - np.array(P.qoi_ref), axis=1)\n",
150+
")\n",
151+
"g_projected = P.lam[closest_fit_index_out, :].ravel()\n",
152+
"lam_true = g_projected\n",
153+
"print(lam_true)"
154+
]
155+
},
156+
{
157+
"cell_type": "markdown",
158+
"id": "aeafe0e4-a777-4aee-a028-9305e4ece09c",
159+
"metadata": {},
160+
"source": [
161+
"# SVD Approach"
162+
]
163+
},
164+
{
165+
"cell_type": "code",
166+
"execution_count": null,
167+
"id": "4766f95c-27b5-4fd4-8d40-0e30bff59471",
168+
"metadata": {},
169+
"outputs": [],
170+
"source": [
171+
"num_obs = 500\n",
172+
"sd = 0.05\n",
173+
"measure_true = P.qoi_ref[0:num_obs]\n",
174+
"measurements = P.qoi[:, 0:num_obs]\n",
175+
"noise = np.random.randn(num_obs) * sd\n",
176+
"data = measure_true + noise\n",
177+
"# q = wme(qoi, data, sd).reshape(-1, 1)"
178+
]
179+
},
180+
{
181+
"cell_type": "code",
182+
"execution_count": null,
183+
"id": "a5147528-1873-471f-a26b-f83ac9118d6e",
184+
"metadata": {},
185+
"outputs": [],
186+
"source": [
187+
"from sklearn.preprocessing import StandardScaler"
188+
]
189+
},
190+
{
191+
"cell_type": "code",
192+
"execution_count": null,
193+
"id": "ee69d7b4-20c1-419c-80f6-2e703d7d33b2",
194+
"metadata": {},
195+
"outputs": [],
196+
"source": [
197+
"scalar = StandardScaler()"
198+
]
199+
},
200+
{
201+
"cell_type": "code",
202+
"execution_count": null,
203+
"id": "ac4493b9-5430-4725-b05d-78f7552c6df6",
204+
"metadata": {},
205+
"outputs": [],
206+
"source": [
207+
"residuals = (measurements - data) / sd\n",
208+
"X = scalar.fit_transform(residuals)"
209+
]
210+
},
211+
{
212+
"cell_type": "code",
213+
"execution_count": null,
214+
"id": "98223fa6-e8fe-4e5e-b603-b7dc8e357da1",
215+
"metadata": {},
216+
"outputs": [],
217+
"source": [
218+
"U, singular_values, singular_vectors = np.linalg.svd(X)"
219+
]
220+
},
221+
{
222+
"cell_type": "code",
223+
"execution_count": null,
224+
"id": "bb4a51d9-6d87-4dee-bd2b-1e79b64fae8a",
225+
"metadata": {},
226+
"outputs": [],
227+
"source": [
228+
"top_singular_values = 5\n",
229+
"inds = np.arange(1, top_singular_values+1)\n",
230+
"plt.plot(inds, singular_values[0:top_singular_values])\n",
231+
"plt.xticks(inds)\n",
232+
"plt.xlabel(\"Index\", fontsize=12)\n",
233+
"plt.ylabel(\"$\\sigma_i$\", fontsize=12)\n",
234+
"plt.yscale('log')\n",
235+
"plt.title(f\"Top {top_singular_values} Singular Values of Scaled Residuals\", fontsize=16)\n",
236+
"plt.show()"
237+
]
238+
},
239+
{
240+
"cell_type": "code",
241+
"execution_count": null,
242+
"id": "0a40a994-8a8e-4e3d-a457-2ecbb7cb574a",
243+
"metadata": {},
244+
"outputs": [],
245+
"source": [
246+
"scalar = StandardScaler()\n",
247+
"S = scalar.fit_transform(singular_vectors)"
248+
]
249+
},
250+
{
251+
"cell_type": "code",
252+
"execution_count": null,
253+
"id": "22bd45ca-0516-4491-b3bf-f3fae16f5707",
254+
"metadata": {},
255+
"outputs": [],
256+
"source": [
257+
"top_singular_vectors = 2\n",
258+
"fig, ax = plt.subplots(1, top_singular_vectors, figsize=(1+5*top_singular_vectors, 5), sharey=True)\n",
259+
"for i in range(top_singular_vectors):\n",
260+
" ax[i].tricontour(P.sensors[0:num_obs, 0], P.sensors[0:num_obs, 1], singular_vectors[i, :], levels=20, alpha=1)\n",
261+
" ax[i].scatter(P.sensors[0:num_obs, 0], P.sensors[0:num_obs, 1], s=100, c=singular_vectors[i, :])\n",
262+
" ax[i].set_title(f\"$\\sigma_{i+1}$ = {singular_values[i]:1.2e}\")\n",
263+
" ax[i].set_xlabel(\"$x_1$\", fontsize=12)\n",
264+
" if i == 0: ax[i].set_ylabel(\"$x_2$\", fontsize=12)\n",
265+
"fig.suptitle(\"Singular Vector Components as Weights for Sensors\", fontsize=16)\n",
266+
"fig.show()"
267+
]
268+
},
269+
{
270+
"cell_type": "code",
271+
"execution_count": null,
272+
"id": "3751ce39-6ed2-4764-ac4d-636e1ed83448",
273+
"metadata": {},
274+
"outputs": [],
275+
"source": [
276+
"fig, ax = plt.subplots(1, 1, figsize=(5,4))\n",
277+
"# plt.plot((singular_vectors[0, :] + singular_vectors[1, :]))\n",
278+
"ax.plot(singular_vectors[0, :], c='xkcd:orange', label=\"$v_1$\", lw=3)\n",
279+
"ax.plot(singular_vectors[1, :], c='xkcd:light blue', label=\"$v_2$\", lw=3, zorder=0)\n",
280+
"ax.legend(fontsize=12)\n",
281+
"ax.set_xlabel(\"Sensor Index\", fontsize=12)\n",
282+
"ax.set_ylabel(\"Singular Vector Component Value\", fontsize=12)\n",
283+
"fig.show()"
284+
]
285+
},
286+
{
287+
"cell_type": "code",
288+
"execution_count": null,
289+
"id": "82a4597d-2202-4952-871d-1f1a73386929",
290+
"metadata": {},
291+
"outputs": [],
292+
"source": [
293+
"# fig, ax = plt.subplots(figsize=(3,3))\n",
294+
"# ax.scatter(P.sensors[0:num_obs, 0], P.sensors[0:num_obs, 1], s=100, c= (singular_vectors[0, :] + singular_vectors[1, :]))\n",
295+
"# fig.show()"
296+
]
297+
},
298+
{
299+
"cell_type": "code",
300+
"execution_count": null,
301+
"id": "6e2d77d7-6512-48e6-b428-b766d16c7820",
302+
"metadata": {},
303+
"outputs": [],
304+
"source": [
305+
"num_qoi = 2\n",
306+
"new_qoi_map = singular_vectors[0:num_qoi, :]\n",
307+
"new_qoi = residuals @ new_qoi_map.T # ok shape, wrong solution, identical to using U@sigma"
308+
]
309+
},
310+
{
311+
"cell_type": "code",
312+
"execution_count": null,
313+
"id": "e85a3cfe-574c-4a9f-8998-7ac86910aeac",
314+
"metadata": {},
315+
"outputs": [],
316+
"source": [
317+
"ax.contour"
318+
]
319+
},
320+
{
321+
"cell_type": "code",
322+
"execution_count": null,
323+
"id": "0f0a725c-1a39-4c53-ac98-3385b0bb17e8",
324+
"metadata": {},
325+
"outputs": [],
326+
"source": [
327+
"top_singular_vectors = 2\n",
328+
"fig, ax = plt.subplots(1, top_singular_vectors, figsize=(1+5*top_singular_vectors, 5), sharey=True)\n",
329+
"for i in range(top_singular_vectors):\n",
330+
" ax[i].scatter(P.lam[:,0], P.lam[:,1], s=50, c=new_qoi[:,i])\n",
331+
" ax[i].tricontour(P.lam[:,0], P.lam[:,1], new_qoi[:,i], levels=20)\n",
332+
"# ax[i].set_title(f\"$\\sigma_{i}$ = {singular_values[i]:1.2e}\")\n",
333+
" ax[i].set_title(f\"$q_{i+1}$\", fontsize=12)\n",
334+
" ax[i].set_xlabel(\"$\\lambda_1$\", fontsize=12)\n",
335+
" if i == 0: ax[i].set_ylabel(\"$\\lambda_2$\", fontsize=12)\n",
336+
"# fig.supylabel(\"$\\lambda_2$\")\n",
337+
"fig.suptitle(\"QoI Component Surface Plot\", fontsize=16)\n",
338+
"fig.show()"
339+
]
340+
},
341+
{
342+
"cell_type": "code",
343+
"execution_count": null,
344+
"id": "66a5be8d-fc6f-47c9-a650-1703d2a49cf2",
345+
"metadata": {},
346+
"outputs": [],
347+
"source": [
348+
"d = DensityProblem(P.lam, new_qoi, P.domain, weights=None)\n",
349+
"print(d.mud_point())"
350+
]
351+
},
352+
{
353+
"cell_type": "code",
354+
"execution_count": null,
355+
"id": "60d1e3af-071e-4820-bff3-38a3b8421db9",
356+
"metadata": {},
357+
"outputs": [],
358+
"source": [
359+
"P.plot()\n",
360+
"plt.plot(np.linspace(0, 1, input_dim + 2), [0] + list(d.mud_point()) + [0], lw=5, c='blue', label='svd')\n",
361+
"plt.plot(np.linspace(0, 1, input_dim + 2), [0] + list(mud_problem.mud_point()) + [0], lw=5, c='purple', label='split')\n",
362+
"plt.legend()\n",
363+
"plt.show()"
364+
]
365+
}
366+
],
367+
"metadata": {
368+
"kernelspec": {
369+
"display_name": "Python (myenv)",
370+
"language": "python",
371+
"name": "myenv"
372+
},
373+
"language_info": {
374+
"codemirror_mode": {
375+
"name": "ipython",
376+
"version": 3
377+
},
378+
"file_extension": ".py",
379+
"mimetype": "text/x-python",
380+
"name": "python",
381+
"nbconvert_exporter": "python",
382+
"pygments_lexer": "ipython3",
383+
"version": "3.9.7"
384+
}
385+
},
386+
"nbformat": 4,
387+
"nbformat_minor": 5
388+
}

0 commit comments

Comments
 (0)