/* 
   mostats5.c 

   compare pareto sets. assumes minimising on each objective 
   By David Corne, Joshua Knowles (C) 1999
 
   Copyright (C) 2000, David Corne and Joshua Knowles

   This program 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 2 of the License, or
   (at your option) any later version. 

   This program 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. 

   The GNU General Public License is available at:
      http://www.gnu.org/copyleft/gpl.html
   or by writing to: 
        The Free Software Foundation, Inc., 
        675 Mass Ave, Cambridge, MA 02139, USA.  



   Updated to include a shift of all points so that highly 
   convex/concave spaces do not result in distorted statistics.

 
   COMPILATION 
 
   Compile with gcc. Eg:

   gcc mostats.c - mostat -lm



   USAGE

   mostat <file> <nobjectives> <gridpoints> <output> <shift>

   Eg: 

   mostat thisfile 5 20 1 1000.0

   does a statistical comparison of the results included in "thisfile"
  which were obtained on a problem where fitness has 5 objectives. The
  comparison is based around sampling the N-dimensional fitness space
  (5 dimensional in this case) on the basis of using 20 gridpoints per 
  co-ordinate. See below for a fuller explanation, but, in general, 
  where g is gridpoints and d is number of objectives, the number of
  lines drawn to sample the tradeoff surfaces is:

     dg^(d-1)

  This is because a line is drawn from the origin to each point on a 
  g by g by ... <d-1 times> g grid for each of the d dimensions.

  ARGUMENTS

  <file>  This is simply the name of the file which contains your results.
          The structure of the file is like this:

           RESULTS FOR ALGORITHM 1
           2000000
           RESULTS FOR ALGORITHM 2
           2000000
           RESULTS FOR ALGORITHM 3
           ...
           2000000

           Ie: it contains results for 1 or more algorithms in sequence
         The output implicitly numbers the algorithms in the order
         they appear in this file. Separating each set is a single line
         containing only the number "2000000". The file must also end with
         such a line.

          Each RESULTS ... set contains results in the following form
 
          SET1
          1000000
          SET2
          1000000
          SET3
          ...
          SETN

           Ie: N (any number, can be different for each algorithm)
          sets of results obtained in different trials of applying the
          algorithm in question to the problem. Sets must have the separator
          "1000000" on a line by itself between them. The sets for an 
          algorthm do not end with this, however,. INstead, after the last 
          set for an algorithm, the 2000000 separator will appear (see above)
          either marking the beginning of sets of results for the next
          algorithm, or marking the end o fthe file.

         Each SET is simply the collection of fitness vectors, one line
         each, returned by a single trial of the algorithm on the problem.
         There can be any number of such points in any set, but each vector
         must of coure contain precisely D numbers, where D is the number of
         objectives.

         Eg: if it was a 4 objective problem, then the following is a
         possible SET:

         2 3.6 199 -0.003 
         6 8 10 -12
         0.4 0.3 0.3 0.2
         11 10 -23.34 22.1
         0 0 1000.0 3

         It doesn't matter if a set includes points which are dominated
        by others within the set. Such points make no difference to
        the statistical analysis.

        EXAMPLE FILE:

        This is a possible file in which three algorithms are compared on
       a 3-objective problem, where, for each algorithm, there were just
       2 trial runs. This is obviously a tiny example, and it will normally
       be better to do at least 20 runs for each algorithm to get conclusive
       results. Different numbers of vectors are returned in different
       trials. This might be the case if you wanted to filter only the
       non-dominate dones from the population, or whatever. ALthough that
       doesn't matter to this code.


       ------------------ cut here ---------------------

       23 10 2
       11 -10 6
       37 2 10
       1000000
       8 5 19
       2000000
       16 -3.06 20
       -11 14 4
       -100 100 2
       3 2 1
       1 2 3 
       1000000
       11 90 -30
       12 70 -90
       2000000
       12 12 -4
       -8 8 9
       1000000
       13 14 10
       6 4 19
       2000000

       ------------------ cut here ----------------------


 <nobjectives>  This is the number of objectives.

 <gridpoints>   The analysis is done by plotting the Pareto Tradeoff
 surfaces found by each algorithm in each trial, and sampling the relative 
positions of the surfaces for each algorithm in a number of places throughout
the fitness space. The number of places, or samples, is controlled by this
parameter. Basically, lines are drawn, each starting from the origin, and
passing through the surfaces in a particualr region, Lines are spread out
equally throughout the space, covering all regions as well as possible,
The accuracy and conclusiveness of the results depends on the number of
lines, and the number of lines depends on this parameter as follows.
Where g is gridpoints, and d is dimensions (nobjectives), the number of
lines drawn is:

   dg^(d-1)

 So, in a 2D problem, making gridpoints=50 will give you 100 lines.
 In a 3D problem, making gridpoints = 10 will give you 300 lines, etc ...

 Recommened for good results are:

  D=2  g = 500
  D=3  g = 19
  D=4  g = 7

 These give around 1000 lines each time


 <output> This parameter shoudl be 0, 1, 2 or 3. Just set it to 0.
          The idea of it is to control whether or not points should
          be output for plotting worst, best, and median surfaces.
          It kind of works, but is not yet at the stage where I want
          to document it. So, keep this param at 0.

 <shift> This parameter shifts all the points by shift away from the origin in
         each dimension. This means that all the points will lie between shift
         and shift+1 in each dimension. Set shift to a large value e.g. 1000.0 to
         ensure that lines pass through the pareto front at roughly the same angle
         to the axes. This means that the front will be sampled more evenly when 
         the front is either highly convex or concave.



  IMPORTANT NOTE: mostate assumes that you are minimising each objective.
 Simply multiply the relevant objectives by -1 if this is not so.

  IMPORTANT NOTE: The fact that 1000000 and 2000000 are used as separators
  means that they cannot appear as values anywhere in your results! If so,
  things will be fouled up. So, just scale your results approprately to
  remove these if necessary.
  
WHAT COMES OUT


 Say your results file had results for N algorithms. The output (maybe
after several lines of other stuff first) will be a 2 by N matrix.
Entry i in the first row gives the percentage of the space on which
algoirthm i's performance is UNBEATEN by any of the other algorithms 
compared. Ie, the percentage of the fitness space for which we cannot
be 95% confident (based on a non-parametric test - The Mann-Whitney 
Rank test) that any other algoirthm beat it. Entry i in the second row gives
the percentage of the space on which we can be 95% confident that algorithm
i BEATS ALL of the other algorithms compared. Eg: you might get this for
three algoirthms:

 77.6 19.0 100.0
 0.0  0.0  11.3

 So, the third algorithm was not beaten anywhere in the space, and it
crushed the other two for sure on 11.3% of the space. The other two algortihms
had regions in which they did OK compared o 3, especially algorithm 1, but
they didn't do well enough to confidently beat 3, or each other, anywhere.

When 2 algorithms are compared the interpretation is obviously simpler,
and you naturally *think to see why) see that entries 1,1 and 2,2 and also
1,2 and 2,1 add to 100. Eg:

 16.2 88.8
 11.2 83.8

 Here, algorithm 2 wins again, although on mor ethan 10% of the space
algorithm 1 is clearly better.


 *
 * The program Uses Peter Ross' getint etc functions.
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#define TRUE 1
#define FALSE 0

#define MAXPOINTS 300000
#define MAXPOINTSETS 100

int print_median_points(int n);
int nextline();
int sort_intersects(int n);
int do_stats(int n, int a1, int a2);

double getdouble(FILE *,double *, int);
double ints;

extern double sqrt(double);
extern double pow(double, double);

int nobjs, nalgs;

int nithis;

int overall_algbeatsall[20];
int overall_algunbeaten[20];

/* lazy */

