amar.c

ADaM is a data mining and image processing toolkit

/*
  Logistic Regression using Truncated Iteratively Re-weighted Least Squares
  (includes several programs)
  Copyright (C) 2005  Paul Komarek

  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

  Author: Paul Komarek, [e-mail]
  Alternate contact: Andrew Moore, [e-mail]

*/



/* !!!!
   File:        amar.c
   Author:      Andrew W. Moore
   Created:     Sat Sep 12 14:53:13 EDT 1992
   Description: Obvious operations on 1-d arrays

   Copyright (C) 1992, Andrew W. Moore
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "amar.h"      /* Obvious operations on 1-d arrays */

/* Private prototypes */
void os_sort_ints(int *farr, int size, int *r_farr);

#ifdef WIN32
int __cdecl double_comp(const void *r1_ptr, const void *r2_ptr);
int __cdecl compare_isf_elt(const void *i_ptr, const void *j_ptr);
#else
int double_comp(const void *r1_ptr, const void *r2_ptr);
int compare_isf_elt(const void *i_ptr, const void *j_ptr);
#endif

void realnums_permute(double *farr,int *indices,double *r_farr,int size);
void ints_permute(int *iarr,int *indices,int *r_iarr,int size);
void my_qs_2(double *farr, int start, int *indices);
void my_qs(double *farr, int size, int start, int end, int *indices);
bool indices_are_sorted(double *farr,int *indices,int size); /* This maybe should be public? */
void spr_indices_sort_realnums(double *farr, int size, int *indices);
void os_sort_realnums(double *farr, int size, double *r_farr);
int ints_argextreme(int *ints, int size, bool lowest);
int doubles_argextreme(double *doubles, int size, bool lowest);
void spr_sort_realnums(double *farr, int size, double *r_farr);
void os_indices_sort_realnums(double *farr, int size, int *indices);
int num_different_values_in_sorted_farr(double *farr,int size); /* Useful. Should be public. */
int frequency_of_most_common_value_in_sorted_farr(double *farr,int size);
bool avoid_os_sort_given_subset(double *subset,int subsize);
bool avoid_os_sort(double *farr,int size);
bool avoid_os_int_sort(int *iarr,int size);
void os_indices_sort_integers(int *iarr, int size, int *indices);
double *mk_farr_from_iarr(int *iarr,int size); /* Should be public */
void spr_indices_sort_integers(int *iarr,int size,int *indices);
void spr_sort_ints(int *iarr, int size, int *r_iarr);



void copy_ints(int *src, int *dest, int size)
{
  int i;
  for ( i = 0 ; i < size ; i++ )
    dest[i] = src[i];
}

void fprintf_ints(FILE *s, char *message1, int *ints, int size, char *message2)
{
  int i;
  fprintf(s,"%s = ",message1);

  if ( ints == NULL && size > 0 )
    fprintf(s,"(int *) NULL");
  else
  {
    fprintf(s,"{ ");
    for ( i = 0 ; i < size ; i++ )
      fprintf(s,"%d %s",ints[i],(i==size-1) ? "" : ", ");
    fprintf(s,"}");
  }

  fprintf(s," %s",message2);
}

void printf_ints(char *message1, int *ints, int size, char *message2)
{
  fprintf_ints(stdout,message1,ints,size,message2);
}

int ints_extreme(int *ints, int size, bool lowest)
{
  int result;
  if ( size <= 0 )
  {
    result = 0;
    my_error("ints_extreme: zero (or -ve) size");
  }
  else
  {
    int i;
    result = ints[0];
    for ( i = 1 ; i < size ; i++ )
      result = (lowest) ? int_min(result,ints[i]) : int_max(result,ints[i]);
  }
  return(result);
}

int ints_min(int *ints, int size)
{
  return(ints_extreme(ints,size,TRUE));
}

int ints_max(int *ints, int size)
{
  return(ints_extreme(ints,size,FALSE));
}

/*
* Frank Dellaert note (6/23/97)
* qsort expects cdecl, but the mblc DLL is compiled using
* __stdcall calling conventions
*/
#ifdef WIN32
int __cdecl int_comp(const void *r1_ptr, const void *r2_ptr)
#else
int int_comp(const void *r1_ptr, const void *r2_ptr)
#endif
{
  int *r1 = (int *) r1_ptr;
  int *r2 = (int *) r2_ptr;
  if ( *r1 < *r2 )
    return(-1);
  else if ( *r1 > *r2 )
    return(1);
  else
    return(0);
}

void os_sort_ints(int *farr, int size, int *r_farr)
/*
   It is fine for farr to be the same memory as r_farr
*/
{
  copy_ints(farr,r_farr,size);
  qsort((char *)r_farr,
        size,
        sizeof(int),
	int_comp
       );
}

void copy_realnums(double *src, double *dest, int size)
{
  int i;
  for ( i = 0 ; i < size ; i++ )
    dest[i] = src[i];
}

