This notebook is intended to introduce and demonstrate some of the features of pyMAISE and examine the performance of machine learning models on a common nuclear engineering application. For further information on the capabilities of the classes and functions shown in this notebook, please refer to the pyMAISE API reference documentation.

MIT Reactor

Inputs: Control rod heights (\(cm\))

Outputs: Pin power (\(W\))

The MIT reactor data set represents the institution’s light-water-cooled 6 MW thermal power reactor. The figure below shows the core contains 22 fuel elements, 5 locations for in-core experiments, and six control blades surrounding them. The data set is used to find a relationship between the six control blade heights and the power produced by the 22 fuel elements in the core. Therefore, the data set is constructed by perturbing the depths of the control blades in the reactor—the corresponding output results in the power levels for each of the 22 fuel elements. The data was generated using MCNP, and the data set size includes 1000 simulations/samples [RPS+23]. The goal is to use pyMAISE to build, tune, and compare the performance of various ML models in predicting the core power distribution based on the control blade insertion depth.

MITR.png

The following are a few standard packages and functions that will prove helpful while using pyMAISE.

[1]:

%load_ext autoreload
%autoreload 2
#import functions as F

import time
import cv2
import numpy as np
import pandas as pd

# Set display option to show all rows and columns
pd.set_option('display.max_rows', None)
pd.set_option('display.max_columns', None)

# Set the width of the columns
pd.set_option('display.width', None)

# See the full content of each column
pd.set_option('display.max_colwidth', None)

import matplotlib.pyplot as plt
from scipy.stats import uniform, randint
from sklearn.preprocessing import MinMaxScaler

# Plot settings
matplotlib_settings = {
    "font.size": 12,
    "legend.fontsize": 11,
    "figure.figsize": (8, 8)
}
plt.rcParams.update(**matplotlib_settings)

We need to import several functions and classes for machine learning tuning and analysis with pyMAISE. Firstly, we need the MITR data set, which we can get using pyMAISE.datasets.load_MITR. We must split the data into training/testing data and scale it. For this, we use the pyMAISE.preprocessing Python module. The remaining classes we can get from pyMAISE directly which we import as mai for convenience.

[2]:

from pyMAISE.datasets import load_MITR
from pyMAISE.preprocessing import correlation_matrix, train_test_split, scale_data
import pyMAISE as mai

pyMAISE Initialization

Starting any pyMAISE job requires initialization. This includes the definition of global settings used throughout pyMAISE. These settings and their defaults include:

  • problem_type: the problem type, either regression or classification, defined by pyMAISE.ProblemType,

  • verbosity=0: the level of output from pyMAISE,

  • random_state=None: the seed for the random number generator, which can be used to get reproducible results from pyMAISE,

  • num_configs_saved=5: the number of top hyperparameter configurations for each model evaluated during tuning,

  • new_nn_architecture=True: a boolean that dictates whether to use the old deprecated pyMAISE neural network tuning architecture,

  • cuda_visible_devices=None: sets the CUDA_VISIBLE_DEVICES environment variable,

  • run_parallel=False: whether to run the NN tuner in parallel (only supported for GPUs),

  • max_models_per_device=numpy.inf: the maximum number of NNs training on a single GPU; if it’s numpy.inf, then pyMAISE will try to submit as many jobs as are available in memory on a given GPU.

The only argument that needs to be specified is problem_type. We also pass "-1" to cuda_visible_devices to ensure we only use tensorflow on the CPU. This is useful for this problem since we will build relatively simple, dense feedforward neural networks with a reasonably small data set. Therefore, running tensorflow on a GPU may hurt our performance. We leave the others default, giving us five hyperparameter configurations for each model, keeping the stochastic nature of some of the algorithms, and using the current neural network hyperparameter tuning architecture.

[3]:

global_settings = mai.init(
    problem_type=mai.ProblemType.REGRESSION,   # Define a regression problem
    cuda_visible_devices="-1"                  # Use CPU only
)

Data Loading and Preprocessing

pyMAISE has several benchmarked data sets, such as this MIT reactor data set. Each data set has an accompanying load function. Most load functions return three xarray.DataArrays: the raw, input, and output data. You can load personal data using the pyMAISE.preprocessing.read_csv function. Refer to the pyMAISE API reference documentation for specifics on the use of this function. Loading the MIT reactor data gives the following data:

[4]:

data, inputs, outputs = load_MITR()

