Changeset 792 for trunk/yat


Ignore:
Timestamp:
Mar 11, 2007, 1:14:24 AM (15 years ago)
Author:
Jari Häkkinen
Message:

Addresses #193. matrix now works as outlined in the ticket
discussion. Added support for const views. Added a clone function that
facilitates resizing of matrices. clone is needed since assignement
operator functionality is changed.

Location:
trunk/yat
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/yat/classifier/Kernel_SEV.cc

    r746 r792  
    4747    : Kernel(data,kf, own)
    4848  {
    49     kernel_matrix_ = utility::matrix(data_->columns(),data_->columns());
     49    kernel_matrix_.clone(utility::matrix(data_->columns(),data_->columns()));
    5050    for (size_t i=0; i<kernel_matrix_.rows(); i++)
    5151      for (size_t j=i; j<kernel_matrix_.columns(); j++)
    … …  
    6666  void Kernel_SEV::build_kernel(void)
    6767  {
    68     kernel_matrix_ = utility::matrix(data_->columns(),data_->columns());
     68    kernel_matrix_.clone(utility::matrix(data_->columns(),data_->columns()));
    6969    for (size_t i=0; i<kernel_matrix_.rows(); i++)
    7070      for (size_t j=i; j<kernel_matrix_.columns(); j++)

  • trunk/yat/classifier/NCC.cc

    r789 r792  
    8787  bool NCC::train()
    8888  {   
    89     centroids_=utility::matrix(data_.rows(), target_.nof_classes());
     89    centroids_.clone(utility::matrix(data_.rows(), target_.nof_classes()));
    9090    utility::matrix nof_in_class(data_.rows(), target_.nof_classes());
    9191    for(size_t i=0; i<data_.rows(); i++) {
    … …  
    101101      }
    102102    }   
    103     centroids_.div_elements(nof_in_class);
     103    centroids_.div(nof_in_class);
    104104    trained_=true;
    105105    return trained_;
    … …  
    118118    // take care of nan's in centroids
    119119    for(size_t j=0; j<centroids_.columns(); j++) {
    120       const utility::vector centroid=utility::vector(centroids_,j,false);
     120      const utility::vector centroid(utility::vector(centroids_,j,false));
    121121      utility::vector wc(centroid.size(),0);
    122122      for(size_t i=0; i<centroid.size(); i++)  {
    … …  
    132132                    utility::matrix& prediction) const
    133133  {   
    134     prediction=utility::matrix(centroids_.columns(), input.columns());   
     134    prediction.clone(utility::matrix(centroids_.columns(), input.columns()));
    135135    try {   
    136136      const MatrixLookupWeighted& data=

  • trunk/yat/classifier/SVM.cc

    r747 r792  
    146146
    147147    assert(input.rows()==alpha_.size());
    148     prediction = utility::matrix(2,input.columns(),0);
     148    prediction.clone(utility::matrix(2,input.columns(),0));
    149149    for (size_t i = 0; i<input.columns(); i++){
    150150      for (size_t j = 0; j<input.rows(); j++){
    … …  
    182182  {
    183183    trained_=false;
    184     alpha_=utility::vector(target_.size(), 0);
     184    alpha_.clone(utility::vector(target_.size(), 0));
    185185  }
    186186
    … …  
    413413  }
    414414
    415 
    416415}}} // of namespace classifier, yat, and theplu

  • trunk/yat/regression/MultiDimensional.cc

    r789 r792  
    5656  {
    5757    assert(x.rows()==y.size());
    58     covariance_=utility::matrix(x.columns(),x.columns());
     58    covariance_.clone(utility::matrix(x.columns(),x.columns()));
    5959    fit_parameters_.clone(utility::vector(x.columns()));
    6060    if (work_)

  • trunk/yat/regression/MultiDimensionalWeighted.cc

    r789 r792  
    5858    assert(x.rows()==y.size());
    5959
    60     covariance_=utility::matrix(x.columns(),x.columns());
     60    covariance_.clone(utility::matrix(x.columns(),x.columns()));
    6161    fit_parameters_.clone(utility::vector(x.columns()));
    6262    if (work_)

  • trunk/yat/utility/PCA.cc

    r789 r792  
    6363    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
    6464    pSVD->decompose();
    65     utility::matrix U = pSVD->U();
    66     utility::matrix V = pSVD->V();
     65    utility::matrix U(pSVD->U());
     66    utility::matrix V(pSVD->V());
    6767
    6868    // Read the eigenvectors and eigenvalues
    69     eigenvectors_ = U;
     69    eigenvectors_.clone(U);
    7070    eigenvectors_ .transpose();
    7171    eigenvalues_.clone(pSVD->s());
    … …  
    103103    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
    104104    pSVD->decompose();
    105     utility::matrix U = pSVD->U();
    106     utility::matrix V = pSVD->V();
     105    utility::matrix U(pSVD->U());
     106    utility::matrix V(pSVD->V());
    107107
    108108    // Read the eigenvectors and eigenvalues
    109     eigenvectors_=V;
     109    eigenvectors_.clone(V);
    110110    eigenvectors_.transpose();
    111111    eigenvalues_.clone(pSVD->s());
    … …  
    197197  {
    198198    size_t DIM = eigenvalues_.size();
    199     explained_intensity_ = utility::vector( DIM );
     199    explained_intensity_.clone(utility::vector( DIM ));
    200200    double sum = 0;
    201201    for( size_t i = 0; i < DIM; ++i )

  • trunk/yat/utility/matrix.cc

    r788 r792  
    2929#include "utility.h"
    3030
     31#include <cassert>
    3132#include <cmath>
    3233#include <sstream>
    … …  
    4142
    4243  matrix::matrix(void)
    43     : blas_result_(NULL), m_(NULL), view_(NULL)
     44    : blas_result_(NULL), m_(NULL), view_(NULL), view_const_(NULL),
     45      proxy_m_(NULL)
    4446  {
    4547  }
    … …  
    4749
    4850  matrix::matrix(const size_t& r, const size_t& c, double init_value)
    49     : blas_result_(NULL), m_(gsl_matrix_alloc(r,c)), view_(NULL)
     51    : blas_result_(NULL), m_(gsl_matrix_alloc(r,c)), view_(NULL),
     52      view_const_(NULL), proxy_m_(m_)
    5053  {
    5154    if (!m_)
    … …  
    5659
    5760  matrix::matrix(const matrix& o)
    58     : blas_result_(NULL), m_(o.create_gsl_matrix_copy()), view_(NULL)
     61    : blas_result_(NULL), m_(o.create_gsl_matrix_copy()), view_(NULL),
     62      view_const_(NULL), proxy_m_(m_)
    5963  {
    6064  }
    … …  
    6367  matrix::matrix(matrix& m, size_t offset_row, size_t offset_column,
    6468                 size_t n_row, size_t n_column)
    65     : blas_result_(NULL)
     69    : blas_result_(NULL), view_const_(NULL)
    6670  {
    6771    view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_,
    … …  
    7074    if (!view_)
    7175      throw utility::GSL_error("matrix::matrix failed to setup view");
    72     m_ = &(view_->matrix);
     76    proxy_m_ = m_ = &(view_->matrix);
    7377  }
    7478
    … …  
    7781  matrix::matrix(std::istream& is, char sep)
    7882    throw (utility::IO_error,std::exception)
    79     : blas_result_(NULL), view_(NULL)
     83    : blas_result_(NULL), view_(NULL), view_const_(NULL)
    8084  {
    8185    // read the data file and store in stl vectors (dynamically
    … …  
    133137    is.clear(std::ios::goodbit);
    134138    // convert the data to a gsl matrix
    135     m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
     139    proxy_m_ = m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
    136140    if (!m_)
    137141      throw utility::GSL_error("matrix::matrix failed to allocate memory");
    … …  
    149153    if (view_)
    150154      delete view_;
    151     else {
    152       if (blas_result_)
    153         gsl_matrix_free(blas_result_);
    154       if (m_)
    155         gsl_matrix_free(m_);
    156     }
     155    else if (view_const_)
     156      delete view_const_;
     157    else if (m_)
     158      gsl_matrix_free(m_);
     159    if (blas_result_)
     160      gsl_matrix_free(blas_result_);
    157161    blas_result_=NULL;
    158162    m_=NULL;
    … …  
    160164
    161165
     166  const matrix& matrix::clone(const matrix& other)
     167  {
     168    if (this!=&other) {
     169
     170      if (view_)
     171        delete view_;
     172      else if (view_const_)
     173        delete view_const_;
     174      else if (m_)
     175        gsl_matrix_free(m_);
     176      if (blas_result_) {
     177        gsl_matrix_free(blas_result_);
     178        blas_result_=NULL;
     179      }
     180      proxy_m_=m_=NULL;
     181
     182      if (other.view_) {
     183        view_ = new gsl_matrix_view(*other.view_);
     184        proxy_m_ = m_ = &(view_->matrix);
     185      }
     186      else if (other.view_const_) {
     187        view_const_ = new gsl_matrix_const_view(*other.view_const_);
     188        proxy_m_ = &(view_const_->matrix);
     189      }
     190      else if (other.m_)
     191        proxy_m_ = m_ = other.create_gsl_matrix_copy();
     192
     193    }
     194    return *this;
     195  }
     196
     197
    162198  size_t matrix::columns(void) const
    163199  {
    164     if (!m_)
     200    if (!proxy_m_)
    165201      return 0;
    166     return m_->size2;
     202    return proxy_m_->size2;
    167203  }
    168204
    … …  
    173209    if (!m)
    174210      throw utility::GSL_error("matrix::create_gsl_matrix_copy failed to allocate memory");
    175     if (gsl_matrix_memcpy(m,m_))
     211    if (gsl_matrix_memcpy(m,proxy_m_))
    176212      throw utility::GSL_error("matrix::create_gsl_matrix_copy dimension mis-match");
    177213    return m;
    … …  
    179215
    180216
    181   void matrix::div_elements(const matrix& b)
    182   {
    183     int status=gsl_matrix_div_elements(m_, b.m_);
     217  void matrix::div(const matrix& other)
     218  {
     219    assert(m_);
     220    int status=gsl_matrix_div_elements(m_, other.gsl_matrix_p());
    184221    if (status)
    185222      throw utility::GSL_error(std::string("matrix::div_elements",status));
    … …  
    205242  const gsl_matrix* matrix::gsl_matrix_p(void) const
    206243  {
     244    return proxy_m_;
     245  }
     246
     247
     248  gsl_matrix* matrix::gsl_matrix_p(void)
     249  {
    207250    return m_;
    208251  }
    209252
    210253
    211   gsl_matrix* matrix::gsl_matrix_p(void)
    212   {
    213     return m_;
    214   }
    215 
    216 
    217254  bool matrix::isview(void) const
    218255  {
    219     return view_;
    220   }
    221 
    222 
    223   void matrix::mul_elements(const matrix& b)
    224   {
    225     int status=gsl_matrix_mul_elements(m_, b.m_);
     256    return view_ || view_const_;
     257  }
     258
     259
     260  void matrix::mul(const matrix& other)
     261  {
     262    assert(m_);
     263    int status=gsl_matrix_mul_elements(m_, other.gsl_matrix_p());
    226264    if (status)
    227265      throw utility::GSL_error(std::string("matrix::mul_elements",status));
    228    
    229266  }
    230267
    … …  
    232269  size_t matrix::rows(void) const
    233270  {
    234     if (!m_)
     271    if (!proxy_m_)
    235272      return 0;
    236     return m_->size1;
    237   }
    238 
    239 
    240   void matrix::set(const matrix& mat)
    241   {
    242     if (gsl_matrix_memcpy(m_, mat.m_))
     273    return proxy_m_->size1;
     274  }
     275
     276
     277  void matrix::set(const matrix& other)
     278  {
     279    assert(m_);
     280    if (gsl_matrix_memcpy(m_, other.gsl_matrix_p()))
    243281      throw utility::GSL_error("matrix::create_gsl_matrix_copy dimension mis-match");
    244282  }
    … …  
    247285  void matrix::set_all(const double value)
    248286  {
     287    assert(m_);
    249288    gsl_matrix_set_all(m_, value);
    250289  }
    … …  
    253292  void matrix::set_column(const size_t column, const vector& vec)
    254293  {
     294    assert(m_);
    255295    int status=gsl_matrix_set_col(m_, column, vec.gsl_vector_p());
    256296    if (status)
    … …  
    261301  void matrix::set_row(const size_t row, const vector& vec)
    262302  {
     303    assert(m_);
    263304    int status=gsl_matrix_set_row(m_, row, vec.gsl_vector_p());
    264305    if (status)
    … …  
    269310  void matrix::swap_columns(const size_t i, const size_t j)
    270311  {
     312    assert(m_);
    271313    int status=gsl_matrix_swap_columns(m_, i, j);
    272314    if (status)
    … …  
    277319  void matrix::swap_rowcol(const size_t i, const size_t j)
    278320  {
     321    assert(m_);
    279322    int status=gsl_matrix_swap_rowcol(m_, i, j);
    280323    if (status)
    … …  
    285328  void matrix::swap_rows(const size_t i, const size_t j)
    286329  {
     330    assert(m_);
    287331    int status=gsl_matrix_swap_rows(m_, i, j);
    288332    if (status)
    … …  
    293337  void matrix::transpose(void)
    294338  {
     339    assert(m_);
    295340    if (columns()==rows())
    296341      gsl_matrix_transpose(m_); // this never fails
    … …  
    302347      gsl_matrix_transpose_memcpy(transposed,m_);
    303348      gsl_matrix_free(m_);
    304       m_=transposed;
     349      proxy_m_ = m_ = transposed;
    305350      if (blas_result_) {
    306351        gsl_matrix_free(blas_result_);
    … …  
    313358  double& matrix::operator()(size_t row, size_t column)
    314359  {
     360    assert(m_);
    315361    double* d=gsl_matrix_ptr(m_, row, column);
    316362    if (!d)
    … …  
    322368  const double& matrix::operator()(size_t row, size_t column) const
    323369  {
    324     const double* d=gsl_matrix_const_ptr(m_, row, column);
     370    const double* d=gsl_matrix_const_ptr(proxy_m_, row, column);
    325371    if (!d)
    326372      throw utility::GSL_error("matrix::operator()",GSL_EINVAL);
    … …  
    343389  const matrix& matrix::operator=( const matrix& other )
    344390  {
    345     if ( this != &other ) {
    346       if (view_)
    347         delete view_;
    348       else if (m_)
    349         gsl_matrix_free(m_);
    350       m_ = other.create_gsl_matrix_copy();
    351     }
    352     return *this;
    353   }
    354 
    355 
    356   const matrix& matrix::operator+=(const matrix& m)
    357   {
    358     int status=gsl_matrix_add(m_, m.m_);
     391    set(other);
     392    return *this;
     393  }
     394
     395
     396  const matrix& matrix::operator+=(const matrix& other)
     397  {
     398    assert(m_);
     399    int status=gsl_matrix_add(m_, other.proxy_m_);
    359400    if (status)
    360401      throw utility::GSL_error(std::string("matrix::operator+=", status));
    … …  
    365406  const matrix& matrix::operator+=(const double d)
    366407  {
     408    assert(m_);
    367409    gsl_matrix_add_constant(m_, d);
    368410    return *this;
    … …  
    370412
    371413
    372   const matrix& matrix::operator-=(const matrix& m)
    373   {
    374     int status=gsl_matrix_sub(m_, m.m_);
     414  const matrix& matrix::operator-=(const matrix& other)
     415  {
     416    assert(m_);
     417    int status=gsl_matrix_sub(m_, other.proxy_m_);
    375418    if (status)
    376419      throw utility::GSL_error(std::string("matrix::operator-=", status));
    … …  
    381424  const matrix& matrix::operator-=(const double d)
    382425  {
     426    assert(m_);
    383427    gsl_matrix_add_constant(m_, -d);
    384428    return *this;
    … …  
    388432  const matrix& matrix::operator*=(const matrix& other)
    389433  {
     434    assert(m_);
    390435    if (!blas_result_) {
    391436      blas_result_ = gsl_matrix_alloc(rows(),other.columns());
    … …  
    393438        throw utility::GSL_error("matrix::operator*= failed to allocate memory");
    394439    }
    395     gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0,
     440    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.proxy_m_, 0.0,
    396441                   blas_result_);
    397442    gsl_matrix* tmp=m_;
    398     m_=blas_result_;
     443    proxy_m_ = m_ = blas_result_;
    399444    blas_result_=tmp;
    400445    return *this;
    … …  
    404449  const matrix& matrix::operator*=(const double d)
    405450  {
     451    assert(m_);
    406452    gsl_matrix_scale(m_, d);
    407453    return *this;
    … …  
    409455
    410456
    411   bool isnull(const matrix& m)
    412   {
    413     return gsl_matrix_isnull(m.gsl_matrix_p());
    414   }
    415 
    416 
    417   double max(const matrix& m)
    418   {
    419     return gsl_matrix_max(m.gsl_matrix_p());
    420   }
    421 
    422 
    423   double min(const matrix& m)
    424   {
    425     return gsl_matrix_min(m.gsl_matrix_p());
    426   }
    427 
    428 
    429   void minmax_index(const matrix& m,
     457  bool isnull(const matrix& other)
     458  {
     459    return gsl_matrix_isnull(other.gsl_matrix_p());
     460  }
     461
     462
     463  double max(const matrix& other)
     464  {
     465    return gsl_matrix_max(other.gsl_matrix_p());
     466  }
     467
     468
     469  double min(const matrix& other)
     470  {
     471    return gsl_matrix_min(other.gsl_matrix_p());
     472  }
     473
     474
     475  void minmax_index(const matrix& other,
    430476                    std::pair<size_t,size_t>& min, std::pair<size_t,size_t>& max)
    431477  {
    432     gsl_matrix_minmax_index(m.gsl_matrix_p(), &min.first, &min.second,
     478    gsl_matrix_minmax_index(other.gsl_matrix_p(), &min.first, &min.second,
    433479                            &max.first, &max.second);
    434480  }
    … …  
    440486    size_t columns=templat.columns();
    441487    if (rows!=flag.rows() && columns!=flag.columns())
    442       flag=matrix(rows,columns,1.0);
     488      flag.clone(matrix(rows,columns,1.0));
    443489    else
    444490      flag.set_all(1.0);
    … …  
    456502  void swap(matrix& a, matrix& b)
    457503  {
     504    assert(v.gsl_matrix_p()); assert(w.gsl_matrix_p());
    458505    int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p());
    459506    if (status)

  • trunk/yat/utility/matrix.h

    r788 r792  
    6464  /// supported are no binary file support and views on arrays,
    6565  /// utility::vectors, gsl_vectors, diagonals, subdiagonals, and
    66   /// superdiagonals. If there is an interest from the user community,
    67   /// the omitted functionality could be included.
     66  /// superdiagonals.
    6867  ///
    6968  class matrix
    … …  
    143142    ~matrix(void);
    144143
     144    /**
     145       \brief Make a copy of \a other.
     146
     147       This function will make a deep copy of \a other. Memory is
     148       resized and view state is changed if needed.
     149    */
     150    const matrix& clone(const matrix& other);
     151
    145152    ///
    146153    /// @return The number of columns in the matrix.
    … …  
    155162       \throw GSL_error if dimensions mis-match.
    156163    */
    157     void div_elements(const matrix& b);
     164    void div(const matrix& b);
    158165
    159166    /**
    … …  
    193200       \throw GSL_error if dimensions mis-match.
    194201    */
    195     void mul_elements(const matrix& b);
     202    void mul(const matrix& b);
    196203
    197204    ///
    … …  
    320327       \brief The assignment operator.
    321328
    322        There is no requirements on dimensions, i.e. the matrix is
    323        remapped in memory if necessary. This implies that in general
    324        views cannot be assigned using this operator. Views will be
    325        mutated into normal matrices. The only exception to this
    326        behaviour on views is when self-assignemnt is done, since
    327        self-assignment is ignored.
     329       Dimensions of the matrices must match. If the LHS matrix is a
     330       view, the underlying data will be changed.
    328331
    329332       \return A const reference to the resulting matrix.
    … …  
    331334       \see void set(const matrix&).
    332335
    333        \throw A GSL_error is indirectly thrown if memory allocation
    334        fails, or if dimensions mis-match.
     336       \throw GSL_error if dimensions mis-match.
    335337    */
    336338    const matrix& operator=(const matrix& other);
    … …  
    419421    gsl_matrix* m_;
    420422    gsl_matrix_view* view_;
     423    const gsl_matrix_const_view* view_const_;
     424    // proxy_m_ is used to access the proper underlying gsl_matrix in
     425    // all const member functions. It is not used by design for
     426    // non-const vector functions and operators. This is to make sure
     427    // that runtime errors occur if a const matrix is used in an
     428    // inappropriate manner such as on left hand side in assignment
     429    // (remember, m_ is null for const matrix views).
     430    const gsl_matrix* proxy_m_;
    421431  };
    422432
Note: See TracChangeset for help on using the changeset viewer.