<jupyter_start><jupyter_text>Build a variational autoencoderBased on the ability to generate UK car registration numbers, one can build a dataset of training and test data. With a suitable vectorizer, one can go ahead and make an autoencoder.The decoder part of the variational encoder can then be used to generate new reg numbers. Referenceshttps://blog.keras.io/building-autoencoders-in-keras.htmlhttps://towardsdatascience.com/understanding-variational-autoencoders-vaes-f70510919f73<jupyter_code>import sys
import math
import numpy as np
import keras
from keras import regularizers
from keras import layers
from keras import backend
from keras import optimizers
import random

if '..' not in sys.path:
    sys.path.append('..')
from car_reg_generator.car_reg_generator.uk_reg import UkRegGenerator
from car_reg_generator.car_reg_generator.uk_reg import UkRegDvlaVectorizer<jupyter_output><empty_output><jupyter_text>Data generation<jupyter_code>n_train = 10000
n_val = 1000
n_test = 1000

random.seed(0)
g = UkRegGenerator()
v = UkRegDvlaVectorizer()

train_strs = [g.get_reg() for _ in range(n_train)]
train_vecs = np.array([v.vectorize(x) for x in train_strs])
val_strs = [g.get_reg() for _ in range(n_val)]
val_strs = [x for x in val_strs if x not in set(train_strs)]  # sanity check
val_vecs = np.array([v.vectorize(x) for x in val_strs])
test_strs = [g.get_reg() for _ in range(n_test)]
test_strs = [x for x in test_strs if x not in set(train_strs) | set(val_strs)]  # sanity check
test_vecs = np.array([v.vectorize(x) for x in test_strs])

vec_length = len(train_vecs[0])
print('input vector length = ' + str(vec_length))
print('actual number of validation samples = ' + str(len(val_strs)))
print('actual number of test samples = ' + str(len(test_strs)))<jupyter_output>input vector length = 150
actual number of validation samples = 1000
actual number of test samples = 1000<jupyter_text>Variational autoencoderDesign the architecture<jupyter_code>latent_dim = 15  # 15 works better for a good model accuracy, 
                # but lower dimensions work better for directly using the latent space

### Define the encoder
inputs = keras.Input(shape=(vec_length,))
h = inputs
h = layers.Dense(vec_length, activation='relu')(h)
h = layers.Dense(50, activation='relu')(h)
h = layers.Dense(25, activation='relu')(h)
z_mean = layers.Dense(latent_dim)(h)
z_log_sigma = layers.Dense(latent_dim)(h)

def sampling(args):
    z_mean, z_log_sigma = args
    epsilon = backend.random_normal(shape=(backend.shape(z_mean)[0], latent_dim),
                                    mean=0., stddev=0.1)
    return z_mean + backend.exp(z_log_sigma) * epsilon

z = layers.Lambda(sampling)([z_mean, z_log_sigma])

encoder = keras.Model(inputs, [z_mean, z_log_sigma, z], name='encoder')

### Define the decoder
latent_inputs = keras.Input(shape=(latent_dim,), name='z_sampling')
x = latent_inputs
x = layers.Dense(25, activation='relu')(x)
x = layers.Dense(50, activation='relu')(x)
x = layers.Dense(vec_length, activation='relu')(x)
x = layers.Dense(vec_length, activation='sigmoid')(x)
outputs = x
decoder = keras.Model(latent_inputs, outputs, name='decoder')

outputs = decoder(encoder(inputs)[2])
vae = keras.Model(inputs, outputs, name='vae_mlp')

### Define the loss function
reconstruction_loss = keras.losses.binary_crossentropy(inputs, outputs)
reconstruction_loss *= vec_length
kl_loss = 1 + z_log_sigma - backend.square(z_mean) - backend.exp(z_log_sigma)
kl_loss = backend.sum(kl_loss, axis=-1)
kl_loss *= -0.5
vae_loss = backend.mean(reconstruction_loss + kl_loss)
vae.add_loss(vae_loss)<jupyter_output><empty_output><jupyter_text>Do some training. The training below has a hard time optimizing when the latent dimension is less than 15, but for some of the techniques used later in the notebook, a much smaller latent dimension is required.<jupyter_code>opt = optimizers.Adamax(learning_rate=0.004)  # larger learning rate for faster descent
vae.compile(optimizer=opt, loss='binary_crossentropy')
vae.fit(train_vecs, train_vecs,
        epochs=200,
        batch_size=8,
        shuffle=True,
        validation_data=(val_vecs, val_vecs))

