<jupyter_start><jupyter_text>Entire workflow for 32 column run In this notebook, I will try to document every step currently required from the raw SP training data to a neural network online run with analysis. The hope is 1) to provide an introduction for people new to this code and (more importantly) 2) show which bottlenecks there currently are in the workflow and where things could/should be improved.Unfortunately, the code and the repository are quite messy at times. Sorry about that.We will use the 32 column SPCAM training run, which has previously crashed, which gives us a good opportunity to do some analysis for a crashed run.I will start working on my local machine in Munich, where I copied the raw data to. There I will do data preprocessing and network training. Then we will switch to Stampede2 to compile CAM, run the neural network version and analyze the output data.<jupyter_code># Imports
import xarray as xr
!pwd
%cd ../../<jupyter_output>/home/s/S.Rasp/repositories/CBRAIN-CAM<jupyter_text>SummaryHere are the commands necessary for the entire workflow. Part 1 - LMU**Preprocessing**Create config file (e.g. `32col_engy_ess_ref.yml`):```inputs : [TBP, QBP, VBP, PS, SOLIN, SHFLX, LHFLX]outputs : [TPHYSTND, PHQ, FSNT, FSNS, FLNT, FLNS, PRECT]in_dir : /project/meteo/w2w/A6/S.Rasp/SP-CAM/sp32fbp_andkua/out_dir: /local/S.Rasp/preprocessed_data/```Run preprocessing scripts:```python cbrain/preprocess_aqua.py \--config pp_config/32col_engy_ess_ref.yml \--aqua_names '*.h1.0000-*-0[1-3]-*' \--out_pref 32_col_engy_ess_3d_trainpython cbrain/preprocess_aqua.py \--config pp_config/32col_engy_ess_ref.yml \--aqua_names '*.h1.0000-*-1[1-3]-*' \--out_pref 32_col_engy_ess_3d_valid --ext_norm Nopepython cbrain/shuffle_ds.py --pref /local/S.Rasp/preprocessed_data/32_col_engy_ess_3d_train```**Neural network training**Create config file (e.g. `001_32cols_engy_ess_3d_pnas.yml`):```exp_name: 001_32cols_engy_ess_3d_pnasdata_dir: /local/S.Rasp/preprocessed_data/train_fn: 32_col_engy_ess_3d_train_shufflevalid_fn: 32_col_engy_ess_3d_validnorm_fn: 32_col_engy_ess_3d_train_norm.ncfsub: feature_meansfdiv: max_rstmult: target_convactivation: LeakyReLUhidden_layers: [256,256,256,256,256,256,256,256,256]loss: mselog_dir: ./logs/epochs: 18lr_step: 3```Run NN script```python run_experiment.py -c nn_config/32cols/001_32cols_engy_ess_3d_pnas.yml```Convert to text files```python save_weights.py -c nn_config/32cols/001_32cols_engy_ess_3d_pnas.yml``` Part 2 - StampedeCopy text files from LMU to Stampede:```weights_from_lmu 001_32cols_engy_ess_3d_pnas```Change network architecture in `cloudbrain_keras_dense.F90`. Additionally, make other changes before compiling.```build```Copy to scratch and create a new run directory```cdsmkdir nncam_001_32cols_engy_ess_3d_pnascd nncam_001_32cols_engy_ess_3d_pnas/```Copy cam binary, network parameters, initial condition file and template namelist and batch files.```cp /home1/05488/tg847872/spcam3-nn-build/run/cam ./cp -r /home1/05488/tg847872/saved_models/001_32cols_engy_ess_3d_pnas ./keras_matricescp /work/05488/tg847872/stampede2/ic_bc/spinup_AndKua_aqua_SPCAM3.0.cam2.i.0000-12-02-00000.nc ./cp /from/somewhere/atm_in /from/somewhere/run.slurm ./```Edit `atm_in`(`caseid`, potentially output, run length) and `run.slurm` (batch ID, queue and **directory**.)Start interactive session:```interactive```Run NNCAM:```bash run.slurm >> out.txt &```Analyze the output data. Would it be useful to write some quick-check functions?Copy the data to somewhere so that they are safe.  Preprocessing: from raw data to neural network inputThe goal of the preprocessing is to convert the raw SPCAM output data to a format which can easily and quickly be used during neural network training.  The raw dataThe raw data is saved in daily files, each approximately 900MB large (h1 files contain the necessary data). For this 32 column version, there is 1 year and 3 months of data. For storage information of raw files, check https://github.com/raspstephan/CBRAIN-CAM/wiki/Raw-SPCAM-runs-(to-be-used-for-network-training)<jupyter_code>#RAW_DATA = '/project/meteo/w2w/A6/S.Rasp/SP-CAM/sp32fbp_andkua/'
RAW_DATA = '/local/S.Rasp/sp32fbp_andkua/'   #SSD drive<jupyter_output><empty_output><jupyter_text>! executes a command line script from within the notebook!<jupyter_code>!ls -hl $RAW_DATA*h1* | head -5
ds = xr.open_dataset(RAW_DATA + 'AndKua_aqua_SPCAM3.0_sp_fbp32.cam2.h1.0000-01-01-00000.nc', 
                     decode_times=False)   # This is necessary because of the way time is stored in SPCAM
