/*************************************************************************************************************************************/
/*                                                                                                                                   */
/*        M A X I M I S A T I O N   D E   L A   P L U S   P E T I T E   D I S T A N C E                                              */
/*        A   L ' I N T E R I E U R   D ' U N   N U A G E   D E   P O I N T S   T R I D I M E N S I O N N E L                        */
/*        E N   C O O R D O N N E E S   S P H E R I Q U E S                                                                          */
/*        P A R   U N E   M E T H O D E   D U   T Y P E   " R E C U I T   S I M U L E "                                              */
/*        A V E C   T E M P E R A T U R E  :                                                                                         */
/*                                                                                                                                   */
/*                                                                                                                                   */
/*        Author of '$xtc/recuit_si.44$c' :                                                                                          */
/*                                                                                                                                   */
/*                    Jean-Francois Colonna (LACTAMME, AAAAMMJJhhmmss).                                                              */
/*                                                                                                                                   */
/*************************************************************************************************************************************/

#include  "INCLUDES.01.I"
                                        /* Introduit le 20051116102508...                                                            */

#define   nPOINTS                                                                                                                       \
                    8                                                                                                                   \
                                        /* Nombre de points du nuage...                                                              */
#define   nITERATIONS                                                                                                                   \
                    90000                                                                                                               \
                                        /* Nombre d'iterations maximum...                                                            */

#define   RHO__MINIMUM                                                                                                                  \
                    (MIN2(dimX,dimY)/2)
#define   RHO__MAXIMUM                                                                                                                  \
                    RHO__MINIMUM
#define   RHO__DELTA_INITIAL                                                                                                            \
                    0.0
#define   RHO__DELTA_FINAL                                                                                                              \
                    0.0
                                        /* Definition de la coordonnee spherique 'rho' initialise tel que rho=constante a priori...  */
#define   THETA_MINIMUM                                                                                                                 \
                    0
#define   THETA_MAXIMUM                                                                                                                 \
                    PI
#define   THETA_DELTA_INITIAL                                                                                                           \
                    0.5
#define   THETA_DELTA_FINAL                                                                                                             \
                    0.1
                                        /* Definition de la coordonnee spherique 'theta'.                                            */
#define   PHI__MINIMUM                                                                                                                  \
                    THETA_MINIMUM
#define   PHI__MAXIMUM                                                                                                                  \
                    (2*THETA_MAXIMUM)
#define   PHI__DELTA_INITIAL                                                                                                            \
                    (2*THETA_DELTA_INITIAL)
#define   PHI__DELTA_FINAL                                                                                                              \
                    (2*THETA_DELTA_FINAL)
                                        /* Definition de la coordonnee spherique 'phi'.                                              */

#define   PROBABILITE_DE_PERTURBER_UN_POINT                                                                                             \
                    0.1                                                                                                                 \
                                        /* Probabilite de perturber le point courant...                                              */
#define   PROBABILITE_D_ACCEPTER_UNE_MAUVAISE_SOLUTION                                                                                  \
                    0.05                                                                                                                \
                                        /* Probabilite d'accepter une mauvaise solution...                                           */

#define   TEMPS                                                                                                                         \
                    (((double)(iterations-1))/((double)(nITERATIONS-1)))                                                                \
                                        /* Simulation d'un temps dans [0,1].                                                         */

extern    double    drand48();
#define   PERTURBATION_ALEATOIRE(delta_initial,delta_final)                                                                             \
                    ((drand48() - 0.5)*((delta_initial*(1-temps))+(delta_final*temps)))                                                 \
                                        /* Plus grande variation possible des coordonnees {x,y}, en rappelant que :                  */ \
                                        /*                                                                                           */ \
                                        /*                  drand48() E [0,1]                                                        */ \
                                        /*                                                                                           */ \
                                        /* d'ou la translation de 0.5 qui permet d'avoir un deplacement algebrique dans n'importe    */ \
                                        /* quelle direction et ce de facon isotrope...                                               */

typedef   struct    point     {
                              double    x;
                              double    y;
                              double    z;
                              double    rho;
                              double    theta;
                              double    phi;
                              }         point;
                                        /* Structure de definition d'un point...                                                     */
