The Fresa plant propagation algorithm for optimization

@misc{fresa,
  author={Eric S. Fraga},
  title={ericsfraga/Fresa.jl: First public release (r2021.06.30)},
  journal={Zenodo},
  note={{doi:10.5281/zenodo.5045812}},
  year={2021},
  doi={10.5281/zenodo.5045812},
  url={https://doi.org/10.5281/zenodo.5045812}
}

@inproceedings{salhi-fraga-2011a,
  author={Abdellah Salhi and Eric S. Fraga},
  title={Nature-inspired optimisation approaches and the new plant propagation algorithm},
  booktitle={Proceedings of ICeMATH 2011, The International Conference on Numerical Analysis and Optimization},
  year={2011},
  pages={K2:1--8},
  address={Yogyakarta},
  month={6}
}

using Fresa

function objective(x)
    # calculate the objective function value
    z = 5*x[1]^2 + 4*x[2]^2 - 60*x[1] - 80*x[2]
    # evaluate the constraints so that feasible points result in a
    # non-positive value, i.e. 0 or less, but infeasible points give a
    # positive value.  We choose the maximum of both constraints as
    # the value to return as an indication of feasibility
    g = maximum( [ 6*x[1] + 5*x[2] - 60
                   10*x[1] + 12*x[2] - 150 ] )
    # return the objective function value along with indication of
    # feasibility
    (z, g)
end

domain = Fresa.Domain( x -> [ 0.0,  0.0 ],  # lower bounds
                       x -> [ 8.0, 12.5 ] ) # upper bounds

initialpopulation = [ Fresa.createpoint( [4.0, 6.25 ], objective) ]

best, population = Fresa.solve( objective, initialpopulation, domain )

# -*- mode: org; -*-
#+startup: show3levels
: Fresa PPA last change [2022-08-06 16:04]
objective (generic function with 1 method)
Fresa.Domain(var"#1#3"(), var"#2#4"())
1-element Vector{Fresa.Point}:
 f([4.0, 6.25])=[-503.75] g=-4.75
** solve objective [2022-08-06 16:04]
#+name: objectivesettings
| variable | value |
|-
| ngen | 100 |
| npop | 10 |
| nrmax | 5 |
| ns | 100 |
| elite | true |
| archive | false |
| fitness | hadamard |
| steepness | 1.0 |
|-
: function evaluations performed sequentially.
*** initial population
#+name: objectiveinitial
|-
| z1 | g | x |
|-
| -503.75 | -4.75 | [4.0, 6.25] |
|-

*** evolution
#+name: objectiveevolution
#+plot: ind:1 deps:(6) with:"points pt 7" set:"logscale x"
|       gen |      npop |        nf |    pruned |     t (s) | z1        |         g |
|-
|         1 |         1 |         1 |         0 |      1.62 | -503.75 | -4.75 |
|         2 |         2 |         2 |         0 |      2.24 | -503.75 | -4.75 |
|         3 |         3 |         4 |         0 |      2.24 | -508.98024613972393 | -3.963481299254326 |
|         4 |        10 |        14 |         1 |      2.24 | -516.708162060443 | -2.25158356765359 |
|         5 |        26 |        40 |         2 |      2.24 | -520.6840859905025 | -1.1924432568181942 |
|         6 |        29 |        68 |         2 |      2.24 | -523.7646499331361 | -0.21175198131829376 |
|         7 |        30 |       103 |         8 |      2.30 | -523.7646499331361 | -0.21175198131829376 |
|         8 |        24 |       129 |        11 |      2.30 | -524.1229955815547 | -0.35088947745457233 |
|         9 |        30 |       160 |        13 |      2.30 | -524.6890696989392 | -0.1186704914001453 |
|        10 |        28 |       188 |        14 |      2.30 | -525.265180854414 | -0.7525770581036753 |
|        20 |        18 |       458 |        67 |      2.30 | -529.2960864238084 | -0.07945091191650988 |
|        30 |        23 |       764 |       125 |      2.31 | -529.6657551053829 | -0.017628945114793737 |
|        40 |        22 |      1028 |       191 |      2.31 | -529.6657551053829 | -0.017628945114793737 |
|        50 |        28 |      1289 |       249 |      2.31 | -529.6657551053829 | -0.017628945114793737 |
|        60 |        15 |      1575 |       364 |      2.32 | -529.6657551053829 | -0.017628945114793737 |
|        70 |        33 |      1894 |       416 |      2.32 | -529.6657551053829 | -0.017628945114793737 |
|        80 |        21 |      2182 |       495 |      2.33 | -529.6657551053829 | -0.017628945114793737 |
|        90 |        22 |      2451 |       545 |      2.34 | -529.6657551053829 | -0.017628945114793737 |
|       100 |        19 |      2741 |       597 |      2.34 | -529.6657551053829 | -0.017628945114793737 |
*** Fresa run finished
: nf=2767 npruned=601
(f([3.7020260052569407, 7.554043004668713])=[-529.6657551053829] g=-0.017628945114793737, |-
| z1 | g | x |
|-
| -529.6657551053829 | -0.017628945114793737 | [3.7020260052569407, 7.554043004668713] |
| -366.54435704507523 | -27.21876533995116 | [1.7692502532766028, 4.433146628077845] |
| -414.2529968479937 | -20.951491873449086 | [1.3966036701257092, 6.133777221159331] |
| -528.8463299506853 | -0.22900672389539523 | [3.6795486014901346, 7.53874033343276] |
| -530.384414390094 | 0.16910579472867937 | [3.7209424919011393, 7.568690168664369] |
| -443.1260111837601 | -16.228285953791016 | [2.86399171858194, 5.317552746943469] |
| -445.6210305308239 | -15.601715891265954 | [3.016294412254267, 5.260103527041689] |
| -414.6744628574651 | -20.448352990168345 | [2.5321364109887003, 4.871765708779892] |
| -429.18276623025895 | -18.390573964068352 | [2.672394747864535, 5.1150115097488875] |
| -527.1298300049249 | -0.5625086104374191 | [3.825709902044894, 7.296646395458643] |
| -524.7657659486813 | -1.1402676069200197 | [3.7873107530029255, 7.227173575012485] |
| -524.5668961935723 | -1.1915319178131085 | [3.780900324334527, 7.224613227235946] |
| -525.8468778251139 | -0.8906651408825468 | [3.7913260301461618, 7.272275735648097] |
| -527.7529757793898 | -0.4269495828330747 | [3.8153534947185843, 7.336185889771085] |
| -528.3084551642806 | -0.347630580465065 | [3.734607564840653, 7.448944806098202] |
| -528.6635168151147 | -0.22202810180591115 | [3.801579393238273, 7.39369910775289] |
| -392.9526971218319 | -22.862738531064636 | [2.5817576982271797, 4.329343055914457] |
| -457.4148735934202 | -13.81074842840058 | [3.112459710647433, 5.5028986615429645] |
| -424.5002108062361 | -12.353125520526376 | [4.648811676286212, 3.9508008843512705] |
| -423.1306398623829 | -15.348331900055584 | [3.973866475752854, 4.161693849085458] |
| -427.5637392336674 | -11.136583646876801 | [4.8401741154869224, 3.964474332040334] |
| -359.63180208155165 | -15.950095856635826 | [5.141747031659222, 2.6398843906817673] |
| -414.0489650018614 | -9.615662219485117 | [5.402713812562485, 3.593610981027995] |
|-
)

# All code copyright © Eric S Fraga. 
# Date of last change in version variable below.
module Fresa
version = "[2022-08-07 13:14]"
using Dates
using Distributed
using Printf
function __init__()
    if myid() == 1
        println("# -*- mode: org; -*-")
        println("#+startup: show3levels")
        println(": Fresa PPA last change $version")
    end
end

function createpoint(x,f,parameters = nothing,ancestor = nothing)
    z = 0
    g = 0
    if ! ( parameters isa Nothing )
        (z, g) = f(x, parameters)
    else
        (z, g) = f(x)
    end
    if g isa Int
        g = float(g)
    end
    p = Nothing
    if rank(z) == 1
        p = Point(x, z, g, ancestor)
    elseif rank(z) == 0
        p = Point(x, [z], g, ancestor)
    else
        error("Fresa can only handle scalar and vector criteria, not $(typeof(z)).")
    end
    return p
end

function createpoint(x,f)
    return createpoint(x,f,nothing,nothing)
end
function createpoint(x,f,parameters)
    return createpoint(x,f,parameters,nothing)
end

function fitness(pop, fitnesstype, steepness, generation, ngen)
    l = length(pop)
    indexfeasible = (1:l)[map(p->p.g,pop) .<= 0]
    indexinfeasible = (1:l)[map(p->p.g,pop) .> 0]
    @debug "Feasible/infeasible breakdown" indexfeasible indexinfeasible maxlog=3
    fit = zeros(l)
    factor = 1              # for placement in fitness interval (0,1)
    if length(indexfeasible) > 0
        feasible = view(pop,indexfeasible)
        # use objective function value(s) for ranking
        feasiblefit = vectorfitness(map(p->p.z,feasible), fitnesstype, steepness, generation, ngen)
        if length(indexinfeasible) > 0
            feasiblefit = feasiblefit./2 .+ 0.5 # upper half of fitness interval
            factor = 2                        # have both feasible & infeasible
        end
        fit[indexfeasible] = (feasiblefit.+factor.-1)./factor
    end
    if length(indexinfeasible) > 0
        # squeeze infeasible fitness values into (0,0.5) or (0,1) depending
        # on factor, i.e. whether there are any feasible solutions as well or not
        infeasible = view(pop,indexinfeasible)
        # use constraint violation for ranking as objective function values
        # may not make any sense given that points are infeasible
        fit[indexinfeasible] = vectorfitness(map(p->p.g, infeasible),
                                             fitnesstype,
                                             steepness,
                                             generation,
                                             ngen
                                             ) / factor;
    end
    fit
end

"""
For single objective problems, the fitness is simply the normalised
objective function value.

For multi-objective cases, there are three alternative measures of
fitness ranking possible.  The first is based on the Hadamard product
of the rank of each member of population accoring to each
criterion.  The second is based on a weighted Borda ranking based on
each criterion ranking.  Finally, a measure based on dominance,
similar to that used by the popular NSGA-II genetic algorithm, is
available.

"""
function vectorfitness(v, fitnesstype, steepness, generation, ngen)
    # determine number of objectives (or pseudo-objectives) to consider in
    # ranking
    l = length(v)
    if l == 1
        # no point in doing much as there is only one solution
        fit = [0.5]
    else
        m = length(v[1])
        # println("VF: v=$v")
        # println("  : of size $(size(v))")
        if m == 1                   # single objective 
            fitness = [v[i][1] for i=1:l]
        else                  # multi-objective
            rank = ones(m,l); #rank of each solution for each objective function 
            if fitnesstype == :hadamard
                for i=1:m
                    rank[i,sortperm([v[j][i] for j=1:l])] = 1:l
                end
                # hadamard product of ranks
                fitness = map(x->prod(x), rank[:,i] for i=1:l)
            elseif fitnesstype == :borda
                for i=1:m
                    rank[i,sortperm([v[j][i] for j=1:l])] = 1:l
                end
                # borda sum of ranks
                fitness = map(x->sum(x), rank[:,i] for i=1:l)
            elseif fitnesstype == :nondominated
                # similar to that used by NSGA-II (Deb 2000)
                fitness = zeros(l)
                maxl = assigndominancefitness!(fitness,v,1)
                # println("Resulting fitness: $fitness")
            else
                throw(ArgumentError("Type of fitness evaluation must be either :borda, :nondominated, or :hadamard, not $(repr(fitnesstype))."))
            end
        end
        # normalise (1=best, 0=worst) while avoiding
        # extreme 0,1 values using the hyperbolic tangent
        fit = adjustfitness(fitness, steepness, generation, ngen)
        # println(":  scaled fitness: $fit")
        @debug "Fitness calculations" v[1][1] v[2][1] v[l][1] fitness[1] fitness[2] fitness[l] fit[1] fit[2] fit[l] maxlog=3
    end
    fit
end

function adjustfitness(fitness, steepness, generation, ngen)
    if (maximum(fitness)-minimum(fitness)) > eps()
        s = steepness
        if steepness isa Tuple
            a = (2*steepness[1]-2*steepness[2])/3
            b = - (3*steepness[1] - 3*steepness[2])/ngen^2
            d = steepness[1]
            s = a*generation^3 + b*generation^2 + c*generation + d
            @debug "Steepness " s "at generation" g
        end  
        fit = 0.5*(tanh.(4*s*(maximum(fitness) .- fitness)
                         / (maximum(fitness)-minimum(fitness))
                         .- 2*s) .+ 1)
    else
        # only one solution (or all solutions the same) in population
        fit = 0.5*ones(length(fitness))
    end
    fit
end

c(g) := a*g^3 + b*g^2 + c*g + d;
define(d(g), diff(c(g),g));
equations: [c(0) = s1,
d(0) = 0,
c(n) = s2,
d(n) = 0];
solution: solve(equations, [a, b, c, d]);
for i: 1 thru length(solution[1]) do print(solution[1][i])$

function assigndominancefitness!(f,v,l)
    # assign value l to all members of v which dominate rest and then
    # recurse on those which are dominated
    (p, d) = paretoindices(v)
    # println("Assigning fitness $l to $p")
    f[p] .= l
    if !isempty(d)
        assigndominancefitness!(view(f,d),v[d],l+1)
    else
        l
    end
end

function neighbour(x :: Array{Float64,1},
                   a :: Array{Float64,1},
                   b :: Array{Float64,1},
                   f :: Float64
                   ) :: Array{Float64,1}
    xnew = x .+ (1.0 .- f) .* 2(rand(length(x)).-0.5) .* (b.-a)
    xnew[xnew.<a] = a[xnew.<a];
    xnew[xnew.>b] = b[xnew.>b];
    return xnew
end

function neighbour(x :: Float64,
                   a :: Float64,
                   b :: Float64,
                   f :: Float64
                   ) :: Float64
    # allow movements both up and down
    # in the domain for this variable
    newx = x + (b-a)*(2*rand()-1)/2.0 * (1-f)
    if newx < a
        newx = a
    elseif newx > b
        newx = b
    end
    newx
end

function printHistoryTrace(p :: Point)
    a = p.ancestor
    while ! (a isa Nothing)
        println("| $(a.generation) | $(a.fitness) |")
        a = a.point.ancestor
    end
end

function prune(pop :: AbstractArray, tolerance)
    npruned = 0
    z = map(p->p.z, pop)
    # @show z[1]
    # println("typeof(z)=$(typeof(z))")
    l = length(z)
    # println("typeof(z[1])=$(typeof(z[1]))")
    n = length(z[1])
    # @show n
    zmin = zeros(n)
    zmax = zeros(n)
    try 
        for i=1:n
            row = [z[j][i] for j=1:l]
            zmin[i] = minimum(row)
            zmax[i] = maximum(row)
            if zmax[i] - zmin[i] < 100*eps()
                zmax[i] = zmin[i]+100*eps()
            end
        end
        pruned = [pop[1]]
        for i=2:l
            similar = false
            for j=1:length(pruned)
                if all(abs.(z[i]-pruned[j].z) .< tolerance*(zmax-zmin))
                    similar = true;
                    break;
                end
            end
            if !similar
                push!(pruned,pop[i])
            else
                npruned += 1
            end
        end
        (pruned, npruned)
    catch e
        # println("prune method error: $e")
        if e isa MethodError
            # probably (possibly) due to objective function type not
            # being a number.  In this case, we try again but looking
            # at the decision variable values instead.
            x = map(p->p.x, pop)
            # println("typeof(z)=$(typeof(z))")
            l = length(x)
            # start building up the population that remains after
            # pruning.  The first entry will always be there as any
            # similar solutions will not be included by the search
            # that follows.
            pruned = [pop[1]]
            try
                for i=2:l
                    similar = false
                    # now check this solution against all those already in
                    # the list we are collating
                    for j=1:length(pruned)
                        if all(Float64.(abs.(x[i]-pruned[j].x)) .< tolerance)
                            similar = true;
                            break;
                        end
                    end
                    if !similar
                        push!(pruned,pop[i])
                    else
                        npruned += 1
                    end
                end
                (pruned, npruned)        
            catch e
                # println("prune method second error: $e")
                if e isa MethodError
                    # this is now probably/possibly due to not being
                    # to find the difference between two decision
                    # points.  In that case, return the whole
                    # original population
                    (pop, 0)
                else
                    @error "Unexpected error in prune method for points" e
                end             # if method error
            end                 # try pruning on points
        else
            @error "Unexpected error in prune method for objective function values" e
        end                     # if is a method error
    end                         # try pruning on objective function values
end                             # function

function randompopulation(n, f, parameters, p0, domain :: Domain)
    p = Point[]                 # population object
    for j in 1:n
        # l = domain.lower(p0.x)
        # @show l
        # u = domain.upper(p0.x)
        # @show u
        # x = randompoint(l,u)
        # push!(p, createpoint(x, f, parameters))
        push!(p, createpoint(randompoint(domain.lower(p0.x), domain.upper(p0.x)), f, parameters))
    end
    p
end

function randompoint(a :: Float64, b :: Float64)
    x = a + rand()*(b-a)
end
function randompoint(a, b)
    x = a + rand(length(a)).*(b-a)
end

function select(f)
    l = length(f)
    ind1 = rand(1:l)
    if ind1 == 0
        ind1 = 1
    end
    ind2 = rand(1:l)
    # println("Comparing $ind1 to $ind2")
    if f[ind1] > f[ind2]
        return ind1
    else
        return ind2
    end
end

""" 

Solve an optimisation problem, defined as the minimization of the
values returned by the objective function, `f`.  `f` returns not only
the objective function values, an array of `Float64` values, but also
a measure of feasibility (≤0) or infeasibility (>0).  The problem is
solved using the Fresa algorithm.  `p0` is the initial population
which has to have at least one member, a `Point`, and `a` and `b` are
*bounds* on the search space.

The return values for the solution of a single criterion problem are
the best point and the full population at the end of the search. 

For a multi-objective problem, the returned values are the set of
indices for the points within the full population (the second returned
value) approximating the *Pareto* front.

The population will consist of an array of `Fresa.Point` objects, each
of which will have the point in the search space, the objective
function value and the feasibility measure.

"""
function solve(f, p0, domain;        # required arguments
               parameters = nothing, # allow parameters for objective function 
               archiveelite = false, # save thinned out elite members
               elite = true,    # elitism by default
               fitnesstype = :hadamard, # how to rank solutions in multi-objective case
               multithreading = false, # use multiple threads for objective function evaluation
               ngen = 100,           # number of generations
               npop = 10,            # population size: fixed (single value) or dynamic (tuple)
               nrmax = 5,            # number of runners maximum
               ns = 100,             # number of stable solutions for stopping
               output = 1,           # how often to output information
               plotvectors = false,  # generate output file for search plot
               populationoutput = false, # output population every generation?
               steepness = 1.0,      # show steep is the adjustment shape for fitness
               tolerance = 0.001,    # tolerance for similarity detection
               usemultiproc = false) # parallel processing by Fresa itself?
    output > 0 && println("** solve $f $(orgtimestamp(now()))")
    tstart = time()
    nf = 1                   # number of function evaluations
    npruned = 0              # number solutions pruned from population
    nz = length(p0[1].z)     # number of criteria
    pop = copy(p0);          # create/initialise the population object
    if archiveelite
        archive = Point[]
    end
    if output > 0
        println("#+name: $(f)settings")
        println("| variable | value |")
        println("|-")
        println("| ngen | $ngen |")
        println("| npop | $npop |")
        println("| nrmax | $nrmax |")
        println("| ns | $ns |")
        println("| elite | $elite |")
        println("| archive | $archiveelite |")
        println("| fitness | $fitnesstype |")
        println("| steepness | $steepness |")
        println("|-")
        # output != 0 && println(": solving with ngen=$ngen npop=$npop nrmax=$nrmax ns=$ns")
        # output != 0 && println(": elite=$elite archive elite=$archiveelite fitness type=$fitnesstype")
    end
    if plotvectors
        plotvectorio = open("fresa-vectors-$(orgtimestamp(now())).data", create=true, write=true)
        output > 0 && println(": output of vectors for subsequent plotting")
    end
    # if npop was given as a tuple, we are to have a dynamic
    # population size.  This only makes sense for multi-objective
    # optimization problems so a warning will be given otherwise.
    npopmin = npop
    npopmax = npop
    if npop isa Tuple
        if nz > 1
            npopmin = npop[1]
            npopmax = npop[2]
            if npopmin > npopmax
                error("Dynamic population sizing requires min <= max; you specified $npop")
            end
            npop = npopmin      # start with minimum possible
        else
            println("*Warning*: you have specified a tuple for population size: $npop")
            println("This only makes sense for multi-objective optimization problems.")
            println("npop will be set to $(npop[1]).")
            npop = npop[1]      # be optimistic and use minimum given
        end
    end
    # we use multithreading if asked for *and* if we have more than
    # one thread available
    multithreading = multithreading && Threads.nthreads() > 1 
    # we use parallel computing if we have more than one processor
    parallel = usemultiproc && nprocs() > 1
    # parallel = false
    if output > 0
        println(": function evaluations performed ",
                parallel
                ? "in parallel with $(nprocs()) processors."
                : (multithreading
                   ? "in parallel with $(Threads.nthreads()) threads."
                   : "sequentially."))
        println("*** initial population")
        println("#+name: $(f)initial")
        println(pop)
    end
    if output > 0
        println("*** evolution")
        println("#+name: $(f)evolution")
        println("#+plot: ind:1 deps:(6) with:\"points pt 7\" set:\"logscale x\"")
        @printf("| %9s | %9s | %9s | %9s | %9s |", "gen", "npop",
                (elite && nz > 1) ? "pareto" : "nf", "pruned", "t (s)")
        for i in 1:nz
            @printf(" z%-8d |", i)
        end
        @printf(" %9s |", "g")
        @printf("\n|-\n")
    end
    # now evolve the population for a predetermined number of generations
    for gen in 1:ngen
        # evaluate fitness which is adjusted depending on value of
        # steepness, a value that may depend on the generation
        fit = fitness(pop, fitnesstype, steepness, gen, ngen)
        if gen == 1
            @debug "Initial fitness" f=fit
        end
        # sort
        index = sortperm(fit)
        if populationoutput
            println("\nGeneration $gen full population is:")
            println(pop)
            println("Fitness vector: $fit")
        end
        # and remember best which really only makes sense in single
        # criterion problems but is best in multi-objective case in
        # the ranking measure used by Fresa
        best = pop[index[end]]
        # if elitism is used
        if elite
            if nz > 1
                # elite set is whole pareto set unless it is too
                # big. Recall that the pareto function returns the set
                # of indices into the population
                wholepareto = pareto(pop)[1]
                # if using dynamic population sizing, adjust the population
                npop = 2 * length(wholepareto)
                if npop < npopmin
                    npop = npopmin
                end
                if npop > npopmax
                    npop = npopmax
                end
                # now check that the pareto is not too big.  if it is, thin it out
                if length(wholepareto) > ceil(npop/2)
                    newpop, removed = thinout(pop, fit, wholepareto, ceil(Int,npop/2))
                    if archiveelite
                        archive = prune(append!(archive, removed), tolerance)[1]
                        archive = archive[pareto(archive)[1]]
                    end
                else
                    newpop = pop[wholepareto]
                end
            else
                # elite set is single element only
                newpop = [best]
            end
            # if plotting vectors for the search, include elitism
            if plotvectors
                for p in newpop
                    write(plotvectorio, "$(gen-1) $(p.x)\n$gen $(p.x)\n\n")
                end
            end
        else
            newpop = Point[]
        end
        if output >= 0
            print(stderr, ": $gen np=$(length(newpop))/$npop",
                  archiveelite ? " na=$(length(archive))" : "",
              " with most fit z=$(best.z)           \r")
            # if output has been requested, check to see if output is
            # required now and then also check to see if the frequency
            # needs to be reduced.
            if output > 0
                if gen%output == 0
                    @printf("| %9d | %9d | %9d | %9d | %9.2f |", gen, length(fit),
                            (elite && nz > 1) ? length(newpop) : nf, npruned, time()-tstart)
                    for i = 1:length(best.z)
                        print(" $(best.z[i]) |")
                    end
                    print(" $(best.g) |")
                    println()
                end
                if 10^(floor(log10(gen))) > output
                    output = 10^(Int(floor(log10(gen))))
                end
            end
        end
        # if we are using any form of multiprocessing, either threads
        # or multiple cores, create an array to store all new points
        # which we evaluate later in parallel.  Ideally, also keep
        # track of the points from which new points are derived to
        # provide the backward link through the evolution but this is
        # currently disabled as the creation of the Ancestor object
        # requires more information than I am currently storing away.
        if multithreading || parallel
            x = Any[] # typeof(newpop[1].x)[]
            # points = Point[]
        end
        # now loop through population, applying selection and then
        # generating neighbours
        l = length(pop)
        for i in 1:min(l,npop)
            s = select(fit)
            # println(": selection $i is $s")
            # println(": size of pop is $(size(pop))")
            selected = pop[s]
            if !elite
                # if no elitism, we ensure selected members remain in population
                push!(newpop, selected)
                if plotvectors
                    write(plotvectorio, "$(gen-1) $(selected.x)\n$gen $(selected.x)\n\n")
                end
            end
            # number of runners to generate, function of fitness
            nr = ceil(fit[s]*nrmax*rand())
            if nr < 1
                nr = 1
            end
            # println(": generating $nr runners")
            for r in 1:nr
                # create a neighbour, also function of fitness,
                # passing through the lower and upper bounds
                # appropriate for the particular solution point
                newx = neighbour(pop[s].x, domain.lower(pop[s].x), domain.upper(pop[s].x), fit[s])
                # for parallel evaluation, we store the neighbours and
                # evaluate them later; otherwise, we evaluate
                # immediately and save the resulting point
                if multithreading || parallel
                    push!(x, newx)
                    # push!(points, pop[s])
                else
                    push!(newpop, createpoint(newx, f, parameters, Ancestor(pop[s],fit[s],gen)))
                    if plotvectors
                        write(plotvectorio, "$(gen-1) $(pop[s].x)\n$gen $newx\n\n")
                    end
                    nf += 1
                end
            end
            # remove selected member from the original population so
            # it is not selected again
            splice!(fit, s)
            splice!(pop, s)
        end
        # if we are making use of parallel computing, we evaluate all
        # points generated in previous loop.  Parallel processing is
        # done either via multithreading or with multiple
        # processors.  The former is easier as it's based on shared
        # memory.
        if multithreading       # using threads and shared memory
            results = Array{Point}(undef,length(x))
            Threads.@threads for i ∈ 1:length(x)
                results[i] = createpoint(x[i],f,parameters)
            end
            append!(newpop, results)
        elseif parallel        # using multiple processors with remote calls
            # will be used to collect results from worker processors
            results = Array{Future,1}(undef, nprocs())
            i = 0;
            while i < length(x)
                # issue remote evaluation call
                for j=1:nprocs()
                    if i+j <= length(x) 
                        # TODO: the information about the ancestor is
                        # not available; this needs to be stored above
                        results[j] = @spawn createpoint(x[i+j],f,parameters)
                        nf += 1
                    end
                end
                # now wait for results
                for j=1:nprocs()
                    if i+j <= length(x)
                        push!(newpop, fetch(results[j]))
                    end
                end
                i += nprocs()
            end
        end
        # and finally, if we have elitism, remove any duplicate points
        # in the new population and make it the current population for
        # the next generation; otherwise, simply copy over
        if elite
            (pop, nn) = prune(newpop, tolerance)
            npruned += nn
        else
            pop = newpop
        end
    end
    output > 0 && println("*** Fresa run finished\n: nf=$nf npruned=$npruned", archiveelite ? " archived=$(length(archive))" : "")
    if plotvectors
        close(plotvectorio)
    end
    if nz == 1
        fit = fitness(pop, fitnesstype, steepness, ngen, ngen)
        index = sortperm(fit)
        best = pop[index[end]]
        return best, pop
    else
        return pareto(archiveelite ? append!(pop,archive) : pop)[1], pop
    end
end

function thinout(pop, fit, pareto, n::Int)
    indices = sortperm(fit[pareto])
    return pop[pareto[indices[end-n+1:end]]], pop[pareto[indices[1:end-n]]]
end

end

# load in the Fresa optimization package
using Fresa
# specify the dimension of the search space
nx = 2
# create an initial point in the search space
x0 = 0.5*ones(nx)
# specify the domain for the search, x ∈ [0,10]ⁿ
domain = Fresa.Domain(x -> zeros(length(x)), x -> 10*ones(length(x)))
# the actual objective function
f = x -> ((x[1]-3)^2+(x[2]-5)^2+8, 0)
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x0,f)]
# now invoke Fresa to solve the problem
@time best, pop = Fresa.solve(f, p0, domain)
# output the results
println("Population at end:\n$pop")
println("Best solution is f($( best.x ))=$( best.z ) with g=$( best.g )")