ds<jupyter_output><empty_output><jupyter_text>Preprocessing scriptThe preprocessing is currently handled by a script called `preprocess_aqua.py` in the `cbrain` directory of the repository. To train the neural network we look at the data column by column. If we want to e.g. use TAP and QAP for the input of the neural network, in the raw data these two variables would be saved at different locations of the disc, which makes reading slow. In order to quickly read the data during neural network training, this script concatenates the input and output vector and saves these vectors to disc. A second step done by the script is to compute variables, that are not directly available in the raw output. For example, `TBP` = `TAP - TPHYSTND * dt`. The exact meaning of the variables for SPCAM can be found here (...eventually): https://github.com/raspstephan/CBRAIN-CAM/wiki/What-exactly-happens-in-a-SPCAM-time-step%3FThe preprocessing script takes a number of arguments from the command line or from a configuration script. Let's write a configuration script.<jupyter_code>!python cbrain/preprocess_aqua.py -h
!cat pp_config/32col_engy_ess_ref.yml
PREPROC_DIR = '/local/S.Rasp/preprocessed_data/'<jupyter_output><empty_output><jupyter_text>The most important arguments are:- inputs: List of input variables, e.g. [TBP, QBP, VBP, PS, SOLIN, SHFLX, LHFLX]. `_BP` variables are computed as: `TAP - TPHYSTND*dt`- outputs: List of output variables, e.g. [TPHYSTND, PHQ, FSNT, FSNS, FLNT, FLNS, PRECT]- in_dir: Directory where SPCAM files are stored. I strongly recommend putting the files on a fast drive (e.g. /scratch/ on Greenplanet)- aqua_names: Here you can put a string with placeholders `*`, e.g. `'*.h1.0001-*-*-*'` for all files from year 1 or `'*.h1.0000-*-1[7-9]-*'` for day 17-19 from each month of year 0.- out_dir: directory where the processed files will be stored- out_pref: Prefix for the output files. Give them some descriptive name, e.g. `my_experiment1`- ext_norm: This is deprecated and ugly. If ext_norm = None (default), a normalization file will be computed from the data. This normalization file contains means, standard deviations, etc. For a full year of data this can take a long time. For this reason I chose to compute the normalization files for a sample of the data. The differences are small. If ext_norm is some string (just pick any jiberish) no normalization file will be computed but the data will still not be normalized. This is controlled by the `norm_features` and `norm_targets` options which are None and should stay so. The normalization happens later on the fly during network training.Let's start by creating a training dataset using the first three days from each month of the first year of available data. This is a total of 36 days, so roughly a tenth of the entire year. I added timestamps to the script, so that we can see how much time each step in the preprocessing takes.<jupyter_code>!python cbrain/preprocess_aqua.py \
--config pp_config/32col_engy_ess_ref.yml \
--aqua_names '*.h1.0000-*-0[1-3]-*' \
--out_pref 32_col_engy_ess_3d_train<jupyter_output>Time checkpoint reading data: 13.15 s
Number of time steps: 1728
Cut time steps: [ 143  287  431  575  719  863 1007 1151 1295 1439 1583]
Cut time steps: [ 143  287  431  575  719  863 1007 1151 1295 1439 1583]
Time checkpoint create datasets: 13.28 s
Time checkpoint reshape and rechunk: 13.75 s
Compute means and stds
Saving normalization file: /local/S.Rasp/preprocessed_data/32_col_engy_ess_3d_train_norm.nc
Time checkpoint normalization arrays: 214.54 s
Time checkpoint rechunk and ds: 215.20 s
Save features: /local/S.Rasp/preprocessed_data/32_col_engy_ess_3d_train_features.nc
Save targets: /local/S.Rasp/preprocessed_data/32_col_engy_ess_3d_train_targets.nc
Total time: 292.96 s<jupyter_text>Computing the normalization file actually takes longest. This involves computing statistics over the entire sample dimension. But it's not crucial to compute the exactly right statistics for normalization. E.g. when using the entire year as training data, the normalization file computed with a subset of the data is totally fine to use.<jupyter_code>!ls $PREPROC_DIR/32_col*<jupyter_output>/local/S.Rasp/preprocessed_data//32_col_engy_ess_3d_train_features.nc
/local/S.Rasp/preprocessed_data//32_col_engy_ess_3d_train_norm.nc
/local/S.Rasp/preprocessed_data//32_col_engy_ess_3d_train_targets.nc<jupyter_text>Now we have a features and a target dataset plus a normalization dataset. Let's check all of them out.<jupyter_code>xr.open_dataset(PREPROC_DIR + '32_col_engy_ess_3d_train_features.nc')
xr.open_dataset(PREPROC_DIR + '32_col_engy_ess_3d_train_targets.nc')<jupyter_output><empty_output><jupyter_text>We have two dimensions, `sample` and `target_lev`, which corresponds to the stacked vectors. Additionally we have the variable `target_names`, which is handy later.<jupyter_code>xr.open_dataset(PREPROC_DIR + '32_col_engy_ess_3d_train_targets.nc').target_names
xr.open_dataset(PREPROC_DIR + '32_col_engy_ess_3d_train_norm.nc')<jupyter_output><empty_output><jupyter_text>The normalization file contains a number of statistics for the features and targets, which can then be used during training for normalization. We would also like to have a separate validation dataset. For this we simply chose a different date range, the months available for the second year. For the validation dataset we do not need to compute a normalization file.<jupyter_code>!python cbrain/preprocess_aqua.py \
--config pp_config/32col_engy_ess_ref.yml \
--aqua_names '*.h1.0001-*-0[1-3]-*' \
--out_pref 32_col_engy_ess_3d_valid --ext_norm Nope
!ls -lh $PREPROC_DIR/32_col*<jupyter_output>-rw-r--r-- 1 S.Rasp ls-craig 5.3G Jan  8 12:49 /local/S.Rasp/preprocessed_data//32_col_engy_ess_3d_train_features.nc
-rw-r--r-- 1 S.Rasp ls-craig  24K Jan  8 12:48 /local/S.Rasp/preprocessed_data//32_col_engy_ess_3d_train_norm.nc
-rw-r--r-- 1 S.Rasp ls-craig 5.3G Jan  8 13:03 /local/S.Rasp/preprocessed_data//32_col_engy_ess_3d_train_shuffle_features.nc
-rw-r--r-- 1 S.Rasp ls-craig 3.8G Jan  8 13:03 /local/S.Rasp/preprocessed_data//32_col_engy_ess_3d_train_shuffle_targets.nc
-rw-r--r-- 1 S.Rasp ls-craig 3.8G Jan  8 12:50 /local/S.Rasp/preprocessed_data//32_col_engy_ess_3d_train_targets.nc
-rw-r--r-- 1 S.Rasp ls-craig 5.3G Jan  8 12:59 /local/S.Rasp/preprocessed_data//32_col_engy_ess_3d_valid_features.nc
-rw-r--r-- 1 S.Rasp ls-craig 3.8G Jan  8 12:59 /local/S.Rasp/preprocessed_data//32_col_engy_ess_3d_valid_targets.nc<jupyter_text>Speed comparison: regular hard-drive vs. SSDI copy all the raw data to an SSD drive, in order to speed up the computation. We can load the data from a regular disc (but still save on SSD), to see the difference.<jupyter_code>!python cbrain/preprocess_aqua.py \
--config pp_config/32col_engy_ess_ref.yml \
--in_dir /project/meteo/w2w/A6/S.Rasp/SP-CAM/sp32fbp_andkua/ --aqua_names '*.h1.0000-*-0[1-3]-*' \
--out_pref test_32_col_engy_ess_3d_train<jupyter_output>Time checkpoint reading data: 40.51 s
Number of time steps: 1728
Cut time steps: [ 143  287  431  575  719  863 1007 1151 1295 1439 1583]
Cut time steps: [ 143  287  431  575  719  863 1007 1151 1295 1439 1583]
Time checkpoint create datasets: 40.64 s
Time checkpoint reshape and rechunk: 42.67 s
Compute means and stds
Saving normalization file: /local/S.Rasp/preprocessed_data/test_32_col_engy_ess_3d_train_norm.nc
Time checkpoint normalization arrays: 1352.22 s
Time checkpoint rechunk and ds: 1352.89 s
Save features: /local/S.Rasp/preprocessed_data/test_32_col_engy_ess_3d_train_features.nc
Save targets: /local/S.Rasp/preprocessed_data/test_32_col_engy_ess_3d_train_targets.nc
Total time: 2062.98 s<jupyter_text>Time to process 36 days:SSD: 5 minutes    Hard-drive: 36 minutes  Entire year of data<jupyter_code>!python cbrain/preprocess_aqua.py --config pp_config/32col_engy_ess_ref.yml \
--aqua_names '*.h1.0000-*' --out_pref 32_col_engy_ess_1y_train --ext_norm Nope<jupyter_output>Time checkpoint reading data: 22.92 s
Number of time steps: 17520
Cut time steps: []
Cut time steps: []
Time checkpoint create datasets: 23.35 s
Time checkpoint reshape and rechunk: 27.97 s
Load external normalization file
Time checkpoint normalization arrays: 27.97 s
Time checkpoint rechunk and ds: 33.70 s
Save features: /local/S.Rasp/preprocessed_data/32_col_engy_ess_1y_train_features.nc
Save targets: /local/S.Rasp/preprocessed_data/32_col_engy_ess_1y_train_targets.nc
Total time: 1725.83 s<jupyter_text>So the preprocessing for an entire year takes around 30 minutes.  Shuffling the training datasetDuring network training it's important that the samples in each batch are random. Currently they are still in order (lat-lon-time only have been flattened). Ideally, you would randomly shuffle during training. Because the datasets exceed RAM however, shuffling means reading random locations on disc for each sample in the batch. This is slow.For this reason, I am pre-shuffling the training dataset once. Then during training, for each batch we will read a continuous slice of the shuffled dataset. Only the order of the batches is random. This should be a good enough approximation. Shuffling is done by the `shuffle_ds.py` script in the `cbrain` directory.<jupyter_code>!python cbrain/shuffle_ds.py -h
!python cbrain/shuffle_ds.py --pref /local/S.Rasp/preprocessed_data/32_col_engy_ess_3d_train<jupyter_output>Reading files: /local/S.Rasp/preprocessed_data/32_col_engy_ess_3d_train_features.nc /local/S.Rasp/preprocessed_data/32_col_engy_ess_3d_train_targets.nc
Creating files: /local/S.Rasp/preprocessed_data/32_col_engy_ess_3d_train_shuffle_features.nc /local/S.Rasp/preprocessed_data/32_col_engy_ess_3d_train_shuffle_targets.nc
100%|█████████████████████████████████████████████| 2/2 [00:51<00:00, 25.74s/it]<jupyter_text>Entire year<jupyter_code>!python cbrain/shuffle_ds.py --pref /local/S.Rasp/preprocessed_data/32_col_engy_ess_1y_train<jupyter_output>Reading files: /local/S.Rasp/preprocessed_data/32_col_engy_ess_1y_train_features.nc /local/S.Rasp/preprocessed_data/32_col_engy_ess_1y_train_targets.nc
Creating files: /local/S.Rasp/preprocessed_data/32_col_engy_ess_1y_train_shuffle_features.nc /local/S.Rasp/preprocessed_data/32_col_engy_ess_1y_train_shuffle_targets.nc
100%|███████████████████████████████████████████| 15/15 [14:27<00:00, 57.80s/it]<jupyter_text>Training the neural network and convert to text filesNow we will use the preprocessed files and train the neural network. Then we will save the parameters of the network in text format to use them inside CAM.For now I will only use the command line script `run_experiment.py` which is in the root directory of the repositories. At some point, I will go deeper and analyze what is actually going on during training. In particular, we need to look at the data_generator function, which reads and modifies the preprocessed data.  Training a neural networkAgain, the script takes a config file. The idea of the config files is to make the runs reproducible.<jupyter_code>!cat nn_config/32cols/001_32cols_engy_ess_3d_pnas.yml<jupyter_output>exp_name: 001_32cols_engy_ess_3d_pnas
data_dir: /local/S.Rasp/preprocessed_data/
train_fn: 32_col_engy_ess_3d_train_shuffle
valid_fn: 32_col_engy_ess_3d_valid
norm_fn: 32_col_engy_ess_3d_train_norm.nc
fsub: feature_means
fdiv: max_rs
tmult: target_conv
activation: LeakyReLU
hidden_layers: [256,256,256,256,256,256,256,256,256]
loss: mse
log_dir: ./logs/
epochs: 18
lr_step: 3<jupyter_text>Some explanation:- log_dir: This simply saves a Tensorboard log. I would not worry about this for now since it doesn't really give more information that the scores.- lr_step: In addition to the simple example above, I aplemented a learning rate scheduler. I basically copied this: https://machinelearningmastery.com/using-learning-rate-schedules-deep-learning-models-python-keras/ lr_step defines after how many epochs the learning rate is dropped.- valid_after: If this is true, the validation score is only computed for the last epoch. This is reasonable once you know that your model behaves ok. Computing the validation score actually takes some time.- model_dir and exp_name: The keras model is saved after training in the directory specified with the experiment name. We need the saved model to implement into CAM later.<jupyter_code>!python run_experiment.py -h
!python run_experiment.py -c nn_config/32cols/001_32cols_engy_ess_3d_pnas.yml<jupyter_output>Using TensorFlow backend.
/home/s/S.Rasp/miniconda3/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: compiletime version 3.5 of module 'tensorflow.python.framework.fast_tensor_util' does not match runtime version 3.6
  return f(*args, **kwds)
