GPUGA.jl

#setup GPUArrays and CUDA backend
#using CUDAnative
using GPUArrays
CUBackend.init()

#overload for simple julia side work
import Base: <,>,==

#constants used in evolution
const absMax = 5.0f0
const D = 14 #Can't go any larger than this! It changes to an allocated version
const PS = 1024#Max batch size 
const gen = 1000
const T_SIZE = 20
const MR = 0.3f0
const MA = 0.2f0
const XR = 0.5f0

#Random number generator
function TausStep(z::Unsigned, S1::Integer, S2::Integer, S3::Integer, M::Unsigned)
    b = (((z << S1) ⊻ z) >> S2)
    return (((z & M) << S3) ⊻ b)
end

LCGStep(z::Unsigned, A::Unsigned, C::Unsigned) = A * z + C

function Random(state::NTuple{4, T}) where T

    s1 = TausStep(state[1], Cint(13), Cint(19), Cint(12), T(4294967294))
    s2 = TausStep(state[2], Cint(2), Cint(25), Cint(4), T(4294967288))
    s3 = TausStep(state[3], Cint(3), Cint(11), Cint(17), T(4294967280))
    s4 = LCGStep(state[4], T(1664525), T(1013904223))
    
    state = (s1,s2,s3,s4)

    return (state,
        Float32(2.3283064f-10 * (state[1] ⊻ state[2] ⊻ state[3] ⊻ state[4]))
    )
end

function gpu_rand(states, tstate)
    threadid = tstate
    states[threadid], result = Random(states[threadid])
    return result
end

#random ntuple generator
#generates ntuples with random float values between 0 and 1
@generated function ntup_rand(::Val{N}, states, tstate) where N
    expr = Expr(:tuple)
    for i = 1:N
        push!(expr.args, :(gpu_rand(states, tstate)))
    end
    expr
end

#generates ntuples which are true < step and false after
@generated function ntup_step(::Val{N}, step) where N
    expr = Expr(:tuple)
    for i = 1:N
        push!(expr.args, :($i <= step))
    end
    expr
end

#immutable individual
struct Individual
    chromosome::NTuple{D, Float32}
    fitness::Float32
    function Individual(size::Int, fitness)
        chromosome = NTuple{D, Float32}(rand(Float32, D) * absMax * 2f0 - absMax)
        return new(chromosome, fitness(new(chromosome,0.0f0)))
    end
    Individual(data, fitness::Float32) = new(data, fitness)
end

<(x::Individual, y::Individual) = x.fitness < y.fitness
>(x::Individual, y::Individual) = x.fitness > y.fitness
==(x::Individual, y::Individual) = x.fitness == y.fitness

#fitness function
function sphereFitness(input::Individual)
    total::Float32 = 0.0f0
    for gene in input.chromosome
        total += gene*gene
    end
    return total
end

#tournament selection
function select(pop, states, tstate)
    pos = round(Int32, gpu_rand(states, tstate)*(PS-1))+1
    bestitem = pop[pos]
    for i = 2:T_SIZE
        pos = round(Int32, gpu_rand(states, tstate)*(PS-1))+1
        item = pop[pos]
        if bestitem > item
            bestitem = item
        end
    end
    return bestitem
end

#search function
function best(pop)
    best = pop[1]
    for item in pop
        if best > item
            best = item
        end
    end
    return best
end

function genIndividuals(state, pop, newpop, states)
    tstate = linear_index(newpop, state)
    #elitism
    if tstate == 1
        newpop[1] = best(pop)
        return
    end

    #crossover
    doxover = gpu_rand(states, tstate) < XR
    newChromosome = select(pop, states, tstate).chromosome
    if doxover
        xpoint = round(Int32, gpu_rand(states, tstate) * D)
        choices = ntup_step(Val{D}(), xpoint) 
        p2 = select(pop, states, tstate).chromosome
        newChromosome = ifelse.(choices, newChromosome, p2)
    end

    #alternate crossover
    # if doxover
    #     choices = ntup_rand(Val{D}(), states, tstate) .< 0.5f0
    #     p2 = select(pop, states, tstate).chromosome
    #     newChromosome = ifelse.(choices, newChromosome, p2)        
    # end

    #mutation
    domutate = ntup_rand(Val{D}(), states, tstate) .< MR
    mutation = ntup_rand(Val{D}(), states, tstate) .* 2 .* MA .- MA
    newChromosome = newChromosome .+ domutate.*mutation

    #output assignment
    newInd = Individual(newChromosome, 0.0f0)
    newInd = Individual(newChromosome, sphereFitness(newInd))
    newpop[tstate] = newInd
    return
end

function main()
    if gen %2 != 0
        println("Even generations only")
        return
    end
    pop = GPUArray(Individual[Individual(D, sphereFitness) for i = 1:PS])
    oldpop = GPUArray(Individual[Individual(D, sphereFitness) for i = 1:PS])
    states = GPUArray(NTuple{4, UInt32}[(rand(UInt32), rand(UInt32), rand(UInt32), rand(UInt32)) for i = 1:PS])
    printpop = [pop;]
    println(best(pop).fitness)
    for i = 1:div(gen,2)
        gpu_call(genIndividuals, pop, (pop, oldpop, states))
        gpu_call(genIndividuals, oldpop, (oldpop, pop, states))
        # println(best(pop).fitness)
    end
    pop = [pop;]
    println(best(pop))    
    println(best(pop).fitness)
    return pop
end

@time main()