amiv.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:        amiv.c
   Author:      Andrew W. Moore
   Created:     Sat Apr  8 18:48:26 EDT 1995
   Updated:     29 March 97
   Description: Integer and Boolean Dynamic vectors

   Copyright 1997, Schenley Park Research
*/

#include "amiv.h"
#include "am_string.h"
#include "am_string_array.h"
#include "amma.h"

#define IVEC_CODE 20541

int Ivecs_mallocked = 0;
int Ivecs_freed = 0;

/* Headers for private functions. Are these all meant to be private? */
int ivec_arg_extreme(const ivec *iv,bool lowest);
int ivec_extreme(const ivec *iv,bool lowest);

ivec *mk_ivec(int size)
{
  ivec *result = AM_MALLOC(ivec);
  if ( size < 0 ) my_error("mk_ivec : size < 0 illegal");
  result -> ivec_code = IVEC_CODE;
  result -> array_size = size;
  result -> size = size;
  result -> iarr = am_malloc_ints(size);
  Ivecs_mallocked += 1;
  return(result);
}

/* Acts as though you created a ivec of size 0,        */
/* but actually allocates an array of size "capacity". */
/* This is very useful if you want to use add_to_ivec  */
/* and happen to have a reasonable upper bound to the  */
/* number of elements you want to add.                 */
ivec *mk_empty_ivec(int capacity) {
  ivec *result = AM_MALLOC(ivec);
  if ( capacity < 0 ) my_error("mk_empty_ivec : capacity < 0 illegal");
  result -> ivec_code = IVEC_CODE;
  result -> array_size = capacity;
  result -> size = 0;
  result -> iarr = am_malloc_ints(capacity);
  Ivecs_mallocked += 1;
  return(result);
}


ivec *mk_ivec_x( int size, ...)
{
  /* Warning: no type checking can be done by the compiler.  You *must*
     send the values as integers for this to work correctly. */
  int i, ival;
  va_list argptr;
  ivec *iv;
  
  iv = mk_ivec(size);

  va_start( argptr, size);
  for (i=0; i<size; ++i) {
    ival = va_arg( argptr, int);
    ivec_set( iv, i, ival);
  }
  va_end(argptr);

  return iv;
}

void free_ivec(ivec *iv)
{
  iv -> ivec_code = 7777;

  am_free_ints(iv->iarr,iv->array_size);
  AM_FREE(iv,ivec);

  Ivecs_freed += 1;
}

void fprintf_ivec(FILE *s,char *m1, const ivec *iv,char *m2)
{
  if ( iv == NULL )
    fprintf(s,"%s = (ivec *)NULL%s",m1,m2);
  else
  {
    char buff[256];

    /* Modified so that we don't cause a seg fault when m1 is bigger than
       the buffer - WKW*/
    if( strlen(m1) <= 255 )
      {
	strncpy(buff,m1,strlen(m1)+1);
      }
    else
      {
	strncpy(buff,m1,sizeof(buff));
	buff[255] = '\0';
      }
    fprintf_ints(s,buff,iv->iarr,iv->size,m2);
  }
}

void copy_iarr_to_ivec(int *iarr,int size,ivec *r_iv)
{
  copy_ints(iarr,r_iv->iarr,size);
}

void copy_ivec_to_iarr(ivec *iv, int *iarr)
{
  copy_ints(iv->iarr,iarr,iv->size);
}

int *mk_iarr_from_ivec(ivec *iv)
{
  int *result;
  result = am_malloc_ints(iv->size);
  copy_ivec_to_iarr(iv,result);
  return(result);
}

/* Makes an ivec of size end - start:
   { start , start+1 , .... end-2 , end-1 } */
ivec *mk_sequence_ivec(int start_value,int end_value)
{
  int size = end_value - start_value;
  ivec *iv = mk_ivec(size);
  int i;
  for ( i = 0 ; i < size ; i++ )
    ivec_set(iv,i,start_value+i);
  return iv;
}

/*
   Allocates and returns an ivec of size size in which ivec[i] = i
   ivec[0] = 0
   ivec[1] = 1
    .
    .
   ivec[size-1] = size-1
*/
ivec *mk_identity_ivec(int size)
{
  return mk_sequence_ivec(0,size);
}

void shuffle_ivec(ivec *iv)
/*
   A random permutation of iv is returned
*/
{
  int size = ivec_size(iv);
  int i;
  for ( i = 0 ; i < size-1 ; i++ )
  {
    int j = int_random(size - i);
    if ( j > 0 )
    {
      int swap_me_1 = ivec_ref(iv,i);
      int swap_me_2 = ivec_ref(iv,i+j);
      ivec_set(iv,i,swap_me_2);
      ivec_set(iv,i+j,swap_me_1);
    }
  }
}