int n1, n2, n[20], algwins[20], alglosses[20];

typedef struct point {
 
 double o[20];
 int set;
 int ndom;
 int alg;
} POINT;

POINT points[MAXPOINTS];

typedef struct is {
 double v;
 int alg;
 double rank;
} IS;

IS iss[10000], tmpiss[1000];
int nis;


typedef struct pointset {

    int p; /* pointer to first point of set */
    int n; /* num ofpoints in this set */
    int alg; /* which algorithm */
    int set;
    POINT *points;
    double intersect;
    int btc; /* ADDED */
    int etc; /* ADDED */

 } POINTSET;

POINTSET pointsets[MAXPOINTSETS];

double max[10], min[10], range[10], v[10];


double getintersect(POINT *, int);

double inc, kk;

double shift;

int print_medians;

/* args:  pointsfile nobjs kk printmedians */
/*
         the data
                    nobjs! 
                          in 2d, 2*kk lines -- in Nd, Nkk^(n-1) lines ...
                              0 = do nothing; 1 = print best surface for each alg; 2 = median ; 3 = worst */

int main(int argc, char **argv)

{ 
 FILE *f1, *f2;
 int i, j, np1, np2, set, alg, n, nps;
 double ind, z, a, b, c, d;
 POINT *p;
 int gd, k;
 POINTSET *ps;
 IS *isp;

 alg = set = 0;

 f1 = fopen(argv[1],"r");

 i= 0;
 n = 0;
 nps = 0;


 nobjs = atoi(argv[2]);


 kk = (double)(atoi(argv[3]));

 print_medians= atoi(argv[4]);

 shift = atof(argv[5]);

 inc = 1.0/(kk+1.0);

 for (i=0;i<20;i++) overall_algbeatsall[i]=overall_algunbeaten[i]=0;

 
/* if(nobjs==2) {inc =  0.0196  0.001996008  ; 
 kk =  500  50 ;} 
 if(nobjs==3) {       inc = 0.1428571; kk = 6;  ( 
                  /* inc = 0.016949153 ; kk = 58; 
                  inc = 0.025 ; kk = 39;
                /*    inc = 0.1; kk = 9; 
                } */

 while (1)
  
   { gd = getdouble(f1,&ind,FALSE);
     if(gd==0) break;

     if(ind>=1000000) {
      
       /* do bookkeeping re pointsets */
       ps = &pointsets[set];
      
       ps->p = i-n; /* printf("%d\n",ps->p);*/
       ps->n = n; /* printf("%d\n",n); */
       ps->alg = alg;
       ps->set = set;
       ps->btc=0;  /* ADDED */
       ps->etc=0;  /* ADDED */
      nps++; n = 0; 
     set++; 

     if(ind==2000000) alg++; 
     continue;
     }

     p = &points[i];

     p->o[0]=ind;

    for(k=1;k<nobjs;k++)
     {getdouble(f1,&ind,FALSE);  p->o[k]=ind;}

     p->set=set;
     p->alg=alg;

    i++;   n++;
   }
 n1 = i;
 nalgs=alg;


 fclose(f1);

/*  printf("n1 = %d  \n",n1); */

 /* gather these into proper pointsets */


for(i=0;i<nps;i++)
  { 
    ps = &pointsets[i];

    ps->points = (POINT *)malloc(ps->n*sizeof(POINT));
    for(j=0;j<ps->n;j++)
      for(k=0;k<nobjs;k++)	 
	ps->points[j].o[k] = points[ps->p+j].o[k];

  }
 


 /*  for(j=0;j<nps;j++)
    { ps = &pointsets[j];
      printf("pointset %d alg %d\n",j, ps->alg);

       for(i=0;i<ps->n;i++)
         {for(k=0;k<nobjs;k++)
           printf(  " %g ", ps->points[i].o[k]);
         printf("\n");}
    }
*/

 /* Normalise the points */

  /* get range for each objective */

 ps = &pointsets[0];

 for(i=0;i<nobjs;i++)
 { min[i] = ps->points[0].o[i];
   max[i] = ps->points[0].o[i];
 }

 for(i=0;i<nps;i++)
   {ps = &pointsets[i];
    for(j=0;j<ps->n;j++)
       for(k=0;k<nobjs;k++)
         {
           ind = ps->points[j].o[k];
           if(ind<min[k]) min[k]=ind;
           if(ind>max[k]) max[k]=ind;
         }
  }   

for(k=0;k<nobjs;k++)
 range[k] = max[k]-min[k];



  ps = &pointsets[0];

  for(i=0;i<nps;i++)
    {ps = &pointsets[i];
        for(j=0;j<ps->n;j++)
          for(k=0;k<nobjs;k++)
	    {
         
            ind = ps->points[j].o[k];
            ps->points[j].o[k] = shift + ( ind - min[k])/range[k];
            } 
   }
 
/*  for(j=0;j<nps;j++)
   { ps = &pointsets[j];
           printf("pointset %d alg %d\n",j, ps->alg);

            for(i=0;i<ps->n;i++)
	      {for(k=0;k<nobjs;k++)
		            printf(  " %g ", ps->points[i].o[k]);
	                printf("\n");}
   }
 */

/* now, get the intersections of lines with the surfaces */

 while(nextline()==1)

   {
    printf("%g %g %g \n", v[0], v[1], v[2]);
      v[0] += shift; v[1] += shift; v[2] += shift;
     nis=0;
     
     for(j=0;j<nps;j++)
       {
         isp = &iss[nis];
         ps = &pointsets[j];
         ps->intersect = 10000000;  
           for(i=0;i<ps->n;i++)
	       for(k=0;k<nobjs;k++)
                 {ints = getintersect(&ps->points[i],k);   
                  if ((ints>-1) && (ints<ps->intersect))
                    ps->intersect = ints; 
		}
        /* printf("%g\n",ps->intersect);*/
        isp->v = ps->intersect;
        isp->alg = ps->alg;
        nis++;
       }


 /* BEGIN ADDED 

 /* if we are comparing a collection of surfaces with a single surface
  (assumed if there are just two algorithms and the second has only one 
 set of points) then we output all sorts of nice detail as follows */

 if ((nalgs==2) && (pointsets[nps-2].alg==0))

   { for(j=0;j<nps;j++)
      {
         ps = &pointsets[j];
         if(ps->intersect<pointsets[nps-1].intersect)
         ps->btc++;
         if(ps->intersect==pointsets[nps-1].intersect) 
         ps->etc++;
       }
   }

/* END ADDED */

 /* we now have the intersects for this line */
 sort_intersects(nis);

 for(i=0;i<nalgs;i++) algwins[i]=alglosses[i]=0;
 for(i=0;i<nalgs-1;i++)
  for(j=i+1;j<nalgs;j++)   
     do_stats(nis,i,j);

 for(i=0;i<nalgs;i++)
 { if(algwins[i]==(nalgs-1)) overall_algbeatsall[i]++;
   if(alglosses[i]==0) overall_algunbeaten[i]++;
 }
 
 /* for(i=0;i<nis;i++)
  printf(" %g %d %g \n", iss[i].v, iss[i].alg, iss[i].rank); */

 /* print medians if required */

 if(print_medians) print_median_points(nis);
 v[0] -= shift; v[1] -= shift; v[2] -= shift;


  }


 /* BEGIN ADDED */

 /* work out comparison details with single input surface if required */

 if ((nalgs==2) && (pointsets[nps-2].alg==0))
   {
     for(j=0;j<nps-1;j++)
      printf("SSC: %g%% better-or-equal  %g%% better\n", 
        100.0*(double)(pointsets[j].btc + pointsets[j].etc)/
       (pow(kk,(double)(nobjs-1))*(double)(nobjs)), 
        100.0*(double)(pointsets[j].btc)/
   (pow(kk,(double)(nobjs-1))*(double)(nobjs)));
   }
	     
/* END ADDED */


/* for(i=0;i<nps;i++)
   {
     j = 0; 

     /* start with the baseline */

/* BEGIN ADDED */
if((nalgs>2) || ((nalgs==2) && (pointsets[nps-2].alg==1)))
  {
/* END ADDED */
 for(i=0;i<nalgs;i++) printf(" %g ",
 (double)(100*overall_algunbeaten[i])/
(pow(kk,(double)(nobjs-1))*(double)(nobjs)));
printf("\n");

  for(i=0;i<nalgs;i++) printf(" %g ",(double)(100*overall_algbeatsall[i])/
(pow(kk,(double)(nobjs-1))*(double)(nobjs)));

 printf("\n");

/* BEGIN ADDED */      
 }
/* END EADDED */

}