#define   X(P)                                                                                                                          \
                    (P.x)
#define   Y(P)                                                                                                                          \
                    (P.y)
#define   Z(P)                                                                                                                          \
                    (P.z)
#define   Rho_(P)                                                                                                                       \
                    (P.rho)
#define   Theta(P)                                                                                                                      \
                    (P.theta)
#define   Phi_(P)                                                                                                                       \
                    (P.phi)
                                        /* Definition de l'acces aux coordonnees d'un point.                                         */

extern    double    cos();
extern    double    sin();
#define   ConversionX(rho,theta,phi)                                                                                                    \
                    ((rho)*cos(phi)*sin(theta))
#define   ConversionY(rho,theta,phi)                                                                                                    \
                    ((rho)*sin(phi)*sin(theta))
#define   ConversionZ(rho,theta,phi)                                                                                                    \
                    ((rho)*cos(theta))
                                        /* Conversion spheriques-cartesiennes.                                                       */

extern    double    sqrt();
#define   Carre(x)                                                                                                                      \
                    ((x)*(x))                                                                                                           \
                                        /* Definition du carre de 'x'.                                                               */
#define   DistanceEuclidienne(A,B)                                                                                                      \
                    sqrt(Carre(X(B)-X(A)) + Carre(Y(B)-Y(A)) + Carre(Z(B)-Z(A)))                                                        \
                                        /* Definition de la distance entre deux points 'A' et 'B'.                                   */

#define   PERTURBATION(coordonnee,delta_initial,delta_final,minimum,maximum,periodiser)                                                 \
                    {                                                                                                                   \
                    int       perturber=VRAI;                                                                                           \
                    double    coordonnee_perturbee;                                                                                     \
                    double    temps=TEMPS;                                                                                              \
                                                                                                                                        \
                    while     (perturber == VRAI)                                                                                       \
                              {                                                                                                         \
                              coordonnee_perturbee = coordonnee + PERTURBATION_ALEATOIRE(delta_initial,delta_final);                    \
                                        /* Perturbation de la coordonnee courante...                                                 */ \
                                                                                                                                        \
                              if        (periodiser == VRAI)                                                                            \
                                        {                                                                                               \
                                        if        (coordonnee_perturbee < minimum)                                                      \
                                                  {                                                                                     \
                                                  coordonnee_perturbee = coordonnee_perturbee - minimum + maximum;                      \
                                                  }                                                                                     \
                                        else                                                                                            \
                                                  {                                                                                     \
                                                  }                                                                                     \
                                                                                                                                        \
                                        if        (coordonnee_perturbee > maximum)                                                      \
                                                  {                                                                                     \
                                                  coordonnee_perturbee = coordonnee_perturbee - maximum + minimum;                      \
                                                  }                                                                                     \
                                        else                                                                                            \
                                                  {                                                                                     \
                                                  }                                                                                     \
                                        }                                                                                               \
                              else                                                                                                      \
                                        {                                                                                               \
                                        }                                                                                               \
                                                                                                                                        \
                              if        ((coordonnee_perturbee < minimum) || (coordonnee_perturbee > maximum))                          \
                                        {                                                                                               \
                                        /* Les coordonnees hors-ecran sont rejetees...                                               */ \
                                        }                                                                                               \
                              else                                                                                                      \
                                        {                                                                                               \
                                        coordonnee = coordonnee_perturbee;                                                              \
                                        perturber=FAUX;                                                                                 \
                                        /* On s'arrete sur la premiere coordonnee perturbee situee dans l'ecran...                   */ \
                                        }                                                                                               \
                              }                                                                                                         \
                    }

