src/linalg.c

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2010 by Matthias Ihrke   *
 *   ihrke@nld.ds.mpg.de
 *                                                                         *
 *   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 "linalg.h"
#include "helper.h"
#include <string.h>
#include <math.h>
#include <gsl/gsl_linalg.h>

gsl_matrix* matrix_to_gsl( Array *in, bool alloc ){
  gsl_matrix *out;
  if( in->ndim != 2 ){
     errprintf("in is not a matrix, ndim=%i\n", in->ndim );
  } 
  if( alloc ){
     out = gsl_matrix_alloc( (size_t)in->size[0], (size_t)in->size[1] );
     if( in->dtype==DOUBLE ){ /* faster */
        memcpy( out->data, in->data, in->nbytes );
     } else { /* slow */
        int i,j;
        double el;
        for( i=0; i<in->size[0]; i++ ){
          for( j=0; j<in->size[1]; j++ ){
             array_dtype_to_double( &el, array_index2( in, i, j ), in->dtype );
             gsl_matrix_set( out, i, j, el );
          }
        }
     }
  } else { /* no mallocing */
     MALLOC( out, 1, gsl_matrix );
     out->size1=(size_t)in->size[0];
     out->size2=(size_t)in->size[1];
     out->tda=(size_t)in->size[1];
     out->data=(double*)in->data;
     out->owner=0; /* GSL is not the owner of the data */
     out->block=NULL;
  }
  
  return out;
}

Array* matrix_mean( Array *a, int dim ){
  Array *v;
  int i,j;
  bool ismatrix;
  matrix_CHECK( ismatrix, a );
  if( !ismatrix ) return NULL;
  if( dim!=0 && dim!=1 ){
     errprintf("wrong dimension dim=%i, must be in {0,1}\n", dim );
     return NULL;
  }
  double *p;
  v = array_new2( DOUBLE, 1, a->size[1-dim] );
  for( i=0; i<a->size[1-dim]; i++ ){
     p=&array_INDEX1( v, double, i );
     *p= 0.0;
     for( j=0; j<a->size[dim]; j++ ){
        array_INDEX1( v, double, i ) += 
          array_INDEX2( a, double, (dim==1)?i:j,(dim==1)?j:i );
     }
     array_INDEX1( v, double, i ) /= (double)a->size[dim];
  }
  return v;
}

Array* matrix_get_row( Array *m, int row, bool alloc ){
  Array *rowvec;
  bool ismatrix;
  matrix_CHECK( ismatrix, m );
  if( !ismatrix ) return NULL;  
  if( row<0 || row>=m->size[0] ){
     errprintf("Index out of bounds: %i/%i\n", row, m->size[0] );
     return NULL;
  }

  if( !alloc ){
     rowvec = array_fromptr2( DOUBLE, 1, array_INDEXMEM2( m, row, 0 ), 
                                      m->size[1] );
  } else {
     rowvec = array_new2( DOUBLE, 1, m->size[1] );
     memcpy( rowvec->data, array_INDEXMEM2( m, row, 0 ), 
                m->size[1]*m->dtype_size );
  }
  return rowvec;
}

Array* matrix_get_col( Array *m, int col ){
  Array *colvec;
  int i;
  bool ismatrix;
  matrix_CHECK( ismatrix, m );
  if( !ismatrix ) return NULL;
  if( col<0 || col>=m->size[1] ){
     errprintf("Index out of bounds: %i/%i\n", col, m->size[1] );
     return NULL;
  }
  colvec = array_new2( DOUBLE, 1, m->size[0] );
  for( i=0; i<m->size[0]; i++ ){
     array_INDEX1( colvec, double, i )=array_INDEX2(m, double, i, col );
  }
  return colvec;
}

