function iterParada = eneReinas_escalada

  ## parámetros
  nRei = 8;           # número de reinas
  nIte = 5000;        # número maximo de iteraciones
  nPob = 100;         # tamaño de la población
  nCan = 5;           # candidatos a reproductores
  pMut = 0.8;         # probabilidad de mutación

  ## contenedores
  media = minimo = maximo = [];
  poblacion = zeros(nPob, nRei);
  adapPobla = zeros(nPob, 1);   # adaptación

  ## prepara población y la evalúa
  for i =1:nPob
    poblacion(i,:) = randperm(nRei);
    adapPobla(i)   = calculaAdaptacion(poblacion(i,:));
  endfor
  hijos = zeros(2, nRei);

  ## bucle principal
  for iter = 1:nIte
    
    ## reproductores
    iCandidatos    = ceil(rand(1,nCan) * nPob);
    [valor indice] = sort(adapPobla(iCandidatos));
    iPadres        = iCandidatos(indice(1:2));
    padres         = poblacion(iPadres,:);

    ## sobrecruce
    pos                 = ceil(rand * (nRei - 1));
    hijos(1,1:pos)      = padres(1,1:pos);
    hijos(1,pos+1:nRei) = quitar(hijos(1,1:pos),
                                 [padres(2,pos+1:nRei) padres(2,1:pos)]);
    hijos(2,1:pos)      = padres(2,1:pos);
    hijos(2,pos+1:nRei) = quitar(hijos(2,1:pos),
                                 [padres(1,pos+1:nRei) padres(1,1:pos)]);
    ## mutación
    if (rand < pMut)
      pos = ceil(rand(1,2) * nRei);
      hijos(1,pos) = hijos(1,fliplr(pos));
    endif
    if (rand < pMut)
      pos = ceil(rand(1,2) * nRei);
      hijos(2,pos) = hijos(2,fliplr(pos));
    endif

    ## escalada
    hijos(1,:) = escalada(hijos(1,:));
    hijos(2,:) = escalada(hijos(2,:));

    ## evaluación
    adapHijos = [calculaAdaptacion(hijos(1,:))
                 calculaAdaptacion(hijos(2,:))];

    ## selección
    [valor indice] = sort(adapPobla);
    poblacion      = [poblacion(indice(1:nPob-2),:)
                      hijos];
    adapPobla      = [adapPobla(indice(1:nPob-2))
                      adapHijos];
    disp([num2str(iter) ": media=" num2str(mean(adapPobla)) ...
          " min=" num2str(min(adapPobla)) ...
          " max=" num2str(max(adapPobla))]);
    
    ## dibujos
    media  = [media mean(adapPobla)];
    minimo = [minimo min(adapPobla)];
    maximo = [maximo max(adapPobla)];
    plot(media, '-r')
    hold on
    plot(minimo, '-b')
    plot(maximo, '-g')
    hold off
    drawnow

    ## criterio de parada
    [m i] = min(adapPobla);
    if m == 0
      solucion = poblacion(i,:)
      tablero  = zeros(nRei);
      tablero(sub2ind([nRei nRei], 1:nRei, solucion)) = 1
      iterParada = iter;
      return
    endif

  endfor

endfunction

## funciones auxiliares

function adaptacion = calculaAdaptacion(disposicion)
  nRei = length(disposicion);
  adaptacion = 0;
  for i = 1:nRei-1
    for j = i+1:nRei
      adaptacion += j-i == abs(disposicion(i)-disposicion(j));
    endfor
  endfor
endfunction

function quedan = quitar (quitados, conjunto)
  ## quita elementos manteniendo el orden original de los que quedan
  [quedanAscendentes indicesEnConjunto] = setdiff(conjunto, quitados);
  quedan = conjunto(sort(indicesEnConjunto));
endfunction

function optimoLocal = escalada(disposicion)
  nRei = length(disposicion);
  do
    vecinos = [];
    for i = 1:nRei-1
      for j = i+1:nRei
        vecino    = disposicion;
        aux       = vecino(i);
        vecino(i) = vecino(j);
        vecino(j) = aux;
        vecinos = [vecinos; vecino];
      endfor
    endfor
    adapVecinos = zeros(size(vecinos,1),1);
    for i = 1:size(vecinos,1)
      adapVecinos(i) = calculaAdaptacion(vecinos(i,:));
    endfor
    [valor indice] = sort(adapVecinos);
    mejora = valor(1) < calculaAdaptacion(disposicion);
    if mejora
      disposicion = vecinos(indice(1),:);
    endif
  until !mejora
  optimoLocal = disposicion;
endfunction