Changeset 42


Ignore:
Timestamp:
Feb 26, 2004, 4:06:20 PM (18 years ago)
Author:
Jari Häkkinen
Message:

Made a major revision of matrix and vector classes. Everything compiles
but the binaries have not been tested.

Location:
trunk
Files:
26 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/GaussianKernelFunction.cc

    r33 r42  
    99#include "vector.h"
    1010
    11 using namespace thep_cpp_tools;
    12 using namespace std;
     11namespace theplu {
     12namespace cpptools { 
    1313
    1414GaussianKernelFunction::GaussianKernelFunction(double sigma)
     
    1717}
    1818
    19 double GaussianKernelFunction::operator()(const thep_gsl_api::vector a,
    20               const thep_gsl_api::vector b) const   
     19double GaussianKernelFunction::operator()(const gslapi::vector& a,
     20                                          const gslapi::vector& b) const   
    2121{
    2222  u_int nof_valid_datapairs = 0;
     
    2424  // check that a size == b size !
    2525  if(a.size() != b.size())
    26     cerr << "vectors must have same length\n";
     26    std::cerr << "vectors must have same length\n";
    2727  for(u_int k=0; k<a.size(); k++)
    2828    // Check that both a and b are Numbers (NaN is NOT equal to NaN)
    29     if (a.get(k) == a.get(k) & b.get(k) == b.get(k)) {
    30      tmp += pow((a.get(k)-b.get(k)),2);
     29    if (a(k) == a(k) & b(k) == b(k)) {
     30     tmp += (a(k)-b(k))*(a(k)-b(k));
    3131     nof_valid_datapairs++;
    3232    }
    3333  //cout << nof_valid_datapairs << endl;
    3434  if (nof_valid_datapairs == 0)
    35     cerr << "data file error: two columns "
    36    << " have no pair of valid data\n";
     35    std::cerr << "data file error: two columns "
     36              << " have no pair of valid data" << std::endl;
    3737 
    3838  tmp/=nof_valid_datapairs / a.size(); //correcting for missing values
     
    4242}
    4343
     44}} // of namespace cpptools and namespace theplu
  • trunk/src/GaussianKernelFunction.h

    r33 r42  
    11// $Id$
    22
    3 #ifndef cs_cpp_tools_gaussian_kernel_h
    4 #define cs_cpp_tools_gaussian_kernel_h
    5 
    6 
     3#ifndef _theplu_cpptools_gaussian_kernel_function_
     4#define _theplu_cpptools_gaussian_kernel_function_
    75
    86#include "KernelFunction.h"
    9 #include "vector.h"
    107
    11 namespace thep_cpp_tools {
    12  
     8namespace theplu {
     9
     10namespace gslapi {
     11  class vector;
     12}
     13
     14namespace cpptools {
     15
    1316  /**
    1417     Class calculating one element in the kernel matrix using the
     
    3841       
    3942     */
    40     double operator()(const thep_gsl_api::vector,
    41           const thep_gsl_api::vector) const;
     43    double operator()(const gslapi::vector&, const gslapi::vector&) const;
    4244   
    4345  private:
     
    4648  }; // class GaussianKernelFunction
    4749
    48 }; // namespace thep_cpp_tools
     50}} // of namespace cpptools and namespace theplu
    4951
    5052#endif
  • trunk/src/Kernel.cc

    r33 r42  
    11// $Id$
    2 
    3 // System
    4 #include <math.h>
    5 
    62
    73// Thep C++ Tools
     
    117#include "vector.h"
    128
    13 using namespace thep_cpp_tools;
     9namespace theplu {
     10namespace cpptools { 
    1411
    15 Kernel::Kernel(const thep_gsl_api::matrix& data, const KernelFunction& kf)
     12Kernel::Kernel(const gslapi::matrix& data, const KernelFunction& kf)
     13  : k_(data.columns(),data.columns())
    1614{
    17   k_=thep_gsl_api::matrix(data.cols(),data.cols(),0);
    18   for(u_int i=0;i<data.cols();i++)
    19     for(u_int j=0;j<i+1;j++) { 
    20      double tmp = kf(data.col_vector(i),data.col_vector(j));
    21      k_.set(i,j,tmp);
    22     }
     15  for(u_int i=0;i<data.columns();i++)
     16    for(u_int j=0;j<i+1;j++)
     17      k_(i,j)=kf(data.TEMP_col_return(i),data.TEMP_col_return(j));
    2318 
    2419  // Copy lower triangle to upper triangle of Kernel matrix
    25   for(u_int i=0;i<data.cols()-1;i++)
    26     for(u_int j=i+1;j<data.cols();j++)
    27       k_.set(i,j,k_.get(j,i));
     20  for(u_int i=0;i<data.columns()-1;i++)
     21    for(u_int j=i+1;j<data.columns();j++)
     22      k_(i,j)=k_(j,i);
    2823}
    2924
     
    3227
    3328
    34 
    35 
     29}} // of namespace cpptools and namespace theplu
  • trunk/src/Kernel.h

    r26 r42  
    11// $Id$
    22
    3 #ifndef cs_cpp_tools_kernel_h
    4 #define cs_cpp_tools_kernel_h
     3#ifndef _theplu_cpptools_kernel_
     4#define _theplu_cpptools_kernel_
    55
    6 // System includes
    7 #include <sys/types.h>
     6#include "matrix.h"
    87
    9 // should be forward declarations!
    10 #include "matrix.h"
    11 #include "KernelFunction.h"
     8namespace theplu {
     9namespace cpptools {
    1210
    13 namespace thep_cpp_tools {
     11  class KernelFunction;
    1412
    1513  /**
     
    2523       Constructor taking the data matrix and KernelFunction as input.
    2624    */
    27     Kernel(const thep_gsl_api::matrix&, const KernelFunction&);
     25    Kernel(const gslapi::matrix&, const KernelFunction&);
    2826    /**
    2927       Destructor
     
    3331       @return Kernel matrix
    3432    */
    35     const thep_gsl_api::matrix& get() const {return k_;}
     33    const gslapi::matrix& get() const {return k_;}
    3634   
    3735  private:
    38     thep_gsl_api::matrix k_;
     36    gslapi::matrix k_;
    3937
    4038  }; // class Kernel
    4139
    42 }; // namespace thep_cpp_tools
    43 
     40}} // of namespace cpptools and namespace theplu
    4441
    4542#endif
  • trunk/src/KernelFunction.h

    r33 r42  
    11// $Id$
    22
    3 #ifndef cs_cpp_tools_kernel_function_h
    4 #define cs_cpp_tools_kernel_function_h
     3#ifndef _theplu_cpptools_kernel_function_
     4#define _theplu_cpptools_kernel_function_
    55
    6 // System includes
    7 #include <sys/types.h>
    8 #include "vector.h"
     6namespace theplu {
    97
    10 // should be forward declarations!
     8namespace gslapi {
     9  class vector;
     10}
    1111
    12 namespace thep_cpp_tools {
     12namespace cpptools {
     13
    1314  /**
    1415     Virtual Class calculating one element in the kernel matrix (i.e. the
     
    2324       Constructor
    2425    */   
    25     KernelFunction() {};
     26    KernelFunction(void) {};
    2627    /**
    2728       Destructor
     
    2930    virtual ~KernelFunction(void) {};
    3031   
    31     virtual double operator()(const thep_gsl_api::vector,const thep_gsl_api::vector) const = 0;
     32    virtual double operator()(const gslapi::vector&,
     33                              const gslapi::vector&) const = 0;
    3234   
    3335  }; // class KernelFunction
    3436
    35 }; // namespace thep_cpp_tools
    36 
     37}} // of namespace cpptools and namespace theplu
    3738
    3839#endif
  • trunk/src/Makefile.am

    r33 r42  
    55lib_LIBRARIES = libc++_tools.a
    66
    7 libc___tools_a_SOURCES = matrix.cc vector.cc SVD.cc PCA.cc \
    8   random_singleton.cc histogram.cc Kernel.cc\
    9   PolynomialKernelFunction.cc SVM.cc GaussianKernelFunction.cc
     7libc___tools_a_SOURCES = \
     8  GaussianKernelFunction.cc histogram.cc Kernel.cc matrix.cc PCA.cc \
     9  PolynomialKernelFunction.cc random_singleton.cc SVD.cc SVM.cc     \
     10  vector.cc
    1011
    1112INCLUDES = -I/usr/local_bio/include -I/usr/local_bio
  • trunk/src/PCA.cc

    r15 r42  
     1// $Id$
     2
     3#include <iostream>
     4#include <memory>
     5
    16#include "PCA.h"
    2 
    3 
    4 using namespace std;
    5 using namespace thep_cpp_tools;
     7#include "SVD.h"
     8
     9
     10// using namespace std;
     11// using namespace thep_cpp_tools;
     12
     13namespace theplu {
     14namespace cpptools { 
    615
    716PCA::PCA() : process_( false ), explained_calc_( false )
     
    918}
    1019
    11 PCA::PCA( const thep_gsl_api::matrix& A ) : A_( A ),
     20PCA::PCA( const gslapi::matrix& A ) : A_( A ),
    1221              process_( false ),
    1322              explained_calc_( false )
     
    1928  process_ = true;
    2029  // Row-center the data matrix
    21   thep_gsl_api::matrix A_center( A_.rows(), A_.cols() );
     30  gslapi::matrix A_center( A_.rows(), A_.columns() );
    2231  this->row_center( A_center );
    2332 
    2433
    2534  // Single value decompose the data matrix
    26   auto_ptr<thep_cpp_tools::SVD> pSVD( new thep_cpp_tools::SVD( A_center ) );
     35  std::auto_ptr<cpptools::SVD> pSVD( new cpptools::SVD( A_center ) );
    2736  pSVD->process();
    28   thep_gsl_api::matrix U = pSVD->get_u();
    29   thep_gsl_api::matrix V = pSVD->get_v();
     37  gslapi::matrix U = pSVD->get_u();
     38  gslapi::matrix V = pSVD->get_v();
    3039
    3140  // Read the eigenvectors and eigenvalues
     
    3948    }
    4049
    41   eigenvalues_ /= A_center.rows(); 
     50  eigenvalues_ *= 1.0/A_center.rows(); 
    4251 
    4352  // Sort the eigenvectors in order of eigenvalues
     
    5160  if( eigenvalues_[ j ] > eigenvalues_[ i ] )
    5261    {
    53      swap( eigenvalues_[ i ], eigenvalues_[ j ] );
    54      eigenvectors_.swap_rows( i, j );
     62      std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
     63     eigenvectors_.row_swap(i,j);
    5564    }
    5665       }
     
    6372  process_ = true;
    6473  // Row-center the data matrix
    65   thep_gsl_api::matrix A_center( A_.rows(), A_.cols() );
     74  gslapi::matrix A_center( A_.rows(), A_.columns() );
    6675  this->row_center( A_center );
    6776
     
    7281
    7382  // Single value decompose the data matrix
    74   auto_ptr<thep_cpp_tools::SVD> pSVD( new thep_cpp_tools::SVD( A_center ) );
     83  std::auto_ptr<cpptools::SVD> pSVD( new cpptools::SVD( A_center ) );
    7584  pSVD->process();
    76   thep_gsl_api::matrix U = pSVD->get_u();
    77   thep_gsl_api::matrix V = pSVD->get_v();
     85  gslapi::matrix U = pSVD->get_u();
     86  gslapi::matrix V = pSVD->get_v();
    7887
    7988  // Read the eigenvectors and eigenvalues
     
    94103    }
    95104
    96   eigenvalues_ /= A_center.rows(); 
     105  eigenvalues_ *= 1.0/A_center.rows(); 
    97106 
    98107  // Sort the eigenvectors in order of eigenvalues
     
    106115  if( eigenvalues_[ j ] > eigenvalues_[ i ] )
    107116    {
    108      swap( eigenvalues_[ i ], eigenvalues_[ j ] );
    109      eigenvectors_.swap_rows( i, j );
     117     std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
     118     eigenvectors_.row_swap(i,j);
    110119    }
    111120       }
     
    117126// that is, A_ = A_ - M, where M is a matrix
    118127// with the meanvalues of each row
    119 void PCA::row_center( thep_gsl_api::matrix& A_center )
    120 {
    121   meanvalues_ = thep_gsl_api::vector( A_.rows() );
    122   thep_gsl_api::vector A_row_sum = A_.row_sum();
     128void PCA::row_center( gslapi::matrix& A_center )
     129{
     130  meanvalues_ = gslapi::vector( A_.rows() );
     131  gslapi::vector A_row_sum(A_.rows());
     132  for (size_t i=0; i<A_row_sum.size(); ++i)
     133    A_row_sum(i)=A_[i].sum();
    123134  for( size_t i = 0; i < A_center.rows(); ++i )
    124135    {
    125      meanvalues_[ i ] = A_row_sum.get( i ) / static_cast<double>( A_.cols() );
    126      for( size_t j = 0; j < A_center.cols(); ++j )
    127        {
    128   A_center.set(  i, j, A_.get( i , j ) - meanvalues_[ i ] );
    129        }
    130     }
    131 }
    132 
    133 
    134 thep_gsl_api::matrix PCA::projection( const thep_gsl_api::matrix&
     136     meanvalues_[i] = A_row_sum(i) / static_cast<double>(A_.columns());
     137     for( size_t j = 0; j < A_center.columns(); ++j )
     138       A_center(i,j) = A_(i,j) - meanvalues_(i);
     139    }
     140}
     141
     142
     143gslapi::matrix PCA::projection( const gslapi::matrix&
    135144              samples ) const
    136145{
    137   assert( process_ );
    138   const size_t Ncol = samples.cols();
     146  const size_t Ncol = samples.columns();
    139147  const size_t Nrow = samples.rows();
    140   thep_gsl_api::matrix projs( Ncol, Ncol );
    141  
     148  gslapi::matrix projs( Ncol, Ncol );
     149 
     150  gslapi::vector temp(samples.rows());
    142151  for( size_t j = 0; j < Ncol; ++j )
    143152    {
    144      thep_gsl_api::matrix temp = samples.col( j );   
    145      thep_gsl_api::vector centered( Nrow );
     153      for (size_t i=0; i<Ncol; ++i )
     154        temp(i) = samples(i,j);   
     155     gslapi::vector centered( Nrow );
    146156
    147157     for( size_t i = 0; i < Nrow; ++i )
    148        {
    149   centered[ i ] = temp.get( i, 0 ) - meanvalues_[ i ];
    150   //cout << temp.get( i, 0 ) << " " << meanvalues_[ i ] << " " << centered[ i ] << endl;
    151        }
    152 
    153      thep_gsl_api::vector proj( eigenvectors_ * centered );
     158       centered(i) = temp(i) - meanvalues_(i);
     159
     160     gslapi::vector proj( eigenvectors_ * centered );
    154161      for( size_t i = 0; i < Ncol; ++i )
    155        {
    156   projs.set( i, j,  proj[ i ] );
    157        }
     162        projs(i,j)=proj(i);
    158163    } 
    159164  return projs;
     
    161166
    162167
    163 thep_gsl_api::matrix PCA::projection_transposed( const thep_gsl_api::matrix&
     168gslapi::matrix PCA::projection_transposed( const gslapi::matrix&
    164169             samples ) const
    165170{
    166   assert( process_ );
    167   const size_t Ncol = samples.cols();
     171  const size_t Ncol = samples.columns();
    168172  const size_t Nrow = samples.rows();
    169   thep_gsl_api::matrix projs( Nrow, Ncol );
    170  
     173  gslapi::matrix projs( Nrow, Ncol );
     174 
     175  gslapi::vector temp(samples.rows());
    171176  for( size_t j = 0; j < Ncol; ++j )
    172177    {
    173      thep_gsl_api::matrix temp = samples.col( j );   
    174      thep_gsl_api::vector centered( Nrow );
     178      for (size_t i=0; i<Ncol; ++i )
     179        temp(i) = samples(i,j);   
     180     gslapi::vector centered( Nrow );
    175181
    176182     for( size_t i = 0; i < Nrow; ++i )
    177        {
    178   centered[ i ] = temp.get( i, 0 ) - meanvalues_[ i ];
    179   //cout << temp.get( i, 0 ) << " " << meanvalues_[ i ] << " " << centered[ i ] << endl;
    180        }
    181 
    182      thep_gsl_api::vector proj( eigenvectors_ * centered );
     183       centered(i)=temp(i)-meanvalues_(i);
     184
     185     gslapi::vector proj( eigenvectors_ * centered );
    183186     for( size_t i = 0; i < Nrow; ++i )
    184        {
    185     projs.set( i, j,  proj[ i ] );
    186        }
     187       projs(i,j)=proj(i);
    187188    }
    188189  return projs;
     
    194195{
    195196  size_t DIM = eigenvalues_.size();
    196   explained_intensity_ = thep_gsl_api::vector( DIM );
     197  explained_intensity_ = gslapi::vector( DIM );
    197198  double sum = 0;
    198199  for( size_t i = 0; i < DIM; ++i )
     
    221222bool PCA::test()
    222223{
    223   cout << "1. Print original matrix A = \n";//Testing row-centering of matrix\n";
    224   thep_gsl_api::matrix A( 3, 4 );
     224  std::cout << "1. Print original matrix A = \n";
     225  gslapi::matrix A( 3, 4 );
    225226  for( size_t i = 0; i < 3; ++i )
    226     {
    227      for( size_t j = 0; j < 4; ++j )
    228        {
    229   A.set( i, j, sin( static_cast<double>( i + j + 2 + (i+1)*(j+1) ) ) );
    230        }
    231     }
    232   A_ = A;
     227    for( size_t j = 0; j < 4; ++j )
     228      A(i,j)= sin( static_cast<double>(i+j+2+(i+1)*(j+1)) );
     229  A_=A;
    233230
    234231// Print only matrix that is row-centered as PCA ( process() )
    235232// will do that before calculation
    236233  this->row_center( A );
    237   cout << A << endl;
    238 
    239   cout << "2. Testing PCA\n";
     234  std::cout << A << std::endl;
     235
     236  std::cout << "2. Testing PCA\n";
    240237  this->process_transposed_problem();
    241238
    242239  for( size_t i = 0; i < eigenvalues_.size(); ++i )
    243240    {
    244      cout << "lambda(" << i << ") = " << eigenvalues_[i] << " "
    245     << "has eigenvector = " << endl;
    246      cout << eigenvectors_.row( i ) << endl;
     241     std::cout << "lambda(" << i << ") = " << eigenvalues_[i] << " "
     242    << "has eigenvector = " << std::endl;
     243     std::cout << eigenvectors_[i] << std::endl;
    247244    }
    248245 
     
    308305//     }
    309306   
    310   cout << "\nThe following projection is obtained:\n";
    311   cout << this->projection_transposed( A_ ).transpose() << endl;
     307  std::cout << "\nThe following projection is obtained:\n";
     308  std::cout << this->projection_transposed( A_ ).transpose() << std::endl;
    312309 
    313310
    314311  return true;
    315312}
     313
     314}} // of namespace cpptools and namespace theplu
  • trunk/src/PCA.h

    r16 r42  
    1 #ifndef GENETICS_PCA_ANALYZER_H
    2 #define GENETICS_PCA_ANALYZER_H
     1// $Id$
    32
    4 // C++ tools include
    5 /////////////////////
     3#ifndef _theplu_cpptools_pca_
     4#define _theplu_cpptools_pca_
     5
    66#include "vector.h"
    77#include "matrix.h"
    8 #include "SVD.h"
    98
    109// Standard C++ includes
    1110////////////////////////
    12 #include <vector>
    13 #include <iostream>
    14 #include <memory>
    15 #include <cstdlib>
     11// #include <vector>
     12// #include <iostream>
     13// #include <memory>
     14// #include <cstdlib>
    1615
    1716
    18 namespace thep_cpp_tools
    19 {
     17namespace theplu {
     18namespace cpptools {
     19
    2020  /**
    2121     Class performing PCA using SVD. This class assumes that
     
    4242       should have been performed and no products.
    4343     */
    44     explicit PCA( const thep_gsl_api::matrix& );
     44    explicit PCA( const gslapi::matrix& );
    4545
    4646   
     
    7474       Returns eigenvector \a i
    7575    */
    76     thep_gsl_api::matrix get_eigenvector( const size_t& i ) const
     76    // Jari, change this to
     77//     const gslapi::vector& get_eigenvector( const size_t& i ) const
     78    const gslapi::vector get_eigenvector( const size_t& i ) const
    7779    {
    78       return eigenvectors_.row( i );
     80      return eigenvectors_[i];
    7981    }
    8082
     
    104106  spanned by the eigenvectors.
    105107    */
    106     thep_gsl_api::matrix projection( const thep_gsl_api::matrix& ) const;
     108    gslapi::matrix projection( const gslapi::matrix& ) const;
    107109   
    108110    /**
    109111  Same as projection() but works when used process_transposed_problem()
    110112    */
    111     thep_gsl_api::matrix projection_transposed( const thep_gsl_api::matrix& ) const;
     113    gslapi::matrix projection_transposed( const gslapi::matrix& ) const;
    112114
    113115   
    114116   
    115117  private:
    116     thep_gsl_api::matrix A_;
    117     thep_gsl_api::matrix  eigenvectors_;
    118     thep_gsl_api::vector  eigenvalues_;
    119     thep_gsl_api::vector  explained_intensity_;
    120     thep_gsl_api::vector  meanvalues_;
     118    gslapi::matrix A_;
     119    gslapi::matrix  eigenvectors_;
     120    gslapi::vector  eigenvalues_;
     121    gslapi::vector  explained_intensity_;
     122    gslapi::vector  meanvalues_;
    121123    bool process_, explained_calc_;
    122124   
     
    126128       with the meanvalues of each row
    127129    */
    128     void row_center( thep_gsl_api::matrix& A_center );
     130    void row_center( gslapi::matrix& A_center );
    129131
    130132    /**
     
    137139  }; // class PCA 
    138140 
    139 }
     141}} // of namespace cpptools and namespace theplu
    140142
    141143#endif
  • trunk/src/PolynomialKernelFunction.cc

    r30 r42  
    44#include <math.h>
    55
    6 
    76// Thep C++ Tools
    87#include "PolynomialKernelFunction.h"
    98#include "vector.h"
    109
    11 using namespace thep_cpp_tools;
    12 using namespace std;
     10
     11namespace theplu {
     12namespace cpptools { 
    1313
    1414PolynomialKernelFunction::PolynomialKernelFunction(int order)
     
    1717}
    1818
    19 double PolynomialKernelFunction::operator()(const thep_gsl_api::vector a,
    20               const thep_gsl_api::vector b) const   
     19double PolynomialKernelFunction::operator()(const gslapi::vector& a,
     20                                            const gslapi::vector& b) const   
    2121{
     22  using namespace std;
     23
    2224  u_int nof_valid_datapairs = 0;
    2325  double tmp = 0;
     
    2729  for(u_int k=0; k<a.size(); k++)
    2830    // Check that both a and b are Numbers (NaN is NOT equal to NaN)
    29     if (a.get(k) == a.get(k) & b.get(k) == b.get(k)) {
    30      tmp += a.get(k)*b.get(k);
     31    if (a(k) == a(k) & b(k) == b(k)) {
     32     tmp += a(k)*b(k);
    3133     nof_valid_datapairs++;
    3234    }
     
    4345}
    4446
     47}} // of namespace cpptools and namespace theplu
  • trunk/src/PolynomialKernelFunction.h

    r30 r42  
    11// $Id$
    22
    3 #ifndef cs_cpp_tools_polynomial_kernel_h
    4 #define cs_cpp_tools_polynomial_kernel_h
    5 
    6 
     3#ifndef _theplu_cpptools_polynomial_kernel_function_
     4#define _theplu_cpptools_polynomial_kernel_function_
    75
    86#include "KernelFunction.h"
    9 #include "vector.h"
    107
    11 namespace thep_cpp_tools {
    12  
     8namespace theplu {
     9
     10namespace gslapi {
     11  class vector;
     12}
     13
     14namespace cpptools {
     15
    1316  /**
    1417     Class calculating one element in the kernel matrix using the
     
    3841       
    3942     */
    40     double operator()(const thep_gsl_api::vector,
    41           const thep_gsl_api::vector) const;
     43    double operator()(const gslapi::vector&, const gslapi::vector&) const;
    4244   
    4345  private:
     
    4648  }; // class PolynomialKernelFunction
    4749
    48 }; // namespace thep_cpp_tools
     50}} // of namespace cpptools and namespace theplu
    4951
    5052#endif
  • trunk/src/SVD.cc

    r11 r42  
    11// $Id$
    22
     3#include <cstdlib>
     4#include <cmath>
     5#include <iostream>
     6
     7#include <gsl/gsl_linalg.h>
     8
    39#include "SVD.h"
    4 using namespace thep_cpp_tools;
     10
     11
     12namespace theplu {
     13namespace cpptools { 
    514
    615//Constructor that can only be used for test-method
     
    918
    1019//Constructor initializing the svd object with matrix A to be decomposed
    11 SVD::SVD( const thep_gsl_api::matrix& Ain ) : A_( Ain )
     20SVD::SVD( const gslapi::matrix& Ain ) : A_( Ain )
    1221{
    1322  //Get dimensions
    1423  ////////////////
    15   size_t n = Ain.cols();
     24  size_t n = Ain.columns();
    1625
    1726  //Instantiate matrix with correct dimenstions
    1827  //and assign all elements the value zero
    1928  /////////////////////////////////////////////
    20   V_ = thep_gsl_api::matrix( n, n, true );
     29  V_ = gslapi::matrix( n, n, true );
    2130 
    2231  //Assign vector s_ with correct dimenstion
    2332  /////////////////////////////////////////////////
    24   s_ = thep_gsl_api::vector( n );
     33  s_ = gslapi::vector( n );
    2534 
    2635  loaded_ = true;
     
    3342//Ax = b
    3443//The x vector will be reached using the get_x() function.
    35 SVD::SVD( const thep_gsl_api::matrix& Ain, const thep_gsl_api::vector& b ) : A_( Ain ), b_( b )
     44SVD::SVD( const gslapi::matrix& Ain, const gslapi::vector& b ) : A_( Ain ), b_( b )
    3645{
    3746  //Get dimensions
    3847  ////////////////
    39   size_t n = Ain.cols();
     48  size_t n = Ain.columns();
    4049 
    4150  //Instantiate matrix with correct dimenstions
    4251  //and assign all elements the value zero
    4352  /////////////////////////////////////////////
    44     V_ = thep_gsl_api::matrix( n, n, true );
     53    V_ = gslapi::matrix( n, n, true );
    4554
    4655
    4756  //Assign vector s_ with correct dimenstion
    4857  /////////////////////////////////////////////////
    49   s_ = thep_gsl_api::vector( n );
     58  s_ = gslapi::vector( n );
    5059 
    5160 
    5261  //An extra vector holding the answer to equation Ax = b is required
    5362  ///////////////////////////////////////////////////////////////////
    54   x_ = thep_gsl_api::vector( n );
     63  x_ = gslapi::vector( n );
    5564
    5665
     
    6877
    6978
     79
     80double SVD::norm(gslapi::matrix& A) const
     81{
     82  double sum=0.0;
     83  for (size_t i=0; i<A.rows(); ++i)
     84    for (size_t j=0; j<A.columns(); ++j)
     85      sum += A(i,j)*A(i,j);
     86  return sqrt(sum);
     87}
     88
     89
     90
    7091//Method to process calculation
    7192bool SVD::process( SVDtypes stype )
     
    113134  //required here
    114135  ///////////////////////////////
    115   thep_gsl_api::vector w( A_.cols() );
    116   if( gsl_linalg_SV_decomp( A_.get_gsl_matrix(), V_.get_gsl_matrix(),
    117           s_.get_gsl_vector(), w.get_gsl_vector()
     136  gslapi::vector w( A_.columns() );
     137  if( gsl_linalg_SV_decomp( A_.gsl_matrix_pointer(), V_.gsl_matrix_pointer(),
     138          s_.gsl_vector_pointer(), w.gsl_vector_pointer()
    118139          ) != 0 )
    119140    {
     
    135156  ///////////////////////////////
    136157  size_t n = A_.rows();
    137   thep_gsl_api::vector w( n );
    138   thep_gsl_api::matrix X( n, n, true );
    139 
    140   if( gsl_linalg_SV_decomp_mod( A_.get_gsl_matrix(), X.get_gsl_matrix(),
    141         V_.get_gsl_matrix(), s_.get_gsl_vector(),
    142         w.get_gsl_vector()
     158  gslapi::vector w( n );
     159  gslapi::matrix X( n, n, true );
     160
     161  if( gsl_linalg_SV_decomp_mod( A_.gsl_matrix_pointer(), X.gsl_matrix_pointer(),
     162        V_.gsl_matrix_pointer(), s_.gsl_vector_pointer(),
     163        w.gsl_vector_pointer()
    143164        ) != 0 )
    144165    {
     
    158179  //Perform SVD using GSL library
    159180  ///////////////////////////////
    160   if( gsl_linalg_SV_decomp_jacobi( A_.get_gsl_matrix(), V_.get_gsl_matrix(),
    161            s_.get_gsl_vector()
     181  if( gsl_linalg_SV_decomp_jacobi( A_.gsl_matrix_pointer(), V_.gsl_matrix_pointer(),
     182           s_.gsl_vector_pointer()
    162183           ) != 0 )
    163184    {
     
    179200 
    180201  //Solve ...
    181   if( gsl_linalg_SV_solve( A_.get_gsl_matrix(), V_.get_gsl_matrix(),
    182          s_.get_gsl_vector(), b_.get_gsl_vector(),
    183          x_.get_gsl_vector()
     202  if( gsl_linalg_SV_solve( A_.gsl_matrix_pointer(), V_.gsl_matrix_pointer(),
     203         s_.gsl_vector_pointer(), b_.gsl_vector_pointer(),
     204         x_.gsl_vector_pointer()
    184205         ) != 0 )
    185206    {
     
    195216
    196217//This method will create a diagonal matrix of vector s_
    197 thep_gsl_api::matrix SVD::get_s() const
     218gslapi::matrix SVD::get_s() const
    198219{
    199220  size_t n = A_.rows();
    200221 
    201   thep_gsl_api::matrix S( n, n, true );
     222  gslapi::matrix S( n, n, true );
    202223
    203224  for( size_t i = 0; i < n; ++i )
    204     {
    205      S.set( i, i, s_.get( i ) );
    206     }
     225     S(i,i)= s_(i);
    207226
    208227  return S;
     
    225244  ///////////////////////////////
    226245    const size_t DIM = 3;
    227     thep_gsl_api::matrix A( DIM, DIM );
     246    gslapi::matrix A( DIM, DIM );
    228247    for( size_t i = 0; i < DIM; ++i )
    229248      {
    230249       for( size_t j = 0; j < DIM; ++j )
    231250   {
    232     A.set( i, j, 10 * (double) rand() / RAND_MAX ); //Use random numbers between 0 and 10
     251     // Jari, change below random numbers to proper RNG
     252    A(i,j)= 10 * (double) rand() / RAND_MAX; //Use random numbers between 0 and 10
    233253   }
    234254      }
     
    239259    A_ = A;
    240260    size_t n = A_.rows();
    241     size_t m = A_.cols();
     261    size_t m = A_.columns();
    242262    n = m = DIM;
    243263    loaded_ = true;
     
    246266    //Set work matrices
    247267    ///////////////////
    248     V_ = thep_gsl_api::matrix( n, n, true );
     268    V_ = gslapi::matrix( n, n, true );
    249269
    250270    //Assign vector s_ with correct dimensions
    251271    /////////////////////////////////////////////////
    252     s_ = thep_gsl_api::vector( n );
     272    s_ = gslapi::vector( n );
    253273
    254274    double error = 1;
    255     thep_gsl_api::matrix E, S;
     275    gslapi::matrix E, S;
    256276
    257277
     
    272292   
    273293
    274     error = sqrt( E.norm( 2 ) );
     294    error = sqrt(norm(E));
    275295    cout << "error(Unmodified method) = " << error << endl;
    276296    cout << "to compare with " << MAXTOL << endl;
     
    285305    S = this->get_s();
    286306    E = A - A_*S*V_.transpose();
    287     error = sqrt( E.norm( 2 ) );
     307    error = sqrt(norm(E));
    288308    cout << "error(Modified method) = " << error << endl;
    289309    cout << "to compare with " << MAXTOL << endl;
     
    298318    S = this->get_s();
    299319    E = A - A_*S*V_.transpose();
    300     error = sqrt( E.norm( 2 ) );
     320    error = sqrt(norm(E));
    301321    cout << "error(Jacobi method) = " << error << endl;
    302322    cout << "to compare with " << MAXTOL << endl;
     
    307327}
    308328
    309 
     329}} // of namespace cpptools and namespace theplu
  • trunk/src/SVD.h

    r11 r42  
    11// $Id$
    22
    3 #ifndef _THEP_CPPTOOLS_SVD_
    4 #define _THEP_CPPTOOLS_SVD_
     3#ifndef _theplu_cpptools_svd_
     4#define _theplu_cpptools_svd_
    55
    6 // C++ tools include
    7 /////////////////////
    86#include "matrix.h"
    97#include "vector.h"
    10 #include <gsl/gsl_vector.h>
    118
    12 
    13 // Standard includes
    14 ////////////////////
    15 #include <iostream>
    16 #include <cmath>   
    17 #include <cstdlib> //To include srand() and rand()
    18 #include <ctime>
    19 
    20 
    21 namespace thep_cpp_tools
    22 
     9namespace theplu {
     10namespace cpptools { 
    2311
    2412  /**
     
    3018     V = Orthogonal matrix, size NxN\n
    3119     A = U S V' = (MxN)(NxN)(NxN) = (MxN)\n   
     20
     21     Shouldn't there be a reference to the algortithm here?
     22
    3223     @version 1.0 $Revision$ $Date$
    3324  */
    3425
    3526  class SVD
    36   {   
    37     static const double MAXTOL = 1E-6;
     27  {
     28    static const double MAXTOL = 1E-6;
    3829
    39     enum SVDtypes
    40       {
    41   Unmodified,
    42   Modified,
    43   Jacobi,
    44   Solver
    45       };
     30    enum SVDtypes {
     31      Unmodified,
     32      Modified,
     33      Jacobi,
     34      Solver
     35    };
    4636
    4737  public:
     
    4939       The default constructor should only be used for testing!!!
    5040    */
    51     SVD();
     41    SVD(void);
    5242
    5343   
     
    5747       object.
    5848    */
    59     SVD( const thep_gsl_api::matrix& );
     49    SVD( const gslapi::matrix& );
    6050   
    6151    /**
     
    6757    */
    6858
    69     SVD( const thep_gsl_api::matrix&, const thep_gsl_api::vector& b);
     59    SVD( const gslapi::matrix&, const gslapi::vector& b);
    7060    ~SVD();
    7161
     
    8676       executed.
    8777    */
    88     thep_gsl_api::matrix get_u() const { return A_; }
     78    gslapi::matrix get_u() const { return A_; }
    8979
    9080    /**
     
    9282       executed.
    9383    */
    94     thep_gsl_api::matrix get_v() const { return V_; }
     84    gslapi::matrix get_v() const { return V_; }
    9585
    9686
     
    10292    */
    10393
    104     thep_gsl_api::vector get_x() const { return x_; }
     94    gslapi::vector get_x() const { return x_; }
    10595
    10696    /**
     
    10898       executed.
    10999    */
    110     thep_gsl_api::matrix get_s() const;
     100    gslapi::matrix get_s() const;
    111101
    112102
     
    114104       Function returning diagonal matrix S in vector form.
    115105    */
    116     thep_gsl_api::vector get_s_vec() const { return s_; }
     106    gslapi::vector get_s_vec() const { return s_; }
    117107
    118108
     
    128118  private:
    129119
    130     thep_gsl_api::matrix A_, V_;   
    131     thep_gsl_api::vector s_, b_, x_;
     120    gslapi::matrix A_, V_;   
     121    gslapi::vector s_, b_, x_;
    132122    bool loaded_, solver_;
    133123
     
    137127    bool process_jacobi();
    138128    bool process_solver();
     129
     130    double norm(gslapi::matrix& A) const; // move into class gslapi::matrix?
    139131  }; 
    140 };
     132
     133}} // of namespace cpptools and namespace theplu
    141134
    142135#endif
  • trunk/src/SVM.cc

    r37 r42  
    1313#include "random_singleton.h"
    1414
    15 using namespace thep_cpp_tools;
    16 using namespace std;
     15namespace theplu {
     16namespace cpptools { 
    1717
    18 SVM::SVM(const thep_gsl_api::matrix& kernel,
    19    const thep_gsl_api::vector& target) : trained_(false),
    20                  kernel_(kernel),
    21                  target_(target),
    22                  alpha_(target.size(),true,true),
    23                  bias_(0)
    24  
     18SVM::SVM(const gslapi::matrix& kernel, const gslapi::vector& target)
     19  : trained_(false),
     20    kernel_(kernel),
     21    target_(target),
     22    alpha_(target.size()),
     23    bias_(0)
    2524{
    2625}
    2726
    28 void SVM::train() //Should be done so one can choose to not train on all the samples
     27void SVM::train(void) //Should be done so one can choose to not train on all the samples
    2928{
    30   thep_gsl_api::vector E = thep_gsl_api::vector(-target_);
     29  using namespace std;
     30  gslapi::vector E = gslapi::vector(-target_);
    3131 
    3232  double upper_bound = pow(10., 32);
     
    3636  double u;
    3737  double v;
    38   thep_gsl_api::vector dalpha; 
     38  gslapi::vector dalpha; 
    3939 
    4040  // Stop criteria 
     
    4747       
    4848     // Choosing a pair of variables to modify (Should be done more clever)
    49      u_long index1 = rnd->get_uniform_int(kernel_.cols());
    50      u_long index2 = rnd->get_uniform_int(kernel_.cols()-1);
     49     u_long index1 = rnd->get_uniform_int(kernel_.columns());
     50     u_long index2 = rnd->get_uniform_int(kernel_.columns()-1);
    5151     if (index2 >= index1)
    5252      index2++;
     
    8080       }
    8181
    82      double k = kernel_.get(index1, index1) + kernel_.get(index2, index2) -
    83        2*kernel_.get(index1, index2);
     82     double k = ( kernel_(index1, index1) + kernel_(index2, index2) -
     83                  2*kernel_(index1, index2) );
    8484     alpha_new = alpha_[index2] + target_[index2]*
    8585       (E[index1]-E[index2])/k;
     
    9393     stop_condition = stop( target_, kernel_, alpha_);
    9494     if (count>10000000){
    95       cerr << "SVM): " << "more than 10,000,000 epochs reached" << endl;
     95      cerr << "SVM): " << "more than 10,000,000 epochs reached" << endl;
    9696       exit(1);
    9797     }
     
    101101  double min_output_positive = 10000;
    102102  double max_output_negative = -10000;
    103   thep_gsl_api::vector output_unbiased = kernel_ * alpha_.mul_elements(target_);
     103  gslapi::vector output_unbiased = kernel_ * alpha_.mul_elements(target_);
    104104  for (u_int i=0; i<target_.size(); i++){
    105105   if (target_[i]==1){
     
    120120}
    121121
    122 bool SVM::stop( const thep_gsl_api::vector& target_,
    123     const thep_gsl_api::matrix& kernel_,
    124     const thep_gsl_api::vector& alpha_)
     122bool SVM::stop( const gslapi::vector& target_,
     123    const gslapi::matrix& kernel_,
     124    const gslapi::vector& alpha_)
    125125{
    126126  double min_output_positive = 10000;
    127127  double max_output_negative = -10000;
    128128  double epsilon = 0.000001; // used twice, should perhaps be two different values
    129   thep_gsl_api::vector output_unbiased = kernel_ * alpha_.mul_elements(target_);
     129  gslapi::vector output_unbiased = kernel_ * alpha_.mul_elements(target_);
    130130  for (u_int i=0; i<target_.size(); i++){
    131131   if (target_[i]==1){
     
    154154}
    155155
    156 
     156}} // of namespace cpptools and namespace theplu
    157157
    158158#endif
  • trunk/src/SVM.h

    r37 r42  
    1414
    1515
    16 namespace thep_cpp_tools
    17 {
     16namespace theplu {
     17namespace cpptools { 
    1818 
    1919  class SVM
     
    2424       Constructor taking the kernel matrix and the target vector as input
    2525    */
    26     SVM(const thep_gsl_api::matrix&, const thep_gsl_api::vector&);
     26    SVM(const gslapi::matrix&, const gslapi::vector&);
    2727           
    2828    /**
     
    3434       Function will return \f$\alpha\f$
    3535    */
    36     inline thep_gsl_api::vector get_alpha() const;
     36    inline gslapi::vector get_alpha() const;
    3737
    3838    /**
    3939       Function will return the output from SVM
    4040    */
    41     inline thep_gsl_api::vector get_output() const;
     41    inline gslapi::vector get_output() const;
    4242
    4343   
     
    4545  private:
    4646    bool trained_;
    47     thep_gsl_api::matrix kernel_;
    48     thep_gsl_api::vector target_;
    49     thep_gsl_api::vector alpha_;
     47    gslapi::matrix kernel_;
     48    gslapi::vector target_;
     49    gslapi::vector alpha_;
    5050    double bias_;
    5151
     
    5858
    5959    */
    60     bool SVM::stop(const thep_gsl_api::vector& target_,
    61        const thep_gsl_api::matrix& kernel_,
    62        const thep_gsl_api::vector& alpha_);
     60    bool SVM::stop(const gslapi::vector& target_,
     61       const gslapi::matrix& kernel_,
     62       const gslapi::vector& alpha_);
    6363  };
    6464 
    6565    // class SVM
    66   thep_gsl_api::vector SVM::get_alpha() const 
     66  gslapi::vector SVM::get_alpha() const 
    6767  {
    6868    return alpha_;
    6969  }
    7070
    71   thep_gsl_api::vector SVM::get_output() const 
     71  gslapi::vector SVM::get_output() const 
    7272  {
    73     thep_gsl_api::vector bias(target_.size(), false, false);
     73    gslapi::vector bias(target_.size());
    7474    bias.set_all(bias_);
    7575    return kernel_ * alpha_.mul_elements(target_) + bias;
    7676   
    7777  }
    78 }; // namespace thep_c++_tools
     78
     79}} // of namespace cpptools and namespace theplu
    7980
    8081#endif
  • trunk/src/matrix.cc

    r30 r42  
     1// $Id$
    12
    23#include <iostream>
     4#include <sstream>
    35#include <string>
    4 #include <sstream>
    56#include <vector>
    67
     8#include <gsl/gsl_blas.h>
     9#include <gsl/gsl_linalg.h>
     10
    711#include "matrix.h"
    8 
    9 using namespace thep_gsl_api;
    10 
    11 
    12 // Constructors and Destructors
    13 ///////////////////////////////
    14 
    15 matrix::matrix() : m_( NULL )
    16 {
    17 }
    18 
    19 
    20 matrix::matrix( const size_t& rows, const size_t& cols,
    21     bool init_to_zero )
    22 {
    23   if( init_to_zero )
    24     {
    25      m_ = gsl_matrix_calloc( rows, cols );
     12#include "vector.h"
     13
     14
     15namespace theplu {
     16namespace gslapi {
     17
     18
     19
     20  matrix::matrix(void)
     21    : m_(NULL)
     22  {
     23  }
     24
     25
     26
     27  matrix::matrix(const size_t& r, const size_t& c, double init_value)
     28  {
     29    m_ = gsl_matrix_alloc(r,c);
     30    set_all(init_value);
     31  }
     32
     33
     34
     35  matrix::matrix(const matrix& other)
     36  {
     37    gsl_matrix_free(m_);
     38    m_ = other.gsl_matrix_copy();
     39  }
     40
     41
     42
     43  // Constructor that gets data from istream
     44  matrix::matrix(std::istream& is)
     45  {
     46    using namespace std;
     47
     48    // read the data file and store in stl vectors (dynamically
     49    // expandable)
     50    std::vector<std::vector<double> > data_matrix;
     51    u_int nof_columns=0;
     52    u_int nof_rows = 0;
     53    string s;
     54    for (nof_rows = 0; getline(is, s, '\n'); nof_rows++) {
     55      istringstream line(s);
     56      std::vector<double> v;
     57      string tmp_string;
     58      while (line >> tmp_string) {
     59        if(!is.good()) {
     60          cerr << "matrix::matrix(std::istream&): "
     61               << "error reading data file!" << endl;
     62          exit(1);
     63        }
     64        double t=atof(tmp_string.c_str());
     65        // Here we should check that it actually was a double!!!!!
     66        v.push_back(t);
     67      }
     68      if(nof_columns==0)
     69        nof_columns=v.size();
     70      else if(v.size()!=nof_columns) {
     71        cerr << "matrix::matrix(std::istream&) data file error: "
     72             << "line" << nof_rows+1 << " has " << v.size()
     73             << " columns; expected " << nof_columns
     74             << " columns"
     75             << endl;
     76        exit(1);
     77      }
     78      data_matrix.push_back(v);
     79    }
     80    // manipulate the state of the stream to be good
     81    is.clear(std::ios::goodbit);
     82    // convert the data to a gsl matrix
     83    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
     84    for(u_int i=0;i<nof_rows;i++)
     85      for(u_int j=0;j<nof_columns;j++)
     86        gsl_matrix_set( m_, i, j, data_matrix[i][j] );
     87  }
     88
     89
     90
     91  matrix::~matrix(void)
     92  {
     93    if (m_) {
     94      gsl_matrix_free(m_);
     95      m_=NULL;
     96    }
     97  }
     98
     99
     100
     101  gsl_matrix* matrix::gsl_matrix_copy(void) const
     102  {
     103    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
     104    gsl_matrix_memcpy(m,m_);
     105    return m;
     106  }
     107
     108
     109
     110  matrix matrix::mul_elements(const matrix& other)
     111  {
     112    matrix res(other);
     113    gsl_matrix_mul_elements( res.m_, other.m_ );
     114    return res;
     115  }
     116
     117
     118
     119  vector matrix::TEMP_col_return(size_t column) const
     120  {
     121    vector vec(rows());
     122    for (size_t i=0; i<rows(); ++i)
     123      vec(i)=operator()(i,column);
     124    return vec;
     125  }
     126
     127
     128
     129  vector matrix::operator[](size_t row) const
     130  {
     131    vector vec(columns());
     132    for (size_t i=0; i<columns(); ++i)
     133      vec(i)=operator()(row,i);
     134    return vec;
     135  }
     136
     137
     138
     139  matrix matrix::operator*( const matrix &other ) const
     140  {
     141    matrix res(rows(), other.columns());
     142    gsl_linalg_matmult(m_,other.m_,res.m_);
     143    return res;
     144  }
     145
     146
     147
     148  vector matrix::operator*( const vector &other ) const
     149  {
     150    gsl_vector* gslvec=gsl_vector_alloc(rows());
     151    gsl_blas_dgemv(CblasNoTrans, 1.0, m_, other.gsl_vector_pointer(), 0.0,
     152                   gslvec );
     153    vector res(gslvec);
     154    return res;
     155  }
     156
     157
     158
     159  matrix matrix::operator+( const matrix &other ) const
     160  {
     161    matrix res( *this );
     162    gsl_matrix_add( res.m_, other.m_ );
     163    return res;
     164  }
     165
     166
     167
     168  matrix matrix::operator-( const matrix &other ) const
     169  {
     170    matrix res( *this );
     171    gsl_matrix_sub( res.m_, other.m_ );
     172    return res;
     173  }
     174
     175
     176 
     177  matrix& matrix::operator=( const matrix& other )
     178  {
     179    if ( this != &other ) {
     180      if ( m_ )
     181        gsl_matrix_free( m_ );
     182      m_ = other.gsl_matrix_copy();
    26183    }
    27   else     
    28     {
    29      m_ = gsl_matrix_alloc ( rows, cols );
     184    return *this;
     185  }
     186
     187
     188
     189  matrix matrix::transpose(void) const
     190  {
     191    matrix res(columns(),rows());
     192    gsl_matrix_transpose_memcpy(res.m_,m_);
     193    return res;
     194  }
     195
     196
     197
     198  std::ostream& operator<<(std::ostream& s, const matrix& a)
     199  {
     200    using namespace std;
     201    s.setf( ios::dec );
     202    s.precision(12);
     203    for(size_t i=0, j=0; i<a.rows(); ++i) {
     204      for (j=0; j<a.columns()-1; ++j)
     205        s << a(i,j) << " ";
     206      s << a(i,j) << endl;
    30207    }
    31 }
    32 
    33 
    34 // Is this the way to do it? No copy ...
    35 // internal data ...
    36 matrix::matrix( gsl_matrix* m ) : m_( m )
    37 {
    38 }
    39 
    40 
    41 // Copy constructor
    42 matrix::matrix( const matrix& other )
    43 {
    44   m_ = new_copy( other.get_gsl_matrix() );
    45 }
    46 
    47 
    48 // Constructor that gets data from istream
    49 matrix::matrix(std::istream& is)
    50 {
    51   using namespace std;
    52 
    53   // read the data file and store in stl vectors (dynamically expandable)
    54   std::vector<std::vector<double> > data_matrix;
    55   u_int nof_columns=0;
    56   u_int nof_rows = 0;
    57   string s;
    58   for (nof_rows = 0; getline(is, s, '\n'); nof_rows++) {
    59    istringstream line(s);
    60    std::vector<double> v;
    61    string tmp_string;
    62    while (line >> tmp_string) {
    63     if(!is.good()) {
    64      cerr << "matrix::matrix(std::istream&): "
    65     << "error reading data file!" << endl;
    66      exit(1);
    67     }
    68     double t=atof(tmp_string.c_str());
    69     // Here we should check that it actually was a double!!!!!
    70     v.push_back(t);
    71    }
    72    if(nof_columns==0)
    73      nof_columns=v.size();
    74    else if(v.size()!=nof_columns) {
    75     cerr << "matrix::matrix(std::istream&) data file error: "
    76    << "line" << nof_rows+1 << " has " << v.size()
    77    << " columns; expected " << nof_columns
    78    << " columns"
    79    << endl;
    80     exit(1);
    81    }
    82    data_matrix.push_back(v);
    83   }
    84   // manipulate the state of the stream to be good
    85   is.clear(std::ios::goodbit);
    86   // convert the data to a gsl matrix
    87   m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
    88   for(u_int i=0;i<nof_rows;i++)
    89     for(u_int j=0;j<nof_columns;j++)
    90       gsl_matrix_set( m_, i, j, data_matrix[i][j] );
    91 }
    92 
    93 matrix::~matrix()
    94 {
    95   if( m_ )
    96     {
    97      gsl_matrix_free( m_ );
    98      m_ = NULL;
    99     }
    100 }
    101 
    102 
    103 gsl_matrix* matrix::new_copy( const gsl_matrix* p_other )
    104 {
    105   // Get dimenstions
    106   size_t rows = p_other->size1;
    107   size_t cols = p_other->size2;
    108  
    109   // Create new empty matrix
    110   gsl_matrix* p_res =  gsl_matrix_alloc ( rows, cols );
    111 
    112   // Copy p_others elements into p_res
    113   gsl_matrix_memcpy( p_res, p_other );
    114 
    115   return p_res;
    116 }
    117 
    118 
    119 // Operators on matrices
    120 ////////////////////////
    121 matrix& matrix::operator=( const matrix& other )
    122 {
    123   if( this != &other )
    124     {
    125      gsl_matrix* v_new = new_copy( other.get_gsl_matrix() );
    126      if( m_ ) gsl_matrix_free( m_ );
    127      m_ = v_new;
    128     }
    129   return *this;
    130 }
    131 
    132 
    133 matrix matrix::operator+( const matrix &other ) const
    134 {
    135   assert( rows() == other.rows() &&
    136     cols() == other.cols() );
    137   matrix res( *this );
    138   gsl_matrix_add( res.get_gsl_matrix(),
    139       other.get_gsl_matrix() );
    140   return res;
    141 }
    142 
    143 
    144 matrix matrix::operator-( const matrix &other ) const
    145 {
    146   assert( rows() == other.rows() &&
    147     cols() == other.cols() );
    148   matrix res( *this );
    149   gsl_matrix_sub( res.get_gsl_matrix(),
    150       other.get_gsl_matrix() );
    151   return res;
    152 }
    153 
    154 
    155 matrix matrix::operator*( const matrix &other ) const
    156 {
    157   assert( rows() == other.cols() );
    158   matrix res( rows(), other.cols() );
    159   gsl_linalg_matmult( m_, other.get_gsl_matrix(),
    160           res.get_gsl_matrix() );
    161   return res;
    162 }
    163 
    164 
    165 // Matrix vector multiplication
    166 vector matrix::operator*( const vector &other ) const
    167 {
    168   assert( cols() == other.size() ); // NxM1 N1x1 ( M1 must equal N1 )
    169   vector res( rows(), other.is_column() );
    170 
    171   gsl_blas_dgemv( CblasNoTrans, 1.0, m_, other.get_gsl_vector(),
    172       0.0, res.get_gsl_vector() );
    173 
    174   return res;
    175 }
    176 
    177 matrix matrix::mul_elements( const matrix& a, const matrix& b )
    178 {
    179   matrix res( a );
    180   gsl_matrix_mul_elements( res.get_gsl_matrix(), b.get_gsl_matrix() );
    181   return res;
    182 }
    183 
    184 
    185 
    186 std::ostream& thep_gsl_api::operator<< ( std::ostream& s_out,
    187            const matrix& a )
    188 {
    189   using namespace std;
    190   s_out.setf( ios::dec );
    191   s_out.precision(12);
    192   for( size_t i = 0, j = 0; i < a.rows(); ++i )
    193     {
    194      for ( j = 0; j < a.cols() - 1; ++j )
    195        {
    196   s_out << a.get( i, j ) << " ";
    197        }
    198      s_out << a.get( i, j ) << endl;
    199     }
    200    return s_out;
    201 }
    202 
    203 
    204 // Public functions (A-Z)
    205 /////////////////////////
    206 matrix matrix::transpose() const
    207 {
    208   matrix res( cols(), rows() );
    209 
    210    for ( size_t i = 0; i < rows(); i++ )
    211      {
    212       for ( size_t j = 0; j < cols(); j++ )
    213   {
    214          gsl_matrix_set( res.get_gsl_matrix(), j, i, gsl_matrix_get( m_, i, j ) );
    215   }
    216      }
    217    
    218    return( res );
    219 }
    220 
    221 
    222 void matrix::swap_cols( const size_t& i, const size_t& j )
    223 {
    224   gsl_matrix_swap_columns( m_, i, j );
    225 }
    226 
    227 
    228 void matrix::swap_rows( const size_t& i, const size_t& j )
    229 {
    230   gsl_matrix_swap_rows( m_, i, j );
    231 }
    232 
    233 
    234 double matrix::norm( double n ) const
    235 {
    236   double sum = 0.0; 
    237   for ( size_t i = 0; i < rows(); ++i )
    238     {
    239      for ( size_t j = 0; j < cols(); ++j )
    240        {
    241   sum += pow( (double)( get( i, j ) ), n );
    242        }
    243     }
    244  
    245   return pow( sum, 1 / n );
    246 }
    247 
    248 
    249 matrix matrix::row( const size_t& i ) const
    250 {
    251   matrix rowmatrix( 1, cols() );
    252   gsl_vector *tempvector = gsl_vector_calloc( cols() );
    253   assert( i >= 0 && i < rows() );
    254   gsl_matrix_get_row( tempvector, m_, i );
    255   gsl_matrix_set_row( rowmatrix.get_gsl_matrix(), 0, tempvector );
    256   gsl_vector_free( tempvector );
    257   return( rowmatrix );
    258 }
    259 
    260 
    261 vector matrix::row_vector( const size_t& i ) const
    262 {
    263   gsl_vector* tmp_vector = gsl_vector_calloc( cols() );
    264   assert( i >= 0 && i < rows() );
    265   gsl_matrix_get_row( tmp_vector, m_, i );
    266   return vector( tmp_vector );
    267 }
    268 
    269 
    270 matrix matrix::col( const size_t& i ) const
    271 {
    272   matrix columnmatrix( rows(), 1 );
    273   gsl_vector *tempvector = gsl_vector_calloc( rows() );
    274   assert( i >= 0 && i < cols() );
    275   gsl_matrix_get_col( tempvector, m_, i );
    276   gsl_matrix_set_col( columnmatrix.get_gsl_matrix(), 0, tempvector );
    277   gsl_vector_free( tempvector );
    278   return( columnmatrix );
    279 }
    280 
    281 
    282 vector matrix::col_vector( const size_t& i ) const
    283 {
    284   gsl_vector* tmp_vector = gsl_vector_calloc( rows() );
    285   assert( i >= 0 && i < cols() );
    286   gsl_matrix_get_col( tmp_vector, m_, i );
    287   return vector( tmp_vector );
    288 }
    289 
    290 
    291 double matrix::sum() const
    292 {
    293   double sum = 0;
    294   for ( size_t i = 0; i < rows(); ++i )
    295     {
    296      for ( size_t j = 0; j < cols(); j++ )
    297        {
    298   sum += gsl_matrix_get( m_, i, j );
    299        }
    300     }
    301   return( sum );
    302 }
    303 
    304 
    305 vector matrix::row_sum() const
    306 {
    307   vector sum( rows() );
    308   for ( size_t i = 0; i < rows(); ++i )
    309     {
    310      sum.set( i, row( i ).sum() );
    311     }
    312   return( sum );
    313 }
    314 
    315 
    316 vector matrix::mean_row_sum() const
    317 {
    318   vector v_sum = row_sum();
    319   for( size_t i = 0; i < rows(); ++i )
    320     {
    321      v_sum[ i ] = v_sum[ i ] / static_cast<double>( cols() );
    322     }
    323   return v_sum;
    324 }
    325 
    326 vector matrix::operator[]( const size_t& i ) const
    327 {
    328   gsl_vector* tmp_vector = gsl_vector_calloc( cols() );
    329   assert( i >= 0 && i < rows() );
    330   gsl_matrix_get_row( tmp_vector, m_, i );
    331   return vector( tmp_vector );
    332 }
    333 
    334 
     208    return s;
     209  }
     210
     211
     212}} // of namespace gslapi and namespace thep
  • trunk/src/matrix.h

    r40 r42  
    1 #ifndef _thep_gsl_api_matrix_
    2 #define _thep_gsl_api_matrix_
     1// $Id$
     2
     3#ifndef _theplu_gslapi_matrix_
     4#define _theplu_gslapi_matrix_
    35
    46#include <iostream>
    5 #include <fstream>
    6 #include <iomanip>
    7 #include <cmath>
    8 #include <cstdlib>
    9 #include <cassert>
    10 
    11 #include <gsl/gsl_math.h>
     7
    128#include <gsl/gsl_matrix.h>
    13 #include <gsl/gsl_linalg.h>
    14 #include <gsl/gsl_blas.h>
    15 
    16 #include "vector.h"
    17 
    18 namespace thep_gsl_api
    19 {
    20 
    21   /**
    22      This is a matrix interface to GSL for c++.
    23      Note that the namespace is "thep_gsl_api".
    24   */
     9
     10namespace theplu {
     11namespace gslapi {
     12
     13  ///
     14  /// This is the C++ tools interface to GSL matrix.
     15  ///
     16
     17  class vector;
    2518
    2619  class matrix
    2720  {
    2821  public:
    29     typedef double Type;
    30  
    31     /**
    32        Default constructor
    33        Does not allocate any memory
    34     */
    35     matrix();
    36 
    37     /**
    38        Constructor taking three argumens:
    39        N = Number of rows
    40        M = Number of columns
    41        init_to_zero, if true matrix will be initialized with zeros.
    42     */
    43     matrix( const size_t& rows, const size_t& cols,
    44       bool init_to_zero = true );
    45 
    46 
    47 
    48 
     22
     23    ///
     24    /// The default constructor.
     25    ///
     26    matrix(void);
     27
     28    ///
     29    /// Constructor. Allocates memory space for \a r times \a c
     30    /// elements, and sets all elements to \a init_value.
     31    ///
     32    matrix(const size_t& r, const size_t& c, const double init_value=0);
    4933 
    50     /**
    51   Copy constructor
    52     */
    53     matrix( const matrix& other );
    54 
    55 
    56     /**
    57        Constructor taking as only argument a matrix of gsl type.
    58        Be causious here! Do not delete this matrix outside of
    59        class: NO COPY IS MADE!!!
    60     */
    61     matrix( gsl_matrix* );
     34    ///
     35    /// The copy constructor.
     36    ///
     37    matrix(const matrix&);
     38
     39    ///
     40    /// The istream constructor.
     41    ///
     42    /// Matrix rows are sepearated with the new line character, and
     43    /// column elements in a row must be separated with white space
     44    /// characters except the new line (\\n) character. Rows, and
     45    /// columns, are read consequently and only rectangular matrices
     46    /// are supported. End of input is at end of file marker.
     47    ///
     48    explicit matrix(std::istream &);
     49
     50    ///
     51    /// The destructor.
     52    ///
     53    ~matrix(void);
     54
     55    ///
     56    /// @return The number of columns in the matrix.
     57    ///
     58    inline size_t columns(void) const { return m_->size2; }
     59
     60    ///
     61    /// Swap column \a i with \j.
     62    ///
     63    inline void column_swap(const size_t& i,const size_t& j)
     64      { gsl_matrix_swap_columns(m_,i,j); }
     65
     66    ///
     67    /// Create a new copy of the internal GSL matrix.
     68    ///
     69    /// Necessary memory for the new GSL matrix is allocated and the
     70    /// caller is responsible for freeing the allocated memory.
     71    ///
     72    /// @return A pointer to a copy of the internal GSL matrix.
     73    ///
     74    gsl_matrix* gsl_matrix_copy(void) const;
     75
     76    ///
     77    /// @return A const pointer to the internal GSL matrix.
     78    ///
     79    inline const gsl_matrix* gsl_matrix_pointer(void) const { return m_; }
     80
     81    ///
     82    /// @return A pointer to the internal GSL matrix.
     83    ///
     84    inline gsl_matrix* gsl_matrix_pointer(void) { return m_; };
     85
     86    ///
     87    /// This function multiplies all elements between two matrices and
     88    /// returns a third matrix containing the result, \f$a_{ij} =
     89    /// b_{ij} * c_{ij} \; \forall i,j\f$.
     90    ///
     91    /// @return The result matrix.
     92    ///
     93    matrix mul_elements(const matrix&);
     94
     95    ///
     96    /// Swap row \a i with \a j.
     97    ///
     98    inline void row_swap( const size_t& i, const size_t& j)
     99      { gsl_matrix_swap_rows(m_,i,j); }
     100
     101    ///
     102    /// @return The number of rows in the matrix.
     103    ///
     104    inline size_t rows(void) const { return m_->size1; }
     105
     106    ///
     107    /// Set all elements to \a value.
     108    ///
     109    inline void set_all(const double value) { gsl_matrix_set_all(m_,value); }
     110
     111    ///
     112    /// Transpose the matrix, the object is changed.
     113    ///
     114    /// @return Transpose of the matrix.
     115    ///
     116    matrix transpose(void) const;
    62117   
    63 
    64     /**
    65        Constructor taking as only argument a an istream that should
    66        contain the matrix where columns are separated with whitespace
    67        (not '\\n') and rows by '\\n'.
    68     */
    69     matrix(std::istream &);
    70 
    71     /**
    72        Destructor will clear matrix's allocated memory
    73     */
    74     ~matrix();
    75 
    76 
    77     /**
    78        The most important function when adding new classes to
    79        c++_tools. This function will return a pointer to the
    80        inner data of the class. Should be forbidden in a strictly
    81        object-orientated environment.
    82        Motivation: Effeciency! ( CAN IT BE MADE CONST? )
    83        To be used when writting new classes using gsl-functions.
    84     */
    85     inline gsl_matrix* get_gsl_matrix() const { return m_; } 
    86 
    87 
    88     /**
    89        Returns the number of rows in the matrix.
    90     */
    91     size_t rows() const { return m_->size1; }
    92 
    93 
    94     /**
    95        Returns the number of rows in the matrix.
    96     */
    97     size_t cols() const { return m_->size2; }
    98 
    99     /**
    100        get( i, j ) will return the Aij element in matrix A
    101        where i = row number and j = column number
    102     */
    103     inline Type get( const size_t& i, const size_t& j ) const
    104       { return gsl_matrix_get(m_,i_row,j_col); }
    105 
    106 
    107     /**
    108        set( i, j, val ) will make Aij = val in matrix A
    109        where i = row number and j = column number
    110     */
    111     inline void set( const size_t& i, const size_t& j, const Type& val )
    112       { gsl_matrix_set( m_, i_row, j_col, val ); }
    113 
    114     /**
    115        set_all( val ) will make Aij = val in matrix A
    116        for all "i"
    117     */
    118     inline void set_all( const Type& val ) { gsl_matrix_set_all(m_,val); }
    119 
    120 
    121 
    122     /**
    123        Assignment-operator. No requirements on dimensions,
    124        i.e. if A = B, then A will be a copy of B with B:s
    125        dimensions.
    126      */
    127     matrix& operator=( const matrix& other );
    128    
    129     /**
    130        operators for comparison are:
    131        == equal
    132        != not equal
    133     */
     118    ///
     119    /// @return Reference to the element position (\a row, \a column).
     120    ///
     121    inline double& operator()(size_t row,size_t column)
     122      { return (*gsl_matrix_ptr(m_,row,column)); }
     123
     124    ///
     125    /// @return Const reference to the element position (\a row, \a
     126    /// column).
     127    ///
     128    inline const double& operator()(size_t row,size_t column) const
     129      { return (*gsl_matrix_const_ptr(m_,row,column)); }
     130
     131    ///
     132    /// Jari, To be implemented. Should return a reference to a vector
     133    /// object. Underlying GSL supports only GSL_vector views, thus
     134    /// some redesign of the vector class is needed.
     135    ///
     136    /// inline vector& operator[](size_t i);
     137    /// So for now, stupid copying is done, and actually a matrix
     138    /// row is returned using this function.
     139    vector operator[](size_t row) const;
     140
     141    ///
     142    /// Jari, to be properly implemented. Should return a reference to
     143    /// a vector object (matrix row). Underlying GSL supports only GSL_vector
     144    /// views, thus some redesign of the vector class is needed.
     145    ///
     146    /// inline const vector& operator[](size_t i) const;
     147    ///
     148    /// So for now, stupid copying is done, and actually a matrix
     149    /// column is returned using this function.
     150    vector TEMP_col_return(size_t column) const;
     151
     152    ///
     153    /// Matrix multiplication.
     154    ///
     155    /// @return The resulting matrix.
     156    ///
     157    matrix operator*(const matrix&) const;
     158
     159    ///
     160    /// Matrix-vector multiplication.
     161    ///
     162    /// @return The resulting vector.
     163    ///
     164    vector operator*(const vector&) const;
     165
     166    ///
     167    /// Matrix addition.
     168    ///
     169    /// @return The resulting matrix.
     170    ///
     171    matrix operator+(const matrix&) const;
     172
     173    ///
     174    /// Matrix subtraction.
     175    ///
     176    /// @return The resulting matrix.
     177    ///
     178    matrix operator-(const matrix&) const;
     179
     180    ///
     181    /// Comparison operator.
     182    ///
     183    /// @return True if ... (to be defined) ... otherwise false.
     184    ///
    134185    bool operator==( const matrix &other ) const;
     186
     187    ///
     188    /// Comparison operator.
     189    ///
     190    /// @return True if ... (to be defined) ... otherwise false.
     191    ///
    135192    bool operator!=( const matrix &other ) const;
    136193
    137     /**
    138        operators for basic arithmetics are: \n
    139        \f$+\f$ matrix-matrix addition \n
    140        \f$-\f$ matrix-matrix subtraction \n
    141        \f$*\f$ matrix-matrix multiplication
    142     */
    143     matrix operator+( const matrix &other ) const;
    144     matrix operator-( const matrix &other ) const;
    145     matrix operator*( const matrix &other ) const;
    146     vector operator*( const vector &other ) const;
    147 
    148     /**
    149   multiplying each element of two matrices
    150   dimension MUST AGREE
    151     */
    152     matrix mul_elements( const matrix& a, const matrix& b );
    153 
    154     /**
    155        Method for transposing a matrix
    156      */
    157     matrix transpose() const;
    158    
    159     /**
    160        Method for calculating the n:th sqrt of the sum of
    161        the n:th power of all elements
    162     */
    163     double norm( double n ) const;
    164 
    165     /**
    166        Calculates sum of all elements in matrix
    167     */   
    168     double matrix::sum() const;
    169      
    170     /**
    171        Exchange row i with row j (vice versa)
    172      */
    173     void swap_rows( const size_t& i, const size_t& j );
    174 
    175     /**
    176        Exchange column i with column j (vice versa)
    177      */
    178     void swap_cols( const size_t& i, const size_t& j );
    179 
    180     /**
    181        Returns row i in matrix
    182     */   
    183     matrix row( const size_t& i ) const;
    184 
    185 
    186     /**
    187        Returns row i in matrix (vector form)
    188     */
    189     vector row_vector( const size_t& i ) const;
    190 
    191 
    192     /**
    193        Returns col i in matrix
    194     */   
    195     matrix col( const size_t& i ) const;
    196 
    197 
    198     /**
    199        Returns col i in matrix (vector form)
    200     */
    201     vector col_vector( const size_t& i ) const;
    202 
    203 
    204     /**
    205        Method for calculating the row sum
    206     */   
    207     vector row_sum() const;
    208 
    209     /**
    210        Method for calculating the mean row sum
    211     */   
    212     vector mean_row_sum() const;
    213 
    214     /**
    215   standard output operator is defined
    216     */
    217     friend std::ostream& operator<< ( std::ostream& s_out, const matrix& );
    218 
    219     /**
    220        Getting vector operator
    221     */
    222     vector operator[]( const size_t& i ) const;
    223  
     194    ///
     195    /// The assignment operator. There is no requirements on
     196    /// dimensions.
     197    ///
     198    matrix& operator=(const matrix& other);
     199
     200    ///
     201    /// Multiply and assign operator.
     202    ///
     203    inline matrix& operator*=(const double d)
     204      { gsl_matrix_scale(m_,d); return *this; }
     205
     206
    224207  private:
    225     gsl_matrix* new_copy( const gsl_matrix* );
    226208
    227209    gsl_matrix* m_;
    228   };
    229 
    230 
    231 };
    232 
    233 
    234 
     210
     211  };
     212
     213  ///
     214  /// The output operator for the vector class.
     215  ///
     216  std::ostream& operator<< (std::ostream& s, const matrix&);
     217
     218
     219
     220}} // of namespace gslapi and namespace thep
    235221
    236222#endif
  • trunk/src/random_singleton.cc

    r23 r42  
    33#include "random_singleton.h"
    44
    5 using namespace thep_cpp_tools;
     5namespace theplu {
     6namespace cpptools { 
    67
    78// Static members
     
    5657}
    5758
    58 
     59}} // of namespace cpptools and namespace theplu
  • trunk/src/random_singleton.h

    r34 r42  
    11// $Id$
    22
    3 #ifndef _THEP_RANDOMSINGLETON_
    4 #define _THEP_RANDOMSINGLETON_
    5 
    6 
    7 #include <cstdlib>
    8 #include <cstdio>
     3#ifndef _theplu_cpptools_random_singleton_
     4#define _theplu_cpptools_random_singleton_
     5
     6
     7// #include <cstdlib>
     8// #include <cstdio>
    99#include <string>
    1010
    11 // GSL INCLUDES
    12 ////////////////////////////
    1311#include <gsl/gsl_rng.h>
    1412#include <gsl/gsl_randist.h>
    1513
    16 
    17 namespace thep_cpp_tools
    18 
     14namespace theplu {
     15namespace cpptools { 
     16
    1917  /**
    2018     Class defining interface for different random
     
    135133       returns the one used.
    136134    */
    137     std::string get_generator_type() const;
     135    std::string get_generator_type(void) const;
    138136
    139137
     
    205203}
    206204
    207 } // namespace thep_c++_tools
     205}} // of namespace cpptools and namespace theplu
    208206
    209207#endif
  • trunk/src/vector.cc

    r37 r42  
    11// $Id$
    22
    3 // #include <iostream>
     3#include <iostream>
     4#include <sstream>
    45#include <string>
    5 #include <sstream>
    66#include <vector>
     7
    78#include "vector.h"
    89
    910
    10 using namespace thep_gsl_api;
    11 
    12 
    13 // Constructors and Destructors
    14 ///////////////////////////////
    15 vector::vector() : v_( NULL )
    16 {
    17   is_col_ = true;
    18 }
    19 
    20 
    21 vector::vector( const size_t& N, bool is_col_vector,
    22     bool init_to_zero ) :
    23   is_col_( is_col_vector )
    24 {
    25   if( init_to_zero )
    26     {
    27      v_ = gsl_vector_calloc( N );
     11namespace theplu {
     12namespace gslapi {
     13
     14
     15
     16  vector::vector(void)
     17    : v_(NULL)
     18  {
     19  }
     20
     21
     22
     23  vector::vector(const size_t n, const double init_value)
     24  {
     25    v_ = gsl_vector_alloc(n);
     26    set_all(init_value);
     27  }
     28
     29
     30
     31  vector::vector(const vector& other)
     32  {
     33    gsl_vector_free(v_);
     34    v_ = other.gsl_vector_copy();
     35  }
     36
     37
     38
     39  vector::vector(gsl_vector* vec)
     40    : v_(vec)
     41  {
     42  }
     43
     44
     45
     46  vector::vector(std::istream& is)
     47  {
     48    using namespace std;
     49 
     50    // read the data file and store in stl vectors (dynamically expandable)
     51    std::vector<std::vector<double> > data_matrix;
     52    u_int nof_columns=0;
     53    u_int nof_rows = 0;
     54    string s;
     55    for (nof_rows = 0; getline(is, s, '\n'); nof_rows++) {
     56      istringstream line(s);
     57      std::vector<double> v;
     58      string tmp_string;
     59      while (line >> tmp_string) {
     60        if(!is.good()) {
     61          // Jari. change to throw exception, remove unnecessary includes
     62          cerr << "vector::vector(std::istream&): "
     63               << "error reading data file!" << endl;
     64          exit(1);
     65        }
     66        double t=atof(tmp_string.c_str());
     67        // Here we should check that it actually was a double!!!!!
     68        v.push_back(t);
     69      }
     70      if(nof_columns==0)
     71        nof_columns=v.size();
     72      else if(v.size()!=nof_columns) {
     73        // Jari. change to throw exception, remove unnecessary includes
     74        cerr << "vector::vector(std::istream&) data file error: "
     75             << "line" << nof_rows+1 << " has " << v.size()
     76             << " columns; expected " << nof_columns
     77             << " columns"
     78             << endl;
     79        exit(1);
     80      }
     81      data_matrix.push_back(v);
     82    }
     83 
     84    // manipulate the state of the stream to be good
     85    is.clear(std::ios::goodbit);
     86 
     87    // convert the data to a gsl vector and check that data file is a
     88    // column vector or a row vector
     89    if(nof_columns==1) {
     90      v_ = gsl_vector_alloc ( nof_rows );
     91      for(u_int i=0;i<nof_rows;i++)
     92        gsl_vector_set( v_, i, data_matrix[i][0] );
     93    }
     94    else if(nof_rows==1){
     95      v_ = gsl_vector_alloc ( nof_columns );
     96      for(u_int i=0;i<nof_columns;i++)
     97        gsl_vector_set( v_, i, data_matrix[0][i] );
     98    }
     99    else {
     100      // Jari. change to throw exception, remove unnecessary includes
     101      cerr << "vector::vector(std::istream&) data file error: "
     102           << "file has " << nof_rows << " rows and " << nof_columns
     103           << " columns; expected a row vector or a column vector"
     104           << endl;
     105      exit(1);
     106    }
     107  }
     108
     109
     110
     111  vector::~vector(void)
     112  {
     113    if (v_) {
     114      gsl_vector_free(v_);
     115      v_=NULL;
    28116    }
    29   else     
    30     {
    31      v_ = gsl_vector_alloc ( N );
     117  }
     118
     119
     120
     121  gsl_vector* vector::gsl_vector_copy(void) const
     122  {
     123    gsl_vector* vec = gsl_vector_alloc(size());
     124    gsl_vector_memcpy(vec,v_);
     125    return vec;
     126  }
     127
     128
     129
     130  vector vector::mul_elements( const vector& other ) const
     131  {
     132    vector res( *this );
     133    gsl_vector_mul( res.v_, other.v_ );
     134    return res;
     135  }
     136
     137
     138
     139  double vector::sum(void) const
     140  {
     141    double sum = 0;
     142    for (u_int i=0; i<size(); i++)
     143      sum += gsl_vector_get( v_, i );
     144    return( sum );
     145  } 
     146
     147
     148
     149  vector vector::operator-(void) const
     150  {
     151    vector res( *this );
     152    gsl_vector_scale (res.v_,-1.);
     153    return res;
     154  }
     155
     156
     157
     158  double vector::operator*( const vector &other ) const
     159  {
     160    // Jari, check for gsl_support
     161    double res = 0.0;;
     162    for ( size_t i = 0; i < size(); ++i )
     163      res += other[i] * (*this)[i];
     164    return res;
     165  }
     166
     167
     168
     169  vector vector::operator+(const vector &other) const
     170  {
     171    vector res(*this);
     172    gsl_vector_add(res.v_,other.v_);
     173    return res;
     174  }
     175
     176
     177
     178  vector vector::operator-( const vector &other ) const
     179  {
     180    vector res( *this );
     181    gsl_vector_sub(res.v_,other.v_);
     182    return res;
     183  }
     184
     185
     186
     187  vector& vector::operator=( const vector& other )
     188  {
     189    if( this != &other ) {
     190      if ( v_ )
     191        gsl_vector_free( v_ );
     192      v_ = other.gsl_vector_copy();
    32193    }
    33 }
    34 
    35 
    36 // Is this the way to do it? No copy ...
    37 // internal data ...
    38 vector::vector( gsl_vector* v, bool is_col ) : v_( v ), is_col_( is_col )
    39 {
    40 }
    41 
    42 
    43 // Copy constructor
    44 vector::vector( const vector& other )
    45 {
    46   v_ = new_copy( other.get_gsl_vector() );
    47   is_col_ = other.is_column();
    48 }
    49 
    50 
    51 // Constructor that gets data from istream
    52 vector::vector(std::istream& is)
    53 {
    54   using namespace std;
    55  
    56   // read the data file and store in stl vectors (dynamically expandable)
    57   std::vector<std::vector<double> > data_matrix;
    58   u_int nof_columns=0;
    59   u_int nof_rows = 0;
    60   string s;
    61   for (nof_rows = 0; getline(is, s, '\n'); nof_rows++) {
    62    istringstream line(s);
    63    std::vector<double> v;
    64    string tmp_string;
    65    while (line >> tmp_string) {
    66     if(!is.good()) {
    67      cerr << "vector::vector(std::istream&): "
    68     << "error reading data file!" << endl;
    69      exit(1);
    70     }
    71     double t=atof(tmp_string.c_str());
    72     // Here we should check that it actually was a double!!!!!
    73     v.push_back(t);
    74    }
    75    if(nof_columns==0)
    76      nof_columns=v.size();
    77    else if(v.size()!=nof_columns) {
    78     cerr << "vector::vector(std::istream&) data file error: "
    79    << "line" << nof_rows+1 << " has " << v.size()
    80    << " columns; expected " << nof_columns
    81    << " columns"
    82    << endl;
    83     exit(1);
    84    }
    85    data_matrix.push_back(v);
    86   }
    87  
    88   // manipulate the state of the stream to be good
    89   is.clear(std::ios::goodbit);
    90  
    91   // convert the data to a gsl vector and check that data file is a
    92   // column vector or a row vector
    93   if(nof_columns==1) {
    94    v_ = gsl_vector_alloc ( nof_rows );
    95    is_col_ = true;
    96    for(u_int i=0;i<nof_rows;i++)
    97      gsl_vector_set( v_, i, data_matrix[i][0] );
    98   }
    99   else if(nof_rows==1){
    100    v_ = gsl_vector_alloc ( nof_columns );
    101    is_col_ = false;
    102    for(u_int i=0;i<nof_columns;i++)
    103      gsl_vector_set( v_, i, data_matrix[0][i] );
    104   }
    105   else {
    106     cerr << "vector::vector(std::istream&) data file error: "
    107    << "file has " << nof_rows << " rows and " << nof_columns
    108    << " columns; expected a row vector or a column vector"
    109    << endl;
    110     exit(1);
    111   }
    112 }
    113 
    114  
    115 vector::~vector()
    116 {
    117   if( v_ )
    118     {
    119      gsl_vector_free( v_ );
    120      v_ = NULL;
    121     }
    122 }
    123 
    124 
    125 gsl_vector* vector::new_copy( const gsl_vector* p_other )
    126 {
    127   // Get dimenstions
    128   size_t N = p_other->size;
    129    
    130   // Create new empty vector
    131   gsl_vector* p_res =  gsl_vector_alloc( N );
    132 
    133   // Copy p_others elements into p_res
    134   gsl_vector_memcpy( p_res, p_other );
    135 
    136   return p_res;
    137 }
    138 
    139 
    140 // Operators on vectors
    141 ////////////////////////
    142 vector& vector::operator=( const vector& other )
    143 {
    144   if( this != &other )
    145     {
    146      gsl_vector* v_new = new_copy( other.get_gsl_vector() );
    147      if( v_ ) gsl_vector_free( v_ );
    148      v_ = v_new;
    149     }
    150   return *this;
    151 }
    152 
    153 
    154 vector vector::operator+( const vector &other ) const
    155 {
    156   assert( size() == other.size() );
    157   vector res( *this );
    158   gsl_vector_add( res.get_gsl_vector(), other.get_gsl_vector() );
    159   return res;
    160 }
    161 
    162 
    163 vector vector::operator-( const vector &other ) const
    164 {
    165   assert( size() == other.size() );
    166   vector res( *this );
    167   gsl_vector_sub( res.get_gsl_vector(), other.get_gsl_vector() );
    168   return res;
    169 }
    170 
    171 vector vector::operator-() const
    172 {
    173   vector res( *this );
    174   gsl_vector_scale (res.get_gsl_vector() , -1.);
    175   return res;
    176 }
    177 double vector::operator*( const vector &other ) const
    178 {
    179   assert( size() == other.size() );
    180   double res = 0.0;;
    181   for ( size_t i = 0; i < other.size(); ++i )
    182     {
    183      res += other.get( i ) * get( i );
    184     }
    185   return res;
    186 }
    187 
    188 
    189 int vector::operator/=( const vector& other ) const
    190 {
    191   return gsl_vector_div( v_, other.get_gsl_vector() );
    192 }
    193 
    194 // Each element is divided by the "constant"
    195 int vector::operator/=( const double& constant )
    196 {
    197   assert( constant != 0 );
    198   return gsl_vector_scale( v_, 1 / constant );
    199 }
    200 
    201 // This is not the dot product, hence a matrix
    202 // is returned
    203 // matrix vector::operator*( const vector &other ) const
    204 // {
    205 //   assert( rows() == other.cols() );
    206 //   matrix res( rows(), other.cols() );
    207 //   gsl_linalg_matmult( v_, other.get_gsl_vector(),
    208 //          res.get_gsl_vector() );
    209 //   return res;
    210 // }
    211 
    212 
    213 
    214 vector vector::mul_elements( const vector& other ) const
    215 {
    216   assert( size() == other.size() );
    217   vector res( *this );
    218   gsl_vector_mul( res.get_gsl_vector(), other.get_gsl_vector() );
    219   return res;
    220 }
    221 
    222 std::ostream& thep_gsl_api::operator<< ( std::ostream& s_out,
    223            const vector& a )
    224 {
    225   using namespace std;
    226   s_out.setf( ios::fixed ); 
    227 
    228   if( a.is_column() )
    229     {
    230      for ( size_t i = 0; i < a.size(); ++i )
    231        {
    232   s_out << a.get( i  ) << endl;
    233        }
    234     }
    235   else
    236     {
    237      for ( size_t j = 0; j < a.size(); ++j )
    238        {
    239   s_out << a.get( j ) << " ";
    240        }
    241     }
    242 
    243   return s_out;
    244 }
    245 
    246 double vector::sum() const
    247 {
    248   double sum = 0;
    249   for (u_int i=0; i<size(); i++)
    250    sum += gsl_vector_get( v_, i );
    251   return( sum );
    252 
    253  
    254  
    255 
    256 // // Public functions (A-Z)
    257 // /////////////////////////
    258 
    259 
     194    return *this;
     195  }
     196
     197
     198
     199  std::ostream& theplu::gslapi::operator<<(std::ostream& s, const vector& a)
     200  {
     201    s.setf(std::ios::fixed);
     202
     203    for (size_t j = 0; j < a.size(); ++j) {
     204      s << a[j];
     205      if ( (j+1)<a.size() )
     206        s << " ";
     207    }
     208
     209    return s;
     210  }
     211
     212
     213}} // of namespace gslapi and namespace thep
  • trunk/src/vector.h

    r40 r42  
    11// $Id$
    22
    3 #ifndef _thep_gsl_api_vector_
    4 #define _thep_gsl_api_vector_
     3#ifndef _theplu_gslapi_vector_
     4#define _theplu_gslapi_vector_
    55
    66#include <iostream>
    7 #include <fstream>
    8 #include <iomanip>
    9 #include <cmath>
    10 #include <cstdlib>
    11 #include <cassert>
    127
    13 #include <gsl/gsl_math.h>
    148#include <gsl/gsl_vector.h>
    15 #include <gsl/gsl_linalg.h>
    16 #include <gsl/gsl_blas.h>
     9
     10namespace theplu {
     11// Jari, this should probably go somewhere else for doxygen to process.
     12///
     13/// All GSL wrapper functionality is to be defined within gslapi
     14/// namespace.
     15///
     16/// There are a number of support operators and functions for the
     17/// wrapper classes that do not belong to the classes themselves, but
     18/// still are defined in this namespace.
     19///
     20namespace gslapi {
     21
     22  ///
     23  /// This is the C++ tools interface to GSL vector.
     24  ///
     25
     26  class vector
     27  {
     28  public:
     29
     30    ///
     31    /// The default constructor.
     32    ///
     33    vector(void);
     34
     35    ///
     36    /// Constructor. Allocates memory space for \a n elements, and
     37    /// sets all elements to \a init_value.
     38    ///
     39    vector(const size_t n, const double init_value=0);
     40
     41    ///
     42    /// The copy constructor.
     43    ///
     44    vector(const vector&);
     45
     46    ///
     47    /// Constructor that imports a GSL vector. The GSL object is owned
     48    /// by the created object.
     49    ///
     50    vector(gsl_vector*);
     51
     52    ///
     53    /// The istream constructor.
     54    ///
     55    /// Elements should be separated with white space characters, the
     56    /// end of input to the vector is at end of file marker.
     57    ///
     58    explicit vector(std::istream &);
     59
     60    ///
     61    /// The destructor.
     62    ///
     63    ~vector(void);
     64
     65    ///
     66    /// Create a new copy of the internal GSL vector.
     67    ///
     68    /// Necessary memory for the new GSL vector is allocated and the
     69    /// caller is responsible for freeing the allocated memory.
     70    ///
     71    /// @return A pointer to a copy of the internal GSL vector.
     72    ///
     73    gsl_vector* gsl_vector_copy(void) const;
     74
     75    ///
     76    /// @return A const pointer to the internal GSL vector,
     77    ///
     78    inline const gsl_vector* gsl_vector_pointer(void) const { return v_; }
     79
     80    ///
     81    /// @return A pointer to the internal GSL vector,
     82    ///
     83    inline gsl_vector* gsl_vector_pointer(void) { return v_; };
     84
     85    ///
     86    /// This function multiplies all elements between two vectors and
     87    /// returns a third vector containing the result, \f$a_i = b_i *
     88    /// c_i \; \forall i\f$.
     89    ///
     90    /// @return The result vector.
     91    ///
     92    vector vector::mul_elements( const vector& other ) const;
     93
     94    ///
     95    /// Set all elements to \a value.
     96    ///
     97    inline void set_all(const double& value) { gsl_vector_set_all(v_,value); }
     98
     99    ///
     100    /// @return the number of elements in the vector.
     101    ///
     102    inline size_t size(void) const { return v_->size; }
     103
     104    ///
     105    /// Calculate the sum of all vector elements.
     106    ///
     107    /// @return The sum.
     108    ///
     109    double vector::sum(void) const;
     110
     111    ///
     112    /// Element access operator.
     113    ///
     114    /// @return Reference to element \a i.
     115    ///
     116    inline double& operator()(size_t i) { return *gsl_vector_ptr(v_,i); }
     117
     118    ///
     119    /// Const element access operator.
     120    ///
     121    /// @return The value of element \a i.
     122    ///
     123    inline const double& operator()(size_t i) const
     124      { return *gsl_vector_ptr(v_,i); }
     125   
     126    ///
     127    /// Element access operator.
     128    ///
     129    /// @return Reference to element \a i.
     130    ///
     131    inline double& operator[](size_t i) { return (*this)(i); }
     132
     133    ///
     134    /// Const element access operator.
     135    ///
     136    /// @return The value of element \a i.
     137    ///
     138    inline const double& operator[](size_t i) const { return (*this)(i); }
     139   
     140    ///
     141    /// @return The negation of the vector.
     142    ///
     143    vector operator-(void) const;
     144
     145    ///
     146    /// @return The dot product.
     147    ///
     148    double operator*( const vector &other ) const;
     149
     150    ///
     151    /// Vector addition.
     152    ///
     153    /// @return The resulting vector.
     154    ///
     155    vector operator+( const vector& other ) const;
     156
     157    ///
     158    /// Vector subtraction.
     159    ///
     160    /// @return The resulting vector.
     161    ///
     162    vector operator-( const vector& other ) const;
     163
     164    ///
     165    /// The assignment operator. There is no requirements on
     166    /// dimensions.
     167    ///
     168    vector& operator=(const vector&);
     169
     170    ///
     171    /// Multiply and assign operator.
     172    ///
     173    vector& operator*=(const double d) { gsl_vector_scale(v_,d); return *this; }
    17174
    18175
    19 namespace thep_gsl_api
    20 {
     176  private:
    21177
    22   /**
    23      This is a vector interface to GSL for c++.
    24      Important: vector has two dimenstions depending on
    25      whether it's a "row" vector or a "column" vector.
    26      A column vector is created by default.
    27      Note that the namespace is "thep_gsl_api".
    28   */
     178    gsl_vector* v_;
    29179
    30   class vector
    31   {
    32   public:
    33     typedef double Type;
    34  
    35     /**
    36        Default constructor
    37        Does not allocate any memory
    38     */
    39     vector();
     180  };
    40181
    41     /**
    42        Constructor taking three argumens:
    43        N = Number of elements
    44        init_to_zero, if true vector will be initialized with zeros.
    45     */
    46     vector( const size_t&, bool is_col_vector = true,
    47       bool init_to_zero = true );
     182  ///
     183  /// The output operator for the vector class.
     184  ///
     185  std::ostream& operator<< ( std::ostream&, const vector& );
    48186
    49187
    50 
    51     /**
    52   Copy constructor
    53     */
    54     vector( const vector& other );
    55 
    56 
    57     /**
    58        Constructor taking as only argument a matrix of gsl type.
    59        Be causious here! Do not delete this matrix outside of
    60        class: NO COPY IS MADE!!!
    61     */
    62     vector( gsl_vector*, bool is_col = true );
    63 
    64     /**
    65        Constructor taking as only argument a an istream that should
    66        contain the vector. Column vector should be sepated with '\\n'
    67        and row vector should be separeted with whitespace (not '\\n').       
    68 
    69     */
    70     vector(std::istream &);
    71 
    72 
    73     /**
    74        Destructor will clear matrix's allocated memory
    75     */
    76     ~vector();
    77 
    78 
    79     /**
    80        The most important function when adding new classes to
    81        c++_tools. This function will return a pointer to the
    82        inner data of the class. Should be forbidden in a strictly
    83        object-orientated environment.
    84        Motivation: Effeciency! ( CAN IT BE MADE CONST? )
    85        To be used when writting new classes using gsl-functions.
    86     */
    87     inline gsl_vector* get_gsl_vector() const { return v_; }
    88 
    89 
    90     /**
    91        Returns the number of rows in the vector.
    92     */
    93     size_t size() const { return v_->size; }
    94 
    95 
    96     /**
    97        Tells whether a vector is column or not
    98      */
    99     bool is_column() const { return is_col_; }
    100 
    101 
    102     /**
    103        get( i ) will return the vi element in vector v
    104        where i = column (if row vector) or row (if column vector)
    105     */
    106     inline Type get( const size_t& i ) const { return gsl_vector_get(v_,i); }
    107 
    108 
    109     /**
    110        set( i, val ) will make vi = val in vector v
    111        where i = column (if row vector) or row (if column vector)
    112     */
    113     inline void set( const size_t& i, const Type& val )
    114       { gsl_vector_set(v_,i,val); }
    115 
    116     /**
    117        set_all( val ) will make vi = val in vector v
    118        for all "i"
    119     */
    120     inline void set_all( const Type& val ) { gsl_vector_set_all(v_,val); }
    121 
    122 
    123     /**
    124        Assignment-operator. No requirements on dimensions,
    125        i.e. if a = b, then a will be a copy of b with b:s
    126        dimensions.
    127      */
    128     vector& operator=( const vector& other );
    129 
    130 
    131     /**
    132        Access element-operator
    133     */
    134     inline double &operator[](size_t i)  { return *gsl_vector_ptr( v_, i ); }
    135     const double &operator[](size_t i) const { return *gsl_vector_ptr( v_, i ); }
    136 
    137    
    138     /**
    139        operators for comparison are: \n
    140        == equal \n
    141        != not equal
    142     */
    143     bool operator==( const vector &other ) const;
    144     bool operator!=( const vector &other ) const;
    145 
    146     /**
    147        operators for basic arithmetics are: \n
    148        \f$+\f$ vector-vector addition \n
    149        \f$-\f$ vector-vector subtraction \n
    150        \f$/=\f$ divide all elements with a constant
    151     */
    152     vector operator+( const vector& other ) const;
    153     vector operator-( const vector& other ) const;
    154     int operator/=( const double& c );
    155     int operator/=( const vector& other ) const;
    156     vector operator-() const;
    157     /**
    158        This operator is implemented as dot-product.
    159        I case of discussions later on; Jari decided
    160        this!
    161     */
    162     double operator*( const vector &other ) const;
    163 
    164     /**
    165        This function multiplies all elements in two vectors
    166        and returns a third vector containing the result.
    167        res = a*b; where size(a) = size(b) = size(res)
    168     */
    169     vector vector::mul_elements( const vector& other ) const;
    170    
    171    
    172     /**
    173        Calculates sum of all elements in vector
    174     */
    175     double vector::sum() const;
    176    
    177     /**
    178         standard output operator is defined
    179     */
    180     friend std::ostream& operator<< ( std::ostream&, const vector& );
    181    
    182 
    183   private:
    184     gsl_vector* new_copy( const gsl_vector* );
    185 
    186     gsl_vector* v_;
    187     bool is_col_;
    188   };
    189 
    190 
    191 };
    192 
    193 
    194 
     188}} // of namespace gslapi and namespace thep
    195189
    196190#endif
  • trunk/test/Makefile.am

    r31 r42  
    33## $Id$
    44
    5 check_PROGRAMS = test_rnd test_pca test_svd test_kernel test_svm
     5check_PROGRAMS = test_kernel test_pca test_rnd test_svd test_svm
    66
    77test_rnd_SOURCES = test_rnd.cc
  • trunk/test/test_kernel.cc

    r41 r42  
    1414// Standard includes
    1515////////////////////
     16#include <fstream>
    1617#include <iostream>
    1718#include <cstdlib>
     
    1920
    2021using namespace std;
    21 using namespace thep_gsl_api;
    22 using namespace thep_cpp_tools;
    2322
    2423int main()
     
    2827  cout << "testing the linear kernel" << endl;   
    2928  ifstream is("nm_data_centralized.txt");
    30   thep_gsl_api::matrix data(is);
     29  theplu::gslapi::matrix data(is);
    3130  is.close();
    3231 
    3332  is.open("nm_kernel.txt");
    34   thep_gsl_api::matrix kernel_matlab(is);
     33  theplu::gslapi::matrix kernel_matlab(is);
    3534  is.close();
    36   thep_cpp_tools::KernelFunction* kf = new
    37     thep_cpp_tools::PolynomialKernelFunction();
    38   thep_cpp_tools::Kernel kernel(data,*kf);
    39   thep_gsl_api::matrix diff;
    40   diff = kernel.get()-kernel_matlab;
     35  theplu::cpptools::KernelFunction* kf =
     36    new theplu::cpptools::PolynomialKernelFunction();
     37  theplu::cpptools::Kernel kernel(data,*kf);
     38  theplu::gslapi::matrix diff(kernel.get()-kernel_matlab);
    4139  double tmp=0;
    4240  for(u_int i=0;i<diff.rows();i++)
    4341    for(u_int j=0;j<diff.rows();j++)
    44       tmp+=abs(diff.get(i,j)/kernel_matlab.get(i,j));
     42      tmp+=abs(diff(i,j)/kernel_matlab(i,j));
    4543  cout << "The mean relative difference is:" << endl; 
    4644  cout << tmp/diff.rows()/diff.rows() << endl;
     
    4947  cout << "testing the polynomial kernel of degree 2" << endl;   
    5048  is.open("nm_kernel2.txt");
    51   thep_gsl_api::matrix kernel_matlab2(is);
     49  theplu::gslapi::matrix kernel_matlab2(is);
    5250  is.close();
    53   kf = new thep_cpp_tools::PolynomialKernelFunction(2);
    54   thep_cpp_tools::Kernel kernel2(data,*kf);
     51  kf = new theplu::cpptools::PolynomialKernelFunction(2);
     52  theplu::cpptools::Kernel kernel2(data,*kf);
    5553  diff = kernel2.get()-kernel_matlab2;
    5654  tmp=0;
    5755  for(u_int i=0;i<diff.rows();i++)
    5856    for(u_int j=0;j<diff.rows();j++)
    59       tmp+=abs(diff.get(i,j)/kernel_matlab2.get(i,j));
     57      tmp+=abs(diff(i,j)/kernel_matlab2(i,j));
    6058  cout << "The mean relative difference is:" << endl; 
    6159  cout << tmp/diff.rows()/diff.rows() << endl;
    6260   
    63   diff.set(1,2,32);
    64   diff.set(2,1,111);
     61  diff(1,2)=32;
     62  diff(2,1)=111;
    6563 
    66   tmp=diff[1][2];
     64  tmp=diff(1,2);
    6765  cout << tmp << endl;
    6866
  • trunk/test/test_pca.cc

    r13 r42  
    55#include <fstream>
    66#include <string>
     7#include <vector>
     8
    79#include "PCA.h"
    810
     
    1517int main()
    1618{
    17   thep_cpp_tools::PCA* mypca = new thep_cpp_tools::PCA;
     19  theplu::cpptools::PCA* mypca = new theplu::cpptools::PCA;
    1820  if( mypca->test() ) { delete mypca; return 0; }
    1921  else { delete mypca; return -1; }
     
    4749
    4850  const size_t cols = data.size();
    49   assert( data.size() > 0 );
    5051  const size_t rows = data[0].size();
    5152 
    52   thep_gsl_api::matrix A( cols, rows );
     53  theplu::gslapi::matrix A( cols, rows );
    5354
    5455  // Copy data into matrix A
    5556  for( size_t i = 0; i < cols; ++i )
    56     {
    57      for( size_t j = 0; j < rows; ++j )
    58        {
    59   A.set( i, j, data[i][j] );
    60        }
    61     }
     57    for( size_t j = 0; j < rows; ++j )
     58      A(i,j)= data[i][j];
    6259
    6360  //  cout << A << endl;
    6461
    65   thep_cpp_tools::PCA* mypca = new thep_cpp_tools::PCA( A.transpose() );
     62  theplu::cpptools::PCA* mypca = new theplu::cpptools::PCA( A.transpose() );
    6663  mypca->process();
    6764
  • trunk/test/test_rnd.cc

    r13 r42  
     1// $Id$
     2
    13#include <iostream>
    24#include "random_singleton.h"
    35
    4 using namespace thep_cpp_tools;
     6using namespace theplu::cpptools;
    57using namespace std;
    68
  • trunk/test/test_svd.cc

    r8 r42  
    77int main()
    88{
    9   thep_cpp_tools::SVD* mysvd = new thep_cpp_tools::SVD;
     9  theplu::cpptools::SVD* mysvd = new theplu::cpptools::SVD;
    1010  if( mysvd->test() ) { delete mysvd; return 0; }
    1111  else { delete mysvd; return -1; }
  • trunk/test/test_svm.cc

    r31 r42  
    11// $Id$
    2 
    3 
    4 #ifndef _THEP_CPPTOOLS_SVM_
    5 #define _THEP_CPPTOOLS_SVM_
    62
    73// C++ tools include
     
    139// Standard includes
    1410////////////////////
     11#include <fstream>
    1512#include <iostream>
    1613#include <cstdlib>
    1714
    1815using namespace std;
    19 using namespace thep_gsl_api;
    20 using namespace thep_cpp_tools;
    2116
    2217int main()
     
    2621  //ifstream is("k.txt");
    2722  ifstream is("nm_kernel.txt");
    28   thep_gsl_api::matrix kernel(is);
     23  theplu::gslapi::matrix kernel(is);
    2924  is.close();
    3025 
     
    3227  is.open("nm_target_bin.txt");
    3328  //is.open("t.txt");
    34   thep_gsl_api::vector target(is);
     29  theplu::gslapi::vector target(is);
    3530  is.close();
    3631
    3732  cout << "nm_alpha_linear_matlab.txt" << endl;   
    3833  is.open("nm_alpha_linear_matlab.txt");
    39   thep_gsl_api::vector alpha_matlab(is);
     34  theplu::gslapi::vector alpha_matlab(is);
    4035  is.close();
    4136
    4237  //training the SVM
    4338  cout << "Training the SVM" << endl;   
    44   thep_cpp_tools::SVM svm(kernel,target);
     39  theplu::cpptools::SVM svm(kernel,target);
    4540  svm.train();
    4641 
    4742  //getting alpha
    4843  cout << "getting alpha" << endl;   
    49   thep_gsl_api::vector alpha;
     44  theplu::gslapi::vector alpha;
    5045  alpha =  svm.get_alpha();
    5146     
    5247  // Comparing alpha to alpha_matlab
    53   thep_gsl_api::vector diff_alpha = alpha - alpha_matlab;
     48  theplu::gslapi::vector diff_alpha = alpha - alpha_matlab;
    5449  cout << "Checking difference" << endl;   
    5550  cout << diff_alpha*diff_alpha << endl;
     
    6156  return 0;
    6257}
    63 
    64 
    65 
    66 #endif
    67 
    68 
    69 
    70 
    71 
Note: See TracChangeset for help on using the changeset viewer.