Solving engineering problems using scikit-learn: Curve fitting with neural network

Image Source

Machine learning is essentially computer doing the learning for you. It doesn't have to just apply to big data. In this series, I hope to capture ways in which machine learning can be used to solve more general engineering problems. In this post, I will show you how to use simple neural network for curve fitting.

This problem is based on a publication in Scientific Reports by Sun et. al.. In their paper, they performed simulations to calculate the bandgap of superlattices (repeated lattices of quantum wells) with various barrier and well thickness:

Source: Sun et. al., Scientific Reports

The result, as extracted from their plots, is shown below:

Obviously since each simulation takes a long time, it is impossible to simulate every single combination. The question then, for a researcher reading this paper, how do he or she find the bandgap of any combination of barrier and well width?

This is essentially a curve fitting problem. There should be a surface, in other words a 2D function that fits this data. The difference is that unlike most curve fitting problems where the form of the function is known, we don't actually know the form of the function here. It is also clear from the plot that the function is non-linear, which means linear interpolation between points will not be accurate.
But we can use machine learning! Machine learning is letting the machine decide what model fits the data the best. The good thing about machine learning is we don't need to know the form of these function. In fact, we can't use it to obtain a mathematical or theoretical model for explaining the phenomena. However, if our goal is simply to obtain the correct value of the unknown function f(X) given a vector of X input, machine learning can be a good solution.

To inject non-linearity, the machine learning algorithm used here will be a simple neural network. I won’t be going in details of what a neural network is at this point, here is a brief introduction:

Neural neworks are typically organized in layers. Layers are made up of a number of interconnected 'nodes' which contain an 'activation function'. Patterns are presented to the network via the 'input layer', which communicates to one or more 'hidden layers' where the actual processing is done via a system of weighted 'connections'. The hidden layers then link to an 'output layer' where the answer is output.

Source: A Basic Introduction To Neural Networks

Simple neural networks can be easily implemented using scikit-learn. The model we are using is neural_network.MLPRegressor. MLP stands for multi-layer Perceptrons, which is just another name for neural networks. Note that only the scikit-learn version v1.80 and above have this implemented, so make sure your scikit-learn is up to date using the command:

pip install -U scikit-learn 

at the command prompt. You should get the following messages:

Next, open your favourite python editor. I have used Spyder for my work. Same as all other machine learning problems, we following the machine learning pipeline:

1. Load data into the data into dataframe using pandas. Here, the data in figure 2 above is entered into a csv with columns named “Barrier”, “Well” and “Bandgap”
2. Separate input and output. Here, the input is the barrier and well-width, and output is the bandgap
3. Create model and train using the input
4. Predict output using the model, and compare the difference between the predicted value and the actual value

Note that since these are clean data, we do not need to worry about training vs testing data. We implement the neural network as follows:

import pandas as pd
from sklearn.neural_network import MLPRegressor
df = pd.read_csv('Bandgap.csv')    
Y = df.Bandgap      
X = df.drop('Bandgap', axis = 1)  
reg = MLPRegressor(hidden_layer_sizes = (4,4), activation = 'tanh', solver = 'lbfgs')  
reg.fit(X, Y)
Y_pred = reg.predict(X)
df['Y_pred'] = Y_pred
df['Difference'] = df.Bandgap - df.Y_pred

Here, I am only using a very small neural network of two hidden layers each with 4 neurons. This is specified with the parameter hidden_layer_size in the form of a tuple. You can add in more hidden layers but increasing the number of elements in the tuple. There is no need for a bigger network as the dataset is quite small and is very clean with well-defined boundaries. Using too big of a network would result in overfitting. Besides, scikit-learn neural network API is not designed for large networks. Another parameters that I have specified are the activation function, which I have specified as activation = “tanh” instead of the default ReLu, as tanh is less abrupt and therefore better for regression problems. The optimisation method was also specified as solver = 'lbfgs', which uses a set of Newtons method, which is more suitable to small dataset with a well defined distribution. A schematic of the neural network is shown below:

The weights of each neuron, which is the weight vector that is applied to the input vector to obtain the output of the neuron can be extracted by extracting the coefs_ parameter of the model reg:

coeff = reg.coefs_
print('W01 =', coeff[0][0])
print('W02 =', coeff[0][1])
print('W11 =', coeff[1][0])
print('W12 =', coeff[1][1])
print('W13 =', coeff[1][2])
print('W14 =', coeff[1][3])
print('W21 =', coeff[2][0])
print('W22 =', coeff[2][1])
print('W23 =', coeff[2][2])
print('W24 =', coeff[2][3])

which gives the following results:

So for example, the input data "barrier_width" would need to multiply elementwise by W01, resulting in four numbers that are the input to the 4 neurons in the 1st hidden layer, etc.

As a sanity check, lets look at what it predicted for the given combination:

It can be seen that the neural network is modelling the original datapoints very accurately. The biggest difference between the prediction and the simulation is 0.025eV, which is sufficient for most applications. We should now examine whether the neural network resulted in a well-defined surface over this parameter space. We can do this using the following script:

import numpy as np
barrier_grid = [i for i in np.arange(0.5, 9.5, 0.025)]
well_grid = [j for j in np.arange(0.5, 4.5, 0.025)]
barrier_list, well_list = np.meshgrid(barrier_grid, well_grid)
df2 = pd.DataFrame()
df2['Barrier'] = barrier_list.reshape((1, barrier_list.size))[0]
df2['Well'] = well_list.reshape((1, well_list.size))[0]
Y_pred2 = reg.predict(df2)
df2['Bandgap'] = Y_pred2

This warrant some explanations here. The meshgrid() function creates a mesh based on the list of barrier and well thickness that we provided (here each grid point is 0.025 in both direction from the next grid point). This saves us from using a double for loop to match every barrier value with every well value. However, the resulting mesh is a 2D array for each of the barrier and well, which cannot be directly used in scikit-learn. We need to convert it back to a vector, and hence the use of the reshape() function.

The resulting surface can then be plotted as a scatterplot, using matplotlib. The code is as follows:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(barrier_list, well_list, Y_pred2, color='r', marker = '.')
ax.scatter(df.Barrier, df.Well, df.Bandgap, c = 'b', marker = 'o')
ax.set_xlabel('Barrier thickness')
ax.set_ylabel('Well thickness')
ax.set_zlabel('Bandgap')
fig.show()

Note the use of axes ax as a handle for the plot, allowing us to do things such as putting two different datasets onto the same axis, or adding labels in the three directions. The final resulting plot is shown below, with the blue points being the original data points:

It can be seen that all the real data points falls nicely on the surface described by the neural network, and the surface is smooth and monotonic as expected. We have therefore successfully fitted a curve using neural network and trained it to give the bandgap for any arbitrary barrier and well width.

---

Posted on Utopian.io - Rewarding Open Source Contributors

---