inexclusion.c

Go to the documentation of this file.

/* 
 * 
 * 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 <strings.h>
#include "point.h"
#include "quickhvolume.h"
#include "macros.h"
#include "subsets.h"
#include "inexclusion.h"

/*** global variables *************************************/
    
/*** file scope macro definitions *************************/

/*** file scope type declarations *************************/
    
/*** file scope variables *********************************/

static point PTS[1<<N] __attribute__ ((aligned (CLS))); 
static int C2P[1<<N] __attribute__ ((aligned (CLS)));   
static int C[N][N+1] __attribute__ ((aligned (CLS)));   
/*** file scope functions declarations ********************/

static double InExClusion(point* zero, point* one, int n, int* idx, point* P);

/*** public functions declarations ************************/

/*** file scope functions *********************************/

/* ------------------------------------------------------ */

/*** public functions *************************************/

void IECinit(void)
{
  int i;
  int j;

  i = 0;
  while(i < N)
    {
      C[i][0] = 1;
      C[i][i] = 1;
      j = 1;
      while(j < i)
        {
          C[i][j] = C[i-1][j-1] + C[i-1][j];
          j++;
        }
      i++;
    }

  i = 0;
  while(i < N)
    {
      j = 1;
      while(j <= i)
        {
          C[i][j] += C[i][j-1];
          j++;
        }
      i++;
    }
}

static double InExClusion(point* zero, point* one, int n, int* idx, point* PS)
{
  assert(n <= N);

  double res;
  int i;
  int j;
  int twoton;
  int ones;                     
  int counts[N+2];              
  res = 0;
  
  i = 0;
  j = 1;
  while(i < n)
    {
      intercept(&(PTS[i+1]), one, &(PS[idx[i]]));
      i++;
      res += HV(zero, &(PTS[i]));
      C2P[j] = i;
      j <<= 1;
    }

  twoton = 1 << n;              
  memcpy(&counts[1], &C[n][0], (n+1)*sizeof(int)); 
  i = 1;
  ones = 1;
  while(i < twoton)             
    {
      if(ones > 1)
        {                       
          C2P[i] = counts[ones];
          counts[ones]++;
          intercept(&(PTS[C2P[i]]), &(PTS[C2P[Rbit(i)]]), &(PTS[C2P[Rpop(i)]]));
          if(ones & 1)          
            res += HV(zero, &(PTS[C2P[i]]));
          else                  
            res -= HV(zero, &(PTS[C2P[i]]));
        }
      i++;
      ones++;
      ones-=fastLog2(Rbit(i));
    }

  return res;
}

/* double exInExClusion(point* zero, point* one, int n, int* idx, point* PS, double* exHV) */
/* { */
/*   assert(n <= N); */

/*   int j;                        /\**< Index value *\/ */
/*   double T;                     /\**< Total volume *\/ */

/*   T = InExClusion(zero, one, n, idx, PS); */

/*   j = 0; */
/*   while(j < n) */
/*     { */
/*       swap(idx, 0, j); */
/*       exHV[idx[0]] += T - InExClusion(zero, one, n-1, &(idx[1]), PS); */
/*       j++; */
/*     } */

/*   return T; */
/* } */

double exInExClusion(point* zero, point* one, int n, int* idx, point* PS, double* exHV)
{
  assert(n <= N);

  int i;
  int j;
  int twoton;
  int ones;                     
  double T;                     
  double H;                     
  T = InExClusion(zero, one, n, idx, PS);

  twoton = 1 << n;              
  i = 1;
  ones = 1;
  while(i < twoton)             
    {
      j = i;
      H = HV(zero, &(PTS[C2P[i]]));
      if(0 == (ones & 1))          
        H = -H;

      while(0 != j)
        {
          exHV[idx[fastLog2(Rbit(j))]] += H;
          j = Rpop(j);
        }
      i++;
      ones++;
      ones-=fastLog2(Rbit(i));
    }

  return T;
}