/* this does permutation test for longitudinal data after Relly 2004, it determines the
smallest number of levels necessary to get the limiting value of the statistic

uses some C++ input/output but otherwise is all in C
*/

/* to use, need to specify (in the lines direclty above the start of main()):

1. nGLob=the total number of subjects,
2. tGlob=the number of observations per subject,
3. numIter=number of permutations to use in permutation test, et at 1000 currently
4. epsGlb=this is a number that is smaller than smallet difference between 2 different
       values of the reponse variable for a given subject-currenlty 1.0e-4

(typically 3 and 4 are fine as is)

then set up the files

1. longVal.dat (holds nGlob by tGlob matrix with data for each
     subject in a single row)
2. longTime.dat (a vector indicating where the times of the
     observations)

in the directory where this program will run-may need to set paths
for these too

*/

/* this assumes there are 2 groups with observations at same time points an no missing
data-this not essential for method, just the way this is programmed up */

#include <math.h>
#pragma hdrstop
#include <iostream.h>
#include <fstream.h>

#pragma argsused
#define FREE_ARG char*

static double dmaxarg1,dmaxarg2;
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
   (dmaxarg1) : (dmaxarg2))

#define SWAP(a,b) itemp=(a);(a)=(b);(b)=itemp;

#define dSWAP(a,b) temp=(a);(a)=(b);(b)=temp;

#define M 7
#define NSTACK 50
#define IAran1 16807
#define IMran1 2147483647
#define AMran1 (1.0/IMran1)
#define IQran1 127773
#define IRran1 2836
#define NDIVran1 (1+(IMran1-1)/32)
#define EPSran1 1.2e-7
#define RNMXran1 (1.0-EPSran1)

void long2surv(int nlevels,double *levels,double **longData,double **survUp,
   double **survDwn,double **upInd,double **dwnInd,double *maxLRup,
   double *maxLRdwn);
void getLevels(double **yGlob,int nGlob,int tGlob,int *nlevels,double *levels);
double logrank(double *times, double *indic);
void rpermute(int n,int *x,long *setSeed);
double *dvector(long nl,long nh);
void free_dvector(double *v, long nl, long nh);
int *ivector(long nl,long nh);
void free_ivector(int *v, long nl, long nh);
unsigned long *lvector(long nl,long nh);
void free_lvector(unsigned long *v, long nl, long nh);
double **dmatrix(long nrl,long nrh,long ncl,long nch);
void free_dmatrix(double **m, long nrl,long nrh,long ncl,long nch);
void sort2(unsigned long n, double *ra, double *rb);
void indexx(unsigned long n, double *arr,unsigned long *indx);
void sort(unsigned long n, double *arr);
double ran1(long *idum);
void nrerror(char error_text[]);

/* nGlob is total number of subjects, not number per group and tGlob is number of time
pts, use MISSINGVAL to specify a value for log rank to return when there are no events */

int nGlob=40,tGlob=5,numIter=1000;
double epsGlob=1.0e-4,MISSINGVAL=1.0e20,LARGEVAL=1.0e200,*timePt;

