Genetic_Algorithm/main.py

+#!/usr/bin/env python
+
+import logging
+from optimizer import Optimizer
+from tqdm import tqdm
+
+# Setup logging.
+logging.basicConfig(
+    format='%(asctime)s - %(levelname)s - %(message)s',
+    datefmt='%m/%d/%Y %I:%M:%S %p',
+    level=logging.DEBUG,
+    filename='log.txt'
+)
+
+def train_networks(networks):
+    """Train each network.
+
+    Args:
+        networks (list): Current population of networks
+    """
+    pbar = tqdm(total=len(networks))
+    for network in networks:
+        network.train()
+        pbar.update(1)
+    pbar.close()
+
+def get_average_accuracy(networks):
+    """Get the average accuracy for a group of networks.
+
+    Args:
+        networks (list): List of networks
+
+    Returns:
+        float: The average accuracy of a population of networks.
+
+    """
+    total_accuracy = 0
+    for network in networks:
+        total_accuracy += network.accuracy
+
+    return total_accuracy / len(networks)
+
+def generate(generations, population, nn_param_choices):
+    """Generate a network with the genetic algorithm.
+
+    Args:
+        generations (int): Number of times to evole the population
+        population (int): Number of networks in each generation
+        nn_param_choices (dict): Parameter choices for networks
+
+    """
+    optimizer = Optimizer(nn_param_choices)
+    networks = optimizer.create_population(population)
+
+    # Evolve the generation.
+    for i in range(generations):
+        logging.info("***Doing generation %d of %d***" %
+                     (i + 1, generations))
+
+        # Train and get accuracy for networks.
+        train_networks(networks)
+
+        # Get the average accuracy for this generation.
+        average_accuracy = get_average_accuracy(networks)
+
+        # Print out the average accuracy each generation.
+        logging.info("Generation average: %.2f%%" % (average_accuracy * 100))
+        logging.info('-'*80)
+
+        # Evolve, except on the last iteration.
+        if i != generations - 1:
+            # Do the evolution.
+            networks = optimizer.evolve(networks)
+
+    # Sort our final population.
+    networks = sorted(networks, key=lambda x: x.accuracy, reverse=True)
+
+    # Print out the top 5 networks.
+    print_networks(networks[:5])
+
+def print_networks(networks):
+    """Print a list of networks.
+
+    Args:
+        networks (list): The population of networks
+
+    """
+    logging.info('-'*80)
+    for network in networks:
+        network.print_network()
+
+def main():
+    """Evolve a network."""
+    generations = 10  # Number of times to evole the population.
+    population = 20  # Number of networks in each generation.
+
+    nn_param_choices = {
+        'nb_neurons': [64, 128, 256, 512, 768, 1024],
+        'nb_layers': [1, 2, 3, 4],
+        'activation': ['relu', 'elu', 'tanh', 'sigmoid'],
+        'optimizer': ['rmsprop', 'adam', 'sgd', 'adagrad',
+                      'adadelta', 'adamax', 'nadam'],
+    }
+
+    logging.info("***Evolving %d generations with population %d***" %
+                 (generations, population))
+
+    generate(generations, population, nn_param_choices)
+
+if __name__ == '__main__':
+    main()

Genetic_Algorithm/network.py

+#!/usr/bin/env python
+
+"""Class that represents the network to be evolved."""
+import random
+import logging
+from train import train_and_score
+
+class Network(object):
+    """Represent a network and let us operate on it.
+
+    Currently only works for an MLP.
+    """
+
+    def __init__(self, nn_param_choices=None):
+        """Initialize our network.
+
+        Args:
+            nn_param_choices (dict): Parameters for the network, includes:
+                nb_neurons (list): [64, 128, 256]
+                nb_layers (list): [1, 2, 3, 4]
+                activation (list): ['relu', 'elu']
+                optimizer (list): ['rmsprop', 'adam']
+        """
+        self.accuracy = 0.
+        self.nn_param_choices = nn_param_choices
+        self.network = {}  # (dic): represents MLP network parameters
+
+    def create_random(self):
+        """Create a random network."""
+        for key in self.nn_param_choices:
+            self.network[key] = random.choice(self.nn_param_choices[key])
+
+    def create_set(self, network):
+        """Set network properties.
+
+        Args:
+            network (dict): The network parameters
+
+        """
+        self.network = network
+
+    def train(self):
+        """Train the network and record the accuracy.
+
+        Args:
+
+        """
+        if self.accuracy == 0.:
+            self.accuracy = train_and_score(self.network)
+
+    def print_network(self):
+        """Print out a network."""
+        logging.info(self.network)
+        logging.info("Network accuracy: %.2f%%" % (self.accuracy * 100))

