Milestone 6 Report

Deliver progress report and preliminary software for the hybrid models/methods.

The preliminary software for this project can be found in a Github repository at https://github.com/jabrams23/UoE-UW-DARPA. Since the delivery of Milestone 5 we have further progressed our analysis. We have also focused on improving documentation for our code. This report serves mainly as a walk through of the current codebase.

This project has several interlinked sub-projects, which is reflected in the folder structure of the Github repo. These folders are: AMOC, Ising Model, Kefi Models, Patterned vegetation and Tom_Models. Here we provide a description of this software, how it works so far and an outline for future developments.

AMOC

Analysis concerning the AMOC tipping point centers around the use of deep learning (DL) to predict the movement towards collapse in a 4-box model as detailed by Gnanadesikan et al.. The model code has been written in two ways, the first two of which 4box_AMOC_null_training_set.R and 4box_AMOC_tip_training_set.R randomly sample the parameters altered in the original paper from a uniform distribution. The code is set up to run an ensemble of 10000 members and we create two ensembles, one where the circulation is not approaching tipping (sampling a freshwater forcing value), and the other where is a movement towards tipping (sampling two freshwater forcing values). The script sw_dens_funcs.R is used in the model to calculate seawater density.

Previous analysis has shown that the DL models created by Bury et al. perform poorly at differentiating between those ensemble members that are forced slowly towards tipping, and those that remain stable. This has led us to create a DL model to determine how far away from tipping the tipping ensemble members are. Alongside each of the parameter combinations in the tipping dataset, a simulation is run with no additive noise, that continues to extrapolate the increase in freshwater forcing to determine when tipping will occur with that specific set of parameters. The script to do this (4box_AMOC_dist_from_tipping.R) takes the saved output of 4box_AMOC_tip_training_set.R as input. A caveat here is that this shows when the equilibrium run would collapse, with additive noise most likely causing an earlier tip in practice.

Generating training data

Each ensemble member runs for 1000 years to determine the equilibrium, then runs for 375 years with stochastic additive noise (of randomly sampled level) and the change in freshwater forcing (for those tipping ensemble members) added. Sometimes a member of the ensemble that uses the set of ‘AMOC off’ equations may tip due to the combination of parameters and tip in the equilibrium phase, or due to a combination of the freshwater forcing and additive noise level. In these cases, the ensemble member is rejected and a new set of parameters are sampled to run a new member.

Examining this tipping dataset showed a strongly exponential distribution of times until the ensemble members would tip, leading to the targets for the DL needing to be logged. To complement this, we also have created code (4box_AMOC_tip_training_set_ver2.R) to produce an ensemble using the standard parameter values. With these, we can determine the critical value of freshwater forcing needed to tip the model, from which we can sample the rate of change in freshwater forcing to create an ensemble member that is a required distance away from tipping. As such, we can create a uniform distribution of targets for the DL.

Training the neural network

A CNN-LSTM neural network coded in Tensorflow has been created to determine the distance from tipping in both cases. In this case we split the data into training (70%), validation (15%) and test (15%) sets. An example of the DL training can be found in CNN_LSTM_eg.R. A different number of CNN and LSTM layers have been tested to determine a best validation fit on the model ensemble created by the first instance detailed above. For the targets that form a uniform distribution, exploration is currently being undertaken. It is hoped that the results from this will inform a better DL model fit across both ensembles.

Ising Model

This folder contains code for a deep learning classifier designed to detect early warning signals in spatiotemporal systems undergoing phase transitions. Training data is generated from implementations of the 2D square Ising model. This data is then processed and fed into a CNN-LSTM neural network coded in Tensorflow. The resulting classifier model is tested on withheld Ising model simulations as well as two other test systems. Scripts to execute this pipeline are as follows:

Generating training data

Ising model data is generated by ising_model_generate_data.py, using functions from ising_model_simulate.py and Ising_model_process_data.py. Ising model lattices are simulated with randomized parameters for coupling strength, temperature bounds, and levels of coarse-graining. Options are given to simulate a system with varied temperature (‘temp’) or varied external magnetic field (‘h’), corresponding to realizations of second- and first-order phase transitions, respectively. In each case, there are some runs simulated in which the actuated quantity passes through its critical value/region and others where it does not. This provides an equal number of positive and null examples for the neural model to learn to distinguish. Functionality is also provided for the Ising lattice to be masked with randomized ellipses (i.e. random elliptical regions of the lattice deleted in order to break symmetry and lend greater complexity to the system). Coupling coefficients J between neighboring lattice sites are also randomized (normally distributed about a mean value). All runs are spatially coarse grained by a factor of at least 4 (up to 14) in order to smooth the binary state values of individual lattice sites. Further smoothing is accomplished by temporal coarse graining during simulation, in which site values are averaged over epochs of 10000 individual spin flip operations performed by the Metropolis algorithm.