The MIT reactor data set has six inputs; one for each control rod position:

[5]:

inputs

[5]:

<xarray.DataArray (index: 1000, variable: 6)>
array([[25.95991738, 22.94937226, 20.85317495, 24.66916763, 20.48104739,
        25.35726597],
       [21.75386762, 25.3606259 , 20.58853009, 20.11087249, 27.46711019,
        25.81658492],
       [27.42919944, 23.5701803 , 27.59630663, 26.39044508, 23.99603712,
        24.61182164],
       ...,
       [24.47247392, 23.21393049, 21.36224456, 25.72137934, 21.72606189,
        20.486756  ],
       [25.14502476, 25.88719097, 22.63841745, 23.11818585, 27.13820033,
        20.12257523],
       [20.80736327, 20.47244047, 26.3283514 , 24.93870275, 22.13608868,
        20.23749172]])
Coordinates:
  * index     (index) int64 0 1 2 3 4 5 6 7 ... 992 993 994 995 996 997 998 999
  * variable  (variable) object 'CR1' 'CR2' 'CR3' 'CR4' 'CR5' 'CR6'

and 22 fuel element power outputs:

[6]:

outputs

[6]:

<xarray.DataArray (index: 1000, variable: 22)>
array([[25930.9161377 , 22958.31494141, 21725.51635742, ...,
        17912.01452637, 18207.04260254, 19089.96942139],
       [25883.078125  , 22856.06195068, 21602.10876465, ...,
        18699.55557251, 18381.29052734, 19052.58782959],
       [25672.20825195, 22584.91094971, 21419.95025635, ...,
        17878.91418457, 17831.21044922, 18702.88995361],
       ...,
       [26011.28759766, 22868.43511963, 21718.16387939, ...,
        17537.16012573, 17521.70654297, 18469.28430176],
       [25867.42346191, 22746.65338135, 21650.78973389, ...,
        17898.8989563 , 17467.56048584, 18369.0256958 ],
       [26047.4777832 , 22729.4630127 , 21553.59606934, ...,
        17622.36529541, 17503.63659668, 18324.41955566]])
Coordinates:
  * index     (index) int64 0 1 2 3 4 5 6 7 ... 992 993 994 995 996 997 998 999
  * variable  (variable) object 'A-2' 'B-1' 'B-2' 'B-4' ... 'C-13' 'C-14' 'C-15'

To get a better understanding of this data set lets plot a correlation matrix of the data using pyMAISE.preprocessing.correlation_matrix.

[7]:

correlation_matrix(data)
plt.show()
../_images/benchmarks_mit_reactor_13_0.png

A negative correlation exists between the control rod positions and the inner fuel elements (A and B elements).

With the data loaded, we can now preprocess it. This includes splitting into training and testing data sets and scaling. For splitting, we use the pyMAISE.preprocessing.train_test_split function and define the fraction of data for testing using test_size. On MITR, we take 30% of the data for testing. The result is a tuple of data: (xtrain, xtest, ytrain, ytest). We can then scale this data using pyMAISE.preprocessing.scale_data with any object that supports fit_transform and transform functions. Common scalers include sklearn.preprocessing.MinMaxScaler and sklearn.preprocessing.StandardScaler. Many machine learning models learn best with scaled data. We use min-max scaling on both the input and output data for this data set. pyMAISE.preprocessing.scale_data returns the scaled data, and the scaler fits the data. These scalers are used in postprocessing to evaluate the performance of the models.

[8]:

xtrain, xtest, ytrain, ytest = train_test_split(data=[inputs, outputs], test_size=0.3)
xtrain, xtest, xscaler = scale_data(xtrain, xtest, scaler=MinMaxScaler())
ytrain, ytest, yscaler = scale_data(ytrain, ytest, scaler=MinMaxScaler())

Model Initialization

For model initialization and tuning, we use the pyMAISE.Tuner. pyMAISE supports both classical ML methods with scikit-learn and sequential neural networks with Keras. We strongly recommend you refer to pyMAISE.Tuner in the pyMAISE API reference documentation for a list of supported models, neural network layers, and other information crucial for optimal tuning. We can define all the models we wish to hyperparameter tune using a dictionary. For the MITR data set, we define the following regression models with their dictionary keys:

  • linear: "Linear",

  • ridge: "RD",

  • elastic net: "EN",

  • lasso: "Lasso",

  • decision tree: "DT",

  • random forest: "RF",

  • extra trees: "ET",

  • k-nearest neighbors: "KN",

  • gradient boosting: GB,

  • stacking: Stacking,

  • multi-output: MultiOutput,

  • feedforward neural network: "FNN".

