amar.c

ADaM is a data mining and image processing toolkit

                 indices[i] .. indices[j-1] point to values == pivot_value
                 indices[k] .. indices[end-1] point to values > pivot_value
           indices[start] .. indices[end-1] are a permutation of the indices on the
           way in */

      while ( j < k && farr[indices[j]] <= pivot_value )
      {
        while ( j < k && farr[indices[j]] < pivot_value )
        {
          if ( i < j )
            my_qs_swap(indices,i,j);
          i += 1;
          j += 1;
        }

        while ( j < k && farr[indices[j]] == pivot_value )
        {
          if ( i < j )
            my_qs_swap(indices,i,j);
          j += 1;
        }

      }

      while ( j < k && farr[indices[k-1]] > pivot_value )
        k -= 1;

      if ( j < k && farr[indices[k-1]] == pivot_value )
      {
        my_qs_swap(indices,j,k-1);
        j += 1;
        if ( j < k )
          k -= 1;
      }
      else if ( j < k )
      {
        if ( i == j )
        {
          my_qs_swap(indices,i,k-1);
          i += 1;
          j += 1;
          if ( j < k )
            k -= 1;
        }
        else
        {
          my_qs_swap(indices,j,k-1);
          my_qs_swap(indices,i,j);
          i += 1;
          j += 1;
          if ( j < k )
            k -= 1;
        }
      }

    }

    my_qs(farr,size,start,i,indices);
    my_qs(farr,size,k,end,indices);
  }
}

bool indices_are_sorted(double *farr,int *indices,int size)
{ /* XXX: This would be more efficient if the && were removed, and the "sorted =" were
     replaced with sorted *=     */
  bool sorted = TRUE;
  int i;
  for ( i = 0 ; sorted && i < size-1 ; i++ )
    sorted = farr[indices[i]] <= farr[indices[i+1]];
  return sorted;
}

void spr_indices_sort_realnums(double *farr, int size, int *indices)
/*
   PRE: farr and "indices" both have "size" elements.
        indices's contents irrelevant.
   POST: indices contains sorted indiceses into farr, so that
         forall j, indices[j] is the j'th smallest element of farr.

         thus forall i,j, (i < j) => farr[indices[i]] <= farr[indices[j]]
          and indices contains a permutation of [0 ...size-1]
*/
{
  int i;
  for ( i = 0 ; i < size ; i++ )
    indices[i] = i;

  /* All the following am_srand stuff is because the above 
     implementation of quicksort randomly selects its pivots and
     (a) I want to do the same pivot selection each time I sort
     the same array and (b) I want the random number generator
     state to be the same after a sort no matter what kind of sort
     I used. */
      
  push_current_am_srand_state();
  am_srand(575297);
  my_qs(farr,size,0,size,indices);
  pop_current_am_srand_state();

#ifndef AMFAST
  if ( !indices_are_sorted(farr,indices,size) )
    my_error("indices_sort_realnums failed");
#endif
}

void spr_sort_realnums(double *farr, int size, double *r_farr)
/*
   It is fine for farr to be the same memory as r_farr
*/
{
  int *indices = AM_MALLOC_ARRAY(int,size);
 
  spr_indices_sort_realnums(farr,size,indices);

  if ( farr == r_farr )
  {
    double *temp = AM_MALLOC_ARRAY(double,size);
    realnums_permute(farr,indices,temp,size);
    copy_realnums(temp,r_farr,size);
    AM_FREE_ARRAY(temp,double,size);
  }
  else
    realnums_permute(farr,indices,r_farr,size);

  AM_FREE_ARRAY(indices,int,size);
}

void os_sort_realnums(double *farr, int size, double *r_farr)
/*
   It is fine for farr to be the same memory as r_farr
*/
{
  copy_realnums(farr,r_farr,size);
  qsort((char *)r_farr,
        size,
        sizeof(double),
	double_comp
       );
}

double doubles_median(double *farr, int size)
{
  double *f2 = am_malloc_realnums(size);
  int med_index = size/2;
  double result;

  sort_realnums(farr,size,f2);

  if ( (size%2) == 0 )
    result = f2[med_index];
  else
    result = (f2[med_index] + f2[med_index+1])/2.0;

  AM_FREE_ARRAY(f2,double,size);

  return(result);
}

double *am_malloc_realnums(
    int size
  )
{
  double *result = AM_MALLOC_ARRAY(double,size);
#ifndef AMFAST
  set_realnums_constant(result,size,-7.7777e27);
#endif
  return(result);
}

void am_free_realnums(
  double *farr,
  int size
  )
{
  AM_FREE_ARRAY(farr,double,size);
}

int *am_malloc_ints(
    int size
  )
{
  int *result = AM_MALLOC_ARRAY(int,size);
  return(result);
}

