#
# Programme developpe par Guy Tremblay, Professeur au dept. d'informatique
# dans le cadre du cours INF7440.
# (Automne 2002, 2003, 2004)
# Modification par Cedric Chauve (Automne 2004).

# Ce programme renvoie un circuit non forcement optimal pour le cas du TSP geometrique
# La technique utilisee est celle du branch-and-bound.
# Le choix de l'extension du chemin courant est l'extension ayant la meilleure borne, 
# si cette borne est meilleure que le meilleur chemin deja rencontre.

# Utilisation : <nom programme> <taille graphe> <germe aleatoire>
#  taille graphe (entier) :
#    nombre de sommets
#  germe aleatoire (entier) : 
#    pour la generation du graphe aleatoire

resource TSP_BB_best()
  const int INFINI = high(int);

  #############################################################
  # Matrices des couts pour le graphe.
  #############################################################

  # Lecture et verification du nombre d'elements a generer et a traiter
  int n; getarg(1, n); 
  if (n <= 0) {
    write("*** Erreur: n <= 0" );
    stop(1);
  }

  # Lecture et verification du germe du generateur aleatoire
  int germe; getarg(2, germe); 
  if (germe <= 0) {
    write("*** Erreur: germe <= 0" );
    stop(1);
  }

  # La liste des coordonnees des sommets du graphe, que l'on suppose
  # dans le quadrant MAX x MAX du plan reel
  real MAX = 10;

  # Generation (aleatoire) des elements du tableau a : graphe geometrique
  # avec distance euclidienne.
  procedure generer( ref real a[*,*], int n ) {
    real X[n], Y[n];
    seed(1.0 * germe);
    for [i = 1 to n] {
      X[i] = real(random(0.0, MAX));
      Y[i] = real(random(0.0, MAX));
    }

    for [i = 1 to n, j = i to n] {
      a[i][j] = sqrt((X[i]-X[j])**2 + (Y[i]-Y[j])**2);
      a[j][i] = a[i][j];
    }
  }

  # Impression des elements d'un tableau.
  procedure imprimerCouts( real a[*][*], int n ) {
    for [i = 1 to n] {
      for [j = 1 to n-1] {
        writes( a[i][j], ", " );
      }
      write( a[i][n] );
    }
  }

  # La matrice d'adjacence du graphe
  real W[n][n]; 


  #############################################################
  # Types et procedures pour les sequences.
  #############################################################

  const int N = 100;
  type Sequence    = ptr SequenceRec;
  type SequenceRec = rec( int elems[N]; int ln; );

  procedure creer() returns Sequence s
  # POSTCONDITION
  #   s = []
  { s = new(SequenceRec); s^.ln = 0; }

  procedure cloner( Sequence s1 ) returns Sequence s2
  # POSTCONDITION
  #   s2 est un nouvel objet qui est une copie (profonde) de s1.
  { s2 = new(SequenceRec); s2^ = s1^; }

  procedure longueur( Sequence s ) returns int l
  # POSTCONDITION
  #   l = length(s)
  { l = s^.ln; }
 
  procedure ajouterEnQueue( Sequence s, int elem )
  # POSTCONDITION
  #   s = s' || [elem]
  { s^.ln += 1; s^.elems[s^.ln] = elem; }
  
  procedure supprimerQueue( Sequence s )
  # PRECONDITION
  #   s != []
  # POSTCONDITION
  #   s = s'[1..length(s')-1]
  { s^.ln -= 1; }
  
  procedure tete( Sequence s ) returns int t
  # PRECONDITION
  #   s != []
  # POSTCONDITION
  #   t = s[1]
  { t = s^.elems[1]; }

  procedure queue( Sequence s ) returns int q
  # PRECONDITION
  #   s != []
  # POSTCONDITION
  #   q = s[length(s)]
  { q = s^.elems[s^.ln]; }

  procedure element( Sequence s, int i ) returns int e
  # PRECONDITION
  #   i IN domain(s)
  # POSTCONDITION
  #   e = s[i]
  { e = s^.elems[i]; }

  procedure estElement( Sequence s, int e ) returns bool r 
  # POSTCONDITION
  #   r <=> e IN s
  { 
    r = false;
    for [i = 1 to s^.ln] {
      if( s^.elems[i] == e ) { r = true; return; }
    }
  }

  procedure imprimer( Sequence s ) 
  # POSTCONDITION
  #   Chacun des elements de s a ete imprime sur stdout.
  {
    writes("[");
    for [i = 1 to s^.ln-1] {
        writes( s^.elems[i], ", " );
    }
    writes( s^.elems[s^.ln], "]" );
  }

  #############################################################
  # Procedures pour le parcours du graphe. 
  #############################################################

  procedure dejaVisite( Sequence chemin, int sommet ) returns bool r
  # ROLE
  #   Determine si le sommet a deja ete visite, c'est-a-dire, 
  #   appartient deja au chemin indique.  
  # POSTCONDITION 
  #   r <=> estElement(chemin, sommmet)
  { r = estElement(chemin, sommet); }


  procedure coutChemin( Sequence chemin ) returns real cout
  # POSTCONDITION 
  #   cout = sommme des couts des aretes reliant les sommets du chemin indique.
  {
    cout = 0;
    for [i = 1 to longueur(chemin)-1] {
      cout += W[element(chemin, i)][element(chemin, i+1)];
    }  
  }

  procedure coutMin( real couts[*], Sequence chemin, int aExclure ) returns real cout
  # POSTCONDITION
  #   cout = minimum des couts indiques, mais en ignorant le sommet aExclure
  #          de meme que les sommets deja visites de chemin (mais a l'exception 
  #          du sommet de depart, pour creer le cycle).
  {
    cout = high(int);
    for [i = 1 to ub(couts) st couts[i] > 0  
                            & i != aExclure 
                            & (~dejaVisite(chemin, i) | element(chemin, 1) == i) ] {
      cout = min(couts[i], cout);
    }
  }

  procedure coutEstime( Sequence chemin, int sommet, int n ) returns real cout
  # ROLE
  #   Determiner un estime (une borne inferieure) pour le cout d'un circuit
  #   obtenu en ajoutant le sommet indique au chemin existant.
  # POSTCONDITION
  #   cout = somme des couts du chemin + somme de l'arete reliant le dernier sommet
  #          du chemin a sommet + somme des couts minimum pour 
  #          aller vers les sommets pas visites puis vers le sommet de depart.
  #
  {
    cout = coutChemin(chemin);
    cout += W[queue(chemin)][sommet];
    cout += coutMin( W[sommet][*], chemin, queue(chemin) );
    for [s = 1 to n st s != sommet & ~dejaVisite(chemin, s)] {
      cout += coutMin( W[s][*], chemin, sommet );
    }
  }

  #########################################################################
  # Procedure pour l'extension du chemin courant.
  #########################################################################

  procedure poursuivreChemin( Sequence chemin, int n, res real cout ) {
    if (longueur(chemin) == n) {
      # On a construit un circuit et on met a jour son cout
      cout = coutChemin(chemin) + W[queue(chemin)][tete(chemin)];
    } 
    else {
      real estime;          # Borne de l'enfant courant
      real opt   = INFINI;  # Borne du meilleur enfant
      int  choix = 0;       # Choix de l'enfant a parcourir

      for [s = 1 to n st ~dejaVisite(chemin, s)] {
        estime = coutEstime(chemin, s, n);
        if (estime < opt) { opt = estime; choix = s; }
      }

      ajouterEnQueue(chemin, choix); 
      poursuivreChemin(chemin, n, cout); 
    }
  }

  #########################################################################
  # Programme principal
  #########################################################################
  
  Sequence chemin = creer(); # Chemin courant
  real     cout;             # Cout du chemin courant
  int      sommetDepart;     # Premier sommet du chemin courant

  # Generation aleatoire d'un graphe et affichage
  generer( W, n );
  write( "W = " );
  imprimerCouts( W, n );
  write();

  # On explore un chemin.
  sommetDepart = 1;
  ajouterEnQueue(chemin, sommetDepart);
  poursuivreChemin(chemin, n, cout);

  # On imprime le resultat.
  imprimer( chemin ); 
  write( " =>", cout );

end