diff --git a/GA.jl b/GA.jl
index d220a80..dacc1f4 100644
--- a/GA.jl
+++ b/GA.jl
@@ -1,3 +1,4 @@
+using StatsBase
 #=
 Parents <— {randomly generated population}
 While not (termination criterion)
@@ -20,7 +21,7 @@ include array of independent variables along with associated range
 
 function initialize_population(populationSize::Int, chromosomeSize::Int, bounds::Array{Tuple{Float64, Float64},1})
     """
-    Initializes the population based on the given population and chromosome size
+    Initializes the population based on the given population size and chromosome size
     Values will be from lb to ub (lower to upper bound)
     """
 
@@ -31,7 +32,7 @@ function initialize_population(populationSize::Int, chromosomeSize::Int, bounds:
     end
 
     # Init population
-    population = Array{Float64}(0,chromosomeSize+1)
+    population = Array{Array{Float64}}(0)
 
     # Foreach member in population
     for i = 1:populationSize
@@ -45,14 +46,13 @@ function initialize_population(populationSize::Int, chromosomeSize::Int, bounds:
         push!(row, 1.0)
 
         # hcat row, append to 2D population matrix
-        population = cat(1, population, hcat(row...))
+        population = push!(population, row)
     end
-
     return population
 
 end
 
-function score_population(population::Array{Float64}, fitness_function::Function, populationSize::Int, chromosomeSize::Int)
+function score_sort_population(population::Array{Array{Float64}}, fitness_function::Function)
     """
     Args
         population: The population to score. Must be 2D array of Float64 elements
@@ -62,17 +62,71 @@ function score_population(population::Array{Float64}, fitness_function::Function
     Returns scored population. Score is the last element of each chromosome (row)
     in population
     """
-    #= 1-3, 4-6, 7-9, 10-12 
-        1    2    3     4
-    =#
-    population = population'
-    for i in 0:populationSize-1
-        println(i)
-        member = population[(chromosomeSize+1)*i+1:(chromosomeSize+1)*(i+1)]
-        println(member)
-        population[(chromosomeSize+1)*(i+1)] = fitness_function(member)
+    # Enumerate through population, set last element of each chromosome to score
+    for (index,member) in enumerate(population)
+        population[index][end] = fitness_function(member)
     end
-    println(population)
+    sort!(population, by = x -> x[end])
+    return population
+end
+
+function create_next_generation(population::Array{Array{Float64}}, recombRate::Float64)
+    """
+    Creates children given a population and recombination rate
+    Uses an elite selection of 0.10 (rounded up)
+    Creates random chance
+    TODO: Implement ranking
+    """
+
+    popSize = length(population)
+    memberSize = length(population[1])
+    # Init children array
+    children = Array{Array{Float64}}(0)
+    # Elite selection
+    for i in 1:Int(ceil(popSize*0.10))
+        push!(children, population[i])
+    end
+    
+    # Generate rest of population
+    while length(children) < popSize
+        # Two random children
+        choices = sample(1:popSize, 2, replace=false)
+        # Check for crossover
+        if rand() < recombRate
+            # Choose xover point
+            # Use -2 since last one is score, and last actual element cant be crossedover since that means there is none
+            xover = rand(1:memberSize-2)
+            child1 = cat(1,population[choices[1]][1:xover],population[choices[2]][xover+1:end])
+            child2 = cat(1,population[choices[2]][1:xover],population[choices[1]][xover+1:end])
+        else
+            # Else no crossover, just copy
+            child1 = population[choices[1]]
+            child2 = population[choices[2]]
+        end
+        # Push child 1 and 2 to children array
+        push!(children, child1, child2)
+    end
+
+    # We might have one extra due to rounding in the elite selection, let's check
+    if length(children) != popSize
+        pop!(children)
+    end
+    return children
+end
+
+function mutate(population::Array{Array{Float64}}, mutateRate::Float64, bounds::Array{Tuple{Float64, Float64},1})
+    """
+    Iterates through each chromosome and mutates if needed
+    """
+
+    for member in population
+        for i = 1:length(member)-1
+            if rand() < mutateRate
+                member[i] = rand(bounds[i][1]:0.10:bounds[i][2])
+            end
+        end
+    end
+
     return population
 end
 
@@ -96,20 +150,28 @@ function GA(populationSize::Int, chromosomeSize::Int, fitness_function::Function
     # Set recombination and mutation rate, lower and upper bound
     recombRate = 0.7
     mutateRate = 0.05
-    maxIterations = 1
+    maxIterations = 1000
     # First initialize the population
     population = initialize_population(populationSize, chromosomeSize, bounds)
+
+    # Then loop for the required iterations
+    bestScores = Array{Float64}(0)
     i = 0
     while i < maxIterations
-        scoredPopulation = score_population(population, fitness_function, populationSize, chromosomeSize)
-
+        # Score and sort the population
+        population = score_sort_population(population, fitness_function)
+        push!(bestScores,population[1][end])
+        population = create_next_generation(population, recombRate)
+        population = mutate(population, mutateRate, bounds)
         i += 1
-
     end
+    println(bestScores)
+    println(population)
 end
 
 function ackley(X)
-    return e - 20*exp(-0.2*sqrt((X[1]^2 + X[2]^2)/2) - exp((cos(2*pi*X[1]) + cos(2*pi*X[2]))/2))
+    return e - 20*exp(-0.2*sqrt((X[1]^2 + X[2]^2)/2)) - exp((cos(2*pi*X[1]) + cos(2*pi*X[2]))/2)
 end
+
 bounds = [(-2.0,2.0), (-2.0,2.0)]
 GA(25,2,ackley,bounds)