Sources de tspgen 0.32 : résolution du problème du voyageur de commerce

/*==========================================================================
Project : tspgen 0.32
File    : TSP.cpp
Purpose : The Traveling Salesman Problem Class : I/O from CSV(Alex) or TSPLIB
    EUC_2D files, distance matrix, neighbour matrix
==========================================================================*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "TSP.h"


/*========================================================================
===================== Traveling Salesman Problem =========================
==========================================================================*/


/*========================================================================
Function : TSP constructor (load a problem from a file)
==========================================================================*/
TSP::TSP(int iInstance0,
         char *sTSPfilename,
         int iNbNear0,
         TspLog * pTspLog0)
{
  iInstance = iInstance0;
  pTspLog = pTspLog0;

  pTspLog->StartFunction ("TSP::TSP(load)", iInstance, LEVEL_DEBUG, 0);

  FILE* fTSP;
  char sLine[MAX_TSP_LENLINE + 1];
  char *sElemX;
  char *sElemY;
  char *sExpX;
  char *sExpY;
  char *sPosX;
  char *sPosY;
  int iExp;
  char sMessage[255];

  int iNbCitiesRead = 0;
  aDistance = NULL;

  // lecture du fichier contenant les coordonnées du TSP (PVC)
  fTSP = fopen(sTSPfilename, "r");
  if (fTSP != NULL)
  {
    //Lecture de l'entête
    //Fichier au format "alex" (csv) ou TSPLIB ?
    fgets(sLine, MAX_TSP_LENLINE, fTSP);

    // si le fichier commence par "NAME=" c'est un TSPLIB
    if (strncmp(sLine, "NAME", 4) == 0)
    {
      pTspLog->DebugMsg("Fichier au format TSPLIB");
      bTSPLIB = true;

      // --- Récupération du nb de lignes ---
      // on recherche le marqueur DIMENSION : <nb de villes>
      do
      {
        fgets(sLine, MAX_TSP_LENLINE, fTSP);
/* DEBUG       sprintf(sMessage, "Ligne = %s", sLine);
        pTspLog->DebugMsg(sMessage);*/
     } while ( !feof(fTSP) && (strncmp(sLine, "DIMENSION", 9) != 0) );
      // le premier élément est la chaine "DIMENSION", qui ne nous intéresse pas
      sPosX = strtok(sLine, ":\n");
      // le deuxième élément est la valeur de DIMENSION, qui nous intéresse
      sPosX = strtok(NULL, ":\n");
      iNbCities = atoi(sPosX);
      sprintf(sMessage, "Nb de villes = %d", iNbCities);
      pTspLog->DebugMsg(sMessage);

      CitiesPosX = new double[iNbCities];
      CitiesPosY = new double[iNbCities];

      // --- Récupération des coordonnées des points ---
      // on recherche le marqueur NODE_COORD_SECTION
      do
      {
        fgets(sLine, MAX_TSP_LENLINE, fTSP);
/* DEBUG       sprintf(sMessage, "Ligne = %s", sLine);
        pTspLog->DebugMsg(sMessage);*/
      } while ( !feof(fTSP) && (strncmp(sLine, "NODE_COORD_SECTION", 17) != 0) );
      // le format est <position> <coordonnée X> <coordonnée Y>
      while ( !feof(fTSP) )
      {
        fgets(sLine, MAX_TSP_LENLINE, fTSP);

        if ( (sLine != NULL) && (*sLine != '\n') &&  (!feof(fTSP)) && (strncmp(sLine, "EOF", 3) != 0) )
        {
          // le premier élément est la position, on ne s'y intéresse pas
          sElemX = strtok(sLine, " \n");
          // donc on l'écrase
          sElemX = strtok(NULL, " \n");
          sElemY = strtok(NULL, " \n");

          // Gestion des chiffres avec exposants (ex : 1.309e+03)
          sPosX = strtok(sElemX, "e");

          CitiesPosX[iNbCitiesRead] = atof(sPosX);
          sExpX = strtok(NULL, "e");
          if (sExpX != NULL)
          {
            iExp=atoi(sExpX);
            CitiesPosX[iNbCitiesRead] = CitiesPosX[iNbCitiesRead] * pow(10, iExp);
/* DEBUG    printf("iExp=%d exp(iExp)=%f\n", iExp, pow(10, iExp));
            fflush(stdout);
            printf("sLine=<<%s>> sElemX=<<%s>> sPosX=<<%s>> sExpX=<<%s>> CitiesPosX[iNbCitiesRead]=%f \n", sLine, sElemX, sPosX, sExpX, CitiesPosX[iNbCitiesRead]);
            fflush(stdout);*/
          }

          sPosY = strtok(sElemY, "e");
          CitiesPosY[iNbCitiesRead] = atof(sPosY);
          sExpY = strtok(NULL, "e");
          if (sExpY != NULL)
          {
            iExp=atoi(sExpY);
            CitiesPosY[iNbCitiesRead] = CitiesPosY[iNbCitiesRead] * pow(10, iExp);
          }

          iNbCitiesRead++;
        }
      }

    }
    else
    {
      pTspLog->DebugMsg("Fichier au format TSP/CSV (Alex)");
      bTSPLIB = false;

      iNbCities = atoi(sLine);
      sprintf(sMessage, "Nb de villes = %d\n", iNbCities);
      pTspLog->DebugMsg(sMessage);

      CitiesPosX = new double[iNbCities];
      CitiesPosY = new double[iNbCities];

      // le format est <coordonnée X>;<coordonnée Y>
      while ( !feof(fTSP) )
      {
        fgets(sLine, MAX_TSP_LENLINE, fTSP);

        if ( (sLine != NULL) &&  (!feof(fTSP)) )
        {
          sPosX = strtok(sLine, ";\n");
          sPosY = strtok(NULL, ";\n");

          CitiesPosX[iNbCitiesRead] = atof(sPosX);
          CitiesPosY[iNbCitiesRead] = atof(sPosY);

          // printf("City %d : X=%f - Y=%f\n", iNbCitiesRead, CitiesPosX[iNbCitiesRead], CitiesPosY[iNbCitiesRead]);

          iNbCitiesRead++;
        }
      }

      bDefaultProblem = ((iNbCities == 250) && (CitiesPosY[iNbCities-1] == 0.306817822158337));

    }

    if (iNbCitiesRead != iNbCities)
    {
      sprintf(sMessage, "*** ERREUR : Nb de villes lu != nb villes ! (%d != %d)", iNbCitiesRead, iNbCities);
      pTspLog->ErrorMsg(sMessage);
    }

    if (iNbNear0 == -1)
    {
      iNbNear = iNbCitiesRead - 1; // nb maximum de voisins
    }
    else
    {
      iNbNear = iNbNear0;
    }

    fclose(fTSP);

    WriteToPosXY();

    FetchDistance();

    FetchNear();

    Stats();
  }
  else
  {
    pTspLog->ErrorMsg("*** ERREUR CRITIQUE : Impossible d'ouvrir le fichier tsp");
  }
  pTspLog->EndFunction ("TSP::TSP(load)", iInstance);
}


