%% Assignment 04
clear all;
clc;

%% Initial parameters
max_generations=300;
generation=1;
target_fitness = 1000;
num_i_nodes = 6;
num_b_nodes = 1;
num_o_nodes = 1;
num_h_nodes = 0;
pop_size = 150;
fixed_topology=1;
rand_h_nodes=2;

stagnation.threshold=1e-2;
stagnation.number_generation=15;
refocus.threshold=1e-2;
refocus.number_generation=20;

initial.kill_percentage=0.2;
initial.number_for_kill=5;
initial.number_copy=5;

selection.pressure=2;

% crossover.percentage=0.8;
crossover.percentage=0;
crossover.rate=0.1;
crossover.probability_interspecies=0.0;
crossover.probability_multipoint=0.6;

% Mutation related parameters
mutation.rate = 0.2;
mutation.perturbation_rate = 0.9;
mutation.perturbation_range=0.1;
mutation.cap = 5;
mutation.range = 2;
mutation.mutate_weight=0.9;
mutation.add_node_rate=0.03;
mutation.add_connection_rate=0.05;
mutation.recurrency_rate=0.0;
mutation.reenabled_rate=0.25;

speciation.c1 = 1;
speciation.c2 = 1;
speciation.c3 = 4;
speciation.threshold=8;

species_record(1).ID=0;
species_record(1).number_individuals=0;
species_record(1).generation_record=[];

if fixed_topology == 1
    num_nodes = num_i_nodes+num_b_nodes+num_o_nodes+num_h_nodes;
    num_conections = num_i_nodes+num_b_nodes;
    connections_from = 1:num_i_nodes+num_b_nodes;
    connections_to = num_nodes*ones(1,num_conections);
    global_innovation_number=num_conections+1;
    innovation_nums = 1:num_conections;
else
    num_h_nodes = randi([0,rand_h_nodes]);
end

innovation_record = [];
max_overall_fitness = [];
%% Initialize population
for ind=1:pop_size

    if fixed_topology ~= 1
        num_nodes = num_i_nodes+num_b_nodes+num_o_nodes+num_h_nodes;
        connections_from = [];
        connections_to = [];
        innovation_nums = [];
        %Initialize random connections for each individual
        conection_matrix=round(rand(num_nodes,num_nodes));
        num_conections=sum(sum(conection_matrix));
        global_innovation_number=num_conections+1;
        for ind1=1:num_nodes
            for ind2=1:num_nodes
                enabled_connection = conection_matrix(ind1,ind2);
                if enabled_connection
                    innovation_nums = [innovation_nums ((ind1-1)*num_nodes + ind2)];
                    connections_from = [connections_from ind1];
                    connections_to = [connections_to ind2];
                end
            end
        end
    end

    %node ID's
    %node type: 1=input 2=output 3=hidden 4=bias
    %node input state
    %node output state
    population(ind).node_genes=[1:(num_i_nodes+num_b_nodes+num_o_nodes+num_h_nodes);
        ones(1,num_i_nodes),4*ones(1,num_b_nodes),2*ones(1,num_o_nodes),3*ones(1,num_h_nodes);
        zeros(1,num_i_nodes),ones(1,num_b_nodes),zeros(1,num_o_nodes),zeros(1,num_h_nodes);
        zeros(1,num_nodes)];

    %innovation number
    %connection from,
    %connection to
    %weights (uniformly distributed in [-2 +2], all connections enabled)
    %enable bit
    population(ind).connection_genes=[innovation_nums;
        connections_from;
        connections_to;
        rand(1,num_conections)*2*mutation.range-mutation.range;
        ones(1,num_conections)];
    population(ind).wMat = weight_matrix(population(ind).node_genes,population(ind).connection_genes);
    population(ind).fitness=0;
    population(ind).species=0;
end
%% Innovation reccord for initial population
%Innovation ID
%Connections from
%Connection to
%New node(Highest)
%Generation this innovation occured
innovation_record=[population(ind).connection_genes(1:3,:);zeros(size(population(ind).connection_genes(1:2,:)))];
innovation_record(4,size(innovation_record,2))=max(population(1).node_genes(1,:));

%% Initial speciation
population(1).species=1;
ref_individuals = 1;
matrix_reference_individuals(ref_individuals).connection_genes=population(1).connection_genes;
species_record(1).ID=1;
species_record(1).number_individuals=1;