pyMAISE uses the dictionary keys for classical models to determine which scikit-learn model you request. If the keys do not match any supported keys, then it is assumed to be a neural network. In the "models" key, we define these models in a list. We can then define the hyperparameters for classical models, which will remain constant throughout tuning. These parameters must only be defined if you want something different from the default. Refer to the model documentation in scikit-learn for hyperparameters and defaults.

Note that this is a multi-output problem. Thus, we need to define how we are going to extend. In pyMAISE, we do so using a dictionary key "multi_output to support models that don’t natively have multi-output capabilities (eg. GB and AB below). We also provided a specific multi-output model, which is demonstrated below for SVM. For Stacking, the user must pass through their estimators using the same dictionary keys as with the other model in the benchmark.

For the neural network model, we need to define the architecture, optimizer, compiling, and fitting parameters under the "structural_params", "optimizer", "compile_params", and "fitting_params" keys. These parameters include those that remain constant and change during tuning. To define hyperparameters for tuning, we use the pyMAISE.Int, pyMAISE.Float, pyMAISE.Choice, pyMAISE.Boolean, and pyMAISE.Fixed classes. These classes are wrappers for KerasTuner.

To define the model architecture, we use the supported Keras layers defined in the pyMAISE API reference documentation under the pyMAISE.Tuner. These layers must differ in name, but use the base layer name. Within the layers, pyMAISE supports additional hyperparameters, which include "sublayer", "wrapper", and "num_layers" to define sublayers, wrappers (such as keras.layers.TimeDistributed), and the number of layers.

We can then pass all this information to pyMAISE.Tuner.

[9]:

model_settings = {
    "models": ["Linear", "Lasso", "Stacking", "MultiOutput", "GB", "DT", "FNN", "RD", "ET", "EN", "AB", "RF", "KN"],
    "FNN": {
        "structural_params": {
            "Dense_hidden": {
                "num_layers": mai.Int(min_value=0, max_value=3),
                "units": mai.Int(min_value=25, max_value=400),
                "activation": "relu",
                "kernel_initializer": "normal",
                "sublayer": mai.Choice(["Dropout_hidden", "None"]),
                "Dropout_hidden": {
                    "rate": mai.Float(min_value=0.4, max_value=0.6),
                },
            },
            "Dense_output": {
                "units": ytrain.shape[-1],
                "activation": "linear",
                "kernel_initializer": "normal",
            },
        },
        "optimizer": "Adam",
        "Adam": {
            "learning_rate": mai.Float(min_value=1e-5, max_value=0.001),
        },
        "compile_params": {
            "loss": "mean_absolute_error",
            "metrics": ["mean_absolute_error"],
        },
        "fitting_params": {
            "batch_size": mai.Choice([8, 16, 32]),
            "epochs": 50,
            "validation_split": 0.15,
        },
    },
    "GB":{
        "multi_output": True
    },
    "AB":{
        "multi_output": True
    },
    "MultiOutput":{
        "estimator": "SVM"
    },
    "Stacking":{
        "estimators": ["EN", "DT"],
        "_multi_output": True,
    },
}
tuner = mai.Tuner(xtrain, ytrain, model_settings=model_settings)

Hyperparameter Tuning

Hyperparameter tuning is split into two options: classical model and neural network tuning. For the classical models, pyMAISE supports grid, random, and Bayesian search. For the neural network models, pyMAISE supports grid, random, Bayesian, and hyperband search. All search methods use cross-validation.

We define the hyperparameter search space for the classical models by defining the arrays, distributions, or skopt.space.space for each hyperparameter we plan to tune. We describe these within subdictionaries within the model keys. For MITR, we use pyMAISE.Tuner.random_search with 200 iterations as the data set is relatively small and the classical models are computationally cheap. For these models, we define the search space using scipy.stats.

We do not hyperparameter tune linear regression because the hyperparameter options are limited, and the default scikit-learn model performs well on linear problems.

Note that for specific models below, we use estimator__ to get access to model parameters. This is because when we use any ensemble methods (eg. multi-output/stacking), you need to enter into that estimator to acess the parameters. For stacking, you must also specific which model you are trying to access with {model key}. Note that here we have estimator__{model key} since we have multi-output ontop of stacking.

