initial commit of turing pattern generator

prepare_data.py

@@ -0,0 +1,9 @@
+#!/usr/bin/env python
+
+import reacdiff.parsing as parsing
+import reacdiff.data.preprocess as preprocess
+
+
+if __name__ == '__main__':
+    args = parsing.parse_dataprep_args()
+    preprocess.preprocess(args.data_path)

reacdiff/data/data.py

@@ -1,3 +1,4 @@
+import h5py
 import numpy as np
 
 import reacdiff.utils as utils
@@ -49,16 +50,21 @@ def __getitem__(self, item):
             return Dataset(self.data[item], self.targets[item], self.data2[item])
 
 
-def load_data(path, targets=False):
-    if targets:
-        return np.random.rand(100, 26)
-    data = np.random.rand(100, 50, 128, 128)
-    data = np.expand_dims(data, axis=-1)
+def load_data(path, max_read=None, key='data'):
+    """
+    :param path: An HDF5 file.
+    :param max_read: Maximum number of samples to read.
+    :param key: Data key in HDF5 file.
+    :return: Numpy array of data
+    """
+    with h5py.File(path, 'r') as f:
+        data = f[key][:] if max_read is None else f[key][:max_read]
     return data
 
 
-def save_data(data, path):
-    raise NotImplementedError
+def save_data(data, path, key='data'):
+    with h5py.File(path, 'w') as f:
+        f.create_dataset(key, data=data)
 
 
 def split_data(data, splits=(0.9, 0.05, 0.05), seed=None):

reacdiff/data/preprocess.py

@@ -0,0 +1,34 @@
+import os
+
+import h5py
+import numpy as np
+
+import reacdiff.data.data as datamod
+
+
+def preprocess(path):
+    """
+    :param path: An HDF5 file.
+    """
+    f = h5py.File(path, 'r')
+    targets = f['A1']
+    states = f['y']
+
+    assert len(targets.shape) == 2
+    assert len(states.shape) == 5
+    assert states.shape[0] == 2
+
+    targets = targets[:]
+    states = states[:]
+    targets = np.swapaxes(targets, 0, 1)
+    states = np.transpose(states, (0, 4, 1, 2, 3))
+    data, data2 = np.expand_dims(states, axis=-1)
+
+    name = os.path.splitext(path)[0]
+    targets_path = name + '_targets.h5'
+    data_path = name + '_states.h5'
+    data_path2 = name + '_states2.h5'
+
+    datamod.save_data(targets, targets_path)
+    datamod.save_data(data, data_path)
+    datamod.save_data(data2, data_path2)

reacdiff/parsing.py

@@ -0,0 +1,89 @@
+import argparse
+import os
+
+
+def parse_dataprep_args():
+    parser = argparse.ArgumentParser(
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter
+    )
+    parser.add_argument('--data_path', type=str, required=True,
+                        help='Path to HDF5 data file')
+    return parser.parse_args()
+
+
+def parse_train_args():
+    parser = argparse.ArgumentParser(
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter
+    )
+
+    # General arguments
+    parser.add_argument('--data_path', type=str, required=True,
+                        help='Path to data containing observable states')
+    parser.add_argument('--targets_path', type=str, required=True,
+                        help='Path to targets')
+    parser.add_argument('--data_path2', type=str,
+                        help='Path to data containing additional observable states')
+    parser.add_argument('--save_dir', type=str, default=os.getcwd(),
+                        help='Directory where model checkpoints will be saved')
+    parser.add_argument('--splits', type=float, nargs=3, default=[0.9, 0.05, 0.05],
+                        help='Split proportions for train/validation/test sets')
+    parser.add_argument('--save_test_data', action='store_true',
+                        help='Save test data and targets to file')
+    parser.add_argument('--test_data_path', type=str,
+                        help='Path to separate test data')
+    parser.add_argument('--test_targets_path', type=str,
+                        help='Path to separate test targets')
+    parser.add_argument('--test_data_path2', type=str,
+                        help='Path to separate test data for additional observable')
+    parser.add_argument('--quiet', action='store_true',
+                        help='Do not print model details and batch information during training')
+    parser.add_argument('--cpu', action='store_true',
+                        help='Run on CPU instead of GPU')
+    parser.add_argument('--seed', type=int, default=7,
+                        help='Random seed for reproducibility')
+
+    # Training arguments
+    parser.add_argument('--batch_size', type=int, default=32,
+                        help='Batch size')
+    parser.add_argument('--epochs', type=int, default=30,
+                        help='Maximum number of epochs to run')
+    parser.add_argument('--lr_start', type=float, default=1e-3,
+                        help='Initial learning rate')
+    parser.add_argument('--lr_end', type=float, default=1e-4,
+                        help='Final learning rate')
+    parser.add_argument('--patience', type=int, default=5,
+                        help='Number of epochs without improvement after which training will be stopped')
+    parser.add_argument('--max_norm', type=float, default=2.0,
+                        help='Maximum gradient norm before clipping occurs')
+
+    # Encoder arguments
+    parser.add_argument('--feat_maps', type=int, default=16,
+                        help='Number of feature maps in first convolutional layer')
+    parser.add_argument('--first_conv_size', type=int, default=7,
+                        help='Filter size in first convolutional layer')
+    parser.add_argument('--first_conv_stride', type=int, default=2,
+                        help='Strides in first convolutional layer')
+    parser.add_argument('--first_pool_size', type=int, default=3,
+                        help='Window size in first pool')
+    parser.add_argument('--first_pool_stride', type=int, default=2,
+                        help='Strides in first pool')
+    parser.add_argument('--growth_rate', type=int, default=12,
+                        help='Growth rate in dense blocks')
+    parser.add_argument('--blocks', type=int, nargs='+', default=[3, 4, 5],
+                        help='Numbers of layers in each dense block')
+    parser.add_argument('--dropout', type=float, default=0.0,
+                        help='Dropout rate after convolutions')
+    parser.add_argument('--reduction', type=float, default=0.5,
+                        help='Compression rate in transition blocks')
+    parser.add_argument('--no_bottleneck', action='store_true',
+                        help='Do not use bottleneck convolution in dense blocks')
+    parser.add_argument('--flatten_last', action='store_true',
+                        help='Flatten instead of global pool before output')
+
+    # RNN arguments
+    parser.add_argument('--rnn_layers', type=int, default=1,
+                        help='Number of RNN layers')
+    parser.add_argument('--rnn_units', type=int, default=100,
+                        help='Number of units in RNN')
+
+    return parser.parse_args()