double getintersect(POINT *p, int d)
{
 /* get intersect of this point with respect to dimension d
    with line v */

 double l;
 POINT in, *inp;
 int i, ok;

 inp = &in;

 l = p->o[d]/v[d];

 for(i=0;i<nobjs;i++)   inp->o[i] = l * v[i];  
 inp->o[d]=p->o[d];

 /* check bounds */

 ok=1;

 for(i=0;i<nobjs;i++)
/*   if((inp->o[i]>1) || (inp->o[i]<p->o[i])) ok=0; */
     if(inp->o[i]<p->o[i]) ok=0;

   

 if(ok==0) return(-1);
 else return(inp->o[0]);
}


/*
int checkdom(int a, int b)
{
 POINT *ap, *bp;
 int i, ones=0, negs=0, zeroes=0, j;
 int vec[100];

 ap = &points[a];
 bp = &points[b];

  for(i=0;i<nobjs;i++) 
    { if (ap->o[j]<bp->o[j]) {vec[j]=-1; negs++;}
      else if(ap->o[j]==bp->o[j]) {vec[j]=0; zeroes++;}
      else {vec[j]=1; ones++;}
    }

/* >0  negs  0 ones  a dominates b
   >0  ones  0 negs b  dominates a */

/*  if((negs>0) && (ones==0))
    bp->dom=1;
  else 
  if((ones>0) && (negs==0))
    ap->dom=1;
}*/

   

    



 

