From b4bace812337c9704febdfdf394c6c1000df9bc8 Mon Sep 17 00:00:00 2001
From: Tiago Freitas Pereira <tiagofrepereira@gmail.com>
Date: Wed, 28 Jan 2015 11:10:19 +0100
Subject: [PATCH] Split the Factor analysis stuff in different files. Change
 the name of some method in order to be more coherent

---
 bob/learn/misc/cpp/FABase.cpp                 | 307 ++++++++++++++++++
 bob/learn/misc/cpp/ISVBase.cpp                |  76 +++++
 bob/learn/misc/cpp/ISVMachine.cpp             | 184 +++++++++++
 bob/learn/misc/cpp/JFABase.cpp                |  75 +++++
 .../misc/include/bob.learn.misc/FABase.h      | 293 +++++++++++++++++
 .../misc/include/bob.learn.misc/ISVBase.h     | 228 +++++++++++++
 .../misc/include/bob.learn.misc/ISVMachine.h  | 230 +++++++++++++
 .../misc/include/bob.learn.misc/JFABase.h     | 253 +++++++++++++++
 8 files changed, 1646 insertions(+)
 create mode 100644 bob/learn/misc/cpp/FABase.cpp
 create mode 100644 bob/learn/misc/cpp/ISVBase.cpp
 create mode 100644 bob/learn/misc/cpp/ISVMachine.cpp
 create mode 100644 bob/learn/misc/cpp/JFABase.cpp
 create mode 100644 bob/learn/misc/include/bob.learn.misc/FABase.h
 create mode 100644 bob/learn/misc/include/bob.learn.misc/ISVBase.h
 create mode 100644 bob/learn/misc/include/bob.learn.misc/ISVMachine.h
 create mode 100644 bob/learn/misc/include/bob.learn.misc/JFABase.h