for index_individual=2:size(population,2)
    assigned_existing_species_flag=0;
    new_species_flag=0;
    ind=1;
    %loops through the existing species, terminates when either the
    %individual is assigned to existing species or there are no more
    %species to test it against, which means it is a new species
    while assigned_existing_species_flag==0 & new_species_flag==0
        distance = compatibility_distance(population(index_individual), ...
            matrix_reference_individuals(ind), speciation);
        if distance<speciation.threshold %If within threshold, assign to the existing species
            population(index_individual).species=ind;
            assigned_existing_species_flag=1;
            species_record(ind).number_individuals=species_record(ind).number_individuals+1;
        end
        ind=ind+1;
        %Outside of species references, must be new species
        if ind>size(matrix_reference_individuals(:),1) & assigned_existing_species_flag==0
            new_species_flag=1;
        end
    end
    %check for new species, if it is, update the species_record and use individual as reference for new species
    if new_species_flag==1
        population(index_individual).species=ind;
        ref_individuals = ref_individuals + 1;
        matrix_reference_individuals(ref_individuals).connection_genes=population(index_individual).connection_genes;
        species_record(ind).ID=ind;
        %if number individuals in a species is zero, that species is extinct
        species_record(ind).number_individuals=1;
    end
end

%% Start generation loop
maximal_fitness=0;
flag_solution=0;
while generation<=max_generations & flag_solution==0
    max_fitnesses_current_generation=zeros(1,size(species_record,2));

    % Fitness of population

    population = get_pop_fitness(population, @RNNet, target_fitness);

    for ind=1:size(species_record,2)
        if species_record(ind).number_individuals>0
            [max_fitness,index_individual_max]=max(([population(:).species]==ind).*[population(:).fitness]);
            mean_fitness=sum(([population(:).species]==ind).*[population(:).fitness])/species_record(ind).number_individuals;
            % Compute stagnation vector (last stagnation.number_generation-1 max fitnesses plus current fitness
            if size(species_record(ind).generation_record,2)>stagnation.number_generation-2
                stagnation_vector=[species_record(ind).generation_record(3,size(species_record(ind).generation_record,2)-stagnation.number_generation+2:size(species_record(ind).generation_record,2)),max_fitness];
                if sum(abs(stagnation_vector-mean(stagnation_vector))<stagnation.threshold)==stagnation.number_generation %Check for stagnation
                    mean_fitness=0.01; %set mean fitness to small value to eliminate species (cannot be set to 0, if only one species is present, we would have divide by zero in fitness sharing. anyways, with only one species present, we have to keep it)
                end
            end
            species_record(ind).generation_record=[species_record(ind).generation_record,[generation;mean_fitness;max_fitness;index_individual_max]];
            max_fitnesses_current_generation(1,ind)=max_fitness;
        end
    end

    c=[];
    for index_species=1:size(species_record,2)
        c=[c,species_record(index_species).generation_record(1:4,size(species_record(index_species).generation_record,2))];
    end
    max_overall_fitness=[max_overall_fitness,[max(c(3,:).*(c(1,:)==generation));generation]];
    maximal_fitness=max(c(3,:).*(c(1,:)==generation));
    if maximal_fitness>=target_fitness
        flag_solution = 1;
    else
        % Function reproduce goes here
        [population,species_record,innovation_record]=...
            reproduce_pop(population, species_record, innovation_record, initial, selection, ...
            crossover, mutation, speciation, generation, pop_size, @RNNet, target_fitness);
    end
    generation
    generation=generation+1;
end

[max_fitness, max_fit_individual] = max([population(:).fitness]);
figure(1);
plot(digraph(population(max_fit_individual).wMat));
% simulate = twoPole_test( population(max_fit_individual).wMat, @RNNet, target_fitness,'vis');


%% Function implementing crossover, mutation, and speciation
function [new_population,updated_species_record,updated_innovation_record]=...
    reproduce_pop(population, species_record, innovation_record, initial, selection, ...
    crossover, mutation, speciation, generation, pop_size, RNNet, target_fitness)

% new_population(1) = population(1);
new_population = rank_based_selection(population,species_record);

% Sum of avg from every species
tot_avg_specie_fitness=0;
for species_ind=1:size(species_record,2)
    tot_avg_specie_fitness=tot_avg_specie_fitness+species_record(species_ind).generation_record(2,size(species_record(species_ind).generation_record,2))...
        *(species_record(species_ind).number_individuals>0);
end

% Matrix containing iformation about species that currently exist and
% will be propagating.
propagating_species=[];
overflow=0;
index_individual=0;
for species_ind = 1:size(species_record,2)
    if species_record(species_ind).number_individuals>0
        % Calculate # offsprings for each species
        avg_specie_fitness = species_record(species_ind).generation_record(2,size(species_record(species_ind).generation_record,2));
        offsprings = (avg_specie_fitness/tot_avg_specie_fitness)*pop_size;
        overflow=overflow+offsprings-floor(offsprings);
        if overflow>=1
            offsprings=ceil(offsprings);
            overflow=overflow-1;
        else
            offsprings=floor(offsprings);
        end

        if offsprings>0
            propagating_species=[propagating_species,[species_record(species_ind).ID;offsprings;0]];
            if species_record(species_ind).number_individuals>=initial.number_copy %check for condition for objective 2
               index_individual=index_individual+1;
               new_population(index_individual)=population(species_record(species_ind).generation_record(4,size(species_record(species_ind).generation_record,2))); % Objective 2
               propagating_species(3,size(propagating_species,2))=1; %Update matrix_existing_and_propagating_species
            end
        end

    end
end
% propagating_species
% Reference individuals from every species
index_ref=0;
for index_species_ref=1:size(species_record,2)
    if sum([population(:).species]==index_species_ref)>0
        index_ref=index_ref+1;
        [discard,index_ref_old]=max([population(:).species]==index_species_ref);
        population_ref(index_ref)=population(index_ref_old);
    end
end

pop_ind = 0;
for index_species=1:size(propagating_species,2)
    count_individuals_species=0;
    species_ID=propagating_species(1,index_species);

    fitnesses_species=[population(find([population(:).species]==species_ID)).fitness];
    index_fitnesses_species=find([population(:).species]==species_ID);
    [discard,sorted_fitnesses]=sort(fitnesses_species);
    ranking=zeros(1,size(fitnesses_species,2));
    ranking(sorted_fitnesses)=1:size(fitnesses_species,2);
    if size(fitnesses_species,2)>1
        FitnV=(2-selection.pressure+2*(selection.pressure-1)/(size(fitnesses_species,2)-1)*(ranking-1))';
    else
        FitnV=2;
    end

    number_overall=propagating_species(2,index_species)-propagating_species(3,index_species);
    number_crossover=round(crossover.percentage*number_overall);
    number_mutate=number_overall-number_crossover;
    Nind=size(fitnesses_species,2);
    Nsel=2*number_crossover+number_mutate;

    if Nsel==0
       Nsel=1;
    end

    cumfit = cumsum(FitnV);
    trials = cumfit(Nind) / Nsel * (rand + (0:Nsel-1)');
    Mf = cumfit(:, ones(1, Nsel));
    Mt = trials(:, ones(1, Nind))';
    [NewChrIx, ans] = find(Mt < Mf & [ zeros(1, Nsel); Mf(1:Nind-1, :) ] <= Mt);
    [ans, shuf] = sort(rand(Nsel, 1));
    NewChrIx = NewChrIx(shuf);
    NewChrIx = index_fitnesses_species(NewChrIx);



%     species_individuals = population(find([population(:).species]==propagating_species(1,index_species)));
    while propagating_species(3,find([propagating_species(1,:)]==species_ID)) < ...
            propagating_species(2,find([propagating_species(1,:)]==species_ID))

        pop_ind=pop_ind+1;
        count_individuals_species=count_individuals_species+1;
        if rand < crossover.rate
            parent1 = population(NewChrIx(randi(size(NewChrIx,2))));
            parent2 = population(NewChrIx(randi(size(NewChrIx,2))));
            parent_nodes = [parent1.node_genes(:,:) parent2.node_genes(:,:)];
            parent_connection_genes = [parent1.connection_genes(:,:) parent2.connection_genes(:,:)];

            [unique_nodes, u_n_ind, u_n_p_ind] = ...
                unique(parent_nodes(1,:));
            [unique_innovation_num, u_i_ind, u_i_p_ind]  = ...
                unique(parent_connection_genes(1,:));
            offspring.node_genes = [];
            offspring.connection_genes = [];
            for node_id = 1:size(unique_nodes,2)
                offspring.node_genes = [offspring.node_genes ...
                    parent_nodes(:,u_n_ind(node_id))];
            end
            for innov_id = 1:size(unique_innovation_num,2)
                offspring.connection_genes = [offspring.connection_genes ...
                    parent_connection_genes(:,u_i_ind(innov_id))];
            end
            offspring.wMat = weight_matrix(offspring.node_genes,offspring.connection_genes);
            offspring.fitness = get_ind_fitness(offspring, RNNet, target_fitness);
        else
            offspring = population(NewChrIx(count_individuals_species));
        end

        % Weight value mutation
        if rand < mutation.mutate_weight
            offspring = weight_value_mutation(offspring,mutation);
        end

        % Add connection mutation
        if rand < mutation.add_node_rate
            [offspring,innovation_record] = ...
                add_connection(offspring,innovation_record,generation,mutation.range);
        end

        % Add node mutation
        if rand < mutation.add_connection_rate
            [offspring,innovation_record] = ...
                add_node(offspring,innovation_record,generation,mutation.range);

        end

        % Assing a species to offspring
        [offspring, species_record] = ...
            speciate_individual(offspring,population_ref,species_record,speciation,generation,pop_ind);

        offspring.fitness = get_ind_fitness(offspring, RNNet, target_fitness);

        new_population(pop_ind) = offspring;
        propagating_species(3,index_species) = propagating_species(3,index_species) + 1;
    end
end

for index_species=1:size(species_record,2)
   species_record(index_species).number_individuals=sum([new_population(:).species]==index_species);
end

updated_species_record = species_record;
updated_innovation_record = innovation_record;
end

%% Rank based selection
function [new_population] = rank_based_selection(population,species_record)
pop_ind = 1;
new_population = population;
for species_ID=1:size(species_record,2)
    species_individuals = population(find([population(:).species]==species_ID));
    if size(species_individuals,2)<0
        num_ind_species = species_record(species_ID).number_individuals;
        competitors = randi([1,size(species_individuals,2)],2,num_ind_species);
        for species_ind=1:num_ind_species
            species_ind
            fitness_1 = species_individuals(competitors(1,species_ind)).fitness
            fitness_2 = species_individuals(competitors(2,species_ind)).fitness;

            if fitness_1 >= fitness_2
                new_population(pop_ind) = species_individuals(competitors(1,species_ind));
            else
                new_population(pop_ind) = species_individuals(competitors(2,species_ind));
            end
            pop_ind = pop_ind + 1;
        end
    end
end
end

%% Weight value mutation
function individual = weight_value_mutation(individual,mutation)
num_connections = size(individual.connection_genes(1,:),2);
mutate_weights = rand(1,num_connections);
perturbate_weight = rand(1,num_connections);
for con_id=1:num_connections
    node_from = find(individual.node_genes(1,:)==individual.connection_genes(2,con_id));
    node_to = find(individual.node_genes(1,:)==individual.connection_genes(3,con_id));
    if mutate_weights(con_id) < mutation.rate
        if perturbate_weight(con_id) < mutation.perturbation_rate
            perturbation = rand(1)*2*mutation.perturbation_range-...
                mutation.perturbation_range;
            new_weight = individual.connection_genes(4,con_id) + perturbation;
            if new_weight < -mutation.cap
                new_weight = -mutation.cap;
            elseif new_weight > mutation.cap
                new_weight = mutation.cap;
            end
            individual.wMat(node_from,node_to) = new_weight;
            individual.connection_genes(4,con_id) = new_weight;
        else
            new_weight = rand(1)*2*mutation.range-mutation.range;
            individual.wMat(node_from,node_to) =  new_weight;
            individual.connection_genes(4,con_id) = new_weight;
        end
    end
end
end

%% Call simulator for every individual in population to get fitness
function population = get_pop_fitness(population, RNNet, target_fitness)
for ind=1:size(population(:),1)
    population(ind).fitness = twoPole_test( population(ind).wMat, RNNet, target_fitness);%,'vis')
end
end

%% Call simulator for one individual to get fitness
function steps = get_ind_fitness(individual, RNNet, target_fitness)
% individual.wMat = weight_matrix(individual.node_genes, ...
%     individual.connection_genes)
steps = twoPole_test(individual.wMat, RNNet, target_fitness);%,'vis')
end