void fprintf_realnums(FILE *s, char *message1, double *doubles, int size, char *message2)
{
  int i;
  fprintf(s,"%s ",message1);

  if ( doubles == NULL && size > 0 )
    fprintf(s,"(double *) NULL");
  else
  {
    fprintf(s,"{ ");
    for ( i = 0 ; i < size ; i++ )
      fprintf(s,"%g %s",doubles[i],(i==size-1) ? "" : ", ");
    fprintf(s,"}");
  }

  fprintf(s," %s",message2);
}

void printf_realnums(char *message1, double *doubles, int size, char *message2)
{
  fprintf_realnums(stdout,message1,doubles,size,message2);
}

void set_realnums_constant(double *doubles, int size, double val)
{
  int i;
  for ( i = 0 ; i < size ; i++ )
    doubles[i] = val;
}

double doubles_extreme(double *doubles, int size, bool lowest)
{
  double result;
  if ( size <= 0 )
  {
    result = 0.0;
    my_error("doubles_extreme: zero (or -ve) size");
  }
  else
  {
    int i;
    result = doubles[0];
    for ( i = 1 ; i < size ; i++ )
      result = (lowest) ? 
               real_min(result,doubles[i]) : real_max(result,doubles[i]);
  }
  return(result);
}

double doubles_min(double *doubles, int size)
{
  return(doubles_extreme(doubles,size,TRUE));
}

double doubles_max(double *doubles, int size)
{
  return(doubles_extreme(doubles,size,FALSE));
}

double doubles_inner_product(double *f1, double *f2, int size)
{
  int i;
  double result = 0.0;
  for ( i = 0 ; i < size ; i++ )
    result += f1[i] * f2[i];
  return(result);
}

double doubles_magnitude_sqd(double *f, int size)
{
  return(doubles_inner_product(f,f,size));
}

void doubles_scalar_multiply(double *f, int size, double alpha, double *res_f)
{
  int i;
  for ( i = 0 ; i < size ; i++ )
   res_f[i] = alpha * f[i];
}

void doubles_add(double *f1, double *f2, int size, double *res_f)
{
  int i;
  for ( i = 0 ; i < size ; i++ )
   res_f[i] = f1[i] + f2[i];
}

void doubles_subtract(double *f1, double *f2, int size, double *res_f)
{
  int i;
  for ( i = 0 ; i < size ; i++ )
   res_f[i] = f1[i] - f2[i];
}

double doubles_scaled_dsqd(double *f1, double *f2, double *scales, int size)
/* Returns scales . ( f1 - f2 ) (. = dot-product) */
{
  double result = 0.0;
  int i;
  for ( i = 0 ; i < size ; i++ )
  {
    double d = scales[i] * (f1[i] - f2[i]);
    result += d * d;
  }
  return(result);
}

double doubles_antiscaled_dsqd(double *f1, double *f2, double *scales, int size)
/* Returns SUM ( f1[i] - f2[i] ) / scales[i] */
{
  double result = 0.0;
  int i;
  for ( i = 0 ; i < size ; i++ )
  {
    double d = (f1[i] - f2[i]) / scales[i];
    result += d * d;
  }
  return(result);
}

void doubles_reciprocal(double *farr, int size, double *res_farr)
{
  int i;
  for ( i = 0 ; i < size ; i++ )
  {
    if ( fabs(farr[i]) < 1e-20 )
      printf("doubles_reciprocal about to hit a near-zero\n");

    res_farr[i] = 1.0 / farr[i];
  }
}

void doubles_scalar_add(double *farr1, double *farr2, int size, double alpha, double *res_farr)
/*
   Any of the input arrays may share memory.
   res_farr := farr1 + alpha * farr2
*/
{
  int i;
  for ( i = 0 ; i < size ; i++ )
    res_farr[i] = farr1[i] + alpha * farr2[i];
}

/* 2d_arrays. The entry at the ith row and jth colum is index as
   x[i][j]. The 2d arrays are implemented as an array of pointers
   to arrays of doubles. indexing is [0..rows-1] [0..columns-1] inclusive
*/

double **am_malloc_2d_realnums(int rows, int cols)
{
  double **result = AM_MALLOC_ARRAY(double_ptr,rows);
  int i;

  for ( i = 0 ; i < rows ; i++ )
    result[i] = am_malloc_realnums(cols);

  return(result);
}

void doubles_2d_scalar_add(double **tdarr1, double **tdarr2, int rows,
 int cols, double alpha, double **res_tdarr)
/*
   Any of the input arrays may share memory.
   res_tdarr := tdarr1 + alpha * tdarr2
*/
{
  int i;
  for ( i = 0 ; i < rows ; i++ )
    doubles_scalar_add(tdarr1[i],tdarr2[i],cols,alpha,res_tdarr[i]);
}

void fprintf_2d_realnums(FILE *s, char *m1, double **tdarr,
 int rows, int cols, char *m2)
{
  char buff[100];
  int i;

  for ( i = 0 ; i < rows ; i++ )
  {
    sprintf(buff,"%s[%d] = ",m1,i);
    fprintf_realnums(s,buff,tdarr[i],cols,m2);
  }
}

