123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161 |
- #include "bz_discr.h"
- #include <stdio.h>
- #define INT_EPS 1e-12
- // define SPAM to get trace output
- #ifdef SPAM
- #define SPIM
- #endif
- // SPIM is more benign
- int comparedim;
- double globaldiscr;
- double fmax(double a,double b)
- {
- return ((a)>(b))?(a):(b);
- }
- // argument is pointer to pointsetmember, i.e. double **.
- int cmpkeyk(double **pt1, double **pt2)
- {
- double a=(*pt1)[comparedim], b=(*pt2)[comparedim];
- if (a<b)
- return -1;
- else if (a>b)
- return 1;
- return 0;
- }
- int intpoints(double **pointset, int dim, int npoints, double *base)
- {
- int n=npoints,i,j;
- for (i=0; i<npoints; i++) {
- for (j=0; j<dim; j++) {
- if (pointset[i][j] >= (base[j]-INT_EPS)) {
- n--;
- break;
- }
- }
- }
- return n;
- }
- // Fixes one dimension at a time, each dimension defined by a point (which
- // must be on the boundary of the box). Maintains list of points that are
- // still possible (first rempoints points) and smallest base compatible with
- // earlier dimension choices (boundary points must not be excluded).
- double int_poly_discr(double **pointset, int dim, int npoints, int rempoints,
- double *base, int cdim)
- {
- double discr,maxdiscr=0.0, basecopy[dim], *ptcopy[rempoints];
- int i,j, resort=0;
- #ifdef SPAM
- fprintf(stderr, "Dim %d points %d/%d base (%g", cdim, rempoints, npoints, base[0]);
- for (i=1; i<dim; i++)
- fprintf(stderr, ", %g", base[i]);
- fprintf(stderr, ")\n");
- #endif
- if (cdim==dim) {
- discr = 1.0;
- for (i=0; i<dim; i++)
- discr *= base[i];
- maxdiscr = fmax((double)rempoints/npoints-discr,
- discr-(double)intpoints(pointset, dim, rempoints, base)/npoints);
- #ifdef SPAM
- fprintf(stderr, "Volume %g point-share %g--%g -> discr. %g\n",
- discr, (double)intpoints(pointset, dim, rempoints, base)/npoints,
- (double)rempoints/npoints, maxdiscr);
- #endif
- if (maxdiscr > globaldiscr) {
- globaldiscr=maxdiscr;
- #ifdef SPIM
- fprintf(stderr, "Improved to %g\n", globaldiscr);
- #endif
- }
- return maxdiscr;
- }
- for (i=cdim; i<dim; i++)
- basecopy[i]=base[i];
- comparedim = cdim;
- qsort(pointset, rempoints, sizeof(double *), cmpkeyk);
- #ifdef SPAM
- for (i=0; i<rempoints; i++) {
- fprintf(stderr, "Point %04d: (%g", i, pointset[i][0]);
- for (j=1; j<dim; j++)
- fprintf(stderr, ", %g", pointset[i][j]);
- fprintf(stderr, ")\n");
- }
- #endif
- for (i=0; i<rempoints; i++) {
- if (!cdim)
- fprintf(stderr, "%d/%d\n", i, rempoints);
- if (pointset[i][cdim] < base[cdim]-INT_EPS) {
- #ifdef SPAM
- fprintf(stderr, "Skip'ng %d: coord too small\n", i);
- #endif
- continue;
- }
-
- base[cdim]=pointset[i][cdim];
- for (j=cdim+1; j<dim; j++) {
- if (pointset[i][j] > base[j])
- base[j]=pointset[i][j];
- }
- if ((i && (pointset[i-1][cdim] == pointset[i][cdim])) ||
- (((i+1)<rempoints) &&
- (pointset[i+1][cdim] == pointset[i][cdim]))) {
- resort=1;
- for (j=0; j<rempoints; j++)
- ptcopy[j]=pointset[j];
- }
- j=i+1;
- while ((j < rempoints) && (pointset[j-1][cdim]==pointset[j][cdim]))
- j++;
-
- #ifdef SPAM
- fprintf(stderr, "Calling %d\n", i);
- #endif
- discr = int_poly_discr(pointset, dim, npoints, j, base, cdim+1);
- if (discr > maxdiscr)
- maxdiscr=discr;
- if (resort) {
- for (j=0; j<rempoints; j++)
- pointset[j]=ptcopy[j];
- resort=0;
- }
- for (j=cdim+1; j<dim; j++)
- base[j]=basecopy[j];
- }
- #ifdef SPAM
- fprintf(stderr, "Calling 1.0\n");
- #endif
- if (!cdim)
- fprintf(stderr, "Trying 1.0\n");
- base[cdim]=1.0;
- discr = int_poly_discr(pointset, dim, npoints, rempoints, base, cdim+1);
- for (j=cdim; j<dim; j++)
- base[j]=basecopy[j];
- if (discr > maxdiscr)
- maxdiscr=discr;
-
- return maxdiscr;
- }
- double poly_discr(double **pointset, int dim, int npoints)
- {
- double *base, discr;
- int i;
- globaldiscr=0;
- base = malloc(dim*sizeof(double));
- for (i=0; i<dim; i++)
- base[i]=0.0;
- discr = int_poly_discr(pointset, dim, npoints, npoints, base, 0);
- free(base);
- return discr;
- }
|