void am_free_ints(
  int *iarr,
  int size
  )
{
  AM_FREE_ARRAY(iarr,int,size);
}

long *am_malloc_longs(
    int size
  )
{
  long *result = AM_MALLOC_ARRAY(long,size);
  return(result);
}

void am_free_longs(
  long *larr,
  int size
  )
{
  AM_FREE_ARRAY(larr,long,size);
}

bool *am_malloc_bools(
    int size
  )
{
  bool *result = AM_MALLOC_ARRAY(bool,size);
  return(result);
}

void am_free_bools(
  bool *barr,
  int size
  )
{
  AM_FREE_ARRAY(barr,bool,size);
}

int ints_argextreme(int *ints, int size, bool lowest)
{
  int i, r=0;
  int result;

  if (size <= 0)
    my_error("int_argextreme: zero (or -ve) size");
 
  result=ints[0];
  for (i=1; i<size; i++) {
    if ( lowest ? (ints[i] < result) : (ints[i] > result)) {
      r=i;
      result=ints[i];
    }
  }
  return r;
}

int ints_argmin(int *ints, int size)
{
  return ints_argextreme(ints,size,TRUE);
}

int ints_argmax(int *ints, int size)
{
  return ints_argextreme(ints,size,FALSE);
}

int doubles_argextreme(double *doubles, int size, bool lowest)
{
  int i, r=0;
  double result;
  if (size <= 0)
    my_error("double_argextreme: zero (or -ve) size");

  result = doubles[0];
  for (i=1; i<size; i++) {
    if ( lowest ? (doubles[i] < result) : (doubles[i] > result)) {
      r=i;
      result=doubles[i];
    }
  }
  return r;
}

int doubles_argmin(double *doubles, int size)
{
  return doubles_argextreme(doubles, size, TRUE);
}

int doubles_argmax(double *doubles, int size)
{
  return doubles_argextreme(doubles, size, FALSE);
}

/* indices_sort_realnums Courtesy of Justin Boyan.... */

typedef struct isf_elt_struct { int ix; double v; } isf_elt;

/* see __cdecl note at int_comp */
#ifdef WIN32
int __cdecl compare_isf_elt(const void *i_ptr, const void *j_ptr)
#else
int compare_isf_elt(const void *i_ptr, const void *j_ptr)
#endif
{
  isf_elt *i = (isf_elt *) i_ptr;
  isf_elt *j = (isf_elt *) j_ptr;
  return (i->v < j->v) ? -1 : (i->v > j->v) ? 1 : 0;
}

void os_indices_sort_realnums(double *farr, int size, int *indices)
/*
   PRE: farr and "indices" both have "size" elements.
        indices's contents irrelevant.
   POST: indices contains sorted indiceses into farr, so that
         forall j, indices[j] is the j'th smallest element of farr.

         thus forall i,j, (i < j) => farr[indices[i]] <= farr[indices[j]]
          and indices contains a permutation of [0 ...size-1]
*/
{
  int i;
  isf_elt *isf = AM_MALLOC_ARRAY(isf_elt, size);

  for (i=0; i<size; i++) {
    isf[i].ix = i;
    isf[i].v = farr[i];
  }
  qsort((char *)isf, size, sizeof(isf_elt), compare_isf_elt);
  for (i=0; i<size; i++) {
    indices[i] = isf[i].ix;
  }
  AM_FREE_ARRAY(isf, isf_elt, size);
}

int num_different_values_in_sorted_farr(double *farr,int size)
{
  int result = 1;
  int i;
  for ( i = 1 ; i < size ; i++ )
  {
    if ( farr[i] > farr[i-1] )
      result += 1;
  }
  return result;
}

int frequency_of_most_common_value_in_sorted_farr(double *farr,int size)
{
  int result = 1;
  int i;
  int count = 1;
  for ( i = 1 ; i < size ; i++ )
  {
    if ( farr[i] == farr[i-1] )
    {
      count += 1;
      result = int_max(result,count);
    }
    else
      count = 1;
  }
  return result;
}

bool avoid_os_sort_given_subset(double *subset,int subsize)
{
  bool result;
  os_sort_realnums(subset,subsize,subset);
  result = num_different_values_in_sorted_farr(subset,subsize) < 40 ||
           frequency_of_most_common_value_in_sorted_farr(subset,subsize) > subsize/2;
  if ( result )
    printf("Avoiding Microsoft's sort because of poor performance with duplicates\n");

  return result;
}

bool avoid_os_sort(double *farr,int size)
{
  bool result;
#ifndef PC_PLATFORM
  result = FALSE;
#else
  if ( size < 5000 )
    result = FALSE;
  else
  {
    int subsize = 200;
    double *subset = AM_MALLOC_ARRAY(double,subsize);
    int i;
    for ( i = 0 ; i < subsize ; i++ )
      subset[i] = farr[int_random(size)];
    result = avoid_os_sort_given_subset(subset,subsize);
    AM_FREE_ARRAY(subset,double,subsize);
  }
#endif

  return result;
}