#define   DISTANCE_MINIMALE(distance_minimale)                                                                                          \
                    {                                                                                                                   \
                    distance_minimale=+1e+300;                                                                                          \
                                                                                                                                        \
                    for       (i=0 ; i<nPOINTS ; i++)                                                                                   \
                              {                                                                                                         \
                              for       (j=0 ; j<nPOINTS ; j++)                                                                         \
                                        {                                                                                               \
                                        if        (i > j)                                                                               \
                                                  {                                                                                     \
                                                  double    distance_courante=DistanceEuclidienne(Lpoints[i],Lpoints[j]);               \
                                                  distance_minimale=MIN2(distance_minimale,distance_courante);                          \
                                        /* Recherche de la distance minimale courante du nuage de points...                          */ \
                                                  }                                                                                     \
                                        else                                                                                            \
                                                  {                                                                                     \
                                                  }                                                                                     \
                                        }                                                                                               \
                              }                                                                                                         \
                    }

extern    void      *malloc();

static    int       dimX=0;
#define   Xmin      0
#define   Xmax      (Xmin + (dimX-1))
                                        /* Definition des abscisses.                                                                 */
static    int       dimY=0;
#define   Ymin      0
#define   Ymax      (Ymin + (dimY-1))
                                        /* Definition des ordonnees.                                                                 */
static    int       dimZ=0;
#define   Zmin      0
#define   Zmax      (Zmin + (dimZ-1))
                                        /* Definition des profondeurs.                                                               */

#define   IMAGE(x,y)                                                                                                                    \
                    (*(image + (((((int)y)-Ymin)*dimX) + (((int)x)-Xmin))))                                                             \
                                        /* Acces a un point de l'image.                                                              */
#define   PROJECT   0.0
#define   store(n,x,y,z)                                                                                                                \
                    {                                                                                                                   \
                    int       xp=(x+(dimX/2))+((int)(PROJECT*((double)(z+(dimZ/2)))));                                                  \
                    int       yp=(y+(dimY/2))+((int)(PROJECT*((double)(z+(dimZ/2)))));                                                  \
                    int       niveau=n;                                                                                                 \
                                        /* "Pseudo-projection" du point {x,y,z}.                                                     */ \
                                                                                                                                        \
                    if        (((xp>=Xmin) && (xp<=Xmax)) && ((yp>=Ymin) && (yp<=Ymax)))                                                \
                              {                                                                                                         \
                              IMAGE(xp,yp) = MAX2(NOIR+1,MIN2(BLANC,n));                                                                \
                              }                                                                                                         \
                    else                                                                                                                \
                              {                                                                                                         \
                              }                                                                                                         \
                    }                                                                                                                   \
                                        /* Rangement d'un point valide d'une image.                                                  */
#define   initialisation(x,y)                                                                                                           \
                    {                                                                                                                   \
                    if        (((x>=Xmin) && (x<=Xmax)) && ((y>=Ymin) && (y<=Ymax)))                                                    \
                              {                                                                                                         \
                              IMAGE(x,y) = NOIR;                                                                                        \
                              }                                                                                                         \
                    else                                                                                                                \
                              {                                                                                                         \
                              }                                                                                                         \
                    }                                                                                                                   \
                                        /* Initialisation de l'image.                                                                */

