diff --git a/Cargo.toml b/Cargo.toml index d3067a3..b76aa49 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -14,6 +14,7 @@ ndarray = { version = "0.15.6", features = ["serde"] } num-traits = "*" ordered-float = "*" rand = "*" +rand_distr = "*" serde = { version = "*", features = ["derive"] } serde_json = "*" tempfile = "*" diff --git a/src/individual.rs b/src/individual.rs index 8aba59f..07831bb 100644 --- a/src/individual.rs +++ b/src/individual.rs @@ -7,6 +7,7 @@ use std::fs::read_to_string; use std::io::Error as IOError; use std::path::Path; use std::ops::Index; +use std::slice::Iter; use log::trace; use num_traits::FromPrimitive; @@ -97,6 +98,16 @@ impl Individual { self.layers.len() } + /// Returns an iterator over the individual's layers. + pub fn iter(&self) -> Iter<'_, Layer> { + self.into_iter() + } + + /// Returns a reference to the individual's cost function. + pub fn get_cost_function(&self) -> &CostFunction { + &self.cost_function + } + /// Passes the `input` through the network and returns the intermediate results of each layer. /// /// # Arguments @@ -339,6 +350,17 @@ impl Index for Individual { } +/// Allows turning a reference to an `Individual` into an iterator over [`Layer`] references. +impl<'a, T: Tensor> IntoIterator for &'a Individual { + type Item = &'a Layer; + type IntoIter = Iter<'a, Layer>; + + fn into_iter(self) -> Self::IntoIter { + self.layers.iter() + } +} + + #[cfg(test)] mod tests { use super::*; diff --git a/src/layer.rs b/src/layer.rs index ad2078c..0e9412a 100644 --- a/src/layer.rs +++ b/src/layer.rs @@ -25,6 +25,10 @@ impl Layer { let activation = (self.activation)(&weighted_input); (weighted_input, activation) } + + pub fn size(&self) -> usize { + self.weights.shape().0 + } } diff --git a/src/lib.rs b/src/lib.rs index 9f41ba4..8d2c733 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -8,7 +8,7 @@ pub mod component; pub mod cost_function; pub mod individual; pub mod layer; -pub mod population; pub mod tensor; pub mod utils; +pub mod world; diff --git a/src/population.rs b/src/population.rs deleted file mode 100644 index 14ee029..0000000 --- a/src/population.rs +++ /dev/null @@ -1,130 +0,0 @@ -use std::collections::HashMap; - -use ordered_float::OrderedFloat; -use rand::distributions::Slice; -use rand::{Rng, thread_rng}; - -use crate::individual::Individual; -use crate::tensor::Tensor; - - -type SelectFunc = fn(&mut Population, &str) -> Vec<(usize, usize)>; -type ProcreateFunc = fn(&Individual, &Individual) -> Individual; -type DetermineSpeciesKey = fn(&Individual) -> String; - - -// TODO: The term "population" refers to a group of organisms of the **same** species. -// We should find a better term. -// It seems the term "community" or even "environment" might be suitable. -pub struct Population { - pub species: HashMap>>, - kill_weak_and_select_parents: SelectFunc, - procreate_pair: ProcreateFunc, // crossover, mutation - determine_species_key: DetermineSpeciesKey, - pub validation_data: (T, T), -} - - -impl Population { - pub fn new( - kill_weak_and_select_parents: SelectFunc, - procreate_pair: ProcreateFunc, - determine_species_key: DetermineSpeciesKey, - validation_data: (T, T), - ) -> Self { - Self { - species: HashMap::new(), - kill_weak_and_select_parents, - procreate_pair, - determine_species_key, - validation_data, - } - } - - // mutates sets of individuals within one species - pub fn apply_selection(&mut self, species_key: String) { - let pairs = (self.kill_weak_and_select_parents)(self, &species_key); - let species = self.species.get(&species_key).unwrap(); - let new_individuals: Vec> = pairs.into_iter().map( - |(index1, index2)| (self.procreate_pair)(&species[index1], &species[index2]) - ).collect(); - for new_individual in new_individuals { - self.add_new_individual(new_individual) - } - } - - pub fn add_new_individual(&mut self, individual: Individual) { - let key = (self.determine_species_key)(&individual); - match self.species.get_mut(&key) { - Some(existing_species) => existing_species.push(individual), - None => { - self.species.insert(key, vec![individual]); - } - } - } -} - - -/// Reference implementation of a `kill_weak_and_select_parents` function. -/// -/// Calculates validation error for all individuals in the specifided species, -/// kills half of them (those with the highest errors), randomly samples parents -/// from the remaining individuals (using a uniform distribution with replacement) -/// and returns their indices as 2-tuples. -/// -/// # Arguments: -/// `population` - The population containing the species select from. -/// `species_key` - The key of the species to select from. -/// -/// # Returns: -/// Vector of 2-tuples of indices of the individuals in the specified species that should procreate. -pub fn select(population: &mut Population, species_key: &str) -> Vec<(usize, usize)> { - // Sort individuals in species by validation error. - let individuals = population.species.get_mut(species_key).unwrap(); - let (input, desired_output) = &population.validation_data; - individuals.sort_by_cached_key( - |individual| OrderedFloat(individual.calculate_error(input, desired_output)) - ); - // Kill the less fit half. - let num_individuals = individuals.len(); - let num_killed = individuals.drain(num_individuals/2..num_individuals).count(); - // Randomly sample (with a uniform distribution, with replacement) twice as many parents. - // If we wanted the probability to become a parent to depend on the validation error, - // we could use `rand::distributions::weighted::WeightedIndex` instead. - let rng = thread_rng(); - let indices = (0..(num_individuals - num_killed)).collect::>(); - let index_distribution = Slice::new(&indices).unwrap(); - rng.sample_iter(&index_distribution) - .take(num_killed * 2) - .array_chunks::<2>() - .map(|chunk| (*chunk[0], *chunk[1])) - .collect() -} - - -pub fn procreate(parent1: &Individual, parent2: &Individual) -> Individual { - parent1.clone() -} - - -pub fn get_species(individual: &Individual) -> String { - individual.num_layers().to_string() -} - - -#[cfg(test)] -mod tests { - use ndarray::array; - - use super::*; - - #[test] - fn test_new() { - let _pop = Population::new( - select, - procreate, - get_species, - (array![[0.]], array![[0.]]), - ); - } -} diff --git a/src/tensor.rs b/src/tensor.rs index 661def1..3c97092 100644 --- a/src/tensor.rs +++ b/src/tensor.rs @@ -3,7 +3,7 @@ use std::fmt::{Debug, Display}; use std::ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign}; -use ndarray::{Array2, Axis}; +use ndarray::{Array2, Axis, ArrayView}; use num_traits::FromPrimitive; use num_traits::real::Real; use serde::Serialize; @@ -30,6 +30,11 @@ pub trait TensorBase: /// Creates a tensor of the specified `shape` with all components equal to `num`. fn from_num(num: Self::Component, shape: (usize, usize)) -> Self; + + fn from_iter(iterator: I) -> Self + where + I: IntoIterator, + J: IntoIterator; /// Returns the shape of the tensor as a tuple of unsigned integers. fn shape(&self) -> (usize, usize); @@ -40,6 +45,12 @@ pub trait TensorBase: /// Returns the transpose of itself as a new tensor. fn transpose(&self) -> Self; + /// Append a row to the tensor. + fn append_row(&mut self, row: &[Self::Component]); + + /// Append a column to the tensor. + fn append_column(&mut self, column: &[Self::Component]); + /// Calls function `f` on each component and returns the result as a new tensor. fn map(&self, f: F) -> Self where F: FnMut(Self::Component) -> Self::Component; @@ -48,6 +59,12 @@ pub trait TensorBase: fn map_inplace(&mut self, f: F) where F: FnMut(Self::Component) -> Self::Component; + fn indexed_iter(&self) -> impl Iterator; + + fn indexed_iter_mut(&mut self) -> impl Iterator; + + fn iter(&self) -> impl Iterator; + /// Returns the sum of all rows (0) or columns (1) as a new tensor. fn sum_axis(&self, axis: usize) -> Self; } @@ -158,6 +175,31 @@ impl TensorBase for Array2 { fn from_num(num: Self::Component, shape: (usize, usize)) -> Self { Self::from_elem(shape, num) } + + fn from_iter(iterator: I) -> Self + where + I: IntoIterator, + J: IntoIterator, + { + let mut data = Vec::new(); + let mut num_rows: usize = 0; + let mut num_columns: Option = None; + for row in iterator { + num_rows += 1; + let row_vec: Vec = row.into_iter().collect(); + match num_columns { + Some(n) => { + if n != row_vec.len() { panic!() } + }, + None => num_columns = Some(row_vec.len()), + } + data.extend_from_slice(&row_vec); + } + Array2::from_shape_vec( + (num_rows, num_columns.unwrap_or(0)), + data, + ).unwrap() + } fn shape(&self) -> (usize, usize) { self.dim() @@ -174,6 +216,14 @@ impl TensorBase for Array2 { self.t().to_owned() } + fn append_row(&mut self, row: &[Self::Component]) { + self.push_row(ArrayView::from(row)).unwrap() + } + + fn append_column(&mut self, column: &[Self::Component]) { + self.push_column(ArrayView::from(column)).unwrap() + } + fn map(&self, f: F) -> Self where F: FnMut(C) -> C { self.mapv(f) @@ -184,6 +234,18 @@ impl TensorBase for Array2 { self.mapv_inplace(f) } + fn indexed_iter(&self) -> impl Iterator { + Array2::::indexed_iter(self) + } + + fn indexed_iter_mut(&mut self) -> impl Iterator { + Array2::::indexed_iter_mut(self) + } + + fn iter(&self) -> impl Iterator { + Array2::::iter(self) + } + fn sum_axis(&self, axis: usize) -> Self { self.sum_axis(Axis(axis)).insert_axis(Axis(axis)) } diff --git a/src/world.rs b/src/world.rs new file mode 100644 index 0000000..d0298cd --- /dev/null +++ b/src/world.rs @@ -0,0 +1,352 @@ +use std::cmp::{min, max}; +use std::collections::HashMap; + +use num_traits::{One, Zero}; +use ordered_float::OrderedFloat; +use rand::distributions::{Slice, WeightedIndex}; +use rand::rngs::ThreadRng; +use rand::seq::SliceRandom; +use rand::{Rng, thread_rng}; +use rand_distr::{Distribution, Poisson}; + +use crate::component::TensorComponent; +use crate::individual::Individual; +use crate::layer::Layer; +use crate::tensor::Tensor; + + +type SelectFunc = fn(&mut World, &str) -> Vec<(usize, usize)>; +type ProcreateFunc = fn(&Individual, &Individual) -> Individual; +type DetermineSpeciesKey = fn(&Individual) -> String; + + +pub struct World { + pub species: HashMap>>, + kill_weak_and_select_parents: SelectFunc, + procreate_pair: ProcreateFunc, // crossover, mutation + determine_species_key: DetermineSpeciesKey, + pub validation_data: (T, T), +} + + +impl World { + pub fn new( + kill_weak_and_select_parents: SelectFunc, + procreate_pair: ProcreateFunc, + determine_species_key: DetermineSpeciesKey, + validation_data: (T, T), + ) -> Self { + Self { + species: HashMap::new(), + kill_weak_and_select_parents, + procreate_pair, + determine_species_key, + validation_data, + } + } + + // mutates sets of individuals within one species + pub fn apply_selection(&mut self, species_key: String) { + let pairs = (self.kill_weak_and_select_parents)(self, &species_key); + let species = self.species.get(&species_key).unwrap(); + let new_individuals: Vec> = pairs.into_iter().map( + |(index1, index2)| (self.procreate_pair)(&species[index1], &species[index2]) + ).collect(); + for new_individual in new_individuals { + self.add_new_individual(new_individual) + } + } + + pub fn add_new_individual(&mut self, individual: Individual) { + let key = (self.determine_species_key)(&individual); + match self.species.get_mut(&key) { + Some(existing_species) => existing_species.push(individual), + None => { + self.species.insert(key, vec![individual]); + } + } + } +} + + +/// Reference implementation of a `kill_weak_and_select_parents` function. +/// +/// Calculates validation error for all individuals in the specifided species, +/// kills half of them (those with the highest errors), randomly samples parents +/// from the remaining individuals (using a uniform distribution with replacement) +/// and returns their indices as 2-tuples. +/// +/// # Arguments: +/// `world` - The world containing the species to select from. +/// `species_key` - The key of the species to select from. +/// +/// # Returns: +/// Vector of 2-tuples of indices of the individuals in the specified species that should procreate. +pub fn select(world: &mut World, species_key: &str) -> Vec<(usize, usize)> { + // Sort individuals in species by validation error. + let individuals = world.species.get_mut(species_key).unwrap(); + let (input, desired_output) = &world.validation_data; + individuals.sort_by_cached_key( + |individual| OrderedFloat(individual.calculate_error(input, desired_output)) + ); + // Kill the less fit half. + let num_individuals = individuals.len(); + let num_killed = individuals.drain(num_individuals/2..num_individuals).count(); + // Randomly sample (with a uniform distribution, with replacement) twice as many parents. + // If we wanted the probability to become a parent to depend on the validation error, + // we could use `rand::distributions::weighted::WeightedIndex` instead. + let rng = thread_rng(); + let indices = (0..(num_individuals - num_killed)).collect::>(); + let index_distribution = Slice::new(&indices).unwrap(); + rng.sample_iter(&index_distribution) + .take(num_killed * 2) + .array_chunks::<2>() + .map(|chunk| (*chunk[0], *chunk[1])) + .collect() +} + + +fn init_weight(rng: &mut ThreadRng, can_zero: bool) -> Option { + if can_zero && *[true, false].choose(rng).unwrap() { + Some(C::zero()) + } else { + C::from_f32(rng.gen_range(0f32..1f32)) + } +} + + +fn crossover_layer_weights( + new_weights: &mut T, + weights_p1: &T, + weights_p2: &T, + rng: &mut ThreadRng, +) { + let (rows_p1, cols_p1) = weights_p1.shape(); + let (rows_p2, cols_p2) = weights_p2.shape(); + for ((row, col), component) in new_weights.indexed_iter_mut() { + // both parents have a component at row|col + if row < rows_p1 && col < cols_p1 && row < rows_p2 && col < cols_p2 { + let weight_p1 = weights_p1[[row, col]]; + let weight_p2 = weights_p2[[row, col]]; + *component = *[weight_p1, weight_p2].choose(rng).unwrap(); + } + // only parent_1 has a component at row|col + else if row < rows_p1 && col < cols_p1 { + *component = weights_p1[[row, col]]; + } + // only parent_2 has a component at row|col + else if row < rows_p2 && col < cols_p2 { + *component = weights_p2[[row, col]]; + } + // neither has a component at row|col + else { + *component = init_weight(rng, true).unwrap() + } + } +} + + +fn crossover_layer_biases( + new_biases: &mut T, + biases_p1: &T, + biases_p2: &T, + rng: &mut ThreadRng, +) { + let rows_p1 = biases_p1.shape().0; + let rows_p2 = biases_p2.shape().0; + for ((row, _), component) in new_biases.indexed_iter_mut() { + // both parents have a component at row + if rows_p1 > row && rows_p2 > row { + let bias_p1 = biases_p1[[row, 0]]; + let bias_p2 = biases_p2[[row, 0]]; + *component = *[bias_p1, bias_p2].choose(rng).unwrap(); + } + // only parent_1 has a component at row + else if rows_p1 > row { + *component = biases_p1[[row, 0]]; + } + // only parent_2 has a component at row + else { + *component = biases_p2[[row, 0]]; + } + } +} + + +pub fn procreate(parent_1: &Individual, parent_2: &Individual) -> Individual { + let num_layers = parent_1.num_layers(); + if num_layers != parent_2.num_layers() { + panic!("Both parents must have the same number of layers!"); + } + let mut layers = Vec::with_capacity(num_layers); + let mut rng = thread_rng(); + let num_cols_p1 = parent_1[0].weights.shape().1; + let num_cols_p2 = parent_2[0].weights.shape().1; + let mut num_cols = rng.gen_range(min(num_cols_p1, num_cols_p2)..=max(num_cols_p1, num_cols_p2)); + for (layer_p1, layer_p2) in parent_1.iter().zip(parent_2) { + let num_rows_p1 = layer_p1.size(); + let num_rows_p2 = layer_p2.size(); + let num_rows = rng.gen_range(min(num_rows_p1, num_rows_p2)..=max(num_rows_p1, num_rows_p2)); + let mut weights = T::zeros((num_rows, num_cols)); + let mut biases = T::zeros((num_rows, 1)); + crossover_layer_weights(&mut weights, &layer_p1.weights, &layer_p2.weights, &mut rng); + crossover_layer_biases(&mut biases, &layer_p1.biases, &layer_p2.biases, &mut rng); + let activation = (*[&layer_p1.activation, &layer_p2.activation].choose(&mut rng).unwrap()).clone(); + layers.push(Layer::new(weights, biases, activation)); + num_cols = num_rows; + } + mutate_add_layer(&mut layers, &mut rng); + mutate_add_connections(&mut layers, &mut rng); + Individual::new(layers, parent_1.get_cost_function().clone()) +} + + +pub fn mutate_add_layer(layers: &mut Vec>, rng: &mut ThreadRng) { + if !should_add_layer(layers.len(), rng) { return; } + + let new_layer_idx = rng.gen_range(1..=layers.len()); + let following_layer = &mut layers[new_layer_idx]; + // columns in new layer's weight matrix = neurons in previous layers + // = columns in weight matrix of following layer BEFORE mutation + let (following_weights_num_rows, following_weights_num_cols) = following_layer.weights.shape(); + // Next layer weight matrix BEFORE mutation: + // 1 0 + // 2 3 + // 0 4 + // New layer weight matrix: + // 1 0 (1) + // 1 0 (2) + // 0 1 (3) + // 0 1 (4) + // Next layer weight matrix AFTER mutation: + // 1 0 0 0 + // 0 2 3 0 + // 0 0 0 4 + + // rows in new layer's weight matrix = neurons = connections between + // = non-zero entries in weight matrix of following layer BEFORE mutation + let zero_component = T::Component::zero(); + let new_weights_num_rows = following_layer.weights.iter().filter(|&weight| *weight != zero_component).count(); + let mut new_weights = T::zeros((new_weights_num_rows, following_weights_num_cols)); + // modified following weight matrix: + // rows unchanged, columns = rows in new layer's weight matrix + let mut following_weights = T::zeros((following_weights_num_rows, new_weights_num_rows)); + let mut connection_idx = 0; + following_layer.weights.indexed_iter().for_each( + |((row_idx, col_idx), weight)| if *weight != zero_component { + new_weights[[connection_idx, col_idx]] = T::Component::one(); + following_weights[[row_idx, connection_idx]] = *weight; + connection_idx += 1; + } + ); + following_layer.weights = following_weights; + let new_biases = T::zeros((new_weights_num_rows, 1)); + let new_layer = Layer::new(new_weights, new_biases, following_layer.activation.clone()); + layers.insert(new_layer_idx, new_layer); +} + + +fn should_add_layer(num_layers: usize, rng: &mut ThreadRng) -> bool { + let threshold = 1.0/2.0_f32.powi((num_layers - 1) as i32); + rng.gen_range(0.0..1.0) < threshold +} + + +/// Randomly adds new connections between existing neurons. +pub fn mutate_add_connections(layers: &mut Vec>, rng: &mut ThreadRng) { + let neuron_index_lookup = get_neuron_index_lookup(layers); + let total_size = neuron_index_lookup.len(); + // Decide how many new connections to create. + let dist = Poisson::new(1.).unwrap(); + let num_new_connections = dist.sample(rng) as usize; + (0..num_new_connections).for_each(|_| add_new_connection(layers, total_size, &neuron_index_lookup, rng)); +} + + +/// Returns a vector that maps "global" neuron index to 2-tuple of layer index and "in-layer" neuron index. +fn get_neuron_index_lookup(layers: &Vec>) -> Vec<(usize, usize)> { + let size = layers.iter().map(|layer| layer.size()).sum(); + let mut neuron_index_lookup = Vec::with_capacity(size); + for (layer_idx, layer) in layers.iter().enumerate() { + for neuron_idx in 0..layer.size() { + neuron_index_lookup.push((layer_idx, neuron_idx)) + } + } + neuron_index_lookup +} + + +/// Randomly chooses a start and end of a new connection. +fn choose_new_connection( + layers: &Vec>, + total_size: usize, + neuron_index_lookup: &Vec<(usize, usize)>, + rng: &mut ThreadRng, +) -> (usize, usize, usize, usize) { + let global_idx = rng.gen_range(0..total_size); + let (mut start_layer_idx, mut start_neuron_idx) = neuron_index_lookup[global_idx]; + let mut distribution_weights: Vec = layers.iter().map(|layer| layer.size() as f32).collect(); + for index in 0..distribution_weights.len() { + let distance = index as isize - start_layer_idx as isize; // cast could theoretically lead to wraparound if there is an absurd amount of layers + distribution_weights[index] *= if distance == 0 { 0. } else { 2f32.powi((1-distance.abs()) as i32) }; + } + let dist = WeightedIndex::new(&distribution_weights).unwrap(); + let mut end_layer_idx = dist.sample(rng); + let mut end_neuron_idx = rng.gen_range(0..layers[end_layer_idx].size()); + if end_layer_idx < start_layer_idx { + (end_layer_idx, end_neuron_idx, start_layer_idx, start_neuron_idx) = (start_layer_idx, start_neuron_idx, end_layer_idx, end_neuron_idx); + } + (start_layer_idx, start_neuron_idx, end_layer_idx, end_neuron_idx) +} + + +/// Creates new connections and possibly neurons in between. +fn add_new_connection( + layers: &mut Vec>, + total_size: usize, + neuron_index_lookup: &Vec<(usize, usize)>, + rng: &mut ThreadRng, +) { + let (start_layer_idx, start_neuron_idx, end_layer_idx, end_neuron_idx) = choose_new_connection(layers, total_size, neuron_index_lookup, rng); + let mut prev_neuron_idx = start_neuron_idx; + for layer_idx in (start_layer_idx + 1)..end_layer_idx { + // Append a row to weight matrix of intermediate layer. + // The row is all zeros except for the column corresponding to the neuron connected to + // it from the previous layer. + let weights = &mut layers[layer_idx].weights; + let mut new_row = vec![T::Component::zero(); weights.shape().1]; + new_row[prev_neuron_idx] = T::Component::one(); + weights.append_row(new_row.as_slice()); + // Remember the index of the new neuron for the next iteration. + prev_neuron_idx = weights.shape().0 - 1; + // Append a zero to the bias vector. + layers[layer_idx].biases.append_row(vec![T::Component::zero()].as_slice()); + // Append a column to the next layer's weight matrix with all zeros. + let weights = &mut layers[layer_idx + 1].weights; + weights.append_column(vec![T::Component::zero(); weights.shape().0].as_slice()); + } + layers[end_layer_idx].weights[[end_neuron_idx, prev_neuron_idx]] = init_weight(rng, false).unwrap(); +} + + +pub fn get_species(individual: &Individual) -> String { + individual.num_layers().to_string() +} + + +#[cfg(test)] +mod tests { + use ndarray::array; + + use super::*; + + #[test] + fn test_new() { + let _pop = World::new( + select, + procreate, + get_species, + (array![[0.]], array![[0.]]), + ); + } +}