Revive old code.

.dockerignore

@@ -0,0 +1,3 @@
+data/
+img/
+**/__pycache__/

.gitignore

@@ -0,0 +1,4 @@
+data/
+img/
+**/__pycache__/
+*.swp

Dockerfile

@@ -0,0 +1,12 @@
+FROM gw000/keras:2.0.8-py3-tf-cpu
+
+RUN pip3 install scikit-image sklearn
+
+RUN mkdir /data
+ADD icenet/ /tmp/icenet
+ADD run_model.py /tmp/run_model.py
+WORKDIR /tmp
+
+ENV PYTHONUNBUFFERED=1
+
+ENTRYPOINT ["python3", "run_model.py"]

README.md

@@ -0,0 +1,62 @@
+# Ship vs iceberg discriminator
+
+TL;DR: Discriminate between ships and icebergs from SAR imagery.
+
+
+## Approach
+
+Data augmentation and parameter sharing. CNN and ResNets.
+
+
+## Data directory
+
+The data directory is expected to have the structure as shown below:
+
+    data/
+    ├── params
+    │   ├── base_cnn-scaling.pkl
+    │   ├── base_cnn-weights-loss.h5
+    │   ├── base_cnn-weights-val_loss.h5
+    │   ├── icenet-weights-loss.h5
+    │   └── icenet-weights-val_loss.h5
+    ├── predictions
+    │   └── icenet-dev.csv
+    ├── sample_submission.csv
+    ├── test.json
+    └── train.json
+
+where `{train,test}.json` is the data from the
+[kaggle website](https://www.kaggle.com/c/statoil-iceberg-classifier-challenge).
+
+
+## Log
+
+### Residual base CNN
+
+Summary:
+
+  * Test loss: 0.5099
+  * Test accuracy: 0.7932
+  * Epochs: 100
+  * Best val loss at epoch 70 (converged until 100, did not overfit)
+
+Comments:
+
+  * Low variance -- training loss is consistently a bit lower than validation
+    loss.
+  * Since images are "artificially" labeled, it is hard to say what the bias is.
+    There should be some bias since this network does not overfit, and it also
+    looks like training converges after 100 epochs (with decaying learning rate).
+  * There may also be labeling noise. It is indeed suspicious that the validation
+    loss converges with very low variance. Perhaps revisit the labeling
+    approach for the base generator.
+  * Conclusion: Check labeling, then bring out the big guns and expand the
+    residual net.
+
+With this model as a basis for the 9 regions, followed by a reshape, conv and
+two dense layers, yields ok performance: Around 0.20 loss after few epochs.
+
+However, validation loss is often lower than training loss. It might be that
+the two distributions are not the same for both networks -- check the random
+seed and verify! Might also be noisy training (because of augmentation).
+

VERSION

@@ -0,0 +1 @@
+0.2

icenet/__init__.py

@@ -0,0 +1,12 @@
+
+try:
+    from icenet import base_cnn, icenet
+except Exception as ex:
+    print("Failed to import keras")
+    print(ex)
+    base_cnn = icenet = None
+
+MODELS = {
+    'base_cnn': base_cnn,
+    'icenet': icenet 
+}

icenet/base_cnn.py

@@ -0,0 +1,280 @@
+"""Base CNN."""
+
+import random
+import numpy as np
+from itertools import islice
+from icenet import util
+from icenet import residual
+from icenet import resnet
+
+from keras.models import Model
+from keras.initializers import glorot_uniform
+from keras.optimizers import Adam, SGD
+from keras.callbacks import ModelCheckpoint, EarlyStopping, ReduceLROnPlateau
+from keras.models import load_model
+from keras.layers import (
+        Input, Add, Dense, Activation, Conv2D, BatchNormalization, Flatten,
+        MaxPooling2D, AveragePooling2D, Dropout, Concatenate, ZeroPadding2D
+)
+
+
+MODEL_NAME = 'base_cnn'
+
+
+def get_model_simple(img_shape):
+    X_img = Input(img_shape)
+    X = ZeroPadding2D((2, 2))(X_img)
+
+    # Conv 1 
+    X = Conv2D(64,
+               kernel_size=(5, 5),
+               strides=(1, 1),
+               padding='valid',
+               activation='relu',
+               kernel_initializer=glorot_uniform(seed=0),
+               name='conv1')(X)
+    X = Conv2D(128,
+               kernel_size=(3, 3),
+               strides=(1, 1),
+               padding='valid',
+               activation='relu',
+               kernel_initializer=glorot_uniform(seed=0),
+               name='conv2')(X)
+    X = MaxPooling2D((3, 3), strides=(2, 2))(X)
+    X = Dropout(0.2)(X)
+
+    X = Conv2D(128,
+               kernel_size=(3, 3),
+               strides=(1, 1),
+               padding='valid',
+               activation='relu',
+               kernel_initializer=glorot_uniform(seed=0),
+               name='conv3')(X)
+    X = MaxPooling2D((2, 2), strides=(2, 2))(X)
+    X = Dropout(0.2)(X)
+
+    X = Conv2D(64,
+               kernel_size=(3, 3),
+               strides=(1, 1),
+               padding='valid',
+               activation='relu',
+               kernel_initializer=glorot_uniform(seed=0),
+               name='conv4')(X)
+    X = MaxPooling2D((3, 3), strides=(1, 1))(X)
+    X = Dropout(0.2)(X)
+
+    X = Flatten()(X)
+    X = Dense(512,
+              activation='relu',
+              kernel_initializer=glorot_uniform(seed=0),
+              name='fc1')(X)
+    X = Dropout(0.2)(X)
+
+    X = Dense(256,
+              activation='relu',
+              kernel_initializer=glorot_uniform(seed=0),
+              name='fc2')(X)
+    X = Dropout(0.2)(X)
+
+    X = Dense(3,
+              activation='sigmoid',
+              kernel_initializer=glorot_uniform(seed=0),
+              name='y_hat')(X)
+
+    return Model(inputs=X_img, outputs=X, name=MODEL_NAME)
+
+
+def get_model_residual(img_shape):
+    X_img = Input(img_shape)
+
+    X = ZeroPadding2D((2, 2))(X_img)
+
+    # Conv 1 
+    X = Conv2D(32,
+               kernel_size=(5, 5),
+               strides=(1, 1),
+               padding='valid',
+               kernel_initializer=glorot_uniform(seed=0),
+               name='conv1')(X)
+    X = BatchNormalization(axis=3, name='bn_conv1')(X)
+    X = Activation('relu')(X)
+    X = MaxPooling2D((2, 2), strides=(2, 2))(X)
+
+    # Conv 2 (residual)
+    X = residual.convolutional_block(X, 4, (32, 32, 128), 'stage2a', s=1)
+    X = residual.identity_block(X, 4, (32, 32, 128), 'stage2b')
+    #X = residual.identity_block(X, 4, (32, 32, 128), 'stage2c')
+
+    # Conv 3 (residual)
+    X = residual.convolutional_block(X, 3, (64, 64, 256), 'stage3a', s=2)
+    X = residual.identity_block(X, 3, (64, 64, 256), 'stage3b')
+    X = residual.identity_block(X, 3, (64, 64, 256), 'stage3c')
+    X = residual.identity_block(X, 3, (64, 64, 256), 'stage3d')
+
+    # Conv 4 (residual)
+    X = residual.convolutional_block(X, 3, (128, 128, 512), 'stage4a', s=2)
+    X = residual.identity_block(X, 3, (128, 128, 512), 'stage4b')
+    X = residual.identity_block(X, 3, (128, 128, 512), 'stage4c')
+    #X = residual.identity_block(X, 3, (128, 128, 512), 'stage4d')
+    #X = residual.identity_block(X, 3, (128, 128, 512), 'stage4e')
+
+    #X = AveragePooling2D(pool_size=(4, 4), name='avg_pool')(X)
+    X = Conv2D(512,
+               kernel_size=(4, 4),
+               strides=(1, 1),
+               padding='valid',
+               kernel_initializer=glorot_uniform(seed=0),
+               name='convend')(X)
+    X = BatchNormalization(axis=3, name='bn_convend')(X)
+    X = Activation('relu')(X)
+
+    # Flatten
+    X = Flatten()(X)
+
+    X = Dense(3,
+              activation='softmax',
+              kernel_initializer=glorot_uniform(seed=0),
+              name='y_hat')(X)
+
+    return Model(inputs=X_img, outputs=X, name=MODEL_NAME)
+
+
+def get_model_res18(img_shape):
+    X_img = Input(img_shape)
+
+    X = ZeroPadding2D((2, 2))(X_img)
+
+    # Conv 1 
+    X = Conv2D(64,
+               kernel_size=(5, 5),
+               strides=(1, 1),
+               padding='valid',
+               kernel_initializer=glorot_uniform(seed=0),
+               name='conv1')(X)
+    X = BatchNormalization(axis=3, name='bn_conv1')(X)
+    X = Activation('relu')(X)
+
+    # Conv 2 (residual)
+    X = residual.basic_block(X, 3, (64, 64), 'stage2a')
+    X = residual.basic_block(X, 3, (64, 64), 'stage2b')
+
+    # Conv 3 (residual)
+    X = residual.basic_block(X, 3, (128, 128), 'stage3a')
+    X = residual.basic_block(X, 3, (128, 128), 'stage3b')
+
+    # Conv 4 (residual)
+    X = residual.basic_block(X, 3, (256, 256), 'stage4a')
+    X = residual.basic_block(X, 3, (256, 256), 'stage4b')
+
+    # Conv 5 (residual)
+    X = residual.basic_block(X, 3, (512, 512), 'stage5a')
+    X = residual.basic_block(X, 3, (512, 512), 'stage5b')
+
+    # AveragePool
+    X = AveragePooling2D(pool_size=(2, 2), name='avg_pool')(X)
+
+    # Flatten
+    X = Flatten()(X)
+
+    X = Dense(3,
+              activation='softmax',
+              kernel_initializer=glorot_uniform(seed=0),
+              name='y_hat')(X)
+
+    return Model(inputs=X_img, outputs=X, name=MODEL_NAME)
+
+
+def train(datadir):
+    print("Load samples from train.json ...")
+    samples = util.load_samples(datadir, 'train.json')
+    m_tot = len(samples)
+    print("Got %d samples" % m_tot)
+
+    # Split the samples with the same seed every time so that the dev set
+    # is "frozen". This makes it easier to monitor and compare different models.
+    # When a good model is found, the dev set can then be removed entirely
+    # in order to fully utilize the available training data.
+    split = 0.90
+    train_samples, dev_samples = util.train_dev_split(samples, split, shuffle=True)
+    m_train = len(train_samples)
+    m_dev = len(dev_samples)
+    print("Split train/test = %.2f" % split)
+    print("Training samples: %d" % m_train)
+    print("Dev samples: %d" % m_dev)
+    print("First 5 dev samples ID's:")
+    print(' '.join([s['id'] for s in dev_samples[0:5]]))
+
+    # The minimum/maximum values of the entire training set will determine
+    # the scaling factors. We store the scaling factors to file so that they
+    # can be retrieved and re-used for predictions.
+    minmax = util.get_minmax(train_samples)
+    print("Write scaling to %s-scaling.csv" % MODEL_NAME)
+    util.save_minmax(datadir, '%s-scaling.pkl' % MODEL_NAME, minmax)
+
+    # Since we make heavy use of augmentation here, we can also augment
+    # the dev set just to smooth out validation losses during training.
+    # Augment the dev set x10.
+    m_dev *= 10
+    print("Dev samples (augmented): %d" % m_dev)
+    dev_generator = util.base_cnn_generator(
+            dev_samples, minmax, m_dev
+    )
+    X_dev, Y_dev = list(islice(dev_generator, 1))[0]
+
+    # Model + optimization parameters.
+    #model = get_model_res18((28, 28, 2))
+    model = resnet.ResnetBuilder.build_resnet_18((2, 28, 28), 3)
+    batch_size = 32
+    #opt = Adam(lr=0.0002, beta_1=0.9, beta_2=0.999, epsilon=1e-08, decay=0.0)
+    opt = SGD(lr=0.005, momentum=0.9, decay=1e-5)
+    model.compile(
+            optimizer=opt,
+            loss="categorical_crossentropy",
+            metrics = ["accuracy"]
+    )
+    model.summary()
+
+    # Callbacks
+    callbacks = [
+        ModelCheckpoint(
+            filepath=util.model_fp(datadir, '%s-weights-val_loss.h5' % MODEL_NAME),
+            verbose=1,
+            monitor='val_loss',
+            save_best_only=True
+        ),
+        ModelCheckpoint(
+            filepath=util.model_fp(datadir, '%s-weights-loss.h5' % MODEL_NAME),
+            verbose=1,
+            monitor='loss',
+            save_best_only=True
+        ),
+        ReduceLROnPlateau(
+            monitor='val_loss',
+            factor=0.1,
+            verbose=1,
+            patience=8,
+            epsilon=0.005,
+            mode='min',
+            min_lr=1e-7
+        ),
+        EarlyStopping(
+            'loss',
+            patience=30,
+            mode="min"
+        )
+    ]
+
+    # TRAIN!
+    model.fit_generator(
+            util.base_cnn_generator(train_samples, minmax, batch_size),
+            steps_per_epoch=int(4*m_train/batch_size),
+            epochs=1000,
+            validation_data=(X_dev, Y_dev),
+            callbacks=callbacks
+    )
+    
+    score = model.evaluate(X_dev, Y_dev)
+    print("")
+    print("Test loss: %.4f" % score[0])
+    print("Test accuracy: %.4f" % score[1])
+

icenet/icenet.py

@@ -0,0 +1,332 @@
+"""Standard CNN with data augmentation."""
+
+import random
+import numpy as np
+from itertools import islice
+from icenet import util
+
+from keras.models import Model
+from keras.initializers import glorot_uniform
+from keras.optimizers import Adam
+from keras import backend as K
+from keras.callbacks import ModelCheckpoint, EarlyStopping, LearningRateScheduler 
+from keras.models import load_model
+from keras.layers import (
+        Input, Dense, Activation, Conv2D, BatchNormalization, Flatten,
+        MaxPooling2D, Dropout, Concatenate, Reshape
+)
+
+
+MODEL_NAME = 'icenet'
+
+
+def get_base(base_model_filepath, postfix, last_layer, freeze=True):
+    base_model = load_model(base_model_filepath)
+    for i,l in enumerate(base_model.layers[0:last_layer+1]):
+        l.trainable = not (i == 0 or freeze)
+        l.name = '%s_%s' % (l.name, postfix)
+    #base_model.layers[last_layer].trainable = True # Should be avrg pool layer
+    return base_model.input, base_model.layers[last_layer].output
+
+
+def get_model(base_model_filepath):
+
+    print("Load base models ...")
+    X_inputs, X_base_outputs = list(zip(*[
+            get_base(base_model_filepath, 'sec%d' % i, -5, freeze=True) for i in range(9)
+    ]))
+    print("Output layer:")
+    print(X_base_outputs[0])
+    
+    #X = Concatenate(axis=2)(list(X_base_outputs))
+    #X = Reshape((nw*3, nh*3, nc))(X)
+    #_, n = K.int_shape(X_base_outputs[0])
+    _, nw, nh, nc = K.int_shape(X_base_outputs[0])
+    assert nw == 1 and nh == 1 and nc == 512
+    X = Concatenate(axis=-1)(list(X_base_outputs))
+    X = Reshape((3, 3, nc))(X)
+
+    X = Conv2D(1024,
+               kernel_size=(3, 3),
+               strides=(1, 1),
+               padding='valid',
+               activation='relu',
+               kernel_initializer=glorot_uniform(seed=0),
+               name='conv_top')(X)
+    #X = MaxPooling2D((3, 3))(X)
+
+    X = Flatten()(X)
+    X = Dropout(0.25)(X)
+    X = Dense(1024,
+              activation='relu',
+              kernel_initializer=glorot_uniform(seed=0),
+              name='fc1')(X)
+    X = Dropout(0.25)(X)
+
+    X = Dense(1,
+              activation='sigmoid',
+              kernel_initializer=glorot_uniform(seed=0),
+              name='y_hat')(X)
+
+    return Model(inputs=list(X_inputs), outputs=X, name=MODEL_NAME)
+
+
+def train(datadir):
+    base_model_name = 'base_cnn'
+    print("Load scaling from %s-scaling.csv" % base_model_name)
+    minmax = util.load_minmax(datadir, '%s-scaling.pkl' % base_model_name)
+
+    print("Load samples from train.json ...")
+    samples = util.load_samples(datadir, 'train.json')
+    m_tot = len(samples)
+    print("Got %d samples" % m_tot)
+
+    split = 0.90
+    train_samples, dev_samples = util.train_dev_split(samples, split)
+    m_train = len(train_samples)
+    m_dev = len(dev_samples)
+    print("Split train/test = %.2f" % split)
+    print("Training samples: %d" % m_train)
+    print("Dev samples: %d" % m_dev)
+    print("First 5 dev samples ID's:")
+    print(' '.join([s['id'] for s in dev_samples[0:5]]))
+
+    # Extract dev_samples
+    dev_generator = util.icenet_generator(
+            dev_samples, minmax, m_dev, crop_offset=3,
+            augment=False
+    )
+    X_dev, Y_dev = list(islice(dev_generator, 1))[0]
+
+    # Model + opt
+    def lr_schedule(epoch):
+        if epoch < 20:
+            return 0.0005
+        elif epoch < 50:
+            return 0.0002
+        elif epoch < 200:
+            return 0.00005
+        else:
+            return 0.00001
+    model = get_model(util.model_fp(datadir, '%s-weights-val_loss.h5' % base_model_name))
+    batch_size = 16
+    opt = Adam(lr=0.001, beta_1=0.9, beta_2=0.999, epsilon=1e-08, decay=0.0)
+    model.compile(
+            optimizer=opt,
+            loss="binary_crossentropy",
+            metrics = ["accuracy"]
+    )
+    _summary = []
+    model.summary(print_fn=lambda x: _summary.append(x))
+    print('\n'.join(_summary[0:8]))
+    print("...")
+    print("...")
+    print("...")
+    print("...")
+    print('\n'.join(_summary[-40:]))
+
+    # Callbacks
+    callbacks = [
+        ModelCheckpoint(
+            filepath=util.model_fp(datadir, '%s-weights-val_loss.h5' % MODEL_NAME),
+            verbose=1,
+            monitor='val_loss',
+            save_best_only=True
+        ),
+        ModelCheckpoint(
+            filepath=util.model_fp(datadir, '%s-weights-loss.h5' % MODEL_NAME),
+            verbose=1,
+            monitor='loss',
+            save_best_only=True
+        ),
+        LearningRateScheduler(lr_schedule),
+        EarlyStopping(
+            'loss',
+            patience=50,
+            mode="min"
+        )
+    ]
+
+    # TRAIN!
+    model.fit_generator(
+            util.icenet_generator(train_samples, minmax, batch_size, crop_offset=3),
+            steps_per_epoch=int(2*m_train/batch_size),
+            epochs=1000,
+            validation_data=(X_dev, Y_dev),
+            callbacks=callbacks
+    )
+    
+    score = model.evaluate(X_dev, Y_dev)
+    print("")
+    print("Test loss: %.4f" % score[0])
+    print("Test accuracy: %.4f" % score[1])
+
+    print("Make dev predictions ...")
+    y_pred = model.predict(X_dev, batch_size=32)
+    print("Write to %s-dev.csv" % MODEL_NAME)
+    util.write_preds(datadir, '%s-dev.csv' % MODEL_NAME, dev_samples, y_pred)
+
+    n_test = 30
+    print("Check invariance on %d random augmented dev samples" % n_test)
+    for i in range(n_test):
+        test_sample = random.choice(dev_samples)
+        dev_generator = util.icenet_generator(
+                [test_sample], minmax, 6*6*4*2, crop_offset=3,
+                augment=True
+        )
+        X_dev, Y_dev = list(islice(dev_generator, 1))[0]
+        y_pred = model.predict(X_dev, batch_size=32)
+        print("Dev sample %s (is_iceberg=%s): Mean = %.4f, std = %.4f, min = %.4f, max = %.4f" % (
+            test_sample.get('id'), test_sample.get('is_iceberg'),
+            np.mean(y_pred), np.std(y_pred),
+            np.min(y_pred), np.max(y_pred)
+        ))
+
+
+
+def predict(datadir):
+    print("Load model from %s-weights-loss.h5 ..." % MODEL_NAME)
+    model = load_model(util.model_fp(datadir, '%s-weights-loss.h5' % MODEL_NAME))
+
+    # Load scaling factors
+    base_model_name = 'base_cnn'
+    print("Load scaling from %s-scaling.pkl" % base_model_name)
+    minmax = util.load_minmax(datadir, '%s-scaling.pkl' % base_model_name)
+
+    target_fp = 'train.json'
+    print("Load samples from %s..." % target_fp)
+    samples = util.load_samples(datadir, target_fp)
+    m_tot = len(samples)
+    print("Got %d samples" % m_tot)
+
+    # Extract samples with generator
+    data_gen = util.icenet_generator(
+            samples, minmax, m_tot, crop_offset=3,
+            augment=False
+    )
+    X, _= list(islice(data_gen, 1))[0]
+    print("X (image) shape:")
+    print(X[0].shape)
+    
+    # Predict ...
+    print("Predict!")
+    y_pred = model.predict(X, batch_size=32)
+
+    filename = '%s-%s.csv' % (MODEL_NAME, target_fp.split('.')[0])
+    print("Write to %s" % filename)
+    util.write_preds(datadir, filename, samples, y_pred)
+
+    n_test = 20
+    print("Check invariance on %d random augmented dev samples" % n_test)
+    for i in range(n_test):
+        test_sample = random.choice(samples)
+        dev_generator = util.icenet_generator(
+                [test_sample], minmax, 6*6*4*2, crop_offset=3,
+                augment=True
+        )
+        X_dev, Y_dev = list(islice(dev_generator, 1))[0]
+        y_pred = model.predict(X_dev, batch_size=32)
+        print("Dev sample %s (is_iceberg=%s): Mean = %.4f, std = %.4f, min = %.4f, max = %.4f" % (
+            test_sample.get('id'), test_sample.get('is_iceberg'),
+            np.mean(y_pred), np.std(y_pred),
+            np.min(y_pred), np.max(y_pred)
+        ))
+
+
+def _sigmoid(z):
+    return 1.0 / (1 + np.exp(-z))
+
+
+def stretch(z, factor=10):
+    return _sigmoid(factor*(z-0.5))
+
+
+def predict_avrg(datadir):
+    print("Load model from %s-weights-loss.h5 ..." % MODEL_NAME)
+    model = load_model(util.model_fp(datadir, '%s-weights-loss.h5' % MODEL_NAME))
+
+    # Load scaling factors
+    base_model_name = 'base_cnn'
+    print("Load scaling from %s-scaling.pkl" % base_model_name)
+    minmax = util.load_minmax(datadir, '%s-scaling.pkl' % base_model_name)
+
+    target_fp = 'train.json'
+    print("Load samples from %s..." % target_fp)
+    samples = util.load_samples(datadir, target_fp)
+    random.shuffle(samples)
+    samples = samples[0:400]
+    random.shuffle(samples)
+    m_tot = len(samples)
+    print("Got %d samples" % m_tot)
+
+    n_test = 6*6*4*2
+    y_pred_first = np.zeros((len(samples), 1))
+    y_pred_avrg = np.zeros((len(samples), 1))
+    y_pred_avrg_sig10 = np.zeros((len(samples), 1))
+    y_pred_avrg_sig20 = np.zeros((len(samples), 1))
+    y_pred_avrg_sig30 = np.zeros((len(samples), 1))
+    y_pred_avrg_sig40 = np.zeros((len(samples), 1))
+    y_pred_avrg_sig50 = np.zeros((len(samples), 1))
+    y_reals = np.zeros((len(samples), 1))
+    print("Average each sample over %d augmented versions" % n_test)
+    for i,s in enumerate(samples):
+        dev_generator = util.icenet_generator(
+                [s], minmax, n_test, crop_offset=3,
+                augment=True
+        )
+        X_dev, Y_dev = list(islice(dev_generator, 1))[0]
+        y_pred = model.predict(X_dev, batch_size=n_test)
+
+        y_pred_first[i,0] = y_pred[0,0]
+        y_pred_avrg[i,0] = np.mean(y_pred)
+        y_pred_avrg_sig10[i,0] = stretch(y_pred_avrg[i,0], factor=10)
+        y_pred_avrg_sig20[i,0] = stretch(y_pred_avrg[i,0], factor=11)
+        y_pred_avrg_sig30[i,0] = stretch(y_pred_avrg[i,0], factor=12)
+        y_pred_avrg_sig40[i,0] = stretch(y_pred_avrg[i,0], factor=13)
+        y_pred_avrg_sig50[i,0] = stretch(y_pred_avrg[i,0], factor=14)
+        y_reals[i,0] = s.get('is_iceberg', 0) * 1.0
+
+        print("Sample %d: %s (iceberg=%s): Mean = %.4f, s(10) = %.4f, s(40) = %.4f std = %.4f, min = %.4f, max = %.4f" % (
+            i, s.get('id'), s.get('is_iceberg'),
+            np.mean(y_pred), 
+            y_pred_avrg_sig10[i,0], y_pred_avrg_sig40[i,0],
+            np.std(y_pred),
+            np.min(y_pred), np.max(y_pred)
+        ))
+
+    print("'First' loss: %.4f" % util.binary_crossentropy(y_reals, y_pred_first))
+    print("'Avrg' loss: %.4f" % util.binary_crossentropy(y_reals, y_pred_avrg))
+    print("'Avrg stretch(10)' loss: %.4f" % util.binary_crossentropy(y_reals, y_pred_avrg_sig10))
+    print("'Avrg stretch(20)' loss: %.4f" % util.binary_crossentropy(y_reals, y_pred_avrg_sig20))
+    print("'Avrg stretch(30)' loss: %.4f" % util.binary_crossentropy(y_reals, y_pred_avrg_sig30))
+    print("'Avrg stretch(40)' loss: %.4f" % util.binary_crossentropy(y_reals, y_pred_avrg_sig40))
+    print("'Avrg stretch(50)' loss: %.4f" % util.binary_crossentropy(y_reals, y_pred_avrg_sig50))
+
+    filename = '%s-%s-avrg.csv' % (MODEL_NAME, target_fp.split('.')[0])
+    print("Write to %s" % filename)
+    util.write_preds(datadir, filename, samples, y_pred_avrg)
+
+    filename = '%s-%s-fst.csv' % (MODEL_NAME, target_fp.split('.')[0])
+    print("Write to %s" % filename)
+    util.write_preds(datadir, filename, samples, y_pred_first)
+
+    filename = '%s-%s-s10.csv' % (MODEL_NAME, target_fp.split('.')[0])
+    print("Write to %s" % filename)
+    util.write_preds(datadir, filename, samples, y_pred_avrg_sig10)
+
+    filename = '%s-%s-s20.csv' % (MODEL_NAME, target_fp.split('.')[0])
+    print("Write to %s" % filename)
+    util.write_preds(datadir, filename, samples, y_pred_avrg_sig20)
+
+    filename = '%s-%s-s30.csv' % (MODEL_NAME, target_fp.split('.')[0])
+    print("Write to %s" % filename)
+    util.write_preds(datadir, filename, samples, y_pred_avrg_sig30)
+
+    filename = '%s-%s-s40.csv' % (MODEL_NAME, target_fp.split('.')[0])
+    print("Write to %s" % filename)
+    util.write_preds(datadir, filename, samples, y_pred_avrg_sig40)
+
+    filename = '%s-%s-s50.csv' % (MODEL_NAME, target_fp.split('.')[0])
+    print("Write to %s" % filename)
+    util.write_preds(datadir, filename, samples, y_pred_avrg_sig50)
+

icenet/residual.py

@@ -0,0 +1,190 @@
+"""Residual blocks: Conv and identity"""
+
+from keras.initializers import glorot_uniform
+from keras.layers import Add, Activation, Conv2D, BatchNormalization
+
+
+def basic_block(X, f, filters, s, label):
+    """
+    Residual net identity block.
+    
+    Arguments:
+    X -- input tensor of shape (m, n_H_prev, n_W_prev, n_C_prev)
+    f -- integer, specifying the shape of the middle CONV's window for the main path
+    filters -- python list of integers, defining the number of filters in the CONV layers of the main path
+    label -- label to use for naming
+    s -- strides
+    
+    Returns:
+    X -- output of the identity block, tensor of shape (n_H, n_W, n_C)
+    """
+    F1, F2 = filters
+    X_shortcut = X
+
+    # First component of main path 
+    X = Conv2D(
+            F1,
+            kernel_size=(f, f),
+            strides=(s, s),
+            padding='same',
+            name='conv_%s_2a' % label,
+            kernel_initializer=glorot_uniform(seed=0)
+    )(X)
+    X = BatchNormalization(axis=3, name='bn_%s_2a' % label)(X)
+    X = Activation('relu')(X)
+
+    # Second component of main path
+    X = Conv2D(
+            F2,
+            kernel_size=(f, f),
+            strides=(1, 1),
+            padding='same',
+            name='conv_%s_2b' % label,
+            kernel_initializer=glorot_uniform(seed=0)
+    )(X)
+    X = BatchNormalization(axis=3, name='bn_%s_2b' % label)(X)
+
+    # Shortcut path
+    X_shortcut = Conv2D(
+            F2,
+            kernel_size=(1, 1),
+            strides=(s, s),
+            padding='valid',
+            name='conv_%s_1' % label,
+            kernel_initializer=glorot_uniform(seed=0)
+    )(X_shortcut)
+    X_shortcut = BatchNormalization(axis=3, name='bn_%s_1' % label)(X_shortcut)
+
+    # Add shortcut value to main path, and pass it through a RELU activation
+    X = Add()([X_shortcut, X])
+    X = Activation('relu')(X)
+    
+    return X
+
+
+def identity_block(X, f, filters, label):
+    """
+    Residual net identity block.
+    
+    Arguments:
+    X -- input tensor of shape (m, n_H_prev, n_W_prev, n_C_prev)
+    f -- integer, specifying the shape of the middle CONV's window for the main path
+    filters -- python list of integers, defining the number of filters in the CONV layers of the main path
+    label -- label to use for naming
+    
+    Returns:
+    X -- output of the identity block, tensor of shape (n_H, n_W, n_C)
+    """
+    F1, F2, F3 = filters
+    X_shortcut = X
+
+    # First component of main path 
+    X = Conv2D(
+            F1,
+            kernel_size=(1, 1),
+            strides=(1,1),
+            padding='valid',
+            name='conv_%s_2a' % label,
+            kernel_initializer=glorot_uniform(seed=0)
+    )(X)
+    X = BatchNormalization(axis=3, name='bn_%s_2a' % label)(X)
+    X = Activation('relu')(X)
+
+    # Second component of main path
+    X = Conv2D(
+            F2,
+            kernel_size=(f, f),
+            strides=(1, 1),
+            padding='same',
+            name='conv_%s_2b' % label,
+            kernel_initializer=glorot_uniform(seed=0)
+    )(X)
+    X = BatchNormalization(axis=3, name='bn_%s_2b' % label)(X)
+    X = Activation('relu')(X)
+
+    # Third component of main path
+    X = Conv2D(
+            F3,
+            kernel_size=(1, 1),
+            strides=(1, 1),
+            padding='valid',
+            name='conv_%s_2c' % label,
+            kernel_initializer=glorot_uniform(seed=0)
+    )(X)
+    X = BatchNormalization(axis=3, name='bn_%s_2c' % label)(X)
+
+    # Add shortcut value to main path, and pass it through a RELU activation
+    X = Add()([X_shortcut, X])
+    X = Activation('relu')(X)
+    
+    return X
+
+
+def convolutional_block(X, f, filters, label, s=2):
+    """
+    Residual net convolutional block.
+    
+    Arguments:
+    X -- input tensor of shape (m, n_H_prev, n_W_prev, n_C_prev)
+    f -- integer, specifying the shape of the middle CONV's window for the main path
+    filters -- python list of integers, defining the number of filters in the CONV layers of the main path
+    label -- label to use for naming
+    s -- Integer, specifying the stride to be used
+    
+    Returns:
+    X -- output of the convolutional block, tensor of shape (n_H, n_W, n_C)
+    """
+    F1, F2, F3 = filters
+    X_shortcut = X
+
+    # First component of main path 
+    X = Conv2D(
+            F1,
+            kernel_size=(1, 1),
+            strides=(s,s),
+            padding='valid',
+            name='conv_%s_2a' % label,
+            kernel_initializer=glorot_uniform(seed=0)
+    )(X)
+    X = BatchNormalization(axis=3, name='bn_%s_2a' % label)(X)
+    X = Activation('relu')(X)
+
+    # Second component of main path
+    X = Conv2D(
+            F2,
+            kernel_size=(f, f),
+            strides=(1, 1),
+            padding='same',
+            name='conv_%s_2b' % label,
+            kernel_initializer=glorot_uniform(seed=0)
+    )(X)
+    X = BatchNormalization(axis=3, name='bn_%s_2b' % label)(X)
+    X = Activation('relu')(X)
+
+    # Third component of main path
+    X = Conv2D(
+            F3,
+            kernel_size=(1, 1),
+            strides=(1, 1),
+            padding='valid',
+            name='conv_%s_2c' % label,
+            kernel_initializer=glorot_uniform(seed=0)
+    )(X)
+    X = BatchNormalization(axis=3, name='bn_%s_2c' % label)(X)
+
+    # Shortcut path
+    X_shortcut = Conv2D(
+            F3,
+            kernel_size=(1, 1),
+            strides=(s, s),
+            padding='valid',
+            name='conv_%s_1' % label,
+            kernel_initializer=glorot_uniform(seed=0)
+    )(X_shortcut)
+    X_shortcut = BatchNormalization(axis=3, name='bn_%s_1' % label)(X_shortcut)
+
+    # Add shortcut value to main path, and pass it through a RELU activation
+    X = Add()([X_shortcut, X])
+    X = Activation('relu')(X)
+
+    return X

icenet/resnet.py

@@ -0,0 +1,252 @@
+from __future__ import division
+
+import six
+from keras.models import Model
+from keras.layers import (
+    Input,
+    Activation,
+    Dense,
+    Flatten
+)
+from keras.layers.convolutional import (
+    Conv2D,
+    MaxPooling2D,
+    AveragePooling2D
+)
+from keras.layers.merge import add
+from keras.layers.normalization import BatchNormalization
+from keras.regularizers import l2
+from keras import backend as K
+
+
+def _bn_relu(input):
+    """Helper to build a BN -> relu block
+    """
+    norm = BatchNormalization(axis=CHANNEL_AXIS)(input)
+    return Activation("relu")(norm)
+
+
+def _conv_bn_relu(**conv_params):
+    """Helper to build a conv -> BN -> relu block
+    """
+    filters = conv_params["filters"]
+    kernel_size = conv_params["kernel_size"]
+    strides = conv_params.setdefault("strides", (1, 1))
+    kernel_initializer = conv_params.setdefault("kernel_initializer", "he_normal")
+    padding = conv_params.setdefault("padding", "same")
+    kernel_regularizer = conv_params.setdefault("kernel_regularizer", l2(1.e-4))
+
+    def f(input):
+        conv = Conv2D(filters=filters, kernel_size=kernel_size,
+                      strides=strides, padding=padding,
+                      kernel_initializer=kernel_initializer,
+                      kernel_regularizer=kernel_regularizer)(input)
+        return _bn_relu(conv)
+
+    return f
+
+
+def _bn_relu_conv(**conv_params):
+    """Helper to build a BN -> relu -> conv block.
+    This is an improved scheme proposed in http://arxiv.org/pdf/1603.05027v2.pdf
+    """
+    filters = conv_params["filters"]
+    kernel_size = conv_params["kernel_size"]
+    strides = conv_params.setdefault("strides", (1, 1))
+    kernel_initializer = conv_params.setdefault("kernel_initializer", "he_normal")
+    padding = conv_params.setdefault("padding", "same")
+    kernel_regularizer = conv_params.setdefault("kernel_regularizer", l2(1.e-4))
+
+    def f(input):
+        activation = _bn_relu(input)
+        return Conv2D(filters=filters, kernel_size=kernel_size,
+                      strides=strides, padding=padding,
+                      kernel_initializer=kernel_initializer,
+                      kernel_regularizer=kernel_regularizer)(activation)
+
+    return f
+
+
+def _shortcut(input, residual):
+    """Adds a shortcut between input and residual block and merges them with "sum"
+    """
+    # Expand channels of shortcut to match residual.
+    # Stride appropriately to match residual (width, height)
+    # Should be int if network architecture is correctly configured.
+    input_shape = K.int_shape(input)
+    residual_shape = K.int_shape(residual)
+    stride_width = int(round(input_shape[ROW_AXIS] / residual_shape[ROW_AXIS]))
+    stride_height = int(round(input_shape[COL_AXIS] / residual_shape[COL_AXIS]))
+    equal_channels = input_shape[CHANNEL_AXIS] == residual_shape[CHANNEL_AXIS]
+
+    shortcut = input
+    # 1 X 1 conv if shape is different. Else identity.
+    if stride_width > 1 or stride_height > 1 or not equal_channels:
+        shortcut = Conv2D(filters=residual_shape[CHANNEL_AXIS],
+                          kernel_size=(1, 1),
+                          strides=(stride_width, stride_height),
+                          padding="valid",
+                          kernel_initializer="he_normal",
+                          kernel_regularizer=l2(0.0001))(input)
+
+    return add([shortcut, residual])
+
+
+def _residual_block(block_function, filters, repetitions, is_first_layer=False):
+    """Builds a residual block with repeating bottleneck blocks.
+    """
+    def f(input):
+        for i in range(repetitions):
+            init_strides = (1, 1)
+            if i == 0 and not is_first_layer:
+                init_strides = (2, 2)
+            input = block_function(filters=filters, init_strides=init_strides,
+                                   is_first_block_of_first_layer=(is_first_layer and i == 0))(input)
+        return input
+
+    return f
+
+
+def basic_block(filters, init_strides=(1, 1), is_first_block_of_first_layer=False):
+    """Basic 3 X 3 convolution blocks for use on resnets with layers <= 34.
+    Follows improved proposed scheme in http://arxiv.org/pdf/1603.05027v2.pdf
+    """
+    def f(input):
+
+        if is_first_block_of_first_layer:
+            # don't repeat bn->relu since we just did bn->relu->maxpool
+            conv1 = Conv2D(filters=filters, kernel_size=(3, 3),
+                           strides=init_strides,
+                           padding="same",
+                           kernel_initializer="he_normal",
+                           kernel_regularizer=l2(1e-4))(input)
+        else:
+            conv1 = _bn_relu_conv(filters=filters, kernel_size=(3, 3),
+                                  strides=init_strides)(input)
+
+        residual = _bn_relu_conv(filters=filters, kernel_size=(3, 3))(conv1)
+        return _shortcut(input, residual)
+
+    return f
+
+
+def bottleneck(filters, init_strides=(1, 1), is_first_block_of_first_layer=False):
+    """Bottleneck architecture for > 34 layer resnet.
+    Follows improved proposed scheme in http://arxiv.org/pdf/1603.05027v2.pdf
+
+    Returns:
+        A final conv layer of filters * 4
+    """
+    def f(input):
+
+        if is_first_block_of_first_layer:
+            # don't repeat bn->relu since we just did bn->relu->maxpool
+            conv_1_1 = Conv2D(filters=filters, kernel_size=(1, 1),
+                              strides=init_strides,
+                              padding="same",
+                              kernel_initializer="he_normal",
+                              kernel_regularizer=l2(1e-4))(input)
+        else:
+            conv_1_1 = _bn_relu_conv(filters=filters, kernel_size=(1, 1),
+                                     strides=init_strides)(input)
+
+        conv_3_3 = _bn_relu_conv(filters=filters, kernel_size=(3, 3))(conv_1_1)
+        residual = _bn_relu_conv(filters=filters * 4, kernel_size=(1, 1))(conv_3_3)
+        return _shortcut(input, residual)
+
+    return f
+
+
+def _handle_dim_ordering():
+    global ROW_AXIS
+    global COL_AXIS
+    global CHANNEL_AXIS
+    if K.image_dim_ordering() == 'tf':
+        ROW_AXIS = 1
+        COL_AXIS = 2
+        CHANNEL_AXIS = 3
+    else:
+        CHANNEL_AXIS = 1
+        ROW_AXIS = 2
+        COL_AXIS = 3
+
+
+def _get_block(identifier):
+    if isinstance(identifier, six.string_types):
+        res = globals().get(identifier)
+        if not res:
+            raise ValueError('Invalid {}'.format(identifier))
+        return res
+    return identifier
+
+
+class ResnetBuilder(object):
+    @staticmethod
+    def build(input_shape, num_outputs, block_fn, repetitions):
+        """Builds a custom ResNet like architecture.
+
+        Args:
+            input_shape: The input shape in the form (nb_channels, nb_rows, nb_cols)
+            num_outputs: The number of outputs at final softmax layer
+            block_fn: The block function to use. This is either `basic_block` or `bottleneck`.
+                The original paper used basic_block for layers < 50
+            repetitions: Number of repetitions of various block units.
+                At each block unit, the number of filters are doubled and the input size is halved
+
+        Returns:
+            The keras `Model`.
+        """
+        _handle_dim_ordering()
+        if len(input_shape) != 3:
+            raise Exception("Input shape should be a tuple (nb_channels, nb_rows, nb_cols)")
+
+        # Permute dimension order if necessary
+        if K.image_dim_ordering() == 'tf':
+            input_shape = (input_shape[1], input_shape[2], input_shape[0])
+
+        # Load function from str if needed.
+        block_fn = _get_block(block_fn)
+
+        input = Input(shape=input_shape)
+        conv1 = _conv_bn_relu(filters=64, kernel_size=(7, 7), strides=(2, 2))(input)
+        pool1 = MaxPooling2D(pool_size=(3, 3), strides=(2, 2), padding="same")(conv1)
+
+        block = pool1
+        filters = 64
+        for i, r in enumerate(repetitions):
+            block = _residual_block(block_fn, filters=filters, repetitions=r, is_first_layer=(i == 0))(block)
+            filters *= 2
+
+        # Last activation
+        block = _bn_relu(block)
+
+        # Classifier block
+        block_shape = K.int_shape(block)
+        pool2 = AveragePooling2D(pool_size=(block_shape[ROW_AXIS], block_shape[COL_AXIS]),
+                                 strides=(1, 1))(block)
+        flatten1 = Flatten()(pool2)
+        dense = Dense(units=num_outputs, kernel_initializer="he_normal",
+                      activation="softmax")(flatten1)
+
+        model = Model(inputs=input, outputs=dense)
+        return model
+
+    @staticmethod
+    def build_resnet_18(input_shape, num_outputs):
+        return ResnetBuilder.build(input_shape, num_outputs, basic_block, [2, 2, 2, 2])
+
+    @staticmethod
+    def build_resnet_34(input_shape, num_outputs):
+        return ResnetBuilder.build(input_shape, num_outputs, basic_block, [3, 4, 6, 3])
+
+    @staticmethod
+    def build_resnet_50(input_shape, num_outputs):
+        return ResnetBuilder.build(input_shape, num_outputs, bottleneck, [3, 4, 6, 3])
+
+    @staticmethod
+    def build_resnet_101(input_shape, num_outputs):
+        return ResnetBuilder.build(input_shape, num_outputs, bottleneck, [3, 4, 23, 3])
+
+    @staticmethod
+    def build_resnet_152(input_shape, num_outputs):
+        return ResnetBuilder.build(input_shape, num_outputs, bottleneck, [3, 8, 36, 3])

icenet/util.py

@@ -0,0 +1,229 @@
+import os
+import json
+import numpy as np
+import random
+
+from itertools import cycle
+from skimage import transform
+from sklearn.externals import joblib
+from sklearn.metrics import log_loss
+
+
+def binary_crossentropy(y_true, y_pred):
+    return log_loss(y_true.flatten(), y_pred.flatten())
+
+
+def get_bands(band1, band2):
+    return np.dstack((
+        np.array(band1).reshape(75, 75),
+        np.array(band2).reshape(75, 75)
+    ))
+
+
+def get_angle(a):
+    if a == 'na':
+        return None
+    else:
+        return a
+
+
+def get_label(l):
+    if l is None:
+        return l
+    else:
+        return int(l)
+
+
+def unpack(samples):
+    images, angles, labels = zip(*[
+        (get_bands(s['band_1'], s['band_2']),
+        get_angle(s['inc_angle']),
+        get_label(s.get('is_iceberg')))
+        for s in samples
+    ])
+    return images, angles, labels
+
+
+def load_samples(datadir, filename):
+    with open(os.path.join(datadir, filename)) as f:
+        samples = json.load(f)
+    return samples
+
+
+def train_dev_split(samples, split, shuffle=True):
+    split_idx = int(split*len(samples))
+    if shuffle:
+        random.seed(0)
+        random.shuffle(samples)
+    return (samples[0:split_idx], samples[split_idx:])
+
+
+def write_preds(datadir, filename, samples, y_pred):
+    with open(os.path.join(datadir, 'predictions', filename), 'w') as f:
+        f.write('id,is_iceberg\n')
+        f.write('\n'.join([
+            '%s,%.6f' % (samples[i].get('id'), y_pred[i][0])
+            for i in range(len(samples))
+        ]))
+        f.write('\n')
+
+
+def model_fp(datadir, filename):
+    return os.path.join(datadir, 'params', filename)
+
+
+def save_minmax(datadir, filename, minmax):
+    joblib.dump(minmax, os.path.join(datadir, 'params', filename))
+
+
+def load_minmax(datadir, filename):
+    return joblib.load(os.path.join(datadir, 'params', filename))
+
+
+def get_minmax(samples):
+    # Extract global mins/max' for normalization
+    images, angles, labels = unpack(samples)
+    b1_min = min([im[:, :, 0].min() for im in images])
+    b1_max = max([im[:, :, 0].max() for im in images])
+    b2_min = min([im[:, :, 1].min() for im in images])
+    b2_max = max([im[:, :, 1].max() for im in images])
+    a_min = min(angles, key=lambda x: x or 1e9)
+    a_max = max(angles, key=lambda x: x or -1e9)
+    print("Band 1 min/max: %.2f, %.2f" % (b1_min, b1_max))
+    print("Band 2 min/max: %.2f, %.2f" % (b2_min, b2_max))
+    print("Angles min/max: %.2f, %.2f" % (a_min, a_max))
+    return (b1_min, b1_max, b2_min, b2_max, a_min, a_max)
+
+
+def base_cnn_generator(samples, minmax, batch_size, verbose=False):
+    img_dim = 75
+    window = 28
+    r = int(window/2)
+    padding = int(window/np.sqrt(2))
+    if verbose:
+        print("window is %d px" % window)
+        print("padding is %d px" % padding)
+
+    # Fuck the angles.
+    images, _, labels = unpack(samples)
+    b1_min, b1_max, b2_min, b2_max, _, _ = minmax
+
+    # Yield batches forever, cycling the samples with augmentation:
+    # Random rotation and crop, and random Bernoulli mirroring.
+    batch_images = []
+    batch_labels = []
+    for image, label in cycle(zip(images, labels)):
+        # Scale each band independently
+        img = np.dstack((
+            (image[:, :, 0] - b1_min) / (b1_max - b1_min),
+            (image[:, :, 1] - b2_min) / (b2_max - b2_min),
+        ))
+
+        # Save the maximum of both bands within this image
+        b1_peak = img[:, :, 0].max()
+        b2_peak = img[:, :, 1].max()
+
+        # Rotate around random midpoint
+        row_mid = random.randint(padding, img_dim-padding-1)
+        col_mid = random.randint(padding, img_dim-padding-1)
+        if verbose:
+            print("\nSample %d" % (len(batch_labels)+1))
+            print("row/col mid: %d/%d" % (row_mid, col_mid))
+            print("Mid values: %.2f, %.2f" % (img[row_mid, col_mid, 0], img[row_mid, col_mid, 1]))
+
+        # NOTE: sklearn.transform.rotate behaves incorrectly when using kwargs 'center'
+        # Clip to mid with padding before rotating.
+        img = img[row_mid-padding:row_mid+padding, col_mid-padding:col_mid+padding, :]
+        rot = random.random() * 360 - 180
+        img = transform.rotate(img, rot, order=1)
+        if verbose:
+            print("Rotation: %.2f deg" % rot)
+            print("Mid values: %.2f, %.2f" % (img[padding, padding, 0], img[padding, padding, 1]))
+        img = img[padding-r:padding+r, padding-r:padding+r, :]
+
+        # Mirror
+        if random.randint(0, 1) and False:
+            img = np.fliplr(img)
+            if verbose:
+                print("Sample was mirrored")
+
+        # Label by checking that subimage contains maximum of either band
+        # If the subimage contains a max value that is within 95% of the image
+        # max, consider it to contain the feature.
+        sub_b1_peak = img[:, :, 0].max()
+        sub_b2_peak = img[:, :, 1].max()
+        l = int(sub_b1_peak/b1_peak > 0.95 and
+                sub_b2_peak/b2_peak > 0.95) * (label + 1)
+        if verbose:
+            print("Label is %d" % l)
+        categorical_label = tuple(int(i==l) for i in range(3))
+
+        batch_images.append(img)
+        batch_labels.append(categorical_label)
+        if len(batch_labels) == batch_size:
+            yield (
+                np.array(batch_images, dtype=np.float32),
+                np.array(batch_labels, dtype=np.float32)
+            )
+            batch_images = []
+            batch_labels = []
+
+
+def augment_image(img, mirror, rots):
+    if mirror:
+        img = np.fliplr(img)
+    img = np.rot90(img, k=rots)
+    return img
+
+
+def icenet_generator(samples, minmax, batch_size, crop_offset=3,
+                   augment=True, verbose=False):
+    img_dim = 75
+    mid_x = mid_y = 75.0 / 2
+    crop_dim = img_dim - 2*crop_offset
+    window = 28
+    r = int(window/2)
+    
+    images, _, labels = unpack(samples)
+    b1_min, b1_max, b2_min, b2_max, _, _ = minmax
+
+    # Yield batches forever
+    # A batch consists of an input-output tuple
+    #   ([X_image_sec1, X_image_sec2, ... , X_image_sec9], Y_labels)
+    # where the input image is split over 9 overlapping sections, starting
+    # from upper left in row-major order
+    batch_image_sections = [[] for _ in range(9)]
+    batch_labels = []
+    for image, label in cycle(zip(images, labels)):
+        # Scale each band independently
+        img = np.dstack((
+            (image[:, :, 0] - b1_min) / (b1_max - b1_min),
+            (image[:, :, 1] - b2_min) / (b2_max - b2_min),
+        ))
+
+        # Crop with random offset from midpoint
+        row_offset = crop_offset * (2*random.random() - 1)
+        col_offset = crop_offset * (2*random.random() - 1)
+        if not augment:
+            row_offset = col_offset = 0
+        row = int(mid_y + row_offset - crop_dim/2)
+        col = int(mid_x + col_offset - crop_dim/2)
+        img = img[row:row+crop_dim, col:col+crop_dim, :]
+
+        if augment:
+            img = augment_image(img, random.randint(0, 1), random.randint(0, 3))
+
+        # Append each section of the image
+        for i, ridx in enumerate([r, int(mid_y), crop_dim-r]):
+            for j, cidx in enumerate([r, int(mid_x), crop_dim-r]):
+                batch_image_sections[i*3+j].append(img[ridx-r:ridx+r, cidx-r:cidx+r, :])
+
+        batch_labels.append(label)
+        if len(batch_labels) == batch_size:
+            yield (
+                [np.array(sec_imgs, dtype=np.float32) for sec_imgs in batch_image_sections],
+                np.array(batch_labels, dtype=np.float32)
+            )
+            batch_image_sections = [[] for _ in range(9)]
+            batch_labels = []
+

plot.py

@@ -0,0 +1,104 @@
+#!/usr/bin/env python3
+"""Plot training/test data. Channels appear side by side."""
+
+import os
+import argparse
+import json
+import random
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from icenet import util
+    
+
+def plot_sample(outdir, sample, minmax, interactive=False):
+    fig, (ax1, ax2) = plt.subplots(1, 2, sharex='col', sharey='row')
+    fig.suptitle('\n'.join([
+        "Id: %s (original)" % sample.get('id'),
+        "Iceberg? %s" % sample.get('is_iceberg', '-'),
+        "Incident angle: %s" % repr(sample.get('inc_angle'))])
+    )
+
+    b1_min, b1_max, b2_min, b2_max, _, _ = minmax
+    b1 = (np.array(sample['band_1']).reshape(75, 75) - b1_min) / (b1_max - b1_min)
+    b2 = (np.array(sample['band_2']).reshape(75, 75) - b2_min) / (b2_max - b2_min)
+    ax1.imshow(b1, vmin=0.2, vmax=0.8, aspect='auto')
+    ax1.set_title("HH")
+    ax2.imshow(b2, vmin=0.2, vmax=0.8, aspect='auto')
+    ax2.set_title("HV")
+
+    #fig.tight_layout()
+    fig.subplots_adjust(top=0.80)
+    if interactive:
+        plt.show()
+    else:
+        fig.savefig(os.path.join(outdir, "%s.png" % sample['id']))
+    plt.close('all')
+
+
+def plot_angle_hist(outdir, samples, interactive=False):
+    angles = [s['inc_angle'] for s in samples if s['inc_angle'] != 'na']
+
+    hist, bins = np.histogram(angles, bins=50)
+    width = 0.7 * (bins[1] - bins[0])
+    center = (bins[:-1] + bins[1:]) / 2
+
+    fig = plt.figure()
+    plt.bar(center, hist, align='center', width=width)
+    plt.title("%d/%d valid angles" % (len(angles), len(samples)))
+    if interactive:
+        plt.show()
+    else:
+        fig.savefig(os.path.join(outdir, "angle_hist.png"))
+    plt.close('all')
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument(
+            '-o',
+            required=False,
+            help="Output directory for samples."
+    )
+    parser.add_argument(
+            '-n',
+            type=int,
+            required=False,
+            help="Number of samples to plot (randomized)"
+    )
+    parser.add_argument(
+            '-i',
+            required=False,
+            action='store_true',
+            help="Show plots (and don't store)"
+    )
+    parser.add_argument(
+            'samples_file',
+            help="JSON file with samples"
+    )
+    args = parser.parse_args()
+
+    outdir = args.o or '.'
+    if not os.path.isdir(outdir):
+        raise Exception("Output directory does not exist")
+
+    print("Loading samples ...")
+    with open(args.samples_file) as f:
+        samples = json.load(f)
+    print("%d samples in set" % len(samples))
+
+    if args.n:
+        print("Pick %d random samples ..." % args.n)
+        random.shuffle(samples)
+        samples = samples[0:args.n]
+
+    minmax = util.get_minmax(samples)
+    for i, s in enumerate(samples):
+        print("Plotting sample %d/%d" % (i+1, len(samples)))
+        plot_sample(outdir, s, minmax, interactive=args.i or False)
+
+    print("Plot angle histogram ..")
+    plot_angle_hist(outdir, samples, interactive=args.i or False)
+
+    print("Done.")

plot_augmented.py

@@ -0,0 +1,133 @@
+#!/usr/bin/env python3
+"""Plot image produced by data generators in icenet/utils.py"""
+
+import os
+import argparse
+import json
+import random
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from icenet import util
+from plot import plot_sample
+from itertools import islice
+
+
+def plot_base_cnn_samples(outdir, samples, interactive=False):
+    m = len(samples)
+    minmax = util.get_minmax(samples)
+    imgs, _ = list(islice(util.base_cnn_generator(
+            samples, minmax, m, shuffle=False, verbose=True
+    ), 1))[0]
+    for i in range(m):
+        print("Plotting sample %d/%d" % (i+1, m))
+        sample = samples[i]
+        img = imgs[i, :, :, :]
+
+        fig, (ax1, ax2) = plt.subplots(1, 2, sharex='col', sharey='row')
+        fig.suptitle('\n'.join([
+            "Id: %s (rotated and cropped)" % sample.get('id'),
+            "Iceberg? %s" % sample.get('is_iceberg', '-'),
+            "Incident angle: %s" % repr(sample.get('inc_angle'))])
+        )
+
+        ax1.imshow(img[:, :, 0], vmin=0.2, vmax=0.8, aspect='auto')
+        ax1.set_title("HH")
+        ax2.imshow(img[:, :, 1], vmin=0.2, vmax=0.8, aspect='auto')
+        ax2.set_title("HV")
+
+        #fig.tight_layout()
+        fig.subplots_adjust(top=0.80)
+        plot_sample(outdir, sample, minmax, interactive=interactive)
+        if not interactive:
+            fig.savefig(os.path.join(outdir, "%s_base.png" % sample['id']))
+        plt.close('all')
+
+
+def plot_icenet_samples(outdir, samples, interactive=False):
+    m = len(samples)
+    minmax = util.get_minmax(samples)
+    sections, y = list(islice(util.icenet_generator(
+            samples, minmax, m, augment=False, verbose=True
+    ), 1))[0]
+    for i in range(m):
+        print("Plotting sample %d/%d" % (i+1, m))
+        sample = samples[i]
+
+        for b in range(2):
+            fig, axes = plt.subplots(3, 3, sharex='col', sharey='row')
+            fig.suptitle('\n'.join([
+                "Id: %s (rotated and cropped, band %d)" % (sample.get('id'), b+1),
+                "Iceberg? %s" % sample.get('is_iceberg', '-'),
+                "Incident angle: %s" % repr(sample.get('inc_angle'))])
+            )
+
+            for ridx in range(3):
+                for cidx in range(3):
+                    axes[ridx, cidx].imshow(
+                            sections[ridx*3+cidx][i, :, :, b], 
+                            vmin=0.2,
+                            vmax=0.8,
+                            aspect='auto'
+                    )
+
+            fig.subplots_adjust(top=0.80)
+            if not interactive:
+                fig.savefig(os.path.join(outdir, "%s_icenet_b%d.png" % (sample['id'], b)))
+
+        plot_sample(outdir, sample, minmax, interactive=interactive)
+        plt.close('all')
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument(
+            '-o',
+            required=False,
+            help="Output directory for samples."
+    )
+    parser.add_argument(
+            '-g',
+            required=True,
+            choices=['base', 'icenet'],
+            help="Generator to use (base|icenet)."
+    )
+    parser.add_argument(
+            '-n',
+            type=int,
+            required=False,
+            help="Number of samples to plot (randomized)"
+    )
+    parser.add_argument(
+            '-i',
+            required=False,
+            action='store_true',
+            help="Show plots (and don't store)"
+    )
+    parser.add_argument(
+            'samples_file',
+            help="JSON file with samples"
+    )
+    args = parser.parse_args()
+
+    outdir = args.o or '.'
+    if not os.path.isdir(outdir):
+        raise Exception("Output directory does not exist")
+
+    print("Loading samples ...")
+    with open(args.samples_file) as f:
+        samples = json.load(f)
+    print("%d samples in set" % len(samples))
+
+    if args.n:
+        print("Pick %d random samples ..." % args.n)
+        random.shuffle(samples)
+        samples = samples[0:args.n]
+
+    if args.g == 'base':
+        plot_base_cnn_samples(outdir, samples, interactive=args.i or False)
+    else:
+        plot_icenet_samples(outdir, samples, interactive=args.i or False)
+
+    print("Done.")

run_model.py

@@ -0,0 +1,29 @@
+#!/usr/bin/env python3
+"""Run a model."""
+
+import os
+import argparse
+from icenet import MODELS
+    
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument(
+            '-d',
+            required=True,
+            help="Directory containing {train,test}.json."
+    )
+    parser.add_argument(
+            'target_fn',
+            help=("Target <module>.<fn> in the icenet package. "
+                  "The function must take the data directory as argument.")
+    )
+    args = parser.parse_args()
+    if not os.path.isdir(args.d):
+        raise Exception("Output directory does not exist")
+
+    model, fn = args.target_fn.split('.')
+    print("Running model %s ..." % model)
+    fn = getattr(MODELS[model], fn)
+    fn(args.d)
+    print("Done.")

run_model_docker.sh

@@ -0,0 +1,5 @@
+#!/usr/bin/env bash
+set -e
+
+docker build -t iceberg:$(cat VERSION) .
+docker run --name iceberg --rm -v $(pwd)/data:/data iceberg:$(cat VERSION) $@

show_predicted.sh

@@ -0,0 +1,11 @@
+#!/usr/bin/env bash
+
+# Args:
+#   path to predictions (CSV)
+#   path to images
+
+tail -n+2 $1 | shuf | while read -r i; do
+  ID=$(echo $i | cut -d ',' -f1)
+  PRED=$(echo $i | cut -d ',' -f2)
+  echo "ID $ID: is_iceberg = $PRED" && gwenview $2/$ID.png 2>/dev/null
+done