Changeset 413 for branches


Ignore:
Timestamp:
Nov 30, 2005, 10:35:19 PM (18 years ago)
Author:
Jari Häkkinen
Message:

All check passes, but more matrix checks needed.

Location:
branches/better_matrix_class
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • branches/better_matrix_class/lib/gslapi/matrix.cc

    r408 r413  
    1111
    1212#include <gsl/gsl_blas.h>
    13 #include <gsl/gsl_linalg.h>
    1413
    1514
    1615namespace theplu {
    1716namespace gslapi {
    18 
    19 
    20 
    21   matrix::matrix(void)
    22     : m_(NULL), view_(NULL)
    23   {
    24   }
    25 
    26 
    27 
    28   matrix::matrix(const size_t& r, const size_t& c, double init_value)
    29     : view_(NULL)
    30   {
    31     m_ = gsl_matrix_alloc(r,c);
    32     set_all(init_value);
    33   }
    34 
    35 
    36 
    37   matrix::matrix(const matrix& other)
    38     : view_(NULL)
    39   {
    40     m_ = other.gsl_matrix_copy();
    41   }
    4217
    4318
    … …  
    133108
    134109  // Jari, checkout GSL transpose support in GSL manual 8.4.9
    135   /*
    136   matrix matrix::transpose(void) const
     110  void matrix::transpose(void)
    137111  {
    138     matrix res(columns(),rows());
    139     gsl_matrix_transpose_memcpy(res.m_,m_);
    140     return res;
     112    if (columns()==rows())
     113      gsl_matrix_transpose(m_);
     114    else {
     115      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
     116      gsl_matrix_transpose_memcpy(transposed,m_);
     117      gsl_matrix_free(m_);
     118      m_=transposed;
     119    }
    141120  }
    142121
    143122
    144123
    145   vector matrix::TEMP_col_return(size_t column) const
    146   {
    147     vector vec(rows());
    148     for (size_t i=0; i<rows(); ++i)
    149       vec(i)=operator()(i,column);
    150     return vec;
    151   }
    152 
    153 
    154 
    155   vector matrix::operator[](size_t row) const
    156   {
    157     vector vec(columns());
    158     for (size_t i=0; i<columns(); ++i)
    159       vec(i)=operator()(row,i);
    160     return vec;
    161   }
    162 
    163 
    164 
    165   matrix matrix::operator*( const matrix &other ) const
    166   {
    167     matrix res(rows(), other.columns());
    168     gsl_linalg_matmult(m_,other.m_,res.m_);
    169     return res;
    170   }
    171 
    172 
    173 
     124  /*
    174125  vector matrix::operator*( const vector &other ) const
    175126  {
    … …  
    178129                   gslvec );
    179130    vector res(gslvec);
    180     return res;
    181   }
    182 
    183 
    184 
    185   matrix matrix::operator+( const matrix &other ) const
    186   {
    187     matrix res( *this );
    188     gsl_matrix_add( res.m_, other.m_ );
    189     return res;
    190   }
    191 
    192 
    193 
    194   matrix matrix::operator-( const matrix &other ) const
    195   {
    196     matrix res( *this );
    197     gsl_matrix_sub( res.m_, other.m_ );
    198131    return res;
    199132  }
    … …  
    208141      m_ = other.gsl_matrix_copy();
    209142    }
     143    return *this;
     144  }
     145
     146
     147
     148  const matrix& matrix::operator*=(const matrix& other)
     149  {
     150    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
     151    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
     152    gsl_matrix_free(m_);
     153    m_=result;
    210154    return *this;
    211155  }

  • branches/better_matrix_class/lib/gslapi/matrix.h

    r409 r413  
    4040  /// the omitted functionality could be included.
    4141  ///
     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  ///
    4247  class matrix
    4348  {
    … …  
    5055    /// structures.
    5156    ///
    52     matrix(void);
     57    inline matrix(void) : m_(NULL), view_(NULL) {}
    5358
    5459    ///
    … …  
    5661    /// elements, and sets all elements to \a init_value.
    5762    ///
    58     matrix(const size_t& r, const size_t& c, const double init_value=0);
     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); }
    5965
    6066    ///
    … …  
    6470    /// of the view will be copied, i.e. the view is not copied.
    6571    ///
    66     matrix(const matrix&);
     72    inline matrix(const matrix& o) : view_(NULL) { m_ = o.gsl_matrix_copy(); }
    6773
    6874    ///

  • branches/better_matrix_class/lib/gslapi/vector.cc

    r408 r413  
    1515
    1616
    17   vector::vector(void)
    18     : v_(NULL), view_(NULL)
    19   {
    20   }
    21 
    22 
    23 
    24   vector::vector(const size_t n,const double init_value)
    25     : view_(NULL)
    26   {
    27     v_ = gsl_vector_alloc(n);
    28     set_all(init_value);
    29   }
    30 
    31 
    32 
    33   vector::vector(const vector& other)
    34     : view_(NULL)
    35   {
    36     v_ = other.gsl_vector_copy();
    37   }
    38 
    39 
    40 
    4117  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
     18    : const_view_(NULL)
    4219  {
    4320    // Jari, exception handling needed here. Failure in setting up a
    … …  
    5128
    5229
     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
    5343  vector::vector(matrix& m, size_t i, bool row)
     44    : const_view_(NULL)
    5445  {
    5546    view_ = new gsl_vector_view(row ?
    … …  
    6152
    6253
    63 
    64   /* Jari, remove this?
    65   vector::vector(gsl_vector* vec)
    66     : v_(vec), view_(NULL)
    67   {
    68   }
    69   */
     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_pointer(),i) :
     59                             gsl_matrix_const_column(m.gsl_matrix_pointer(),i) );
     60    v_ = const_cast<gsl_vector*>(&(const_view_->vector));
     61  }
    7062
    7163
    7264
    7365  vector::vector(std::istream& is) throw (utility::IO_error,std::exception)
    74     : view_(NULL)
     66    : view_(NULL), const_view_(NULL)
    7567  {
    7668    // read the data file and store in stl vectors (dynamically
    … …  
    122114    if (view_)
    123115      delete view_;
     116    else if (const_view_)
     117      delete const_view_;
    124118    else if (v_)
    125119      gsl_vector_free(v_);

  • branches/better_matrix_class/lib/gslapi/vector.h

    r408 r413  
    4545    /// The default constructor.
    4646    ///
    47     vector(void);
     47    inline vector(void) : v_(NULL), view_(NULL), const_view_(NULL) {}
    4848
    4949    ///
    … …  
    5151    /// sets all elements to \a init_value.
    5252    ///
    53     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); }
    5456
    5557    ///
    … …  
    5961    /// of the view will be copied, i.e. the view is not copied.
    6062    ///
    61     vector(const vector&);
     63    inline vector(const vector& other) : view_(NULL), const_view_(NULL)
     64      { v_ = other.gsl_vector_copy(); }
    6265
    6366    ///
    … …  
    8184
    8285    ///
     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    ///
    83103    /// Matrix row/column view constructor.
    84104    ///
    … …  
    103123
    104124    ///
    105     /// Constructor that imports a GSL vector. The GSL object is owned
    106     /// by the created object.
    107     ///
    108     // Jari, remove this?
    109 //    vector(gsl_vector*);
     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);
    110144
    111145    ///
    … …  
    172206    /// @return True if the object is a view, false othwerwise.
    173207    ///
    174     inline bool isview(void) const { return view_; }
     208    inline bool isview(void) const { return view_ || const_view_; }
    175209
    176210    ///
    … …  
    414448    gsl_vector* v_;
    415449    gsl_vector_view* view_;
     450    gsl_vector_const_view* const_view_;
    416451  };
    417452

  • branches/better_matrix_class/lib/svm/InputRanker.cc

    r301 r413  
    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

  • branches/better_matrix_class/lib/svm/InputRanker.h

    r295 r413  
    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 

  • branches/better_matrix_class/lib/svm/Kernel_MEV.cc

    r336 r413  
    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

  • branches/better_matrix_class/lib/svm/Kernel_MEV.h

    r350 r413  
    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    ///

  • branches/better_matrix_class/lib/svm/Kernel_SEV.cc

    r336 r413  
    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

  • branches/better_matrix_class/lib/svm/Kernel_SEV.h

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

  • branches/better_matrix_class/lib/utility/PCA.cc

    r301 r413  
    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

  • branches/better_matrix_class/lib/utility/PCA.h

    r334 r413  
    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 

  • branches/better_matrix_class/lib/utility/SVD.cc

    r301 r413  
    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

  • branches/better_matrix_class/lib/utility/SVD.h

    r303 r413  
    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    ///
    … …  
    9899    /// is undefined.
    99100    ///
    100     gslapi::matrix U(void) const { return U_; }
     101    inline const gslapi::matrix& U(void) const { return U_; }
    101102
    102103    ///
    … …  
    108109    /// is undefined.
    109110    ///
    110     gslapi::matrix V(void) const { return V_; }
     111    inline const gslapi::matrix& V(void) const { return V_; }
    111112
    112113  private:

  • branches/better_matrix_class/lib/utility/stl_utility.cc

    r341 r413  
    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

  • branches/better_matrix_class/test/roc_test.cc

    r342 r413  
    7070  const double tol = 0.001;
    7171  for (size_t i=0; i<data.rows(); i++){
    72     gslapi::vector vec = data[i];
     72    gslapi::vector vec(data,i);
    7373    area = roc.score(target2,vec);
    7474    if (area<correct_area(i)-tol || area>correct_area(i)+tol){
    … …  
    8383  gslapi::vector weight(target2.size(),1);
    8484  for (size_t i=0; i<data.rows(); i++){
    85     gslapi::vector vec = data[i];
     85    gslapi::vector vec(data,i);
    8686    area = roc.score(target2, vec, weight);
    8787    if (area<correct_area(i)-tol || area>correct_area(i)+tol){

  • branches/better_matrix_class/test/svd_test.cc

    r377 r413  
    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
Note: See TracChangeset for help on using the changeset viewer.