/***************************************************************/
/* Get the next number from the input: put it in the location  */
/* addressed by second argument. This function returns 0 on    */
/* EOF. If stopateol is true, it returns -1 when it hits \n    */
/* (after which some other procedure has to read past the \n), */
/* otherwise it continues looking for the next number.         */
/* A number has an optional sign, perhaps followed by digits,  */
/* perhaps followed by a decimal point, perhaps followed by    */
/* more digits. There must be a digit somewhere for it to count*/
/* as a number. So it would read any of:                       */
/*  -.5                                                        */
/*  -0.5                                                       */
/*  -.5.7                                                      */
/* as minus-a-half. In the last case, it would read .7 next    */
/* time around.                                                */
/*   There doesn't seem to be a neat and reliable way to do    */
/* all this, including stopateol, using scanf?                 */
/***************************************************************/

double 
getdouble (FILE *file, double *valaddr, int stopateol)
{
  int c;
  int found = FALSE, indecimal = FALSE;
  int sign = +1;
  double n = 0.0, p = 1.0;

  /* First find what looks like start of a number - the first digit. */
  /* And note any sign and whether we just passed a decimal point.   */  
  do {
    c = fgetc(file);
    if(c == EOF) return (0);
    else if(stopateol && c =='\n') return(-1);
    else if(c == '+' || c == '-') {
      sign = (c == '+')? +1 : -1;
      c = fgetc(file);
      if(c == EOF) return (0);
      else if(stopateol && c =='\n') return(-1);
    }
    if(c == '.') {
      indecimal = TRUE;
      c = fgetc(file);
      if(c == EOF) return (0);
      else if(stopateol && c =='\n') return(-1);
    }
    if(c >= '0' && c <= '9') {
      found = TRUE;
    } else {
      sign = +1;
      indecimal =  FALSE;
    }
  } while(!found);

  /* Now we've got digit(s) ... */
  do {
    n = 10.0*n + c - '0';
    p = 10.0*p;
    c = fgetc(file);

    if((c < '0') || (c > '9')) {
      found = FALSE;
      /* We've run out. If we already saw a decimal point, return now */
      if(indecimal) {
      	if(c != EOF) ungetc(c,file);
        *valaddr =  sign * n/p;
      	return(1);
      } else p = 1.0;
    }
  } while(found);

  /* We ran out and we didn't see a decimal point, so is this a decimal? */
  if(c != '.') {
    /* No, give it back to caller */
    if(c != EOF) ungetc(c,file);
    *valaddr = sign * n;
    return(1);
  } else {
    /* It is. Step past it, carry on hoping for more digits */
    c = fgetc(file);
    while(c >= '0' && c <= '9') {
      n = 10.0*n + c - '0';
      p = p*10.0;
      c =  fgetc(file);
    }
    /* We've run out of digits but we have a number to give */
    if(c != EOF) ungetc(c,file);
    *valaddr = sign * n/p;
    return(1);
  }
}

