<jupyter_start><jupyter_text>Visualise the data<jupyter_code>with open('data/label_dict.pkl', 'rb') as f:
    label_dict = pkl.load(f)
data = tf.data.TFRecordDataset(data_dir).map(parse_example, num_parallel_calls=tf.data.experimental.AUTOTUNE)

data_viz = iter(data.batch(1))
for i in range(50):
    img, label = next(data_viz)
plt.imshow(img.numpy()[0]) #Add 0.5 if it's zero-centered. (it's not quite 0.5 strictly but it's close enough)
print("Label: {}".format(label_dict[label.numpy()[0]]))
print(img.numpy()[0].mean())<jupyter_output>Label: denseresidential
0.518327<jupyter_text>Build the modelLi 2016 - Li 2020 - 1e6 patches from 1e5 images  Build the custom layer in keras<jupyter_code># class CustomKernelInit(tf.keras.initializers.Initializer):
#     def __init__(self, visual_words, w):
#         K = visual_words.shape[0]
#         words_kernel = -2*visual_words.reshape((K, w, w, -1)) # ||x-ck||^2 = x t x - 2ck t x + ck t ck 
#         self.full_kernel = np.append(words_kernel, np.ones(words_kernel.shape), axis=3)
#         self.full_kernel = np.moveaxis(self.full_kernel, 0, -1) #(input rows, input columns, input channels, output channels)
        
#     def __call__(self, shape, dtype=None):
#         return self.full_kernel
def kernel_init(visual_words, w):
        K = visual_words.shape[0]
        words_kernel = -2*visual_words.reshape((K, w, w, -1)) # ||x-ck||^2 = x t x - 2ck t x + ck t ck 
        full_kernel = np.append(words_kernel, np.ones(words_kernel.shape), axis=3)
        full_kernel = np.moveaxis(full_kernel, 0, -1) #(input rows, input columns, input channels, output channels)
        
        return full_kernel
class FeatureMapLayer(tf.keras.layers.Layer):
    def __init__(self, K, w, visual_words):
        super(FeatureMapLayer, self).__init__()
        self.visual_words = visual_words
        self.K = K
        self.w = w
        
        
        
    def build(self, input_shape):
        kernel_initializer = kernel_init(self.visual_words, self.w)
        augmented_shape = (input_shape[0], input_shape[1], input_shape[2], input_shape[3]*2)
        
        self.conv2 = tf.keras.layers.Conv2D(filters=self.K, kernel_size=self.w, use_bias=False,
                                           input_shape=augmented_shape[1:])
        
        self.conv2.build(augmented_shape) 
        self.conv2.set_weights([kernel_init(self.visual_words, self.w)])        

        
    def call(self, input):
        #Augment N image layers w/ N additional layers that are element-wise squared
        input_augmented = tf.concat([input, tf.math.square(input)], axis = 3) 
        #print(input_augmented.shape)
        conv_out = self.conv2(input_augmented) #Returns x t x - 2ck t x
        
        square_word_norms = np.linalg.norm(self.visual_words, axis=1)**2 #Squared norms of the visual words - ck t ck
        
        z2 = conv_out + square_word_norms #Note to self - this might break if K is equal to one of the image dimensions
        z = tf.math.sqrt(z2)
        
        u = tf.math.reduce_mean(z, axis=3) #Take mean over the axis of features (which is of length K)
        u = tf.expand_dims(u, 3)
        u = tf.tile(u, [1,1,1, self.K])
        
        f = tf.math.maximum(0.0, u-z)
        
        return f
def zca_whitening_matrix(X):
    """
    Function to compute ZCA whitening matrix (aka Mahalanobis whitening).
    INPUT:  X: [M x N] matrix.
        Rows: Variables
        Columns: Observations
    OUTPUT: ZCAMatrix: [M x M] matrix
    """
    # Covariance matrix [column-wise variables]: Sigma = (X-mu)' * (X-mu) / N
    sigma = np.cov(X, rowvar=True) # [M x M]
    # Singular Value Decomposition. X = U * np.diag(S) * V
    U,S,V = np.linalg.svd(sigma)
        # U: [M x M] eigenvectors of sigma.
        # S: [M x 1] eigenvalues of sigma.
        # V: [M x M] transpose of U
    # Whitening constant: prevents division by zero
    epsilon = 1e-5
    # ZCA Whitening matrix: U * Lambda * U'
    ZCAMatrix = np.dot(U, np.dot(np.diag(1.0/np.sqrt(S + epsilon)), U.T)) # [M x M]
    return ZCAMatrix