/*========================================================================
Function : TSP constructor (random generation)
==========================================================================*/
TSP::TSP(int iInstance0,
         int iNbCities0,
         int iNbNear0,
         TspLog * pTspLog0)
{
  int i;

  iInstance = iInstance0;
  pTspLog = pTspLog0;

  pTspLog->StartFunction ("TSP::TSP(rand)", iInstance, LEVEL_DEBUG, 0);

  iNbCities = iNbCities0;
  iNbNear = iNbNear0;

  CitiesPosX = new double[iNbCities];
  CitiesPosY = new double[iNbCities];

  for (i=0;i<iNbCities;i++)
  {
    CitiesPosX[i] = (rand()/(RAND_MAX + 1.0));
    CitiesPosY[i] = (rand()/(RAND_MAX + 1.0));
  }
  WriteToPosXY();

  FetchDistance();

  FetchNear();

  pTspLog->EndFunction ("TSP::TSP(rand)", iInstance);
}

/*========================================================================
Function : TSP destructor
==========================================================================*/
TSP::~TSP()
{
  pTspLog->StartFunction ("TSP::~TSP", iInstance, LEVEL_DEBUG, 0);

  delete[] CitiesPosX;
  delete[] CitiesPosY;

  delete[] aDistance;

  pTspLog->EndFunction ("TSP::~TSP", iInstance);
}