println("\nHistory trace, by generation number, of fitness value of solution selected for propagation which results in a new best solution:")
println("#+plot: ind:1 deps:(2) with:\"linespoints pt 7 ps 0.25\" set:nokey set:\"yrange [0:1]\" set:\"xrange [0:*]\" set:\"xlabel 'Generation'\" set:\"ylabel 'fitness'\"")
Fresa.printHistoryTrace(best)

using Fresa
nx = 2
x0 = 0.5*ones(nx)
# specify the domain for the search, x ∈ [0,10]ⁿ
domain = Fresa.Domain(x -> zeros(length(x)), x -> 10*ones(length(x)))
rosenbrock(x) = ([(1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2], 0)
# f = x -> ((x[1]-3)^2+(x[2]-5)^2+8, 0)
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x0,rosenbrock)]
# now invoke Fresa to solve the problem
best, pop = Fresa.solve(rosenbrock, p0, domain; ngen=1000, tolerance=1e-8)
println("Population at end: $pop")
println("Best solution is f($( best.x ))=$( best.z ) with g=$( best.g )")

using Fresa
nx = 20
x0 = 0.5*ones(nx)
# specify the domain for the search, x ∈ [0,10]ⁿ
domain = Fresa.Domain(x -> zeros(length(x)), x -> 10*ones(length(x)))
rosenbrock(x) = (sum([100 * (x[i+1]-x[i]^2)^2 + (1-x[i])^2 for i ∈ 1:length(x)-1]), 0)
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x0,rosenbrock)]
# now invoke Fresa to solve the problem
best, pop = Fresa.solve(rosenbrock, p0, domain; npop=100, ngen=1000, tolerance=1e-8, multithreading=true)
println("Best solution is f($( best.x ))=$( best.z ) with g=$( best.g )")