void constant_ivec(ivec *iv,int value)
{
  int i;
  for ( i = 0 ; i < ivec_size(iv) ; i++ )
    ivec_set(iv,i,value);
}

ivec *mk_constant_ivec(int size,int value)
{
  ivec *iv = mk_ivec(size);
  constant_ivec(iv,value);
  return(iv);
}

void zero_ivec(ivec *iv)
{
  constant_ivec(iv,0);
}

ivec *mk_zero_ivec(int size)
{
  return(mk_constant_ivec(size,0));
}

void ivec_scalar_plus(const ivec *iv, int addme, ivec *riv)
{
  int i;
  for ( i = 0 ; i < ivec_size(iv) ; i++ )
    ivec_set(riv,i,ivec_ref(iv,i) + addme);
}

ivec *mk_ivec_scalar_plus(ivec *iv,int delta)
{
  ivec *result = mk_ivec(ivec_size(iv));
  ivec_scalar_plus(iv,delta,result);
  return(result);
}

void ivec_scalar_mult(const ivec *iv,int scale,ivec *riv)
{
  int i;
  for ( i = 0 ; i < ivec_size(iv) ; i++ )
    ivec_set(riv,i,ivec_ref(iv,i) * scale);
}

ivec *mk_ivec_scalar_mult(ivec *iv,int scale)
{
  ivec *result = mk_ivec(ivec_size(iv));
  ivec_scalar_mult(iv,scale,result);
  return(result);
}

void copy_ivec(const ivec *iv,ivec *r_iv)
{
  int i, size;
  size = iv->size;
  for (i=0; i<size; ++i) ivec_set(r_iv, i, ivec_ref(iv, i));
  return;
  /* ivec_scalar_mult(iv,1,r_iv); */
}

ivec *mk_copy_ivec(const ivec *iv)
{
  ivec *result = mk_ivec(ivec_size(iv));
  copy_ivec(iv,result);
  return(result);
}

void copy_ivec_subset(ivec *iv, int start, int size, ivec *r_iv)
{
  int i;
  for (i = start; i < start + size; i++)
    ivec_set(r_iv, i - start, ivec_ref(iv, i));
}

ivec *mk_copy_ivec_subset(ivec *iv, int start, int size)
{
  ivec *result = mk_ivec(size);
  copy_ivec_subset(iv, start, size, result);
  return (result);
}

int num_of_given_value(ivec *iv,int value)
{
  int i;
  int result = 0;
  for ( i = 0 ; i < ivec_size(iv) ; i++ )
    if ( ivec_ref(iv,i) == value ) result += 1;

  return(result);
}

int num_zero_entries(ivec *iv)
{
  return(num_of_given_value(iv,0));
}

int num_nonzero_entries(ivec *iv)
{
  return(ivec_size(iv) - num_zero_entries(iv));
}

int ivec_arg_extreme(const ivec *iv,bool lowest)
{
  int extreme_index;
  int extremest_value;
  int i;

  if ( ivec_size(iv) <= 0 )
    my_error("Can't find min or max of empty ivec");
  
  extreme_index = 0;
  extremest_value = ivec_ref(iv,0);

  for ( i = 1 ; i < ivec_size(iv) ; i++ )
  {
    int v = ivec_ref(iv,i);
    if ( (lowest && v < extremest_value) || (!lowest && v > extremest_value) ) 
    {
      extremest_value = v;
      extreme_index = i;
    }
  }

  return(extreme_index);
}

int ivec_extreme(const ivec *iv,bool lowest)
{
  return ivec_ref(iv,ivec_arg_extreme(iv,lowest));
}

int ivec_min(const ivec *iv)
{
  return ivec_ref(iv,ivec_arg_extreme(iv,TRUE));
}

int ivec_max(const ivec *iv)
{
  return ivec_ref(iv,ivec_arg_extreme(iv,FALSE));
}

int ivec_argmin(const ivec *iv)
{
  return ivec_arg_extreme(iv,TRUE);
}

int ivec_argmax(const ivec *iv)
{
  return ivec_arg_extreme(iv,FALSE);
}

/* Sensible if args are NULL. False if different size */
bool ivec_equal(const ivec *x1,const ivec *x2)
{
  bool result = TRUE;

  if ( EQ_PTR(x1,x2) )
    result = TRUE;
  else if ( x1 == NULL || x2 == NULL )
    result = FALSE;
  else if ( ivec_size(x1) != ivec_size(x2) )
    result = FALSE;
  else
  {
    int i;
    for ( i = 0 ; result && i < ivec_size(x1) ; i++ ) 
      result = result && ivec_ref(x1,i) == ivec_ref(x2,i);
  }
  return(result);
}

