From 31a7c6e7607d6e55b81a15af3fd44752d65caf84 Mon Sep 17 00:00:00 2001 From: mhjensen Date: Thu, 2 Jan 2020 10:23:52 +0100 Subject: [PATCH] more pca --- doc/pub/DimRed/html/._DimRed-bs000.html | 2 +- doc/pub/DimRed/html/._DimRed-bs018.html | 49 ++++++----- doc/pub/DimRed/html/DimRed-bs.html | 2 +- doc/pub/DimRed/html/DimRed-reveal.html | 51 +++++++----- doc/pub/DimRed/html/DimRed-solarized.html | 51 +++++++----- doc/pub/DimRed/html/DimRed.html | 51 +++++++----- doc/pub/DimRed/ipynb/DimRed.ipynb | 81 ++++++++++++------- doc/pub/DimRed/ipynb/ipynb-DimRed-src.tar.gz | Bin 191 -> 190 bytes doc/pub/DimRed/pdf/DimRed-minted.pdf | Bin 263267 -> 263343 bytes doc/src/DimRed/DimRed.do.txt | 46 +++++++---- doc/src/DimRed/Programs/PCAsimple.py | 41 +++++++--- 11 files changed, 241 insertions(+), 133 deletions(-) diff --git a/doc/pub/DimRed/html/._DimRed-bs000.html b/doc/pub/DimRed/html/._DimRed-bs000.html index 59e9e22ca..2fdf5d407 100644 --- a/doc/pub/DimRed/html/._DimRed-bs000.html +++ b/doc/pub/DimRed/html/._DimRed-bs000.html @@ -224,7 +224,7 @@ MathJax.Hub.Config({
[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

-

Jan 1, 2020

+

Jan 2, 2020


diff --git a/doc/pub/DimRed/html/._DimRed-bs018.html b/doc/pub/DimRed/html/._DimRed-bs018.html index 1dc23b946..26757403b 100644 --- a/doc/pub/DimRed/html/._DimRed-bs018.html +++ b/doc/pub/DimRed/html/._DimRed-bs018.html @@ -233,12 +233,7 @@ n = 100 mean = (-1, 2) cov = [[4, 2], [2, 2]] X = np.random.multivariate_normal(mean, cov, n) -# Print the X-matrix -print(X) -

-Try to add to this code your own calculation of the covariance matrix. -

Now we are going to implement the PCA algorithm. We will break it down into various substeps. @@ -258,31 +253,31 @@ $$ When you are done with these steps, print out \( \mu_n \) to verify it is close to \( \mu \) and plot your mean centered data to verify it is centered at the origin! Compare your code with the functionality from Scikit-Learn discussed above. -The following code elements perform these operations using pandas or using our own functionality for doing so. The latter, using numpy is rather simply through the mean() function. +The following code elements perform these operations using pandas or using our own functionality for doing so. The latter, using numpy is rather simple through the mean() function. 

df = pd.DataFrame(X)
 # Pandas does the centering for us
 df = df -df.mean()
-display(df)
-
 # we center it ourselves
 X_centered = X - X.mean(axis=0)
-# test that we get the same as Pandas
-print(X_centered-df)

-Alternatively, you could also have used the functions we discussed earlier for scaling the data set. -That is, we could have used the StandardScaler function in Scikit-Learn, a function which -ensures that for each feature/predictor we study the mean value is -zero and the variance is one (every column in the design/feature -matrix). You would then not get the same results, since we divide by the variance. For diagonal covariance matrix elements will then be one, while the non-diagonal ones will be divided by \( 2\sqrt{2} \) for our specific case. +Alternatively, we could use the functions we discussed +earlier for scaling the data set. That is, we could have used the +StandardScaler function in Scikit-Learn, a function which ensures +that for each feature/predictor we study the mean value is zero and +the variance is one (every column in the design/feature matrix). You +would then not get the same results, since we divide by the +variance. The diagonal covariance matrix elements will then be one, +while the non-diagonal ones need to be divided by \( 2\sqrt{2} \) for our +specific case.

Compute the sample covariance

-Now we are going to use the mean centered data to compute the sample covariance of the data. Recall it is given by: +Now we are going to use the mean centered data to compute the sample covariance of the data. $$ \begin{equation*} \Sigma_n = \frac{1}{n-1} \sum_{i=1}^n \bar{x}_i^T \bar{x}_i = \frac{1}{n-1} \sum_{i=1}^n (x_i - \mu_n)^T (x_i - \mu_n) @@ -294,11 +289,29 @@ We can write our own code or simply use either the functionaly of numpy o

-

print(np.cov(X.T))
-print(df.cov())
+
print(df.cov())
 print(np.cov(X_centered.T))
 

 

+Note that the way we define the covariance matrix here has a factor \( n-1 \) instead of \( n \).
+Our own code here is not very elegant and asks for improvements. 
+

+
+
+
# extract the relevant columns from the centered design matrix of dim n x 2                                                      x = X_centered[:,[0]]
+y = X_centered[:,[1]]
+Cov = np.zeros((2,2))
+Cov[0,1] = np.sum(x.T@y)/(n-1.0)
+Cov[0,0] = np.sum(x.T@x)/(n-1.0)
+Cov[1,1] = np.sum(y.T@y)/(n-1.0)
+Cov[1,0]= Cov[0,1]
+print("Centered covariance using own code")
+print(Cov)
+plt.plot(x, y, 'x')
+plt.axis('equal')
+plt.show()
+

+

 Depending on the number of points \( n \), we will get results that are close to the covariance values defined above.
 
 
Diagonalize the sample covariance matrix to obtain the principal components 

diff --git a/doc/pub/DimRed/html/DimRed-bs.html b/doc/pub/DimRed/html/DimRed-bs.html
index 59e9e22ca..2fdf5d407 100644
--- a/doc/pub/DimRed/html/DimRed-bs.html
+++ b/doc/pub/DimRed/html/DimRed-bs.html
@@ -224,7 +224,7 @@ MathJax.Hub.Config({
 
[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

 

 

-
Jan 1, 2020
 
+
Jan 2, 2020
 
 

 

 
diff --git a/doc/pub/DimRed/html/DimRed-reveal.html b/doc/pub/DimRed/html/DimRed-reveal.html
index d40e41dbd..5a9a90a8f 100644
--- a/doc/pub/DimRed/html/DimRed-reveal.html
+++ b/doc/pub/DimRed/html/DimRed-reveal.html
@@ -148,7 +148,7 @@ MathJax.Hub.Config({
 
[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

 

 
 

-
Jan 1, 2020
 
+
Jan 2, 2020
 
 

 

 
@@ -1039,12 +1039,7 @@ n = 100
 mean = (-1, 2)
 cov = [[4, 2], [2, 2]]
 X = np.random.multivariate_normal(mean, cov, n)
-# Print the X-matrix
-print(X)
-

-Try to add to this code your own calculation of the covariance matrix. -

Now we are going to implement the PCA algorithm. We will break it down into various substeps. @@ -1068,31 +1063,31 @@ $$ When you are done with these steps, print out \( \mu_n \) to verify it is close to \( \mu \) and plot your mean centered data to verify it is centered at the origin! Compare your code with the functionality from Scikit-Learn discussed above. -The following code elements perform these operations using pandas or using our own functionality for doing so. The latter, using numpy is rather simply through the mean() function. +The following code elements perform these operations using pandas or using our own functionality for doing so. The latter, using numpy is rather simple through the mean() function. 

df = pd.DataFrame(X)
 # Pandas does the centering for us
 df = df -df.mean()
-display(df)
-
 # we center it ourselves
 X_centered = X - X.mean(axis=0)
-# test that we get the same as Pandas
-print(X_centered-df)

-Alternatively, you could also have used the functions we discussed earlier for scaling the data set. -That is, we could have used the StandardScaler function in Scikit-Learn, a function which -ensures that for each feature/predictor we study the mean value is -zero and the variance is one (every column in the design/feature -matrix). You would then not get the same results, since we divide by the variance. For diagonal covariance matrix elements will then be one, while the non-diagonal ones will be divided by \( 2\sqrt{2} \) for our specific case. +Alternatively, we could use the functions we discussed +earlier for scaling the data set. That is, we could have used the +StandardScaler function in Scikit-Learn, a function which ensures +that for each feature/predictor we study the mean value is zero and +the variance is one (every column in the design/feature matrix). You +would then not get the same results, since we divide by the +variance. The diagonal covariance matrix elements will then be one, +while the non-diagonal ones need to be divided by \( 2\sqrt{2} \) for our +specific case.

Compute the sample covariance

-Now we are going to use the mean centered data to compute the sample covariance of the data. Recall it is given by: +Now we are going to use the mean centered data to compute the sample covariance of the data.

 
$$ \begin{equation*} @@ -1106,11 +1101,29 @@ We can write our own code or simply use either the functionaly of numpy o

-

print(np.cov(X.T))
-print(df.cov())
+
print(df.cov())
 print(np.cov(X_centered.T))
 

 

+Note that the way we define the covariance matrix here has a factor \( n-1 \) instead of \( n \).
+Our own code here is not very elegant and asks for improvements. 
+

+
+
+
# extract the relevant columns from the centered design matrix of dim n x 2                                                      x = X_centered[:,[0]]
+y = X_centered[:,[1]]
+Cov = np.zeros((2,2))
+Cov[0,1] = np.sum(x.T@y)/(n-1.0)
+Cov[0,0] = np.sum(x.T@x)/(n-1.0)
+Cov[1,1] = np.sum(y.T@y)/(n-1.0)
+Cov[1,0]= Cov[0,1]
+print("Centered covariance using own code")
+print(Cov)
+plt.plot(x, y, 'x')
+plt.axis('equal')
+plt.show()
+

+

 Depending on the number of points \( n \), we will get results that are close to the covariance values defined above.
 
 
Diagonalize the sample covariance matrix to obtain the principal components 

diff --git a/doc/pub/DimRed/html/DimRed-solarized.html b/doc/pub/DimRed/html/DimRed-solarized.html
index a15038ff7..0c1ae2da3 100644
--- a/doc/pub/DimRed/html/DimRed-solarized.html
+++ b/doc/pub/DimRed/html/DimRed-solarized.html
@@ -179,7 +179,7 @@ MathJax.Hub.Config({
 
[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

 

 

-
Jan 1, 2020
 
+
Jan 2, 2020
 
 

 

 










@@ -1019,12 +1019,7 @@ n = 100
 mean = (-1, 2)
 cov = [[4, 2], [2, 2]]
 X = np.random.multivariate_normal(mean, cov, n)
-# Print the X-matrix
-print(X)
-

-Try to add to this code your own calculation of the covariance matrix. -

Now we are going to implement the PCA algorithm. We will break it down into various substeps. @@ -1044,31 +1039,31 @@ $$ When you are done with these steps, print out \( \mu_n \) to verify it is close to \( \mu \) and plot your mean centered data to verify it is centered at the origin! Compare your code with the functionality from Scikit-Learn discussed above. -The following code elements perform these operations using pandas or using our own functionality for doing so. The latter, using numpy is rather simply through the mean() function. +The following code elements perform these operations using pandas or using our own functionality for doing so. The latter, using numpy is rather simple through the mean() function. 

df = pd.DataFrame(X)
 # Pandas does the centering for us
 df = df -df.mean()
-display(df)
-
 # we center it ourselves
 X_centered = X - X.mean(axis=0)
-# test that we get the same as Pandas
-print(X_centered-df)

-Alternatively, you could also have used the functions we discussed earlier for scaling the data set. -That is, we could have used the StandardScaler function in Scikit-Learn, a function which -ensures that for each feature/predictor we study the mean value is -zero and the variance is one (every column in the design/feature -matrix). You would then not get the same results, since we divide by the variance. For diagonal covariance matrix elements will then be one, while the non-diagonal ones will be divided by \( 2\sqrt{2} \) for our specific case. +Alternatively, we could use the functions we discussed +earlier for scaling the data set. That is, we could have used the +StandardScaler function in Scikit-Learn, a function which ensures +that for each feature/predictor we study the mean value is zero and +the variance is one (every column in the design/feature matrix). You +would then not get the same results, since we divide by the +variance. The diagonal covariance matrix elements will then be one, +while the non-diagonal ones need to be divided by \( 2\sqrt{2} \) for our +specific case.

Compute the sample covariance

-Now we are going to use the mean centered data to compute the sample covariance of the data. Recall it is given by: +Now we are going to use the mean centered data to compute the sample covariance of the data. $$ \begin{equation*} \Sigma_n = \frac{1}{n-1} \sum_{i=1}^n \bar{x}_i^T \bar{x}_i = \frac{1}{n-1} \sum_{i=1}^n (x_i - \mu_n)^T (x_i - \mu_n) @@ -1080,11 +1075,29 @@ We can write our own code or simply use either the functionaly of numpy o

-

print(np.cov(X.T))
-print(df.cov())
+
print(df.cov())
 print(np.cov(X_centered.T))
 

 

+Note that the way we define the covariance matrix here has a factor \( n-1 \) instead of \( n \).
+Our own code here is not very elegant and asks for improvements. 
+

+
+
+
# extract the relevant columns from the centered design matrix of dim n x 2                                                      x = X_centered[:,[0]]
+y = X_centered[:,[1]]
+Cov = np.zeros((2,2))
+Cov[0,1] = np.sum(x.T@y)/(n-1.0)
+Cov[0,0] = np.sum(x.T@x)/(n-1.0)
+Cov[1,1] = np.sum(y.T@y)/(n-1.0)
+Cov[1,0]= Cov[0,1]
+print("Centered covariance using own code")
+print(Cov)
+plt.plot(x, y, 'x')
+plt.axis('equal')
+plt.show()
+

+

 Depending on the number of points \( n \), we will get results that are close to the covariance values defined above.
 
 
Diagonalize the sample covariance matrix to obtain the principal components 

diff --git a/doc/pub/DimRed/html/DimRed.html b/doc/pub/DimRed/html/DimRed.html
index 424c95f51..479e67b4f 100644
--- a/doc/pub/DimRed/html/DimRed.html
+++ b/doc/pub/DimRed/html/DimRed.html
@@ -184,7 +184,7 @@ MathJax.Hub.Config({
 
[2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

 

 

-
Jan 1, 2020
 
+
Jan 2, 2020
 
 

 

 










@@ -1024,12 +1024,7 @@ n = 100
 mean = (-1, 2)
 cov = [[4, 2], [2, 2]]
 X = np.random.multivariate_normal(mean, cov, n)
-# Print the X-matrix
-print(X)
-

-Try to add to this code your own calculation of the covariance matrix. -

Now we are going to implement the PCA algorithm. We will break it down into various substeps. @@ -1049,31 +1044,31 @@ $$ When you are done with these steps, print out \( \mu_n \) to verify it is close to \( \mu \) and plot your mean centered data to verify it is centered at the origin! Compare your code with the functionality from Scikit-Learn discussed above. -The following code elements perform these operations using pandas or using our own functionality for doing so. The latter, using numpy is rather simply through the mean() function. +The following code elements perform these operations using pandas or using our own functionality for doing so. The latter, using numpy is rather simple through the mean() function. 

df = pd.DataFrame(X)
 # Pandas does the centering for us
 df = df -df.mean()
-display(df)
-
 # we center it ourselves
 X_centered = X - X.mean(axis=0)
-# test that we get the same as Pandas
-print(X_centered-df)

-Alternatively, you could also have used the functions we discussed earlier for scaling the data set. -That is, we could have used the StandardScaler function in Scikit-Learn, a function which -ensures that for each feature/predictor we study the mean value is -zero and the variance is one (every column in the design/feature -matrix). You would then not get the same results, since we divide by the variance. For diagonal covariance matrix elements will then be one, while the non-diagonal ones will be divided by \( 2\sqrt{2} \) for our specific case. +Alternatively, we could use the functions we discussed +earlier for scaling the data set. That is, we could have used the +StandardScaler function in Scikit-Learn, a function which ensures +that for each feature/predictor we study the mean value is zero and +the variance is one (every column in the design/feature matrix). You +would then not get the same results, since we divide by the +variance. The diagonal covariance matrix elements will then be one, +while the non-diagonal ones need to be divided by \( 2\sqrt{2} \) for our +specific case.

Compute the sample covariance

-Now we are going to use the mean centered data to compute the sample covariance of the data. Recall it is given by: +Now we are going to use the mean centered data to compute the sample covariance of the data. $$ \begin{equation*} \Sigma_n = \frac{1}{n-1} \sum_{i=1}^n \bar{x}_i^T \bar{x}_i = \frac{1}{n-1} \sum_{i=1}^n (x_i - \mu_n)^T (x_i - \mu_n) @@ -1085,11 +1080,29 @@ We can write our own code or simply use either the functionaly of numpy o

-

print(np.cov(X.T))
-print(df.cov())
+
print(df.cov())
 print(np.cov(X_centered.T))
 

 

+Note that the way we define the covariance matrix here has a factor \( n-1 \) instead of \( n \).
+Our own code here is not very elegant and asks for improvements. 
+

+
+
+
# extract the relevant columns from the centered design matrix of dim n x 2                                                      x = X_centered[:,[0]]
+y = X_centered[:,[1]]
+Cov = np.zeros((2,2))
+Cov[0,1] = np.sum(x.T@y)/(n-1.0)
+Cov[0,0] = np.sum(x.T@x)/(n-1.0)
+Cov[1,1] = np.sum(y.T@y)/(n-1.0)
+Cov[1,0]= Cov[0,1]
+print("Centered covariance using own code")
+print(Cov)
+plt.plot(x, y, 'x')
+plt.axis('equal')
+plt.show()
+

+

 Depending on the number of points \( n \), we will get results that are close to the covariance values defined above.
 
 
Diagonalize the sample covariance matrix to obtain the principal components 

diff --git a/doc/pub/DimRed/ipynb/DimRed.ipynb b/doc/pub/DimRed/ipynb/DimRed.ipynb
index 42b7a9f71..bded76346 100644
--- a/doc/pub/DimRed/ipynb/DimRed.ipynb
+++ b/doc/pub/DimRed/ipynb/DimRed.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: **Jan 1, 2020**\n",
+    "Date: **Jan 2, 2020**\n",
     "\n",
     "Copyright 1999-2020, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n",
     "\n",
@@ -1111,17 +1111,13 @@
     "n = 100\n",
     "mean = (-1, 2)\n",
     "cov = [[4, 2], [2, 2]]\n",
-    "X = np.random.multivariate_normal(mean, cov, n)\n",
-    "# Print the X-matrix\n",
-    "print(X)"
+    "X = np.random.multivariate_normal(mean, cov, n)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Try to add to this code your own calculation of the covariance matrix.\n",
-    "\n",
     "Now we are going to implement the PCA algorithm. We will break it down into various substeps.\n",
     "\n",
     "### Compute the sample mean and center the data\n",
@@ -1161,7 +1157,7 @@
     "When you are done with these steps, print out $\\mu_n$ to verify it is\n",
     "close to $\\mu$ and plot your mean centered data to verify it is\n",
     "centered at the origin! Compare your code with the functionality from **Scikit-Learn** discussed above.\n",
-    "The following code elements perform these operations using **pandas** or using our own functionality for doing so. The latter, using **numpy** is rather simply through the **mean()** function."
+    "The following code elements perform these operations using **pandas** or using our own functionality for doing so. The latter, using **numpy** is rather simple through the **mean()** function."
    ]
   },
   {
@@ -1175,27 +1171,27 @@
     "df = pd.DataFrame(X)\n",
     "# Pandas does the centering for us\n",
     "df = df -df.mean()\n",
-    "display(df)\n",
-    "\n",
     "# we center it ourselves\n",
-    "X_centered = X - X.mean(axis=0)\n",
-    "# test that we get the same as Pandas\n",
-    "print(X_centered-df)"
+    "X_centered = X - X.mean(axis=0)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Alternatively, you could also have used the functions we discussed earlier for scaling the data set. \n",
-    "That is, we could have used the **StandardScaler** function in **Scikit-Learn**, a function which \n",
-    "ensures that for each feature/predictor we study the mean value is\n",
-    "zero and the variance is one (every column in the design/feature\n",
-    "matrix).  You would then not get the same results, since we divide by the variance. For diagonal covariance matrix elements will then be one, while the non-diagonal ones will be divided by $2\\sqrt{2}$ for our specific case.\n",
+    "Alternatively, we could use the functions we discussed\n",
+    "earlier for scaling the data set.  That is, we could have used the\n",
+    "**StandardScaler** function in **Scikit-Learn**, a function which ensures\n",
+    "that for each feature/predictor we study the mean value is zero and\n",
+    "the variance is one (every column in the design/feature matrix).  You\n",
+    "would then not get the same results, since we divide by the\n",
+    "variance. The diagonal covariance matrix elements will then be one,\n",
+    "while the non-diagonal ones need to be divided by $2\\sqrt{2}$ for our\n",
+    "specific case.\n",
     "\n",
     "### Compute the sample covariance\n",
     "\n",
-    "Now we are going to use the mean centered data to compute the sample covariance of the data. Recall it is given by:"
+    "Now we are going to use the mean centered data to compute the sample covariance of the data."
    ]
   },
   {
@@ -1223,11 +1219,40 @@
    },
    "outputs": [],
    "source": [
-    "print(np.cov(X.T))\n",
     "print(df.cov())\n",
     "print(np.cov(X_centered.T))"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note that the way we define the covariance matrix here has a factor $n-1$ instead of $n$.\n",
+    "Our own code here is not very elegant and asks for improvements."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# extract the relevant columns from the centered design matrix of dim n x 2                                                      x = X_centered[:,[0]]\n",
+    "y = X_centered[:,[1]]\n",
+    "Cov = np.zeros((2,2))\n",
+    "Cov[0,1] = np.sum(x.T@y)/(n-1.0)\n",
+    "Cov[0,0] = np.sum(x.T@x)/(n-1.0)\n",
+    "Cov[1,1] = np.sum(y.T@y)/(n-1.0)\n",
+    "Cov[1,0]= Cov[0,1]\n",
+    "print(\"Centered covariance using own code\")\n",
+    "print(Cov)\n",
+    "plt.plot(x, y, 'x')\n",
+    "plt.axis('equal')\n",
+    "plt.show()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -1275,7 +1300,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 15,
    "metadata": {
     "collapsed": false
    },
@@ -1586,7 +1611,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 16,
    "metadata": {
     "collapsed": false
    },
@@ -1633,7 +1658,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 17,
    "metadata": {
     "collapsed": false
    },
@@ -1657,7 +1682,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 18,
    "metadata": {
     "collapsed": false
    },
@@ -1681,7 +1706,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 19,
    "metadata": {
     "collapsed": false
    },
@@ -1705,7 +1730,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": 20,
    "metadata": {
     "collapsed": false
    },
@@ -1759,7 +1784,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 21,
    "metadata": {
     "collapsed": false
    },
@@ -1782,7 +1807,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 21,
+   "execution_count": 22,
    "metadata": {
     "collapsed": false
    },
@@ -1829,7 +1854,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": 23,
    "metadata": {
     "collapsed": false
    },
diff --git a/doc/pub/DimRed/ipynb/ipynb-DimRed-src.tar.gz b/doc/pub/DimRed/ipynb/ipynb-DimRed-src.tar.gz
index 3dcbbe279d3f4a50048718bdc8bc698ae30ff2c9..603419a5eb62e8043df1425cde6a9db54ab251ac 100644
GIT binary patch
literal 190
zcmV;v073sBiwFRlwhdkY1MSbv3c@f92k@Qu6nTP?uG>5d+rfh%!x!ji=c=xqZ3ph{
z-3RDN@iIi{@A4-kgrqFC`M!wU-36;b#4Jf*%1xHeiSbg;2uXp^4CyGQ6p&DwqygZ1
zC%yF6F*BahoQ1MNy`3A!>iWZ;kM811MSaC3c@fD2H>uHia9|^Op~sKcHu&h;ssKY+Ne!xl7hXx
zeSoeMH${Yeo1bBZVWup$`MyZ}-AAiI2(dJUDKkdrM4P3aVN3xNh6rYqW+@;{#wr4|
z-bpXL^W2W7G-s)u)Nkj;vAX`SXL$vl`6rHza7@F%F%4Q&P{!=4#sEqB;)>1B?-uP2polWX8T
z&axkh#~FsY=2b8RCd0IQOO`-`V(jV3A48n=6mI$fY$xj}?BO}IQ<wvx##Q1CG
z$ngw2_!`oquVIa&Ip#}NgIhr4GB;LN-LjhDBGGeur1%^QktUEqCwv|Rue~^8FQy`H6=Hdn1-@c;c`vDxvsgJ5)g(+DDdvd=7s*pefy4t-
z4cx{nG~!9f=>f61*}e~5aKYxjd(T6}7W90`r%mb-0-tNJtD;5I#56N&BhN_C&F;`}
z&$suW(z&yZ<&6-XZKTfw;Cxmw+`p)j^nbonz2~(LyFv4L&oYNW@;8CQY}eof^UP-A29}g*8p2AIj$0Svi^HN=_+FMns@w}NZ#y?e<=cj%`yriZe)sRCJKd=jt7fhR
z{)=WWWnJO4%HTP)hVRX@geMfCj>d&3yw~MYMZ2LaXo_;dKW8?&fj3DbfP;60669?h
z#Z{``Fp1_7BhZ}#_=!TXQK;B#0~dtP2t0Jb)zKUboTtDF!9|^>8|~ssj*}(KCv6H}
z7ZEPv0x@KQlfFtsFH==3i)C)`ZD+Or$|y3UgkggD&vHx`tr4toW7-uOk7Y6ekJZbI
zbo0){8V8nk^W9a%6-Fxm97}L{Iwu;sx5rT0Bk_VDNEnL(^26C73b0FsyOiRPfX|a}
zA#%Wf-zZL}T~2wYxA!R)xY@VYA^EnYCT#OWT%Ut%1A^eafpmT!ePz}PmNw3_{)&}l
z0-0Kv7XxPOnYUU%c-+~K2s&Kb)N$hwhXf~zyhI&Ozhhh(0>dcl9>OB{+){Y+Skvj}
z#+jkgdbfRmkQ7g{87+vr_wvYO3mw14KzYA65tv(I3lh;>rW4W%BWAl;xL5#Rn|oTF
zkJ~4F(qABK&Y*fa$}iSq`>nv66UJUa00?20&d>llk7$^UCEDF(U()y%1~B+;d&GNM
z=8f_}`QzFI!F>}K=O04cnOQtPIQr#KU#E$iN6-gA9B2Lyn=hVN!U0!J<&Ms~Zjfjp
zR8C(_MWv%?ya3cm7C;V$zhWtng&!sspee(;*!(9+L~3bBAI`j+5*GwUFVdX9Dy91K
z8LqHH!xQ#dvZleL4io67u~IEF1fSju@q5AaYqb)CZj4KH9k)oopOP;lf*LyP;{0ix
z-*8R)2@4P|S`vD6kq{wL@Tn6etDgNDd*wL1f~Lig2$sW0!i~(#{Py7#hQj)Boxa<>
zrR>z>F>DgJrD`*ftU!6mFGH)oQUSW*&A97pmjA--{0*OgRx{h81p6Td%9D}w51baz
zt0$kb)q&dgp-sbGWha4-3>MxVu5y(Aawfr|qJFk7+#tx*h^tx{Ve
zo<9-!?%Q7+MG17!+{R@DvGxI7m?ycTezx>?oOtLf!|-(Iylm1mEb>w_no{4^)sIIU
zKV{gKwDPS-yRp>`55)jmv->iyGTq`lr?1|vj2)}#->yQbct(=H6ULuliDr-kWm=Z2
z-5f$zYy+-phc=VUB3b$xvOnHY<+Cno@m{v_Uw-AZ^?$gnB@HvYcNkhkx$9ZQdnbBM
z0ITM%8sHaAQ=HMr8FDSbl=#|BXpm}I4z
zZ&W3C)9^K<-Io>$q|Xu+1S
z1Ac6!uSO&erA^Gg$m-8=hhUH_H+WMGvtHfJ;KVs*I>fawXddn8flZnG*IU1nF6!~N
zAg;i&%RLTIH7rqy-6nB5Ll8}(gjo*3&
z?~9yp+>6@d(77GY=+n%jcSX~3i-SG<7XF|Car$i{$^8@2-+OSJQgaT*JPevoaRxG&F4*
zr36-DL3C(=i7kq*$DIFcop!uOii6$4O;J
z`02)h_H%|Hj0jbPJ_^UFer~u$2^0CpHH?MD$}9j17qEp>B}Cxuq)X56>BHb?It;=u
z>uc*6WOA0kgR}yg3=ttZ!P}3rKT{Tg_N;M?d&i#O=UZvsmZx7&^Idz*>q%C`+i`~a
z1@XV}i(F*;Enp_C=p}>NB%1N9qismZQVFXbEYpU}a&>}O-3Oi(U@_*L{$y3)jz~I+Q!j&gp%Aqm=u=joX%6@stVAlH*wnVs!=Ide*
zGj#YXUG+;OHGqA8ge5gxeWJ~!x}hFy85fIkOoOQa0AEu)Y|97>ycp8fgms?5X^I_B
zXH@&)brn6Xy(t$LA&By8(Db)8H7A!+4ZDEi)K?j7Zoq2$nWS9!lK}%yOeL6q1QcUw5{G9yw2|JqiIwRo$5UBxu
z8?NLX$ggtopu*MP@#28s2PxVC_6@ciO!1W)r9kI>
z9<-c;6!(}@zmOZ@wccEXkuE$y;65NWQB@XLP5ZU2ECF1Y-PG>a{-6cilPl+;Mdssi
z+M1%XvImQj=Z7GP(q>9hY_H5ZL@MwbZz4e^
z(<soSPhgi4m0k2W_YRcC&YxAoj1w8Iazdm{pcKTF_Gz6=9&8wno6roP@CRx#U<2pnqs^WI@{isL2A7RM?&r+iDKrrtRd=vDr
z&L5ft5~PH|yjlgV261z-itV9c2a2I`EdH+X2B5Vw2&V0ZT&iVatxsqQd0UXR{{uDs
z+|3cn8_voY6*(j2HFjfBpaA^xHFiK0FISRF9uyc)Q9B*vE>9`GQjY2LYIx(g`DQr^
z4`#}EEHrIECU+`XmWy;6m$Q}~@#F3Mw(cyMVf9;}4bfuGXz5kPbUyG=sK4bc$6rO4
zU`cK5EU62O`iyhbgpi}zmXD&yw#gm@+rHn-3*EYfOywCDh=oSUjRUY5AXTOfi5mEV
z+&v5Uf5Rwos)jQEnVHK%%omnO7A@C!ogB6W{oa#H+;2p1uxAo3FRjNLN4?D|(2@ip
zY4`liiXIyxEF*75P6C}bi6Wa_Bhws%gI&)BE!5`aSXyvKKBnLt-}ZOaq0W!2Ytfo}
zET)OBNdeRc3B%E)*a0ALtKfYP9HpM5#m$+GBlXDTQ6A-$|2@!wfw9bU$dvs|tB=z?
zg|ui%`h7zqsMZQyz5FSB@Z3c;V&?uUTwehCK5U(z1WNVn0(;7B&1-zZwhfL7Bp{eW
zNK_#3S>b-`uZCm+m3rY_S--E|cMk0=w__WM3D(4o*dNvb?mfVz!h=2c^En^W6DRLq
zi)4IcsNKy|PbANPQaKHOH)5=Dmo!io>BciVJ>xDq=eB+pRk7~o6=kH}5w>msB$7Iia
zhzy^UHF<>a3Zqa4$$|EO!!Q@D2Bx68d*~b;U~H%4c-7TvkMwk}J#EO4Y-%jrrky
z94o)#DG+iZ8<(#cVQ+%Eovi$`uCm57J#W2hC@qy~+5sSktpaGN399*ZwKr94pYL0U
z3G)LdE8qCNF*9BKf-teXc3r34s_P}B=^v3IcWwWVURrp>T_J~pP`>?`H}E$4C@jrR
z_QfQRiEq{>f7!t>FRt3}+Y~gfA|I!Nfi)?c0z0;sxUf(XAlh(j;em&e9ZA3C)K_P@
zCCG{b5CT{?5ubig{=qBfaDS^cIWG1{j;DO%?~eFY{81|Z;JJ4I$Ho}w7ChL6Q#ixT
zgzu{aWgt(70KK)|dJORyS1^9jKn#Y%c)IVFqZW+3JvvzIf6#mn8%`y7v$^zT`bpE^
zLR}`4Y1Al96fF1Lra=Vi=TXPT*;jw{qi&u0#tXPQj&03YFq*o#3q(Zaf?N@ksXLdR
zQgZ$TA5Vo^%7GLFuJpy7Cg%ZuWnM!(>mXo@9c{xl?o(&Scfb
zA%I@z52OW=s(2bmDG^vgAyn4udmYaLW#hVb$R5a&I{Hc1dCj?B#tWD&W_+kwu)2Qi
zi2y^-4{jb5CESK;CDvJ*oDK+#^6GT?!|WQn5=IhaUYGiLITv}}#u?p8a_yd@H_otp
z@gmB&EHgbw`KjfB^q`D8V@fq&Hui@>DF4>k-hghhED)ahPfc`w8)?S-cG-d{G1s!{V1~T&RL#mcch#9IBNqrs~9%=
zR#UO_ORIvf#?9{PP9u3*?SgvGrfPtfv~+@}@(;kvFqXWjI-5}u@x`gDtft0_eQSj)
zo-1&|^IFIdzQt`VU`258Lb03wdbrm!C6`W1e{Icsc1vFyas030kIzFSgLPZ3dl3h^
zrAaYf*E_qx=*E*CZ1?+zXO;HVogZV2*^>*SoX`5;mwqnmH(j6IS&=QYo;TohXh1zQ
ztZy??_4NJF>1=O6eNsDc#{mc;4L$JYg}Lf@q5Q-ZCw$PG&_XV1K@}6ncvtSDH6TnI
z8tU}#^1LZdC)6N4dDLX{%F^6)qe>P*G^P+XZ5NdH#(2TR#F_#B-t?
zvnKMx?{w@Bhs~yypjzlrfi55esDM{Y7!UV6&x8GYTtDUByIA!V8ckmZ^$04mKalxo
zIaBhlPYCv|!@h^8JVgP6Dv3gVT4QagB(#CY)9mbua06LcAv;wEgd8Kg$hj5ZX*j{R
zEh3{55%%2YR{mwbLtx8t&-g2Ao4)AG1A|qNl-_kXdD1|l%Xy25XA;mv{_L-^7~;`b
zC48jj7@DIYhzn~#Pc>fGc`!qhQWs%W#(1{fu=%ymQz$WgrdTF|FLXIXfrI$ER&gOl
zN!5YE?vR)^179BE;mD~jm8942s49*uV;;Y%YzN_9%O{=U?YOhOFn@u+OZaii^vat_
zw@D>Q+-SX4i+v@VkqoeP`gbzuxggiJqEd$S2luUDtTG8e?3K9wP`R>8m6Ix*%cxTb
z<%*3q;*_N^I3#rBi<-V!YqC0lmVq}`2c^gg`SUqswJ*}Bk)ycU^S}ne&+U_*lwST8
zis|ly461SwiV4I})@e$Nn|8ArJ6?`E)SY6%z3R|)VW!ZQ^($w$6C+~;a
zrntiI`+KEjZ+3t^Lq_2lu*mE?!$*dcV$>f1BTaP)(=k&lL({gUig)kIfM8;wut@P`
zo)#(+vr1)u{^ClU`K+Hp7IV
zb@eP8+XQKeq5?>b#y<~q8!UtF(`J;vk!P8$kGJ1F`U?Fge4^1eCE9
zM|6$uPojC1oFUk!j$4Y5&|u!z9;JFqi7iwLc;~W?;ciQwELMIJI*tJP2G$#TWltqo
zgLbTzLfhR;7^@X<{xf;d#u1`zV2fWfH!}5a`rhvua0@TVb1{D-AN})6?_bL`l
zV$YY1SYEB>&GTM*1E|_Bf|A@A95-CtLk#(=L>uJ%X_DSW)l*g
zNHgDPU~DUlp&6UPvq=i=?0LhQicw@NcV9fD8c~;D)X3{S(WLmP9}IjU)TkH$zvH$S
zMBs%JXM=o}ket6B%Bg#uGPJ5m*(Cfh>*b0Tqw^nHm8g*P
z(+YDOas>qAgQ2{qe1eA{tJ+U<21qA=;RI7g@VlDN
z!vr>W;_=e+jRLlioRs!x1d*{cirmM!gAT&QX135F>!UnN`h71NeIMo{homA|zzCcL
z0FX;31t;d+Ol5`hTkdUd1?16z3@DU
zu?xBq*%%`vQLXb>KCOdGkAF(mPqP`HHgOch;&9DNl`qBxzY%@p5oFsRe~#ab|4W)8
zHL~pEGuko1h&jEK_ROpO01USt?Zhh$JSi()4tIMeH@Q{0DmcdbM#z{rpa
zR4qS1wSr~*(ncqaTfD#_`__A
z`8ZC;O${>P12SA*8LgI4`7+4{e-ck8w?_D63mobQ_1lybI2ISnnhDCoE@tb>v
zq}oM0hnDHtUv>=6cz>R;uf);f)yPT&R^lK*&_MEwHmR=wb1NXBRXOFd(Cs)>5&uH3
zFr@pIK+`A{xyX$|j@UeMxV)%|jzY-J+C8u-FU^Z#OZp*)3!7i*RFD`0`jiF-n<|b{
zT?pkPsogw17hI*Th|vpZ(DnM6(jfKCoog3Rm`!AITtK
ze!H6Tcp=Uv9}iG$pJ%5e2drH5;$W#ltXdji3^T3_e$Qt9W}-nI)~2-*5H69h}HwYvppk`
zsO%jt-bzOwSd%gj7L|$J`K4f1c|8PQ@K8=4>2*6emacw-B#F-hKtiek4+`I0tWNL&
z9$sQ7pcx96+QpR8+xW-mEGQl1D2qS5KBqFQg)`LheqxViOsfTd(PT}b*_RF@(c{4o
zlHv_u1mcZU(a4GI-4T+)y#D>bRuY8}DR^-2#;Sv`9=0{ZNGU+fA%W~tls7!-yhK5u
zpT?f52Q=M2PKR>?z8J7+LIW=Me_jpii!cZSVaIhANDY62U&Yu-MIT6%2z*^PX_RtF
zT@NZLy*4*aRQ+CoWszYdy)N#l^T{j`pnva4W+4Cd!fM(vXipNbsKpr~#6|NTcE)lp
zyjC4?GIED3?3&w!Shc0ra6hm>QxJ<`es%r|za(35`cA6^=+QaV1(Pb1FK2{>!&KB@
zG^fPsFqc9Rc7Ci;2K-E`Ok)E-9>9dHv#}dOb@52H6mC9_LsQWt-oC*sYWm~e9&sko
z&4L1ZPas!ePG5T>xCSYxPhyEPG?N1AY6`$}Gp^NCV^7euWgFXom@p^vg76aU;*5oy+Z;&Ox;1vB4k?8wz#o2}>y~_i8
zRPeg0q03Sxcq=qS%ZwY>@t=v*Y^NB;KpzfNhdUUGN#C
z{-c3>A>381Urhopqv@vc{gxmX+UP6GIP_1O5(Z0oU&bd+271aF2AA;k!;4sJzTJ)0
zRBIktam+ZYbaRkCzWXYwJ}fTji+A(-z&WmGfas&km??mF=tO{KU+@ItR7}-niPJm=lbL
z=l@#sh7-X!W?(@><@nLy&L8Fw6Om@}aS7<8EW1gpUH3PK%R;ksx5dlj0W>8&b^){mRSwW;YxGax3t(5%sS?jo!
zd?!9Yd~0cJ6c9VZl|bD#O!OvP*%D4Ap)m!tP6*S(%1o0JT@-y$8_mW?&3u`P`_)&?
z0cUpPtLUx;4qB*R>TGmSmMySw#at!oBWkf;bY|3J>5`##3?3k=)(f;$mn!Bz`G79B
z`7pNysI+;a1z3IeNu4x;dL4qcYq^jEhPK>8GV_{WF%KPLG?>SAvMLlBAhXt$R3K4t
z&%ugVm1L>4aBW1rn!G?bQ-pu9O|fE_q92s|aK{*x65<^g4upqbF+fIAuK6@8NJdtj
zp)k02d=F&UV3X{`W@DV+qVSvzMF5&BSW?a~PFz*W!BM=DL8p-Fm|zf0n^~zg5k(|#
zC>}}d_2fdE*rcXm995h~dZ}R3!mUgWOAsZpufh*74Gjih$RQ+rW6w+K3!N>>d@unl
zzY$I|tErM|P3@}fS&Ba;b)5=g?acgdP))fIxVXbH!haG|rp~plxPQ0oFaU|aW3A_w
zU-f<)*Yz9L+{t^Qn({U*zusH9R|)+PBw!!5_ir<;yw0v{9h-fc^=>j~=lE;VwBr~v
z-mKl+)Mvaj64<_O-TznLW1r^E#g?EZQ_YyM&w?U@;1tlBPHR^1N(=@>c=wLY*fInR
zx^<*UqbA>QaYU`*WRYcd-LZgqtLGxJ@$SP>-(`Y+*$c*SpFx=g`
z8zXi4!luW!6T24O-8vDc&aa2LoCYK3dM^3kv0tOFXrHP_+wuwyfb4a7c`2t0VBOy8
z;NtxGEvJ>dQnr(wr-p))Cl6saA$1s*Cjo=^=ZRC`mmlny$#3I6jl=fo0`8C6E%T1q
z4~+IZah^W)O2rmiHq5#ylsx=ex#(rrXRynWP&eKi
z8b`rZ%;jl6DrDIC!v_fQcVvzjcse&SH#a^<3<}vv8TZo8{!6K;_ugN{FG5*_7xD}i
zc^e$!8}b7qSZp1drsDwwE{ID+eZk5kYG`9N#nT9xGF`VBKs99lajOq14EmAk`6S{M
z&SUFIEq0&MYuox^A{|JyzA!V_F_V-!4W625JqeZZD|MPdw_J_r{j07v;8%BE_z)2g
zT7IYHJDvc2I*4Sx(}cKVmuK#YsvG$^U=2TERQv+^N7a4RD5E0ICO#>jZ3b2-kncf5
zVKJm#@mfFuuo&Q+6Yr~E0Zn{5$%kXoM$GB)40j!p4dPE%>J_Zq^AA>>56R?4cMgI*
z2b^Wv*ThN?f2PXjp_~sf;Dt?k6=4buCsluu3T1Jh&7Eh3(KmV$Mt&e2jz`Hb>{>IeyYTf|i
zgWyp>ad*jB4P6XeOxK+Rwh1
zexe|}Q$p*YzFJtC+Zf@${rdc4g7gr#A^dpz64P0DH7RuclK>Z)3Wa_$>1Pg=VO>9z
zbS~ZDuCeZR#N%$E?!q+ptv*sD8^@dGA@FpS#tMCt|2z7$d3*-&)kFQ^lkKem;IRRw
zA5`$b7GTqQMW`KKZQVHE>;_dZ5v{TZb$-sxKHaZqGatP)#Nr$#7FyvX$~DGq%Vj*M
zye?UEH@J-scawXq4?1mqy(XEgSD!(|0lffmxyXpZLG)-mWMy!#r1;Z!ih}xanBU{vk{sn-bfi
zjIGy;Ca3?&`6gy;F=4H6xOcpa^wULg
zkr-@2^;p8{JY6
zN3w_WiVA_qpfl@;J1(1u5mp_dEPkW}rlEGky-p*X&1PH_uSI>fbzvZ5$
zW4sRyZ6=u|K9$}JYv6qX;D#AuX${p5dCPZdi#E6|){RbgyhVr3jd^p79Ajg}0wktz
z1uc_&yGWS_OiFN|SK6Fc*s|PRmjNbfHjiDy}ElXYO?2
z$e>l~Ou&9yd;rH|e5map^wOGMIk=p6GWhE>3gkILq8TaDQcy*S?Xo$Zv*@3qaU}h=
z?UYFZ)5L|0*{qdxc|uvknMK{_cBK#SY_N!q#2=$LG7=%=lN3x;P@xJXd4{tlvmWWQ
zZ%P&q%N)Gn&pw?1>`-Y3J{Mknf=t)3PeQMim5N2O`L~}6k$MQJ&2kf3_5#<%_Y_yr(>XY_R)zZgnoBTTulU_HbdE^W_LLy4fH
z%$5zns#dRv+=QP#U07lgQO%9CuKStGic4uegb6v-I7hitzm*jTp
z$byCSD022kOx4D*A^7quyQFJ9BJIjSw%Xw=Y)XzPeU6XW{n0~PugGs%ftf-Msynm}
z?Z;sxBwvf9ICwk;Y^3+;Q66n()o64Cg5VN&;7M%`m2kcK1rWoesarBLPD
z5KyG&c#~1%2vlHZfkqZrbKLy-Mq>s<2%mh8w%jnFf%xMWO;U@tRE>6e)Qp5CXV)@A
zV-3G6Z?A&7mtz8+rMlXQmY&tpq&9BC1ahQBOsW%qY#~j#vs#3Gi%jT4JME;&Ha+SXtx5|YYSPv8R
zpy5jlI2Y0Lx$56!UJPE8MO6!NOw7~K1S?R07!q`6n|8LO@4|q&C8^8?;uE-ajLqZ{
z*@hGXFDZ6jNoG>Iy6(MM&)dF^t;ao^&|Iggnv2lbZb_O%fHofa4g2QknrPgaWGoTF
zgK&9lG#{vvtOM7dnfjnYef3`y&a1gTVA>Ca4muA3s3Wcrdb7s|fNO$qIyg%w3wMvO
znX<^4GSL4HtbL{o+*k>Ujg^Isojrr~5uEV^J0
z7^|ohC!3TQs}#E=kC+&b5Q)J5YcYglRuD*y3&Vj&1AaY3tYrUtYa2||q*K~V!4;s$6p)OKGG
zhiyAxtS#oGdR>tI4nyKZvE+K4GJ&}kO`FA!Cim{u(>L$Q|VtnTP6nTKc%MF
z(dQSp3fNp$DZRfgs`ow~{`Z4#eP2-*-B6oU4{Lr=7Q{@eQ}&~yo#}mEuDsod%$=9U
b*aX-kGrG>;gMc80<6>cfqokBnmV)~qI)42w

delta 10814
zcmai)Q*fV;6R%^dvDNq++qP}nw%*uI8Z@@e#&hJ#kucRa>XTx@bgOINXPVE|5#n}Y%6w}XG@li!I51V1$tK6&$1_+Gu1FYYk}IXwF-lS89yvI#|JOo6`2aUp1NEJOkh84bvL
zMi1q4EHXKRCoADggNfQO6#f|*+d&4!VUDf6lE+6^HwF??R(aUFaMY
zD!U&PT7WQH3w{=U-!Xw0^fb%Bgp2Iuf>FP_?!2JJb)AJ?dCqq8UX}inTt$P7%979Z
zC9^bIj^pH)q*Vx3ep{v*?h;bhVl*(2KS>UdK;KQ(I~DEv@DXawn?c`u|H>u)xsu23
z+>~U{u*0>U1&=*dvk!e;2hX&|!IGNv;f&Oh_I6nXyr=@sCbc6TQp@7v!YY+^^S467
zy}bXV$Gni;ndintscwx6-~Nyj9Dd{A#p^PTZIF1!Nt
zFX>NU8Umcn3v(|AvKN0Ci~IV+bh4y@iKLG?>Vi>ziOqb##GYnJau=Q$8yicl=V-0{~PAUcwK8j5I5jrc9!NE*9F#1?V55|jj*oIU+%B{;N
z!=4@o6T>vqFuNud428Qcgg>l^7fu)R(3`yCe-yX$0KuY7*03tyOz1n~5};NH^J)?I
z7s61`0{`c7G$Rxu|K%kQEXS~ToUWdtx!N<#I)_)acGGbZWWVqO)*&FnEq3E0@>+a8
zdgsm3M{u$FuHy#Rr8+NG81y{+jr|z!9h*FHMF>=pnb*<~PfW+OPF|N7W;q
zR=QWibSFb?G#XLpIYZ%{pnH=Rx{WGQk|llKvBsZHB^TkNHOjT}f`=Evy+=2{b@t12
zGPmKj*`rR;%b(HVTmfcb-)Yk!n7yPFa`>y8#Gvh8{Z7b8VR_E?ZsYzQD{rg1fW&Er
z%gy+-)?=v0`q2I&?&m7Xre-4M3^oWDhis%`qGC
zrTIrVP__+!yL@9dun))0jPqW~!y4Rpw_%$@8~$m9V}57QBY-z&XXq|vZqrKv4=V8B
zZcF8@oOG%Oi;>gFFg9v=MOe`VuanZKjiKl^UYDE4C0cQcmn}54#|Y2;*7>Mcwx2Nmla8!;9X{DVg9EN^EeXut)ja4H6KPg~!KbQL02e?1(^WUkSZdQVS6U|^*
zh6Dk+a)}3{-+Hg>l^rpS7d|Y1N1!H{z)C}mMZ-1(8<1E>vx#Y=QMc_%C~pUlKg$3i{#>O?-|waLiiwvlwkl#+n}H1Ww{N=H
zMR=XZ9S%afPeb!p*=nni1ChtT@O5$tc@3-Wfl*1#WM%BK!{AOt>$^KW
zFV6%Cf^GXC8w*-GXKQs@+pZ1bo$inDTF-y~R(_ik7Q~S3ktA$5`%azD)AQ
zRiS_&$2Rlp^XCpmo-~M&&RxUSUk1Sbpc?v{gKl!xWa`2DH^9p`=Nnp@4$jQyR-F%v
zr8&Qm8XdY)=iw&vm?d1sWSXbQGar?Erl}=tom}kw&W)->WApT2m3o=s3{1eU9cMKI
zy|fnTm5D0;KHDEU4~V|FUlscYa=Sd6&UI}?bdBTwedT*_!}e)(JlU9cL}mcJ7_YjX
zyxBj{Nl^zKN7l0bd|!7FscggIT6h(wikj7HaK2|(;j$*vl72t6^Z7~MD_OblS_Dfr
zcWZbaV$KJ<@ibjzg^(-{_SipHP~5Mm#|_OpGj^)#NG)}$p(bIPij1%&F8G(e%0)!S
zPbsEJJSGsO;7n{&Q=8nSX%&!0Zhb5iVZc|zRaDo5Q!1&v$SZ#MfhS1s?>olsHkp&X
z*Mje)SKf3b*DOK$Uos;<*vedY*NHxw`!0~BF7vuP`I4sWN4qUfbxG-kV({j}5eU;&vmmm;g%TS2s>_r^Ugx
zoGJ?gCQ5!`re#2=n8XIK@h24sNX#s>Dy3Y
zZ+)}XoS)rQke+rJvrmmg8>~0cWQ4cMh>5l8Izw6xC5#{f_m^^y+L%JV7mNz`h5X2Q
zZ8^hkq1OzO-j@xopgykC-PZP_;lip@b&mJQS^*QG+FIwDB>>Af$*h-i#K%6M_||Cj
z>IN$P*~Iv*kcajA6ZZns5P-tk|yumQ+~5*+{gr{YnK=z
z#C$ii09bk7kX=@Ok$QGwfjUG6tIc^oPLM|K89)`0tbcRrT}RQVwJ($WS^s=;^n4{r
z^p$A>?5`fVBzNm6KWIlk1PLTt6@xy8F27@pt8(nS(IpJ{=$?S{Ww>iCC=DX=X}x4!
z=lbu$xv#wnYKvrU6EOuo9MjT;0R3kt~i|9J9Tv`I)l8S$q_{J*lrZ4^6N(RnU9
zn2#E0@uE_)_hxNN>tGC`en^IyzzdAKD*Xzjc|**wM#}a-NG2ki_uRb72@hp2%Y379
z+oTUhXGVr}O{e>8nEi0Xtziub4!Cb1fZ^ss>augUYdcLMoQa7jHK=OEQgzsDbV=Bs
z!Ds+gMkdor+(B-ljMc%5P{{#S{p;_p+V{$(j$&!x5#pLc*Ep96{L%mVnTa_bJ~+og};sHe-D=cS+FHDmD~K%br(
zUto)6ISX?5#u|_9ky0$VWTp1n|8$5)D$sr58;sDuwA(aHRWn;*#0wH%+{c4pc@+&P
zR`qW%&`Bb)94#|WCd&;=jFV*5m-93e)tRI9-w%^*Q^Wz3LE^bSQH;D%RW#cyvj#xCxjIaK~U+C;rn7KosHziaLXA?!?FX$1~iryBnC++#Y)7uM!Ne1&-uD
zDuJvp`4ax#^?S?xygWuO@kN@93y-YiCN#&TCEP7kuClp_g)A;Ho_m;P+lKch11<-t
zI1ruwgAoyYg8sJzEYw82MAW%Qc
zz6<|@PpXc>^-G1p<(3pq=O9GW)J*~%+i}ckL>lA$*pq(K;z?NqZJ4&&`dw%iBY8~1
z?g$eDm@WjZv}hf3ly*^F0D)kjP;;R$a9q=n8m%Vf68&IsL08zXlWlOi-vZ8^t%s+B
z{Ey4Q5=KLiVI=Opsf!o?InP9qOOuxC**(IAF&_Y@R<_S
z6xey+TxLxRc#*<~ndY89V0C*s6R$o(O_f-Z^o@ePo$)8O3?Q4}2@A}qo#ho|&Uo|b
zXbB&>|H_cRhO1b<(anewaKRopkir|+X%o$>hTneyZEcQ(+j?D`(fj7uzx
z5H|eMPxwM`O%u-fNsUF+JvjNCi_wJ*rR0Ad#Xbo__U#cXuH2C;va?q*)^>xD_>Te9
zqDUDQIoBC(q}^0!D^B=TzRfyx-u6O$>(q>(T;KSDj_YWgP%aVu0`HQEO=a=JNPW@9
zg5=_1`){#93*b+rY>qw(8U=6l$E=Au345vXQZ)UD0sR>HSw&&-0+yPIH%!_ORv>vG
z)+=!Rcs?rW%$@v0ozc8FW7TOlh@vww<ObQ|2e?$`rh*^FhfhC>*ez)d-d9zQ(y$W4gfbrzm;~LK&O?2gAH=_s`3h0
zU&^k<)X`soG7EI$&Znbd$j=A*Kmh&QB}q(oZT%KgxTWga%ZrwGW3*c*5ECCwLNd7x
z^GR~rSka}{^~Ysl!R;&QnY}Y~eKa`{hYI{pu!Ss}$!RY+iLfl*E|w2q8@nn6Hjn-0
zx(bO4wBJQ_A>ey`*q{>RomK)jzat^2?tc;sI<_R-pN0KfkY8W;QsQ3bd;r@Y48|?&
zA7G!|WUcpV4kdRy3z5bj|3-X+2IAA)TM2mu5`c8DK0J>qn1|@U@x~7-x{jP&VzZD=<6t*LJjj9tkIDb
zbEaY&>}BXwTlT7e+8Q>H;ombg7s|Npeb>JFct_>_62EsRkcU8Q?iZB`PS|ybOUR`+
z2{1Z-R>;+Iw_mjnGib_%Wt`6|qC*j;x;XDD&@>f|gaAp?zx&Ygr(=mKc^E5Ifh{ut
z*6t#|N#!V5)hd=WST0IMwwE^HQ=D<1yjnOz^1#&dXx|zj5S4b^kG9?+HV|rU^w&20OuwxqL$u8xY^)72toh&;0Q|Q=y`wDRDUiYThX-vc<8WvykaJ
zWziZo-X5fI-hp#H9Y~r9Jn0n&n3zjE$4eOFJH-J%-yc-*X$1Uu$UX#JbhQ1sp
zkRnB*apO8yrT*K55Y?6VZE^j@w*eTX3{|g6AR;8&*4H001-6>vgj#h6(
zQpkL0D$`i9**M#cb}Tof(w=hxDdwPP6zcMlrmz+bih2-M?Ur^@`pN%|O3&m&-rF7J%8MF<4~LX69|HmHz7
zY`~suplN4RK?6=Mj5fa+!dtBoIe(VAv6oAO!6niE!+nup&ZOIY^?;|jwX(_ANQ#eU
zDh>c&$ygM832#yQL~Ugp^m-xq);N>m1PBKv?q}Rsm
z^NLswK~fonq6k5S=_xDF5}^ZPcx1orQQi^6mMtP6Fla6+bvmg|rq%tEw%c*O*O=Q#
zzb0tmQm8_%&DaQ!;Q6p?3KpV)1XUDWvpABp*&J~Eh3Zb0A?a8F$iNVkIg}m|7xN6)
z)!Uj0T<>Q}yd34{%3OCWqCV)xN|uYF6;XZ>OF;%6H+FB$pEsRQ10fKD@yq=$?Lv71kISmVPHgd
z!Kiy^8W%RvVI3AQF9<0Jmh1n3k+%7E6tH&q>`rODh7qq6oKpm4yP14N+T>2#j5JMp
zbl)p)_^WCsPF}}u#+k;dqsa}81Dm*f=-XEwBQp(A4urL!`Vkx05@KWFf-5q7@<5JG
zKeuUURhMhmA%a@k`hb?Z*UrL>mx4wz6tC{Tsv8Tvo@WmLwLl9t4Dn~h?eO2AKB4)m
zC)I&s)fw@9k0A{z*)LebA6B&Hx&A0+fxYfWlEUFV+@bCuzI1O}Yw-RkUIylycnsg0
zJCzMTB-{%+D+R|^_%44Z*}skCYN;{uv_fjR2&g_r#xnXiZ_BcBHWQ*Gbkv``-;x2}
zI^QRhL(nd8PXoA?k%#Evf~*PiwdT3!A|vAY*5-kGvHXljh_yzBa|55@9fScsdogTO
zpP$QO)9b3MGkA)-1{J?w&`RPxWkUAR_Q&W1+%!5tD*>o>6^nmf>~buR
z6&7&ydXb4TAcw+#{)>X4JrhzM-C#Vdg*y
zWFf$u5HfPZg{aW!EtBUwhxoJ+QDRx)@9!bSB0fH&<^3>`1Rugk+lZByx2=#~bPq)K
z{+!l*cqgB~;1@zsB+{9PJ1u1Rh=YCNBx5)f9eq|f`(7<;d>{atd;|4@@ppvN@bg`>
zi~9jOXGrGo6?SN6*M7P;dEn_AMMA9{Gb#K16kEDW#c>sU+Y3}@;zRvTu9n&1ErYd4
zJD;j@Ea4^k^wz!PVn1ANZfVU0!S2FE)%WyTb#5YRJ1L~iNfUf{zHF>36Sw`bgbs@n
z@p)f!jdI#Pyb9Z|Ar%nQeYCqFZ%p9rcVGcB?WryjzEv0eonv~wni!pV33^Tju0oo<
z_5EoHjfJE)Q447F2>t2E_=4#rC%yA=6?)JHQKfkP?%alRYHq;qlio*87Dp0+hJ#3f
z;qSddHwq%VK%UsZ_=4w$zC5@G^u%Ejx%^hPoM*<>C6A$FrScGG9@gu<54;I&Dk5Mg
z>K~hk)qAA85s8t*i#e;)uWn4)D&-Q7z^d=|6t1=hK7YKWzg&>MR})DE!}&`GfZI1}
zv^)FYz2>gy^r5}(#9?C!z`HMl?;VrYzJc7>ysPyC54q)8Lr+L#pq%~!`x0dn&vz&x}w0{S#w~e6suV8Ega(=ZBLSMa1<0&d|&nW-w
zJh)V#pSQD=_p{KZFa&~HNx%TLh;O#bH~4vhJm2~UY;>d>kYJZ^5EgPo3G|kM5g{Lb
zf?!bOdw8TH42l8Q5VHQ_KG*xprIR)D%ctnw-gC#0pi}H}T&&&@)HB}dn!O+WTcn#7
zk0(dN=&9I75Q)H2*O@p|%l2>Gs@xFK1Wuo$sZ>meAB16Zv-lQy>%dk?!I5W_!5@s6
z!k*B?!s@8^^8kLYd>(e#zM&)Kx=q|DtRX2VMg*pkP%Wf{&&?B(CY@kujD8OcDgY|=rj`kpi6+fW(
z<|IU|j7F=dN!KUi2gjsUagD}+?+=YA88Dv14AFcVyRL=et|YSM
zUc^#~8HLizEOC>xY7J>%IYm}P8qScv5&UR`VLnbWZS0xt=&d_eO-n-LvH4KPE%Z5T1+cK)rblt6sj>j`o^zB
z$F!s}<9YR(cFwSM-T<_3dhb8wCez67$1O}6=0Uwk&y*m|%1cP9)
zYe^9h6TydMV$__+;EG^PtJi&-K}Y{qvzHSO?ycs`gKu9Rk1c!fjS>Oik>l&@X2{ab
z@2lKJzb@iy=iS>)jx6GFdo#Ij^kf$Rnsq0)fu!V)Q)2cCXtx#86GT4?kVy5~v$d|B_!>h@a$-YZU?i3`~_WiEZA3|T!
z{(++G42F|$>hFcNtIs?4k8zIkGpz-$FN`c(r-X4pa-2P1^t$IB!edh~CDjX-Gjjv}
zM=GimA)yPfITlLlT=wxE{PCXkbbo(&l86i4N}2QA?iQig*mLJs{D6>;dz2=kB)+#J
z+;7yN#G(xM27_?L4+$-_DXg}5Y&<@=K2z*rh)S8R(?UgV`*x#;7ifB+wzEvo1UW2I
z&zGE%PpDA5&E>gtg2(ss1uKI`)zt@=@>r*3pp*=CFvUen|n-~{9wMl0dAUtmyd5UB@)J9CAr2bG>8pc
zz=J6TU!^YP_DDTOoFcoQkuFMY6;nSaqR}$PufS?fkeL41cAGTKe<)bId4YG;DVWk0
zjAaso5hX`Dmw&mT{vot&?|7a5{0zyP3Rgf3Ku)P1hgK6-M-QPT@;C8Io~z;dyEa=p
z400SL;>+;vr!_!
zk&pVoAB$=yyWK%1htk0{`7Mc4`AU~?0~H|Z9@=n@X-)aqrFMr$u1el-WqGVnFRlMx5mVPL5EEj
zJNc+?*PX=4AfMrPsG2an2xe~P#mzdPniyPn`8fG23?!_hJY&g{P*(L*7W!@-SbWzA
zN13fhXhH^C9Tsk(7m2N*P+Cb=m>vf6+}8+$HI~?-q?5)2(YU*iQ~n`H9lH#Je~%0L
z9uQklZWronw#S9$-eAlz!kQ3UKm(%=c^x4w1EH>Zh-IJVRaT=VO~L_@D&LP$GpKfC>g
zj5J{4wk?y4oyN7Us7I1F!JKcrJH`4&BQ+3Wjrb;03A@`Y&pV9jux>rw35H_!4P_--}|Ys1(flkiggkw2;80*(~}
zisH$?J3dke^)43XF=c_u>2bkP2ju3{*59PNU^>HL|0M~KZ&Ok!
z9F~8@nt9gwd#gFTjqm>10^rG)@(6YsIHAVAnFk;V%G#<~r%{s&O%K|0Ka7cpZzQ;s1F3MVuL;O_Fv@}Z=r@ooY
zOCYBFn=FD^evA#*nsTm+P6&opVss90Tt3p>L%a`pR(}w`{svP%0j~A+RnbjWtIrE|
zQv%Dc1Pe73>ty{XGMVLDC|cR(hU#M4AdoAHD~LZ@B|dGZk~*Gzs|=Us2c)h%Hu<<-
zAM@6%mdc1Fg$)#Fij!5BI^cL49~e#ybr#V^VcJf
zhE}Th!Y{=Ir6C4u046-!zUI6*MSEK>epe%ZL1irjYUm=+U_E|G*JjcZ773sHLm{N`
za3yLXIuT<2zEqC9OREvajoQX!dE2_A3~5(7?9Z*&N+=mu`7ih%p-v-ol)*lwhJ}Kn
z#qKHoA*k}sH(nH=UQEdlvj9~ae5*h6uj|g~6v3UcdRn?^fP^J08l~IQr1eJ(dSO;Y
zm1Tx5v}XV{S$Vn5WZi*Ev4iX+W~l}zQ9$YK_&c0*ne=&EbrjZ-i@wANlv1K(Dv>Qw
zRo6S<@j3MV%;)u)3+7-t-)D7F9&6PZd+#2#Q7`K&tM{_R`GcIX)h@vDaYr>`gv*
zHcu)$ExF?eFie9KU1^mh&6Rs3G&bk7(+6o(Y~NMhzhZ?#VNKg@fpS~8xhKhb2x11i
zwZzV~%n)Y8WsEQEXuk^(Q|8=mZa2R`n`JP<~q+
z^fpikG1~r)^3QwiT#85TSou(RTsRq7Lp>*h`{xK6;L1Jj(LG2%845!TLX1z_R_=46
zicaTcKzC>%!#ON-!{l0YrIRzv;SwSQ#VZ>AyU!sgg*tHjI)n0*-EXec@O*sR(q*OaX*#Ar@6f##HlDecNb3J~%8?QK1#V3n`t?72
zcXQ8>dW1j@r7nx4E(8D1Alj+RAXOD%SXg*CI5{$a2S}p-r`wAjAj6SW#W=Y|B}6&I
zctkjvStPk6*xA`d#Qq~-HVz(URt^?UK~nzzuZsZ!le~qamAf@5J2wyM|Na93qQx#H
zbrH<5+xxC&NcwsSWtDH362dpZ`)5`B{-7v?|867_H>d@pxS%Q8nz+(vJsn&ObUrit
z!mzMeZEbC9)1{TMXZHN}j^xz=DsYQZkW|SFqs!4MafWjhZdSQR^YB*GT9~=
z?t`3b(0V3_Z{+I2l4#aQ+5kj5a&iVE-o@M{$sr8`)})H78Xd;6!hB2pTt7r322)IlAqX!7v=u
zWR*oi5DA|dmVIG{=@N4o+6mpR!`kExTc-1>emy9V9$)gt42f!teO&Q#uzt>;)kR5-
z3RM5XI-Z{IHz`v}omT{?e_mCJ>JmPv)7B>VRi*UuKc#7Le(k8WHm&C7*`|FY5

diff --git a/doc/src/DimRed/DimRed.do.txt b/doc/src/DimRed/DimRed.do.txt
index 17ddc720e..0ad06f689 100644
--- a/doc/src/DimRed/DimRed.do.txt
+++ b/doc/src/DimRed/DimRed.do.txt
@@ -786,13 +786,8 @@ n = 100
 mean = (-1, 2)
 cov = [[4, 2], [2, 2]]
 X = np.random.multivariate_normal(mean, cov, n)
-# Print the X-matrix
-print(X)
 !ec
 
-
-Try to add to this code your own calculation of the covariance matrix.
-
 Now we are going to implement the PCA algorithm. We will break it down into various substeps.
 
 === Compute the sample mean and center the data ===
@@ -812,27 +807,28 @@ and the mean-centered data $\bar{X} = \{ \bar{x}_1, \ldots, \bar{x}_n \}$ takes
 When you are done with these steps, print out $\mu_n$ to verify it is
 close to $\mu$ and plot your mean centered data to verify it is
 centered at the origin! Compare your code with the functionality from _Scikit-Learn_ discussed above.
-The following code elements perform these operations using _pandas_ or using our own functionality for doing so. The latter, using _numpy_ is rather simply through the _mean()_ function. 
+The following code elements perform these operations using _pandas_ or using our own functionality for doing so. The latter, using _numpy_ is rather simple through the _mean()_ function. 
 !bc pycod
 df = pd.DataFrame(X)
 # Pandas does the centering for us
 df = df -df.mean()
-display(df)
-
 # we center it ourselves
 X_centered = X - X.mean(axis=0)
-# test that we get the same as Pandas
-print(X_centered-df)
 !ec
-Alternatively, you could also have used the functions we discussed earlier for scaling the data set. 
-That is, we could have used the _StandardScaler_ function in _Scikit-Learn_, a function which 
-ensures that for each feature/predictor we study the mean value is
-zero and the variance is one (every column in the design/feature
-matrix).  You would then not get the same results, since we divide by the variance. For diagonal covariance matrix elements will then be one, while the non-diagonal ones will be divided by $2\sqrt{2}$ for our specific case.
+
+Alternatively, we could use the functions we discussed
+earlier for scaling the data set.  That is, we could have used the
+_StandardScaler_ function in _Scikit-Learn_, a function which ensures
+that for each feature/predictor we study the mean value is zero and
+the variance is one (every column in the design/feature matrix).  You
+would then not get the same results, since we divide by the
+variance. The diagonal covariance matrix elements will then be one,
+while the non-diagonal ones need to be divided by $2\sqrt{2}$ for our
+specific case.
 
 === Compute the sample covariance ===
 
-Now we are going to use the mean centered data to compute the sample covariance of the data. Recall it is given by:
+Now we are going to use the mean centered data to compute the sample covariance of the data. 
 !bt
 \begin{equation*}
 \Sigma_n = \frac{1}{n-1} \sum_{i=1}^n \bar{x}_i^T \bar{x}_i = \frac{1}{n-1} \sum_{i=1}^n (x_i - \mu_n)^T (x_i - \mu_n)
@@ -841,10 +837,26 @@ Now we are going to use the mean centered data to compute the sample covariance
 where the data points $x_i \in \mathbb{R}^p$ (here in this example $p = 2$) are column vectors and $x^T$ is the transpose of $x$.
 We can write our own code or simply use either the functionaly of _numpy_ or that of _pandas_, as follows
 !bc pycod
-print(np.cov(X.T))
 print(df.cov())
 print(np.cov(X_centered.T))
 !ec
+Note that the way we define the covariance matrix here has a factor $n-1$ instead of $n$.
+Our own code here is not very elegant and asks for improvements. 
+!bc pycod
+# extract the relevant columns from the centered design matrix of dim n x 2                                                      x = X_centered[:,[0]]
+y = X_centered[:,[1]]
+Cov = np.zeros((2,2))
+Cov[0,1] = np.sum(x.T@y)/(n-1.0)
+Cov[0,0] = np.sum(x.T@x)/(n-1.0)
+Cov[1,1] = np.sum(y.T@y)/(n-1.0)
+Cov[1,0]= Cov[0,1]
+print("Centered covariance using own code")
+print(Cov)
+plt.plot(x, y, 'x')
+plt.axis('equal')
+plt.show()
+!ec
+
 Depending on the number of points $n$, we will get results that are close to the covariance values defined above.
 
 
diff --git a/doc/src/DimRed/Programs/PCAsimple.py b/doc/src/DimRed/Programs/PCAsimple.py
index bfbccb78e..9e7716348 100644
--- a/doc/src/DimRed/Programs/PCAsimple.py
+++ b/doc/src/DimRed/Programs/PCAsimple.py
@@ -2,28 +2,47 @@ import numpy as np
 import pandas as pd
 from IPython.display import display
 import matplotlib.pyplot as plt
-from sklearn.model_selection import  train_test_split 
-from sklearn.preprocessing import MinMaxScaler, StandardScaler
 
-
-n = 10000
+n = 1000
 mean = (-1, 2)
 cov = [[4, 2], [2, 2]]
-print(cov)
 X = np.random.multivariate_normal(mean, cov, n)
-print(np.cov(X.T))
 
 df = pd.DataFrame(X)
 # Pandas does the centering for us
 df = df -df.mean()
-correlation_matrix = df.cov()
-print(correlation_matrix)
+print("Centered covariance with Pandas")
+covarianceX = df.cov()
+
+print(covarianceX)
 
 # we center it ourselves
 X_centered = X - X.mean(axis=0)
+print("Centered covariance using numpy")
 print(np.cov(X_centered.T))
-#print("test that we get the same as Pandas")
-#print(X_centered-X_train_scaled)
+# extract the relevant columns from the centered design matrix
+x = X_centered[:,[0]]
+y = X_centered[:,[1]]
+Cov = np.zeros((2,2))
+cov_xy = np.sum(x.T@y)/(n-1.0)
+cov_xx = np.sum(x.T@x)/(n-1.0)
+cov_yy = np.sum(y.T@y)/(n-1.0)
+
+Cov[0,0]= cov_xx
+Cov[1,1]= cov_yy
+Cov[0,1]= cov_xy
+Cov[1,0]= Cov[0,1]
+print("Centered covariance using own code")
+print(Cov)
+
+plt.plot(x, y, 'x')
+plt.axis('equal')
+plt.show()
+
+
+
+
+"""
 #Now we do an SVD
 U, s, V = np.linalg.svd(X_centered)
 c1 = V.T[:, 0]
@@ -38,6 +57,6 @@ print("Check that we get the same")
 print(X2D-X2Dsl)
 
 print(pca.components_.T[:, 0])
-
+"""