## a proximity search for creating groups in a network
## aplied for the Portuguese 2020/21 3rd level male football 11 competition

clear all
##TODO:20200701:ricardo-andre:generalise, generalise, generalise!
##TODO:20200701:ricardo-andre:add a finetuning 2-3-4way systematic switch


## teamData
## name, coordinates, region, series
teamData = ...
{
 "1º Dezembro", [38.7882, -9.3753], "C", "D";
 "Aljustrelense", [37.8858, -8.1620], "C", "D";
 "Alverca", [38.8979, -9.0350], "C", "D";
 "Amarante", [41.2821, -8.0757], "C", "B";
 "Amora", [38.6331, -9.1166], "C", "D";
 "Anadia", [40.4364, -8.4326], "C", "C";
 "Armacenenses", [37.1083, -8.3653], "C", "D";
 "Arouca", [40.9329, -8.2503], "C", "B";
 "Beira-Mar", [40.6178, -8.5937], "C", "C";
 "Benfica Castelo Branco", [39.8353, -7.4951], "C", "C";
 "Berço", [41.4301, -8.3210], "C", "A";
 "Bragança", [41.8034, -6.7703], "C", "A";
 "Caldas", [39.4033, -9.1265], "C", "C";
 "Canelas 2010", [41.0818, -8.6057], "C", "B";
 "Castro Daire", [40.9018, -7.9284], "C", "B";
 "Cerveira", [41.9489, -8.7423], "C", "A";
 "Chaves Satélite", [41.7506, -7.4649], "C", "A";
 "Coimbrões", [41.1276, -8.6304], "C", "B";
 "Condeixa", [40.1179, -8.4919], "C", "C";
 "Câmara de Lobos", [41.2421, -8.6786], "M", "A";
 "Esp. Lagos", [37.1166, -8.6778], "C", "D";
 "Fabril Barreiro", [38.6577, -9.0539], "C", "D";
 "Fafe", [41.4471, -8.1691], "C", "A";
 "Felgueiras 1932", [41.3560, -8.1906], "C", "B";
 "Fontinhas", [38.7756, -9.1354], "A", "C";
 "Fátima", [39.6035, -8.6633], "C", "C";
 "Gin. Figueirense (C.Rodrigo)", [40.9023, -6.9609], "C", "B";
 "Gondomar", [41.1359, -8.5418], "C", "B";
 "Ideal", [38.7756, -9.1354], "A", "C";
 "Juv. Pedras Salgadas", [41.5463, -7.6138], "C", "A";
 "Leça", [41.2033, -8.6888], "C", "B";
 "Louletano", [37.0883, -7.9747], "C", "D";
 "Loures", [38.8277, -9.1649], "C", "D";
 "Lusit. Évora", [38.5686, -7.8908], "C", "D";
 "Lusitano Vildemoinhos", [40.6582, -7.9257], "C", "B";
 "Lusitânia de Lourosa", [40.9878, -8.5436], "C", "B";
 "Maria da Fonte", [41.5695, -8.2642], "C", "A";
 "Marinhense", [39.7366, -8.9312], "C", "C";
 "Marítimo B", [41.2421, -8.6786], "M", "A";
 "Merelinense", [41.5767, -8.4581], "C", "A";
 "Mirandela", [41.4901, -7.1767], "C", "A";
 "Montalegre", [41.8244, -7.8051], "C", "A";
 "Oleiros", [39.9231,-7.9124], "C", "C";
 "Olhanense", [37.0292, -7.8486], "C", "D";
 "Oliv. Hospital", [40.3229, -7.8557], "C", "C";
 "Oliveirense", [41.4042, -8.4162], "C", "A";
 "Olímpico Montijo", [38.7134, -8.9676], "C", "D";
 "Oriental", [38.7439, -9.1098], "C", "D";
 "Paredes", [41.2119, -8.3499], "C", "B";
 "Pedras Rubras", [41.2571, -8.6780], "C", "B";
 "Pinhalnovense", [38.6347, -8.9155], "C", "D";
 "Praiense", [38.7756, -9.1354], "A", "C";
 "Real", [38.7520, -9.2697], "C", "D";
 "Rec. Águeda", [40.5624, -8.4417], "C", "C";
 "Sacavenense", [38.7898, -9.1030], "C", "D";
 "Sanjoanense", [40.8982, -8.4999], "C", "B";
 "Sertanense", [39.7975, -8.0967], "C", "C";
 "Sintra Football", [38.7006, -9.3039], "C", "D";
 "Sintrense", [38.7994, -9.3769], "C", "D";
 "Sp. Braga B", [41.5031, -8.7690], "C", "A";
 "Sp. Espinho", [40.9581, -8.5563], "C", "B";
 "São Martinho", [41.3604, -8.3682], "C", "A";
 "Torreense", [39.0949, -9.2567], "C", "C";
 "Trofense", [41.3329, -8.5644], "C", "B";
 "Un. Leiria", [39.7488, -8.8129], "C", "C";
 "Un. Madeira", [41.2421, -8.6786], "M", "A";
 "Un. Santarém", [39.2296, -8.6927], "C", "C";
 "Valadares Gaia", [41.0977, -8.6398], "C", "B";
 "Vila Real", [41.2876, -7.7557], "C", "B";
 "Vit. Guimarães B", [41.4488, -8.2798], "C", "A";
 "Vit. Sernache",[39.8069, -8.1932], "C", "C";
 "Vizela", [41.3888, -8.3074], "C", "A"
};
## {"Loures", [38.8277, -9.1649], "C", "D"},
##   has played at Zambujal:  {"Loures", [38.8697, -9.1211], "C", "D"},
## {"Sintra Football", [38.7006, -9.3039], "C", "D"},
##   may move to Amadora:  {"Sintra Football", [38.7520, -9.2279], "C", "C"},
## Açores inserted at Lisboa (airport) 38.7756, -9.1354
##   "Fontinhas", [38.7409, -27.1045], "A", "C";
##   "Ideal", [37.8186, -25.5301], "A", "C";
##   "Praiense", [38.7321, -27.0700], "A", "C";
## Madeira inserted at Porto (airport) 41.2421, -8.6786
##   "Câmara de Lobos", [32.6551, -16.9823], "M", "A";
##   "Marítimo B", [32.6706, -16.9362], "M", "A";
##   "Un. Madeira", [32.6494, -16.9004], "M", "A";

