9 #ifndef  MRPT_DATA_UTILS_MATH_H    10 #define  MRPT_DATA_UTILS_MATH_H    33                 template<
class VECTORLIKE1,
class VECTORLIKE2,
class MAT>
    36                         const VECTORLIKE2 &MU,
    40                         #if defined(_DEBUG) || (MRPT_ALWAYS_CHECKS_DEBUG_MATRICES)    45                         const size_t N = X.size();
    46                         Eigen::Matrix<typename MAT::Scalar,Eigen::Dynamic,1> X_MU(N);
    47                         for (
size_t i=0;i<N;i++) X_MU[i]=X[i]-MU[i];
    48                         const Eigen::Matrix<typename MAT::Scalar,Eigen::Dynamic,1>  
z = COV.llt().solve(X_MU);
    57                 template<
class VECTORLIKE1,
class VECTORLIKE2,
class MAT>
    60                         const VECTORLIKE2 &MU,
    70                 template<
class VECTORLIKE,
class MAT1,
class MAT2,
class MAT3>
    73                         const VECTORLIKE &mean_diffs,
    76                         const MAT3 &CROSS_COV12 )
    79                         #if defined(_DEBUG) || (MRPT_ALWAYS_CHECKS_DEBUG_MATRICES)    82                                 ASSERT_( COV1.isSquare() && COV2.isSquare() );
    87                         COV.substract_An(CROSS_COV12,2);
    89                         COV.inv_fast(COV_inv);
    97                 template<
class VECTORLIKE,
class MAT1,
class MAT2,
class MAT3> 
inline typename VECTORLIKE::Scalar    99                         const VECTORLIKE &mean_diffs,
   102                         const MAT3 &CROSS_COV12 )
   110                 template<
class VECTORLIKE,
class MATRIXLIKE>
   115                         ASSERTDEB_(
size_t(
cov.getColCount())==
size_t(delta_mu.size()))
   122                 template<
class VECTORLIKE,
class MATRIXLIKE>
   132                 template <
typename T>
   134                         const std::vector<T> &mean_diffs,
   139                         const size_t vector_dim = mean_diffs.size();
   144                         const T cov_det = C.det();
   148                         return std::pow( 
M_2PI, -0.5*vector_dim ) * (1.0/std::sqrt( cov_det ))
   149                                 * exp( -0.5 * mean_diffs.multiply_HCHt_scalar(C_inv) );
   155                 template <
typename T, 
size_t DIM>
   157                         const std::vector<T> &mean_diffs,
   162                         ASSERT_(mean_diffs.size()==DIM);
   166                         const T cov_det = C.det();
   170                         return std::pow( 
M_2PI, -0.5*DIM ) * (1.0/std::sqrt( cov_det ))
   171                                 * exp( -0.5 * mean_diffs.multiply_HCHt_scalar(C_inv) );
   177                 template <
typename T, 
class VECLIKE,
class MATLIKE1, 
class MATLIKE2>
   179                         const VECLIKE   &mean_diffs,
   180                         const MATLIKE1  &COV1,
   181                         const MATLIKE2  &COV2,
   184                         const MATLIKE1  *CROSS_COV12=NULL
   187                         const size_t vector_dim = mean_diffs.size();
   192                         if (CROSS_COV12) { C-=*CROSS_COV12; C-=*CROSS_COV12; }
   193                         const T cov_det = C.det();
   197                         maha2_out = mean_diffs.multiply_HCHt_scalar(C_inv);
   198                         intprod_out = std::pow( 
M_2PI, -0.5*vector_dim ) * (1.0/std::sqrt( cov_det ))*exp(-0.5*maha2_out);
   204                 template <
typename T, 
class VECLIKE,
class MATRIXLIKE>
   206                         const VECLIKE           &diff_mean,
   207                         const MATRIXLIKE        &
cov,
   213                         ASSERTDEB_(
size_t(
cov.getColCount())==
size_t(diff_mean.size()))
   228                 template <
typename T, 
class VECLIKE,
class MATRIXLIKE>
   230                         const VECLIKE           &diff_mean,
   231                         const MATRIXLIKE        &
cov,
   236                         pdf_out = std::exp(pdf_out); 
   250                 template<