int main()
{
   /* this file holds data for both groups*/
   ifstream filer1("longVal.dat");
   if(!filer1){
      cerr << "Failure to open filer1\n";
      exit(EXIT_FAILURE);
   }
   /* this file holds time points of meaurements */
   ifstream filer2("longTime.dat");
   if(!filer2){
      cerr << "Failure to open filer2\n";
      exit(EXIT_FAILURE);
   }
   /* the p-value */
   ofstream filew1("longOut1.dat");
   if(!filew1){
      cerr << "Failure to open filew1\n";
      exit(EXIT_FAILURE);
   }
   /* numIter controls number of random permutations */
   int i,j,k,nlevels,*ints;
   long setSeed=-1876;
   double maxUp,maxDwn,obsStat,pval,maxstat,stat,upstat,dwnstat,minLevel=0.8,
      maxLevel=6.1,*vec1,*vec2,*vec3,*vec4,*vec5,*vec6,*vec7,*vec8,*uptimes,
      *upIndv,*InitLevels,*dwntimes,*dwnIndv,*times,*indic,*levels,**yGlob,**survUp,
      **survDwn,**upInd,**dwnInd;
   ints=ivector(1,nGlob);
   timePt=dvector(1,tGlob);
   times=dvector(1,nGlob);
   indic=dvector(1,nGlob);
   InitLevels=dvector(1,2*tGlob*nGlob+2*nGlob);
   vec1=dvector(1,nGlob/2);
   vec2=dvector(1,nGlob/2);
   vec3=dvector(1,nGlob/2);
   vec4=dvector(1,nGlob/2);
   vec5=dvector(1,nGlob/2);
   vec6=dvector(1,nGlob/2);
   vec7=dvector(1,nGlob/2);
   vec8=dvector(1,nGlob/2);
   uptimes=dvector(1,nGlob);
   upIndv=dvector(1,nGlob);
   dwntimes=dvector(1,nGlob);
   dwnIndv=dvector(1,nGlob);
   yGlob=dmatrix(1,nGlob,1,tGlob);
   /* read in data */
   filer1.seekg(0);
   for(i=1;i<=nGlob;i++){
      for(j=1;j<=tGlob;j++){
         filer1 >> yGlob[i][j];
      }
   }
   /* read in data */
   filer2.seekg(0);
   for(i=1;i<=tGlob;i++)
     filer2 >> timePt[i];
   /* set up levels */
   getLevels(yGlob,nGlob,tGlob,&nlevels,InitLevels);
   levels=dvector(1,nlevels);
   for(i=1;i<=nlevels;i++){
      levels[i]=InitLevels[i];
   }
   survUp=dmatrix(1,nGlob,1,nlevels);
   upInd=dmatrix(1,nGlob,1,nlevels);
   survDwn=dmatrix(1,nGlob,1,nlevels);
   dwnInd=dmatrix(1,nGlob,1,nlevels);
   /* intialize variable for random permutations */
   for(i=1;i<=nGlob;i++) ints[i]=i;
   /* convert to survival data sets and get minlog rank test stat */
   long2surv(nlevels,levels,yGlob,survUp,survDwn,upInd,dwnInd,&maxUp,&maxDwn);
   obsStat=DMAX(fabs(maxUp),fabs(maxDwn));
   cout << "obs stat=" << obsStat << endl;
   /* randomly permute the rows of the survival data set and
      recalculate the largest log rank statistic */
   pval=0.0;
   for(i=1;i<=numIter;i++){
      maxstat= -1.0e30;
      rpermute(nGlob,ints,&setSeed);
      for(j=1;j<=nlevels;j++){
         for(k=1;k<=nGlob/2;k++){
            vec1[k]=survUp[ints[k]][j];
            vec2[k]=upInd[ints[k]][j];
            vec3[k]=survDwn[ints[k]][j];
            vec4[k]=dwnInd[ints[k]][j];
         }
         for(k=nGlob/2+1;k<=nGlob;k++){
            vec5[k-nGlob/2]=survUp[ints[k]][j];
            vec6[k-nGlob/2]=upInd[ints[k]][j];
            vec7[k-nGlob/2]=survDwn[ints[k]][j];
            vec8[k-nGlob/2]=dwnInd[ints[k]][j];
         }
         /* data needs to be sorted within groups for log rank function-but
         should be a quicker way since I sort the data in that routine */
         sort2(nGlob/2,vec1,vec2);
         sort2(nGlob/2,vec3,vec4);
         sort2(nGlob/2,vec5,vec6);
         sort2(nGlob/2,vec7,vec8);
         for(k=1;k<=nGlob/2;k++){
           uptimes[k]=vec1[k];
           upIndv[k]=vec2[k];
           dwntimes[k]=vec3[k];
           dwnIndv[k]=vec4[k];
         }
         for(k=nGlob/2+1;k<=nGlob;k++){
           uptimes[k]=vec5[k-nGlob/2];
           upIndv[k]=vec6[k-nGlob/2];
           dwntimes[k]=vec7[k-nGlob/2];
           dwnIndv[k]=vec8[k-nGlob/2];
         }
         upstat=fabs(logrank(uptimes,upIndv));
         dwnstat=fabs(logrank(dwntimes,dwnIndv));
         if(upstat<MISSINGVAL && dwnstat<MISSINGVAL)
            stat=DMAX(upstat,dwnstat);
         if(upstat < dwnstat && dwnstat==MISSINGVAL) stat=upstat;
         if(upstat > dwnstat && upstat==MISSINGVAL) stat=dwnstat;
         if(stat > maxstat && stat<MISSINGVAL) maxstat=stat;
      }
      if(maxstat > obsStat) pval+=1.0;
   }
   pval/=numIter;
   cout << "pval=" << pval << endl;
   filew1 << pval << endl;
   free_ivector(ints,1,nGlob);
   free_dvector(timePt,1,tGlob);
   free_dvector(times,1,nGlob);
   free_dvector(indic,1,nGlob);
   free_dvector(InitLevels,1,2*tGlob*nGlob+2*nGlob);
   free_dvector(levels,1,nlevels);
   free_dvector(vec1,1,nGlob/2);
   free_dvector(vec2,1,nGlob/2);
   free_dvector(vec3,1,nGlob/2);
   free_dvector(vec4,1,nGlob/2);
   free_dvector(vec5,1,nGlob/2);
   free_dvector(vec6,1,nGlob/2);
   free_dvector(vec7,1,nGlob/2);
   free_dvector(vec8,1,nGlob/2);
   free_dvector(uptimes,1,nGlob);
   free_dvector(upIndv,1,nGlob);
   free_dvector(dwntimes,1,nGlob);
   free_dvector(dwnIndv,1,nGlob);
   free_dmatrix(yGlob,1,nGlob,1,tGlob);
   free_dmatrix(survUp,1,nGlob,1,nlevels);
   free_dmatrix(upInd,1,nGlob,1,nlevels);
   free_dmatrix(survDwn,1,nGlob,1,nlevels);
   free_dmatrix(dwnInd,1,nGlob,1,nlevels);
}

