Ignore:
Timestamp:
Jul 17, 2008, 12:36:23 AM (13 years ago)
Author:
Peter
Message:

working on #363

File:
1 copied

Legend:

Unmodified
Added
Removed
  • trunk/yat/utility/MatrixWeighted.h

    r1360 r1380  
    1 #ifndef _theplu_yat_utility_matrix_
    2 #define _theplu_yat_utility_matrix_
     1#ifndef _theplu_yat_utility_matrix_weighted_
     2#define _theplu_yat_utility_matrix_weighted_
    33
    44// $Id$
     
    2828*/
    2929
    30 #include "Exception.h"
     30#include "DataWeight.h"
    3131#include "StrideIterator.h"
    32 #include "Vector.h"
    33 #include "VectorConstView.h"
    34 #include "VectorView.h"
    35 
    36 #include <gsl/gsl_matrix.h>
    37 #include <iostream>
    38 #include <utility>
     32
     33#include <vector>
    3934
    4035namespace theplu {
     
    4237namespace utility {
    4338
    44   class VectorBase;
    45 
    4639  ///
    47   /// @brief Interface to GSL matrix.
     40  /// @brief Weighted Matrix
    4841  ///
    49   /// For the time being 'double' is the only type supported.
    50   ///
    51   /// \par[File streams] Reading and writing vectors to file streams
    52   /// are of course supported. These are implemented without using GSL
    53   /// functionality, and thus binary read and write to streams are not
    54   /// supported.
    55   ///
    56   /// @note All GSL matrix related functions are not implement but
    57   /// most functionality defined for GSL matrices can be achieved with
    58   /// this interface class. Most notable GSL functionality not
    59   /// supported are no binary file support and views on arrays,
    60   /// utility::vectors, gsl_vectors, diagonals, subdiagonals, and
    61   /// superdiagonals.
    62   ///
    63   class Matrix
     42  class MatrixWeighted
    6443  {
    6544  public:
     
    6746       Mutable iterator that iterates over all elements
    6847     */
    69     typedef StrideIterator<double*> iterator;
     48    typedef StrideIterator<std::vector<DataWeight>::iterator> iterator;
    7049
    7150    /**
    7251       Read-only iterator that iterates over all elements
    7352     */
    74     typedef StrideIterator<const double*> const_iterator;
     53    typedef StrideIterator<std::vector<DataWeight>::const_iterator>
     54    const_iterator;
    7555
    7656    /**
    7757       Mutable iterator that iterates over one column
    7858     */
    79     typedef StrideIterator<double*> column_iterator;
     59    typedef StrideIterator<std::vector<DataWeight>::iterator> column_iterator;
    8060
    8161    /**
    8262       Read-only iterator that iterates over one column
    8363     */
    84     typedef StrideIterator<const double*> const_column_iterator;
     64    typedef StrideIterator<std::vector<DataWeight>::const_iterator>
     65    const_column_iterator;
    8566
    8667    /**
    8768       Mutable iterator that iterates over one row
    8869     */
    89     typedef StrideIterator<double*> row_iterator;
     70    typedef StrideIterator<std::vector<DataWeight>::iterator> row_iterator;
    9071
    9172    /**
    9273       Read-only iterator that iterates over one row
    9374     */
    94     typedef StrideIterator<const double*> const_row_iterator;
     75    typedef StrideIterator<std::vector<DataWeight>::const_iterator>
     76    const_row_iterator;
    9577
    9678    /**
     
    10082       structures.
    10183    */
    102     Matrix(void);
     84    MatrixWeighted(void);
    10385
    10486    /**
    10587       \brief Constructor allocating memory space for \a r times \a c
    106        elements, and sets all elements to \a init_value.
    107 
    108        \throw GSL_error if memory allocation fails.
    109     */
    110     Matrix(const size_t& r, const size_t& c, double init_value=0);
     88       elements, and sets all elements to
     89       DataWeight(\a init_value, \a init_weight)
     90    */
     91    MatrixWeighted(const size_t& r, const size_t& c, double init_value=0,
     92                 double init_weight=1.0);
    11193
    11294    /**
    11395       \brief The copy constructor.
    114 
    115        \throw A GSL_error is indirectly thrown if memory allocation
    116        fails.
    117     */
    118     Matrix(const Matrix&);
     96    */
     97    MatrixWeighted(const MatrixWeighted&);
    11998
    12099    /**
    121100       \brief The istream constructor.
    122101
    123        Matrix rows are sepearated with the new line character, and
    124        column elements in a row must be separated either with white
    125        space characters except the new line (\\n) character (default),
    126        or by the delimiter \a sep.  Rows, and columns, are read
    127        consecutively and only rectangular matrices are
    128        supported. Empty lines are ignored. End of input is at end of
    129        file marker.
    130 
    131        \throw GSL_error if memory allocation fails, IO_error if
    132        unexpected input is found in the input stream.
    133     */
    134     explicit Matrix(std::istream &, char sep='\0')
    135       throw(utility::IO_error, std::exception);
    136 
    137     /**
    138        \brief The destructor.
    139     */
    140     ~Matrix(void);
    141 
    142     ///
    143     /// Set all elements to \a value.
    144     ///
    145     void all(const double value);
     102       Constructor is reading data values using corresponding
     103       constructor in Matrix class. The weights are calculate dusing
     104       the nan function.
     105    */
     106    explicit MatrixWeighted(std::istream &, char sep='\0');
    146107
    147108    /**
     
    190151
    191152    /**
    192        \return Vector view of column \a i
    193      */
    194     VectorView column_view(size_t i);
    195 
    196     /**
    197        \return const vector view of column \a i
    198      */
    199     const VectorConstView column_const_view(size_t) const;
    200 
    201     ///
    202     /// @return The number of columns in the matrix.
    203     ///
     153       \return The number of columns in the matrix.
     154    */
    204155    size_t columns(void) const;
    205156
    206157    /**
    207        Elementwise division of the elemnts of the calling matrix by
    208        the elements of matrix \a b, \f$ a_{ij} = a_{ij} / b_{ij} \;
    209        \forall i,j \f$. The result is stored into the calling matrix.
    210 
    211        \throw GSL_error if dimensions mis-match.
    212     */
    213     void div(const Matrix& b);
    214 
    215     /**
    216158       \return iterator pointing to end of matrix
    217159     */
     
    244186
    245187    /**
    246        \brief Check whether matrices are equal within a user defined
    247        precision, set by \a precision.
    248 
    249        \return True if each element deviates less or equal than \a
    250        d. If any matrix contain a NaN, false is always returned.
    251 
    252        \see operator== and operator!=
    253     */
    254     bool equal(const Matrix&, const double precision=0) const;
    255 
    256     ///
    257     /// @return A const pointer to the internal GSL matrix.
    258     ///
    259     const gsl_matrix* gsl_matrix_p(void) const;
    260 
    261     ///
    262     /// @return A pointer to the internal GSL matrix.
    263     ///
    264     gsl_matrix* gsl_matrix_p(void);
    265 
    266     /**
    267        Multiply the elements of Matrix \a b with the elements of the
    268        calling Matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall i,j
    269        \f$. The result is stored into the calling Matrix.
    270 
    271        \throw GSL_error if dimensions mis-match.
    272     */
    273     void mul(const Matrix& b);
    274 
    275     /**
    276188       \brief Resize Matrix
    277        
    278        All elements are set to @a init_value.
    279 
    280        \note underlying GSL matrix is destroyed and views into this
    281        Matrix becomes invalid.
    282     */
    283     void resize(size_t, size_t, double init_value=0);
    284 
    285     ///
    286     /// @return The number of rows in the matrix.
    287     ///
     189
     190       \note this function may invalidate iterators.
     191    */
     192    void resize(size_t, size_t);
     193
     194    /**
     195       \return The number of rows in the matrix.
     196    */
    288197    size_t rows(void) const;
    289198
    290199    /**
    291        \return Vector view of \a i
    292      */
    293     VectorView row_view(size_t);
    294 
    295     /**
    296        \return const vector view of \a i
    297      */
    298     const VectorConstView row_const_view(size_t) const;
    299 
    300     /**
    301        \brief Swap columns \a i and \a j.
    302 
    303        \throw GSL_error if either index is out of bounds.
    304     */
    305     void swap_columns(const size_t i,const size_t j);
    306 
    307     /**
    308        \brief Swap row \a i and column \a j.
    309 
    310        The Matrix must be square.
    311 
    312        \throw GSL_error if either index is out of bounds, or if matrix
    313        is not square.
    314     */
    315     void swap_rowcol(const size_t i,const size_t j);
    316 
    317     /**
    318        \brief Swap rows \a i and \a j.
    319 
    320        \throw GSL_error if either index is out of bounds.
    321     */
    322     void swap_rows(const size_t i, const size_t j);
    323 
    324     /**
    325200       \brief Transpose the matrix.
    326201
    327        \throw GSL_error if memory allocation fails for the new
    328        transposed matrix.
     202       \note This function may invalidate iterators.
    329203    */
    330204    void transpose(void);
     
    334208
    335209       \return Reference to the element position (\a row, \a column).
    336 
    337        \throw If GSL range checks are enabled in the underlying GSL
    338        library a GSL_error exception is thrown if either index is out
    339        of range.
    340     */
    341     double& operator()(size_t row,size_t column);
     210    */
     211    DataWeight& operator()(size_t row,size_t column);
    342212
    343213    /**
     
    346216       \return Const reference to the element position (\a row, \a
    347217       column).
    348 
    349        \throw If GSL range checks are enabled in the underlying GSL
    350        library a GSL_error exception is thrown if either index is out
    351        of range.
    352     */
    353     const double& operator()(size_t row,size_t column) const;
    354 
    355     /**
    356        \brief Comparison operator. Takes squared time.
    357 
    358        Checks are performed with exact matching, i.e., rounding off
    359        effects may destroy comparison. Use the equal function for
    360        comparing elements within a user defined precision.
    361 
    362        \return True if all elements are equal otherwise false.
    363 
    364        \see equal
    365     */
    366     bool operator==(const Matrix& other) const;
    367 
    368     /**
    369        \brief Comparison operator. Takes squared time.
    370 
    371        Checks are performed with exact matching, i.e., rounding off
    372        effects may destroy comparison. Use the equal function for
    373        comparing elements within a user defined precision.
    374 
    375        \return False if all elements are equal otherwise true.
    376 
    377        \see equal
    378     */
    379     bool operator!=(const Matrix& other) const;
     218    */
     219    const DataWeight& operator()(size_t row,size_t column) const;
     220
    380221
    381222    /**
     
    384225       \return A const reference to the resulting Matrix.
    385226    */
    386     const Matrix& operator=(const Matrix& other);
    387 
    388     /**
    389        \brief Add and assign operator.
    390 
    391        Elementwise addition of the elements of Matrix \a b to the left
    392        hand side Matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall i,j
    393        \f$.
    394 
    395        \return A const reference to the resulting Matrix.
    396 
    397        \throw GSL_error if dimensions mis-match.
    398     */
    399     const Matrix& operator+=(const Matrix& b);
    400 
    401     /**
    402        \brief Add and assign operator
    403 
    404        Add the scalar value \a d to the left hand side Matrix, \f$
    405        a_{ij} = a_{ij} + d \; \forall i,j \f$.
    406     */
    407     const Matrix& operator+=(const double d);
    408 
    409     /**
    410        \brief Subtract and assign operator.
    411 
    412        Elementwise subtraction of the elements of Matrix \a b to the
    413        left hand side Matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall
    414        i,j \f$.
    415 
    416        \return A const reference to the resulting Matrix.
    417 
    418        \throw GSL_error if dimensions mis-match.
    419     */
    420     const Matrix& operator-=(const Matrix&);
    421 
    422     /**
    423        \brief Subtract and assign operator
    424 
    425        Subtract the scalar value \a d to the left hand side Matrix,
    426        \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$.
    427     */
    428     const Matrix& operator-=(const double d);
    429 
    430     /**
    431        \brief Multiply and assigment operator.
    432 
    433        \return Const reference to the resulting Matrix.
    434 
    435        \throw GSL_error if memory allocation fails.
    436     */
    437     const Matrix& operator*=(const Matrix&);
    438 
    439     /**
    440        \brief Multiply and assignment operator
    441 
    442        Multiply the elements of the left hand side Matrix with a
    443        scalar \a d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$.
    444 
    445        \throw GSL_error if memory allocation fails.
    446     */
    447     const Matrix& operator*=(double d);
     227    const MatrixWeighted& operator=(const MatrixWeighted& other);
    448228
    449229  private:
    450 
    451     /**
    452        \brief Create a new copy of the internal GSL matrix.
    453 
    454        Necessary memory for the new GSL matrix is allocated and the
    455        caller is responsible for freeing the allocated memory.
    456 
    457        \return A pointer to a copy of the internal GSL matrix.
    458 
    459        \throw GSL_error if memory cannot be allocated for the new
    460        copy, or if dimensions mis-match.
    461     */
    462     gsl_matrix* create_gsl_matrix_copy(void) const;
    463 
    464     /**
    465        \brief Clear all dynamically allocated memory.
    466 
    467        Internal utility function.
    468     */
    469     void delete_allocated_memory(void);
    470 
    471     // blas_result_ is used to temporarily store result in BLAS calls
    472     // and when not NULL it should always have the same size as
    473     // m_. Memory is not allocated for blas_result_ until it is
    474     // needed.
    475     gsl_matrix* blas_result_;
    476     gsl_matrix* m_;
     230    std::vector<DataWeight> vec_;
     231    size_t columns_;
     232
    477233  };
    478234
    479235  /**
    480      \brief Check if all elements of the Matrix are zero.
    481 
    482      \return True if all elements in the Matrix is zero, false
    483      othwerwise.
     236     \brief Exchange all elements.
     237
     238     Takes constant time.
     239
     240     \throw std::runtime_error if the two matrices disagree in size.
    484241  */
    485   bool isnull(const Matrix&);
    486 
    487   /**
    488      \brief Get the maximum value of the Matrix.
    489 
    490      \return The maximum value of the Matrix.
    491   */
    492   double max(const Matrix&);
    493 
    494   /**
    495      \brief Get the minimum value of the Matrix.
    496 
    497      \return The minimum value of the Matrix.
    498   */
    499   double min(const Matrix&);
    500 
    501   /**
    502      \brief Locate the maximum and minumum element in the Matrix.
    503 
    504      \return The indecies to the element with the minimum and maximum
    505      values of the Matrix, respectively.
    506 
    507      \note Lower index has precedence (searching in row-major order).
    508   */
    509   void minmax_index(const Matrix&,
    510                     std::pair<size_t,size_t>& min,
    511                     std::pair<size_t,size_t>& max);
    512 
    513   /**
    514      \brief Create a Matrix \a flag indicating NaN's in another Matrix
    515      \a templat.
    516 
    517      The \a flag Matrix is changed to contain 1's and 0's only. A 1
    518      means that the corresponding element in the \a templat Matrix is
    519      valid and a zero means that the corresponding element is a NaN.
    520 
    521      \note Space for Matrix \a flag is reallocated to fit the size of
    522      Matrix \a templat if sizes mismatch.
    523 
    524      \return True if the \a templat Matrix contains at least one NaN.
    525   */
    526   bool nan(const Matrix& templat, Matrix& flag);
    527 
    528   /**
    529      \brief Exchange all elements between the matrices by copying.
    530 
    531      The two matrices must have the same size.
    532 
    533      \throw GSL_error if either index is out of bounds.
    534   */
    535   void swap(Matrix&, Matrix&);
    536 
    537   /**
    538      \brief The output operator for the Matrix class.
    539   */
    540   std::ostream& operator<< (std::ostream& s, const Matrix&);
    541 
    542   /**
    543      \brief Vector Matrix multiplication
    544    */
    545   Vector operator*(const Matrix&, const VectorBase&);
    546 
    547   /**
    548      \brief Matrix Vector multiplication
    549    */
    550   Vector operator*(const VectorBase&, const Matrix&);
     242  void swap(MatrixWeighted&, MatrixWeighted&);
    551243
    552244}}} // of namespace utility, yat, and theplu
Note: See TracChangeset for help on using the changeset viewer.