Genetic_Algorithm/optimizer.py

+#!/usr/bin/env python
+
+"""
+Class that holds a genetic algorithm for evolving a network.
+
+Credit:
+    A lot of those code was originally inspired by:
+    http://lethain.com/genetic-algorithms-cool-name-damn-simple/
+"""
+from functools import reduce
+from operator import add
+import random
+from network import Network
+
+class Optimizer(object):
+    """Class that implements genetic algorithm for MLP optimization."""
+
+    def __init__(self, nn_param_choices, retain=0.4,
+                 random_select=0.1, mutate_chance=0.2):
+        """Create an optimizer.
+
+        Args:
+            nn_param_choices (dict): Possible network paremters
+            retain (float): Percentage of population to retain after
+                each generation
+            random_select (float): Probability of a rejected network
+                remaining in the population
+            mutate_chance (float): Probability a network will be
+                randomly mutated
+
+        """
+        self.mutate_chance = mutate_chance
+        self.random_select = random_select
+        self.retain = retain
+        self.nn_param_choices = nn_param_choices
+
+    def create_population(self, count):
+        """Create a population of random networks.
+
+        Args:
+            count (int): Number of networks to generate, aka the
+                size of the population
+
+        Returns:
+            (list): Population of network objects
+
+        """
+        pop = []
+        for _ in range(0, count):
+            # Create a random network.
+            network = Network(self.nn_param_choices)
+            network.create_random()
+
+            # Add the network to our population.
+            pop.append(network)
+
+        return pop
+
+    @staticmethod
+    def fitness(network):
+        """Return the accuracy, which is our fitness function."""
+        return network.accuracy
+
+    def grade(self, pop):
+        """Find average fitness for a population.
+
+        Args:
+            pop (list): The population of networks
+
+        Returns:
+            (float): The average accuracy of the population
+
+        """
+        summed = reduce(add, (self.fitness(network) for network in pop))
+        return summed / float((len(pop)))
+
+    def breed(self, mother, father):
+        """Make two children as parts of their parents.
+
+        Args:
+            mother (dict): Network parameters
+            father (dict): Network parameters
+
+        Returns:
+            (list): Two network objects
+
+        """
+        children = []
+        for _ in range(2):
+
+            child = {}
+
+            # Loop through the parameters and pick params for the kid.
+            for param in self.nn_param_choices:
+                child[param] = random.choice(
+                    [mother.network[param], father.network[param]]
+                )
+
+            # Now create a network object.
+            network = Network(self.nn_param_choices)
+            network.create_set(child)
+
+            # Randomly mutate some of the children.
+            if self.mutate_chance > random.random():
+                network = self.mutate(network)
+
+            children.append(network)
+
+        return children
+
+    def mutate(self, network):
+        """Randomly mutate one part of the network.
+
+        Args:
+            network (dict): The network parameters to mutate
+
+        Returns:
+            (Network): A randomly mutated network object
+
+        """
+        # Choose a random key.
+        mutation = random.choice(list(self.nn_param_choices.keys()))
+
+        # Mutate one of the params.
+        network.network[mutation] = random.choice(self.nn_param_choices[mutation])
+
+        return network
+
+    def evolve(self, pop):
+        """Evolve a population of networks.
+
+        Args:
+            pop (list): A list of network parameters
+
+        Returns:
+            (list): The evolved population of networks
+
+        """
+        # Get scores for each network.
+        graded = [(self.fitness(network), network) for network in pop]
+
+        # Sort on the scores.
+        graded = [x[1] for x in sorted(graded, key=lambda x: x[0], reverse=True)]
+
+        # Get the number we want to keep for the next gen.
+        retain_length = int(len(graded)*self.retain)
+
+        # The parents are every network we want to keep.
+        parents = graded[:retain_length]
+
+        # For those we aren't keeping, randomly keep some anyway.
+        for individual in graded[retain_length:]:
+            if self.random_select > random.random():
+                parents.append(individual)
+
+        # Now find out how many spots we have left to fill.
+        parents_length = len(parents)
+        desired_length = len(pop) - parents_length
+        children = []
+
+        # Add children, which are bred from two remaining networks.
+        while len(children) < desired_length:
+
+            # Get a random mom and dad.
+            male = random.randint(0, parents_length-1)
+            female = random.randint(0, parents_length-1)
+
+            # Assuming they aren't the same network...
+            if male != female:
+                male = parents[male]
+                female = parents[female]
+
+                # Breed them.
+                babies = self.breed(male, female)
+
+                # Add the children one at a time.
+                for baby in babies:
+                    # Don't grow larger than desired length.
+                    if len(children) < desired_length:
+                        children.append(baby)
+
+        parents.extend(children)
+
+        return parents