2019-01-08 13:50:03.328552: I tensorflow/core/platform/cpu_feature_guard.cc:140] Your CPU supports instructions that this TensorFlow binary was not compiled to use: AVX512F FMA
2019-01-08 13:50:03.483059: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1356] Found device 0 with properties: 
name: GeForce GTX 1080 major: 6 minor: 1 memoryClockRate(GHz): 1.7335
pciBusID: 0000:b3:00.0
totalMemory: 7.93GiB freeMemory: 6.47GiB
2019-01-08 13:50:03.483092: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1435] Adding visible gpu devices: 0
2019-01-08 13:50:03.736182: I tensorflow/core/common_runtime/gpu/gpu_device.cc:923] Device interconnect StreamExecutor with strength 1 edge matrix:
2019-01-08 13:50:03.736232: I tensorflow/core/common_ru[...]<jupyter_text>Here we can see that we trained for way too long. How long to train depends on the diversity of the training data and the data amount. Improving this is something to work on, see: https://github.com/raspstephan/CBRAIN-CAM/issues/1We can also see that we do not get any noticeable overfitting in the traditional sense, i.e. difference between training and validation loss.The neural network matrices and architecture are saved in the `./saved_models` directory as h5 files.<jupyter_code>!ls saved_models/001_*<jupyter_output>../../saved_models/001_32cols_engy_ess_3d_pnas.h5<jupyter_text>Entire year<jupyter_code>!cat nn_config/32cols/002_32cols_engy_ess_1y_pnas.yml
!python run_experiment.py -c nn_config/32cols/002_32cols_engy_ess_1y_pnas.yml<jupyter_output>Using TensorFlow backend.
/home/s/S.Rasp/miniconda3/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: compiletime version 3.5 of module 'tensorflow.python.framework.fast_tensor_util' does not match runtime version 3.6
  return f(*args, **kwds)