opt = optimizers.Adamax(learning_rate=0.001)  # slower learning rate for approaching local min
vae.compile(optimizer=opt, loss='binary_crossentropy')
vae.fit(train_vecs, train_vecs,
        epochs=200,
        batch_size=8,
        shuffle=True,
        validation_data=(val_vecs, val_vecs))<jupyter_output>Epoch 1/200
1250/1250 [==============================] - 1s 945us/step - loss: 28.4330 - val_loss: 26.7790
Epoch 2/200
1250/1250 [==============================] - 1s 824us/step - loss: 25.5594 - val_loss: 24.5344
Epoch 3/200
1250/1250 [==============================] - 1s 887us/step - loss: 23.4502 - val_loss: 22.9401
Epoch 4/200
1250/1250 [==============================] - 1s 913us/step - loss: 22.2339 - val_loss: 22.0361
Epoch 5/200
1250/1250 [==============================] - 1s 976us/step - loss: 21.5440 - val_loss: 21.5660
Epoch 6/200
1250/1250 [==============================] - 1s 847us/step - loss: 21.0753 - val_loss: 21.1477
Epoch 7/200
1250/1250 [==============================] - 1s 851us/step - loss: 20.6930 - val_loss: 20.8693
Epoch 8/200
1250/1250 [==============================] - 1s 858us/step - loss: 20.3107 - val_loss: 20.6205
Epoch 9/200
1250/1250 [==============================] - 1s 879us/step - loss: 19.8597 - val_loss: 20.1434
Epoch 10/200
1250/1250 [=============[...]<jupyter_text>Evaluation<jupyter_code>decoded_regs = vae.predict(test_vecs)
print(test_strs[:10])
recovered = [v.recover(x) for x in decoded_regs]
print(recovered[:10])
acc = np.sum([x == y for x, y in zip(recovered, test_strs)]) / len(test_strs)
print('overall accuracy = ' + str(acc))
print('by digit accuracy:')
for i_digit in range(len(test_strs[0])):
    acc = np.sum([x[i_digit] == y[i_digit] for x, y in zip(recovered, test_strs)]) / len(test_strs)
    print('    ' + str(acc))<jupyter_output>['MY33WSZ', 'WE83GXF', 'LU19JLV', 'KX41GVE', 'WP91NCV', 'SG47LFF', 'HB17WIY', 'CA89AXN', 'GT57AKA', 'CD77SJT']
['MY33WSZ', 'WE83GXF', 'LU19JLV', 'KX41GVE', 'WP91NCX', 'SG47KFF', 'HB17WIY', 'CA89OXN', 'GT57YKA', 'CD77SJT']
overall accuracy = 0.517
by digit accuracy:
    0.991
    0.967
    0.992
    0.993
    0.688
    0.903
    0.882<jupyter_text>Use decoder as generatorBy directly sampling in the latent space, and running through the decoder, one effectively has a data generator.<jupyter_code>[v.recover(x) for x in decoder.predict(np.random.random((5, latent_dim)))]<jupyter_output><empty_output><jupyter_text>One can do some sampling and look at some stats to see what kind of data comes out. For instance, the number of reg numbers with pair of integers below 50 and above 50<jupyter_code>random_vectors = np.random.random((1000, latent_dim)) - 0.5
generated_regs = [v.recover(x) for x in decoder.predict(random_vectors)]
n_reg_below_50 = len([x for x in generated_regs if int(x[2:4]) < 50])
n_reg_above_50 = len([x for x in generated_regs if int(x[2:4]) >= 50])
print('reg numbers below 50 = ' + str(n_reg_below_50))
print('reg numbers above 50 = ' + str(n_reg_above_50))<jupyter_output>reg numbers below 50 = 442
reg numbers above 50 = 558<jupyter_text>Avoid generating any training dataThere is a chance that a random generation of data might produce one of the actual registration numbers used in training. However, any input registrations are represented as a probability distribution in the latent space; so one could require that any sample vector used in the latent space to be 'sufficiently far' down the tails of all the recorded distributions.The distributions of the training data in the latent space are readily available:<jupyter_code># Get the distributions for each reg number in the data
all_strs = train_strs + val_strs + test_strs
all_vecs = np.array([v.vectorize(x) for x in all_strs])
means, log_sigmas, _ = encoder.predict(all_vecs)
distn_dict = {x: (y,z) for x, y, z in zip(all_strs, means, np.exp(log_sigmas))}  # dict<reg no: (mean, sigma)>
example_reg = list(distn_dict.keys())[0]
print('distribution for input data registration ' + str(example_reg))
print('    mu = ' + str(distn_dict[example_reg][0]))
print('    diagonal sigmas = ' + str(distn_dict[example_reg][1]))
# for a multivariate distribution ~ N(mean, Var) valued in the latent space, 
# mean is the first vector below, and Var is a diagonal matrix, the square root
# of the diagonal entries being given in the second vector.
# If the KL loss function pulls the distributions to be close to N(0,1)<jupyter_output>distribution for input data registration YK66BIQ
    mu = [ 6.05341792e-01  1.34610564e-01  4.82149050e-03  1.87374011e-01
  6.64951280e-04 -2.46015787e-02 -1.18131995e-01  8.47651958e-01
  1.24904223e-01 -4.60249305e-01  4.73329067e-01 -5.29569853e-03
  1.75410032e-01  7.85022378e-01 -9.21758711e-02]
    diagonal sigmas = [0.5313385  0.44662854 0.98994035 0.6331824  0.994979   0.3902884
 0.9549835  0.4150101  0.39563507 0.33123204 0.4422034  0.99358207
 0.48427805 0.4332095  0.41187736]<jupyter_text>One can work out a rough area in the latent space that the data vectors occupy. This could be important; certainly any point in the latent space can be used (since there is a distribution over the space), but how valid are samples from far outside the region covered by the 'heavy' parts of the obtained distributions?<jupyter_code>all_means_mtr = np.array([x[0] for x in distn_dict.values()])
print('minimum of means over each latent dimension:')
print(np.min(all_means_mtr, axis=0))
print('maximum of means over each latent dimension:')
print(np.max(all_means_mtr, axis=0))<jupyter_output>minimum of means over each latent dimension:
[-1.690712   -1.7782323  -0.03395703 -1.4910543  -0.03753981 -1.7040544
 -0.5108407  -1.7530599  -1.7487404  -1.5404701  -1.1142322  -0.0425394
 -1.4865733  -1.6822727  -1.7714614 ]
maximum of means over each latent dimension:
[1.6918695  1.9303555  0.02382918 1.2332906  0.01244401 1.698646
 0.74276054 1.8050513  1.8302792  1.6911187  1.1751331  0.01992849
 1.6836305  1.8916475  1.789688  ]<jupyter_text>Generating subsets of data with certain characteristicsAs a more general example of 'avoid generating any training data', consider the case of 'avoid generating data with certain characteristics'. Since the latent space has a joint probability density function over it, this can be answered in a statistical manner. Consider the distributions of all data, that represent numbers below 50 (chosen since it is a characteristic that is to set up an experiment for). Recall all 7 digits are alpha numeric, with only the third and fourth digits being numeric<jupyter_code>distns_below_50 = {k: v for k, v in distn_dict.items() if int(k[2:4]) < 50}
print('number of distributions with number below 50 = ' + str(len(distns_below_50)))<jupyter_output>number of distributions with number below 50 = 6103<jupyter_text>and use these to find only random vectors that are 'sufficiently far away' from the means of those reg numbers with numbers below 50. So the 'characteristic' chosen, will be for the number to be above 50; the reg numbers are of the form AAXXAAA, so the characteristic is that XX is greater than 50.Employing a cheap-and-cheerful example of rejection sampling, and a cheap-and-cheerful example of statistical hypothesis testing, one can obtain the a set of random vectors that are 'sufficiently far' in the tails of the distributions belonging to those data vectors whose registration numbers are below 50.<jupyter_code># magnitude tuned so that it is roughly the volume covered by the means of the distributions
random_vectors = 2 * (np.random.random((1000, latent_dim)) - 0.5)
generated_regs = [v.recover(x) for x in decoder.predict(random_vectors)]
print('reg numbers below 50 = ' + str(len([x for x in generated_regs if int(x[2:4]) < 50])))
print('reg numbers above 50 = ' + str(len([x for x in generated_regs if int(x[2:4]) >= 50])))
# for each of the random vectors, find the smallest Mahalanobis distance from all the distributions with
# registration numeric part less than 50
means_mtr = np.array([mu for mu, sigma in distns_below_50.values()])
sigma_mtr = np.array([sigma for mu, sigma in distns_below_50.values()])
min_mahalanobis_dists = [np.min(np.linalg.norm((x - means_mtr) / sigma_mtr, axis=1))
                         for x in random_vectors]

# set a threshold on the Mahalanobis distance (which roughly corresponds to a hypothesis
# test significance level) to perform some crude rejection sampling
mahalanobis_threshold = 3.0
# I'm not convinced the latent_dim should be in the denominator here....
filtered_random_vectors = [x for x, y in zip(random_vectors, np.array(min_mahalanobis_dists)) 
                           if y > mahalanobis_threshold]
print('theoretical probability of sample from normal distribution having mahalanobis distance greater')
print('than mahalanobis threshold = ' +str(2*(1 - math.erf(0.5 * mahalanobis_threshold))))
print('number of filtered random vectors = ' + str(len(filtered_random_vectors)))
if len(filtered_random_vectors) > 0:
    filtered_generated_regs = [v.recover(x) for x in decoder.predict(np.array(filtered_random_vectors))]
    filtered_n_reg_above_50 = len([x for x in filtered_generated_regs if int(x[2:4]) > 50])
    print('number of filtered reg numbers with numeric part above 50 = ' + str(filtered_n_reg_above_50))
    print('proportion of filtered reg numbers obeying chosen characteristic = '
          + str(filtered_n_reg_above_50/len(filtered_random_vectors)))
else:
    print('no filtered random vectors to work with')<jupyter_output>number of filtered reg numbers with numeric part above 50 = 93
proportion of filtered reg numbers obeying chosen characteristic = 0.8773584905660378<jupyter_text>Using the distributions to choose characteristics in the generated data works best when- the latent dimension is low- the random samples used in the latent space are not 'too large, i.e. not far away from the means of the distributionsThe higher the dimension of the latent space, the greater care must be taken that the random samples used are in a suitable region of the latent space.As a sanity check, consider the opposite chosen characterstic, and generating reg numbers where the numerics are below 50. This should yield similar results.<jupyter_code>distns_above_50 = {k: v for k, v in distn_dict.items() if int(k[2:4]) >= 50}

means_mtr = np.array([mu for mu, sigma in distns_above_50.values()])
sigma_mtr = np.array([sigma for mu, sigma in distns_above_50.values()])
min_mahalanobis_dists = [np.min(np.linalg.norm((x - means_mtr) / sigma_mtr, axis=1))
                         for x in random_vectors]

mahalanobis_threshold = 3.0  
# filter on random samples that are 'far away' from test vectors that violate characteristic choice
filtered_random_vectors = [x for x, y in zip(random_vectors, np.array(min_mahalanobis_dists)) 
                           if y > mahalanobis_threshold]
filtered_generated_regs = [v.recover(x) for x in decoder.predict(np.array(filtered_random_vectors))]
filtered_n_reg_below_50 = len([x for x in filtered_generated_regs if int(x[2:4]) < 50])
print('number of filtered reg numbers with numeric part below 50 = ' + str(filtered_n_reg_below_50))
print('proportion of filtered reg numbers obeying chosen characteristic = '
        + str(filtered_n_reg_below_50/len(filtered_random_vectors)))<jupyter_output>number of filtered reg numbers with numeric part below 50 = 108
proportion of filtered reg numbers obeying chosen characteristic = 0.8780487804878049
{'fraction_non_alphanumeric': 0.11568816169393648, 'fraction_numerical': 0.07340391401989092, 'mean_word_length': 4.930745814307458, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 48, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '580741', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5958, 'n_tokens_neox': 5071, 'n_words': 1629, 'path': 'python/notebooks/variational_autoencoding.ipynb', 'repo_name': 'nagi49000/uk-car-reg-numbers', 'size': '15585'}