diff --git a/GA.jl b/GA.jl
new file mode 100644
index 0000000..d220a80
--- /dev/null
+++ b/GA.jl
@@ -0,0 +1,115 @@
+#=
+Parents <— {randomly generated population}
+While not (termination criterion)
+    Calculate the fitness of each parent in the population
+    Children <- 0
+    While | Children | < | Parents |
+        Use fitnesses to probabilistically select a pair of parents for mating
+        Mate the parents to create children c\ and c<i
+        Children <— Children U {ci,C2}
+    Loop
+    Randomly mutate some of the children
+    Parents <— Children
+Next generation
+
+Use continuous GA, not binary
+
+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
+    Values will be from lb to ub (lower to upper bound)
+    """
+
+    # Check to see if we dont have enough bounds or have too many
+    if length(bounds) != chromosomeSize
+        println("Length of bounds ($(length(bounds))) not equal to specified chromosome size ($chromosomeSize)")
+        exit()
+    end
+
+    # Init population
+    population = Array{Float64}(0,chromosomeSize+1)
+
+    # Foreach member in population
+    for i = 1:populationSize
+        # Init the member
+        row = []
+        # Create member elements based on bounds
+        for bound in bounds
+            push!(row, rand(bound[1]:0.10:bound[2]))
+        end
+        # Add another element to the end to keep score
+        push!(row, 1.0)
+
+        # hcat row, append to 2D population matrix
+        population = cat(1, population, hcat(row...))
+    end
+
+    return population
+
+end
+
+function score_population(population::Array{Float64}, fitness_function::Function, populationSize::Int, chromosomeSize::Int)
+    """
+    Args
+        population: The population to score. Must be 2D array of Float64 elements
+        fitness_function: Function to score the population
+                          Must be able to take array as argument
+                          For example f(x,y) should be f(A) where A[1]=x, A[2]=y
+    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)
+    end
+    println(population)
+    return population
+end
+
+function GA(populationSize::Int, chromosomeSize::Int, fitness_function::Function, bounds::Array{Tuple{Float64, Float64},1})
+    """
+    Args
+        populationSize: Total number of chromosomes
+        chromosomeSize: How long each chromosome should be (variables in fitness function)
+        fitness_function: Function to determine fitness of each solution
+                          Should take an array as an arg
+                          Example: f(x,y,z) = x+y+z
+                           should be f(X) = X[1]+X[2]+X[3]
+                          in order to allow GA to take in any function with any
+                          amount of args to the fitness function
+        bounds: 1D array of tuples, each tuple is the (lower, upper) bounds
+                for each variable. Both lower and upper need to be Float64 types
+                Length should match chromosomeSize
+                e.g: [(1.0,2.0), (3.5,4.5)]
+    """
+
+    # Set recombination and mutation rate, lower and upper bound
+    recombRate = 0.7
+    mutateRate = 0.05
+    maxIterations = 1
+    # First initialize the population
+    population = initialize_population(populationSize, chromosomeSize, bounds)
+    i = 0
+    while i < maxIterations
+        scoredPopulation = score_population(population, fitness_function, populationSize, chromosomeSize)
+
+        i += 1
+
+    end
+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))
+end
+bounds = [(-2.0,2.0), (-2.0,2.0)]
+GA(25,2,ackley,bounds)