%% Built the weight matrix from connected genes
function wMat = weight_matrix(node_genes,connection_genes)
wMat = zeros(size(node_genes(1,:),2),size(node_genes(1,:),2));
for c = 1:size(connection_genes(1,:),2)
    enabled = connection_genes(5,c);
    from_node = find(node_genes(1,:)==connection_genes(2,c));
    to_node = find(node_genes(1,:)==connection_genes(3,c));
    weight=connection_genes(4,c);
    if enabled
        wMat(from_node, to_node) = weight;
    else
        wMat(from_node, to_node) = 0;
    end
end
end

%% Add a connection to the network
function [individual,innovation_record] = ...
    add_connection(individual,innovation_record,generation, mutation_range)

% individual.connection_genes
%get all possible connection from enabled ones
inn_old=max((innovation_record(5,:)<generation).*innovation_record(1,:));

pos_conn=individual.connection_genes(2:3,(individual.connection_genes(5,:)==1)  & (individual.connection_genes(1,:)<=inn_old));

%choose a random 'from node' and 'to node' connection

rand_node_connection=pos_conn(:,round(rand*size(pos_conn,2)+0.5));
%
% rand_from_node = rand_node_connection(1);
% rand_to_node = rand_node_connection(2);
rand_from_node = individual.node_genes(1,round(rand*size(individual.node_genes,2)+0.5));
rand_to_node = individual.node_genes(1,round(rand*size(individual.node_genes,2)+0.5));

