
--# Notes
--[[
The GA class provides a fairly simple genetic algorithm

The key variables are
# Size of population - more is generally better, but slower
# Number of iterations - again, more is often better, but GA will eventually get stuck, so there's no point making this unnecessarily high
# Crossover rate - the % of the population replaced in each iteration. If it is high, eg 75%, it will produce quicker results initially, by focussing on the most successful solutions so far, but these may be a dead end, and the downside is that a high crossover drives out a lot of apparently bad solutions that may come good later on, ie it destroys variety, which is crucial to evolution. Conversely, a low crossover means slow convergence towards the final answer.
# Mutation - when a new solution is created, its random numbers come from its parents (see below), but for each of them, an extra random number is generated, and if it is lower than the mutation %, then that random number is regenerated randomly instead of coming from a parent.

Replacement is done by sorting the population by fitness, and replacing the worst X%, where X% is the crossover rate. For each replacement, two 'parents' are chosen using a common technique called 'tournament', choosing three solutions at random, and the best of them is used as one parent, then doing the same to choose the other parent. The reason for this is that if you just used the very best solutions to replace bad ones, you will lose genetic variability, so tournament is a way of accessing a wide range of solutions. while ensuring they are not too bad to start with.

The random numbers from the two parents are combined in one of two ways
# 2 slices - the replacement array is made up of the first N random numbers from one parent, and the remaining random numbers from the other parent,where N is chosen randomly. This technique works best when it matters that random numbers are adjacent to each other (as in this demo, where pairs of random numbers are used together)
# Individual selection (called EitherOr in the code), where each random number in the replacement is chosen from one of the parents randomly.
As implemented, GA chooses a replacement method randomly for each replacement, so it mixes these two methods.

The critical choices are mainly with
# how to use the random numbers produced by GA, to generate practical solutions, and
# how to calculate fitness.


--]]
--# Main
--Main

--NB The Notes tab has explanations and some advice on various aspects!

--This class is set up with a very simple exanple, solving the coefficients of a quadratic equation
-- ie solve for a,b,c in ax^2 + bx + c, given the correct answer for x=1,2,3,4,5
--set a,b,c
local a,b,c=-2,33,4

local ga=nil --is the instance of the GA class
local Best= 0 --this stores the best score so far
local BestSol = {} -- this array stores the best solution so far
local nIter=0 --number of iterations,set below

function setup()
    SetupGA()
end