Resulting spatiotemporal time series data is then processed to prepare it for model training. Randomized subregions of the lattice (at different levels of spatial coarse-graining) are selected and detrended with a Gaussian filter (along the time dimension). Secondary statistics of variance and correlation along the temporal and spatial axes are then computed and saved. Data is prepared in a Tensorflow-friendly format by ising_model_prepare_train_data.py.

Training neural network

Training of the neural model is handled by Ising_train_CNN_LSTM.py. This script reads in the prepared training data derived from the Ising model (containing phase transitions of either or both orders) and instantiates a CNN-LSTM neural network with any of a preset variety of combinations of hyperparameters (defined by a hardcoded table of CL_par_* variables). The model is trained using the standard Keras Tensorflow functionality, with the Adam optimization algorithm on categorical cross-entropy loss. Tests are carried out to evaluate the resulting classifier on a withheld test set of Ising data, and plots are produced to visualize both the training process and the performance on the test set.

Additional test systems

Two other spatiotemporal systems exhibiting phase transition behavior are provided as out-of-sample test cases for the model. A coupled water-vegetation model, based on one presented by Kefi et. al., is implemented in wv_generate_data.py. A model of our own design, intended to provide a simplified realization of percolation phase transitions in sea ice dynamics, is presented in perc_generate_data.py. Each of these files has a corresponding wv_process_EWS.py and perc_process_EWS.py file which applies the same processing pipeline used for the Ising model training data. These data sets are supplied to the trained neural model in spatiotemporal_test_cases.py, with classification results saved to the disk. These results are visualized in spatiotemporal_test_cases_plots.py, which produces plots of different ROC curves.

Patterned Vegetation

This folder contains the code to train a deep learning model to identify patterned vegetation and to distinguish it from non-patterned vegetation. It currently consists of three main components which include the code to model patterned vegetation which is dependent on precipitation (GeneratePattern.R), the implementation of this to generate training datasets of our 3 pattern clases (CreatePatternDatasets.R), and then a script to train and test a DL model (TrainDLModel.R). This code is written in the R coding language. Further details are provided below.

GeneratePattern.R

This script contains code to generate a 2-dimensional model of patterned vegetation, which first appeared in Konings et al. and is an adaptation of the Kefi et al.1 model. To generate a pattern, this model requires an initial, pre-defined JSON configuration file (patternConfig.json), which defines numerous variables, such as rainfall, bare soil infiltration, plant loss due to grazing, among others. In order to simulate across a range of rainfall values, we enable the user to define the rainfall level in the function. The grid size of the pattern can also be altered.

CreatePatternDatasets.R

This script contains code to create a dataset for our DL model which will later be split into a training, testing and validation dataset. This initial dataset contains 40,000 samples, which are divided evenly between four pattern classes: gaps, labyrinths, spots and no patterns.

Figure 1. Examples of gaps, labyrinths and spot vegetation as simulated from the model. The patterns are precipitation dependents and are generated with rainfall levels of 2.3 mm/day, 1.6 mm/day and 0.8 mm/day respectively.

Precipitation levels for each pattern are sampled from a normal distribution with a standard deviation of 0.1 and a mean of 0.8 mm/day for spots, 1.6 mm/day for labyrinths and 2.3 mm/day for gaps. These mean values were chosen by inspection, and examples of these patterns are given in Figure 1. We generate 10,000 iterations of each pattern on a 150 x 150 pixel grid. These are then converted to a binary image, where the value of 1 is given for vegetation and 0 is given for bare soil, this is in keeping with the pattern vegetation analysis in Buxton et al.. A set of non-patterned vegetation is also constructed, which aims to represent bare soil, complete vegetation coverage, patchy vegetation or other pattern structures which are not vegetation, such as roads.

TrainDLModel.R

This script shuffles the pattern vegetation datasets and the corresponding labels and splits them into training (80% of data), validation (10%) and test sets (10%). A 2-dimensional convolutional neural network is then trained to identify patterned vegetation using these datasets.

Preliminary results suggest that this approach works well for identifying vegetation patterning from models (>99% accuracy), however the model does not perform as well at identifying and classifying out of dataset satellite data. Further work is underway to diversify the training set, i.e. through introducing noise, constructing mixed images of a variety of patterns and non-patterns and with borderline/transitional patterns.

Tom_Models

This directory contains the models used in the analysis of Bury et al.. These are used in the above analysis for the AMOC models. In addition to this, these models are used to analyze the CMIP5 abrupt shifts from the Drijfhout et al. database, as a comparison exercise between conventional statistical indicators and DL techniques.