/*========================================================================
Function : Displays the position of one city (text mode)
==========================================================================*/
void TSP::Display()
{
  int i;

  for (i=0;i<iNbCities;i++)
  {
        printf("City %d : X=%.15f - Y=%.15f\n", i, CitiesPosX[i], CitiesPosY[i]);
  }
}

/*========================================================================
Function : Write the current Traveling Salesman Problem to a file
==========================================================================*/
void TSP::WriteToCSV(char *sTSPfilename)
{
  int i;
  FILE* fTSP;

  pTspLog->StartFunction ("TSP::WriteToCSV", iInstance, LEVEL_DEBUG, 0);

  fTSP = fopen(sTSPfilename, "w");
  if (fTSP != NULL)
  {
    fprintf(fTSP, "%d\n", iNbCities);
    for (i=0;i<iNbCities;i++)
    {
      fprintf(fTSP, "%.15f;%.15f\n", CitiesPosX[i], CitiesPosY[i]);
    }

    fclose(fTSP);
  }
  pTspLog->EndFunction ("TSP::WriteToCSV", iInstance);
}

/*========================================================================
Function : Write the current Traveling Salesman Problem *****
==========================================================================*/
void TSP::WriteToPosXY()
{
  int i;
  char sTemp[17];
  strcpy(sCitiesPosX, "");
  strcpy(sCitiesPosY, "");
  for (i=0;i<iNbCities;i++)
  {
    sprintf(sTemp, "%f;", CitiesPosX[i]);
    // A améliorer... enlever les 0 inutiles
    strcat(sCitiesPosX, sTemp);
    sprintf(sTemp, "%f;", CitiesPosY[i]);
    strcat(sCitiesPosY, sTemp);
  }
}

/*========================================================================
Function :
==========================================================================*/
char *TSP::getCitiesPosX()
{
  return sCitiesPosX;
}

/*========================================================================
Function :
==========================================================================*/
char *TSP::getCitiesPosY()
{
  return sCitiesPosY;
}

/*========================================================================
Function : computes the distance array by calculating the distance
              between all cities
==========================================================================*/
void TSP::FetchDistance()
{
  int i, j;

  pTspLog->StartFunction ("TSP::FetchDistance", iInstance, LEVEL_DEBUG, 0);

  aDistance = new double[iNbCities*iNbCities];

  // Calcul TSPLIB EUCL_2D
  if (bTSPLIB)
  {
    for (i=0; i<iNbCities;i++)
    {
      for (j=0; j<iNbCities;j++)
      {
        aDistance[i*iNbCities+j] = (double)
          ( (long)
            (sqrt(  (CitiesPosX[i]-CitiesPosX[j])*(CitiesPosX[i]-CitiesPosX[j])
                +  (CitiesPosY[i]-CitiesPosY[j])*(CitiesPosY[i]-CitiesPosY[j]) )
            + 0.5)
          );
      }
    }
  }
  else
  {
    for (i=0; i<iNbCities;i++)
    {
      for (j=0; j<iNbCities;j++)
      {
        aDistance[i*iNbCities+j] = sqrt(
                (CitiesPosX[i]-CitiesPosX[j])*(CitiesPosX[i]-CitiesPosX[j])
            +  (CitiesPosY[i]-CitiesPosY[j])*(CitiesPosY[i]-CitiesPosY[j])
                                        );
      }
    }
  }
  pTspLog->EndFunction ("TSP::FetchDistance", iInstance);
}

/*========================================================================
Function : Returns the distance between cities using the distance array
==========================================================================*/
double TSP::Distance(int i, int j)
{
  return aDistance[i*iNbCities+j];
}