Since we have already defined the feedforward neural network search space in the previous section, all we need to do is call the search method function. For these NNs, we can call pyMAISE.Tuner.nn_bayesian_search with 50 iterations as NNs tend to be more computationally expensive to train over the classical models, and the Bayesian optimization may converge on the optimal model for fewer iterations.

We set cv for both methods to 5, which gives us five cross-validation splits for each model. This will avoid overfitting. Therefore, for each classical model, we run 1000 hyperparameter configurations, and for the neural networks, we run 250 hyperparameter configurations.

[10]:

random_search_spaces = {
    "Lasso": {
        "alpha": uniform(loc=0.0001, scale=0.0099), # 0.0001 - 0.01
    },
    "DT": {
        "max_depth": randint(low=5, high=50), # 5 - 50
        "max_features": [None, "sqrt", "log2", 2, 4, 6],
        "min_samples_split": randint(low=2, high=20), # 2 - 20
        "min_samples_leaf": randint(low=1, high=20), # 1 - 20
    },
    "RF": {
        "n_estimators": randint(low=50, high=200), # 50 - 200
        "criterion": ["squared_error", "absolute_error", "poisson"],
        "min_samples_split": randint(low=2, high=20), # 2 - 20
        "min_samples_leaf": randint(low=1, high=20), # 1 - 20
        "max_features": [None, "sqrt", "log2", 2, 4, 6],
    },
    "KN": {
        "n_neighbors": randint(low=1, high=20), # 1 - 20
        "weights": ["uniform", "distance"],
        "leaf_size": randint(low=1, high=30), # 1 - 30
        "p": randint(low=1, high=10), # 1 - 10
    },
    "RD":{
        "alpha": uniform(loc=1e-11, scale=1e-9)
    },
    "EN": {
        "alpha": uniform(loc=1e-11, scale=1e-9)
    },
    "ET":{
        "max_depth": randint(low=5, high=50), # 5 - 50
        "min_samples_split": randint(low=2, high=20), # 2 - 20
        "min_samples_leaf": randint(low=1, high=20), # 1 - 20
    },
    "GB": {
        "estimator__learning_rate": uniform(loc=0.1, scale=1e-2)
    },
    "AB":{
        "estimator__learning_rate": uniform(loc=1, scale=0.5)
    },
    "MultiOutput":{
        "estimator__C": uniform(loc=0.01, scale=99.99), # 0.01 - 100
    },
    "Stacking":{
        "estimator__EN__alpha": uniform(loc=1e-11, scale=1e-9),
        "estimator__DT__max_depth": randint(low=5, high=50), # 5 - 50
        "estimator__DT__max_features": [None, "sqrt", "log2", 2, 4, 6],
        "estimator__DT__min_samples_split": randint(low=2, high=20), # 2 - 20
        "estimator__DT__min_samples_leaf": randint(low=1, high=20), # 1 - 20
    },
}
start = time.time()
random_search_configs = tuner.random_search(
    param_spaces=random_search_spaces,
    n_iter=200,
    n_jobs=6,
    cv=5,
)
bayesian_search_configs = tuner.nn_bayesian_search(
    objective="r2_score",
    max_trials=50,
    cv=5,
)
print("Hyperparameter tuning took " + str((time.time() - start) / 60) + " minutes to process.")

Hyperparameter tuning took 23.69085413614909 minutes to process.

After the training, we can see training results for each iteration using the pyMAISE.Tuner.convergence_plot function. For the neural networks, we have the following:

[20]:

ax = tuner.convergence_plot(model_types="FNN")
ax.set_ylim([0, 1])
plt.show()
../_images/benchmarks_mit_reactor_21_0.png

The Bayesian search explores the search space but finds several suitable hyperparameter configurations.

Model Postprocessing

With the models tuned and the top pyMAISE.Settings.num_configs_saved saved, we can now pass these models to the pyMAISE.PostProcessor for model comparison and analysis. Additionally, we can update some hyperparameters after tuning. For MITR, we changed the number of FNN epochs to 200 since we are only running five models. We give the pyMAISE.PostProcessor all the scaled data, the model configurations, the new model settings, and the output scaler.

[21]:

postprocessor = mai.PostProcessor(
    data=(xtrain, xtest, ytrain, ytest),
    model_configs=[random_search_configs, bayesian_search_configs],
    new_model_settings={
        "FNN": {"fitting_params": {"epochs": 200}},
    },
    yscaler=yscaler,
)

