{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Setup","metadata":{}},{"cell_type":"code","source":"!pip install scikit-image\n!pip install -q segmentation_models --no-index --find-links=file:///kaggle/input/segmentation-models/segmentation_models\n!pip install -q albumentations \n%env SM_FRAMEWORK=tf.keras\n    \nimport os\nimport numpy as np\nfrom matplotlib import animation\nimport matplotlib.pyplot as plt\nfrom IPython import display\nimport math\nimport pickle\nfrom scipy import ndimage\nfrom tqdm import tqdm \nfrom contextlib import nullcontext\nimport segmentation_models\nfrom tqdm import tqdm\nimport gc \nimport albumentations as A\n\nimport tensorflow as tf\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras.layers import Conv2D, Activation,Conv2DTranspose,MaxPool2D,Concatenate, BatchNormalization\nfrom tensorflow.keras.models import Model\nfrom tensorflow.python.keras import losses\nfrom tensorflow.keras.applications.efficientnet import EfficientNetB0\nfrom tensorflow.keras.callbacks import ModelCheckpoint,ReduceLROnPlateau\nfrom tensorflow.keras.optimizers import Adam\nfrom tensorflow.keras.optimizers.schedules import ExponentialDecay\n\n\n\nBASE_DIR = '/kaggle/input/google-research-identify-contrails-reduce-global-warming'\nSUBMISSION_DIR = \"/kaggle/working\"\nN_TIMES_BEFORE = 4\nIMG_W=256\nIMG_H=256\nN_TIMES = 8\nN_CHANNELS = 3\n\nBANDS = ['band_10.npy', 'band_14.npy', 'band_15.npy', 'band_16.npy',\n'band_08.npy', 'band_09.npy', 'band_13.npy', 'band_11.npy', 'band_12.npy']\n\nOUTPUT_DIR = \"/kaggle/working\"\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-06-06T03:06:59.844566Z","iopub.execute_input":"2023-06-06T03:06:59.844955Z","iopub.status.idle":"2023-06-06T03:08:12.667563Z","shell.execute_reply.started":"2023-06-06T03:06:59.844926Z","shell.execute_reply":"2023-06-06T03:08:12.666196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DO_TRAINING = True\nPREPROCESS_DATA = True","metadata":{"execution":{"iopub.status.busy":"2023-06-06T03:08:12.670321Z","iopub.execute_input":"2023-06-06T03:08:12.670728Z","iopub.status.idle":"2023-06-06T03:08:12.677238Z","shell.execute_reply.started":"2023-06-06T03:08:12.670677Z","shell.execute_reply":"2023-06-06T03:08:12.676332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dataset visualization","metadata":{}},{"cell_type":"code","source":"sample_record_id = '1704010292581573769'\n\nsample_band = np.load(os.path.join(BASE_DIR,\"train\",sample_record_id,\"band_08.npy\"))\nsample_ground_truth = np.load(os.path.join(BASE_DIR,\"train\",sample_record_id,\"human_pixel_masks.npy\"))\n\nprint(\"band shape:\")\nprint(sample_band.shape)\nprint(\"ground truth shape:\")\nprint(sample_ground_truth.shape)\n\n    \nfile_names = os.listdir(os.path.join(BASE_DIR,\"train\",sample_record_id))\n\nband_files = [file_name for file_name in file_names if file_name.startswith(\"band\")]\n\n\nn_band_files = len(band_files)\n\nplt.figure(figsize=(16, 6))\n\nfor i in range(n_band_files):\n        file_path = os.path.join(BASE_DIR,\"train\",sample_record_id, band_files[i])\n        with open(file_path, 'rb') as file:\n            band = np.load(file)\n            band = band[...,N_TIMES_BEFORE]\n            plt.subplot(2, n_band_files, i+1)\n            plt.imshow(band, interpolation='none')\n            \n\nfor i in range(n_band_files):\n        file_path = os.path.join(BASE_DIR,\"train\",sample_record_id, band_files[i])\n        with open(file_path, 'rb') as file:\n            band = np.load(file)\n            band = band[...,N_TIMES_BEFORE-1]\n            plt.subplot(2, n_band_files, i+1)\n            plt.imshow(band, interpolation='none')\n\nplt.subplot(2, n_band_files, n_band_files+1)\nplt.imshow(sample_ground_truth[...,0], interpolation='none')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-06-06T03:08:12.678896Z","iopub.execute_input":"2023-06-06T03:08:12.679219Z","iopub.status.idle":"2023-06-06T03:08:14.001022Z","shell.execute_reply.started":"2023-06-06T03:08:12.679189Z","shell.execute_reply":"2023-06-06T03:08:14.000187Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preprocessing","metadata":{}},{"cell_type":"code","source":"_T11_BOUNDS = (243, 303)\n_CLOUD_TOP_TDIFF_BOUNDS = (-4, 5)\n_TDIFF_BOUNDS = (-4, 2)\n\ndef normalize_range(data, bounds):\n    \"\"\"Maps data to the range [0, 1].\"\"\"\n    return (data - bounds[0]) / (bounds[1] - bounds[0])\n\ndef process_record(subset,record_id):\n    with open(os.path.join(BASE_DIR,subset, record_id, 'band_11.npy'), 'rb') as f:\n        band11 = np.load(f)\n    with open(os.path.join(BASE_DIR,subset, record_id, 'band_14.npy'), 'rb') as f:\n        band14 = np.load(f)\n    with open(os.path.join(BASE_DIR,subset, record_id, 'band_15.npy'), 'rb') as f:\n        band15 = np.load(f)\n\n    r = normalize_range(band15 - band14, _TDIFF_BOUNDS)\n    g = normalize_range(band14 - band11, _CLOUD_TOP_TDIFF_BOUNDS)\n    b = normalize_range(band14, _T11_BOUNDS)\n    false_color = np.clip(np.stack([r, g, b], axis=2), 0, 1)\n    false_color = false_color[...,N_TIMES_BEFORE]\n    with open(os.path.join(BASE_DIR,subset,record_id,\"human_pixel_masks.npy\"), 'rb') as file:\n        ground_truth = np.load(file)\n        ground_truth = ground_truth[...,0]\n    return false_color,ground_truth\n    \n\ndef process_batch(subset,batch_size,records):\n    gc.collect()\n    n_records = min(len(records),batch_size)\n    data = np.zeros(shape=(n_records,IMG_W,IMG_H,N_CHANNELS),dtype=\"float32\")\n    data_y = np.zeros(shape=(n_records,IMG_W,IMG_H),dtype=\"bool\")\n    for i,record_id in tqdm(enumerate(records[:n_records]),total = n_records,desc=f\"Processing {subset} data\"):\n        data[i],data_y[i] = process_record(subset,record_id)\n        \n    return n_records,data,data_y","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-06-06T03:08:14.003529Z","iopub.execute_input":"2023-06-06T03:08:14.004459Z","iopub.status.idle":"2023-06-06T03:08:14.33753Z","shell.execute_reply.started":"2023-06-06T03:08:14.004428Z","shell.execute_reply":"2023-06-06T03:08:14.336447Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if (PREPROCESS_DATA):\n    batch_size = 256\n    subset = 'train'\n    records = os.listdir(os.path.join(BASE_DIR,subset))\n\n    n_batches = math.ceil(len(records)/batch_size)\n\n    for batch in range(n_batches):\n        print(f\"Processing batch {batch+1}/{n_batches}...\")\n        processed,data,data_y = process_batch(subset,batch_size,records)\n        \n        if not os.path.exists(os.path.join(OUTPUT_DIR,subset)):\n            os.makedirs(os.path.join(OUTPUT_DIR,subset))\n        if not os.path.exists(os.path.join(OUTPUT_DIR,f\"{subset}_y\")):\n            os.makedirs(os.path.join(OUTPUT_DIR,f\"{subset}_y\"))\n        np.save(os.path.join(OUTPUT_DIR,subset,f\"{batch}.npy\"), data)\n        np.save(os.path.join(OUTPUT_DIR,f\"{subset}_y\",f\"{batch}.npy\"), data_y)\n        \n        records = records[processed:]\n        \n    batch_size = 256\n    subset = 'validation'\n    records = os.listdir(os.path.join(BASE_DIR,subset))\n\n    n_batches = math.ceil(len(records)/batch_size)\n\n    for batch in range(n_batches):\n        print(f\"Processing batch {batch+1}/{n_batches}...\")\n        processed,data,data_y = process_batch(subset,batch_size,records)\n        if not os.path.exists(os.path.join(OUTPUT_DIR,subset)):\n            os.makedirs(os.path.join(OUTPUT_DIR,subset))\n        if not os.path.exists(os.path.join(OUTPUT_DIR,f\"{subset}_y\")):\n            os.makedirs(os.path.join(OUTPUT_DIR,f\"{subset}_y\"))\n        np.save(os.path.join(OUTPUT_DIR,subset,f\"{batch}.npy\"), data)\n        np.save(os.path.join(OUTPUT_DIR,f\"{subset}_y\",f\"{batch}.npy\"), data_y)\n        \n        \n        records = records[processed:]","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-06-06T04:29:16.960894Z","iopub.execute_input":"2023-06-06T04:29:16.961338Z","iopub.status.idle":"2023-06-06T04:32:02.395897Z","shell.execute_reply.started":"2023-06-06T04:29:16.961284Z","shell.execute_reply":"2023-06-06T04:32:02.393886Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Utility Functions","metadata":{}},{"cell_type":"code","source":"def count_train_steps():\n    batch_files = os.listdir(os.path.join(OUTPUT_DIR,\"train\"))\n    train_size = 0\n    for i,batch_file in enumerate(batch_files):\n        with open(os.path.join(OUTPUT_DIR,\"train\",batch_file),\"rb\") as file:\n            train_X = np.load(file)\n        train_size += len(train_X)\n    return train_size\n\ndef count_validation_steps():\n    batch_files = os.listdir(os.path.join(OUTPUT_DIR,\"validation\"))\n    validation_size = 0\n    for i,batch_file in enumerate(batch_files):\n        with open(os.path.join(OUTPUT_DIR,\"validation\",batch_file),\"rb\") as file:\n            validation_X = np.load(file) \n        validation_size += len(validation_X)\n    return validation_size","metadata":{"execution":{"iopub.status.busy":"2023-06-06T04:35:36.645482Z","iopub.execute_input":"2023-06-06T04:35:36.645888Z","iopub.status.idle":"2023-06-06T04:35:36.655849Z","shell.execute_reply.started":"2023-06-06T04:35:36.645858Z","shell.execute_reply":"2023-06-06T04:35:36.654965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model Definition","metadata":{}},{"cell_type":"code","source":"EPOCHS = 20\n\nbatch_size = 8\n\nN_AUG = 1 #Number of augmentations per image\n\nLR = 0.016 #initial lr\n\nbackbone_name = 'efficientnetb0'\n    \nencoder_weights='/kaggle/input/contrail-model/efficientnet-b0_weights_tf_dim_ordering_tf_kernels_autoaugment_notop.h5'\n\nsteps_per_epoch = (N_AUG+1)*math.ceil(count_train_steps()/batch_size)\nvalidation_steps = math.ceil(count_validation_steps()/batch_size)\n\n\nmodel = segmentation_models.Unet(backbone_name=backbone_name,input_shape=(256, 256, 3),\n                                 decoder_block_type='transpose',classes = 1,\n                                 activation='sigmoid',encoder_weights = encoder_weights)\n\ndice_loss = segmentation_models.losses.DiceLoss()\nfocal_loss = segmentation_models.losses.BinaryFocalLoss()\ntotal_loss = dice_loss + (1 * focal_loss)\n\n\nlr_decayed_fn = ExponentialDecay(\n    initial_learning_rate=LR, decay_steps=steps_per_epoch,decay_rate=0.7\n)\n\noptimizer = Adam(lr_decayed_fn)\n\nmetrics = [segmentation_models.metrics.IOUScore(threshold=0.5), segmentation_models.metrics.FScore(threshold=0.5)]\n\nmodel.compile(loss=total_loss, optimizer=optimizer, metrics=[metrics])","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-06-06T04:35:38.156341Z","iopub.execute_input":"2023-06-06T04:35:38.156708Z","iopub.status.idle":"2023-06-06T04:37:12.400599Z","shell.execute_reply.started":"2023-06-06T04:35:38.156678Z","shell.execute_reply":"2023-06-06T04:37:12.399364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#model.summary()","metadata":{"execution":{"iopub.status.busy":"2023-06-06T03:08:37.435805Z","iopub.status.idle":"2023-06-06T03:08:37.436271Z","shell.execute_reply.started":"2023-06-06T03:08:37.436035Z","shell.execute_reply":"2023-06-06T03:08:37.436057Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training","metadata":{}},{"cell_type":"code","source":"\n\ntransform = A.Compose([\n    A.ToFloat(max_value=1),\n    A.OneOf([\n        A.HorizontalFlip(p=0.3),\n        A.VerticalFlip(p=0.3),\n    ]),\n    A.ShiftScaleRotate(p=1),\n    A.Resize(256,256),    \n],p=1,\nadditional_targets={\n        'image1': 'image'\n})\n\n\n\ndef data_gen(batch_size):\n    batch_files = os.listdir(os.path.join(OUTPUT_DIR, \"train\"))\n    while True:\n        for batch_file in batch_files:\n            with open(os.path.join(OUTPUT_DIR, \"train\", batch_file), \"rb\") as file:\n                train_X = np.load(file)\n            with open(os.path.join(OUTPUT_DIR, \"train_y\", batch_file), \"rb\") as file:\n                train_y = np.load(file)\n            n_records = len(train_X)\n            train_X = train_X.astype(\"float32\")\n            train_y = train_y.astype(\"float32\")\n            for i in range(N_AUG):\n                augmentations_X = np.zeros_like(train_X)\n                augmentations_y = np.zeros_like(train_y)\n                for j in range(n_records):\n                    transformations = transform(image=train_X[j],image1=train_y[j])\n                    augmentations_X[j],augmentations_y[j] = transformations['image'],transformations['image1']\n                train_X= np.concatenate((train_X, augmentations_X))\n                train_y=np.concatenate((train_y, augmentations_y))\n            p = np.random.permutation(len(train_X))\n            train_X = train_X[p]\n            train_y = train_y[p]\n            if (len(train_X)>batch_size):\n                train_X = np.split(train_X, len(train_X) / batch_size)\n                train_y = np.split(train_y, len(train_y) / batch_size)\n                for X, y in zip(train_X, train_y):\n                    yield X, y\n            else:\n                yield train_X, train_y\n            gc.collect()\n      \n    \ndef validation_data_gen(batch_size):\n    batch_files = os.listdir(os.path.join(OUTPUT_DIR, \"validation\"))\n    while True:\n        for i, batch_file in enumerate(batch_files):\n            with open(os.path.join(OUTPUT_DIR, \"validation\", batch_file), \"rb\") as file:\n                val_X = np.load(file)\n            with open(os.path.join(OUTPUT_DIR, \"validation_y\", batch_file), \"rb\") as file:\n                val_y = np.load(file)\n            val_X = val_X.astype(\"float32\")\n            val_y = val_y.astype(\"float32\")\n            if (len(val_X)>batch_size):\n                val_X = np.split(val_X, len(val_X) / batch_size)\n                val_y = np.split(val_y, len(val_y) / batch_size)\n                for X, y in zip(val_X, val_y):\n                    yield X, y\n            else:\n                yield val_X, val_y\n            gc.collect()\n\n\ntrain_gen = data_gen(batch_size)\nval_gen = validation_data_gen(batch_size)","metadata":{"execution":{"iopub.status.busy":"2023-06-06T04:39:22.332627Z","iopub.execute_input":"2023-06-06T04:39:22.333051Z","iopub.status.idle":"2023-06-06T04:39:22.352521Z","shell.execute_reply.started":"2023-06-06T04:39:22.333016Z","shell.execute_reply":"2023-06-06T04:39:22.351349Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ngc.collect()\n\n            \nCHECKPOINT_DIR = os.path.join(OUTPUT_DIR,\"best_checkpoint.h5\")\n\nclass LearningRateLogger(tf.keras.callbacks.Callback):\n    def on_epoch_end(self, epoch, logs=None):\n        current_learning_rate = self.model.optimizer.lr.numpy()\n        print(f\"\\nLearning Rate at the end of epoch {epoch + 1}: {current_learning_rate}\")\n\n        \ncallbacks = [\n    ModelCheckpoint(CHECKPOINT_DIR, save_weights_only=True, save_best_only=True, mode='min'),\n    LearningRateLogger()\n]\n\nif (DO_TRAINING):\n    history = model.fit(x=train_gen,batch_size=batch_size,epochs=EPOCHS,\n                        steps_per_epoch=steps_per_epoch,\n                        validation_steps=validation_steps,\n                        validation_data=val_gen,\n                        callbacks=callbacks)\n\n\n    \n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.load_weights(CHECKPOINT_DIR)\nloss = model.evaluate(validation_data_gen(batch_size),steps=validation_steps)\nprint(\"Validation Loss:\", loss[2])","metadata":{"execution":{"iopub.status.busy":"2023-06-06T06:21:14.461264Z","iopub.execute_input":"2023-06-06T06:21:14.46241Z","iopub.status.idle":"2023-06-06T06:21:35.607696Z","shell.execute_reply.started":"2023-06-06T06:21:14.462364Z","shell.execute_reply":"2023-06-06T06:21:35.606455Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Predictions on test dataset and Submission","metadata":{}},{"cell_type":"code","source":"import csv\n\ndef process_test_record(subset,record_id):\n    with open(os.path.join(BASE_DIR,subset, record_id, 'band_11.npy'), 'rb') as f:\n        band11 = np.load(f)\n    with open(os.path.join(BASE_DIR,subset, record_id, 'band_14.npy'), 'rb') as f:\n        band14 = np.load(f)\n    with open(os.path.join(BASE_DIR,subset, record_id, 'band_15.npy'), 'rb') as f:\n        band15 = np.load(f)\n\n    r = normalize_range(band15 - band14, _TDIFF_BOUNDS)\n    g = normalize_range(band14 - band11, _CLOUD_TOP_TDIFF_BOUNDS)\n    b = normalize_range(band14, _T11_BOUNDS)\n    false_color = np.clip(np.stack([r, g, b], axis=2), 0, 1)\n    false_color = false_color[...,N_TIMES_BEFORE]\n    return false_color\n\ndef process_test_batch(subset,cur_batch,batch_size,records):\n    \n    n_records = min(len(records),batch_size)\n    data = np.zeros(shape=(n_records,IMG_W,IMG_H,N_CHANNELS),dtype=\"float32\")\n    for i,record_id in tqdm(enumerate(records[:n_records]),total = n_records,desc=f\"Processing {subset} data\"):\n        data[i] = process_test_record(subset,record_id)\n\n    \n    return n_records,data\n\ndef rle_encode(x):\n    mask = np.where(x.T.flatten() == 1)[0]\n    rle = []\n    prev = -2\n    for i in mask:\n        if i==prev+1:\n            rle[-1]+=1\n        else:\n            rle.extend((i+1,1))\n        prev = i\n    if not rle:\n        text = '-'\n    else:\n        text = str(rle).replace(\"[\", \"\").replace(\"]\", \"\").replace(\",\", \"\")\n    return text\n\n\ndef clean(binary_prediction): \n    #removes single dots in prediction(improvement is not significant)\n    k = np.ones((10,10,1))\n    conv = ndimage.convolve(binary_prediction,k,mode='constant',cval=0.25 ) > 2\n    return np.logical_and(binary_prediction,conv)\n\n\ncsvfile = open(os.path.join(SUBMISSION_DIR,'submission.csv'), 'w')\n\nfield_names= ['record_id', 'encoded_pixels']\nwriter = csv.DictWriter(csvfile, fieldnames=field_names)\nwriter.writeheader()\n\nrecords = os.listdir(os.path.join(BASE_DIR,'test'))\nbatch_size = 256\nbatches = math.ceil(len(records)/batch_size)\n\nfor batch in range(batches):\n    processed,test_X = process_test_batch('test', 1, batch_size,records)\n    predictions = model.predict(test_X)\n    \n    submission = []                 \n    for i,record in enumerate(records[:processed]):\n\n        binary_prediction = np.vectorize(lambda x: 1 if x>=0.05 else 0)(predictions[i])\n        binary_prediction = clean(binary_prediction)\n        submission.append({})\n        submission[-1]['record_id']=record\n        submission[-1]['encoded_pixels']=rle_encode(binary_prediction)\n\n    writer.writerows(submission)\n    records = records[processed:]\n    gc.collect()\ncsvfile.close()","metadata":{"execution":{"iopub.status.busy":"2023-06-06T06:21:37.448951Z","iopub.execute_input":"2023-06-06T06:21:37.449368Z","iopub.status.idle":"2023-06-06T06:21:40.28834Z","shell.execute_reply.started":"2023-06-06T06:21:37.449332Z","shell.execute_reply":"2023-06-06T06:21:40.287172Z"},"trusted":true},"execution_count":null,"outputs":[]}]}