Array* matrix_mult( const Array *m1, const Array *m2 ){
  Array *out=NULL;
  int i,j,k;
  bool ismatrix;
  matrix_CHECK( ismatrix, m1 ); if( !ismatrix ) return NULL;
  matrix_CHECK( ismatrix, m2 ); if( !ismatrix ) return NULL;
  if( m1->size[1] != m2->size[0] ){
     errprintf("Mismatch in size, cannot calculate matrix product: %i!=%i\n",
                  m1->size[1], m2->size[0] );
     return NULL;
  }

  out = array_new2( DOUBLE, 2, m1->size[0], m2->size[1] );

  double *field;
  for( i=0; i<out->size[0]; i++ ){
     for( j=0; j<out->size[1]; j++ ){
        field=(double*)array_INDEXMEM2( out, i, j );
        *field=0.0;
        for( k=0; k<m1->size[1]; k++ ){
          *field += array_INDEX2( m1, double, i, k )*
             array_INDEX2( m2, double, k, j );
        }
     }
  }
  
  return out;
}
Array* matrix_transpose( Array *m, bool alloc ){
  int i,j;
  bool ismatrix;
  matrix_CHECK( ismatrix, m );
  if( !ismatrix ) return NULL;

  Array *out;
  out = array_new2( DOUBLE, 2, m->size[1], m->size[0] );
  for( i=0; i<m->size[0]; i++ ){
     for( j=0; j<m->size[1]; j++ ){
        array_INDEX2( out, double, j,i) = array_INDEX2( m, double, i,j);
     }
  }
  if( !alloc ){
     memcpy( m->data, out->data, out->nbytes );
     array_free( out );
     out = m;
     double tmp=out->size[0];
     out->size[0]=out->size[1];
     out->size[1]=tmp;
  }
  
  return out;
}

Array* matrix_pca( Array *X, Array **var, bool alloc ){
  Array *out=NULL;
  int i,j;
  bool ismatrix;
  matrix_CHECK( ismatrix, X );
  if( !ismatrix ) return NULL;

  int N,K; /* N observations, K variables */
  N = X->size[0];
  K = X->size[1];

  Array *tmp = array_copy( X, TRUE );

  /* subtract mean from observations */
  Array *mean=matrix_mean( X, 0 ); 
  for( i=0; i<N; i++ ){ 
     for( j=0; j<K; j++ ){
        array_INDEX2( tmp, double, i, j ) -=
          array_INDEX1( mean, double, j );
     }
  }

  array_scale( tmp, 1.0/sqrt((double) N-1 ) );

  gsl_matrix *A=matrix_to_gsl( tmp, TRUE ); /* copy */
  gsl_matrix *V=gsl_matrix_alloc( K, K );
  gsl_vector *S=gsl_vector_alloc( K );
  gsl_vector *workspace=gsl_vector_alloc( K );

  /* A->U, V->V, S->S */
  gsl_linalg_SV_decomp( A, V, S, workspace);
  gsl_matrix_transpose( V );

  if( var ){
     (*var)=array_fromptr2( DOUBLE, 1, S->data, S->size );
     S->owner=0; /* transfer ownership to array */
     (*var)->free_data=1;
     for( i=0; i<array_NUMEL( *var ); i++ ){ 
        array_INDEX1( *var, double, i ) = SQR( array_INDEX1( *var, double, i ) );
     }
  }

  Array *Vp=array_fromptr2( DOUBLE, 2, V->data, V->size1, V->size2 );
  matrix_transpose( tmp, FALSE );
  out=matrix_mult( Vp, tmp ); /* PCA'd data */
  matrix_transpose( out, FALSE );

  if( out->size[0]!=X->size[0] || out->size[1]!=X->size[1] ){
     errprintf("Input/Output dimension mismatch: (%i,%i) vs. (%i,%i)\n",
                  X->size[0], X->size[1], out->size[0], out->size[1] );
  }

  if( !alloc ){ /* write back out->X */
     memcpy( X->data, out->data, out->nbytes );
     array_free( out );
     out = X;
  }

  /* clean up */
  gsl_matrix_free( A );
  gsl_matrix_free( V );
  gsl_vector_free( S );
  gsl_vector_free( workspace );
  array_free( mean );
  array_free( Vp );
  array_free( tmp );

  return out;
}