void getLevels(double **yGlob,int nGlob,int tGlob,int *nlevels,double *levels)
{
/* levels hould be a vector of length 2*tGlob*nGlob+2*nGlob which should
   be trimmed down in main-the last elements of this vector set to the global
   variable value of LARGEVAL */
   int i,j,k,k1;
/* eps should be smaller than the minimum difference between any 2 values of
the reponse variable-so be careful */
   double eps=epsGlob,miny,maxy,locM,*cc,*lowy,*hiy,*tempVec;
   cc=dvector(1,2*tGlob*nGlob);
   lowy=dvector(1,nGlob);
   hiy=dvector(1,nGlob);
   /* first do upcrosing levels, but do miny and maxy here */
   k1=1;
   for(i=1;i<=nGlob;i++){
      miny=maxy=yGlob[i][1];
      for(j=1;j<=tGlob-1;j++){
         if(yGlob[i][j+1]>maxy) maxy=yGlob[i][j+1];
         if(yGlob[i][j+1]<miny) miny=yGlob[i][j+1];
         if(yGlob[i][j+1]>yGlob[i][j]){
            locM=yGlob[i][j];
            for(k=1;k<=j;k++) if(yGlob[i][k]>locM) locM=yGlob[i][k];
            if(yGlob[i][j+1]>locM){
               cc[k1]=locM;
               k1+=1;
            }
         }
      }
      cc[k1]=maxy;
      k1+=1;
      cc[k1]=miny;
      k1+=1;
      lowy[i]=miny;
      hiy[i]=maxy;
   }
   /* then do downcrosing levels */
   for(i=1;i<=nGlob;i++){
      for(j=1;j<=tGlob-1;j++){
         if(yGlob[i][j+1]<yGlob[i][j]){
            locM=yGlob[i][j];
            for(k=1;k<=j;k++) if(yGlob[i][k]<locM) locM=yGlob[i][k];
            if(yGlob[i][j+1]<locM){
               cc[k1]=locM;
               k1+=1;
            }
         }
      }
   }
   *nlevels=k1-1+2*nGlob;
   tempVec=dvector(1,k1-1+2*nGlob);
   for(i=1;i<=k1-1;i++){
      tempVec[i]=cc[i];
   }
   for(i=1;i<=nGlob;i++) tempVec[k1-1+i]=lowy[i];
   for(i=1;i<=nGlob;i++) tempVec[k1-1+nGlob+i]=hiy[i];
   sort(k1-1+2*nGlob,tempVec);
   for(i=1;i<=k1-1;i++) levels[i]=cc[i]+eps;
   for(i=1;i<=nGlob;i++) levels[k1-1+i]=lowy[i]-eps;
   for(i=1;i<=nGlob;i++) levels[k1-1+nGlob+i]=hiy[i]-eps;
   for(i=k1+2*nGlob;i<=2*tGlob*nGlob+2*nGlob;i++) levels[i]=LARGEVAL;
   sort(2*tGlob*nGlob+2*nGlob,levels);
   free_dvector(tempVec,1,k1-1+2*nGlob);
   free_dvector(cc,1,2*tGlob*nGlob);
   free_dvector(lowy,1,nGlob);
   free_dvector(hiy,1,nGlob);
}

void long2surv(int nlevels,double *levels,double **longData,double **survUp,
   double **survDwn,double **upInd,double **dwnInd,double *maxLRup,
   double *maxLRdwn)
/* this takes an nGlob by tGlob matrix of longitudinal data
and converts it to 2 survival data sets-these are matrices of times
of first upcrossing and times of first downcrossing and 2 matrices
that indicate if uncensored-also gets log rank stat along the way */
{
   int i,j,k;
   double *vec1,*vec2,*vec3,*vec4,*vec5,*vec6,*vec7,*vec8,*survUpV,*survDwnV,
      *upIndV,*dwnIndV,*tup,*tdwn;
   vec1=dvector(1,nGlob/2);
   vec2=dvector(1,nGlob/2);
   vec3=dvector(1,nGlob/2);
   vec4=dvector(1,nGlob/2);
   vec5=dvector(1,nGlob/2);
   vec6=dvector(1,nGlob/2);
   vec7=dvector(1,nGlob/2);
   vec8=dvector(1,nGlob/2);
   survUpV=dvector(1,nGlob);
   survDwnV=dvector(1,nGlob);
   upIndV=dvector(1,nGlob);
   dwnIndV=dvector(1,nGlob);
   tup=dvector(1,nlevels);
   tdwn=dvector(1,nlevels);
   for(i=1;i<=nlevels;i++){
      for(j=1;j<=nGlob;j++){
         /* next line just a trick to easily test if subject ever reaches the level */
         survUp[j][i]= -1.0;
         survDwn[j][i]= -1.0;
         /* see if subject censored drom start, i.e. already above/below the level */
         if(longData[j][1]>levels[i]){
            survUpV[j]=survUp[j][i]=timePt[1];
            upIndV[j]=upInd[j][i]=0;
         }
         else{
            for(k=2;k<=tGlob;k++){
               if(longData[j][k]>levels[i]){
                  survUpV[j]=survUp[j][i]=timePt[k];
                  upIndV[j]=upInd[j][i]=1;
                  break;
               }
            }
         }
         if(longData[j][1]<levels[i]){
            survDwnV[j]=survDwn[j][i]=timePt[1];
            dwnIndV[j]=dwnInd[j][i]=0;
         }
         /* otherwise find first timePt that it occurs */
         else{
            for(k=2;k<=tGlob;k++){
               if(longData[j][k]<levels[i]){
                  survDwnV[j]=survDwn[j][i]=timePt[k];
                  dwnIndV[j]=dwnInd[j][i]=1;
                  break;
               }
            }
         }
         /* and if never reach the level they are censored at the end */
         if(survUp[j][i] < 0.0){
            survUpV[j]=survUp[j][i]=timePt[tGlob];
            upIndV[j]=upInd[j][i]=0;
         }
         if(survDwn[j][i] < 0.0){
            survDwnV[j]=survDwn[j][i]=timePt[tGlob];
            dwnIndV[j]=dwnInd[j][i]=0;
         }
      }
      /* now find log rank statistic at this level-need to sort within groups first */
      for(j=1;j<=nGlob/2;j++){
         vec1[j]=survUpV[j];
         vec2[j]=upIndV[j];
         vec3[j]=survDwnV[j];
         vec4[j]=dwnIndV[j];
      }
      for(j=nGlob/2+1;j<=nGlob;j++){
         vec5[j-nGlob/2]=survUpV[j];
         vec6[j-nGlob/2]=upIndV[j];
         vec7[j-nGlob/2]=survDwnV[j];
         vec8[j-nGlob/2]=dwnIndV[j];
      }
      /* data needs to be sorted within groups for log rank function-but
         should be a quicker way since I sort the data in that routine */
      sort2(nGlob/2,vec1,vec2);
      sort2(nGlob/2,vec3,vec4);
      sort2(nGlob/2,vec5,vec6);
      sort2(nGlob/2,vec7,vec8);
      for(j=1;j<=nGlob/2;j++){
         survUpV[j]=vec1[j];
         upIndV[j]=vec2[j];
         survDwnV[j]=vec3[j];
         dwnIndV[j]=vec4[j];
      }
      for(j=nGlob/2+1;j<=nGlob;j++){
         survUpV[j]=vec5[j-nGlob/2];
         upIndV[j]=vec6[j-nGlob/2];
         survDwnV[j]=vec7[j-nGlob/2];
         dwnIndV[j]=vec8[j-nGlob/2];
      }
      tdwn[i]=logrank(survDwnV,dwnIndV);
      tup[i]=logrank(survUpV,upIndV);
   }
   *maxLRup=*maxLRdwn=0.0;
   for(i=1;i<=nlevels;i++)
      if(fabs(*maxLRup) < fabs(tup[i]) && tup[i]<MISSINGVAL) *maxLRup=fabs(tup[i]);
   for(i=1;i<=nlevels;i++)
      if(fabs(*maxLRdwn) < fabs(tdwn[i]) && tdwn[i]<MISSINGVAL) *maxLRdwn=fabs(tdwn[i]);
   free_dvector(vec1,1,nGlob/2);
   free_dvector(vec2,1,nGlob/2);
   free_dvector(vec3,1,nGlob/2);
   free_dvector(vec4,1,nGlob/2);
   free_dvector(vec5,1,nGlob/2);
   free_dvector(vec6,1,nGlob/2);
   free_dvector(vec7,1,nGlob/2);
   free_dvector(vec8,1,nGlob/2);
   free_dvector(survUpV,1,nGlob);
   free_dvector(survDwnV,1,nGlob);
   free_dvector(upIndV,1,nGlob);
   free_dvector(dwnIndV,1,nGlob);
   free_dvector(tup,1,nlevels);
   free_dvector(tdwn,1,nlevels);
}