## coordinates and sort north to south
coords = cell2mat(teamData(:,2));
[~, i] = sort(-coords(:,1));
teamData = teamData(i,:);
coords = coords(i,:);

## point to point distances
## stable expression for "small" distances
## https://en.wikipedia.org/wiki/Great-circle_distance
phi = coords(:,2)*pi/180;
lambda = coords(:,1)*pi/180;
dPhi = phi-phi';
dLambda = lambda-lambda';
R = 6378;
global dist = (R*2*asin(sqrt(sin(dPhi/2).^2+cos(phi).*cos(phi').*sin(dLambda/2).^2))).^1;
##TODO:20200629:ricardo-andré:(rough) corretion for Açores and Madeira
## get distance through Porto/Lisbon(/Faro?) whatever is shortest (both ways)
##  (or only from mainland to islands, but this distance is not metric - right?)


function d = totalDist (list)
  ## total distance
  global dist;
  d = 0;
  for g = 1:4
    ## groups are 1:18,19:36,39:54,55:72, i.e. (18*g-17):(18*g), g=1:4
    grp = list((18*g-17):(18*g));
    d += sum(dist(grp,grp)(:));
  endfor
endfunction

function output (bestGrp, coords, teamData)
  M=[];
  for g = 1:4
    grp = sort(bestGrp((18*g-17):(18*g)));
    M=[M;grp];
    yy{g}=coords(grp,1);
    xx{g}=coords(grp,2);
  endfor
  ## print groups
  teamData(M')
  ## plot groups
  plot(xx{1}, yy{1}, ":*k",
       xx{2}, yy{2}, ":or",
       xx{3}, yy{3}, ":vb",
       xx{4}, yy{4}, ":sg")
  axis("image");
  axis([-9.5, -6.6, 36.8, 42.2]);
  drawnow () # force update
  ## save plot
  fname = ["map" num2str(time()) ".png"]
  print (fname)
endfunction


## initial population (~n*log2(n))
pMax = 72 * 6;
bestDist = Inf;
bestGrp = zeros(1,72);
pop = {};
tdist = [];
for p = 1:pMax
  pop{p} = randperm(72);
  tdist(p) = totalDist(pop{p});
  if (tdist(p) < bestDist)
    bestGrp = pop{p}
    output (bestGrp, coords, teamData)
    bestDist = tdist(p)
  endif
endfor


## new generations
kNoImprovement = 0;
while (kNoImprovement<10*pMax)
  ## reproduction
  newPop = {};
  newTdist = [];
  for p=1:pMax

    ## parents
    for r=1:2
      candidates = randperm(pMax,10);
      [~, idx] = min(tdist(candidates));
      idx=idx(1); #there can be only one!
      parent{r}=pop{candidates(idx)};
    endfor

    ## cross
    ## inspired by https://www.youtube.com/watch?v=a1DUUnhk3uE?t=348
    ## indices yet to fill
    ii = randperm(72);
    ## initial parents to copy (p0) and check (p1)
    r0 = randi(2,1,1);
    r1 = 3-r0;
    newPop{p}=zeros(1,72);
    while (size(ii))
      i=ii(1);
      ii(1)=[];
      newPop{p}(i) = parent{r0}(i);
      while (isempty(find(parent{r1}(i)==newPop{p})))
	## the team in the other parent is not yet used
	## so, keep it in the location of the parent p0
	## (so that it stays in one of the two previous locations)
	## index, on p0, of the value in p1
	i = find(parent{r0}==parent{r1}(i));
	## clear the index from the todo list
	ii(find(i==ii))=[];
	newPop{p}(i) = parent{r0}(i);
      endwhile
      ## switch parents, and repeat
      r0=3-r0;
      r1=3-r1;
    endwhile

    ## mutate
    ##TODO:20200701:ricardo-andre:threeway mutation
    for i = 1:min(randp(2,1,1),72)
      i1 = randi(72,1,1);
      i2 = randi(72,1,1);
      x = newPop{p}(i1);
      newPop{p}(i1) = newPop{p}(i2);
      newPop{p}(i2) = x;
    endfor

    ## distance
    newTdist(p) = totalDist(newPop{p});
    if (newTdist(p) < bestDist)
      bestGrp = newPop{p}
      output (bestGrp, coords, teamData)
      bestDist = newTdist(p)
      kNoImprovement = 0;
    else
      ++kNoImprovement;
    endif
  endfor
  ## replace the population
  pop=newPop;
  tdist=newTdist;
endwhile


## summary (org-mode format)
printf("| name | d0 | d1 | d2 | gain | g0 | g1 | g2 | coord |\n")
for t=1:72
  ## official group
  g0 = teamData{t,4};
  grp0=find([teamData{:,4}]==g0);
  d0 = sum(dist(t,grp0));
  ## my north to south group (note that data is sorted north to south
  g1 = 1 + floor((t-1)/18);
  grp1 = (18*g1-17):(18*g1);
  d1 = sum(dist(t,grp1));
  ## ideal group
  g2 = 1 + floor((find(bestGrp==t)-1) / 18);
  grp2 = bestGrp((18*g2-17):(18*g2));
  d2 = sum(dist(t,grp2));
  ## print it
  printf("| %s | %7.1f | %7.1f | %7.1f | %7.1f | %c | %d | %d | %f |\n",
	 teamData{t,1}, d0, d1, d2, d1-d2, g0, g1, g2, teamData{t,2}(1) )
endfor