Changeset 420


Ignore:
Timestamp:
Dec 2, 2005, 1:15:50 AM (18 years ago)
Author:
Jari Häkkinen
Message:

Merged better_matrix_class changes r402:419 into the trunk.

Location:
trunk
Files:
30 edited

Legend:

Unmodified
Added
Removed

  • trunk/lib/gslapi/matrix.cc

    r357 r420  
    77
    88#include <cmath>
    9 #include <iostream>
    109#include <sstream>
    11 #include <string>
    1210#include <vector>
    1311
    1412#include <gsl/gsl_blas.h>
    15 #include <gsl/gsl_linalg.h>
    1613
    1714
    … …  
    2118
    2219
    23   matrix::matrix(void)
    24     : m_(NULL)
     20  matrix::matrix(matrix& m, size_t offset_row, size_t offset_column,
     21                 size_t n_row, size_t n_column)
    2522  {
    26   }
    27 
    28 
    29 
    30   matrix::matrix(const size_t& r, const size_t& c, double init_value)
    31   {
    32     m_ = gsl_matrix_alloc(r,c);
    33     set_all(init_value);
    34   }
    35 
    36 
    37 
    38   matrix::matrix(const matrix& other)
    39   {
    40     m_ = other.gsl_matrix_copy();
     23    // Jari, exception handling needed here. Failure in setting up a
     24    // proper gsl_matrix_view is communicated by NULL pointer in the
     25    // view structure (cf. GSL manual). How about GSL error state?
     26    view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_,
     27                                                     offset_row,offset_column,
     28                                                     n_row,n_column));
     29    m_ = &(view_->matrix);
    4130  }
    4231
    … …  
    4433
    4534  // Constructor that gets data from istream
    46   matrix::matrix(std::istream& is)
     35  matrix::matrix(std::istream& is) throw (utility::IO_error,std::exception)
     36    : view_(NULL)
    4737  {
    4838    // read the data file and store in stl vectors (dynamically
    … …  
    5343    std::vector<double> v;
    5444    for (nof_rows = 0; utility::read_to_double(is, v); nof_rows++) {
    55    
    5645      // Ignoring empty lines
    57       if(!v.size()){
     46      if (!v.size()) {
    5847        nof_rows--;
    5948        continue;
    6049      }
    61      
    62       if(nof_columns==0)
     50      if (!nof_columns)
    6351        nof_columns=v.size();
    64       else if(v.size()!=nof_columns) {
    65         std::cerr << "matrix::matrix(std::istream&) data file error: "
    66                   << "line" << nof_rows+1 << " has " << v.size()
    67                   << " columns; expected " << nof_columns
    68                   << " columns"
    69                   << std::endl;
    70         exit(1);
     52      else if (v.size()!=nof_columns) {
     53        std::ostringstream s;
     54        s << "matrix::matrix(std::istream&) data file error: "
     55          << "line" << nof_rows+1 << " has " << v.size()
     56          << " columns; expected " << nof_columns << " columns.";
     57        throw utility::IO_error(s.str());
    7158      }
    7259      data_matrix.push_back(v);
    … …  
    7764    // convert the data to a gsl matrix
    7865    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
    79     for(u_int i=0;i<nof_rows;i++) 
    80       for(u_int j=0;j<nof_columns;j++) 
     66    for(u_int i=0;i<nof_rows;i++)
     67      for(u_int j=0;j<nof_columns;j++)
    8168        gsl_matrix_set( m_, i, j, data_matrix[i][j] );
    8269  }
    … …  
    8673  matrix::~matrix(void)
    8774  {
    88     if (m_) {
     75    if (view_)
     76      delete view_;
     77    else if (m_)
    8978      gsl_matrix_free(m_);
    90       m_=NULL;
     79    m_=NULL;
     80  }
     81
     82
     83
     84  bool matrix::equal(const matrix& other, const double d) const
     85  {
     86    if (columns()!=other.columns() || rows()!=other.rows())
     87      return false;
     88    for (size_t i=0; i<rows(); i++)
     89      for (size_t j=0; j<columns(); j++)
     90        // The two last condition checks are needed for NaN detection
     91        if (fabs( (*this)(i,j)-other(i,j) ) > d ||
     92            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
     93          return false;
     94    return true;
     95  }
     96
     97
     98
     99  gsl_matrix* matrix::create_gsl_matrix_copy(void) const
     100  {
     101    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
     102    gsl_matrix_memcpy(m,m_);  // Jari, a GSL return value is ignored here
     103    return m;
     104  }
     105
     106
     107
     108  // Jari, checkout GSL transpose support in GSL manual 8.4.9
     109  void matrix::transpose(void)
     110  {
     111    if (columns()==rows())
     112      gsl_matrix_transpose(m_);
     113    else {
     114      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
     115      gsl_matrix_transpose_memcpy(transposed,m_);
     116      gsl_matrix_free(m_);
     117      m_=transposed;
    91118    }
    92119  }
    … …  
    94121
    95122
    96   gsl_matrix* matrix::gsl_matrix_copy(void) const
    97   {
    98     gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
    99     gsl_matrix_memcpy(m,m_);
    100     return m;
    101   }
    102 
    103   bool matrix::equal(const matrix& other, const double d) const
    104   {
    105     if (columns()!=other.columns() || rows()!=other.rows())
    106       return false;
    107     for (size_t i=0; i<rows(); i++)
    108       for (size_t j=0; j<columns(); j++)
    109         // two last cond. are to check for nans
    110         if (fabs( (*this)(i,j)-other(i,j) ) > d ||
    111             (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
    112           return false;
    113 
    114     return true;
    115   }
    116 
    117   vector matrix::TEMP_col_return(size_t column) const
    118   {
    119     vector vec(rows());
    120     for (size_t i=0; i<rows(); ++i)
    121       vec(i)=operator()(i,column);
    122     return vec;
    123   }
    124 
    125 
    126 
    127   vector matrix::operator[](size_t row) const
    128   {
    129     vector vec(columns());
    130     for (size_t i=0; i<columns(); ++i)
    131       vec(i)=operator()(row,i);
    132     return vec;
    133   }
    134 
    135 
    136 
    137   matrix matrix::operator*( const matrix &other ) const
    138   {
    139     matrix res(rows(), other.columns());
    140     gsl_linalg_matmult(m_,other.m_,res.m_);
    141     return res;
    142   }
    143 
    144 
    145 
    146   vector matrix::operator*( const vector &other ) const
    147   {
    148     gsl_vector* gslvec=gsl_vector_alloc(rows());
    149     gsl_blas_dgemv(CblasNoTrans, 1.0, m_, other.gsl_vector_pointer(), 0.0,
    150                    gslvec );
    151     vector res(gslvec);
    152     return res;
    153   }
    154 
    155 
    156 
    157   matrix matrix::operator+( const matrix &other ) const
    158   {
    159     matrix res( *this );
    160     gsl_matrix_add( res.m_, other.m_ );
    161     return res;
    162   }
    163 
    164 
    165 
    166   matrix matrix::operator-( const matrix &other ) const
    167   {
    168     matrix res( *this );
    169     gsl_matrix_sub( res.m_, other.m_ );
    170     return res;
    171   }
    172 
    173 
    174  
    175   matrix& matrix::operator=( const matrix& other )
     123  const matrix& matrix::operator=( const matrix& other )
    176124  {
    177125    if ( this != &other ) {
    178       if ( m_ )
    179         gsl_matrix_free( m_ );
    180       m_ = other.gsl_matrix_copy();
    181     }
     126      if (view_)
     127        delete view_;
     128      else if (m_)
     129        gsl_matrix_free(m_);
     130      m_ = other.create_gsl_matrix_copy();
     131    }
    182132    return *this;
    183133  }
    … …  
    185135
    186136
    187   matrix matrix::transpose(void) const
     137  const matrix& matrix::operator*=(const matrix& other)
    188138  {
    189     matrix res(columns(),rows());
    190     gsl_matrix_transpose_memcpy(res.m_,m_);
    191     return res;
     139    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
     140    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
     141    gsl_matrix_free(m_);
     142    m_=result;
     143    return *this;
    192144  }
    193145
    … …  
    196148  std::ostream& operator<<(std::ostream& s, const matrix& m)
    197149  {
    198     using namespace std;
    199     s.setf( ios::dec );
     150    s.setf(std::ios::dec);
    200151    s.precision(12);
    201     for(size_t i=0, j=0; i<m.rows(); i++) 
     152    for(size_t i=0, j=0; i<m.rows(); i++)
    202153      for (j=0; j<m.columns(); j++) {
    203154        s << m(i,j);
    … …  
    207158          s << "\n";
    208159      }
    209     return s;
     160    return s;
    210161  }
    211162

  • trunk/lib/gslapi/matrix.h

    r390 r420  
    11// $Id$
    22
    3 #ifndef _theplu_gslapi_matrix_ 
    4 #define _theplu_gslapi_matrix_ 
     3#ifndef _theplu_gslapi_matrix_
     4#define _theplu_gslapi_matrix_
    55
    66#include <c++_tools/gslapi/vector.h>
     7#include <c++_tools/utility/Exception.h>
    78
    89#include <gsl/gsl_matrix.h>
    910#include <iostream>
    10 
    1111
    1212namespace theplu {
    … …  
    2020  /// well in the future.
    2121  ///
    22   class matrix
    23   {
    24   public:
     22  /// \par[File streams] Reading and writing vectors to file streams
     23  /// are of course supported. These are implemented without using GSL
     24  /// functionality, and thus binary read and write to streams are not
     25  /// supported.
     26  ///
     27  /// \par[Matrix views] GSL matrix views are supported and these are
     28  /// disguised as ordinary gslapi::matrix'es. A support function is
     29  /// added, gslapi::matrix::isview(), that can be used to check if a
     30  /// matrix object is a view. Note that view matricess do not own the
     31  /// undelying data, and a view is not valid if the matrix owning the
     32  /// data is deallocated.
     33  ///
     34  /// @note All GSL matrix related functions are not implement but
     35  /// most functionality defined for GSL matrices can be achieved with
     36  /// this interface class. Most notable GSL functionality not
     37  /// supported are no binary file support and views on arrays,
     38  /// gslapi::vectors, gsl_vectors, diagonals, subdiagonals, and
     39  /// superdiagonals. If there is an interest from the user community,
     40  /// the omitted functionality could be included.
     41  ///
     42  /// @todo Maybe it would be smart to create temporary objects need
     43  /// for BLAS in constructors. As it is now temporary objects are
     44  /// called before BLAS functionality iss used, cf. const matrix&
     45  /// matrix::operator*=(const matrix& other)
     46  ///
     47  class matrix
     48  {
     49  public:
    2550
    2651    ///
    2752    /// The default constructor.
    28     ///
    29     matrix(void);
     53    ///
     54    /// This contructor does not initialize underlying (essential)
     55    /// structures.
     56    ///
     57    inline matrix(void) : m_(NULL), view_(NULL) {}
    3058
    3159    ///
    3260    /// Constructor. Allocates memory space for \a r times \a c
    3361    /// elements, and sets all elements to \a init_value.
    34     ///
    35     matrix(const size_t& r, const size_t& c, const double init_value=0);
    36  
    37     ///
     62    ///
     63    inline matrix(const size_t& r, const size_t& c, double init_value=0)
     64      : view_(NULL) { m_ = gsl_matrix_alloc(r,c); set_all(init_value); }
     65
     66    ///
    3867    /// The copy constructor.
    39     ///
    40     matrix(const matrix&);
     68    ///
     69    /// @note If the object to be copied is a matrix view, the values
     70    /// of the view will be copied, i.e. the view is not copied.
     71    ///
     72    inline matrix(const matrix& o) : view_(NULL) {m_=o.create_gsl_matrix_copy();}
     73
     74    ///
     75    /// The matrix view constructor.
     76    ///
     77    /// Create a view of matrix \a m, with starting row \a offset_row,
     78    /// starting column \a offset_column, row size \a n_row, and
     79    /// column size \a n_column.
     80    ///
     81    /// A matrix view can be used as any matrix with the difference
     82    /// that changes made to the view will also change the object that
     83    /// is viewed. Also, using the copy constructor will create a new
     84    /// matrix object that is a copy of whatever is viewed. If a copy
     85    /// of the view is needed then you should use this constructor to
     86    /// obtain a copy.
     87    ///
     88    /// @note If the object viewed by the view goes out of scope or is
     89    /// deleted, the view becomes invalid and the result of further
     90    /// use is undefined.
     91    ///
     92    matrix(matrix& m, size_t offset_row, size_t offset_column, size_t n_row,
     93           size_t n_column);
    4194
    4295    ///
    … …  
    4699    /// column elements in a row must be separated with white space
    47100    /// characters except the new line (\\n) character. Rows, and
    48     /// columns, are read consequently and only rectangular matrices
     101    /// columns, are read consecutively and only rectangular matrices
    49102    /// are supported. Empty lines are ignored. End of input is at end
    50103    /// of file marker.
    51104    ///
    52     explicit matrix(std::istream &);
    53 
    54     ///
     105    explicit matrix(std::istream &) throw (utility::IO_error,std::exception);
     106
     107    ///
    55108    /// The destructor.
    56109    ///
    57     ~matrix(void);
    58 
    59     ///
     110    ~matrix(void);
     111
     112    ///
     113    /// Elementwise addition of the elements of matrix \a b to the
     114    /// elements of the calling matrix ,\f$a_{ij} = a_{ij} + b_{ij} \;
     115    /// \forall i,j\f$. The result is stored into the calling matrix.
     116    ///
     117    /// @return Whatever GSL returns.
     118    ///
     119    inline int add(const matrix& b) { return gsl_matrix_add(m_,b.m_); }
     120
     121    ///
     122    /// Add the scalar value \a d to the elements of the calling
     123    /// matrix, \f$a_{ij} = a_{ij} + d \; \forall i,j\f$. The result
     124    /// is stored into the calling matrix.
     125    ///
     126    /// @return Whatever GSL returns.
     127    ///
     128    inline int
     129    add_constant(const double d) { return gsl_matrix_add_constant(m_,d); }
     130
     131    ///
    60132    /// @return The number of columns in the matrix.
    61133    ///
    62     inline size_t columns(void) const { return m_->size2; }
    63 
    64     ///
    65     /// Swap column \a i with \a j.
    66     ///
    67     inline void column_swap(const size_t& i,const size_t& j)
    68       { gsl_matrix_swap_columns(m_,i,j); }
    69 
    70     ///
    71     /// @return true if each element deviates less or equal than \a
    72     /// d. If any matrix contain a nan false is always returned.
    73     ///
    74     bool equal(const matrix&, const double d=0) const;
     134    inline size_t columns(void) const { return m_->size2; }
     135
     136    ///
     137    /// Elementwise division of the elemnts of the calling matrix by
     138    /// the elements of matrix \a b, \f$a_{ij} = a_{ij} / b_{ij} \;
     139    /// \forall i,j\f$. The result is stored into the calling matrix.
     140    ///
     141    /// @return Whatever GSL returns.
     142    ///
     143    inline int
     144    div_elements(const matrix& b) { return gsl_matrix_div_elements(m_,b.m_); }
     145
     146    ///
     147    /// Check whether matrices are equal within a user defined
     148    /// precision, set by \a precision.
     149    ///
     150    /// @return True if each element deviates less or equal than \a
     151    /// d. If any matrix contain a NaN, false is always returned.
     152    ///
     153    bool equal(const matrix&, const double precision=0) const;
     154
     155    ///
     156    /// @return A const pointer to the internal GSL matrix.
     157    ///
     158    inline const gsl_matrix* gsl_matrix_p(void) const { return m_; }
     159
     160    ///
     161    /// @return A pointer to the internal GSL matrix.
     162    ///
     163    // Jari, is this needed?
     164    inline gsl_matrix* gsl_matrix_p(void) { return m_; };
     165
     166    ///
     167    /// @return True if all elements in the matrix is zero, false
     168    /// othwerwise;
     169    ///
     170    inline bool isnull(void) const { return gsl_matrix_isnull(m_); }
     171
     172    ///
     173    /// Check if the matrix object is a view (sub-matrix) to another
     174    /// matrix.
     175    ///
     176    /// @return True if the object is a view, false othwerwise.
     177    ///
     178    inline bool isview(void) const { return view_; }
     179
     180    ///
     181    /// @return The maximum value of the matrix.
     182    ///
     183    inline double max(void) const { return gsl_matrix_max(m_); }
     184
     185    ///
     186    /// @return The index to the element with the maximum value of the
     187    /// matrix. The lowest index has precedence (searching in
     188    /// row-major order).
     189    ///
     190    inline std::pair<size_t,size_t> max_index(void) const;
     191
     192    ///
     193    /// @return The minimum value of the matrix.
     194    ///
     195    inline double min(void) const { return gsl_matrix_min(m_); }
     196
     197    ///
     198    /// @return The index to the element with the minimum value of the
     199    /// matrix. The lowest index has precedence (searching in
     200    /// row-major order).
     201    ///
     202    inline std::pair<size_t,size_t> min_index(void) const;
     203
     204    ///
     205    /// @return The indecies to the element with the minimum and
     206    /// maximum values of the matrix, respectively. The lowest index
     207    /// has precedence (searching in row-major order).
     208    ///
     209    inline void minmax_index(std::pair<size_t,size_t>& min,
     210                             std::pair<size_t,size_t>& max) const
     211      { gsl_matrix_minmax_index(m_,&min.first,&min.second,
     212                                &max.first,&max.second); }
     213
     214    ///
     215    /// @return The minimum and maximum values of the matrix, as the
     216    /// \a first and \a second member of the returned \a pair,
     217    /// respectively.
     218    ///
     219    std::pair<double,double> minmax(void) const;
     220
     221    ///
     222    /// Multiply the elements of matrix \a b with the elements of the
     223    /// calling matrix ,\f$a_{ij} = a_{ij} * b_{ij} \; \forall
     224    /// i,j\f$. The result is stored into the calling matrix.
     225    ///
     226    /// @return Whatever GSL returns.
     227    ///
     228    inline int
     229    mul_elements(const matrix& b) { return gsl_matrix_mul_elements(m_,b.m_); }
     230
     231    ///
     232    /// @return The number of rows in the matrix.
     233    ///
     234    inline size_t rows(void) const { return m_->size1; }
     235
     236    ///
     237    /// Multiply the elements of the calling matrix with a scalar \a
     238    /// d, \f$a_{ij} = d * a_{ij} \; \forall i,j\f$. The result is
     239    /// stored into the calling matrix.
     240    ///
     241    /// @return Whatever GSL returns.
     242    ///
     243    inline int scale(const double d) { return gsl_matrix_scale(m_,d); }
     244
     245    ///
     246    /// Set element values to values in \a mat. This function is
     247    /// needed for assignment of viewed elements.
     248    ///
     249    /// @return Whatever GSL returns.
     250    ///
     251    /// @note No check on size is done.
     252    ///
     253    /// @see const matrix& operator=(const matrix&)
     254    ///
     255    inline int set(const matrix& mat) { return gsl_matrix_memcpy(m_,mat.m_); }
     256
     257    ///
     258    /// Set all elements to \a value.
     259    ///
     260    inline void set_all(const double value) { gsl_matrix_set_all(m_,value); }
     261
     262    ///
     263    /// Set \a column values to values in \a vec.
     264    ///
     265    /// @return Whatever GSL returns.
     266    ///
     267    /// @note No check on size is done.
     268    ///
     269    inline int set_column(const size_t column, const vector& vec)
     270      { return gsl_matrix_set_col(m_, column, vec.gsl_vector_p()); }
     271
     272    ///
     273    /// Set \a row values to values in \a vec.
     274    ///
     275    /// @return Whatever GSL returns.
     276    ///
     277    /// @note No check on size is done.
     278    ///
     279    inline int set_row(const size_t row, const vector& vec)
     280      { return gsl_matrix_set_row(m_, row, vec.gsl_vector_p()); }
     281
     282    ///
     283    /// Subtract the elements of matrix \a b from the elements of the
     284    /// calling matrix ,\f$a_{ij} = a_{ij} - b_{ij} \; \forall
     285    /// i,j\f$. The result is stored into the calling matrix.
     286    ///
     287    /// @return Whatever GSL returns.
     288    ///
     289    inline int sub(const matrix& b) { return gsl_matrix_sub(m_,b.m_); }
     290
     291    ///
     292    /// Exchange the elements of the this and \a other by copying. The
     293    /// two matrices must have the same size.
     294    ///
     295    /// @return Whatever GSL returns.
     296    ///
     297    inline int swap(matrix& other) { return gsl_matrix_swap(m_,other.m_); }
     298
     299    ///
     300    /// Swap columns \a i and \a j.
     301    ///
     302    inline int swap_columns(const size_t i,const size_t j)
     303      { return gsl_matrix_swap_columns(m_,i,j); }
     304
     305    ///
     306    /// Swap row \a i and column \a j.
     307    ///
     308    inline int swap_rowcol(const size_t i,const size_t j)
     309      { return gsl_matrix_swap_rowcol(m_,i,j); }
     310
     311    ///
     312    /// Swap rows \a i and \a j.
     313    ///
     314    inline int swap_rows(const size_t i, const size_t j)
     315      { return gsl_matrix_swap_rows(m_,i,j); }
     316
     317    ///
     318    /// Transpose the matrix.
     319    ///
     320    void transpose(void);
     321
     322    ///
     323    /// @return Reference to the element position (\a row, \a column).
     324    ///
     325    inline double& operator()(size_t row,size_t column)
     326      { return (*gsl_matrix_ptr(m_,row,column)); }
     327
     328    ///
     329    /// @return Const reference to the element position (\a row, \a
     330    /// column).
     331    ///
     332    inline const double& operator()(size_t row,size_t column) const
     333      { return (*gsl_matrix_const_ptr(m_,row,column)); }
     334
     335    ///
     336    /// Matrix-vector multiplication.
     337    ///
     338    /// @return The resulting vector.
     339    ///
     340    // Jari, where should this go?
     341    const vector operator*(const vector&) const;
     342
     343    ///
     344    /// Comparison operator.
     345    ///
     346    /// @return True if all elements are equal otherwise False.
     347    ///
     348    /// @see equal
     349    ///
     350    inline bool operator==(const matrix& other) const { return equal(other); }
     351
     352    ///
     353    /// Comparison operator.
     354    ///
     355    /// @return False if all elements are equal otherwise True.
     356    ///
     357    /// @see equal
     358    ///
     359    inline bool operator!=(const matrix& other) const { return !equal(other); }
     360
     361    ///
     362    /// The assignment operator. There is no requirements on
     363    /// dimensions, i.e. the matrix is remapped in memory if
     364    /// necessary. This implies that in general views cannot be
     365    /// assigned using this operator. Views will be mutated into
     366    /// normal matrices. The only exception to this behaviour on views
     367    /// is when self-assignemnt is done, since self-assignment is
     368    /// ignored.
     369    ///
     370    /// @return A const reference to the resulting matrix.
     371    ///
     372    /// @see int set(const matrix&)
     373    ///
     374    const matrix& operator=(const matrix& other);
     375
     376    ///
     377    /// Add and assign operator.
     378    ///
     379    inline const matrix& operator+=(const matrix& m) { add(m); return *this; }
     380
     381    ///
     382    /// Add and assign operator.
     383    ///
     384    inline const matrix&
     385    operator+=(const double d) { add_constant(d); return *this; }
     386
     387    ///
     388    /// Subtract and assign operator.
     389    ///
     390    inline const matrix& operator-=(const matrix& m) { sub(m); return *this; }
     391
     392    ///
     393    /// Multiply and assigment operator.
     394    ///
     395    /// @return Const reference to the resulting matrix.
     396    ///
     397    const matrix& operator*=(const matrix&);
     398
     399    ///
     400    /// Multiply and assign operator.
     401    ///
     402    inline const matrix& operator*=(const double d) { scale(d); return *this; }
     403
     404
     405  private:
    75406
    76407    ///
    … …  
    82413    /// @return A pointer to a copy of the internal GSL matrix.
    83414    ///
    84     gsl_matrix* gsl_matrix_copy(void) const;
    85 
    86     ///
    87     /// @return A const pointer to the internal GSL matrix.
    88     ///
    89     inline const gsl_matrix* gsl_matrix_pointer(void) const { return m_; }
    90 
    91     ///
    92     /// @return A pointer to the internal GSL matrix.
    93     ///
    94     inline gsl_matrix* gsl_matrix_pointer(void) { return m_; };
    95 
    96     ///
    97     /// This function multiplies all elements between two matrices and
    98     /// the result is stored back into the calling object, \f$a_{ij} =
    99     /// a_{ij} * b_{ij} \; \forall i,j\f$.
    100     ///
    101     inline void mul_elements(const matrix& b)
    102       { gsl_matrix_mul_elements(m_,b.m_); }
    103 
    104     ///
    105     /// Swap row \a i with \a j.
    106     ///
    107     inline void row_swap( const size_t& i, const size_t& j)
    108       { gsl_matrix_swap_rows(m_,i,j); }
    109 
    110     ///
    111     /// @return The number of rows in the matrix.
    112     ///
    113     inline size_t rows(void) const { return m_->size1; }
    114 
    115     ///
    116     /// Set all elements to \a value.
    117     ///
    118     inline void set_all(const double value) { gsl_matrix_set_all(m_,value); }
    119 
    120     ///
    121     /// Set column \a i to \a vec.
    122     ///
    123     inline void set_column(const size_t i, const vector vec)
    124     {gsl_matrix_set_col(m_, i, vec.gsl_vector_pointer());}
    125 
    126     ///
    127     /// Set row \a i to \a vec.
    128     ///
    129     inline void set_row(const size_t i, const vector vec)
    130     {gsl_matrix_set_row(m_, i, vec.gsl_vector_pointer());}
    131 
    132     ///
    133     /// Set column \a i to \a vec
    134     ///
    135     //void matrix set_column(const size_t i, const vector vec) const;
    136 
    137     ///
    138     /// Transpose the matrix.
    139     ///
    140     /// @return Transpose of the matrix.
    141     ///
    142     matrix transpose(void) const;
    143    
    144     ///
    145     /// @return Reference to the element position (\a row, \a column).
    146     ///
    147     inline double& operator()(size_t row,size_t column)
    148       { return (*gsl_matrix_ptr(m_,row,column)); }
    149 
    150     ///
    151     /// @return Const reference to the element position (\a row, \a
    152     /// column).
    153     ///
    154     inline const double& operator()(size_t row,size_t column) const
    155       { return (*gsl_matrix_const_ptr(m_,row,column)); }
    156 
    157     ///
    158     /// @todo to be properly implemented. Should return a reference to a vector
    159     /// object. Underlying GSL supports only GSL_vector views, thus
    160     /// some redesign of the vector class is needed.
    161     ///
    162     /// inline vector& operator[](size_t i);
    163     /// So for now, stupid copying is done, and actually a matrix
    164     /// row is returned using this function.
    165     vector operator[](size_t row) const;
    166 
    167     ///
    168     /// @todo to be properly implemented. Should return a reference to
    169     /// a vector object (matrix row). Underlying GSL supports only GSL_vector
    170     /// views, thus some redesign of the vector class is needed.
    171     ///
    172     /// inline const vector& operator[](size_t i) const;
    173     ///
    174     /// So for now, stupid copying is done, and actually a matrix
    175     /// column is returned using this function.
    176     vector TEMP_col_return(size_t column) const;
    177 
    178     ///
    179     /// Matrix multiplication.
    180     ///
    181     /// @return The resulting matrix.
    182     ///
    183     matrix operator*(const matrix&) const;
    184 
    185     ///
    186     /// Matrix-vector multiplication.
    187     ///
    188     /// @return The resulting vector.
    189     ///
    190     vector operator*(const vector&) const;
    191 
    192     ///
    193     /// Matrix addition.
    194     ///
    195     /// @return The resulting matrix.
    196     ///
    197     matrix operator+(const matrix&) const;
    198 
    199     ///
    200     /// @todo Matrix addition.
    201     ///
    202     /// @return The resulting matrix.
    203     ///
    204     matrix operator+=(const matrix&);
    205 
    206     ///
    207     /// Matrix subtraction.
    208     ///
    209     /// @return The resulting matrix.
    210     ///
    211     matrix operator-(const matrix&) const;
    212 
    213     ///
    214     /// @todo Matrix subtraction.
    215     ///
    216     /// @return The resulting matrix.
    217     ///
    218     matrix operator-=(const matrix&) const;
    219 
    220     ///
    221     /// Comparison operator.
    222     ///
    223     /// @return True if all elements are equal otherwise False.
    224     ///
    225     /// @see equal
    226     ///
    227     inline bool operator==( const matrix& other ) const { return equal(other);}
    228 
    229     ///
    230     /// Comparison operator.
    231     ///
    232     /// @return False if all elements are equal otherwise True.
    233     ///
    234     /// @see equal
    235     ///
    236     inline bool operator!=( const matrix& other ) const {return !equal(other);}
    237 
    238     ///
    239     /// The assignment operator. There is no requirements on
    240     /// dimensions.
    241     ///
    242     matrix& operator=(const matrix& other);
    243 
    244     ///
    245     /// Subtract and assign operator.
    246     ///
    247     inline matrix& operator-=(const matrix& m)
    248       { gsl_matrix_sub(m_,m.m_); return *this; }
    249 
    250     ///
    251     /// Multiply and assign operator.
    252     ///
    253     inline matrix& operator*=(const double d)
    254       { gsl_matrix_scale(m_,d); return *this; }
    255 
    256 
    257   private:
    258 
    259     gsl_matrix* m_;
    260 
     415    gsl_matrix* create_gsl_matrix_copy(void) const;
     416
     417    gsl_matrix* m_;
     418    gsl_matrix_view* view_;
    261419  };
    262420
    263421  ///
    264422  /// The output operator for the matrix class.
    265   ///
    266   std::ostream& operator<< (std::ostream& s, const matrix&);
    267 
     423  ///
     424  std::ostream& operator<< (std::ostream& s, const matrix&);
    268425
    269426
    270427}} // of namespace gslapi and namespace theplu
    271428
    272 #endif 
     429#endif

  • trunk/lib/gslapi/vector.cc

    r357 r420  
    22
    33#include <c++_tools/gslapi/vector.h>
     4#include <c++_tools/gslapi/matrix.h>
    45#include <c++_tools/utility/stl_utility.h>
    56
    6 #include <iostream>
    77#include <sstream>
    8 #include <string>
    98#include <vector>
    109#include <utility>
    … …  
    1615
    1716
    18   vector::vector(void)
    19     : v_(NULL), view_(NULL)
    20   {
    21   }
    22 
    23 
    24 
    25   vector::vector(const size_t n,const double init_value)
    26     : view_(NULL)
    27   {
    28     v_ = gsl_vector_alloc(n);
    29     set_all(init_value);
    30   }
    31 
    32 
    33 
    34   vector::vector(const vector& other)
    35     : view_(NULL)
    36   {
    37     v_ = other.TEMP_gsl_vector_copy();
    38   }
    39 
    40 
    41 
    4217  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
    43   {
     18    : const_view_(NULL)
     19  {
     20    // Jari, exception handling needed here. Failure in setting up a
     21    // proper gsl_vector_view is communicated by NULL pointer in the
     22    // view structure (cf. GSL manual). How about GSL error state?
    4423    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
    4524                                                                 stride,n));
    … …  
    4928
    5029
    51   vector::vector(gsl_vector* vec)
    52     : v_(vec), view_(NULL)
    53   {
     30  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
     31    : view_(NULL)
     32  {
     33    // Jari, exception handling needed here. Failure in setting up a
     34    // proper gsl_vector_view is communicated by NULL pointer in the
     35    // view structure (cf. GSL manual). How about GSL error state?
     36    const_view_ = new gsl_vector_const_view(
     37                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
     38    v_ = const_cast<gsl_vector*>(&(const_view_->vector));
     39  }
     40
     41
     42
     43  vector::vector(matrix& m, size_t i, bool row)
     44    : const_view_(NULL)
     45  {
     46    view_=new gsl_vector_view(row ?
     47                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
     48                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
     49    v_ = &(view_->vector);
     50  }
     51
     52
     53
     54  vector::vector(const matrix& m, size_t i, bool row)
     55    : view_(NULL)
     56  {
     57    const_view_ = new gsl_vector_const_view(row ?
     58                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
     59                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
     60    v_ = const_cast<gsl_vector*>(&(const_view_->vector));
    5461  }
    5562
    … …  
    5764
    5865  vector::vector(std::istream& is) throw (utility::IO_error,std::exception)
    59     : view_(NULL)
     66    : view_(NULL), const_view_(NULL)
    6067  {
    6168    // read the data file and store in stl vectors (dynamically
    … …  
    7986          << ") and columns (" << nof_columns
    8087          << ").\n    Expected a row or a column vector.";
    81         std::string m=s.str();
    82         throw utility::IO_error(m);
     88        throw utility::IO_error(s.str());
    8389      }
    84       else if(v.size()!=nof_columns) {
     90      else if (v.size()!=nof_columns) {
    8591        std::ostringstream s;
    8692        s << "vector::vector(std::istream&) data file error:\n"
    8793          << "    Line " << (nof_rows+1) << " has " << v.size()
    8894          << " columns; expected " << nof_columns << " column.";
    89         std::string m=s.str();
    90         throw utility::IO_error(m);
     95        throw utility::IO_error(s.str());
    9196      }
    9297      data_matrix.push_back(v);
    … …  
    109114    if (view_)
    110115      delete view_;
     116    else if (const_view_)
     117      delete const_view_;
    111118    else if (v_)
    112119      gsl_vector_free(v_);
    … …  
    116123
    117124
    118   gsl_vector* vector::TEMP_gsl_vector_copy(void) const
     125  gsl_vector* vector::create_gsl_vector_copy(void) const
    119126  {
    120127    gsl_vector* vec = gsl_vector_alloc(size());
    … …  
    159166    for ( size_t i = 0; i < size(); ++i )
    160167      res += other[i] * (*this)[i];
    161     return res;
    162   }
    163 
    164 
    165 
    166   vector vector::operator+(const vector &other) const
    167   {
    168     vector res(*this);
    169     gsl_vector_add(res.v_,other.v_);
    170     return res;
    171   }
    172 
    173 
    174 
    175   vector vector::operator-( const vector &other ) const
    176   {
    177     vector res( *this );
    178     gsl_vector_sub(res.v_,other.v_);
    179168    return res;
    180169  }
    … …  
    197186  {
    198187    if( this != &other ) {
    199       if ( v_ )
     188      if (view_)
     189        delete view_;
     190      else if (const_view_)
     191        delete const_view_;
     192      else if ( v_ )
    200193        gsl_vector_free( v_ );
    201       v_ = other.TEMP_gsl_vector_copy();
     194      v_ = other.create_gsl_vector_copy();
    202195    }
    203196    return *this;
    … …  
    208201  std::ostream& operator<<(std::ostream& s, const vector& a)
    209202  {
    210     s.setf(std::ios::fixed);
     203    s.setf(std::ios::dec);
    211204    s.precision(12);
    212205    for (size_t j = 0; j < a.size(); ++j) {

  • trunk/lib/gslapi/vector.h

    r357 r420  
    11// $Id$
    22
    3 #ifndef _theplu_gslapi_vector_ 
     3#ifndef _theplu_gslapi_vector_
    44#define _theplu_gslapi_vector_
     5
     6#include <c++_tools/utility/Exception.h>
    57
    68#include <iostream>
    … …  
    911
    1012#include <gsl/gsl_vector.h>
    11 #include <c++_tools/utility/Exception.h>
    1213
    1314namespace theplu {
    1415namespace gslapi {
     16
     17  class matrix;
     18
    1519  ///
    1620  /// This is the C++ tools interface to GSL vector. 'double' is the
    … …  
    2529  /// \par[Vector views] GSL vector views are supported and these are
    2630  /// disguised as ordinary gslapi::vectors. A support function is
    27   /// added, gslapi::vector::is_view(), that can be used to check if
    28   /// the vector object is a view. Note that view vectors do not own
    29   /// the undelying data, and is not valid if the original vector is
    30   /// deallocated.
     31  /// added, gslapi::vector::isview(), that can be used to check if a
     32  /// vector object is a view. Note that view vectors do not own the
     33  /// undelying data, and a view is not valid if the vector owning the
     34  /// data is deallocated.
     35  ///
     36  /// @note Missing support to underlying GSL vector features can in
     37  /// principle be included to this class if requested by the user
     38  /// community.
    3139  ///
    3240  class vector
    … …  
    3745    /// The default constructor.
    3846    ///
    39     vector(void);
     47    inline vector(void) : v_(NULL), view_(NULL), const_view_(NULL) {}
    4048
    4149    ///
    … …  
    4351    /// sets all elements to \a init_value.
    4452    ///
    45     vector(const size_t n, const double init_value=0);
     53    inline vector(const size_t n,const double init_value=0)
     54      : view_(NULL), const_view_(NULL)
     55      { v_ = gsl_vector_alloc(n); set_all(init_value); }
    4656
    4757    ///
    … …  
    5161    /// of the view will be copied, i.e. the view is not copied.
    5262    ///
    53     vector(const vector&);
    54 
    55     ///
    56     /// The vector view constructor.
     63    inline vector(const vector& other) : view_(NULL), const_view_(NULL)
     64      { v_ = other.create_gsl_vector_copy(); }
     65
     66    ///
     67    /// Vector view constructor.
    5768    ///
    5869    /// Create a view of vector \a v, with starting index \a offset,
    … …  
    6374    /// is viewed. Also, using the copy constructor will create a new
    6475    /// vector object that is a copy of whatever is viewed. If a copy
    65     /// of the view is needed the you should use this consturctor to
    66     /// get one.
     76    /// of the view is needed then you should use this constructor to
     77    /// obtain a copy.
    6778    ///
    6879    /// @note If the object viewed by the view goes out of scope or is
    69     /// deleted, the view becomes invalid.
    70     ///
    71     vector(vector& v, size_t offset, size_t n, size_t stride=1);
    72 
    73     ///
    74     /// Constructor that imports a GSL vector. The GSL object is owned
    75     /// by the created object.
    76     ///
    77     vector(gsl_vector*);
     80    /// deleted, the view becomes invalid and the result of further
     81    /// use is undefined.
     82    ///
     83    vector(vector& v, size_t offset, size_t n, size_t stride=1);
     84
     85    ///
     86    /// Vector const view constructor.
     87    ///
     88    /// Create a view of vector \a v, with starting index \a offset,
     89    /// size \a n, and an optional \a stride.
     90    ///
     91    /// A vector view can be used as any const vector. Using the copy
     92    /// constructor will create a new vector object that is a copy of
     93    /// whatever is viewed. If a copy of the view is needed then you
     94    /// should use this constructor to obtain a copy.
     95    ///
     96    /// @note If the object viewed by the view goes out of scope or is
     97    /// deleted, the view becomes invalid and the result of further
     98    /// use is undefined.
     99    ///
     100    vector(const vector& v, size_t offset, size_t n, size_t stride=1);
     101
     102    ///
     103    /// Matrix row/column view constructor.
     104    ///
     105    /// Create a row/column vector view of matrix \a m, pointing at
     106    /// row/column \a i. The parameter \a row is used to set whether
     107    /// the view should be a row or column view. If \a row is set to
     108    /// true, the view will be a row view (default behaviour), and,
     109    /// naturally, a column view otherwise.
     110    ///
     111    /// A vector view can be used as any vector with the difference
     112    /// that changes made to the view will also change the object that
     113    /// is viewed. Also, using the copy constructor will create a new
     114    /// vector object that is a copy of whatever is viewed. If a copy
     115    /// of the view is needed then you should use the vector view
     116    /// constructor to obtain a copy.
     117    ///
     118    /// @note If the object viewed by the view goes out of scope or is
     119    /// deleted, the view becomes invalid and the result of further
     120    /// use is undefined.
     121    ///
     122    vector(matrix& m, size_t i, bool row=true);
     123
     124    ///
     125    /// Matrix row/column const view constructor.
     126    ///
     127    /// Create a row/column vector view of matrix \a m, pointing at
     128    /// row/column \a i. The parameter \a row is used to set whether
     129    /// the view should be a row or column view. If \a row is set to
     130    /// true, the view will be a row view (default behaviour), and,
     131    /// naturally, a column view otherwise.
     132    ///
     133    /// A const vector view can be used as any const vector. Using the
     134    /// copy constructor will create a new vector object that is a
     135    /// copy of whatever is viewed. If a copy of the view is needed
     136    /// then you should use the vector view constructor to obtain a
     137    /// copy.
     138    ///
     139    /// @note If the object viewed by the view goes out of scope or is
     140    /// deleted, the view becomes invalid and the result of further
     141    /// use is undefined.
     142    ///
     143    vector(const matrix& m, size_t i, bool row=true);
    78144
    79145    ///
    … …  
    117183
    118184    ///
     185    /// @return A const pointer to the internal GSL vector,
     186    ///
     187    inline const gsl_vector* gsl_vector_p(void) const { return v_; }
     188
     189    ///
     190    /// @return A pointer to the internal GSL vector,
     191    ///
     192    inline gsl_vector* gsl_vector_p(void) { return v_; }
     193
     194    ///
    119195    /// @return True if all elements in the vector is zero, false
    120196    /// othwerwise;
    … …  
    123199
    124200    ///
    125     /// Check if the vector object is a view (sub vector) to another
     201    /// Check if the vector object is a view (sub-vector) to another
    126202    /// vector.
    127203    ///
    128     /// @return True if the object is a view, false othwerwise;
    129     ///
    130     inline bool isview(void) const { return view_; }
    131 
    132     ///
    133     /// @return A const pointer to the internal GSL vector,
    134     ///
    135     inline const gsl_vector* gsl_vector_pointer(void) const { return v_; }
    136 
    137     ///
    138     /// @return A pointer to the internal GSL vector,
    139     ///
    140     inline gsl_vector* TEMP_gsl_vector_pointer(void) { return v_; }
     204    /// @return True if the object is a view, false othwerwise.
     205    ///
     206    inline bool isview(void) const { return view_ || const_view_; }
    141207
    142208    ///
    … …  
    208274
    209275    ///
     276    /// Set element values to values in \a vec. This function is
     277    /// needed for assignment of viewed elements.
     278    ///
     279    /// @return Whatever GSL returns.
     280    ///
     281    /// @note No check on size is done.
     282    ///
     283    /// @see const vector& operator=(const vector&)
     284    ///
     285    inline int set(const vector& vec) { return gsl_vector_memcpy(v_,vec.v_); }
     286
     287    ///
    210288    /// Set all elements to \a value.
    211289    ///
    … …  
    279357    ///
    280358    // Jari, doxygen group as Accessing vector elements
    281     inline const double& operator()(size_t i) const
    282       { return *gsl_vector_const_ptr(v_,i); }
    283    
     359    inline const double&
     360    operator()(size_t i) const { return *gsl_vector_const_ptr(v_,i); }
     361
    284362    ///
    285363    /// Element access operator.
    … …  
    296374    ///
    297375    // Jari, doxygen group as Accessing vector elements
    298     inline const double& operator[](size_t i) const
    299       { return *gsl_vector_const_ptr(v_,i); }
    300    
     376    inline const double&
     377    operator[](size_t i) const { return *gsl_vector_const_ptr(v_,i); }
     378
    301379    ///
    302380    /// @return The dot product.
    303381    ///
    304382    double operator*( const vector &other ) const;
    305 
    306     ///
    307     /// Vector with scalar multiplication.
    308     ///
    309     /// @note This operator is not implemented!
    310     ///
    311     vector operator*(const double d) const;
    312 
    313     ///
    314     /// Vector addition.
    315     ///
    316     /// @return The resulting vector.
    317     ///
    318     vector operator+( const vector& other ) const;
    319 
    320     ///
    321     /// Vector subtraction.
    322     ///
    323     /// @return The resulting vector.
    324     ///
    325     vector operator-( const vector& other ) const;
    326383
    327384    ///
    … …  
    335392    ///
    336393    /// The assignment operator. There is no requirements on
    337     /// dimensions.
     394    /// dimensions, i.e. the vector is remapped in memory if
     395    /// necessary. This implies that in general views cannot be
     396    /// assigned using this operator. Views will be mutated into
     397    /// normal vectors. The only exception to this behaviour on views
     398    /// is when self-assignemnt is done, since self-assignment is
     399    /// ignored.
    338400    ///
    339401    /// @return A const reference to the resulting vector.
    340402    ///
     403    /// @see int set(const vector&)
     404    ///
    341405    const vector& operator=(const vector&);
    342406
    … …  
    346410    /// @return A const reference to the resulting vector.
    347411    ///
    348     inline const vector& operator+=( const vector& other )
    349     { gsl_vector_add(v_,other.v_); return *this; }
     412    inline const vector&
     413    operator+=(const vector& other) { gsl_vector_add(v_,other.v_); return *this;}
    350414
    351415    ///
    … …  
    354418    /// @return A const reference to the resulting vector.
    355419    ///
    356     inline const vector& operator-=( const vector& other )
    357     { gsl_vector_sub(v_,other.v_); return *this; }
     420    inline const vector&
     421    operator-=(const vector& other) { gsl_vector_sub(v_,other.v_); return *this;}
    358422
    359423    ///
    … …  
    362426    /// @return A const reference to the resulting vector.
    363427    ///
    364     inline const vector& operator*=(const double d)
    365     { gsl_vector_scale(v_,d); return *this; }
    366 
    367 
    368   private:
     428    inline const vector&
     429    operator*=(const double d) { gsl_vector_scale(v_,d); return *this; }
     430
     431
     432  private:
    369433
    370434    ///
    … …  
    376440    /// @return A pointer to a copy of the internal GSL vector.
    377441    ///
    378     gsl_vector* TEMP_gsl_vector_copy(void) const;
    379 
    380     gsl_vector* v_;
     442    gsl_vector* create_gsl_vector_copy(void) const;
     443
     444    gsl_vector* v_;
    381445    gsl_vector_view* view_;
    382   };
     446    gsl_vector_const_view* const_view_;
     447  };
    383448
    384449  ///
    385450  /// The output operator for the vector class.
    386   /// 
     451  ///
    387452  std::ostream& operator<<(std::ostream&, const vector& );
    388453
    … …  
    390455}} // of namespace gslapi and namespace theplu
    391456
    392 #endif 
     457#endif

  • trunk/lib/statistics/Averager.cc

    r295 r420  
    1111 
    1212
    13   Averager::Averager(void)
    14     : n_(0), x_(0), xx_(0)
    15   {
    16   }
    17 
    18   Averager::Averager(const double x,const double xx,const long n)
    19     : n_(n), x_(x), xx_(xx)
    20   {
    21   }
    22 
    23   Averager::Averager(const Averager& a)
    24     : n_(a.n_), x_(a.x_), xx_(a.xx_)
    25   {
    26   }
    27 
    2813  const Averager& Averager::operator+=(const Averager& a)
    2914  {
    … …  
    3419  }
    3520
     21
    3622}} // of namespace statistics and namespace theplu

  • trunk/lib/statistics/Averager.h

    r349 r420  
    2626    /// Default constructor
    2727    ///
    28     Averager(void);
     28    inline Averager(void) : n_(0), x_(0), xx_(0) {}
    2929   
    3030    ///
    … …  
    3232    /// number of samples \a n.
    3333    ///
    34     Averager(const double x, const double xx, const long n);
     34    inline Averager(const double x,const double xx,const long n)
     35      : n_(n), x_(x), xx_(xx) {}
    3536
    3637    ///
    3738    /// Copy constructor
    3839    ///
    39     Averager(const Averager&);
    40    
     40    inline Averager(const Averager& a) : n_(a.n_), x_(a.x_), xx_(a.xx_) {}
     41
    4142    ///
    4243    /// Adding \a n (default=1) number of data point(s) with value \a d.
    4344    ///
    44     inline void add(const double d,const long n=1)
    45     {n_+=n; x_+=n*d; xx_+=n*d*d;}
     45    inline void add(const double d,const long n=1) { n_+=n; x_+=n*d; xx_+=n*d*d;}
    4646
    4747    ///

  • trunk/lib/statistics/AveragerPair.cc

    r295 r420  
    1010 
    1111
    12   AveragerPair::AveragerPair(void)
    13     : x_(Averager()), y_(Averager()), xy_(0.0)
    14   {
    15   }
    16 
    17   AveragerPair::AveragerPair(const double x, const double xx, const double y,
    18                      const double yy, const double xy, const long n)
    19     : x_(Averager(x,xx,n)), y_(Averager(y,yy,n)), xy_(xy)
    20   {
    21   }
    22 
    23   AveragerPair::AveragerPair(const AveragerPair& a)
    24     : x_(a.x_averager()), y_(a.y_averager()), xy_(a.sum_xy())
    25   {
    26   }
    27 
    2812  const AveragerPair& AveragerPair::operator+=(const AveragerPair& a)
    2913  {
    … …  
    3418  }
    3519
     20
    3621}} // of namespace statistics and namespace theplu

  • trunk/lib/statistics/AveragerPair.h

    r355 r420  
    2525    /// Default constructor
    2626    ///
    27     AveragerPair(void);
     27    inline AveragerPair(void) : x_(Averager()), y_(Averager()), xy_(0.0) {}
    2828   
    2929    ///
    … …  
    3131    /// number of pair of values \a n
    3232    ///
    33     AveragerPair(const double x, const double xx, const double y,
    34                  const double yy, const double xy, const long);
     33    inline AveragerPair(const double x, const double xx, const double y,
     34                        const double yy, const double xy, const long n)
     35      : x_(Averager(x,xx,n)), y_(Averager(y,yy,n)), xy_(xy) {}
    3536
    3637    ///
    3738    /// Copy constructor
    3839    ///
    39     AveragerPair(const AveragerPair&);
     40    inline AveragerPair(const AveragerPair& a)
     41      : x_(a.x_averager()), y_(a.y_averager()), xy_(a.sum_xy()) {}
    4042   
    4143    ///

  • trunk/lib/statistics/AveragerWeighted.cc

    r397 r420  
    1111 
    1212
    13   AveragerWeighted::AveragerWeighted(void)
    14     : w_(Averager()), wx_(Averager()), wwx_(0), wxx_(0)
    15   {
    16   }
    17 
    18   AveragerWeighted::AveragerWeighted(const AveragerWeighted& a)
    19     : w_(Averager(a.sum_w(),a.sum_ww(),1)),
    20       wx_(Averager(a.sum_wx(),a.sum_wwxx(),1)),
    21       wwx_(a.sum_wwx()),
    22       wxx_(a.sum_wxx())
    23   {
    24   }
    25 
    2613  AveragerWeighted AveragerWeighted::operator+=(const AveragerWeighted& a)
    2714  {

  • trunk/lib/statistics/AveragerWeighted.h

    r414 r420  
    2424    ///
    2525    /// Default constructor
    26     ///
    27     AveragerWeighted(void);
     26    ///
     27    inline AveragerWeighted(void)
     28      : w_(Averager()), wx_(Averager()), wwx_(0), wxx_(0) {}
    2829
    2930    ///
    3031    /// Copy constructor
    3132    ///
    32     AveragerWeighted(const AveragerWeighted&);
     33    inline AveragerWeighted(const AveragerWeighted& a)
     34      : w_(Averager(a.sum_w(),a.sum_ww(),1)),
     35        wx_(Averager(a.sum_wx(),a.sum_wwxx(),1)), wwx_(a.sum_wwx()),
     36        wxx_(a.sum_wxx()) {}
    3337
    3438    ///

  • trunk/lib/statistics/Linear.h

    r389 r420  
    3232    /// Destructor
    3333    ///
    34     virtual ~Linear(void) {};
     34    inline virtual ~Linear(void) {};
    3535         
    3636    ///

  • trunk/lib/statistics/MultiDimensional.cc

    r386 r420  
    1010
    1111
    12   MultiDimensional::~MultiDimensional(void)
    13   {
    14     if (work_)
    15       gsl_multifit_linear_free(work_);
    16   }
    17 
    18 
    19 
    2012  void MultiDimensional::fit(const gslapi::matrix& x, const gslapi::vector& y)
    2113  {
    … …  
    2517      gsl_multifit_linear_free(work_);
    2618    work_=gsl_multifit_linear_alloc(x.rows(),fit_parameters_.size());
    27     gsl_multifit_linear(x.gsl_matrix_pointer(),y.gsl_vector_pointer(),
    28                         fit_parameters_.TEMP_gsl_vector_pointer(),
    29                         covariance_.gsl_matrix_pointer(),&chisquare_,work_);
     19    gsl_multifit_linear(x.gsl_matrix_p(),y.gsl_vector_p(),
     20                        fit_parameters_.gsl_vector_p(),
     21                        covariance_.gsl_matrix_p(),&chisquare_,work_);
    3022  }
    3123

  • trunk/lib/statistics/MultiDimensional.h

    r389 r420  
    2929    /// @brief Destructor
    3030    ///
    31     ~MultiDimensional(void);
     31    inline ~MultiDimensional(void) { if (work_) gsl_multifit_linear_free(work_);}
    3232
    3333    ///

  • trunk/lib/svm/InputRanker.cc

    r301 r420  
    1414namespace svm { 
    1515
    16   InputRanker::InputRanker()
    17     :train_set_(std::vector<size_t>()),
    18      id_(std::vector<size_t>()),
    19      rank_(std::vector<size_t>())
    20  
    21   {
    22   }
    23  
     16
    2417  InputRanker::InputRanker(const gslapi::matrix& data,
    2518                           const gslapi::vector& target,
    … …  
    4134    //scoring each input
    4235    std::vector<std::pair<size_t, double> > score;
    43     for (size_t i=0; i<nof_genes; i++){
    44       double area = score_object.score(target, data[i], train_set_); 
     36    for (size_t i=0; i<nof_genes; i++) {
     37      double area = score_object.score(target,gslapi::vector(data,i),train_set_);
    4538      std::pair<size_t, double> tmp(i,area);
    4639      score.push_back(tmp);
    … …  
    5649    }
    5750  }
     51
     52
    5853
    5954  InputRanker::InputRanker(const gslapi::matrix& data,
    … …  
    7671    //scoring each input
    7772    std::vector<std::pair<size_t, double> > score;
    78     for (size_t i=0; i<nof_genes; i++){
    79       double area = score_object.score(target, data[i],
    80                                        weight[i], train_set_); 
     73    for (size_t i=0; i<nof_genes; i++) {
     74      double area = score_object.score(target, gslapi::vector(data,i),
     75                                       gslapi::vector(weight,i), train_set_);
    8176      std::pair<size_t, double> tmp(i,area);
    8277      score.push_back(tmp);
    … …  
    9186    }
    9287  }
    93  
     88
     89
    9490}} // of namespace svm and namespace theplu

  • trunk/lib/svm/InputRanker.h

    r295 r420  
    2121  class InputRanker
    2222  {
    23  
     23
    2424  public:
    2525    ///
    26     /// Default constructor. Not implemented.
     26    /// Default constructor.
    2727    ///
    28     InputRanker();
     28    inline InputRanker(void) : train_set_(std::vector<size_t>()),
     29                               id_(std::vector<size_t>()),
     30                               rank_(std::vector<size_t>()) {}
    2931
    3032    ///
    … …  
    5456    ///
    5557    inline size_t id(const size_t i) const {return id_[i];}
    56    
     58
    5759    ///
    5860    /// highest ranked gene is ranked as number zero @return rank for
    … …  
    6668    std::vector<size_t> id_;
    6769    std::vector<size_t> rank_;
     70  };
    6871
    69          
    70   };
    7172
    7273}} // of namespace svm and namespace theplu
    7374
    7475#endif
    75 

  • trunk/lib/svm/Kernel_MEV.cc

    r336 r420  
    33#include <c++_tools/svm/Kernel_MEV.h>
    44
    5 #include <c++_tools/svm/KernelFunction.h>
    6 #include <c++_tools/gslapi/matrix.h>
    7 #include <c++_tools/gslapi/vector.h>
    8 
    95namespace theplu {
    106namespace svm { 
    117
    12   Kernel_MEV::Kernel_MEV(const gslapi::matrix& data, const KernelFunction& kf)
    13     : Kernel(data, kf), data_(data), kf_(&kf)
    14   {
    15   }
    16 
    17 
    18 
    19   Kernel_MEV::~Kernel_MEV(void)
    20   {
    21   }
    22 
    23 
    24 
    258
    269}} // of namespace svm and namespace theplu

  • trunk/lib/svm/Kernel_MEV.h

    r350 r420  
    3030    ///  Default constructor (not implemented)
    3131    ///
    32     Kernel_MEV();
     32    Kernel_MEV(void);
    3333
    3434    ///
    … …  
    3737    ///   sample. @note Can not handle NaNs.
    3838    ///
    39     Kernel_MEV(const gslapi::matrix&, const KernelFunction&);
    40    
     39    inline Kernel_MEV(const gslapi::matrix& data, const KernelFunction& kf)
     40      : Kernel(data,kf), data_(data), kf_(&kf) {}
     41
    4142    ///
    4243    ///   @todo Constructor taking the \a data matrix, the KernelFunction and a
    … …  
    5455    ///   Destructor
    5556    ///
    56     virtual ~Kernel_MEV(void);
     57    inline virtual ~Kernel_MEV(void) {};
    5758
    5859    ///
    59     /// @return element at position (\a row, \a column) of the Kernel
     60    /// @return Element at position (\a row, \a column) of the Kernel
    6061    /// matrix
    6162    ///
    62     double operator()(const size_t row,const size_t column) const
    63     { return (*kf_)(data_.TEMP_col_return(row),data_.TEMP_col_return(column)); }
     63    inline double operator()(const size_t row, const size_t column) const
     64      { return (*kf_)(gslapi::vector(data_,row,false),
     65                      gslapi::vector(data_,column,false)); }
    6466
    6567    ///

  • trunk/lib/svm/Kernel_SEV.cc

    r336 r420  
    1111namespace svm { 
    1212
     13
    1314Kernel_SEV::Kernel_SEV(const gslapi::matrix& data, const KernelFunction& kf)
    1415  : Kernel(data,kf), data_(data), kf_(&kf)
    … …  
    1718  for (size_t i=0; i<kernel_matrix_.rows(); i++)
    1819    for (size_t j=i; j<kernel_matrix_.columns(); j++)
    19       kernel_matrix_(i,j) = kernel_matrix_(j,i) =
    20         (*kf_)(data_.TEMP_col_return(i),data_.TEMP_col_return(j));
    21 
     20      kernel_matrix_(i,j) = kernel_matrix_(j,i) =
     21        (*kf_)(gslapi::vector(data_,i,false),gslapi::vector(data_,j,false));
    2222}
    2323
    24 Kernel_SEV::~Kernel_SEV(void)
    25 {
    26 } 
    2724
    2825}} // of namespace svm and namespace theplu

  • trunk/lib/svm/Kernel_SEV.h

    r336 r420  
    4848    ///   Destructor
    4949    ///
    50     virtual ~Kernel_SEV(void);
     50    inline virtual ~Kernel_SEV(void) {};
    5151
    5252    ///

  • trunk/lib/utility/Exception.h

    r411 r420  
    1717  public:
    1818    IO_error(void) throw() : std::runtime_error("IO_error:") {}
    19     IO_error(std::string& message) throw()
     19    IO_error(std::string message) throw()
    2020      : std::runtime_error("IO_error: " + message) {}
    2121  };

  • trunk/lib/utility/PCA.cc

    r301 r420  
    1212namespace utility { 
    1313
    14   //PCA::PCA() : process_( false ), explained_calc_( false ){}
    15 
    16 PCA::PCA( const gslapi::matrix& A ) : A_( A ),
    17               process_( false ),
    18               explained_calc_( false )
    19 {
    20 }
    2114
    2215void PCA::process()
    … …  
    3528
    3629  // Read the eigenvectors and eigenvalues
    37   eigenvectors_ = U.transpose();
    38   eigenvalues_ = pSVD->s();
     30  eigenvectors_ = U;
     31  eigenvectors_ .transpose();
     32  eigenvalues_ = pSVD->s();
    3933
    4034  // T
    … …  
    5044  // make sure that the i:th element is in its correct
    5145  // position (N element --> Ordo( N*N ))
    52   for( size_t i = 0; i < eigenvalues_.size(); ++i )
    53     {
    54      for( size_t j = i + 1; j < eigenvalues_.size(); ++j )
    55        {
    56   if( eigenvalues_[ j ] > eigenvalues_[ i ] )
    57     {
    58       std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
    59      eigenvectors_.row_swap(i,j);
    60     }
    61        }
    62     }
    63 }
     46  for ( size_t i = 0; i < eigenvalues_.size(); ++i )
     47    for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
     48      if ( eigenvalues_[ j ] > eigenvalues_[ i ] ) {
     49        std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
     50        eigenvectors_.swap_rows(i,j);
     51      }
     52}
     53
    6454
    6555
    … …  
    7363  // Transform into SVD friendly dimension
    7464  ////////////////////////////////////////
    75   A_ = A_.transpose();
    76   A_center = A_center.transpose();
     65  A_.transpose();
     66  A_center.transpose();
    7767
    7868  // Single value decompose the data matrix
    … …  
    8373
    8474  // Read the eigenvectors and eigenvalues
    85   eigenvectors_ = V.transpose();//U.transpose();
     75  eigenvectors_=V;
     76  eigenvectors_.transpose();
    8677  eigenvalues_ = pSVD->s();
    8778
    … …  
    8980  // (used V insted of U now for eigenvectors)
    9081  ////////////////////////////////////////////
    91   A_ = A_.transpose();
    92   A_center = A_center.transpose();
    93 
     82  A_.transpose();
     83  A_center.transpose();
    9484
    9585  // T
    … …  
    10595  // make sure that the i:th element is in its correct
    10696  // position (N element --> Ordo( N*N ))
    107   for( size_t i = 0; i < eigenvalues_.size(); ++i )
    108     {
    109      for( size_t j = i + 1; j < eigenvalues_.size(); ++j )
    110        {
    111   if( eigenvalues_[ j ] > eigenvalues_[ i ] )
    112     {
    113      std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
    114      eigenvectors_.row_swap(i,j);
    115     }
    116        }
    117     }
    118 }
     97  for ( size_t i = 0; i < eigenvalues_.size(); ++i )
     98    for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
     99      if ( eigenvalues_[ j ] > eigenvalues_[ i ] ) {
     100        std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
     101        eigenvectors_.swap_rows(i,j);
     102      }
     103}
     104
    119105
    120106
    … …  
    127113  gslapi::vector A_row_sum(A_.rows());
    128114  for (size_t i=0; i<A_row_sum.size(); ++i)
    129     A_row_sum(i)=A_[i].sum();
     115    A_row_sum(i)=gslapi::vector(A_,i).sum();
    130116  for( size_t i = 0; i < A_center.rows(); ++i )
    131117    {
    … …  
    135121    }
    136122}
     123
    137124
    138125
    … …  
    162149
    163150
     151
    164152gslapi::matrix PCA::projection_transposed( const gslapi::matrix&
    165153             samples ) const
    … …  
    205193}
    206194
     195
     196
    207197double PCA::get_explained_intensity( const size_t& k )
    208198{
    … …  
    215205}
    216206
     207
    217208}} // of namespace utility and namespace theplu

  • trunk/lib/utility/PCA.h

    r334 r420  
    66#include <c++_tools/gslapi/matrix.h>
    77#include <c++_tools/gslapi/vector.h>
    8 
    9 // Standard C++ includes
    10 ////////////////////////
    11 // #include <vector>
    12 // #include <iostream>
    13 // #include <memory>
    14 // #include <cstdlib>
    158
    169
    … …  
    2922     projection_transposed()...
    3023  */
    31  
    3224  class PCA
    3325  {
    … …  
    3628       Default constructor (not implemented)
    3729    */
    38     PCA();
     30    PCA(void);
    3931
    4032    /**
    … …  
    4234       should have been performed and no products.
    4335     */
    44     explicit PCA( const gslapi::matrix& );
     36    inline explicit PCA(const gslapi::matrix& A)
     37      : A_(A), process_(false), explained_calc_(false) {}
    4538
    46    
    4739    /**
    4840       Will perform PCA according to the following scheme: \n
    … …  
    5143       3: Calculate eigenvalues according to \n
    5244          \f$ \lambda_{ii} = s_{ii}/N_{rows} \f$ \n
    53        4: Sort eigenvectors (from matrix V) according to descending eigenvalues \n
     45       4: Sort eigenvectors (from matrix V) according to descending eigenvalues\n
    5446    */
    55     void process();
     47    void process(void);
    5648
    5749    /**
    … …  
    6254       projection_transposed() method.
    6355     */
    64     void process_transposed_problem();
    65    
     56    void process_transposed_problem(void);
     57
    6658    /**
    67        Returns eigenvector \a i
     59       @return Eigenvector \a i.
    6860    */
    69     // Jari, change this to
    70 //     const gslapi::vector& get_eigenvector( const size_t& i ) const
    71     const gslapi::vector get_eigenvector( const size_t& i ) const
    72     {
    73       return eigenvectors_[i];
    74     }
     61    inline gslapi::vector get_eigenvector(const size_t& i) const
     62      { return gslapi::vector(eigenvectors_,i); }
    7563
    7664    /**
    … …  
    7866       \f$ C = \frac{1}{N^2}A^TA \f$
    7967    */
    80     double get_eigenvalue( const size_t& i ) const
    81     {
    82       return eigenvalues_[ i ];
    83     }
     68    inline double
     69    get_eigenvalue(const size_t& i) const { return eigenvalues_[i]; }
    8470
    8571    /**
    … …  
    9076    double get_explained_intensity( const size_t& k );
    9177
    92 
    93 
    94     /**
    95         This function will project data onto the new coordinate-system
    96   where the axes are the calculated eigenvectors. This means that
    97   PCA must have been run before this function can be used!
    98   Output is presented as coordinates in the N-dimensional room
    99   spanned by the eigenvectors.
     78    /**
     79       This function will project data onto the new coordinate-system
     80       where the axes are the calculated eigenvectors. This means that
     81       PCA must have been run before this function can be used!
     82       Output is presented as coordinates in the N-dimensional room
     83       spanned by the eigenvectors.
    10084    */
    10185    gslapi::matrix projection( const gslapi::matrix& ) const;
    102    
    103     /**
    104   Same as projection() but works when used process_transposed_problem()
     86
     87    /**
     88       Same as projection() but works when used
     89       process_transposed_problem().
    10590    */
    10691    gslapi::matrix projection_transposed( const gslapi::matrix& ) const;
    10792
    108    
    109    
     93
    11094  private:
    111     gslapi::matrix A_;
     95    gslapi::matrix  A_;
    11296    gslapi::matrix  eigenvectors_;
    11397    gslapi::vector  eigenvalues_;
    … …  
    128112    */
    129113    void calculate_explained_intensity();
     114  }; // class PCA 
    130115
    131    
    132   }; // class PCA 
    133  
     116
    134117}} // of namespace utility and namespace theplu
    135118
    136119#endif
    137 

  • trunk/lib/utility/SVD.cc

    r301 r420  
    66namespace theplu {
    77namespace utility {
    8 
    9 
    10   SVD::SVD(const gslapi::matrix& Ain)
    11     : U_(Ain), V_(Ain.columns(),Ain.columns()), s_(Ain.columns())
    12   {
    13   }
    14 
    15 
    16 
    17   SVD::~SVD(void)
    18   {
    19   }
    20 
    218
    229
    … …  
    4229  {
    4330    gslapi::vector w(U_.columns());
    44     return gsl_linalg_SV_decomp(U_.gsl_matrix_pointer(),
    45                                 V_.gsl_matrix_pointer(),
    46                                 s_.TEMP_gsl_vector_pointer(),
    47                                 w.TEMP_gsl_vector_pointer());
     31    return gsl_linalg_SV_decomp(U_.gsl_matrix_p(), V_.gsl_matrix_p(),
     32                                s_.gsl_vector_p(), w.gsl_vector_p());
    4833  }
    4934
    … …  
    5439    gslapi::vector w(U_.columns());
    5540    gslapi::matrix X(U_.columns(),U_.columns());
    56     return gsl_linalg_SV_decomp_mod(U_.gsl_matrix_pointer(),
    57                                     X.gsl_matrix_pointer(),
    58                                     V_.gsl_matrix_pointer(),
    59                                     s_.TEMP_gsl_vector_pointer(),
    60                                     w.TEMP_gsl_vector_pointer());
     41    return gsl_linalg_SV_decomp_mod(U_.gsl_matrix_p(), X.gsl_matrix_p(),
     42                                    V_.gsl_matrix_p(), s_.gsl_vector_p(),
     43                                    w.gsl_vector_p());
    6144  }
    6245

  • trunk/lib/utility/SVD.h

    r303 r420  
    5252    /// input matrix is copied for further use in the object.
    5353    ///
    54     SVD(const gslapi::matrix& );
     54    inline SVD(const gslapi::matrix& Ain)
     55      : U_(Ain), V_(Ain.columns(),Ain.columns()), s_(Ain.columns()) {}
    5556
    56     ~SVD(void);
     57    inline ~SVD(void) {}
    5758
    5859    ///
    … …  
    7273    /// is undefined.
    7374    ///
    74     gslapi::vector s(void) const { return s_; }
     75    inline const gslapi::vector& s(void) const { return s_; }
    7576
    7677    ///
    … …  
    8384    ///
    8485    inline int solve(gslapi::vector b, gslapi::vector x)
    85       { return gsl_linalg_SV_solve(U_.gsl_matrix_pointer(),
    86                                    V_.gsl_matrix_pointer(),
    87                                    s_.gsl_vector_pointer(),
    88                                    b.gsl_vector_pointer(),
    89                                    x.TEMP_gsl_vector_pointer());
    90       }
     86      { return gsl_linalg_SV_solve(U_.gsl_matrix_p(), V_.gsl_matrix_p(),
     87                                   s_.gsl_vector_p(), b.gsl_vector_p(),
     88                                   x.gsl_vector_p()); }
    9189
    9290    ///
    … …  
    9896    /// is undefined.
    9997    ///
    100     gslapi::matrix U(void) const { return U_; }
     98    inline const gslapi::matrix& U(void) const { return U_; }
    10199
    102100    ///
    … …  
    108106    /// is undefined.
    109107    ///
    110     gslapi::matrix V(void) const { return V_; }
     108    inline const gslapi::matrix& V(void) const { return V_; }
    111109
    112110  private:
    113111    inline int jacobi(void)
    114       { return gsl_linalg_SV_decomp_jacobi(U_.gsl_matrix_pointer(),
    115                                            V_.gsl_matrix_pointer(),
    116                                            s_.TEMP_gsl_vector_pointer());
    117       }
     112      { return gsl_linalg_SV_decomp_jacobi(U_.gsl_matrix_p(), V_.gsl_matrix_p(),
     113                                           s_.gsl_vector_p()); }
    118114    int golub_reinsch(void);
    119115    int modified_golub_reinsch(void);

  • trunk/lib/utility/stl_utility.cc

    r341 r420  
    2020      for (size_t i=0; i<vec_str.size(); i++) {
    2121        if (!is_float(vec_str[i])){
    22           std::cerr << "Warning: '" << vec_str[i]
    23                     << "' is not a number." << std::endl;
     22          // Jari, this should be communicated with as an exception.
     23//          std::cerr << "Warning: '" << vec_str[i]
     24//                    << "' is not a number." << std::endl;
    2425        }
    2526        else
    … …  
    4142      for (size_t i=0; i<vec_str.size(); i++) {
    4243        if (!is_int(vec_str[i])){
    43           std::cerr << "Warning: '" << vec_str[i]
    44                     << "' is not an integer." << std::endl;
     44          // Jari, this should be communicated with as an exception.
     45//          std::cerr << "Warning: '" << vec_str[i]
     46//                    << "' is not an integer." << std::endl;
    4547        }
    4648        else

  • trunk/test/data/knni_matrix.data

    r175 r420  
     1
    12.2348 .234773 .098934 .21627  .32782  .2328 .347454 .25784  .345376 .56582   .1234
    23.2346 .786345 .234734 .34935  .23733  .2347 .21347  .286758 .456454 .97732   .1234
    34.8964 .471294 .947234 .38272  .48328  .2874 .489892 .23478  .45645  .23899   .1234
    4 .4372 .382144 .763342 .32924  .27823  .3827 .234764 .384734 .398795 .23498   .1234
     5  .4372 .382144 .763342 .32924  .27823  .3827 .234764 .384734 .398795 .23498   .1234
     6
    57.8624 .286482 .83248  .72846  .28468  .2347 .23234  .98973  .29874  .32894   .1234
    68.8624 .286482 .83248  .72846  .28468  .2347 .23234  .98973  .29874  .32894   .1234
    79.2342 .23878  .97437  .76236  .32764  .3474 .43276  .34555  .29948  .34598   .1234
    810.2347 .56783  .78435  .35763  .23468  .7384 .23578  .345533 .235948 .32458   .1234
    9 .2234 .77639  .38783  .28953  .62376  .5445 .32489  .22388  .238746 .23468   .1234
     11.2234 .77639  .38783  .28953  .62376  .5445 .32489  .22388  .238746 .23468   .1234 
    1012.2388 .234688 .43687  .86456  .65654  .6864 .879645 .98454  .98731  .164897  .1234
    1113.1968 .153987 .13684  .14685  .316135 .5286 .69435  .13634  .18351  .1984531 .1234
    1214.6865 .213687 .12368  .18564  .897654 .4945 .789654 .45547  .98466  .66136   .1234
    13 .9835 .963154 .964664 .132369 .983316 .3276 .376513 .34941  .387433 .31685   .1234
     15.9835 .963154 .964664 .132369     .983316 .3276 .376513 .34941  .387433 .31685   .1234
    1416.9783 .973154 .897432 .643344 .64343  .6731 .69436  .765313 .547865 .89713   .1234
    1517.9783 .973154 .897432 .643344 .64343  .6731 .69436  .765313 .547865 .89713   .1234

  • trunk/test/matrix_test.cc

    r399 r420  
    33#include <c++_tools/gslapi/matrix.h>
    44
     5#include <unistd.h>
    56#include <fstream>
    67#include <iostream>
    78
     9class matrixwrapper
     10{
     11public:
     12  matrixwrapper(size_t i,size_t j,double value=7)
     13    : m_(i,j,value) {}
     14
     15  inline theplu::gslapi::vector
     16  row(const size_t& i) { return theplu::gslapi::vector(m_,i); }
     17
     18  inline const theplu::gslapi::matrix& matrix(void) const { return m_; }
     19
     20private:
     21  theplu::gslapi::matrix m_;
     22};
     23
     24
    825int main(const int argc,const char* argv[])
    9 { 
     26{
    1027  using namespace theplu;
    1128  std::ostream* error;
    … …  
    1734      std::cout << "matrix_test -v : for printing extra information\n";
    1835  }
    19   *error << "testing matrix" << std::endl;
     36
     37  *error << "Testing matrix class" << std::endl;
    2038  bool ok = true;
    2139
     40  *error << "\tcopy constructor and operator!=" << std::endl;
    2241  gslapi::matrix m(3,3,9);
    2342  gslapi::matrix m2(m);
    2443  if (m2!=m)
    2544    ok=false;
     45
     46  *error << "\toutput operator and istream constructor" << std::endl;
     47  // Checking that the matrix output operator writes a file that the
     48  // input operator can read.
    2649  std::ofstream my_out("data/tmp_test_matrix.txt");
    2750  my_out << m2;
    2851  my_out.close();
    29 
    3052  std::ifstream is("data/tmp_test_matrix.txt");
    3153  gslapi::matrix m3(is);
    … …  
    3355  if (m3!=m2)
    3456    ok=false;
    35 
     57  unlink("data/tmp_test_matrix.txt");
     58
     59  *error << "\toperator*(double)" << std::endl;
    3660  gslapi::matrix m4(3,3,1);
    3761  m4 *= 9;
    38  
    39   if (m4!=m){
     62  if (m4!=m) {
    4063    ok=false;
    4164    *error << "error operator*=(double)" << std::endl;
    4265  }
    4366
     67  *error << "\tsub-matrix" << std::endl;
     68  // Checking that sub-matrices work, i.e. mutation to the view are
     69  // reflected in the viewed object.
     70  is.open("data/knni_matrix.data");
     71  // The stream input is a proper matrix file, with some stray empty
     72  // lines and other whitespaces. The file is not expected to break
     73  // things.
     74  gslapi::matrix m5(is);
     75  is.close();
     76  double m5_sum=0;
     77  for (size_t i=0; i<m5.rows(); ++i)
     78    for (size_t j=0; j<m5.columns(); ++j)
     79      m5_sum+=m5(i,j);
     80  gslapi::matrix* m5sub=new gslapi::matrix(m5,3,3,3,3);
     81  double m5sub_sum=0;
     82  for (size_t i=0; i<m5sub->rows(); ++i)
     83    for (size_t j=0; j<m5sub->columns(); ++j) {
     84      m5sub_sum+=(*m5sub)(i,j);
     85      (*m5sub)(i,j)=1;
     86    }
     87  delete m5sub;
     88  double m5_sum2=0;
     89  for (size_t i=0; i<m5.rows(); ++i)
     90    for (size_t j=0; j<m5.columns(); ++j)
     91      m5_sum2+=m5(i,j);
     92  if (m5_sum2-m5_sum-9+m5sub_sum>1e-13) {
     93    ok=false;
     94    *error << "error sub-matrix test" << std::endl;
     95  }
     96
     97  *error << "\tsub-(row)vector" << std::endl;
     98  // Checking that the row view works, i.e. mutation to the view are
     99  // reflected in the viewed object.
     100  gslapi::vector v5subrow(m5,3);
     101  double v5subrow_sum=0;
     102  for (size_t i=0; i<v5subrow.size(); ++i) {
     103    v5subrow_sum+=v5subrow(i);
     104    v5subrow[i]=0;
     105  }
     106  double m5_sum3=0;
     107  for (size_t i=0; i<m5.rows(); ++i)
     108    for (size_t j=0; j<m5.columns(); ++j)
     109      m5_sum3+=m5(i,j);
     110  if (m5_sum3-m5_sum2+v5subrow_sum>1e-13) {
     111    ok=false;
     112    *error << "error sub-vector test 1" << std::endl;
     113  }
     114
     115  *error << "\tsub-(column)vector" << std::endl;
     116  // Checking that the column view works, i.e. mutation to the view
     117  // are reflected in the viewed object.
     118  gslapi::vector v5subcolumn(m5,0,false);
     119  double v5subcolumn_sum=0;
     120  for (size_t i=0; i<v5subcolumn.size(); ++i) {
     121    v5subcolumn_sum+=v5subcolumn(i);
     122    v5subcolumn(i)=1;
     123  }
     124  double m5_sum4=0;
     125  for (size_t i=0; i<m5.rows(); ++i)
     126    for (size_t j=0; j<m5.columns(); ++j)
     127      m5_sum4+=m5(i,j);
     128  if (m5_sum4-m5_sum3-v5subcolumn.size()+v5subcolumn_sum>1e-13) {
     129    ok=false;
     130    *error << "error sub-vector test 2" << std::endl;
     131  }
     132  // Checking that the column view above mutates the values in the row
     133  // view.
     134  double v5subrow_sum2=0;
     135  for (size_t i=0; i<v5subrow.size(); ++i)
     136    v5subrow_sum2+=v5subrow(i);
     137  if (v5subrow_sum2-v5subcolumn(3)>1e-13) {
     138    ok=false;
     139    *error << "error sub-vector test 3" << std::endl;
     140  }
     141
     142  *error << "\tsub-vector and vector copying" << std::endl;
     143  // Checking that a view is not inherited through the copy
     144  // contructor.
     145  gslapi::vector v6(v5subrow);
     146  v6.set_all(2);
     147  double v5subrow_sum3=0;
     148  for (size_t i=0; i<v5subrow.size(); ++i)
     149    v5subrow_sum3+=v5subrow(i);
     150  if (v5subrow_sum3-v5subrow_sum2>1e-13) {
     151    ok=false;
     152    *error << "error sub-vector test 4" << std::endl;
     153  }
     154  // Checking that values in a vector is copied into a viewed matrix.
     155  v5subrow.set(v6);
     156  double m5_sum5=0;
     157  for (size_t i=0; i<m5.rows(); ++i)
     158    for (size_t j=0; j<m5.columns(); ++j)
     159      m5_sum5+=m5(i,j);
     160  if (m5_sum5-m5_sum4-v5subrow.size()*2+v5subrow_sum3>1e-13) {
     161    ok=false;
     162    *error << "error sub-vector test 5" << std::endl;
     163  }
     164  // Checking that vector::operator= indeed makes a view to become a
     165  // "normal" vector.
     166  v5subrow=gslapi::vector(23,22);
     167  double m5_sum6=0;
     168  for (size_t i=0; i<m5.rows(); ++i)
     169    for (size_t j=0; j<m5.columns(); ++j)
     170      m5_sum6+=m5(i,j);
     171  if (m5_sum6-m5_sum5>1e-13) {
     172    ok=false;
     173    *error << "error sub-vector test 6" << std::endl;
     174  }
     175
     176  // Checking that the memberfunction matrixwrapper::row() returns a
     177  // view
     178  *error << "\tthat class member returns a view" << std::endl;
     179  matrixwrapper mw(5,2);
     180  gslapi::vector mwrow=mw.row(2);
     181  if (mwrow.gsl_vector_p()->data != &(mw.matrix()(2,0))) {
     182    ok=false;
     183    *error << "error sub-vector test 7" << std::endl;
     184  }
     185
     186  *error << "\tsub-matrix of a sub-matrix" << std::endl;
     187  // Checking that a sub-matrix of a sub-matrix can be created.
     188  gslapi::matrix sub(m5,3,3,3,3);
     189  gslapi::matrix subsub(sub,2,1,1,2);
     190  subsub(0,0)=23221121;
     191  if (&sub(2,1)!=&subsub(0,0) || subsub(0,0)!=m5(5,4) || &subsub(0,0)!=&m5(5,4)){
     192    ok=false;
     193    *error << "error sub-matrix test 1" << std::endl;
     194  }
    44195
    45196  if (error!=&std::cerr)
    … …  
    47198
    48199  return (ok ? 0 : -1);
    49 
    50200}

  • trunk/test/score_test.cc

    r414 r420  
    7979  const double tol = 0.001;
    8080  for (size_t i=0; i<data.rows(); i++){
    81     gslapi::vector vec = data[i];
     81    gslapi::vector vec(data,i);
    8282    area = roc.score(target2,vec);
    8383    if (area<correct_area(i)-tol || area>correct_area(i)+tol){
    … …  
    9090  gslapi::vector weight(target2.size(),1);
    9191  for (size_t i=0; i<data.rows(); i++){
    92     gslapi::vector vec = data[i];
     92    gslapi::vector vec(data,i);
    9393    area = roc.score(target2, vec, weight);
    9494    if (area<correct_area(i)-tol || area>correct_area(i)+tol){

  • trunk/test/svd_test.cc

    r377 r420  
    3737  for (size_t i=0; i<s.size(); ++i)
    3838    S(i,i)=s[i];
    39   gslapi::matrix V=svd.V();
    40   gslapi::matrix U=svd.U();
    41   double error = this_norm(A-U*S*V.transpose());
     39  gslapi::matrix Vtranspose=svd.V();
     40  Vtranspose.transpose();
     41  // Reconstructing A = U*S*Vtranspose
     42  gslapi::matrix Areconstruct=svd.U();
     43  Areconstruct*=S;
     44  Areconstruct*=Vtranspose;
     45  Areconstruct-=A;  // Expect null matrix
     46  double error = this_norm(Areconstruct);
    4247  bool testerror=false;
    4348  if (error>MAXTOL) {
    … …  
    4954  }
    5055
    51   error=this_norm(V.transpose()*V)-n;
     56  Vtranspose*=svd.V();  // Expect unity matrix
     57  error=this_norm(Vtranspose)-n;
    5258  if (error>MAXTOL) {
    5359    std::cerr << "test_svd: FAILED, algorithm " << algo
    … …  
    5763  }
    5864
    59   error=this_norm(U.transpose()*U)-n;
     65  gslapi::matrix Utranspose=svd.U();
     66  Utranspose.transpose();
     67  Utranspose*=svd.U();  // Expect unity matrix
     68  error=this_norm(Utranspose)-n;
    6069  if (error>MAXTOL) {
    6170    std::cerr << "test_svd: FAILED, algorithm " << algo

  • trunk/test/svm_test.cc

    r337 r420  
    8282     
    8383  // Comparing alpha to alpha_matlab
    84   theplu::gslapi::vector diff_alpha = alpha - alpha_matlab;
     84  theplu::gslapi::vector diff_alpha(alpha);
     85  diff_alpha-=alpha_matlab;
    8586  if (diff_alpha*diff_alpha> 1e-10 ){
    8687    if (print)
    … …  
    9091
    9192  // Comparing output to target
    92   theplu::gslapi::vector output = svm.output();
     93  theplu::gslapi::vector output(svm.output());
    9394  double slack = 0;
    9495  for (unsigned int i=0; i<target.size(); i++){
Note: See TracChangeset for help on using the changeset viewer.