Commit 6069544c authored by André Anjos's avatar André Anjos 💬

Remove all traces of Bob dependencies

parent e958429e
This diff is collapsed.
This diff is collapsed.
/**
* @date Fri Jan 27 14:10:23 2012 +0100
* @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
*
* Copyright (C) Idiap Research Institute, Martigny, Switzerland
*/
#include "det.h"
#include "linear.h"
#include "lu.h"
#include <bob.core/assert.h>
#include <limits>
double bob::math::det(const blitz::Array<double,2>& A)
{
bob::core::array::assertSameDimensionLength(A.extent(0),A.extent(1));
return bob::math::det_(A);
}
double bob::math::det_(const blitz::Array<double,2>& A)
{
// Size variable
int N = A.extent(0);
// Perform an LU decomposition
blitz::Array<double,2> L(N,N);
blitz::Array<double,2> U(N,N);
blitz::Array<double,2> P(N,N);
math::lu(A, L, U, P);
// Compute the determinant of A = det(P*L)*PI(diag(U))
// where det(P*L) = +- 1 (Number of permutation in P)
// and PI(diag(U)) is the product of the diagonal elements of U
blitz::Array<double,2> Lperm(N,N);
math::prod(P,L,Lperm);
int s = 1;
double Udiag=1.;
for (int i=0; i<N; ++i)
{
for (int j=i+1; j<N; ++j)
if (P(i,j) > 0)
{
s = -s;
break;
}
Udiag *= U(i,i);
}
return s*Udiag;
}
double bob::math::slogdet(const blitz::Array<double,2>& A, int& sign)
{
bob::core::array::assertSameDimensionLength(A.extent(0),A.extent(1));
return bob::math::slogdet_(A, sign);
}
double bob::math::slogdet_(const blitz::Array<double,2>& A, int& sign)
{
// Size variable
int N = A.extent(0);
// Perform an LU decomposition
blitz::Array<double,2> L(N,N);
blitz::Array<double,2> U(N,N);
blitz::Array<double,2> P(N,N);
math::lu(A, L, U, P);
// Compute the determinant of A = det(P*L)*SI(diag(U))
// where det(P*L) = +- 1 (Number of permutation in P)
// and SI(diag(log|U|)) is the sum of the logarithm of the
// diagonal elements of U
blitz::Array<double,2> Lperm(N,N);
math::prod(P,L,Lperm);
sign = 1;
double Udiag=0.;
for (int i=0; i<N; ++i)
{
for (int j=i+1; j<N; ++j)
if (P(i,j) > 0)
{
sign = -sign;
break;
}
Udiag += log(fabs(U(i,i)));
}
// Check for infinity
if ((Udiag*-1) == std::numeric_limits<double>::infinity())
sign = 0;
return Udiag;
}
/**
* @date Fri Jan 27 14:10:23 2012 +0100
* @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
*
* @brief This file defines function to compute the determinant of
* a 2D blitz array matrix.
*
* Copyright (C) Idiap Research Institute, Martigny, Switzerland
*/
#ifndef BOB_MATH_DET_H
#define BOB_MATH_DET_H
#include <blitz/array.h>
namespace bob { namespace math {
/**
* @brief Function which computes the determinant of a square matrix
* @param A The A matrix to consider (size NxN)
*/
double det(const blitz::Array<double,2>& A);
/**
* @brief Function which computes the determinant of a square matrix
* @param A The A matrix to consider (size NxN)
* @warning Does not check the input matrix
*/
double det_(const blitz::Array<double,2>& A);
/**
* @brief Function which computes the sign and (natural) logarithm of
* the determinant of a square matrix
* @param A The A matrix to consider (size NxN)
* @param sign The (output) sign of the determinant
* (-1 if negative, 0 if zero, +1 if positive)
*/
double slogdet(const blitz::Array<double,2>& A, int& sign);
/**
* @brief Function which computes the sign and (natural) logarithm of
* the determinant of a square matrix
* @param A The A matrix to consider (size NxN)
* @param sign The (output) sign of the determinant
* (-1 if negative, 0 if zero, +1 if positive)
* @warning Does not check the input matrix
*/
double slogdet_(const blitz::Array<double,2>& A, int& sign);
}}
#endif /* BOB_MATH_DET_H */
This diff is collapsed.
/**
* @date Mon May 16 21:45:27 2011 +0200
* @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
* @author Andre Anjos <andre.anjos@idiap.ch>
*
* @brief This file defines a function to determine an eigenvalue decomposition
* a 2D blitz array using LAPACK.
*
* Copyright (C) Idiap Research Institute, Martigny, Switzerland
*/
#ifndef BOB_MATH_EIG_H
#define BOB_MATH_EIG_H
#include <blitz/array.h>
#include <complex>
namespace bob { namespace math {
/**
* @brief Function which performs an eigenvalue decomposition of a real
* matrix, using the LAPACK function <code>dgeev</code>.
*
* @param A The A matrix to decompose (size NxN)
*
* @param V The V matrix of eigenvectors (size NxN) stored in columns
*
* @param D The vector of eigenvalues (size N) (in ascending order). Note
* that this is a 1D array rather than a 2D diagonal matrix!
*/
void eig(const blitz::Array<double,2>& A,
blitz::Array<std::complex<double>,2>& V,
blitz::Array<std::complex<double>,1>& D);
/**
* @brief Function which performs an eigenvalue decomposition of a real
* matrix, using the LAPACK function <code>dgeev</code>. This version does
* <b>NOT</b> perform any checks on the input data and should be, therefore,
* faster.
*
* @param A The A matrix to decompose (size NxN)
*
* @param V The V matrix of eigenvectors (size NxN) stored in columns
*
* @param D The vector of eigenvalues (size N) (in ascending order). Note
* that this is a 1D array rather than a 2D diagonal matrix!
*/
void eig_(const blitz::Array<double,2>& A,
blitz::Array<std::complex<double>,2>& V,
blitz::Array<std::complex<double>,1>& D);
/**
* @brief Function which performs an eigenvalue decomposition of a real
* symmetric matrix, using the LAPACK function <code>dsyevd</code>.
*
* @warning The input matrix should be symmetric.
*
* @param A The A matrix to decompose (size NxN)
*
* @param V The V matrix of eigenvectors (size NxN) stored in columns
*
* @param D The vector of eigenvalues (size N) (in ascending order). Note
* that this is a 1D array rather than a 2D diagonal matrix!
*/
void eigSym(const blitz::Array<double,2>& A, blitz::Array<double,2>& V,
blitz::Array<double,1>& D);
/**
* @brief Function which performs an eigenvalue decomposition of a real
* symmetric matrix, using the LAPACK function <code>dsyevd</code>. This
* version does <b>NOT</b> perform any checks on the input data and should
* be, therefore, faster.
*
* @warning The input matrix should be symmetric.
*
* @param A The A matrix to decompose (size NxN)
*
* @param V The V matrix of eigenvectors (size NxN) stored in columns
*
* @param D The vector of eigenvalues (size N) (in ascending order). Note
* that this is a 1D array rather than a 2D diagonal matrix!
*/
void eigSym_(const blitz::Array<double,2>& A, blitz::Array<double,2>& V,
blitz::Array<double,1>& D);
/**
* @brief Computes all the eigenvalues and the eigenvectors of a real
* generalized symmetric-definite eigenproblem, of the form:
* A*x=(lambda)*B*x, using the LAPACK function <code>dsygvd</code>.
*
* @warning The input matrices A and B are assumed to be symmetric and B is
* also positive definite.
*
* @param A The A input matrix (size NxN) of the problem
*
* @param B The B input matrix (size NxN) of the problem
*
* @param V The V matrix of eigenvectors (size NxN) stored in columns
*
* @param D The vector of eigenvalues (size N) (in ascending order). Note
* that this is a 1D array rather than a 2D diagonal matrix!
*/
void eigSym(const blitz::Array<double,2>& A, const blitz::Array<double,2>& B,
blitz::Array<double,2>& V, blitz::Array<double,1>& D);
/**
* @brief Computes all the eigenvalues and the eigenvectors of a real
* generalized symmetric-definite eigenproblem, of the form:
* A*x=(lambda)*B*x, using the LAPACK function <code>dsygvd</code>. This
* version does <b>NOT</b> perform any checks on the input data and should
* be, therefore, faster.
*
* @warning The input matrices A and B are assumed to be symmetric and B is
* also positive definite.
*
* @param A The A input matrix (size NxN) of the problem
*
* @param B The B input matrix (size NxN) of the problem
*
* @param V The V matrix of eigenvectors (size NxN) stored in columns
*
* @param D The vector of eigenvalues (size N) (in ascending order). Note
* that this is a 1D array rather than a 2D diagonal matrix!
*/
void eigSym_(const blitz::Array<double,2>& A, const blitz::Array<double,2>& B,
blitz::Array<double,2>& V, blitz::Array<double,1>& D);
}}
#endif /* BOB_MATH_EIG_H */
This diff is collapsed.
/**
* @date Mon Apr 16
* @author Manuel Guenther <Manuel.Guenther@idiap.ch>
*
* @brief Implements fast versions of some histogram measures
*
* Copyright (C) Idiap Research Institute, Martigny, Switzerland
*/
#ifndef BOB_MATH_HISTOGRAM_H
#define BOB_MATH_HISTOGRAM_H
#include <blitz/array.h>
#include <numeric>
#include <functional>
#include <cmath>
#include <bob.core/assert.h>
namespace bob{ namespace math{
template <class T>
// wrapper for the std::min function (required to be compilable under OsX)
static inline T minimum(const T& v1, const T& v2){
return std::min<T>(v1,v2);
}
template <class T>
// helper function to compute the chi_square distance between the given values
static inline T chi_square_distance(const T& v1, const T& v2){
return v1 != v2 ? (v1 - v2) * (v1 - v2) / (v1 + v2) : T(0);
}
template <class T>
// helper function to compute the symmetric Kullback-Leibler divergence
static inline double kullback_leibler_divergence(const T& v1, const T& v2){
// This implementation is inspired by http://www.informatik.uni-freiburg.de/~tipaldi/FLIRTLib/HistogramDistances_8hpp_source.html
// and http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/WAHL1/node5.html
double a1 = std::max((double)v1, 1e-5);
double a2 = std::max((double)v2, 1e-5);
return (a1 - a2) * std::log(a1/a2);
}
template <class T>
//! Fast implementation of the histogram intersection measure
inline T histogram_intersection(const blitz::Array<T,1>& h1, const blitz::Array<T,1>& h2){
bob::core::array::assertCContiguous(h1);
bob::core::array::assertCContiguous(h2);
bob::core::array::assertSameShape(h1,h2);
// we use the std::inner_product function (using blitz iterators!),
// but instead of computing the element-wise multiplication,
// we use the std::min of the two elements
return std::inner_product(
h1.begin(),
h1.end(),
h2.begin(),
T(0),
std::plus<T>(),
bob::math::minimum<T>
);
}
template <class T1, class T2>
//! Fast implementation of the sparse histogram intersection measure
inline T2 histogram_intersection(
const blitz::Array<T1,1>& index_1, const blitz::Array<T2,1>& values_1,
const blitz::Array<T1,1>& index_2, const blitz::Array<T2,1>& values_2
){
bob::core::array::assertSameShape(index_1,values_1);
bob::core::array::assertSameShape(index_2,values_2);
int i1 = 0, i2 = 0, i1_end = index_1.shape()[0], i2_end = index_2.shape()[0];
T1 p1 = index_1(i1), p2 = index_2(i2);
T2 sum = T2(0);
while (i1 < i1_end && i2 < i2_end){
p1 = index_1(i1);
p2 = index_2(i2);
if (p1 == p2) sum += bob::math::minimum(values_1(i1++), values_2(i2++));
else if (p1 < p2) ++i1;
else ++i2;
}
return sum;
}
template <class T>
//! Fast implementation of the chi square histogram distance measure
inline T chi_square(const blitz::Array<T,1>& h1, const blitz::Array<T,1>& h2){
bob::core::array::assertCContiguous(h1);
bob::core::array::assertCContiguous(h2);
bob::core::array::assertSameShape(h1,h2);
// we use the std::inner_product function (using blitz iterators!),
// but instead of computing the element-wise multiplication,
// we use our own chi_square_distance function
return std::inner_product(
h1.begin(),
h1.end(),
h2.begin(),
T(0),
std::plus<T>(),
bob::math::chi_square_distance<T>
);
}
template <class T1, class T2>
//! Fast implementation of the sparse chi square measure
inline T2 chi_square(
const blitz::Array<T1,1>& index_1, const blitz::Array<T2,1>& values_1,
const blitz::Array<T1,1>& index_2, const blitz::Array<T2,1>& values_2
){
bob::core::array::assertSameShape(index_1,values_1);
bob::core::array::assertSameShape(index_2,values_2);
int i1 = 0, i2 = 0, i1_end = index_1.shape()[0], i2_end = index_2.shape()[0];
T1 p1 = index_1(i1), p2 = index_2(i2);
T2 sum = T2(0);
while (i1 < i1_end && i2 < i2_end){
p1 = index_1(i1);
p2 = index_2(i2);
if (p1 == p2){
sum += bob::math::chi_square_distance(values_1(i1++), values_2(i2++));
} else if (p1 < p2) {
sum += bob::math::chi_square_distance(values_1(i1++), T2(0));
} else{
sum += bob::math::chi_square_distance(T2(0), values_2(i2++));
}
}
// roll up the remaining parts of the histograms
while (i1 < i1_end) sum += chi_square_distance(values_1(i1++), T2(0));
while (i2 < i2_end) sum += chi_square_distance(T2(0), values_2(i2++));
return sum;
}
template <class T>
//! Fast implementation of the symmetric Kullback-Leibler histogram divergence measure (a distance measure)
inline double kullback_leibler(const blitz::Array<T,1>& h1, const blitz::Array<T,1>& h2){
bob::core::array::assertCContiguous(h1);
bob::core::array::assertCContiguous(h2);
bob::core::array::assertSameShape(h1,h2);
// we use the std::inner_product function (using blitz iterators!)
// calling out implementation of the Kullback-Leibler divergence of two values
return std::inner_product(
h1.begin(),
h1.end(),
h2.begin(),
0.,
std::plus<T>(),
bob::math::kullback_leibler_divergence<T>
);
}
template <class T1, class T2>
//! Fast implementation of the sparse symmetric Kullback-Leibler histogram divergence measure (a distance measure)
inline double kullback_leibler(
const blitz::Array<T1,1>& index_1, const blitz::Array<T2,1>& values_1,
const blitz::Array<T1,1>& index_2, const blitz::Array<T2,1>& values_2
){
bob::core::array::assertSameShape(index_1,values_1);
bob::core::array::assertSameShape(index_2,values_2);
int i1 = 0, i2 = 0, i1_end = index_1.shape()[0], i2_end = index_2.shape()[0];
T1 p1 = index_1(i1), p2 = index_2(i2);
double sum = 0.;
while (i1 < i1_end && i2 < i2_end){
p1 = index_1(i1);
p2 = index_2(i2);
if (p1 == p2){
sum += bob::math::kullback_leibler_divergence(values_1(i1++), values_2(i2++));
} else if (p1 < p2) {
sum += bob::math::kullback_leibler_divergence(values_1(i1++), T2(0));
} else{
sum += bob::math::kullback_leibler_divergence(T2(0), values_2(i2++));
}
}
// roll up the remaining parts of the histograms
while (i1 < i1_end) sum += bob::math::kullback_leibler_divergence(values_1(i1++), T2(0));
while (i2 < i2_end) sum += bob::math::kullback_leibler_divergence(T2(0), values_2(i2++));
return sum;
}
}}
#endif // BOB_MATH_HISTOGRAM_H
/**
* @date Fri Jan 27 14:10:23 2012 +0100
* @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
*
* Copyright (C) Idiap Research Institute, Martigny, Switzerland
*/
#include <stdexcept>
#include <boost/shared_array.hpp>
#include "inv.h"
#include "linear.h"
#include <bob.core/assert.h>
#include <bob.core/check.h>
#include <bob.core/array_copy.h>
// Declaration of the external LAPACK function
// LU decomposition of a general matrix (dgetrf)
extern "C" void dgetrf_( const int *M, const int *N, double *A, const int *lda,
int *ipiv, int *info);
// Inverse of a general matrix (dgetri)
extern "C" void dgetri_( const int *N, double *A, const int *lda,
const int *ipiv, double *work, const int *lwork, int *info);
void bob::math::inv(const blitz::Array<double,2>& A, blitz::Array<double,2>& B)
{
// Size variable
const int N = A.extent(0);
const blitz::TinyVector<int,2> shapeA(N,N);
bob::core::array::assertZeroBase(A);
bob::core::array::assertZeroBase(B);
bob::core::array::assertSameShape(A,shapeA);
bob::core::array::assertSameShape(B,shapeA);
bob::math::inv_(A, B);
}
void bob::math::inv_(const blitz::Array<double,2>& A, blitz::Array<double,2>& B)
{
// Size variable
const int N = A.extent(0);
//////////////////////////////////////
// Prepares to call LAPACK functions
// Initializes LAPACK variables
int info = 0;
const int lda = N;
// Initializes LAPACK arrays
boost::shared_array<int> ipiv(new int[N]);
// Tries to use B directly if possible
// Input and output arrays are both column-major order.
// Hence, we can ignore the problem of column- and row-major order
// conversions.
bool B_direct_use = bob::core::array::isCZeroBaseContiguous(B);
blitz::Array<double,2> A_blitz_lapack;
if (B_direct_use)
{
A_blitz_lapack.reference(B);
A_blitz_lapack = A;
}
else
A_blitz_lapack.reference(bob::core::array::ccopy(A));
double *A_lapack = A_blitz_lapack.data();
// Calls the LAPACK functions
// 1/ Computes the LU decomposition
dgetrf_( &N, &N, A_lapack, &lda, ipiv.get(), &info);
// Checks the info variable
if (info != 0)
throw std::runtime_error("The LAPACK dgetrf function returned a non-zero value.");
// TODO: We might consider adding a real invertibility test as described in
// this thread (Btw, this is what matlab does):
// http://icl.cs.utk.edu/lapack-forum/archives/lapack/msg00778.html
// 2/ Computes the inverse matrix
// 2/A/ Queries the optimal size of the working array
const int lwork_query = -1;
double work_query;
dgetri_( &N, A_lapack, &lda, ipiv.get(), &work_query, &lwork_query, &info);
// 2/B/ Computes the inverse
const int lwork = static_cast<int>(work_query);
boost::shared_array<double> work(new double[lwork]);
dgetri_( &N, A_lapack, &lda, ipiv.get(), work.get(), &lwork, &info);
// Checks info variable
if (info != 0)
throw std::runtime_error("The LAPACK dgetri function returned a non-zero value. The matrix might not be invertible.");
// Copy back content to B if required
if (!B_direct_use)
B = A_blitz_lapack;
}
/**
* @date Fri Jan 27 14:10:23 2012 +0100
* @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
*
* @brief This file defines function to inverse a 2D blitz array matrix.
*
* Copyright (C) Idiap Research Institute, Martigny, Switzerland
*/
#ifndef BOB_MATH_INV_H
#define BOB_MATH_INV_H
#include <blitz/array.h>
namespace bob { namespace math {
/**
* @brief Function which computes the inverse of a matrix,
* using the dgetrf and dgetri LAPACK functions.
* @param A The A matrix to decompose (size NxN)
* @param B The B=inverse(A) matrix (size NxN)
*/
void inv(const blitz::Array<double,2>& A, blitz::Array<double,2>& B);
void inv_(const blitz::Array<double,2>& A, blitz::Array<double,2>& B);
}}
#endif /* BOB_MATH_INV_H */
This diff is collapsed.
This diff is collapsed.
/**
* @date Sat Mar 19 19:49:51 2011 +0100
* @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
*
* @brief This file defines functions to solve linear systems
* A*x=b using LAPACK or a conjugate gradient implementation.
*
* Copyright (C) Idiap Research Institute, Martigny, Switzerland
*/
#ifndef BOB_MATH_LINSOLVE_H
#define BOB_MATH_LINSOLVE_H
#include <blitz/array.h>
namespace bob { namespace math {
/**
* @brief Function which solves a linear system of equation using the
* 'generic' dgsev LAPACK function.
* @param A The A squared-matrix of the system A*x=b (size NxN)
* @param b The b vector of the system A*x=b (size N)
* @param x The x vector of the system A*x=b which will be updated
* at the end of the function.
*/
void linsolve(const blitz::Array<double,2>& A, blitz::Array<double,1>& x,
const blitz::Array<double,1>& b);
void linsolve_(const blitz::Array<double,2>& A, blitz::Array<double,1>& x,
const blitz::Array<double,1>& b);
/**
* @brief Function which solves a linear system of equation using the
* 'generic' dgsev LAPACK function.
* @param A The A squared-matrix of the system A*X=B (size NxN)
* @param B The B matrix of the system A*X=B (size NxP)
* @param X The X matrix of the system A*X=B which will be updated
* at the end of the function (size NxP).
*/
void linsolve(const blitz::Array<double,2>& A, blitz::Array<double,2>& X,
const blitz::Array<double,2>& B);
void linsolve_(const blitz::Array<double,2>& A, blitz::Array<double,2>& X,
const blitz::Array<double,2>& B);
/**
* @brief Function which solves a symmetric positive definite linear
* system of equation using the dposv LAPACK function.
* @warning No check is performed wrt. to the fact that A should be
* symmetric positive definite.
* @param A The A squared-matrix, symmetric definite positive, of the
* system A*x=b (size NxN)
* @param b The b vector of the system A*x=b (size N)
* @param x The x vector of the system A*x=b which will be updated
* at the end of the function (size N)
*/
void linsolveSympos(const blitz::Array<double,2>& A,
blitz::Array<double,1>& x, const blitz::Array<double,1>& b);
void linsolveSympos_(const blitz::Array<double,2>& A,
blitz::Array<double,1>& x, const blitz::Array<double,1>& b);
/**
* @brief Function which solves a symmetric positive definite linear
* system of equation using the dposv LAPACK function.
* @warning No check is performed wrt. to the fact that A should be
* symmetric positive definite.
* @param A The A squared-matrix, symmetric definite positive, of the
* system A*X=B (size NxN)
* @param B The B vector of the system A*x=b (size NxP)
* @param X The X vector of the system A*x=b which will be updated
* at the end of the function (size NxP)
*/
void linsolveSympos(const blitz::Array<double,2>& A,
blitz::Array<double,2>& X, const blitz::Array<double,2>& B);
void linsolveSympos_(const blitz::Array<double,2>& A,
blitz::Array<double,2>& X, const blitz::Array<double,2>& B);
/**
* @brief Function which solves a symmetric positive-definite linear
* system of equation via conjugate gradients.
* @param A The A symmetric positive-definite squared-matrix of the
* system A*x=b (size NxN)
* @param b The b vector of the system A*x=b (size N)
* @param x The x vector of the system A*x=b which will be updated
* at the end of the function.
* @param acc The desired accuracy. The algorithm terminates when
* norm(Ax-b)/norm(b) < acc
* @param max_iter The maximum number of iterations
*/
void linsolveCGSympos(const blitz::Array<double,2>& A, blitz::Array<double,1>& x,
const blitz::Array<double,1>& b, const double acc, const int max_iter);
void linsolveCGSympos_(const blitz::Array<double,2>& A, blitz::Array<double,1>& x,
const blitz::Array<double,1>& b, const double acc, const int max_iter);
}}
#endif /* BOB_MATH_LINSOLVE_H */
/**
* @date Fri Feb 10 20:02:07 2012 +0200
* @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
*
* Copyright (C) Idiap Research Institute, Martigny, Switzerland
*/
#include <stdexcept>
#include <boost/format.hpp>