diff --git a/bob/learn/misc/cpp/FABase.cpp b/bob/learn/misc/cpp/FABase.cpp
new file mode 100644
index 0000000..dab2aeb
--- /dev/null
+++ b/bob/learn/misc/cpp/FABase.cpp
@@ -0,0 +1,307 @@
+/**
+ * @date Tue Jan 27 15:51:15 2015 +0200
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+
+#include <bob.learn.misc/FABase.h>
+#include <bob.core/array_copy.h>
+#include <bob.math/linear.h>
+#include <bob.math/inv.h>
+#include <limits>
+
+
+//////////////////// FABase ////////////////////
+bob::learn::misc::FABase::FABase():
+  m_ubm(boost::shared_ptr<bob::learn::misc::GMMMachine>()), m_ru(1), m_rv(1),
+  m_U(0,1), m_V(0,1), m_d(0)
+{}
+
+bob::learn::misc::FABase::FABase(const boost::shared_ptr<bob::learn::misc::GMMMachine> ubm,
+    const size_t ru, const size_t rv):
+  m_ubm(ubm), m_ru(ru), m_rv(rv),
+  m_U(getSupervectorLength(),ru), m_V(getSupervectorLength(),rv), m_d(getSupervectorLength())
+{
+  if (ru < 1) {
+    boost::format m("value for parameter `ru' (%lu) cannot be smaller than 1");
+    m % ru;
+    throw std::runtime_error(m.str());
+  }
+  if (rv < 1) {
+    boost::format m("value for parameter `rv' (%lu) cannot be smaller than 1");
+    m % ru;
+    throw std::runtime_error(m.str());
+  }
+  updateCache();
+}
+
+bob::learn::misc::FABase::FABase(const bob::learn::misc::FABase& other):
+  m_ubm(other.m_ubm), m_ru(other.m_ru), m_rv(other.m_rv),
+  m_U(bob::core::array::ccopy(other.m_U)),
+  m_V(bob::core::array::ccopy(other.m_V)),
+  m_d(bob::core::array::ccopy(other.m_d))
+{
+  updateCache();
+}
+
+bob::learn::misc::FABase::~FABase() {
+}
+
+bob::learn::misc::FABase& bob::learn::misc::FABase::operator=
+(const bob::learn::misc::FABase& other)
+{
+  if (this != &other)
+  {
+    m_ubm = other.m_ubm;
+    m_ru = other.m_ru;
+    m_rv = other.m_rv;
+    m_U.reference(bob::core::array::ccopy(other.m_U));
+    m_V.reference(bob::core::array::ccopy(other.m_V));
+    m_d.reference(bob::core::array::ccopy(other.m_d));
+
+    updateCache();
+  }
+  return *this;
+}
+
+bool bob::learn::misc::FABase::operator==(const bob::learn::misc::FABase& b) const
+{
+  return ( (((m_ubm && b.m_ubm) && *m_ubm == *(b.m_ubm)) || (!m_ubm && !b.m_ubm)) &&
+          m_ru == b.m_ru && m_rv == b.m_rv &&
+          bob::core::array::isEqual(m_U, b.m_U) &&
+          bob::core::array::isEqual(m_V, b.m_V) &&
+          bob::core::array::isEqual(m_d, b.m_d));
+}
+
+bool bob::learn::misc::FABase::operator!=(const bob::learn::misc::FABase& b) const
+{
+  return !(this->operator==(b));
+}
+
+bool bob::learn::misc::FABase::is_similar_to(const bob::learn::misc::FABase& b,
+    const double r_epsilon, const double a_epsilon) const
+{
+  // TODO: update is_similar_to of the GMMMachine with the 2 epsilon's
+  return (( ((m_ubm && b.m_ubm) && m_ubm->is_similar_to(*(b.m_ubm), a_epsilon)) ||
+            (!m_ubm && !b.m_ubm) ) &&
+          m_ru == b.m_ru && m_rv == b.m_rv &&
+          bob::core::array::isClose(m_U, b.m_U, r_epsilon, a_epsilon) &&
+          bob::core::array::isClose(m_V, b.m_V, r_epsilon, a_epsilon) &&
+          bob::core::array::isClose(m_d, b.m_d, r_epsilon, a_epsilon));
+}
+
+void bob::learn::misc::FABase::resize(const size_t ru, const size_t rv)
+{
+  if (ru < 1) {
+    boost::format m("value for parameter `ru' (%lu) cannot be smaller than 1");
+    m % ru;
+    throw std::runtime_error(m.str());
+  }
+  if (rv < 1) {
+    boost::format m("value for parameter `rv' (%lu) cannot be smaller than 1");
+    m % ru;
+    throw std::runtime_error(m.str());
+  }
+
+  m_ru = ru;
+  m_rv = rv;
+  m_U.resizeAndPreserve(m_U.extent(0), ru);
+  m_V.resizeAndPreserve(m_V.extent(0), rv);
+
+  updateCacheUbmUVD();
+}
+
+void bob::learn::misc::FABase::resize(const size_t ru, const size_t rv, const size_t cd)
+{
+  if (ru < 1) {
+    boost::format m("value for parameter `ru' (%lu) cannot be smaller than 1");
+    m % ru;
+    throw std::runtime_error(m.str());
+  }
+  if (rv < 1) {
+    boost::format m("value for parameter `rv' (%lu) cannot be smaller than 1");
+    m % ru;
+    throw std::runtime_error(m.str());
+  }
+
+  if (!m_ubm || (m_ubm && getSupervectorLength() == cd))
+  {
+    m_ru = ru;
+    m_rv = rv;
+    m_U.resizeAndPreserve(cd, ru);
+    m_V.resizeAndPreserve(cd, rv);
+    m_d.resizeAndPreserve(cd);
+
+    updateCacheUbmUVD();
+  }
+  else {
+    boost::format m("value for parameter `cd' (%lu) is not set to %lu");
+    m % cd % getSupervectorLength();
+    throw std::runtime_error(m.str());
+  }
+}
+
+void bob::learn::misc::FABase::setUbm(const boost::shared_ptr<bob::learn::misc::GMMMachine> ubm)
+{
+  m_ubm = ubm;
+  m_U.resizeAndPreserve(getSupervectorLength(), m_ru);
+  m_V.resizeAndPreserve(getSupervectorLength(), m_rv);
+  m_d.resizeAndPreserve(getSupervectorLength());
+
+  updateCache();
+}
+
+void bob::learn::misc::FABase::setU(const blitz::Array<double,2>& U)
+{
+  if(U.extent(0) != m_U.extent(0)) { //checks dimension
+    boost::format m("number of rows in parameter `U' (%d) does not match the expected size (%d)");
+    m % U.extent(0) % m_U.extent(0);
+    throw std::runtime_error(m.str());
+  }
+  if(U.extent(1) != m_U.extent(1)) { //checks dimension
+    boost::format m("number of columns in parameter `U' (%d) does not match the expected size (%d)");
+    m % U.extent(1) % m_U.extent(1);
+    throw std::runtime_error(m.str());
+  }
+  m_U.reference(bob::core::array::ccopy(U));
+
+  // update cache
+  updateCacheUbmUVD();
+}
+
+void bob::learn::misc::FABase::setV(const blitz::Array<double,2>& V)
+{
+  if(V.extent(0) != m_V.extent(0)) { //checks dimension
+    boost::format m("number of rows in parameter `V' (%d) does not match the expected size (%d)");
+    m % V.extent(0) % m_V.extent(0);
+    throw std::runtime_error(m.str());
+  }
+  if(V.extent(1) != m_V.extent(1)) { //checks dimension
+    boost::format m("number of columns in parameter `V' (%d) does not match the expected size (%d)");
+    m % V.extent(1) % m_V.extent(1);
+    throw std::runtime_error(m.str());
+  }
+  m_V.reference(bob::core::array::ccopy(V));
+}
+
+void bob::learn::misc::FABase::setD(const blitz::Array<double,1>& d)
+{
+  if(d.extent(0) != m_d.extent(0)) { //checks dimension
+    boost::format m("size of input vector `d' (%d) does not match the expected size (%d)");
+    m % d.extent(0) % m_d.extent(0);
+    throw std::runtime_error(m.str());
+  }
+  m_d.reference(bob::core::array::ccopy(d));
+}
+
+
+void bob::learn::misc::FABase::updateCache()
+{
+  updateCacheUbm();
+  updateCacheUbmUVD();
+  resizeTmp();
+}
+
+void bob::learn::misc::FABase::resizeTmp()
+{
+  m_tmp_IdPlusUSProdInv.resize(getDimRu(),getDimRu());
+  m_tmp_Fn_x.resize(getSupervectorLength());
+  m_tmp_ru.resize(getDimRu());
+  m_tmp_ruD.resize(getDimRu(), getNInputs());
+  m_tmp_ruru.resize(getDimRu(), getDimRu());
+}
+
+void bob::learn::misc::FABase::updateCacheUbm()
+{
+  // Put supervectors in cache
+  if (m_ubm)
+  {
+    m_cache_mean.resize(getSupervectorLength());
+    m_cache_sigma.resize(getSupervectorLength());
+    m_cache_mean  = m_ubm->getMeanSupervector();
+    m_cache_sigma = m_ubm->getVarianceSupervector();
+  }
+}
+
+void bob::learn::misc::FABase::updateCacheUbmUVD()
+{
+  // Compute and put  U^{T}.diag(sigma)^{-1} in cache
+  if (m_ubm)
+  {
+    blitz::firstIndex i;
+    blitz::secondIndex j;
+    m_cache_UtSigmaInv.resize(getDimRu(), getSupervectorLength());
+    m_cache_UtSigmaInv = m_U(j,i) / m_cache_sigma(j); // Ut * diag(sigma)^-1
+  }
+}
+
+void bob::learn::misc::FABase::computeIdPlusUSProdInv(const bob::learn::misc::GMMStats& gmm_stats,
+  blitz::Array<double,2>& output) const
+{
+  // Computes (Id + U^T.Sigma^-1.U.N_{i,h}.U)^-1 =
+  // (Id + sum_{c=1..C} N_{i,h}.U_{c}^T.Sigma_{c}^-1.U_{c})^-1
+
+  // Blitz compatibility: ugly fix (const_cast, as old blitz version does not
+  // provide a non-const version of transpose())
+  blitz::Array<double,2> Ut = const_cast<blitz::Array<double,2>&>(m_U).transpose(1,0);
+
+  blitz::firstIndex i;
+  blitz::secondIndex j;
+  blitz::Range rall = blitz::Range::all();
+
+  bob::math::eye(m_tmp_ruru); // m_tmp_ruru = Id
+  // Loop and add N_{i,h}.U_{c}^T.Sigma_{c}^-1.U_{c} to m_tmp_ruru at each iteration
+  const size_t dim_c = getNGaussians();
+  const size_t dim_d = getNInputs();
+  for(size_t c=0; c<dim_c; ++c) {
+    blitz::Range rc(c*dim_d,(c+1)*dim_d-1);
+    blitz::Array<double,2> Ut_c = Ut(rall,rc);
+    blitz::Array<double,1> sigma_c = m_cache_sigma(rc);
+    m_tmp_ruD = Ut_c(i,j) / sigma_c(j); // U_{c}^T.Sigma_{c}^-1
+    blitz::Array<double,2> U_c = m_U(rc,rall);
+    // Use m_cache_IdPlusUSProdInv as an intermediate array
+    bob::math::prod(m_tmp_ruD, U_c, output); // U_{c}^T.Sigma_{c}^-1.U_{c}
+    // Finally, add N_{i,h}.U_{c}^T.Sigma_{c}^-1.U_{c} to m_tmp_ruru
+    m_tmp_ruru += output * gmm_stats.n(c);
+  }
+  // Computes the inverse
+  bob::math::inv(m_tmp_ruru, output);
+}
+
+
+void bob::learn::misc::FABase::computeFn_x(const bob::learn::misc::GMMStats& gmm_stats,
+  blitz::Array<double,1>& output) const
+{
+  // Compute Fn_x = sum_{sessions h}(N*(o - m) (Normalised first order statistics)
+  blitz::Range rall = blitz::Range::all();
+  const size_t dim_c = getNGaussians();
+  const size_t dim_d = getNInputs();
+  for(size_t c=0; c<dim_c; ++c) {
+    blitz::Range rc(c*dim_d,(c+1)*dim_d-1);
+    blitz::Array<double,1> Fn_x_c = output(rc);
+    blitz::Array<double,1> mean_c = m_cache_mean(rc);
+    Fn_x_c = gmm_stats.sumPx(c,rall) - mean_c*gmm_stats.n(c);
+  }
+}
+
+void bob::learn::misc::FABase::estimateX(const blitz::Array<double,2>& IdPlusUSProdInv,
+  const blitz::Array<double,1>& Fn_x, blitz::Array<double,1>& x) const
+{
+  // m_tmp_ru = UtSigmaInv * Fn_x = Ut*diag(sigma)^-1 * N*(o - m)
+  bob::math::prod(m_cache_UtSigmaInv, Fn_x, m_tmp_ru);
+  // x = IdPlusUSProdInv * m_cache_UtSigmaInv * Fn_x
+  bob::math::prod(IdPlusUSProdInv, m_tmp_ru, x);
+}
+
+
+void bob::learn::misc::FABase::estimateX(const bob::learn::misc::GMMStats& gmm_stats, blitz::Array<double,1>& x) const
+{
+  if (!m_ubm) throw std::runtime_error("No UBM was set in the JFA machine.");
+  computeIdPlusUSProdInv(gmm_stats, m_tmp_IdPlusUSProdInv); // Computes first term
+  computeFn_x(gmm_stats, m_tmp_Fn_x); // Computes last term
+  estimateX(m_tmp_IdPlusUSProdInv, m_tmp_Fn_x, x); // Estimates the value of x
+}
+
diff --git a/bob/learn/misc/cpp/ISVBase.cpp b/bob/learn/misc/cpp/ISVBase.cpp
new file mode 100644
index 0000000..7884786
--- /dev/null
+++ b/bob/learn/misc/cpp/ISVBase.cpp
@@ -0,0 +1,76 @@
+/**
+ * @date Tue Jan 27 16:02:00 2015 +0200
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+
+#include <bob.learn.misc/ISVBase.h>
+#include <bob.core/array_copy.h>
+#include <bob.math/linear.h>
+#include <bob.math/inv.h>
+#include <bob.learn.misc/LinearScoring.h>
+#include <limits>
+
+
+//////////////////// ISVBase ////////////////////
+bob::learn::misc::ISVBase::ISVBase()
+{
+}
+
+bob::learn::misc::ISVBase::ISVBase(const boost::shared_ptr<bob::learn::misc::GMMMachine> ubm,
+    const size_t ru):
+  m_base(ubm, ru, 1)
+{
+  blitz::Array<double,2>& V = m_base.updateV();
+  V = 0;
+}
+
+bob::learn::misc::ISVBase::ISVBase(const bob::learn::misc::ISVBase& other):
+  m_base(other.m_base)
+{
+}
+
+
+bob::learn::misc::ISVBase::ISVBase(bob::io::base::HDF5File& config)
+{
+  load(config);
+}
+
+bob::learn::misc::ISVBase::~ISVBase() {
+}
+
+void bob::learn::misc::ISVBase::save(bob::io::base::HDF5File& config) const
+{
+  config.setArray("U", m_base.getU());
+  config.setArray("d", m_base.getD());
+}
+
+void bob::learn::misc::ISVBase::load(bob::io::base::HDF5File& config)
+{
+  //reads all data directly into the member variables
+  blitz::Array<double,2> U = config.readArray<double,2>("U");
+  blitz::Array<double,1> d = config.readArray<double,1>("d");
+  const int ru = U.extent(1);
+  if (!m_base.getUbm())
+    m_base.resize(ru, 1, U.extent(0));
+  else
+    m_base.resize(ru, 1);
+  m_base.setU(U);
+  m_base.setD(d);
+  blitz::Array<double,2>& V = m_base.updateV();
+  V = 0;
+}
+
+bob::learn::misc::ISVBase&
+bob::learn::misc::ISVBase::operator=(const bob::learn::misc::ISVBase& other)
+{
+  if (this != &other)
+  {
+    m_base = other.m_base;
+  }
+  return *this;
+}
+
diff --git a/bob/learn/misc/cpp/ISVMachine.cpp b/bob/learn/misc/cpp/ISVMachine.cpp
new file mode 100644
index 0000000..b205337
--- /dev/null
+++ b/bob/learn/misc/cpp/ISVMachine.cpp
@@ -0,0 +1,184 @@
+/**
+ * @date Tue Jan 27 16:06:00 2015 +0200
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+
+#include <bob.learn.misc/ISVMachine.h>
+#include <bob.core/array_copy.h>
+#include <bob.math/linear.h>
+#include <bob.math/inv.h>
+#include <bob.learn.misc/LinearScoring.h>
+#include <limits>
+
+
+//////////////////// ISVMachine ////////////////////
+bob::learn::misc::ISVMachine::ISVMachine():
+  m_z(1)
+{
+  resizeTmp();
+}
+
+bob::learn::misc::ISVMachine::ISVMachine(const boost::shared_ptr<bob::learn::misc::ISVBase> isv_base):
+  m_isv_base(isv_base),
+  m_z(isv_base->getSupervectorLength())
+{
+  if (!m_isv_base->getUbm())
+    throw std::runtime_error("No UBM was set in the JFA machine.");
+  updateCache();
+  resizeTmp();
+}
+
+
+bob::learn::misc::ISVMachine::ISVMachine(const bob::learn::misc::ISVMachine& other):
+  m_isv_base(other.m_isv_base),
+  m_z(bob::core::array::ccopy(other.m_z))
+{
+  updateCache();
+  resizeTmp();
+}
+
+bob::learn::misc::ISVMachine::ISVMachine(bob::io::base::HDF5File& config)
+{
+  load(config);
+}
+
+bob::learn::misc::ISVMachine::~ISVMachine() {
+}
+
+bob::learn::misc::ISVMachine&
+bob::learn::misc::ISVMachine::operator=(const bob::learn::misc::ISVMachine& other)
+{
+  if (this != &other)
+  {
+    m_isv_base = other.m_isv_base;
+    m_z.reference(bob::core::array::ccopy(other.m_z));
+  }
+  return *this;
+}
+
+bool bob::learn::misc::ISVMachine::operator==(const bob::learn::misc::ISVMachine& other) const
+{
+  return (*m_isv_base == *(other.m_isv_base) &&
+          bob::core::array::isEqual(m_z, other.m_z));
+}
+
+bool bob::learn::misc::ISVMachine::operator!=(const bob::learn::misc::ISVMachine& b) const
+{
+  return !(this->operator==(b));
+}
+
+
+bool bob::learn::misc::ISVMachine::is_similar_to(const bob::learn::misc::ISVMachine& b,
+    const double r_epsilon, const double a_epsilon) const
+{
+  return (m_isv_base->is_similar_to(*(b.m_isv_base), r_epsilon, a_epsilon) &&
+          bob::core::array::isClose(m_z, b.m_z, r_epsilon, a_epsilon));
+}
+
+void bob::learn::misc::ISVMachine::save(bob::io::base::HDF5File& config) const
+{
+  config.setArray("z", m_z);
+}
+
+void bob::learn::misc::ISVMachine::load(bob::io::base::HDF5File& config)
+{
+  //reads all data directly into the member variables
+  blitz::Array<double,1> z = config.readArray<double,1>("z");
+  if (!m_isv_base)
+    m_z.resize(z.extent(0));
+  setZ(z);
+  // update cache
+  updateCache();
+  resizeTmp();
+}
+
+void bob::learn::misc::ISVMachine::setZ(const blitz::Array<double,1>& z)
+{
+  if(z.extent(0) != m_z.extent(0)) { //checks dimension
+    boost::format m("size of input vector `z' (%d) does not match the expected size (%d)");
+    m % z.extent(0) % m_z.extent(0);
+    throw std::runtime_error(m.str());
+  }
+  m_z.reference(bob::core::array::ccopy(z));
+  // update cache
+  updateCache();
+}
+
+void bob::learn::misc::ISVMachine::setISVBase(const boost::shared_ptr<bob::learn::misc::ISVBase> isv_base)
+{
+  if (!isv_base->getUbm())
+    throw std::runtime_error("No UBM was set in the JFA machine.");
+  m_isv_base = isv_base;
+  // Resize variables
+  resize();
+}
+
+void bob::learn::misc::ISVMachine::resize()
+{
+  m_z.resizeAndPreserve(getSupervectorLength());
+  updateCache();
+  resizeTmp();
+}
+
+void bob::learn::misc::ISVMachine::resizeTmp()
+{
+  if (m_isv_base)
+  {
+    m_tmp_Ux.resize(getSupervectorLength());
+  }
+}
+
+void bob::learn::misc::ISVMachine::updateCache()
+{
+  if (m_isv_base)
+  {
+    // m + Dz
+    m_cache_mDz.resize(getSupervectorLength());
+    m_cache_mDz = m_isv_base->getD()*m_z + m_isv_base->getUbm()->getMeanSupervector();
+    m_cache_x.resize(getDimRu());
+  }
+}
+
+void bob::learn::misc::ISVMachine::estimateUx(const bob::learn::misc::GMMStats& gmm_stats,
+  blitz::Array<double,1>& Ux)
+{
+  estimateX(gmm_stats, m_cache_x);
+  bob::math::prod(m_isv_base->getU(), m_cache_x, Ux);
+}
+
+void bob::learn::misc::ISVMachine::forward(const bob::learn::misc::GMMStats& input,
+  double& score) const
+{
+  forward_(input, score);
+}
+
+void bob::learn::misc::ISVMachine::forward(const bob::learn::misc::GMMStats& gmm_stats,
+  const blitz::Array<double,1>& Ux, double& score) const
+{
+  // Checks that a Base machine has been set
+  if (!m_isv_base) throw std::runtime_error("No UBM was set in the JFA machine.");
+
+  score = bob::learn::misc::linearScoring(m_cache_mDz,
+            m_isv_base->getUbm()->getMeanSupervector(), m_isv_base->getUbm()->getVarianceSupervector(),
+            gmm_stats, Ux, true);
+}
+
+void bob::learn::misc::ISVMachine::forward_(const bob::learn::misc::GMMStats& input,
+  double& score) const
+{
+  // Checks that a Base machine has been set
+  if(!m_isv_base) throw std::runtime_error("No UBM was set in the JFA machine.");
+
+  // Ux and GMMStats
+  estimateX(input, m_cache_x);
+  bob::math::prod(m_isv_base->getU(), m_cache_x, m_tmp_Ux);
+
+  score = bob::learn::misc::linearScoring(m_cache_mDz,
+            m_isv_base->getUbm()->getMeanSupervector(), m_isv_base->getUbm()->getVarianceSupervector(),
+            input, m_tmp_Ux, true);
+}
+
diff --git a/bob/learn/misc/cpp/JFABase.cpp b/bob/learn/misc/cpp/JFABase.cpp
new file mode 100644
index 0000000..a87e57f
--- /dev/null
+++ b/bob/learn/misc/cpp/JFABase.cpp
@@ -0,0 +1,75 @@
+/**
+ * @date Tue Jan 27 15:54:00 2015 +0200
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+
+#include <bob.learn.misc/JFABase.h>
+#include <bob.core/array_copy.h>
+#include <bob.math/linear.h>
+#include <bob.math/inv.h>
+#include <bob.learn.misc/LinearScoring.h>
+#include <limits>
+
+
+//////////////////// JFABase ////////////////////
+bob::learn::misc::JFABase::JFABase()
+{
+}
+
+bob::learn::misc::JFABase::JFABase(const boost::shared_ptr<bob::learn::misc::GMMMachine> ubm,
+    const size_t ru, const size_t rv):
+  m_base(ubm, ru, rv)
+{
+}
+
+bob::learn::misc::JFABase::JFABase(const bob::learn::misc::JFABase& other):
+  m_base(other.m_base)
+{
+}
+
+
+bob::learn::misc::JFABase::JFABase(bob::io::base::HDF5File& config)
+{
+  load(config);
+}
+
+bob::learn::misc::JFABase::~JFABase() {
+}
+
+void bob::learn::misc::JFABase::save(bob::io::base::HDF5File& config) const
+{
+  config.setArray("U", m_base.getU());
+  config.setArray("V", m_base.getV());
+  config.setArray("d", m_base.getD());
+}
+
+void bob::learn::misc::JFABase::load(bob::io::base::HDF5File& config)
+{
+  //reads all data directly into the member variables
+  blitz::Array<double,2> U = config.readArray<double,2>("U");
+  blitz::Array<double,2> V = config.readArray<double,2>("V");
+  blitz::Array<double,1> d = config.readArray<double,1>("d");
+  const int ru = U.extent(1);
+  const int rv = V.extent(1);
+  if (!m_base.getUbm())
+    m_base.resize(ru, rv, U.extent(0));
+  else
+    m_base.resize(ru, rv);
+  m_base.setU(U);
+  m_base.setV(V);
+  m_base.setD(d);
+}
+
+bob::learn::misc::JFABase&
+bob::learn::misc::JFABase::operator=(const bob::learn::misc::JFABase& other)
+{
+  if (this != &other)
+  {
+    m_base = other.m_base;
+  }
+  return *this;
+}
diff --git a/bob/learn/misc/include/bob.learn.misc/FABase.h b/bob/learn/misc/include/bob.learn.misc/FABase.h
new file mode 100644
index 0000000..35b1667
--- /dev/null
+++ b/bob/learn/misc/include/bob.learn.misc/FABase.h
@@ -0,0 +1,293 @@
+/**
+ * @date Tue Jan 27 15:51:15 2015 +0200
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * @brief A base class for Factor Analysis
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#ifndef BOB_LEARN_MISC_FABASE_H
+#define BOB_LEARN_MISC_FABASE_H
+
+#include <stdexcept>
+
+#include <bob.learn.misc/GMMMachine.h>
+#include <boost/shared_ptr.hpp>
+
+namespace bob { namespace learn { namespace misc {
+
+/**
+ * @brief A FA Base class which contains U, V and D matrices
+ * TODO: add a reference to the journal articles
+ */
+class FABase
+{
+  public:
+    /**
+     * @brief Default constructor. Builds an otherwise invalid 0 x 0 FABase
+     * The Universal Background Model and the matrices U, V and diag(d) are
+     * not initialized.
+     */
+    FABase();
+
+    /**
+     * @brief Constructor. Builds a new FABase.
+     * The Universal Background Model and the matrices U, V and diag(d) are
+     * not initialized.
+     *
+     * @param ubm The Universal Background Model
+     * @param ru size of U (CD x ru)
+     * @param rv size of U (CD x rv)
+     * @warning ru and rv SHOULD BE  >= 1. Just set U/V/D to zero if you want
+     *   to ignore one subspace. This is the case for ISV.
+     */
+    FABase(const boost::shared_ptr<bob::learn::misc::GMMMachine> ubm, const size_t ru=1, const size_t rv=1);
+
+    /**
+     * @brief Copy constructor
+     */
+    FABase(const FABase& other);
+
+    /**
+     * @brief Just to virtualise the destructor
+     */
+    virtual ~FABase();
+
+    /**
+     * @brief Assigns from a different JFA machine
+     */
+    FABase& operator=(const FABase &other);
+
+    /**
+     * @brief Equal to
+     */
+    bool operator==(const FABase& b) const;
+
+    /**
+     * @brief Not equal to
+     */
+    bool operator!=(const FABase& b) const;
+
+    /**
+     * @brief Similar to
+     */
+    bool is_similar_to(const FABase& b, const double r_epsilon=1e-5,
+      const double a_epsilon=1e-8) const;
+
+    /**
+     * @brief Returns the UBM
+     */
+    const boost::shared_ptr<bob::learn::misc::GMMMachine> getUbm() const
+    { return m_ubm; }
+
+    /**
+     * @brief Returns the U matrix
+     */
+    const blitz::Array<double,2>& getU() const
+    { return m_U; }
+
+    /**
+     * @brief Returns the V matrix
+     */
+    const blitz::Array<double,2>& getV() const
+    { return m_V; }
+
+    /**
+     * @brief Returns the diagonal matrix diag(d) (as a 1D vector)
+     */
+    const blitz::Array<double,1>& getD() const
+    { return m_d; }
+
+    /**
+     * @brief Returns the UBM mean supervector (as a 1D vector)
+     */
+    const blitz::Array<double,1>& getUbmMean() const
+    { return m_cache_mean; }
+
+    /**
+     * @brief Returns the UBM variance supervector (as a 1D vector)
+     */
+    const blitz::Array<double,1>& getUbmVariance() const
+    { return m_cache_sigma; }
+
+    /**
+     * @brief Returns the number of Gaussian components C
+     * @warning An exception is thrown if no Universal Background Model has
+     *   been set yet.
+     */
+    const size_t getNGaussians() const
+    { if(!m_ubm) throw std::runtime_error("No UBM was set in the JFA machine.");
+      return m_ubm->getNGaussians(); }
+
+    /**
+     * @brief Returns the feature dimensionality D
+     * @warning An exception is thrown if no Universal Background Model has
+     *   been set yet.
+     */
+    const size_t getNInputs() const
+    { if(!m_ubm) throw std::runtime_error("No UBM was set in the JFA machine.");
+      return m_ubm->getNInputs(); }
+
+    /**
+     * @brief Returns the supervector length CD
+     * (CxD: Number of Gaussian components by the feature dimensionality)
+     * @warning An exception is thrown if no Universal Background Model has
+     *   been set yet.
+     */
+    const size_t getSupervectorLength() const
+    { if(!m_ubm) throw std::runtime_error("No UBM was set in the JFA machine.");
+      return m_ubm->getNInputs()*m_ubm->getNGaussians(); }
+
+    /**
+     * @brief Returns the size/rank ru of the U matrix
+     */
+    const size_t getDimRu() const
+    { return m_ru; }
+
+    /**
+     * @brief Returns the size/rank rv of the V matrix
+     */
+    const size_t getDimRv() const
+    { return m_rv; }
+
+    /**
+     * @brief Resets the dimensionality of the subspace U and V
+     * U and V are hence uninitialized.
+     */
+    void resize(const size_t ru, const size_t rv);
+
+    /**
+     * @brief Resets the dimensionality of the subspace U and V,
+     * assuming that no UBM has yet been set
+     * U and V are hence uninitialized.
+     */
+    void resize(const size_t ru, const size_t rv, const size_t cd);
+
+    /**
+     * @brief Returns the U matrix in order to update it
+     * @warning Should only be used by the trainer for efficiency reason,
+     *   or for testing purpose.
+     */
+    blitz::Array<double,2>& updateU()
+    { return m_U; }
+
+    /**
+     * @brief Returns the V matrix in order to update it
+     * @warning Should only be used by the trainer for efficiency reason,
+     *   or for testing purpose.
+     */
+    blitz::Array<double,2>& updateV()
+    { return m_V; }
+
+    /**
+     * @brief Returns the diagonal matrix diag(d) (as a 1D vector) in order
+     * to update it
+     * @warning Should only be used by the trainer for efficiency reason,
+     *   or for testing purpose.
+     */
+    blitz::Array<double,1>& updateD()
+    { return m_d; }
+
+
+    /**
+     * @brief Sets (the mean supervector of) the Universal Background Model
+     * U, V and d are uninitialized in case of dimensions update (C or D)
+     */
+    void setUbm(const boost::shared_ptr<bob::learn::misc::GMMMachine> ubm);
+
+    /**
+     * @brief Sets the U matrix
+     */
+    void setU(const blitz::Array<double,2>& U);
+
+    /**
+     * @brief Sets the V matrix
+     */
+    void setV(const blitz::Array<double,2>& V);
+
+    /**
+     * @brief Sets the diagonal matrix diag(d)
+     * (a 1D vector is expected as an argument)
+     */
+    void setD(const blitz::Array<double,1>& d);
+
+
+    /**
+     * @brief Estimates x from the GMM statistics considering the LPT
+     * assumption, that is the latent session variable x is approximated
+     * using the UBM
+     */
+    void estimateX(const bob::learn::misc::GMMStats& gmm_stats, blitz::Array<double,1>& x) const;
+
+    /**
+     * @brief Compute and put U^{T}.Sigma^{-1} matrix in cache
+     * @warning Should only be used by the trainer for efficiency reason,
+     *   or for testing purpose.
+     */
+    void updateCacheUbmUVD();
+
+
+  private:
+    /**
+     * @brief Update cache arrays/variables
+     */
+    void updateCache();
+    /**
+     * @brief Put GMM mean/variance supervector in cache
+     */
+    void updateCacheUbm();
+    /**
+     * @brief Resize working arrays
+     */
+    void resizeTmp();
+    /**
+     * @brief Computes (Id + U^T.Sigma^-1.U.N_{i,h}.U)^-1 =
+     *   (Id + sum_{c=1..C} N_{i,h}.U_{c}^T.Sigma_{c}^-1.U_{c})^-1
+     */
+    void computeIdPlusUSProdInv(const bob::learn::misc::GMMStats& gmm_stats,
+      blitz::Array<double,2>& out) const;
+    /**
+     * @brief Computes Fn_x = sum_{sessions h}(N*(o - m))
+     * (Normalised first order statistics)
+     */
+    void computeFn_x(const bob::learn::misc::GMMStats& gmm_stats,
+      blitz::Array<double,1>& out) const;
+    /**
+     * @brief Estimates the value of x from the passed arguments
+     * (IdPlusUSProdInv and Fn_x), considering the LPT assumption
+     */
+    void estimateX(const blitz::Array<double,2>& IdPlusUSProdInv,
+      const blitz::Array<double,1>& Fn_x, blitz::Array<double,1>& x) const;
+
+
+    // UBM
+    boost::shared_ptr<bob::learn::misc::GMMMachine> m_ubm;
+
+    // dimensionality
+    size_t m_ru; // size of U (CD x ru)
+    size_t m_rv; // size of V (CD x rv)
+
+    // U, V, D matrices
+    // D is assumed to be diagonal, and only the diagonal is stored
+    blitz::Array<double,2> m_U;
+    blitz::Array<double,2> m_V;
+    blitz::Array<double,1> m_d;
+
+    // Vectors/Matrices precomputed in cache
+    blitz::Array<double,1> m_cache_mean;
+    blitz::Array<double,1> m_cache_sigma;
+    blitz::Array<double,2> m_cache_UtSigmaInv;
+
+    mutable blitz::Array<double,2> m_tmp_IdPlusUSProdInv;
+    mutable blitz::Array<double,1> m_tmp_Fn_x;
+    mutable blitz::Array<double,1> m_tmp_ru;
+    mutable blitz::Array<double,2> m_tmp_ruD;
+    mutable blitz::Array<double,2> m_tmp_ruru;
+};
+
+
+} } } // namespaces
+
+#endif // BOB_LEARN_MISC_FABASE_H
diff --git a/bob/learn/misc/include/bob.learn.misc/ISVBase.h b/bob/learn/misc/include/bob.learn.misc/ISVBase.h
new file mode 100644
index 0000000..895a81d
--- /dev/null
+++ b/bob/learn/misc/include/bob.learn.misc/ISVBase.h
@@ -0,0 +1,228 @@
+/**
+ * @date Tue Jan 27 16:02:00 2015 +0200
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * @brief A base class for Joint Factor Analysis-like machines
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#ifndef BOB_LEARN_MISC_ISVBASE_H
+#define BOB_LEARN_MISC_ISVBASE_H
+
+#include <stdexcept>
+
+#include <bob.learn.misc/GMMMachine.h>
+#include <bob.learn.misc/FABase.h>
+
+#include <bob.io.base/HDF5File.h>
+#include <boost/shared_ptr.hpp>
+
+namespace bob { namespace learn { namespace misc {
+
+
+/**
+ * @brief An ISV Base class which contains U and D matrices
+ * TODO: add a reference to the journal articles
+ */
+class ISVBase
+{
+  public:
+    /**
+     * @brief Default constructor. Builds an otherwise invalid 0 x 0 ISVBase
+     * The Universal Background Model and the matrices U, V and diag(d) are
+     * not initialized.
+     */
+    ISVBase();
+
+    /**
+     * @brief Constructor. Builds a new ISVBase.
+     * The Universal Background Model and the matrices U, V and diag(d) are
+     * not initialized.
+     *
+     * @param ubm The Universal Background Model
+     * @param ru size of U (CD x ru)
+     * @warning ru SHOULD BE >= 1.
+     */
+    ISVBase(const boost::shared_ptr<bob::learn::misc::GMMMachine> ubm, const size_t ru=1);
+
+    /**
+     * @brief Copy constructor
+     */
+    ISVBase(const ISVBase& other);
+
+    /**
+     * @deprecated Starts a new JFAMachine from an existing Configuration object.
+     */
+    ISVBase(bob::io::base::HDF5File& config);
+
+    /**
+     * @brief Just to virtualise the destructor
+     */
+    virtual ~ISVBase();
+
+    /**
+     * @brief Assigns from a different JFA machine
+     */
+    ISVBase& operator=(const ISVBase &other);
+
+    /**
+     * @brief Equal to
+     */
+    bool operator==(const ISVBase& b) const
+    { return m_base.operator==(b.m_base); }
+
+    /**
+     * @brief Not equal to
+     */
+    bool operator!=(const ISVBase& b) const
+    { return m_base.operator!=(b.m_base); }
+
+    /**
+     * @brief Similar to
+     */
+    bool is_similar_to(const ISVBase& b, const double r_epsilon=1e-5,
+      const double a_epsilon=1e-8) const
+    { return m_base.is_similar_to(b.m_base, r_epsilon, a_epsilon); }
+
+    /**
+     * @brief Saves machine to an HDF5 file
+     */
+    void save(bob::io::base::HDF5File& config) const;
+
+    /**
+     * @brief Loads data from an existing configuration object. Resets
+     * the current state.
+     */
+    void load(bob::io::base::HDF5File& config);
+
+    /**
+     * @brief Returns the UBM
+     */
+    const boost::shared_ptr<bob::learn::misc::GMMMachine> getUbm() const
+    { return m_base.getUbm(); }
+
+    /**
+     * @brief Returns the U matrix
+     */
+    const blitz::Array<double,2>& getU() const
+    { return m_base.getU(); }
+
+    /**
+     * @brief Returns the diagonal matrix diag(d) (as a 1D vector)
+     */
+    const blitz::Array<double,1>& getD() const
+    { return m_base.getD(); }
+
+    /**
+     * @brief Returns the number of Gaussian components C
+     * @warning An exception is thrown if no Universal Background Model has
+     *   been set yet.
+     */
+    const size_t getNGaussians() const
+    { return m_base.getNGaussians(); }
+
+    /**
+     * @brief Returns the feature dimensionality D
+     * @warning An exception is thrown if no Universal Background Model has
+     *   been set yet.
+     */
+    const size_t getNInputs() const
+    { return m_base.getNInputs(); }
+
+    /**
+     * @brief Returns the supervector length CD
+     * (CxD: Number of Gaussian components by the feature dimensionality)
+     * @warning An exception is thrown if no Universal Background Model has
+     *   been set yet.
+     */
+    const size_t getSupervectorLength() const
+    { return m_base.getSupervectorLength(); }
+
+    /**
+     * @brief Returns the size/rank ru of the U matrix
+     */
+    const size_t getDimRu() const
+    { return m_base.getDimRu(); }
+
+    /**
+     * @brief Resets the dimensionality of the subspace U and V
+     * U and V are hence uninitialized.
+     */
+    void resize(const size_t ru)
+    { m_base.resize(ru, 1);
+      blitz::Array<double,2>& V = m_base.updateV();
+      V = 0;
+     }
+
+    /**
+     * @brief Returns the U matrix in order to update it
+     * @warning Should only be used by the trainer for efficiency reason,
+     *   or for testing purpose.
+     */
+    blitz::Array<double,2>& updateU()
+    { return m_base.updateU(); }
+
+    /**
+     * @brief Returns the diagonal matrix diag(d) (as a 1D vector) in order
+     * to update it
+     * @warning Should only be used by the trainer for efficiency reason,
+     *   or for testing purpose.
+     */
+    blitz::Array<double,1>& updateD()
+    { return m_base.updateD(); }
+
+
+    /**
+     * @brief Sets (the mean supervector of) the Universal Background Model
+     * U, V and d are uninitialized in case of dimensions update (C or D)
+     */
+    void setUbm(const boost::shared_ptr<bob::learn::misc::GMMMachine> ubm)
+    { m_base.setUbm(ubm); }
+
+    /**
+     * @brief Sets the U matrix
+     */
+    void setU(const blitz::Array<double,2>& U)
+    { m_base.setU(U); }
+
+    /**
+     * @brief Sets the diagonal matrix diag(d)
+     * (a 1D vector is expected as an argument)
+     */
+    void setD(const blitz::Array<double,1>& d)
+    { m_base.setD(d); }
+
+    /**
+     * @brief Estimates x from the GMM statistics considering the LPT
+     * assumption, that is the latent session variable x is approximated
+     * using the UBM
+     */
+    void estimateX(const bob::learn::misc::GMMStats& gmm_stats, blitz::Array<double,1>& x) const
+    { m_base.estimateX(gmm_stats, x); }
+
+    /**
+     * @brief Precompute (put U^{T}.Sigma^{-1} matrix in cache)
+     * @warning Should only be used by the trainer for efficiency reason,
+     *   or for testing purpose.
+     */
+    void precompute()
+    { m_base.updateCacheUbmUVD(); }
+
+    /**
+     * @brief Returns the FABase member
+     */
+    const bob::learn::misc::FABase& getBase() const
+    { return m_base; }
+
+
+  private:
+    // FABase
+    bob::learn::misc::FABase m_base;
+};
+
+
+} } } // namespaces
+
+#endif // BOB_LEARN_MISC_JFABASE_H
diff --git a/bob/learn/misc/include/bob.learn.misc/ISVMachine.h b/bob/learn/misc/include/bob.learn.misc/ISVMachine.h
new file mode 100644
index 0000000..83b8795
--- /dev/null
+++ b/bob/learn/misc/include/bob.learn.misc/ISVMachine.h
@@ -0,0 +1,230 @@
+/**
+ * @date Tue Jan 27 16:06:00 2015 +0200
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * @brief A base class for Joint Factor Analysis-like machines
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#ifndef BOB_LEARN_MISC_ISVMACHINE_H
+#define BOB_LEARN_MISC_ISVMACHINE_H
+
+#include <stdexcept>
+
+#include <bob.learn.misc/ISVBase.h>
+#include <bob.learn.misc/GMMMachine.h>
+#include <bob.learn.misc/LinearScoring.h>
+
+#include <bob.io.base/HDF5File.h>
+#include <boost/shared_ptr.hpp>
+
+namespace bob { namespace learn { namespace misc {
+
+
+/**
+ * @brief A ISVMachine which is associated to a ISVBase that contains
+ *   U D matrices.
+ * TODO: add a reference to the journal articles
+ */
+class ISVMachine
+{
+  public:
+    /**
+     * @brief Default constructor. Builds an otherwise invalid 0 x 0 ISVMachine
+     * The Universal Background Model and the matrices U, V and diag(d) are
+     * not initialized.
+     */
+    ISVMachine();
+
+    /**
+     * @brief Constructor. Builds a new ISVMachine.
+     *
+     * @param isv_base The ISVBase associated with this machine
+     */
+    ISVMachine(const boost::shared_ptr<bob::learn::misc::ISVBase> isv_base);
+
+    /**
+     * @brief Copy constructor
+     */
+    ISVMachine(const ISVMachine& other);
+
+    /**
+     * @brief Starts a new ISVMachine from an existing Configuration object.
+     */
+    ISVMachine(bob::io::base::HDF5File& config);
+
+    /**
+     * @brief Just to virtualise the destructor
+     */
+    virtual ~ISVMachine();
+
+    /**
+     * @brief Assigns from a different ISV machine
+     */
+    ISVMachine& operator=(const ISVMachine &other);
+
+    /**
+     * @brief Equal to
+     */
+    bool operator==(const ISVMachine& b) const;
+
+    /**
+     * @brief Not equal to
+     */
+    bool operator!=(const ISVMachine& b) const;
+
+    /**
+     * @brief Similar to
+     */
+    bool is_similar_to(const ISVMachine& b, const double r_epsilon=1e-5,
+      const double a_epsilon=1e-8) const;
+
+    /**
+     * @brief Saves machine to an HDF5 file
+     */
+    void save(bob::io::base::HDF5File& config) const;
+
+    /**
+     * @brief Loads data from an existing configuration object. Resets
+     * the current state.
+     */
+    void load(bob::io::base::HDF5File& config);
+
+
+    /**
+     * @brief Returns the number of Gaussian components C
+     * @warning An exception is thrown if no Universal Background Model has
+     *   been set yet.
+     */
+    const size_t getNGaussians() const
+    { return m_isv_base->getNGaussians(); }
+
+    /**
+     * @brief Returns the feature dimensionality D
+     * @warning An exception is thrown if no Universal Background Model has
+     *   been set yet.
+     */
+    const size_t getNInputs() const
+    { return m_isv_base->getNInputs(); }
+
+    /**
+     * @brief Returns the supervector length CD
+     * (CxD: Number of Gaussian components by the feature dimensionality)
+     * @warning An exception is thrown if no Universal Background Model has
+     *   been set yet.
+     */
+    const size_t getSupervectorLength() const
+    { return m_isv_base->getSupervectorLength(); }
+
+    /**
+     * @brief Returns the size/rank ru of the U matrix
+     */
+    const size_t getDimRu() const
+    { return m_isv_base->getDimRu(); }
+
+    /**
+     * @brief Returns the x session factor
+     */
+    const blitz::Array<double,1>& getX() const
+    { return m_cache_x; }
+
+    /**
+     * @brief Returns the z speaker factor
+     */
+    const blitz::Array<double,1>& getZ() const
+    { return m_z; }
+
+    /**
+     * @brief Returns the z speaker factors in order to update it
+     */
+    blitz::Array<double,1>& updateZ()
+    { return m_z; }
+
+    /**
+     * @brief Returns the V matrix
+     */
+    void setZ(const blitz::Array<double,1>& z);
+
+    /**
+     * @brief Returns the ISVBase
+     */
+    const boost::shared_ptr<bob::learn::misc::ISVBase> getISVBase() const
+    { return m_isv_base; }
+
+    /**
+     * @brief Sets the ISVBase
+     */
+    void setISVBase(const boost::shared_ptr<bob::learn::misc::ISVBase> isv_base);
+
+
+    /**
+     * @brief Estimates x from the GMM statistics considering the LPT
+     * assumption, that is the latent session variable x is approximated
+     * using the UBM
+     */
+    void estimateX(const bob::learn::misc::GMMStats& gmm_stats, blitz::Array<double,1>& x) const
+    { m_isv_base->estimateX(gmm_stats, x); }
+    /**
+     * @brief Estimates Ux from the GMM statistics considering the LPT
+     * assumption, that is the latent session variable x is approximated
+     * using the UBM
+     */
+    void estimateUx(const bob::learn::misc::GMMStats& gmm_stats, blitz::Array<double,1>& Ux);
+
+   /**
+    * @brief Execute the machine
+    *
+    * @param input input data used by the machine
+    * @param score value computed by the machine
+    * @warning Inputs are checked
+    */
+    void forward(const bob::learn::misc::GMMStats& input, double& score) const;
+    /**
+     * @brief Computes a score for the given UBM statistics and given the
+     * Ux vector
+     */
+    void forward(const bob::learn::misc::GMMStats& gmm_stats,
+      const blitz::Array<double,1>& Ux, double& score) const;
+
+    /**
+     * @brief Execute the machine
+     *
+     * @param input input data used by the machine
+     * @param score value computed by the machine
+     * @warning Inputs are NOT checked
+     */
+    void forward_(const bob::learn::misc::GMMStats& input, double& score) const;
+
+  private:
+    /**
+     * @brief Resize latent variable according to the ISVBase
+     */
+    void resize();
+    /**
+     * @ Update cache
+     */
+    void updateCache();
+    /**
+     * @brief Resize working arrays
+     */
+    void resizeTmp();
+
+    // UBM
+    boost::shared_ptr<bob::learn::misc::ISVBase> m_isv_base;
+
+    // y and z vectors/factors learned during the enrolment procedure
+    blitz::Array<double,1> m_z;
+
+    // cache
+    blitz::Array<double,1> m_cache_mDz;
+    mutable blitz::Array<double,1> m_cache_x;
+
+    // x vector/factor in cache when computing scores
+    mutable blitz::Array<double,1> m_tmp_Ux;
+};
+
+} } } // namespaces
+
+#endif // BOB_LEARN_MISC_ISVMACHINE_H
diff --git a/bob/learn/misc/include/bob.learn.misc/JFABase.h b/bob/learn/misc/include/bob.learn.misc/JFABase.h
new file mode 100644
index 0000000..7fbc669
--- /dev/null
+++ b/bob/learn/misc/include/bob.learn.misc/JFABase.h
@@ -0,0 +1,253 @@
+/**
+ * @date Tue Jan 27 15:54:00 2015 +0200
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * @brief A base class for Joint Factor Analysis-like machines
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#ifndef BOB_LEARN_MISC_JFABASE_H
+#define BOB_LEARN_MISC_JFABASE_H
+
+#include <stdexcept>
+
+#include <bob.learn.misc/GMMMachine.h>
+#include <bob.learn.misc/FABase.h>
+//#include <bob.learn.misc/LinearScoring.h>
+
+#include <bob.io.base/HDF5File.h>
+#include <boost/shared_ptr.hpp>
+
+namespace bob { namespace learn { namespace misc {
+
+
+/**
+ * @brief A JFA Base class which contains U, V and D matrices
+ * TODO: add a reference to the journal articles
+ */
+class JFABase
+{
+  public:
+    /**
+     * @brief Default constructor. Builds a 1 x 1 JFABase
+     * The Universal Background Model and the matrices U, V and diag(d) are
+     * not initialized.
+     */
+    JFABase();
+
+    /**
+     * @brief Constructor. Builds a new JFABase.
+     * The Universal Background Model and the matrices U, V and diag(d) are
+     * not initialized.
+     *
+     * @param ubm The Universal Background Model
+     * @param ru size of U (CD x ru)
+     * @param rv size of U (CD x rv)
+     * @warning ru and rv SHOULD BE  >= 1.
+     */
+    JFABase(const boost::shared_ptr<bob::learn::misc::GMMMachine> ubm, const size_t ru=1, const size_t rv=1);
+
+    /**
+     * @brief Copy constructor
+     */
+    JFABase(const JFABase& other);
+
+    /**
+     * @deprecated Starts a new JFAMachine from an existing Configuration object.
+     */
+    JFABase(bob::io::base::HDF5File& config);
+
+    /**
+     * @brief Just to virtualise the destructor
+     */
+    virtual ~JFABase();
+
+    /**
+     * @brief Assigns from a different JFA machine
+     */
+    JFABase& operator=(const JFABase &other);
+
+    /**
+     * @brief Equal to
+     */
+    bool operator==(const JFABase& b) const
+    { return m_base.operator==(b.m_base); }
+
+    /**
+     * @brief Not equal to
+     */
+    bool operator!=(const JFABase& b) const
+    { return m_base.operator!=(b.m_base); }
+
+    /**
+     * @brief Similar to
+     */
+    bool is_similar_to(const JFABase& b, const double r_epsilon=1e-5,
+      const double a_epsilon=1e-8) const
+    { return m_base.is_similar_to(b.m_base, r_epsilon, a_epsilon); }
+
+    /**
+     * @brief Saves model to an HDF5 file
+     */
+    void save(bob::io::base::HDF5File& config) const;
+
+    /**
+     * @brief Loads data from an existing configuration object. Resets
+     * the current state.
+     */
+    void load(bob::io::base::HDF5File& config);
+
+    /**
+     * @brief Returns the UBM
+     */
+    const boost::shared_ptr<bob::learn::misc::GMMMachine> getUbm() const
+    { return m_base.getUbm(); }
+
+    /**
+     * @brief Returns the U matrix
+     */
+    const blitz::Array<double,2>& getU() const
+    { return m_base.getU(); }
+
+    /**
+     * @brief Returns the V matrix
+     */
+    const blitz::Array<double,2>& getV() const
+    { return m_base.getV(); }
+
+    /**
+     * @brief Returns the diagonal matrix diag(d) (as a 1D vector)
+     */
+    const blitz::Array<double,1>& getD() const
+    { return m_base.getD(); }
+
+    /**
+     * @brief Returns the number of Gaussian components
+     * @warning An exception is thrown if no Universal Background Model has
+     *   been set yet.
+     */
+    const size_t getNGaussians() const
+    { return m_base.getNGaussians();}
+
+    /**
+     * @brief Returns the feature dimensionality D
+     * @warning An exception is thrown if no Universal Background Model has
+     *   been set yet.
+     */
+    const size_t getNInputs() const
+    { return m_base.getNInputs(); }
+
+    /**
+     * @brief Returns the supervector length CD
+     * (CxD: Number of Gaussian components by the feature dimensionality)
+     * @warning An exception is thrown if no Universal Background Model has
+     *   been set yet.
+     */
+    const size_t getSupervectorLength() const
+    { return m_base.getSupervectorLength(); }
+
+    /**
+     * @brief Returns the size/rank ru of the U matrix
+     */
+    const size_t getDimRu() const
+    { return m_base.getDimRu(); }
+
+    /**
+     * @brief Returns the size/rank rv of the V matrix
+     */
+    const size_t getDimRv() const
+    { return m_base.getDimRv(); }
+
+    /**
+     * @brief Resets the dimensionality of the subspace U and V
+     * U and V are hence uninitialized.
+     */
+    void resize(const size_t ru, const size_t rv)
+    { m_base.resize(ru, rv); }
+
+    /**
+     * @brief Returns the U matrix in order to update it
+     * @warning Should only be used by the trainer for efficiency reason,
+     *   or for testing purpose.
+     */
+    blitz::Array<double,2>& updateU()
+    { return m_base.updateU(); }
+
+    /**
+     * @brief Returns the V matrix in order to update it
+     * @warning Should only be used by the trainer for efficiency reason,
+     *   or for testing purpose.
+     */
+    blitz::Array<double,2>& updateV()
+    { return m_base.updateV(); }
+
+    /**
+     * @brief Returns the diagonal matrix diag(d) (as a 1D vector) in order
+     * to update it
+     * @warning Should only be used by the trainer for efficiency reason,
+     *   or for testing purpose.
+     */
+    blitz::Array<double,1>& updateD()
+    { return m_base.updateD(); }
+
+
+    /**
+     * @brief Sets (the mean supervector of) the Universal Background Model
+     * U, V and d are uninitialized in case of dimensions update (C or D)
+     */
+    void setUbm(const boost::shared_ptr<bob::learn::misc::GMMMachine> ubm)
+    { m_base.setUbm(ubm); }
+
+    /**
+     * @brief Sets the U matrix
+     */
+    void setU(const blitz::Array<double,2>& U)
+    { m_base.setU(U); }
+
+    /**
+     * @brief Sets the V matrix
+     */
+    void setV(const blitz::Array<double,2>& V)
+    { m_base.setV(V); }
+
+    /**
+     * @brief Sets the diagonal matrix diag(d)
+     * (a 1D vector is expected as an argument)
+     */
+    void setD(const blitz::Array<double,1>& d)
+    { m_base.setD(d); }
+
+    /**
+     * @brief Estimates x from the GMM statistics considering the LPT
+     * assumption, that is the latent session variable x is approximated
+     * using the UBM
+     */
+    void estimateX(const bob::learn::misc::GMMStats& gmm_stats, blitz::Array<double,1>& x) const
+    { m_base.estimateX(gmm_stats, x); }
+
+    /**
+     * @brief Precompute (put U^{T}.Sigma^{-1} matrix in cache)
+     * @warning Should only be used by the trainer for efficiency reason,
+     *   or for testing purpose.
+     */
+    void precompute()
+    { m_base.updateCacheUbmUVD(); }
+
+    /**
+     * @brief Returns the FABase member
+     */
+    const bob::learn::misc::FABase& getBase() const
+    { return m_base; }
+
+
+  private:
+    // FABase
+    bob::learn::misc::FABase m_base;
+};
+
+
+} } } // namespaces
+
+#endif // BOB_LEARN_MISC_JFABASE_H
-- 
GitLab