Strawberry plant propagation algorithm

function [x z nf ninf bestgen] = strawberry(x0, a, b, f, ngen, npop, nrmax, ns, population_strategy, output)

global esfdebug
global mosd_numberimproved
global npruned
global nsimilar
global strawberry_archive
global strawberry_fitness_method
global strawberry_diversity_tolerance
global strawberry_neighbourfn 
global strawberry_octave
global strawberry_output_population
mosd_numberimproved = 0;
npruned = 0;
nsimilar = 0;
stderr = 2;                     # for MATLAB compatibility
warning ('error', 'Octave:broadcast');
## process arguments
if nargin<9, population_strategy = 4; end # population selection strategy
if nargin<8, ns = 100; end    # number of stable generations
if nargin<7, nrmax = 5; end   # number of runners
if nargin<6, npop = 100; end  # population size
if nargin<5, ngen = 1000; end # number of generations
if nargin<4
  disp('Error: need at least 4 arguments to Strawberry method: x0, a, b & f.')
  return
end
if nargin<10, 
  output = ngen/50;
end             # how often to output status of population
## process global variables
if isempty (strawberry_archive)
  strawberry_archive = 0;
endif
if strawberry_archive > 0
  printf(': pareto solutions are archived\n');
endif
if isempty (strawberry_fitness_method)
  strawberry_fitness_method = 1;  #choose MOGA fitness for now
end
strawberry_diversity_tolerance = norm(b-a) / 1000;

## determine the number of processors to use in parallel.  If the
## variable has not been set, we use the system function to determine
## how many processors there are
global strawberry_numberofprocessors
if isempty(strawberry_numberofprocessors)
  strawberry_numberofprocessors = 1 # not currently working for np>1, nproc();
endif
if strawberry_numberofprocessors > 1
  printf (': using %d processors for parallel processing.\n', strawberry_numberofprocessors);
else
  printf (': not using any parallel processing on a single processor system.\n');
endif
if strawberry_octave
  if strawberry_numberofprocessors > 1
    starttime = time;           # use wall clock time instead of CPU time
  else
    starttime = cputime;
  endif
else
  tic;                        # start wall clock timer
  starttime = 0;              # for later use in output statistics
endif
if isempty(strawberry_output_population)
  strawberry_output_population = 0
endif
strawberry_initialise
printf(': ngen=%d npop=%d nrmax=%d ns=%d population strategy=%d and tolerance=%f fitness method=%d\n', ngen, npop, nrmax, ns, population_strategy, strawberry_diversity_tolerance, strawberry_fitness_method)

constrained = 1;
try
  [z g] = f(x0);
catch err
  printf('Error identifier is %s\n', err.identifier);
  printf('Error message is %s\n', err.message);
  printf('\nThis could be due to the change in Strawberry that\n')
  printf('now requires the objective function to return a\n')
  printf('feasibility indication as well as the objective function\n')
  printf('value(s).\n\n')
  if (strcmp (err.identifier, 'Octave:undefined-function') || strcmp(err.identifier,'MATLAB:UndefinedFunction'))
    printf(': objective function is unconstrained.\n')
    constrained = 0;
  else
    ## not an error we expected to pass this on
    rethrow(err);
  endif
end_try_catch
if constrained
     printf(': objective function is constrained.\n');
endif

[phi, nx, nz, nf, ninf] = strawberry_newpopulation (x0,a,b,npop,f,constrained);
pop = strawberry_sort (phi,nx,nz);
gen = 0;
best = [];                      #will be defined after fitness evaluation
bestgen = gen;
nfunctionevaluations = npop; # initial population size