reacdiff/train/run_training.py

@@ -0,0 +1,77 @@
+import os
+
+import reacdiff.data.data as datamod
+import reacdiff.models.crnn as models
+import reacdiff.train.train as train
+
+
+def run_training(args):
+    # Load data
+    print('Loading data')
+    data = datamod.Dataset(
+        datamod.load_data(args.data_path),
+        datamod.load_data(args.targets_path),
+        None if args.data_path2 is None else datamod.load_data(args.data_path2)
+    )
+
+    # Split data
+    print(f'Splitting data with seed {args.seed}')
+    train_data, val_data, test_data = datamod.split_data(data, splits=args.splits, seed=args.seed)
+
+    if args.test_data_path is not None and args.test_targets_path is not None:
+        test_data = datamod.Dataset(
+            datamod.load_data(args.test_data_path),
+            datamod.load_data(args.test_targets_path),
+            None if args.test_data_path2 is None else datamod.load_data(args.test_data_path2)
+        )
+
+    os.makedirs(args.save_dir, exist_ok=True)
+
+    if args.save_test_data:
+        test_data.save(os.path.join(args.save_dir, 'test'))
+
+    print(f'Train size = {len(train_data):,} | val size = {len(val_data):,} | test size = {len(test_data):,}')
+
+    # Build model
+    crnn = models.CRNN(
+        time_steps=train_data.data.shape[1],
+        input_shape=train_data.data.shape[2:],
+        output_dim=train_data.targets.shape[1],
+        observables=train_data.get_num_observables(),
+        rnn_layers=args.rnn_layers,
+        rnn_units=args.rnn_units,
+        use_gpu=not args.cpu
+    )
+    crnn.build(
+        feat_maps=args.feat_maps,
+        first_conv_size=args.first_conv_size,
+        first_conv_stride=args.first_conv_stride,
+        first_pool_size=args.first_pool_size,
+        first_pool_stride=args.first_pool_stride,
+        growth_rate=args.growth_rate,
+        blocks=args.blocks,
+        dropout=args.dropout,
+        reduction=args.reduction,
+        bottleneck=not args.no_bottleneck,
+        flatten_last=args.flatten_last
+    )
+
+    model_dir = os.path.join(args.save_dir, 'model')
+    os.makedirs(model_dir, exist_ok=True)
+
+    # Train model
+    train.train(crnn, train_data, model_dir,
+                batch_size=args.batch_size,
+                epochs=args.epochs,
+                val_data=val_data,
+                patience=args.patience,
+                lr_start=args.lr_start,
+                lr_end=args.lr_end,
+                max_norm=args.max_norm,
+                quiet=args.quiet)
+
+    # Evaluate on test set
+    print('Evaluating test data')
+    metrics = crnn.model.evaluate(test_data.get_data(), test_data.targets,
+                                  batch_size=args.batch_size)
+    print(*(f'{name}: {metric:.4f}' for name, metric in zip(crnn.model.metrics_names, metrics)), sep='; ')

reacdiff/train/train.py

