HSO.c

/* 
 * 
 * Copyright (c) Year(s), 2013, Luis M. S. Russo and Alexandre 
 * P. Francisco / KDBIO / INESC-ID, <qhv@kdbio.inesc-id.pt> 
 *  
 * Any published media that is related with to use of the distributed 
 * software, or derived software, must contain a reference to "Extending 
 * quick hypervolume. Luís M. S. Russo, Alexandre P. Francisco: 
 * J. Heuristics 22(3): 245-271 (2016)".
 *  
 * Permission to use, copy, modify, and/or distribute this software for 
 * any purpose with or without fee is hereby granted, provided that the 
 * above copyright notice and this permission notice appear in all 
 * copies. 
 *  
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL 
 * WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED 
 * WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE 
 * AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL 
 * DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR 
 * PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER 
 * TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR 
 * PERFORMANCE OF THIS SOFTWARE. 
 * 
 */

#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <string.h>
#include "constants.h"
#include "HSO.h"
#include "macros.h"

/* global variables */

/* file scope macro definitions */

/* file scope type declarations */

/* file scope variables */
static int i;                 
static unsigned char p;       
static unsigned char c[D] __attribute__ ((aligned (CLS)));
static double V[D-1] __attribute__ ((aligned (CLS)));
static unsigned char S[D][N+2] __attribute__ ((aligned (CLS)));
static point P[N+2] __attribute__ ((aligned (CLS)));
static double tX[N+2] __attribute__ ((aligned (CLS)));
    
/* file scope functions declarations */

static sureInline(void) showPoint(int a);

static sureInline(void) showS();

static sureInline(void) get2D();

static sureInline(void) insert();

/* file scope functions */

static sureInline(void) showPoint(int a)
{
  int t;
  
  t = 0;
  printf("(");
  while(t < D)      {
    printf("%.1f ", P[a].x[t]);
    t++;
  }
  printf("), ");    
}

static sureInline(void) showS()
{                             
  int t;
  int cc;
    
  t = 0;
  while(t < D-1)
    {
      printf("V[%d]=%.1f ", t, V[t]);
      cc = c[t];
      while(cc <= N+1)
        {
          showPoint(S[t][cc]);
          cc++;
        }
      printf("/\n");
      t++;
    }
}

static sureInline(void) insert()
{
  unsigned char* s;
  int k, l;

  if (i < D-2) {
    s = &S[i][c[i]-1];          
    do
      {
        s++;
        *s = *(s+1);
      }while(P[p].x[i] > P[*s].x[i]); 
    *s = p;
    c[i]--;                       
    S[i][c[i]] = 0;

    V[i] = 0;
    c[i+1] = N;
    S[i+1][c[i+1]] = 0;
    S[i+1][N+1] = 1;
    if (i == D-3) {
      V[i+1] = 0;
      c[i+2] = N;
      S[i+2][c[i+2]] = 1;
      S[i+2][N+1] = 1;
    }
  } else {

    tX[p] = P[p].x[i+1];
    tX[1] = P[1].x[i+1];

    k = c[i];
    S[i][k] = S[i][k+1];
    S[i+1][k] = S[i+1][k+1];

    c[i]--;
    l = c[i];
    S[i][c[i]] = 0;
    c[i+1]--;
    S[i+1][c[i+1]] = 1;

    while (S[i+1][k] == 0 || P[p].x[i] > P[S[i][k]].x[i]) {

      if (S[i+1][k] == 1 && tX[p] > tX[S[i][k]]) {
        V[D-2] += (tX[p] - tX[S[i][k]])
                * (P[S[i][k]].x[i] - P[S[i][l]].x[i]);
    
        tX[S[i][k]] = tX[p];
      }

      if (S[i+1][k] == 1) l = k;
      k++;
      S[i][k] = S[i][k+1];
      S[i+1][k] = S[i+1][k+1];
    }

    if (S[i+1][k] == 1 && tX[p] > tX[S[i][k]]) {
      V[D-2] += (tX[p] - tX[S[i][k]])
              * (P[p].x[i] - P[S[i][l]].x[i]);

      if (P[p].x[i] == P[S[i][k]].x[i]) { 
        tX[S[i][k]] = tX[p];
        S[i+1][k+1] = 0;
      }    
    }
    
    if (P[p].x[i] <= P[S[i][k]].x[i] && tX[p] < tX[S[i][k]]) {
      tX[p] = tX[S[i][k]];
      S[i+1][k] = 0;
    } else {
      S[i+1][k] = 1;
    }
    
    S[i][k] = p;
  } 
}

static sureInline(void) get2D()
{
  double m;

  V[D-2] = 0;

  c[D-1] = N+1;
  m = P[S[D-2][c[D-1]]].x[D-1];
  while(c[D-1] > c[D-2])
    {
      c[D-1]--;
      V[D-2] += (m-P[0].x[D-1])*(P[S[D-2][c[D-1]+1]].x[D-2] - P[S[D-2][c[D-1]]].x[D-2]);
      if(P[S[D-2][c[D-1]]].x[D-1] > m)
        m = P[S[D-2][c[D-1]]].x[D-1];
    }
}

/* public functions */

void resetHSO()
{
  i = 0;
  while(i < D)
    {
      S[i][N+1] = 1;          
      c[i]      = N;
      S[i][N]   = 0;
      i++;
    }
  S[D-1][N] = 1;
}

double HSO(point* zero,  point* one, int n, int* idx, point* PS)
{
  int t;

  assert(n <= N);

  i = 0;                        
  V[i] = 0;
  c[i] = N;
  S[i][c[i]] = 0;               
  V[D-2] = 0;

  P[0] = *zero;
  P[1] = *one;
  P[1].x[D-1] = P[0].x[D-1];    
  tX[0] = P[0].x[D-1];
  tX[1] = P[1].x[D-1];

  t = n;
  while(t > 0)             
    {
      t--;
      p = t+2;
      intercept(&(P[p]), one, &(PS[idx[t]]));
      /* project(&(P[p]), to0, &one, 0, &(PS[idx[t]])); */
      insert();
    }

  do {   
    /* printf("Before Insert\n"); */
    /* showS(); */
      if(c[i] < c[i+1])         
        {                       
          while(i+2 < D)       /* Propagate elements down */
            {
              i++;
              do{
                p = S[i-1][c[i]];   
                insert();
              }while(c[i-1] < c[i] && equal(&P[p], &P[S[i-1][c[i]]], i-1) );
            }

          // get2D(); /// Calculate a 2D volume
        }
      else
        {                       
          i--;
          if(i >= 0)
            V[i] += V[i+1] * (P[S[i][c[i+1]+1]].x[i] - P[S[i][c[i+1]]].x[i]);
        }
    } while(i >= 0);     
    /* printf("End\n"); */
    /* showS(); */

  return V[0];
}
/* ------------------------------------------------------ */