Deep Learning#
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import numpy as np
import pandas as pd
import os.path
import subprocess
import matplotlib.collections
import scipy.signal
from sklearn import model_selection
import warnings
warnings.filterwarnings('ignore')
import tensorflow.compat.v1 as tf
tf.disable_v2_behavior()
print(tf.__version__)
Helpers for Getting, Loading and Locating Data
def wget_data(url: str):
local_path = './tmp_data'
p = subprocess.Popen(["wget", "-nc", "-P", local_path, url], stderr=subprocess.PIPE, encoding='UTF-8')
rc = None
while rc is None:
line = p.stderr.readline().strip('\n')
if len(line) > 0:
print(line)
rc = p.poll()
def locate_data(name, check_exists=True):
local_path='./tmp_data'
path = os.path.join(local_path, name)
if check_exists and not os.path.exists(path):
raise RuxntimeError('No such data file: {}'.format(path))
return path
Get Data#
wget_data('https://raw.githubusercontent.com/illinois-ipaml/MachineLearningForPhysics/main/data/circles_data.hf5')
wget_data('https://raw.githubusercontent.com/illinois-ipaml/MachineLearningForPhysics/main/data/circles_targets.hf5')
wget_data('https://raw.githubusercontent.com/illinois-ipaml/MachineLearningForPhysics/main/data/spectra_data.hf5')
Neural Network Architectures for Deep Learning#
We previously took a bottom-up look at how a neural network is composed of basic building blocks. Now, we take a top-down look at some of the novel network architectures that are enabling the current deep-learning revolution:
Convolutional networks
Unsupervised learning networks
Recurrent networks
Reinforcement learning
We conclude with some reflections on where “deep learning” is headed.
The examples below use higher-level tensorflow APIs than we have seen before, so we start with a brief introduction to them.
High-Level Tensorflow APIs#
In our earlier examples, we built our networks using low-level tensorflow primitives. For more complex networks composed of standard building blocks, there are convenient higher-level application programming interfaces (APIs) that abstract aways the low-level graphs and sessions.
Reading Data#
The tf.data API handles data used to train and test a network, replacing the low-level placeholders we used earlier. For a small dataset that fits in memory, use:
dataset = tf.data.Dataset.from_tensor_slices((dict(X), y))
Creating a Dataset adds nodes to a graph so you should normally wrap your code to create a Dataset in a function that tensorflow will call in the appropriate context. For example, to split the 300 circles samples above into train (200) and test (100) datasets:
X = pd.read_hdf(locate_data('circles_data.hf5'))
y = pd.read_hdf(locate_data('circles_targets.hf5'))
X_train, X_test, y_train, y_test = model_selection.train_test_split(
X, y, test_size=100, random_state=123)
def get_train_data(batch_size=50):
dataset = tf.data.Dataset.from_tensor_slices((dict(X_train), y_train))
return dataset.shuffle(len(X_train)).repeat().batch(batch_size)
def get_test_data(batch_size=50):
dataset = tf.data.Dataset.from_tensor_slices((dict(X_test), y_test))
return dataset.batch(batch_size)
While from_tensor_slices is convenient, it is not very efficient since the whole dataset is added to the graph with constant nodes (and potentially copied multiple times). Alternatively, convert your data to tensorflow’s binary file format so it can be read as a TFRecordDataset.
Building a Model#
The tf.estimator API builds and runs a graph for training, evaluation and prediction. This API generates a lot of INFO log messages, which can be suppressed using:
tf.logging.set_verbosity(tf.logging.WARN)
First specify the names and types (but not values) of the features that feed the network’s input layer:
inputs = [tf.feature_column.numeric_column(key=key) for key in X]
Next, build the network graph. There are pre-made estimators for standard architectures that are easy to use. For example, to recreate our earlier architecture of a single 4-node hidden layer with sigmoid activation:
config = tf.estimator.RunConfig(
model_dir='tfs/circle',
tf_random_seed=123
)
classifier = tf.estimator.DNNClassifier(
config=config,
feature_columns=inputs,
hidden_units=[4],
activation_fn=tf.sigmoid,
n_classes=2
)
There are only a limited number of pre-defined models available so you often have to build a custom estimator using the intermediate-level layers API. See convolutional-network example below.
Training a Model#
An estimator remembers any previous training (using files saved to its model_dir) so if you really want to start from scratch you will need to clear this history:
!rm -rf tfs/circle/*
The train method runs a specified number of steps (each learning from one batch of training data):
classifier.train(input_fn=get_train_data, steps=5000);
After training, you can list the model parameters and access their values:
classifier.get_variable_names()
classifier.get_variable_value('dnn/hiddenlayer_0/kernel')
Testing a Model#
results = classifier.evaluate(input_fn=get_test_data)
results
Convolutional Networks#
A Convolutional Neural Network (CNN) is a special architecture that:
Assumes that input features measure some property on a grid. The grid is usually spatial or temporal, but this is not required. For example, a 1D spectrum or time series, a 2D monochrome image, or a 3D stack of 2D images in different filters (RGB, etc).
Performs translation-invariant learning efficiently. For example, identifying a galaxy wherever it appears in an image, or a transient pulse wherever it appears in a time series. The main efficiency is a much reduced number of parameters compared to the number of input features, relative to the dense fully connected networks we have seen so far.
As we saw in the previous lecture, Neural Networks receive an input (a single vector), and transform it through a series of hidden layers. Each hidden layer is made up of a set of neurons, where each neuron is fully connected to all neurons in the previous layer, and where neurons in a single layer function completely independently and do not share any connections. The last fully-connected layer is called the “output layer” and in classification settings it represents the class scores.
The fully-connected, feed-forward neural networks we have studied thus far do not scale well to large image data. For example, a modest 200 \(\times\) 200 \(\times\) 3 (x-pixels, y-pixels, 3 colors) image would lead to neurons that have 200 \(\times\) 200 \(\times\) 3 = 120,000 weights. Moreover, we would almost certainly want to have several such neurons, so the parameters would add up quickly! Clearly, this full connectivity is wasteful and the huge number of parameters would quickly lead to overfitting.
Convolutional Neural Networks take advantage of the fact that the input consists of images and they constrain the architecture in a more sensible way to reduce the number of parameters. In particular, unlike a regular Neural Network, the layers of a CNN have neurons arranged in 3 dimensions: width, height, depth.
(Note that the word “depth” here refers to the third dimension of an activation volume, not to the depth of a full Neural Network, which can refer to the total number of layers in a network…)
The neurons in a CNN layer will only be connected to a small region of the layer before it, instead of all of the neurons in a fully-connected manner.
A CNN is made up of layers of different types (convolutions, pooling, fully-connected), in general. Every layer has a simple API: It transforms an input 3D volume to an output 3D volume with some differentiable function that may or may not have parameters.
We will use the following problem to motivate and demonstration a CNN:
The input data consists of triplets of digitized waveforms.
Each waveform has a slowly varying level with some narrow pulses superimposed.
Each triplet has a single pulse that is synchronized (coincident) in all three waveforms.
Waveforms also contain a random number of unsynchronized “background” pulses.
Synchronized and unsynchronized pulses can overlap in time and between traces.
The goal is to identify the location of the synchronized pulses in each triplet. This is a simplified version of a common task in data acquisition trigger systems and transient analysis pipelines.
def generate(N=10000, ntrace=3, nt=100, nbg=1., A=5., nsmooth=3, T=1., seed=123):
gen = np.random.RandomState(seed=seed)
t_grid = np.linspace(0., T, nt)
# Generate the smooth background shapes as superpositions of random cosines.
wlen = 2 * T * gen.lognormal(mean=0., sigma=0.2, size=(nsmooth, N, ntrace, 1))
phase = gen.uniform(size=wlen.shape)
X = np.cos(2 * np.pi * (t_grid + phase * wlen) / wlen).sum(axis=0)
# Superimpose short pulses.
sigma = 0.02 * T
tsig = T * gen.uniform(0.05, 0.95, size=N)
y = np.empty(N, dtype=int)
nbg = gen.poisson(lam=nbg, size=(N, ntrace))
for i in range(N):
# Add a coincident pulse to all traces.
xsig = A * np.exp(-0.5 * (t_grid - tsig[i]) ** 2 / sigma ** 2)
y[i] = np.argmax(xsig)
X[i] += xsig
# Add non-coincident background pulses to each trace.
for j in range(ntrace):
if nbg[i, j] > 0:
t0 = T * gen.uniform(size=(nbg[i, j], 1))
X[i, j] += (A * np.exp(-0.5 * (t_grid - t0) ** 2 / sigma ** 2)).sum(axis=0)
return X.astype(np.float32), y
X, y = generate()
def plot_traces(X, y):
Nsample, Ntrace, D = X.shape
_, ax = plt.subplots(Nsample, 1, figsize=(9, 1.5 * Nsample))
t = np.linspace(0., 1., 100)
dt = t[1] - t[0]
for i in range(Nsample):
for j in range(Ntrace):
ax[i].plot(t, X[i, j], lw=1)
ax[i].axvline(t[y[i]], c='k', ls=':')
ax[i].set_yticks([])
ax[i].set_xticks([])
ax[i].set_xlim(-0.5 * dt, 1 + 0.5 * dt)
plt.subplots_adjust(left=0.01, right=0.99, bottom=0.01, top=0.99, hspace=0.1)
plot_traces(X[:5], y[:5])
The derivative of \(f(x)\) can be approximated as
for small \(\delta\). We can use this approximation to convert an array of \(f(n \Delta x)\) values into an array of estimated \(f'(n \Delta x)\) values using:
K = np.array([-1, 0, +1]) / ( 2 * dx)
fp[0] = K.dot(f[[0,1,2]])
fp[1] = K.dot(f[[1,2,3]])
...
fp[N-2] = K.dot(f[[N-3,N-2,N-1]]
The numpy convolve function automates this process of sliding an arbitrary kernel \(K\) along an input array like this. The result only estimates a first (or higher-order) derivative when the kernel contains special values (and you should normally use the numpy gradient function for this), but any convolution is a valid and potentially useful transformation.
A clarifying word about terminology: In the context of convolutional networks, kernel is a simple group of weights shared all over the input space that is engineered to determine what specific features are to be detected. The kernel is also sometimes referred to as a “feature map” or “filter” in this context.
See for example the application of a kernel in convolution over a simple black-and-white image: here.
The kernel needs to completely overlap the input array it is being convolved with, which means that the output array is smaller and offset. Alternatively, you can pad the input array with zeros to extend the output array. There are three different conventions for handling these edge effects via the mode parameter to np.convolve:
valid: no zero padding, so output length is \(N - K + 1\) and offset is \((K-1)/2\).
same: apply zero padding and trim so output length equals input length \(N\), and offset is zero.
full: apply zero padding without trimming, so output length is \(N + K - 1\) and offset is \(-(K-1)/2\).
(Here \(N\) and \(K\) are the input and kernel lengths, respectively).
We can use a convolution to identify features in our input data:
def plot_convolved(x, kernel, smax=50):
t = np.arange(len(x))
plt.plot(t, x, lw=1, c='gray')
z = np.convolve(x, kernel, mode='same')
for sel, c in zip(((z > 0), (z < 0)), 'rb'):
plt.scatter(t[sel], x[sel], c=c, s=smax * np.abs(z[sel]), lw=0)
plt.gca()
plt.grid(False)
First, let’s pick out regions of large positive (red) or negative slope (notice how the edge padding causes some artifacts):
plot_convolved(X[1, 1], [0.5,0,-0.5])
We can also pick out regions of large curvature (using the finite-difference coefficients for a second derivative):
plot_convolved(X[1, 1], [1.,-2.,1.])
We can apply both of these convolutions to transform our input data to a new representation that highlights regions of large first or second derivative. Use a tanh activation to accentuate the effect:
def apply_convolutions(X, *kernels):
N1, N2, D = X.shape
out = []
for i in range(N1):
sample = []
for j in range(N2):
for K in kernels:
sample.append(np.tanh(np.convolve(X[i, j], K, mode='valid')))
out.append(sample)
return np.asarray(out)
out = apply_convolutions(X, [0.5,0,-0.5], [1.,-2.,1.])
The resulting array can be viewed as a synthetic image and offers an easy way to visually identify individual narrow peaks and their correlations between traces:
def plot_synthetic(Z):
_, ax = plt.subplots(len(Z), 1, figsize=(9, len(Z)))
for i, z in enumerate(Z):
ax[i].imshow(z, aspect='auto', origin='upper', interpolation='none',
cmap='coolwarm', vmin=-1, vmax=+1);
ax[i].grid(False)
ax[i].axis('off')
plt.subplots_adjust(left=0.01, right=0.99, bottom=0.01, top=0.99, hspace=0.1)
plot_synthetic(out[:5])
The patterns that identify individual and coincident peaks are all translation invariant so can be identified in this array using a new convolution, but now in the 2D space of these synthetic images.
Since matrix convolution is a linear operation, it is a special case of our general neural network unit,
but with the matrix \(W\) now having many repeated elements so its effective number of dimensions is greatly reduced in typical applications.
A convolutional layer takes an arbitrary input array and applies a number of filters with the same shape in parallel. By default, the filter kernels march with single-element steps through the input array, but you can also specify larger stride vector.
In the general case, the input array, kernels and stride vector are all multidimensional, but with the same dimension. Tensorflow provides convenience functions for 1D, 2D and 3D convolutional layers, for example:
hidden = tf.layers.Conv2D(
filters=3, kernel_size=[4, 5], strides=[2, 1],
padding='same', activation=tf.nn.relu)
Note that padding specifies how edges effects are handled, but only same and valid are supported (and valid is the default). You can also implement higher-dimensional convolutional layers using the lower-level APIs.
A convolutional neural network (CNN) is a network containing convolutional layers. A typical architecture starts with convolutional layers, processing the input, then finishes with some fully connected dense layers to calculate the output. Since one of the goals of a CNN is reduce the number of parameters, a CNN often also incorporates pooling layers to reduce the size of the array fed to to later layers by “downsampling” (typically using a maximum or mean value). See these Stanford CS231n notes for more details in the context of image classification.
def pulse_model(features, labels, mode, params):
"""Build a graph to TRAIN/TEST/PREDICT a pulse coincidence detection model.
"""
D = params['time_steps']
M = params['number_of_traces']
n1 = params['conv1_width']
n2 = params['conv2_width']
eta = params['learning_rate']
assert n1 % 2 == 1 and n2 % 2 == 1
# Build the input layer.
inputs = tf.reshape(features['X'], [-1, M, D, 1])
# Add the first convolutional layer.
conv1 = tf.layers.conv2d(
inputs=inputs, filters=2, kernel_size=[1, n1],
padding='same', activation=tf.tanh, name='conv1')
# Add the second convolutional (and output) layer.
logits = tf.layers.conv2d(
inputs=conv1, filters=1, kernel_size=[M, n2],
padding='valid', activation=None, name='conv2')
# Flatten the outputs.
logits = tf.reshape(logits, [-1, D - n2 + 1])
# Calculate the offset between input labels and the output-layer node index
# that is introduced by using padding='valid' for the output layer below.
offset = (n2 - 1) // 2
# Calculate the network's predicted best label.
predicted_labels = tf.argmax(logits, axis=1) + offset
# Calculate the network's predicted probability of each label.
probs = tf.nn.softmax(logits)
# Calculate the network's predicted mean label.
bins = tf.range(0., D - n2 + 1., dtype=np.float32) + offset
mean_labels = tf.reduce_sum(bins * probs, axis=-1)
# Return predicted labels and probabilities in PREDICT mode.
if mode == tf.estimator.ModeKeys.PREDICT:
return tf.estimator.EstimatorSpec(mode, predictions={
'label': predicted_labels,
'probs': tf.nn.softmax(logits)
})
# Calculate the loss for TRAIN and EVAL modes. We need to offset the labels
# used here so they correspond to output-layer node indices.
loss = tf.losses.sparse_softmax_cross_entropy(labels=labels - offset, logits=logits)
# Compute evaluation metrics.
if mode == tf.estimator.ModeKeys.EVAL:
accuracy = tf.metrics.accuracy(labels=labels, predictions=predicted_labels)
rmse = tf.metrics.root_mean_squared_error(
labels=tf.cast(labels, np.float32), predictions=mean_labels)
return tf.estimator.EstimatorSpec(
mode, loss=loss, eval_metric_ops={'accuracy': accuracy, 'rmse': rmse})
# Create optimizer.
assert mode == tf.estimator.ModeKeys.TRAIN
optimizer = tf.train.AdamOptimizer(learning_rate=eta)
step = optimizer.minimize(loss, global_step=tf.train.get_global_step())
return tf.estimator.EstimatorSpec(mode, loss=loss, train_op=step)
tf.logging.set_verbosity(tf.logging.WARN)
!rm -rf tfs/pulses
config = tf.estimator.RunConfig(
model_dir='tfs/pulses',
tf_random_seed=123
)
pulse = tf.estimator.Estimator(
config=config,
model_fn=pulse_model,
params = dict(
time_steps=100,
number_of_traces=3,
conv1_width=3,
conv2_width=7,
learning_rate=0.01))
X_train, X_test, y_train, y_test = model_selection.train_test_split(
X, y, test_size=0.4, random_state=123)
pulse.train(
input_fn=tf.estimator.inputs.numpy_input_fn(
x={'X': X_train}, y=y_train,
batch_size=500, num_epochs=None, shuffle=True),
steps=500);
Compare the kernels learned during training with the derivative kernels we used above. We find that they are qualitatively similar:
The “odd” kernel correlates most strongly with a rising slope, so approximately measures \(+f'(t)\).
The “even” kernel correlates most strongly with a local maximum, so approximately measures \(-f''(t)\).
The odd-numbered rows of the image are correlated with the odd kernel, and correlate with a pulse that rises (red) on the left and falls on the right (blue).
The even-numbered rows of the image are correlated with the even kernel, and correlate with a pulse that peaks (dark red) at the center.
Note that nothing in the network architecture requires that the three traces be processed the same way in the second convolutional layer (right-hand image), and we do find some variations. A more detailed analysis of these weights would take into account the additional bias parameters and the influence of the activations.
def plot_kernels():
M = pulse.params['number_of_traces']
n1 = pulse.params['conv1_width']
n2 = pulse.params['conv2_width']
K1 = pulse.get_variable_value('conv1/kernel')
K2 = pulse.get_variable_value('conv2/kernel')
assert K1.shape == (1, n1, 1, 2)
assert K2.shape == (M, n2, 2, 1)
_, ax = plt.subplots(1, 2, figsize=(10, 3))
# Plot the two 1D kernels used in the first layer.
dt = np.arange(n1) - 0.5 * (n1 - 1)
ax[0].plot(dt, K1[0, :, 0, 0], 'o:', label='even')
ax[0].plot(dt, K1[0, :, 0, 1], 'o:', label='odd')
ax[0].legend(fontsize='x-large')
# Assemble an image of the second-layer kernel that can be compared with plot_synthetic().
K2img = np.empty((M, 2, n2))
K2img[:, 0] = K2[:, :, 0, 0]
K2img[:, 1] = K2[:, :, 1, 0]
vlim = np.max(np.abs(K2))
ax[1].imshow(K2img.reshape(2 * M, n2), aspect='auto', origin='upper',
interpolation='none', cmap='coolwarm', vmin=-vlim, vmax=+vlim)
ax[1].axis('off')
ax[1].grid(False)
plt.tight_layout()
plot_kernels()
Evaluate how well the trained network performs on the test data:
results = pulse.evaluate(
input_fn=tf.estimator.inputs.numpy_input_fn(
x={'X': X_test}, y=y_test,
num_epochs=1, shuffle=False))
We find that about 95% of test samples are classified “correctly”, defined as the network predicting the bin containing the the coincidence maximum exactly. However, The RMS error between the predicted and true bins is only 0.4 bins, indicating that the network usually predicts a neighboring bin in the 5% of “incorrect” test cases.
results
Finally, compare the predicted (gray histogram) and true (dotted line) coincidence locations for a few test samples:
def plot_predictions(X, y):
# Calculate predicted labels and PDFs over labels.
predictions = pulse.predict(
input_fn=tf.estimator.inputs.numpy_input_fn(
x={'X': X}, y=None, num_epochs=1, shuffle=False))
Nsample, Ntrace, D = X.shape
t = np.linspace(0., 1., 100)
dt = t[1] - t[0]
bins = np.linspace(-0.5 * dt, 1 + 0.5 * dt, len(t) + 1)
probs = np.zeros(D)
# Plot input data, truth, and predictions.
_, ax = plt.subplots(Nsample, 1, figsize=(9, 1.5 * Nsample))
for i, pred in enumerate(predictions):
label = pred['label']
# Plot the input traces.
for x in X[i]:
ax[i].plot(t, x, lw=1)
# Indicate the true coincidence position.
ax[i].axvline(t[y[i]], c='k', ls=':')
# Indicate the predicted probability distribution.
n2 = D - len(pred['probs']) + 1
offset = (n2 - 1) // 2
probs[offset:-offset] = pred['probs']
rhs = ax[i].twinx()
rhs.hist(t, weights=probs, bins=bins, histtype='stepfilled', alpha=0.25, color='k')
rhs.set_ylim(0., 1.)
rhs.set_xlim(bins[0], bins[-1])
rhs.set_yticks([])
ax[i].set_xticks([])
ax[i].set_yticks([])
ax[i].grid(False)
ax[i].set_xlim(bins[0], bins[-1])
plt.subplots_adjust(left=0.01, right=0.99, bottom=0.01, top=0.99, hspace=0.1)
plot_predictions(X_test[:5], y_test[:5])
Note that our loss function does not know that consecutive labels are close and being off by one is almost as good as getting the right label. We could change this by treating this as a regression problem, but a nice feature of our multi-category approach is that we can predict a a full probability density over labels (the gray histograms above) which is often useful.
Recurrent Networks#
All the architectures we have seen so far are feed-foward networks, with input data always from left (input layer) to right (output layer). A recurrent neural network (RNN) adds links that feed back into a previous layer. This simple modification adds significant complexity but also expressive power (comparable to the electronics revolution associated with the idea of transistor feedback).
Architectures with feedback are still maturing but some useful building blocks have emerged, such as the long short-term memory unit, which allows a network to remember some internal state but also forget it based on new input.
Some practical considerations for RNN designs:
The order of training data is now significant and defines a “model time”, but the network can be reset whenever needed.
Input data can be packaged into variable-length messages that generate variable (and different) length output messages. This is exactly what language translation needs.
Optimization of the weights using gradients is still possible but requires “unrolling” the network by cloning it enough times to process the longest allowed messages.
A feed-foward network implements a universal approximating function. Since the internal state of an RNN acts like local variables, you can think of an RNN as a universal approximating program.
See this blog post for an example based on natural language synthesis.
Reinforcement Learning#
The architectures we have seen so far all have target output values associated with each input sample, which are necessary to update the network parameters during the learning (loss optimization) phase:
However, we can relax this requirement of being able to calculate a loss after each new input as long as we eventually get some feedback on how well our input-to-output mapping is doing. This is the key idea of reinforcement learning (RL):
A RL network watches some external “reality” (which is often simulated) and learns a policy for how to take actions. A sequence of actions eventually leads to some feedback, which is then used to take a single step in optimizing the policy network’s parameters:
See this blog post for an example based on image generation.
See this blog post for an example based on image generation.
Deep Learning Outlook#
The depth of “deep learning” comes primarily from network architectures that stack many layers. In another sense, deep learning is very shallow since it often performs well using little to no specific knowledge about the problem it is solving, using generic building blocks.
The field of modern deep learning started around 2012 when the architectures described above were first used successfully, and the necessary large-scale computing and datasets were available. Massive neural networks are now the state of the art for many benchmark problems, including image classification, speech recognition and language translation.
However, less than a decade into the field, there are signs that deep learning is reaching its limits. Some of the pioneers are focusing on new directions such as capsule networks and causal inference. Others are taking a critical look at the current state of the field:
Deep learning does not use data efficiently.
Deep learning does not integrate prior knowledge.
Deep learning often give correct answers but without associated uncertainties.
Deep learning applications are hard to interpret and transfer to related problems.
Deep learning is excellent at learning stable input-output mappings but does not cope well with varying conditions.
Deep learning cannot distinguish between correlation and causation.
These are mostly concerns for the future of neural networks as a general model for artificial intelligence, but they also limit the potential of scientific applications.
However, there are many challenges in scientific data analysis and interpretation that could benefit from deep learning approaches, so I encourage you to follow the field and experiment. Through this course, you now have a pretty solid foundation in data science and machine learning to further your studies toward more advanced and current topics!