Genetic_Algorithm/train.py

+#!/usr/bin/env python
+
+import os,sys,inspect
+currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
+parentdir = os.path.dirname(currentdir)
+sys.path.insert(0,parentdir) 
+import toolkit
+from toolkit import KerasROOTClassification
+
+def init_model(network):
+    
+    nb_layers = network['nb_layers']
+    nb_neurons = network['nb_neurons']
+    activation = network['activation']
+    optimizer = network['optimizer']
+
+    filename = "/project/etp4/nhartmann/trees/allTrees_m1.8_NoSys.root"
+    
+    c = KerasROOTClassification("",
+                                signal_trees = [(filename, "GG_oneStep_1705_1105_505_NoSys")],
+                                bkg_trees = [(filename, "ttbar_NoSys"),
+                                             (filename, "wjets_Sherpa221_NoSys")
+                                ],
+                                dumping_enabled=False,
+                                optimizer=optimizer,
+                                layers=nb_layers,
+                                nodes=nb_neurons,
+                                activation_function=activation,
+                                #optimizer_opts=dict(lr=100., decay=1e-6, momentum=0.9),
+                                earlystopping_opts=dict(monitor='val_loss',
+                                    min_delta=0, patience=2, verbose=0, mode='auto'),
+                                # optimizer="Adam",
+                                selection="lep1Pt<5000", # cut out a few very weird outliers
+                                branches = ["met", "mt"],
+                                weight_expr = "eventWeight*genWeight",
+                                identifiers = ["DatasetNumber", "EventNumber"],
+                                step_bkg = 100)
+    return c
+
+def train_and_score(network):
+    model = init_model(network)
+
+    model.train()
+
+    score = model.score
+
+    return score[1]  # 1 is accuracy. 0 is loss.

toolkit.py

@@ -44,19 +44,21 @@ class KerasROOTClassification(object):
 
 
     # Datasets that are stored to (and dynamically loaded from) hdf5
-    dataset_names = ["x_train", "x_test", "y_train", "y_test", "w_train", "w_test", "scores_train", "scores_test"]
+    dataset_names = ["x_train", "x_test", "y_train", "y_test", "w_train", "w_test", "pred_train", "pred_test"]
 
     # Datasets that are retrieved from ROOT trees the first time
     dataset_names_tree = ["x_train", "x_test", "y_train", "y_test", "w_train", "w_test"]
 
     def __init__(self, name, *args, **kwargs):
         self._init_from_args(name, *args, **kwargs)
-        with open(os.path.join(self.project_dir, "options.json"), "w") as of:
-            json.dump(dict(args=args, kwargs=kwargs), of)
+        if self.dumping_enabled:
+            with open(os.path.join(self.project_dir, "options.json"), "w") as of:
+                json.dump(dict(args=args, kwargs=kwargs), of)
 
 
     def _init_from_args(self, name,
                         signal_trees, bkg_trees, branches, weight_expr, identifiers,
+                        dumping_enabled=True,
                         selection=None,
                         layers=3,
                         nodes=64,
@@ -72,6 +74,7 @@ class KerasROOTClassification(object):
                         earlystopping_opts=None):
 
         self.name = name
+        self.dumping_enabled = dumping_enabled
         self.signal_trees = signal_trees
         self.bkg_trees = bkg_trees
         self.branches = branches
@@ -114,8 +117,9 @@ class KerasROOTClassification(object):
         self._y_test = None
         self._w_train = None
         self._w_test = None
-        self._scores_train = None
-        self._scores_test = None
+        self.pred_train = None
+        self.pred_test = None
+        self.score = None
 
         self.s_eventlist_train = None
         self.b_eventlist_train = None