class VECTOR_OF_VECTORS, 
class MATRIXLIKE,
class VECTORLIKE,
class VECTORLIKE2,
class VECTORLIKE3>
   252                         const VECTOR_OF_VECTORS &elements,
   253                         MATRIXLIKE &covariances,
   255                         const VECTORLIKE2 *weights_mean,
   256                         const VECTORLIKE3 *weights_cov,
   257                         const bool *elem_do_wrap2pi = NULL
   260                         ASSERTMSG_(elements.size()!=0,
"No samples provided, so there is no way to deduce the output size.")
   262                         const size_t DIM = elements[0].size();
   264                         covariances.setSize(DIM,DIM);
   265                         const size_t nElms=elements.size();
   266                         const T NORM=1.0/nElms;
   267                         if (weights_mean) { 
ASSERTDEB_(
size_t(weights_mean->size())==
size_t(nElms)) }
   269                         for (
size_t i=0;i<DIM;i++)
   272                                 if (!elem_do_wrap2pi || !elem_do_wrap2pi[i])
   276                                                 for (
size_t j=0;j<nElms;j++)
   277                                                         accum+= (*weights_mean)[j] * elements[j][i];
   281                                                 for (
size_t j=0;j<nElms;j++) accum+=elements[j][i];
   287                                         double accum_L=0,accum_R=0;
   288                                         double Waccum_L=0,Waccum_R=0;
   289                                         for (
size_t j=0;j<nElms;j++)
   291                                                 double ang = elements[j][i];
   292                                                 const double w   = weights_mean!=NULL ? (*weights_mean)[j] : NORM;
   293                                                 if (fabs( ang )>0.5*
M_PI)
   295                                                         if (ang<0) ang = (
M_2PI + ang);
   305                                         if (Waccum_L>0) accum_L /= Waccum_L;  
   306                                         if (Waccum_R>0) accum_R /= Waccum_R;  
   308                                         accum = (accum_L* Waccum_L + accum_R * Waccum_R );      
   313                         for (
size_t i=0;i<DIM;i++)
   314                                 for (
size_t j=0;j<=i;j++)       
   319                                                 ASSERTDEB_(
size_t(weights_cov->size())==
size_t(nElms))
   320                                                 for (
size_t k=0;k<nElms;k++)
   322                                                         const T Ai = (elements[k][i]-means[i]);
   323                                                         const T Aj = (elements[k][j]-means[j]);
   324                                                         if (!elem_do_wrap2pi || !elem_do_wrap2pi[i])
   325                                                                         elem+= (*weights_cov)[k] * Ai * Aj;
   331                                                 for (
size_t k=0;k<nElms;k++)
   333                                                         const T Ai = (elements[k][i]-means[i]);
   334                                                         const T Aj = (elements[k][j]-means[j]);
   335                                                         if (!elem_do_wrap2pi || !elem_do_wrap2pi[i])
   341                                         covariances.get_unsafe(i,j) = elem;
   342                                         if (i!=j) covariances.get_unsafe(j,i)=elem;
   353                 template<
class VECTOR_OF_VECTORS, 
class MATRIXLIKE,
class VECTORLIKE>
   354                 void covariancesAndMean(
const VECTOR_OF_VECTORS &elements,MATRIXLIKE &covariances,VECTORLIKE &means, 
const bool *elem_do_wrap2pi = NULL)
   356                         covariancesAndMeanWeighted<VECTOR_OF_VECTORS,MATRIXLIKE,VECTORLIKE,CVectorDouble,CVectorDouble>(elements,covariances,means,NULL,NULL,elem_do_wrap2pi);
   368                 template<
class VECTORLIKE1,
class VECTORLIKE2>
   370                                 const VECTORLIKE1       &
values,
   373                                 VECTORLIKE2     &out_binCenters,
   374                                 VECTORLIKE2     &out_binValues )
   382                                 unsigned int    nBins = 
static_cast<unsigned>(ceil((
maximum( 
values )-minBin) / binWidth));
   385                                 out_binCenters.resize(nBins);
   386                                 out_binValues.clear(); out_binValues.resize(nBins,0);
   387                                 TNum halfBin = TNum(0.5)*binWidth;;
   388                                 VECTORLIKE2   binBorders(nBins+1,minBin-halfBin);
   389                                 for (
unsigned int i=0;i<nBins;i++)
   391                                         binBorders[i+1] = binBorders[i]+binWidth;
   392                                         out_binCenters[i] = binBorders[i]+halfBin;
   400                                         int idx = 
round(((*itVal)-minBin)/binWidth);
   401                                         if (idx>=
int(nBins)) idx=nBins-1;
   403                                         out_binValues[idx] += *itW;
   408                                         out_binValues /= totalSum;
   422                 template<
class VECTORLIKE1,
class VECTORLIKE2>
   424                         const VECTORLIKE1       &
values,
   425                         const VECTORLIKE1       &log_weights,
   427                         VECTORLIKE2     &out_binCenters,
   428                         VECTORLIKE2     &out_binValues )
   436                         unsigned int    nBins = 
static_cast<unsigned>(ceil((
maximum( 
values )-minBin) / binWidth));
   439                         out_binCenters.resize(nBins);
   440                         out_binValues.clear(); out_binValues.resize(nBins,0);
   441                         TNum halfBin = TNum(0.5)*binWidth;;
   442                         VECTORLIKE2   binBorders(nBins+1,minBin-halfBin);
   443                         for (
unsigned int i=0;i<nBins;i++)
   445                                 binBorders[i+1] = binBorders[i]+binWidth;
   446                                 out_binCenters[i] = binBorders[i]+halfBin;
   450                         const TNum max_log_weight = 
maximum(log_weights);
   453                         for (itVal = 
values.begin(), itW = log_weights.begin(); itVal!=
values.end(); ++itVal, ++itW )
   455                                 int idx = 