if output > 0
  printf('In the following table, $n_f$ is number of feasible points,\n')
  printf('$n_i$ for infeasible, $n_p$ number in Pareto set\n');
  printf('$c$ best constraint value and $g_b$ best generation.\n');
  printf('| %9s |', 'gen')
  if nz > 1
    for i=1:nz, printf('  $z_{%02d}$ |', i); endfor
    if constrained; printf(' %9s |', 'c'); endif
  else
    index = nx;
    if nx > 4
      index = 4;
    endif
    for i=1:index; printf(' %6sx%02d |', '', i); endfor
    for i=1:nz; printf(' %6sy%02d |', '', i); endfor
    if constrained; printf(' %9s |', 'c'); endif
  endif
  printf(' %5s | %5s |', '$n_f$', '$n_i$')
  if nz > 1
    printf(' %5s |', '$n_p$')
  else
    printf(' %5s |', '$g_b$')
  endif
  printf('\n');
  printf('|-\n');
endif

while gen < ngen && (nz > 1 || gen < bestgen+ns)
  gen = gen + 1;

endwhile
if gen < ngen
  printf('\n: strawberry terminating condition satisfied due to lack of improvement at gen=%d\n', gen);
endif
if output
  printf('\n: strawberry at end, nf=%d ninf=%d best solution:\n', nf, ninf)
  disp(best');
endif
## ensure the best solution is returned
pop = [best, pop];
x = pop(1:nx,:);
z = pop(nx+1:end,:);
if strawberry_octave
  if strawberry_numberofprocessors > 1
    endtime = time;            # wall clock time
  else
    endtime = cputime;
  endif
else
  endtime = toc;
endif
printf('\n: strawberry method elapsed time: %.1f with %d functions evaluated and %d solutions pruned with %d similar.\n', ...
       endtime-starttime, nf, npruned, nsimilar)
printf('MOSD improvements: %d\n', mosd_numberimproved);
endfunction

function d = strawberry_distance(f,n)
  d = (1-f) * 2*(rand(n,1)-0.5);
endfunction

function [fit, best] = strawberry_fitness(p,nx,nz,npop)
  global strawberry_fitness_method
  [n m] = size(p);
  debugprint (1, 'fitness, nx nz m n npop', [nx, nz, m, n, npop])
  if m <= 0
    disp('No solutions in the population');
  endif
  fit = zeros(m,1);
  factor = 1;                # for placement in fitness interval (0,1)
  indexfeasible = p(n,:) <= 0;
  nfeas = sum(indexfeasible);
  indexinfeasible = p(n,:) > 0;
  ninfeas = sum(indexinfeasible);
  if nfeas > 0
    feasiblefit = strawberry_vectorfitness(nz, p(nx+1:nx+nz,indexfeasible));
    ## the ranking of the feasible solutions is possibly adjusted to
    ## fit in only the upper half of the [0,1] interval if there are
    ## infeasible solutions as well.  The latter will have fitness
    ## values in the bottom half of the interval.
    if ninfeas > 0
      factor = 2;
    end
    fit(indexfeasible) = (feasiblefit+factor-1)/factor;
  end
  if ninfeas > 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
    fit(indexinfeasible) = strawberry_vectorfitness(1, p(n,indexinfeasible))/factor;
  end
  if nz < 2
    [bestfit index] = max(fit);
    best = p(:,index);
  else
    ## printf('There are %d feasible solutions when calculating fitness\n',nfeas)
    if nfeas > 0
      allindices = 1:m;
      feasibleindices = allindices(indexfeasible);
      ## printf('feasible population:\n')
      ## disp(p(:,indexfeasible))
      ## printf('and pareto from that population\n')
      ## disp(p(:,feasibleindices(findpareto(p(nx+1:nx+nz,indexfeasible)))))
      best = p(:,feasibleindices(findpareto(p(nx+1:nx+nz,indexfeasible))));
    else
      best = [];
    endif
  endif
endfunction

function fit = strawberry_vectorfitness(m,v)
  ## println("VF: v=$v")
  ## println("  : of size $(size(v))")
  if m == 1                   # single objective 
    l = length(v);
    [sorted indices] = sort(v);
    zmin = v(indices(1));
    zmax = v(indices(l));
    if l == 1 || abs(zmax-zmin) < eps()
      fit = 0.5*ones(1,l);
    else
      ## avoid extreme 0,1 values
      fit = tanh((zmax - v) / (zmax - zmin) - 0.5)+0.5;
    end
  else                  # multi-objective
    [m, l] = size(v);
    rank = ones(m,l); #rank of each solution for each objective function 
    for i=1:m
      [sorted indices] = sort(v(i,:));
      rank(i,indices) = 1:l;
    end
    ## hadamard product of ranks
    fitness = prod(rank);
    ## normalise and reverse meaning (1=best, 0=worst)
    if l == 1
      fit = 0.5;
    else
      ## avoid extreme 0,1 values
      fit = tanh(0.5 - fitness / max(fitness)) + 0.5;
    end
  end
  ## println("VF: fit=$fit")
end

function strawberry_initialise
  global strawberry_version
  global strawberry_numberofprocessors
  global strawberry_octave
  global moga_version
  global mosd_version
  if isempty (strawberry_version)
    strawberry_version = '2019.06.24 12:21:04';
    mosd_initialise
    moga_initialise
    printf (': STRAWBERRY %s, main function\n', strawberry_version)
    printf (': using MOSD %s, multi-objective steepest descent method.\n', mosd_version)
    printf (': and MOGA %s, multi-objective genetic algorithm fitness function.\n', moga_version)
    strawberry_octave = exist ('OCTAVE_VERSION', 'builtin') > 0;
    if strawberry_octave
      printf(': running in Octave\n')
      if strawberry_numberofprocessors > 1
        pkg load parallel             #need parcellfun in particular
      endif
    endif
  endif

global strawberry_neighbourfn
global strawberry_randomsolutionfn
if isempty(strawberry_randomsolutionfn)
  strawberry_randomsolutionfn = @strawberry_randomsolution;
  printf (': using default random solution generator function.\n');
endif
if isempty(strawberry_neighbourfn)
  strawberry_neighbourfn = @strawberry_neighbour;
  printf (': using default neighbouring solution generator function.\n')
endif

function diverse = strawberry_isdiverse (pop, x)
  global nsimilar
  global strawberry_diversity_tolerance
  i = 0;
  diverse = 1;
  while i < length(pop) && diverse
    diverse = norm(pop{i}-x) > strawberry_diversity_tolerance;
    i = i + 1;
  endwhile
  if ~ diverse, nsimilar = nsimilar + 1; endif
endfunction

function n = strawberry_neighbour(x,a,b,f)
  nx = length(x);
  dx = (b-a)/2 .* strawberry_distance (f,nx); # how far to run
  ## fprintf(stderr, 'Selected %d with fitness %g to move %g\n', s, N(s), norm(dx))
  debugprint(3, ': dx', dx(1:nx))
  n = x+dx;
  n(n<a) = a(n<a);        # reset to boundary
  n(n>b) = b(n>b);        # reset to boundary
endfunction

function [pop, nx, nz, nf, ninf] = strawberry_newpopulation (x0,a,b,npop,f,constrained)
  global strawberry_randomsolutionfn
  global strawberry_numberofprocessors
  nx = length(x0);
  printf(': %d decision variables\n', nx);
  nf = 0;
  ninf = 0;
  pop = [];                     # will be array of x+y values
  n = 0;
  ntries = 0;
  xcell{1} = x0;                        # start with initial guess

  for i=2:npop
    xcell{i} = strawberry_randomsolutionfn(a,b);
  endfor
  if strawberry_numberofprocessors > 1
    if constrained
      pop = parcellfun (strawberry_numberofprocessors, ...
                        @(x) [x; f(x); g(x)], xcell, 'VerboseLevel', 0 ...
                       );
    else
      pop = parcellfun (strawberry_numberofprocessors, ...
                        @(x) [x; f(x); 0], xcell, 'VerboseLevel', 0);
    endif
  else                          #sequential operation
    x = xcell{1};
    if constrained
      [z g] = f(x);
    else
      z = f(x);
      g = 0;
    endif
    y = [x; z; g];
    pop = zeros(length(y),npop);
    pop(:,1) = y;
    for i=2:npop
      x = xcell{i};
      nf = nf+1;
      if constrained
        [z g] = f(x);
      else
        z = f(x);
        g = 0;
      endif
      pop(:,i) = [x; z; g];
    endfor
    ## pop = cellfun (@(x) [x; f(x); 0], xcell);
  endif
  mtotal = size(pop,2);
  [n m] = size(pop);
  ninf = ninf + mtotal - m;
  ## each row consists of nx x values, nz y values and a feasibility
  ## indication
  nz = n-nx-1;
  printf(': %d objective function values\n', nz);
endfunction

function populationstatistics(s,p,nx,nz)
  [m n] = size(p);
  feasible = p(end,:) <= 0;
  nfeasible = sum(feasible);
  infeasible = p(end,:) >= 0;
  ninfeasible = sum(infeasible);
  printf('Population statistics for %s:\n', s)
  printf(': there are %d members of which %d are feasible.\n', n, nfeasible)
  printf(': min objective function values: ')
  printf('%g ', min(p(nx+1:nx+nz,:),[],2))
  printf('\n');
  if nfeasible > 0
    printf(': constraint values for feasible solutions in [%g,%g]\n', ...
           min(p(end,feasible)), max(p(end,feasible)));
  endif
  if ninfeasible > 0
    printf(': constraint values for infeasible solutions in [%g,%g]\n', ...
           min(p(end,infeasible)), max(p(end,infeasible)));
  endif
endfunction

function pruned = strawberry_prune (pop, nx, nz)
  global npruned
  global strawberry_diversity_tolerance
  npop = size(pop,2);
  if npop > 0
    pruned = pop(:,1);
    np = 1;                       # size of pruned set
    for i=2:npop
      diverse = 1;
      j = 0;
      while diverse && j < np
        j = j + 1;
        ## diversity based on x values alone
        diverse = norm (pop(1:nx,i)-pruned(1:nx,j)) > strawberry_diversity_tolerance;
        ## printf('div(%d,%d)=%g %d\r', i, j, norm (pop(1:nx,i)-pruned(1:nx,j)), diverse);
      endwhile
      if diverse
        pruned = [pruned, pop(:,i)];
        np = np + 1;            #keep track of how many copied over
      else
        ## prune one or the other.  ideally, we keep the better
        ## one.  If we cannot distinguish (in a multi-objective
        ## sense), we simply keep the one that is already there. 
        npruned = npruned + 1;
        if ~ dominates(pruned(nx+1:nx+nz,j),pop(nx+1:nx+nz,i))
          if dominates(pop(nx+1:nx+nz,i),pruned(nx+1:nx+nz,j))
            ## current solution is better, i.e. dominates, than previous
            ## one so replace
            pruned(:,j) = pop(:,i);
          endif
        endif
      endif
    endfor
  endif
endfunction

function x = strawberry_randomsolution(a, b)
  r = rand(length(a),1);
  x = a + r.*(b-a);
endfunction

function nr = strawberry_runners(N, nrmax)
  r = rand();
  nr = ceil (nrmax*N*r); # number of runners
  if nr < 1, nr = 1; endif
  debugprint (4, ':runners N r runners max', [N, r, nr, nrmax])
endfunction

function s = strawberry_select (N)
  ## select only from those solutions that have not been selected
  ## before, i.e. those with non-negative fitness values
  n = length(N);
  available = N >= 0;
  debugprint (4, 'select, available = ', available)
  na = sum(available);
  allpossibleindices = 1:n;
  indices = allpossibleindices(available);
  switch na
    case 0
      s = 0;                  #no solutions available.
    case 1
      s = indices(1);           #only one available so select it
    case 2
      if N(indices(1)) > N(indices(2))
        s = indices(1);
      else
        s = indices(2);
      endif
    otherwise
      i = ceil(na * rand(2,1));
      i(i>na) = na;
      i(i<1) = 1;
      if N(indices(i(1))) > N(indices(i(2)))
        s = indices(i(1));
      else
        s = indices(i(2));
      endif
  endswitch
  if s>0
    debugprint(5, 'select: s N(s)', [s, N(s)])
  else
    debugprint(5, 'select: none available for selection')
  endif
endfunction

function pop = strawberry_sort (pop,nx,nz)
  [n m] = size(pop);
  ## sorting only makes sense in single criterion problems
  if nz == 1                    # single objective optimisation
    if n-nx == 2                # y + constraint violation
                                # find feasible and infeasible sets
      feasible = pop(:, pop(n,:) <= 0);
      [s, i] = sort(feasible(nx+1,:)); # sort on objective function value
      feasible = feasible(:,i);        # sorted
      ## now check for infeasible solutions
      nfeasible = size(feasible,2);
      if nfeasible < size(pop,2)
        infeasible = pop(:, pop(n,:) > 0);
        [s, i] = sort(infeasible(n,:));
        infeasible = infeasible(:,i);
        pop = [feasible, infeasible];
      else
        pop = feasible;
      endif
    else
      error('Number of objective functions does not seem to match');
    endif
  endif
endfunction

strawberry_test

function [z g] = strawberry_testobjective(x)
  z = [ sum((x-0.5).^2+1)
        sum(cos(x))];
  g = 0;
endfunction

function [x z pareto fitness] = strawberry_test
  global esfdebug
  esfdebug = 0
  global strawberry_numberofprocessors
  strawberry_numberofprocessors = 1;
  x0 = [0 0 0 0 0]';
  a = zeros(5,1);
  b = ones(5,1);
  printf('Testing bicriteria problem with %d decision variables\n', length(x0))
  printf('starting with initial point '), x0
  z0 = strawberry_testobjective(x0)
  printf('which has objective %d function values ', length(z0))
  for i=1:length(z0)
    printf('%f ', z0(i))
  endfor
  printf('\n')
  ## [x y nf ninf bestgen] = strawberry(x0, a, b, f, ngen, npop, nrmax, ns, population_strategy, output, g)
  [x z nf ninf bestgen] = strawberry(x0, a, b, @strawberry_testobjective, 200, 50);
  paretoset = z(1:2, findpareto(z(1:2,:)));
  [zz indices] = sort(paretoset(1,:));
  paretoset = paretoset(:,indices)
  plot(paretoset(1,:),paretoset(2,:),'- r',z(1,:),z(2,:),' *g')

function [z g] = strawberry_test3objective(x)
  z = [ sum((x-0.5).^2+1)
        sum(cos(x))
        sum(sin(x))];
  g = 0;
endfunction

function [x z paretoset] = strawberry_test3
  global esfdebug
  esfdebug = 0
  global strawberry_numberofprocessors
  strawberry_numberofprocessors = 1;
  global strawberry_fitness_method
  strawberry_fitness_method = 0
  global strawberry_output_population
  strawberry_output_population = 100
  x0 = [0 0 0 0 0]';
  a = zeros(5,1);
  b = ones(5,1);
  printf('Testing three criteria problem with %d decision variables\n', length(x0))
  printf('starting with initial point '), x0
  z0 = strawberry_test3objective(x0);
  printf('which has objective %d function values ', length(z0))
  for i=1:length(z0)
    printf('%f ', z0(i))
  endfor
  printf('\n')
  ## [x y nf ninf bestgen] = strawberry(x0, a, b, f, ngen, npop, nrmax, ns, population_strategy, output, g)
  [x z nf ninf bestgen] = strawberry(x0, a, b, @strawberry_test3objective, 1000, 20)
  paretoset = z(1:3, findpareto(z(1:3,:)));
  [zz indices] = sort(paretoset(1,:));
  paretoset = paretoset(:,indices)
  plot3(paretoset(1,:),paretoset(2,:), paretoset(3,:),'- r',z(1,:),z(2,:), z(3,:),' *g')