/**** Removal functions on ivecs ****/

/* Reduces the size of iv by one.
   ivec_ref(iv,index) disappears.
   Everything to the right of ivec_ref(iv,index) is copied one to the left.

   Formally: Let ivold be the ivec value beore calling this function
             Let ivnew be the ivec value after calling this function.

PRE: ivec_size(ivold) > 0
     0 <= index < ivec_size(ivold)

POST: ivec_size(ivnew) = ivec_size(ivold)-1
      for j = 0 , 1, 2 ... index-1  : 
         ivec_ref(ivnew,j) == ivec_ref(ivold,j)

      for j = index , index+1 , ... ivec_size(ivnew)-1:
         ivec_ref(ivnew,j) == ivec_ref(ivold,j+1)
*/
void ivec_remove(ivec *iv,int idx)
{
  int i;
  int isize = ivec_size(iv);

#ifndef AMFAST
  if ( isize <= 0 ) my_error("ivec_remove: empty ivec");
  if ( idx < 0 || idx >= isize )
    my_error("ivec_remove: bad index");
#endif /* #ifndef AMFAST */

  for ( i = idx ; i < isize - 1 ; i++ )
    ivec_set(iv,i,ivec_ref(iv,i+1));
  iv -> size -= 1;
}

int ivec_sum(const ivec *iv)
{
  int result = 0;
  int i;
  for ( i = 0 ; i < ivec_size(iv) ; i++ )
    result += ivec_ref(iv,i);
  return(result);
}

bool equal_ivecs(const ivec *iv1, const ivec *iv2)
{
  bool result = TRUE;
  int i;

  if ( ivec_size(iv1) != ivec_size(iv2) )
    my_error("qual_ivecs wrong sizes");

  for ( i = 0 ; i < ivec_size(iv1) && result ; i++ )
    result = ivec_ref(iv1,i) == ivec_ref(iv2,i);
  return(result);
}

void add_to_ivec(ivec *iv,int new_val)
{
  int size;


  if ( iv->array_size < 0 )
    my_error("An ivec should never have a negative array_size if\n"
	     "it was constructed using legal API calls");

  if ( iv->size > iv->array_size )
    my_error("An ivec should never have size greater than array_size if\n"
    "it was constructed using legal API calls");

  size = iv->size;

  if ( size == iv->array_size )
  {
    int new_array_size = 2 * size + 10;
    int *iarr_new = AM_MALLOC_ARRAY(int,new_array_size);
    int i;

    for ( i = 0 ; i < size ; i++ )
      iarr_new[i] = iv->iarr[i];
    AM_FREE_ARRAY(iv->iarr,int,size);
    iv -> iarr = iarr_new;
    iv -> array_size = new_array_size;
  }
  iv->iarr[size] = new_val;
  iv -> size += 1;
}


void add_to_ivec_unique(ivec *iv,int val) {
  printf("AWM comments: You should probably be using sivecs\n");
  if (!is_in_ivec(iv, val)) add_to_ivec(iv, val);
}

/* Increases iv in length by 1 and shifts all elements
   with original index greater or equal to index one to the
   right and inserts val at index. */
void ivec_insert(ivec *iv,int idx,int val)
{
  int i;
  add_to_ivec(iv,0);
  for ( i = ivec_size(iv)-1 ; i > idx ; i-- )
    ivec_set(iv,i,ivec_ref(iv,i-1));
  ivec_set(iv,idx,val);
}

/* Creates a new ivec with the given element inserted.
   This is more efficient than copying and inserting. */
ivec *mk_ivec_insert( ivec *iv, int idx, int val)
{
  int i, newsize;
  ivec *result;

  newsize = ivec_size( iv) + 1;
  result = mk_ivec( newsize);

#ifndef AMFAST
  if (idx >= newsize || idx < 0) {
    my_errorf( "mk_ivec_insert: idx=%d is out-of-bounds [0,%d]",
               idx, newsize-1);
  }
#endif

  for (i=0; i<idx; ++i) ivec_set( result, i, ivec_ref( iv, i));
  ivec_set( result, idx, val);
  for (i=idx+1; i<newsize; ++i) ivec_set( result, i, ivec_ref( iv, i-1));
  
  return result;
}



/* Find least index i such that value = ivec_ref(iv,i).
  If not found, returns -1
*/
int find_index_in_ivec(const ivec *iv, int value)
{
  int result = -1;
  int i;

  for ( i = 0 ; i < ivec_size(iv) && result < 0 ; i++ )
    if (value == ivec_ref(iv,i)) 
      result = i;
  return(result);
}