/*========================================================================
Function : Computes the nearest city array
==========================================================================*/
void TSP::FetchNear()
{
  bool bSorted;
  int iNewPos, iNbSorted;
  int i, j, k, temp_pos;
  double d, temp_d;
  char sMessage[255];

  pTspLog->StartFunction ("TSP::FetchNear", iInstance, LEVEL_DEBUG, 0);

  double *NearBestDistance;
  NearBestDistance = new double[iNbNear];

  aNear = new int[iNbCities*iNbNear];

  // pour chaque ville
  for (i=0; i<iNbCities; i++)
  {
    if (i%100==0)
    {
      sprintf(sMessage, "Ville %d / %d", i+1, iNbCities);
      pTspLog->DebugMsg(sMessage);
    }

    iNbSorted = 0;
    // recherche des villes les plus proches
    for (j=0; j<iNbCities; j++)
    {
      // ne pas prendre la même ville
      if (i != j)
      {
        d = Distance(i, j);
        bSorted = false;
        iNewPos = 0;

        // si le tableau des villes proches était vide jusqu'à présent
        if (iNbSorted == 0)
        {
          NearBestDistance[0] = d;
          aNear[i*iNbNear+0] = j;

          iNbSorted++;
        }
        // si déjà rempli et qu'on peut l'insérer dans le groupe déjà existant
        else if (d < NearBestDistance[iNbSorted-1])
        {
          // on parcours le tableau des villes proches pour trouvé la place de l'insertion
          do
          {
            if (d < NearBestDistance[iNewPos])
            {
              if (iNbSorted < iNbNear)
              {
                iNbSorted ++;
              }

              for (k=iNbSorted - 1; k>iNewPos; k--)
              {
                temp_d = NearBestDistance[k];
                temp_pos = aNear[i*iNbNear+k];

                NearBestDistance[k] = NearBestDistance[k-1];
                aNear[i*iNbNear+k] = aNear[i*iNbNear+k-1];

                NearBestDistance[k-1] = temp_d;
                aNear[i*iNbNear+k-1] = temp_pos;
              }

              NearBestDistance[iNewPos] = d;
              aNear[i*iNbNear+iNewPos] = j;
              bSorted = true;
            }
            else
            {
              iNewPos++;
            }
          } while ( (!bSorted) && (iNewPos < iNbSorted) );
        }
        // si le groupe n'est pas complètement rempli
        else if (iNbSorted < iNbNear)
        {
          iNbSorted++;
          iNewPos =  iNbSorted - 1;

          NearBestDistance[iNewPos] = d;
          aNear[i*iNbNear+iNewPos] = j;
        }
      }
    }
  }
  delete[] NearBestDistance;

  pTspLog->EndFunction ("TSP::FetchNear", iInstance);
}

/*========================================================================
Function : Returns the j-ieth nearer city of city i
==========================================================================*/
int TSP::Near(int i, int j)
{
  if (j<=iNbNear)
  {
    return aNear[i*iNbNear+j];
  }
  else
  {
    //printf("j > iNbNear !\n");
    //fflush(stdout);
    return aNear[i*iNbNear+iNbNear];
  }
}


/*========================================================================
Function : Returns the distance (int terms of nb of cities) between 2 cities
==========================================================================*/
int TSP::IndexNear(int iCity1, int iCity2)
{
  int i, iIndex;
  bool bFound;

  i=0;
  iIndex=-1;
  bFound=false;
  while ( (i<iNbNear) && (!bFound) )
  {
    if (Near(iCity1, i)==iCity2)
    {
      iIndex=i;
      bFound = true;
    }
    i++;
  }

  return iIndex;
}

/*========================================================================
Function : guess the optimum distance of the TSP instance
==========================================================================*/
void TSP::Stats()
{
  double fDistance, fDistance1, fDistance2;
  int i;
  char sMessage[255];

  pTspLog->InfoMsg("Stats of the problem :");
  fDistance = 0.0;
  for (i=0; i<iNbCities; i++)
  {
    fDistance += Distance(i, Near(i,0));
  }
  sprintf(sMessage, "NN0 distance (nearest neigbours level 0)= %f", fDistance);
  pTspLog->InfoMsg(sMessage);

  fDistance1 = 0.0;
  for (i=0; i<iNbCities; i++)
  {
    fDistance1 += Distance(i, Near(i,1));
  }
  sprintf(sMessage, "NN1 distance (nearest neigbours level 1)= %f", fDistance1);
  pTspLog->InfoMsg(sMessage);

  fDistance2 = 0.0;
  for (i=0; i<iNbCities; i++)
  {
    fDistance2 += Distance(i, Near(i,2));
  }
  sprintf(sMessage, "NN2 distance (nearest neigbours level 2)= %f", fDistance2);
  pTspLog->InfoMsg(sMessage);

  sprintf(sMessage, "Optimum guessed distance (1/3 * NN0 + 1/3 * NN1)= %f", (fDistance/3)+(fDistance1/3)+(fDistance2/3));
  pTspLog->InfoMsg(sMessage);
}