round(((*itVal)-minBin)/binWidth);
   456                                 if (idx>=
int(nBins)) idx=nBins-1;
   458                                 const TNum 
w = exp(*itW-max_log_weight);
   459                                 out_binValues[idx] += 
w;
   464                                 out_binValues /= totalSum;
 void mahalanobisDistance2AndPDF(const VECLIKE &diff_mean, const MATRIXLIKE &cov, T &maha2_out, T &pdf_out)
Computes both, the PDF and the square Mahalanobis distance between a point (given by its difference w...
 
double BASE_IMPEXP averageWrap2Pi(const CVectorDouble &angles)
Computes the average of a sequence of angles in radians taking into account the correct wrapping in t...
 
GLboolean GLenum GLenum GLvoid * values
 
size_t size(const MATRIXLIKE &m, const int dim)
 
This file implements miscelaneous matrix and matrix/vector operations, and internal functions in mrpt...
 
MAT_C::Scalar multiply_HCHt_scalar(const VECTOR_H &H, const MAT_C &C)
r (a scalar) = H * C * H^t (with a vector H and a symmetric matrix C) 
 
const Scalar * const_iterator
 
GLubyte GLubyte GLubyte GLubyte w
 
CONTAINER::Scalar minimum(const CONTAINER &v)
 
A numeric matrix of compile-time fixed size. 
 
void mahalanobisDistance2AndLogPDF(const VECLIKE &diff_mean, const MATRIXLIKE &cov, T &maha2_out, T &log_pdf_out)
Computes both, the logarithm of the PDF and the square Mahalanobis distance between a point (given by...
 
Eigen::Matrix< typename MATRIX::Scalar, MATRIX::ColsAtCompileTime, MATRIX::ColsAtCompileTime > cov(const MATRIX &v)
Computes the covariance matrix from a list of samples in an NxM matrix, where each row is a sample...
 
double BASE_IMPEXP averageLogLikelihood(const CVectorDouble &logLikelihoods)
A numerically-stable method to compute average likelihood values with strongly different ranges (unwe...
 
void productIntegralAndMahalanobisTwoGaussians(const VECLIKE &mean_diffs, const MATLIKE1 &COV1, const MATLIKE2 &COV2, T &maha2_out, T &intprod_out, const MATLIKE1 *CROSS_COV12=NULL)
Computes both, the integral of the product of two Gaussians and their square Mahalanobis distance...
 
CONTAINER::Scalar maximum(const CONTAINER &v)
 
T wrapToPi(T a)
Modifies the given angle to translate it into the ]-pi,pi] range. 
 
MAT::Scalar mahalanobisDistance2(const VECTORLIKE1 &X, const VECTORLIKE2 &MU, const MAT &COV)
Computes the squared mahalanobis distance of a vector X given the mean MU and the covariance inverse ...
 
void weightedHistogram(const VECTORLIKE1 &values, const VECTORLIKE1 &weights, float binWidth, VECTORLIKE2 &out_binCenters, VECTORLIKE2 &out_binValues)
Computes the weighted histogram for a vector of values and their corresponding weights. 
 
This is the global namespace for all Mobile Robot Programming Toolkit (MRPT) libraries. 
 
#define ASSERTDEB_(f)
Defines an assertion mechanism - only when compiled in debug. 
 
VECTORLIKE1::Scalar mahalanobisDistance(const VECTORLIKE1 &X, const VECTORLIKE2 &MU, const MAT &COV)
Computes the mahalanobis distance of a vector X given the mean MU and the covariance inverse COV_inv ...
 
void covariancesAndMeanWeighted(const VECTOR_OF_VECTORS &elements, MATRIXLIKE &covariances, VECTORLIKE &means, const VECTORLIKE2 *weights_mean, const VECTORLIKE3 *weights_cov, const bool *elem_do_wrap2pi=NULL)
Computes covariances and mean of any vector of containers, given optional weights for the different s...
 
T productIntegralTwoGaussians(const std::vector< T > &mean_diffs, const CMatrixTemplateNumeric< T > &COV1, const CMatrixTemplateNumeric< T > &COV2)
Computes the integral of the product of two Gaussians, with means separated by "mean_diffs" and covar...
 
A matrix of dynamic size. 
 
int round(const T value)
Returns the closer integer (int) to x. 
 
CONTAINER::value_type element_t
 
dynamic_vector< double > CVectorDouble
Column vector, like Eigen::MatrixXd, but automatically initialized to zeros since construction...
 
void covariancesAndMean(const VECTOR_OF_VECTORS &elements, MATRIXLIKE &covariances, VECTORLIKE &means, const bool *elem_do_wrap2pi=NULL)
Computes covariances and mean of any vector of containers. 
 
void weightedHistogramLog(const VECTORLIKE1 &values, const VECTORLIKE1 &log_weights, float binWidth, VECTORLIKE2 &out_binCenters, VECTORLIKE2 &out_binValues)
Computes the weighted histogram for a vector of values and their corresponding log-weights. 
 
#define ASSERTMSG_(f, __ERROR_MSG)