%check innovation duplicacy
current_gen_changes = innovation_record(:,innovation_record(5,:)==generation);
dup = 0;
if (~isempty(current_gen_changes))

    for cc=1:size(current_gen_changes,2)
        if current_gen_changes(2,cc) == rand_from_node & current_gen_changes(3,cc) == rand_to_node
            next_inn_num = current_gen_changes(1,cc);

            dup = 1;
        else
            next_inn_num = max(innovation_record(1,:))+1;
        end
    end
else
    next_inn_num = max(innovation_record(1,:))+1;
end

%choose a random 'from node' and 'to node'

weight = rand(1)*2*mutation_range-mutation_range;
new_connection = [next_inn_num; rand_from_node; rand_to_node; weight; 1];
individual.connection_genes = [individual.connection_genes new_connection];
if dup==0
    innovation_record = [innovation_record [next_inn_num rand_from_node rand_to_node 0 generation]'];
end
% innovation_record

individual.wMat = weight_matrix(individual.node_genes,individual.connection_genes);

end

%% Add node to network
function [individual,innovation_record] = ...
    add_node(individual,innovation_record,generation,mutation_range)

num_nodes = size(individual.node_genes(1,:),2);
num_connections = size(individual.connection_genes(1,:),2);
new_node_ID =max(innovation_record(4,:))+1;

%get all possible connection from enabled ones
inn_old=max((innovation_record(5,:)<generation).*innovation_record(1,:));

pos_conn=individual.connection_genes(2:3,(individual.connection_genes(5,:)==1)  & (individual.connection_genes(1,:)<=inn_old));

%choose a random 'from node' and 'to node' connection

rand_node_connection=pos_conn(:,round(rand*size(pos_conn,2)+0.5));
%
rand_from_node = rand_node_connection(1);
rand_to_node = rand_node_connection(2);

%sanity check
% rand_from_node = innovation_record(2,8);
% rand_to_node = innovation_record(3,9);

%check innovation duplicacy
current_gen_changes = innovation_record(:,innovation_record(5,:)==generation);
% current_gen_changes;
dup = 0;
current_gen_changes;
if (~isempty(current_gen_changes))
    next_inn_num = max(innovation_record(1,:))+1;
    for cc=1:size(current_gen_changes,2)-1
        if current_gen_changes(2,cc) == rand_from_node & current_gen_changes(3,cc+1) == rand_to_node
            next_inn_num = current_gen_changes(1,cc);
            new_node_ID =current_gen_changes(3,cc);
            dup = 1;

        end
    end
else
    next_inn_num = max(innovation_record(1,:))+1;
end
% new node gene's properties [assign node ID, type=hidden,
% inp_state = 0, o_state = 0
new_node_gene = [ new_node_ID ; 3; 0; 0];
individual.node_genes(:,num_nodes+1) = new_node_gene;



%check if connection it already exists
exist_from_con = find(individual.connection_genes(2,:)==rand_from_node);

if(isempty(exist_from_con))
    weight1=rand(1)*2*mutation_range-mutation_range;
    weight2=rand(1)*2*mutation_range-mutation_range;
else
    for i=1:length(exist_from_con)
        if(individual.connection_genes(3,exist_from_con(i)) == rand_to_node)
            individual.connection_genes(5,exist_from_con(i))=0;

            weight1=1;
            weight2=individual.connection_genes(4,exist_from_con(i));

            break;
        else
            weight1=rand(1)*2*mutation_range-mutation_range;
            weight2=rand(1)*2*mutation_range-mutation_range;
        end
    end
end
new_connection_gene_1 = [ next_inn_num; rand_from_node;new_node_ID;weight1;1];
new_connection_gene_2 = [ next_inn_num+1; new_node_ID;rand_to_node;weight2;1];

%return incremented global innovation number
individual.connection_genes(:,num_connections+1) = new_connection_gene_1;
individual.connection_genes(:,num_connections+2) = new_connection_gene_2;

individual.wMat = weight_matrix(individual.node_genes,individual.connection_genes);

%update innovation
if dup==0
    innovation_record = [innovation_record  [next_inn_num;rand_from_node;new_node_ID;new_node_ID;generation] [next_inn_num+1;new_node_ID;rand_to_node;0;generation]];
end
innovation_record;
% individual
end

%% Computes compatibility distance, only average weight distance considered
function [distance] = compatibility_distance(individual, reference_individual,speciation)
    weight_diff = 0;
    max_inn_ind=max(individual.connection_genes(1,:));
    max_inn_ref=max(reference_individual.connection_genes(1,:));

    if max_inn_ind < max_inn_ref
        disjoint_ind = ismember(individual.connection_genes(1,:),reference_individual.connection_genes(1,:));
        disjoint_ref = ismember(find(reference_individual.connection_genes(1,:)<max_inn_ind),individual.connection_genes(1,:));
        disjoint = size(find([disjoint_ind disjoint_ref] == 0),2);

        excess = size(find(reference_individual.connection_genes(1,:)>max_inn_ind),2);
    else
        disjoint_ref = ismember(reference_individual.connection_genes(1,:),individual.connection_genes(1,:));
        disjoint_ind = ismember(find(individual.connection_genes(1,:)<max_inn_ref),reference_individual.connection_genes(1,:));
        disjoint = size(find([disjoint_ind disjoint_ref] == 0),2);

        excess = size(find(individual.connection_genes(1,:)>max_inn_ref),2);
    end

    matching_genes = intersect(individual.connection_genes(1,:), reference_individual.connection_genes(1,:));

    for weight_ind = 1:size(matching_genes,2)
        weight1 = individual.connection_genes(4,(individual.connection_genes(1,:)==matching_genes(weight_ind)));
        weight2 = reference_individual.connection_genes(4,(reference_individual.connection_genes(1,:)==matching_genes(weight_ind)));
        weight_diff = weight_diff + abs(weight1 - weight2);
    end
    num_matching_genes = size(matching_genes,2);

    average_weight_difference=weight_diff/num_matching_genes;

    distance=speciation.c1*excess/num_matching_genes+...
        speciation.c2*disjoint/num_matching_genes+...
        speciation.c3*average_weight_difference;
end

%% Function to implement speciation
function [offspring, species_record] = ...
    speciate_individual(offspring,population_ref,species_record,speciation,generation,pop_ind)
    species_assigned=0;
    index_population_ref=0;
    while species_assigned==0 & index_population_ref<size(population_ref,2)
        index_population_ref=index_population_ref+1;
        ref_individual=population_ref(index_population_ref);

        distance = compatibility_distance(offspring, ...
            ref_individual, speciation);

        if distance<speciation.threshold
            % assign individual to same species as current reference individual
            offspring.species=ref_individual.species;
            species_assigned=1; %set flag indicating new_individual has been assigned to species
        end
    end
    % Assign a new species to offspring if it doesn't match any other species
    if species_assigned==0
        new_species_ID=size(species_record,2)+1;
        % assign individual to new species
        offspring.species=new_species_ID;
        % update species_record
        species_record(new_species_ID).ID=new_species_ID;
        species_record(new_species_ID).number_individuals=1;
        species_record(new_species_ID).generation_record=[generation;offspring.fitness;offspring.fitness;pop_ind];
        % update population reference with new species individual
        population_ref(size(population_ref,2)+1)=offspring;
    end
end