double logrank(double *times, double *indic)
{
/* times should be sorted within each of the 2 groups */
/* this can return answers like 1.0e-300 when e1=dd[1], so test stat is zero */
/* returns the global variable MISSINGVAL if  < 1 event */
   int i,j,noTimes,nDie=0;
   double yy,dd,e1,v1,qnum=0.0,qden=0.0,*temp,*loctimes,*timO,*yy1,*yy2,*dd1,*dd2,
      yy1init,yy2init;
   temp=dvector(1,nGlob);
/* first order and pick out the distinct failures times */
   for(i=1;i<=nGlob;i++){
      nDie+=indic[i];
      if(indic[i]==1) temp[nDie]=times[i];
   }
   if(nDie < 1){
      /* use large value to indicate that no events took place */
      qnum=MISSINGVAL;
      qden=1.0;
   }
   else{
      loctimes=dvector(1,nDie);
      timO=dvector(1,nDie);
      for(i=1;i<=nDie;i++) loctimes[i]=temp[i];
      sort(nDie,loctimes);
      timO[1]=loctimes[1];
      noTimes=1;
      for(i=2;i<=nDie;i++){
         if(loctimes[i]!=timO[noTimes]){
            noTimes+=1;
            timO[noTimes]=loctimes[i];
         }
      }
      yy1=dvector(1,noTimes+1);
      yy2=dvector(1,noTimes+1);
      dd1=dvector(1,noTimes);
      dd2=dvector(1,noTimes);
      for(i=1;i<=noTimes;i++)
         dd1[i]=dd2[i]=0.0;
      /* now determine the number failures and no. at risk in each group at
      each of these distinct times */
      yy1[1]=nGlob/2;
      yy1init=yy1[1];
      for(i=1;i<=noTimes;i++){
         yy1[i]=yy1init;
         for(j=1;j<=nGlob/2;j++){
            if(times[j]>timO[i]) break;
            if((times[j] == timO[i]) && (indic[j]==1)) dd1[i]+=1.0;
            if(times[j] < timO[i]) yy1[i]-=1.0;
         }
      }
      yy2[1]=nGlob/2;
      yy2init=yy2[1];
      for(i=1;i<=noTimes;i++){
         yy2[i]=yy2init;
         for(j=nGlob/2+1;j<=nGlob;j++){
            if(times[j]>timO[i]) break;
            if((times[j] == timO[i]) && (indic[j]==1)) dd2[i]+=1.0;
            if(times[j] < timO[i]) yy2[i]-=1.0;
         }
      }
      for(i=1;i<=noTimes;i++){
         yy=yy1[i]+yy2[i];
         if(yy > 1.0){
            dd=dd1[i]+dd2[i];
            e1=dd*yy1[i]/yy;
            v1=e1*yy2[i]*(yy-dd)/(yy*(yy-1.0));
            qnum+=(dd1[i]-e1);
            qden+=v1;
         }
      }
      free_dvector(loctimes,1,nDie);
      free_dvector(timO,1,nDie);
      free_dvector(yy1,1,noTimes+1);
      free_dvector(yy2,1,noTimes+1);
      free_dvector(dd1,1,noTimes);
      free_dvector(dd2,1,noTimes);
   }
   free_dvector(temp,1,nGlob);
   if(qden==0.0 && qnum==0.0) return MISSINGVAL;
   else return qnum/sqrt(qden);
}