/* Use getdouble() above but convert result to int. */
int 
getint (FILE *f, int *valaddr, int stopateol)
{
  int r;
  double x;
  r = getdouble(f,&x,stopateol);
  *valaddr = (int)x;
  return(r);
}

/* Use getdouble above but convert result to long. */
int 
getlong (FILE *f, long *valaddr, int stopateol)
{
  int r;
  double x;
  r = getdouble(f,&x,stopateol);
  *valaddr = (long)x;
  return(r);
}


 

int nextline()
{

 int i, a, b, r;
 static int started = 0;

 static int d;

 r=1;

 if (started==0) 

   { started=1; d=0;

     for(i=0;i<nobjs;i++) v[i]=inc;

     v[d]=1;
   }
 else
   {

     /* increment the line */

     /* a will be the first nonmax dimension from the right */

    a = -1;

    for(i=nobjs-1;i>-1;i--)
     if(i!=d)
         if (v[i]<((kk-1)*inc + inc/2)) {a = i; break;}
      
 if(a==-1) r=0; 

 else {
  v[a]+=inc;
  
  for(i=a+1;i<nobjs;i++) v[i]=inc;
  v[d]=1;

  } }


 if(r==0)
   {
     /* we had all max value things */

 if(d<(nobjs-1))

   {d++; for(i=0;i<nobjs;i++) v[i]=inc; v[d]=1; r=1;}

}

return(r);

}
  
                 
   
     
     



   
  
