source: branches/HeuristicLab.Analysis.AlgorithmBehavior/liblrs/liblrs/extfunc.c @ 10200

/* HeuristicLab
* Copyright (C) 2002-2013 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
*
* This file is part of HeuristicLab.
*
* HeuristicLab is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* HeuristicLab is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
*/

/* Source based on lrslib's chdemo.c 
Copyright: David Avis 2001, avis@cs.mcgill.ca 
Licensed under the GNU GPL
*/

#include <stdio.h>
#include <string.h>
#include "extfunc.h"

 __declspec(dllexport) double calculate_volume(int numPoints, int dimension, long *num, long *den) {
  double result = 0.0;

  lrs_dic *P; /* structure for holding current dictionary and indices  */
  lrs_dat *Q; /* structure for holding static problem data             */
  lrs_mp_vector output; /* one line of output:ray,vertex,facet,linearity */
  lrs_mp_matrix Lin;    /* holds input linearities if any are found      */

  long i;
  long col;     /* output column index for dictionary            */

  printf("\nnumPoints: %d", numPoints);
  printf("\ndimension: %d", dimension);
  printf("\nnum: %d", num[11]);
  printf("\nden: %d", den[11]);
  printf("\nsizeofint: %d", sizeof(int));
  printf("\nsizeoflong: %d", sizeof(long));    

  /* Global initialization - done once */   
  if ( !lrs_init ("\nHeuristicLab lrslib wrapper"))
    return -1.0;

  /* allocate and init structure for static problem data */
  Q = lrs_alloc_dat ("LRS globals");
  if (Q == NULL)
    return -1.0;

  /* now flags in lrs_dat can be set */     
  Q->m=numPoints;           /* number of input rows = number of vertices   */ 
  Q->n=dimension;         /* number of input columns  */
  Q->hull = TRUE;     /* convex hull problem: facet enumeration      */
  Q->polytope= TRUE;  /* input is a polytope                         */
  Q->getvolume= TRUE; /* compute the volume                          */

  output = lrs_alloc_mp_vector (Q->n);

  
  P = lrs_alloc_dic (Q);   /* allocate and initialize lrs_dic      */
  if (P == NULL)
    return -1.0;
       
  /* Build polytope */   
  fill_data(P,Q, num, den);
    
  /* code from here is borrowed from lrs_main */

  /* Pivot to a starting dictionary                      */    
  if (!lrs_getfirstbasis (&P, Q, &Lin, FALSE))
    return -1.0;

  /* We initiate reverse search from this dictionary       */
  /* getting new dictionaries until the search is complete */
  /* User can access each output line from output which is */
  /* vertex/ray/facet from the lrs_mp_vector output        */
  do
  {
    for (col = 0; col <= P->d; col++)
      if (lrs_getsolution (P, Q, output, col))
        lrs_printoutput (Q, output);
  }
  while (lrs_getnextbasis (&P, Q, FALSE));
  result = get_volume(P,Q);

  lrs_printtotals (P, Q);    /* print final totals */

  /* free space : do not change order of next 3 lines! */
  lrs_clear_mp_vector (output, Q->n);
  lrs_free_dic (P,Q);           /* deallocate lrs_dic */
  lrs_free_dat (Q);             /* deallocate lrs_dat */
  lrs_close ("HeuristicLab lrslib wrapper");        

  return result;
}

double get_volume (lrs_dic *P, lrs_dat *Q) {
  double vol;

  if (Q->getvolume)
  {
    rescalevolume (P, Q, Q->Nvolume, Q->Dvolume);

    if (Q->polytope){
      prat ("\n*Volume=", Q->Nvolume, Q->Dvolume);
    } else {
      prat ("\n*Pseudovolume=", Q->Nvolume, Q->Dvolume); 
    }          

    rattodouble (Q->Nvolume, Q->Dvolume, &vol);
    return vol;
  }
  return -1.0;
}

void fill_data (lrs_dic *P, lrs_dat *Q, long *num, long *den)   
{
  int i, j;
  long *tmpNum = (long*)malloc(sizeof(long)* Q->n);
  long *tmpDen = (long*)malloc(sizeof(long)* Q->n);

  for (i=0;i<Q->m;i++)
  {     
    for(j=0;j<Q->n;j++) {
      tmpNum[j] = num[i*Q->n+j];
      tmpDen[j] = den[i*Q->n+j];
    }
    lrs_set_row(P,Q,i+1,tmpNum,tmpDen,GE);          
  }
  free(tmpNum);
  free(tmpDen);
}