bool avoid_os_int_sort(int *iarr,int size)
{
  bool result;
#ifndef PC_PLATFORM
  result = FALSE;
#else
  if ( size < 5000 )
    result = FALSE;
  else
  {
    int subsize = 200;
    double *subset = AM_MALLOC_ARRAY(double,subsize);
    int i;
    for ( i = 0 ; i < subsize ; i++ )
      subset[i] = (double) iarr[int_random(size)];
    result = avoid_os_sort_given_subset(subset,subsize);
    AM_FREE_ARRAY(subset,double,subsize);
  }
#endif
  return result;
}

/* by Artur ----------------------------------------------------------------- */
void os_indices_sort_integers(int *iarr, int size, int *indices)
/*
   PRE: iarr and "indices" both have "size" elements.
        indices's contents irrelevant.
   POST: indices contains sorted indiceses into iarr, so that
         forall j, indices[j] is the j'th smallest element of iarr.

         thus forall i,j, (i < j) => iarr[indices[i]] <= iarr[indices[j]]
          and indices contains a permutation of [0 ...size-1]
*/
{
  int i;
  isf_elt *isf = AM_MALLOC_ARRAY(isf_elt, size);
  for (i=0; i<size; i++) {
    isf[i].ix = i;
    isf[i].v = iarr[i];
  }
  qsort((char *)isf, size, sizeof(isf_elt), compare_isf_elt);
  for (i=0; i<size; i++) {
    indices[i] = isf[i].ix;
  }
  AM_FREE_ARRAY(isf, isf_elt, size);
}

double *mk_farr_from_iarr(int *iarr,int size)
{
  double *farr = AM_MALLOC_ARRAY(double,size);
  int i;
  for ( i = 0 ; i < size ; i++ )
    farr[i] = (double) iarr[i];
  return farr;
}

void spr_indices_sort_integers(int *iarr,int size,int *indices)
{
  double *farr = mk_farr_from_iarr(iarr,size);
  spr_indices_sort_realnums(farr,size,indices);
  AM_FREE_ARRAY(farr,double,size);
}

void spr_sort_ints(int *iarr, int size, int *r_iarr)
/*
   It is fine for iarr to be the same memory as r_iarr
*/
{
  int *indices = AM_MALLOC_ARRAY(int,size);
  double *farr = mk_farr_from_iarr(iarr,size);
  spr_indices_sort_realnums(farr,size,indices);

  if ( iarr == r_iarr )
  {
    int *temp = AM_MALLOC_ARRAY(int,size);
    ints_permute(iarr,indices,temp,size);
    copy_ints(temp,r_iarr,size);
    AM_FREE_ARRAY(temp,int,size);
  }
  else
    ints_permute(iarr,indices,r_iarr,size);

  AM_FREE_ARRAY(indices,int,size);
  AM_FREE_ARRAY(farr,double,size);
}

/*
   PRE: farr and "indices" both have "size" elements.
        indices's contents irrelevant.
   POST: indices contains sorted indiceses into farr, so that
         forall j, indices[j] is the j'th smallest element of farr.

         thus forall i,j, (i < j) => farr[indices[i]] <= farr[indices[j]]
          and indices contains a permutation of [0 ...size-1]
*/
void indices_sort_realnums(double *farr, int size, int *indices)
{
  if ( avoid_os_sort(farr,size) )
    spr_indices_sort_realnums(farr,size,indices);
  else
    os_indices_sort_realnums(farr,size,indices);
}

void sort_realnums(double *farr, int size, double *r_farr)
{
  if ( avoid_os_sort(farr,size) )
    spr_sort_realnums(farr,size,r_farr);
  else
    os_sort_realnums(farr,size,r_farr);
}

/*
   PRE: farr and "indices" both have "size" elements.
        indices's contents irrelevant.
   POST: indices contains sorted indiceses into farr, so that
         forall j, indices[j] is the j'th smallest element of farr.

         thus forall i,j, (i < j) => farr[indices[i]] <= farr[indices[j]]
          and indices contains a permutation of [0 ...size-1]
*/
void indices_sort_integers(int *iarr, int size, int *indices)
{
  if ( avoid_os_int_sort(iarr,size) )
    spr_indices_sort_integers(iarr,size,indices);
  else
    os_indices_sort_integers(iarr,size,indices);
}

void sort_ints(int *iarr, int size, int *r_iarr)
{
  if ( avoid_os_int_sort(iarr,size) )
    spr_sort_ints(iarr,size,r_iarr);
  else
    os_sort_ints(iarr,size,r_iarr);
}