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)