2019-01-10 17:14:27.929078: I tensorflow/core/platform/cpu_feature_guard.cc:140] Your CPU supports instructions that this TensorFlow binary was not compiled to use: AVX512F FMA
2019-01-10 17:14:28.176396: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1356] Found device 0 with properties: 
name: GeForce GTX 1080 major: 6 minor: 1 memoryClockRate(GHz): 1.7335
pciBusID: 0000:b3:00.0
totalMemory: 7.93GiB freeMemory: 6.85GiB
2019-01-10 17:14:28.176438: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1435] Adding visible gpu devices: 0
2019-01-10 17:14:29.185087: I tensorflow/core/common_runtime/gpu/gpu_device.cc:923] Device interconnect StreamExecutor with strength 1 edge matrix:
2019-01-10 17:14:29.185123: I tensorflow/core/common_ru[...]<jupyter_text>During the training above, my GPU usage is around 30%.  Checking the offline performance of the neural networkSometimes it can be handy to check out how well the neural network actually does when run offline for the validation set. For this I wrote a `ModelDiagnostics` class. This probably should be revised at some point, but for now I will go through some basic functionality.One useful thing you can do with this class is to test the neural network with a validation set from a different reference climate and see how well it does. This gives you a first sense of generalization.If someone has questions about this, try understanding the source code which is somewhat documented or ask me!<jupyter_code>from cbrain.model_diagnostics import *
from cbrain.utils import *
# Do this if you are using the GPU!
limit_mem()
# LMU
MODEL_DIR = '/home/s/S.Rasp/repositories/CBRAIN-CAM/saved_models/'
# Create a ModelDiagnostics object
norms = ('feature_means', 'max_rs', None, 'target_conv')