We can now evaluate the performance of the models using the pyMAISE.PostProcessor.metrics function. This returns an ordered table with training and testing performance metrics. By default, pyMAISE.PostProcessor.metrics evaluates the training and testing \(R^2\), mean absolute error (MAE), mean absolute percentage error (MAPE), root mean squared error (RMSE), and root mean squared percentage error (RMSPE) for regression problems. The table is sorted by descending testing \(R^2\). We can also add additional metrics, sort by a different metric, or choose the output pyMAISE computes the metrics with. Since all outputs are in watts, we do not need to evaluate each output individually. For further functionality on pyMAISE.PostProcessor.metrics refer to the pyMAISE API reference documentation.

[22]:

postprocessor.metrics().drop("Parameter Configurations", axis=1)

[22]:
Model Types Train R2 Train MAE Train MAPE Train RMSE Train RMSPE Test R2 Test MAE Test MAPE Test RMSE Test RMSPE
57 FNN 0.998522 7.651470 0.040188 10.990297 0.052270 0.998321 8.178382 0.042891 11.580021 0.055818
60 FNN 0.998508 7.014305 0.036585 10.486746 0.051757 0.998316 7.574305 0.039428 11.020411 0.054630
56 FNN 0.998393 7.634138 0.039536 10.204548 0.049374 0.998201 8.314544 0.043119 11.126426 0.054253
59 FNN 0.998328 7.755986 0.040237 10.704594 0.051773 0.998101 8.419915 0.043704 11.552209 0.056252
58 FNN 0.998166 8.086968 0.041680 10.607666 0.052269 0.997801 9.016867 0.046574 11.829436 0.058491
0 Linear 0.994589 12.895922 0.067056 18.689971 0.093390 0.995613 12.226755 0.063559 17.099414 0.085183
30 RD 0.994589 12.895922 0.067056 18.689971 0.093390 0.995613 12.226755 0.063559 17.099414 0.085183
29 RD 0.994589 12.895922 0.067056 18.689971 0.093390 0.995613 12.226755 0.063559 17.099414 0.085183
27 RD 0.994589 12.895922 0.067056 18.689971 0.093390 0.995613 12.226755 0.063559 17.099414 0.085183
26 RD 0.994589 12.895922 0.067056 18.689971 0.093390 0.995613 12.226755 0.063559 17.099414 0.085183
28 RD 0.994589 12.895922 0.067056 18.689971 0.093390 0.995613 12.226755 0.063559 17.099414 0.085183
39 EN 0.994589 12.895922 0.067056 18.689971 0.093390 0.995613 12.226755 0.063559 17.099414 0.085183
40 EN 0.994589 12.895922 0.067056 18.689971 0.093390 0.995613 12.226755 0.063559 17.099414 0.085183
38 EN 0.994589 12.895922 0.067056 18.689971 0.093390 0.995613 12.226755 0.063559 17.099414 0.085183
37 EN 0.994589 12.895922 0.067056 18.689971 0.093390 0.995613 12.226755 0.063559 17.099414 0.085183
36 EN 0.994589 12.895922 0.067056 18.689971 0.093390 0.995613 12.226755 0.063559 17.099414 0.085183
1 Lasso 0.994566 13.064965 0.067923 18.734206 0.093580 0.995543 12.419661 0.064541 17.230611 0.085796
2 Lasso 0.994524 13.215825 0.068701 18.814318 0.093962 0.995469 12.593155 0.065434 17.374703 0.086490
3 Lasso 0.994367 13.629701 0.070847 19.112596 0.095419 0.995243 13.044244 0.067763 17.817935 0.088656
4 Lasso 0.994208 13.972201 0.072627 19.408889 0.096881 0.995038 13.417249 0.069694 18.217272 0.090624
5 Lasso 0.994113 14.162937 0.073619 19.585710 0.097756 0.994918 13.624632 0.070770 18.446925 0.091759
9 Stacking 0.993333 15.427020 0.080379 21.221407 0.105473 0.993089 16.241389 0.084547 21.709552 0.107421
7 Stacking 0.993185 15.580625 0.080994 21.191958 0.105642 0.993064 16.311708 0.084765 21.520878 0.107160
10 Stacking 0.993168 15.612716 0.081260 21.368638 0.106510 0.992609 16.860849 0.087705 22.254331 0.110713
8 Stacking 0.993084 15.895874 0.082805 21.712463 0.107736 0.992456 17.094955 0.088988 22.847127 0.112987
6 Stacking 0.992452 16.769559 0.087348 22.670303 0.112780 0.991845 17.900181 0.093148 23.654289 0.117097
18 GB 0.996414 11.652970 0.060569 15.410644 0.076915 0.979806 28.037931 0.145550 36.429789 0.181778
17 GB 0.996419 11.641918 0.060503 15.389189 0.076796 0.979779 28.053830 0.145631 36.425126 0.181777
19 GB 0.996428 11.633960 0.060466 15.364124 0.076670 0.979773 28.055968 0.145653 36.466874 0.181903
16 GB 0.996420 11.639010 0.060487 15.382631 0.076799 0.979767 28.056177 0.145653 36.473334 0.181956
20 GB 0.996433 11.662542 0.060624 15.426311 0.076829 0.979521 28.268099 0.146710 36.747966 0.182866
53 KN 1.000000 0.000000 0.000000 0.000000 0.000000 0.941764 50.280813 0.263252 66.749644 0.330709
51 KN 1.000000 0.000000 0.000000 0.000000 0.000000 0.941513 50.674117 0.265386 67.026931 0.331908
52 KN 1.000000 0.000000 0.000000 0.000000 0.000000 0.941218 50.897468 0.266673 67.375957 0.333413
55 KN 1.000000 0.000000 0.000000 0.000000 0.000000 0.941099 50.342355 0.263431 66.963304 0.331805
13 MultiOutput 0.943199 49.330086 0.257756 64.754166 0.321625 0.938272 49.116152 0.256356 66.778377 0.332670
12 MultiOutput 0.943199 49.330086 0.257756 64.754166 0.321625 0.938272 49.116152 0.256356 66.778377 0.332670
15 MultiOutput 0.943199 49.330086 0.257756 64.754166 0.321625 0.938272 49.116152 0.256356 66.778377 0.332670
14 MultiOutput 0.943199 49.330086 0.257756 64.754166 0.321625 0.938272 49.116152 0.256356 66.778377 0.332670
11 MultiOutput 0.943199 49.330086 0.257756 64.754166 0.321625 0.938272 49.116152 0.256356 66.778377 0.332670
54 KN 1.000000 0.000000 0.000000 0.000000 0.000000 0.937977 52.551795 0.275380 69.226016 0.342367
31 ET 0.995062 14.486439 0.075838 19.467468 0.096187 0.934138 51.232586 0.267875 70.578198 0.349732
32 ET 0.995230 14.237860 0.074581 19.132383 0.094537 0.933101 51.899850 0.271556 71.262238 0.353028
34 ET 0.984064 25.068874 0.131177 34.774141 0.172325 0.925610 54.291129 0.284014 75.000859 0.372073
33 ET 0.983971 25.175986 0.131659 34.884440 0.172725 0.925002 54.750377 0.286241 75.314955 0.372986
35 ET 0.980793 28.152816 0.147287 38.419301 0.189901 0.923701 55.579070 0.290643 76.039966 0.376509
43 AB 0.930691 53.138599 0.275749 66.690694 0.333620 0.907799 60.139175 0.312135 77.689743 0.387559
44 AB 0.929857 53.569285 0.278074 67.381991 0.336612 0.906902 60.698030 0.315230 78.364714 0.390649
45 AB 0.930191 53.577110 0.278190 67.384974 0.336485 0.906518 60.794817 0.315599 78.486415 0.391348
41 AB 0.930494 53.165917 0.275953 66.736017 0.333969 0.906339 60.375576 0.313350 77.993015 0.389640
42 AB 0.929634 53.591065 0.278197 67.321103 0.336413 0.905960 60.561209 0.314278 78.411401 0.391092
47 RF 0.962164 40.484389 0.211932 53.899722 0.266339 0.884928 71.136447 0.372161 93.766269 0.463894
46 RF 0.981332 27.860081 0.145888 37.839969 0.186820 0.884815 70.002958 0.365806 93.259719 0.460897
48 RF 0.932460 54.016528 0.282515 71.968080 0.355049 0.868691 75.083263 0.392446 99.518885 0.492844
49 RF 0.930966 54.842553 0.287011 72.954149 0.360228 0.865970 76.409399 0.399578 101.213943 0.501018
50 RF 0.931283 54.857252 0.287130 72.814569 0.359563 0.864165 77.082358 0.403241 101.910463 0.504550
22 DT 0.903526 64.734217 0.339298 86.170114 0.426064 0.720725 110.630982 0.579576 146.601384 0.724877
21 DT 0.879990 72.471361 0.379819 96.091326 0.475046 0.715738 112.062329 0.587186 147.475184 0.730209
25 DT 0.875137 74.045449 0.387996 97.909219 0.483931 0.714681 112.972956 0.592017 148.488063 0.734188
24 DT 0.875137 74.045449 0.387996 97.909219 0.483931 0.713980 113.097348 0.592596 148.724465 0.735041
23 DT 0.850537 81.763843 0.428436 107.331722 0.530457 0.712929 114.007295 0.597593 148.918379 0.736472