function SetupGA()
    nIter=500 --number of iterations
    --create GA object
    local numVars = 6 --numVars = number of random numbers to be used (reason below)
    local numPopulation = 5000 --numPopulation = the number of sets of random numbers to create
    --the last parameter is a callback for the function that will calculate fitness
    --(GA doesn't do this, because it is application specific)
    --so we will use a function called Fitness in this tab to do the work
    ga=GA(numVars,numPopulation,Fitness)
    --NB if you are calling from within a class, add an extra parameter, self, like so
    --ga=GA(numVars,numPopulation,Fitness,self)
    --set up a coroutine with the function that will do the iterations
    --coroutines are used so that you get a chance to draw the screen and other stuff between iterations
    --if you are calling from within a class, prefix DoRuns with the classname, eg MyClass.DoRuns
    cos=coroutine.create(DoRuns)
    coroutine.resume(cos)
end

--this function manages the iterations
function DoRuns()
    for i=1,nIter do
        --call the GA to process a generation of solutions
        --the first one just creates random numbers for all of them
        --from then on, GA uses crossover and mutation to replace bad solutions with combinations of good ones

         --Parameters are
        --Iteration number (GA only really needs to know if it's the first one or not)
        --Crossover % (>0,<=100) the percentage of solutions to be replaced
        local Crossover=0.4 --generally 20-80%
        --Mutation % when good solutions are combined to replace bad ones, the random numbers
        --  are taken from both "parents", but a small percentage is chosen randomly - that's what this is
        local Mutation=0.01 --generally 1-2%
         --Sample isn't a classic genetic algorithm component, just something I was trying.
         -- It lets you choose the candidates from a subset of the population. If you
        --set it to 100% it samples the entire population, but if you set it to 60% it samples the best
        --60% of solutions
        local Sample=1.0
        --NB if these variables are not going to change, you could pass them into GA just once, at the top. I
        --put them here because the user can play with the parameters as it runs, so they can change

        --call GA
        --it returns two things
        --tblRand is an array of the random numbers in the best solution so far
        --Best is the fitness value of that solution
        BestSol,Best=GA:DoRun(i,Crossover,Mutation,Sample)
        currIter=i
        --pause so we can draw the screen, use the results whatever
        coroutine.yield()
    end
end

--this function takes a set of random numbers from GA and calculates its fitness
--what this is depends entirely on what you're using this for, and it's worth thinking hard about it
--in this case, we're trying to solve for coefficients that are real numbers of any size, plus and minus
--so we need three candidates that can have a very wide range
--but GA actually only gives us numbers in the range 0-1
--how on earth do we produce such a huge range of numbers, using what GA gives us?
--What I've done below is to combine two random numbers to produce each candidate, like so
-- candidate = (rand1 - 0.5) / (rand2 - 0.5)
--there are probably better ways, but this seems to work quite well
--That's why I asked for 6 random numbers up above, rather than three
function Fitness(t) --t is an array of 6 random numbers
    --create the numbers we need, using the random numbers
    local aa=(t[1] - 0.5)/(t[2] - 0.5)
    local bb=(t[3] - 0.5)/(t[4] - 0.5)
    local cc=(t[5] - 0.5)/(t[6] - 0.5)
    --now test them
    --I've done this very inefficently, calculating the same test results everytime
    --but this keeps it all in one place and makes it easy for you to replace it with your stuff
    --If I had to solve this function over a wider range of x, the variation in results could bias
    --the GA. For example, if I just calculate the squared error as I do below, and one of the
    --test answers is 3 and another is 33333, the first error is likely to be dominated by the second one
    --so the GA will do all it can to fix the second one, whereas the first one might be just as important
    --So in my own implementation of this, I dampened the results by dividing the answers by a carefully
    --chosen set of factors that reduced the differences between them, so the GA didn't ignore any of them
    local error=0
    for x=1,5 do
        local e=Fn(a,b,c,x)-Fn(aa,bb,cc,x)
        error=error+e*e
    end
    return error
end

function draw()
    --we'll print error results after each iteration, when GA yields and the screen refreshes
    --to see the coefficients, we'll have to transform the random numbers like we did in Fitness
    --these results should all trend towards zero
    local aa=(BestSol[1] - 0.5)/(BestSol[2] - 0.5)
    local bb=(BestSol[3] - 0.5)/(BestSol[4] - 0.5)
    local cc=(BestSol[5] - 0.5)/(BestSol[6] - 0.5)
    print(currIter,"Err",string.format("%.1f",Best),"CoeffErr",
            string.format("%.1f",aa-a),
            string.format("%.1f",bb-b),
            string.format("%.1f",cc-c))
    coroutine.resume(cos) --important, this has to go somewhere to keep the iterations going
end

--this function just calculates the quadratic
function Fn(a,b,c,x)
    return a*x*x + b*x + c
end
--# GA
--GA class
--You can treat this like a black box

GA = class()

local nVar
local nPop
local pctCrossover
local pctMutation
local pctSample
local replaceMethod
local fnFitness
local arrScore={}
local arrVar={}
local arrSort={}
local bestPop
local bestFitness

function GA:init(nVar1,nPop1,fnFitness1,fnClass)
    nVar=nVar1
    nPop=nPop1
    replaceMethod="Random"
    fnFitness=fnFitness1
    fnClass=fnClass1
end

function GA:DoRun(r,c,m,s)
    pctCrossover=c
    pctMutation=m
    pctSample=s
    nCross = math.floor(pctCrossover * nPop)
    --if first iteration, create array ofrandom number
    if r==1 then
        arrVar={}
        for i=1,nPop do
            arrVar[i]={}
            for j=1,nVar do
                arrVar[i][j]=math.random()
            end
            --get score
            if fnClass==nil then
                arrScore[i]=fnFitness(arrVar[i])
            else
                arrScore[i]=fnFitness(fnClass,arrVar[i])
            end
        end
        arrSort=GA:SortLinked(arrScore)
        bestPop = arrSort[1]
        bestFitness = arrScore[bestPop]
        return arrVar[bestPop],bestFitness
    else
        for j = 1, nCross do
            GA:ReplacePop(arrSort[nPop + 1 - j])
        end
        arrSort=GA:SortLinked(arrScore)
        bestPop = arrSort[1]
        bestFitness = arrScore[bestPop]
        return arrVar[bestPop],bestFitness
    end
end

function GA:ReplacePop(p)
    local p1
    local p2
    local x1
    local x2
    local x3
  x1 = math.floor(nPop * pctSample * math.random() + 1)
  x2 = math.floor(nPop * pctSample * math.random() + 1)
  x3 = math.floor(nPop * pctSample * math.random() + 1)
  if arrScore[x1] < arrScore[x2] then
    if arrScore[x1] < arrScore[x3] then p1 = x1 else p1 = x3 end
  else
    if arrScore[x2] < arrScore[x3] then p1 = x2 else p1 = x3 end
  end
  x1 = math.floor(nPop * pctSample * math.random() + 1)
  x2 = math.floor(nPop * pctSample * math.random() + 1)
  x3 = math.floor(nPop * pctSample * math.random() + 1)
  if arrScore[x1] < arrScore[x2] then
    if arrScore[x1] < arrScore[x3] then p2 = x1 else p2 = x3 end
  else
    if arrScore[x2] < arrScore[x3] then p2 = x2 else p2 = x3 end
  end

  local j
  local Tmp1
  local Tmp2
  local rMethod
  if replaceMethod == "Random" then
    j = math.floor(math.random() * 2 + 1)
    if j==1 then
      rMethod = "EitherOr"
    else
      rMethod="2 Slices"
    end
  else
    rMethod = replaceMethod
  end
  if rMethod == "2 Slices" then
    j = math.floor(math.random() * nVar) + 1
    for i = 1, nVar do
      if math.random() < pctMutation then
        arrVar[p][i] = math.random()
      end
    end
    k=math.floor(math.random()*nVar+1)
    if k<1 then k=1 elseif k>=nVar then k=nVar end
    for i=1,k do
          arrVar[p][i] = arrVar[p1][i]
    end
    for i=k+1,nVar do
          arrVar[p][i] = arrVar[p2][i]
    end
  elseif rMethod=="EitherOr" then
    for i = 1, nVar do
      if math.random() < pctMutation then
        arrVar[p][i] = math.random()
      else
        if math.random()< 0.5 then
          arrVar[p][i] = arrVar[p1][i]
        else
          arrVar[p][i] = arrVar[p2][i]
        end
      end
    end
  end

  for i = 1,nVar do
    if arrVar[p][i] < 0.000001 then
        arrVar[p][i] = 0.000001
    elseif arrVar[p][i]> 0.999999 then
        arrVar[p][i] = 0.999999
    end
  end

  if fnClass==nil then
        arrScore[p]=fnFitness(arrVar[p])
  else
        arrScore[p]=fnFitness(fnClass,arrVar[p])
  end
  arrScore[p]=fnFitness(arrVar[p])

end

function GA:SortLinked(A)
    local switches
    local gap
    local Low
    local Num
    local hold
    Low = 1
    Num = #A
    sList={}

    for ii = 1,Num do
      sList[ii] = ii
    end

    gap = Num - Low + 1
    repeat
      gap = math.floor((gap / 13) * 10)
      if gap < 1 then gap = 1 end
      if gap == 9 or gap == 10 then gap = 11 end
      switches = 0
      for ii = Low, Num - gap do
        if A[sList[ii]] > A[sList[ii + gap]] then
          sList[ii],sList[ii + gap]=sList[ii + gap],sList[ii]
          switches = 1
        end
      end
    until switches + gap == 1
    return sList
end