diag = ModelDiagnostics(
    MODEL_DIR + '001_32cols_engy_ess_3d_pnas.h5',
    fpath=PREPROC_DIR + '32_col_engy_ess_3d_valid_features.nc',
    tpath=PREPROC_DIR + '32_col_engy_ess_3d_valid_targets.nc',
    npath=PREPROC_DIR + '32_col_engy_ess_3d_train_norm.nc',
    norms=norms,
)<jupyter_output><empty_output><jupyter_text>With this object, you can essentially do two things. First, make plots that compare the neural network predictions for one time step to the truth. And second, compute some average statistics. Let's do the latter one first.<jupyter_code># First do this to compute the statistics over the entire dataset
diag.compute_stats()
# The print the mean_statistics ignoring the upper 10 levels 
# (because the tendencies are small, some metrics are nonsense there)
diag.mean_stats(9)<jupyter_output><empty_output><jupyter_text>A lot is printed here. Things you might be interested in are the r2 (cf. Gentine et al 2018). Check out the hor_r2 for a more useful measure. Here we can see, e.g. that moisture isn't so good but the radiative fluxes are really good.Now let's do some plotting.<jupyter_code>diag.plot_double_xy(itime=0, ilev=0, var='PRECT');
diag.plot_double_yz(itime=0, ilon=0, var='TPHYSTND', vmin=-1e-4, vmax=1e-4, cmap='RdBu_r');<jupyter_output><empty_output><jupyter_text>In these figures you can see the lack of variability in the CBRAIN output.  Converting .h5 file to txt fileThe next step here is now to convert the h5 file to files that can be read by the neural network version of CAM. For this we will use another script, called `save_weights.py`, which takes the same config file as `run_experiment.py`.<jupyter_code>!python save_weights.py -c nn_config/32cols/001_32cols_engy_ess_3d_pnas.yml
!ls saved_models/001_32cols_engy_ess_3d_pnas/<jupyter_output>inp_max_rs.txt	     layer1_bias.txt	layer5_bias.txt    layer9_bias.txt
inp_maxs.txt	     layer1_kernel.txt	layer5_kernel.txt  layer9_kernel.txt
inp_means.txt	     layer2_bias.txt	layer6_bias.txt    outp_maxs.txt
inp_mins.txt	     layer2_kernel.txt	layer6_kernel.txt  outp_mins.txt
inp_stds_by_var.txt  layer3_bias.txt	layer7_bias.txt    weights.h5
inp_stds.txt	     layer3_kernel.txt	layer7_kernel.txt
layer10_bias.txt     layer4_bias.txt	layer8_bias.txt
layer10_kernel.txt   layer4_kernel.txt	layer8_kernel.txt<jupyter_text>The weights and biases for all the layers are saved in separate text files, plus the relevant normalization arrays.<jupyter_code>!cat saved_models/001_32cols_engy_ess_3d_pnas/layer1_kernel.txt | head -1<jupyter_output>-2.247076e-01,2.444012e-02,-6.513681e-02,1.084918e-01,-1.971805e-02,1.679797e-02,6.645612e-02,1.205308e-01,9.274448e-02,-6.410713e-02,-5.922029e-02,-8.371486e-02,-6.580098e-02,-7.378058e-02,-1.398824e-02,-1.351725e-01,-1.762785e-01,7.083420e-02,2.030234e-03,2.106156e-01,3.582359e-01,3.041640e-01,2.500775e-01,-5.012782e-02,-1.125699e-01,4.082102e-02,-1.728023e-01,-6.359089e-03,-1.227781e-01,-1.985432e-02,1.931775e-01,1.839796e-02,-7.847431e-02,-6.119738e-02,-1.337386e-01,1.668242e-01,3.829718e-02,5.141954e-02,5.827349e-01,1.015465e+00,8.116463e-01,6.864675e-01,5.228361e-01,1.373079e-01,-8.082748e-02,-4.842604e-02,8.416072e-03,7.908968e-02,-4.162578e-02,5.653955e-02,2.271484e-02,9.763756e-02,4.991648e-02,-3.098556e-02,-6.799575e-03,1.183797e-01,-1.354993e-01,9.173439e-02,2.124729e-02,-4.702663e-02,-6.526815e-02,1.570719e-01,5.800200e-02,-9.230268e-02,-6.970107e-02,-3.576361e-02,8.898021e-03,7.295761e-02,-7.373397e-02,-4.990729e-02,7.852868e-02,-1.300494e-02,-1.297906e-01,8.923770e-02,5.7[...]<jupyter_text>Switch to StampedeNow I am on Stampede2. Here we will compile NNCAM, run it and analyze the output. At this point, my workflow becomes very manual and non-reproducible. Improvement is desparately needed!The NNCAM code is here: https://gitlab.com/mspritch/spcam3.0-neural-net<jupyter_code>%pwd
%cd
!ls<jupyter_output>bin			       miniconda3    saved_weights
cloudbrain_ctrl_aquaplanet_03  netcdf-3.6.3  share
include			       nn_weights    spcam3-nn-build
intel			       repositories  spcam3-precip_output-build
lib			       saved_models  tmp<jupyter_text>Copy parameters and compile CAM  Copy files to StampedeFirst we need to copy the text files containing the weights and parameters to Stampede. For this I have the following command in my `.bashrc`:```weights_from_lmu(){   rsync -av <EMAIL>:/home/s/S.Rasp/repositories/CBRAIN-CAM/saved_models/"$1"* /home1/05488/tg847872/saved_models/}```Now I can simply type:```weights_from_lmu 001_32cols_engy_ess_3d_pnas```to copy the files to my local `saved_models` directory in my home directory.<jupyter_code>!ls saved_models/001*<jupyter_output>saved_models/001_32cols_engy_ess_3d_pnas.h5

