Binded the Generalized SVD from LAPACK (issue #6) swaping properly its output (issue #7)

parent 5010cd86
Pipeline #6541 passed with stages
in 4 minutes and 17 seconds
This diff is collapsed.
/**
* @date Tue Jan 10 21:14 2016 +0100
* @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
*
* @brief Binds the Generalized SVD
*/
#include "gsvd.h"
#include <bob.blitz/cppapi.h>
#include <bob.blitz/cleanup.h>
#include <bob.math/gsvd.h>
PyObject* py_gsvd (PyObject*, PyObject* args, PyObject* kwds) {
/* Parses input arguments in a single shot */
static const char* const_kwlist[] = { "A", "B", 0 /* Sentinel */ };
static char** kwlist = const_cast<char**>(const_kwlist);
PyBlitzArrayObject* A = 0;
PyBlitzArrayObject* B = 0;
if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&O&", kwlist,
&PyBlitzArray_Converter, &A,
&PyBlitzArray_Converter, &B
))
return 0;
auto A_ = make_safe(A);
auto B_ = make_safe(B);
if (A->ndim != 2 || A->type_num != NPY_FLOAT64) {
PyErr_Format(PyExc_TypeError, "`A` matrix only supports 2D 64-bit float array");
return 0;
}
if (B->ndim != 2 || B->type_num != NPY_FLOAT64) {
PyErr_Format(PyExc_TypeError, "`B` matrix only supports 2D 64-bit float array");
return 0;
}
auto A_bz = PyBlitzArrayCxx_AsBlitz<double,2>(A);
auto B_bz = PyBlitzArrayCxx_AsBlitz<double,2>(B);
const int M = A_bz->extent(0);
const int N = A_bz->extent(1);
const int P = B_bz->extent(0);
// Creating the output matrices
blitz::Array<double,2> U(M, M); U=0;
blitz::Array<double,2> V(P, P); V=0;
blitz::Array<double,2> Q(N, N); Q=0;
blitz::Array<double,2> zeroR(N, N); zeroR=0;
blitz::Array<double,2> X(N, N); X=0;
blitz::Array<double,2> C(N, N); C=0;
blitz::Array<double,2> S(N, N); S=0;
try {
bob::math::gsvd(*A_bz,*B_bz,U,V,zeroR,Q,X,C,S);
return Py_BuildValue("NNNNN",
PyBlitzArrayCxx_AsConstNumpy(U),
PyBlitzArrayCxx_AsConstNumpy(V),
PyBlitzArrayCxx_AsConstNumpy(X),
PyBlitzArrayCxx_AsConstNumpy(C),
PyBlitzArrayCxx_AsConstNumpy(S));
}
catch (std::exception& e) {
PyErr_SetString(PyExc_RuntimeError, e.what());
}
catch (...) {
PyErr_SetString(PyExc_RuntimeError, "norminv failed: unknown exception caught");
}
return 0;
}
/**
* @author Andre Anjos <andre.anjos@idiap.ch>
* @date Thu 5 Dec 12:01:57 2013
*
* @brief Declaration of components relevant for main.cpp
*/
#include <Python.h>
PyObject* py_gsvd(PyObject*, PyObject* args, PyObject* kwds);
/**
* @date Tue Jan 10 21:14 2016 +0100
* @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
*
* @brief This file defines a function to determine the SVD decomposition
* a 2D blitz array using LAPACK.
*
* Copyright (C) Idiap Research Institute, Martigny, Switzerland
*/
#ifndef BOB_MATH_GSVD_H
#define BOB_MATH_GSVD_H
#include <blitz/array.h>
namespace bob { namespace math {
/**
* @ingroup MATH
* @{
*/
/**
* @brief Function which performs the Generalized Singular Value Decomposition
* using the 'simple' driver routine dggsvd3 of LAPACK.
* Just remembering that given A and B as input we have:
*
* A = U C X and B = V S X
* where X = [0,R]Q^T
*
*
*
* @warning The output blitz::array sigma should have the correct
* size, with zero base index. Checks are performed.
* @param A The A matrix to decompose (size MxP)
* @param B The B matrix to decompose (size NxP)
* @warning A and B must have the same number of columns, but may have different numbers of rows. If A is m-by-p and B is n-by-p, then U is m-by-m, V is n-by-n, X is p-by-q, C is m-by-q
* and S is n-by-q, where q = min(m+n,p).
* @param U The U matrix (MxM)
* @param V The V matrix (NxN)
* @param zeroR The X matrix (rxN)
* @param Q The Q matrix (Q)
* @param C The C matrix (MxQ)
* @param S The S matrix (NxQ)
*/
void gsvd(blitz::Array<double,2>& A,
blitz::Array<double,2>& B,
blitz::Array<double,2>& U,
blitz::Array<double,2>& V,
blitz::Array<double,2>& zeroR,
blitz::Array<double,2>& Q,
blitz::Array<double,2>& X,
blitz::Array<double,2>& C,
blitz::Array<double,2>& S
);
/**
* @brief Swaping using the LAPACK variable iWork
* http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga4a187519e5c71da3b3f67c85e9baf0f2.html#ga4a187519e5c71da3b3f67c85e9baf0f2
*
* @param A 1D Matrix
* @param indexes Sorting information
* @param n number of operations
*/
template<typename T>
void swap_(blitz::Array<T,1>& A, int* indexes, int begin, int end) {
T aux = 0;
int fortran_index = 0;
for (int i=0; i<A.extent(0); i++){
fortran_index = indexes[i]-1;
aux = A(i);
A(i) = A(fortran_index);
A(fortran_index) = aux;
}
}
/**
* @brief Swaping using the LAPACK variable iWork
* http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga4a187519e5c71da3b3f67c85e9baf0f2.html#ga4a187519e5c71da3b3f67c85e9baf0f2
*
* @param A D Matrix
* @param indexes Sorting information
* @param n number of operations
*/
template<typename T>
void swap_(blitz::Array<T,2>& A, int* indexes, int begin, int end) {
blitz::Array<T,1> aux(A.extent(1)); aux = 0;
int fortran_index = 0;
for (int i=begin; i<end; i++){
fortran_index = indexes[i]-1;
aux = A(blitz::Range::all(), i);
A(blitz::Range::all(), i) = A(blitz::Range::all(), fortran_index);
A(blitz::Range::all(), fortran_index) = aux;
}
}
/**
* @}
*/
}}
#endif /* BOB_MATH_GSVD_H */
......@@ -20,6 +20,7 @@
#include "norminv.h"
#include "scatter.h"
#include "lp_interior_point.h"
#include "gsvd.h"
static bob::extension::FunctionDoc s_histogram_intersection = bob::extension::FunctionDoc(
"histogram_intersection",
......@@ -355,6 +356,26 @@ static bob::extension::FunctionDoc s_scatters_nocheck = bob::extension::Function
;
static bob::extension::FunctionDoc s_gsvd = bob::extension::FunctionDoc(
"gsvd",
"Computes the Generalized SVD",
"Computes the Generalized SVD. The output of this function is similar with the one found in Matlab\n"
"[U,V,X,C,S] = gsvd(A,B) returns unitary matrices :math:`U` and :math:`V`, the square matrix :math:`X` (which is :math:`[0 R] Q^{T}`), and nonnegative diagonal matrices :math:`C` and :math:`S` such that:\n\n"
".. math:: C^{T}C + S^{T}S = I \n"
".. math:: A = (XC^{T}U^{T})^{T}\n"
".. math:: B = (XS^{T}V^{T})^{T}\n"
)
.add_prototype("A, B", "")
.add_parameter("A", "[array_like (float, 2D)]", "Must be :math:`m \\times n`")
.add_parameter("B", "[array_like (float, 2D)]", "Must be :math:`p \\times n`")
.add_return("U", "[array_like (float, 2D)]", "Contains a :math:`m \\times m` orthogonal matrix.")
.add_return("V", "[array_like (float, 2D)]", "Contains a :math:`n \\times n` orthogonal matrix.")
.add_return("X", "[array_like (float, 2D)]", "Contains a :math:`p \\times q` matrix, where :math:`p=\\min(m+n,p)` and :math:`X=[0, R] Q^{T}` (Check LAPACK documentation).")
.add_return("C", "[array_like (float, 2D)]", "Contains a :math:`p \\times q` matrix.")
.add_return("S", "[array_like (float, 2D)]", "Contains a :math:`n \\times q` matrix.")
;
static PyMethodDef module_methods[] = {
{
s_histogram_intersection.name(),
......@@ -470,6 +491,13 @@ static PyMethodDef module_methods[] = {
METH_VARARGS|METH_KEYWORDS,
s_scatters_nocheck.doc()
},
{
s_gsvd.name(),
(PyCFunction)py_gsvd,
METH_VARARGS|METH_KEYWORDS,
s_gsvd.doc()
},
{0} /* Sentinel */
};
......
#!/usr/bin/env python
# Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
# Sun Jan 15 19:12:43 CET 2017
#
"""
Tests GSVD
Basically these tests test the GSVD relation.
Given 2 matrices A and B GSVD(A,B) = [U,V,X,C,S] where,
A= (X * C.T * U^T)^T and
B= (X * S.T * V^T)^T and
C**2 + S**2 = 1
"""
import bob.math
import numpy
import nose.tools
numpy.random.seed(10)
def gsvd_relations(A,B):
[U,V,X,C,S] = bob.math.gsvd(A, B)
# Cheking the relation C**2 + S**2 = 1
I = numpy.eye(A.shape[1])
I_check = numpy.dot(C.T, C) + numpy.dot(S.T, S)
nose.tools.eq_( (abs(I-I_check) < 1e-10).all(), True )
# Cheking the relation A= (X * C.T * U^T)^T
A_check = numpy.dot(numpy.dot(X,C.T),U.T).T
nose.tools.eq_( (abs(A-A_check) < 1e-10).all(), True )
# Cheking the relation B= (X * S.T * V^T)^T
B_check = numpy.dot(numpy.dot(X,S.T),V.T).T
nose.tools.eq_( (abs(B-B_check) < 1e-10).all(), True )
del U,V,X,C,S
def test_first_case():
"""
Testing the first scenario:
M-K-L >= 0 (check http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga4a187519e5c71da3b3f67c85e9baf0f2.html#ga4a187519e5c71da3b3f67c85e9baf0f2)
"""
A = numpy.random.rand(10,10)
B = numpy.random.rand(790,10)
gsvd_relations(A, B)
def test_second_case():
"""
Testing the second scenario:
M-K-L < 0 (check http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga4a187519e5c71da3b3f67c85e9baf0f2.html#ga4a187519e5c71da3b3f67c85e9baf0f2)
"""
A = numpy.random.rand(4,5)
B = numpy.random.rand(11,5)
gsvd_relations(A, B)
......@@ -15,7 +15,19 @@ Summary
bob.math.LPInteriorPoint
bob.math.LPInteriorPointShortstep
bob.math.LPInteriorPointLongstep
bob.math.chi_square
bob.math.gsvd
bob.math.histogram_intersection
bob.math.kullback_leibler
bob.math.linsolve
bob.math.linsolve_cg_sympos
bob.math.linsolve_sympos
bob.math.norminv
bob.math.pavx
bob.math.pavxWidth
bob.math.pavxWidthHeight
bob.math.scatter
bob.math.scatters
Details
.......
......
......@@ -179,6 +179,7 @@ setup(
"bob/math/cpp/pavx.cpp",
"bob/math/cpp/pinv.cpp",
"bob/math/cpp/svd.cpp",
"bob/math/cpp/gsvd.cpp",
"bob/math/cpp/sqrtm.cpp",
],
version = version,
......@@ -195,6 +196,7 @@ setup(
"bob/math/linsolve.cpp",
"bob/math/pavx.cpp",
"bob/math/norminv.cpp",
"bob/math/gsvd.cpp",
"bob/math/scatter.cpp",
"bob/math/lp_interior_point.cpp",
"bob/math/main.cpp",
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment