/* truecluster/src/tieorder.c
 * A version of R's order() function breaking ties at random
 * (C) 2000-2007 Dr. Jens Oehlschlägel
 *
 * 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.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include <R_ext/Arith.h>        /* NA handling */
#include <Rmath.h>              /* fmin2 */
#include <R_ext/Random.h>       /* ..RNGstate */

#include <R_ext/Applic.h>       /* NA handling */

#include <R_ext/Utils.h>      /* sort routines et al */

#define RELTOL 1e-5

void R_tc_tieorder(
  double *x
, int *indx
, int *n_
)
{
  int tie, ntie, i, n, n2;
  double tol;
  Rboolean used_random=FALSE, thisNA, lastNA;
  double *randbuffer;
  int *posbuffer;
  n = *n_;
  // sort down with randomly breaking ties
  rsort_with_index(x, indx, n);
  tol = RELTOL * x[0];
  lastNA = ISNA(x[0]);
  ntie = 0;
  for (i=1; i<n; i++){
    thisNA = ISNA(x[i]);
    if ( (lastNA && thisNA) || ((x[i]-x[i-1])<=tol) ){
      if (!ntie){
        // start of equal series, let's break ties
        if (!used_random) {
          GetRNGstate(); used_random = TRUE;
          randbuffer = Calloc(sizeof(double)*n, double);
          posbuffer = Calloc(sizeof(int)*n, int);
        }
        posbuffer[ntie] = i-1;
        randbuffer[ntie] = unif_rand();
        //Rprintf("i=%d ntieon=%d rand=%f\n", i, ntie, randbuffer[ntie]);
        ntie++;
      }
      posbuffer[ntie] = i;
      randbuffer[ntie] = unif_rand();
      //Rprintf("i=%d ntie++=%d rand=%f\n", i, ntie, randbuffer[ntie]);
      ntie++;
    }else{
      if (ntie){
        //Rprintf("i=%d ntieoff=%d\n", i, ntie);
        // end of equal series, let's break ties
        rsort_with_index(randbuffer, posbuffer, ntie);
        // now copy x to randbuffer and indx to posbuffer and back
        for (tie=0; tie<ntie; tie++){
          //Rprintf("tie=%d rand=%f pos=%d\n", tie, randbuffer[tie], posbuffer[tie]);
          randbuffer[tie] = x[posbuffer[tie]];
          posbuffer[tie] = indx[posbuffer[tie]];
        }
        n2 = i - ntie;
        for (tie=0; tie<ntie; tie++){
          x[n2+tie] = randbuffer[tie];
          indx[n2+tie] = posbuffer[tie];
        }
        ntie = 0;
      }
    }
    lastNA = thisNA;
  }
  // don't forget the last one
  if (ntie){
    //Rprintf("i=%d lastntieoff=%d\n", i, ntie);
    // end of equal series, let's break ties
    rsort_with_index(randbuffer, posbuffer, ntie);
    // now copy x to randbuffer and indx to posbuffer and back
    for (tie=0; tie<ntie; tie++){
      //Rprintf("tie=%d rand=%f pos=%d\n", tie, randbuffer[tie], posbuffer[tie]);
      randbuffer[tie] = x[posbuffer[tie]];
      posbuffer[tie] = indx[posbuffer[tie]];
    }
    n2 = i - ntie;
    for (tie=0; tie<ntie; tie++){
      x[n2+tie] = randbuffer[tie];
      indx[n2+tie] = posbuffer[tie];
    }
    ntie = 0;
  }

  if(used_random) {
    PutRNGstate();
    Free(randbuffer);
    Free(posbuffer);
  }

}

