The MRPT project: ops_containers.h Source File

Go to the documentation of this file.
00001 /* +---------------------------------------------------------------------------+
00002    |          The Mobile Robot Programming Toolkit (MRPT) C++ library          |
00003    |                                                                           |
00004    |                   http://mrpt.sourceforge.net/                            |
00005    |                                                                           |
00006    |   Copyright (C) 2005-2010  University of Malaga                           |
00007    |                                                                           |
00008    |    This software was written by the Machine Perception and Intelligent    |
00009    |      Robotics Lab, University of Malaga (Spain).                          |
00010    |    Contact: Jose-Luis Blanco  <jlblanco@ctima.uma.es>                     |
00011    |                                                                           |
00012    |  This file is part of the MRPT project.                                   |
00013    |                                                                           |
00014    |     MRPT is free software: you can redistribute it and/or modify          |
00015    |     it under the terms of the GNU General Public License as published by  |
00016    |     the Free Software Foundation, either version 3 of the License, or     |
00017    |     (at your option) any later version.                                   |
00018    |                                                                           |
00019    |   MRPT is distributed in the hope that it will be useful,                 |
00020    |     but WITHOUT ANY WARRANTY; without even the implied warranty of        |
00021    |     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         |
00022    |     GNU General Public License for more details.                          |
00023    |                                                                           |
00024    |     You should have received a copy of the GNU General Public License     |
00025    |     along with MRPT.  If not, see <http://www.gnu.org/licenses/>.         |
00026    |                                                                           |
00027    +---------------------------------------------------------------------------+ */
00028 #ifndef  mrpt_math_container_ops_H
00029 #define  mrpt_math_container_ops_H
00030 
00031 #include <mrpt/math/math_frwds.h>  // Fordward declarations
00032 
00033 #include <functional>
00034 #include <algorithm>
00035 #include <cmath>
00036 
00037 /** \file ops_containers.h
00038   * This file implements several operations that operate element-wise on individual or pairs of containers.
00039   *  Containers here means any of: mrpt::math::CVectorTemplace, mrpt::math::CArray, mrpt::math::CMatrixFixedNumeric, mrpt::math::CMatrixTemplate.
00040   *
00041   *  In general, any container having a type "mrpt_autotype" self-referencing to the type itself, and a dummy struct mrpt_container<>
00042   *   which is only used as a way to force the compiler to assure that BOTH containers are valid ones in binary operators.
00043   *   This restrictions
00044   *   have been designed as a way to provide "polymorphism" at a template level, so the "+,-,..." operators do not
00045   *   generate ambiguities for ANY type, and limiting them to MRPT containers.
00046   *
00047   *   In some cases, the containers provide specializations of some operations, for increased performance.
00048   */
00049 
00050 #include <algorithm>
00051 #include <numeric>
00052 #include <functional>
00053 
00054 #include <mrpt/math/CHistogram.h>  // Used in ::histogram()
00055 
00056 namespace mrpt
00057 {
00058         namespace math
00059         {
00060 
00061                 /** \name Generic container element-wise operations - Basic arithmetic
00062                   * @{
00063                   */
00064 
00065                 /** a+=b */
00066                 template <class CONTAINER1,class CONTAINER2>
00067                 inline RET_CONT1_ASSERT_MRPTCONTAINERS(CONTAINER1,CONTAINER2) &
00068                 operator +=(CONTAINER1 &a, const CONTAINER2 &b)
00069                 {
00070                         ASSERT_(a.size()==b.size());
00071                         typename CONTAINER1::iterator ita = a.begin();
00072                         typename CONTAINER2::const_iterator itb = b.begin();
00073                         const typename CONTAINER1::iterator last=a.end();
00074                         while (ita!=last) { *(ita++)+=*(itb++); }
00075                         return a;
00076                 }
00077 
00078                 /** return a+b */
00079                 template <class CONTAINER1,class CONTAINER2>
00080                 inline RET_CONT1_ASSERT_MRPTCONTAINERS(CONTAINER1,CONTAINER2)
00081                 operator +(const CONTAINER1 &a, const CONTAINER2 &b)
00082                 {
00083                         CONTAINER1 ret = a;
00084                         ret+=b;
00085                         return ret;
00086                 }
00087 
00088                 /** a-=b */
00089                 template <class CONTAINER1,class CONTAINER2>
00090                 inline RET_CONT1_ASSERT_MRPTCONTAINERS(CONTAINER1,CONTAINER2) &
00091                 operator -=(CONTAINER1 &a, const CONTAINER2 &b)
00092                 {
00093                         ASSERT_(a.size()==b.size());
00094                         typename CONTAINER1::iterator ita = a.begin();
00095                         typename CONTAINER2::const_iterator itb = b.begin();
00096                         const typename CONTAINER1::iterator last=a.end();
00097                         while (ita!=last) { *(ita++)-=*(itb++); }
00098                         return a;
00099                 }
00100 
00101                 /** return a-b */
00102                 template <class CONTAINER1,class CONTAINER2>
00103                 inline RET_CONT1_ASSERT_MRPTCONTAINERS(CONTAINER1,CONTAINER2)
00104                 operator -(const CONTAINER1 &a, const CONTAINER2 &b)
00105                 {
00106                         CONTAINER1 ret = a;
00107                         ret-=b;
00108                         return ret;
00109                 }
00110 
00111                 /** a/=b (element-wise division) */
00112                 template<class CONTAINER1,class CONTAINER2> CONTAINER1 &divide_elementwise(CONTAINER1 &a, const CONTAINER2 &b)
00113                 {
00114                         ASSERT_(a.size()==b.size());
00115                         typename CONTAINER1::iterator ita = a.begin();
00116                         typename CONTAINER2::const_iterator itb = b.begin();
00117                         const typename CONTAINER1::iterator last=a.end();
00118                         while (ita!=last) { *(ita++)/=*(itb++); }
00119                         return a;
00120                 }
00121 
00122                 /** Sum a scalar to all the elements */
00123                 template <class CONTAINER>
00124                 inline RET_CONT_ASSERT_MRPTCONTAINER(CONTAINER) & operator+=(CONTAINER &c, typename CONTAINER::value_type  val)
00125                 {
00126                         const typename CONTAINER::iterator last=c.end();
00127                         for (typename CONTAINER::iterator it=c.begin();it!=last;++it)  *it += val;
00128                         return c;
00129                 }
00130 
00131                 /** Sum a scalar to all the elements */
00132                 template <class CONTAINER>
00133                 inline RET_CONT_ASSERT_MRPTCONTAINER(CONTAINER) operator +(CONTAINER &c, typename CONTAINER::value_type  val)
00134                 {
00135                         CONTAINER ret = c;
00136                         ret+=val;
00137                         return ret;
00138                 }
00139 
00140                 /** Substract a scalar to all the elements */
00141                 template <class CONTAINER>
00142                 inline RET_CONT_ASSERT_MRPTCONTAINER(CONTAINER) & operator -= (CONTAINER &c, typename CONTAINER::value_type  val)
00143                 {
00144                         const typename CONTAINER::iterator last=c.end();
00145                         for (typename CONTAINER::iterator it=c.begin();it!=last;++it)  *it -= val;
00146                         return c;
00147                 }
00148 
00149                 /** Substract a scalar to all the elements */
00150                 template <class CONTAINER>
00151                 inline RET_CONT_ASSERT_MRPTCONTAINER(CONTAINER) operator -(CONTAINER &c, typename CONTAINER::value_type  val)
00152                 {
00153                         CONTAINER ret = c;
00154                         ret-=val;
00155                         return ret;
00156                 }
00157 
00158                 /** Multiplies all the elements by a scalar C *= s */
00159                 template <class CONTAINER>
00160                 inline RET_CONT_ASSERT_MRPTCONTAINER(CONTAINER) & operator *= (CONTAINER &c, typename CONTAINER::value_type  val)
00161                 {
00162                         const typename CONTAINER::iterator last=c.end();
00163                         for (typename CONTAINER::iterator it=c.begin();it!=last;++it)  *it *= val;
00164                         return c;
00165                 }
00166 
00167                 /** Multiplies all the elements by a scalar C2 = C*s */
00168                 template <class CONTAINER>
00169                 inline RET_CONT_ASSERT_MRPTCONTAINER(CONTAINER) operator *(const CONTAINER &c, const typename CONTAINER::value_type  val)
00170                 {
00171                         CONTAINER ret = c;
00172                         ret*=val;
00173                         return ret;
00174                 }
00175                 /** Multiplies all the elements by a scalar C2= s*C */
00176                 template <class CONTAINER>
00177                 inline RET_CONT_ASSERT_MRPTCONTAINER(CONTAINER) operator *(const typename CONTAINER::value_type  val, const CONTAINER &c)
00178                 {
00179                         CONTAINER ret = c;
00180                         ret*=val;
00181                         return ret;
00182                 }
00183 
00184                 /** Divides all the elements by a scalar */
00185                 template <class CONTAINER>
00186                 inline RET_CONT_ASSERT_MRPTCONTAINER(CONTAINER) & operator /= (CONTAINER &c, typename CONTAINER::value_type  val)
00187                 {
00188                         typename CONTAINER::value_type k = 1/val; // Multiply is cheaper..
00189                         const typename CONTAINER::iterator last=c.end();
00190                         for (typename CONTAINER::iterator it=c.begin();it!=last;++it)  *it *= k;
00191                         return c;
00192                 }
00193 
00194                 /** Divides all the elements by a scalar */
00195                 template <class CONTAINER>
00196                 inline RET_CONT_ASSERT_MRPTCONTAINER(CONTAINER) operator /(CONTAINER &c, typename CONTAINER::value_type  val)
00197                 {
00198                         CONTAINER ret = c;
00199                         ret/=val;
00200                         return ret;
00201                 }
00202 
00203 
00204                 /** @} end Arithmetic  */
00205 
00206 
00207                 /** \name Generic container element-wise operations - Miscelaneous
00208                   * @{
00209                   */
00210 
00211                 /** A template for counting how many elements in a container are non-Zero. */
00212                 template <class CONTAINER>
00213                 size_t countNonZero(const CONTAINER &a)
00214                 {
00215                         typename CONTAINER::const_iterator      it_a;
00216                         size_t count=0;
00217                         for (it_a=a.begin(); it_a!=a.end(); it_a++) if (*it_a) count++;
00218                         return count;
00219                 }
00220 
00221                 /** Finds the maximum value (and the corresponding zero-based index) from a given container.
00222                   */
00223                 template <class CONTAINER>
00224                 typename CONTAINER::value_type
00225                 maximum(const CONTAINER &v, size_t *maxIndex)
00226                 {
00227                         typename CONTAINER::const_iterator maxIt = std::max_element(v.begin(),v.end());
00228                         if (maxIndex) *maxIndex = std::distance(v.begin(),maxIt);
00229                         return *maxIt;
00230                 }
00231 
00232                 /** Finds the maximum value (and the corresponding zero-based index) from a given vector.
00233                   * \sa maximum, minimum_maximum
00234                   */
00235                 template <class CONTAINER>
00236                 typename CONTAINER::value_type
00237                 minimum(const CONTAINER &v, size_t *minIndex)
00238                 {
00239                         typename CONTAINER::const_iterator minIt = std::min_element(v.begin(),v.end());
00240                         if (minIndex) *minIndex = std::distance(v.begin(),minIt);
00241                         return *minIt;
00242                 }
00243 
00244                 /** Compute the minimum and maximum of a vector at once
00245                   * \sa maximum, minimum
00246                   */
00247                 template <class CONTAINER>
00248                 void minimum_maximum(
00249                         const CONTAINER &v,
00250                         typename CONTAINER::mrpt_autotype::value_type & out_min,
00251                         typename CONTAINER::mrpt_autotype::value_type& out_max,
00252                         size_t *minIndex,
00253                         size_t *maxIndex)
00254                 {
00255                         const size_t N = v.size();
00256                         if (N)
00257                         {
00258                                 typename CONTAINER::const_iterator it = v.begin();
00259                                 size_t min_idx=0,max_idx=0;
00260                                 out_max = out_min = *it;
00261                                 size_t i;
00262                                 for (i=0;i<N;++it,++i)
00263                                 {
00264                                         if (*it<out_min)
00265                                         {
00266                                                 out_min=*it;
00267                                                 min_idx = i;
00268                                         }
00269                                         if (*it>out_max)
00270                                         {
00271                                                 out_max=*it;
00272                                                 max_idx = i;
00273                                         }
00274                                 }
00275                                 if (minIndex) *minIndex = min_idx;
00276                                 if (maxIndex) *maxIndex = max_idx;
00277                         }
00278                         else {
00279                                 out_min=out_max=0;
00280                                 if (minIndex) *minIndex = 0;
00281                                 if (maxIndex) *maxIndex = 0;
00282                         }
00283                 }
00284                 /** Compute the norm-infinite of a vector ($f[ ||\mathbf{v}||_\infnty $f]), ie the maximum absolute value of the elements.
00285                   */
00286                 template <class CONTAINER>
00287                 typename CONTAINER::value_type
00288                 norm_inf(const CONTAINER &v, size_t *maxIndex)
00289                 {
00290                         double  M=0;
00291                         int             i,M_idx=-1;
00292                         typename CONTAINER::const_iterator      it;
00293                         for (i=0, it=v.begin(); it!=v.end();it++,i++)
00294                         {
00295                                 double  it_abs = fabs( static_cast<double>(*it));
00296                                 if (it_abs>M || M_idx==-1)
00297                                 {
00298                                         M = it_abs;
00299                                         M_idx = i;
00300                                 }
00301                         }
00302                         if (maxIndex) *maxIndex = M_idx;
00303                         return static_cast<typename CONTAINER::mrpt_autotype::value_type>(M);
00304                 }
00305 
00306 
00307                 /** Compute the square norm of a vector/array/matrix (the Euclidean distance to the origin, taking all the elements as a single vector).
00308                   \sa norm */
00309                 template <class CONTAINER>
00310                 typename CONTAINER::value_type
00311                 squareNorm(const CONTAINER &v)
00312                 {
00313                         typename CONTAINER::value_type  total=0;
00314                         typename CONTAINER::const_iterator      it;
00315                         for (it=v.begin(); it!=v.end();it++)
00316                                 total += square(*it);
00317                         return total;
00318                 }
00319 
00320                 /** Compute the square norm of anything implementing [].
00321                   \sa norm */
00322                 template<size_t N,class T,class U>
00323                 inline T squareNorm(const U &v) {
00324                         T res=0;
00325                         for (size_t i=0;i<N;i++) res+=square(v[i]);
00326                         return res;
00327                 }
00328 
00329                 /** Compute the norm of a vector/array/matrix (the Euclidean distance to the origin, taking all the elements as a single vector).
00330                   \sa squareNorm */
00331                 template <class CONTAINER>
00332                 inline typename CONTAINER::value_type
00333                 norm(const CONTAINER &v)
00334                 {
00335                         return ::sqrt(squareNorm(v));
00336                 }
00337 
00338                 /** v1·v2: The dot product of two containers (vectors/arrays/matrices) */
00339                 template <class CONTAINER1,class CONTAINER2>
00340                 inline typename CONTAINER1::value_type
00341                 dotProduct(const CONTAINER1 &v1,const CONTAINER1 &v2)
00342                 {
00343                         ASSERT_(v1.size()==v2.size());
00344                         const size_t N = v1.size();
00345                         typename CONTAINER1::value_type res=0;
00346                         typename CONTAINER1::const_iterator     it1;
00347                         typename CONTAINER2::const_iterator     it2;
00348                         size_t i;
00349                         for (i=0,it1=v1.begin(),it2=v2.begin();i<N;++i,++it1,++it2) res+=(*it1)*(*it2);
00350                         return res;
00351                 }
00352 
00353                 /** v1·v2: The dot product of any two objects supporting []  */
00354                 template<size_t N,class T,class U,class V>
00355                 inline T dotProduct(const U &v1,const V &v2)    {
00356                         T res=0;
00357                         for (size_t i=0;i<N;i++) res+=v1[i]*v2[i];
00358                         return res;
00359                 }
00360 
00361                 /** Computes the mean value of a vector  \return The mean, as a double number.
00362                   * \sa math::stddev,math::meanAndStd  */
00363                 template <class CONTAINER>
00364                 inline double
00365                 mean(const CONTAINER &v)
00366                 {
00367                         if (v.empty()) return 0;
00368                         else    return static_cast<typename CONTAINER::mrpt_autotype::value_type>( std::accumulate(v.begin(),v.end(), static_cast<typename CONTAINER::mrpt_autotype::value_type>(0) ) ) / double(v.size());
00369                 }
00370 
00371                 /** Computes the sum of all the elements  \sa cumsum  */
00372                 template <class CONTAINER>
00373                 inline typename CONTAINER::value_type
00374                 sum(const CONTAINER &v)
00375                 {
00376                         return std::accumulate(v.begin(),v.end(), static_cast<typename CONTAINER::mrpt_autotype::value_type>(0) );
00377                 }
00378 
00379                 /** Computes the cumulative sum of all the elements, saving the result in another container.
00380                   *  This works for both matrices (even mixing their types) and vectores/arrays (even mixing types),
00381                   *  and even to store the cumsum of any matrix into any vector/array, but not in opposite direction.
00382                   * \sa sum */
00383                 template <class CONTAINER1,class CONTAINER2>
00384                 void cumsum(const CONTAINER1 &in_data, CONTAINER2 &out_cumsum)
00385                 {
00386                         out_cumsum.resize(in_data.size());
00387                         typename CONTAINER1::value_type last=0;
00388                         const size_t N = in_data.size();
00389                         typename CONTAINER1::const_iterator     it1;
00390                         typename CONTAINER2::iterator   it2;
00391                         size_t i;
00392                         for (i=0,it1=in_data.begin(),it2=out_cumsum.begin();i<N;++i,++it1,++it2)
00393                                 last = (*it2) = last + (*it1);
00394                 }
00395 
00396                 /** Computes the cumulative sum of all the elements
00397                   * \sa sum  */
00398                 template<class CONTAINER>
00399                 inline CONTAINER
00400                 cumsum(const CONTAINER &in_data)
00401                 {
00402                         CONTAINER ret;
00403                         ret.resize(in_data.size());
00404                         cumsum(in_data,ret);
00405                         return ret;
00406                 }
00407 
00408                 /** Counts the number of elements that appear in both containers (comparison through the == operator)
00409                   *  It is assumed that no repeated elements appear within each of the containers.  */
00410                 template <class CONTAINER1,class CONTAINER2>
00411                 size_t  countCommonElements(const CONTAINER1 &a,const CONTAINER2 &b)
00412                 {
00413                     size_t ret=0;
00414                         for (typename CONTAINER1::const_iterator it1 = a.begin();it1!=a.end();it1++)
00415                             for (typename CONTAINER2::const_iterator it2 = b.begin();it2!=b.end();it2++)
00416                     if ( (*it1) == (*it2) )
00417                          ret++;
00418                         return ret;
00419                 }
00420 
00421                 /** Computes the normalized or normal histogram of a sequence of numbers given the number of bins and the limits.
00422                   *  In any case this is a "linear" histogram, i.e. for matrices, all the elements are taken as if they were a plain sequence, not taking into account they were in columns or rows.
00423                   */
00424                 template<class CONTAINER>
00425                 std::vector<double> histogram(
00426                         const CONTAINER &v,
00427                         double limit_min,
00428                         double limit_max,
00429                         size_t number_bins,
00430                         bool do_normalization)
00431                 {
00432                         mrpt::math::CHistogram  H( limit_min, limit_max, number_bins );
00433                         std::vector<double> ret(number_bins);
00434                         for (typename CONTAINER::const_iterator i=v.begin();i!=v.end();++i)
00435                                 H.add(static_cast<double>(*i));
00436                         for (size_t i=0;i<number_bins;i++)
00437                                 ret[i] = do_normalization ? H.getBinRatio(i) : H.getBinCount(i);
00438                         return ret;
00439                 }
00440 
00441                 /** Adjusts the range of all the elements such as the minimum and maximum values being those supplied by the user.  */
00442                 template <class CONTAINER>
00443                 void  adjustRange(CONTAINER &m, const typename CONTAINER::value_type minVal,const typename CONTAINER::value_type maxVal)
00444                 {
00445                         if (size_t(m.size())==0) return;
00446                         typename CONTAINER::value_type curMin,curMax;
00447                         minimum_maximum(m,curMin,curMax);
00448                         const typename CONTAINER::value_type curRan = curMax-curMin;
00449                         m -= (curMin+minVal);
00450                         if (curRan!=0) m *= (maxVal-minVal)/curRan;
00451                 }
00452 
00453                 /** @} Misc ops */
00454 
00455                 /** \name Generic container element-wise operations - Mathematical functions
00456                   * @{
00457                   */
00458                 /** out = abs(in), element by element */
00459                 template <class CONTAINER1,class CONTAINER2>
00460                 void Abs(const CONTAINER1 &dat_in, CONTAINER2 &dat_out) {
00461                         dat_out.resize(dat_in.size());
00462                         std::transform<typename CONTAINER1::const_iterator,typename CONTAINER2::iterator,typename CONTAINER1::value_type (*)(typename CONTAINER1::value_type)>
00463                                 (dat_in.begin(),dat_in.end(), dat_out.begin(), &::abs );
00464                 }
00465                 /** return abs(in), element by element */
00466                 template <class CONTAINER>  inline CONTAINER Abs(const CONTAINER &dat_in)
00467                 { CONTAINER ret; Abs(dat_in,ret); return ret; }
00468 
00469                 /** out = log(in), element by element */
00470                 template <class CONTAINER1,class CONTAINER2>
00471                 void Log(const CONTAINER1 &dat_in, CONTAINER2 &dat_out) {
00472                         dat_out.resize(dat_in.size());
00473                         std::transform<typename CONTAINER1::const_iterator,typename CONTAINER2::iterator,typename CONTAINER1::value_type (*)(typename CONTAINER1::value_type)>
00474                                 (dat_in.begin(),dat_in.end(), dat_out.begin(), &::log );
00475                 }
00476                 /** return log(in), element by element */
00477                 template <class CONTAINER>  inline CONTAINER Log(const CONTAINER &dat_in)
00478                 { CONTAINER ret; Log(dat_in,ret); return ret; }
00479 
00480                 /** out = exp(in), element by element */
00481                 template <class CONTAINER1,class CONTAINER2>
00482                 void Exp(const CONTAINER1 &dat_in, CONTAINER2 &dat_out) {
00483                         dat_out.resize(dat_in.size());
00484                         std::transform<typename CONTAINER1::const_iterator,typename CONTAINER2::iterator,typename CONTAINER1::value_type (*)(typename CONTAINER1::value_type)>
00485                                 (dat_in.begin(),dat_in.end(), dat_out.begin(), &::exp );
00486                 }
00487                 /** return exp(in), element by element */
00488                 template <class CONTAINER>  inline CONTAINER Exp(const CONTAINER &dat_in)
00489                 { CONTAINER ret; Exp(dat_in,ret); return ret; }
00490 
00491 
00492                 /** Logical equal-to == operator between any pair of matrix or vectors */
00493                 template <class CONTAINER1,class CONTAINER2>
00494                 RET_TYPE_ASSERT_MRPTCONTAINER(CONTAINER1, bool)
00495                 operator == (const CONTAINER1 &m1, const CONTAINER2 &m2)
00496                 {
00497                    if (m1.size()!=m2.size()) return false;
00498                    typename CONTAINER1::const_iterator it1=m1.begin();
00499                    typename CONTAINER2::const_iterator it2=m2.begin();
00500                    const size_t N = m1.size();
00501                    for(size_t i=0;i<N;++i, ++it1,++it2)
00502                                 if (*it1!=*it2)
00503                                          return false;
00504                    return true;
00505                 }
00506 
00507                 /** Logical not-equal-to != operator between any pair of matrix or vectors  */
00508                 template <class CONTAINER1,class CONTAINER2>
00509                 inline RET_TYPE_ASSERT_MRPTCONTAINER(CONTAINER1, bool)
00510                 operator != (const CONTAINER1 &m1, const CONTAINER2 &m2)
00511                 {
00512                         return !(m1 == m2);
00513                 }
00514 
00515                 /** Computes the standard deviation of a vector
00516                   * \param v The set of data
00517                   * \param unbiased If set to true or false the std is normalized by "N-1" or "N", respectively.
00518                   * \sa math::mean,math::meanAndStd
00519                   */
00520                 template<class VECTORLIKE>
00521                 double  stddev(const VECTORLIKE &v, bool unbiased)
00522                 {
00523                         if (v.size()<2)
00524                                 return 0;
00525                         else
00526                         {
00527                                 // Compute the mean:
00528                                 double   vector_std=0,vector_mean = 0;
00529                                 for (typename VECTORLIKE::const_iterator it = v.begin();it!=v.end();++it) vector_mean += (*it);
00530                                 vector_mean /= static_cast<double>(v.size());
00531                                 // Compute the std:
00532                                 for (typename VECTORLIKE::const_iterator it = v.begin();it!=v.end();++it) vector_std += square((*it)-vector_mean);
00533                                 vector_std = sqrt(vector_std  / static_cast<double>(v.size() - (unbiased ? 1:0)) );
00534                                 return vector_std;
00535                         }
00536                 }
00537 
00538                 /** Computes the mean vector and covariance from a list of values given as a vector of vectors, where each row is a sample.
00539                   * \param v The set of data, as a vector of N vectors of M elements.
00540                   * \param out_mean The output M-vector for the estimated mean.
00541                   * \param out_cov The output MxM matrix for the estimated covariance matrix.
00542                   * \sa math::mean,math::stddev, math::cov
00543                   */
00544                 template<class VECTOR_OF_VECTOR, class VECTORLIKE, class MATRIXLIKE>
00545                 void  meanAndCov(
00546                         const VECTOR_OF_VECTOR &v,
00547                         VECTORLIKE &out_mean,
00548                         MATRIXLIKE &out_cov
00549                         )
00550                 {
00551                         const size_t N = v.size();
00552                         ASSERTMSG_(N>0,"The input vector contains no elements");
00553                         const double N_inv = 1.0/N;
00554 
00555                         const size_t M = v[0].size();
00556                         ASSERTMSG_(M>0,"The input vector contains rows of length 0");
00557 
00558                         // First: Compute the mean
00559                         out_mean.assign(M,0);
00560                         for (size_t i=0;i<N;i++)
00561                                 for (size_t j=0;j<M;j++)
00562                                         out_mean[j]+=v[i][j];
00563                         out_mean*=N_inv;
00564 
00565                         // Second: Compute the covariance
00566                         //  Save only the above-diagonal part, then after averaging
00567                         //  duplicate that part to the other half.
00568                         out_cov.zeros(M,M);
00569                         for (size_t i=0;i<N;i++)
00570                         {
00571                                 for (size_t j=0;j<M;j++)
00572                                         out_cov.get_unsafe(j,j)+=square(v[i][j]-out_mean[j]);
00573 
00574                                 for (size_t j=0;j<M;j++)
00575                                         for (size_t k=j+1;k<M;k++)
00576                                                 out_cov.get_unsafe(j,k)+=(v[i][j]-out_mean[j])*(v[i][k]-out_mean[k]);
00577                         }
00578                         for (size_t j=0;j<M;j++)
00579                                 for (size_t k=j+1;k<M;k++)
00580                                         out_cov.get_unsafe(k,j) = out_cov.get_unsafe(j,k);
00581                         out_cov*=N_inv;
00582                 }
00583 
00584 
00585                 /** Computes the standard deviation of a vector
00586                   * \param v The set of data
00587                   * \param out_mean The output for the estimated mean
00588                   * \param out_std The output for the estimated standard deviation
00589                   * \param unbiased If set to true or false the std is normalized by "N-1" or "N", respectively.
00590                   * \sa math::mean,math::stddev
00591                   */
00592                 template<class VECTORLIKE>
00593                         void  meanAndStd(
00594                                 const VECTORLIKE &v,
00595                                 double                  &out_mean,
00596                                 double                  &out_std,
00597                                 bool                    unbiased)
00598                 {
00599                         if (v.size()<2)
00600                         {
00601                                 out_std = 0;
00602                                 out_mean = (v.size()==1) ? *v.begin() : 0;
00603                         }
00604                         else
00605                         {
00606                                 // Compute the mean:
00607                                 typename VECTORLIKE::const_iterator it;
00608                                 out_std=0,out_mean = 0;
00609                                 for (it = v.begin();it!=v.end();it++) out_mean += (*it);
00610                                 out_mean /= static_cast<double>(v.size());
00611 
00612                                 // Compute the std:
00613                                 for (it = v.begin();it!=v.end();it++) out_std += mrpt::utils::square(static_cast<double>(*it)-out_mean);
00614                                 out_std = std::sqrt(out_std / static_cast<double>((v.size() - (unbiased ? 1:0)) ));
00615                         }
00616                 }
00617 
00618                 /** @} Math ops */
00619 
00620         } // End of math namespace
00621 } // End of mrpt namespace
00622 
00623 
00624 #endif