Something went wrong. Try again.
[READ-ONLY] Mirror of https://github.com/emilyhunt/masters_project. My Master's project code, using a mixture density network to infer photometric redshifts to galaxies with spectroscopic redshifts as training data.
Something went wrong. Try again.
40 kB · 801 lines
Python
at master
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802"""Placeholder for a potential future class, depending on how I go about things.""" # todo: docstrings
import tensorflow as tfimport numpy as npimport matplotlib.pyplot as pltimport timeimport gcimport pandas as pdimport pymc3from matplotlib import cmfrom scripts import loss_funcsfrom scripts.twitter import calc_local_timefrom typing import Optional, Unionfrom sklearn.model_selection import train_test_splitfrom scipy.optimize import minimize as scipy_minimizefrom sklearn.preprocessing import RobustScalerfrom sklearn.preprocessing import StandardScalerfrom sklearn.preprocessing import MinMaxScaler
class MixtureDensityNetwork:
def __init__(self, loss_function, summary_directory: str, regularization: Optional[str]=None, regularization_scale=0.1, x_features: int=1, y_features: int=1, x_scaling: Optional[str]=None, y_scaling: Optional[str]=None, layer_sizes=15, mixture_components: int=5, learning_rate=0.001, convolution_layer: bool=False, convolution_window_size: int=5, convolution_stride: int=1) -> None: """Initialises a mixture density network in tensorflow given the specified (or default) parameters.
Args: loss_function (loss_funcs class): an instance of the desired loss function to use. summary_directory (str): location of where we want TensorBoard to write our summaries to. regularization (str): controls the type of weight regularisation. Accepts 'none' (default), 'L1' or 'L2'. x_features (int): number of input x data points. Default is 1. y_features (int): number of input y data points layer_sizes (int, list-like): size of one layer (int) or a list of different sizes of each layer. mixture_components (int): number of mixtures to try to use todo: fnstring out of date
Returns: None """ # Cast layer_sizes as a 1D numpy array and use it to calculate how many hidden layers there should be layer_sizes = np.array([layer_sizes]).flatten() hidden_layers = layer_sizes.size
# Create a tensorflow graph for this class with tf.name_scope('graph'): self.graph = tf.Graph() self.graph_output_names = loss_function.coefficient_names
# Setup our graph with self.graph.as_default(): # Set seeds for reproducibility self.set_seeds()
# Placeholders for input data with tf.variable_scope('data'):
if convolution_layer: self.x_placeholder = tf.placeholder(tf.float32, [None, x_features, 1], name='x') else: self.x_placeholder = tf.placeholder(tf.float32, [None, x_features], name='x') self.y_placeholder = tf.placeholder(tf.float32, [None, y_features], name='y')
# Decide on the type of weight co-efficient regularisation to use based on what the user specified with tf.variable_scope('regularization'): if regularization is None: self.regularisation_function = None #self.loss_from_regularisation = np.float32(0) elif regularization is 'L1': self.regularisation_function = tf.contrib.layers.l1_regularizer(regularization_scale) #self.loss_from_regularisation = tf.losses.get_regularization_loss() elif regularization is 'L2': self.regularisation_function = tf.contrib.layers.l2_regularizer(regularization_scale) #self.loss_from_regularisation = tf.losses.get_regularization_loss() else: raise ValueError('specified regularisation type is invalid or unsupported.')
# Setup of the requisite number of layers, kept in a list - this lets us use easy numerical indexes # (including just -1 to get to the last one) to access different hidden layers. i = 0 with tf.variable_scope('hidden_layers'): self.graph_layers = []
# Join layers to x data, with a convolution layer if specified if convolution_layer: # We 'squeeze' the conv layer to get rid of the last dimension self.graph_layers.append(tf.layers.flatten( tf.layers.conv1d(self.x_placeholder, layer_sizes[i], [convolution_window_size], strides=convolution_stride, activation=tf.nn.tanh, kernel_regularizer=self.regularisation_function, name='hidden_layer_1_convolution'), name='convolution_layer_flattening')) else: self.graph_layers.append(tf.layers.dense(self.x_placeholder, layer_sizes[i], activation=tf.nn.tanh, kernel_regularizer=self.regularisation_function, name='hidden_layer_1'))
# Join layers to each other from here on out i += 1 while i < hidden_layers: self.graph_layers.append(tf.layers.dense(self.graph_layers[i - 1], layer_sizes[i], activation=tf.nn.tanh, kernel_regularizer=self.regularisation_function, name='hidden_layer_' + str(i+1))) i += 1
# Setup of the outputs as a dictionary of output layers, by cycling over the names of outputs and the with tf.variable_scope('output_layers'): self.graph_output = {} for output_name, activation_function, bias_initializer, kernel_initializer in zip( loss_function.coefficient_names, loss_function.activation_functions, loss_function.bias_initializers, loss_function.kernel_initializers): self.graph_output[output_name] = tf.layers.dense(self.graph_layers[-1], mixture_components, activation=activation_function, kernel_regularizer=self.regularisation_function, name='output_' + output_name, bias_initializer=bias_initializer, kernel_initializer=kernel_initializer)
# Collect all the regularization losses with tf.variable_scope('regularization'): if regularization is None: self.loss_from_regularisation = np.float32(0) else: self.loss_from_regularisation = tf.losses.get_regularization_loss()
# Initialise the loss function (storing the user-specified one with the class) and training scheme with tf.variable_scope('loss_calculation'): self.loss_function = loss_function loss_function_residuals = self.loss_function.tensor_evaluate(self.y_placeholder, self.graph_output) self.loss_from_function = loss_function_residuals['total_residual']
# Get all the other different residuals from loss_from_function, either using or summary writing them residuals = list(loss_function_residuals.keys()) residuals.remove('total_residual')
# Record summaries of all the individual residuals self.summary_residuals = [] for a_residual in residuals: self.summary_residuals.append(tf.summary.scalar(a_residual, loss_function_residuals[a_residual]))
self.loss_total = tf.add(self.loss_from_function, self.loss_from_regularisation) self.train_function = tf.train.AdamOptimizer(learning_rate=learning_rate, name='optimizer').minimize(self.loss_total)
# Initialise a tensorflow session object using our lovely graph we just made, and initialise the variables self.session = tf.Session() self.session.run(tf.global_variables_initializer())
# Initialise some performance evaluation tools with tf.variable_scope('summaries'): self.summary_loss_function = tf.summary.scalar('loss_function', self.loss_from_function) self.summary_loss_regularisation = tf.summary.scalar('loss_regularisation', self.loss_from_regularisation) self.summary_loss_total = tf.summary.scalar('loss_total', self.loss_total)
# Create a summary of the tf.layers variables self.summary_of_variables = {} for a_var in tf.trainable_variables(): self.summary_of_variables[str(a_var)] = tf.summary.histogram(a_var.name, a_var)
# Make us a cheeky summary writer to write all of our summaries to file self.summary_writer = tf.summary.FileWriter(summary_directory, self.graph)
# Create a blank loss & time array we can append more stuff to later for recording the loss function evolution self.recording_steps = np.array([]) self.loss = np.array([]) self.epoch_times = np.array([])
# Initialise blank feed dictionaries and a data range list we use to constrain later MAP value guessing self.y_data_range = None self.training_data = None self.validation_data = None
# An exit reason to keep class-wide and a training_success bool that we can use to prevent plotting self.exit_reason = None self.training_success = True
# Store some deets about the conv layer self.convolution_layer = convolution_layer
# Initialise scalers for the x and y data. We keep these with the class because it means that if x_scaling is 'min_max': self.x_scaler = MinMaxScaler(feature_range=(0.1, 0.9)) # todo: feature ranges are hard coded! elif x_scaling is 'robust': self.x_scaler = RobustScaler() elif x_scaling is 'standard': self.x_scaler = StandardScaler() elif x_scaling is None: self.x_scaler = None else: raise NotImplementedError('selected x scaling type has not been implemented!')
if y_scaling is 'min_max': self.y_scaler = MinMaxScaler(feature_range=(0.1, 0.9)) elif y_scaling is 'robust': self.y_scaler = RobustScaler() elif y_scaling is 'standard': self.y_scaler = StandardScaler() elif y_scaling is None: self.y_scaler = None else: raise NotImplementedError('selected y scaling type has not been implemented!')
print('An MDN has been initialised!')
@staticmethod def set_seeds(seed=42): """Sets seed in numpy and tensorflow. Allows for repeatability!
Args: seed (int, float): universal seed to set in tf and numpy. Default is 42. """ np.random.seed(seed) tf.set_random_seed(seed)
def __del__(self): """Closes the tensorflow session.""" self.session.close()
def set_training_data(self, x_data, y_data): """Modifies the class-unique training feed dictionary to give to the network.
Args: x_data (any): independent variables to feed to the network. y_data (any): dependent variables to feed to the network. """ # Scale the x data if requested if self.x_scaler is not None: x_data = self.x_scaler.fit_transform(x_data)
# Scale the y data if requested if self.y_scaler is not None: y_data = self.y_scaler.fit_transform(y_data)
# Keep an idea of what the maximum and minimum training data range is. This is helpful for plotting later. self.y_data_range = [y_data.min(), y_data.max()]
# Set an extra dimension in the data if we're using a convolution layer if self.convolution_layer: x_data = np.expand_dims(x_data, axis=2)
# Add the new x_data, y_data self.training_data = {self.x_placeholder: x_data, self.y_placeholder: y_data}
def set_validation_data(self, x_data, y_data): """Modifies the class-unique validation feed dictionary to give to the network.
Args: x_data (any): independent variables to feed to the network. y_data (any): dependent variables to feed to the network. """ # Scale the x data if requested, using the same scaling parameters as the training data if self.x_scaler is not None: x_data = self.x_scaler.transform(x_data)
# Scale the y data if requested, using the same scaling parameters as the training data #if self.y_scaler is not None: # y_data = self.y_scaler.transform(y_data) todo: does this need ydata??!?!?
# Set an extra dimension in the data if we're using a convolution layer if self.convolution_layer: x_data = np.expand_dims(x_data, axis=2)
# Add the new x_data, y_data self.validation_data = {self.x_placeholder: x_data} #, self.y_placeholder: y_data}
def train(self, max_epochs: int=50, max_runtime: float=1., reporting_time: float=10., max_epochs_without_report: int=1000): """Trains the tensorflow graph for the specified amount of time.
Args: max_epochs (int): number of training epochs to run the network for. max_runtime (int): maximum number of hours the code will run for before exiting. reporting_time (float): how often, in seconds, we print to the console. max_epochs_without_report (int): longest time to wait without letting the user know how everything is going and storing diagnostic info.
Returns: A list of: 0. The exit reason 1. The epoch reached 2. A bool of whether or not the training was a success """ # optimise: a way to record loss in-place would make this a little bit faster # todo: exit automatically when loss function converges.
# Get a start time and calculate the maximum amount of time we can run for. start_time = time.time() cutoff_time = time.time() + max_runtime * 60**2
# Let the user know about our grand plans print('\n=== BEGINNING TRAINING ===') print('start_time = {}'.format(calc_local_time(start_time))) print('max_epochs = {}'.format(max_epochs)) print('max_runtime = {:.2f} hour'.format(max_runtime)) print('reporting_time = {:.1f} seconds'.format(reporting_time)) print('==========================')
# Correct for if training has been ran before. if self.recording_steps.size != 0: start_epoch = self.recording_steps[-1] else: start_epoch = 0 epoch = start_epoch max_epochs = start_epoch + max_epochs
with self.graph.as_default(): # Set some defaults so that we don't fail on the first run self.exit_reason = 'max_epochs reached' # Cheeky string that says why we finished training have_done_more_than_one_step = False # Stop stupid predictions being made too early epochs_per_report = 1 step_start_time = time.time()
# Cycle over, doing some running while epoch < max_epochs:
# Run for requisite number of epochs until refresh (this stops print from being spammed on fast code) epochs_in_this_report = 0 while epochs_in_this_report < epochs_per_report: # optimise: make logging toggleable # Train the network in a new epoch self.session.run(self.train_function, feed_dict=self.training_data) epoch += 1 epochs_in_this_report += 1
# Calculate the new epoch time and ETA step_end_time = time.time() epoch_time = (step_end_time - step_start_time) / epochs_per_report step_start_time = time.time() finish_time = (max_epochs - epoch) * epoch_time + step_start_time
# Record the loss self.recording_steps = np.append(self.recording_steps, epoch) self.loss = np.append(self.loss, self.session.run(self.loss_from_function, feed_dict=self.training_data)) self.epoch_times = np.append(self.epoch_times, epoch_time)
# Check for nan loss, and if not, then write a summary if np.isnan(self.loss[-1]): # todo: could replace this with an exception catch instead for when nan loss happens self.exit_reason = 'nan loss encountered' self.training_success = False else: # Write out the summaries summary_merge = tf.summary.merge_all() summary = self.session.run(summary_merge, feed_dict=self.training_data) self.summary_writer.add_summary(summary, epoch + start_epoch)
# Force garbage collection (sometimes helps to prevent memory issues... sometimes.) gc.collect()
# Work out how many epochs we can do before our next check-in epochs_per_report = int(np.ceil(reporting_time / epoch_time))
# Make sure we aren't gonna go over max_epochs if epoch + epochs_per_report > max_epochs: epochs_per_report = max_epochs - epoch if epochs_per_report > max_epochs_without_report: epochs_per_report = max_epochs_without_report
# Output some details on the last few epochs print('--------------------------') print('CURRENT TIME: {}'.format(calc_local_time(step_start_time))) print('epoch = {} ({:.1f}% done)'.format(epoch, epoch / float(max_epochs) * 100)) print('epoch_time = {:.3f} seconds'.format(epoch_time)) print('loss = {:.5f}'.format(self.loss[-1]))
if have_done_more_than_one_step: print('finish_time = {}'.format(calc_local_time(finish_time))) # todo: evaluate whether or not we're gonna finish before the cutoff time
# Decide if we need to end if step_start_time > cutoff_time: self.exit_reason = 'time limit reached' self.training_success = True break if self.exit_reason is 'nan loss encountered': self.training_success = False break
have_done_more_than_one_step = True
# Make sure to update the training success bool if it did indeed work if self.exit_reason == 'max_epochs reached': self.training_success = True
# It's all over! =( (we let the user know, unsurprisingly) print('=== ENDING TRAINING ======') print('reason = {}'.format(self.exit_reason)) print('epochs done = {}'.format(epoch)) print('==========================')
return [self.exit_reason, epoch, self.training_success]
def save_graph(self, location: str): """Saves a complete copy of the current network to a specified location. Args: location (str): place where you want it saved, including the filename. Returns: None """ # Make a tf.train.Saver object and use it to save the graph. with self.graph.as_default(): saver = tf.train.Saver() save_path = saver.save(self.session, location)
print("A copy of the model has been saved to {}".format(save_path))
def open_graph(self, location: str): """Saves a complete copy of the current network to a specified location. Args: location (str): place where you want it saved, including the filename. Returns: None """ # todo: make this work, as it requires doing some extra stuff with loss functions and writing shit to the class pass
def validate(self): """Returns mixture parameters for the code given the verification data.
Returns: Validation data mixture co-efficients in a dictionary that sorts them by name. """ print('Validating the graph on the validation data...') # Run the session on validation data with the aim of returning mixture co-efficients result = {}
# Cycle over the different output constants, writing their validation run to 'result' for a_constant in self.graph_output_names: result[a_constant] = self.session.run(self.graph_output[a_constant], feed_dict=self.validation_data)
return result
def plot_loss_function_evolution(self, start: int=0, end: int=-1, y_log: bool=False, figure_directory: Optional=None) -> None: """Returns a plot of the change of the loss function over time.
Args: start (int): start epoch to plot. end (int): end epoch to plot. Default: -1, which sets the end to the last training step. y_log (bool): if True, sets y axis to be logarithmic. """ print('Plotting the evolution of the loss function...')
# See if the user wants the end to be the most recent loss function evaluation. if end == -1: end = self.loss.size
# Plot some stuff! fig = plt.figure(figsize=(4, 6)) ax_loss = fig.add_subplot(2, 1, 1) ax_time = fig.add_subplot(2, 1, 2, sharex=ax_loss) ax_loss.plot(self.recording_steps, self.loss[start:end], 'r-', label='loss') ax_time.plot(self.recording_steps, self.epoch_times[start:end], 'b-', label='epoch time')
ax_loss.set_title('Loss function evolution') ax_loss.set_ylabel('Loss')
ax_time.set_xlabel('Epoch') ax_time.set_ylabel('Epoch time')
# Set the plot to be log if desired if y_log: plt.yscale('log')
# Save if desired if figure_directory is not None: plt.savefig(figure_directory + '_lossf.png')
plt.show()
def calculate_validation_stats(self, validation_mixtures, data_range: Union[tuple, list], reporting_interval: int=100, start_resolution: int=100, uncertainty_integration_resolution: int=2000, uncertainty_sigma_level: float=1.): """Calculates the MAP (maximum a posteriori), uncertainty and modality of a given set of mixture distributions.
Args: validation_mixtures (dict): the data from a .validate call. data_range (list-like of length 2): the allowed range of data. reporting_interval (int): how often to let the user know which objects we're working on. Default: 100. start_resolution (int): number of points to test against when finding the initial guess. uncertainty_integration_resolution (int): number of random variables to draw when integrating the pdfs to find the uncertainty levels. uncertainty_sigma_level (float): todo: implement this! calculate with erf(x / sqrt(2))
Returns: A pandas.DataFrame() containing map values, limits and a pdf modality count, addressed with keys 'map', 'lower', 'upper' and 'modality'.
""" # Only do this next step if prior training actually worked! if self.training_success is False: print('Prior training failed! Unable to calculate MAP values. Exiting calculate_MAP.') return pd.DataFrame()
print('Attempting to calculate the MAP values of all distributions...')
# Define a function to minimise, multiplied by -1 to make sure we're minimising not maximising def function_to_minimise(x_data, my_object_dictionary, my_loss_function): return -1 * my_loss_function.pdf_single_point(x_data, my_object_dictionary)
# Create blank arrays of np.nan values to populate with hopefully successful minimisations n_objects = validation_mixtures[self.graph_output_names[0]].shape[0] map_values = np.empty(n_objects) map_values[:] = np.nan upper_limits = map_values.copy() lower_limits = map_values.copy() cdf_constant = map_values.copy() cdf_multiplier = map_values.copy() modality = np.zeros(n_objects, dtype=int) valid_map_values = np.zeros(n_objects, dtype=int)
# A bit of setup for our initial guess of MAP values guess_x_range = np.linspace(data_range[0], data_range[1], num=start_resolution)
# Loop over each object and work out the stats for each i = 0 successes = 0 mean_number_of_iterations = 0.0 while i < n_objects:
# Make a dictionary that only has data on this specific galaxy object_dictionary = {} for a_name in self.graph_output_names: object_dictionary[a_name] = validation_mixtures[a_name][i]
# Make a sensible starting guess and look at the y values in parameter space guess_y_range = self.loss_function.pdf_multiple_points(guess_x_range, object_dictionary, sum_mixtures=True) starting_guess = guess_x_range[np.argmax(guess_y_range)]
# Estimate the modality of the distribution by looking for sign changes, using a method from: # https://stackoverflow.com/questions/2652368/how-to-detect-a-sign-change-for-elements-in-a-numpy-array # First, we find the difference between consecutive values and hence estimate the gradient guess_difference = guess_y_range[1:] - guess_y_range[0:-1]
# Then, we find the sign of the gradient (np.sign returns +1 or -1) and then look for signchanges, where # consecutive guess_sign values being added or subtracted make +2 or -2 (not zero.) We also only care about # sign changes that occur for points that have a y value at least 1% of the MAP (as we don't care about # tiny modes or accidentally detecting the flat tails either side of the distribution.) guess_sign = np.sign(guess_difference) guess_signchange = np.logical_and(guess_sign[1:] - guess_sign[0:-1] != 0, guess_y_range[1:-1] > guess_y_range.max() * 0.01)
# Count the number of turning points in the pdf turning_points = np.count_nonzero(guess_signchange)
# Convert this number of turning points into a measure of modality with some fancyness if turning_points != 0: modality[i] = int(np.ceil(0.5 + turning_points / 2)) else: modality[i] = 0
# Attempt to minimise and find the MAP value result = scipy_minimize(function_to_minimise, np.array([starting_guess]), args=(object_dictionary, self.loss_function), method='Nelder-Mead', options={'disp': False})
# Store the result only if we're able to if result.success: map_values[i] = result.x valid_map_values[i] = 1 successes += 1 mean_number_of_iterations += result.nit
# Now, use random deviates to find the uncertainty on the MAP random_deviates = self.loss_function.draw_random_variables(uncertainty_integration_resolution, object_dictionary)
# pymc3 docs for the following functions: https://docs.pymc.io/api/stats.html # Use the highest probability density region if the pdf is unimodal if modality[i] == 1: limits = pymc3.stats.hpd(random_deviates, alpha=0.3173) # Corresponds to 1 sigma of error. lower_limits[i] = limits[0] upper_limits[i] = limits[1]
# If not, then we find quantiles corresponding to +- 1 sigma. This isn't theoretically as # good, as the IQR is not guaranteed to contain the MAP value we found earlier. else: limits = pymc3.stats.quantiles(random_deviates, qlist=[15.865, 84.135]) # Again, 1 sigma error. lower_limits[i] = limits[15.865] upper_limits[i] = limits[84.135]
# Set the MAP to the median if the MAP we found isn't within the limits if map_values[i] > upper_limits[i] or map_values[i] < lower_limits[i]: map_values[i] = np.median(random_deviates) valid_map_values[i] = 2
# Calculate constants to let us normalise the CDFs later cdf_stuff = self.loss_function.cdf_multiple_points(data_range, object_dictionary) cdf_constant[i] = -1. * cdf_stuff[0] cdf_multiplier[i] = 1. / cdf_stuff[1]
else: print('Failed to find MAP for object {}!'.format(i))
# Keep the user updated on what interval of objects we're working on (prevents panic if this takes ages) if i % reporting_interval == 0: print('Working on objects {} to {}...'.format(i, i+reporting_interval))
i += 1
# Scale the data in reverse if we're using a y-data scaler if self.y_scaler is not None: finite_map_values = np.isfinite(map_values) map_values[finite_map_values] = self.y_scaler.inverse_transform(map_values[finite_map_values] .reshape(-1, 1)).flatten() finite_lower_limits = np.isfinite(lower_limits) lower_limits[finite_lower_limits] = self.y_scaler.inverse_transform(lower_limits[finite_lower_limits] .reshape(-1, 1)).flatten() finite_upper_limits = np.isfinite(upper_limits) upper_limits[finite_upper_limits] = self.y_scaler.inverse_transform(upper_limits[finite_upper_limits] .reshape(-1, 1)).flatten()
# Calculate the mean number of iterations of the minimizer mean_number_of_iterations /= float(successes)
print('Found MAP values for {:.2f}% of objects.'.format(100 * successes / float(n_objects))) print('Mean number of MAP minimisation iterations = {}'.format(mean_number_of_iterations))
return pd.DataFrame({'map': map_values, 'lower': lower_limits, 'upper': upper_limits, 'modality': modality, 'valid': valid_map_values, 'cdf_constant': cdf_constant, 'cdf_multiplier': cdf_multiplier})
def plot_pdf(self, validation_data: dict, values_to_highlight, data_range=None, resolution: int=100, map_values=None, true_values=None, figure_directory: Optional[str]=None, show_fig: bool=False): """Plots the mixture pdf of a given set of parameters.
# todo: true_values should get flattened first as there can be issues with it being an odd shape
Args: validation_data (dict): as returned by network.validate, this is the validation data to plot with. values_to_highlight (int, list-like of ints): IDs of the objects to plot pdfs for. Default: None. data_range (list-like of floats): The range resolution (int): how many points to evaluate the pdf at. map_values ( #todo FINISH YOUR FUCKING DOCSTRINGS EMILY FFS
Returns: pretty graphs """ # Only do this next step if prior training actually worked! if self.training_success is False: print('Prior training failed! Unable to plot PDFs. Exiting calculate_MAP.') return 0
# See if the user specified their own data range if data_range is None: data_range = self.y_data_range
# Typecast the list of values to highlight as a numpy 1D array values_to_highlight = np.array([values_to_highlight]).flatten() n_mixtures = validation_data[self.graph_output_names[0]][0, :].size
# Cycle over intrinsic_start and intrinsic_end. We scale back into actual value space later. y_range = np.linspace(data_range[0], data_range[1], num=resolution) if self.y_scaler is not None: actual_y_range = self.y_scaler.inverse_transform(y_range.reshape(-1, 1)).flatten() else: actual_y_range = y_range
# Setup a list of colours to plot each mixture with colors = cm.viridis(np.linspace(0, 1, n_mixtures))
# Evaluate all the different pdfs. for an_object in values_to_highlight:
# Make a dictionary that only has data on this specific galaxy object_dictionary = {} for a_name in self.graph_output_names: object_dictionary[a_name] = validation_data[a_name][an_object]
# Get the pdf data, including a total probability mixture_pdfs = self.loss_function.pdf_multiple_points(y_range, object_dictionary, sum_mixtures=False) total_pdf = np.sum(mixture_pdfs, axis=0)
# Make some cool plots plt.figure()
# Plot the total mixture plt.plot(actual_y_range, total_pdf, 'k-', lw=2, label='Total pdf', alpha=1)
# Plot each mixture individually for mixture_number, a_color in enumerate(colors): plt.plot(actual_y_range, mixture_pdfs[mixture_number, :], '--', lw=1, color=a_color, label='Mixture ' + str(mixture_number))
# If specified, plot the MAP value if map_values is not None: plt.plot([map_values[an_object], map_values[an_object]], [total_pdf.min(), total_pdf.max()], 'r--', label='MAP value')
# If specified, plot the MAP value if true_values is not None: plt.plot([true_values[an_object], true_values[an_object]], [total_pdf.min(), total_pdf.max()], 'b--', label='True value')
plt.legend(edgecolor='k', facecolor='w', fancybox=True) plt.title('PDF of object ' + str(an_object))
# Save if desired if figure_directory is not None: plt.savefig(figure_directory + '_pdf_' + str(an_object) + '.png')
if show_fig: plt.show() else: plt.close() # todo make this stuff work
# Unit tests: implements the class on the toy_mdn_emily example data, using data from the following blog post:# http://blog.otoro.net/2015/11/24/mixture-density-networks-with-tensorflow/if __name__ == '__main__': print('Commencing mdn.py unit tests!')
# Create some data to play wit1h def build_toy_dataset(dataset_size): y_data = np.random.uniform(5, 10.5, dataset_size) # DEFAULT: -10.5, +10.5 r_data = np.random.normal(size=dataset_size) # random noise x_data = np.sin(0.75 * y_data) * 7.0 + y_data * 0.5 + r_data * 0.1 x_data = x_data.reshape((dataset_size, 1)) y_data = y_data.reshape((dataset_size, 1)) return train_test_split(x_data, y_data, random_state=42)
points = 5000
x_train, x_test, y_train, y_test = build_toy_dataset(points) print("Size of features in training data: {}".format(x_train.shape)) print("Size of output in training data: {}".format(y_train.shape)) print("Size of features in test data: {}".format(x_test.shape)) print("Size of output in test data: {}".format(y_test.shape)) # plt.plot(x_train, y_train, 'or', mew=0, ms=3, alpha=0.2) # plt.title('Training data') # plt.show()
# Initialise the network network = MixtureDensityNetwork(loss_funcs.NormalPDFLoss(), './logs/mdn_tests_tensorboard/' + str(time.strftime('%H-%M-%S', time.localtime(time.time()))), regularization=None, x_features=1, y_features=1, layer_sizes=[20, 20], mixture_components=15, learning_rate=5e-3, y_scaling=None, x_scaling=None)
# Set the data network.set_training_data(x_train, y_train) network.set_validation_data(x_test, y_test)
# Train the network for max_epochs epochs network.train(max_epochs=250)
# Plot the loss function network.plot_loss_function_evolution()
# Validate the network validation_results = network.validate()
# Calculate some MAP values validation_stats = network.calculate_validation_stats(validation_results, reporting_interval=500)
# Plot some pdfs network.plot_pdf(validation_results, [0, 100, 200], map_values=validation_stats['map'], show_fig=True)
# Some code to plot a validation plot. Firstly, we have to generate points to pull from: def generate_points(my_x_test, my_weights, my_means, my_std_deviations): """Generates points randomly given a loada points. Uses uniform deviates to guess a mixture coefficient to use. Then draws a point randomly from said selected distribution. We do this instead of picking the mode because we're fitting a model to data with intrinsic scatter!""" n_test_points = my_x_test.size
# Pick out some uniform deviates mixtures_to_use = np.random.rand(n_test_points).reshape(n_test_points, 1)
# Create a cumulative sum of weights my_weights_sum = np.cumsum(my_weights, axis=1)
# Find the first argument that's greater than the uniform deviate (since np.argmax stops at the first instance) random_weights_indexes = np.argmax(np.greater(my_weights_sum, mixtures_to_use), axis=1)
# Grab the random means and standard deviations random_means = my_means[np.arange(0, n_test_points), random_weights_indexes] random_std_deviations = my_std_deviations[np.arange(0, n_test_points), random_weights_indexes]
# Use these parameters to make some random numbers that are normal distributed return np.random.normal(loc=random_means, scale=random_std_deviations)
# Make some points y_test_random = generate_points(x_test, validation_results[network.graph_output_names[0]], validation_results[network.graph_output_names[1]], validation_results[network.graph_output_names[2]])
# Plot some stuff plt.figure() plt.plot(x_test, y_test, 'or', mew=0, ms=3, alpha=0.5, label='Training data') plt.plot(x_test, y_test_random, 'ob', mew=0, ms=3, alpha=0.5, label='Predictions') plt.title('Network prediction vs training data') plt.legend(fancybox=True) plt.ylim(-25, 25) plt.show()