int sort_intersects(int n)
{

 int i;
 IS *isp1, *isp2;
 int nochange=1;
 double tv;
 int ta, change;
 

 while(nochange)

   {
    change=0;
    for(i=0;i<n-1;i++)
    { isp1 = &iss[i];
      isp2 = &iss[i+1];
      
      if(isp1->v>isp2->v)
	{tv = isp1->v; ta = isp1->alg;
         isp1->v = isp2->v;
         isp1->alg = isp2->alg;
         isp2->v = tv;
         isp2->alg = ta;
         change=1;
       }
    }
    if(change==0) nochange=0;
  }
}
         
   
int do_stats(int n, int a1, int a2)
{

 int n1, n2, i, j, k, ntmpiss;
 double last, x, sum, u, t1, t2, u2, z, t;
 double dn1, dn2;
 IS *isp;



/* first, get into tmpiss the results concerning these algs only */

 ntmpiss=0;

  for(i=0;i<n;i++)

     if ((iss[i].alg==a1) || (iss[i].alg==a2))
         { isp = &tmpiss[ntmpiss++];
            isp->v = iss[i].v;
            isp->alg = iss[i].alg;
	 }


 last=-1;

 for(i=0;i<ntmpiss;i++)

   {
     isp = &tmpiss[i];
     isp->rank = (double)(i) + 1.0;
   }


for(i=0;i<ntmpiss;i++)
  {
    x = tmpiss[i].v;
    for(j=i+1;j<ntmpiss;j++)
      if(tmpiss[j].v!=tmpiss[i].v) break;

 
   /* started at i; j is first found which is different */

   /* now average the ranks between i and j-1 inclusive. */

  
    sum=0.0;
    for(k=i;k<j;k++) sum+=tmpiss[k].rank;
    for(k=i;k<j;k++)
    { isp = &tmpiss[k];
      isp->rank = sum/((double)(j-i));}
      
 /* now start again from j; */

    i = j-1;
  }

 n1 = n2 = 0;

t1 = t2  = 0.0; 
for(i=0;i<ntmpiss;i++)

 if(tmpiss[i].alg==a1) {n1++; t1 +=tmpiss[i].rank;}  
  else {n2++; t2+=tmpiss[i].rank;}


/* if(t1<t2) printf(" ALG %d BEATS ALG %d\n", a1, a2);
else  printf("ALG %d BEATS ALG %d\n", a2, a1); */

 dn1 = (double)(n1);
 dn2 = (double)(n2);

  u =  dn1*dn2 + (dn1*(dn1+1.0)/2.0)   - t1;

  u2 = (dn1*dn2-u);

if(u2<u) u=u2;

 z = ( u - (dn1*dn2/2.0) ) / ( sqrt(dn1*dn2*(dn1+dn2+1.0)/12.0));


  if( (t1<t2) && (z<-1.645)) {algwins[a1]++; alglosses[a2]++;}
  else if ((t2<t1) && (z<-1.645)) {algwins[a2]++; alglosses[a1]++;}
/*  printf(" Z = %g\n", z) ; */

}

  
int print_median_points(int n)
{

 int i, j, k;
 double med, lastv, sum;
 int num, num2;

 for(i=0;i<nalgs;i++)
   {
     /* find the median point for this alg. */

    /* first count the points with reasonable intersects. */

    sum = 0.0;
    num=0;
    for(j=0;j<n;j++)
     if(iss[j].alg==i) num++;
       
  
   /* now get the median or best or worst. */


  num2=0; 

 if(print_medians==2)
       for(j=0;j<n;j++)
         if(iss[j].alg==i)
               { num2++; 
                if(((double)(num2))>(((double)(num))/2.0)) break;
                else lastv = iss[j].v;
	       }
 if(print_medians==1)
       for(j=0;j<n;j++)
              if(iss[j].alg==i) break;
 if(print_medians==3)
    for(j=0;j<n;j++)
               if(iss[j].alg==i)
		 { num2++;
	              if(num2==num) break;
		 }
    
 

   if(print_medians==2) {
        if(num==(2*(num2-1))) med =  (iss[j].v +  lastv)/2.0;
        else med = iss[j].v;}

   else med = iss[j].v;


   
/*  printf(" num %d  num2 %d iss[j].v %g lastv %g med %g\n", num, num2, iss[j].v, lastv, med); */

  for(k=0;k<nobjs;k++) printf("%g ", 
    min[k] + (med*(v[k]/v[0])*range[k])) ; /* ADDED (CHANGED) FOR UNNORMALIZATION */
  /*  printf("%g %g ", med, (v[1]/v[0])*med); */
   }
  printf("\n");
}