def open_image(image_loc, expected_dim):
    pixel_data = tifffile.imread(image_loc)
    
    height = pixel_data.shape[0]
    width = pixel_data.shape[1]
    channels = pixel_data.shape[2]

    if width != height :
        print('Image {} is not square, returning empty variable'.format(image_loc))
        return None
    
    if width != expected_dim:
        print('Image {} has dimension {}, expected {}, returning empty variable'.format(image_loc, width, expected_dim))
        return None
    
    return pixel_data / 10000 #Normalised

def load_to_numpy(data_dir):
    data = []
    for image_loc in os.listdir(data_dir):
        data.append(open_image(data_dir + image_loc, 64))
    return np.stack(data)<jupyter_output><empty_output><jupyter_text>Need data as numpy array:<jupyter_code>data_source = 'sentinel' #{'sentinel', UCM}
if data_source == 'UCM':
    data_numpy = np.stack([x[0] for x in list(data.as_numpy_iterator())]) #Could change this to not load whole dataset into memory if too big
elif data_source == 'sentinel':
    train_images_dir = 'data/sentinel_unlabelled/'
    data_numpy = load_to_numpy(train_images_dir)<jupyter_output><empty_output><jupyter_text>Define hyperparameters:<jupyter_code>num_images = data_numpy.shape[0]
image_size = data_numpy.shape[1]
no_channels = data_numpy.shape[3]
seed = 0

datagen_args = dict(
    featurewise_center=True, samplewise_center=False,
    featurewise_std_normalization=False, samplewise_std_normalization=False,
    zca_whitening=True, zca_epsilon=0.1, rotation_range=0, width_shift_range=0.0,
    height_shift_range=0.0, brightness_range=None, shear_range=0.0, zoom_range=0.0,
    channel_shift_range=0.0, fill_mode='nearest', cval=0.0,
    horizontal_flip=False, vertical_flip=False, rescale=None,
    preprocessing_function=None, data_format=None, validation_split=0.0, dtype=None
    )

K = [100, 900]
featureset_sizes = [411, 3600]
w = 6 #Size of receptive field for feature mapping
s = 2 #Size of local max pooling window

normalise_feature_set = True<jupyter_output><empty_output><jupyter_text>Unsupervised training for first layer (K Means)<jupyter_code>%%time
#Unsupervised training for first layer

feature_sets = []
visual_dicts = []

#Randomly sample image patches - for the input layer
image_indices = choices(range(num_images), k = featureset_sizes[0])
row_indices = choices(range(image_size - w + 1), k = featureset_sizes[0])
col_indices = choices(range(image_size - w + 1), k = featureset_sizes[0])

feature_set = []

for i in range(len(image_indices)):
    ii, ri, ci = image_indices[i], row_indices[i], col_indices[i]
    sample_patch = data_numpy[ii][ri:ri+w,ci:ci+w]
    feature_set.append(sample_patch.flatten())
    
feature_set = np.array(feature_set)

if normalise_feature_set:
    feature_set = feature_set - feature_set.mean(axis=0)

z = zca_whitening_matrix(feature_set.T)
feature_set_whitened = np.dot(z, feature_set.T).T
feature_sets.append(feature_set_whitened)

kmeans = KMeans(n_clusters=K[0], random_state = seed).fit(feature_sets[0])    

visual_dicts.append(kmeans.cluster_centers_)<jupyter_output>CPU times: user 2.46 s, sys: 6.97 s, total: 9.42 s
Wall time: 422 ms<jupyter_text>Build model for first layer:<jupyter_code>inputs = tf.keras.Input(shape=(image_size,image_size,no_channels))

feature_map_0 = FeatureMapLayer(K[0], w, np.array(visual_dicts[0]))
f0 = feature_map_0(inputs)

local_pool_0 = tf.keras.layers.MaxPool2D(pool_size=(2,2))
p0 = local_pool_0(f0)