This data set is very linear, with fantastic performance from linear, elastic-net, ridge, and lasso regression. The FNNs proved to be marginally better-performing models over the linear classical models; however, this is at the cost of complexity. K-nearest neighbor, extra-trees, adaboost, random forest, and decision tree struggled on this data set with all three overfittings.

Using the pyMAISE.PostProcessor.print_model function, we can see all the top performing (based on test \(R^2\)) models’ hyperparameters.

[14]:

for model in ["Lasso", "DT", "FNN", "RD", "ET", "EN", "GB", "AB", "MultiOutput", "Stacking", "RF", "KN"]:
    postprocessor.print_model(model_type=model)
    print()

Model Type: Lasso
  alpha: 0.0001133351011862626

Model Type: DT
  max_depth: 43
  max_features: 6
  min_samples_leaf: 3
  min_samples_split: 6

Model Type: FNN
  Structural Hyperparameters
    Layer: Dense_hidden_0
      units: 309
      sublayer: None
    Layer: Dense_output_0
  Compile/Fitting Hyperparameters
    Adam_learning_rate: 0.0008321972582830564
    batch_size: 8
Model: "FNN"
_________________________________________________________________
 Layer (type)                Output Shape              Param #
=================================================================
 Dense_hidden_0 (Dense)      (None, 309)               2163

 Dense_output_0 (Dense)      (None, 22)                6820

