amdym.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]

*/


#include<errno.h>
#include "amma.h"
#include "ambs.h"
#include "amiv.h"
#include "amdyv.h"
#include "amdyv_array.h"
#include "am_string_array.h"
#include "amdym.h"


/*private header */
static int is_comment(const char* line);
static int find_number_of_columns(const char* line);
static int find_number_of_rows_and_columns(PFILE *fp , int *nrows, int *ncols);
static int assign_dym_values(PFILE* fp, dym* d);
static int assign_dym_values_for_csv(PFILE* fp, dym *factors, dyv *outputs);
static char *mk_readline(PFILE* fp );

#define MAX_LINE_SIZE (1024*8)
static char* delim = ",";



/** Allocate a matrix nrows X ncolumns */
dym* __stdcall mk_dym(int nrows, int ncolumns)
{
  int i;
  dym *d;
  double **array;

  i = 0;
  d = (dym *) malloc(sizeof(dym));

  array = (double **) malloc(nrows * sizeof(double *));
  array[0] = (double *) malloc(nrows * ncolumns * sizeof(double));

  if( d && array && array[0] ) {
    for(i = 1; i < nrows; i++) array[i] = array[0] + (i * ncolumns);

    d->rows = nrows;
    d->cols = ncolumns;
    d->rows_allocated = nrows;
    d->tdarr = array;

  }
  else my_error( "mk_dym: Failed to allocate memory");

  return d;
}

dym *mk_copy_dym(const dym *d)
{
  return(mk_dym_scalar_mult(d,1.0));
}

void copy_dym(dym *d, dym *r_d)
{
 
  dym_scalar_mult(d,1.0,r_d);
}
 
void __stdcall free_dym(dym* d)
{
  if(d != NULL) {
    if(d->tdarr != NULL) {
      if(d->tdarr[0]) free( d->tdarr[0]);
      free(d->tdarr);
    }
    free( d);
  }
  return;
}

int dym_rows(const dym *d)
{
  return(d->rows);
}

int dym_cols(const dym *d)
{
  return(d->cols);
}



/**
 *
 * 1. open the file to read
 * 2. Find the max rows and columns 
 * 3. Allocate memory to hold the matrix
 * 4. Read the values from the file an assign them to the matrix
 *
 */
dym* mk_read_dym(const char* filename)
{
  int nrows =0 , ncols = 0;
  dym* d = NULL;
  PFILE* fp;

  fp = safe_pfopen(filename,"r");

  if( find_number_of_rows_and_columns(fp,&nrows,&ncols ) != -1 ) {
    d = mk_dym(nrows,ncols);
    if (d == NULL) {
      my_error( "mk_read_dym: unable to allocate memory for dym.");
    }

    prewind(fp); /* start from the beginning of the file */
    
    if( assign_dym_values(fp,d) == -1) {
      free_dym(d); /* error reading, free the allocated memory */
      my_error( "mk_read_dym: Unhandled error while reading.  Freeing dym.");
    }
  }

  pfclose(fp);
  return d;
}

/**
 *
 * 1. open the file to read
 * 2. Find the max rows and columns 
 * 3. Allocate memory to hold the matrix
 * 4. Read the values from the file an assign them to the matrix
 *
 */
void mk_read_dym_for_csv(const char* filename, dym **factors, dyv **outputs)
{
  int nrows =0 , ncols = 0;
  PFILE* fp;

  if (factors == NULL) my_error("mk_read_dym_for_csv: Error: factors is NULL");
  if (outputs == NULL) my_error("mk_read_dym_for_csv: Error: outputs is NULL");
  
  fp = safe_pfopen(filename,"r");

  if( find_number_of_rows_and_columns(fp,&nrows,&ncols ) != -1 ) {
    *factors = mk_dym( nrows, ncols-1);
    if (*factors == NULL) {
      my_error( "mk_read_dym: unable to allocate memory for dym.");
    }
    *outputs = mk_dyv( nrows);

    prewind(fp); /* Jump back to the beginning of the file */
    
    if( assign_dym_values_for_csv( fp, *factors, *outputs) == -1) {
      free_dym( *factors); /* error reading, free the allocated memory */
      free_dyv( *outputs);
      my_error( "mk_read_dym: Unhandled error while reading.  Freeing dym.");
    }
  }

  pfclose(fp);
  return;
}

/** Returns 0 on success , -1 on failure */
int save_dym(PFILE* s, dym* d)
{
  int rc = 0 ;
  
  if( s && d ){
    int i; 
    int j; 

    for(i = 0 ; i < d->rows; i++){
      for(j =0 ; j < d->cols ; j++){
	rc = pfprintf(s,"%g",d->tdarr[i][j]);
	if( rc == -1) return rc;
	if( j == (d->cols - 1))/* last column */
	  rc = pfprintf(s,"\n");
	else
	  rc = pfprintf(s,delim);
	if( rc == -1) return rc;
      }
    }

  }
  pfflush(s);
  return rc;
}