saved_models/001_32cols_engy_ess_3d_pnas:
inp_max_rs.txt	     layer1_bias.txt	layer5_bias.txt    layer9_bias.txt
inp_maxs.txt	     layer1_kernel.txt	layer5_kernel.txt  layer9_kernel.txt
inp_means.txt	     layer2_bias.txt	layer6_bias.txt    outp_maxs.txt
inp_mins.txt	     layer2_kernel.txt	layer6_kernel.txt  outp_mins.txt
inp_stds_by_var.txt  layer3_bias.txt	layer7_bias.txt    weights.h5
inp_stds.txt	     layer3_kernel.txt	layer7_kernel.txt
layer10_bias.txt     layer4_bias.txt	layer8_bias.txt
layer10_kernel.txt   layer4_kernel.txt	layer8_kernel.txt<jupyter_text>Compiling CAMHere is where it gets ugly. The CAM repository is not very well structures and has a lot of branches for different options that should instead be namelist switches.The "correct" NNCAM branch is `nn_fbp_engy_ess`.<jupyter_code>%cd repositories/spcam3.0-neural-net/
!git status 
!git log -1<jupyter_output>[33mcommit af6d10bc2566982dec7344f49080cb0d7c0c4582[m
Author: raspstephan <<EMAIL>>
Date:   Tue Jul 24 11:03:25 2018 -0500

    Not entirely sure what is happening here, new modes for everything?<jupyter_text>The neural network is hardcoded in `spcam3.0-neural-net/models/atm/cam/src/physics/cam1/cloudbrain_keras_dense.F90`.This this means that that every time the network architecture is changed, the parameters in this file also have to be changed accordingly and the model has to be recompiled. All these options should really be namelist switches, see this issue: https://gitlab.com/mspritch/spcam3.0-neural-net/issues/1To compile cam this script has to be called `spcam3.0-neural-net/scripts/build-spsld-stampede2-intel.csh`The compiled file will be placed in `spcam3-nn-build/run/`.Again, I have an alias: `alias build='/home1/05488/tg847872/repositories/spcam3.0-neural-net/scripts/build-spsld-stampede2-intel.csh && ll /home1/05488/tg847872/spcam3-nn-build/run/`To compile I simply type:```build```<jupyter_code>!ls -l /home1/05488/tg847872/spcam3-nn-build/run/<jupyter_output>total 18776
-rwx------ 1 tg847872 G-818279 19220288 Jan  9 06:06 cam<jupyter_text>Running NNCAMNow we have a compiles version of CAM and the weights and biases in a format CAM can read.To run cam we first go to the scratch directory. On Stampede you can simply use `cds`.<jupyter_code>%cd /scratch/05488/tg847872<jupyter_output>/scratch/05488/tg847872<jupyter_text>Then we create a new directory for the run.<jupyter_code>%mkdir nncam_001_32cols_engy_ess_3d_pnas
%cd nncam_001_32cols_engy_ess_3d_pnas/<jupyter_output>/scratch/05488/tg847872/nncam_001_32cols_engy_ess_3d_pnas<jupyter_text>Here we need to do a number of things.First, we need to copy the cam binary. Again I have an alias: `alias copy='cp /home1/05488/tg847872/spcam3-nn-build/run/cam ./'`<jupyter_code>!cp /home1/05488/tg847872/spcam3-nn-build/run/cam ./<jupyter_output><empty_output><jupyter_text>Second, we need to copy the neural network weights and biases. These must be saved in a directory called `keras_matrices`.<jupyter_code>!cp -r /home1/05488/tg847872/saved_models/001_32cols_engy_ess_3d_pnas ./keras_matrices<jupyter_output><empty_output><jupyter_text>Third, we need to copy the initial conditions from `/work/05488/tg847872/stampede2/ic_bc/`.<jupyter_code>!cp /work/05488/tg847872/stampede2/ic_bc/spinup_AndKua_aqua_SPCAM3.0.cam2.i.0000-12-02-00000.nc ./<jupyter_output><empty_output><jupyter_text>Finally, we need to copy and modify the namelist `atm_in` and the batch script `run.slurm`.In `atm_in`, change the `caseid`. Note that it cannot exceed 32 characters.In `run.slurm`, change the jobname, potentially the queue (`normal` or `flat-quadrant`), and most importantly the directory.i<jupyter_code>!cat atm_in
!cat run.slurm<jupyter_output>#!/bin/bash
#SBATCH -J 001
#SBATCH -o myjob.o%j       
#SBATCH -e myjob.e%j       
#SBATCH -p normal          
#SBATCH -N 4              
#SBATCH -n 64              
#SBATCH -t 2-00:00:00        
#SBATCH --mail-user=<EMAIL>
#SBATCH --mail-type=all    
#SBATCH -A TG-ATM170029       
export I_MPI_COMPATIBILITY=4
export OMP_NUM_THREADS=16
cd /scratch/05488/tg847872/nncam_001_32cols_engy_ess_3d_pnas
ibrun -n 64 ./cam < atm_in<jupyter_text>Interactive testingIn order to check whether everything runs, its better to start a quick interactive job than to wait for the job to start in the queue. The maximum interactive time is 120 minutes.For this type: `idev -N 4 -n 64 -m 120 -A TG-ATM170029`Or with an alias: `alias interactive='idev -N 4 -n 64 -m 120 -A TG-ATM170029'`Then to run cam:```bash run.slurm```There will be a lot of output, so it might be best to write the output in a text file ` >> out.txt &`.If the model crashes, it is likely it will crash within two hours, so there is no need to queue a job unless it survives the two hour testing.<jupyter_code>!ls<jupyter_output>atm_in
cam
cam2.nn_001_32cols_engy_ess_3d_pnas.rpointer
keras_matrices
nn_001_32cols_engy_ess_3d_pnas.cam2.h0.0000-01.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h0.0000-02.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h0.0000-03.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h1.0000-01-01-00000.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h1.0000-01-11-00000.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h1.0000-01-21-00000.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h1.0000-01-31-00000.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h1.0000-02-10-00000.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h1.0000-02-20-00000.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h1.0000-03-02-00000.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h1.0000-03-12-00000.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h1.0000-03-22-00000.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h1.0000-04-01-00000.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h1.0000-04-11-00000.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h2.0000-01-01-00000.nc
nn_001_32cols_engy_ess_3d_pnas.cam2.h2.0000-01-11-00[...]<jupyter_text>Analyze model outputNow, we would like to actually check out the model output. **IMPORTANT**: The scratch directory on Stampede is not safe. If you have a run that you want to save, it must be copied somewhere.For a first sanity check, we will load the `h0` output files. These are monthly aggregate statistics of the most important quantities. As a comparison, we will use the reference training run.<jupyter_code>import xarray as xr
import matplotlib.pyplot as plt
%matplotlib inline
REF_DIR = '/work/05488/tg847872/stampede2/sp32fbp_andkua/'
h0s = xr.open_mfdataset('*h0*', decode_times=False, concat_dim='time')
# Unfortunately, this old run does not yet have TE (total energy) and TW (total water) in the h0 output.
# I only added this to later model version.
h0s_ref = xr.open_mfdataset(REF_DIR + '*h0*', decode_times=False, concat_dim='time')
def gw_avg(ds, var=None, da=None):
    """Computes the area weighted mean over lat-lon"""
    da = ds[var] if da is None else da
    return (da * ds['gw']).sum('lat').mean('lon') / 2.
# Mean surface temperature drift
gw_avg(h0s.isel(lev=29), 'TAP').plot(label='NNCAM')
gw_avg(h0s_ref.isel(lev=29), 'TAP').plot(label='SPCAM')
plt.legend();
# Temperature difference for month 12
(h0s.TAP.isel(time=11).mean(dim='lon') - h0s_ref.TAP.isel(time=11).
 mean(dim='lon')).plot(yincrease=False, vmin=-5, vmax=5, cmap='RdBu_r')
plt.title('NNCAM - SPCAM');
# Humidity difference for month 12
(h0s.QAP.isel(time=11).mean(dim='lon') - h0s_ref.QAP.isel(time=11).
 mean(dim='lon')).plot(yincrease=False, cmap='RdBu_r')
plt.title('NNCAM - SPCAM');
# U wind difference for month 12
(h0s.UAP.isel(time=11).mean(dim='lon') - h0s_ref.UAP.isel(time=11).
 mean(dim='lon')).plot(yincrease=False, cmap='RdBu_r')
plt.title('NNCAM - SPCAM');
# Physics tendency teperature difference for month 12
(h0s.TPHYSTND.isel(time=11).mean(dim='lon') - h0s_ref.TPHYSTND.isel(time=11).
 mean(dim='lon')).plot(yincrease=False, cmap='RdBu_r')
plt.title('NNCAM - SPCAM');
# Physics tendency humidity difference for month 12
(h0s.PHQ.isel(time=11).mean(dim='lon') - h0s_ref.PHQ.isel(time=11).
 mean(dim='lon')).plot(yincrease=False, cmap='RdBu_r')
plt.title('NNCAM - SPCAM');<jupyter_output><empty_output><jupyter_text>So we can see that for this particular NNCAM version we have quite large deviations with a drift in NNCAM that seems to saturate at some point. The model seems to be stable, however. This is interesting since I remember a previous attempt using 32 columns crashed.We can also look at the time-step level output which is saved in the `h1` files.<jupyter_code>%%time
# Reading can take a while. Let's only read the files for one month.
h1s = xr.open_mfdataset('*h1.0000-11*', decode_times=False, concat_dim='time')<jupyter_output>CPU times: user 108 ms, sys: 10.7 ms, total: 119 ms
Wall time: 1.57 s<jupyter_text>We can now check whether the immediate outputs by the neural network are reasonable in the variables `TPHYSTND = NNDT` and `PHQ = NNDQ`.<jupyter_code># Plot a horizontal slice at level 10
h1s.TPHYSTND.isel(time=10, lev=10).plot();
# Plot a vertical slice
h1s.TPHYSTND.isel(time=10, lon=3).plot(yincrease=False);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09666009509451468, 'fraction_numerical': 0.07903281920445321, 'mean_word_length': 5.201546206400575, 'pattern_counts': {'":': 0, '<': 129, '<?xml version=': 0, '>': 130, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '36', 'dir': 'jupyter-structured-clean-dedup', 'id': '141918', 'license': 'MIT', 'n_tokens_mistral': 14251, 'n_tokens_neox': 12375, 'n_words': 3399, 'path': 'notebooks/stephans-devlog/1.0-Entire-workflow-for-32-column-run.ipynb', 'repo_name': 'ankitesh97/CBRAIN-CAM', 'size': '34546'}