[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
-
Aug 21, 2019
+
Aug 27, 2019
diff --git a/doc/pub/Regression/html/._Regression-bs028.html b/doc/pub/Regression/html/._Regression-bs028.html
index d66b2671d..8f1870b01 100644
--- a/doc/pub/Regression/html/._Regression-bs028.html
+++ b/doc/pub/Regression/html/._Regression-bs028.html
@@ -398,10 +398,15 @@ The examples we have looked at so far are cases where we normally can
invert the matrix \( \boldsymbol{X}^T\boldsymbol{X} \). Using a polynomial expansion as we
did both for the masses and the fitting of the equation of state,
leads to row vectors of the design matrix which are essentially
-orthogonal due to the polynomial character of our model. This may
+orthogonal due to the polynomial character of our model. Obtaining the inverse of the design matrix is then often done via a so-called LU, QR or Cholesky decomposition.
+More material to come here.
+
+
+This may
however not the be case in general and a standard matrix inversion
-algorithm based on say LU decomposition may lead to singularities. We will see an example of this below when we try to fit
-the coupling constant of the widely used Ising model.
+algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below.
+
+
There is however a way to partially circumvent this problem and also gain some insight about the ordinary least squares approach.
@@ -415,6 +420,9 @@ inversion algorithm. Thereafter we dive into the math of the SVD.
+
[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
-
Aug 21, 2019
+
Aug 27, 2019
@@ -1300,10 +1300,15 @@ The examples we have looked at so far are cases where we normally can
invert the matrix \( \boldsymbol{X}^T\boldsymbol{X} \). Using a polynomial expansion as we
did both for the masses and the fitting of the equation of state,
leads to row vectors of the design matrix which are essentially
-orthogonal due to the polynomial character of our model. This may
+orthogonal due to the polynomial character of our model. Obtaining the inverse of the design matrix is then often done via a so-called LU, QR or Cholesky decomposition.
+More material to come here.
+
+
+This may
however not the be case in general and a standard matrix inversion
-algorithm based on say LU decomposition may lead to singularities. We will see an example of this below when we try to fit
-the coupling constant of the widely used Ising model.
+algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below.
+
+
There is however a way to partially circumvent this problem and also gain some insight about the ordinary least squares approach.
@@ -1314,6 +1319,9 @@ inversion algorithm. Thereafter we dive into the math of the SVD.
+
+
[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
-
Aug 21, 2019
+
Aug 27, 2019
@@ -1368,10 +1368,15 @@ The examples we have looked at so far are cases where we normally can
invert the matrix \( \boldsymbol{X}^T\boldsymbol{X} \). Using a polynomial expansion as we
did both for the masses and the fitting of the equation of state,
leads to row vectors of the design matrix which are essentially
-orthogonal due to the polynomial character of our model. This may
+orthogonal due to the polynomial character of our model. Obtaining the inverse of the design matrix is then often done via a so-called LU, QR or Cholesky decomposition.
+More material to come here.
+
+
+This may
however not the be case in general and a standard matrix inversion
-algorithm based on say LU decomposition may lead to singularities. We will see an example of this below when we try to fit
-the coupling constant of the widely used Ising model.
+algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below.
+
+
There is however a way to partially circumvent this problem and also gain some insight about the ordinary least squares approach.
@@ -1384,6 +1389,9 @@ inversion algorithm. Thereafter we dive into the math of the SVD.
+
[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
-
Aug 21, 2019
+
Aug 27, 2019
@@ -1373,10 +1373,15 @@ The examples we have looked at so far are cases where we normally can
invert the matrix \( \boldsymbol{X}^T\boldsymbol{X} \). Using a polynomial expansion as we
did both for the masses and the fitting of the equation of state,
leads to row vectors of the design matrix which are essentially
-orthogonal due to the polynomial character of our model. This may
+orthogonal due to the polynomial character of our model. Obtaining the inverse of the design matrix is then often done via a so-called LU, QR or Cholesky decomposition.
+More material to come here.
+
+
+This may
however not the be case in general and a standard matrix inversion
-algorithm based on say LU decomposition may lead to singularities. We will see an example of this below when we try to fit
-the coupling constant of the widely used Ising model.
+algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below.
+
+
There is however a way to partially circumvent this problem and also gain some insight about the ordinary least squares approach.
@@ -1389,6 +1394,9 @@ inversion algorithm. Thereafter we dive into the math of the SVD.
+
+
+
diff --git a/doc/pub/Regression/ipynb/Regression.ipynb b/doc/pub/Regression/ipynb/Regression.ipynb
index 3517bc4d6..d63d6e1c6 100644
--- a/doc/pub/Regression/ipynb/Regression.ipynb
+++ b/doc/pub/Regression/ipynb/Regression.ipynb
@@ -10,7 +10,7 @@
" \n",
"**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n",
"\n",
- "Date: **Aug 21, 2019**\n",
+ "Date: **Aug 27, 2019**\n",
"\n",
"Copyright 1999-2019, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n",
"\n",
@@ -379,610 +379,10 @@
{
"cell_type": "code",
"execution_count": 1,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/html": [
- "
\n",
- "\n",
- "
\n",
- " \n",
- "
\n",
- "
\n",
- "
1
\n",
- "
A
\n",
- "
A^(2/3)
\n",
- "
A^(-1/3)
\n",
- "
1/A
\n",
- "
\n",
- "
\n",
- "
A
\n",
- "
\n",
- "
\n",
- "
\n",
- "
\n",
- "
\n",
- "
\n",
- " \n",
- " \n",
- "
\n",
- "
1
\n",
- "
1.0
\n",
- "
1.0
\n",
- "
1.000000
\n",
- "
1.000000
\n",
- "
1.000000
\n",
- "
\n",
- "
\n",
- "
2
\n",
- "
1.0
\n",
- "
2.0
\n",
- "
1.587401
\n",
- "
0.793701
\n",
- "
0.500000
\n",
- "
\n",
- "
\n",
- "
3
\n",
- "
1.0
\n",
- "
3.0
\n",
- "
2.080084
\n",
- "
0.693361
\n",
- "
0.333333
\n",
- "
\n",
- "
\n",
- "
4
\n",
- "
1.0
\n",
- "
4.0
\n",
- "
2.519842
\n",
- "
0.629961
\n",
- "
0.250000
\n",
- "
\n",
- "
\n",
- "
5
\n",
- "
1.0
\n",
- "
5.0
\n",
- "
2.924018
\n",
- "
0.584804
\n",
- "
0.200000
\n",
- "
\n",
- "
\n",
- "
6
\n",
- "
1.0
\n",
- "
6.0
\n",
- "
3.301927
\n",
- "
0.550321
\n",
- "
0.166667
\n",
- "
\n",
- "
\n",
- "
7
\n",
- "
1.0
\n",
- "
7.0
\n",
- "
3.659306
\n",
- "
0.522758
\n",
- "
0.142857
\n",
- "
\n",
- "
\n",
- "
8
\n",
- "
1.0
\n",
- "
8.0
\n",
- "
4.000000
\n",
- "
0.500000
\n",
- "
0.125000
\n",
- "
\n",
- "
\n",
- "
9
\n",
- "
1.0
\n",
- "
9.0
\n",
- "
4.326749
\n",
- "
0.480750
\n",
- "
0.111111
\n",
- "
\n",
- "
\n",
- "
10
\n",
- "
1.0
\n",
- "
10.0
\n",
- "
4.641589
\n",
- "
0.464159
\n",
- "
0.100000
\n",
- "
\n",
- "
\n",
- "
11
\n",
- "
1.0
\n",
- "
11.0
\n",
- "
4.946087
\n",
- "
0.449644
\n",
- "
0.090909
\n",
- "
\n",
- "
\n",
- "
12
\n",
- "
1.0
\n",
- "
12.0
\n",
- "
5.241483
\n",
- "
0.436790
\n",
- "
0.083333
\n",
- "
\n",
- "
\n",
- "
13
\n",
- "
1.0
\n",
- "
13.0
\n",
- "
5.528775
\n",
- "
0.425290
\n",
- "
0.076923
\n",
- "
\n",
- "
\n",
- "
14
\n",
- "
1.0
\n",
- "
14.0
\n",
- "
5.808786
\n",
- "
0.414913
\n",
- "
0.071429
\n",
- "
\n",
- "
\n",
- "
15
\n",
- "
1.0
\n",
- "
15.0
\n",
- "
6.082202
\n",
- "
0.405480
\n",
- "
0.066667
\n",
- "
\n",
- "
\n",
- "
16
\n",
- "
1.0
\n",
- "
16.0
\n",
- "
6.349604
\n",
- "
0.396850
\n",
- "
0.062500
\n",
- "
\n",
- "
\n",
- "
17
\n",
- "
1.0
\n",
- "
17.0
\n",
- "
6.611489
\n",
- "
0.388911
\n",
- "
0.058824
\n",
- "
\n",
- "
\n",
- "
18
\n",
- "
1.0
\n",
- "
18.0
\n",
- "
6.868285
\n",
- "
0.381571
\n",
- "
0.055556
\n",
- "
\n",
- "
\n",
- "
19
\n",
- "
1.0
\n",
- "
19.0
\n",
- "
7.120367
\n",
- "
0.374756
\n",
- "
0.052632
\n",
- "
\n",
- "
\n",
- "
20
\n",
- "
1.0
\n",
- "
20.0
\n",
- "
7.368063
\n",
- "
0.368403
\n",
- "
0.050000
\n",
- "
\n",
- "
\n",
- "
21
\n",
- "
1.0
\n",
- "
21.0
\n",
- "
7.611663
\n",
- "
0.362460
\n",
- "
0.047619
\n",
- "
\n",
- "
\n",
- "
22
\n",
- "
1.0
\n",
- "
22.0
\n",
- "
7.851424
\n",
- "
0.356883
\n",
- "
0.045455
\n",
- "
\n",
- "
\n",
- "
23
\n",
- "
1.0
\n",
- "
23.0
\n",
- "
8.087579
\n",
- "
0.351634
\n",
- "
0.043478
\n",
- "
\n",
- "
\n",
- "
24
\n",
- "
1.0
\n",
- "
24.0
\n",
- "
8.320335
\n",
- "
0.346681
\n",
- "
0.041667
\n",
- "
\n",
- "
\n",
- "
25
\n",
- "
1.0
\n",
- "
25.0
\n",
- "
8.549880
\n",
- "
0.341995
\n",
- "
0.040000
\n",
- "
\n",
- "
\n",
- "
26
\n",
- "
1.0
\n",
- "
26.0
\n",
- "
8.776383
\n",
- "
0.337553
\n",
- "
0.038462
\n",
- "
\n",
- "
\n",
- "
27
\n",
- "
1.0
\n",
- "
27.0
\n",
- "
9.000000
\n",
- "
0.333333
\n",
- "
0.037037
\n",
- "
\n",
- "
\n",
- "
28
\n",
- "
1.0
\n",
- "
28.0
\n",
- "
9.220873
\n",
- "
0.329317
\n",
- "
0.035714
\n",
- "
\n",
- "
\n",
- "
29
\n",
- "
1.0
\n",
- "
29.0
\n",
- "
9.439131
\n",
- "
0.325487
\n",
- "
0.034483
\n",
- "
\n",
- "
\n",
- "
30
\n",
- "
1.0
\n",
- "
30.0
\n",
- "
9.654894
\n",
- "
0.321830
\n",
- "
0.033333
\n",
- "
\n",
- "
\n",
- "
...
\n",
- "
...
\n",
- "
...
\n",
- "
...
\n",
- "
...
\n",
- "
...
\n",
- "
\n",
- "
\n",
- "
238
\n",
- "
1.0
\n",
- "
238.0
\n",
- "
38.404723
\n",
- "
0.161364
\n",
- "
0.004202
\n",
- "
\n",
- "
\n",
- "
239
\n",
- "
1.0
\n",
- "
239.0
\n",
- "
38.512224
\n",
- "
0.161139
\n",
- "
0.004184
\n",
- "
\n",
- "
\n",
- "
240
\n",
- "
1.0
\n",
- "
240.0
\n",
- "
38.619575
\n",
- "
0.160915
\n",
- "
0.004167
\n",
- "
\n",
- "
\n",
- "
241
\n",
- "
1.0
\n",
- "
241.0
\n",
- "
38.726778
\n",
- "
0.160692
\n",
- "
0.004149
\n",
- "
\n",
- "
\n",
- "
242
\n",
- "
1.0
\n",
- "
242.0
\n",
- "
38.833832
\n",
- "
0.160470
\n",
- "
0.004132
\n",
- "
\n",
- "
\n",
- "
243
\n",
- "
1.0
\n",
- "
243.0
\n",
- "
38.940738
\n",
- "
0.160250
\n",
- "
0.004115
\n",
- "
\n",
- "
\n",
- "
244
\n",
- "
1.0
\n",
- "
244.0
\n",
- "
39.047499
\n",
- "
0.160031
\n",
- "
0.004098
\n",
- "
\n",
- "
\n",
- "
245
\n",
- "
1.0
\n",
- "
245.0
\n",
- "
39.154113
\n",
- "
0.159813
\n",
- "
0.004082
\n",
- "
\n",
- "
\n",
- "
246
\n",
- "
1.0
\n",
- "
246.0
\n",
- "
39.260582
\n",
- "
0.159596
\n",
- "
0.004065
\n",
- "
\n",
- "
\n",
- "
247
\n",
- "
1.0
\n",
- "
247.0
\n",
- "
39.366908
\n",
- "
0.159380
\n",
- "
0.004049
\n",
- "
\n",
- "
\n",
- "
248
\n",
- "
1.0
\n",
- "
248.0
\n",
- "
39.473090
\n",
- "
0.159166
\n",
- "
0.004032
\n",
- "
\n",
- "
\n",
- "
249
\n",
- "
1.0
\n",
- "
249.0
\n",
- "
39.579129
\n",
- "
0.158952
\n",
- "
0.004016
\n",
- "
\n",
- "
\n",
- "
250
\n",
- "
1.0
\n",
- "
250.0
\n",
- "
39.685026
\n",
- "
0.158740
\n",
- "
0.004000
\n",
- "
\n",
- "
\n",
- "
251
\n",
- "
1.0
\n",
- "
251.0
\n",
- "
39.790783
\n",
- "
0.158529
\n",
- "
0.003984
\n",
- "
\n",
- "
\n",
- "
252
\n",
- "
1.0
\n",
- "
252.0
\n",
- "
39.896399
\n",
- "
0.158319
\n",
- "
0.003968
\n",
- "
\n",
- "
\n",
- "
253
\n",
- "
1.0
\n",
- "
253.0
\n",
- "
40.001875
\n",
- "
0.158110
\n",
- "
0.003953
\n",
- "
\n",
- "
\n",
- "
254
\n",
- "
1.0
\n",
- "
254.0
\n",
- "
40.107212
\n",
- "
0.157902
\n",
- "
0.003937
\n",
- "
\n",
- "
\n",
- "
255
\n",
- "
1.0
\n",
- "
255.0
\n",
- "
40.212412
\n",
- "
0.157696
\n",
- "
0.003922
\n",
- "
\n",
- "
\n",
- "
256
\n",
- "
1.0
\n",
- "
256.0
\n",
- "
40.317474
\n",
- "
0.157490
\n",
- "
0.003906
\n",
- "
\n",
- "
\n",
- "
257
\n",
- "
1.0
\n",
- "
257.0
\n",
- "
40.422399
\n",
- "
0.157286
\n",
- "
0.003891
\n",
- "
\n",
- "
\n",
- "
258
\n",
- "
1.0
\n",
- "
258.0
\n",
- "
40.527188
\n",
- "
0.157082
\n",
- "
0.003876
\n",
- "
\n",
- "
\n",
- "
259
\n",
- "
1.0
\n",
- "
259.0
\n",
- "
40.631842
\n",
- "
0.156880
\n",
- "
0.003861
\n",
- "
\n",
- "
\n",
- "
260
\n",
- "
1.0
\n",
- "
260.0
\n",
- "
40.736361
\n",
- "
0.156678
\n",
- "
0.003846
\n",
- "
\n",
- "
\n",
- "
261
\n",
- "
1.0
\n",
- "
261.0
\n",
- "
40.840746
\n",
- "
0.156478
\n",
- "
0.003831
\n",
- "
\n",
- "
\n",
- "
262
\n",
- "
1.0
\n",
- "
262.0
\n",
- "
40.944999
\n",
- "
0.156279
\n",
- "
0.003817
\n",
- "
\n",
- "
\n",
- "
264
\n",
- "
1.0
\n",
- "
264.0
\n",
- "
41.153106
\n",
- "
0.155883
\n",
- "
0.003788
\n",
- "
\n",
- "
\n",
- "
265
\n",
- "
1.0
\n",
- "
265.0
\n",
- "
41.256962
\n",
- "
0.155687
\n",
- "
0.003774
\n",
- "
\n",
- "
\n",
- "
266
\n",
- "
1.0
\n",
- "
266.0
\n",
- "
41.360688
\n",
- "
0.155491
\n",
- "
0.003759
\n",
- "
\n",
- "
\n",
- "
269
\n",
- "
1.0
\n",
- "
269.0
\n",
- "
41.671089
\n",
- "
0.154911
\n",
- "
0.003717
\n",
- "
\n",
- "
\n",
- "
270
\n",
- "
1.0
\n",
- "
270.0
\n",
- "
41.774300
\n",
- "
0.154720
\n",
- "
0.003704
\n",
- "
\n",
- " \n",
- "
\n",
- "
267 rows × 5 columns
\n",
- "
"
- ],
- "text/plain": [
- " 1 A A^(2/3) A^(-1/3) 1/A\n",
- "A \n",
- "1 1.0 1.0 1.000000 1.000000 1.000000\n",
- "2 1.0 2.0 1.587401 0.793701 0.500000\n",
- "3 1.0 3.0 2.080084 0.693361 0.333333\n",
- "4 1.0 4.0 2.519842 0.629961 0.250000\n",
- "5 1.0 5.0 2.924018 0.584804 0.200000\n",
- "6 1.0 6.0 3.301927 0.550321 0.166667\n",
- "7 1.0 7.0 3.659306 0.522758 0.142857\n",
- "8 1.0 8.0 4.000000 0.500000 0.125000\n",
- "9 1.0 9.0 4.326749 0.480750 0.111111\n",
- "10 1.0 10.0 4.641589 0.464159 0.100000\n",
- "11 1.0 11.0 4.946087 0.449644 0.090909\n",
- "12 1.0 12.0 5.241483 0.436790 0.083333\n",
- "13 1.0 13.0 5.528775 0.425290 0.076923\n",
- "14 1.0 14.0 5.808786 0.414913 0.071429\n",
- "15 1.0 15.0 6.082202 0.405480 0.066667\n",
- "16 1.0 16.0 6.349604 0.396850 0.062500\n",
- "17 1.0 17.0 6.611489 0.388911 0.058824\n",
- "18 1.0 18.0 6.868285 0.381571 0.055556\n",
- "19 1.0 19.0 7.120367 0.374756 0.052632\n",
- "20 1.0 20.0 7.368063 0.368403 0.050000\n",
- "21 1.0 21.0 7.611663 0.362460 0.047619\n",
- "22 1.0 22.0 7.851424 0.356883 0.045455\n",
- "23 1.0 23.0 8.087579 0.351634 0.043478\n",
- "24 1.0 24.0 8.320335 0.346681 0.041667\n",
- "25 1.0 25.0 8.549880 0.341995 0.040000\n",
- "26 1.0 26.0 8.776383 0.337553 0.038462\n",
- "27 1.0 27.0 9.000000 0.333333 0.037037\n",
- "28 1.0 28.0 9.220873 0.329317 0.035714\n",
- "29 1.0 29.0 9.439131 0.325487 0.034483\n",
- "30 1.0 30.0 9.654894 0.321830 0.033333\n",
- ".. ... ... ... ... ...\n",
- "238 1.0 238.0 38.404723 0.161364 0.004202\n",
- "239 1.0 239.0 38.512224 0.161139 0.004184\n",
- "240 1.0 240.0 38.619575 0.160915 0.004167\n",
- "241 1.0 241.0 38.726778 0.160692 0.004149\n",
- "242 1.0 242.0 38.833832 0.160470 0.004132\n",
- "243 1.0 243.0 38.940738 0.160250 0.004115\n",
- "244 1.0 244.0 39.047499 0.160031 0.004098\n",
- "245 1.0 245.0 39.154113 0.159813 0.004082\n",
- "246 1.0 246.0 39.260582 0.159596 0.004065\n",
- "247 1.0 247.0 39.366908 0.159380 0.004049\n",
- "248 1.0 248.0 39.473090 0.159166 0.004032\n",
- "249 1.0 249.0 39.579129 0.158952 0.004016\n",
- "250 1.0 250.0 39.685026 0.158740 0.004000\n",
- "251 1.0 251.0 39.790783 0.158529 0.003984\n",
- "252 1.0 252.0 39.896399 0.158319 0.003968\n",
- "253 1.0 253.0 40.001875 0.158110 0.003953\n",
- "254 1.0 254.0 40.107212 0.157902 0.003937\n",
- "255 1.0 255.0 40.212412 0.157696 0.003922\n",
- "256 1.0 256.0 40.317474 0.157490 0.003906\n",
- "257 1.0 257.0 40.422399 0.157286 0.003891\n",
- "258 1.0 258.0 40.527188 0.157082 0.003876\n",
- "259 1.0 259.0 40.631842 0.156880 0.003861\n",
- "260 1.0 260.0 40.736361 0.156678 0.003846\n",
- "261 1.0 261.0 40.840746 0.156478 0.003831\n",
- "262 1.0 262.0 40.944999 0.156279 0.003817\n",
- "264 1.0 264.0 41.153106 0.155883 0.003788\n",
- "265 1.0 265.0 41.256962 0.155687 0.003774\n",
- "266 1.0 266.0 41.360688 0.155491 0.003759\n",
- "269 1.0 269.0 41.671089 0.154911 0.003717\n",
- "270 1.0 270.0 41.774300 0.154720 0.003704\n",
- "\n",
- "[267 rows x 5 columns]"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
@@ -1486,7 +886,9 @@
{
"cell_type": "code",
"execution_count": 2,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"# matrix inversion to find beta\n",
@@ -1505,7 +907,9 @@
{
"cell_type": "code",
"execution_count": 3,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"fit = np.linalg.lstsq(X, Energies, rcond =None)[0]\n",
@@ -1522,19 +926,10 @@
{
"cell_type": "code",
"execution_count": 4,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEMCAYAAAAs8rYIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xl0XOd93//3d3bsALFwAUiCFElZG03ZlCxLsqXYqSXbkhXLjY+cxLFyqqqpk1r1qU9/TlvHSdpfj9PTNvGWOPpVqRtboVq78ZLWchzZkmvZtWTS1kpqpUAS4AaCAEEAs977/f0xAwgAQRLEzjuf1zk4mJk7c+/zYMjPPPO99z7X3B0REakuseVugIiILD2Fv4hIFVL4i4hUIYW/iEgVUviLiFQhhb+ISBVS+IuIVCGFv4hIFVL4i4hUocRyN+Bs2travLu7e7mbISJyUdmzZ88Jd28/3/NWbPh3d3eze/fu5W6GiMhFxcwOzOZ5KvuIiFQhhb+ISBVS+IuIVCGFv4hIFVL4i4hUIYW/iEgVUviLiFShFXucv5zbcK7IoZNjtNalOTAwSkMmyXOHT1GXSpAtBhRLAVtX1TEyUqI2neDJ/QNc2dVEY22SxkySTe11EINYKoaZLXd3RGSJKfwvAmEp5NiBUfJHCzzyeB+5E0VyQ0UYCUlmnWTOSeYgVoSxkmMBxEvwwqR1NAO9HJ+4/1ODRMyoSSXIE9LakiZdkyBdF6doTrouQW1DEs8Yja1p4vVx4nVxyBjJxkT5fn2ceEOcRFP5fhA6iXj5y2QYOvuODlMKnMvWNpJK6EumyEqi8F9BwlJIoa9AridHridH9rUs2QM5Bg6NceDEGACZyg9AOhmjFDi1qTjF0KlNxvGEY0nDgJyHxFMxcqWQxnSC09kisZiRKwaUAicMnVK2CMCR7NhZ25VKxAhCJ52IMVYISCVijH9ZiJtRDJx40jgVD0m2JBhNOdmMk8tAsdb4Wp2RakmyYXMDdR1pLt/azCUd9WSS8cX7Y4rIOSn8l1FYChnbO8boc6OMPj/K2ItjeN45MZJnOFskGY8xOFYgBPINRrEpRu26NFdc3kK6OUFnVz2JxgTxxjiJhgSxmhiWNGLp8u+zlXPcnWLgHDw5yuHBLKnQ2P3ySWIlZ2AwR10swbGBLOFYSE3JCEcDknlI5J1EPl75DYmck8o5yayTyEEaYLRIA9AAJBMx4jHIFUpAERgjD+xOHOQn9UayNUmmPUVsVYLOTfXUdqS4ZGsTmY4UyZYkFlc5SmSxmLsvdxtmtHPnTo/S3D4j+RLJuBHPOcNPDvPcw8c4/tMh1qTThO7EzBgYLTDa4BxrDBltjTHSZoy0GcXmGO/avpY739S5ZPX5QikkCB0zOHIqR106zuGhHBtX1XLkVI5UotyOXDGkqSbJicEc3ekMfYdGaCrGSeeA4YDgVEBhqMjg4RyDfVkKA0WGh/IUSiEz/dOLxQxwalIJaI6TXp2ivbuOdEeKjk110BqnZX0tyfYksaRKSSLTmdked9953ucp/Bff4cEx/uzPn2fNcwHrDjhjo6WJZSNtxuD6GKfWxRjqjFGsLYfqDVvaOJ0rcdsb17K5rS5SO2WDsYDc8Twvv3SKXH+Bwb4cg4ezZPsLFPoLpEecZBbsLP8008kYBqTbUqQ6UqRXp6hdkybWlqB7WyPx9iRN62uIpfThINVntuGvss8i+MXBQf7u+WP8gy3t8PgoL361jzcMlAN/1GBofYyhrQmarmvk6bERVjdlMHeuaqnlxi1tJOMxLl/XuMy9WDzx2jh13bXs6K6d8ri7cypbxDBe6BvGTxZ5ed8pRo/kKBwrUugvUDsMqVMh6dNO7kgOjuSmrGN8J3cqESNfb8RXJwlb48Q7kqzqrsPaE1y+fRXp5gRt9elIfaiKXAiN/BfQcK7Ii0dP8+Xv76fj50XW/7xEKltZ2J6g8ZeaGNqe4o53bqQuFSceM44N52mrT00cJSNnN76v4sDAKITw4ouDBMeLnOzNEvaXYKBE9miezLCTPuVn/eYAUEpD0BrHWxPUdKZJr02RWpNm46WNJNqSbGyvIxEz7ZSWi47KPkssDJ3f/9az+E/G2PyTEulRJxE3SuuTtHygjbf9w/U01qSWu5mR5+4EodM3kCU+FHDk1RG8v8ipA1myR/KUjhXJHc4TKzilYOZ/+2Eccg1GttlIdCQptcap6cqwalMtmTUpOttqScRirG3OkIgZTTVJfYOQFUNlnyX28939rP2LYZqOOm31abp/qYXOj6yh/o31CoYlZGYk4sbGjjrogK5tZ5bP3J3SqRJHXh0h25vnyKunKR0rMnY4T/FoATsVkBgOqDvleE++8qoxcpwka3CkyRhrNsZayh8QsTVJ6tdnWLepgTdvWsX6VTXUpvRfS1Y2/Qudp2IpYNdn9uLfHKKpBOs21nP1x7tpvqlZob9CmRnJ5iQb3twCb4ZLWTNleZgPKRwtkD2c5+grI5SOFjj52hjFI3mKJ0oU8iHhESd/KCB0CMMSkKUUH+SxpkNkW2LUr8+QXpdi0+XNrN5SR/eWRpIJlZBk5VD4z8PR42Ps+0wPsR8MATD4xiTv/+wV1DWrvHMxi6VjZDZmyGzM0PLWpinLwkL5gyHfl6dwuPx74LVRsr15Bo9kGRsJyJ4M4NVRYJQeBukBnkwZqXVpajrTtF1SR+OGGja8oZF0V5pEU0IDBVlyCv852rP7OLv/5UvUnXRKaWj87TXc8aFu6jLJ5W6aLKJYKkZmQ4bMhszEY12V30EuoHC4wFDPGEdfHmGoZ4wT+0cpHi0QngrI9+TI9+QY+vEpAJ6OGbGYUduQIL4mRe36DHXrMzR117B2WwO1GzLEa/RtQRaHwn8Ocgdy9Pzea9SddEZbjcN31fP7H9lCPKbRWzWLZ+LUbK6hZnMNa9/ROvG4u9N3eJQjr4xwcN8pRg7mGO3NwfEitYMhxcECDBY4uW9kyvpiMSPdmiS+Nln+1rA+Tfsl9azeWk/LhloSOslN5kHhf4HGXhrj5X+zn2x/gcH1MTb+m43cfmmrgl/Oyszo6qynq7Oea24q718YP6dhOFvk1deGGevLM7h/jOBIgVxfHo4VqRl0wv4C9BfgmVEA+sZXmjDiHUlia5I0bMyQ7kyzZlsDnW9opLEtRSymDwY5tyUNfzP7OHAP4MCzwG+5e+7cr1o5codyvPKp1zh4aIQTm2PE/nE773tL1/lfKDKNmdFcm6K5NsWG1jqYdmCeu3NypEDPS6cY7skyeihHrjfPyKEcdqxIfCgkOFwgOFzg5M/LHwxHgF8AYY1BR4JUZ5r6DRlq12dYf1kTnZc2UFunsqSULVn4m1kn8DHgcnfPmtn/AO4CvrxUbZiP4wdGOfSvXqPv4Aj71wTsuz3Fx69ec/4XisyBmdHakKb1zR3w5jOXnxrKc+Tl05zYP8rQ/jEKRwqMHsjhx4rEsg4HipQOFBn6yQhDwGF6AQhb4sTWJEl1pqjfUMOqS+rY8IZG2jfW6mikKrPUZZ8EUGNmRaAWOLzE25+1UhASj5VnxgwLIQ9/9FnifUVOrY3x0u1p/vDOK1nTlDn/ikQWQVNzmqZr0nBN25THwzBkuL/A8ZdOc/yVEYYPZMn25hk7lIMTJWKDAQwGFPblOMkwJ4FXgDABQVuifM7CxgxNG2tZtbmWddsaaOnIkNYHQ+QsWfi7e5+Z/UfgIJAFvufu35v8HDO7F7gXYMOGDUvVtDPsPTzMf/rei9x9Qzdv29rOK587RLyvSLbZeOZ9SX7pyg4Fv6xIsViM5tUZmldn2Pa29inLgkLI8Z4R+l8Z5eRrY4wcyDJyMEfxcIHYcEDsaAmOlsg+lSXLIEeBvUCxFrwjSdAep2FjDanONKs217FqYy0tjSla61Mk4zGSmqLkorKUZZ8W4A5gEzAEfM3MfsPdvzr+HHe/H7gfytM7LFXbpvvFoUEAfrp/gCt7Exz/7gBhHJ57b5K2NTXceqXKPXLxiadirN3WyNoZznrODxc5sX+UE6+OcvTlEXKHcuT68oRHi5ANoacIPUVKP8tRAsaAQwa5ytnO2VVGqjMNHUnqNqRJrkrSUp+iIZ2kIZOgvSFNTSpOU02SmmRc5zWsAEtZ9vll4DV37wcws78Brge+es5XLYODA+WrWh1+9TS9PygxVijx8s0J3v6OTj7wZu3glehJNybp3NFM545m3jjpcXdn7FiB46+cZuRgjhOvjFI8kid7KI8PlEiPhNQNO2FPCD9//WpwQRKONhv7VxnZ5hhjLa9PiWG1MeorHwqb2up404YWtq7Wld2W2lKG/0HgOjOrpVz2eSew4mZuc3cODY6BO93fKzKWTXB8a5zDV8KvtNUtd/NElpSZUbcmzaY16TOWTZzt3Jtn9GCOk/tHyfXlyR7KEYwEFPIhQZ9TPOgUg/LFgUqhk8s4Yy15xlqM11qG2dt8lOwqo3VTHdvWNbJtdT1bOhpoqtGRSYtpKWv+T5jZ14GfAyXKR6Xdv1Tbn63jp/PkiyHtr4S09oT0Nud58toEmLFJ4S8yYfLZzk3XN7Fu0rLS6VJ5Coy+8hQYE9NhHM5TypU/CIr9IcMHi4zkSowVA5wCI01D/KjF+F5zuYy0ems93Zc3sXVbEx2NGZWLFtCSHu3j7p8GPr2U27xQB0+OEc87W39YvvjKU1c7o2ln/apaWmo1EhGZjURDgsQbEtS9YeqAyd0p9henfiBUjkY61ZdjNF9i9EiJ0Z4S/osxYIwejvNqEoL2OLXra6jtSrPu0gbWbWugY2s9iTqdqzoX+qsBj75wnNFCidu2r+PAwBgb9pTotBT5yxJc+aFWPriumUvXNGjUITJPZkaqo3z5zYarG6YsCwshhSPlbwrZ3hzHXhqhf/8oIweyZIeKBIcDSodHGH5ihGEGeAEwg3hTgnRXmoYNGRq7y4entl1SR2ptilhCRyCdTdWHv7vz1Z8eAODa7lXse3mQ9U8FNDRleNvvX0rdZSr1iCyFWOr12VSbaGINqyeWFYeLHH7pNH37TnOyZ5ST+8coHSkQPxHgQyVKQyVGnxvlKAO8RPlDIZmMYW0JaroyNGysoXFjDa2ba2nfUk+mPVX1g7mqD/98KZy4/cRrJ7EfnCZZhNVvaVbwi6wQycYkG3euYuPOVVMeH80XOfzqCIdeHK6cu5Aj25sj0R+QOR1ihwvkDxcYenL49RcZJGvi0FH+xlD+UKhjzdZ6WjfXkaivjlisjl6ew0i+NHH7f/2kl7c+FdCQSbL211ef41UishLUpZNsvbyFrZe3THk8VwwYGMxxbP8oR148zfDBLPm+PMGRArHjJRgLoCeg2JNn5PFhDlOebMwMrCmOrU6S6UrTuKGGls11rN5Sx+pLGkikolNGqvrwH8sHE7fX/yIgUYDGG+o16he5iGWScTo76ujsqIPrOqYsyxUDjvaNMXhglP5XRzn12hijh3IER4okBgLiQwEMBYy+mGOUUxyhfKYzMfDWBLYmSe36DE0ba2i7pI512xpo7aq56GZSrfrwHx/5x0pO5zMlUokY2+9dvqklRGRxZZJxursb6O5ugJumLhvLl+g/OMrA/lEGXh1j+ECWXG+O0pEiscEA6y9Bf4mxZ7OMMcgRyt8YPF2eSTW9Lk3Txhqau2tYvaV8RFJyhV6prerDf7RQDv+Ol0I64kk2XbuKlu0N53mViERRbTrBxq1NbNzadMaysdEiA6+N0f/qKAP7Rzhd2b8QHCkSGw3hUJHCoSL9T4zQD7xM+YI8maYEmy5rpnlTLbVdadKdaVLrUqTXpYnXLt9ZzQr/ysj/LQdTbFiVZO0d7SvyU1pElldtXZLaK5tYf+WZHwxDJ3KceHWUwy+d5uRro4z15sn15Yn1lwgHi+z76Qn8/zrpRIy1TTU0ZMrfBhLNiYkPglRn5XflfmyR9y8o/PMBtQMhmcNObF2apred+caKiJxLc1uG5rYMW97SOuXxodECf/W9/fTsPUXtkFMzGHJoKEvLSIw1+TixU2AHjfp0grrUpAnvDK74+hWL+gGg8M+XWLs3IG5Jmt/WvOiftiJSPZrrUvyT27fy+BtOsLG1jleOj/DoC8d5biQP7qRPQ+1QQM1giebRGN1BitpBx925SiP/xTWaK7L6hZB42mh+Z/NyN0dEIiadiPPOy8qHjm/pqOeWK1bz/OFhDlRmDz6dK/J07yn2DufYS7kMnYwZ71nkdlV9+Acv58mMOsn1SWovrV3u5ohIxJkZV3Y2cWXn6yXmu66Fo6dyPNM7RCl0alNx3H1R9z9WffjHnsoCUPOWRu3oFZFls6Ypw5qmpbtQVOQL3I/sPcajLxyfuJ8rBjyxf4BcMcDdqXkuD0DTDWde3UhEJKoiPfIPQmfXkwcBuPnS8iGcDz15kB+9fII3d7dwd9c67FRIrsFovrx+mVsrIrJ0Ij3yD8LXLwNcqtzefaB8fd49PYOM7BkhCEMGumPUpSP9OSgiMkWkwz/018O/GJRn75x8ndChnw3jDkOb4qQ177eIVJFIJ96k7KdQGg//cpcTOef0vlHCOBS3am5vEakukQ7/qSP/8u2aysi/5VBIsRgyvCZGc8uZF6cWEYmyKgr/qWWfpsMhhVLIYFeM1nqFv4hUl2iH/+sX6ZoI/0Rlzu3m3pBcKWCoK0ZrXWo5micismyiHf4zjPwLQUAi79SfcMaCkOE1Rmu9wl9EqkvVhH+h5JXfIU19IeZwbFVImDSVfUSk6kQ8/F+/PT7yLwZO47HygoHK1d1U9hGRahPp8PfJI/+Jsk9Iw7Hy7eHV5e63aeQvIlUm0uE/08i/UAxoPFq+fXpNjJpUnJrU8l1KTURkOUQ6/IMpNf9y4NvJgGQOijWQa9CoX0SqU6TDP5w8t0/lJK9MX/liCcOrY2DG7W9ctyxtExFZTpGezWzK9A5BWJ7C+WgAwI63dfChO9fT0ZhZptaJiCyfSId/MO04/1Lo1J5wzGDHdR00KfhFpEpFu+wzbW6fYhBSP1C+NFp6o2r9IlK9Ih3+Pm2Hb264SGbYsYSRXqvwF5HqFenwDybN7VMKQ0YP5ADIt8WwuKZwFpHqFemafzht5J89XA7/wmod1y8i1W1JR/5m1mxmXzezF8xsn5m9dTG3F047wzd3sHyx9mBNpD/zRETOa6lT8LPAd939H5pZCqhdzI1NPtSzWHIKfZXw71D4i0h1W7IUNLMm4O3A3QDuXgAKi7nN6VM6F4+WN+ftCn8RqW5LWfbZBPQD/9XMfmFm/8XM6iY/wczuNbPdZra7v79/3hucMrdPISA4Xixvp03hLyLVbSnDPwG8Cfhzd78aGAU+OfkJ7n6/u+90953t7e3z3mAwKf19MCAMnFyDkazRDl8RqW5LGf69QK+7P1G5/3XKHwaLZvJx/vETJUJ3sk1GKhHpI1xFRM5ryVLQ3Y8Ch8zs0spD7wT2LuY2J5d94icCwhByTUYyrvAXkeq21MXvfwY8WDnSZz/wW4u5sck7fBMDAe4xxpqNlMJfRKrckoa/uz8F7Fyq7U0J/5MhITGyzSr7iIhEOgXDSdM7pIZCwtDJNarsIyIS6RScPPJPng4JvXy0j0b+IlLtIp2C4+EfKzrJHARWvnxjUpO6iUiVi3T4jw/80yPlG8OZEMxYVZdaxlaJiCy/SIf/+Ele6ZHy/ZOp8k6ADasWdUohEZEVL9LhP172yZwu/87XG401SZprNfIXkeoW8fAv/x4v++QaNOoXEYHIh3859LsT5Qu1F+qN9Qp/EZHqCP/WQnkit1yD0Vavko+ISKTnNh4v+8RPBWxorSV3aS3XbW5d3kaJiKwAkQ7/8Vk9Y6dCWmpTfOQ9W0gmNZ2ziMh5yz5m9kUzu2EpGrPQQndwx8bKh3jGGxX8IiIwu5r/S8B/NLMeM/sPZnb1YjdqoQQhJHNgDvH6ODFN6yAiAswi/N39s+7+VuAmYAD4SzN7wcw+bWbbFr2F8xC6kxxzDIg3adQvIjJu1kNhdz/g7n9cuQTjh4BfAfYtWssWgLuTygIGiaZI794QEbkgsw5/M0uY2e1m9iDwMPAicOeitWwBhA7JbHnkr/AXEXndeRPRzP4B5ZH+e4AngYeAe919dJHbNm/jZR8whb+IyCSzScTfA/4a+BfuPrjI7VlQoUMqC2aq+YuITHbe8Hf3dwBY2W8Am939j8xsA7DG3Z9c7EbOVRi6yj4iIjO4kGMf/wx4K+USEMBp4IsL3qIFFLqTGnPt8BURmeZCEvEt7v4mM/sFgLsPmtmKniinvMMXTDV/EZEpLmTkXzSzOOAAZtYOhOd+yfLy8ZE/qvmLiEx2IeH/OeAbQIeZ/b/A48C/X5RWLZBANX8RkRnNOhHd/UEz2wO8EzDgV9x9RZ/kFYblC7dTC/EGjfxFRMZd0HDY3V8AXliktiw4z4aYQywT07w+IiKTzOYkr2+fa7m7v2/hmrPAKrN5UqvgFxGZbDYj/7cCh4BdwBOUSz4XBctW5vNX+IuITDGb8F8DjE/x8GvA/wZ2ufvzi9mwBZEtj/ytRvV+EZHJZjOlc+Du33X3jwDXAa8Aj5nZ7y566+bJcuXwj9Vp5C8iMtmsdviaWRp4L+XRfzevH/a5ollOZR8RkZnMZofvXwFXAt8B/tDdn1v0Vi2U8ZF/rco+IiKTzWbk/xvAKHAfcJ+ZeeVxA9zdGxercfM1scO3TuEvIjLZbML/RuD/uruf95krTKxS9omr5i8iMsVsUvHDwB4ze8jM7jazNYvdqIUyscNXNX8RkSlmM5//PwUwszcA7wa+bGZNwKPAd4Efu3sw2w1WJofbDfS5+21zavVst1Up+8RV8xcRmeJCLuD+grv/ibvfCryD8sRuv0r5xK8LcR9LdOH3WF7hLyIykznVQ9w96+7fcfd/5u47Z/s6M+uifMjof5nLdi+UTdT8Ff4iIpMtdTH8T4F/yRJdByBWqfnH6xX+IiKTXfAk92b2NaAOSAEB5cM9b53F624Djrv7HjO7+SzPuRe4F2DDhg0X2rQzxHPlT5l4vXb4iohMdsGp6O6/SnmH7S3ArcDfz/KlNwDvM7Me4CHgHWb21Wnrvt/dd7r7zvb29gtt2hli+fLIP6Gav4jIFHMdEm8DOilP+rZ5Ni9w999z9y537wbuAn7g7r8xx+2ff3uBYwVw06GeIiLTzfXahp8GPkb5er6fX7jmLJxgLACHIAWJmMJfRGSyOYW/u78IfGKuG3X3x4DH5vr62QhGAxynlDJidtFcgkBEZEnMKfzN7F8BV1AuG7m7/9qCtmoBhJUjfYIkCn8RkWnmWvYxd//1BW3JAvOC414Of1PVR0RkirmG/yVm9quUZ/vE3b+zcE1aGGE+xIEwobKPiMh0cw3/HwK1lZ8VOdunF8rNChIQU/aLiEwxl5O8bgD6ef1C7isy/MN8iLsTJjXyFxGZbi4j/1VAC6+fI7AkUzVcqLByglegso+IyBnmcobv31I+ueu9wHuA1QvdqIVQHvlDqLKPiMgZ5nocTLu7/6q7f5AVGv5emc45TEBc6S8iMsVcav7vBdaa2YcrD63IK3sFufL1ZYIEmMo+IiJTnHfkb2ZXTHuojfJkbrHKz2wndltSQaXmHyYV/CIi081m5P8V4E0AZnaPu09ciMXMat19bLEaNx9B5QxfUgp/EZHpZlPzn5yeH5227EcL2JYFNT7yd438RUTOMJvwn3wc//QkXbETJ4QKfxGRs5pN2WeNmd0NPM2Z4b8iT/ACCCpH+3hymRsiIrICzSb8/wB4M/BbQJeZ7QX2AS9Q3vm7IoX58tE+GvmLiJzpvOHv7vdPvm9mXcBVwHbg/yxSu+YtzFW+lCj8RUTOcMHH+bt7L9ALPLzwzVk4Ezt8dbSPiMgZVuwO2/kKCzrUU0TkbKIb/jraR0TkrCIb/l4Jf0sr/EVEpots+I8f6mnpyHZRRGTOIpuM4xdwjyv8RUTOEMlkdPeJHb5Jhb+IyBkimYxecEJ3ggQkEpHsoojIvEQyGSdfxSsRi2QXRUTmJZLJGBZCQpwgaaQ08hcROUMkk9HzPmnkr0M9RUSmi2T4h/lwouafjEeyiyIi8xLJZHy95m8kVfYRETlDJJNxfOQfJiCpso+IyBkiGf5eqtT845BQ2UdE5AzRTMYA3MFjkIxr5C8iMl0kw99L5ZO8yuEfyS6KiMxLJJPRg9fLPgp/EZEzLVkymtl6M3vUzPaa2fNmdt9ibctLjuN4zFT2ERGZwQVfxnEeSsC/cPefm1kDsMfM/t7d9y70hjxQ2UdE5FyWLBnd/Yi7/7xy+zSwD+hclG2VnNAhjOkMXxGRmSzLsNjMuoGrgScWY/3jh3p6DJ3kJSIygyVPRjOrB/4n8M/dfXjasnvNbLeZ7e7v75/zNso7fCtlH83qKSJyhiVNRjNLUg7+B939b6Yvd/f73X2nu+9sb2+f83bGyz7lkb/KPiIi0y3l0T4GPADsc/f/vKgbCypX89KhniIiM1rKZLwB+DDwDjN7qvLznsXYkAdOCLiZyj4iIjNYskM93f1xYElqMBM1/7iR0HH+IiJniOSweErNX2UfEZEzRDIZx0f+oSZ2ExGZUTTDXyN/EZFzimYyVo72cR3tIyIyo0gmY1gMK2f4amI3EZGZRDP8Sw6AxaF8eoGIiEwWyfAvFUIAYslIdk9EZN4imY5BUB75xzW1g4jIjCIZ/hr5i4icWyTTMSiWwz+u6ZxFRGYUyXQc3+GbSKrsIyIyk0iG//jIX2UfEZGZRTIdg8rIP66Rv4jIjCIZ/mFl5J9QzV9EZEaRTMfxkX8yFcnuiYjMWyTTcfxQz1QqvswtERFZmSIZ/uMj/1Qmkt0TEZm3SKbj+NE+aY38RURmtGSXcVxK4+GfSkfys03kolUsFunt7SWXyy13Uy56mUyGrq4uksnknF4fyfAfP8krlY5k90QuWr29vTQ0NNDd3a0Zd+fB3RkYGKC3t5dNmzbNaR2RHBqPh39GNX+RFSWXy9Ha2qrgnyczo7W1dV7foCKZjmFp/GgfjfxFVhoF/8KY798xmuFf1MhfRM7um9/8JmbGCy+8sCDre/DBB9m+fTtXXXUV119/PU8//fTEsu9+97tceumlbNmyhc985jMTj3977GFPAAAMUUlEQVThC19gy5YtmBknTpyYsr7HHnuMHTt2cMUVV3DTTTctSBuni2Q6hpX5/NOq+YvIDHbt2sWNN97Irl27FmR9mzZt4oc//CHPPvssn/rUp7j33nsBCIKA3/md3+Hhhx9m79697Nq1i7179wJwww038Mgjj7Bx48Yp6xoaGuKjH/0o3/72t3n++ef52te+tiBtnC5y4e/ueKXmn9bIX0SmGRkZ4fHHH+eBBx7goYceAsoj7Ztuuok77riDzZs388lPfpIHH3yQa6+9lquuuopXX30VgP7+fj7wgQ9wzTXXcM011/DjH/8YgOuvv56WlhYArrvuOnp7ewF48skn2bJlC5s3byaVSnHXXXfxrW99C4Crr76a7u7uM9r313/919x5551s2LABgI6OjkX5O0RuaOyBE7jjBjWq+YusWP/oyz9blPU+cPc151z+rW99i1tvvZVt27bR2trKnj17AHj66afZt28fq1atYvPmzdxzzz08+eSTfPazn+Xzn/88f/qnf8p9993Hxz/+cW688UYOHjzILbfcwr59+6Zu/4EHePe73w1AX18f69evn1jW1dXFE088cc72vfTSSxSLRW6++WZOnz7Nfffdx2/+5m/O5U9xTtFLxwDCEMI4pDWxm4hMs2vXLu677z4A7rrrLnbt2sVtt93GNddcw9q1awG45JJLeNe73gXAVVddxaOPPgrAI488MlG2ARgeHmZkZIT6+noAHn30UR544AEef/zxObevVCqxZ88evv/975PNZnnrW9/Kddddx7Zt2+a8zplELvy95ITueBwySZ3hK7JSnW+EvhhOnjzJD37wA5599lnMjCAIMDPe+973kk6nJ54Xi8Um7sdiMUqlEgBhGPLTn/6UTCZzxrqfeeYZ7rnnHh5++GFaW1sB6Ozs5NChQxPP6e3tpbOz85xt7OrqorW1lbq6Ourq6nj729/O008/veDhH7mhsQeV8I9BRiN/EZnk61//Oh/+8Ic5cOAAPT09HDp0iE2bNvGjH/1oVq9/17vexec///mJ+0899RQABw8e5M477+QrX/nKlJC+5pprePnll3nttdcoFAo89NBDvO997zvnNu644w4ef/xxSqUSY2NjPPHEE1x22WVz6O25RS4di/kAdyBuJOKR656IzMOuXbt4//vfP+WxD3zgA7M+6udzn/scu3fvZvv27Vx++eV86UtfAuCP/uiPGBgY4KMf/Sg7duxg586dACQSCb7whS9wyy23cNlll/HBD36QK664YmJdXV1d9Pb2sn37du655x4ALrvsMm699Va2b9/Otddeyz333MOVV165UH+CCebuC77ShbBz507fvXv3Bb9usC/L/37vbkpNMe7+4Q2L0DIRmat9+/Ytyii2Ws309zSzPe6+83yvjdzQOJ8r1+ZiCZ1FKCJyNtEL/3zl4u0KfxGRs4pe+I+P/FXvFxE5qyVNSDO71cxeNLNXzOyTi7GNfOUSjvGkRv4iImezZOFvZnHgi8C7gcuBD5nZ5Qu9nXwuACCm8BcROaulHPlfC7zi7vvdvQA8BNyx0BvJ58vhH9cx/iIiZ7WUCdkJHJp0v7fy2IIqjod/UuEvImeKx+Ps2LFj4qenp4fdu3fzsY99DChP8vaTn/xkmVu5+FbU9A5mdi9wLzAxo92FakwnaKxJ0lA/t+taiki01dTUTJyZO667u3vixKzHHnuM+vp6rr/++uVo3pJZyuFxH7B+0v2uymMT3P1+d9/p7jvb29vntJEtrfVsaquje3Xd3FsqIlXlscce47bbbqOnp4cvfelL/Mmf/Ak7duyY9bQPF6OlHPn/DNhqZpsoh/5dwK8t+FbKVR8srh2+IivZs7c/uyjrvepvrzrn8mw2y44dO4DyRVi+8Y1vTCzr7u7mt3/7t6mvr+cTn/jEorRvpViy8Hf3kpn9LvB3QBz4S3d/fsG3U7mKl8JfRGYyU9mnGi1pzd/dvwN8Z1G3UVL4i1wMzjdCl8UVuUNiJkb+mt5BROagoaGB06dPL3czFl30wr8y8kfXcRGRObj99tv5xje+oR2+FxuN/EXkXEZGRs547Oabb+bmm28GYNu2bTzzzDNL3KqlF7mRvyWMeFOceK2G/iIiZxO5kf+qX17Fql9etdzNEBFZ0SI38hcRkfNT+IvIklqpl4692Mz376jwF5Elk8lkGBgY0AfAPLk7AwMDZDKZOa8jcjV/EVm5urq66O3tpb+/f7mbctHLZDJ0dXXN+fUKfxFZMslkkk2bNi13MwSVfUREqpLCX0SkCin8RUSqkK3Uve5m1g8cmMNL24ATC9yclUZ9vPhFvX8Q/T6u1P5tdPfzXg1rxYb/XJnZbnffudztWEzq48Uv6v2D6PfxYu+fyj4iIlVI4S8iUoWiGP73L3cDloD6ePGLev8g+n28qPsXuZq/iIicXxRH/iIich6RCn8zu9XMXjSzV8zsk8vdnoViZj1m9qyZPWVmuyuPrTKzvzezlyu/W5a7nbNlZn9pZsfN7LlJj83YHyv7XOU9fcbM3rR8LZ+9s/TxD8ysr/I+PmVm75m07PcqfXzRzG5ZnlbPnpmtN7NHzWyvmT1vZvdVHo/M+3iOPkbjfXT3SPxQvmrvq8BmIAU8DVy+3O1aoL71AG3THvsPwCcrtz8J/PFyt/MC+vN24E3Ac+frD/Ae4GHAgOuAJ5a7/fPo4x8An5jhuZdX/r2mgU2Vf8fx5e7Defq3FnhT5XYD8FKlH5F5H8/Rx0i8j1Ea+V8LvOLu+929ADwE3LHMbVpMdwD/rXL7vwG/soxtuSDu/n+Ak9MePlt/7gD+yst+CjSb2dqlaencnaWPZ3MH8JC75939NeAVyv+eVyx3P+LuP6/cPg3sAzqJ0Pt4jj6ezUX1PkYp/DuBQ5Pu93LuN+pi4sD3zGyPmd1beWy1ux+p3D4KrF6epi2Ys/Unau/r71bKHn85qVR3UffRzLqBq4EniOj7OK2PEIH3MUrhH2U3uvubgHcDv2Nmb5+80MvfOSNz2FbU+jPJnwOXADuAI8B/Wt7mzJ+Z1QP/E/jn7j48eVlU3scZ+hiJ9zFK4d8HrJ90v6vy2EXP3fsqv48D36D8VfLY+Nfmyu/jy9fCBXG2/kTmfXX3Y+4euHsI/H+8XhK4KPtoZknKofigu/9N5eFIvY8z9TEq72OUwv9nwFYz22RmKeAu4NvL3KZ5M7M6M2sYvw28C3iOct8+UnnaR4BvLU8LF8zZ+vNt4DcrR4tcB5yaVFa4qEyrcb+f8vsI5T7eZWZpM9sEbAWeXOr2XQgzM+ABYJ+7/+dJiyLzPp6tj5F5H5d7j/NC/lA+ouAlynvZ//Vyt2eB+rSZ8hEETwPPj/cLaAW+D7wMPAKsWu62XkCfdlH+ulykXBf9R2frD+WjQ75YeU+fBXYud/vn0cevVPrwDOWgWDvp+f+60scXgXcvd/tn0b8bKZd0ngGeqvy8J0rv4zn6GIn3UWf4iohUoSiVfUREZJYU/iIiVUjhLyJShRT+IiJVSOEvIlKFFP4iIlVI4S8iUoUU/lIVzOzzZnZgnuv4J2bmZnbZpMf2Vc7mnMv6/uukOeGfMrOjZjbbmUBF5kXhL5FXmZHxl4DU+FQZc3QV5bM831tZb4byrJU9Z9nuzWb25bOtzN1/y913uPsOytMElIC759E+kVlT+Es1+EPg3wF7gSvmsZ7twB9TCX/KF+94wed5mryZtQHfBf6tu1/081HJxUHhL5FmZlcAVwL/nfLFOK6c4Tk/mlZ+Gf/55WlPvZzyRGUdZtZE+ZvAM/NsXy3wt8D/cPe/mM+6RC5EYrkbILLI/h3w++7uZraPGUb+7v62863EzNYDA+6eNbO/B26h/E3g2Rme+wTlS/nVA6vM7KnKov/H3f9u0vPilK8494K7f+rCuyYydwp/iSwzewtwK3C1mX0RyDBzWP+I8jVap/uEuz9SuX3VpNd+B/h1ytd4/eb0F7n7WyrrvRm4293vPksT/wxIAv94dj0SWTgKf4myfw/cPh7gZrYa+MX0J81m5M/UUf4Pgb8Aapjhw2Q2zOzTwJuBm929NJd1iMyHav4SSZV6fWrSyB13PwbUm9mqOaxyYuTv7nnKtf6Cuw/NoW3dwB9Qnvv+8Un7GP77HNolMieaz19EpApp5C8iUoUU/iIiVUjhLyJShRT+IiJVSOEvIlKFFP4iIlVI4S8iUoUU/iIiVej/B43ExyNcthy2AAAAAElFTkSuQmCC\n",
- "text/plain": [
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
"source": [
"Masses['Eapprox'] = ytilde\n",
"# Generate a plot comparing the experimental with the fitted values values.\n",
@@ -1563,7 +958,9 @@
{
"cell_type": "code",
"execution_count": 5,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"def R2(y_data, y_model):\n",
@@ -1580,16 +977,10 @@
{
"cell_type": "code",
"execution_count": 6,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "0.9547578478889096\n"
- ]
- }
- ],
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
"source": [
"print(R2(Energies,ytilde))"
]
@@ -1604,16 +995,10 @@
{
"cell_type": "code",
"execution_count": 7,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "0.037875961483052376\n"
- ]
- }
- ],
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
"source": [
"def MSE(y_data,y_model):\n",
" n = np.size(y_model)\n",
@@ -1632,78 +1017,10 @@
{
"cell_type": "code",
"execution_count": 8,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "A \n",
- "1 0 inf\n",
- "2 1 1.123190\n",
- "3 2 0.327631\n",
- "4 6 0.344172\n",
- "5 9 0.044402\n",
- "6 14 0.076899\n",
- "7 19 0.091110\n",
- "8 24 0.090653\n",
- "9 29 0.033010\n",
- "10 34 0.060536\n",
- "11 40 0.021348\n",
- "12 46 0.057821\n",
- "13 52 0.012456\n",
- "14 57 0.002538\n",
- "15 64 0.011360\n",
- "16 72 0.033230\n",
- "17 78 0.006406\n",
- "18 85 0.014667\n",
- "19 93 0.022515\n",
- "20 102 0.001432\n",
- "21 110 0.013773\n",
- "22 118 0.007012\n",
- "23 128 0.009449\n",
- "24 137 0.003110\n",
- "25 146 0.006650\n",
- "26 154 0.001906\n",
- "27 164 0.002782\n",
- "28 174 0.006990\n",
- "29 183 0.003376\n",
- "30 192 0.008386\n",
- " ... \n",
- "238 3089 0.000277\n",
- "239 3099 0.000507\n",
- "240 3109 0.000072\n",
- "241 3118 0.000309\n",
- "242 3127 0.000048\n",
- "243 3136 0.000287\n",
- "244 3144 0.000074\n",
- "245 3154 0.000242\n",
- "246 3162 0.000222\n",
- "247 3170 0.000012\n",
- "248 3177 0.000331\n",
- "249 3186 0.000059\n",
- "250 3194 0.000165\n",
- "251 3201 0.000063\n",
- "252 3209 0.000285\n",
- "253 3216 0.000087\n",
- "254 3224 0.000199\n",
- "255 3232 0.000294\n",
- "256 3241 0.000068\n",
- "257 3248 0.000319\n",
- "258 3256 0.000969\n",
- "259 3264 0.001239\n",
- "260 3275 0.008019\n",
- "261 3280 0.003034\n",
- "262 3289 0.006113\n",
- "264 3304 0.009911\n",
- "265 3310 0.009154\n",
- "266 3317 0.007824\n",
- "269 3338 0.011347\n",
- "270 3344 0.009790\n",
- "Name: Ebinding, Length: 267, dtype: float64\n"
- ]
- }
- ],
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
"source": [
"def RelativeError(y_data,y_model):\n",
" return abs((y_data-y_model)/y_data)\n",
@@ -2079,7 +1396,9 @@
{
"cell_type": "code",
"execution_count": 9,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"# Common imports\n",
@@ -2198,7 +1517,9 @@
{
"cell_type": "code",
"execution_count": 10,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import os\n",
@@ -2278,10 +1599,14 @@
"invert the matrix $\\boldsymbol{X}^T\\boldsymbol{X}$. Using a polynomial expansion as we\n",
"did both for the masses and the fitting of the equation of state,\n",
"leads to row vectors of the design matrix which are essentially\n",
- "orthogonal due to the polynomial character of our model. This may\n",
+ "orthogonal due to the polynomial character of our model. Obtaining the inverse of the design matrix is then often done via a so-called LU, QR or Cholesky decomposition. \n",
+ "More material to come here. \n",
+ "\n",
+ "\n",
+ "This may\n",
"however not the be case in general and a standard matrix inversion\n",
- "algorithm based on say LU decomposition may lead to singularities. We will see an example of this below when we try to fit\n",
- "the coupling constant of the widely used Ising model. \n",
+ "algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below.\n",
+ "\n",
"There is however a way to partially circumvent this problem and also gain some insight about the ordinary least squares approach. \n",
"\n",
"This is given by the **Singular Value Decomposition** algorithm, perhaps\n",
@@ -2291,6 +1616,9 @@
"\n",
"\n",
"\n",
+ "\n",
+ "\n",
+ "\n",
"## The Ising model\n",
"\n",
"The one-dimensional Ising model with nearest neighbor interaction, no\n",
@@ -2328,7 +1656,9 @@
{
"cell_type": "code",
"execution_count": 11,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import numpy as np\n",
@@ -2441,7 +1771,9 @@
{
"cell_type": "code",
"execution_count": 12,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"X = np.zeros((n, L ** 2))\n",
@@ -2505,7 +1837,9 @@
{
"cell_type": "code",
"execution_count": 13,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"X_train_own = np.concatenate(\n",
@@ -2521,7 +1855,9 @@
{
"cell_type": "code",
"execution_count": 14,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"def ols_inv(x: np.ndarray, y: np.ndarray) -> np.ndarray:\n",
@@ -2606,7 +1942,9 @@
{
"cell_type": "code",
"execution_count": 15,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"def ols_svd(x: np.ndarray, y: np.ndarray) -> np.ndarray:\n",
@@ -2617,7 +1955,9 @@
{
"cell_type": "code",
"execution_count": 16,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"beta = ols_svd(X_train_own,y_train)"
@@ -2633,7 +1973,9 @@
{
"cell_type": "code",
"execution_count": 17,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"J = beta[1:].reshape(L, L)"
@@ -2649,7 +1991,9 @@
{
"cell_type": "code",
"execution_count": 18,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"fig = plt.figure(figsize=(20, 14))\n",
@@ -5038,7 +4382,9 @@
{
"cell_type": "code",
"execution_count": 19,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"from numpy import *\n",
@@ -5179,7 +4525,9 @@
{
"cell_type": "code",
"execution_count": 20,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"from numpy import *\n",
@@ -5237,7 +4585,9 @@
{
"cell_type": "code",
"execution_count": 21,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import numpy as np\n",
@@ -5466,7 +4816,9 @@
{
"cell_type": "code",
"execution_count": 22,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
@@ -5535,7 +4887,9 @@
{
"cell_type": "code",
"execution_count": 23,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
@@ -5632,7 +4986,9 @@
{
"cell_type": "code",
"execution_count": 24,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"\"\"\"\n",
@@ -5749,7 +5105,9 @@
{
"cell_type": "code",
"execution_count": 25,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import numpy as np\n",
@@ -5855,7 +5213,9 @@
{
"cell_type": "code",
"execution_count": 26,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"X = np.zeros((n, L ** 2))\n",
@@ -5885,7 +5245,9 @@
{
"cell_type": "code",
"execution_count": 27,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"clf = skl.LinearRegression().fit(X_train, y_train)"
@@ -5901,7 +5263,9 @@
{
"cell_type": "code",
"execution_count": 28,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"J_sk = clf.coef_.reshape(L, L)"
@@ -5917,7 +5281,9 @@
{
"cell_type": "code",
"execution_count": 29,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"fig = plt.figure(figsize=(20, 14))\n",
@@ -5973,7 +5339,9 @@
{
"cell_type": "code",
"execution_count": 30,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"_lambda = 0.1\n",
@@ -6024,7 +5392,9 @@
{
"cell_type": "code",
"execution_count": 31,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"clf_lasso = skl.Lasso(alpha=_lambda).fit(X_train, y_train)\n",
@@ -6058,7 +5428,9 @@
{
"cell_type": "code",
"execution_count": 32,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"lambdas = np.logspace(-4, 5, 10)\n",
@@ -6121,7 +5493,9 @@
{
"cell_type": "code",
"execution_count": 33,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"fig = plt.figure(figsize=(20, 14))\n",
@@ -6177,7 +5551,9 @@
{
"cell_type": "code",
"execution_count": 34,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"x = np.random.rand(100,1)\n",
@@ -6308,7 +5684,9 @@
{
"cell_type": "code",
"execution_count": 35,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"x = np.random.rand(100,1)\n",
@@ -6432,7 +5810,9 @@
{
"cell_type": "code",
"execution_count": 36,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"from mpl_toolkits.mplot3d import Axes3D\n",
@@ -6554,25 +5934,7 @@
]
}
],
- "metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "language": "python",
- "name": "python3"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 3
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython3",
- "version": "3.7.3"
- }
- },
+ "metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}
diff --git a/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz b/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz
index d919eb94b..8f021628a 100644
Binary files a/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz and b/doc/pub/Regression/ipynb/ipynb-Regression-src.tar.gz differ
diff --git a/doc/pub/Regression/pdf/Regression-minted.pdf b/doc/pub/Regression/pdf/Regression-minted.pdf
index 33987bff5..3ad30e25a 100644
Binary files a/doc/pub/Regression/pdf/Regression-minted.pdf and b/doc/pub/Regression/pdf/Regression-minted.pdf differ
diff --git a/doc/src/Regression/Regression.do.txt b/doc/src/Regression/Regression.do.txt
index c377ef3e0..b77215506 100644
--- a/doc/src/Regression/Regression.do.txt
+++ b/doc/src/Regression/Regression.do.txt
@@ -974,10 +974,14 @@ The examples we have looked at so far are cases where we normally can
invert the matrix $\bm{X}^T\bm{X}$. Using a polynomial expansion as we
did both for the masses and the fitting of the equation of state,
leads to row vectors of the design matrix which are essentially
-orthogonal due to the polynomial character of our model. This may
+orthogonal due to the polynomial character of our model. Obtaining the inverse of the design matrix is then often done via a so-called LU, QR or Cholesky decomposition.
+More material to come here.
+
+
+This may
however not the be case in general and a standard matrix inversion
-algorithm based on say LU decomposition may lead to singularities. We will see an example of this below when we try to fit
-the coupling constant of the widely used Ising model.
+algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below.
+
There is however a way to partially circumvent this problem and also gain some insight about the ordinary least squares approach.
This is given by the _Singular Value Decomposition_ algorithm, perhaps
@@ -987,6 +991,9 @@ inversion algorithm. Thereafter we dive into the math of the SVD.
!eblock
+
+# todo: change model here.
+
!split
===== The Ising model =====