Looking for the lost harmony

by Mithrandir

In fact, this is only a successful attempt at trying the Harmony Search for a toy problem.

I began by creating a library for random numbers, taking into account the things gathered from the Monte Carlo articles and other sources of inspiration. Though it’s not fully tested yet, I can say for sure that it is random enough for this article. And the interface gives no detail about the implementation, which is just a good thing. Here are the function headers:

 /*! initializes the random number generator with a specific seed */
 inline void init_seed(int seed);
 /*! initializes the random number generator with a time specific seed */
 inline void init();

 /*! gets a random double value in [0, 1)*/
 double random_double_unit();
 /*! gets a random double value in [0, M)*/
 double random_double(double M);
 /*! gets a random int value in [0, M)*/
 int random_int(double M);
 /*! gets a random double in range [a, b)*/
 double random_range_double(double a, double b);
 /*! gets a random int in range [a, b)*/
 int random_range_int(int a, int b);
 

Then, I took the toy problem into consideration. Not wanting to use one of those hardcore examples (for debugging purposes), I've chosen the following: let f:\{1,2,3,4,5\}^3 -> N where f(a, b, c) = |a-3| + |b-4| + |c-1| be one function. We are interested in it's minimum value.

As usual, we have to define some parameters and some datatypes and functions, before dwelling into code.

 #define K 10
 #define VS 3
 #define PHMCR 0.9
 #define PPAR 0.1
 #define ITCount 10
 #define DBGE 1
 #define DBGCount 5

 typedef struct HVector{
   int vector[VS];
   int value;
 } HVector;
 HVector memory[K];

 int value(int vector[])
 {
   return abs(vector[0]-3) + abs(vector[1] - 4) + abs(vector[2] - 1);
 }
 

Now, we can successfully write the rest of the code.

 void initialize_memory()
 {
   int i, j;
   for (i = 0; i < K; i++){//each vector
     for (j = 0; j < VS; j++)
       memory[i].vector[j] = random_range_int(1, 5);
     memory[i].value = value(memory[i].vector);
   }
   qsort(memory, K, sizeof(memory[0]), cmp);
 }

 void print_memory()
 {
   int i, j;
   printf("Memory:\n");
   for (i = 0; i < K; i++){
     printf("\t");
     for (j = 0; j < VS; j++)
       printf("%d ", memory[i].vector[j]);
     printf(" [%d]\n", memory[i].value);
   }
 }

 void evolve_once()
 {
   int newvector[VS];
   int i, val, j, k;
   float p;
   for (i = 0; i < VS; i++){
     p = random_double_unit();
     if (p < PHMCR){
       newvector[i] = memory[random_int(K)].vector[i];
       p = random_double_unit();
       if (p < PPAR){
         p = random_double_unit();
         if (p < 0.5 && newvector[i] > 1)
           newvector[i]--;
         if (p > 0.5 && newvector[i] < 4)
           newvector[i]++;
       }
     } else {
       newvector[i] = random_range_int(1, 5);
     }
   }
   val = value(newvector);
   if (val < memory[K-1].value){
     for (j = K-1; j >= 0; j--)
       if (val > memory[j].value)
         break;
     for (k = K-2; k > j; k--){
       for (i = 0; i < VS; i++)
         memory[k+1].vector[i] = memory[k].vector[i];
       memory[k+1].value = memory[k].value;
     }
     for (i = 0; i < VS; i++)
       memory[j+1].vector[i] = newvector[i];
     memory[j+1].value = val;
   }
 }

 void evolve()
 {
   int i;
   for (i = 0; i < ITCount; i++){
     evolve_once();
     if (DBGE && i && i % DBGCount == 0){
       printf("Current status %d/%d\n", i/DBGCount, ITCount/DBGCount);
       print_memory();
     }
   }
 }

 int main()
 {
   init();
   initialize_memory();
   print_memory();
   printf("Starting evolution..\n");
   evolve();
   printf("Done\n");
   print_memory();
   return 0;
 }
 

Here’s the result:

 $ make test 
 time  ./main
 Memory:
 3 4 2  [1]
 2 4 1  [1]
 2 3 1  [2]
 3 3 4  [4]
 1 2 1  [4]
 3 2 4  [5]
 1 2 2  [5]
 3 1 4  [6]
 1 2 3  [6]
 1 2 4  [7]
 Starting evolution..
 Current status 1/2
 Memory:
 2 4 1  [1]
 3 4 2  [1]
 2 4 1  [1]
 3 2 1  [2]
 2 3 1  [2]
 3 1 1  [3]
 3 2 2  [3]
 3 2 3  [4]
 3 3 4  [4]
 1 2 1  [4]
 Done
 Memory:
 3 4 1  [0]
 2 4 1  [1]
 3 4 2  [1]
 2 4 1  [1]
 2 4 2  [2]
 3 2 1  [2]
 2 3 1  [2]
 2 3 2  [3]
 3 1 1  [3]
 3 2 2  [3]
 0.00user 0.00system 0:00.00elapsed 0%CPU (0avgtext+0avgdata 0maxresident)k
 0inputs+0outputs (0major+139minor)pagefaults 0swaps
 

I didn’t implement a library for the Harmony Search because there is a need that some parameters of the problem be interleaved with some procedures of the search. Of course, we can do that via some function calls but sometimes there would be too many function calls. And I was lazy right now.

Now, seeing that the code snippets in wordpress are a headache I will turn my head to something new. Expect another article soon.

About these ads