mtx.h

00001 
00004 
00005 //
00006 // Copyright (c) 2001-2004 California Institute of Technology
00007 // Copyright (c) 2004-2007 University of Southern California
00008 // Rob Peters <rjpeters at usc dot edu>
00009 //
00010 // created: Mon Mar 12 12:23:11 2001
00011 // commit: $Id: mtx.h 10065 2007-04-12 05:54:56Z rjpeters $
00012 // $HeadURL: file:///lab/rjpeters/svnrepo/code/trunk/groovx/src/mtx/mtx.h $
00013 //
00014 // --------------------------------------------------------------------
00015 //
00016 // This file is part of GroovX.
00017 //   [http://ilab.usc.edu/rjpeters/groovx/]
00018 //
00019 // GroovX is free software; you can redistribute it and/or modify it
00020 // under the terms of the GNU General Public License as published by
00021 // the Free Software Foundation; either version 2 of the License, or
00022 // (at your option) any later version.
00023 //
00024 // GroovX is distributed in the hope that it will be useful, but
00025 // WITHOUT ANY WARRANTY; without even the implied warranty of
00026 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00027 // General Public License for more details.
00028 //
00029 // You should have received a copy of the GNU General Public License
00030 // along with GroovX; if not, write to the Free Software Foundation,
00031 // Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
00032 //
00034 
00035 #ifndef GROOVX_PKGS_MTX_MTX_H_UTC20050626084022_DEFINED
00036 #define GROOVX_PKGS_MTX_MTX_H_UTC20050626084022_DEFINED
00037 
00038 #include "datablock.h"
00039 
00040 #include "arithfunctor.h"
00041 
00042 #include <iosfwd>
00043 #include <iterator>
00044 
00045 namespace rutz
00046 {
00047   class fstring;
00048 }
00049 
00050 //  #######################################################
00051 //  =======================================================
00053 
00054 namespace range_checking
00055 {
00056   void geq(const void* x, const void* lim, const char* f, int ln);
00057   void lt(const void* x, const void* lim, const char* f, int ln);
00058   void leq(const void* x, const void* lim, const char* f, int ln);
00059   void in_half_open(const void* x, const void* llim, const void* ulim,
00060                     const char* f, int ln);
00061   void in_full_open(const void* x, const void* llim, const void* ulim,
00062                     const char* f, int ln);
00063 
00064   void geq(int x, int lim, const char* f, int ln);
00065   void lt(int x, int lim, const char* f, int ln);
00066   void leq(int x, int lim, const char* f, int ln);
00067   void in_half_open(int x, int llim, int ulim, const char* f, int ln);
00068   void in_full_open(int x, int llim, int ulim, const char* f, int ln);
00069 }
00070 
00071 #ifdef RANGE_CHECK
00072 // Range check
00073 #  define RC_geq(x,lim) range_checking::geq((x),(lim),__FILE__,__LINE__)
00074 #  define RC_less(x,lim) range_checking::lt((x),(lim),__FILE__,__LINE__)
00075 #  define RC_leq(x,lim) range_checking::leq((x),(lim),__FILE__,__LINE__)
00076 #  define RC_in_half_open(x,llim,ulim) range_checking::in_half_open((x),(llim),(ulim),__FILE__,__LINE__)
00077 #  define RC_in_full_open(x,llim,ulim) range_checking::in_full_open((x),(llim),(ulim),__FILE__,__LINE__)
00078 
00079 // Range check, and return the checked value
00080 #  define RCR_geq(x,lim) (range_checking::geq((x),(lim),__FILE__,__LINE__), x)
00081 #  define RCR_less(x,lim) (range_checking::lt((x),(lim),__FILE__,__LINE__), x)
00082 #  define RCR_leq(x,lim) (range_checking::leq((x),(lim),__FILE__,__LINE__), x)
00083 #  define RCR_in_half_open(x,llim,ulim) (range_checking::in_half_open((x),(llim),(ulim),__FILE__,__LINE__), x)
00084 #  define RCR_in_full_open(x,llim,ulim) (range_checking::in_full_open((x),(llim),(ulim),__FILE__,__LINE__), x)
00085 
00086 #else // !RANGE_CHECK
00087 
00088 #  define RC_geq(x,lim)
00089 #  define RC_less(x,lim)
00090 #  define RC_leq(x,lim)
00091 #  define RC_in_half_open(x,llim,ulim)
00092 #  define RC_in_full_open(x,llim,ulim)
00093 
00094 #  define RCR_geq(x,lim) (x)
00095 #  define RCR_less(x,lim) (x)
00096 #  define RCR_leq(x,lim) (x)
00097 #  define RCR_in_half_open(x,llim,ulim) (x)
00098 #  define RCR_in_full_open(x,llim,ulim) (x)
00099 #endif
00100 
00101 
00102 //  #######################################################
00103 //  #######################################################
00104 //  =======================================================
00105 //  classes for representing ranges of indices
00106 
00108 
00109 class index_range
00110 {
00111 private:
00112   const int m_first;
00113   const int m_count;
00114 
00115 public:
00116   index_range(int first, int count) : m_first(first), m_count(count) {}
00117 
00118   int begin() const { return m_first; }
00119   int end() const { return m_first+m_count; }
00120   int count() const { return m_count; }
00121 };
00122 
00124 class row_index_range : public index_range
00125 {
00126 public:
00128   row_index_range(int first, int count) : index_range(first, count) {}
00129 };
00130 
00132 class col_index_range : public index_range
00133 {
00134 public:
00136   col_index_range(int first, int count) : index_range(first, count) {}
00137 };
00138 
00139 inline index_range range(int i) { return index_range(i, 1); }
00140 inline index_range range(int begin, int end) { return index_range(begin, end-begin); }
00141 inline index_range range_n(int begin, int count) { return index_range(begin, count); }
00142 
00143 inline row_index_range row_range(int r) { return row_index_range(r, 1); }
00144 inline row_index_range row_range(int begin, int end) { return row_index_range(begin, end-begin); }
00145 inline row_index_range row_range_n(int begin, int count) { return row_index_range(begin, count); }
00146 
00147 inline col_index_range col_range(int c) { return col_index_range(c, 1); }
00148 inline col_index_range col_range(int begin, int end) { return col_index_range(begin, end-begin); }
00149 inline col_index_range col_range_n(int begin, int count) { return col_index_range(begin, count); }
00150 
00151 
00152 
00153 
00154 //  #######################################################
00155 //  #######################################################
00156 //  =======================================================
00157 //  iterators
00158 
00159 
00160 class mtx;
00161 
00163 template <class T>
00164 class stride_iterator_base
00165 {
00166 private:
00167   T* data;
00168   int stride;
00169   T* stop;
00170 
00171 protected:
00172   // these friend declarations are to allow access to the protected
00173   // "raw" constructor that starts from a raw pointer plus
00174   // stride+length info
00175   friend class slice;
00176   friend class mtx;
00177   template <class U> friend class stride_iterator_base;
00178 
00179   stride_iterator_base(T* d, int str, int n) :
00180     data(d), stride(str), stop(data+str*n) {}
00181 
00182 public:
00183   typedef std::random_access_iterator_tag iterator_category;
00184 
00185   stride_iterator_base(const stride_iterator_base& other) :
00186     data(other.data), stride(other.stride), stop(other.stop) {}
00187 
00188   template <class U>
00189   explicit stride_iterator_base(const stride_iterator_base<U>& other) :
00190     data(other.data), stride(other.stride), stop(other.stop) {}
00191 
00192   typedef T            value_type;
00193   typedef ptrdiff_t    difference_type;
00194   typedef T*           pointer;
00195   typedef T&           reference;
00196 
00197   stride_iterator_base end() const { stride_iterator_base e(*this); e.data = stop; return e; }
00198 
00199   bool has_more() const { return data < stop; }
00200 
00201   // Dereference
00202 
00203   reference operator*() const { RC_less(data, stop); return *data; }
00204 
00205   // Comparison
00206 
00207   bool operator==(const stride_iterator_base& other) const { return data == other.data; }
00208 
00209   bool operator!=(const stride_iterator_base& other) const { return data != other.data; }
00210 
00211   bool operator<(const stride_iterator_base& other) const { return data < other.data; }
00212 
00213   difference_type operator-(const stride_iterator_base& other) const
00214     { return (data - other.data) / stride; }
00215 
00216   // Increment/Decrement
00217 
00218   stride_iterator_base& operator++() { data += stride; return *this; }
00219   stride_iterator_base& operator--() { data -= stride; return *this; }
00220 
00221   stride_iterator_base operator++(int) { stride_iterator_base cpy(*this); data += stride; return cpy; }
00222   stride_iterator_base operator--(int) { stride_iterator_base cpy(*this); data -= stride; return cpy; }
00223 
00224   stride_iterator_base& operator+=(int x) { data += x*stride; return *this; }
00225   stride_iterator_base& operator-=(int x) { data -= x*stride; return *this; }
00226 
00227   stride_iterator_base operator+(int x) const { stride_iterator_base res(*this); res += x; return res; }
00228   stride_iterator_base operator-(int x) const { stride_iterator_base res(*this); res -= x; return res; }
00229 };
00230 
00231 typedef stride_iterator_base<double> mtx_iter;
00232 typedef stride_iterator_base<const double> mtx_const_iter;
00233 
00234 
00235 //  #######################################################
00236 //  =======================================================
00237 
00239 class slice
00240 {
00241 private:
00242   data_holder&          m_data_source;
00243   ptrdiff_t    const    m_offset;
00244   int          const    m_stride;
00245   int          const    m_nelems;
00246 
00247   const double* data_start() const
00248     { return m_data_source.storage() + m_offset; }
00249 
00250   double* data_start_nc()
00251     { return m_data_source.storage_nc() + m_offset; }
00252 
00253   ptrdiff_t data_offset(int i) const { return m_stride*i; }
00254   const double* address(int i) const { return data_start() + data_offset(i); }
00255 
00256   ptrdiff_t storage_offset(int i) const { return m_offset + m_stride*i; }
00257 
00258   slice(const data_holder& src, ptrdiff_t offset, int s, int n)
00259     :
00260     // by doing this const_cast, slice is promising to uphold the const
00261     // correctness of the data_holder
00262     m_data_source(const_cast<data_holder&>(src)),
00263     m_offset(offset),
00264     m_stride(s),
00265     m_nelems(n)
00266   {}
00267 
00268   friend class mtx; // so that mtx can get at slice's constructor
00269 
00270 public:
00271 
00272   // Default copy constructor OK
00273 
00274   // No "default" assignment operator, since that would be assignment
00275   // of reference; instead we have an operator=() that does assignment
00276   // of value; see below
00277 
00278   //
00279   // Data access
00280   //
00281 
00282   double operator[](int i) const
00283   {
00284     RC_in_half_open(i, 0, m_nelems);
00285     return *(address(i));
00286   }
00287 
00288   int nelems() const { return m_nelems; }
00289 
00290   slice operator()(const index_range& rng) const;
00291 
00293   void print(std::ostream& s) const;
00294 
00296   void print_stdout() const;
00297 
00298   //
00299   // Iteration
00300   //
00301 
00302   mtx_iter begin_nc()
00303     { return mtx_iter(data_start_nc(), m_stride, m_nelems); }
00304 
00305   mtx_iter end_nc()
00306     { return begin_nc().end(); }
00307 
00308   mtx_const_iter begin() const
00309     { return mtx_const_iter(data_start(), m_stride, m_nelems); }
00310 
00311   mtx_const_iter end() const
00312     { return begin().end(); }
00313 
00314   //
00315   // Const functions
00316   //
00317 
00318   double sum() const;
00319   double min() const;
00320   double max() const;
00321 
00322   double mean() const { return sum()/m_nelems; }
00323 
00324 
00325   // Returns an index matrix so that this->reorder(mtx) will put the
00326   // slice in sorted order
00327   mtx get_sort_order() const;
00328 
00329   bool operator==(const slice& other) const;
00330 
00331   bool operator!=(const slice& other) const
00332   { return !(operator==(other)); }
00333 
00334   //
00335   // Non-const functions
00336   //
00337 
00338   template <class F>
00339   void apply(F func)
00340   {
00341     for (mtx_iter i = begin_nc(); i.has_more(); ++i)
00342       *i = func(*i);
00343   }
00344 
00345   void sort();
00346 
00347   // Reorders the slice by applying *this[i] = *this[index[i]]
00348   void reorder(const mtx& index);
00349 
00350   void operator+=(double val) { apply(dash::add(val)); }
00351   void operator-=(double val) { apply(dash::sub(val)); }
00352 
00353   void operator*=(double factor) { apply(dash::mul(factor)); }
00354   void operator/=(double div) { apply(dash::div(div)); }
00355 
00356   slice& operator+=(const slice& other);
00357   slice& operator-=(const slice& other);
00358 
00360   slice& operator=(double val);
00361 
00363   slice& operator=(const slice& other);
00364 
00366   slice& operator=(const mtx& other);
00367 };
00368 
00369 
00370 
00371 //  #######################################################
00372 //  =======================================================
00374 
00375 class mtx_shape
00376 {
00377 public:
00378   mtx_shape(int mr, int nc) : m_mrows(mr), m_ncols(nc) {}
00379 
00380   int mrows() const { return m_mrows; }
00381   int ncols() const { return m_ncols; }
00382 
00383   int max_dim() const { return (m_mrows > m_ncols) ? m_mrows : m_ncols; }
00384 
00385   int nelems() const { return m_mrows*m_ncols; }
00386 
00387 private:
00388   int m_mrows;
00389   int m_ncols;
00390 };
00391 
00392 
00393 //  #######################################################
00394 //  =======================================================
00396 
00397 class mtx_specs
00398 {
00399 public:
00400   mtx_specs() : m_shape(0, 0), m_rowstride(0), m_offset(0) {}
00401 
00402   mtx_specs(const mtx_shape& shape) :
00403     m_shape(shape),
00404     m_rowstride(shape.mrows()),
00405     m_offset(0)
00406   {}
00407 
00408   mtx_specs(int mrows, int ncols) :
00409     m_shape(mrows, ncols),
00410     m_rowstride(mrows),
00411     m_offset(0)
00412   {}
00413 
00414   mtx_specs(const mtx_specs& other) :
00415     m_shape(other.m_shape),
00416     m_rowstride(other.m_rowstride),
00417     m_offset(other.m_offset)
00418   {}
00419 
00420   void swap(mtx_specs& other);
00421 
00422   mtx_specs as_shape(const mtx_shape& s) const;
00423 
00424   mtx_specs as_shape(int mr, int nc) const { return as_shape(mtx_shape(mr,nc)); }
00425 
00426   void select_rows(const row_index_range& rng);
00427   void select_cols(const col_index_range& rng);
00428 
00429   mtx_specs sub_rows(const row_index_range& rng) const
00430   {
00431     mtx_specs result(*this); result.select_rows(rng); return result;
00432   }
00433 
00434   mtx_specs sub_cols(const col_index_range& rng) const
00435   {
00436     mtx_specs result(*this); result.select_cols(rng); return result;
00437   }
00438 
00439   const mtx_shape& shape() const { return m_shape; }
00440 
00441   ptrdiff_t offset() const { return m_offset; }
00442 
00443   int max_dim() const { return m_shape.max_dim(); }
00444   int nelems() const { return mrows()*ncols(); }
00445 
00446   int mrows() const { return m_shape.mrows(); }
00447   int rowstride() const { return m_rowstride; }
00448 
00449   int ncols() const { return m_shape.ncols(); }
00450   int colstride() const { return m_colstride; }
00451 
00452   unsigned int rowgap() const { return m_rowstride - mrows(); }
00453 
00454   int offset_from_start(int row, int col) const
00455   {
00456     // Strictly speaking, the only valid offsets are those that pass
00457     // RC_in_half_open(), but in order to allow "one-past-the-end"
00458     // indexing for iterators etc., we use the more permissive
00459     // RC_in_full_open() instead:
00460     RC_in_full_open(row, 0, mrows());
00461     RC_in_full_open(col, 0, ncols());
00462 
00463     return row + (col*m_rowstride);
00464   }
00465 
00466   int offset_from_start(int elem) const
00467   { return offset_from_start(elem%mrows(), elem/mrows()); }
00468 
00469   int offset_from_storage(int row, int col) const
00470   { return m_offset + offset_from_start(row, col); }
00471 
00472   int offset_from_storage(int elem) const
00473   { return m_offset + offset_from_start(elem); }
00474 
00475 private:
00476   mtx_shape m_shape;
00477   int m_rowstride;
00478   static const int m_colstride = 1;
00479   ptrdiff_t m_offset;
00480 };
00481 
00482 
00483 //  #######################################################
00484 //  #######################################################
00485 //  =======================================================
00486 //  Iterator templates for iterating over two-dimensional data
00487 
00488 
00489 //  #######################################################
00490 //  =======================================================
00492 
00493 template <class M, class T>
00494 class index_iterator_base
00495 {
00496   M* m_src;
00497   int m_index;
00498 
00499 public:
00500   index_iterator_base(M* m, int e) : m_src(m), m_index(e) {}
00501 
00502   typedef std::random_access_iterator_tag iterator_category;
00503 
00504   typedef T            value_type;
00505   typedef ptrdiff_t    difference_type;
00506   typedef T*           pointer;
00507   typedef T&           reference;
00508 
00509   index_iterator_base end() const { return index_iterator_base(m_src, m_src->nelems()); }
00510 
00511   bool has_more() const { return m_index < m_src->nelems(); }
00512 
00513 private:
00514   // Need this pair of overloads to distinguish between const and
00515   // non-const M types, so that we can call either at() or at_nc()
00516   // as appropriate.
00517   template <class MM>
00518   static reference get_at(MM* m, int e) { return m->at_nc(e); }
00519 
00520   template <class MM>
00521   static reference get_at(const MM* m, int e) { return m->at(e); }
00522 
00523 public:
00524   reference operator*() const { return get_at(m_src, m_index); }
00525 
00526   // Comparison
00527 
00528   bool operator==(const index_iterator_base& other) const { return m_index == other.m_index; }
00529 
00530   bool operator!=(const index_iterator_base& other) const { return m_index != other.m_index; }
00531 
00532   bool operator<(const index_iterator_base& other) const { return m_index < other.m_index; }
00533 
00534   difference_type operator-(const index_iterator_base& other) const
00535   { return m_index - other.m_index; }
00536 
00537   // Increment/Decrement
00538 
00539   index_iterator_base& operator++() { ++m_index; return *this; }
00540   index_iterator_base& operator--() { --m_index; return *this; }
00541 
00542   index_iterator_base operator++(int) { return index_iterator_base(m_src, m_index++); }
00543   index_iterator_base operator--(int) { return index_iterator_base(m_src, m_index--); }
00544 
00545   index_iterator_base& operator+=(int x) { m_index += x; return *this; }
00546   index_iterator_base& operator-=(int x) { m_index -= x; return *this; }
00547 
00548   index_iterator_base operator+(int x) const { return index_iterator_base(m_src, m_index+x); }
00549   index_iterator_base operator-(int x) const { return index_iterator_base(m_src, m_index-x); }
00550 };
00551 
00552 
00553 //  #######################################################
00554 //  =======================================================
00556 
00557 template <class T>
00558 class colmaj_iterator_base
00559 {
00560   int m_rowgap;
00561   int m_rowstride;
00562   T* m_ptr;
00563   T* m_current_end;
00564 
00565 public:
00566 
00567   typedef std::forward_iterator_tag iterator_category;
00568 
00569   typedef T            value_type;
00570   typedef ptrdiff_t    difference_type;
00571   typedef T*           pointer;
00572   typedef T&           reference;
00573 
00574   colmaj_iterator_base(int rg, int rs, T* ptr) :
00575     m_rowgap(rg),
00576     m_rowstride(rs),
00577     m_ptr(ptr),
00578     m_current_end(m_ptr+(m_rowstride-m_rowgap))
00579   {}
00580 
00581   T& operator*() const { return *m_ptr; }
00582 
00583   colmaj_iterator_base& operator++()
00584   {
00585     if (++m_ptr == m_current_end)
00586       {
00587         m_ptr += m_rowgap;
00588         m_current_end += m_rowstride;
00589       }
00590     return *this;
00591   }
00592 
00593   bool operator==(const colmaj_iterator_base& other) const
00594   { return m_ptr == other.m_ptr; }
00595 
00596   bool operator!=(const colmaj_iterator_base& other) const
00597   { return m_ptr != other.m_ptr; }
00598 };
00599 
00600 
00601 //  #######################################################
00602 //  =======================================================
00604 
00605 template <class T>
00606 class rowmaj_iterator_base
00607 {
00608   int m_stride;
00609   T* m_current_start;
00610   T* m_ptr;
00611   T* m_current_end;
00612 
00613 public:
00614 
00615   typedef std::forward_iterator_tag iterator_category;
00616 
00617   typedef T            value_type;
00618   typedef ptrdiff_t    difference_type;
00619   typedef T*           pointer;
00620   typedef T&           reference;
00621 
00622   rowmaj_iterator_base(int s, int ncols, T* ptr) :
00623     m_stride(s),
00624     m_current_start(ptr),
00625     m_ptr(ptr),
00626     m_current_end(m_ptr+(ncols*s))
00627   {}
00628 
00629   T& operator*() const { return *m_ptr; }
00630 
00631   rowmaj_iterator_base& operator++()
00632   {
00633     m_ptr += m_stride;
00634     if (m_ptr == m_current_end)
00635       {
00636         ++m_current_start;
00637         ++m_current_end;
00638         m_ptr = m_current_start;
00639       }
00640     return *this;
00641   }
00642 
00643   bool operator==(const rowmaj_iterator_base& other) const
00644   { return m_ptr == other.m_ptr; }
00645 
00646   bool operator!=(const rowmaj_iterator_base& other) const
00647   { return m_ptr != other.m_ptr; }
00648 };
00649 
00650 
00651 //  #######################################################
00652 //  =======================================================
00654 
00655 template <class Data>
00656 class mtx_base : public mtx_specs, public mtx_policies
00657 {
00658 private:
00659   mtx_base& operator=(const mtx_base& other); // not allowed
00660 
00661 protected:
00662   Data m_data;
00663 
00664   void swap(mtx_base& other);
00665 
00666   mtx_base(const mtx_base& other);
00667 
00668   mtx_base(int mrows, int ncols, const Data& data);
00669 
00670   mtx_base(const mtx_specs& specs, const Data& data);
00671 
00672   ~mtx_base();
00673 
00674   ptrdiff_t offset_from_storage(int r, int c) const
00675   { return RCR_leq(mtx_specs::offset_from_storage(r, c), m_data.storage_length()); }
00676 
00677   double* address_nc(int row, int col)
00678   { return m_data.storage_nc() + offset_from_storage(row, col); }
00679 
00680   const double* address(int row, int col) const
00681   { return m_data.storage() + offset_from_storage(row, col); }
00682 
00683   // Does the same thing as address_nc(), but doesn't range-check, since
00684   // this result is assumed to be some kind of "one-past-the-end" offset.
00685   ptrdiff_t end_offset_from_storage(int r, int c) const
00686   { return mtx_specs::offset_from_storage(r, c); }
00687 
00688   // Does the same thing as address_nc(), but doesn't range-check, since
00689   // this result is assumed to be some kind of "one-past-the-end" address.
00690   double* end_address_nc(int row, int col)
00691   { return m_data.storage_nc() + end_offset_from_storage(row, col); }
00692 
00693   // Does the same thing as address(), but doesn't range-check, since
00694   // this result is assumed to be some kind of "one-past-the-end" address.
00695   const double* end_address(int row, int col) const
00696   { return m_data.storage() + end_offset_from_storage(row, col); }
00697 
00698   const double* storage() const { return m_data.storage(); }
00699   double* storage_nc() { return m_data.storage_nc(); }
00700 
00701 public:
00702 
00703   const mtx_shape& shape() const { return specs().shape(); }
00704   const mtx_specs& specs() const { return *this; }
00705 
00706   const double& at(int i) const
00707   {
00708     RC_less(i+offset(), m_data.storage_length());
00709     return m_data.storage()[i+offset()];
00710   }
00711 
00712   double& at_nc(int i)
00713   {
00714     RC_less(i+offset(), m_data.storage_length());
00715     return m_data.storage_nc()[i+offset()];
00716   }
00717 
00718   template <class F>
00719   void apply(F func)
00720   {
00721     double* p = m_data.storage_nc() + offset();
00722     unsigned int gap = rowgap();
00723 
00724     if (gap == 0)
00725       {
00726         double* end = p + nelems();
00727         for (; p < end; ++p)
00728           *p = func(*p);
00729       }
00730     else
00731       {
00732         for (int c = 0; c < ncols(); ++c)
00733           {
00734             for (int r = 0; r < mrows(); ++r)
00735               {
00736                 *p = func(*p);
00737                 ++p;
00738               }
00739             p += gap;
00740           }
00741       }
00742   }
00743 
00744   //
00745   // Iterators
00746   //
00747 
00748   // Index-based iteration
00749 
00750   typedef index_iterator_base<mtx_base, double> iterator;
00751   typedef index_iterator_base<const mtx_base, const double> const_iterator;
00752 
00753   iterator begin_nc() { return iterator(this, 0); }
00754   iterator end_nc() { return iterator(this, nelems()); }
00755 
00756   const_iterator begin() const { return const_iterator(this, 0); }
00757   const_iterator end() const { return const_iterator(this, nelems()); }
00758 
00759   // Column-major iteration
00760 
00761   typedef colmaj_iterator_base<double> colmaj_iter;
00762   typedef colmaj_iterator_base<const double> const_colmaj_iter;
00763 
00764   colmaj_iter colmaj_begin_nc()
00765   { return colmaj_iter(rowgap(), rowstride(), address_nc(0,0)); }
00766 
00767   colmaj_iter colmaj_end_nc()
00768   { return colmaj_iter(rowgap(), rowstride(), end_address_nc(0,ncols())); }
00769 
00770   const_colmaj_iter colmaj_begin() const
00771   { return const_colmaj_iter(rowgap(), rowstride(), address(0,0)); }
00772 
00773   const_colmaj_iter colmaj_end() const
00774   { return const_colmaj_iter(rowgap(), rowstride(), end_address(0,ncols())); }
00775 
00776   // Row-major iteration
00777 
00778   typedef rowmaj_iterator_base<double> rowmaj_iter;
00779   typedef rowmaj_iterator_base<const double> const_rowmaj_iter;
00780 
00781   rowmaj_iter rowmaj_begin_nc()
00782   { return rowmaj_iter(rowstride(), ncols(), address_nc(0,0)); }
00783 
00784   rowmaj_iter rowmaj_end_nc()
00785   { return rowmaj_iter(rowstride(), ncols(), address_nc(mrows(),0)); }
00786 
00787   const_rowmaj_iter rowmaj_begin() const
00788   { return const_rowmaj_iter(rowstride(), ncols(), address(0,0)); }
00789 
00790   const_rowmaj_iter rowmaj_end() const
00791   { return const_rowmaj_iter(rowstride(), ncols(), address(mrows(),0)); }
00792 
00793 
00794   //
00795   // Functions
00796   //
00797 
00798   void clear(double x = 0.0) { apply(dash::setter(x)); }
00799 };
00800 
00801 
00802 //  #######################################################
00803 //  =======================================================
00805 
00806 class sub_mtx_ref : public mtx_base<data_ref_holder>
00807 {
00808 public:
00809   typedef mtx_base<data_ref_holder> Base;
00810 
00811   sub_mtx_ref(const mtx_specs& specs, data_holder& ref) :
00812     Base(specs, data_ref_holder(&ref))
00813   {}
00814 
00815   sub_mtx_ref& operator=(const sub_mtx_ref& other);
00816   sub_mtx_ref& operator=(const mtx& other);
00817 };
00818 
00819 
00820 //  #######################################################
00821 //  =======================================================
00823 
00824 class mtx : public mtx_base<data_holder>
00825 {
00826 public:
00827 
00828   typedef mtx_base<data_holder> Base;
00829 
00830   //
00831   // Constructors + Conversion
00832   //
00833 
00834   static const mtx& empty_mtx();
00835 
00836   mtx(const mtx_specs& specs, const data_holder& data) :
00837     Base(specs, data)
00838   {}
00839 
00841   static mtx colmaj_copy_of(const double* data, int mrows, int ncols);
00842 
00844   static mtx colmaj_borrow_from(double* data, int mrows, int ncols);
00845 
00847   static mtx colmaj_refer_to(double* data, int mrows, int ncols);
00848 
00849   static mtx zeros(const mtx_shape& s);
00850 
00851   static mtx uninitialized(const mtx_shape& s);
00852 
00853   static mtx zeros(int mrows, int ncols)
00854   { return zeros(mtx_shape(mrows, ncols)); }
00855 
00856   static mtx uninitialized(int mrows, int ncols)
00857   { return uninitialized(mtx_shape(mrows, ncols)); }
00858 
00859   static mtx from_stream(std::istream& s);
00860 
00861   static mtx from_string(const char* s);
00862 
00863   mtx(const slice& s);
00864 
00865   mtx(const mtx& other) : Base(other) {}
00866 
00867   virtual ~mtx();
00868 
00869   mtx& operator=(const mtx& other)
00870   {
00871     mtx temp(other);
00872     Base::swap(temp);
00873     return *this;
00874   }
00875 
00876   // This will destroy any data in the process of changing the size of
00877   // the mtx to the specified dimensions; its only advantage over just
00878   // declaring a new mtx is that it will avoid a deallocate/allocate
00879   // cycle if the new dimensions are the same as the current dimensions.
00880   void resize(int mrows_new, int ncols_new);
00881 
00886   mtx contig() const;
00887 
00888   data_holder& get_data_holder() { return m_data; }
00889 
00890 
00891   //
00892   // I/O
00893   //
00894 
00896   void print(std::ostream& s, const char* mtx_name = "") const;
00897 
00899   void print_stdout() const;
00900 
00902   void print_stdout_named(const char* mtx_name) const;
00903 
00905   rutz::fstring as_string() const;
00906 
00908   void scan(std::istream& s);
00909 
00911   void scan_string(const char* s);
00912 
00913   //
00914   // Iteration
00915   //
00916 
00917 
00918   //
00919   // Data access
00920   //
00921 
00922   double& at(int row, int col)
00923     { return Base::at_nc(offset_from_start(row, col)); }
00924 
00925   const double& at(int row, int col) const
00926     { return Base::at(offset_from_start(row, col)); }
00927 
00928   double& at(int elem)
00929     { return Base::at_nc(offset_from_start(elem)); }
00930 
00931   const double& at(int elem) const
00932     { return Base::at(offset_from_start(elem)); }
00933 
00934   bool same_size(const mtx& x) const
00935   { return (mrows() == x.mrows()) && (ncols() == x.ncols()); }
00936 
00937   //
00938   // Slices, submatrices
00939   //
00940 
00941   sub_mtx_ref sub(const row_index_range& rng)
00942   {
00943     return sub_mtx_ref(this->specs().sub_rows(rng), this->m_data);
00944   }
00945 
00946   sub_mtx_ref sub(const col_index_range& rng)
00947   {
00948     return sub_mtx_ref(this->specs().sub_cols(rng), this->m_data);
00949   }
00950 
00951   sub_mtx_ref sub(const row_index_range& rr, const col_index_range& cc)
00952   {
00953     return sub_mtx_ref(this->specs().sub_rows(rr).sub_cols(cc), this->m_data);
00954   }
00955 
00956 
00957   mtx sub(const row_index_range& rng) const
00958   {
00959     return mtx(this->specs().sub_rows(rng), this->m_data);
00960   }
00961 
00962   mtx sub(const col_index_range& rng) const
00963   {
00964     return mtx(this->specs().sub_cols(rng), this->m_data);
00965   }
00966 
00967   mtx sub(const row_index_range& rr, const col_index_range& cc) const
00968   {
00969     return mtx(this->specs().sub_rows(rr).sub_cols(cc), this->m_data);
00970   }
00971 
00972 
00973   mtx operator()(const row_index_range& rng) const { return sub(rng); }
00974 
00975   mtx operator()(const col_index_range& rng) const { return sub(rng); }
00976 
00977   mtx operator()(const row_index_range& rr, const col_index_range& cc) const
00978   { return sub(rr, cc); }
00979 
00980 
00981   mtx as_shape(const mtx_shape& s) const
00982   { return mtx(this->specs().as_shape(s), this->m_data); }
00983 
00984   mtx as_shape(int mr, int nc) const
00985   { return mtx(this->specs().as_shape(mr, nc), this->m_data); }
00986 
00987   slice row(int r) const
00988     { return slice(this->m_data, offset_from_storage(r,0), rowstride(), ncols()); }
00989 
00990   mtx_iter row_iter(int r)
00991     { return mtx_iter(address_nc(r,0), rowstride(), ncols()); }
00992 
00993   mtx_const_iter row_iter(int r) const
00994     { return mtx_const_iter(address(r,0), rowstride(), ncols()); }
00995 
00996   mtx as_row() const { return as_shape(1, nelems()); }
00997 
00998   void reorder_rows(const mtx& index);
00999 
01000 
01001 
01002   slice column(int c) const
01003     { return slice(this->m_data, offset_from_storage(0,c), colstride(), mrows()); }
01004 
01005   mtx_iter column_iter(int c)
01006     { return mtx_iter(address_nc(0,c), colstride(), mrows()); }
01007 
01008   mtx_const_iter column_iter(int c) const
01009     { return mtx_const_iter(address(0,c), colstride(), mrows()); }
01010 
01011   mtx as_column() const { return as_shape(nelems(), 1); }
01012 
01013   void reorder_columns(const mtx& index);
01014 
01015   void swap_columns(int c1, int c2);
01016 
01017 
01018   //
01019   // Functions
01020   //
01021 
01022   mtx mean_row() const;
01023 
01024   mtx mean_column() const;
01025 
01026   const_iterator find_min() const;
01027   const_iterator find_max() const;
01028 
01029   double min() const;
01030   double max() const;
01031 
01032   double sum() const;
01033 
01034   double mean() const { return sum() / nelems(); }
01035 
01036   mtx& operator+=(double x) { apply(dash::add(x)); return *this; }
01037   mtx& operator-=(double x) { apply(dash::sub(x)); return *this; }
01038   mtx& operator*=(double fac) { apply(dash::mul(fac)); return *this; }
01039   mtx& operator/=(double div) { apply(dash::div(div)); return *this; }
01040 
01041   mtx& operator+=(const mtx& other);
01042   mtx& operator-=(const mtx& other);
01043 
01044   bool operator==(const mtx& other) const;
01045 
01046   bool operator!=(const mtx& other) const
01047   { return !(operator==(other)); }
01048 
01049   // result = vec * mtx;
01050   static void VMmul_assign(const slice& vec, const mtx& mtx,
01051                            slice& result);
01052 
01053   // this = m1 * m2;
01054   void assign_MMmul(const mtx& m1, const mtx& m2);
01055 
01056   void make_unique() { Base::m_data.make_unique(); }
01057 
01058 private:
01059   friend class slice;
01060 };
01061 
01062 
01063 
01064 //  #######################################################
01065 //  #######################################################
01066 //  #######################################################
01067 //  =======================================================
01068 //  Inline free functions
01069 
01070 
01071 inline double inner_product(mtx_const_iter s1, mtx_const_iter s2)
01072 {
01073   double result = 0.0;
01074 
01075   for (; s1.has_more(); ++s1, ++s2)
01076     result += (*s1) * (*s2);
01077 
01078   return result;
01079 }
01080 
01081 mtx operator+(const mtx& m, double x);
01082 mtx operator-(const mtx& m, double x);
01083 mtx operator*(const mtx& m, double x);
01084 mtx operator/(const mtx& m, double x);
01085 
01086 mtx operator+(const mtx& m1, const mtx& m2);
01087 mtx operator-(const mtx& m1, const mtx& m2);
01088 
01089 // Simple element-by-element multiplication (i.e. NOT matrix multiplication)
01090 mtx arr_mul(const mtx& m1, const mtx& m2);
01091 
01092 // Simple element-by-element division (i.e. NOT matrix division)
01093 mtx arr_div(const mtx& m1, const mtx& m2);
01094 
01095 mtx min(const mtx& m1, const mtx& m2);
01096 mtx max(const mtx& m1, const mtx& m2);
01097 
01098 static const char __attribute__((used)) vcid_groovx_pkgs_mtx_mtx_h_utc20050626084022[] = "$Id: mtx.h 10065 2007-04-12 05:54:56Z rjpeters $ $HeadURL: file:
01099 #endif // !GROOVX_PKGS_MTX_MTX_H_UTC20050626084022_DEFINED

The software described here is Copyright (c) 1998-2005, Rob Peters.
This page was generated Wed Dec 3 06:49:39 2008 by Doxygen version 1.5.5.