main()
     {
     int                 generer_les_positions_finales=VRAI;
                                        /* Pour controler la generation des positions finales uniquement (=1) ou des positions       */
                                        /* finales et intermediaires (=0).                                                           */
     int                 sortir_l_image=VRAI;
                                        /* Pour controler la sortie (=1) ou pas (=0) de l'image...                                   */
     int                 x,y;
                                        /* Definition des coordonnees.                                                               */
     unsigned  char      *image;
                                        /* Definition de l'image a generer...                                                        */

     int       i,j,n;
     int       iterations;
                                        /* Index divers...                                                                           */
     point     Lpoints[nPOINTS];
     point     SLpoints[nPOINTS];
                                        /* Liste des points et sa Sauvegarde.                                                        */
     double    distance_minimale_avant;
     double    distance_minimale_apres;
     int       distance_minimale_apres_valide=FAUX;

     Get(dimX,"dimX");
     Get(dimY,"dimY");
     Get(dimZ,"dimZ");
     dimZ = 4*dimZ;
                                        /* Recuperation des dimensions en 'X', en 'Y' et en 'Z' de l'image a generer. On notera la   */
                                        /* multiplication par 4 de 'dimZ' car il est en general 4 fois plus petit que les autres...  */

     image=malloc(dimY*dimX);
                                        /* Definition de l'image a generer...                                                        */

     for  (y=Ymin ; y<=Ymax ; y++)
          {
          for  (x=Xmin ; x<=Xmax ; x++)
               {
               initialisation(x,y);
                                        /* Initialisation de l'image a generer...                                                    */
               }
          }

     for       (n=0 ; n<nPOINTS ; n++)
               {
               double    aveugle1=drand48();
               double    aveugle2=drand48();
               double    aveugle3=drand48();
               double    aveugle4=drand48();
               double    aveugle5=drand48();
               double    aveugle6=drand48();
                                        /* Necessaire (c'est l'experience qui le montre...).                                         */

               Rho_(Lpoints[n]) = RHO__MINIMUM;
               Theta(Lpoints[n]) = THETA_MINIMUM;
               Phi_(Lpoints[n]) = PHI__MINIMUM;
                                        /* Initialisation arbitraire des points. L'experience montre qu'une configuration            */
                                        /* initiale tres fortement symetrique est preferable a une situation aleatoire...            */
               X(Lpoints[n]) = ConversionX(Rho_(Lpoints[n]),Theta(Lpoints[n]),Phi_(Lpoints[n]));
               Y(Lpoints[n]) = ConversionY(Rho_(Lpoints[n]),Theta(Lpoints[n]),Phi_(Lpoints[n]));
               Z(Lpoints[n]) = ConversionZ(Rho_(Lpoints[n]),Theta(Lpoints[n]),Phi_(Lpoints[n]));
                                        /* Et conversion...                                                                          */
               }

     for       (iterations=1 ; iterations<=nITERATIONS ; iterations++)
               {
               for       (n=0 ; n<nPOINTS ; n++)
                         {
                         Rho_(SLpoints[n]) = Rho_(Lpoints[n]);
                         Theta(SLpoints[n]) = Theta(Lpoints[n]);
                         Phi_(SLpoints[n]) = Phi_(Lpoints[n]);
                         X(SLpoints[n]) = X(Lpoints[n]);
                         Y(SLpoints[n]) = Y(Lpoints[n]);
                         Z(SLpoints[n]) = Z(Lpoints[n]);
                                        /* Sauvegarde de l'etat courant des points...                                                */
                         }

               if        (distance_minimale_apres_valide == FAUX)
                         {
                         if        (iterations == 1)
                                   {
                                   DISTANCE_MINIMALE(distance_minimale_avant);
                                        /* Recherche de la distance minimale avant la perturbation.                                  */
                                   }
                         else
                                   {
                                        /* Dans ce cas 'distance_minimale_avant' est encore valide (calculee lors d'une iteration    */
                                        /* anterieure...).                                                                           */
                                   }
                         }
               else
                         {
                         distance_minimale_avant = distance_minimale_apres;
                         }

               if        (iterations == 1)
                         {
                         fprintf(stderr,"distance minimale initiale=%f\n",distance_minimale_avant);
                         }
               else
                         {
                         }

               for       (n=0 ; n<nPOINTS ; n++)
                         {
                         if        (drand48() < PROBABILITE_DE_PERTURBER_UN_POINT)
                                   {
                                   double    Rho__perturbe=Rho_(Lpoints[n]);
                                   double    Theta_perturbe=Theta(Lpoints[n]);
                                   double    Phi__perturbe=Phi_(Lpoints[n]);

                                   PERTURBATION(Rho__perturbe,RHO__DELTA_INITIAL,RHO__DELTA_FINAL,RHO__MINIMUM,RHO__MAXIMUM,FAUX);
                                   PERTURBATION(Theta_perturbe,THETA_DELTA_INITIAL,THETA_DELTA_FINAL,THETA_MINIMUM,THETA_MAXIMUM,FAUX);
                                   PERTURBATION(Phi__perturbe,PHI__DELTA_INITIAL,PHI__DELTA_FINAL,PHI__MINIMUM,PHI__MAXIMUM,FAUX);

                                   Rho_(Lpoints[n]) = Rho__perturbe;
                                   Theta(Lpoints[n]) = Theta_perturbe;
                                   Phi_(Lpoints[n]) = Phi__perturbe;
                                        /* Perturbation aleatoire des points.                                                        */
                                   X(Lpoints[n]) = ConversionX(Rho_(Lpoints[n]),Theta(Lpoints[n]),Phi_(Lpoints[n]));
                                   Y(Lpoints[n]) = ConversionY(Rho_(Lpoints[n]),Theta(Lpoints[n]),Phi_(Lpoints[n]));
                                   Z(Lpoints[n]) = ConversionZ(Rho_(Lpoints[n]),Theta(Lpoints[n]),Phi_(Lpoints[n]));
                                        /* Et conversion...                                                                          */
                                   }
                         else
                                   {
                                   }
                         }

               DISTANCE_MINIMALE(distance_minimale_apres);
                                        /* Recherche de la distance minimale apres la perturbation.                                  */

               if        (    (distance_minimale_apres < distance_minimale_avant)
                         ||   (drand48() < PROBABILITE_D_ACCEPTER_UNE_MAUVAISE_SOLUTION)
                          )
                                        /* On cherche a maximiser la distance minimale. On notera le test '<' (alors que la          */
                                        /* logique voudrait '<=') a cause du cas ou les points occuperaient initialement tous        */
                                        /* la meme position (d'ou : distance_minimale_avant=0) : dans ce cas, comme un certain       */
                                        /* nombre de points ne sont pas perturbes, on a donc distance_minimale_apres=0, d'ou avec    */
                                        /* avec le test '<=' le rejet de cette solution ; ainsi, on risque d'etre bloque             */
                                        /* indefiniment puisque l'on restaure alors l'etat precedent...                              */
                         {
                         distance_minimale_apres_valide = FAUX;
                                        /* L'etat anterieur est inchange...                                                          */

                         for       (n=0 ; n<nPOINTS ; n++)
                                   {
                                   Rho_(Lpoints[n]) = Rho_(SLpoints[n]);
                                   Theta(Lpoints[n]) = Theta(SLpoints[n]);
                                   Phi_(Lpoints[n]) = Phi_(SLpoints[n]);
                                   X(Lpoints[n]) = X(SLpoints[n]);
                                   Y(Lpoints[n]) = Y(SLpoints[n]);
                                   Z(Lpoints[n]) = Z(SLpoints[n]);
                                        /* Restauration de l'etat anterieur lorsqu'il y a degradation des performances...            */
                                   }
                         }
               else
                         {
                         distance_minimale_apres_valide = VRAI;
                                        /* On conserve un etat perturbe qui maximise la distance minimale...                         */
                         }

               if        (generer_les_positions_finales == 1)
                         {
                         }
               else
                         {
                         for       (n=0 ; n<nPOINTS ; n++)
                                   {
                                   store((BLANC*iterations)/nITERATIONS,X(Lpoints[n]),Y(Lpoints[n]),Z(Lpoints[n]));
                                        /* Marquage des points...                                                                    */
                                   }
                         }

               if        (iterations == nITERATIONS)
                         {
                         fprintf(stderr,"distance minimale finale..=%f\n",distance_minimale_avant);
                                        /* ATTENTION : c'est la distance minimale a l'avant-derniere iteration qui est editee car,   */
                                        /* en effet, il n'est pas sur que 'distance_minimale_apres' etait maximale...                */
                         }
               else
                         {
                         }
               }

     if        (generer_les_positions_finales == 1)
               {
               for       (n=0 ; n<nPOINTS ; n++)
                         {
                         store((BLANC*(Z(Lpoints[n])+(dimZ/2))/dimZ),X(Lpoints[n]),Y(Lpoints[n]),Z(Lpoints[n]));
                                        /* Marquage des points...                                                                    */
                         }
               }
     else
               {
               }

     if        (sortir_l_image == 1)
               {
               write(1,image,dimX*dimY);
                                        /* Sortie de l'image...                                                                      */
               }
     else
               {
               }
     }



Copyright © Jean-François Colonna, 2021-2023.
Copyright © CMAP (Centre de Mathématiques APpliquées) UMR CNRS 7641 / École polytechnique, Institut Polytechnique de Paris, 2021-2023.