{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport xml.etree.ElementTree as ET\nimport matplotlib.pyplot as plt\nfrom PIL import Image\nimport os\nimport glob\nimport time\nimport zipfile\nimport operator\nimport collections\nfrom skimage import transform as tform\nfrom sklearn import decomposition, mixture, cluster\nfrom scipy.spatial.distance import pdist, squareform","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"root_images = \"../input/all-dogs/all-dogs/\"\nroot_annots = \"../input/annotation/Annotation/\"\n\nall_images = os.listdir(\"../input/all-dogs/all-dogs/\")\nbreeds = glob.glob('../input/annotation/Annotation/*')\n\nannotation=[]\nfor b in breeds:\n    annotation+=glob.glob(b+\"/*\")\n\nbreed_map={}\nfor annot in annotation:\n    breed=annot.split(\"/\")[-2]\n    index=breed.split(\"-\")[0]\n    breed_map.setdefault(index,breed)\n\nbreed_folders = glob.glob('../input/annotation/Annotation/*')\nbreed_index_name_map = {}\nfor breed_folder in breed_folders:\n    full_breed_name = breed_folder.split(\"/\")[-1]\n    breed_index  = full_breed_name.split(\"-\")[ 0]\n    breed_name   = full_breed_name.split(\"-\")[-1]\n    breed_index_name_map[breed_index] = breed_name\n\nnum_dog_images = len(all_images)\nnum_dog_breeds = len(breed_map)\nnum_avg_images_per_breed = num_dog_images / num_dog_breeds\nprint('Total %d dog images of %d different breeds (on average %.1f images per breed)' %(num_dog_images,num_dog_breeds,num_avg_images_per_breed))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def bounding_box(image_filename):\n    bpath=root_annots+str(breed_map[image_filename.split(\"_\")[0]])+\"/\"+str(image_filename.split(\".\")[0])\n    tree = ET.parse(bpath)\n    root = tree.getroot()\n    objects = root.findall('object')\n    for o in objects:\n        bndbox = o.find('bndbox') # reading bound box\n        xmin = int(bndbox.find('xmin').text)\n        ymin = int(bndbox.find('ymin').text)\n        xmax = int(bndbox.find('xmax').text)\n        ymax = int(bndbox.find('ymax').text)\n        \n    return (xmin,ymin,xmax,ymax)\n\ndef expand_bounding_box(bbox, orig_image, expand_margin_fraction=0.1):\n    im_width, im_height = orig_image.size\n\n    bbox_w = bbox[2]-bbox[0]\n    bbox_h = bbox[3]-bbox[1]\n    \n    xmin = max(0, bbox[0] - 0.5*expand_margin_fraction*bbox_w)\n    ymin = max(0, bbox[1] - 0.5*expand_margin_fraction*bbox_h)\n    \n    xmax = min(im_width,  bbox[2] + 0.5*expand_margin_fraction*bbox_w)\n    ymax = min(im_height, bbox[3] + 0.5*expand_margin_fraction*bbox_h)\n    \n    return [xmin,ymin,xmax,ymax]\n\nnum_rows = 5\nnum_cols = 9\nnum_images_to_show = num_rows*num_cols\nselected_images = np.random.choice(all_images, size=num_images_to_show, replace=False)\n\nexpand_margin_fraction = 0.125\n\nfig = plt.figure(figsize=(30,16)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.94, bottom=0.02, hspace=0.05, wspace=0.05); \nplt.suptitle('figure 1: raw samples from dataset', fontsize=30)\nfor k, image_filename in enumerate(selected_images):\n    bbox = bounding_box(image_filename)\n    orig_image = Image.open(os.path.join(root_images, image_filename))\n    \n    bbox_expanded = expand_bounding_box(bbox, orig_image, expand_margin_fraction=expand_margin_fraction)\n    cropped_image = orig_image.crop(bbox_expanded)\n    \n    plt.subplot(num_rows,num_cols,k+1); plt.imshow(cropped_image); plt.axis(\"off\")\nfig.savefig('figure_1.png')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# create a small dataset of resized images\nnum_images_in_dataset = num_dog_images\n#num_images_in_dataset = 16384\n#num_images_in_dataset = 8192\n\nexpantion_factor = 20\n\nstart_time = time.time()\n\nimage_dimention = 64 + expantion_factor\nresize_shape = (image_dimention, image_dimention)\n\nselected_images = np.random.choice(all_images, size=num_images_in_dataset, replace=False)\n\n# create a matrix to hold all images\nimage_dataset_4D_matrix = np.zeros((image_dimention,image_dimention,3,num_images_in_dataset), dtype=np.uint8)\nimage_breed_label_list = []\n\n# fill up the matrix with images\nfor k, image_filename in enumerate(selected_images):\n    bbox = bounding_box(image_filename)\n    orig_image = Image.open(os.path.join(root_images, image_filename))\n    bbox_expanded = expand_bounding_box(bbox, orig_image, expand_margin_fraction=expand_margin_fraction)\n    cropped_image = orig_image.crop(bbox_expanded)\n    resized_image = tform.resize(np.array(cropped_image), resize_shape, preserve_range=True).astype(np.uint8)\n\n    image_dataset_4D_matrix[:,:,:,k] = resized_image\n    image_breed_label_list.append(breed_index_name_map[image_filename.split('_')[0]])\n    \ntraining_duration_min = (time.time()-start_time)/60\nprint('finished collecting dataset. took %.1f minutes' %(training_duration_min))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Show collected rescaled dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"num_rows = 6\nnum_cols = 9\nnum_images_to_show = num_rows*num_cols\nselected_inds = np.random.choice(num_images_in_dataset,size=num_images_to_show,replace=False)\n\nfig = plt.figure(figsize=(30,22)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.93, bottom=0.02, hspace=0.13, wspace=0.05); \nplt.suptitle('figure 2: rescaled loosly cropped images', fontsize=30)\nfor k, image_ind in enumerate(selected_inds):\n    dog_image = image_dataset_4D_matrix[:,:,:,image_ind]\n    dog_breed = image_breed_label_list[image_ind]\n    \n    plt.subplot(num_rows,num_cols,k+1); plt.imshow(dog_image); plt.title(dog_breed, fontsize=16); plt.axis(\"off\")\nfig.savefig('figure_2.png')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn import preprocessing\n\nbreed_encoder = preprocessing.LabelEncoder()\nbreed_encoder.fit(image_breed_label_list)\n\nnum_breeds = len(list(breed_encoder.classes_))\nprint('total num breeds in dataset is %d' %(num_breeds))\n\n# short test:\ntest_labels = image_breed_label_list[:7]\nbreed_index = breed_encoder.transform(test_labels)\nprint(breed_index, breed_index.shape)\ntest_labels_hat = list(breed_encoder.inverse_transform(breed_index))\nprint(test_labels)\nprint(test_labels_hat)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Define Autoencoder Architecture and Learning params"},{"metadata":{"trusted":true},"cell_type":"code","source":"import keras\nfrom keras.models import Model\nfrom keras.layers import GlobalMaxPooling2D, GlobalAveragePooling2D\nfrom keras.layers import Input, Dropout, Dense, Conv2D, MaxPooling2D, AveragePooling2D, UpSampling2D, Reshape, LeakyReLU, Flatten\nfrom keras.preprocessing.image import ImageDataGenerator\nfrom keras.callbacks import LearningRateScheduler\nfrom keras.layers.normalization import BatchNormalization\nfrom keras.layers.noise import GaussianNoise, GaussianDropout\nfrom keras.regularizers import l1,l2,l1_l2\nfrom keras import optimizers","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# from keras.applications.vgg19 import VGG19\n\n# vgg_encoder = VGG19(include_top=False, weights='imagenet', input_shape=(64,64,3))\n\n# # Creating dictionary that maps layer names to the layers\n# layer_dict = dict([(layer.name, layer) for layer in vgg_encoder.layers])\n\n# # Getting output tensor of the last VGG layer that we want to include\n# desired_output = layer_dict['block4_conv2'].output\n\n# # Creating new model. Please note that this is NOT a Sequential() model.\n# from keras.models import Model\n# custom_model = Model(input=vgg_model.input, output=desired_output)\n\n# # Make sure that the pre-trained bottom layers are not trainable\n# for layer in custom_model.layers[:7]:\n#     layer.trainable = False","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# hyperparams\nkernel_reg = 1e-7\nactivity_reg = 1e-6\n\nencoder_output_channel_size = 256\n\nmultiplicative_noise_sigma = 0.3\nadditive_noise_sigma = 0.04\ndropout_noise_rate = 0.04\n\nleaky_relu_slope = 0.33\nleaky_relu = lambda x: LeakyReLU(alpha=leaky_relu_slope)(x)\n\nmultiplicative_gaussian_noise_level = (multiplicative_noise_sigma**2)/(1+multiplicative_noise_sigma**2)\nadditive_gaussian_noise_level = additive_noise_sigma\n\n# encoder\ninput_image = Input(shape=(64, 64, 3), name='input_image')\nx = BatchNormalization(name='encoder_BN_1_1')(input_image)\nx = Conv2D( 64, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='encoder_conv1_1')(x)\nx = BatchNormalization(name='encoder_BN_1_2')(x)\nx = Conv2D( 64, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='encoder_conv1_2')(x)\nx = MaxPooling2D((2, 2), padding='same', name='encoder_pool1')(x)\n\nx = BatchNormalization(name='encoder_BN_2_1')(x)\nx = Conv2D(128, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='encoder_conv2_1')(x)\nx = BatchNormalization(name='encoder_BN_2_2')(x)\nx = Conv2D(128, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='encoder_conv2_2')(x)\nx = MaxPooling2D((2, 2), padding='same', name='encoder_pool2')(x)\n\nx = BatchNormalization(name='encoder_BN_3_1')(x)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='encoder_conv3_1')(x)\nx = BatchNormalization(name='encoder_BN_3_2')(x)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='encoder_conv3_2')(x)\nx = MaxPooling2D((2, 2), padding='same', name='encoder_pool3')(x)\n\nx = BatchNormalization(name='encoder_BN_4_1')(x)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='encoder_conv4_1')(x)\nx = BatchNormalization(name='encoder_BN_4_2')(x)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='encoder_conv4_2')(x)\nx = MaxPooling2D((2, 2), padding='same', name='encoder_pool4')(x)\n\n# bottleneck layer. Try to make it sparse using L1 activity regularization\nencoder_output = Conv2D(encoder_output_channel_size, (1, 1), activation='linear', kernel_regularizer=l2(kernel_reg), activity_regularizer=l1(activity_reg), name='encoder_output')(x)\n\n# add some noise during training to force interpulation smoothness in the latent space (use both additive and multiplicative noise and dropout noise)\nencoder_output = GaussianDropout(multiplicative_gaussian_noise_level, name='multiplicative_noise')(encoder_output)\nencoder_output = GaussianNoise(additive_gaussian_noise_level        , name='additive_noise'      )(encoder_output)\nencoder_output = Dropout(dropout_noise_rate                         , name='dropout_noise'       )(encoder_output)\n\nencoder = Model(input_image, encoder_output, name='encoder')\nprint('---------------------------------------------------------------------------------------------------')\nprint(' Encoder:')\nprint('---------------------------------------------------------------------------------------------------')\nencoder.summary()\nprint('---------------------------------------------------------------------------------------------------')\n\n\n# latent auxilary classifier\nexternal_encoder_output = Input(shape=(4,4,encoder_output_channel_size), name='extenral_image_rep')\n\n# CNN\nx = BatchNormalization(name='latent_aux_BN_1')(external_encoder_output)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='latent_aux_conv1')(x)\nx = Dropout(0.85, name='latent_aux_dropout_1')(x)\nx = BatchNormalization(name='latent_aux_BN_2')(x)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='latent_aux_conv2')(x)\nx = GlobalMaxPooling2D(name='global_pool')(x)\n\n# FCN\n'''\nx = Flatten(name='flatten_encoder_latent')(external_encoder_output)\nx = BatchNormalization(name='latent_aux_BN_1')(x)\nx = Dense(768, activation=leaky_relu, name='latent_aux_h1', kernel_regularizer=l2(kernel_reg))(x)\nx = Dropout(0.7, name='dropout_1')(x)\nx = BatchNormalization(name='latent_aux_BN_2')(x)\nx = Dense(768, activation=leaky_relu, name='latent_aux_h2', kernel_regularizer=l2(kernel_reg))(x)\nx = Dropout(0.7, name='dropout_2')(x)\nx = BatchNormalization(name='latent_aux_BN_3')(x)\nx = Dense(768, activation=leaky_relu, name='latent_aux_h3', kernel_regularizer=l2(kernel_reg))(x)\nx = Dropout(0.7, name='dropout_3')(x)\n'''\n\nlatent_aux_classifier_output = Dense(num_breeds, activation='softmax', name='latent_aux_pred')(x)\nlatent_aux_classifier = Model(external_encoder_output, latent_aux_classifier_output, name='latent_aux_classifier')\nprint('---------------------------------------------------------------------------------------------------')\nprint(' Latent Auxilary Classifier:')\nprint('---------------------------------------------------------------------------------------------------')\nlatent_aux_classifier.summary()\nprint('---------------------------------------------------------------------------------------------------')\n\n\n# decoder\n#external_encoder_output = Input(shape=(4,4,encoder_output_channel_size), name='extenral_image_rep')\n\nx = BatchNormalization(name='decoder_BN_1_1')(external_encoder_output)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='decoder_conv1_1')(x)\nx = BatchNormalization(name='decoder_BN_1_2')(external_encoder_output)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='decoder_conv1_2')(x)\nx = UpSampling2D((2, 2), interpolation='bilinear', name='decoder_upsample_1')(x)\n\nx = BatchNormalization(name='decoder_BN_2_1')(x)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='decoder_conv2_1')(x)\nx = BatchNormalization(name='decoder_BN_2_2')(x)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='decoder_conv2_2')(x)\nx = UpSampling2D((2, 2), interpolation='bilinear', name='decoder_upsample_2')(x)\n\nx = BatchNormalization(name='decoder_BN_3_1')(x)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='decoder_conv3_1')(x)\nx = BatchNormalization(name='decoder_BN_3_2')(x)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='decoder_conv3_2')(x)\nx = UpSampling2D((2, 2), interpolation='bilinear', name='decoder_upsample_3')(x)\n\nx = BatchNormalization(name='decoder_BN_4_1')(x)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='decoder_conv4_1')(x)\nx = BatchNormalization(name='decoder_BN_4_2')(x)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='decoder_conv4_2')(x)\nx = UpSampling2D((2, 2), interpolation='bilinear', name='decoder_upsample_4')(x)\n\nx = BatchNormalization(name='decoder_BN_5_1')(x)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='decoder_conv5_1')(x)\nx = BatchNormalization(name='decoder_BN_5_2')(x)\nx = Conv2D(256, (3, 3), activation=leaky_relu, padding='same', kernel_regularizer=l2(kernel_reg), name='decoder_conv5_2')(x)\n\noutput_image = Conv2D(3, (3, 3), activation='linear', padding='same', name='generated_image')(x)\n\ndecoder = Model(external_encoder_output, output_image, name='decoder')\nprint('---------------------------------------------------------------------------------------------------')\nprint(' Decoder:')\nprint('---------------------------------------------------------------------------------------------------')\ndecoder.summary()\nprint('---------------------------------------------------------------------------------------------------')\n\n# autoencoder\nencoder_latent = encoder(input_image)\nautoencoder = Model(input_image, decoder(encoder_latent))\nprint('---------------------------------------------------------------------------------------------------')\nprint(' Autoencoder:')\nprint('---------------------------------------------------------------------------------------------------')\nautoencoder.summary()\nprint('---------------------------------------------------------------------------------------------------')\n\n# autoencoder with latent aux classifier\nautoencoder_with_latent_aux_classifier = Model(input_image, outputs=[decoder(encoder_latent), latent_aux_classifier(encoder_latent)])\n#autoencoder_with_latent_aux_classifier = Model(input_image, outputs=[decoder(encoder(input_image)), latent_aux_classifier(encoder(input_image))])\nprint('---------------------------------------------------------------------------------------------------')\nprint(' Autoencoder with latent auxilary classifer:')\nprint('---------------------------------------------------------------------------------------------------')\nautoencoder_with_latent_aux_classifier.summary()\nprint('---------------------------------------------------------------------------------------------------')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Generator function"},{"metadata":{"trusted":true},"cell_type":"code","source":"def normalize_image(orig_scale_images):\n    # map [0,255] range to [16/256,240/256]\n    normlized_images = (orig_scale_images.astype(np.float32) * (224/255) + 16) / 256\n    return normlized_images\n\ndef unnormalize_image(normlized_images):\n    # map from the range [16/256,240/256] back to [0,255]\n    orig_scale_images = 255 * ((normlized_images - 16/256) / (224/256))\n    return orig_scale_images\n\ndef generate_batches(possible_inds_list, batch_size=64, random_crops=True, random_flips=True):\n    num_possible_images = len(possible_inds_list)\n    assert(num_possible_images >= batch_size)\n    \n    while True:\n        curr_batch = np.zeros((batch_size,64,64,3))\n        curr_batch_lables_list = []\n        selected_images_for_batch = np.random.choice(possible_inds_list,size=batch_size,replace=False)\n        for k, selected_image_ind in enumerate(selected_images_for_batch):\n            if random_crops:\n                h_start = np.random.randint(expantion_factor)\n                w_start = np.random.randint(expantion_factor)\n            else:\n                h_start = int(expantion_factor/2)\n                w_start = int(expantion_factor/2)\n    \n            h_end = h_start + 64\n            w_end = w_start + 64\n                        \n            if random_flips:\n                image = np.fliplr(image_dataset_4D_matrix[h_start:h_end,w_start:w_end,:,selected_image_ind])\n            else:\n                image = image_dataset_4D_matrix[h_start:h_end,w_start:w_end,:,selected_image_ind]\n            \n            curr_batch[k,:,:,:] = normalize_image(image)\n            curr_batch_lables_list.append(image_breed_label_list[selected_image_ind])\n            \n        curr_batch_lables_vec = breed_encoder.transform(curr_batch_lables_list)\n        curr_batch_lables = keras.utils.to_categorical(curr_batch_lables_vec, num_breeds)\n\n        yield (curr_batch, [curr_batch, curr_batch_lables])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Train autoencoder"},{"metadata":{"trusted":true},"cell_type":"code","source":"num_epochs = 2200\nbatch_size = 16\nlearning_rate = 0.0001\n\nvalid_data_fraction = 0.125\ntrain_steps_per_epoch = 96\nvalid_steps_per_epoch = 16\n\nvalid_cutoff = int((1-valid_data_fraction) * num_images_in_dataset)\ntrain_inds = [x for x in range(valid_cutoff)]\nvalid_inds = [x for x in range(valid_cutoff,num_images_in_dataset)]\n\ntrain_data_generator = generate_batches(train_inds, batch_size=batch_size, random_crops=True, random_flips=True)\nvalid_data_generator = generate_batches(valid_inds, batch_size=batch_size, random_crops=True, random_flips=True)\n\nlosses_to_use = ['mae','categorical_crossentropy']\nloss_weights_to_use = [1.0,0.01]\nmetrics_to_use = ['accuracy']\noptimizer_to_use = optimizers.Nadam(lr=learning_rate)\nautoencoder_with_latent_aux_classifier.compile(optimizer=optimizer_to_use, loss=losses_to_use, loss_weights=loss_weights_to_use, metrics=metrics_to_use)\n\nstart_time = time.time()\nhistory = autoencoder_with_latent_aux_classifier.fit_generator(generator=train_data_generator, epochs=num_epochs, steps_per_epoch=train_steps_per_epoch,\n                                                               validation_data=valid_data_generator, validation_steps=valid_steps_per_epoch)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"training_duration_sec = time.time()-start_time\ntraining_duration_hours = training_duration_sec / 3600\ntraining_duration_remaining_minutes = 60 * (training_duration_hours - int(training_duration_hours))\nprint('finished training Autoencoder. took in total %d hours and %d minutes' %(training_duration_hours, training_duration_remaining_minutes))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"losses_to_show = ['loss','decoder_loss','latent_aux_classifier_loss','latent_aux_classifier_acc']\nylim_ranges = {}\nylim_ranges['loss'] = [0.01,0.27]\nylim_ranges['decoder_loss'] = [0.01,0.17]\nylim_ranges['latent_aux_classifier_loss'] = [0.3,5.8]\nylim_ranges['latent_aux_classifier_acc'] = [0.02,1.01]\nnum_rows_in_subplot = len(losses_to_show)\n\n# show learning curves\nepoch_number = np.arange(1,num_epochs+1)\nfig = plt.figure(figsize=(min(30,int(5+0.5*num_epochs)),25)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.92, bottom=0.02, hspace=0.25, wspace=0.05);\nplt.suptitle('figure 3: Auto-Encoder learning curves', fontsize=24);\nfor k, loss_string in enumerate(losses_to_show):\n    plt.subplot(num_rows_in_subplot,1,k+1);\n    \n    final_train_loss = np.array(history.history[loss_string][-25:]).mean()\n    final_valid_loss = np.array(history.history['val_'+loss_string][-25:]).mean()\n    plt.title('final (train,valid) %s = (%.4f,%.4f)' %(loss_string, final_train_loss, final_valid_loss),fontsize=22)\n    plt.plot(epoch_number, history.history[loss_string],'b')\n    plt.plot(epoch_number, history.history['val_'+loss_string],'g')\n    plt.legend(['train','valid'], fontsize=18);\n    plt.xlabel('num iterations', fontsize=16); plt.ylabel(loss_string, fontsize=16)\n    plt.ylim(ylim_ranges[loss_string][0],ylim_ranges[loss_string][1])\n    \nfig.savefig('figure_3.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Gather data for GMM learning"},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data_generator = generate_batches(list(range(num_images_in_dataset)), batch_size=2048, random_crops=True, random_flips=True)\n\nh_start = int(expantion_factor/2); h_end = h_start + 64;\nw_start = int(expantion_factor/2); w_end = w_start + 64;\n\nX_normlized_center_crop = np.transpose(normalize_image(image_dataset_4D_matrix[h_start:h_end,w_start:w_end,:,:]),[3,0,1,2])\nX_normlized_center_crop_flipped = np.flip(np.transpose(normalize_image(image_dataset_4D_matrix[h_start:h_end,w_start:w_end,:,:]),[3,0,1,2]),axis=2)\n#X_normlized_random_crop_random_flip = next(all_data_generator)[0]\n\nX_normlized = np.concatenate((X_normlized_center_crop,X_normlized_center_crop_flipped,next(all_data_generator)[0],\n                                                                                      next(all_data_generator)[0],\n                                                                                      next(all_data_generator)[0],\n                                                                                      next(all_data_generator)[0],\n                                                                                      next(all_data_generator)[0],\n                                                                                      next(all_data_generator)[0],\n                                                                                      next(all_data_generator)[0],\n                                                                                      next(all_data_generator)[0],\n                                                                                      next(all_data_generator)[0],\n                                                                                      next(all_data_generator)[0]),axis=0)\n\nprint(X_normlized.shape, X_normlized.min(), X_normlized.max())\ndel X_normlized_center_crop, X_normlized_center_crop_flipped","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# quick verification of flipping and data that will be used for GMM fitting\nfig = plt.figure(figsize=(30,12)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.94, bottom=0.02, hspace=0.08, wspace=0.05); \nplt.suptitle('figure 4: flipping and crops quick check', fontsize=30)\nfor k in range(8):\n    rand_ind_center_crop = np.random.randint(num_images_in_dataset)\n    rand_ind_center_crop_flipped = rand_ind_center_crop + num_images_in_dataset\n    rand_ind_random_crop = np.random.randint(2*num_images_in_dataset, X_normlized.shape[0])\n    plt.subplot(3,8,k+1+ 0);  plt.imshow(unnormalize_image(X_normlized[rand_ind_center_crop,:,:,:]).astype(np.uint8)); plt.axis('off')\n    plt.subplot(3,8,k+1+ 8);  plt.imshow(unnormalize_image(X_normlized[rand_ind_center_crop_flipped,:,:,:]).astype(np.uint8)); plt.axis('off')\n    plt.subplot(3,8,k+1+ 16); plt.imshow(unnormalize_image(X_normlized[rand_ind_random_crop,:,:,:]).astype(np.uint8)); plt.axis('off')\nfig.savefig('figure_4.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Show some Autoencoder reconstructions"},{"metadata":{"trusted":true},"cell_type":"code","source":"# show several model reconstructions\nnum_rows = 6\nnum_cols = 9\nnum_images_to_show = num_rows*num_cols\n\nselected_inds = np.random.choice(X_normlized.shape[0], size=num_images_to_show, replace=False)\n\nX_rec_autuencoder = np.transpose(unnormalize_image(autoencoder.predict(X_normlized[selected_inds])), [1,2,3,0])\nprint(X_rec_autuencoder.shape, X_rec_autuencoder.mean(), X_rec_autuencoder.std())\n\nfig = plt.figure(figsize=(30,20)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.95, bottom=0.02, hspace=0.05, wspace=0.05); \nplt.suptitle('figure 5: AE reconstuctions', fontsize=30)\nfor k in range(num_images_to_show):\n    doglike_image = X_rec_autuencoder[:,:,:,k]\n    doglike_image[doglike_image > 255] = 255\n    doglike_image[doglike_image <   0] =   0\n    \n    plt.subplot(num_rows,num_cols,k+1); plt.imshow(doglike_image.astype(np.uint8)); plt.axis(\"off\")\nfig.savefig('figure_5.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Show histograms of encoder latent space units"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_rep_autoencoder = encoder.predict(X_normlized).reshape((X_normlized.shape[0],-1))\nprint(X_rep_autoencoder.shape)\n\nfig = plt.figure(figsize=(30,12)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.92, bottom=0.02, hspace=0.08, wspace=0.12); \nplt.suptitle('figure 6: AE unit activations', fontsize=30)\nfor k in range(36):\n    selected_ind = np.random.randint(X_rep_autoencoder.shape[-1])\n    unit_activations = X_rep_autoencoder[:,selected_ind]\n    range_limit = max(abs(unit_activations.min()), abs(unit_activations.max()))\n    activation_range = np.linspace(-range_limit,range_limit,100)\n    plt.subplot(4,9,k+1); plt.hist(unit_activations, bins=activation_range, log=True);\nfig.savefig('figure_6.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Add noise pertubations in AE latent space"},{"metadata":{"trusted":true},"cell_type":"code","source":"# calc std for each latent direction\nae_latent_rep = encoder.predict(X_normlized)\nae_latent_std = ae_latent_rep.std(axis=0,keepdims=True)\n\nprint(ae_latent_std.shape, ae_latent_std.mean(), ae_latent_std.std())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Show distribution of euclidean distances between samples in latent space"},{"metadata":{"trusted":true},"cell_type":"code","source":"# show the distribution of distnaces between different samples in the latent space\nsubset_size = 4500 # to limit compute complexity\nsubset_inds = np.random.choice(ae_latent_rep.shape[0], size=subset_size, replace=False)\n\nae_latent_rep_table = ae_latent_rep.reshape((ae_latent_rep.shape[0],-1))\nae_latent_space_distances = pdist(ae_latent_rep_table[subset_inds,:], 'euclidean')\n\nfig = plt.figure(figsize=(20,10)); \nplt.subplots_adjust(left=0.05, right=0.95, top=0.93, bottom=0.05, hspace=0.15, wspace=0.05); \nplt.suptitle('figure 7: euclidean distance distribution in latent space', fontsize=30)\nplt.subplot(2,1,1); plt.hist(ae_latent_space_distances, bins=200);\nplt.subplot(2,1,2); plt.hist(ae_latent_space_distances, bins=200, log=True); plt.xlabel('euclidean distance');\nfig.savefig('figure_7.png')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"d_min = np.percentile(ae_latent_space_distances, 1)\nd_max = np.percentile(ae_latent_space_distances,99)\nd_mean = ae_latent_space_distances.mean()\nd_std = ae_latent_space_distances.std()\nprint('98%s of euclidian distnaces range between %.4f to %.4f (mean = %.4f, std = %.4f)' %('%', d_min,d_max,d_mean,d_std))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Show noisy images around an image with various noise levels"},{"metadata":{"trusted":true},"cell_type":"code","source":"noise_levels_to_show = [0.2,0.5,0.75,1.0,1.25,1.5,2.0] # these are in units of std per feature\nnum_pertubation_per_noise_level = 4\n\nnum_cols = len(noise_levels_to_show)\nnum_rows = num_pertubation_per_noise_level\n\nselected_image_rep = ae_latent_rep[np.random.randint(ae_latent_rep.shape[0])][np.newaxis,:,:,:]\n\nfig = plt.figure(figsize=(30,18)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.92, bottom=0.02, hspace=0.15, wspace=0.05); \nplt.suptitle('figure 8: image + random noise in latent space (a)', fontsize=30)\nfor col, noise_level in enumerate(noise_levels_to_show):\n    for row in range(num_pertubation_per_noise_level):\n        noise_to_add = noise_level * np.random.normal(scale=ae_latent_std, size=ae_latent_std.shape)\n        noisy_image_rep = selected_image_rep + noise_to_add\n        noisy_image = unnormalize_image(decoder.predict(noisy_image_rep)[0])\n        noisy_image[noisy_image > 255] = 255\n        noisy_image[noisy_image <   0] =   0\n        \n        src_to_noisy_dist = pdist(np.concatenate((selected_image_rep.reshape((1,-1)), noisy_image_rep.reshape((1,-1))), axis=0))[0]\n        plt.subplot(num_rows,num_cols,col+1+row*num_cols); plt.imshow(noisy_image.astype(np.uint8)); plt.axis(\"off\"); \n        plt.title('L2 distance = %.3f' %(src_to_noisy_dist), fontsize=18)\nfig.savefig('figure_8.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# another image with more noise samples"},{"metadata":{"trusted":true},"cell_type":"code","source":"noise_levels_to_show = [0.2,0.5,0.75,1.0,1.25,1.5,2.0,2.5] # these are in units of std per feature\nnum_pertubation_per_noise_level = 6\n\nnum_cols = len(noise_levels_to_show)\nnum_rows = num_pertubation_per_noise_level\n\nselected_image_rep = ae_latent_rep[np.random.randint(ae_latent_rep.shape[0])][np.newaxis,:,:,:]\n\nfig = plt.figure(figsize=(30,24)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.93, bottom=0.02, hspace=0.15, wspace=0.05); \nplt.suptitle('figure 9: image + random noise in latent space (b)', fontsize=30)\nfor col, noise_level in enumerate(noise_levels_to_show):\n    for row in range(num_pertubation_per_noise_level):        \n        noise_to_add = noise_level * np.random.normal(scale=ae_latent_std, size=ae_latent_std.shape)\n        noisy_image_rep = selected_image_rep + noise_to_add\n        noisy_image = unnormalize_image(decoder.predict(noisy_image_rep)[0])\n        noisy_image[noisy_image > 255] = 255\n        noisy_image[noisy_image <   0] =   0\n        \n        src_to_noisy_dist = pdist(np.concatenate((selected_image_rep.reshape((1,-1)), noisy_image_rep.reshape((1,-1))), axis=0))[0]\n        plt.subplot(num_rows,num_cols,col+1+row*num_cols); plt.imshow(noisy_image.astype(np.uint8)); plt.axis(\"off\"); \n        plt.title('L2 distance = %.3f' %(src_to_noisy_dist), fontsize=16)\nfig.savefig('figure_9.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# another image with more noise samples"},{"metadata":{"trusted":true},"cell_type":"code","source":"noise_levels_to_show = [0.2,0.5,0.75,1.0,1.25,1.5,2.0,2.5] # these are in units of std per feature\nnum_pertubation_per_noise_level = 6\n\nnum_cols = len(noise_levels_to_show)\nnum_rows = num_pertubation_per_noise_level\n\nselected_image_rep = ae_latent_rep[np.random.randint(ae_latent_rep.shape[0])][np.newaxis,:,:,:]\n\nfig = plt.figure(figsize=(30,24)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.93, bottom=0.02, hspace=0.15, wspace=0.05); \nplt.suptitle('figure 10: image + random noise in latent space (c)', fontsize=30)\nfor col, noise_level in enumerate(noise_levels_to_show):\n    for row in range(num_pertubation_per_noise_level):        \n        noise_to_add = noise_level * np.random.normal(scale=ae_latent_std, size=ae_latent_std.shape)\n        noisy_image_rep = selected_image_rep + noise_to_add\n        noisy_image = unnormalize_image(decoder.predict(noisy_image_rep)[0])\n        noisy_image[noisy_image > 255] = 255\n        noisy_image[noisy_image <   0] =   0\n        \n        src_to_noisy_dist = pdist(np.concatenate((selected_image_rep.reshape((1,-1)), noisy_image_rep.reshape((1,-1))), axis=0))[0]\n        plt.subplot(num_rows,num_cols,col+1+row*num_cols); plt.imshow(noisy_image.astype(np.uint8)); plt.axis(\"off\"); \n        plt.title('L2 distance = %.3f' %(src_to_noisy_dist), fontsize=16)\nfig.savefig('figure_10.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Show latent interpulations between two images"},{"metadata":{"trusted":true},"cell_type":"code","source":"interpulation_weights = [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]\nnum_interpulations = len(interpulation_weights)\nnum_pairs = 5\n\nnum_cols = len(interpulation_weights) + 2\nnum_rows = num_pairs\n\nfig = plt.figure(figsize=(30,20)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.95, bottom=0.02, hspace=0.05, wspace=0.05); \nplt.suptitle('figure 11: latent space interpulations', fontsize=30)\nfor row in range(num_pairs):\n    # randomly select two pairs of images \n    selected_image_1_rep = ae_latent_rep[np.random.randint(ae_latent_rep.shape[0])][np.newaxis,:,:,:]\n    selected_image_2_rep = ae_latent_rep[np.random.randint(ae_latent_rep.shape[0])][np.newaxis,:,:,:]\n\n    # show left and right images\n    selected_image_1 = unnormalize_image(decoder.predict(selected_image_1_rep)[0])\n    selected_image_1[selected_image_1 > 255] = 255\n    selected_image_1[selected_image_1 <   0] =   0\n    \n    selected_image_2 = unnormalize_image(decoder.predict(selected_image_2_rep)[0])\n    selected_image_2[selected_image_2 > 255] = 255\n    selected_image_2[selected_image_2 <   0] =   0\n\n    image_ind = 1+row*num_cols\n    plt.subplot(num_rows,num_cols,1 + row*num_cols); plt.imshow(selected_image_1.astype(np.uint8)); plt.axis(\"off\"); \n    plt.subplot(num_rows,num_cols,(row+1)*num_cols); plt.imshow(selected_image_2.astype(np.uint8)); plt.axis(\"off\"); \n    \n    for col, weight in enumerate(interpulation_weights):\n\n        # create latent interpulations between them\n        interpulated_image_rep = (1-weight)*selected_image_1_rep + weight*selected_image_2_rep\n    \n        interpulated_image = unnormalize_image(decoder.predict(interpulated_image_rep)[0])\n        interpulated_image[interpulated_image > 255] = 255\n        interpulated_image[interpulated_image <   0] =   0\n        \n        plt.subplot(num_rows,num_cols,col+2+row*num_cols); plt.imshow(interpulated_image.astype(np.uint8)); plt.axis(\"off\"); \nfig.savefig('figure_11.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Show interpulations side by side along with linear blend of the pixels"},{"metadata":{"trusted":true},"cell_type":"code","source":"interpulation_weights = [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]\nnum_interpulations = len(interpulation_weights)\nnum_pairs = 10\n\nnum_cols = len(interpulation_weights) + 2\nnum_rows = 2*num_pairs\n\nfig = plt.figure(figsize=(30,50)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.96, bottom=0.02, hspace=0.06, wspace=0.05); \nplt.suptitle('figure 12: AE interpulations vs pixel space interpulations', fontsize=30)\nfor pair in range(num_pairs):\n    # randomly select two pairs of images \n    selected_image_1_rep = ae_latent_rep[np.random.randint(ae_latent_rep.shape[0])][np.newaxis,:,:,:]\n    selected_image_2_rep = ae_latent_rep[np.random.randint(ae_latent_rep.shape[0])][np.newaxis,:,:,:]\n\n    # show left and right images\n    selected_image_1 = unnormalize_image(decoder.predict(selected_image_1_rep)[0])\n    selected_image_1[selected_image_1 > 255] = 255\n    selected_image_1[selected_image_1 <   0] =   0\n    \n    selected_image_2 = unnormalize_image(decoder.predict(selected_image_2_rep)[0])\n    selected_image_2[selected_image_2 > 255] = 255\n    selected_image_2[selected_image_2 <   0] =   0\n\n    # plot odd rows (latent space interpulations)\n    row = 2*pair\n    plt.subplot(num_rows,num_cols,1 + row*num_cols); plt.imshow(selected_image_1.astype(np.uint8)); plt.axis(\"off\"); \n    plt.subplot(num_rows,num_cols,(row+1)*num_cols); plt.imshow(selected_image_2.astype(np.uint8)); plt.axis(\"off\"); \n\n    for col, weight in enumerate(interpulation_weights):\n        # create latent interpulations between them\n        interpulated_image_rep = (1-weight)*selected_image_1_rep + weight*selected_image_2_rep\n        interpulated_image = unnormalize_image(decoder.predict(interpulated_image_rep)[0])\n        interpulated_image[interpulated_image > 255] = 255\n        interpulated_image[interpulated_image <   0] =   0\n        plt.subplot(num_rows,num_cols,col+2+row*num_cols); plt.imshow(interpulated_image.astype(np.uint8)); plt.axis(\"off\"); \n        \n    # plot even rows (pixel space interpulations)\n    row = 2*pair + 1\n    plt.subplot(num_rows,num_cols,1 + row*num_cols); plt.imshow(selected_image_1.astype(np.uint8)); plt.axis(\"off\"); \n    plt.subplot(num_rows,num_cols,(row+1)*num_cols); plt.imshow(selected_image_2.astype(np.uint8)); plt.axis(\"off\"); \n\n    for col, weight in enumerate(interpulation_weights):\n        interpulated_image = (1-weight)*selected_image_1 + weight*selected_image_2      \n        plt.subplot(num_rows,num_cols,col+2+row*num_cols); plt.imshow(interpulated_image.astype(np.uint8)); plt.axis(\"off\"); \nfig.savefig('figure_12.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Apply Kmeans on the latents and show the cluster centers"},{"metadata":{},"cell_type":"markdown","source":"### first train a PCA to reduce dimentionality so that kmeans will ever finish running"},{"metadata":{"trusted":true},"cell_type":"code","source":"# create PCA model of the data for more efficient kmeans\nnum_components_for_kmeans = 384\n\ndog_PCA_for_kmeans = decomposition.PCA(n_components=num_components_for_kmeans, whiten=True)\ndog_PCA_for_kmeans.fit(ae_latent_rep_table)\n\nprint('total explained percent by %d components - %.1f%s' %(num_components_for_kmeans, 100*dog_PCA_for_kmeans.explained_variance_ratio_.sum(),'%'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num_clusters = 180\n\nX_for_kmeans = dog_PCA_for_kmeans.transform(ae_latent_rep_table)\n\ndog_Kmeans = cluster.KMeans(n_clusters=num_clusters)\n\nstart_time = time.time()\ncluster_inds = dog_Kmeans.fit_predict(X_for_kmeans)\nprint('finished training Kmeans model. took %.1f minutes' %((time.time()-start_time)/60))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# sort the clusters according to their frequency\ncluster_counter_dict = collections.Counter(cluster_inds)\nsorted_cluster_count = sorted(cluster_counter_dict.items(), key=operator.itemgetter(1))\nsorted_cluster_inds = [x[0] for x in sorted_cluster_count]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"num_rows = 20\nnum_cols = 9\n\nfig = plt.figure(figsize=(30,70)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.96, bottom=0.02, hspace=0.17, wspace=0.05); \nplt.suptitle('figure 13: Kmeans cluster centers', fontsize=30)\nfor k in range(num_clusters):\n    cluster_ind = sorted_cluster_inds[-k-1]\n    cluster_count = cluster_counter_dict[cluster_ind]\n    cluster_center_rep_row = dog_PCA_for_kmeans.inverse_transform(dog_Kmeans.cluster_centers_[cluster_ind,:][np.newaxis,:])\n    cluster_center_rep = np.reshape(cluster_center_rep_row, (1,4,4,encoder_output_channel_size))\n    cluster_doglike_image = unnormalize_image(decoder.predict(cluster_center_rep)[0])\n    cluster_doglike_image[cluster_doglike_image > 255] = 255\n    cluster_doglike_image[cluster_doglike_image <   0] =   0\n\n    plt.subplot(num_rows,num_cols,k+1); plt.imshow(cluster_doglike_image.astype(np.uint8)); \n    plt.title('(%d,%d)' %(k+1,cluster_count), fontsize=16); plt.axis(\"off\"); \nfig.savefig('figure_13.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Interpulate between two images from the same cluster"},{"metadata":{"trusted":true},"cell_type":"code","source":"interpulation_weights = [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]\nnum_interpulations = len(interpulation_weights)\nnum_pairs = 12\n\nmost_frequent_clusters_cutoff = min(60,num_clusters)\n\nnum_cols = len(interpulation_weights) + 2\nnum_rows = num_pairs\n\nfig = plt.figure(figsize=(30,42)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.95, bottom=0.02, hspace=0.05, wspace=0.05); \nplt.suptitle('figure 14: same cluster interpulations', fontsize=30)\nfor row in range(num_pairs):\n    \n    # randomly select a cluster amoungs the most frequent clusters\n    selected_cluster_ind = np.random.choice(sorted_cluster_inds[-most_frequent_clusters_cutoff:], size=1)\n    \n    # randomly select two pairs of images from the same cluster\n    possible_candidates = np.nonzero(cluster_inds == selected_cluster_ind)[0]\n    chosen_pair_inds = np.random.choice(possible_candidates, size=2, replace=False)\n    selected_image_1_rep = ae_latent_rep[chosen_pair_inds[0]][np.newaxis,:,:,:]\n    selected_image_2_rep = ae_latent_rep[chosen_pair_inds[1]][np.newaxis,:,:,:]\n\n    # show left and right images\n    selected_image_1 = unnormalize_image(decoder.predict(selected_image_1_rep)[0])\n    selected_image_1[selected_image_1 > 255] = 255\n    selected_image_1[selected_image_1 <   0] =   0\n    \n    selected_image_2 = unnormalize_image(decoder.predict(selected_image_2_rep)[0])\n    selected_image_2[selected_image_2 > 255] = 255\n    selected_image_2[selected_image_2 <   0] =   0\n\n    image_ind = 1+row*num_cols\n    plt.subplot(num_rows,num_cols,1 + row*num_cols); plt.imshow(selected_image_1.astype(np.uint8)); plt.axis(\"off\"); \n    plt.subplot(num_rows,num_cols,(row+1)*num_cols); plt.imshow(selected_image_2.astype(np.uint8)); plt.axis(\"off\"); \n    \n    for col, weight in enumerate(interpulation_weights):\n\n        # create latent interpulations between them\n        interpulated_image_rep = (1-weight)*selected_image_1_rep + weight*selected_image_2_rep\n    \n        interpulated_image = unnormalize_image(decoder.predict(interpulated_image_rep)[0])\n        interpulated_image[interpulated_image > 255] = 255\n        interpulated_image[interpulated_image <   0] =   0\n        \n        plt.subplot(num_rows,num_cols,col+2+row*num_cols); plt.imshow(interpulated_image.astype(np.uint8)); plt.axis(\"off\"); \nfig.savefig('figure_14.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Train PCA on encoder represnetation to sample from"},{"metadata":{"trusted":true},"cell_type":"code","source":"# create PCA model of the data\nnum_components = 512\n\ndog_PCA = decomposition.PCA(n_components=num_components, whiten=True)\ndog_PCA.fit(X_rep_autoencoder)\n\nprint('finished training PCA model')\nX_pca = dog_PCA.transform(X_rep_autoencoder)\n\n# show cumulative variance explained\nfig = plt.figure(figsize=(16,10)); plt.title('figure 15: PCA variance explained', fontsize=26)\nplt.plot(100*np.concatenate((np.array([0]),np.cumsum(dog_PCA.explained_variance_ratio_))))\nplt.xlabel('num components', fontsize=16); plt.ylabel('% variance explained', fontsize=16); plt.ylim(-1,101); plt.xlim(-1,num_components+1);\nfig.savefig('figure_15.png')\nprint('total explained percent by %d components - %.1f%s' %(num_components, 100*dog_PCA.explained_variance_ratio_.sum(),'%'))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Show some Autoencoder + PCA reconstructions (maximal expected performace)"},{"metadata":{"trusted":true},"cell_type":"code","source":"num_rows = 6\nnum_cols = 9\nnum_images_to_show = num_rows*num_cols\n\nX_rep_autoencoder_rec = dog_PCA.inverse_transform(X_pca)\nselected_inds = np.random.choice(X_rep_autoencoder_rec.shape[0], size=num_images_to_show, replace=False)\n\nfig = plt.figure(figsize=(30,20)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.95, bottom=0.02, hspace=0.05, wspace=0.05); \nplt.suptitle('figure 16: AE + PCA reconstrctions', fontsize=30)\nfor k, selected_ind in enumerate(selected_inds):\n    decoder_input = np.reshape(X_rep_autoencoder_rec[selected_ind,:], (1,4,4,encoder_output_channel_size))\n    doglike_image = unnormalize_image(decoder.predict(decoder_input)[0])\n    doglike_image[doglike_image > 255] = 255\n    doglike_image[doglike_image <   0] =   0\n\n    plt.subplot(num_rows,num_cols,k+1); plt.imshow(doglike_image.astype(np.uint8)); plt.axis(\"off\")\nfig.savefig('figure_16.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Show some PCA unit activation histograms (should appear approx. gaussian)"},{"metadata":{"trusted":true},"cell_type":"code","source":"selected_inds = np.random.choice(num_components, size=36, replace=False)\n\nfig = plt.figure(figsize=(30,12)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.93, bottom=0.02, hspace=0.13, wspace=0.13); \nplt.suptitle('figure 17: PCA unit activations', fontsize=30)\nfor k, selected_ind in enumerate(selected_inds):\n    unit_activations = X_pca[:,selected_ind]\n    range_limit = max(abs(unit_activations.min()), abs(unit_activations.max()))\n    activation_range = np.linspace(-range_limit,range_limit,100)\n    plt.subplot(4,9,k+1); plt.hist(unit_activations, bins=activation_range, log=True);\nfig.savefig('figure_17.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Decompose data with ICA (optional)\nthe goal is to find statistically independent directions in the latent space for more efficient sampling"},{"metadata":{"trusted":true},"cell_type":"code","source":"#decomposition_method = 'ICA'\ndecomposition_method = 'PCA'\n\nif decomposition_method == 'ICA':\n    num_ICA_components = num_components\n    dog_ICA = decomposition.FastICA(n_components=num_ICA_components, algorithm='parallel', whiten=True)\n\n    start_time = time.time()\n    X_pca_ica = dog_ICA.fit_transform(X_pca)\n    print('finished training ICA model. took %.1f minutes' %((time.time()-start_time)/60))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Show unit activation histograms for ICA model (should be non gaussians)"},{"metadata":{"trusted":true},"cell_type":"code","source":"if decomposition_method == 'ICA':\n    selected_inds = np.random.choice(num_ICA_components, size=36, replace=False)\n\n    fig = plt.figure(figsize=(30,12)); \n    plt.subplots_adjust(left=0.02, right=0.98, top=0.93, bottom=0.02, hspace=0.13, wspace=0.13); \n    plt.suptitle('figure 18: ICA unit activations', fontsize=30)\n    for k, selected_ind in enumerate(selected_inds):\n        unit_activations = X_pca_ica[:,selected_ind]\n        range_limit = max(abs(unit_activations.min()), abs(unit_activations.max()))\n        activation_range = np.linspace(-range_limit,range_limit,100)\n        plt.subplot(4,9,k+1); plt.hist(unit_activations, bins=activation_range, log=True);\n    fig.savefig('figure_18.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Sample from single gaussian in the selected decomposition methond and present samples"},{"metadata":{"trusted":true},"cell_type":"code","source":"if decomposition_method == 'ICA':\n    covariance_matrix_regularization = 1e-5\nelse:\n    covariance_matrix_regularization = 1e-3\n\n# generate several random samples from the gaussian model and present them\ndog_single_gaussian_model = mixture.GaussianMixture(n_components=1, covariance_type='diag', reg_covar=covariance_matrix_regularization, n_init=5)\n\nif decomposition_method == 'ICA':\n    dog_single_gaussian_model.fit(X_pca_ica)\nelse:\n    dog_single_gaussian_model.fit(X_pca)\n\nnum_rows = 6\nnum_cols = 9\nnum_images_to_show = num_rows*num_cols\n\nrandom_latents = dog_single_gaussian_model.sample(num_images_to_show)[0]\nprint(random_latents.shape, random_latents.mean(), random_latents.std())\n\nif decomposition_method == 'ICA':\n    random_doglike_vectors = dog_PCA.inverse_transform(dog_ICA.inverse_transform(random_latents))\nelse:\n    random_doglike_vectors = dog_PCA.inverse_transform(random_latents)\nprint(random_doglike_vectors.shape, random_doglike_vectors.mean(), random_doglike_vectors.std())\n\nfig = plt.figure(figsize=(30,20)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.95, bottom=0.02, hspace=0.05, wspace=0.05); \nplt.suptitle('figure 19: single gaussian samples', fontsize=30)\nfor k in range(num_images_to_show):\n    # convert to \n    decoder_input = np.reshape(random_doglike_vectors[k,:], (1,4,4,encoder_output_channel_size))\n    doglike_image = unnormalize_image(decoder.predict(decoder_input)[0])\n    doglike_image[doglike_image > 255] = 255\n    doglike_image[doglike_image <   0] =   0\n\n    plt.subplot(num_rows,num_cols,k+1); plt.imshow(doglike_image.astype(np.uint8)); plt.axis(\"off\")\nfig.savefig('figure_19.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Find best number gaussians for GMM"},{"metadata":{"trusted":true},"cell_type":"code","source":"num_gaussians_to_try = [1,2,3,4,5,7,10,13,17,22,30,40,50,60,90,120,180,300]\ncovariance_matrix_type = 'diag'\n\nif decomposition_method == 'ICA':\n    covariance_matrix_regularization = 1e-5\nelse:\n    covariance_matrix_regularization = 1e-3\n\nvalid_fraction = 0.3\nvalid_cutoff = int((1-valid_fraction)*X_pca.shape[0])\nrand_perm = np.random.permutation(X_pca.shape[0])\n\nX_pca_train = X_pca[rand_perm[:valid_cutoff],:]\nX_pca_valid = X_pca[rand_perm[valid_cutoff:],:]\nif decomposition_method == 'ICA':\n    X_pca_ica_train = X_pca_ica[rand_perm[:valid_cutoff],:]\n    X_pca_ica_valid = X_pca_ica[rand_perm[valid_cutoff:],:]\n\ntrain_LogLikelihood = []\nvalid_LogLikelihood = []\nfor num_gaussians in num_gaussians_to_try:\n    \n    curr_dog_GMM = mixture.GaussianMixture(n_components=num_gaussians, covariance_type=covariance_matrix_type, n_init=2,\n                                           reg_covar=covariance_matrix_regularization, verbose=0, verbose_interval=1)\n    \n    if decomposition_method == 'ICA':\n        curr_dog_GMM.fit(X_pca_ica_train)\n        train_LL = curr_dog_GMM.score_samples(X_pca_ica_train).mean()\n        valid_LL = curr_dog_GMM.score_samples(X_pca_ica_valid).mean()\n    else:\n        curr_dog_GMM.fit(X_pca_train)\n        train_LL = curr_dog_GMM.score_samples(X_pca_train).mean()\n        valid_LL = curr_dog_GMM.score_samples(X_pca_valid).mean()\n    \n    print('for %d gaussians: (train,valid) LogLikelihood = (%.5f,%.5f)' %(num_gaussians, train_LL, valid_LL))\n    \n    train_LogLikelihood.append(train_LL)\n    valid_LogLikelihood.append(valid_LL)\n\nfig = plt.figure(figsize=(20,10)); plt.title('figure 20: GMM LL vs number of gaussians', fontsize=26)\nplt.plot(num_gaussians_to_try, train_LogLikelihood, color='b')\nplt.plot(num_gaussians_to_try, valid_LogLikelihood, color='g')\nplt.legend(['train','valid'], fontsize=16)\nplt.ylabel('Log Likelihood', fontsize=16); plt.xlabel('num gaussians', fontsize=16)\nfig.savefig('figure_20.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Train final GMM"},{"metadata":{"trusted":true},"cell_type":"code","source":"# train a mixture of gaussians model in the PCA space\nnum_gaussians = num_gaussians_to_try[np.argmax(valid_LogLikelihood)]\n\nprint('selected number of gaussians is %d' %(num_gaussians))\n\ncovariance_matrix_type = 'diag'\nif decomposition_method == 'ICA':\n    covariance_matrix_regularization = 1e-5\nelse:\n    covariance_matrix_regularization = 1e-3\ndog_gaussian_mixture_model = mixture.GaussianMixture(n_components=num_gaussians, covariance_type=covariance_matrix_type, n_init=3, \n                                                     reg_covar=covariance_matrix_regularization, verbose=2, verbose_interval=1)\n\nif decomposition_method == 'ICA':\n    dog_gaussian_mixture_model.fit(X_pca_ica_train)\nelse:\n    dog_gaussian_mixture_model.fit(X_pca_train)\n\nprint('finished training GMM')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Show GMM samples"},{"metadata":{"trusted":true},"cell_type":"code","source":"# generate several random samples from the mixture model and present them\nnum_rows = 6\nnum_cols = 9\nnum_images_to_show = num_rows*num_cols\n\nrandom_latents = dog_gaussian_mixture_model.sample(num_images_to_show)[0]\n\nif decomposition_method == 'ICA':\n    random_doglike_vectors = dog_PCA.inverse_transform(dog_ICA.inverse_transform(random_latents))\nelse:\n    random_doglike_vectors = dog_PCA.inverse_transform(random_latents)\n\nfig = plt.figure(figsize=(30,20)); \nplt.subplots_adjust(left=0.02, right=0.98, top=0.95, bottom=0.02, hspace=0.05, wspace=0.05); \nplt.suptitle('figure 21: GMM samples', fontsize=30)\nfor k in range(num_images_to_show):\n    decoder_input = np.reshape(random_doglike_vectors[k,:], (1,4,4,encoder_output_channel_size))\n    doglike_image = unnormalize_image(decoder.predict(decoder_input)[0])\n    doglike_image[doglike_image > 255] = 255\n    doglike_image[doglike_image <   0] =   0\n\n    plt.subplot(num_rows,num_cols,k+1); plt.imshow(doglike_image.astype(np.uint8)); plt.axis(\"off\")\nfig.savefig('figure_21.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Create a submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"num_images_to_submit = 10000\n\nsampled_latents = dog_gaussian_mixture_model.sample(num_images_to_submit)[0]\n\nif decomposition_method == 'ICA':\n    random_doglike_vectors = dog_PCA.inverse_transform(dog_ICA.inverse_transform(sampled_latents))\nelse:\n    random_doglike_vectors = dog_PCA.inverse_transform(sampled_latents)\n\nz = zipfile.PyZipFile('images.zip', mode='w')\nfor k in range(num_images_to_submit):\n    decoder_input = np.reshape(random_doglike_vectors[k,:], (1,4,4,encoder_output_channel_size))\n    doglike_image = unnormalize_image(decoder.predict(decoder_input)[0])\n    doglike_image[doglike_image > 255] = 255\n    doglike_image[doglike_image <   0] =   0\n    image_to_save = Image.fromarray(doglike_image.astype(np.uint8))\n\n    image_filename = '%d.png' %(k)\n    image_to_save.save(image_filename,'PNG'); z.write(image_filename); os.remove(image_filename)\nprint('finished writing \"image.zip\"')","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}