using Fresa
nx = 2
# specify the domain for the search, x ∈ [0,10]ⁿ
domain = Fresa.Domain(x -> zeros(length(x)), x -> ones(length(x)))
# initial point in domain
x = rand(nx)
# objective function 
f = x -> ( [sin(x[1]-x[2]); cos(x[1]+x[2])], 0)
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x,f)]
# now invoke Fresa to solve the problem
pareto, population = Fresa.solve(f, p0, domain;
                                 #fitnesstype = :hadamard,
                                 #fitnesstype = :borda,
                                 fitnesstype = :nondominated,
                                 ngen=200,
                                 npop=(20,40),
                                 plotvectors=true,
                                 tolerance=0.01)

println("**** Pareto front:")
println("#+plot: ind:1 deps:(2) with:points")
println(population[pareto])
#using BenchmarkTools
#@benchmark

using Fresa
using Profile
nx = 5
# specify the domain for the search, x ∈ [0,1]ⁿ
domain = Fresa.Domain(x -> zeros(length(x)), x -> ones(length(x)))
x = zeros(nx)
f = x -> ([ sum((x.-0.5).^2 .+ 1)
            sum(cos.(x))
            sum(sin.(x))],
          0)
# create the initial population consisting of this single point
p0 = [Fresa.createpoint(x,f)]
# now invoke Fresa to solve the problem
@profile for i=1
    pareto, population = Fresa.solve(f, p0, domain;
                                     archiveelite = false,
                                     npop=20, ngen=300,
                                     #output=100,
                                     tolerance=0.01)

    println("*** Pareto front:")
    println(population[pareto])