/* make and return an ivec containing the indices of iv that have value */
ivec *mk_indices_in_ivec(const ivec *iv, int value)
{
  int i;
  ivec *res = mk_ivec(0);
  for (i=0;i<ivec_size(iv);i++) if (ivec_ref(iv,i)==value) add_to_ivec(res,i);
  return res;
}

int find_index_in_ivec_hint(const ivec *iv, int value, int hint)
{
  int i;
  int sign = -1;
  int sz = ivec_size(iv);

  for (i=1; i<=sz; i++) {
	int disp = i/2;
	int idx = (hint + sign*disp);
    /* make sure idx is in bounds. Adding sz will take care of negative idx.*/
    idx = (idx + sz) % sz;
    if (value == ivec_ref(iv, idx)) 
      return idx;
    sign *= -1;
  }
  return -1;
}

/* Finds leftmost index such that siv[index] > val. If none
   such, returns ivec_size(siv) */
int index_in_sorted_ivec(ivec *siv,int val)
{
  int result = -1;
  int i;
  for ( i = 0 ; i < ivec_size(siv) && result < 0 ; i++ )
    if ( ivec_ref(siv,i) > val )
      result = i;
  if ( result < 0 ) result = ivec_size(siv);
  return(result);
}

int find_in_sorted_ivec(ivec *siv, int val)
{
  int begin = 0;
  int end = ivec_size(siv);
  int cur = (begin+end)/2;
  while (end > begin)
  {
    if (ivec_ref(siv,cur) == val) return cur;
    else if (ivec_ref(siv,cur) < val) begin = cur+1;
    else                              end = cur;
    cur = (begin+end)/2;
  }
  return -1;
}

void add_to_sorted_ivec(ivec *siv, int val)
{
  int i;
  add_to_ivec(siv, 0);
  for (i = ivec_size(siv) - 2; i >= 0 && val < ivec_ref(siv, i); i--)
    ivec_set(siv, i + 1, ivec_ref(siv, i));
  ivec_set(siv, i + 1, val);
}

bool is_in_ivec( const ivec *iv,int value)
{
  return(find_index_in_ivec(iv,value) >= 0);
}

/* Returns list of indices at which value appears. */
ivec *mk_is_in_ivec( ivec *iv, int value)
{
  int size, num, i;
  ivec *where, *r_iv;

  /* Find locations of value. */
  size = ivec_size( iv);
  where = mk_ivec( size);
  num = 0;
  for (i=0; i<size; ++i) {
    if (value == ivec_ref( iv, i)) {
      ivec_set( where, num, i);
      num += 1;
    }
  }

  /* Copy useful part of where into r_iv. */
  r_iv = mk_ivec( num);
  for (i=0; i<num; ++i) ivec_set( r_iv, i, ivec_ref( where, i));
  free_ivec( where);

  return r_iv;
}

/* return the number of elements appearing in both ivecs */
int count_ivec_overlap(ivec *a, ivec *b)
{
  int i;
  int result = 0;

  for (i=0;i<ivec_size(a);i++) if (is_in_ivec(b,ivec_ref(a,i))) result++;

  return result;
}

/* It is fine for iv and r_iv to occupy the same memory
   (see amar.c:sort_ints). */
void ivec_sort(const ivec *iv,ivec *r_iv)
{
  copy_ivec(iv,r_iv);
  sort_ints(r_iv->iarr,ivec_size(r_iv),r_iv->iarr);
}

ivec *mk_ivec_sort(const ivec *iv)
{
  ivec *result = mk_ivec(ivec_size(iv));
  ivec_sort(iv,result);
  return result;
}

/*
 Creates a ivec of indices such that indices[i] is the origional
 location (in the unsorted iv) of the ith smallest value.
 Used when you want the location of the sorted values instead of 
 the sorted vector.
*/
ivec *mk_sorted_ivec_indices(ivec *iv)
{
  ivec* indices;
  int size = ivec_size(iv);
  int i;
  int* iarr_ind;
  int* iarr;

  iarr = mk_iarr_from_ivec(iv);
  iarr_ind = am_malloc_ints(size);

  indices_sort_integers(iarr,size,iarr_ind);
  indices = mk_ivec(size);

  for(i=0;i<size;i++)
  {
    ivec_set(indices,i,iarr_ind[i]);
  }

  am_free_ints(iarr_ind,size);
  am_free_ints(iarr,size);

  return indices;
}

string_array *mk_string_array_from_ivec(ivec *iv)
{
  string_array *sa;
  
  if ( ivec_size(iv) == 0 )
  {
    sa = mk_string_array(1);