void rpermute(int n,int *x,long *setSeed)
/* this takes a vector and returns a random permutation of it using
the algorithm from Knuth Art of Scient. Comp. vol. 2-there it is not
clear what m should be, so just set it equal to n */
{
   int i,j,temp,m=n;
   for(i=m;i>1;i--){
      j=int(i*ran1(setSeed))+1;
      temp=x[i];
      x[i]=x[j];
      x[j]=temp;
   }
}

/* these routines from Numerical Recipes in C by Press, Teukolsky, Vetterling, and Flannery
1992 Cambridge University Press */

double *dvector(long nl,long nh)
/* allocate a vector with index range v[nl...nh] */
{
   void nrerror(char error_text[]);
   double *v;
   v=(double *)malloc((size_t) ((nh-nl+2)*sizeof(double)));
   if (!v) nrerror("allocation failure in dvector()");
   return v-nl+1;
}

void free_dvector(double *v, long nl, long nh)
{
   int n;
   free((FREE_ARG) (v+nl-1));
}

int *ivector(long nl,long nh)
/* allocate a vector with index range v[nl...nh] */
{
   void nrerror(char error_text[]);
   int *v;
   v=(int *)malloc((size_t) ((nh-nl+2)*sizeof(int)));
   if (!v) nrerror("allocation failure in dvector()");
   return v-nl+1;
}