end
println("*** profile data")
println(": this may take some time so please wait")
Profile.print(format=:flat, sortedby=:count)

using Distributed
using Printf
@everywhere using Fresa
# define new type for mixed integer problems
# in general, this would be vectors of real and integer values
@everywhere struct MI
    x :: Float64
    y :: Int32
end
import Base
Base.show(io::IO, m::MI) = print(io, m.x, " ", m.y)
f = s -> (3s.y - 5s.x,
          max(2s.y + 3s.x - 24,
              3s.x - 2s.y - 8,
              2s.y^2 - 2*√s.y + 11s.y + 8s.x - 39 - 2*√s.x*s.y^2))
# bounds
domain = Fresa.Domain(x -> MI(1.0, 1),
                      x -> MI(6.0, 6))
# function to find a neighbouring solution for MI type decision points
function Fresa.neighbour(s :: MI,
                         a :: MI,
                         b :: MI,
                         f :: Float64) :: MI
    x = s.x + (b.x-a.x)*(1-f)*2*(rand()-0.5)
    x = x < a.x ? a.x : (x > b.x ? b.x : x)
    # for the integer variable, we move in one direction or the other
    # a random number of places depending on fitness
    positive = rand(Bool)
    r = rand()
    # @printf(": neighbour: f=%g r=%g\n", f, r)
    inc = ceil(f*r*(b.y-a.y)/2)
    # @printf(": neighbour: positive=%s inc=%d\n", positive, inc)
    y = s.y + (positive ? inc : -inc)
    y = y < a.y ? a.y : (y > b.y ? b.y : y)
    return MI(x,y)
end
# create the initial population consisting of a single MI point
p0 = [Fresa.createpoint(MI(1.0, 1),f)]
# now invoke Fresa to solve the problem
best, pop = Fresa.solve(f, p0, domain; ngen=100)
println("Population: $pop")
println("Best: f($(best.x)) = $(best.z), $(best.g)")

println("#+plot: ind:3 deps:(2) with:\"linespoints pt 7\" set:nokey set:\"yrange [0:1]\"")
ancestor = best.ancestor;
while ancestor != Some(nothing) && ! (ancestor isa Nothing)
    global ancestor
    println("| $(ancestor.point.z) | $(ancestor.fitness) | $(ancestor.generation) |")
    ancestor = ancestor.point.ancestor
end