@@ -169,8 +173,8 @@ class KerasROOTClassification(object):
                                      branches=self.branches+[self.weight_expr],
                                      selection=self.selection,
                                      start=1, step=self.step_bkg)
-
-            self._dump_training_list()
+            if self.dumping_enabled:
+                self._dump_training_list()
             self.s_eventlist_train = self.s_train[self.identifiers]
             self.b_eventlist_train = self.b_train[self.identifiers]
 
@@ -196,7 +200,8 @@ class KerasROOTClassification(object):
             self.y_test[:len(self.s_test)] = 1
             self.y_test[len(self.s_test):] = 0
 
-            self._dump_to_hdf5(*self.dataset_names_tree)
+            if self.dumping_enabled:
+                self._dump_to_hdf5(*self.dataset_names_tree)
 
         self.data_loaded = True
 
@@ -261,7 +266,8 @@ class KerasROOTClassification(object):
                 # probably we either want to fit only training data or training and test data together
                 # logger.info("Fitting StandardScaler to test data")
                 # self._scaler.fit(self.x_test)
-                joblib.dump(self._scaler, filename)
+                if self.dumping_enabled:
+                    joblib.dump(self._scaler, filename)
         return self._scaler
 
 
@@ -354,8 +360,9 @@ class KerasROOTClassification(object):
                 logger.info("No weights found, starting completely new model")
 
             # dump to json for documentation
-            with open(os.path.join(self.project_dir, "model.json"), "w") as of:
-                of.write(self._model.to_json())
+            if self.dumping_enabled:
+                with open(os.path.join(self.project_dir, "model.json"), "w") as of:
+                    of.write(self._model.to_json())
 
         return self._model
 
@@ -386,9 +393,9 @@ class KerasROOTClassification(object):
         np.random.shuffle(self.y_train)
         np.random.set_state(rn_state)
         np.random.shuffle(self.w_train)
-        if self._scores_test is not None:
+        if self.pred_test is not None:
             np.random.set_state(rn_state)
-            np.random.shuffle(self._scores_test)
+            np.random.shuffle(self.pred_test)
 
 
     def train(self, epochs=10):
@@ -418,20 +425,25 @@ class KerasROOTClassification(object):
         except KeyboardInterrupt:
             logger.info("Interrupt training - continue with rest")
 
-        logger.info("Save history")
-        self._dump_history()
+        if self.dumping_enabled:
+            logger.info("Save history")
+            self._dump_history()
+
+            logger.info("Save weights")
+            self.model.save_weights(os.path.join(self.project_dir, "weights.h5"))
 
-        logger.info("Save weights")
-        self.model.save_weights(os.path.join(self.project_dir, "weights.h5"))
+            self.total_epochs += epochs
+            self._write_info("epochs", self.total_epochs)
 
-        self.total_epochs += epochs
-        self._write_info("epochs", self.total_epochs)
+        logger.info("Create/Update predictions for ROC curve")
+        self.pred_test = self.model.predict(self.x_test)
+        self.pred_train = self.model.predict(self.x_train)
 
-        logger.info("Create/Update scores for ROC curve")
-        self.scores_test = self.model.predict(self.x_test)
-        self.scores_train = self.model.predict(self.x_train)
+        logger.info("Get test loss and metrics of the model")
+        self.score = self.model.evaluate(self.x_test, self.y_test, verbose=0, sample_weight=None)
 
-        self._dump_to_hdf5("scores_train", "scores_test")
+        if self.dumping_enabled:
+            self._dump_to_hdf5("pred_train", "pred_test")
 
 
 
@@ -521,7 +533,7 @@ class KerasROOTClassification(object):
     def plot_ROC(self):
 
         logger.info("Plot ROC curve")
-        fpr, tpr, threshold = roc_curve(self.y_test, self.scores_test, sample_weight = self.w_test)
+        fpr, tpr, threshold = roc_curve(self.y_test, self.pred_test, sample_weight = self.w_test)
 
         fpr = 1.0 - fpr
         roc_auc = auc(tpr, fpr)
@@ -571,7 +583,10 @@ class KerasROOTClassification(object):
 def create_getter(dataset_name):
     def getx(self):
         if getattr(self, "_"+dataset_name) is None:
-            self._load_from_hdf5(dataset_name)
+            try:
+                self._load_from_hdf5(dataset_name)
+            except KeyError:
+                logger.info("KeyError")
         return getattr(self, "_"+dataset_name)
     return getx