void free_ivector(int *v, long nl, long nh)
{
   int n;
   free((FREE_ARG) (v+nl-1));
}

unsigned long *lvector(long nl,long nh)
/* allocate a vector with index range v[nl...nh] */
{
   void nrerror(char error_text[]);
   unsigned long *v;
   v=(unsigned long *)malloc((size_t) ((nh-nl+2)*sizeof(long)));
   if (!v) nrerror("allocation failure in dvector()");
   return v-nl+1;
}

void free_lvector(unsigned long *v, long nl, long nh)
{
   int n;
   free((FREE_ARG) (v+nl-1));
}


double **dmatrix(long nrl,long nrh,long ncl,long nch)
/* allocate a double matrix with subscript range m[nrl...nrh][ncl...nch] */
{
   long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
   double **m;
   void nrerror(char error_text[]);
   /* allocate pointers to rows */
   m=(double **) malloc((size_t)((nrow+1)*sizeof(double*)));
   if (!m) nrerror("allocation failure 1 in dmatrix()");
   m += 1;
   m -= nrl;

   /* allocate rows and set pointers to them */
   m[nrl]=(double *) malloc((size_t)((nrow*ncol+1)*sizeof(double)));
   if (!m[nrl]) nrerror("allocation failure 2 in dmatrix()");
   m[nrl] += 1;
   m[nrl] -= ncl;

   for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

   /* return pointer to array of pointers to rows */
   return m;
}

void free_dmatrix(double **m, long nrl,long nrh,long ncl,long nch)
/* free a double matrix allocated by dmatrix() */
{
   free((FREE_ARG) (m[nrl]+ncl-1));
   free((FREE_ARG) (m+nrl-1));
}

void sort2(unsigned long n, double *ra, double *rb)
/* this sorts ra and rearranges rb so that its order matches ra */
{
   void indexx(unsigned long n, double *arr,unsigned long *indx);
   unsigned long j,*iwksp;
   double *wksp;
   iwksp=lvector(1,n);
   wksp=dvector(1,n);
   indexx(n,ra,iwksp);
   for(j=1;j<=n;j++) wksp[j]=ra[j];
   for(j=1;j<=n;j++) ra[j]=wksp[iwksp[j]];
   for(j=1;j<=n;j++) wksp[j]=rb[j];
   for(j=1;j<=n;j++) rb[j]=wksp[iwksp[j]];
   free_dvector(wksp,1,n);
   free_lvector(iwksp,1,n);
}


void indexx(unsigned long n, double *arr,unsigned long *indx)
/* this indexes an array */
{
   unsigned long i,indxt,ir=n,itemp,j,k,l=1;
   int jstack=0,*istack;
   double a;
   istack=ivector(1,NSTACK);
   for(j=1;j<=n;j++) indx[j]=j;
   for(;;){
      if(ir-l < M){
         for(j=l+1;j<=ir;j++){
            indxt=indx[j];
            a=arr[indxt];
            for(i=j-1;i>=l;i--){
               if(arr[indx[i]] <= a) break;
               indx[i+1]=indx[i];
            }
            indx[i+1]=indxt;
         }
         if(jstack==0) break;
         ir=istack[jstack--];
         l=istack[jstack--];
      }
      else{
         k=(l+ir) >> 1;
         SWAP(indx[k],indx[l+1]);
         if(arr[indx[l]] > arr[indx[ir]]){
            SWAP(indx[l],indx[ir])
         }
         if(arr[indx[l+1]]>arr[indx[ir]]){
            SWAP(indx[l+1],indx[ir])
         }
         if(arr[indx[l]]>arr[indx[l+1]]){
            SWAP(indx[l],indx[l+1])
         }
         i=l+1;
         j=ir;
         indxt=indx[l+1];
         a=arr[indxt];
         for(;;){
            do i++; while (arr[indx[i]]<a);
            do j--; while (arr[indx[j]]>a);
            if(j<i) break;
            SWAP(indx[i],indx[j])
         }
         indx[l+1]=indx[j];
         indx[j]=indxt;
         jstack += 2;
         if(jstack > NSTACK) nrerror("NSTACK too small in indexx.");
         if(ir-i+1 >= j-l){
            istack[jstack]=ir;
            istack[jstack-1]=i;
            ir=j-1;
         }
         else{
            istack[jstack]=j-1;
            istack[jstack-1]=l;
            l=i;
         }
      }
   }
   free_ivector(istack,1,NSTACK);
}

void sort(unsigned long n, double *arr)
/* does quicksort */
{
   unsigned long i,ir=n,j,k,l=1;
   int jstack=0,*istack;
   double a, temp;
   istack=ivector(1,NSTACK);
   for(;;){
      if(ir-l<M){
         for(j=l+1;j<=ir;j++){
            a=arr[j];
            for(i=j-1;i>=l;i--){
               if(arr[i] <= a) break;
               arr[i+1]=arr[i];
            }
            arr[i+1]=a;
         }
         if(jstack==0) break;
         ir=istack[jstack--];
         l=istack[jstack--];
      }
      else{
         k=(l+ir) >> 1;
         dSWAP(arr[k],arr[l+1])
         if(arr[l] > arr[ir]){
            dSWAP(arr[l],arr[ir])
         }
         if(arr[l+1] > arr[ir]){
            dSWAP(arr[l+1],arr[ir])
         }
         if(arr[l] > arr[l+1]){
            dSWAP(arr[l],arr[l+1])
         }
         i=l+1;
         j=ir;
         a=arr[l+1];
         for(;;){
            do i++; while(arr[i] < a);
            do j--; while(arr[j] > a);
            if(j < i) break;
            dSWAP(arr[i],arr[j])
         }
         arr[l+1]=arr[j];
         arr[j]=a;
         jstack+=2;
         if(jstack>NSTACK) nrerror("NSTACK too small in sort");
         if(ir-i+1>= j-l){
            istack[jstack]=ir;
            istack[jstack-1]=i;
            ir=j-1;
         }
         else{
            istack[jstack]=j-1;
            istack[jstack-1]=l;
            l=i;
         }
      }
   }
   free_ivector(istack,1,NSTACK);
}

double ran1(long *idum)
{
/* generates U(0,1) deviate */
   const int NTABran1=32;
   int j;
   long k;
   static long iy=0;
   static long iv[NTABran1];
   double temp;

   if (*idum <= 0 || !iy) {
      if (-(*idum)<1) *idum=1;
      else *idum = -(*idum);
      for(j=NTABran1+7;j>=0;j--) {
         k=(*idum)/IQran1;
         *idum=IAran1*(*idum-k*IQran1)-IRran1*k;
         if (*idum<0) *idum += IMran1;
         if (j<NTABran1) iv[j] = *idum;
      }
      iy=iv[0];
   }
   k=(*idum)/IQran1;
   *idum=IAran1*(*idum-k*IQran1)-IRran1*k;
   if (*idum<0) *idum += IMran1;
   j=iy/NDIVran1;
   iy=iv[j];
   iv[j] = *idum;
   if ((temp=AMran1*iy)>RNMXran1) return RNMXran1;
   else return temp;
}

void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
   fprintf(stderr,"Run-time error...\n");
   fprintf(stderr,"%s\n",error_text);
   fprintf(stderr,"...now exiting to system...\n");
   exit(1);
}