volume.c

#include <stdlib.h>
#include <ctype.h>
#include <math.h>
/*--------------------------------------------------------------

    This file contains six recursive C routines for calculating
    the volume and derivatives of a convex polyhedron in n 
    dimensions. The routines volume, volumef, volumeb and dvda are 
    fortran callable interfaces to their C counterparts cvolume,
    cvolumeb and cdvda.

    ROUTINE      COMMENT			CALLS 

    cvolume	 calculates volume of region    cvolume 
                 which satisfies Ax <= b.

    cvolumef	 calculates volume of region    cvolumef 
                 which satisfies Ax <= b.
                 (faster version, see below)

    cvolumeb     calculates volume and the      cvolume
                 derivative d(vol)/db(0).

    cdvda        calculates derivatives		cdvda, cvolumeb, 
                 d(vol)/da(0,j) (j=1,n).        & cvolumebj

    cdvdaf        calculates derivatives	cdvdaf, cvolumeb, 
                 d(vol)/da(0,j) (j=1,n).        & cvolumebj
                 (faster version, see below)

    cvolumebj    calculates partial volumes
                 and derivatives d(vol)/d(b(j)
                 (j=1,n). (called by cdvda.)

    Here b(j) is the jth element of vector b, and a(i,j) is the
    (i,j)th coefficient of the matrix A corresponding to the
    ith constraint and the jth dimension.

    The two faster versions cvolumef and cdvdaf do not continue down
    the recursive tree if b(j) = 0 for any j at any time. This avoids
    extra work but restricts the applicability of the routines to
    cases where the origin does not pass through any inconsistent
    constraints.
     
				Malcolm Sambridge, April 1995.

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

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

			ROUTINE: cvolume

    This routine calculates the `volume' in dimension n of the region
    bounded by a set of m linear inequality constraints of the form
    A x <= b, where a has m rows and n columns and is given by a(m,n),
    b is the n-vector and is contained in b(m). The recursive formula
    of Lasserre (1983) is used. Redundant constraints are allowed 
    If the inequality constraints are inconsistent then the volume 
    is returned as zero.  If the original polyhedron is unbounded 
    then a warning is issued and the volume is return as -1.0 
    (even though it is undefined). 

			Jean Braun and Malcolm Sambridge, Jan 1995.

    Lasserre, J. B., 1983. "An analytical expression and an Algorithm
    for the Volume of a Convex Polyhedron in R^n.", JOTA, vol. 39, 
    no. 3, 363-377.

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

float cvolume (a,b,m,n,mmax,nmax)

int   *n, *m, *mmax, *nmax;
float *a, *b; 

{

float v,amax,pivot;
int   i,j,t,k,l;
int   jj,kk;
float   *ai, *aj, *ajt, *apjj, *bi;
int   kmm,tmm;
int   nn,mm,nn_max,mm_max;
int   nm1,mm1;
float  *ap, *bp;
int   firstmin,firstmax;
float bmin,bmax,bb;
float partialv;

nn= *n;
mm= *m;
nn_max= *nmax;
mm_max= *mmax;

/* one-dimensional case (full reduction) */

if (nn == 1) 
  {
  firstmin=0;
  firstmax=0;
  for (l=0;l<mm;l++)
    {
    if ( *(a+l) > 0.) 
      {
      bb= *(b+l)/ *(a+l);
      if (firstmin==0) {firstmin=1;bmin=bb;}
      else if (bb<bmin) bmin=bb;
      }
    else if ( *(a+l) < 0.)
      {
      bb= *(b+l)/ *(a+l);
      if (firstmax==0) {firstmax=1;bmax=bb;}
      else if (bb>bmax) bmax=bb;
      }
    else if ( *(b+l) < 0.)         /* Constraints are inconsistent */
      { 
       printf("Inconsistent constraints found after reduction to n = 1 \n");
       return(0.); 
      }
    }
  v=0.;
  if (firstmin*firstmax == 1) v=bmin-bmax;
  else 
    {
/*     printf("Volume is unbounded at 1-D; volume returned as -1\n");*/
     return(-1.0);
    }

  if (v<0.) {v=0.;} 

  return(v);
  }

nm1=nn-1;
mm1=mm-1;
v=0.;

for (i=0;i<mm;i++)
  {
   ai=a+i;
   bi=b+i;

/* find largest pivot */

  amax=0.;
  for (j=0;j<nn;j++) 
    if (fabs( *(ai+j*mm_max)) >= amax) {amax= fabs( *(ai+j*mm_max)); t=j;}
  tmm=t*mm_max;
  pivot=*(ai+tmm);

/* finds contribution to v from this pivot (if not nil) */

  if (amax == 0.)
  {
                                 /* Constraint is inconsistent */

    if ( *(bi) < 0.) 
      { printf("Constraint %d is inconsistent\n",i+1); return(0.); }

                                 /* otherwise constraint is redundant */

    printf("Constraint %d is redundant\n",i+1); 
  }
  else
  {

/* allocate memory */

  ap = (float *) malloc((sizeof *a)*nm1*mm1);
  bp = (float *) malloc((sizeof *b)*mm1);

/* reduce a and b into ap and bp eliminating variable t and constraint i */

  jj=-1;
    for (j=0;j<mm;j++)
      if (j != i)
      {
      jj=jj+1;
      aj=a+j;
      ajt=aj+tmm;
      *(bp+jj)= *(b+j) - *(bi) * *(ajt) / pivot;
      apjj=ap+jj;
      kk=-1;
      for (k=0;k<nn;k++)
        if (k != t)
        {
        kk=kk+1;
        kmm=k*mm_max;
        *(apjj+kk*mm1)= *(aj+kmm)- *(ajt) * *(ai+kmm)/ pivot;
        }
      }
  
/* add contribution to volume from volume calculated in smaller dimension */

  partialv=cvolume(ap,bp, &mm1, &nm1, &mm1, &nm1);
  if(partialv == -1.0)return(-1.0);
  v=v+ *(bi)/amax*partialv/nn; 

  free(ap);
  free(bp);