global_pool_size_0 = (p0.shape[1]//2, p0.shape[2]//2)
global_pool_strides_0 = (p0.shape[1]//4, p0.shape[2]//4)
global_pool_0 = tf.keras.layers.AveragePooling2D(pool_size=global_pool_size_0, strides=global_pool_strides_0)
g0 = global_pool_0(p0)

flat_0 = tf.keras.layers.Flatten()
out_0 = flat_0(g0)


model_firstlayer = tf.keras.Model(inputs=inputs, outputs=[out_0, p0], name="MultilayerFeatureLearning")<jupyter_output><empty_output><jupyter_text>Unsupervised training for second layer<jupyter_code>%%time
featuremap_0_size = model_firstlayer.predict(data_numpy[0:1])[1].shape[1]

#Randomly sample image patches - for the input layer
image_indices = choices(range(num_images), k = featureset_sizes[1])
row_indices = choices(range(featuremap_0_size - w + 1), k = featureset_sizes[1])
col_indices = choices(range(featuremap_0_size - w + 1), k = featureset_sizes[1])

feature_set = []

for i in range(len(image_indices)):
    ii, ri, ci = image_indices[i], row_indices[i], col_indices[i]
    first_layer_image_out = model_firstlayer.predict(data_numpy[ii:ii+1])
    sample_patch = first_layer_image_out[1][0][ri:ri+w,ci:ci+w]
    feature_set.append(sample_patch.flatten())
    
feature_set = np.array(feature_set)

if normalise_feature_set:
    feature_set = feature_set - feature_set.mean(axis=0)
    
z = zca_whitening_matrix(feature_set.T)
feature_set_whitened = np.dot(z, feature_set.T).T


if len(feature_sets) == 1:
    feature_sets.append(feature_set_whitened)
else:
    feature_sets[1] = feature_set_whitened
%%time
kmeans = KMeans(n_clusters=K[1], random_state = seed).fit(feature_sets[1])    

if len(visual_dicts) == 1:
    visual_dicts.append(kmeans.cluster_centers_)
else:
    visual_dicts[1] = kmeans.cluster_centers_
    
save_file('visual_dicts.pkl', visual_dicts)<jupyter_output>CPU times: user 6min 23s, sys: 10min 27s, total: 16min 51s
Wall time: 48.3 s<jupyter_text>Build rest of model<jupyter_code>feature_map_1 = FeatureMapLayer(K[1], w, np.array(visual_dicts[1]))
f1 = feature_map_1(p0)

local_pool_1 = tf.keras.layers.MaxPool2D(pool_size=(2,2))
p1 = local_pool_1(f1)

global_pool_size_1 = (p1.shape[1]//2, p1.shape[2]//2)
global_pool_strides_1 = (p1.shape[1]//4, p1.shape[2]//4)
global_pool_1 = tf.keras.layers.AveragePooling2D(pool_size=global_pool_size_1, strides=global_pool_strides_1)
g1 = global_pool_1(p1)

flat_1 = tf.keras.layers.Flatten()
out_1 = flat_1(g1)

concat = tf.keras.layers.Concatenate(axis=1)
output = concat([out_0, out_1])

model_twolayers = tf.keras.Model(inputs=inputs, outputs = output)
model_twolayers(data_numpy[0:1])
timestr = time.strftime("%Y%m%d-%H%M%S")

model_twolayers.save("{}+UnsupervisedFeatureLearning{}+w{}+K{}K{}+M{}M{}_ZCANorm{}_Seed{}".format(data_source, timestr, w, K[0], K[1], featureset_sizes[0], featureset_sizes[1], normalise_feature_set, seed))<jupyter_output>INFO:tensorflow:Assets written to: sentinel+UnsupervisedFeatureLearning20210317-162239+w3+K100K900+M411M3600_ZCANormTrue_Seed0/assets
{'fraction_non_alphanumeric': 0.12642911585365854, 'fraction_numerical': 0.027820121951219513, 'mean_word_length': 4.261654135338346, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '607547', 'license': 'RSA-MD', 'n_tokens_mistral': 4061, 'n_tokens_neox': 3879, 'n_words': 767, 'path': 'settlement_segmentation/liunsupervised/deep_feature_learning.ipynb', 'repo_name': 'ai4er-cdt/gtc-exposure', 'size': '10496'}