@@ -0,0 +1,59 @@
+import os
+
+import keras
+
+import reacdiff.utils as utils
+
+
+def train(crnn, data, save_dir,
+          batch_size=32,
+          epochs=30,
+          val_data=None,
+          patience=5,
+          lr_start=1e-3,
+          lr_end=1e-4,
+          max_norm=2.0,
+          quiet=False):
+    """
+    :param crnn: CRNN instance.
+    :param data: Training data.
+    :param save_dir: directory to save models to.
+    :param batch_size: int, batch size.
+    :param epochs: maximum number of epochs.
+    :param val_data: Validation data.
+    :param patience: patience.
+    :param lr_start: float, initial learning rate.
+    :param lr_end: float, final learning rate.
+    :param max_norm: float, maximum gradient norm.
+    :param quiet: bool, print less information
+    """
+    optimizer = keras.optimizers.Adam(lr=lr_start, clipnorm=max_norm)
+    crnn.model.compile(optimizer=optimizer, loss='mse', metrics=[utils.rmse, utils.mae])
+    if not quiet:
+        crnn.encoder.summary()
+        crnn.rnn.summary()
+        crnn.model.summary()
+
+    model_name = 'model.{epoch:03d}.h5'
+    model_path = os.path.join(save_dir, model_name)
+
+    def lr_schedule(epoch, _):
+        return lr_start * ((lr_end / lr_start) ** (epoch / epochs))
+
+    checkpoint = keras.callbacks.ModelCheckpoint(filepath=model_path,
+                                                 verbose=1,
+                                                 save_best_only=True)
+    early_stopping = keras.callbacks.EarlyStopping(patience=patience,
+                                                   verbose=1,
+                                                   restore_best_weights=True)
+    lr_scheduler = keras.callbacks.LearningRateScheduler(lr_schedule,
+                                                         verbose=1)
+    callbacks = [checkpoint, early_stopping, lr_scheduler]
+
+    crnn.model.fit(data.get_data(), data.targets,
+                   batch_size=batch_size,
+                   epochs=epochs,
+                   validation_data=(val_data.get_data(), val_data.targets),
+                   shuffle=True,
+                   verbose=2 if quiet else 1,
+                   callbacks=callbacks)

solver/DDFT.m

@@ -0,0 +1,51 @@
+%parameter search
+L = [5,5]*2;
+N = [256,256]*2;
+n = prod(N);
+n0 = 0.07;
+sigma = 0.01;
+
+rng(1);
+n0 = n0 + sigma*randn(N);
+y0 = fftn(n0);
+y0 = y0(:);
+% y00 = y0(1);
+% y0(1) = [];
+
+% %use image?
+% im = imread('cameraman.tif');
+% im = imresize(im,N);
+% C2 = fft(im);
+% C2 = C2(:);
+% C2(1) = 0;
+% C2 = C2/max(abs(C2))/1.5;
+
+[k2,k] = formk(N,L);
+k2 = k2(:);
+% k2(1) = [];
+
+kind1 = k{1}; kind2 = k{2};
+[kind1,kind2] = ndgrid(kind1(:),kind2(:));
+theta = angle(kind1+i*kind2);
+
+theta = theta(:);
+
+k0 = 10;
+alpha = 2;
+C2 = exp(-(sqrt(k2)-k0).^2/(2*alpha^2));
+C22 = C2;
+% C22 = C2 .* (1+cos(theta))/2;
+% C22 = C2 .* (exp(-(theta-pi/2).^2/(2*(50/180*pi)^2))+exp(-(theta+pi/2).^2/(2*(50/180*pi)^2)))/2;
+J = -k2.*(-C22);
+
+tspan = linspace(0,0.14,100);
+tspan = [0,1.7];
+tspan = linspace(0,1.07,100);
+% tspan = [0,0.14];
+options = odeset('AbsTol',1e-3);
+[tout,yout] = odeimex(@(t,y) DDFT_nlin(t,y,k2,N),J,tspan,y0);
+% Nt = 10000;
+% h = (tspan(end)-tspan(1))/Nt;
+% [tout,yout] = odeimexez(@(t,y) DDFT_nlin(t,y,k2,N),J,h,Nt,y0,[],1:100:Nt);
+
+figure; k2real(yout(:,end),N);

solver/DDFT_nlin.m

@@ -0,0 +1,17 @@
+function dy = DDFT_nlin(t,y,k2,N,y00)
+  if nargin > 4 && ~isempty(y00)
+    %k=0 component is omitted
+    k0flag = true;
+    y = [y00;y(:)];
+  else
+    k0flag = false;
+  end
+  y = real(ifftn(reshape(y,N)));
+  % mu = y - y.^2/2 + y.^3/3;
+  mu = log(1+y);
+  fmu = reshape(fftn(mu,N),[],1);
+  if k0flag
+    fmu(1) = [];
+  end
+  dy = -k2 .* fmu;
+end