void free_2d_realnums(double **tdarr, int rows, int cols)
{
  int i;
  for ( i = 0 ; i < rows ; i++ )
    AM_FREE_ARRAY(tdarr[i],double,cols);
  AM_FREE_ARRAY(tdarr,double_ptr,rows);
}

void am_free_2d_realnums(double **tdarr, int rows, int cols)
/* Exactly the same as "free_2d_realnums". Added for naming consistency. */
{
  free_2d_realnums(tdarr,rows,cols);
}

void set_2d_realnums_constant(double **tdarr, int rows, int cols, double val)
{
  int i;
  for ( i = 0 ; i < rows ; i++ )
    set_realnums_constant(tdarr[i],cols,val);
}

void copy_2d_realnums(double **tdarr1, double **tdarr2, int rows, int cols)
{
  int i;
  for ( i = 0 ; i < rows ; i++ )
    copy_realnums(tdarr1[i],tdarr2[i],cols);
}

double doubles_sum(double *farr, int size)
{
  double result = 0.0;
  int i;
  for ( i = 0 ; i < size ; i++ )
    result += farr[i];
  return(result);
}

double doubles_mean(double *farr, int size)
{
  return(doubles_sum(farr,size) / int_max(size,1));
}

double doubles_sdev(double *farr, int size)
{
  double mean = doubles_mean(farr,size);
  int i;
  double sum_sqs = 0.0;
  for ( i = 0 ; i < size ; i++ )
    sum_sqs += real_square(farr[i] - mean);

  return(sqrt(sum_sqs / int_max(size,1)));
}

void random_shuffle_ints(int *iarr, int size)
{
  int i;
  for ( i = 0 ; i < size-1 ; i++ )
  {
    int j = i + int_random(size-i);
    int temp = iarr[i];
    iarr[i] = iarr[j];
    iarr[j] = temp;
  }
}


/* see __cdecl note at int_comp */
#ifdef WIN32
int __cdecl double_comp(const void *r1_ptr, const void *r2_ptr)
#else
int double_comp(const void *r1_ptr, const void *r2_ptr)
#endif
{
  double *r1 = (double *) r1_ptr;
  double *r2 = (double *) r2_ptr;
  if ( *r1 < *r2 )
    return(-1);
  else if ( *r1 > *r2 )
    return(1);
  else
    return(0);
}

/* PRE:  farr and r_farr must be different bits of memory.
         indices must be a permutation of 0,1,...size-1 
   POST: r_farr[i] = farr[indices[i]] */
void realnums_permute(double *farr,int *indices,double *r_farr,int size)
{
  int i;
  for ( i = 0 ; i < size ; i++ )
    r_farr[i] = farr[indices[i]];
}

/* PRE:  iarr and r_iarr must be different bits of memory.
         indices must be a permutation of 0,1,...size-1 
   POST: r_iarr[i] = iarr[indices[i]] */
void ints_permute(int *iarr,int *indices,int *r_iarr,int size)
{
  int i;
  for ( i = 0 ; i < size ; i++ )
    r_iarr[i] = iarr[indices[i]];
}

/*********** SPR's own QuickSort implementation.

    AWM: It would usually be insane to use anything except the OS's
    built in sort but it so happens that microsoft's qsort goes
    slowly when there are lots of duplicates. The following implementation
    behaves very well with many duplicates. 

    Later on you'll see that our bottom level sorting chooses whether or
    not to use the OS sort according to whether there are many duplicates
    or not
*/

#define my_qs_swap(indices,i,j) \
{ \
  int temp = indices[i]; \
  indices[i] = indices[j]; \
  indices[j] = temp; \
}

void my_qs_2(double *farr, int start, int *indices)
{
  double x0 = farr[indices[start+0]];
  double x1 = farr[indices[start+1]];
  if ( x1 < x0 )
    my_qs_swap(indices,start,start+1);
}
  
void my_qs(double *farr, int size, int start, int end, int *indices)
{
  if ( end <= start+1 )
  {
    /* skip */
  }
  else if ( end == start + 2 )
  {
    my_qs_2(farr,start,indices);
  }
  else if ( end == start + 3 )
  {
    double x0 = farr[indices[start+0]];
    double x1 = farr[indices[start+1]];
    double x2 = farr[indices[start+2]];
    if ( x0 > x1 || x0 > x2 )
    {
      if ( x1 < x2 )
      {
        my_qs_swap(indices,start,start+1);
        my_qs_2(farr,start+1,indices);
      }
      else
      {
        my_qs_swap(indices,start,start+2);
        my_qs_2(farr,start+1,indices);
      }
    }
    else
      my_qs_2(farr,start+1,indices);
  }
  else
  {
    int pivot_index = start + int_random(end-start);
    double pivot_value = farr[indices[pivot_index]];
    int i = start;
    int j = start;
    int k = end;

    while ( j < k )
    {
      /* Assert: indices[start] ... indices[i-1] point to values < pivot_value