void dym_scalar_mult(const dym *d, double alpha, dym *r_d)
{
  int i,j;
 
  for ( i = 0 ; i < r_d -> rows ; i++ )
    for ( j = 0 ; j < r_d -> cols ; j++ )
      r_d -> tdarr[i][j] = d->tdarr[i][j] * alpha;
}

dym *mk_dym_scalar_mult(const dym *d,double alpha)
{
  dym *result;

  result = mk_dym(d->rows,d->cols);
  dym_scalar_mult(d,alpha,result);
  return(result);
}

void dym_times_dyv(const dym *a, const dyv *b, dyv *result)
{
  int i;
  dyv *temp = mk_dyv(a->rows); 
             /* We need a copy in case b and result share memory */

  if ( a->cols != b -> size )
    my_error("dym_times_dyv: sizes wrong");

  for ( i = 0 ; i < a->rows ; i++ )
  {
    double sum = 0.0;
    int j;
    for ( j = 0 ; j < a->cols ; j++ )
      sum += dym_ref( a, i, j) * dyv_ref( b, j);
    dyv_set( temp, i, sum);
  }

  copy_dyv(temp,result);
  free_dyv(temp);
}

dyv *mk_dym_times_dyv(const dym *a, const dyv *b)
{
  dyv *result = mk_dyv(a->rows);
  dym_times_dyv(a,b,result);
  return(result);
}

void dym_transpose( const dym *d, dym *r_d)
{
  dym *a = mk_dym(d->cols,d->rows);
             /* Note we have to first do the transpose to the result
                a, in case the routine was called with d's memory
                = r_d's memory */
  int i,j;

  for ( i = 0 ; i < d -> rows ; i++ )
    for ( j = 0 ; j < d -> cols ; j++ )
      a->tdarr[j][i] = d->tdarr[i][j];

  copy_dym(a,r_d);
  free_dym(a);
}

dym *mk_dym_transpose( const dym *a)
{
  dym *result = mk_dym(a->cols,a->rows);
  dym_transpose(a,result);
  return(result);
}

void dym_transpose_times_dyv(dym *p,dyv *v, dyv *r_dyv)
{
  int v_size, r_size, i, k;
  double factor, val;
  dyv *temp;

  v_size = dyv_size( v);
  r_size = dyv_size( r_dyv);

  temp = mk_zero_dyv( r_size);

  for (k=0; k<v_size; ++k) {
    factor = dyv_ref( v, k);
    for (i=0; i<r_size; ++i) {
      val = dym_ref( p, k, i);
      dyv_increment( temp, i, factor * val);
    }
  }

  copy_dyv(temp,r_dyv);
  free_dyv(temp);
  return;
}

dyv *mk_dym_transpose_times_dyv(dym *p,dyv *v)
/* Returns p^T v */
{
  dyv *result = mk_dyv(dym_cols(p));
  dym_transpose_times_dyv(p,v,result);
  return(result);
}

/*
  mk_dym_subset - builds a dym out of a subset of rows
  and columns from a matrix x such that if:
    r_given = rows(i)
    c_given = cols(j)
    val = x(r_given,c_given)
  then:
    new_dym(i,j) = val;
  Note: this can also be used to permute rows and columns.

  rows may be NULL meaning "use all rows"
  cols may be NULL meaning "use all columns"
*/
dym *mk_dym_subset( const dym *x, const ivec *rows, const ivec *cols)
{
  int Rn = (rows==NULL) ? dym_rows(x) : ivec_size(rows);
  int Cn = (cols==NULL) ? dym_cols(x) : ivec_size(cols);
  int r, c;
  dym* nu;

  nu = mk_dym(Rn,Cn);

  for(r=0; r<Rn; r++) {
    int row = (rows==NULL) ? r : ivec_ref(rows,r);
    for(c=0;c<Cn;c++) {
      int col = (cols==NULL) ? c : ivec_ref(cols,c);
      dym_set(nu,r,c,dym_ref(x,row,col));
    }
  }

  return nu;
}

dyv *mk_dyv_from_dym_column( const dym *dm, int col)
{
  int numrows, i;
  double val;
  dyv *dv;

  numrows = dym_rows( dm);
  dv = mk_dyv( numrows);

  for (i=0; i<numrows; ++i) {
    val = dym_ref( dm, i, col);
    dyv_set( dv, i, (int) val);
  }

  return dv;
}

dyv *mk_dyv_from_dym_row( const dym *dm, int row)
{
  int numcols;
  dyv *dv;
  numcols = dym_cols( dm);
  dv = mk_dyv_from_farr( dm->tdarr[row], numcols);
  return dv;
}

ivec *mk_ivec_from_binary_dym_column( const dym *dm, int col)
{
  int numrows, i;
  double val;
  ivec *iv;

  numrows = dym_rows( dm);
  iv = mk_ivec( numrows);

  for (i=0; i<numrows; ++i) {
    val = dym_ref( dm, i, col);