diff --git a/doc/src/Clustering/.Clustering-reveal_html_file_collection b/doc/src/Clustering/.Clustering-reveal_html_file_collection
new file mode 100644
index 000000000..146687edb
--- /dev/null
+++ b/doc/src/Clustering/.Clustering-reveal_html_file_collection
@@ -0,0 +1,2 @@
+Clustering-reveal.html
+clustering_example_images/simple_clustering.gif
diff --git a/doc/src/Clustering/Clustering.dlog b/doc/src/Clustering/Clustering.dlog
new file mode 100644
index 000000000..147687a4f
--- /dev/null
+++ b/doc/src/Clustering/Clustering.dlog
@@ -0,0 +1,96 @@
+Translating doconce text in Clustering.do.txt to html
+*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax)
+*** error: figure file "clustering_example_images/some_image.jpg" does not exist!
+Translating doconce text in Clustering.do.txt to html
+*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax)
+*** error: syntax error in table!
+ missing three horizontal rules and heading
+ in the right places
+
+ lines surrounding the table:
+
+
+- For a given cluster assignment \( C \), and \( k \) cluster means
+ \( \{m_1, \cdots, m_k\} \). We minimize the total cluster variance with respect to
+ the cluster means \( \{m_k\} \) yielding the means of the currently assigned
+ clusters.
+
- Given a current set of \( k \) means \( \{m_k\} \) the total cluster variance is
+ minimized by assigning each observation to the closest (current) cluster mean.
+ That is $$C(i) = \underset{1\leq k\leq K}{\mathrm{argmin}}
+ ||\boldsymbol{x_i} - \boldsymbol{m_k}||^2$$
+
- Steps 1 and 2 are repeated until the assignments do not change.
+
+
+As previously stated the above formulation can be a bit difficult to understand,
+at least the first time, due to the dense notation used. But all in all the
+
+ here are the recorded table rows:
+NOTE: do not use pipes in horizontal rule of this type:
+(write instead |-\boldsymbol{x_i} - \boldsymbol{m_k}-|)
+| | \boldsymbol{x_i} - \boldsymbol{m_k} | |
+
+possible trouble:
+1. Not a table, just an opening pipe symbol at the beginning of the line?
+2. Something wrong with the syntax in a preceding table?
+Translating doconce text in Clustering.do.txt to html
+*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax)
+*** error: syntax error in table!
+ missing three horizontal rules and heading
+ in the right places
+
+ lines surrounding the table:
+
+
+- For a given cluster assignment \( C \), and \( k \) cluster means
+ \( \{m_1, \cdots, m_k\} \). We minimize the total cluster variance with respect to
+ the cluster means \( \{m_k\} \) yielding the means of the currently assigned
+ clusters.
+
- Given a current set of \( k \) means \( \{m_k\} \) the total cluster variance is
+ minimized by assigning each observation to the closest (current) cluster mean.
+ That is $$C(i) = \underset{1\leq k\leq K}{\mathrm{argmin}}
+ ||\boldsymbol{x_i} - \boldsymbol{m_k}||^2$$
+
- Steps 1 and 2 are repeated until the assignments do not change.
+
+
+As previously stated the above formulation can be a bit difficult to understand,
+at least the first time, due to the dense notation used. But all in all the
+
+ here are the recorded table rows:
+NOTE: do not use pipes in horizontal rule of this type:
+(write instead |-\boldsymbol{x_i} - \boldsymbol{m_k}-|)
+| | \boldsymbol{x_i} - \boldsymbol{m_k} | |
+
+possible trouble:
+1. Not a table, just an opening pipe symbol at the beginning of the line?
+2. Something wrong with the syntax in a preceding table?
+Translating doconce text in Clustering.do.txt to html
+*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax)
+*** error: syntax error in table!
+ missing three horizontal rules and heading
+ in the right places
+
+ lines surrounding the table:
+
+
+- For a given cluster assignment \( C \), and \( k \) cluster means
+ \( \left{m_1, \cdots, m_k\right} \). We minimize the total cluster variance with respect to
+ the cluster means \( \{m_k\} \) yielding the means of the currently assigned
+ clusters.
+
- Given a current set of \( k \) means \( \{m_k\} \) the total cluster variance is
+ minimized by assigning each observation to the closest (current) cluster mean.
+ That is $$C(i) = \underset{1\leq k\leq K}{\mathrm{argmin}}
+ ||\boldsymbol{x_i} - \boldsymbol{m_k}||^2$$
+
- Steps 1 and 2 are repeated until the assignments do not change.
+
+
+As previously stated the above formulation can be a bit difficult to understand,
+at least the first time, due to the dense notation used. But all in all the
+
+ here are the recorded table rows:
+NOTE: do not use pipes in horizontal rule of this type:
+(write instead |-\boldsymbol{x_i} - \boldsymbol{m_k}-|)
+| | \boldsymbol{x_i} - \boldsymbol{m_k} | |
+
+possible trouble:
+1. Not a table, just an opening pipe symbol at the beginning of the line?
+2. Something wrong with the syntax in a preceding table?
diff --git a/doc/src/How2ReadData/Git.dlog b/doc/src/How2ReadData/Git.dlog
new file mode 100644
index 000000000..ae49b77d9
--- /dev/null
+++ b/doc/src/How2ReadData/Git.dlog
@@ -0,0 +1,4 @@
+Translating doconce text in Git.do.txt to ipynb
+Failed to remove ans_at_end environment
+Failed to remove sol_at_end environment
+output in Git.ipynb
diff --git a/doc/src/How2ReadData/Git.do.txt~ b/doc/src/How2ReadData/Git.do.txt~
new file mode 100644
index 000000000..08f6118ad
--- /dev/null
+++ b/doc/src/How2ReadData/Git.do.txt~
@@ -0,0 +1,505 @@
+
+===== Simple linear regression model using _scikit-learn_ =====
+
+We start with perhaps our simplest possible example, using _Scikit-Learn_ to perform linear regression analysis on a data set produced by us.
+
+What follows is a simple Python code where we have defined a function
+$y$ in terms of the variable $x$. Both are defined as vectors with $100$ entries.
+The numbers in the vector $\hat{x}$ are given
+by random numbers generated with a uniform distribution with entries
+$x_i \in [0,1]$ (more about probability distribution functions
+later). These values are then used to define a function $y(x)$
+(tabulated again as a vector) with a linear dependence on $x$ plus a
+random noise added via the normal distribution.
+
+
+The Numpy functions are imported used the _import numpy as np_
+statement and the random number generator for the uniform distribution
+is called using the function _np.random.rand()_, where we specificy
+that we want $100$ random variables. Using Numpy we define
+automatically an array with the specified number of elements, $100$ in
+our case. With the Numpy function _randn()_ we can compute random
+numbers with the normal distribution (mean value $\mu$ equal to zero and
+variance $\sigma^2$ set to one) and produce the values of $y$ assuming a linear
+dependence as function of $x$
+
+!bt
+\[
+y = 2x+N(0,1),
+\]
+!et
+
+where $N(0,1)$ represents random numbers generated by the normal
+distribution. From _Scikit-Learn_ we import then the
+_LinearRegression_ functionality and make a prediction $\tilde{y} =
+\alpha + \beta x$ using the function _fit(x,y)_. We call the set of
+data $(\hat{x},\hat{y})$ for our training data. The Python package
+_scikit-learn_ has also a functionality which extracts the above
+fitting parameters $\alpha$ and $\beta$ (see below). Later we will
+distinguish between training data and test data.
+
+For plotting we use the Python package
+"matplotlib":"https://matplotlib.org/" which produces publication
+quality figures. Feel free to explore the extensive
+"gallery":"https://matplotlib.org/gallery/index.html" of examples. In
+this example we plot our original values of $x$ and $y$ as well as the
+prediction _ypredict_ ($\tilde{y}$), which attempts at fitting our
+data with a straight line.
+
+The Python code follows here.
+!bc pycod
+# Importing various packages
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import LinearRegression
+
+x = np.random.rand(100,1)
+y = 2*x+np.random.randn(100,1)
+linreg = LinearRegression()
+linreg.fit(x,y)
+xnew = np.array([[0],[1]])
+ypredict = linreg.predict(xnew)
+
+plt.plot(xnew, ypredict, "r-")
+plt.plot(x, y ,'ro')
+plt.axis([0,1.0,0, 5.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Simple Linear Regression')
+plt.show()
+!ec
+
+This example serves several aims. It allows us to demonstrate several
+aspects of data analysis and later machine learning algorithms. The
+immediate visualization shows that our linear fit is not
+impressive. It goes through the data points, but there are many
+outliers which are not reproduced by our linear regression. We could
+now play around with this small program and change for example the
+factor in front of $x$ and the normal distribution. Try to change the
+function $y$ to
+
+!bt
+\[
+y = 10x+0.01 \times N(0,1),
+\]
+!et
+
+where $x$ is defined as before. Does the fit look better? Indeed, by
+reducing the role of the noise given by the normal distribution we see immediately that
+our linear prediction seemingly reproduces better the training
+set. However, this testing 'by the eye' is obviouly not satisfactory in the
+long run. Here we have only defined the training data and our model, and
+have not discussed a more rigorous approach to the _cost_ function.
+
+We need more rigorous criteria in defining whether we have succeeded or
+not in modeling our training data. You will be surprised to see that
+many scientists seldomly venture beyond this 'by the eye' approach. A
+standard approach for the *cost* function is the so-called $\chi^2$
+function (a variant of the mean-squared error (MSE))
+
+!bt
+\[ \chi^2 = \frac{1}{n}
+\sum_{i=0}^{n-1}\frac{(y_i-\tilde{y}_i)^2}{\sigma_i^2},
+\]
+!et
+
+where $\sigma_i^2$ is the variance (to be defined later) of the entry
+$y_i$. We may not know the explicit value of $\sigma_i^2$, it serves
+however the aim of scaling the equations and make the cost function
+dimensionless.
+
+Minimizing the cost function is a central aspect of
+our discussions to come. Finding its minima as function of the model
+parameters ($\alpha$ and $\beta$ in our case) will be a recurring
+theme in these series of lectures. Essentially all machine learning
+algorithms we will discuss center around the minimization of the
+chosen cost function. This depends in turn on our specific
+model for describing the data, a typical situation in supervised
+learning. Automatizing the search for the minima of the cost function is a
+central ingredient in all algorithms. Typical methods which are
+employed are various variants of _gradient_ methods. These will be
+discussed in more detail later. Again, you'll be surprised to hear that
+many practitioners minimize the above function ''by the eye', popularly dubbed as
+'chi by the eye'. That is, change a parameter and see (visually and numerically) that
+the $\chi^2$ function becomes smaller.
+
+There are many ways to define the cost function. A simpler approach is to look at the relative difference between the training data and the predicted data, that is we define
+the relative error (why would we prefer the MSE instead of the relative error?) as
+
+!bt
+\[
+\epsilon_{\mathrm{relative}}= \frac{\vert \hat{y} -\hat{\tilde{y}}\vert}{\vert \hat{y}\vert}.
+\]
+!et
+
+The squared cost function results in an arithmetic mean-unbiased
+estimator, and the absolute-value cost function results in a
+median-unbiased estimator (in the one-dimensional case, and a
+geometric median-unbiased estimator for the multi-dimensional
+case). The squared cost function has the disadvantage that it has the tendency
+to be dominated by outliers.
+
+We can modify easily the above Python code and plot the relative error instead
+!bc pycod
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import LinearRegression
+
+x = np.random.rand(100,1)
+y = 5*x+0.01*np.random.randn(100,1)
+linreg = LinearRegression()
+linreg.fit(x,y)
+ypredict = linreg.predict(x)
+
+plt.plot(x, np.abs(ypredict-y)/abs(y), "ro")
+plt.axis([0,1.0,0.0, 0.5])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$\epsilon_{\mathrm{relative}}$')
+plt.title(r'Relative error')
+plt.show()
+!ec
+
+Depending on the parameter in front of the normal distribution, we may
+have a small or larger relative error. Try to play around with
+different training data sets and study (graphically) the value of the
+relative error.
+
+As mentioned above, _Scikit-Learn_ has an impressive functionality.
+We can for example extract the values of $\alpha$ and $\beta$ and
+their error estimates, or the variance and standard deviation and many
+other properties from the statistical data analysis.
+
+Here we show an
+example of the functionality of _Scikit-Learn_.
+!bc pycod
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import LinearRegression
+from sklearn.metrics import mean_squared_error, r2_score, mean_squared_log_error, mean_absolute_error
+
+x = np.random.rand(100,1)
+y = 2.0+ 5*x+0.5*np.random.randn(100,1)
+linreg = LinearRegression()
+linreg.fit(x,y)
+ypredict = linreg.predict(x)
+print('The intercept alpha: \n', linreg.intercept_)
+print('Coefficient beta : \n', linreg.coef_)
+# The mean squared error
+print("Mean squared error: %.2f" % mean_squared_error(y, ypredict))
+# Explained variance score: 1 is perfect prediction
+print('Variance score: %.2f' % r2_score(y, ypredict))
+# Mean squared log error
+print('Mean squared log error: %.2f' % mean_squared_log_error(y, ypredict) )
+# Mean absolute error
+print('Mean absolute error: %.2f' % mean_absolute_error(y, ypredict))
+plt.plot(x, ypredict, "r-")
+plt.plot(x, y ,'ro')
+plt.axis([0.0,1.0,1.5, 7.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Linear Regression fit ')
+plt.show()
+
+!ec
+The function _coef_ gives us the parameter $\beta$ of our fit while _intercept_ yields
+$\alpha$. Depending on the constant in front of the normal distribution, we get values near or far from $alpha =2$ and $\beta =5$. Try to play around with different parameters in front of the normal distribution. The function _meansquarederror_ gives us the mean square error, a risk metric corresponding to the expected value of the squared (quadratic) error or loss defined as
+!bt
+\[ MSE(\hat{y},\hat{\tilde{y}}) = \frac{1}{n}
+\sum_{i=0}^{n-1}(y_i-\tilde{y}_i)^2,
+\]
+!et
+
+The smaller the value, the better the fit. Ideally we would like to
+have an MSE equal zero. The attentive reader has probably recognized
+this function as being similar to the $\chi^2$ function defined above.
+
+The _r2score_ function computes $R^2$, the coefficient of
+determination. It provides a measure of how well future samples are
+likely to be predicted by the model. Best possible score is 1.0 and it
+can be negative (because the model can be arbitrarily worse). A
+constant model that always predicts the expected value of $\hat{y}$,
+disregarding the input features, would get a $R^2$ score of $0.0$.
+
+If $\tilde{\hat{y}}_i$ is the predicted value of the $i-th$ sample and $y_i$ is the corresponding true value, then the score $R^2$ is defined as
+!bt
+\[
+R^2(\hat{y}, \tilde{\hat{y}}) = 1 - \frac{\sum_{i=0}^{n - 1} (y_i - \tilde{y}_i)^2}{\sum_{i=0}^{n - 1} (y_i - \bar{y})^2},
+\]
+!et
+where we have defined the mean value of $\hat{y}$ as
+!bt
+\[
+\bar{y} = \frac{1}{n} \sum_{i=0}^{n - 1} y_i.
+\]
+!et
+Another quantity taht we will meet again in our discussions of regression analysis is
+ the mean absolute error (MAE), a risk metric corresponding to the expected value of the absolute error loss or what we call the $l1$-norm loss. In our discussion above we presented the relative error.
+The MAE is defined as follows
+!bt
+\[
+\text{MAE}(\hat{y}, \hat{\tilde{y}}) = \frac{1}{n} \sum_{i=0}^{n-1} \left| y_i - \tilde{y}_i \right|.
+\]
+!et
+We present the
+squared logarithmic (quadratic) error
+!bt
+\[
+\text{MSLE}(\hat{y}, \hat{\tilde{y}}) = \frac{1}{n} \sum_{i=0}^{n - 1} (\log_e (1 + y_i) - \log_e (1 + \tilde{y}_i) )^2,
+\]
+!et
+
+where $\log_e (x)$ stands for the natural logarithm of $x$. This error
+estimate is best to use when targets having exponential growth, such
+as population counts, average sales of a commodity over a span of
+years etc.
+
+
+
+We conclude this part with another example. Instead of
+a linear $x$-dependence we study now a cubic polynomial and use the polynomial regression analysis tools of scikit-learn.
+
+!bc pycod
+import matplotlib.pyplot as plt
+import numpy as np
+import random
+from sklearn.linear_model import Ridge
+from sklearn.preprocessing import PolynomialFeatures
+from sklearn.pipeline import make_pipeline
+from sklearn.linear_model import LinearRegression
+
+x=np.linspace(0.02,0.98,200)
+noise = np.asarray(random.sample((range(200)),200))
+y=x**3*noise
+yn=x**3*100
+poly3 = PolynomialFeatures(degree=3)
+X = poly3.fit_transform(x[:,np.newaxis])
+clf3 = LinearRegression()
+clf3.fit(X,y)
+
+Xplot=poly3.fit_transform(x[:,np.newaxis])
+poly3_plot=plt.plot(x, clf3.predict(Xplot), label='Cubic Fit')
+plt.plot(x,yn, color='red', label="True Cubic")
+plt.scatter(x, y, label='Data', color='orange', s=15)
+plt.legend()
+plt.show()
+
+def error(a):
+ for i in y:
+ err=(y-yn)/yn
+ return abs(np.sum(err))/len(err)
+
+print (error(y))
+!ec
+
+
+
+===== The Boston housing data example =====
+
+The Boston housing
+data set was originally a part of UCI Machine Learning Repository
+and has been removed now. The data set is now included in _Scikit-Learn_'s
+library. There are 506 samples and 13 feature (predictor) variables
+in this data set. The objective is to predict the value of prices of
+the house using the features (predictors) listed here.
+
+The features/predictors are
+ o CRIM: Per capita crime rate by town
+ o ZN: Proportion of residential land zoned for lots over 25000 square feet
+ o INDUS: Proportion of non-retail business acres per town
+ o CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
+ o NOX: Nitric oxide concentration (parts per 10 million)
+ o RM: Average number of rooms per dwelling
+ o AGE: Proportion of owner-occupied units built prior to 1940
+ o DIS: Weighted distances to five Boston employment centers
+ o RAD: Index of accessibility to radial highways
+ o TAX: Full-value property tax rate per USD10000
+ o B: $1000(Bk - 0.63)^2$, where $Bk$ is the proportion of [people of African American descent] by town
+ o LSTAT: Percentage of lower status of the population
+ o MEDV: Median value of owner-occupied homes in USD 1000s
+
+
+We start by importing the libraries
+!bc pycod
+import numpy as np
+import matplotlib.pyplot as plt
+
+import pandas as pd
+import seaborn as sns
+!ec
+and load the Boston Housing DataSet from _Scikit-Learn_
+
+
+!bc pycod
+from sklearn.datasets import load_boston
+
+boston_dataset = load_boston()
+
+# boston_dataset is a dictionary
+# let's check what it contains
+boston_dataset.keys()
+!ec
+Then we invoke Pandas
+!bc pycod
+boston = pd.DataFrame(boston_dataset.data, columns=boston_dataset.feature_names)
+boston.head()
+boston['MEDV'] = boston_dataset.target
+!ec
+and preprocess the data
+!bc pycod
+# check for missing values in all the columns
+boston.isnull().sum()
+!ec
+We can then visualize the data
+!bc pycod
+# set the size of the figure
+sns.set(rc={'figure.figsize':(11.7,8.27)})
+
+# plot a histogram showing the distribution of the target values
+sns.distplot(boston['MEDV'], bins=30)
+plt.show()
+!ec
+
+It is now useful to look at the correlation matrix
+!bc pycod
+# compute the pair wise correlation for all columns
+correlation_matrix = boston.corr().round(2)
+# use the heatmap function from seaborn to plot the correlation matrix
+# annot = True to print the values inside the square
+sns.heatmap(data=correlation_matrix, annot=True)
+!ec
+From the above coorelation plot we can see that _MEDV_ is strongly correlated to _LSTAT_ and _RM_. We see also that _RAD_ and _TAX_ are stronly correlated, but we don't include this in our features together to avoid multi-colinearity
+
+!bc pycod
+plt.figure(figsize=(20, 5))
+
+features = ['LSTAT', 'RM']
+target = boston['MEDV']
+
+for i, col in enumerate(features):
+ plt.subplot(1, len(features) , i+1)
+ x = boston[col]
+ y = target
+ plt.scatter(x, y, marker='o')
+ plt.title(col)
+ plt.xlabel(col)
+ plt.ylabel('MEDV')
+!ec
+Now we start training our model
+!bc pycod
+X = pd.DataFrame(np.c_[boston['LSTAT'], boston['RM']], columns = ['LSTAT','RM'])
+Y = boston['MEDV']
+!ec
+We split the data into training and test sets
+
+!bc pycod
+from sklearn.model_selection import train_test_split
+
+# splits the training and test data set in 80% : 20%
+# assign random_state to any value.This ensures consistency.
+X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.2, random_state=5)
+print(X_train.shape)
+print(X_test.shape)
+print(Y_train.shape)
+print(Y_test.shape)
+!ec
+Then we use the linear regression functionality from _Scikit-Learn_
+!bc pycod
+from sklearn.linear_model import LinearRegression
+from sklearn.metrics import mean_squared_error, r2_score
+
+lin_model = LinearRegression()
+lin_model.fit(X_train, Y_train)
+
+# model evaluation for training set
+
+y_train_predict = lin_model.predict(X_train)
+rmse = (np.sqrt(mean_squared_error(Y_train, y_train_predict)))
+r2 = r2_score(Y_train, y_train_predict)
+
+print("The model performance for training set")
+print("--------------------------------------")
+print('RMSE is {}'.format(rmse))
+print('R2 score is {}'.format(r2))
+print("\n")
+
+# model evaluation for testing set
+
+y_test_predict = lin_model.predict(X_test)
+# root mean square error of the model
+rmse = (np.sqrt(mean_squared_error(Y_test, y_test_predict)))
+
+# r-squared score of the model
+r2 = r2_score(Y_test, y_test_predict)
+
+print("The model performance for testing set")
+print("--------------------------------------")
+print('RMSE is {}'.format(rmse))
+print('R2 score is {}'.format(r2))
+!ec
+
+!bc pycod
+# plotting the y_test vs y_pred
+# ideally should have been a straight line
+plt.scatter(Y_test, y_test_predict)
+plt.show()
+!ec
+
+
+
+
+Many Machine Learning problems involve thousands or even millions of
+features for each training instance. Not only does this make training
+extremely slow, it can also make it much harder to find a good
+solution, as we will see. This problem is often referred to as the
+curse of dimensionality. Fortunately, in real-world problems, it is
+often possible to reduce the number of features considerably, turning
+an intractable problem into a tractable one.
+
+Later we will discuss some of the most popular dimensionality reduction
+techniques: the principal component analysis (PCA), Kernel PCA, and
+Locally Linear Embedding (LLE).
+
+
+Principal component analysis and its various variants deal with the
+problem of fitting a low-dimensional "affine
+subspace":"https://en.wikipedia.org/wiki/Affine_space" to a set of of
+data points in a high-dimensional space. With its family of methods it
+is one of the most used tools in data modeling, compression and
+visualization.
+
+
+Before we proceed however, we will discuss how to preprocess our
+data. Till now and in connection with our previous examples we have
+not met so many cases where we are too sensitive to the scaling of our
+data. Normally the data may need a rescaling and/or may be sensitive
+to extreme values. Scaling the data renders our inputs much more
+suitable for the algorithms we want to employ.
+
+_Scikit-Learn_ has several functions which allow us to rescale the
+data, normally resulting in much better results in terms of various
+accuracy scores. The _StandardScaler_ function in _Scikit-Learn_
+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). This scaling has the drawback that it does not ensure that
+we have a particular maximum or minimum in our data set. Another
+function included in _Scikit-Learn_ is the _MinMaxScaler_ which
+ensures that all features are exactly between $0$ and $1$. The
+
+
+The _Normalizer_ scales each data
+point such that the feature vector has a euclidean length of one. In other words, it
+projects a data point on the circle (or sphere in the case of higher dimensions) with a
+radius of 1. This means every data point is scaled by a different number (by the
+inverse of it’s length).
+This normalization is often used when only the direction (or angle) of the data matters,
+not the length of the feature vector.
+
+The _RobustScaler_ works similarly to the StandardScaler in that it
+ensures statistical properties for each feature that guarantee that
+they are on the same scale. However, the RobustScaler uses the median
+and quartiles, instead of mean and variance. This makes the
+RobustScaler ignore data points that are very different from the rest
+(like measurement errors). These odd data points are also called
+outliers, and might often lead to trouble for other scaling
+techniques.
+
diff --git a/doc/src/How2ReadData/add.dlog b/doc/src/How2ReadData/add.dlog
new file mode 100644
index 000000000..8838deedd
--- /dev/null
+++ b/doc/src/How2ReadData/add.dlog
@@ -0,0 +1,11 @@
+Translating doconce text in add.do.txt to ipynb
+*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax)
+
+*** warning: latex envir \begin{cases} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
+Failed to remove ans_at_end environment
+Failed to remove sol_at_end environment
+output in add.ipynb
+Translating doconce text in add.do.txt to ipynb
+Failed to remove ans_at_end environment
+Failed to remove sol_at_end environment
+output in add.ipynb
diff --git a/doc/src/How2ReadData/add.do.txt~ b/doc/src/How2ReadData/add.do.txt~
new file mode 100644
index 000000000..c14a7b2ff
--- /dev/null
+++ b/doc/src/How2ReadData/add.do.txt~
@@ -0,0 +1,519 @@
+
+===== Simple linear regression model using _scikit-learn_ =====
+
+We start with perhaps our simplest possible example, using _Scikit-Learn_ to perform linear regression analysis on a data set produced by us.
+
+What follows is a simple Python code where we have defined a function
+$y$ in terms of the variable $x$. Both are defined as vectors with $100$ entries.
+The numbers in the vector $\hat{x}$ are given
+by random numbers generated with a uniform distribution with entries
+$x_i \in [0,1]$ (more about probability distribution functions
+later). These values are then used to define a function $y(x)$
+(tabulated again as a vector) with a linear dependence on $x$ plus a
+random noise added via the normal distribution.
+
+
+The Numpy functions are imported used the _import numpy as np_
+statement and the random number generator for the uniform distribution
+is called using the function _np.random.rand()_, where we specificy
+that we want $100$ random variables. Using Numpy we define
+automatically an array with the specified number of elements, $100$ in
+our case. With the Numpy function _randn()_ we can compute random
+numbers with the normal distribution (mean value $\mu$ equal to zero and
+variance $\sigma^2$ set to one) and produce the values of $y$ assuming a linear
+dependence as function of $x$
+
+!bt
+\[
+y = 2x+N(0,1),
+\]
+!et
+
+where $N(0,1)$ represents random numbers generated by the normal
+distribution. From _Scikit-Learn_ we import then the
+_LinearRegression_ functionality and make a prediction $\tilde{y} =
+\alpha + \beta x$ using the function _fit(x,y)_. We call the set of
+data $(\hat{x},\hat{y})$ for our training data. The Python package
+_scikit-learn_ has also a functionality which extracts the above
+fitting parameters $\alpha$ and $\beta$ (see below). Later we will
+distinguish between training data and test data.
+
+For plotting we use the Python package
+"matplotlib":"https://matplotlib.org/" which produces publication
+quality figures. Feel free to explore the extensive
+"gallery":"https://matplotlib.org/gallery/index.html" of examples. In
+this example we plot our original values of $x$ and $y$ as well as the
+prediction _ypredict_ ($\tilde{y}$), which attempts at fitting our
+data with a straight line.
+
+The Python code follows here.
+!bc pycod
+# Importing various packages
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import LinearRegression
+
+x = np.random.rand(100,1)
+y = 2*x+np.random.randn(100,1)
+linreg = LinearRegression()
+linreg.fit(x,y)
+xnew = np.array([[0],[1]])
+ypredict = linreg.predict(xnew)
+
+plt.plot(xnew, ypredict, "r-")
+plt.plot(x, y ,'ro')
+plt.axis([0,1.0,0, 5.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Simple Linear Regression')
+plt.show()
+!ec
+
+This example serves several aims. It allows us to demonstrate several
+aspects of data analysis and later machine learning algorithms. The
+immediate visualization shows that our linear fit is not
+impressive. It goes through the data points, but there are many
+outliers which are not reproduced by our linear regression. We could
+now play around with this small program and change for example the
+factor in front of $x$ and the normal distribution. Try to change the
+function $y$ to
+
+!bt
+\[
+y = 10x+0.01 \times N(0,1),
+\]
+!et
+
+where $x$ is defined as before. Does the fit look better? Indeed, by
+reducing the role of the noise given by the normal distribution we see immediately that
+our linear prediction seemingly reproduces better the training
+set. However, this testing 'by the eye' is obviouly not satisfactory in the
+long run. Here we have only defined the training data and our model, and
+have not discussed a more rigorous approach to the _cost_ function.
+
+We need more rigorous criteria in defining whether we have succeeded or
+not in modeling our training data. You will be surprised to see that
+many scientists seldomly venture beyond this 'by the eye' approach. A
+standard approach for the *cost* function is the so-called $\chi^2$
+function (a variant of the mean-squared error (MSE))
+
+!bt
+\[ \chi^2 = \frac{1}{n}
+\sum_{i=0}^{n-1}\frac{(y_i-\tilde{y}_i)^2}{\sigma_i^2},
+\]
+!et
+
+where $\sigma_i^2$ is the variance (to be defined later) of the entry
+$y_i$. We may not know the explicit value of $\sigma_i^2$, it serves
+however the aim of scaling the equations and make the cost function
+dimensionless.
+
+Minimizing the cost function is a central aspect of
+our discussions to come. Finding its minima as function of the model
+parameters ($\alpha$ and $\beta$ in our case) will be a recurring
+theme in these series of lectures. Essentially all machine learning
+algorithms we will discuss center around the minimization of the
+chosen cost function. This depends in turn on our specific
+model for describing the data, a typical situation in supervised
+learning. Automatizing the search for the minima of the cost function is a
+central ingredient in all algorithms. Typical methods which are
+employed are various variants of _gradient_ methods. These will be
+discussed in more detail later. Again, you'll be surprised to hear that
+many practitioners minimize the above function ''by the eye', popularly dubbed as
+'chi by the eye'. That is, change a parameter and see (visually and numerically) that
+the $\chi^2$ function becomes smaller.
+
+There are many ways to define the cost function. A simpler approach is to look at the relative difference between the training data and the predicted data, that is we define
+the relative error (why would we prefer the MSE instead of the relative error?) as
+
+!bt
+\[
+\epsilon_{\mathrm{relative}}= \frac{\vert \hat{y} -\hat{\tilde{y}}\vert}{\vert \hat{y}\vert}.
+\]
+!et
+
+The squared cost function results in an arithmetic mean-unbiased
+estimator, and the absolute-value cost function results in a
+median-unbiased estimator (in the one-dimensional case, and a
+geometric median-unbiased estimator for the multi-dimensional
+case). The squared cost function has the disadvantage that it has the tendency
+to be dominated by outliers.
+
+We can modify easily the above Python code and plot the relative error instead
+!bc pycod
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import LinearRegression
+
+x = np.random.rand(100,1)
+y = 5*x+0.01*np.random.randn(100,1)
+linreg = LinearRegression()
+linreg.fit(x,y)
+ypredict = linreg.predict(x)
+
+plt.plot(x, np.abs(ypredict-y)/abs(y), "ro")
+plt.axis([0,1.0,0.0, 0.5])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$\epsilon_{\mathrm{relative}}$')
+plt.title(r'Relative error')
+plt.show()
+!ec
+
+Depending on the parameter in front of the normal distribution, we may
+have a small or larger relative error. Try to play around with
+different training data sets and study (graphically) the value of the
+relative error.
+
+As mentioned above, _Scikit-Learn_ has an impressive functionality.
+We can for example extract the values of $\alpha$ and $\beta$ and
+their error estimates, or the variance and standard deviation and many
+other properties from the statistical data analysis.
+
+Here we show an
+example of the functionality of _Scikit-Learn_.
+!bc pycod
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import LinearRegression
+from sklearn.metrics import mean_squared_error, r2_score, mean_squared_log_error, mean_absolute_error
+
+x = np.random.rand(100,1)
+y = 2.0+ 5*x+0.5*np.random.randn(100,1)
+linreg = LinearRegression()
+linreg.fit(x,y)
+ypredict = linreg.predict(x)
+print('The intercept alpha: \n', linreg.intercept_)
+print('Coefficient beta : \n', linreg.coef_)
+# The mean squared error
+print("Mean squared error: %.2f" % mean_squared_error(y, ypredict))
+# Explained variance score: 1 is perfect prediction
+print('Variance score: %.2f' % r2_score(y, ypredict))
+# Mean squared log error
+print('Mean squared log error: %.2f' % mean_squared_log_error(y, ypredict) )
+# Mean absolute error
+print('Mean absolute error: %.2f' % mean_absolute_error(y, ypredict))
+plt.plot(x, ypredict, "r-")
+plt.plot(x, y ,'ro')
+plt.axis([0.0,1.0,1.5, 7.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Linear Regression fit ')
+plt.show()
+
+!ec
+The function _coef_ gives us the parameter $\beta$ of our fit while _intercept_ yields
+$\alpha$. Depending on the constant in front of the normal distribution, we get values near or far from $alpha =2$ and $\beta =5$. Try to play around with different parameters in front of the normal distribution. The function _meansquarederror_ gives us the mean square error, a risk metric corresponding to the expected value of the squared (quadratic) error or loss defined as
+!bt
+\[ MSE(\hat{y},\hat{\tilde{y}}) = \frac{1}{n}
+\sum_{i=0}^{n-1}(y_i-\tilde{y}_i)^2,
+\]
+!et
+
+The smaller the value, the better the fit. Ideally we would like to
+have an MSE equal zero. The attentive reader has probably recognized
+this function as being similar to the $\chi^2$ function defined above.
+
+The _r2score_ function computes $R^2$, the coefficient of
+determination. It provides a measure of how well future samples are
+likely to be predicted by the model. Best possible score is 1.0 and it
+can be negative (because the model can be arbitrarily worse). A
+constant model that always predicts the expected value of $\hat{y}$,
+disregarding the input features, would get a $R^2$ score of $0.0$.
+
+If $\tilde{\hat{y}}_i$ is the predicted value of the $i-th$ sample and $y_i$ is the corresponding true value, then the score $R^2$ is defined as
+!bt
+\[
+R^2(\hat{y}, \tilde{\hat{y}}) = 1 - \frac{\sum_{i=0}^{n - 1} (y_i - \tilde{y}_i)^2}{\sum_{i=0}^{n - 1} (y_i - \bar{y})^2},
+\]
+!et
+where we have defined the mean value of $\hat{y}$ as
+!bt
+\[
+\bar{y} = \frac{1}{n} \sum_{i=0}^{n - 1} y_i.
+\]
+!et
+Another quantity taht we will meet again in our discussions of regression analysis is
+ the mean absolute error (MAE), a risk metric corresponding to the expected value of the absolute error loss or what we call the $l1$-norm loss. In our discussion above we presented the relative error.
+The MAE is defined as follows
+!bt
+\[
+\text{MAE}(\hat{y}, \hat{\tilde{y}}) = \frac{1}{n} \sum_{i=0}^{n-1} \left| y_i - \tilde{y}_i \right|.
+\]
+!et
+We present the
+squared logarithmic (quadratic) error
+!bt
+\[
+\text{MSLE}(\hat{y}, \hat{\tilde{y}}) = \frac{1}{n} \sum_{i=0}^{n - 1} (\log_e (1 + y_i) - \log_e (1 + \tilde{y}_i) )^2,
+\]
+!et
+
+where $\log_e (x)$ stands for the natural logarithm of $x$. This error
+estimate is best to use when targets having exponential growth, such
+as population counts, average sales of a commodity over a span of
+years etc.
+
+
+Finally, another cost function is the Huber cost function used in robust regression.
+
+The rationale behind this possible cost function is its reduced
+sensitivity to outliers in the data set. In our discussions on
+dimensionality reduction and normalization of data we will meet other
+ways of dealing with outliers.
+
+The Huber cost function is defined as
+!bt
+\[
+H_{\delta}(a)={\begin{cases}{\frac {1}{2}}{a^{2}}&{\text{for }}|a|\leq \delta ,\\\delta (|a|-{\frac {1}{2}}\delta ),&{\text{otherwise.}}\end{cases}}}.
+\]
+!et
+Here $a=\bm{y} - \bm{\tilde{y}}$.
+We will discuss in more
+detail these and other functions in the various lectures. We conclude this part with another example. Instead of
+a linear $x$-dependence we study now a cubic polynomial and use the polynomial regression analysis tools of scikit-learn.
+
+!bc pycod
+import matplotlib.pyplot as plt
+import numpy as np
+import random
+from sklearn.linear_model import Ridge
+from sklearn.preprocessing import PolynomialFeatures
+from sklearn.pipeline import make_pipeline
+from sklearn.linear_model import LinearRegression
+
+x=np.linspace(0.02,0.98,200)
+noise = np.asarray(random.sample((range(200)),200))
+y=x**3*noise
+yn=x**3*100
+poly3 = PolynomialFeatures(degree=3)
+X = poly3.fit_transform(x[:,np.newaxis])
+clf3 = LinearRegression()
+clf3.fit(X,y)
+
+Xplot=poly3.fit_transform(x[:,np.newaxis])
+poly3_plot=plt.plot(x, clf3.predict(Xplot), label='Cubic Fit')
+plt.plot(x,yn, color='red', label="True Cubic")
+plt.scatter(x, y, label='Data', color='orange', s=15)
+plt.legend()
+plt.show()
+
+def error(a):
+ for i in y:
+ err=(y-yn)/yn
+ return abs(np.sum(err))/len(err)
+
+print (error(y))
+!ec
+
+
+
+===== The Boston housing data example =====
+
+The Boston housing
+data set was originally a part of UCI Machine Learning Repository
+and has been removed now. The data set is now included in _Scikit-Learn_'s
+library. There are 506 samples and 13 feature (predictor) variables
+in this data set. The objective is to predict the value of prices of
+the house using the features (predictors) listed here.
+
+The features/predictors are
+ o CRIM: Per capita crime rate by town
+ o ZN: Proportion of residential land zoned for lots over 25000 square feet
+ o INDUS: Proportion of non-retail business acres per town
+ o CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
+ o NOX: Nitric oxide concentration (parts per 10 million)
+ o RM: Average number of rooms per dwelling
+ o AGE: Proportion of owner-occupied units built prior to 1940
+ o DIS: Weighted distances to five Boston employment centers
+ o RAD: Index of accessibility to radial highways
+ o TAX: Full-value property tax rate per USD10000
+ o B: $1000(Bk - 0.63)^2$, where $Bk$ is the proportion of [people of African American descent] by town
+ o LSTAT: Percentage of lower status of the population
+ o MEDV: Median value of owner-occupied homes in USD 1000s
+
+
+We start by importing the libraries
+!bc pycod
+import numpy as np
+import matplotlib.pyplot as plt
+
+import pandas as pd
+import seaborn as sns
+!ec
+and load the Boston Housing DataSet from _Scikit-Learn_
+
+
+!bc pycod
+from sklearn.datasets import load_boston
+
+boston_dataset = load_boston()
+
+# boston_dataset is a dictionary
+# let's check what it contains
+boston_dataset.keys()
+!ec
+Then we invoke Pandas
+!bc pycod
+boston = pd.DataFrame(boston_dataset.data, columns=boston_dataset.feature_names)
+boston.head()
+boston['MEDV'] = boston_dataset.target
+!ec
+and preprocess the data
+!bc pycod
+# check for missing values in all the columns
+boston.isnull().sum()
+!ec
+We can then visualize the data
+!bc pycod
+# set the size of the figure
+sns.set(rc={'figure.figsize':(11.7,8.27)})
+
+# plot a histogram showing the distribution of the target values
+sns.distplot(boston['MEDV'], bins=30)
+plt.show()
+!ec
+
+It is now useful to look at the correlation matrix
+!bc pycod
+# compute the pair wise correlation for all columns
+correlation_matrix = boston.corr().round(2)
+# use the heatmap function from seaborn to plot the correlation matrix
+# annot = True to print the values inside the square
+sns.heatmap(data=correlation_matrix, annot=True)
+!ec
+From the above coorelation plot we can see that _MEDV_ is strongly correlated to _LSTAT_ and _RM_. We see also that _RAD_ and _TAX_ are stronly correlated, but we don't include this in our features together to avoid multi-colinearity
+
+!bc pycod
+plt.figure(figsize=(20, 5))
+
+features = ['LSTAT', 'RM']
+target = boston['MEDV']
+
+for i, col in enumerate(features):
+ plt.subplot(1, len(features) , i+1)
+ x = boston[col]
+ y = target
+ plt.scatter(x, y, marker='o')
+ plt.title(col)
+ plt.xlabel(col)
+ plt.ylabel('MEDV')
+!ec
+Now we start training our model
+!bc pycod
+X = pd.DataFrame(np.c_[boston['LSTAT'], boston['RM']], columns = ['LSTAT','RM'])
+Y = boston['MEDV']
+!ec
+We split the data into training and test sets
+
+!bc pycod
+from sklearn.model_selection import train_test_split
+
+# splits the training and test data set in 80% : 20%
+# assign random_state to any value.This ensures consistency.
+X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size = 0.2, random_state=5)
+print(X_train.shape)
+print(X_test.shape)
+print(Y_train.shape)
+print(Y_test.shape)
+!ec
+Then we use the linear regression functionality from _Scikit-Learn_
+!bc pycod
+from sklearn.linear_model import LinearRegression
+from sklearn.metrics import mean_squared_error, r2_score
+
+lin_model = LinearRegression()
+lin_model.fit(X_train, Y_train)
+
+# model evaluation for training set
+
+y_train_predict = lin_model.predict(X_train)
+rmse = (np.sqrt(mean_squared_error(Y_train, y_train_predict)))
+r2 = r2_score(Y_train, y_train_predict)
+
+print("The model performance for training set")
+print("--------------------------------------")
+print('RMSE is {}'.format(rmse))
+print('R2 score is {}'.format(r2))
+print("\n")
+
+# model evaluation for testing set
+
+y_test_predict = lin_model.predict(X_test)
+# root mean square error of the model
+rmse = (np.sqrt(mean_squared_error(Y_test, y_test_predict)))
+
+# r-squared score of the model
+r2 = r2_score(Y_test, y_test_predict)
+
+print("The model performance for testing set")
+print("--------------------------------------")
+print('RMSE is {}'.format(rmse))
+print('R2 score is {}'.format(r2))
+!ec
+
+!bc pycod
+# plotting the y_test vs y_pred
+# ideally should have been a straight line
+plt.scatter(Y_test, y_test_predict)
+plt.show()
+!ec
+
+
+
+
+Many Machine Learning problems involve thousands or even millions of
+features for each training instance. Not only does this make training
+extremely slow, it can also make it much harder to find a good
+solution, as we will see. This problem is often referred to as the
+curse of dimensionality. Fortunately, in real-world problems, it is
+often possible to reduce the number of features considerably, turning
+an intractable problem into a tractable one.
+
+Later we will discuss some of the most popular dimensionality reduction
+techniques: the principal component analysis (PCA), Kernel PCA, and
+Locally Linear Embedding (LLE).
+
+
+Principal component analysis and its various variants deal with the
+problem of fitting a low-dimensional "affine
+subspace":"https://en.wikipedia.org/wiki/Affine_space" to a set of of
+data points in a high-dimensional space. With its family of methods it
+is one of the most used tools in data modeling, compression and
+visualization.
+
+
+Before we proceed however, we will discuss how to preprocess our
+data. Till now and in connection with our previous examples we have
+not met so many cases where we are too sensitive to the scaling of our
+data. Normally the data may need a rescaling and/or may be sensitive
+to extreme values. Scaling the data renders our inputs much more
+suitable for the algorithms we want to employ.
+
+_Scikit-Learn_ has several functions which allow us to rescale the
+data, normally resulting in much better results in terms of various
+accuracy scores. The _StandardScaler_ function in _Scikit-Learn_
+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). This scaling has the drawback that it does not ensure that
+we have a particular maximum or minimum in our data set. Another
+function included in _Scikit-Learn_ is the _MinMaxScaler_ which
+ensures that all features are exactly between $0$ and $1$. The
+
+
+The _Normalizer_ scales each data
+point such that the feature vector has a euclidean length of one. In other words, it
+projects a data point on the circle (or sphere in the case of higher dimensions) with a
+radius of 1. This means every data point is scaled by a different number (by the
+inverse of it’s length).
+This normalization is often used when only the direction (or angle) of the data matters,
+not the length of the feature vector.
+
+The _RobustScaler_ works similarly to the StandardScaler in that it
+ensures statistical properties for each feature that guarantee that
+they are on the same scale. However, the RobustScaler uses the median
+and quartiles, instead of mean and variance. This makes the
+RobustScaler ignore data points that are very different from the rest
+(like measurement errors). These odd data points are also called
+outliers, and might often lead to trouble for other scaling
+techniques.
+
diff --git a/doc/src/Splines/.ipynb_checkpoints/GradientMethods-checkpoint.ipynb b/doc/src/Optimization/.ipynb_checkpoints/GradientMethods-checkpoint.ipynb
similarity index 100%
rename from doc/src/Splines/.ipynb_checkpoints/GradientMethods-checkpoint.ipynb
rename to doc/src/Optimization/.ipynb_checkpoints/GradientMethods-checkpoint.ipynb
diff --git a/doc/src/Splines/GradientMethods.ipynb b/doc/src/Optimization/GradientMethods.ipynb
similarity index 100%
rename from doc/src/Splines/GradientMethods.ipynb
rename to doc/src/Optimization/GradientMethods.ipynb
diff --git a/doc/src/Splines/Splines.do.txt b/doc/src/Optimization/Splines.do.txt
similarity index 100%
rename from doc/src/Splines/Splines.do.txt
rename to doc/src/Optimization/Splines.do.txt
diff --git a/doc/src/Splines/autodiff/.ipynb_checkpoints/example_plot-checkpoint.ipynb b/doc/src/Optimization/autodiff/.ipynb_checkpoints/example_plot-checkpoint.ipynb
similarity index 100%
rename from doc/src/Splines/autodiff/.ipynb_checkpoints/example_plot-checkpoint.ipynb
rename to doc/src/Optimization/autodiff/.ipynb_checkpoints/example_plot-checkpoint.ipynb
diff --git a/doc/src/Splines/autodiff/.ipynb_checkpoints/examples_allowed_functions-Copy1-checkpoint.ipynb b/doc/src/Optimization/autodiff/.ipynb_checkpoints/examples_allowed_functions-Copy1-checkpoint.ipynb
similarity index 100%
rename from doc/src/Splines/autodiff/.ipynb_checkpoints/examples_allowed_functions-Copy1-checkpoint.ipynb
rename to doc/src/Optimization/autodiff/.ipynb_checkpoints/examples_allowed_functions-Copy1-checkpoint.ipynb
diff --git a/doc/src/Splines/autodiff/.ipynb_checkpoints/examples_allowed_functions-checkpoint.ipynb b/doc/src/Optimization/autodiff/.ipynb_checkpoints/examples_allowed_functions-checkpoint.ipynb
similarity index 100%
rename from doc/src/Splines/autodiff/.ipynb_checkpoints/examples_allowed_functions-checkpoint.ipynb
rename to doc/src/Optimization/autodiff/.ipynb_checkpoints/examples_allowed_functions-checkpoint.ipynb
diff --git a/doc/src/Splines/autodiff/autograd.pdf b/doc/src/Optimization/autodiff/autograd.pdf
similarity index 100%
rename from doc/src/Splines/autodiff/autograd.pdf
rename to doc/src/Optimization/autodiff/autograd.pdf
diff --git a/doc/src/Splines/autodiff/example_plot.ipynb b/doc/src/Optimization/autodiff/example_plot.ipynb
similarity index 100%
rename from doc/src/Splines/autodiff/example_plot.ipynb
rename to doc/src/Optimization/autodiff/example_plot.ipynb
diff --git a/doc/src/Splines/autodiff/examples_allowed_functions-Copy1.ipynb b/doc/src/Optimization/autodiff/examples_allowed_functions-Copy1.ipynb
similarity index 100%
rename from doc/src/Splines/autodiff/examples_allowed_functions-Copy1.ipynb
rename to doc/src/Optimization/autodiff/examples_allowed_functions-Copy1.ipynb
diff --git a/doc/src/Splines/autodiff/examples_allowed_functions.ipynb b/doc/src/Optimization/autodiff/examples_allowed_functions.ipynb
similarity index 100%
rename from doc/src/Splines/autodiff/examples_allowed_functions.ipynb
rename to doc/src/Optimization/autodiff/examples_allowed_functions.ipynb
diff --git a/doc/src/Optimization/backup/lecture4.pdf b/doc/src/Optimization/backup/lecture4.pdf
new file mode 100644
index 000000000..229abfdc7
Binary files /dev/null and b/doc/src/Optimization/backup/lecture4.pdf differ
diff --git a/doc/src/Optimization/backup/lecture5.pdf b/doc/src/Optimization/backup/lecture5.pdf
new file mode 100644
index 000000000..260c8749a
Binary files /dev/null and b/doc/src/Optimization/backup/lecture5.pdf differ
diff --git a/doc/src/Optimization/backup/lecture7.pdf b/doc/src/Optimization/backup/lecture7.pdf
new file mode 100644
index 000000000..b2433bbc0
Binary files /dev/null and b/doc/src/Optimization/backup/lecture7.pdf differ
diff --git a/doc/src/Splines/beamerthemered_plain.sty b/doc/src/Optimization/beamerthemered_plain.sty
similarity index 100%
rename from doc/src/Splines/beamerthemered_plain.sty
rename to doc/src/Optimization/beamerthemered_plain.sty
diff --git a/doc/src/Splines/beamerthemered_shadow.sty b/doc/src/Optimization/beamerthemered_shadow.sty
similarity index 100%
rename from doc/src/Splines/beamerthemered_shadow.sty
rename to doc/src/Optimization/beamerthemered_shadow.sty
diff --git a/doc/src/Splines/clean.sh b/doc/src/Optimization/clean.sh
similarity index 100%
rename from doc/src/Splines/clean.sh
rename to doc/src/Optimization/clean.sh
diff --git a/doc/src/Splines/logsgd.py b/doc/src/Optimization/logsgd.py
similarity index 100%
rename from doc/src/Splines/logsgd.py
rename to doc/src/Optimization/logsgd.py
diff --git a/doc/src/Splines/make.sh b/doc/src/Optimization/make.sh
similarity index 100%
rename from doc/src/Splines/make.sh
rename to doc/src/Optimization/make.sh
diff --git a/doc/src/Splines/ridgesgd.py b/doc/src/Optimization/ridgesgd.py
similarity index 100%
rename from doc/src/Splines/ridgesgd.py
rename to doc/src/Optimization/ridgesgd.py
diff --git a/doc/src/week37/olscode.py~ b/doc/src/week37/olscode.py~
new file mode 100644
index 000000000..a6094b632
--- /dev/null
+++ b/doc/src/week37/olscode.py~
@@ -0,0 +1,47 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.model_selection import KFold
+from sklearn.linear_model import Ridge, LinearRegression
+from sklearn.model_selection import cross_val_score
+from sklearn.preprocessing import PolynomialFeatures
+
+# A seed just to ensure that the random numbers are the same for every run.
+# Useful for eventual debugging.
+np.random.seed(3155)
+
+# Generate the data.
+nsamples = 1000
+x = np.random.randn(nsamples)
+y = 3*x**2 + np.random.randn(nsamples)
+
+## Cross-validation on Ridge regression using KFold only
+
+# Decide degree on polynomial to fit
+poly = PolynomialFeatures(degree = 6)
+
+
+# Initialize a KFold instance
+k = 10
+kfold = KFold(n_splits = k)
+
+# Perform the cross-validation to estimate MSE using OLS
+scores_KFold = np.zeros((k))
+model = LinearRegression()
+j = 0
+for train_inds, test_inds in kfold.split(x):
+ xtrain = x[train_inds]
+ ytrain = y[train_inds]
+ xtest = x[test_inds]
+ ytest = y[test_inds]
+ Xtrain = poly.fit_transform(xtrain[:, np.newaxis])
+ model.fit(Xtrain, ytrain[:, np.newaxis])
+ Xtest = poly.fit_transform(xtest[:, np.newaxis])
+ ypred = model.predict(Xtest)
+ scores_KFold[j] = np.sum((ypred - ytest[:, np.newaxis])**2)/np.size(ypred)
+ print(f"Score for each fold:{scores_KFold[j]}")
+ j += 1
+
+
+estimated_mse_KFold = np.mean(scores_KFold)
+print(f"Average OLS score:{estimated_mse_KFold}")
+
diff --git a/doc/src/week39/codes/test10.py~ b/doc/src/week39/codes/test10.py~
new file mode 100644
index 000000000..fbeb00c02
--- /dev/null
+++ b/doc/src/week39/codes/test10.py~
@@ -0,0 +1,44 @@
+
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import cm
+from matplotlib.ticker import LinearLocator, FormatStrFormatter
+import sys
+
+# the number of datapoints
+n = 100
+x = 2*np.random.rand(n,1)
+y = 4+3*x+np.random.randn(n,1)
+
+X = np.c_[np.ones((n,1)), x]
+XT_X = X.T @ X
+
+#Ridge parameter lambda
+lmbda = 0.001
+Id = lmbda* np.eye(XT_X.shape[0])
+
+beta_linreg = np.linalg.inv(XT_X+Id) @ X.T @ y
+print(beta_linreg)
+# Start plain gradient descent
+beta = np.random.randn(2,1)
+
+eta = 0.1
+Niterations = 100
+
+for iter in range(Niterations):
+ gradients = 2.0/n*X.T @ (X @ (beta)-y)+2*lmbda*beta
+ beta -= eta*gradients
+
+print(beta)
+ypredict = X @ beta
+ypredict2 = X @ beta_linreg
+plt.plot(x, ypredict, "r-")
+plt.plot(x, ypredict2, "b-")
+plt.plot(x, y ,'ro')
+plt.axis([0,2.0,0, 15.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Gradient descent example for Ridge')
+plt.show()
diff --git a/doc/src/week39/codes/test12.py~ b/doc/src/week39/codes/test12.py~
new file mode 100644
index 000000000..7abc500c5
--- /dev/null
+++ b/doc/src/week39/codes/test12.py~
@@ -0,0 +1,1170 @@
+import os
+import sys
+import pytest
+import numba
+from matplotlib import cm
+from matplotlib.ticker import LinearLocator, FormatStrFormatter
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+import numpy as np
+import functools
+import time
+from numba import jit, njit
+from PIL import Image
+import pandas as pd
+import seaborn as sns
+sns.set()
+import math
+
+from sklearn.model_selection import train_test_split, cross_val_score, KFold
+from sklearn.preprocessing import StandardScaler
+from sklearn.linear_model import SGDRegressor, LinearRegression, LogisticRegression, Lasso, Ridge
+from sklearn import datasets
+from sklearn.metrics import confusion_matrix, mean_squared_error, r2_score
+from sklearn.utils import resample
+
+
+# Bootstrap
+
+def Bootstrap(x1,x2, y, N_boot=500, method = 'ols', degrees = 5, random_state = 42):
+ """
+ Computes bias^2, variance and the mean squared error using bootstrap resampling method
+ for the provided data and the method.
+
+ Arguments:
+ x1: 1D numpy array, covariate
+ x2: 1D numpy array, covariate
+ N_boot: integer type, the number of bootstrap samples
+ method: string type, accepts 'ols', 'ridge' or 'lasso' as arguments
+ degree: integer type, polynomial degree for generating the design matrix
+ random_state: integer, ensures the same split when using the train_test_split functionality
+
+ Returns: Bias_vec, Var_vec, MSE_vec, betaVariance_vec
+ numpy arrays. Bias, Variance, MSE and the variance of beta for the predicted model
+ """
+ ##split x1, x2 and y arrays as a train and test data and generate design matrix
+ x1_train, x1_test,x2_train, x2_test, y_train, y_test = train_test_split(x1,x2, y, test_size=0.2, random_state = random_state)
+ y_pred_test = np.zeros((y_test.shape[0], N_boot))
+ X_test = designMatrix(x1_test, x2_test, degrees)
+
+ betaMatrix = np.zeros((X_test.shape[1], N_boot))
+
+ ##resample and fit the corresponding method on the train data
+ for i in range(N_boot):
+ x1_,x2_, y_ = resample(x1_train, x2_train, y_train)
+ X_train = designMatrix(x1_, x2_, degrees)
+ scaler = StandardScaler()
+ scaler.fit(X_train)
+ X_train = scaler.transform(X_train)
+ X_train[:, 0] = 1
+ X_test = designMatrix(x1_test, x2_test, degrees)
+ X_test = scaler.transform(X_test)
+ X_test[:, 0] = 1
+
+ if method == 'ols':
+ manual_regression = linregOwn(method = 'ols')
+ beta = manual_regression.fit(X_train, y_)
+ if method == 'ridge':
+ manual_regression = linregOwn(method = 'ridge')
+ beta = manual_regression.fit(X_train, y_, lambda_ = 0.05)
+ if method == 'lasso':
+ manual_regression = linregOwn(method = 'lasso')
+ beta = manual_regression.fit(X_train, y_, lambda_ = 0.05)
+
+ ##predict on the same test data
+ y_pred_test[:, i] = np.dot(X_test, beta)
+ betaMatrix[:, i] = beta
+ y_test = y_test.reshape(len(y_test),1)
+
+ Bias_vec = []
+ Var_vec = []
+ MSE_vec = []
+ betaVariance_vec = []
+ R2_score = []
+ y_test = y_test.reshape(len(y_test),1)
+ MSE = np.mean( np.mean((y_test - y_pred_test)**2, axis=1, keepdims=True) )
+ bias = np.mean( (y_test - np.mean(y_pred_test, axis=1, keepdims=True))**2 )
+ variance = np.mean( np.var(y_pred_test, axis=1, keepdims=True) )
+ betaVariance = np.var(betaMatrix, axis=1)
+ print("-------------------------------------------------------------")
+ print("Degree: %d" % degrees)
+ print('MSE:', np.round(MSE, 3))
+ print('Bias^2:', np.round(bias, 3))
+ print('Var:', np.round(variance,3))
+ print('{} >= {} + {} = {}'.format(MSE, bias, variance, bias+variance))
+ print("-------------------------------------------------------------")
+
+ Bias_vec.append(bias)
+ Var_vec.append(variance)
+ MSE_vec.append(MSE)
+ betaVariance_vec.append(betaVariance)
+ return Bias_vec, Var_vec, MSE_vec, betaVariance_vec
+
+
+
+class CrossValidation:
+ """
+ A class of cross-validation technique. Performs cross-validation with shuffling.
+ """
+ def __init__(self, LinearRegression, DesignMatrix):
+ """
+ Initialization
+
+ Arguments:
+ LinearRegression: Instance from the class created by either linregOwn or linregSKl
+ DesignMatrix: Function that generates design matrix
+ """
+ self.LinearRegression = LinearRegression
+ self.DesignMatrix = DesignMatrix
+
+ def kFoldCV(self, x1, x2, y, k = 10, lambda_ = 0, degree = 5):
+ """
+ Performs shuffling of the data, holds a split of the data as a test set at each split and evaluates the model
+ on the rest of the data.
+ Calculates the MSE , R2_score, variance, bias on the test data and MSE on the train data.
+
+ Arguments:
+ x1: 1D numpy array
+ x2: 1D numpy array
+ y: 1D numpy array
+ k: integer, the number of splits
+ lambda_: float type, shrinkage parameter for ridge and lasso methods.
+ degree: integer type, the number of polynomials, complexity parameter
+
+ """
+ self.lambda_ = lambda_
+ M = x1.shape[0]//k ## Split input data x in k folds of size M
+
+
+ ##save the statistic in the list
+ MSE_train = []
+ MSE_k = []
+ R2_k = []
+ var_k = []
+ bias_k = []
+
+ ##shuffle the data randomly
+ shf = np.random.permutation(x1.size)
+ x1_shuff = x1[shf]
+ x2_shuff = x2[shf]
+ y_shuff = y[shf]
+
+ for i in range(k):
+ # x_k and y_k are the hold out data for fold k
+ x1_k = x1_shuff[i*M:(i+1)*M]
+ x2_k = x2_shuff[i*M:(i+1)*M]
+ y_k = y_shuff[i*M:(i+1)*M]
+
+ ## Generate train data and then scale both train and test
+ index_true = np.array([True for i in range(x1.shape[0])])
+ index_true[i*M:(i+1)*M] = False
+ X_train = self.DesignMatrix(x1_shuff[index_true], x2_shuff[index_true], degree)
+ y_train = y_shuff[index_true]
+ scaler = StandardScaler()
+ scaler.fit(X_train)
+ X_train = scaler.transform(X_train)
+ X_train[:, 0] = 1
+
+ ### Fit the regression on the train data
+ beta = self.LinearRegression.fit(X_train, y_train, lambda_)
+ y_predict_train = np.dot(X_train, beta)
+ MSE_train.append(np.sum( (y_train-y_predict_train)**2)/len(y_train))
+
+ ## Predict on the hold out data and calculate statistic of interest
+ X_k = self.DesignMatrix(x1_k, x2_k, degree)
+ X_k = scaler.transform(X_k)
+ X_k[:, 0] = 1
+ y_predict = np.dot(X_k,beta)
+ MSE_k.append(np.sum((y_k-y_predict)**2, axis=0, keepdims=True)/len(y_predict))
+ R2_k.append(1.0 - np.sum((y_k - y_predict)**2, axis=0, keepdims=True) / np.sum((y_k - np.mean(y_k))**2, axis=0, keepdims=True) )
+ var_k.append(np.var(y_predict,axis=0, keepdims=True))
+ bias_k.append((y_k - np.mean(y_predict, axis=0, keepdims=True))**2 )
+
+ means = [np.mean(MSE_k), np.mean(R2_k), np.mean(var_k),
+ np.mean(bias_k),np.mean(MSE_train)]
+ #print('MSE_test: {}' .format(np.round(np.mean(MSE_k),3)))
+ #print('R2: {}' .format(np.round(np.mean(R2_k),3)))
+ #print('Variance of the predicted outcome: {}' .format(np.round(np.mean(var_k),3)))
+ #print('Bias: {}' .format(np.round(np.mean(bias_k),3)))
+ #print('MSE_train {}' .format(np.round(np.mean(MSE_train),3)))
+ return means
+
+
+# Franke Function
+
+def franke(x, y):
+ """
+ Computes Franke function.
+ Franke's function has two Gaussian peaks of different heights, and a smaller dip.
+ It is used as a test function in interpolation problems.
+
+ Franke's function is normally defined on the grid [0, 1] for each x, y.
+
+ Arguments of the function:
+ x : numpy array
+ y : numpy array
+
+ Output of the function:
+ f : Franke function values at specific coordinate points of x and y
+ """
+ f = (0.75 * np.exp(-((9*x - 2)**2)/4 - ((9*y - 2)**2)/4 )
+ + 0.75 * np.exp(-((9*x + 1)**2)/49 - (9*y + 1) /10)
+ + 0.5 * np.exp(-((9*x - 7)**2)/4 - ((9*y - 3)**2)/4 )
+ - 0.2 * np.exp(-((9*x - 4)**2) - ((9*y - 7)**2) ))
+ return f
+
+
+
+
+class linregOwn:
+ """
+ A class of linear regressions. Perform ordinarly least squares (OLS) and Ridge regression manually. Lasso
+ is performed using scikit-learn functionality.
+ """
+ def __init__(self, method = 'ols'):
+ """
+ Constructor
+
+ Determines the method used in the fitting
+
+ Arguments:
+ method: string type. Accepts either 'ols', 'ridge' or 'lasso'.
+
+ """
+ self.method = method
+ self.yHat = None
+ self.X = None
+ self.y = None
+ self.beta = None
+
+ self._MSE = None
+ self._R2 = None
+ self._betaVariance = None
+ self.lambda_ = None
+
+ def fit(self, X_train, y_train, lambda_ = 0):
+ """
+ Performs the fit of OLS, Ridge or Lasso, depending on the argument provided initially.
+
+ Arguments:
+ X_train: Covariate matrix of the train data set, i.e. design matrix of
+ the shape m x p where m is the number of rows and p is the number of columns
+ (i.e. p is the complexity parameter).
+ y_train: Outcome variable, 1D numpy array
+ lambda_: float type. Shrinkage parameter for ridge and lasso methods. The higher value, higher shrinkage.
+ lambda_ is set to 0 for the OLS regression
+
+ """
+ self.X_train = X_train
+ self.y_train = y_train
+ self.lambda_ = lambda_
+ if self.method == 'ols':
+ self._olsFit(X_train, y_train)
+ if self.method == 'ridge':
+ self._ridgeFit(X_train, y_train, lambda_)
+ if self.method == 'lasso':
+ self._lassoFitSKL(X_train, y_train, lambda_)
+ return self.beta
+
+ def _olsFit(self, X_train, y_train):
+ """
+ Performs the ordinary least squares (OLS) fit on the provided data using singular value decomposition(SVD).
+
+
+ Arguments:
+
+ X_train: Covariate matrix of the train data set, i.e. design matrix of
+ the shape m x p where m is the number of rows and p is the number of columns
+ (i.e. p is the complexity parameter).
+ y_train: Outcome variable, 1D numpy array
+
+ Returns:
+ beta : numpy.array
+ The beta parameters from the performed fit
+ """
+ self.X_train = X_train
+ self.y_test = y_train
+ U, S, VT = np.linalg.svd(self.X_train, full_matrices=True)
+ S_inverse = np.zeros(shape=self.X_train.shape)
+ ##S is a vector, with shape of the number of columns
+ S_inverse[:S.shape[0], :S.shape[0]] = np.diag(1/S)
+ self.beta = np.dot(VT.T, np.dot(S_inverse.T, np.dot(U.T, self.y_train)))
+ #self.beta = np.linalg.inv(np.dot(X.T,X)).dot(X.T, y)
+
+ def _ridgeFit(self, X_train, y_train, lambda_):
+ """
+ Performs the ridge regression fit
+
+ Arguments:
+ X_train: Covariate matrix of the train data set, design matrix of
+ the shape m x p (m_train_rows, p_columns).
+ y_train: Outcome variable, 1D numpy array, dimension m x 1
+ lambda_: Integer type. The shrinkage parameter
+
+ Returns:
+ beta : numpy.array
+ The beta parameters from the performed fit
+ """
+ self.X_train = X_train
+ self.y_train = y_train
+ self.lambda_ = lambda_
+ self.beta = np.dot(np.linalg.inv(np.dot(X_train.T,X_train) + self.lambda_ * np.eye(X_train.shape[1])), np.dot(X_train.T,y_train))
+
+ def _lassoFitSKL(self, X_train, y_train, lambda_):
+ """
+ Performs lasso fit using scikit-learn functionality.
+
+ Arguments:
+ X_train: Covariate matrix of the train data set, design matrix of
+ the shape m x p (m_train_datapoints, p_parameters).
+ y_train: Outcome variable, 1D numpy array, dimension m x 1
+ lambda_: Integer type. The shrinkage parameter
+
+ Returns:
+ self.beta : numpy.array
+ The beta parameters from the performed fit
+ """
+ self.regression = Lasso(fit_intercept=True, max_iter=1000000, alpha=self.lambda_)
+ self.regression.fit(X_train,y_train)
+ self.beta = self.regression.coef_
+ self.beta[0] = self.regression.intercept_
+
+ def predict(self, X_test):
+ """
+ Performs prediction of the fitted model on the provided test data set.
+
+ Arguments:
+ X_test: Design matrix, covariate matrix, dimension k x p (k_test_rows, p_columns)
+
+ Returns: self.yHat
+ numpy 1D array, prediction values of dimension k x p
+ """
+ self.X_test = X_test
+ self._predictOwntest(X_test)
+ return self.yHat
+
+ def _predictOwntest(self, X_test):
+ """
+ Performs manual prediction of the given model on the train data.
+ """
+ self.X_test = X_test
+ self.yHat = np.dot(self.X_test, self.beta)
+
+ def MSE(self, y_test):
+ """
+ Calculates the mean squared error (MSE) manually after the fit and prediction have been implemented.
+
+ Arguments:
+ y_test: Outcome variable, 1D numpy array, dimension k x 1 (k_test_rows, 1_column)
+
+ Returns: self._MSE
+ The mean squared error of the predicted model
+ """
+ self.y_test = y_test
+ if self.yHat is None :
+ self._predictOwntest(X_test)
+ N = self.yHat.size
+ self._MSE = (np.sum((self.y_test - self.yHat)**2))/N
+ return self._MSE
+
+ def R2(self, y_test):
+ """
+ Calculates R2 score manually after the fit and prediction have been implemented.
+
+ Arguments:
+ y_test: Outcome variable, 1D numpy array, dimension k x 1 (k_test_rows, 1_column)
+
+ Returns: self._R2
+ The R2 score of the predicted model
+ """
+ self.y_test = y_test
+ if self.yHat is None:
+ self._predictOwntest(X_test)
+ yMean = (1.0 / self.y_test.size) * np.sum(self.y_test)
+ self._R2 = 1.0 - np.sum((self.y_test - self.yHat)**2) / np.sum((self.y_test - yMean)**2)
+ return self._R2
+
+ def CI(self, y_test):
+ """
+ Calculates confidence intervals manually after the fit and prediction have been implemented.
+
+ Arguments:
+ y_test: Outcome variable, 1D numpy array, dimension k x 1 (k_test_rows, 1_column)
+
+ Returns: var, Lower, Upper
+ Variance, Lower and Upper bounds of the confidence intervals for the parameter self.beta
+ """
+ self.y_test = y_test
+ if self.yHat is None:
+ self._predictOwntest(X_test)
+ sigma2 = np.sum(((self.y_test - self.yHat)**2))/(self.y_test.size - self.beta.size)
+ var = np.diag(np.linalg.inv(np.dot(self.X_test.T, self.X_test))) * sigma2
+ Lower = self.beta - 1.96*np.sqrt(var)
+ Upper = self.beta + 1.96*np.sqrt(var)
+ return var, Lower, Upper
+
+
+ ###Implementation through scikitlearn
+class linregSKL:
+ def __init__(self, method = 'ols'):
+ """
+ A class of linear regressions. Perform ordinarly least squares (OLS) and Ridge and Lasso
+ using scikit-learn functionality.
+
+ """
+ self.method = method
+ self.yHat = None
+ self.X = None
+ self.y = None
+ self.beta = None
+
+ self._MSE = None
+ self._R2 = None
+ self._betaVariance = None
+
+
+ def fit(self, X_train, y_train, lambda_ = 0):
+ self.X_train = X_train
+ self.y_train = y_train
+ if self.method == 'ols':
+ self._olsSKLfit(X_train, y_train)
+ if self.method == 'ridge':
+ self._sklRidgeFit(X_train, y_train, lambda_)
+ if self.method == 'lasso':
+ self._SKLlassoFit(X_train, y_train, lambda_)
+ return self.beta
+
+ def _olsSKLfit(self, X_train, y_train):
+ self.X_train = X_train
+ self.y_train = y_train
+ ##We already have standardized data from design matrix
+ self.ols = LinearRegression().fit(self.X_train, self.y_train)
+ self.beta = self.ols.coef_
+ self.beta[0] = self.ols.intercept_
+
+ def _SKLlassoFit(self, X_train, y_train, lambda_):
+ self.regression = Lasso(fit_intercept=True, max_iter=100000, alpha=self.lambda_)
+ self.regression.fit(X_train,y_train)
+ self.beta = self.regression.coef_
+ self.beta[0] = self.regression.intercept_
+
+ def _sklRidgeFit(self, X_train, y_train, lambda_):
+ self.regression = Ridge(fit_intercept=True, alpha=self.lambda_)
+ self.regression.fit(X,y)
+ self.beta = self.regression.coef_
+ self.beta[0] = self.regression.intercept_
+
+ def predict(self, X_test):
+ self.X_test = X_test
+ if self.method == 'ols':
+ self._sklPredict(X_test)
+ return self.yHat
+
+ def _sklPredict(self, X_test):
+ self.X_test = X_test
+ ## Since our data contains 1-s, we should subtract intercept, since scikit learn additionally
+ ##generates the 1-s
+ self.yHat = self.ols.predict(self.X_test) - self.beta[0]
+
+ def MSE(self, y_test):
+ self.y_test = y_test
+ if self.yHat is None :
+ self._sklPredict(X_test)
+ self._MSE = mean_squared_error(self.y_test, self.yHat)
+ return self._MSE
+
+ def R2(self, y_test):
+ self.y_test = y_test
+ if self.yHat is None :
+ self._sklPredict()
+ self._R2 = r2_score(self.y_test, self.yHat)
+ return self._R2
+
+
+
+def designMatrix(x, y, k=5):
+ """
+ Generates the design matrix (covariates of polynomial degree k).
+ Intercept is included in the design matrix.
+ Scaling does not apply to the intercept term.
+ if k = 2, generated column vectors: 1, x, y, x^2, xy, y^2
+ if k = 3, generated column vectors: 1, x, y, x^2, xy, y^2, x^3, x^2y, xy^2, y^3
+ ...
+
+ Arguments:
+ x: 1D numpy array
+ y: 1D numpy array
+ k: integer type. complexity parameter (i.e polynomial degree)
+ """
+
+ xb = np.ones((x.size, 1))
+
+ for i in range(1, k+1):
+ for j in range(i+1):
+ xb = np.c_[xb, (x**(i-j))*(y**j)]
+
+ xb[:, 0] = 1
+ return xb
+
+
+# # Functions for Project 2
+
+# In[7]:
+
+
+# MLP
+
+class Layer:
+ """
+ Represents a layer (hidden or output) in our neural network.
+ """
+
+ def __init__(self, n_input, n_neurons, activation=None, alpha = 0.01, lam = 0.1):
+ """
+ :param int n_input: The input size (coming from the input layer or a previous hidden layer)
+ :param int n_neurons: The number of neurons in this layer.
+ :param str activation: The activation function to use (if any).
+ :param weights: The layer's weights.
+ :param bias: The layer's bias.
+ """
+ self.activation = activation
+ self.alpha = alpha
+ self.lam = lam
+ self.last_activation = None
+ self.error = None
+ self.delta = None
+ # self.weights = weights if weights is not None else np.random.rand(n_input, n_neurons)
+ # self.activation = activation
+ # self.bias = bias if bias is not None else np.random.rand(n_neurons)
+
+ np.random.seed(1921)
+ # Xavier initializations (http://proceedings.mlr.press/v9/glorot10a/glorot10a.pdf).
+ if self.activation == 'sigmoid' :
+ r_inputs = np.sqrt(6.0 / (n_input + n_neurons))
+ self.weights = np.random.uniform(-r_inputs, r_inputs, size=(n_input, n_neurons))
+ #self.weights = np.random.rand(n_input, n_neurons)
+ self.bias = np.random.rand(n_neurons)
+ #self.v_reg = np.random.rand(n_neurons, 1)
+ self.v_reg = np.random.uniform(-(6/(n_neurons+1)), (6/(n_neurons+1)), size=(n_neurons, 1))
+ self.bias_reg = np.random.rand(1)
+
+ elif self.activation == 'tanh' :
+ r_inputs = 4.0 * np.sqrt(6.0 / (n_input + n_neurons))
+ self.weights = np.random.uniform(-r_inputs, r_inputs, size=(n_input, n_neurons))
+ #self.weights = np.random.rand(n_input, n_neurons)
+ self.v_reg = np.random.uniform(-4.0 * np.sqrt(6.0 / (n_neurons+1)), 4.0 * np.sqrt(6.0 / (n_neurons+1)), size=(n_neurons, 1))
+ #self.v_reg = np.random.rand(n_neurons, 1)
+ #self.bias = np.zeros(shape=(n_neurons))
+ self.bias = np.random.rand(n_neurons)
+ self.bias_reg = np.random.rand(1) ##output bias
+
+ # He initializations (https://arxiv.org/pdf/1502.01852.pdf).
+ elif self.activation == 'relu' or self.activation == 'leaky_relu' or self.activation == 'elu' :
+ self.weights = np.random.normal(size=(n_input, n_neurons)) * np.sqrt(2.0 / n_input)
+ #self.weights = np.random.rand(n_input, n_neurons)
+ #self.v_reg = np.random.rand(n_neurons, 1)
+ self.v_reg = np.random.normal(size=(n_neurons,1)) * np.sqrt(2.0 / n_neurons)
+ self.bias = np.random.rand(n_neurons)
+ self.bias_reg = np.random.rand(1)
+
+ else :
+ self.weights = np.random.normal(size=(n_input, n_neurons))
+ #self.weights = np.random.rand(n_input, n_neurons)
+ #self.v_reg = np.random.rand(n_neurons, 1)
+ self.v_reg = np.random.normal(size=(n_neurons, 1))
+ self.bias = np.random.rand(n_neurons)
+ self.bias_reg = np.random.rand(1)
+
+ def activate(self, x):
+ """
+ Calculates the dot product of this layer.
+ :param x: The input.
+ :return: The result.
+ """
+
+ r = np.dot(x, self.weights) + self.bias
+ self.last_activation = self._apply_activation(r)
+ return self.last_activation
+
+ def _apply_activation(self, r):
+ """
+ Applies the chosen activation function (if any).
+ :param r: The normal value.
+ :return: The activated value.
+ """
+
+ # In case no activation function was chosen
+
+ if self.activation is None:
+ return r
+
+ if self.activation == 'tanh':
+ return np.tanh(r)
+
+ if self.activation == 'sigmoid':
+ return self._sigmoid(r)
+
+ if self.activation == 'relu':
+ return self._relu(r)
+
+ if self.activation == 'leaky_relu':
+ return self._leakyrelu(r)
+
+ if self.activation == 'elu':
+ return self._elu(r)
+
+ if self.activation == 'softmax':
+ return self._softmax(r)
+
+ if self.activation == 'identity':
+ return self._identity(r)
+
+ return r
+
+ def apply_activation_derivative(self, r):
+ """
+ Applies the derivative of the activation function (if any).
+ :param r: The normal value.
+ :return: The "derived" value.
+ """
+
+ # We use 'r' directly here because its already activated, the only values that
+
+ # are used in this function are the last activations that were saved.
+
+ if self.activation is None:
+ return r
+
+ if self.activation == 'tanh':
+ return 1 - r ** 2
+
+ if self.activation == 'sigmoid':
+ return r * (1 - r)
+
+ if self.activation == 'relu':
+
+ r[r > 0] = self.lam
+ return r
+
+ if self.activation == 'leaky_relu':
+
+ r[r > 0] = self.lam
+ r[r <= 0] = self.lam * self.alpha
+ return r
+
+ if self.activation == 'identity':
+ return 1
+
+ if self.activation == 'elu':
+
+ r[r > 0] = 1
+ r[r <= 0] = r[r <= 0] + self.alpha
+
+
+ return r
+
+
+ def _sigmoid(self, x):
+ return 1.0/ (1.0 + np.exp(-x))
+
+
+ def _tanh(self, x):
+
+ return np.tanh(x)
+
+
+ def _relu(self, x):
+
+ x = self.lam * x
+ x[x <= 0] = 0
+
+ return x
+
+
+ def _leakyrelu(self, x) :
+
+ x = self.lam * x
+ x[x <= 0] = self.alpha * x[x <= 0]
+
+ return x
+
+ def _identity(self, x) : ##linear activation function
+ return x
+
+
+ def _elu(self, x) :
+ neg = x<0.0
+ x[neg] = self.alpha * (np.exp(x[neg]) - 1.0)
+
+ return x
+
+ def _softmax(self, x) :
+ exps = np.exp(x - np.max(x))
+ return exps / np.sum(exps, axis=0, keepdims=True) ## sum along the column
+
+
+
+class NeuralNetwork:
+ """
+ Represents a neural network.
+ """
+
+ def __init__(self):
+ self._layers = []
+
+ def add_layer(self, layer):
+ """
+ Adds a layer to the neural network.
+ :param Layer layer: The layer to add.
+ """
+
+ self._layers.append(layer)
+
+
+ def feed_forward(self, X):
+ """
+ Feed forward the input through the layers.
+ :param X: The input values.
+ :return: The result.
+ """
+
+ for layer in self._layers:
+ X = layer.activate(X)
+
+ return X
+
+ """
+ N.B: Having a sigmoid activation in the output layer can be interpreted
+ as expecting probabilities as outputs.
+ W'll need to choose a winning class, this is usually done by choosing the
+ index of the biggest probability.
+ """
+ def predict(self, X, net_type = 'regression', n_neurons = 3):
+ """
+ Predicts a class (or classes).
+ :param X: The input values.
+ :return: The predictions.
+ """
+
+ ff = self.feed_forward(X)
+
+ if net_type == 'classification':
+
+ # One row
+
+ if ff.ndim == 1:
+ pred = np.argmax(ff)
+ else:
+ pred = np.argmax(ff, axis = 1)
+
+ if net_type == 'regression':
+
+ pred = ff
+
+ return pred
+
+
+ # # Multiple rows
+
+ # return np.argmax(ff, axis=1)
+
+
+ def backpropagation(self, X, y, learning_rate, lmbd, net_type = 'classification'):
+ """
+ Performs the backward propagation algorithm and updates the layers weights.
+ :param X: The input values.
+ :param y: The target values.
+ :param float learning_rate: The learning rate (between 0 and 1).
+ """
+ ntarget = y.size
+
+ # Feed forward for the output
+
+ output = self.feed_forward(X)
+
+ # Loop over the layers backward
+
+ for i in reversed(range(len(self._layers))):
+ layer = self._layers[i]
+
+ # If this is the output layer
+ if layer == self._layers[-1]:
+ layer.error = y - output
+
+ # The output = layer.last_activation in this case
+ layer.delta = layer.error * layer.apply_activation_derivative(output)
+
+ else:
+ next_layer = self._layers[i + 1]
+ layer.error = np.dot(next_layer.weights, next_layer.delta)
+ layer.delta = layer.error * layer.apply_activation_derivative(layer.last_activation)
+
+
+ # Update the weights
+
+ for i in range(len(self._layers)):
+ layer = self._layers[i]
+ # The input is either the previous layers output or X itself (for the first hidden layer)
+
+ input_to_use = np.atleast_2d(X if i == 0 else self._layers[i - 1].last_activation)
+
+ layer.weights = layer.weights + layer.delta * input_to_use.T * learning_rate
+
+ if lmbd > 0: ###adding L2 regularization
+
+ layer.weights = layer.weights*(1-lmbd*learning_rate) + layer.delta * input_to_use.T * learning_rate
+
+ layer.bias = layer.bias + layer.delta * learning_rate
+
+
+
+ def train(self, X, y, learning_rate, max_epochs, net_type = 'regression', lmbd = 0):
+ """
+ Trains the neural network using backpropagation.
+ :param X: The input values.
+ :param y: The target values.
+ :param float learning_rate: The learning rate (between 0 and 1).
+ :param int max_epochs: The maximum number of epochs (cycles).
+ :return: The list of calculated MSE errors.
+ """
+
+ mses = []
+
+ for i in range(max_epochs):
+ for j in range(len(X)): ##len(X) is rows
+ self.backpropagation(X[j], y[j], learning_rate, lmbd)
+
+ # if i % 10 == 0: #At every 10th epoch, we will print out the Mean Squared Error and save it in mses which we will return at the end.
+ # nn = NeuralNetwork()
+ # mse = np.mean(np.square(y - nn.feed_forward(X)))
+ # mses.append(mse)
+ # print('Epoch: #%s, MSE: %f' % (i, float(mse)))
+ # return mses
+
+
+
+ def MSE(self, y_pred, y_true):
+ return (1/len(y_true))*np.sum((y_pred - y_true)**2)
+
+
+ def accuracy(self, y_pred, y_true):
+ """
+ Calculates the accuracy between the predicted labels and true labels.
+ :param y_pred: The predicted labels.
+ :param y_true: The true labels.
+ :return: The calculated accuracy.
+ """
+
+ return (y_pred == y_true).mean()
+
+ def cal_err(self, y_pred, y_true, costf):
+
+ if costf == "squared-error":
+ err = np.sum((y_pred - y_true)**2)
+ elif costf == "MSE":
+ err = (1/len(y_true))*np.sum((y_pred - y_true)**2)
+
+ return err
+
+ def confusion_table(self, y_pred, y_true):
+
+ conf = confusion_matrix(y_true, y_pred, labels=[0, 1])
+
+ return conf
+
+ def cal_r2(self, y_pred, y_true):
+
+ mu = np.mean(y_true)
+ SS_tot = np.sum((y_true - mu)**2)
+ SS_res = np.sum((y_true - y_pred)**2)
+
+ r2 = 1 - (SS_res/SS_tot)
+
+ return r2
+
+
+# In[8]:
+
+
+# Multiclass Logistig Regression
+
+class multiclassLogistic:
+
+ def __init__(self, X, y, y_onehot, learning_rate, lambda_):
+ self.X = X
+ self.y = y
+ self.learning_rate = learning_rate
+ self.lambda_ = lambda_
+ self.y_onehot = y_onehot
+
+
+ def softmax(self, X) :
+ exps = np.exp(X - np.max(X))
+ return exps / np.sum(exps, axis=1, keepdims=True)
+
+
+ def sgd(self, X, y, y_onehot, iterations, lambda_, learning_rate):
+ w = np.zeros([X.shape[1],len(np.unique(y))])
+ losses = []
+
+ for i in range(0,iterations):
+ loss,grad = self._getLoss(w,X,y_onehot, lambda_)
+ losses.append(loss)
+ w = w - (learning_rate * grad)
+ return w
+
+
+ def _getLoss(self, w, X, y_onehot, lambda_):
+ m = X.shape[0] #First we get the number of training examples
+ scores = np.dot(X,w) #Then we compute raw class scores given our input and current weights
+ prob = self.softmax(scores) #Next we perform a softmax on these scores to get their probabilities
+ loss = (-1 / m) * np.sum(y_onehot * np.log(prob)) + (lambda_/2)*np.sum(w*w) #We then find the loss of the probabilities
+ grad = (-1 / m) * np.dot(X.T,(y_onehot - prob)) + lambda_*w #And compute the gradient for that loss
+ return loss,grad
+
+
+ def accuracy(self, someX,someY, w):
+ prob,prede = self._getProbsAndPreds(someX, w)
+ accuracy = sum(prede == someY)/(float(len(someY)))
+ return accuracy
+
+
+ def _getProbsAndPreds(self, someX, w):
+ probs = self.softmax(np.dot(someX,w))
+ preds = np.argmax(probs,axis=1) ##returns the highest value index (i.e. class) along the row
+ return probs,preds
+
+
+# In[9]:
+
+
+# Stochastic Gradient Descent
+
+from matplotlib.ticker import LinearLocator, FormatStrFormatter
+
+def compute_square_loss(X, y, theta):
+ loss = 0 #Initialize the average square loss
+
+ m = len(y)
+ loss = (1.0/m)*(np.linalg.norm((X.dot(theta) - y)) ** 2)
+ return loss
+
+
+def gradient_ridge(X, y, beta, lambda_):
+ return 2*(np.dot(X.T, (X.dot(beta) - y))) + 2*lambda_*beta
+
+def gradient_ols(X, y, beta):
+ m = X.shape[0]
+
+ grad = 2/m * X.T.dot(X.dot(beta) - y)
+
+ return grad
+
+def learning_schedule(t):
+ t0, t1 = 5, 50
+ return t0/(t+t1)
+
+
+def iterate_minibatches(inputs, targets, batchsize, shuffle=True):
+ assert inputs.shape[0] == targets.shape[0]
+ if shuffle:
+ indices = np.random.permutation(inputs.shape[0])
+ for start_idx in range(0, inputs.shape[0], batchsize):
+ end_idx = min(start_idx + batchsize, inputs.shape[0])
+ if shuffle:
+ excerpt = indices[start_idx:end_idx]
+ else:
+ excerpt = slice(start_idx, end_idx)
+ yield inputs[excerpt], targets[excerpt]
+
+
+###sgd
+def SGD(X, y, learning_rate = 0.02, n_epochs = 100, lambda_ = 0.01, batch_size = 20, method = 'ols'):
+ num_instances, num_features = X.shape[0], X.shape[1]
+ beta = np.random.randn(num_features) ##initialize beta
+
+ for epoch in range(n_epochs+1):
+
+ for batch in iterate_minibatches(X, y, batch_size, shuffle=True):
+
+ X_batch, y_batch = batch
+
+ # for i in range(batch_size):
+ # learning_rate = learning_schedule(n_epochs*epoch + i)
+
+ if method == 'ols':
+ gradient = gradient_ols(X_batch, y_batch, beta)
+ beta = beta - learning_rate*gradient
+ if method == 'ridge':
+ gradient = gradient_ridge(X_batch, y_batch, beta, lambda_ = lambda_)
+ beta = beta - learning_rate*gradient
+
+ mse_ols_train = compute_square_loss(X, y, beta)
+ mse_ridge_train = compute_square_loss(X, y, beta) + lambda_*np.dot(beta.T, beta)
+
+ return beta
+
+def compute_test_mse(X_test, y_test, beta, lambda_ = 0.01):
+ mse_ols_test = compute_square_loss(X_test, y_test, beta)
+ mse_ridge_test = compute_square_loss(X_test, y_test, beta) + lambda_*np.dot(beta.T, beta)
+ return mse_ols_test, mse_ridge_test
+
+
+# # Part A
+
+# In[10]:
+
+
+# a
+
+##Make synthetic data
+n = 1000
+np.random.seed(20)
+x1 = np.random.rand(n)
+x2 = np.random.rand(n)
+X = designMatrix(x1, x2, 4)
+y = franke(x1, x2)
+
+##Train-validation-test samples.
+# We choose / play with hyper-parameters on the validation data and then test predictions on the test data
+
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)
+X_train, X_val, y_train, y_val = train_test_split(X_train, y_train, test_size=0.25, random_state=1) # 0.25 x 0.8 = 0.2
+
+scaler = StandardScaler()
+scaler.fit(X_train)
+X_train = scaler.transform(X_train)
+X_test = scaler.transform(X_test)
+X_val = scaler.transform(X_val)
+
+X_train[:, 0] = 1
+X_test[:, 0] = 1
+X_val[:, 0] = 1
+
+
+
+linreg = linregOwn(method='ols')
+#print('Invert OLS:', linreg.fit(X_train, y_train))
+beta = SGD(X_train, y_train, learning_rate=0.07)
+#print('SGD OLS:', beta)
+
+
+linreg = linregOwn(method='ridge')
+#print('Invert Ridge:', linreg.fit(X_train, y_train, lambda_= 0.01))
+beta = SGD(X_train, y_train, learning_rate=0.0004, method='ridge')
+#print('SGD Ridge:', beta)
+
+
+sgdreg = SGDRegressor(max_iter = 100, penalty=None, eta0=0.1)
+sgdreg.fit(X_train[:, 1:],y_train.ravel())
+#print('sklearn:', sgdreg.coef_)
+#print('sklearn intercept:', sgdreg.intercept_)
+
+
+def plot_MSE(method = 'ridge', scheme = None):
+ eta = np.logspace(-5, -3, 10)
+ lambda_ = np.logspace(-5, -1, 10)
+ MSE_ols = []
+ MSE_ridge = []
+
+ if scheme == 'joint':
+
+ if method == 'ridge':
+
+ for lmbd in lambda_:
+
+ for i in eta:
+
+ beta = SGD(X_train, y_train, learning_rate=i, lambda_ = lmbd, method = method)
+ mse_ols_test, mse_ridge_test = compute_test_mse(X_val, y_val, lambda_ = lmbd, beta = beta)
+ MSE_ridge.append(mse_ridge_test)
+
+ fig = plt.figure()
+ ax = fig.gca(projection='3d') ##get current axis
+ lambda_ = np.ravel(lambda_)
+ eta = np.ravel(eta)
+ ax.zaxis.set_major_locator(LinearLocator(5))
+ ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
+ ax.xaxis.set_major_formatter(FormatStrFormatter('%.02f'))
+ ax.yaxis.set_major_formatter(FormatStrFormatter('%.03f'))
+ ax.plot_trisurf(lambda_, eta, MSE_ridge, cmap='viridis', edgecolor='none')
+ ax.set_xlabel(r'$\lambda$')
+ ax.set_ylabel(r'$\eta$')
+ ax.set_title(r'MSE Ridge')
+ ax.view_init(30, 60)
+ plt.show()
+
+ if scheme == 'separate':
+
+ if method == 'ols':
+
+ eta = np.logspace(-5, 0, 10)
+
+ for i in eta:
+
+ beta = SGD(X_train, y_train, learning_rate=i, lambda_ = 0.01, method = method)
+ mse_ols_test, mse_ridge_test = compute_test_mse(X_val, y_val, beta = beta)
+ MSE_ols.append(mse_ols_test)
+
+ print('The learning rate {} performs best for the OLS' .format(eta[MSE_ols.index(min(MSE_ols))]))
+ print('Corresponding minimum MSE for OLS: {}'.format(min(MSE_ols)))
+ plt.semilogx(eta, MSE_ols)
+ plt.xlabel(r'Learning rate, $\eta$')
+ plt.ylabel('MSE OLS')
+ plt.title('Stochastic Gradient Descent')
+ plt.show()
+ if scheme == 'separate':
+
+ if method == 'ridge':
+
+ eta = np.logspace(-5, 0, 10)
+
+ for i in eta:
+
+ beta = SGD(X_train, y_train, learning_rate=i, lambda_ = 0.01, method = method)
+ mse_ols_test, mse_ridge_test = compute_test_mse(X_val, y_val, beta = beta)
+ MSE_ols.append(mse_ridge_test)
+
+ print('The learning rate {} performs best for Ridge' .format(eta[MSE_ols.index(min(MSE_ols))]))
+ print('Corresponding minimum MSE for Ridge: {}'.format(min(MSE_ols)))
+
+ plt.plot(eta, MSE_ols)
+ plt.xlabel(r'Learning rate, $\eta$')
+ plt.ylabel('MSE Ridge')
+ plt.title('Stochastic Gradient Descent')
+ plt.show()
+
+
+
+
+####Predict OLS, Ridge on test data after tuning learning rate and lambda on validation data
+
+def plot_scatter(y_true, method = 'ols'):
+ if method == 'ols':
+ beta = SGD(X_train, y_train, learning_rate=0.07, lambda_ = 0, method = method, n_epochs=300)
+ if method == 'ridge':
+ beta = SGD(X_train, y_train, learning_rate=0.0001, lambda_ = 0, method = method, n_epochs=300)
+ y_pred = np.dot(X_test, beta)
+ mse_ols_test, mse_ridge_test = compute_test_mse(X_test, y_true, beta = beta)
+ print('Test MSE OLS: {}' .format(mse_ols_test))
+ print('Test MSE Ridge: {}' .format(mse_ridge_test))
+ a = plt.axes(aspect='equal')
+ plt.scatter(y_pred, y_pred, color= 'blue', label = "True values")
+ plt.scatter(y_pred, y_true, color = 'red', label = "Predicted values")
+ plt.xlabel('True y values')
+ plt.ylabel('Predicted y')
+ plt.title(f"Prediction - {method}")
+ plt.legend()
+ # if method == 'ols':
+ # plt.savefig(os.path.join(os.path.dirname(__file__), 'Plots', 'ols_reg_pred.png'), transparent=True, bbox_inches='tight')
+ # if method == 'ridge':
+ # plt.savefig(os.path.join(os.path.dirname(__file__), 'Plots', 'ridge_reg_pred.png'), transparent=True, bbox_inches='tight')
+
+ plt.show()
+
+plot_scatter(y_test, method='ols')
+
+plot_scatter(y_test, method='ridge')
+
+
+
diff --git a/doc/src/week39/codes/test8.py~ b/doc/src/week39/codes/test8.py~
new file mode 100644
index 000000000..7061dcae6
--- /dev/null
+++ b/doc/src/week39/codes/test8.py~
@@ -0,0 +1,77 @@
+"""
+Code to test Ridge and NNs using Scikit-Learn only
+"""
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from sklearn.model_selection import train_test_split
+from sklearn import linear_model
+from sklearn.neural_network import MLPRegressor
+from sklearn.metrics import accuracy_score
+import seaborn as sns
+
+
+def MSE(y_data,y_model):
+ n = np.size(y_model)
+ return np.sum((y_data-y_model)**2)/n
+# A seed just to ensure that the random numbers are the same for every run.
+# Useful for eventual debugging.
+np.random.seed(315)
+
+n = 100
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 5
+X = np.zeros((n,Maxpolydegree-1))
+
+for degree in range(1,Maxpolydegree): #No intercept column
+ X[:,degree-1] = x**(degree)
+
+# We split the data in test and training data
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
+
+# Decide which values of lambda to use
+
+nlambdas = 10
+lmbd_vals = np.logspace(-4, 0, nlambdas)
+MSERidgePredict = np.zeros(nlambdas)
+for i in range(nlambdas):
+ lmb = lmbd_vals[i]
+ RegRidge = linear_model.Ridge(lmb)
+ RegRidge.fit(X_train,y_train)
+ ypredictRidge = RegRidge.predict(X_test)
+ MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+
+plt.figure()
+plt.plot(np.log10(lmbd_vals), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+
+# Neural Network part
+
+n_hidden_neurons = 100
+epochs = 100
+# store models for later use
+eta_vals = np.logspace(-4, 0, 4)
+# store the models for later use
+DNN_scikit = np.zeros((len(eta_vals), len(lmbd_vals)), dtype=object)
+test_accuracy = np.zeros((len(eta_vals), len(lmbd_vals)))
+sns.set()
+for i, eta in enumerate(eta_vals):
+ for j, lmbd in enumerate(lmbd_vals):
+ dnn = MLPRegressor(hidden_layer_sizes=(n_hidden_neurons), activation='logistic',
+ alpha=lmbd, learning_rate_init=eta, max_iter=epochs)
+ dnn.fit(X_train, y_train)
+ ypredictMLP = dnn.predict(X_test)
+ test_accuracy[i][j] = MSE(ypredictMLP, y_test)
+
+fig, ax = plt.subplots(figsize = (10, 10))
+sns.heatmap(test_accuracy, annot=True, ax=ax, cmap="viridis")
+ax.set_title("Training Accuracy")
+ax.set_ylabel("$\eta$")
+ax.set_xlabel("$\lambda$")
+plt.show()
diff --git a/doc/src/week39/codes/test9.py~ b/doc/src/week39/codes/test9.py~
new file mode 100644
index 000000000..e10d5b982
--- /dev/null
+++ b/doc/src/week39/codes/test9.py~
@@ -0,0 +1,85 @@
+"""
+Code to test Ridge with own gradient descent and SGD
+"""
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from sklearn.model_selection import train_test_split
+from sklearn import linear_model
+from sklearn.neural_network import MLPRegressor
+from sklearn.metrics import accuracy_score
+import seaborn as sns
+import autograd.numpy as np
+from autograd import grad
+
+
+def MSE(y_data,y_model):
+ n = np.size(y_model)
+ return np.sum((y_data-y_model)**2)/n
+# A seed just to ensure that the random numbers are the same for every run.
+# Useful for eventual debugging.
+np.random.seed(315)
+
+n = 100
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 5
+X = np.zeros((n,Maxpolydegree-1))
+
+for degree in range(1,Maxpolydegree): #No intercept column
+ X[:,degree-1] = x**(degree)
+
+# We split the data in test and training data
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
+
+
+nlambdas = 10
+lmbd_vals = np.logspace(-4, 0, nlambdas)
+MSERidgePredict = np.zeros(nlambdas)
+for i in range(nlambdas):
+ lmb = lmbd_vals[i]
+ RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
+ RegRidge.fit(X_train,y_train)
+ ypredictRidge = RegRidge.predict(X_test)
+ MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+
+beta = np.random.randn(X_train.shape[1],1)
+loss = np.mean((y_train.reshape(-1,1) - X_train@beta)**2)
+print(loss)
+get_grad = grad(loss,argnum=2)
+print(get_grad)
+#grad_beta = get_grad(X_train,y_train,beta)
+#print(grad_beta)
+
+"""
+print(beta)
+print( (X_train.T @ y_train).T)
+# Make own gradient descent and define precalculated quantities, saves cycles
+XT_X = X_train.T @ X_train
+XTy = X_train.T @ y_train
+MSERidgeGDPredict = np.zeros(nlambdas)
+for i in range(nlambdas):
+ lmb = lmbd_vals[i]
+ Id = lmb* np.eye(XT_X.shape[0])
+ beta = np.random.randn(X_train.shape[1],1)
+ eta = 0.01
+ Niterations = 2
+# beta_linreg = np.linalg.pinv(XT_X+Id) @ X_train.T @ y_train
+ for iter in range(Niterations):
+ XX = XT_X @ beta-XTy
+ gradients = (2.0/n)*XX *lmb*beta
+ beta -= eta*gradients
+ ypredictRidgeGD = X_test @ beta
+ MSERidgeGDPredict[i] = MSE(y_test,ypredictRidgeGD)
+
+plt.figure()
+plt.plot(np.log10(lmbd_vals), MSERidgePredict, 'g--', label = 'MSE Sklearn Ridge Test')
+plt.plot(np.log10(lmbd_vals), MSERidgeGDPredict, 'r', label = 'MSE GD Ridge Test')
+plt.xlabel('log10(lambda)')
+plt.ylabel('MSE')
+plt.legend()
+plt.show()
+"""
+
diff --git a/doc/src/week41/programs/breast.py~ b/doc/src/week41/programs/breast.py~
new file mode 100644
index 000000000..205822b14
--- /dev/null
+++ b/doc/src/week41/programs/breast.py~
@@ -0,0 +1,164 @@
+import tensorflow as tf
+from tensorflow.keras.layers import Input
+from tensorflow.keras.models import Sequential #This allows appending layers to existing models
+from tensorflow.keras.layers import Dense #This allows defining the characteristics of a particular layer
+from tensorflow.keras import optimizers #This allows using whichever optimiser we want (sgd,adam,RMSprop)
+from tensorflow.keras import regularizers #This allows using whichever regularizer we want (l1,l2,l1_l2)
+from tensorflow.keras.utils import to_categorical #This allows using categorical cross entropy as the cost function
+import numpy as np
+import matplotlib.pyplot as plt
+import seaborn as sns
+from sklearn.model_selection import train_test_split as splitter
+from sklearn.datasets import load_breast_cancer
+import pickle
+import os
+
+
+"""Load breast cancer dataset"""
+
+np.random.seed(0) #create same seed for random number every time
+
+cancer=load_breast_cancer() #Download breast cancer dataset
+
+inputs=cancer.data #Feature matrix of 569 rows (samples) and 30 columns (parameters)
+outputs=cancer.target #Label array of 569 rows (0 for benign and 1 for malignant)
+labels=cancer.feature_names[0:30]
+
+print('The content of the breast cancer dataset is:') #Print information about the datasets
+print(labels)
+print('-------------------------')
+print("inputs = " + str(inputs.shape))
+print("outputs = " + str(outputs.shape))
+print("labels = "+ str(labels.shape))
+
+x=inputs #Reassign the Feature and Label matrices to other variables
+y=outputs
+
+# Visualisation of dataset (for correlation analysis)
+
+plt.figure()
+plt.scatter(x[:,0],x[:,2],s=40,c=y,cmap=plt.cm.Spectral)
+plt.xlabel('Mean radius',fontweight='bold')
+plt.ylabel('Mean perimeter',fontweight='bold')
+plt.show()
+
+plt.figure()
+plt.scatter(x[:,5],x[:,6],s=40,c=y, cmap=plt.cm.Spectral)
+plt.xlabel('Mean compactness',fontweight='bold')
+plt.ylabel('Mean concavity',fontweight='bold')
+plt.show()
+
+
+plt.figure()
+plt.scatter(x[:,0],x[:,1],s=40,c=y,cmap=plt.cm.Spectral)
+plt.xlabel('Mean radius',fontweight='bold')
+plt.ylabel('Mean texture',fontweight='bold')
+plt.show()
+
+plt.figure()
+plt.scatter(x[:,2],x[:,1],s=40,c=y,cmap=plt.cm.Spectral)
+plt.xlabel('Mean perimeter',fontweight='bold')
+plt.ylabel('Mean compactness',fontweight='bold')
+plt.show()
+
+
+# Generate training and testing datasets
+
+#Select features relevant to classification (texture,perimeter,compactness and symmetery)
+#and add to input matrix
+
+temp1=np.reshape(x[:,1],(len(x[:,1]),1))
+temp2=np.reshape(x[:,2],(len(x[:,2]),1))
+X=np.hstack((temp1,temp2))
+temp=np.reshape(x[:,5],(len(x[:,5]),1))
+X=np.hstack((X,temp))
+temp=np.reshape(x[:,8],(len(x[:,8]),1))
+X=np.hstack((X,temp))
+
+X_train,X_test,y_train,y_test=splitter(X,y,test_size=0.1) #Split datasets into training and testing
+
+y_train=to_categorical(y_train) #Convert labels to categorical when using categorical cross entropy
+y_test=to_categorical(y_test)
+
+del temp1,temp2,temp
+
+# Define tunable parameters"
+
+eta=np.logspace(-3,-1,3) #Define vector of learning rates (parameter to SGD optimiser)
+lamda=0.01 #Define hyperparameter
+n_layers=2 #Define number of hidden layers in the model
+n_neuron=np.logspace(0,3,4,dtype=int) #Define number of neurons per layer
+epochs=100 #Number of reiterations over the input data
+batch_size=100 #Number of samples per gradient update
+
+"""Define function to return Deep Neural Network model"""
+
+def NN_model(inputsize,n_layers,n_neuron,eta,lamda):
+ model=Sequential()
+ for i in range(n_layers): #Run loop to add hidden layers to the model
+ if (i==0): #First layer requires input dimensions
+ model.add(Dense(n_neuron,activation='relu',kernel_regularizer=regularizers.l2(lamda),input_dim=inputsize))
+ else: #Subsequent layers are capable of automatic shape inferencing
+ model.add(Dense(n_neuron,activation='relu',kernel_regularizer=regularizers.l2(lamda)))
+ model.add(Dense(2,activation='softmax')) #2 outputs - ordered and disordered (softmax for prob)
+ sgd=optimizers.SGD(lr=eta)
+ model.compile(loss='categorical_crossentropy',optimizer=sgd,metrics=['accuracy'])
+ return model
+
+
+Train_accuracy=np.zeros((len(n_neuron),len(eta))) #Define matrices to store accuracy scores as a function
+Test_accuracy=np.zeros((len(n_neuron),len(eta))) #of learning rate and number of hidden neurons for
+
+for i in range(len(n_neuron)): #run loops over hidden neurons and learning rates to calculate
+ for j in range(len(eta)): #accuracy scores
+ DNN_model=NN_model(X_train.shape[1],n_layers,n_neuron[i],eta[j],lamda)
+ DNN_model.fit(X_train,y_train,epochs=epochs,batch_size=batch_size,verbose=1)
+ Train_accuracy[i,j]=DNN_model.evaluate(X_train,y_train)[1]
+ Test_accuracy[i,j]=DNN_model.evaluate(X_test,y_test)[1]
+
+
+def plot_data(x,y,data,title=None):
+
+ # plot results
+ fontsize=16
+
+
+ fig = plt.figure()
+ ax = fig.add_subplot(111)
+ cax = ax.matshow(data, interpolation='nearest', vmin=0, vmax=1)
+
+ cbar=fig.colorbar(cax)
+ cbar.ax.set_ylabel('accuracy (%)',rotation=90,fontsize=fontsize)
+ cbar.set_ticks([0,.2,.4,0.6,0.8,1.0])
+ cbar.set_ticklabels(['0%','20%','40%','60%','80%','100%'])
+
+ # put text on matrix elements
+ for i, x_val in enumerate(np.arange(len(x))):
+ for j, y_val in enumerate(np.arange(len(y))):
+ c = "${0:.1f}\\%$".format( 100*data[j,i])
+ ax.text(x_val, y_val, c, va='center', ha='center')
+
+ # convert axis vaues to to string labels
+ x=[str(i) for i in x]
+ y=[str(i) for i in y]
+
+
+ ax.set_xticklabels(['']+x)
+ ax.set_yticklabels(['']+y)
+
+ ax.set_xlabel('$\\mathrm{learning\\ rate}$',fontsize=fontsize)
+ ax.set_ylabel('$\\mathrm{hidden\\ neurons}$',fontsize=fontsize)
+ if title is not None:
+ ax.set_title(title)
+
+ plt.tight_layout()
+
+ plt.show()
+
+plot_data(eta,n_neuron,Train_accuracy, 'training')
+plot_data(eta,n_neuron,Test_accuracy, 'testing')
+
+
+
+
+
diff --git a/doc/src/week41/programs/nnown.py~ b/doc/src/week41/programs/nnown.py~
new file mode 100644
index 000000000..7875ade5f
--- /dev/null
+++ b/doc/src/week41/programs/nnown.py~
@@ -0,0 +1,171 @@
+
+# import necessary packages
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.neural_network import MLPClassifier
+from sklearn.metrics import accuracy_score
+import seaborn as sns
+
+# ensure the same random numbers appear every time
+np.random.seed(0)
+
+# Design matrix
+X = np.array([ [0, 0], [0, 1], [1, 0],[1, 1]],dtype=np.float64)
+
+# The XOR gate
+yXOR = np.array( [ 0, 1 ,1, 0])
+# The OR gate
+yOR = np.array( [ 0, 1 ,1, 1])
+# The AND gate
+yAND = np.array( [ 0, 0 ,0, 1])
+
+
+# Defining the neural network
+n_inputs, n_features = X.shape
+n_hidden_neurons = 2
+n_categories = 2
+n_features = 2
+
+def sigmoid(x):
+ return 1/(1 + np.exp(-x))
+
+
+class NeuralNetwork:
+ def __init__(
+ self,
+ X_data,
+ Y_data,
+ n_hidden_neurons=2,
+ n_categories=2,
+ epochs=10,
+ batch_size=100,
+ eta=0.1,
+ lmbd=0.0):
+
+ self.X_data_full = X_data
+ self.Y_data_full = Y_data
+
+ self.n_inputs = X_data.shape[0]
+ self.n_features = X_data.shape[1]
+ self.n_hidden_neurons = n_hidden_neurons
+ self.n_categories = n_categories
+
+ self.epochs = epochs
+ self.batch_size = batch_size
+ self.iterations = self.n_inputs // self.batch_size
+ self.eta = eta
+ self.lmbd = lmbd
+
+ self.create_biases_and_weights()
+
+
+ def create_biases_and_weights(self):
+ self.hidden_weights = np.random.randn(self.n_features, self.n_hidden_neurons)
+ self.hidden_bias = np.zeros(self.n_hidden_neurons) + 0.01
+
+ self.output_weights = np.random.randn(self.n_hidden_neurons, self.n_categories)
+ self.output_bias = np.zeros(self.n_categories) + 0.01
+
+ def feed_forward(self):
+ # feed-forward for training
+ self.z_h = np.matmul(self.X_data, self.hidden_weights) + self.hidden_bias
+ self.a_h = sigmoid(self.z_h)
+
+ self.z_o = np.matmul(self.a_h, self.output_weights) + self.output_bias
+
+ exp_term = np.exp(self.z_o)
+ self.probabilities = exp_term / np.sum(exp_term, axis=1, keepdims=True)
+
+ def feed_forward_out(self, X):
+ # feed-forward for output
+ z_h = np.matmul(X, self.hidden_weights) + self.hidden_bias
+ a_h = sigmoid(z_h)
+
+ z_o = np.matmul(a_h, self.output_weights) + self.output_bias
+
+ exp_term = np.exp(z_o)
+ probabilities = exp_term / np.sum(exp_term, axis=1, keepdims=True)
+ return probabilities
+
+ def backpropagation(self):
+ error_output = self.probabilities - self.Y_data
+ error_hidden = np.matmul(error_output, self.output_weights.T) * self.a_h * (1 - self.a_h)
+
+ self.output_weights_gradient = np.matmul(self.a_h.T, error_output)
+ self.output_bias_gradient = np.sum(error_output, axis=0)
+
+ self.hidden_weights_gradient = np.matmul(self.X_data.T, error_hidden)
+ self.hidden_bias_gradient = np.sum(error_hidden, axis=0)
+
+ if self.lmbd > 0.0:
+ self.output_weights_gradient += self.lmbd * self.output_weights
+ self.hidden_weights_gradient += self.lmbd * self.hidden_weights
+
+ self.output_weights -= self.eta * self.output_weights_gradient
+ self.output_bias -= self.eta * self.output_bias_gradient
+ self.hidden_weights -= self.eta * self.hidden_weights_gradient
+ self.hidden_bias -= self.eta * self.hidden_bias_gradient
+
+ def predict(self, X):
+ probabilities = self.feed_forward_out(X)
+ return np.argmax(probabilities, axis=1)
+
+ def predict_probabilities(self, X):
+ probabilities = self.feed_forward_out(X)
+ return probabilities
+
+ def train(self):
+ data_indices = np.arange(self.n_inputs)
+
+ for i in range(self.epochs):
+ for j in range(self.iterations):
+ # pick datapoints with replacement
+ chosen_datapoints = np.random.choice(
+ data_indices, size=self.batch_size, replace=False
+ )
+
+ # minibatch training data
+ self.X_data = self.X_data_full[chosen_datapoints]
+ self.Y_data = self.Y_data_full[chosen_datapoints]
+
+ self.feed_forward()
+ self.backpropagation()
+
+epochs = 100
+batch_size = 100
+
+eta_vals = np.logspace(-5, 1, 7)
+lmbd_vals = np.logspace(-5, 1, 7)
+# store the models for later use
+DNN_numpy = np.zeros((len(eta_vals), len(lmbd_vals)), dtype=object)
+
+# grid search
+for i, eta in enumerate(eta_vals):
+ for j, lmbd in enumerate(lmbd_vals):
+ dnn = NeuralNetwork(X, yXOR, eta=eta, lmbd=lmbd, epochs=epochs, batch_size=batch_size,
+ n_hidden_neurons=n_hidden_neurons, n_categories=n_categories)
+ dnn.train()
+ DNN_numpy[i][j] = dnn
+ test_predict = dnn.predict(X)
+ print("Learning rate = ", eta)
+ print("Lambda = ", lmbd)
+ print("Accuracy score on test set: ", accuracy_score(yXOR, test_predict))
+ print()
+
+sns.set()
+test_accuracy = np.zeros((len(eta_vals), len(lmbd_vals)))
+
+for i in range(len(eta_vals)):
+ for j in range(len(lmbd_vals)):
+ dnn = DNN_numpy[i][j]
+ test_pred = dnn.predict(X)
+ test_accuracy[i][j] = accuracy_score(yXOR, test_pred)
+
+fig, ax = plt.subplots(figsize = (10, 10))
+sns.heatmap(test_accuracy, annot=True, ax=ax, cmap="viridis")
+ax.set_title("Test Accuracy")
+ax.set_ylabel("$\eta$")
+ax.set_xlabel("$\lambda$")
+plt.show()
+
+
diff --git a/doc/src/week41/programs/test.py~ b/doc/src/week41/programs/test.py~
new file mode 100644
index 000000000..c150131e3
--- /dev/null
+++ b/doc/src/week41/programs/test.py~
@@ -0,0 +1,139 @@
+# import necessary packages
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn import datasets
+
+
+# ensure the same random numbers appear every time
+np.random.seed(0)
+
+# display images in notebook
+%matplotlib inline
+plt.rcParams['figure.figsize'] = (12,12)
+
+
+# download MNIST dataset
+digits = datasets.load_digits()
+
+# define inputs and labels
+inputs = digits.images
+labels = digits.target
+
+# RGB images have a depth of 3
+# our images are grayscale so they should have a depth of 1
+inputs = inputs[:,:,:,np.newaxis]
+
+print("inputs = (n_inputs, pixel_width, pixel_height, depth) = " + str(inputs.shape))
+print("labels = (n_inputs) = " + str(labels.shape))
+
+
+# choose some random images to display
+n_inputs = len(inputs)
+indices = np.arange(n_inputs)
+random_indices = np.random.choice(indices, size=5)
+
+for i, image in enumerate(digits.images[random_indices]):
+ plt.subplot(1, 5, i+1)
+ plt.axis('off')
+ plt.imshow(image, cmap=plt.cm.gray_r, interpolation='nearest')
+ plt.title("Label: %d" % digits.target[random_indices[i]])
+plt.show()
+
+from tensorflow.keras import datasets, layers, models
+from tensorflow.keras.layers import Input
+from tensorflow.keras.models import Sequential #This allows appending layers to existing models
+from tensorflow.keras.layers import Dense #This allows defining the characteristics of a particular layer
+from tensorflow.keras import optimizers #This allows using whichever optimiser we want (sgd,adam,RMSprop)
+from tensorflow.keras import regularizers #This allows using whichever regularizer we want (l1,l2,l1_l2)
+from tensorflow.keras.utils import to_categorical #This allows using categorical cross entropy as the cost function
+#from tensorflow.keras import Conv2D
+#from tensorflow.keras import MaxPooling2D
+#from tensorflow.keras import Flatten
+
+from sklearn.model_selection import train_test_split
+
+# representation of labels
+labels = to_categorical(labels)
+
+# split into train and test data
+# one-liner from scikit-learn library
+train_size = 0.8
+test_size = 1 - train_size
+X_train, X_test, Y_train, Y_test = train_test_split(inputs, labels, train_size=train_size,
+ test_size=test_size)
+
+def create_convolutional_neural_network_keras(input_shape, receptive_field,
+ n_filters, n_neurons_connected, n_categories,
+ eta, lmbd):
+ model = Sequential()
+ model.add(layers.Conv2D(n_filters, (receptive_field, receptive_field), input_shape=input_shape, padding='same',
+ activation='relu', kernel_regularizer=regularizers.l2(lmbd)))
+ model.add(layers.MaxPooling2D(pool_size=(2, 2)))
+ model.add(layers.Flatten())
+ model.add(layers.Dense(n_neurons_connected, activation='relu', kernel_regularizer=regularizers.l2(lmbd)))
+ model.add(layers.Dense(n_categories, activation='softmax', kernel_regularizer=regularizers.l2(lmbd)))
+
+ sgd = optimizers.SGD(lr=eta)
+ model.compile(loss='categorical_crossentropy', optimizer=sgd, metrics=['accuracy'])
+
+ return model
+
+epochs = 100
+batch_size = 100
+input_shape = X_train.shape[1:4]
+receptive_field = 3
+n_filters = 10
+n_neurons_connected = 50
+n_categories = 10
+
+eta_vals = np.logspace(-5, 1, 7)
+lmbd_vals = np.logspace(-5, 1, 7)
+
+CNN_keras = np.zeros((len(eta_vals), len(lmbd_vals)), dtype=object)
+
+for i, eta in enumerate(eta_vals):
+ for j, lmbd in enumerate(lmbd_vals):
+ CNN = create_convolutional_neural_network_keras(input_shape, receptive_field,
+ n_filters, n_neurons_connected, n_categories,
+ eta, lmbd)
+ CNN.fit(X_train, Y_train, epochs=epochs, batch_size=batch_size, verbose=0)
+ scores = CNN.evaluate(X_test, Y_test)
+
+ CNN_keras[i][j] = CNN
+
+ print("Learning rate = ", eta)
+ print("Lambda = ", lmbd)
+ print("Test accuracy: %.3f" % scores[1])
+ print()
+
+# visual representation of grid search
+# uses seaborn heatmap, could probably do this in matplotlib
+import seaborn as sns
+
+sns.set()
+
+train_accuracy = np.zeros((len(eta_vals), len(lmbd_vals)))
+test_accuracy = np.zeros((len(eta_vals), len(lmbd_vals)))
+
+for i in range(len(eta_vals)):
+ for j in range(len(lmbd_vals)):
+ CNN = CNN_keras[i][j]
+
+ train_accuracy[i][j] = CNN.evaluate(X_train, Y_train)[1]
+ test_accuracy[i][j] = CNN.evaluate(X_test, Y_test)[1]
+
+
+fig, ax = plt.subplots(figsize = (10, 10))
+sns.heatmap(train_accuracy, annot=True, ax=ax, cmap="viridis")
+ax.set_title("Training Accuracy")
+ax.set_ylabel("$\eta$")
+ax.set_xlabel("$\lambda$")
+plt.show()
+
+fig, ax = plt.subplots(figsize = (10, 10))
+sns.heatmap(test_accuracy, annot=True, ax=ax, cmap="viridis")
+ax.set_title("Test Accuracy")
+ax.set_ylabel("$\eta$")
+ax.set_xlabel("$\lambda$")
+plt.show()
+
diff --git a/doc/src/week43/LatexFigures/fig1.aux b/doc/src/week43/LatexFigures/fig1.aux
new file mode 100644
index 000000000..b64012178
--- /dev/null
+++ b/doc/src/week43/LatexFigures/fig1.aux
@@ -0,0 +1,2 @@
+\relax
+\gdef \@abspage@last{1}
diff --git a/doc/src/week43/LatexFigures/fig1.log b/doc/src/week43/LatexFigures/fig1.log
new file mode 100644
index 000000000..3ef45abce
--- /dev/null
+++ b/doc/src/week43/LatexFigures/fig1.log
@@ -0,0 +1,458 @@
+This is pdfTeX, Version 3.141592653-2.6-1.40.22 (TeX Live 2021) (preloaded format=pdflatex 2021.4.23) 27 OCT 2021 09:29
+entering extended mode
+ restricted \write18 enabled.
+ %&-line parsing enabled.
+**fig1.tex
+(./fig1.tex
+LaTeX2e <2020-10-01> patch level 4
+L3 programming layer <2021-02-18>
+(/usr/local/texlive/2021/texmf-dist/tex/latex/standalone/standalone.cls
+Document Class: standalone 2018/03/26 v1.3a Class to compile TeX sub-files stan
+dalone
+(/usr/local/texlive/2021/texmf-dist/tex/latex/tools/shellesc.sty
+Package: shellesc 2019/11/08 v1.0c unified shell escape interface for LaTeX
+Package shellesc Info: Restricted shell escape enabled on input line 77.
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/iftex/ifluatex.sty
+Package: ifluatex 2019/10/25 v1.5 ifluatex legacy package. Use iftex instead.
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/iftex/iftex.sty
+Package: iftex 2020/03/06 v1.0d TeX engine tests
+))
+(/usr/local/texlive/2021/texmf-dist/tex/latex/xkeyval/xkeyval.sty
+Package: xkeyval 2020/11/20 v2.8 package option processing (HA)
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/xkeyval/xkeyval.tex
+(/usr/local/texlive/2021/texmf-dist/tex/generic/xkeyval/xkvutils.tex
+\XKV@toks=\toks15
+\XKV@tempa@toks=\toks16
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/xkeyval/keyval.tex))
+\XKV@depth=\count179
+File: xkeyval.tex 2014/12/03 v2.7a key=value parser (HA)
+))
+\sa@internal=\count180
+\c@sapage=\count181
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/standalone/standalone.cfg
+File: standalone.cfg 2018/03/26 v1.3a Default configuration file for 'standalon
+e' class
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/base/article.cls
+Document Class: article 2020/04/10 v1.4m Standard LaTeX document class
+(/usr/local/texlive/2021/texmf-dist/tex/latex/base/size10.clo
+File: size10.clo 2020/04/10 v1.4m Standard LaTeX file (size option)
+)
+\c@part=\count182
+\c@section=\count183
+\c@subsection=\count184
+\c@subsubsection=\count185
+\c@paragraph=\count186
+\c@subparagraph=\count187
+\c@figure=\count188
+\c@table=\count189
+\abovecaptionskip=\skip47
+\belowcaptionskip=\skip48
+\bibindent=\dimen138
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/frontendlayer/tikz.sty
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/basiclayer/pgf.sty
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/utilities/pgfrcs.sty
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfutil-common.te
+x
+\pgfutil@everybye=\toks17
+\pgfutil@tempdima=\dimen139
+\pgfutil@tempdimb=\dimen140
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfutil-common-li
+sts.tex))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfutil-latex.def
+\pgfutil@abb=\box47
+) (/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfrcs.code.tex
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/pgf.revision.tex)
+Package: pgfrcs 2020/12/27 v3.1.8b (3.1.8b)
+))
+Package: pgf 2020/12/27 v3.1.8b (3.1.8b)
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/basiclayer/pgfcore.sty
+(/usr/local/texlive/2021/texmf-dist/tex/latex/graphics/graphicx.sty
+Package: graphicx 2020/09/09 v1.2b Enhanced LaTeX Graphics (DPC,SPQR)
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/graphics/graphics.sty
+Package: graphics 2020/08/30 v1.4c Standard LaTeX Graphics (DPC,SPQR)
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/graphics/trig.sty
+Package: trig 2016/01/03 v1.10 sin cos tan (DPC)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/graphics-cfg/graphics.cfg
+File: graphics.cfg 2016/06/04 v1.11 sample graphics configuration
+)
+Package graphics Info: Driver file: pdftex.def on input line 105.
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/graphics-def/pdftex.def
+File: pdftex.def 2020/10/05 v1.2a Graphics/color driver for pdftex
+))
+\Gin@req@height=\dimen141
+\Gin@req@width=\dimen142
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/systemlayer/pgfsys.sty
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/systemlayer/pgfsys.code.tex
+Package: pgfsys 2020/12/27 v3.1.8b (3.1.8b)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfkeys.code.tex
+\pgfkeys@pathtoks=\toks18
+\pgfkeys@temptoks=\toks19
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfkeysfiltered.c
+ode.tex
+\pgfkeys@tmptoks=\toks20
+))
+\pgf@x=\dimen143
+\pgf@y=\dimen144
+\pgf@xa=\dimen145
+\pgf@ya=\dimen146
+\pgf@xb=\dimen147
+\pgf@yb=\dimen148
+\pgf@xc=\dimen149
+\pgf@yc=\dimen150
+\pgf@xd=\dimen151
+\pgf@yd=\dimen152
+\w@pgf@writea=\write3
+\r@pgf@reada=\read2
+\c@pgf@counta=\count190
+\c@pgf@countb=\count191
+\c@pgf@countc=\count192
+\c@pgf@countd=\count193
+\t@pgf@toka=\toks21
+\t@pgf@tokb=\toks22
+\t@pgf@tokc=\toks23
+\pgf@sys@id@count=\count194
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/systemlayer/pgf.cfg
+File: pgf.cfg 2020/12/27 v3.1.8b (3.1.8b)
+)
+Driver file for pgf: pgfsys-pdftex.def
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/systemlayer/pgfsys-pdftex.d
+ef
+File: pgfsys-pdftex.def 2020/12/27 v3.1.8b (3.1.8b)
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/systemlayer/pgfsys-common-p
+df.def
+File: pgfsys-common-pdf.def 2020/12/27 v3.1.8b (3.1.8b)
+)))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/systemlayer/pgfsyssoftpath.
+code.tex
+File: pgfsyssoftpath.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgfsyssoftpath@smallbuffer@items=\count195
+\pgfsyssoftpath@bigbuffer@items=\count196
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/systemlayer/pgfsysprotocol.
+code.tex
+File: pgfsysprotocol.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)) (/usr/local/texlive/2021/texmf-dist/tex/latex/xcolor/xcolor.sty
+Package: xcolor 2016/05/11 v2.12 LaTeX color extensions (UK)
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/graphics-cfg/color.cfg
+File: color.cfg 2016/01/02 v1.6 sample color configuration
+)
+Package xcolor Info: Driver file: pdftex.def on input line 225.
+Package xcolor Info: Model `cmy' substituted by `cmy0' on input line 1348.
+Package xcolor Info: Model `hsb' substituted by `rgb' on input line 1352.
+Package xcolor Info: Model `RGB' extended on input line 1364.
+Package xcolor Info: Model `HTML' substituted by `rgb' on input line 1366.
+Package xcolor Info: Model `Hsb' substituted by `hsb' on input line 1367.
+Package xcolor Info: Model `tHsb' substituted by `hsb' on input line 1368.
+Package xcolor Info: Model `HSB' substituted by `hsb' on input line 1369.
+Package xcolor Info: Model `Gray' substituted by `gray' on input line 1370.
+Package xcolor Info: Model `wave' substituted by `hsb' on input line 1371.
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcore.code.tex
+Package: pgfcore 2020/12/27 v3.1.8b (3.1.8b)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmath.code.tex
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathcalc.code.tex
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathutil.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathparser.code.tex
+\pgfmath@dimen=\dimen153
+\pgfmath@count=\count197
+\pgfmath@box=\box48
+\pgfmath@toks=\toks24
+\pgfmath@stack@operand=\toks25
+\pgfmath@stack@operation=\toks26
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.code.
+tex
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.basic
+.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.trigo
+nometric.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.rando
+m.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.compa
+rison.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.base.
+code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.round
+.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.misc.
+code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.integ
+erarithmetics.code.tex)))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfloat.code.tex
+\c@pgfmathroundto@lastzeros=\count198
+)) (/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfint.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorepoints.co
+de.tex
+File: pgfcorepoints.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@picminx=\dimen154
+\pgf@picmaxx=\dimen155
+\pgf@picminy=\dimen156
+\pgf@picmaxy=\dimen157
+\pgf@pathminx=\dimen158
+\pgf@pathmaxx=\dimen159
+\pgf@pathminy=\dimen160
+\pgf@pathmaxy=\dimen161
+\pgf@xx=\dimen162
+\pgf@xy=\dimen163
+\pgf@yx=\dimen164
+\pgf@yy=\dimen165
+\pgf@zx=\dimen166
+\pgf@zy=\dimen167
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorepathconst
+ruct.code.tex
+File: pgfcorepathconstruct.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@path@lastx=\dimen168
+\pgf@path@lasty=\dimen169
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorepathusage
+.code.tex
+File: pgfcorepathusage.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@shorten@end@additional=\dimen170
+\pgf@shorten@start@additional=\dimen171
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorescopes.co
+de.tex
+File: pgfcorescopes.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgfpic=\box49
+\pgf@hbox=\box50
+\pgf@layerbox@main=\box51
+\pgf@picture@serial@count=\count199
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoregraphicst
+ate.code.tex
+File: pgfcoregraphicstate.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgflinewidth=\dimen172
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoretransform
+ations.code.tex
+File: pgfcoretransformations.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@pt@x=\dimen173
+\pgf@pt@y=\dimen174
+\pgf@pt@temp=\dimen175
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorequick.cod
+e.tex
+File: pgfcorequick.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoreobjects.c
+ode.tex
+File: pgfcoreobjects.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorepathproce
+ssing.code.tex
+File: pgfcorepathprocessing.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorearrows.co
+de.tex
+File: pgfcorearrows.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgfarrowsep=\dimen176
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoreshade.cod
+e.tex
+File: pgfcoreshade.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@max=\dimen177
+\pgf@sys@shading@range@num=\count266
+\pgf@shadingcount=\count267
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoreimage.cod
+e.tex
+File: pgfcoreimage.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoreexternal.
+code.tex
+File: pgfcoreexternal.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgfexternal@startupbox=\box52
+))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorelayers.co
+de.tex
+File: pgfcorelayers.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoretranspare
+ncy.code.tex
+File: pgfcoretransparency.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorepatterns.
+code.tex
+File: pgfcorepatterns.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorerdf.code.
+tex
+File: pgfcorerdf.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/modules/pgfmoduleshapes.cod
+e.tex
+File: pgfmoduleshapes.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgfnodeparttextbox=\box53
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/modules/pgfmoduleplot.code.
+tex
+File: pgfmoduleplot.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/compatibility/pgfcomp-version
+-0-65.sty
+Package: pgfcomp-version-0-65 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@nodesepstart=\dimen178
+\pgf@nodesepend=\dimen179
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/compatibility/pgfcomp-version
+-1-18.sty
+Package: pgfcomp-version-1-18 2020/12/27 v3.1.8b (3.1.8b)
+))
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/utilities/pgffor.sty
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/utilities/pgfkeys.sty
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfkeys.code.tex)
+) (/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/math/pgfmath.sty
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmath.code.tex))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgffor.code.tex
+Package: pgffor 2020/12/27 v3.1.8b (3.1.8b)
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmath.code.tex)
+\pgffor@iter=\dimen180
+\pgffor@skip=\dimen181
+\pgffor@stack=\toks27
+\pgffor@toks=\toks28
+))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/frontendlayer/tikz/tikz.cod
+e.tex
+Package: tikz 2020/12/27 v3.1.8b (3.1.8b)
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/libraries/pgflibraryplothan
+dlers.code.tex
+File: pgflibraryplothandlers.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@plot@mark@count=\count268
+\pgfplotmarksize=\dimen182
+)
+\tikz@lastx=\dimen183
+\tikz@lasty=\dimen184
+\tikz@lastxsaved=\dimen185
+\tikz@lastysaved=\dimen186
+\tikz@lastmovetox=\dimen187
+\tikz@lastmovetoy=\dimen188
+\tikzleveldistance=\dimen189
+\tikzsiblingdistance=\dimen190
+\tikz@figbox=\box54
+\tikz@figbox@bg=\box55
+\tikz@tempbox=\box56
+\tikz@tempbox@bg=\box57
+\tikztreelevel=\count269
+\tikznumberofchildren=\count270
+\tikznumberofcurrentchild=\count271
+\tikz@fig@count=\count272
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/modules/pgfmodulematrix.cod
+e.tex
+File: pgfmodulematrix.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgfmatrixcurrentrow=\count273
+\pgfmatrixcurrentcolumn=\count274
+\pgf@matrix@numberofcolumns=\count275
+)
+\tikz@expandcount=\count276
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/frontendlayer/tikz/librarie
+s/tikzlibrarytopaths.code.tex
+File: tikzlibrarytopaths.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)))
+\sa@box=\box58
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/frontendlayer/tikz/librarie
+s/tikzlibrarypositioning.code.tex
+File: tikzlibrarypositioning.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/frontendlayer/tikz/librarie
+s/tikzlibrarychains.code.tex
+File: tikzlibrarychains.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/l3backend/l3backend-pdftex.def
+File: l3backend-pdftex.def 2021-03-18 L3 backend support: PDF output (pdfTeX)
+\l__color_backend_stack_int=\count277
+\l__pdf_internal_box=\box59
+)
+No file fig1.aux.
+\openout1 = `fig1.aux'.
+
+LaTeX Font Info: Checking defaults for OML/cmm/m/it on input line 3.
+LaTeX Font Info: ... okay on input line 3.
+LaTeX Font Info: Checking defaults for OMS/cmsy/m/n on input line 3.
+LaTeX Font Info: ... okay on input line 3.
+LaTeX Font Info: Checking defaults for OT1/cmr/m/n on input line 3.
+LaTeX Font Info: ... okay on input line 3.
+LaTeX Font Info: Checking defaults for T1/cmr/m/n on input line 3.
+LaTeX Font Info: ... okay on input line 3.
+LaTeX Font Info: Checking defaults for TS1/cmr/m/n on input line 3.
+LaTeX Font Info: ... okay on input line 3.
+LaTeX Font Info: Checking defaults for OMX/cmex/m/n on input line 3.
+LaTeX Font Info: ... okay on input line 3.
+LaTeX Font Info: Checking defaults for U/cmr/m/n on input line 3.
+LaTeX Font Info: ... okay on input line 3.
+(/usr/local/texlive/2021/texmf-dist/tex/context/base/mkii/supp-pdf.mkii
+[Loading MPS to PDF converter (version 2006.09.02).]
+\scratchcounter=\count278
+\scratchdimen=\dimen191
+\scratchbox=\box60
+\nofMPsegments=\count279
+\nofMParguments=\count280
+\everyMPshowfont=\toks29
+\MPscratchCnt=\count281
+\MPscratchDim=\dimen192
+\MPnumerator=\count282
+\makeMPintoPDFobject=\count283
+\everyMPtoPDFconversion=\toks30
+) (/usr/local/texlive/2021/texmf-dist/tex/latex/epstopdf-pkg/epstopdf-base.sty
+Package: epstopdf-base 2020-01-24 v2.11 Base part for package epstopdf
+Package epstopdf-base Info: Redefining graphics rule for `.eps' on input line 4
+85.
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/latexconfig/epstopdf-sys.cfg
+File: epstopdf-sys.cfg 2010/07/13 v1.3 Configuration of (r)epstopdf for TeX Liv
+e
+))
+LaTeX Font Info: External font `cmex10' loaded for size
+(Font) <7> on input line 8.
+LaTeX Font Info: External font `cmex10' loaded for size
+(Font) <5> on input line 8.
+
+[1
+
+{/usr/local/texlive/2021/texmf-var/fonts/map/pdftex/updmap/pdftex.map}]
+(./fig1.aux) )
+Here is how much of TeX's memory you used:
+ 12295 strings out of 478994
+ 254218 string characters out of 5858184
+ 526219 words of memory out of 5000000
+ 29583 multiletter control sequences out of 15000+600000
+ 403738 words of font info for 28 fonts, out of 8000000 for 9000
+ 1141 hyphenation exceptions out of 8191
+ 128i,7n,130p,424b,832s stack positions out of 5000i,500n,10000p,200000b,80000s
+
+Output written on fig1.pdf (1 page, 41515 bytes).
+PDF statistics:
+ 31 PDF objects out of 1000 (max. 8388607)
+ 22 compressed objects within 1 object stream
+ 0 named destinations out of 1000 (max. 500000)
+ 13 words of extra memory for PDF output out of 10000 (max. 10000000)
+
diff --git a/doc/src/week43/LatexFigures/fig2.aux b/doc/src/week43/LatexFigures/fig2.aux
new file mode 100644
index 000000000..b64012178
--- /dev/null
+++ b/doc/src/week43/LatexFigures/fig2.aux
@@ -0,0 +1,2 @@
+\relax
+\gdef \@abspage@last{1}
diff --git a/doc/src/week43/LatexFigures/fig2.log b/doc/src/week43/LatexFigures/fig2.log
new file mode 100644
index 000000000..5a7f8bbe2
--- /dev/null
+++ b/doc/src/week43/LatexFigures/fig2.log
@@ -0,0 +1,454 @@
+This is pdfTeX, Version 3.141592653-2.6-1.40.22 (TeX Live 2021) (preloaded format=pdflatex 2021.4.23) 27 OCT 2021 10:35
+entering extended mode
+ restricted \write18 enabled.
+ %&-line parsing enabled.
+**fig2
+(./fig2.tex
+LaTeX2e <2020-10-01> patch level 4
+L3 programming layer <2021-02-18>
+(/usr/local/texlive/2021/texmf-dist/tex/latex/standalone/standalone.cls
+Document Class: standalone 2018/03/26 v1.3a Class to compile TeX sub-files stan
+dalone
+(/usr/local/texlive/2021/texmf-dist/tex/latex/tools/shellesc.sty
+Package: shellesc 2019/11/08 v1.0c unified shell escape interface for LaTeX
+Package shellesc Info: Restricted shell escape enabled on input line 77.
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/iftex/ifluatex.sty
+Package: ifluatex 2019/10/25 v1.5 ifluatex legacy package. Use iftex instead.
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/iftex/iftex.sty
+Package: iftex 2020/03/06 v1.0d TeX engine tests
+))
+(/usr/local/texlive/2021/texmf-dist/tex/latex/xkeyval/xkeyval.sty
+Package: xkeyval 2020/11/20 v2.8 package option processing (HA)
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/xkeyval/xkeyval.tex
+(/usr/local/texlive/2021/texmf-dist/tex/generic/xkeyval/xkvutils.tex
+\XKV@toks=\toks15
+\XKV@tempa@toks=\toks16
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/xkeyval/keyval.tex))
+\XKV@depth=\count179
+File: xkeyval.tex 2014/12/03 v2.7a key=value parser (HA)
+))
+\sa@internal=\count180
+\c@sapage=\count181
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/standalone/standalone.cfg
+File: standalone.cfg 2018/03/26 v1.3a Default configuration file for 'standalon
+e' class
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/base/article.cls
+Document Class: article 2020/04/10 v1.4m Standard LaTeX document class
+(/usr/local/texlive/2021/texmf-dist/tex/latex/base/size10.clo
+File: size10.clo 2020/04/10 v1.4m Standard LaTeX file (size option)
+)
+\c@part=\count182
+\c@section=\count183
+\c@subsection=\count184
+\c@subsubsection=\count185
+\c@paragraph=\count186
+\c@subparagraph=\count187
+\c@figure=\count188
+\c@table=\count189
+\abovecaptionskip=\skip47
+\belowcaptionskip=\skip48
+\bibindent=\dimen138
+)
+\sa@box=\box47
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/frontendlayer/tikz.sty
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/basiclayer/pgf.sty
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/utilities/pgfrcs.sty
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfutil-common.te
+x
+\pgfutil@everybye=\toks17
+\pgfutil@tempdima=\dimen139
+\pgfutil@tempdimb=\dimen140
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfutil-common-li
+sts.tex))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfutil-latex.def
+\pgfutil@abb=\box48
+) (/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfrcs.code.tex
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/pgf.revision.tex)
+Package: pgfrcs 2020/12/27 v3.1.8b (3.1.8b)
+))
+Package: pgf 2020/12/27 v3.1.8b (3.1.8b)
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/basiclayer/pgfcore.sty
+(/usr/local/texlive/2021/texmf-dist/tex/latex/graphics/graphicx.sty
+Package: graphicx 2020/09/09 v1.2b Enhanced LaTeX Graphics (DPC,SPQR)
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/graphics/graphics.sty
+Package: graphics 2020/08/30 v1.4c Standard LaTeX Graphics (DPC,SPQR)
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/graphics/trig.sty
+Package: trig 2016/01/03 v1.10 sin cos tan (DPC)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/graphics-cfg/graphics.cfg
+File: graphics.cfg 2016/06/04 v1.11 sample graphics configuration
+)
+Package graphics Info: Driver file: pdftex.def on input line 105.
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/graphics-def/pdftex.def
+File: pdftex.def 2020/10/05 v1.2a Graphics/color driver for pdftex
+))
+\Gin@req@height=\dimen141
+\Gin@req@width=\dimen142
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/systemlayer/pgfsys.sty
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/systemlayer/pgfsys.code.tex
+Package: pgfsys 2020/12/27 v3.1.8b (3.1.8b)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfkeys.code.tex
+\pgfkeys@pathtoks=\toks18
+\pgfkeys@temptoks=\toks19
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfkeysfiltered.c
+ode.tex
+\pgfkeys@tmptoks=\toks20
+))
+\pgf@x=\dimen143
+\pgf@y=\dimen144
+\pgf@xa=\dimen145
+\pgf@ya=\dimen146
+\pgf@xb=\dimen147
+\pgf@yb=\dimen148
+\pgf@xc=\dimen149
+\pgf@yc=\dimen150
+\pgf@xd=\dimen151
+\pgf@yd=\dimen152
+\w@pgf@writea=\write3
+\r@pgf@reada=\read2
+\c@pgf@counta=\count190
+\c@pgf@countb=\count191
+\c@pgf@countc=\count192
+\c@pgf@countd=\count193
+\t@pgf@toka=\toks21
+\t@pgf@tokb=\toks22
+\t@pgf@tokc=\toks23
+\pgf@sys@id@count=\count194
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/systemlayer/pgf.cfg
+File: pgf.cfg 2020/12/27 v3.1.8b (3.1.8b)
+)
+Driver file for pgf: pgfsys-pdftex.def
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/systemlayer/pgfsys-pdftex.d
+ef
+File: pgfsys-pdftex.def 2020/12/27 v3.1.8b (3.1.8b)
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/systemlayer/pgfsys-common-p
+df.def
+File: pgfsys-common-pdf.def 2020/12/27 v3.1.8b (3.1.8b)
+)))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/systemlayer/pgfsyssoftpath.
+code.tex
+File: pgfsyssoftpath.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgfsyssoftpath@smallbuffer@items=\count195
+\pgfsyssoftpath@bigbuffer@items=\count196
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/systemlayer/pgfsysprotocol.
+code.tex
+File: pgfsysprotocol.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)) (/usr/local/texlive/2021/texmf-dist/tex/latex/xcolor/xcolor.sty
+Package: xcolor 2016/05/11 v2.12 LaTeX color extensions (UK)
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/graphics-cfg/color.cfg
+File: color.cfg 2016/01/02 v1.6 sample color configuration
+)
+Package xcolor Info: Driver file: pdftex.def on input line 225.
+Package xcolor Info: Model `cmy' substituted by `cmy0' on input line 1348.
+Package xcolor Info: Model `hsb' substituted by `rgb' on input line 1352.
+Package xcolor Info: Model `RGB' extended on input line 1364.
+Package xcolor Info: Model `HTML' substituted by `rgb' on input line 1366.
+Package xcolor Info: Model `Hsb' substituted by `hsb' on input line 1367.
+Package xcolor Info: Model `tHsb' substituted by `hsb' on input line 1368.
+Package xcolor Info: Model `HSB' substituted by `hsb' on input line 1369.
+Package xcolor Info: Model `Gray' substituted by `gray' on input line 1370.
+Package xcolor Info: Model `wave' substituted by `hsb' on input line 1371.
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcore.code.tex
+Package: pgfcore 2020/12/27 v3.1.8b (3.1.8b)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmath.code.tex
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathcalc.code.tex
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathutil.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathparser.code.tex
+\pgfmath@dimen=\dimen153
+\pgfmath@count=\count197
+\pgfmath@box=\box49
+\pgfmath@toks=\toks24
+\pgfmath@stack@operand=\toks25
+\pgfmath@stack@operation=\toks26
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.code.
+tex
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.basic
+.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.trigo
+nometric.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.rando
+m.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.compa
+rison.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.base.
+code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.round
+.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.misc.
+code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfunctions.integ
+erarithmetics.code.tex)))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmathfloat.code.tex
+\c@pgfmathroundto@lastzeros=\count198
+)) (/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfint.code.tex)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorepoints.co
+de.tex
+File: pgfcorepoints.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@picminx=\dimen154
+\pgf@picmaxx=\dimen155
+\pgf@picminy=\dimen156
+\pgf@picmaxy=\dimen157
+\pgf@pathminx=\dimen158
+\pgf@pathmaxx=\dimen159
+\pgf@pathminy=\dimen160
+\pgf@pathmaxy=\dimen161
+\pgf@xx=\dimen162
+\pgf@xy=\dimen163
+\pgf@yx=\dimen164
+\pgf@yy=\dimen165
+\pgf@zx=\dimen166
+\pgf@zy=\dimen167
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorepathconst
+ruct.code.tex
+File: pgfcorepathconstruct.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@path@lastx=\dimen168
+\pgf@path@lasty=\dimen169
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorepathusage
+.code.tex
+File: pgfcorepathusage.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@shorten@end@additional=\dimen170
+\pgf@shorten@start@additional=\dimen171
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorescopes.co
+de.tex
+File: pgfcorescopes.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgfpic=\box50
+\pgf@hbox=\box51
+\pgf@layerbox@main=\box52
+\pgf@picture@serial@count=\count199
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoregraphicst
+ate.code.tex
+File: pgfcoregraphicstate.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgflinewidth=\dimen172
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoretransform
+ations.code.tex
+File: pgfcoretransformations.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@pt@x=\dimen173
+\pgf@pt@y=\dimen174
+\pgf@pt@temp=\dimen175
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorequick.cod
+e.tex
+File: pgfcorequick.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoreobjects.c
+ode.tex
+File: pgfcoreobjects.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorepathproce
+ssing.code.tex
+File: pgfcorepathprocessing.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorearrows.co
+de.tex
+File: pgfcorearrows.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgfarrowsep=\dimen176
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoreshade.cod
+e.tex
+File: pgfcoreshade.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@max=\dimen177
+\pgf@sys@shading@range@num=\count266
+\pgf@shadingcount=\count267
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoreimage.cod
+e.tex
+File: pgfcoreimage.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoreexternal.
+code.tex
+File: pgfcoreexternal.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgfexternal@startupbox=\box53
+))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorelayers.co
+de.tex
+File: pgfcorelayers.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcoretranspare
+ncy.code.tex
+File: pgfcoretransparency.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorepatterns.
+code.tex
+File: pgfcorepatterns.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/basiclayer/pgfcorerdf.code.
+tex
+File: pgfcorerdf.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/modules/pgfmoduleshapes.cod
+e.tex
+File: pgfmoduleshapes.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgfnodeparttextbox=\box54
+)
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/modules/pgfmoduleplot.code.
+tex
+File: pgfmoduleplot.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/compatibility/pgfcomp-version
+-0-65.sty
+Package: pgfcomp-version-0-65 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@nodesepstart=\dimen178
+\pgf@nodesepend=\dimen179
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/compatibility/pgfcomp-version
+-1-18.sty
+Package: pgfcomp-version-1-18 2020/12/27 v3.1.8b (3.1.8b)
+))
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/utilities/pgffor.sty
+(/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/utilities/pgfkeys.sty
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgfkeys.code.tex)
+) (/usr/local/texlive/2021/texmf-dist/tex/latex/pgf/math/pgfmath.sty
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmath.code.tex))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/utilities/pgffor.code.tex
+Package: pgffor 2020/12/27 v3.1.8b (3.1.8b)
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/math/pgfmath.code.tex)
+\pgffor@iter=\dimen180
+\pgffor@skip=\dimen181
+\pgffor@stack=\toks27
+\pgffor@toks=\toks28
+))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/frontendlayer/tikz/tikz.cod
+e.tex
+Package: tikz 2020/12/27 v3.1.8b (3.1.8b)
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/libraries/pgflibraryplothan
+dlers.code.tex
+File: pgflibraryplothandlers.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgf@plot@mark@count=\count268
+\pgfplotmarksize=\dimen182
+)
+\tikz@lastx=\dimen183
+\tikz@lasty=\dimen184
+\tikz@lastxsaved=\dimen185
+\tikz@lastysaved=\dimen186
+\tikz@lastmovetox=\dimen187
+\tikz@lastmovetoy=\dimen188
+\tikzleveldistance=\dimen189
+\tikzsiblingdistance=\dimen190
+\tikz@figbox=\box55
+\tikz@figbox@bg=\box56
+\tikz@tempbox=\box57
+\tikz@tempbox@bg=\box58
+\tikztreelevel=\count269
+\tikznumberofchildren=\count270
+\tikznumberofcurrentchild=\count271
+\tikz@fig@count=\count272
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/modules/pgfmodulematrix.cod
+e.tex
+File: pgfmodulematrix.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+\pgfmatrixcurrentrow=\count273
+\pgfmatrixcurrentcolumn=\count274
+\pgf@matrix@numberofcolumns=\count275
+)
+\tikz@expandcount=\count276
+
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/frontendlayer/tikz/librarie
+s/tikzlibrarytopaths.code.tex
+File: tikzlibrarytopaths.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)))
+(/usr/local/texlive/2021/texmf-dist/tex/generic/pgf/frontendlayer/tikz/librarie
+s/tikzlibrarypositioning.code.tex
+File: tikzlibrarypositioning.code.tex 2020/12/27 v3.1.8b (3.1.8b)
+)
+(/usr/local/texlive/2021/texmf-dist/tex/latex/l3backend/l3backend-pdftex.def
+File: l3backend-pdftex.def 2021-03-18 L3 backend support: PDF output (pdfTeX)
+\l__color_backend_stack_int=\count277
+\l__pdf_internal_box=\box59
+)
+(./fig2.aux)
+\openout1 = `fig2.aux'.
+
+LaTeX Font Info: Checking defaults for OML/cmm/m/it on input line 4.
+LaTeX Font Info: ... okay on input line 4.
+LaTeX Font Info: Checking defaults for OMS/cmsy/m/n on input line 4.
+LaTeX Font Info: ... okay on input line 4.
+LaTeX Font Info: Checking defaults for OT1/cmr/m/n on input line 4.
+LaTeX Font Info: ... okay on input line 4.
+LaTeX Font Info: Checking defaults for T1/cmr/m/n on input line 4.
+LaTeX Font Info: ... okay on input line 4.
+LaTeX Font Info: Checking defaults for TS1/cmr/m/n on input line 4.
+LaTeX Font Info: ... okay on input line 4.
+LaTeX Font Info: Checking defaults for OMX/cmex/m/n on input line 4.
+LaTeX Font Info: ... okay on input line 4.
+LaTeX Font Info: Checking defaults for U/cmr/m/n on input line 4.
+LaTeX Font Info: ... okay on input line 4.
+
+(/usr/local/texlive/2021/texmf-dist/tex/context/base/mkii/supp-pdf.mkii
+[Loading MPS to PDF converter (version 2006.09.02).]
+\scratchcounter=\count278
+\scratchdimen=\dimen191
+\scratchbox=\box60
+\nofMPsegments=\count279
+\nofMParguments=\count280
+\everyMPshowfont=\toks29
+\MPscratchCnt=\count281
+\MPscratchDim=\dimen192
+\MPnumerator=\count282
+\makeMPintoPDFobject=\count283
+\everyMPtoPDFconversion=\toks30
+) (/usr/local/texlive/2021/texmf-dist/tex/latex/epstopdf-pkg/epstopdf-base.sty
+Package: epstopdf-base 2020-01-24 v2.11 Base part for package epstopdf
+Package epstopdf-base Info: Redefining graphics rule for `.eps' on input line 4
+85.
+
+(/usr/local/texlive/2021/texmf-dist/tex/latex/latexconfig/epstopdf-sys.cfg
+File: epstopdf-sys.cfg 2010/07/13 v1.3 Configuration of (r)epstopdf for TeX Liv
+e
+))
+LaTeX Font Info: External font `cmex10' loaded for size
+(Font) <7> on input line 23.
+LaTeX Font Info: External font `cmex10' loaded for size
+(Font) <5> on input line 23.
+
+[1
+
+{/usr/local/texlive/2021/texmf-var/fonts/map/pdftex/updmap/pdftex.map}]
+(./fig2.aux) )
+Here is how much of TeX's memory you used:
+ 12246 strings out of 478994
+ 253742 string characters out of 5858184
+ 524465 words of memory out of 5000000
+ 29536 multiletter control sequences out of 15000+600000
+ 403430 words of font info for 27 fonts, out of 8000000 for 9000
+ 1141 hyphenation exceptions out of 8191
+ 113i,7n,116p,409b,758s stack positions out of 5000i,500n,10000p,200000b,80000s
+
+Output written on fig2.pdf (1 page, 39214 bytes).
+PDF statistics:
+ 27 PDF objects out of 1000 (max. 8388607)
+ 19 compressed objects within 1 object stream
+ 0 named destinations out of 1000 (max. 500000)
+ 13 words of extra memory for PDF output out of 10000 (max. 10000000)
+
diff --git a/doc/src/week43/LatexFigures/fig2.tex~ b/doc/src/week43/LatexFigures/fig2.tex~
new file mode 100644
index 000000000..22c777fcd
--- /dev/null
+++ b/doc/src/week43/LatexFigures/fig2.tex~
@@ -0,0 +1,55 @@
+\documentclass[border=0.125cm]{standalone}
+\usepackage{tikz}
+\usetikzlibrary{positioning}
+\begin{document}
+
+\tikzset{%
+ every neuron/.style={
+ circle,
+ draw,
+ minimum size=1cm
+ },
+ neuron missing/.style={
+ draw=none,
+ scale=4,
+ text height=0.333cm,
+ execute at begin node=\color{black}$\vdots$
+ },
+}
+
+\begin{tikzpicture}[x=1.5cm, y=1.5cm, >=stealth]
+
+\foreach \m/\l [count=\y] in {1,2,3,missing,4}
+ \node [every neuron/.try, neuron \m/.try] (input-\m) at (0,2.5-\y) {};
+
+\foreach \m [count=\y] in {1,missing,2}
+ \node [every neuron/.try, neuron \m/.try ] (hidden-\m) at (2,2-\y*1.25) {};
+
+\foreach \m [count=\y] in {1,missing,2}
+ \node [every neuron/.try, neuron \m/.try ] (output-\m) at (4,1.5-\y) {};
+
+\foreach \l [count=\i] in {1,2,3,n}
+ \draw [<-] (input-\i) -- ++(-1,0)
+ node [above, midway] {$I_\l$};
+
+\foreach \l [count=\i] in {1,n}
+ \node [above] at (hidden-\i.north) {$H_\l$};
+
+\foreach \l [count=\i] in {1,n}
+ \draw [->] (output-\i) -- ++(1,0)
+ node [above, midway] {$O_\l$};
+
+\foreach \i in {1,...,4}
+ \foreach \j in {1,...,2}
+ \draw [->] (input-\i) -- (hidden-\j);
+
+\foreach \i in {1,...,2}
+ \foreach \j in {1,...,2}
+ \draw [->] (hidden-\i) -- (output-\j);
+
+\foreach \l [count=\x from 0] in {Input, Hidden, Ouput}
+ \node [align=center, above] at (\x*2,2) {\l \\ layer};
+
+\end{tikzpicture}
+
+\end{document}
diff --git a/doc/src/week43/LatexFigures/nn.py~ b/doc/src/week43/LatexFigures/nn.py~
new file mode 100644
index 000000000..d41673ce2
--- /dev/null
+++ b/doc/src/week43/LatexFigures/nn.py~
@@ -0,0 +1,51 @@
+import matplotlib.pyplot as plt
+
+def draw_neural_net(ax, left, right, bottom, top, layer_sizes):
+ '''
+ Draw a neural network cartoon using matplotilb.
+
+ :usage:
+ >>> fig = plt.figure(figsize=(12, 12))
+ >>> draw_neural_net(fig.gca(), .1, .9, .1, .9, [4, 7, 2])
+
+ :parameters:
+ - ax : matplotlib.axes.AxesSubplot
+ The axes on which to plot the cartoon (get e.g. by plt.gca())
+ - left : float
+ The center of the leftmost node(s) will be placed here
+ - right : float
+ The center of the rightmost node(s) will be placed here
+ - bottom : float
+ The center of the bottommost node(s) will be placed here
+ - top : float
+ The center of the topmost node(s) will be placed here
+ - layer_sizes : list of int
+ List of layer sizes, including input and output dimensionality
+ '''
+ n_layers = len(layer_sizes)
+ v_spacing = (top - bottom)/float(max(layer_sizes))
+ h_spacing = (right - left)/float(len(layer_sizes) - 1)
+ # Nodes
+ for n, layer_size in enumerate(layer_sizes):
+ layer_top = v_spacing*(layer_size - 1)/2. + (top + bottom)/2.
+ for m in range(layer_size):
+ circle = plt.Circle((n*h_spacing + left, layer_top - m*v_spacing), v_spacing/4.,
+ color='w', ec='k', zorder=4)
+ ax.add_artist(circle)
+ # Edges
+ for n, (layer_size_a, layer_size_b) in enumerate(zip(layer_sizes[:-1], layer_sizes[1:])):
+ layer_top_a = v_spacing*(layer_size_a - 1)/2. + (top + bottom)/2.
+ layer_top_b = v_spacing*(layer_size_b - 1)/2. + (top + bottom)/2.
+ for m in range(layer_size_a):
+ for o in range(layer_size_b):
+ line = plt.Line2D([n*h_spacing + left, (n + 1)*h_spacing + left],
+ [layer_top_a - m*v_spacing, layer_top_b - o*v_spacing], c='k')
+ ax.add_artist(line)
+
+
+fig = plt.figure(figsize=(12, 12))
+ax = fig.gca()
+ax.axis('off')
+draw_neural_net(ax, .1, .9, .1, .9, [4, 7, 2])
+fig.savefig('nn.png')
+fig.show()
diff --git a/doc/src/week45/Datafiles/ride.csv~ b/doc/src/week45/Datafiles/ride.csv~
new file mode 100644
index 000000000..7a908c5f7
--- /dev/null
+++ b/doc/src/week45/Datafiles/ride.csv~
@@ -0,0 +1,15 @@
+Outlook,Temperature,Humidity,Wind,Ride
+Sunny,Hot,High,Weak,0
+Sunny,Hot,High,Strong,1
+Overcast,Hot,High,Weak,1
+Rain,Mild,High,Weak,1
+Rain,Cool,Normal,Weak,1
+Rain,Cool,Normal,Strong,0
+Overcast,Cool,Normal,Strong,1
+Sunny,Mild,High,Weak,0
+Sunny,Cool,Normal,Weak,1
+Rain,Mild,Normal,Weak,1
+Sunny,Mild,Normal,Strong,1
+Overcast,Mild,High,Strong,1
+Overcast,Hot,Normal,Weak,1
+Rain,Mild,High,Strong,0