=================================================================
Total params: 8983 (35.09 KB)
Trainable params: 8983 (35.09 KB)
Non-trainable params: 0 (0.00 Byte)
_________________________________________________________________

Model Type: RD
  alpha: 9.643960473209805e-10

Model Type: ET
  max_depth: 33
  min_samples_leaf: 1
  min_samples_split: 3

Model Type: EN
  alpha: 3.724600700139202e-11

Model Type: GB
  estimator__learning_rate: 0.10929699916511056

Model Type: AB
  estimator__learning_rate: 1.47118944569141

Model Type: MultiOutput
  estimator__C: 46.16762541225824

Model Type: Stacking
  estimator__DT__max_depth: 9
  estimator__DT__max_features: 2
  estimator__DT__min_samples_leaf: 16
  estimator__DT__min_samples_split: 7
  estimator__EN__alpha: 4.850855576743397e-10

Model Type: RF
  criterion: absolute_error
  max_features: None
  min_samples_leaf: 1
  min_samples_split: 4
  n_estimators: 186

Model Type: KN
  leaf_size: 4
  n_neighbors: 8
  p: 2
  weights: distance

Additionally, we can visualize the structure of the FNN using pyMAISE.PostProcessor.nn_network_plot.

[15]:

postprocessor.nn_network_plot(
    to_file="./supporting/mit_reactor_network.png",
    show_shapes=True,
    show_layer_names=True,
    expand_nested=True,
    show_layer_activations=True,
)

[15]:
../_images/benchmarks_mit_reactor_29_0.png

We can use the pyMAISE.PostProcessor.diagonal_validation_plot and pyMAISE.PostProcessor.validation_plot functions to visualize the performance of these models. Diagonal validation plots show the actual versus predicted outcome. A well-performing model follows \(y=x\) on this plot. You can choose what outputs to display on the diagonal validation plot; however, since all outputs have the same units, we can display all of them. Both functions use the top-performing model on test \(R^2\) by default.

[16]:

models = np.array([["Linear", "Lasso"], ["RD", "FNN"], ["EN", "ET"], ["AB", "MultiOutput"], ["Stacking", "GB"], ["RF", "KN"], ["RD", "RD"]] )
fig, axarr = plt.subplots(models.shape[0], models.shape[1], sharex=True, sharey=True, figsize=(15,40))
for i in range(models.shape[0]):
    for j in range(models.shape[1]):
        plt.sca(axarr[i, j])
        axarr[i, j] = postprocessor.diagonal_validation_plot(model_type=models[i, j])
        axarr[i, j].set_title(models[i, j])

plt.show()
../_images/benchmarks_mit_reactor_31_0.png

All models display a close spread to \(y=x\); however, the lack of performance of the decision tree, extra-trees, random forest, and k-nearest neighbor is apparent, given the more extensive spread.

Validation plots show the absolute error of the predicted output relative to the actual outputs of the testing data set. The function can evaluate all output or show just what is given in a list. This list can include column positions in the data set or the output names. Only the A-2, B-8, and C-8 fuel elements are shown for the validation plot below.

[17]:

fig, axarr = plt.subplots(models.shape[0], models.shape[1], figsize=(17,22))

y = ["A-2", "B-8", "C-8"]

for i in range(models.shape[0]):
    for j in range(models.shape[1]):
        plt.sca(axarr[i, j])
        axarr[i, j] = postprocessor.validation_plot(model_type=models[i, j], y=y)
        axarr[i, j].set_title(models[i, j])
        axarr[i, j].get_legend().remove()

plt.subplots_adjust(hspace=0.5)  # Adjust hspace to increase spacing between rows
fig.legend(y, loc="upper center", ncol=4)
plt.show()
../_images/benchmarks_mit_reactor_33_0.png

The diagonal validation and validation plots agree with the performance metrics. Interestingly, all models increase prediction error as we move radially from the middle of the core.

To further understand the behavior of the top neural network configurations, we can plot the learning curve. Here, the top neural network learning curve is shown but, similar to the diagonal and validation plot functions, pyMAISE.PostProcessor.nn_learning_plot shows the neural network based on the index in pyMAISE.PostProcessor.metrics or, if no index is provided, the one with the best test \(R^2\).

[18]:

postprocessor.nn_learning_plot()
plt.show()
../_images/benchmarks_mit_reactor_35_0.png

The FNN is not overfitting, as the validation curve closely follows the training curve.

Finally, using the best FNN model, we can generate a predicted power for each fuel element given the control blade heights. This code block retrieves the top-performing model using the pyMAISE.PostProcessor.get_model function.

[19]:

def plot_mitr(model, x, yscaler=None):
    pos=[(400,330), (465,230), (500,300), (395,400), (320,400), (285,255), (318,193),
         (500,165), (535,230), (575,295), (535,360), (500,430), (430,465), (355,465),
         (280,465), (240,400), (205,335), (210,255), (243,193), (280,130), (355,130),
         (430,130), (393,193), (460,360), (280,335), (400,265), (340,295)]

    Ynn=model.predict(np.array([x,]))
    if yscaler:
        Ynn=yscaler.inverse_transform(Ynn)
    Ynn=Ynn.flatten().tolist()

    image = cv2.imread("./supporting/mitr.png")
    for i in range(len(pos)):
        if i==0:
            image=cv2.putText(
                img=np.copy(image),
                text=str(int(Ynn[i])),
                org=pos[i],
                fontFace=1,
                fontScale=1.1,
                color=(0,0,0)
            )
        if i in [1,2,3,4,5,6]:
            image=cv2.putText(
                img=np.copy(image),
                text=str(int(Ynn[i])),
                org=pos[i],
                fontFace=1,
                fontScale=1.1,
                color=(0,0,255)
            )
        if i in list(range(7,22)):
            image=cv2.putText(
                img=np.copy(image),
                text=str(int(Ynn[i])),
                org=pos[i],
                fontFace=1,
                fontScale=1.1,
                color=(0, 100, 0)
            )
        if i in list(range(22,28)):
            image=cv2.putText(
                img=np.copy(image),
                text=str('Empty'),
                org=pos[i],
                fontFace=1,
                fontScale=1.1,
                color=(233, 150, 122)
            )
    plt.figure(figsize=(12, 12))
    plt.imshow(image)
    plt.axis('off')

x=[0.75266553, 0.90280633, 0.00539489, 0.25308624, 0.57678792, 0.77792903]

plot_mitr(postprocessor.get_model(), x, yscaler)

1/1 [==============================] - 0s 20ms/step
../_images/benchmarks_mit_reactor_37_1.png

pyMAISElogo.png