/* One-Way layouts:
   If treatment groups are significantly different by the Kruskal-Wallis
     test, which PAIRS of treatments are significantly different?
     If the first treatment group is a control, which INDIVIDUAL
     treatments are significantly better than the control?
   Recall that if you do e.g. 50 tests, then even if there are no 
     differences, approximately 0.05*50 = 2 1/2 tests will show up
     as false positives by chance. This is called the
     MULTIPLE COMPARISON problem. The goal is to make inferences
     that are not artifacts of multiple comparisons. */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* Prototypes for all functions in  statlib.c */
#include "statlib.h"

#define MAXK 10       /* Maximum possible number of treatments */
#define MAXN 30       /* Maximum number of observations per treatment */
#define MAXTOT (MAXK*MAXN)

/* Description of data set and treatment groups: */

char *title="Numbers of glucocorticoid receptor sites per Leukocyte cell";
char *source="Hollander & Wolfe (Table 6.4, p201)";

/* Number of treatment groups and the numbers of observations */
/*   in each treatment group */

const int kk=5;
const int nn[MAXK] = { 14, 5, 6, 4, 8 };

/* A double array for k treatment groups: */
/* See OneWayLayout.c for a discussion of double arrays in C */
/*   and how to initialize them. */
/* Here  xx[i][j]  means the j-th element of the  */
/*   i-th row or treatment group.              */

/* Treatment-group names */

const char *treatnames[MAXK] =
 { "Normal", "HCell Anemia", "CLL", "CML", "Acute" };

/* Treatment groups are rows: */

const double xx[MAXK][MAXN] = {
 { 3500, 3500, 3500, 4000, 4000, 4000, 4300, 4500, 4500, 4900,
      5200, 6000, 6750, 8000 },
 { 5710, 6110, 8060, 8080, 11400 },
 { 2930, 3330, 3580, 3880, 4280, 5120 },
 { 6320, 6860, 11400, 14000 },
 { 3230, 3880, 7640, 7890, 8280, 16200, 18250, 29900 } };

/* The starting random-numbr seed. */

unsigned long wseed=123456;

void do_multcomp(const double xx[MAXK][MAXN], int kk, const int nn[MAXK],
     const char *treatnames[MAXK]);


int main(int argc, char *argv[])
 { 
   printf ("Multiple Comparison procedures for One-Way layouts:\n");
   /* Helpful messages on segmentation faults and some other fatal errors */
   /* This function is in statlib.c */
   set_newsignals();

   printf ("Kruskal-Wallis and multiple comparison tests.\n");

   /* If a number is entered on the command line, use it as the */
   /*   starting seed instead of wseed. If the starting seed is 0, */
   /*   initialize the seed from the system clock. */
   /* `atol' stands for `ASCII (text) to long (integer)' */
   if (argc>=2) wseed=atol(argv[1]);
   wseed=setrand(wseed);  /* If wseed=0, replace by true starting seed */
   printf ("\nStarting random-number seed: %lu\n",wseed);
   printf ("  (Enter a number on the command line "
     "to specify a different seed.)\n");

   /* See above for dataset title and source */
   printf ("\nData: %s\n%s\n", title, source);
   do_multcomp(xx,kk,nn,treatnames);
   show_randnums_used();

   return 0; }  /* Return with no errors */



double get_hobs (double *mrank, int *nfirst, int *nmax, int kk);
double get_wilcoxon(const double *wboth,
   int nn1, int nn2, double *wstatptr);


void do_multcomp(const double xx[MAXK][MAXN], int kk, const int nn[MAXK],
     const char *treatnames[MAXK])
 { int i,j,r, ns,nbig,nsims, ntot,ntg, ncon;
   int nfirst[MAXK], nmax[MAXK], tgsize[MAXTOT], nbigg[MAXK][MAXK];
   double hsum, hh, hobs, tiesum,tiecorr, maxzstar, pp,pwid;
   double wwdat[MAXTOT], mrank[MAXTOT], rankav[MAXK];
   double wmcbuff[2*MAXN], wdobs[MAXK][MAXK];
   double zdobs[MAXK][MAXK], wzdobs[MAXK][MAXK];

   /* Set `ntot' to be the total number of observations */
   ntot=0;
   for (i=0; i<kk; i++)  ntot+=nn[i];
   printf ("\nNumber of treatments: %d    "
     "Total number of observations: %d\n", kk,ntot);
   /* Display the data */
   for (i=0; i<kk; i++)
    { printf ("Group #%d (n=%d): ", i+1, nn[i]);
      for (j=0; j<nn[i]; j++) 
       { if (j>0 && (j%10)==0) printf("\n");
         printf (" %5g", xx[i][j]); }
      printf("\n"); }

   /* Find the Kruskal-Wallis midranks */
   /* See OneWayLayout.c for comments about method */
   for (r=i=0; i<kk; i++)       /* Sum over treatment groups */
    { nfirst[i]=r;              /* First offset for i-th group */
      for (j=0; j<nn[i]; j++)   /* Sum over i-th treatment group */
        wwdat[r++]=xx[i][j];    /* Copy values to wwdat[] */
      nmax[i]=r; }              /* Last-plus-one offset for i-th group */

   /* Put midranks of the array  wwdat[]  in  mrank[] */
   /* This function is in statlib.c */
   makeranks(mrank, wwdat,ntot, tgsize,&ntg);

   /* Display data with midranks */
   printf ("\nData with midranks (RAv=RankSum/n_i):\n");
   /* Sum over treatment groups */
   for (r=i=0; i<kk; i++)
    { int nf=nfirst[i], len=nn[i];  double rsum;
      /* Rank average of i-th treatment group */
      for (rsum=0.0, j=0; j<len; j++) rsum += mrank[nf+j];
      rankav[i]=rsum/len;
      printf (" #%d (RAv=%5.2f, n=%d, %s):",
        i+1, rankav[i], nn[i], treatnames[i]);
      /* Display the row */
      for (j=0; j<len; j++)
       { if ( (j%5)==0)  printf("\n     ");
         printf ("  %4g(%4.1f)", wwdat[nf+j], mrank[nf+j]); }
      printf("\n"); }

   /* Find the large-sample P-value of the Kruskal-Wallis test */
   printf ("\nThe large-sample Kruskal-Wallis test:\n");
   hsum=0.0;
   for (i=0; i<kk; i++)
    { double rdif=rankav[i]-(double)(ntot+1)/2;
      hsum+=nn[i]*rdif*rdif; }
   hh=(12.0/(ntot*(ntot+1)))*hsum;
   /* Use  (tgsize[],ntg)  to find the Kruskal-Wallis test Tiesum */
   /*   and the tie-corrected statistic  HstatPrime  */
   tiesum=0.0;
   for (i=0; i<ntg; i++)
    { int ts=tgsize[i];  tiesum += (ts-1)*ts*(ts+1); }
   tiecorr = tiesum/((double)(ntot-1)*ntot*(ntot+1));
   hh /= (1 - tiecorr);

   printf ("  HP=%g   Tiegroups=%d   Tiesum=%g   Tiecorr=%g\n", 
     hh,ntg,tiesum,tiecorr);
   printf ("  P=%6.4f (Chi-square approx, df=%d)\n",
     pvchisq(kk-1,hh), kk-1);

   /* Simulate the true Kruskal-Wallis P-value: */
   /* Try with nsims=1,000,000 */
   nsims=1000000;
   /* commnum(nn)  writes   */
   /*    nn:  100  1000  10000   100000   1000000    */
   /*    as:  100  1000  10,000  100,000  1,000,000  */
   /* Warning: commnum() writes to a char buffer, so replace %d by %s */
   /* commnum() is in statlib.c */
   printf ("\nSimulating the true KW P-value with "
     "n=%s permutations:\n",commnum(nsims));
   hobs=get_hobs(mrank,nfirst,nmax,kk);
   printf ("  Hobs = Sum RankSum_i^2/n_i = %g\n", hobs);
   /* Carry out permutation test with `nsims' permutations of mrank[] */
   nbig=0;
   for (ns=0; ns<nsims; ns++)
    { double hscore;
      shuffle_dub(mrank,ntot);  /* Shuffle `mrank' array */
      hscore=get_hobs(mrank,nfirst,nmax,kk);
      /* Allow for small roundoff error */
      if (hscore>=hobs-1e-8) nbig++; }
   pp=(double)nbig/nsims;
   pwid = 1.960*sqrt(pp*(1-pp)/nsims);
   printf ("  P = %d/%s = %7.5f   95%% CI: (%7.5f, %7.5f, %7.5f)\n",
     nbig,commnum(nsims), pp, pp-pwid, pp, pp+pwid);

   /* Two-way multiple comparisons:
      Given that the Kruskal-Wallis test is significant,
        which PAIRS OF TREATMENTS are significantly different,
        allowing for all possible multiple comparisons? */
   printf (
   "\nWhich of the %d (= k*(k-1)/2) pairwise comparisons between the\n"
     "  %d (=k) treatment groups are significant, allowing for\n"
     "  multiple comparisons?\n",  kk*(kk-1)/2, kk);
   printf (
     "Using the Steel-Dwass-Crichlow-Fligner multiple-comparison\n"
     "  procedure described in Section 6.5 of Hollander & Wolfe.\n"
     "That is, let ObsZ(i,j) be the pairwise Wilcoxon rank-sum score\n"
     "  between all treatment-group pairs, normalized to have mean zero\n"
     "  and variance one. For a permutation of the data, let mstar =\n"
     "  max_{a,b}Z(a,b) where Z(i,j) is the corresponding normalized\n"
     "  Wilcoxon score between pairs of permuted data. Finally, let\n"
     "  P(i,j) be the proportion of permutations for which\n"
     "  mstar >= |ObsZ(i,j)|.\n"
     "Then P(i,j) are multiple-comparison-corrected P-values for all\n"
     "  treatment pairs (i,j).\n");
   /* First, find the observed pairwise Wilcoxon rank-sum statistics
        (Wij=wdobs[i][j]) and the corresponding large-sample tie-corrected
        Z-statistics  (zdobs[i][j])
      Loop over pairs of treatments: */
   maxzstar=0.0;  /* For display */
   for (i=0; i<kk; i++)  for (j=i+1; j<kk; j++)
    { int s, ni=nn[i], nj=nn[j];  double zstar;
      /* For use by get_wilcoxon(), copy the i-th treatment group to
           wmcbuff[]  followed by the j-th treatment group */
      for (r=s=0; s<ni; s++)  wmcbuff[r++]=xx[i][s];
      for (s=0; s<nj; s++)  wmcbuff[r++]=xx[j][s];
      /* Compute the Wilcoxon ranksum statistic for groups #i and #j
           based on group #j (call it W) and the large-sample tie-corrected
           normal approximation Zstar.  The function sets  wd[i][j]=W
           and returns the normal score Zstar. */
      zstar=get_wilcoxon(wmcbuff, ni,nj, &(wdobs[i][j]));
      /* Save the absolute value and keep track of the maximum |Zstar| */
      zdobs[i][j]=zstar=fabs(zstar);
      if (zstar>maxzstar) maxzstar=zstar; }
   printf ("\nFor comparison, the following are "
      "multiple-comparison-UNCORRECTED P-values.\n");
   printf ("Pairwise Wilcoxon W and tie-corrected Z values:\n");
   for (i=0; i<kk; i++)  for (j=i+1; j<kk; j++)
     printf ("  #%d vs. #%d:  W=%5.1f  E(W)=%5.1f  "
       "Z=%.3f   P=%.4f (2-sided, UNCORR)\n",
       i+1,j+1, wdobs[i][j], (double)nn[j]*(nn[i]+nn[j]+1)/2,
       zdobs[i][j], normpval(zdobs[i][j]));
   printf ("Maximum observed absolute Z-value:  %.4f\n", maxzstar);

   nsims=10000;
   printf ("\nSimulating P-values for MULTIPLE-COMPARISON-CORRECTED\n"
     "  comparisons for all pairs (using %s permutations of data):\n",
       commnum(nsims));

   /* Compare each observed Z(i,j) with the max Z(*,*) */
   /*   for random permutations */
   /* For safety, copy the observed data back to wwdat[] */
   for (r=i=0; i<kk; i++)
     for (j=0; j<nn[i]; j++)  wwdat[r++]=xx[i][j];
   /* Initialize the `nbig' array  nbigg[][]  */
   for (i=0; i<kk; i++)  for (j=0; j<kk; j++)  nbigg[i][j]=0;
   /* Carry out `nsims' permutations for all pairs */
   for (ns=0; ns<nsims; ns++)
    { double mstar=0.0;  /* Initialize mstar at 0 */
      /* Permute data */
      shuffle_dub(wwdat,ntot);
      /* Calculate  mstar = max |Z(*,*)| for the permuted data */
      for (i=0; i<kk; i++)  for (j=i+1; j<kk; j++)
       { int s, ni=nn[i], nj=nn[j];  double zstar;
         /* Copy treatment groups #i,#j to wmcbuff[] */
         for (r=s=0; s<ni; s++)  wmcbuff[r++]=wwdat[nfirst[i]+s];
         for (s=0; s<nj; s++)  wmcbuff[r++]=wwdat[nfirst[j]+s];
         zstar=get_wilcoxon(wmcbuff, ni,nj, NULL);
         if (fabs(zstar)>mstar) mstar=fabs(zstar); }
      /* mstar>=zdobs[i][j] is a `success' for pair (i,j) */
      for (i=0; i<kk; i++)  for (j=i+1; j<kk; j++)
       { /* Allow a little leeway for roundoff error */
         if (mstar > zdobs[i][j]-(1e-12))  nbigg[i][j]++; }}
   printf ("Pairwise MULTIPLE-COMPARISON-CORRECTED P-values:\n");
   /* starp(p) returns (*) for P<0.05 and (**) for P<0.01 */
   for (i=0; i<kk; i++)  for (j=i+1; j<kk; j++)
    { printf ("  #%d vs. #%d:  ObsZ=%.2f", i+1,j+1, zdobs[i][j]);
      pp=(double)nbigg[i][j]/nsims;
      printf ("   P = %5d/%s = %6.4f %s\n",
        nbigg[i][j],commnum(nsims), pp, starp(pp)); }
   printf ("\n(**) P<0.01    (*) 0.01<=P<0.05\n");

   printf (
   "\nNOTE: The text (p240-242) uses Wstar=Sqrt(2)*Z, not Z.\n"
   "Then Max Wstar(i,j) has a normal-range distribution for large ni\n"
   "Thus large-sample multiple-comparison-corrected Wilcoxon P-values\n"
   "  can be found by the probability that a normal-range statistic\n"
   "  is larger than Wstar(i,j). By definition, the normal-range\n"
   "  statistic is the maximum minus the minimum of k independent\n"
   "  N(0,1) variables, which can be simulated:\n"
   "Note that these P-values are similar to the exact MC-corrected\n"
   "  P-values above, but are more conservative.\n"
   "Pairwise MULTIPLE-COMPARISON-CORRECTED P-values via the\n"
   "  normal-range statistic:\n");
   /* Define Wdzobs=root(2)*zdobs and reset the counts nbigg[][] */
   for (i=0; i<kk; i++)  for (j=i+1; j<kk; j++)
    { wzdobs[i][j]=zdobs[i][j]*sqrt(2);  nbigg[i][j]=0; }
   /* Carry out `nsims' simulations of the normal-range P-values */
   /*   for wzdobs[][]: */
   for (ns=0; ns<nsims; ns++)
    { double rmin,rmax, zrange;
      /* Simulate one value of the normal-range statistic: */
      /* Rmin=min_k N_k,  Rmax=max_k N_k,  Range=Rmax-Rmin  */
      rmin=rmax=normrand();
      for (i=1; i<kk; i++)
       { double zz=normrand();
         if (zz<rmin) rmin=zz;  if (zz>rmax) rmax=zz; }
      zrange=rmax-rmin;
      /* zrange>=wzdobs[i][j] is a `success' for pair (i,j) */
      for (i=0; i<kk; i++)  for (j=i+1; j<kk; j++)
        if (zrange>wzdobs[i][j]-(1e-12))  nbigg[i][j]++; }
   for (i=0; i<kk; i++)  for (j=i+1; j<kk; j++)
    { pp=(double)nbigg[i][j]/nsims;
      printf ("  #%d vs. #%d:  W=%5.1f  Wstar=%.3f  P=%6d/%s = %6.4f %s\n",
        i+1,j+1, wdobs[i][j], wzdobs[i][j], nbigg[i][j],commnum(nsims),
        pp, starp(pp)); }

   /* Multiple comparisons WITH A CONTROL */
   printf ("\n\nYou can also do multiple comparisons WITH A CONTROL,\n"
   "  by assuming one of the treatment groups is a `null group'\n"
   "  and comparing all other treatment groups to it.\n"
   "\nThis can be done either TWO-SIDED (looking for absolute differences)\n"
   "  or (more commonly) ONE-SIDED (looking for larger or smaller),\n"
   "  in both cases in comparison with the control (or null) group.\n"
   "\nHere we consider treatment group #1 AS THE CONTROL and use a\n"
   "  ONE-SIDED procedure to look for treatment groups that are\n"
   "  SIGNIFICANTLY LARGER than the control. (See Section 6.7 in the text.)\n"
   "As before, we use normalized Wilcoxon rank-sum differences, and base\n"
   "  a multiple-comparison-corrected P-value for treatment #i on the\n"
   "  number of times that the maximum of treatment #i vs. control group\n"
   "  normalized comparisons is larger than the observe treatment #1 vs\n"
   "  control comparison. Since this is a one-sided procedure, absolute\n"
   "  values are not used.\n"
   "NOTE THAT this is the method of Steel (1959) (text p259) and NOT\n"
   "  the Nemenyi-Damico-Wolfe method discussed in most of Section 6.7.\n");

   /* We proceed as before, except that we compare all treatment groups */
   /*   with the control, and consider only larger than the control. */
   /* For clarity, we assume that the control group comes first, */
   /*   and in particular put it first in wmcbuff[] */
   ncon=nn[0];
   for (i=0; i<ncon; i++)  wmcbuff[i]=xx[0][i];
   /* maxzstar will be initialized at i=1 below */
   /* Note that maxzstar = max zd[0][i] can be negative */
   for (i=1; i<kk; i++)   /* All but the control group*/
    { int s, ni=nn[i];  double zstar;
      /* Copy the i-th treatment group after the control in wmcbuff[] */
      for (s=0; s<ni; s++)  wmcbuff[ncon+s]=xx[i][s];
      /* Compute the Wilcoxon ranksum statistic as before */
      zstar=get_wilcoxon(wmcbuff, ncon,ni, &(wdobs[0][i]));
      /* Save zstar and keep track of the minimum (more negative) zstar */
      zdobs[0][i]=zstar;
      /* The `i==1' allows maxzstar to be negative */
      if (i==1 || zstar>maxzstar) maxzstar=zstar; }
   printf (
   "\nFor comparison, the following are "
       "multiple-comparison-UNCORRECTED P-values.\n"
   "Pairwise Wilcoxon W and Z values:\n");
   for (i=1; i<kk; i++)
     printf ("  #1 vs. #%d:  W=%5.1f  E(W)=%5.1f  "
       "Z=%6.3f   P=%.4f (1-sided, UNCORR)\n",
       i+1, wdobs[0][i], (double)nn[i]*(nn[i]+ncon+1)/2,
       zdobs[0][i], 1-normint(zdobs[0][i]));
   printf ("Maximum observed Z-value:  %.4f\n", maxzstar);

   printf (
   "\nSimulating P-values for multiple-comparison-corrected comparisons\n"
     "  for control vs groups #2-#%d, using %s permutations:\n",
       kk, commnum(nsims));

   /* Compare each observed Z(1,i) with the max Z(1,*) */
   /*   for random permutations */
   /* Copy the observed data to wwdat[] */
   for (r=i=0; i<kk; i++)
     for (j=0; j<nn[i]; j++)  wwdat[r++]=xx[i][j];
   /* Clear the counts  nbigg[0][i] */
   for (i=0; i<kk; i++) nbigg[0][i]=0;
   /* Carry out `nsims' permutations */
   for (ns=0; ns<nsims; ns++)
    { int s;  double mstar;  /* mstar will be initialized at i=1 */
      /* Permute data */
      shuffle_dub(wwdat,ntot);
      /* Put the new control group first in wmcbuff[] */
      for (s=0; s<ncon; s++)  wmcbuff[s]=wwdat[s];
      /* Calculate  mstar = max Z(1,i) for the permuted data */
      for (i=1; i<kk; i++)
       { int s, ni=nn[i];  double zstar;
         /* Copy treatment #i to wmcbuff[] after #1 */
         for (s=0; s<ni; s++)  wmcbuff[ncon+s]=wwdat[nfirst[i]+s];
         zstar=get_wilcoxon(wmcbuff, ncon,ni, NULL);
         /* The `i==1' initializes mstar at zstar[1] */
         if (i==1 || zstar>mstar) mstar=zstar; }
      /* Add to counts of mstar >= Z(1,i) */
      for (i=1; i<kk; i++)
       { /* Allow a little leeway for roundoff error */
         if (mstar > zdobs[0][i]-(1e-12))  nbigg[0][i]++; }}
   printf ("Pairwise MULTIPLE-COMPARISON-CORRECTED P-values:\n");
   for (i=1; i<kk; i++)
    { printf ("  #1 vs. #%d:  ObsZ=%5.2f", i+1, zdobs[0][i]);
      pp=(double)nbigg[0][i]/nsims;
      printf ("   P = %6d/%s = %6.4f (1-sided) %s\n",
        nbigg[0][i],commnum(nsims), pp,starp(pp)); }
   printf ("\n(**) P<0.01    (*) 0.01<=P<0.05\n");
        }


/* Compute Hobs = Sum RankSum_i^2/nn_i  given offsets and nsizes */

double get_hobs (double *mrank, int *nfirst, int *nmax, int kk)
 { int i,j;  double hobs=0.0;
   for (i=0; i<kk; i++) if (nmax[i]>nfirst[i])
    { double rsum=0.0;
      for (j=nfirst[i]; j<nmax[i]; j++)  rsum+=mrank[j];
      hobs+=rsum*rsum/(nmax[i]-nfirst[i]); }
   return hobs; }


/* For samples x1[] of length n1  and  x2[] of length n2, */
/*   calculate the Wilcoxon rank-sum statistic W based on sample x2[] */
/*   and the large-sample normal approximation Z of W.  */
/* If  wstatptr!=NULL, return W at wstatptr, and return  Z */

double get_wilcoxon(const double *wboth, int nn1, int nn2, double *wstatptr)
 { int i, nboth=nn1+nn2, ntg, tgsize[MAXTOT];
   double wval, wmean,wvar, tsum,tcorr, zstar;
   double mrank[2*MAXN];
   /* Put midranks in mrank[] and tie-group sizes in tgsize[] */
   makeranks(mrank, wboth,nn1+nn2, tgsize,&ntg);
   /* Compute and save the Wilcoxon rank-sum statistic Wval */
   /*   based on sample #j */
   wval=0.0;  for (i=0; i<nn2; i++) wval+=mrank[nn1+i];
   if (wstatptr!=NULL) *wstatptr=wval;
   /* Find the large-sample Z-approximation Wstar of W: */
   /* The mean and variance of wval (without tie correction) */
   wmean=(double)(nn2*(nboth+1))/2;
   wvar=(double)(nn1*nn2*(nboth+1))/12;
   /* Now include the tie correction */
   tsum=0.0;
   for (i=0; i<ntg; i++)
    { int ts=tgsize[i];  tsum+=(ts-1)*ts*(ts+1); }
   tcorr = 1 - tsum/((double)(nboth-1)*nboth*(nboth+1));
   /* Tie-corrected zstar */
   zstar=(wval-wmean)/sqrt(wvar*tcorr);
   return zstar; }


/* Put midranks of array  wwdat[]  of length ntot */
/*   into the array  mrank[]  at corresponding offsets  */

double wtemp[MAXTOT];
double *mranptr[MAXTOT];

void makeranks_oneway(double *mrank, const double *wwdat, int ntot,
  int *tgsize, int *ntgptr)
 { int i,j,r, ntg;
   /* Make sure that the buffers  wtemp[],mranptr[]  are long enough */
   if (ntot>MAXTOT)    /* Then buffers will overflow */
    { printf ("\nERROR IN SUBROUTINE  MakeRanks!\n");
      printf ("Array length: %d    Buffer sizes: %d\n", ntot,MAXTOT);
      exit(1); }
   /* Copy wwdat[] to buffer  wtemp[]  so that it won't be corrupted */
   /* Also, initialize mranptr[] array to point to values in mrank[] */
   /* The next step will sort (wtemp[r],mranptr[r]) together. */
   /* The  mranptr[r]  value means that this subroutine will always know */
   /*   where to put the midrank of  wwdat[r]  */
   for (r=0; r<ntot; r++)
    { wtemp[r]=wwdat[r];  mranptr[r]=&(mrank[r]); }
   /* Bubble-sort the values (wtemp[r],mranptr[r]) together */
   for (i=0; i<ntot; i++) for (j=0; j<i; j++)
     if (wtemp[i]<wtemp[j])    /* Out of order, since j<i */
      { /* Switch values of  (wtemp[i],mranptr[i]) */
        /*   and  (wtemp[j],mranptr[j]) */
        double *wptr, wval=wtemp[i];  wtemp[i]=wtemp[j];  wtemp[j]=wval;
        wptr=mranptr[i];  mranptr[i]=mranptr[j];  mranptr[j]=wptr; }

   /* Assign midranks */
   ntg=0;   /* Number of tie-groups of size > 1 */
   for (i=0; i<ntot; i=j)
    { int k, ts;  double mr, firstval=wtemp[i];  /* First of tie group */
      /* The tie-group will be from i to j-1, inclusive */
      /* The first statement after the `do' is always executed */
      /*   at least once */
      /* The loop repeats as long as the `while' condition is TRUE */
      /* When the loop ends, j will point to the next value */
      j=i;  do j++; while(j<ntot && wtemp[j]==firstval);
      mr=(i+j+1)/2.0;  /* The average rank for the tie group */
      for (k=i; k<j; k++) *(mranptr[k])=mr;  /* Fill in the midranks */
      /* Accumulate tiegroup sums */
      ts=j-i;  /* The tie-group size */
      if (ts>1 && tgsize!=NULL)
        tgsize[ntg++]=ts; }  /* Keep track of tie-group sizes */
   /* If the user doesn't want the value ntg, don't insist. */
   if (ntgptr!=NULL) *ntgptr=ntg;
        }