TooN 2.1
helpers.h
00001 // -*- c++ -*-
00002 
00003 // Copyright (C) 2005,2009 Tom Drummond (twd20@cam.ac.uk),
00004 // Ed Rosten (er258@cam.ac.uk), Gerhard Reitmayr (gr281@cam.ac.uk)
00005 //
00006 // This file is part of the TooN Library.  This library is free
00007 // software; you can redistribute it and/or modify it under the
00008 // terms of the GNU General Public License as published by the
00009 // Free Software Foundation; either version 2, or (at your option)
00010 // any later version.
00011 
00012 // This library is distributed in the hope that it will be useful,
00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 // GNU General Public License for more details.
00016 
00017 // You should have received a copy of the GNU General Public License along
00018 // with this library; see the file COPYING.  If not, write to the Free
00019 // Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
00020 // USA.
00021 
00022 // As a special exception, you may use this file as part of a free software
00023 // library without restriction.  Specifically, if other files instantiate
00024 // templates or use macros or inline functions from this file, or you compile
00025 // this file and link it with other files to produce an executable, this
00026 // file does not by itself cause the resulting executable to be covered by
00027 // the GNU General Public License.  This exception does not however
00028 // invalidate any other reasons why the executable file might be covered by
00029 // the GNU General Public License.
00030 
00031 
00032 #ifndef TOON_INCLUDE_HELPERS_H
00033 #define TOON_INCLUDE_HELPERS_H
00034 
00035 #include <TooN/TooN.h>
00036 #include <TooN/gaussian_elimination.h>
00037 #include <cmath>
00038 #include <functional>
00039 #include <utility>
00040 
00041 #ifndef M_PI
00042 #define M_PI 3.14159265358979323846
00043 #endif
00044 
00045 #ifndef M_SQRT1_2
00046 #define M_SQRT1_2 0.707106781186547524401
00047 #endif
00048 
00049 #ifdef WIN32
00050 namespace std {
00051     inline int isfinite( const double & v ){ return _finite(v); }
00052     inline int isfinite( const float & v ){ return _finite(v); }
00053     inline int isnan( const double & v ){ return _isnan(v); }
00054     inline int isnan( const float & v ){ return _isnan(v); }
00055 }
00056 #endif
00057 
00058 namespace TooN {
00059     
00060     ///Invert a matrix.
00061     ///
00062     ///For sizes other than 2x2, @link gDecomps decompositions @endlink provide a suitable solition.
00063     ///
00064     ///@param m Matrix to invert.
00065     ///@ingroup gDecomps
00066     inline Matrix<2> inv(const Matrix<2>& m)
00067     {
00068         double d = 1./(m[0][0]*m[1][1] - m[1][0]*m[0][1]);
00069 
00070         return Data(
00071              m[1][1]*d, -m[0][1]*d,
00072             -m[1][0]*d,  m[0][0]*d
00073         );
00074     }
00075 
00076 
00077 
00078     ///\deprecated
00079     ///@ingroup gLinAlg
00080     template<int Size, class Precision, class Base> TOON_DEPRECATED void Fill(Vector<Size, Precision, Base>& v, const Precision& p)
00081     {
00082         for(int i=0; i < v.size(); i++)
00083             v[i]= p;
00084     }
00085 
00086     ///\deprecated
00087     ///@ingroup gLinAlg
00088     template<int Rows, int Cols, class Precision, class Base> TOON_DEPRECATED void Fill(Matrix<Rows, Cols, Precision, Base>& m, const Precision& p)
00089     {
00090         for(int i=0; i < m.num_rows(); i++)
00091             for(int j=0; j < m.num_cols(); j++)
00092                 m[i][j] = p;
00093     }
00094 
00095     ///Compute the \f$L_2\f$ norm of \e v
00096     ///@param v \e v
00097     ///@ingroup gLinAlg
00098     template<int Size, class Precision, class Base> inline Precision norm(const Vector<Size, Precision, Base>& v)
00099     {
00100         using std::sqrt;
00101         return sqrt(v*v);
00102     }
00103 
00104     ///Compute the \f$L_2^2\f$ norm of \e v
00105     ///@param v \e v
00106     ///@ingroup gLinAlg
00107     template<int Size, class Precision, class Base> inline Precision norm_sq(const Vector<Size, Precision, Base>& v)
00108     {
00109         return v*v;
00110     }
00111     
00112     ///Compute the \f$L_1\f$ norm of \e v
00113     ///@param v \e v
00114     ///@ingroup gLinAlg
00115     template<int Size, class Precision, class Base> inline Precision norm_1(const Vector<Size, Precision, Base>& v)
00116     {
00117         using std::abs;
00118         Precision n = 0;
00119         for(int i=0; i < v.size(); i++)
00120             n += abs(v[i]);
00121         return n;
00122     }
00123 
00124     ///Compute the \f$L_\infty\f$ norm of \e v
00125     ///@param v \e v
00126     ///@ingroup gLinAlg
00127     template<int Size, class Precision, class Base> inline Precision norm_inf(const Vector<Size, Precision, Base>& v)
00128     {
00129         using std::abs;
00130         using std::max;
00131         Precision n = 0;
00132         n = abs(v[0]);
00133 
00134         for(int i=1; i < v.size(); i++)
00135             n = max(n, abs(v[i]));
00136         return n;
00137     }
00138     
00139     ///Compute the \f$L_2\f$ norm of \e v.
00140     ///Synonym for norm()
00141     ///@param v \e v
00142     ///@ingroup gLinAlg
00143     template<int Size, class Precision, class Base> inline Precision norm_2(const Vector<Size, Precision, Base>& v)
00144     {
00145         return norm(v);
00146     }
00147     
00148 
00149 
00150 
00151     ///Compute a the unit vector \f$\hat{v}\f$.
00152     ///@param v \e v
00153     ///@ingroup gLinAlg
00154     template<int Size, class Precision, class Base> inline Vector<Size, Precision> unit(const Vector<Size, Precision, Base> & v)
00155     {
00156         using std::sqrt;
00157         return TooN::operator*(v,(1/sqrt(v*v)));
00158     }
00159     
00160     //Note because of the overload later, this function will ONLY receive sliced vectors. Therefore
00161     //a copy can be made, which is still a slice, so operating on the copy operates on the original
00162     //data.
00163     ///Normalize a vector in place
00164     ///@param v Vector to normalize
00165     ///@ingroup gLinAlg
00166     template<int Size, class Precision, class Base> inline void normalize(Vector<Size, Precision, Base> v)
00167     {
00168         using std::sqrt;
00169         v /= sqrt(v*v);
00170     }
00171     
00172     //This overload is required to operate on non-slice vectors
00173     /**
00174         \overload
00175     */  
00176     template<int Size, class Precision> inline void normalize(Vector<Size, Precision> & v)
00177     {
00178         normalize(v.as_slice());
00179     }
00180 
00181     ///For a vector \e v of length \e i, return \f$[v_1, v_2, \cdots, v_{i-1}] / v_i \f$
00182     ///@param v \e v
00183     ///@ingroup gLinAlg
00184     // Don't remove the +0 in the first template parameter of the return type. It is a work around for a Visual Studio 2010 bug:
00185     // https://connect.microsoft.com/VisualStudio/feedback/details/735283/vc-2010-parse-error-in-template-parameters-using-ternary-operator
00186     template<int Size, typename Precision, typename Base> inline Vector<(Size==Dynamic?Dynamic:Size-1)+0, Precision> project( const Vector<Size, Precision, Base> & v){
00187         static const int Len = (Size==Dynamic?Dynamic:Size-1);
00188         return TooN::operator/(v.template slice<0, Len>(0 , v.size() - 1) , v[v.size() - 1]);
00189     }
00190 
00191     //This should probably be done with an operator to prevent an extra new[] for dynamic vectors.
00192     ///For a vector \e v of length \e i, return \f$[v_1, v_2, \cdots, v_{i}, 1]\f$
00193     ///@param v \e v
00194     ///@ingroup gLinAlg
00195     // Don't remove the +0 in the first template parameter of the return type. It is a work around for a Visual Studio 2010 bug:
00196     // https://connect.microsoft.com/VisualStudio/feedback/details/735283/vc-2010-parse-error-in-template-parameters-using-ternary-operator
00197     template<int Size, typename Precision, typename Base> inline Vector<(Size==Dynamic?Dynamic:Size+1)+0, Precision> unproject( const Vector<Size, Precision, Base> & v){
00198         Vector<(Size==Dynamic?Dynamic:Size+1), Precision> result(v.size()+1);
00199         static const int Len = (Size==Dynamic?Dynamic:Size);
00200         result.template slice<0, Len>(0, v.size()) = v;
00201         result[v.size()] = 1;
00202         return result;
00203     }
00204     
00205     /**
00206        \overload
00207     */
00208     template<int R, int C, typename Precision, typename Base> inline Matrix<R-1, C, Precision> project( const Matrix<R,C, Precision, Base> & m){
00209         Matrix<R-1, C, Precision> result = m.slice(0,0,R-1,m.num_cols());
00210         for( int c = 0; c < m.num_cols(); ++c ) {
00211             result.slice(0,c,R-1,1) /= m[R-1][c];
00212         }
00213         return result;
00214     }
00215 
00216     template<int C, typename Precision, typename Base> inline Matrix<-1, C, Precision> project( const Matrix<-1,C, Precision, Base> & m){
00217         Matrix<-1, C, Precision> result = m.slice(0,0,m.num_rows()-1,m.num_cols());
00218         for( int c = 0; c < m.num_cols(); ++c ) {
00219             result.slice(0,c,m.num_rows()-1,1) /= m[m.num_rows()-1][c];
00220         }
00221         return result;
00222     }
00223 
00224     /**
00225        \overload
00226     */
00227     template<int R, int C, typename Precision, typename Base> inline Matrix<R+1, C, Precision> unproject( const Matrix<R, C, Precision, Base> & m){
00228         Matrix<R+1, C, Precision> result;
00229         result.template slice<0,0,R,C>() = m;
00230         result[R] = Ones;
00231         return result;
00232     }
00233 
00234     template<int C, typename Precision, typename Base> inline Matrix<-1, C, Precision> unproject( const Matrix<-1, C, Precision, Base> & m){
00235         Matrix<-1, C, Precision> result( m.num_rows()+1, m.num_cols() );
00236         result.slice(0,0,m.num_rows(),m.num_cols()) = m;
00237         result[m.num_rows()] = Ones;
00238         return result;
00239     }
00240 
00241     /// Frobenius (root of sum of squares) norm of input matrix \e m
00242     ///@param m \e m
00243     ///@ingroup gLinAlg
00244     template <int R, int C, typename P, typename B>
00245     P inline norm_fro( const Matrix<R,C,P,B> & m ){
00246         using std::sqrt;
00247         P n = 0;
00248         for(int r = 0; r < m.num_rows(); ++r)
00249             for(int c = 0; c < m.num_cols(); ++c)
00250                 n += m[r][c] * m[r][c];
00251 
00252         return sqrt(n);
00253     }
00254 
00255     /// \e L<sub>&#8734;</sub> (row sum) norm of input matrix m
00256     /// computes the maximum of the sums of absolute values over rows
00257     ///@ingroup gLinAlg
00258     template <int R, int C, typename P, typename B>
00259     P inline norm_inf( const Matrix<R,C,P,B> & m ){
00260         using std::abs;
00261         using std::max;
00262         P n = 0;
00263         for(int r = 0; r < m.num_rows(); ++r){
00264             P s = 0;
00265             for(int c = 0; c < m.num_cols(); ++c)
00266                 s += abs(m(r,c));
00267             n = max(n,s);
00268         }
00269         return n;
00270     }
00271     
00272     /// \e L<sub>1</sub> (col sum) norm of input matrix m
00273     /// computes the maximum of the sums of absolute values over columns
00274     ///@ingroup gLinAlg
00275     template <int R, int C, typename P, typename B>
00276     P inline norm_1( const Matrix<R,C,P,B> & m ){
00277         using std::abs;
00278         using std::max;
00279         P n = 0;
00280         for(int c = 0; c < m.num_cols(); ++c){
00281             P s = 0;
00282             for(int r = 0; r < m.num_rows(); ++r)
00283                 s += abs(m(r,c));
00284             n = max(n,s);
00285         }
00286         return n;
00287     }
00288 
00289     namespace Internal {
00290         ///@internal
00291         ///@brief Exponentiate a matrix using a the Taylor series
00292         ///This will not work if the norm of the matrix is too large.
00293         template <int R, int C, typename P, typename B>
00294         inline Matrix<R, C, P> exp_taylor( const Matrix<R,C,P,B> & m ){
00295             TooN::SizeMismatch<R, C>::test(m.num_rows(), m.num_cols());
00296             Matrix<R,C,P> result = TooN::Zeros(m.num_rows(), m.num_cols());
00297             Matrix<R,C,P> f = TooN::Identity(m.num_rows());
00298             P k = 1;
00299             while(norm_inf((result+f)-result) > 0){
00300                 result += f;
00301                 f = (m * f) / k;
00302                 k += 1;
00303             }
00304             return result;
00305         }
00306 
00307         ///@internal
00308         ///@brief Logarithm of a matrix using a the Taylor series
00309         ///This will not work if the norm of the matrix is too large.
00310         template <int R, int C, typename P, typename B>
00311         inline Matrix<R, C, P> log_taylor( const Matrix<R,C,P,B> & m ){
00312             TooN::SizeMismatch<R, C>::test(m.num_rows(), m.num_cols());
00313             Matrix<R,C,P> X = m - TooN::Identity * 1.0;
00314             Matrix<R,C,P> F = X;
00315             Matrix<R,C,P> sum = TooN::Zeros(m.num_rows(), m.num_cols());
00316             P k = 1;
00317             while(norm_inf((sum+F/k)-sum) > 0){
00318                 sum += F/k;
00319                 F = -F*X;
00320                 k += 1;
00321             }
00322             return sum;
00323         }
00324 
00325     };
00326     
00327     /// computes the matrix exponential of a matrix m by 
00328     /// scaling m by 1/(powers of 2), using Taylor series and 
00329     /// squaring again.
00330     /// @param m input matrix, must be square
00331     /// @return result matrix of the same size/type as input
00332     /// @ingroup gLinAlg
00333     template <int R, int C, typename P, typename B>
00334     inline Matrix<R, C, P> exp( const Matrix<R,C,P,B> & m ){
00335         using std::max;
00336         SizeMismatch<R, C>::test(m.num_rows(), m.num_cols());
00337         const P l = log2(norm_inf(m));
00338         const int s = max(0,(int)ceil(l));
00339         Matrix<R,C,P> result = Internal::exp_taylor(m/(1<<s));
00340         for(int i = 0; i < s; ++i)
00341             result = result * result;
00342         return result;
00343     }
00344     
00345     /// computes a matrix square root of a matrix m by
00346     /// the product form of the Denman and Beavers iteration
00347     /// as given in Chen et al. 'Approximating the logarithm of a matrix to specified accuracy', 
00348     /// J. Matrix Anal Appl, 2001. This is used for the matrix
00349     /// logarithm function, but is useable by on its own.
00350     /// @param m input matrix, must be square
00351     /// @return a square root of m of the same size/type as input
00352     /// @ingroup gLinAlg
00353     template <int R, int C, typename P, typename B>
00354     inline Matrix<R, C, P> sqrt( const Matrix<R,C,P,B> & m){
00355         SizeMismatch<R, C>::test(m.num_rows(), m.num_cols());
00356         Matrix<R,C,P> M = m;
00357         Matrix<R,C,P> Y = m;
00358         Matrix<R,C,P> M_inv(m.num_rows(), m.num_cols());
00359         const Matrix<R,C,P> id = Identity(m.num_rows());
00360         do {
00361             M_inv = gaussian_elimination(M, id);
00362             Y = Y * (id + M_inv) * 0.5;
00363             M = 0.5 * (id + (M + M_inv) * 0.5);
00364         } while(norm_inf(M - M_inv) > 0);
00365         return Y;
00366     }
00367     
00368     /// computes the matrix logarithm of a matrix m using the inverse scaling and 
00369     /// squaring method. The overall approach is described in
00370     /// Chen et al. 'Approximating the logarithm of a matrix to specified accuracy', 
00371     /// J. Matrix Anal Appl, 2001, but this implementation only uses a simple
00372     /// taylor series after the repeated square root operation.
00373     /// @param m input matrix, must be square
00374     /// @return the log of m of the same size/type as input
00375     /// @ingroup gLinAlg
00376     template <int R, int C, typename P, typename B>
00377     inline Matrix<R, C, P> log( const Matrix<R,C,P,B> & m){
00378         int counter = 0;
00379         Matrix<R,C,P> A = m;
00380         while(norm_inf(A - Identity*1.0) > 0.5){
00381             ++counter;
00382             A = sqrt(A);
00383         }
00384         return Internal::log_taylor(A) * pow(2.0, counter);
00385     }
00386     
00387     /// Returns true if every element is finite
00388     ///@ingroup gLinAlg
00389     template<int S, class P, class B> bool isfinite(const Vector<S, P, B>& v)
00390     { 
00391         using std::isfinite;
00392         for(int i=0; i < v.size(); i++)
00393             if(!isfinite(v[i]))
00394                 return 0;
00395         return 1;
00396     }
00397 
00398     /// Returns true if any element is NaN
00399     ///@ingroup gLinAlg
00400     template<int S, class P, class B> bool isnan(const Vector<S, P, B>& v)
00401     { 
00402         using std::isnan;
00403         for(int i=0; i < v.size(); i++)
00404             if(isnan(v[i]))
00405                 return 1;
00406         return 0;
00407     }
00408 
00409     /// Symmetrize a matrix 
00410     ///@param m \e m
00411     ///@return \f$ \frac{m + m^{\mathsf T}}{2} \f$
00412     ///@ingroup gLinAlg
00413     template<int Rows, int Cols, typename Precision, typename Base>
00414     void Symmetrize(Matrix<Rows,Cols,Precision,Base>& m){
00415         SizeMismatch<Rows,Cols>::test(m.num_rows(), m.num_cols());
00416         for(int r=0; r<m.num_rows()-1; r++){
00417             for(int c=r+1; c<m.num_cols(); c++){
00418                 const Precision temp=(m(r,c)+m(c,r))/2;
00419                 m(r,c)=temp;
00420                 m(c,r)=temp;
00421             }
00422         }
00423     }
00424     
00425     /// computes the trace of a square matrix
00426     ///@ingroup gLinAlg
00427     template<int Rows, int Cols, typename Precision, typename Base>
00428     Precision trace(const Matrix<Rows, Cols, Precision, Base> & m ){
00429         SizeMismatch<Rows,Cols>::test(m.num_rows(), m.num_cols());
00430         Precision tr = 0;
00431         for(int i = 0; i < m.num_rows(); ++i)
00432             tr += m(i,i);
00433         return tr;
00434     }
00435 
00436     /// creates an returns a cross product matrix M from a 3 vector v, such that for all vectors w, the following holds: v ^ w = M * w
00437     /// @param vec the 3 vector input
00438     /// @return the 3x3 matrix to set to the cross product matrix
00439     ///@ingroup gLinAlg
00440     template<int Size, class P, class B> inline TooN::Matrix<3, 3, P> cross_product_matrix(const Vector<Size, P, B>& vec)
00441     {
00442         SizeMismatch<Size,3>::test(vec.size(), 3);
00443 
00444         TooN::Matrix<3, 3, P> result;
00445 
00446         result(0,0) = 0; 
00447         result(0,1) = -vec[2]; 
00448         result(0,2) = vec[1];
00449         result(1,0) = vec[2]; 
00450         result(1,1) = 0; 
00451         result(1,2) = -vec[0];
00452         result(2,0) = -vec[1]; 
00453         result(2,1) = vec[0]; 
00454         result(2,2) = 0;
00455 
00456         return result;
00457     }
00458 
00459     namespace Internal {
00460         template<int Size, typename Precision, typename Base, typename Func, typename Ret> inline Ret accumulate( const Vector<Size, Precision, Base> & v )  {
00461             Func func;
00462             if( v.size() == 0 ) {
00463                 return func.null(); // What should we return, exception?
00464             }
00465             func.initialise( v[0], 0 );
00466             for( int ii = 1; ii < v.size(); ii++ ) {
00467                 func( v[ii], ii );
00468             }
00469             return func.ret();
00470         }
00471 
00472         template<int R, int C, typename Precision, typename Base, typename Func, typename Ret> inline Ret accumulate( const Matrix<R, C, Precision, Base> & m )  {
00473             Func func;
00474             if( m.num_rows() == 0 || m.num_cols() == 0) {
00475                 return func.null(); // What should we return, exception?
00476             }
00477             func.initialise( m[0][0], 0, 0 );
00478             for(int r=0; r<m.num_rows(); r++){
00479                 for(int c=0; c<m.num_cols(); c++){
00480                     func( m[r][c], r, c );
00481                 }
00482             }
00483             return func.ret();
00484         }
00485         template<int R, int C, typename Precision, typename Base, typename Func, typename Ret> inline Ret accumulate_horizontal( const Matrix<R, C, Precision, Base> & m ) {
00486             Func func( m.num_rows() );
00487             if( m.num_cols() == 0 || m.num_rows() == 0 ) {
00488                 func.null(); // What should we return, exception?
00489             }
00490             for(int r=0; r<m.num_rows(); r++){
00491                 func.initialise( m[r][0], r, 0 );
00492                 for(int c=1; c<m.num_cols(); c++){
00493                     func( m[r][c], r, c );
00494                 }
00495             }
00496             return func.ret();
00497         }
00498         template<int R, int C, typename Precision, typename Base, typename Func, typename Ret> inline Ret accumulate_vertical( const Matrix<R, C, Precision, Base> & m ) {
00499             Func func( m.num_cols() );
00500             if( m.num_cols() == 0 || m.num_rows() == 0 ) {
00501                 func.null(); // What should we return, exception?
00502             }
00503             for(int c=0; c<m.num_cols(); c++){
00504                 func.initialise( m[0][c], 0, c );
00505                 for(int r=1; r<m.num_rows(); r++){
00506                     func( m[r][c], r, c );
00507                 }
00508             }
00509             return func.ret();
00510         }        
00511 
00512         template<typename Precision, typename ComparisonFunctor>
00513         class accumulate_functor_vector {
00514             Precision bestVal;
00515         public:
00516             Precision null() {
00517                 return 0;
00518             }
00519             void initialise( Precision initialVal, int ) {
00520                 bestVal = initialVal;
00521             }
00522             void operator()( Precision curVal, int ) {
00523                 if( ComparisonFunctor()( curVal, bestVal ) ) {
00524                     bestVal = curVal;
00525                 }
00526             }
00527             Precision ret() { return bestVal; }            
00528         };
00529         template<typename Precision, typename ComparisonFunctor>
00530         class accumulate_element_functor_vector {
00531             Precision bestVal;
00532             int nBestIndex;
00533         public:
00534             std::pair<Precision,int> null() {
00535                 return std::pair<Precision,int>( 0, 0 );
00536             }
00537             void initialise( Precision initialVal, int nIndex ) {
00538                 bestVal = initialVal;
00539                 nBestIndex = nIndex;
00540             }
00541             void operator()( Precision curVal, int nIndex ) {
00542                 if( ComparisonFunctor()( curVal, bestVal ) ) {
00543                     bestVal = curVal;
00544                     nBestIndex = nIndex;
00545                 }
00546             }
00547             std::pair<Precision,int> ret() {
00548                 return std::pair<Precision,int>( bestVal, nBestIndex );
00549             }            
00550         };
00551         template<typename Precision, typename ComparisonFunctor>
00552         class accumulate_functor_matrix {
00553             Precision bestVal;
00554         public:
00555             Precision null() {
00556                 return 0;
00557             }
00558             void initialise( Precision initialVal, int, int ) {
00559                 bestVal = initialVal;
00560             }
00561             void operator()( Precision curVal, int, int ) {
00562                 if( ComparisonFunctor()( curVal, bestVal ) ) {
00563                     bestVal = curVal;
00564                 }
00565             }
00566             Precision ret() { return bestVal; }            
00567         };
00568         template<typename Precision, typename ComparisonFunctor>
00569         class accumulate_element_functor_matrix {
00570             Precision bestVal;
00571             int nBestRow;
00572             int nBestCol;
00573         public:
00574             std::pair<Precision,std::pair<int,int> > null() {
00575                 return std::pair<Precision,std::pair<int,int> >( 0, std::pair<int,int>( 0, 0 ) );
00576             }
00577             void initialise( Precision initialVal, int nRow, int nCol ) {
00578                 bestVal = initialVal;
00579                 nBestRow = nRow;
00580                 nBestCol = nCol;
00581             }
00582             void operator()( Precision curVal, int nRow, int nCol ) {
00583                 if( ComparisonFunctor()( curVal, bestVal ) ) {
00584                     bestVal = curVal;
00585                     nBestRow = nRow;
00586                     nBestCol = nCol;
00587                 }
00588             }
00589             std::pair<Precision,std::pair<int,int> > ret() {
00590                 return std::pair<Precision,std::pair<int,int> >( bestVal, 
00591                                                                  std::pair<int,int>( nBestRow, nBestCol ) );
00592             }            
00593         };
00594         template<typename Precision, typename ComparisonFunctor>
00595         class accumulate_vertical_functor {
00596             Vector<Dynamic,Precision>* bestVal;
00597         public:
00598             accumulate_vertical_functor() {
00599                 bestVal = NULL;
00600             }
00601             accumulate_vertical_functor( int nNumCols ) {
00602                 bestVal = new Vector<Dynamic,Precision>( nNumCols );
00603             }
00604             Vector<Dynamic,Precision> null() {
00605                 return Vector<Dynamic,Precision>( 0 );
00606             }
00607             void initialise( Precision initialVal, int, int nCol ) {
00608                 (*bestVal)[nCol] = initialVal;
00609             }
00610             void operator()( Precision curVal, int, int nCol ) {
00611                 if( ComparisonFunctor()( curVal, (*bestVal)[nCol] ) ) {
00612                     (*bestVal)[nCol] = curVal;
00613                 }
00614             }
00615             Vector<Dynamic,Precision> ret() {
00616                 if( bestVal == NULL ) {
00617                     return null();
00618                 }
00619                 Vector<Dynamic,Precision> vRet = *bestVal; 
00620                 delete bestVal;
00621                 return vRet;
00622             }            
00623         };
00624         template<typename Precision, typename ComparisonFunctor>
00625         class accumulate_element_vertical_functor {
00626             Vector<Dynamic,Precision>* bestVal;
00627             Vector<Dynamic,Precision>* bestIndices;
00628         public:
00629             accumulate_element_vertical_functor() {
00630                 bestVal = NULL;
00631                 bestIndices = NULL;
00632             }
00633             accumulate_element_vertical_functor( int nNumCols ) {
00634                 bestVal = new Vector<Dynamic,Precision>( nNumCols );
00635                 bestIndices = new Vector<Dynamic,Precision>( nNumCols );
00636             }
00637             std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > null() {
00638                 Vector<Dynamic,Precision> vEmpty( 0 );
00639                 return std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> >( vEmpty, vEmpty );
00640             }
00641             void initialise( Precision initialVal, int nRow, int nCol ) {
00642                 (*bestVal)[nCol] = initialVal;
00643                 (*bestIndices)[nCol] = nRow;
00644             }
00645             void operator()( Precision curVal, int nRow, int nCol ) {
00646                 if( ComparisonFunctor()( curVal, (*bestVal)[nCol] ) ) {
00647                     (*bestVal)[nCol] = curVal;
00648                     (*bestIndices)[nCol] = nRow;
00649                 }
00650             }
00651             std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > ret() {
00652                 if( bestVal == NULL ) {
00653                     return null();
00654                 }
00655                 std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > vRet = 
00656                     std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > (*bestVal, *bestIndices );
00657                 delete bestVal; bestVal = NULL;
00658                 delete bestIndices; bestIndices = NULL;
00659                 return vRet;
00660             }            
00661         };
00662         template<typename Precision, typename ComparisonFunctor>
00663         class accumulate_horizontal_functor {
00664             Vector<Dynamic,Precision>* bestVal;
00665         public: 
00666             accumulate_horizontal_functor() {
00667                 bestVal = NULL;
00668             }
00669             accumulate_horizontal_functor( int nNumRows ) {
00670                 bestVal = new Vector<Dynamic,Precision>( nNumRows );
00671             }
00672             Vector<Dynamic,Precision> null() {
00673                 return Vector<Dynamic,Precision>( 0 );
00674             }
00675             void initialise( Precision initialVal, int nRow, int ) {
00676                 (*bestVal)[nRow] = initialVal;
00677             }
00678             void operator()( Precision curVal, int nRow, int ) {
00679                 if( ComparisonFunctor()( curVal, (*bestVal)[nRow] ) ) {
00680                     (*bestVal)[nRow] = curVal;
00681                 }
00682             }
00683             Vector<Dynamic,Precision> ret() { 
00684                 if( bestVal == NULL ) {
00685                     return null();
00686                 }
00687                 Vector<Dynamic,Precision> vRet = *bestVal;
00688                 delete bestVal; bestVal = NULL;
00689                 return vRet; 
00690             }            
00691         };
00692         template<typename Precision, typename ComparisonFunctor>
00693         class accumulate_element_horizontal_functor {
00694             Vector<Dynamic,Precision>* bestVal;
00695             Vector<Dynamic,Precision>* bestIndices;
00696         public:
00697             accumulate_element_horizontal_functor() {
00698                 bestVal = NULL;
00699                 bestIndices = NULL;
00700             }
00701             accumulate_element_horizontal_functor( int nNumRows ) {
00702                 bestVal = new Vector<Dynamic,Precision>( nNumRows );
00703                 bestIndices = new Vector<Dynamic,Precision>( nNumRows );
00704             }
00705             std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > null() {
00706                 Vector<Dynamic,Precision> vEmpty( 0 );
00707                 return std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> >( vEmpty, vEmpty );
00708             }
00709             void initialise( Precision initialVal, int nRow, int nCol ) {
00710                 (*bestVal)[nRow] = initialVal;
00711                 (*bestIndices)[nRow] = nCol;
00712             }
00713             void operator()( Precision curVal, int nRow, int nCol ) {
00714                 if( ComparisonFunctor()( curVal, (*bestVal)[nRow] ) ) {
00715                     (*bestVal)[nRow] = curVal;
00716                     (*bestIndices)[nRow] = nCol;
00717                 }
00718             }
00719             std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > ret() {
00720                 if( bestVal == NULL ) {
00721                     return null();
00722                 }
00723                 std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > vRet = 
00724                     std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> >( *bestVal, *bestIndices );
00725                 delete bestVal; bestVal = NULL;
00726                 delete bestIndices; bestIndices = NULL;
00727                 return vRet;
00728             }            
00729         };
00730     }
00731 
00732 
00733     /// Finds the minimal value of a vector.
00734     /// @param v a vector
00735     /// @return the smallest value of v
00736     template<int Size, typename Precision, typename Base> inline Precision min_value( const Vector<Size, Precision, Base>& v) {
00737         typedef Internal::accumulate_functor_vector<Precision, std::less<Precision> > vector_accumulate_functor;
00738         return Internal::accumulate<Size,Precision,Base,
00739             vector_accumulate_functor, Precision >( v ); 
00740     }
00741     /// Finds the largest value of a vector.
00742     /// @param v a vector
00743     /// @return the largest value of v
00744     template<int Size, typename Precision, typename Base> inline Precision max_value( const Vector<Size, Precision, Base>& v) {
00745         typedef Internal::accumulate_functor_vector<Precision, std::greater<Precision> > vector_accumulate_functor;
00746         return Internal::accumulate<Size,Precision,Base,
00747             vector_accumulate_functor, Precision >( v ); 
00748     }
00749 
00750     /// Finds the smallest value of a matrix.
00751     /// @param m a matrix
00752     /// @return the smallest value of m
00753     template<int R, int C, typename Precision, typename Base> inline Precision min_value( const Matrix<R, C, Precision, Base> & m) {
00754         typedef Internal::accumulate_functor_matrix<Precision, std::less<Precision> > matrix_accumulate_functor;
00755         return Internal::accumulate<R,C,Precision,Base,
00756             matrix_accumulate_functor, Precision>( m );
00757     }
00758     /// Finds the largest value of a matrix.
00759     /// @param m a matrix
00760     /// @return the largest value of m
00761     template<int R, int C, typename Precision, typename Base> inline Precision max_value( const Matrix<R, C, Precision, Base> & m) {
00762         typedef Internal::accumulate_functor_matrix<Precision, std::greater<Precision> > matrix_accumulate_functor;
00763         return Internal::accumulate<R,C,Precision,Base,
00764             matrix_accumulate_functor, Precision>( m );
00765     }
00766     /// Finds the smallest values of each column of a matrix.
00767     /// @param m a matrix
00768     /// @return a vector of size C
00769     template<int R, int C, typename Precision, typename Base> inline Vector<Dynamic,Precision> min_value_vertical( const Matrix<R, C, Precision, Base> & m) {
00770         typedef Internal::accumulate_vertical_functor<Precision,std::less<Precision> > matrix_accumulate_vertical_functor;
00771         return Internal::accumulate_vertical<R,C,Precision,Base,
00772             matrix_accumulate_vertical_functor, Vector<Dynamic,Precision> >( m );
00773     }
00774     /// Finds the largest values of each column of a matrix.
00775     /// @param m a matrix
00776     /// @return a vector of size C
00777     template<int R, int C, typename Precision, typename Base> inline Vector<Dynamic,Precision> max_value_vertical( const Matrix<R, C, Precision, Base> & m) {
00778         typedef Internal::accumulate_vertical_functor<Precision,std::greater<Precision> > matrix_accumulate_vertical_functor;
00779         return Internal::accumulate_vertical<R,C,Precision,Base,
00780             matrix_accumulate_vertical_functor, Vector<Dynamic,Precision> >( m );
00781     }
00782     /// Finds the smallest values of each row of a matrix.
00783     /// @param m a matrix
00784     /// @return a vector of size R
00785     template<int R, int C, typename Precision, typename Base> inline Vector<Dynamic,Precision> min_value_horizontal( const Matrix<R, C, Precision, Base> & m) {
00786         typedef Internal::accumulate_horizontal_functor<Precision,std::less<Precision> > matrix_accumulate_horizontal_functor;
00787         return Internal::accumulate_horizontal<R,C,Precision,Base,
00788             matrix_accumulate_horizontal_functor, Vector<Dynamic,Precision> >( m );
00789     }
00790     /// Finds the largest values of each row of a matrix.
00791     /// @param m a matrix
00792     /// @return a vector of size R
00793     template<int R, int C, typename Precision, typename Base> inline Vector<Dynamic,Precision> max_value_horizontal( const Matrix<R, C, Precision, Base> & m) {
00794         typedef Internal::accumulate_horizontal_functor<Precision,std::greater<Precision> > matrix_accumulate_horizontal_functor;
00795         return Internal::accumulate_horizontal<R,C,Precision,Base,
00796             matrix_accumulate_horizontal_functor, Vector<Dynamic,Precision> >( m );
00797     }
00798     /// Finds the smallest value of a vector and its index.
00799     /// @param v a vector
00800     /// @return a pair containing the smallest value and its index
00801     template<int Size, typename Precision, typename Base> inline std::pair<Precision,int> min_element( const Vector<Size, Precision, Base>& v) {
00802         typedef Internal::accumulate_element_functor_vector<Precision, std::less<Precision> > vector_accumulate_functor;
00803         return Internal::accumulate<Size,Precision,Base,
00804             vector_accumulate_functor, std::pair<Precision,int> >( v ); 
00805     }
00806     /// Finds the largest value of a vector and its index.
00807     /// @param v a vector
00808     /// @return a pair containing the largest value and its index
00809     template<int Size, typename Precision, typename Base> inline std::pair<Precision,int> max_element( const Vector<Size, Precision, Base>& v) {
00810         typedef Internal::accumulate_element_functor_vector<Precision, std::greater<Precision> > vector_accumulate_functor;
00811         return Internal::accumulate<Size,Precision,Base,
00812             vector_accumulate_functor, std::pair<Precision,int> >( v ); 
00813     }    
00814     /// Finds the smallest value of a matrix and its row and column.
00815     /// @param m a matrix
00816     /// @return a pair containing the smallest value and a pair
00817     /// containing its row and column
00818     template<int R, int C, typename Precision, typename Base> inline std::pair<Precision,std::pair<int,int> > min_element( const Matrix<R, C, Precision, Base> & m) {
00819         typedef Internal::accumulate_element_functor_matrix<Precision, std::less<Precision> > matrix_accumulate_functor;
00820         typedef std::pair<Precision,std::pair<int,int> > Ret;
00821         return Internal::accumulate<R,C,Precision,Base,
00822             matrix_accumulate_functor, Ret>( m );
00823     }
00824     /// Finds the largest value of a matrix and its row and column.
00825     /// @param m a matrix
00826     /// @return a pair containing the largest value and a pair
00827     /// containing its row and column
00828     template<int R, int C, typename Precision, typename Base> inline std::pair<Precision,std::pair<int,int> > max_element( const Matrix<R, C, Precision, Base> & m) {
00829         typedef Internal::accumulate_element_functor_matrix<Precision, std::greater<Precision> > matrix_accumulate_functor;
00830         typedef std::pair<Precision,std::pair<int,int> > Ret;
00831         return Internal::accumulate<R,C,Precision,Base,
00832             matrix_accumulate_functor, Ret>( m );
00833     }   
00834     /// Finds the smallest values of each column of a matrix and their
00835     /// indices.
00836     /// @param m a matrix
00837     /// @return a pair of vectors of size C containg the values and
00838     /// their indices
00839     template<int R, int C, typename Precision, typename Base> inline std::pair<Vector<Dynamic,Precision>,Vector<Dynamic,Precision> > min_element_vertical( const Matrix<R, C, Precision, Base> & m) {
00840         typedef Internal::accumulate_element_vertical_functor<Precision,std::less<Precision> > matrix_accumulate_vertical_functor;
00841         typedef std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > Ret;
00842         return Internal::accumulate_vertical<R,C,Precision,Base,
00843             matrix_accumulate_vertical_functor, Ret >( m );
00844     }
00845     /// Finds the largest values of each column of a matrix and their
00846     /// indices.
00847     /// @param m a matrix
00848     /// @return a pair of vectors of size C containg the values and
00849     /// their indices
00850     template<int R, int C, typename Precision, typename Base> inline std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > max_element_vertical( const Matrix<R, C, Precision, Base> & m) {
00851         typedef Internal::accumulate_element_vertical_functor<Precision,std::greater<Precision> > matrix_accumulate_vertical_functor;
00852         typedef std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > Ret;
00853         return Internal::accumulate_vertical<R,C,Precision,Base,
00854             matrix_accumulate_vertical_functor, Ret >( m );
00855     }
00856     /// Finds the smallest values of each row of a matrix and their
00857     /// indices.
00858     /// @param m a matrix
00859     /// @return a pair of vectors of size R containg the values and
00860     /// their indices
00861     template<int R, int C, typename Precision, typename Base> inline std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > min_element_horizontal( const Matrix<R, C, Precision, Base> & m) {
00862         typedef Internal::accumulate_element_horizontal_functor<Precision,std::less<Precision> > matrix_accumulate_vertical_functor;
00863         typedef std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > Ret;
00864         return Internal::accumulate_horizontal<R,C,Precision,Base,
00865             matrix_accumulate_vertical_functor, Ret >( m );
00866     }
00867     /// Finds the largest values of each row of a matrix and their
00868     /// indices.
00869     /// @param m a matrix
00870     /// @return a pair of vectors of size R containg the values and
00871     /// their indices
00872     template<int R, int C, typename Precision, typename Base> inline std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > max_element_horizontal( const Matrix<R, C, Precision, Base> & m) {
00873         typedef Internal::accumulate_element_horizontal_functor<Precision,std::greater<Precision> > matrix_accumulate_vertical_functor;
00874         typedef std::pair<Vector< Dynamic, Precision >,Vector< Dynamic, Precision > > Ret;
00875         return Internal::accumulate_horizontal<R,C,Precision,Base,
00876             matrix_accumulate_vertical_functor, Ret >( m );
00877     }
00878 }
00879 #endif