diff --git a/bob/learn/misc/cpp/ISVMachine.cpp b/bob/learn/misc/cpp/ISVMachine.cpp
index b205337fde165e6374c4f994dd8e9f667b9524e6..ea5ad28ba224c7cbaa5359b5c7e6bcb8b9268098 100644
--- a/bob/learn/misc/cpp/ISVMachine.cpp
+++ b/bob/learn/misc/cpp/ISVMachine.cpp
@@ -150,25 +150,23 @@ void bob::learn::misc::ISVMachine::estimateUx(const bob::learn::misc::GMMStats&
   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
+double bob::learn::misc::ISVMachine::forward(const bob::learn::misc::GMMStats& input)
 {
-  forward_(input, score);
+  return forward_(input);
 }
 
-void bob::learn::misc::ISVMachine::forward(const bob::learn::misc::GMMStats& gmm_stats,
-  const blitz::Array<double,1>& Ux, double& score) const
+double bob::learn::misc::ISVMachine::forward(const bob::learn::misc::GMMStats& gmm_stats,
+  const blitz::Array<double,1>& Ux)
 {
   // 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,
+  return 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
+double bob::learn::misc::ISVMachine::forward_(const bob::learn::misc::GMMStats& input)
 {
   // Checks that a Base machine has been set
   if(!m_isv_base) throw std::runtime_error("No UBM was set in the JFA machine.");
@@ -177,7 +175,7 @@ void bob::learn::misc::ISVMachine::forward_(const bob::learn::misc::GMMStats& in
   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,
+  return 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/JFAMachine.cpp b/bob/learn/misc/cpp/JFAMachine.cpp
index fb521bfe5b915f7ebb587ba0798bcc281c1f1fcf..ae4d1d2f9a6ea848b619151c8a643b77cffe99c4 100644
--- a/bob/learn/misc/cpp/JFAMachine.cpp
+++ b/bob/learn/misc/cpp/JFAMachine.cpp
@@ -174,25 +174,23 @@ void bob::learn::misc::JFAMachine::estimateUx(const bob::learn::misc::GMMStats&
   bob::math::prod(m_jfa_base->getU(), m_cache_x, Ux);
 }
 
-void bob::learn::misc::JFAMachine::forward(const bob::learn::misc::GMMStats& input,
-  double& score) const
+double bob::learn::misc::JFAMachine::forward(const bob::learn::misc::GMMStats& input)
 {
-  forward_(input, score);
+  return forward_(input);
 }
 
-void bob::learn::misc::JFAMachine::forward(const bob::learn::misc::GMMStats& gmm_stats,
-  const blitz::Array<double,1>& Ux, double& score) const
+double bob::learn::misc::JFAMachine::forward(const bob::learn::misc::GMMStats& gmm_stats,
+  const blitz::Array<double,1>& Ux)
 {
   // Checks that a Base machine has been set
   if (!m_jfa_base) throw std::runtime_error("No UBM was set in the JFA machine.");
 
-  score = bob::learn::misc::linearScoring(m_cache_mVyDz,
+  return bob::learn::misc::linearScoring(m_cache_mVyDz,
             m_jfa_base->getUbm()->getMeanSupervector(), m_jfa_base->getUbm()->getVarianceSupervector(),
             gmm_stats, Ux, true);
 }
 
-void bob::learn::misc::JFAMachine::forward_(const bob::learn::misc::GMMStats& input,
-  double& score) const
+double bob::learn::misc::JFAMachine::forward_(const bob::learn::misc::GMMStats& input)
 {
   // Checks that a Base machine has been set
   if (!m_jfa_base) throw std::runtime_error("No UBM was set in the JFA machine.");
@@ -201,7 +199,7 @@ void bob::learn::misc::JFAMachine::forward_(const bob::learn::misc::GMMStats& in
   estimateX(input, m_cache_x);
   bob::math::prod(m_jfa_base->getU(), m_cache_x, m_tmp_Ux);
 
-  score = bob::learn::misc::linearScoring(m_cache_mVyDz,
+  return bob::learn::misc::linearScoring(m_cache_mVyDz,
             m_jfa_base->getUbm()->getMeanSupervector(), m_jfa_base->getUbm()->getVarianceSupervector(),
             input, m_tmp_Ux, true);
 }
diff --git a/bob/learn/misc/include/bob.learn.misc/ISVBase.h b/bob/learn/misc/include/bob.learn.misc/ISVBase.h
index 895a81d4dcc438e14adc97aa9363e29d32ce8ba7..c062076362582d878b2827b8aea0e924ba4811c8 100644
--- a/bob/learn/misc/include/bob.learn.misc/ISVBase.h
+++ b/bob/learn/misc/include/bob.learn.misc/ISVBase.h
@@ -147,8 +147,8 @@ class ISVBase
     { return m_base.getDimRu(); }
 
     /**
-     * @brief Resets the dimensionality of the subspace U and V
-     * U and V are hence uninitialized.
+     * @brief Resets the dimensionality of the subspace U
+     * U is hence uninitialized.
      */
     void resize(const size_t ru)
     { m_base.resize(ru, 1);
diff --git a/bob/learn/misc/include/bob.learn.misc/ISVMachine.h b/bob/learn/misc/include/bob.learn.misc/ISVMachine.h
index 83b8795551efd547c95a6adbda1f265696b6668f..16adfb92f5c22bdb08dc3f596ae7fecfafaab54e 100644
--- a/bob/learn/misc/include/bob.learn.misc/ISVMachine.h
+++ b/bob/learn/misc/include/bob.learn.misc/ISVMachine.h
@@ -177,25 +177,25 @@ class ISVMachine
     * @brief Execute the machine
     *
     * @param input input data used by the machine
-    * @param score value computed by the machine
     * @warning Inputs are checked
+    * @return score value computed by the machine    
     */
-    void forward(const bob::learn::misc::GMMStats& input, double& score) const;
+    double forward(const bob::learn::misc::GMMStats& input);
     /**
      * @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;
+    double forward(const bob::learn::misc::GMMStats& gmm_stats,
+      const 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 NOT checked
+     * @return score value computed by the machine     
      */
-    void forward_(const bob::learn::misc::GMMStats& input, double& score) const;
+    double forward_(const bob::learn::misc::GMMStats& input);
 
   private:
     /**
diff --git a/bob/learn/misc/include/bob.learn.misc/JFAMachine.h b/bob/learn/misc/include/bob.learn.misc/JFAMachine.h
index 48d1f4a48ea10d958ad4dbd642b98c7f832f38fa..54c93b9ef4e85b99218f730341e908de00876837 100644
--- a/bob/learn/misc/include/bob.learn.misc/JFAMachine.h
+++ b/bob/learn/misc/include/bob.learn.misc/JFAMachine.h
@@ -200,16 +200,16 @@ class JFAMachine
     * @brief Execute the machine
     *
     * @param input input data used by the machine
-    * @param score value computed by the machine
     * @warning Inputs are checked
+    * @return score value computed by the machine
     */
-    void forward(const bob::learn::misc::GMMStats& input, double& score) const;
+    double forward(const bob::learn::misc::GMMStats& input);
     /**
      * @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;
+    double forward(const bob::learn::misc::GMMStats& gmm_stats,
+      const blitz::Array<double,1>& Ux);
 
     /**
      * @brief Execute the machine
@@ -218,7 +218,7 @@ class JFAMachine
      * @param score value computed by the machine
      * @warning Inputs are NOT checked
      */
-    void forward_(const bob::learn::misc::GMMStats& input, double& score) const;
+    double forward_(const bob::learn::misc::GMMStats& input);
 
   private:
     /**
diff --git a/bob/learn/misc/isv_base.cpp b/bob/learn/misc/isv_base.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ca169377dc62fc6d353e9057f26197686fac63c3
--- /dev/null
+++ b/bob/learn/misc/isv_base.cpp
@@ -0,0 +1,523 @@
+/**
+ * @date Wed Jan 28 11:13:15 2015 +0200
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * @brief Python API for bob::learn::em
+ *
+ * Copyright (C) 2011-2014 Idiap Research Institute, Martigny, Switzerland
+ */
+
+#include "main.h"
+
+/******************************************************************/
+/************ Constructor Section *********************************/
+/******************************************************************/
+
+static auto ISVBase_doc = bob::extension::ClassDoc(
+  BOB_EXT_MODULE_PREFIX ".ISVBase",
+
+  "A ISVBase instance can be seen as a container for U and D when performing Joint Factor Analysis (JFA)."
+  "References: [Vogt2008,McCool2013]",
+  ""
+).add_constructor(
+  bob::extension::FunctionDoc(
+    "__init__",
+    "Creates a ISVBase",
+    "",
+    true
+  )
+  .add_prototype("gmm,ru","")
+  .add_prototype("other","")
+  .add_prototype("hdf5","")
+  .add_prototype("","")
+
+  .add_parameter("gmm", ":py:class:`bob.learn.misc.GMMMachine`", "The Universal Background Model.")
+  .add_parameter("ru", "int", "Size of U (Within client variation matrix). In the end the U matrix will have (number_of_gaussians * feature_dimension x ru)")
+  .add_parameter("other", ":py:class:`bob.learn.misc.ISVBase`", "A ISVBase object to be copied.")
+  .add_parameter("hdf5", ":py:class:`bob.io.base.HDF5File`", "An HDF5 file open for reading")
+
+);
+
+
+static int PyBobLearnMiscISVBase_init_copy(PyBobLearnMiscISVBaseObject* self, PyObject* args, PyObject* kwargs) {
+
+  char** kwlist = ISVBase_doc.kwlist(1);
+  PyBobLearnMiscISVBaseObject* o;
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!", kwlist, &PyBobLearnMiscISVBase_Type, &o)){
+    ISVBase_doc.print_usage();
+    return -1;
+  }
+
+  self->cxx.reset(new bob::learn::misc::ISVBase(*o->cxx));
+  return 0;
+}
+
+
+static int PyBobLearnMiscISVBase_init_hdf5(PyBobLearnMiscISVBaseObject* self, PyObject* args, PyObject* kwargs) {
+
+  char** kwlist = ISVBase_doc.kwlist(2);
+
+  PyBobIoHDF5FileObject* config = 0;
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&", kwlist, &PyBobIoHDF5File_Converter, &config)){
+    ISVBase_doc.print_usage();
+    return -1;
+  }
+
+  self->cxx.reset(new bob::learn::misc::ISVBase(*(config->f)));
+
+  return 0;
+}
+
+
+static int PyBobLearnMiscISVBase_init_ubm(PyBobLearnMiscISVBaseObject* self, PyObject* args, PyObject* kwargs) {
+
+  char** kwlist = ISVBase_doc.kwlist(0);
+  
+  PyBobLearnMiscGMMMachineObject* ubm;
+  int ru = 1;
+
+  //Here we have to select which keyword argument to read  
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!i", kwlist, &PyBobLearnMiscGMMMachine_Type, &ubm, &ru)){
+    ISVBase_doc.print_usage();
+    return -1;
+  }
+  
+  self->cxx.reset(new bob::learn::misc::ISVBase(ubm->cxx, ru));
+  return 0;
+}
+
+
+static int PyBobLearnMiscISVBase_init(PyBobLearnMiscISVBaseObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+
+  // get the number of command line arguments
+  int nargs = (args?PyTuple_Size(args):0) + (kwargs?PyDict_Size(kwargs):0);
+    
+  switch (nargs) {
+
+    case 1:{
+      //Reading the input argument
+      PyObject* arg = 0;
+      if (PyTuple_Size(args))
+        arg = PyTuple_GET_ITEM(args, 0);
+      else {
+        PyObject* tmp = PyDict_Values(kwargs);
+        auto tmp_ = make_safe(tmp);
+        arg = PyList_GET_ITEM(tmp, 0);
+      }
+
+      // If the constructor input is Gaussian object
+     if (PyBobLearnMiscISVBase_Check(arg))
+       return PyBobLearnMiscISVBase_init_copy(self, args, kwargs);
+      // If the constructor input is a HDF5
+     else if (PyBobIoHDF5File_Check(arg))
+       return PyBobLearnMiscISVBase_init_hdf5(self, args, kwargs);
+    }
+    case 2:
+      return PyBobLearnMiscISVBase_init_ubm(self, args, kwargs);
+    default:
+      PyErr_Format(PyExc_RuntimeError, "number of arguments mismatch - %s requires 1 or 2 arguments, but you provided %d (see help)", Py_TYPE(self)->tp_name, nargs);
+      ISVBase_doc.print_usage();
+      return -1;
+  }
+  BOB_CATCH_MEMBER("cannot create ISVBase", 0)
+  return 0;
+}
+
+
+static void PyBobLearnMiscISVBase_delete(PyBobLearnMiscISVBaseObject* self) {
+  self->cxx.reset();
+  Py_TYPE(self)->tp_free((PyObject*)self);
+}
+
+static PyObject* PyBobLearnMiscISVBase_RichCompare(PyBobLearnMiscISVBaseObject* self, PyObject* other, int op) {
+  BOB_TRY
+
+  if (!PyBobLearnMiscISVBase_Check(other)) {
+    PyErr_Format(PyExc_TypeError, "cannot compare `%s' with `%s'", Py_TYPE(self)->tp_name, Py_TYPE(other)->tp_name);
+    return 0;
+  }
+  auto other_ = reinterpret_cast<PyBobLearnMiscISVBaseObject*>(other);
+  switch (op) {
+    case Py_EQ:
+      if (*self->cxx==*other_->cxx) Py_RETURN_TRUE; else Py_RETURN_FALSE;
+    case Py_NE:
+      if (*self->cxx==*other_->cxx) Py_RETURN_FALSE; else Py_RETURN_TRUE;
+    default:
+      Py_INCREF(Py_NotImplemented);
+      return Py_NotImplemented;
+  }
+  BOB_CATCH_MEMBER("cannot compare ISVBase objects", 0)
+}
+
+int PyBobLearnMiscISVBase_Check(PyObject* o) {
+  return PyObject_IsInstance(o, reinterpret_cast<PyObject*>(&PyBobLearnMiscISVBase_Type));
+}
+
+
+/******************************************************************/
+/************ Variables Section ***********************************/
+/******************************************************************/
+
+/***** shape *****/
+static auto shape = bob::extension::VariableDoc(
+  "shape",
+  "(int,int, int)",
+  "A tuple that represents the number of gaussians, dimensionality of each Gaussian, dimensionality of the rU (within client variability matrix) `(#Gaussians, #Inputs, #rU)`.",
+  ""
+);
+PyObject* PyBobLearnMiscISVBase_getShape(PyBobLearnMiscISVBaseObject* self, void*) {
+  BOB_TRY
+  return Py_BuildValue("(i,i,i)", self->cxx->getNGaussians(), self->cxx->getNInputs(), self->cxx->getDimRu());
+  BOB_CATCH_MEMBER("shape could not be read", 0)
+}
+
+/***** supervector_length *****/
+static auto supervector_length = bob::extension::VariableDoc(
+  "supervector_length",
+  "int",
+
+  "Returns the supervector length."
+  "NGaussians x NInputs: Number of Gaussian components by the feature dimensionality",
+  
+  "@warning An exception is thrown if no Universal Background Model has been set yet."
+);
+PyObject* PyBobLearnMiscISVBase_getSupervectorLength(PyBobLearnMiscISVBaseObject* self, void*) {
+  BOB_TRY
+  return Py_BuildValue("i", self->cxx->getSupervectorLength());
+  BOB_CATCH_MEMBER("supervector_length could not be read", 0)
+}
+
+
+/***** u *****/
+static auto U = bob::extension::VariableDoc(
+  "u",
+  "array_like <float, 2D>",
+  "Returns the U matrix (within client variability matrix)",
+  ""
+);
+PyObject* PyBobLearnMiscISVBase_getU(PyBobLearnMiscISVBaseObject* self, void*){
+  BOB_TRY
+  return PyBlitzArrayCxx_AsConstNumpy(self->cxx->getU());
+  BOB_CATCH_MEMBER("``u`` could not be read", 0)
+}
+int PyBobLearnMiscISVBase_setU(PyBobLearnMiscISVBaseObject* self, PyObject* value, void*){
+  BOB_TRY
+  PyBlitzArrayObject* o;
+  if (!PyBlitzArray_Converter(value, &o)){
+    PyErr_Format(PyExc_RuntimeError, "%s %s expects a 2D array of floats", Py_TYPE(self)->tp_name, U.name());
+    return -1;
+  }
+  auto o_ = make_safe(o);
+  auto b = PyBlitzArrayCxx_AsBlitz<double,2>(o, "u");
+  if (!b) return -1;
+  self->cxx->setU(*b);
+  return 0;
+  BOB_CATCH_MEMBER("``u`` matrix could not be set", -1)
+}
+
+
+/***** d *****/
+static auto D = bob::extension::VariableDoc(
+  "d",
+  "array_like <float, 1D>",
+  "Returns the diagonal matrix diag(d) (as a 1D vector)",
+  ""
+);
+PyObject* PyBobLearnMiscISVBase_getD(PyBobLearnMiscISVBaseObject* self, void*){
+  BOB_TRY
+  return PyBlitzArrayCxx_AsConstNumpy(self->cxx->getD());
+  BOB_CATCH_MEMBER("``d`` could not be read", 0)
+}
+int PyBobLearnMiscISVBase_setD(PyBobLearnMiscISVBaseObject* self, PyObject* value, void*){
+  BOB_TRY
+  PyBlitzArrayObject* o;
+  if (!PyBlitzArray_Converter(value, &o)){
+    PyErr_Format(PyExc_RuntimeError, "%s %s expects a 1D array of floats", Py_TYPE(self)->tp_name, D.name());
+    return -1;
+  }
+  auto o_ = make_safe(o);
+  auto b = PyBlitzArrayCxx_AsBlitz<double,1>(o, "d");
+  if (!b) return -1;
+  self->cxx->setD(*b);
+  return 0;
+  BOB_CATCH_MEMBER("``d`` matrix could not be set", -1)
+}
+
+
+/***** ubm *****/
+static auto ubm = bob::extension::VariableDoc(
+  "ubm",
+  ":py:class:`bob.learn.misc.GMMMachine`",
+  "Returns the UBM (Universal Background Model",
+  ""
+);
+PyObject* PyBobLearnMiscISVBase_getUBM(PyBobLearnMiscISVBaseObject* self, void*){
+  BOB_TRY
+
+  boost::shared_ptr<bob::learn::misc::GMMMachine> ubm_gmmMachine = self->cxx->getUbm();
+
+  //Allocating the correspondent python object
+  PyBobLearnMiscGMMMachineObject* retval =
+    (PyBobLearnMiscGMMMachineObject*)PyBobLearnMiscGMMMachine_Type.tp_alloc(&PyBobLearnMiscGMMMachine_Type, 0);
+  retval->cxx = ubm_gmmMachine;
+
+  return Py_BuildValue("O",retval);
+  BOB_CATCH_MEMBER("ubm could not be read", 0)
+}
+int PyBobLearnMiscISVBase_setUBM(PyBobLearnMiscISVBaseObject* self, PyObject* value, void*){
+  BOB_TRY
+
+  if (!PyBobLearnMiscGMMMachine_Check(value)){
+    PyErr_Format(PyExc_RuntimeError, "%s %s expects a :py:class:`bob.learn.misc.GMMMachine`", Py_TYPE(self)->tp_name, ubm.name());
+    return -1;
+  }
+
+  PyBobLearnMiscGMMMachineObject* ubm_gmmMachine = 0;
+  PyArg_Parse(value, "O!", &PyBobLearnMiscGMMMachine_Type,&ubm_gmmMachine);
+
+  self->cxx->setUbm(ubm_gmmMachine->cxx);
+
+  return 0;
+  BOB_CATCH_MEMBER("ubm could not be set", -1)  
+}
+
+
+
+static PyGetSetDef PyBobLearnMiscISVBase_getseters[] = { 
+  {
+   shape.name(),
+   (getter)PyBobLearnMiscISVBase_getShape,
+   0,
+   shape.doc(),
+   0
+  },
+  
+  {
+   supervector_length.name(),
+   (getter)PyBobLearnMiscISVBase_getSupervectorLength,
+   0,
+   supervector_length.doc(),
+   0
+  },
+  
+  {
+   U.name(),
+   (getter)PyBobLearnMiscISVBase_getU,
+   (setter)PyBobLearnMiscISVBase_setU,
+   U.doc(),
+   0
+  },
+  
+  {
+   D.name(),
+   (getter)PyBobLearnMiscISVBase_getD,
+   (setter)PyBobLearnMiscISVBase_setD,
+   D.doc(),
+   0
+  },
+
+  {
+   ubm.name(),
+   (getter)PyBobLearnMiscISVBase_getUBM,
+   (setter)PyBobLearnMiscISVBase_setUBM,
+   ubm.doc(),
+   0
+  },
+
+
+  {0}  // Sentinel
+};
+
+
+/******************************************************************/
+/************ Functions Section ***********************************/
+/******************************************************************/
+
+
+/*** save ***/
+static auto save = bob::extension::FunctionDoc(
+  "save",
+  "Save the configuration of the ISVBase to a given HDF5 file"
+)
+.add_prototype("hdf5")
+.add_parameter("hdf5", ":py:class:`bob.io.base.HDF5File`", "An HDF5 file open for writing");
+static PyObject* PyBobLearnMiscISVBase_Save(PyBobLearnMiscISVBaseObject* self,  PyObject* args, PyObject* kwargs) {
+
+  BOB_TRY
+  
+  // get list of arguments
+  char** kwlist = save.kwlist(0);  
+  PyBobIoHDF5FileObject* hdf5;
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&", kwlist, PyBobIoHDF5File_Converter, &hdf5)) return 0;
+
+  auto hdf5_ = make_safe(hdf5);
+  self->cxx->save(*hdf5->f);
+
+  BOB_CATCH_MEMBER("cannot save the data", 0)
+  Py_RETURN_NONE;
+}
+
+/*** load ***/
+static auto load = bob::extension::FunctionDoc(
+  "load",
+  "Load the configuration of the ISVBase to a given HDF5 file"
+)
+.add_prototype("hdf5")
+.add_parameter("hdf5", ":py:class:`bob.io.base.HDF5File`", "An HDF5 file open for reading");
+static PyObject* PyBobLearnMiscISVBase_Load(PyBobLearnMiscISVBaseObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+  
+  char** kwlist = load.kwlist(0);  
+  PyBobIoHDF5FileObject* hdf5;
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&", kwlist, PyBobIoHDF5File_Converter, &hdf5)) return 0;
+  
+  auto hdf5_ = make_safe(hdf5);  
+  self->cxx->load(*hdf5->f);
+
+  BOB_CATCH_MEMBER("cannot load the data", 0)
+  Py_RETURN_NONE;
+}
+
+
+/*** is_similar_to ***/
+static auto is_similar_to = bob::extension::FunctionDoc(
+  "is_similar_to",
+  
+  "Compares this ISVBase with the ``other`` one to be approximately the same.",
+  "The optional values ``r_epsilon`` and ``a_epsilon`` refer to the "
+  "relative and absolute precision for the ``weights``, ``biases`` "
+  "and any other values internal to this machine."
+)
+.add_prototype("other, [r_epsilon], [a_epsilon]","output")
+.add_parameter("other", ":py:class:`bob.learn.misc.ISVBase`", "A ISVBase object to be compared.")
+.add_parameter("r_epsilon", "float", "Relative precision.")
+.add_parameter("a_epsilon", "float", "Absolute precision.")
+.add_return("output","bool","True if it is similar, otherwise false.");
+static PyObject* PyBobLearnMiscISVBase_IsSimilarTo(PyBobLearnMiscISVBaseObject* self, PyObject* args, PyObject* kwds) {
+
+  /* Parses input arguments in a single shot */
+  char** kwlist = is_similar_to.kwlist(0);
+
+  //PyObject* other = 0;
+  PyBobLearnMiscISVBaseObject* other = 0;
+  double r_epsilon = 1.e-5;
+  double a_epsilon = 1.e-8;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwds, "O!|dd", kwlist,
+        &PyBobLearnMiscISVBase_Type, &other,
+        &r_epsilon, &a_epsilon)){
+
+        is_similar_to.print_usage(); 
+        return 0;        
+  }
+
+  if (self->cxx->is_similar_to(*other->cxx, r_epsilon, a_epsilon))
+    Py_RETURN_TRUE;
+  else
+    Py_RETURN_FALSE;
+}
+
+
+/*** resize ***/
+static auto resize = bob::extension::FunctionDoc(
+  "resize",
+  "Resets the dimensionality of the subspace U. "
+  "U is hence uninitialized.",
+  0,
+  true
+)
+.add_prototype("rU")
+.add_parameter("rU", "int", "Size of U (Within client variation matrix)");
+static PyObject* PyBobLearnMiscISVBase_resize(PyBobLearnMiscISVBaseObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+
+  /* Parses input arguments in a single shot */
+  char** kwlist = resize.kwlist(0);
+
+  int rU = 0;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "i", kwlist, &rU)) Py_RETURN_NONE;
+
+  if (rU <= 0){
+    PyErr_Format(PyExc_TypeError, "rU must be greater than zero");
+    resize.print_usage();
+    return 0;
+  }
+
+  self->cxx->resize(rU);
+
+  BOB_CATCH_MEMBER("cannot perform the resize method", 0)
+
+  Py_RETURN_NONE;
+}
+
+
+
+
+static PyMethodDef PyBobLearnMiscISVBase_methods[] = {
+  {
+    save.name(),
+    (PyCFunction)PyBobLearnMiscISVBase_Save,
+    METH_VARARGS|METH_KEYWORDS,
+    save.doc()
+  },
+  {
+    load.name(),
+    (PyCFunction)PyBobLearnMiscISVBase_Load,
+    METH_VARARGS|METH_KEYWORDS,
+    load.doc()
+  },
+  {
+    is_similar_to.name(),
+    (PyCFunction)PyBobLearnMiscISVBase_IsSimilarTo,
+    METH_VARARGS|METH_KEYWORDS,
+    is_similar_to.doc()
+  },
+  {
+    resize.name(),
+    (PyCFunction)PyBobLearnMiscISVBase_resize,
+    METH_VARARGS|METH_KEYWORDS,
+    resize.doc()
+  },
+  
+  {0} /* Sentinel */
+};
+
+
+/******************************************************************/
+/************ Module Section **************************************/
+/******************************************************************/
+
+// Define the ISV type struct; will be initialized later
+PyTypeObject PyBobLearnMiscISVBase_Type = {
+  PyVarObject_HEAD_INIT(0,0)
+  0
+};
+
+bool init_BobLearnMiscISVBase(PyObject* module)
+{
+  // initialize the type struct
+  PyBobLearnMiscISVBase_Type.tp_name      = ISVBase_doc.name();
+  PyBobLearnMiscISVBase_Type.tp_basicsize = sizeof(PyBobLearnMiscISVBaseObject);
+  PyBobLearnMiscISVBase_Type.tp_flags     = Py_TPFLAGS_DEFAULT;
+  PyBobLearnMiscISVBase_Type.tp_doc       = ISVBase_doc.doc();
+
+  // set the functions
+  PyBobLearnMiscISVBase_Type.tp_new         = PyType_GenericNew;
+  PyBobLearnMiscISVBase_Type.tp_init        = reinterpret_cast<initproc>(PyBobLearnMiscISVBase_init);
+  PyBobLearnMiscISVBase_Type.tp_dealloc     = reinterpret_cast<destructor>(PyBobLearnMiscISVBase_delete);
+  PyBobLearnMiscISVBase_Type.tp_richcompare = reinterpret_cast<richcmpfunc>(PyBobLearnMiscISVBase_RichCompare);
+  PyBobLearnMiscISVBase_Type.tp_methods     = PyBobLearnMiscISVBase_methods;
+  PyBobLearnMiscISVBase_Type.tp_getset      = PyBobLearnMiscISVBase_getseters;
+  //PyBobLearnMiscISVBase_Type.tp_call = reinterpret_cast<ternaryfunc>(PyBobLearnMiscISVBase_forward);
+
+
+  // check that everything is fine
+  if (PyType_Ready(&PyBobLearnMiscISVBase_Type) < 0) return false;
+
+  // add the type to the module
+  Py_INCREF(&PyBobLearnMiscISVBase_Type);
+  return PyModule_AddObject(module, "ISVBase", (PyObject*)&PyBobLearnMiscISVBase_Type) >= 0;
+}
+
diff --git a/bob/learn/misc/isv_machine.cpp b/bob/learn/misc/isv_machine.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0d7f0da6e450406b5d0b86edad6221c7342108e4
--- /dev/null
+++ b/bob/learn/misc/isv_machine.cpp
@@ -0,0 +1,604 @@
+/**
+ * @date Wed Jan 28 13:03:15 2015 +0200
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * @brief Python API for bob::learn::em
+ *
+ * Copyright (C) 2011-2014 Idiap Research Institute, Martigny, Switzerland
+ */
+
+#include "main.h"
+
+/******************************************************************/
+/************ Constructor Section *********************************/
+/******************************************************************/
+
+static auto ISVMachine_doc = bob::extension::ClassDoc(
+  BOB_EXT_MODULE_PREFIX ".ISVMachine",
+  "A ISVMachine. An attached :py:class:`bob.learn.misc.ISVBase` should be provided for Joint Factor Analysis. The :py:class:`bob.learn.misc.ISVMachine` carries information about the speaker factors y and z, whereas a :py:class:`bob.learn.misc.JFABase` carries information about the matrices U, V and D."
+  "References: [Vogt2008,McCool2013]",
+  ""
+).add_constructor(
+  bob::extension::FunctionDoc(
+    "__init__",
+    "Constructor. Builds a new ISVMachine",
+    "",
+    true
+  )
+  .add_prototype("isv_base","")
+  .add_prototype("other","")
+  .add_prototype("hdf5","")
+
+  .add_parameter("isv", ":py:class:`bob.learn.misc.ISVBase`", "The ISVBase associated with this machine")
+  .add_parameter("other", ":py:class:`bob.learn.misc.ISVMachine`", "A ISVMachine object to be copied.")
+  .add_parameter("hdf5", ":py:class:`bob.io.base.HDF5File`", "An HDF5 file open for reading")
+
+);
+
+
+static int PyBobLearnMiscISVMachine_init_copy(PyBobLearnMiscISVMachineObject* self, PyObject* args, PyObject* kwargs) {
+
+  char** kwlist = ISVMachine_doc.kwlist(1);
+  PyBobLearnMiscISVMachineObject* o;
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!", kwlist, &PyBobLearnMiscISVMachine_Type, &o)){
+    ISVMachine_doc.print_usage();
+    return -1;
+  }
+
+  self->cxx.reset(new bob::learn::misc::ISVMachine(*o->cxx));
+  return 0;
+}
+
+
+static int PyBobLearnMiscISVMachine_init_hdf5(PyBobLearnMiscISVMachineObject* self, PyObject* args, PyObject* kwargs) {
+
+  char** kwlist = ISVMachine_doc.kwlist(2);
+
+  PyBobIoHDF5FileObject* config = 0;
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&", kwlist, &PyBobIoHDF5File_Converter, &config)){
+    ISVMachine_doc.print_usage();
+    return -1;
+  }
+
+  self->cxx.reset(new bob::learn::misc::ISVMachine(*(config->f)));
+
+  return 0;
+}
+
+
+static int PyBobLearnMiscISVMachine_init_isvbase(PyBobLearnMiscISVMachineObject* self, PyObject* args, PyObject* kwargs) {
+
+  char** kwlist = ISVMachine_doc.kwlist(0);
+  
+  PyBobLearnMiscISVBaseObject* isv_base;
+
+  //Here we have to select which keyword argument to read  
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!", kwlist, &PyBobLearnMiscISVBase_Type, &isv_base)){
+    ISVMachine_doc.print_usage();
+    return -1;
+  }
+  
+  self->cxx.reset(new bob::learn::misc::ISVMachine(isv_base->cxx));
+  return 0;
+}
+
+
+static int PyBobLearnMiscISVMachine_init(PyBobLearnMiscISVMachineObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+
+  // get the number of command line arguments
+  int nargs = (args?PyTuple_Size(args):0) + (kwargs?PyDict_Size(kwargs):0);
+
+  if(nargs == 1){
+    //Reading the input argument
+    PyObject* arg = 0;
+    if (PyTuple_Size(args))
+      arg = PyTuple_GET_ITEM(args, 0);
+    else {
+      PyObject* tmp = PyDict_Values(kwargs);
+      auto tmp_ = make_safe(tmp);
+      arg = PyList_GET_ITEM(tmp, 0);
+    }
+
+    // If the constructor input is Gaussian object
+    if (PyBobLearnMiscISVMachine_Check(arg))
+      return PyBobLearnMiscISVMachine_init_copy(self, args, kwargs);
+    // If the constructor input is a HDF5
+    else if (PyBobIoHDF5File_Check(arg))
+      return PyBobLearnMiscISVMachine_init_hdf5(self, args, kwargs);
+    // If the constructor input is a JFABase Object
+    else
+      return PyBobLearnMiscISVMachine_init_isvbase(self, args, kwargs);
+  }
+  else{
+    PyErr_Format(PyExc_RuntimeError, "number of arguments mismatch - %s requires only 1 argument, but you provided %d (see help)", Py_TYPE(self)->tp_name, nargs);
+    ISVMachine_doc.print_usage();
+    return -1;
+  }
+  
+  BOB_CATCH_MEMBER("cannot create ISVMachine", 0)
+  return 0;
+}
+
+static void PyBobLearnMiscISVMachine_delete(PyBobLearnMiscISVMachineObject* self) {
+  self->cxx.reset();
+  Py_TYPE(self)->tp_free((PyObject*)self);
+}
+
+static PyObject* PyBobLearnMiscISVMachine_RichCompare(PyBobLearnMiscISVMachineObject* self, PyObject* other, int op) {
+  BOB_TRY
+
+  if (!PyBobLearnMiscISVMachine_Check(other)) {
+    PyErr_Format(PyExc_TypeError, "cannot compare `%s' with `%s'", Py_TYPE(self)->tp_name, Py_TYPE(other)->tp_name);
+    return 0;
+  }
+  auto other_ = reinterpret_cast<PyBobLearnMiscISVMachineObject*>(other);
+  switch (op) {
+    case Py_EQ:
+      if (*self->cxx==*other_->cxx) Py_RETURN_TRUE; else Py_RETURN_FALSE;
+    case Py_NE:
+      if (*self->cxx==*other_->cxx) Py_RETURN_FALSE; else Py_RETURN_TRUE;
+    default:
+      Py_INCREF(Py_NotImplemented);
+      return Py_NotImplemented;
+  }
+  BOB_CATCH_MEMBER("cannot compare ISVMachine objects", 0)
+}
+
+int PyBobLearnMiscISVMachine_Check(PyObject* o) {
+  return PyObject_IsInstance(o, reinterpret_cast<PyObject*>(&PyBobLearnMiscISVMachine_Type));
+}
+
+
+/******************************************************************/
+/************ Variables Section ***********************************/
+/******************************************************************/
+
+/***** shape *****/
+static auto shape = bob::extension::VariableDoc(
+  "shape",
+  "(int,int, int, int)",
+  "A tuple that represents the number of gaussians, dimensionality of each Gaussian and dimensionality of the rU (within client variability matrix)) ``(#Gaussians, #Inputs, #rU)``.",
+  ""
+);
+PyObject* PyBobLearnMiscISVMachine_getShape(PyBobLearnMiscISVMachineObject* self, void*) {
+  BOB_TRY
+  return Py_BuildValue("(i,i,i)", self->cxx->getNGaussians(), self->cxx->getNInputs(), self->cxx->getDimRu());
+  BOB_CATCH_MEMBER("shape could not be read", 0)
+}
+
+/***** supervector_length *****/
+static auto supervector_length = bob::extension::VariableDoc(
+  "supervector_length",
+  "int",
+
+  "Returns the supervector length."
+  "NGaussians x NInputs: Number of Gaussian components by the feature dimensionality",
+  
+  "@warning An exception is thrown if no Universal Background Model has been set yet."
+);
+PyObject* PyBobLearnMiscISVMachine_getSupervectorLength(PyBobLearnMiscISVMachineObject* self, void*) {
+  BOB_TRY
+  return Py_BuildValue("i", self->cxx->getSupervectorLength());
+  BOB_CATCH_MEMBER("supervector_length could not be read", 0)
+}
+
+/***** z *****/
+static auto Z = bob::extension::VariableDoc(
+  "z",
+  "array_like <float, 1D>",
+  "Returns the z speaker factor. Eq (31) from [McCool2013]",
+  ""
+);
+PyObject* PyBobLearnMiscISVMachine_getZ(PyBobLearnMiscISVMachineObject* self, void*){
+  BOB_TRY
+  return PyBlitzArrayCxx_AsConstNumpy(self->cxx->getZ());
+  BOB_CATCH_MEMBER("`z` could not be read", 0)
+}
+int PyBobLearnMiscISVMachine_setZ(PyBobLearnMiscISVMachineObject* self, PyObject* value, void*){
+  BOB_TRY
+  PyBlitzArrayObject* o;
+  if (!PyBlitzArray_Converter(value, &o)){
+    PyErr_Format(PyExc_RuntimeError, "%s %s expects a 1D array of floats", Py_TYPE(self)->tp_name, Z.name());
+    return -1;
+  }
+  auto o_ = make_safe(o);
+  auto b = PyBlitzArrayCxx_AsBlitz<double,1>(o, "z");
+  if (!b) return -1;
+  self->cxx->setZ(*b);
+  return 0;
+  BOB_CATCH_MEMBER("`z` vector could not be set", -1)
+}
+
+
+/***** x *****/
+static auto X = bob::extension::VariableDoc(
+  "x",
+  "array_like <float, 1D>",
+  "Returns the X session factor. Eq (29) from [McCool2013]",
+  "The latent variable x (last one computed). This is a feature provided for convenience, but this attribute is not 'part' of the machine. The session latent variable x is indeed not class-specific, but depends on the sample considered. Furthermore, it is not saved into the machine or used when comparing machines."
+);
+PyObject* PyBobLearnMiscISVMachine_getX(PyBobLearnMiscISVMachineObject* self, void*){
+  BOB_TRY
+  return PyBlitzArrayCxx_AsConstNumpy(self->cxx->getX());
+  BOB_CATCH_MEMBER("`x` could not be read", 0)
+}
+
+
+/***** isv_base *****/
+static auto isv_base = bob::extension::VariableDoc(
+  "isv_base",
+  ":py:class:`bob.learn.misc.ISVBase`",
+  "The ISVBase attached to this machine",
+  ""
+);
+PyObject* PyBobLearnMiscISVMachine_getISVBase(PyBobLearnMiscISVMachineObject* self, void*){
+  BOB_TRY
+
+  boost::shared_ptr<bob::learn::misc::ISVBase> isv_base_o = self->cxx->getISVBase();
+
+  //Allocating the correspondent python object
+  PyBobLearnMiscISVBaseObject* retval =
+    (PyBobLearnMiscISVBaseObject*)PyBobLearnMiscISVBase_Type.tp_alloc(&PyBobLearnMiscISVBase_Type, 0);
+  retval->cxx = isv_base_o;
+
+  return Py_BuildValue("O",retval);
+  BOB_CATCH_MEMBER("isv_base could not be read", 0)
+}
+int PyBobLearnMiscISVMachine_setISVBase(PyBobLearnMiscISVMachineObject* self, PyObject* value, void*){
+  BOB_TRY
+
+  if (!PyBobLearnMiscISVBase_Check(value)){
+    PyErr_Format(PyExc_RuntimeError, "%s %s expects a :py:class:`bob.learn.misc.ISVBase`", Py_TYPE(self)->tp_name, isv_base.name());
+    return -1;
+  }
+
+  PyBobLearnMiscISVBaseObject* isv_base_o = 0;
+  PyArg_Parse(value, "O!", &PyBobLearnMiscISVBase_Type,&isv_base_o);
+
+  self->cxx->setISVBase(isv_base_o->cxx);
+
+  return 0;
+  BOB_CATCH_MEMBER("isv_base could not be set", -1)  
+}
+
+
+
+
+static PyGetSetDef PyBobLearnMiscISVMachine_getseters[] = { 
+  {
+   shape.name(),
+   (getter)PyBobLearnMiscISVMachine_getShape,
+   0,
+   shape.doc(),
+   0
+  },
+  
+  {
+   supervector_length.name(),
+   (getter)PyBobLearnMiscISVMachine_getSupervectorLength,
+   0,
+   supervector_length.doc(),
+   0
+  },
+  
+  {
+   isv_base.name(),
+   (getter)PyBobLearnMiscISVMachine_getISVBase,
+   (setter)PyBobLearnMiscISVMachine_setISVBase,
+   isv_base.doc(),
+   0
+  },
+
+  {
+   Z.name(),
+   (getter)PyBobLearnMiscISVMachine_getZ,
+   (setter)PyBobLearnMiscISVMachine_setZ,
+   Z.doc(),
+   0
+  },
+
+  {
+   X.name(),
+   (getter)PyBobLearnMiscISVMachine_getX,
+   0,
+   X.doc(),
+   0
+  },
+
+
+  {0}  // Sentinel
+};
+
+
+/******************************************************************/
+/************ Functions Section ***********************************/
+/******************************************************************/
+
+
+/*** save ***/
+static auto save = bob::extension::FunctionDoc(
+  "save",
+  "Save the configuration of the ISVMachine to a given HDF5 file"
+)
+.add_prototype("hdf5")
+.add_parameter("hdf5", ":py:class:`bob.io.base.HDF5File`", "An HDF5 file open for writing");
+static PyObject* PyBobLearnMiscISVMachine_Save(PyBobLearnMiscISVMachineObject* self,  PyObject* args, PyObject* kwargs) {
+
+  BOB_TRY
+  
+  // get list of arguments
+  char** kwlist = save.kwlist(0);  
+  PyBobIoHDF5FileObject* hdf5;
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&", kwlist, PyBobIoHDF5File_Converter, &hdf5)) return 0;
+
+  auto hdf5_ = make_safe(hdf5);
+  self->cxx->save(*hdf5->f);
+
+  BOB_CATCH_MEMBER("cannot save the data", 0)
+  Py_RETURN_NONE;
+}
+
+/*** load ***/
+static auto load = bob::extension::FunctionDoc(
+  "load",
+  "Load the configuration of the ISVMachine to a given HDF5 file"
+)
+.add_prototype("hdf5")
+.add_parameter("hdf5", ":py:class:`bob.io.base.HDF5File`", "An HDF5 file open for reading");
+static PyObject* PyBobLearnMiscISVMachine_Load(PyBobLearnMiscISVMachineObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+  
+  char** kwlist = load.kwlist(0);  
+  PyBobIoHDF5FileObject* hdf5;
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&", kwlist, PyBobIoHDF5File_Converter, &hdf5)) return 0;
+  
+  auto hdf5_ = make_safe(hdf5);  
+  self->cxx->load(*hdf5->f);
+
+  BOB_CATCH_MEMBER("cannot load the data", 0)
+  Py_RETURN_NONE;
+}
+
+
+/*** is_similar_to ***/
+static auto is_similar_to = bob::extension::FunctionDoc(
+  "is_similar_to",
+  
+  "Compares this ISVMachine with the ``other`` one to be approximately the same.",
+  "The optional values ``r_epsilon`` and ``a_epsilon`` refer to the "
+  "relative and absolute precision for the ``weights``, ``biases`` "
+  "and any other values internal to this machine."
+)
+.add_prototype("other, [r_epsilon], [a_epsilon]","output")
+.add_parameter("other", ":py:class:`bob.learn.misc.ISVMachine`", "A ISVMachine object to be compared.")
+.add_parameter("r_epsilon", "float", "Relative precision.")
+.add_parameter("a_epsilon", "float", "Absolute precision.")
+.add_return("output","bool","True if it is similar, otherwise false.");
+static PyObject* PyBobLearnMiscISVMachine_IsSimilarTo(PyBobLearnMiscISVMachineObject* self, PyObject* args, PyObject* kwds) {
+
+  /* Parses input arguments in a single shot */
+  char** kwlist = is_similar_to.kwlist(0);
+
+  //PyObject* other = 0;
+  PyBobLearnMiscISVMachineObject* other = 0;
+  double r_epsilon = 1.e-5;
+  double a_epsilon = 1.e-8;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwds, "O!|dd", kwlist,
+        &PyBobLearnMiscISVMachine_Type, &other,
+        &r_epsilon, &a_epsilon)){
+
+        is_similar_to.print_usage(); 
+        return 0;        
+  }
+
+  if (self->cxx->is_similar_to(*other->cxx, r_epsilon, a_epsilon))
+    Py_RETURN_TRUE;
+  else
+    Py_RETURN_FALSE;
+}
+
+
+/*** estimate_x ***/
+static auto estimate_x = bob::extension::FunctionDoc(
+  "estimate_x",
+  "Estimates the session offset x (LPT assumption) given GMM statistics.",
+  "Estimates x from the GMM statistics considering the LPT assumption, that is the latent session variable x is approximated using the UBM", 
+  true
+)
+.add_prototype("stats,input")
+.add_parameter("stats", ":py:class:`bob.learn.misc.GMMStats`", "Statistics of the GMM")
+.add_parameter("input", "array_like <float, 1D>", "Input vector");
+static PyObject* PyBobLearnMiscISVMachine_estimateX(PyBobLearnMiscISVMachineObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+
+  char** kwlist = estimate_x.kwlist(0);
+
+  PyBobLearnMiscGMMStatsObject* stats = 0;
+  PyBlitzArrayObject* input           = 0;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!O&", kwlist, &PyBobLearnMiscGMMStats_Type, &stats, 
+                                                                 &PyBlitzArray_Converter,&input))
+    Py_RETURN_NONE;
+
+  //protects acquired resources through this scope
+  auto input_ = make_safe(input);
+  self->cxx->estimateX(*stats->cxx, *PyBlitzArrayCxx_AsBlitz<double,1>(input));
+
+  BOB_CATCH_MEMBER("cannot estimate X", 0)
+  Py_RETURN_NONE;
+}
+
+
+/*** estimate_ux ***/
+static auto estimate_ux = bob::extension::FunctionDoc(
+  "estimate_ux",
+  "Estimates Ux (LPT assumption) given GMM statistics.",
+  "Estimates Ux from the GMM statistics considering the LPT assumption, that is the latent session variable x is approximated using the UBM.", 
+  true
+)
+.add_prototype("stats,input")
+.add_parameter("stats", ":py:class:`bob.learn.misc.GMMStats`", "Statistics of the GMM")
+.add_parameter("input", "array_like <float, 1D>", "Input vector");
+static PyObject* PyBobLearnMiscISVMachine_estimateUx(PyBobLearnMiscISVMachineObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+
+  char** kwlist = estimate_ux.kwlist(0);
+
+  PyBobLearnMiscGMMStatsObject* stats = 0;
+  PyBlitzArrayObject* input           = 0;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!O&", kwlist, &PyBobLearnMiscGMMStats_Type, &stats, 
+                                                                 &PyBlitzArray_Converter,&input))
+    Py_RETURN_NONE;
+
+  //protects acquired resources through this scope
+  auto input_ = make_safe(input);
+  self->cxx->estimateUx(*stats->cxx, *PyBlitzArrayCxx_AsBlitz<double,1>(input));
+
+  BOB_CATCH_MEMBER("cannot estimate Ux", 0)
+  Py_RETURN_NONE;
+}
+
+
+/*** forward_ux ***/
+static auto forward_ux = bob::extension::FunctionDoc(
+  "forward_ux",
+  "Computes a score for the given UBM statistics and given the Ux vector",
+  "", 
+  true
+)
+.add_prototype("stats,ux")
+.add_parameter("stats", ":py:class:`bob.learn.misc.GMMStats`", "Statistics as input")
+.add_parameter("ux", "array_like <float, 1D>", "Input vector");
+static PyObject* PyBobLearnMiscISVMachine_ForwardUx(PyBobLearnMiscISVMachineObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+
+  char** kwlist = forward_ux.kwlist(0);
+
+  PyBobLearnMiscGMMStatsObject* stats = 0;
+  PyBlitzArrayObject* ux_input        = 0;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!O&", kwlist, &PyBobLearnMiscGMMStats_Type, &stats, 
+                                                                 &PyBlitzArray_Converter,&ux_input))
+    Py_RETURN_NONE;
+
+  //protects acquired resources through this scope
+  auto ux_input_ = make_safe(ux_input);
+  double score = self->cxx->forward(*stats->cxx, *PyBlitzArrayCxx_AsBlitz<double,1>(ux_input));
+  
+  return Py_BuildValue("d", score);
+  BOB_CATCH_MEMBER("cannot forward_ux", 0)
+}
+
+
+/*** forward ***/
+static auto forward = bob::extension::FunctionDoc(
+  "forward",
+  "Execute the machine",
+  "", 
+  true
+)
+.add_prototype("stats")
+.add_parameter("stats", ":py:class:`bob.learn.misc.GMMStats`", "Statistics as input");
+static PyObject* PyBobLearnMiscISVMachine_Forward(PyBobLearnMiscISVMachineObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+
+  char** kwlist = forward.kwlist(0);
+
+  PyBobLearnMiscGMMStatsObject* stats = 0;
+  
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!", kwlist, &PyBobLearnMiscGMMStats_Type, &stats))
+    Py_RETURN_NONE;
+
+  //protects acquired resources through this scope
+  double score = self->cxx->forward(*stats->cxx);
+
+  return Py_BuildValue("d", score);
+  BOB_CATCH_MEMBER("cannot forward", 0)
+
+}
+
+
+static PyMethodDef PyBobLearnMiscISVMachine_methods[] = {
+  {
+    save.name(),
+    (PyCFunction)PyBobLearnMiscISVMachine_Save,
+    METH_VARARGS|METH_KEYWORDS,
+    save.doc()
+  },
+  {
+    load.name(),
+    (PyCFunction)PyBobLearnMiscISVMachine_Load,
+    METH_VARARGS|METH_KEYWORDS,
+    load.doc()
+  },
+  {
+    is_similar_to.name(),
+    (PyCFunction)PyBobLearnMiscISVMachine_IsSimilarTo,
+    METH_VARARGS|METH_KEYWORDS,
+    is_similar_to.doc()
+  },
+  
+  {
+    estimate_x.name(),
+    (PyCFunction)PyBobLearnMiscISVMachine_estimateX,
+    METH_VARARGS|METH_KEYWORDS,
+    estimate_x.doc()
+  },
+  
+  {
+    estimate_ux.name(),
+    (PyCFunction)PyBobLearnMiscISVMachine_estimateUx,
+    METH_VARARGS|METH_KEYWORDS,
+    estimate_ux.doc()
+  },
+
+  {
+    forward_ux.name(),
+    (PyCFunction)PyBobLearnMiscISVMachine_ForwardUx,
+    METH_VARARGS|METH_KEYWORDS,
+    forward_ux.doc()
+  },
+
+  {0} /* Sentinel */
+};
+
+
+/******************************************************************/
+/************ Module Section **************************************/
+/******************************************************************/
+
+// Define the JFA type struct; will be initialized later
+PyTypeObject PyBobLearnMiscISVMachine_Type = {
+  PyVarObject_HEAD_INIT(0,0)
+  0
+};
+
+bool init_BobLearnMiscISVMachine(PyObject* module)
+{
+  // initialize the type struct
+  PyBobLearnMiscISVMachine_Type.tp_name      = ISVMachine_doc.name();
+  PyBobLearnMiscISVMachine_Type.tp_basicsize = sizeof(PyBobLearnMiscISVMachineObject);
+  PyBobLearnMiscISVMachine_Type.tp_flags     = Py_TPFLAGS_DEFAULT;
+  PyBobLearnMiscISVMachine_Type.tp_doc       = ISVMachine_doc.doc();
+
+  // set the functions
+  PyBobLearnMiscISVMachine_Type.tp_new         = PyType_GenericNew;
+  PyBobLearnMiscISVMachine_Type.tp_init        = reinterpret_cast<initproc>(PyBobLearnMiscISVMachine_init);
+  PyBobLearnMiscISVMachine_Type.tp_dealloc     = reinterpret_cast<destructor>(PyBobLearnMiscISVMachine_delete);
+  PyBobLearnMiscISVMachine_Type.tp_richcompare = reinterpret_cast<richcmpfunc>(PyBobLearnMiscISVMachine_RichCompare);
+  PyBobLearnMiscISVMachine_Type.tp_methods     = PyBobLearnMiscISVMachine_methods;
+  PyBobLearnMiscISVMachine_Type.tp_getset      = PyBobLearnMiscISVMachine_getseters;
+  PyBobLearnMiscISVMachine_Type.tp_call = reinterpret_cast<ternaryfunc>(PyBobLearnMiscISVMachine_Forward);
+
+
+  // check that everything is fine
+  if (PyType_Ready(&PyBobLearnMiscISVMachine_Type) < 0) return false;
+
+  // add the type to the module
+  Py_INCREF(&PyBobLearnMiscISVMachine_Type);
+  return PyModule_AddObject(module, "ISVMachine", (PyObject*)&PyBobLearnMiscISVMachine_Type) >= 0;
+}
+
diff --git a/bob/learn/misc/jfa_base.cpp b/bob/learn/misc/jfa_base.cpp
index 7fa60b5dcc9684359acabb837a61facee30a4fb3..d995114e65cee760d96697e5645903243a644e05 100644
--- a/bob/learn/misc/jfa_base.cpp
+++ b/bob/learn/misc/jfa_base.cpp
@@ -1,5 +1,5 @@
 /**
- * @date Tue Jan 27 17:03:15 2015 +0200
+ * @date Wed Jan 27 17:03:15 2015 +0200
  * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
  *
  * @brief Python API for bob::learn::em
@@ -15,14 +15,13 @@
 
 static auto JFABase_doc = bob::extension::ClassDoc(
   BOB_EXT_MODULE_PREFIX ".JFABase",
-
-  "Constructor. Builds a new JFABase. "
-  "TODO: add a reference to the journal articles",
+  "A JFABase instance can be seen as a container for U, V and D when performing Joint Factor Analysis (JFA)."
+  "References: [Vogt2008,McCool2013]",
   ""
 ).add_constructor(
   bob::extension::FunctionDoc(
     "__init__",
-    "Creates a FABase",
+    "Constructor. Builds a new JFABase",
     "",
     true
   )
@@ -462,35 +461,36 @@ static PyObject* PyBobLearnMiscJFABase_IsSimilarTo(PyBobLearnMiscJFABaseObject*
 /*** resize ***/
 static auto resize = bob::extension::FunctionDoc(
   "resize",
-  "Allocates space for the statistics and resets to zero.",
+  "Resets the dimensionality of the subspace U and V. "
+  "U and V are hence uninitialized",
   0,
   true
 )
-.add_prototype("n_gaussians,n_inputs")
-.add_parameter("n_gaussians", "int", "Number of gaussians")
-.add_parameter("n_inputs", "int", "Dimensionality of the feature vector");
+.add_prototype("rU,rV")
+.add_parameter("rU", "int", "Size of U (Within client variation matrix)")
+.add_parameter("rV", "int", "Size of V (Between client variation matrix)");
 static PyObject* PyBobLearnMiscJFABase_resize(PyBobLearnMiscJFABaseObject* self, PyObject* args, PyObject* kwargs) {
   BOB_TRY
 
   /* Parses input arguments in a single shot */
   char** kwlist = resize.kwlist(0);
 
-  int n_gaussians = 0;
-  int n_inputs = 0;
+  int rU = 0;
+  int rV = 0;
 
-  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "ii", kwlist, &n_gaussians, &n_inputs)) Py_RETURN_NONE;
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "ii", kwlist, &rU, &rV)) Py_RETURN_NONE;
 
-  if (n_gaussians <= 0){
-    PyErr_Format(PyExc_TypeError, "n_gaussians must be greater than zero");
+  if (rU <= 0){
+    PyErr_Format(PyExc_TypeError, "rU must be greater than zero");
     resize.print_usage();
     return 0;
   }
-  if (n_inputs <= 0){
-    PyErr_Format(PyExc_TypeError, "n_inputs must be greater than zero");
+  if (rV <= 0){
+    PyErr_Format(PyExc_TypeError, "rV must be greater than zero");
     resize.print_usage();
     return 0;
   }
-  self->cxx->resize(n_gaussians, n_inputs);
+  self->cxx->resize(rU, rV);
 
   BOB_CATCH_MEMBER("cannot perform the resize method", 0)
 
diff --git a/bob/learn/misc/jfa_machine.cpp b/bob/learn/misc/jfa_machine.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..af714125e44421599934534422f4eca07bab120e
--- /dev/null
+++ b/bob/learn/misc/jfa_machine.cpp
@@ -0,0 +1,650 @@
+/**
+ * @date Wed Jan 28 17:03:15 2015 +0200
+ * @author Tiago de Freitas Pereira <tiago.pereira@idiap.ch>
+ *
+ * @brief Python API for bob::learn::em
+ *
+ * Copyright (C) 2011-2014 Idiap Research Institute, Martigny, Switzerland
+ */
+
+#include "main.h"
+
+/******************************************************************/
+/************ Constructor Section *********************************/
+/******************************************************************/
+
+static auto JFAMachine_doc = bob::extension::ClassDoc(
+  BOB_EXT_MODULE_PREFIX ".JFAMachine",
+  "A JFAMachine. An attached :py:class:`bob.learn.misc.JFABase` should be provided for Joint Factor Analysis. The :py:class:`bob.learn.misc.JFAMachine` carries information about the speaker factors y and z, whereas a :py:class:`bob.learn.misc.JFABase` carries information about the matrices U, V and D."
+  "References: [Vogt2008,McCool2013]",
+  ""
+).add_constructor(
+  bob::extension::FunctionDoc(
+    "__init__",
+    "Constructor. Builds a new JFAMachine",
+    "",
+    true
+  )
+  .add_prototype("jfa_base","")
+  .add_prototype("other","")
+  .add_prototype("hdf5","")
+
+  .add_parameter("jfa", ":py:class:`bob.learn.misc.JFABase`", "The JFABase associated with this machine")
+  .add_parameter("other", ":py:class:`bob.learn.misc.JFAMachine`", "A JFAMachine object to be copied.")
+  .add_parameter("hdf5", ":py:class:`bob.io.base.HDF5File`", "An HDF5 file open for reading")
+
+);
+
+
+static int PyBobLearnMiscJFAMachine_init_copy(PyBobLearnMiscJFAMachineObject* self, PyObject* args, PyObject* kwargs) {
+
+  char** kwlist = JFAMachine_doc.kwlist(1);
+  PyBobLearnMiscJFAMachineObject* o;
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!", kwlist, &PyBobLearnMiscJFAMachine_Type, &o)){
+    JFAMachine_doc.print_usage();
+    return -1;
+  }
+
+  self->cxx.reset(new bob::learn::misc::JFAMachine(*o->cxx));
+  return 0;
+}
+
+
+static int PyBobLearnMiscJFAMachine_init_hdf5(PyBobLearnMiscJFAMachineObject* self, PyObject* args, PyObject* kwargs) {
+
+  char** kwlist = JFAMachine_doc.kwlist(2);
+
+  PyBobIoHDF5FileObject* config = 0;
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&", kwlist, &PyBobIoHDF5File_Converter, &config)){
+    JFAMachine_doc.print_usage();
+    return -1;
+  }
+
+  self->cxx.reset(new bob::learn::misc::JFAMachine(*(config->f)));
+
+  return 0;
+}
+
+
+static int PyBobLearnMiscJFAMachine_init_jfabase(PyBobLearnMiscJFAMachineObject* self, PyObject* args, PyObject* kwargs) {
+
+  char** kwlist = JFAMachine_doc.kwlist(0);
+  
+  PyBobLearnMiscJFABaseObject* jfa_base;
+
+  //Here we have to select which keyword argument to read  
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!", kwlist, &PyBobLearnMiscJFABase_Type, &jfa_base)){
+    JFAMachine_doc.print_usage();
+    return -1;
+  }
+  
+  self->cxx.reset(new bob::learn::misc::JFAMachine(jfa_base->cxx));
+  return 0;
+}
+
+
+static int PyBobLearnMiscJFAMachine_init(PyBobLearnMiscJFAMachineObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+
+  // get the number of command line arguments
+  int nargs = (args?PyTuple_Size(args):0) + (kwargs?PyDict_Size(kwargs):0);
+
+  if(nargs == 1){
+    //Reading the input argument
+    PyObject* arg = 0;
+    if (PyTuple_Size(args))
+      arg = PyTuple_GET_ITEM(args, 0);
+    else {
+      PyObject* tmp = PyDict_Values(kwargs);
+      auto tmp_ = make_safe(tmp);
+      arg = PyList_GET_ITEM(tmp, 0);
+    }
+
+    // If the constructor input is Gaussian object
+    if (PyBobLearnMiscJFAMachine_Check(arg))
+      return PyBobLearnMiscJFAMachine_init_copy(self, args, kwargs);
+    // If the constructor input is a HDF5
+    else if (PyBobIoHDF5File_Check(arg))
+      return PyBobLearnMiscJFAMachine_init_hdf5(self, args, kwargs);
+    // If the constructor input is a JFABase Object
+    else
+      return PyBobLearnMiscJFAMachine_init_jfabase(self, args, kwargs);
+  }
+  else{
+    PyErr_Format(PyExc_RuntimeError, "number of arguments mismatch - %s requires only 1 argument, but you provided %d (see help)", Py_TYPE(self)->tp_name, nargs);
+    JFAMachine_doc.print_usage();
+    return -1;
+  }
+  
+  BOB_CATCH_MEMBER("cannot create JFAMachine", 0)
+  return 0;
+}
+
+static void PyBobLearnMiscJFAMachine_delete(PyBobLearnMiscJFAMachineObject* self) {
+  self->cxx.reset();
+  Py_TYPE(self)->tp_free((PyObject*)self);
+}
+
+static PyObject* PyBobLearnMiscJFAMachine_RichCompare(PyBobLearnMiscJFAMachineObject* self, PyObject* other, int op) {
+  BOB_TRY
+
+  if (!PyBobLearnMiscJFAMachine_Check(other)) {
+    PyErr_Format(PyExc_TypeError, "cannot compare `%s' with `%s'", Py_TYPE(self)->tp_name, Py_TYPE(other)->tp_name);
+    return 0;
+  }
+  auto other_ = reinterpret_cast<PyBobLearnMiscJFAMachineObject*>(other);
+  switch (op) {
+    case Py_EQ:
+      if (*self->cxx==*other_->cxx) Py_RETURN_TRUE; else Py_RETURN_FALSE;
+    case Py_NE:
+      if (*self->cxx==*other_->cxx) Py_RETURN_FALSE; else Py_RETURN_TRUE;
+    default:
+      Py_INCREF(Py_NotImplemented);
+      return Py_NotImplemented;
+  }
+  BOB_CATCH_MEMBER("cannot compare JFAMachine objects", 0)
+}
+
+int PyBobLearnMiscJFAMachine_Check(PyObject* o) {
+  return PyObject_IsInstance(o, reinterpret_cast<PyObject*>(&PyBobLearnMiscJFAMachine_Type));
+}
+
+
+/******************************************************************/
+/************ Variables Section ***********************************/
+/******************************************************************/
+
+/***** shape *****/
+static auto shape = bob::extension::VariableDoc(
+  "shape",
+  "(int,int, int, int)",
+  "A tuple that represents the number of gaussians, dimensionality of each Gaussian, dimensionality of the rU (within client variability matrix) and dimensionality of the rV (between client variability matrix) ``(#Gaussians, #Inputs, #rU, #rV)``.",
+  ""
+);
+PyObject* PyBobLearnMiscJFAMachine_getShape(PyBobLearnMiscJFAMachineObject* self, void*) {
+  BOB_TRY
+  return Py_BuildValue("(i,i,i,i)", self->cxx->getNGaussians(), self->cxx->getNInputs(), self->cxx->getDimRu(), self->cxx->getDimRv());
+  BOB_CATCH_MEMBER("shape could not be read", 0)
+}
+
+/***** supervector_length *****/
+static auto supervector_length = bob::extension::VariableDoc(
+  "supervector_length",
+  "int",
+
+  "Returns the supervector length."
+  "NGaussians x NInputs: Number of Gaussian components by the feature dimensionality",
+  
+  "@warning An exception is thrown if no Universal Background Model has been set yet."
+);
+PyObject* PyBobLearnMiscJFAMachine_getSupervectorLength(PyBobLearnMiscJFAMachineObject* self, void*) {
+  BOB_TRY
+  return Py_BuildValue("i", self->cxx->getSupervectorLength());
+  BOB_CATCH_MEMBER("supervector_length could not be read", 0)
+}
+
+
+/***** y *****/
+static auto Y = bob::extension::VariableDoc(
+  "y",
+  "array_like <float, 1D>",
+  "Returns the y speaker factor. Eq (30) from [McCool2013]",
+  ""
+);
+PyObject* PyBobLearnMiscJFAMachine_getY(PyBobLearnMiscJFAMachineObject* self, void*){
+  BOB_TRY
+  return PyBlitzArrayCxx_AsConstNumpy(self->cxx->getY());
+  BOB_CATCH_MEMBER("`y` could not be read", 0)
+}
+int PyBobLearnMiscJFAMachine_setY(PyBobLearnMiscJFAMachineObject* self, PyObject* value, void*){
+  BOB_TRY
+  PyBlitzArrayObject* o;
+  if (!PyBlitzArray_Converter(value, &o)){
+    PyErr_Format(PyExc_RuntimeError, "%s %s expects a 1D array of floats", Py_TYPE(self)->tp_name, Y.name());
+    return -1;
+  }
+  auto o_ = make_safe(o);
+  auto b = PyBlitzArrayCxx_AsBlitz<double,1>(o, "y");
+  if (!b) return -1;
+  self->cxx->setY(*b);
+  return 0;
+  BOB_CATCH_MEMBER("`y` vector could not be set", -1)
+}
+
+
+/***** z *****/
+static auto Z = bob::extension::VariableDoc(
+  "z",
+  "array_like <float, 1D>",
+  "Returns the z speaker factor. Eq (31) from [McCool2013]",
+  ""
+);
+PyObject* PyBobLearnMiscJFAMachine_getZ(PyBobLearnMiscJFAMachineObject* self, void*){
+  BOB_TRY
+  return PyBlitzArrayCxx_AsConstNumpy(self->cxx->getZ());
+  BOB_CATCH_MEMBER("`z` could not be read", 0)
+}
+int PyBobLearnMiscJFAMachine_setZ(PyBobLearnMiscJFAMachineObject* self, PyObject* value, void*){
+  BOB_TRY
+  PyBlitzArrayObject* o;
+  if (!PyBlitzArray_Converter(value, &o)){
+    PyErr_Format(PyExc_RuntimeError, "%s %s expects a 1D array of floats", Py_TYPE(self)->tp_name, Z.name());
+    return -1;
+  }
+  auto o_ = make_safe(o);
+  auto b = PyBlitzArrayCxx_AsBlitz<double,1>(o, "z");
+  if (!b) return -1;
+  self->cxx->setZ(*b);
+  return 0;
+  BOB_CATCH_MEMBER("`z` vector could not be set", -1)
+}
+
+
+/***** x *****/
+static auto X = bob::extension::VariableDoc(
+  "x",
+  "array_like <float, 1D>",
+  "Returns the X session factor. Eq (29) from [McCool2013]",
+  "The latent variable x (last one computed). This is a feature provided for convenience, but this attribute is not 'part' of the machine. The session latent variable x is indeed not class-specific, but depends on the sample considered. Furthermore, it is not saved into the machine or used when comparing machines."
+);
+PyObject* PyBobLearnMiscJFAMachine_getX(PyBobLearnMiscJFAMachineObject* self, void*){
+  BOB_TRY
+  return PyBlitzArrayCxx_AsConstNumpy(self->cxx->getX());
+  BOB_CATCH_MEMBER("`x` could not be read", 0)
+}
+
+
+/***** jfa_base *****/
+static auto jfa_base = bob::extension::VariableDoc(
+  "jfa_base",
+  ":py:class:`bob.learn.misc.JFABase`",
+  "The JFABase attached to this machine",
+  ""
+);
+PyObject* PyBobLearnMiscJFAMachine_getJFABase(PyBobLearnMiscJFAMachineObject* self, void*){
+  BOB_TRY
+
+  boost::shared_ptr<bob::learn::misc::JFABase> jfa_base_o = self->cxx->getJFABase();
+
+  //Allocating the correspondent python object
+  PyBobLearnMiscJFABaseObject* retval =
+    (PyBobLearnMiscJFABaseObject*)PyBobLearnMiscJFABase_Type.tp_alloc(&PyBobLearnMiscJFABase_Type, 0);
+  retval->cxx = jfa_base_o;
+
+  return Py_BuildValue("O",retval);
+  BOB_CATCH_MEMBER("jfa_base could not be read", 0)
+}
+int PyBobLearnMiscJFAMachine_setJFABase(PyBobLearnMiscJFAMachineObject* self, PyObject* value, void*){
+  BOB_TRY
+
+  if (!PyBobLearnMiscJFABase_Check(value)){
+    PyErr_Format(PyExc_RuntimeError, "%s %s expects a :py:class:`bob.learn.misc.JFABase`", Py_TYPE(self)->tp_name, jfa_base.name());
+    return -1;
+  }
+
+  PyBobLearnMiscJFABaseObject* jfa_base_o = 0;
+  PyArg_Parse(value, "O!", &PyBobLearnMiscJFABase_Type,&jfa_base_o);
+
+  self->cxx->setJFABase(jfa_base_o->cxx);
+
+  return 0;
+  BOB_CATCH_MEMBER("jfa_base could not be set", -1)  
+}
+
+
+
+
+static PyGetSetDef PyBobLearnMiscJFAMachine_getseters[] = { 
+  {
+   shape.name(),
+   (getter)PyBobLearnMiscJFAMachine_getShape,
+   0,
+   shape.doc(),
+   0
+  },
+  
+  {
+   supervector_length.name(),
+   (getter)PyBobLearnMiscJFAMachine_getSupervectorLength,
+   0,
+   supervector_length.doc(),
+   0
+  },
+  
+  {
+   jfa_base.name(),
+   (getter)PyBobLearnMiscJFAMachine_getJFABase,
+   (setter)PyBobLearnMiscJFAMachine_setJFABase,
+   jfa_base.doc(),
+   0
+  },
+
+  {
+   Y.name(),
+   (getter)PyBobLearnMiscJFAMachine_getY,
+   (setter)PyBobLearnMiscJFAMachine_setY,
+   Y.doc(),
+   0
+  },
+
+  {
+   Z.name(),
+   (getter)PyBobLearnMiscJFAMachine_getZ,
+   (setter)PyBobLearnMiscJFAMachine_setZ,
+   Z.doc(),
+   0
+  },
+
+  {
+   X.name(),
+   (getter)PyBobLearnMiscJFAMachine_getX,
+   0,
+   X.doc(),
+   0
+  },
+
+
+  {0}  // Sentinel
+};
+
+
+/******************************************************************/
+/************ Functions Section ***********************************/
+/******************************************************************/
+
+
+/*** save ***/
+static auto save = bob::extension::FunctionDoc(
+  "save",
+  "Save the configuration of the JFAMachine to a given HDF5 file"
+)
+.add_prototype("hdf5")
+.add_parameter("hdf5", ":py:class:`bob.io.base.HDF5File`", "An HDF5 file open for writing");
+static PyObject* PyBobLearnMiscJFAMachine_Save(PyBobLearnMiscJFAMachineObject* self,  PyObject* args, PyObject* kwargs) {
+
+  BOB_TRY
+  
+  // get list of arguments
+  char** kwlist = save.kwlist(0);  
+  PyBobIoHDF5FileObject* hdf5;
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&", kwlist, PyBobIoHDF5File_Converter, &hdf5)) return 0;
+
+  auto hdf5_ = make_safe(hdf5);
+  self->cxx->save(*hdf5->f);
+
+  BOB_CATCH_MEMBER("cannot save the data", 0)
+  Py_RETURN_NONE;
+}
+
+/*** load ***/
+static auto load = bob::extension::FunctionDoc(
+  "load",
+  "Load the configuration of the JFAMachine to a given HDF5 file"
+)
+.add_prototype("hdf5")
+.add_parameter("hdf5", ":py:class:`bob.io.base.HDF5File`", "An HDF5 file open for reading");
+static PyObject* PyBobLearnMiscJFAMachine_Load(PyBobLearnMiscJFAMachineObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+  
+  char** kwlist = load.kwlist(0);  
+  PyBobIoHDF5FileObject* hdf5;
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&", kwlist, PyBobIoHDF5File_Converter, &hdf5)) return 0;
+  
+  auto hdf5_ = make_safe(hdf5);  
+  self->cxx->load(*hdf5->f);
+
+  BOB_CATCH_MEMBER("cannot load the data", 0)
+  Py_RETURN_NONE;
+}
+
+
+/*** is_similar_to ***/
+static auto is_similar_to = bob::extension::FunctionDoc(
+  "is_similar_to",
+  
+  "Compares this JFAMachine with the ``other`` one to be approximately the same.",
+  "The optional values ``r_epsilon`` and ``a_epsilon`` refer to the "
+  "relative and absolute precision for the ``weights``, ``biases`` "
+  "and any other values internal to this machine."
+)
+.add_prototype("other, [r_epsilon], [a_epsilon]","output")
+.add_parameter("other", ":py:class:`bob.learn.misc.JFAMachine`", "A JFAMachine object to be compared.")
+.add_parameter("r_epsilon", "float", "Relative precision.")
+.add_parameter("a_epsilon", "float", "Absolute precision.")
+.add_return("output","bool","True if it is similar, otherwise false.");
+static PyObject* PyBobLearnMiscJFAMachine_IsSimilarTo(PyBobLearnMiscJFAMachineObject* self, PyObject* args, PyObject* kwds) {
+
+  /* Parses input arguments in a single shot */
+  char** kwlist = is_similar_to.kwlist(0);
+
+  //PyObject* other = 0;
+  PyBobLearnMiscJFAMachineObject* other = 0;
+  double r_epsilon = 1.e-5;
+  double a_epsilon = 1.e-8;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwds, "O!|dd", kwlist,
+        &PyBobLearnMiscJFAMachine_Type, &other,
+        &r_epsilon, &a_epsilon)){
+
+        is_similar_to.print_usage(); 
+        return 0;        
+  }
+
+  if (self->cxx->is_similar_to(*other->cxx, r_epsilon, a_epsilon))
+    Py_RETURN_TRUE;
+  else
+    Py_RETURN_FALSE;
+}
+
+
+/*** estimate_x ***/
+static auto estimate_x = bob::extension::FunctionDoc(
+  "estimate_x",
+  "Estimates the session offset x (LPT assumption) given GMM statistics.",
+  "Estimates x from the GMM statistics considering the LPT assumption, that is the latent session variable x is approximated using the UBM", 
+  true
+)
+.add_prototype("stats,input")
+.add_parameter("stats", ":py:class:`bob.learn.misc.GMMStats`", "Statistics of the GMM")
+.add_parameter("input", "array_like <float, 1D>", "Input vector");
+static PyObject* PyBobLearnMiscJFAMachine_estimateX(PyBobLearnMiscJFAMachineObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+
+  char** kwlist = estimate_x.kwlist(0);
+
+  PyBobLearnMiscGMMStatsObject* stats = 0;
+  PyBlitzArrayObject* input           = 0;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!O&", kwlist, &PyBobLearnMiscGMMStats_Type, &stats, 
+                                                                 &PyBlitzArray_Converter,&input))
+    Py_RETURN_NONE;
+
+  //protects acquired resources through this scope
+  auto input_ = make_safe(input);
+  self->cxx->estimateX(*stats->cxx, *PyBlitzArrayCxx_AsBlitz<double,1>(input));
+
+  BOB_CATCH_MEMBER("cannot estimate X", 0)
+  Py_RETURN_NONE;
+}
+
+
+/*** estimate_ux ***/
+static auto estimate_ux = bob::extension::FunctionDoc(
+  "estimate_ux",
+  "Estimates Ux (LPT assumption) given GMM statistics.",
+  "Estimates Ux from the GMM statistics considering the LPT assumption, that is the latent session variable x is approximated using the UBM.", 
+  true
+)
+.add_prototype("stats,input")
+.add_parameter("stats", ":py:class:`bob.learn.misc.GMMStats`", "Statistics of the GMM")
+.add_parameter("input", "array_like <float, 1D>", "Input vector");
+static PyObject* PyBobLearnMiscJFAMachine_estimateUx(PyBobLearnMiscJFAMachineObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+
+  char** kwlist = estimate_ux.kwlist(0);
+
+  PyBobLearnMiscGMMStatsObject* stats = 0;
+  PyBlitzArrayObject* input           = 0;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!O&", kwlist, &PyBobLearnMiscGMMStats_Type, &stats, 
+                                                                 &PyBlitzArray_Converter,&input))
+    Py_RETURN_NONE;
+
+  //protects acquired resources through this scope
+  auto input_ = make_safe(input);
+  self->cxx->estimateUx(*stats->cxx, *PyBlitzArrayCxx_AsBlitz<double,1>(input));
+
+  BOB_CATCH_MEMBER("cannot estimate Ux", 0)
+  Py_RETURN_NONE;
+}
+
+
+/*** forward_ux ***/
+static auto forward_ux = bob::extension::FunctionDoc(
+  "forward_ux",
+  "Computes a score for the given UBM statistics and given the Ux vector",
+  "", 
+  true
+)
+.add_prototype("stats,ux")
+.add_parameter("stats", ":py:class:`bob.learn.misc.GMMStats`", "Statistics as input")
+.add_parameter("ux", "array_like <float, 1D>", "Input vector");
+static PyObject* PyBobLearnMiscJFAMachine_ForwardUx(PyBobLearnMiscJFAMachineObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+
+  char** kwlist = forward_ux.kwlist(0);
+
+  PyBobLearnMiscGMMStatsObject* stats = 0;
+  PyBlitzArrayObject* ux_input        = 0;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!O&", kwlist, &PyBobLearnMiscGMMStats_Type, &stats, 
+                                                                 &PyBlitzArray_Converter,&ux_input))
+    Py_RETURN_NONE;
+
+  //protects acquired resources through this scope
+  auto ux_input_ = make_safe(ux_input);
+  double score = self->cxx->forward(*stats->cxx, *PyBlitzArrayCxx_AsBlitz<double,1>(ux_input));
+  
+  return Py_BuildValue("d", score);
+  BOB_CATCH_MEMBER("cannot forward_ux", 0)
+
+}
+
+
+/*** forward ***/
+static auto forward = bob::extension::FunctionDoc(
+  "forward",
+  "Execute the machine",
+  "", 
+  true
+)
+.add_prototype("stats")
+.add_parameter("stats", ":py:class:`bob.learn.misc.GMMStats`", "Statistics as input");
+static PyObject* PyBobLearnMiscJFAMachine_Forward(PyBobLearnMiscJFAMachineObject* self, PyObject* args, PyObject* kwargs) {
+  BOB_TRY
+
+  char** kwlist = forward.kwlist(0);
+
+  PyBobLearnMiscGMMStatsObject* stats = 0;
+  
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O!", kwlist, &PyBobLearnMiscGMMStats_Type, &stats))
+    Py_RETURN_NONE;
+
+  //protects acquired resources through this scope
+  double score = self->cxx->forward(*stats->cxx);
+
+  return Py_BuildValue("d", score);
+  BOB_CATCH_MEMBER("cannot forward", 0)
+
+}
+
+
+static PyMethodDef PyBobLearnMiscJFAMachine_methods[] = {
+  {
+    save.name(),
+    (PyCFunction)PyBobLearnMiscJFAMachine_Save,
+    METH_VARARGS|METH_KEYWORDS,
+    save.doc()
+  },
+  {
+    load.name(),
+    (PyCFunction)PyBobLearnMiscJFAMachine_Load,
+    METH_VARARGS|METH_KEYWORDS,
+    load.doc()
+  },
+  {
+    is_similar_to.name(),
+    (PyCFunction)PyBobLearnMiscJFAMachine_IsSimilarTo,
+    METH_VARARGS|METH_KEYWORDS,
+    is_similar_to.doc()
+  },
+  
+  {
+    estimate_x.name(),
+    (PyCFunction)PyBobLearnMiscJFAMachine_estimateX,
+    METH_VARARGS|METH_KEYWORDS,
+    estimate_x.doc()
+  },
+  
+  {
+    estimate_ux.name(),
+    (PyCFunction)PyBobLearnMiscJFAMachine_estimateUx,
+    METH_VARARGS|METH_KEYWORDS,
+    estimate_ux.doc()
+  },
+
+  {
+    forward_ux.name(),
+    (PyCFunction)PyBobLearnMiscJFAMachine_ForwardUx,
+    METH_VARARGS|METH_KEYWORDS,
+    forward_ux.doc()
+  },
+/*
+  {
+    forward.name(),
+    (PyCFunction)PyBobLearnMiscJFAMachine_Forward,
+    METH_VARARGS|METH_KEYWORDS,
+    forward.doc()
+  },*/
+
+
+  {0} /* Sentinel */
+};
+
+
+/******************************************************************/
+/************ Module Section **************************************/
+/******************************************************************/
+
+// Define the JFA type struct; will be initialized later
+PyTypeObject PyBobLearnMiscJFAMachine_Type = {
+  PyVarObject_HEAD_INIT(0,0)
+  0
+};
+
+bool init_BobLearnMiscJFAMachine(PyObject* module)
+{
+  // initialize the type struct
+  PyBobLearnMiscJFAMachine_Type.tp_name      = JFAMachine_doc.name();
+  PyBobLearnMiscJFAMachine_Type.tp_basicsize = sizeof(PyBobLearnMiscJFAMachineObject);
+  PyBobLearnMiscJFAMachine_Type.tp_flags     = Py_TPFLAGS_DEFAULT;
+  PyBobLearnMiscJFAMachine_Type.tp_doc       = JFAMachine_doc.doc();
+
+  // set the functions
+  PyBobLearnMiscJFAMachine_Type.tp_new         = PyType_GenericNew;
+  PyBobLearnMiscJFAMachine_Type.tp_init        = reinterpret_cast<initproc>(PyBobLearnMiscJFAMachine_init);
+  PyBobLearnMiscJFAMachine_Type.tp_dealloc     = reinterpret_cast<destructor>(PyBobLearnMiscJFAMachine_delete);
+  PyBobLearnMiscJFAMachine_Type.tp_richcompare = reinterpret_cast<richcmpfunc>(PyBobLearnMiscJFAMachine_RichCompare);
+  PyBobLearnMiscJFAMachine_Type.tp_methods     = PyBobLearnMiscJFAMachine_methods;
+  PyBobLearnMiscJFAMachine_Type.tp_getset      = PyBobLearnMiscJFAMachine_getseters;
+  PyBobLearnMiscJFAMachine_Type.tp_call        = reinterpret_cast<ternaryfunc>(PyBobLearnMiscJFAMachine_Forward);
+
+
+  // check that everything is fine
+  if (PyType_Ready(&PyBobLearnMiscJFAMachine_Type) < 0) return false;
+
+  // add the type to the module
+  Py_INCREF(&PyBobLearnMiscJFAMachine_Type);
+  return PyModule_AddObject(module, "JFAMachine", (PyObject*)&PyBobLearnMiscJFAMachine_Type) >= 0;
+}
+
diff --git a/bob/learn/misc/main.cpp b/bob/learn/misc/main.cpp
index 6456cdc96b1611441399e7399b9919d31dc7f2b9..f57fdfe3f96beb72d0f3eaeb9870793f3c5ac4f2 100644
--- a/bob/learn/misc/main.cpp
+++ b/bob/learn/misc/main.cpp
@@ -49,6 +49,10 @@ static PyObject* create_module (void) {
   if (!init_BobLearnMiscMLGMMTrainer(module)) return 0;  
   if (!init_BobLearnMiscMAPGMMTrainer(module)) return 0;
   if (!init_BobLearnMiscJFABase(module)) return 0;
+  if (!init_BobLearnMiscISVBase(module)) return 0;
+
+  if (!init_BobLearnMiscJFAMachine(module)) return 0;
+  if (!init_BobLearnMiscISVMachine(module)) return 0;  
 
 
   static void* PyBobLearnMisc_API[PyBobLearnMisc_API_pointers];
diff --git a/bob/learn/misc/main.h b/bob/learn/misc/main.h
index 64b501e7a8278b0af2661e8e201139c43232d726..aa36669e70d44455100f17211043fa698fad586a 100644
--- a/bob/learn/misc/main.h
+++ b/bob/learn/misc/main.h
@@ -28,6 +28,10 @@
 #include <bob.learn.misc/MAP_GMMTrainer.h>
 
 #include <bob.learn.misc/JFABase.h>
+#include <bob.learn.misc/ISVBase.h>
+
+#include <bob.learn.misc/JFAMachine.h>
+#include <bob.learn.misc/ISVMachine.h>
 
 
 #if PY_VERSION_HEX >= 0x03000000
@@ -167,4 +171,38 @@ bool init_BobLearnMiscJFABase(PyObject* module);
 int PyBobLearnMiscJFABase_Check(PyObject* o);
 
 
+// ISVBase
+typedef struct {
+  PyObject_HEAD
+  boost::shared_ptr<bob::learn::misc::ISVBase> cxx;
+} PyBobLearnMiscISVBaseObject;
+
+extern PyTypeObject PyBobLearnMiscISVBase_Type;
+bool init_BobLearnMiscISVBase(PyObject* module);
+int PyBobLearnMiscISVBase_Check(PyObject* o);
+
+
+// JFAMachine
+typedef struct {
+  PyObject_HEAD
+  boost::shared_ptr<bob::learn::misc::JFAMachine> cxx;
+} PyBobLearnMiscJFAMachineObject;
+
+extern PyTypeObject PyBobLearnMiscJFAMachine_Type;
+bool init_BobLearnMiscJFAMachine(PyObject* module);
+int PyBobLearnMiscJFAMachine_Check(PyObject* o);
+
+
+// ISVMachine
+typedef struct {
+  PyObject_HEAD
+  boost::shared_ptr<bob::learn::misc::ISVMachine> cxx;
+} PyBobLearnMiscISVMachineObject;
+
+extern PyTypeObject PyBobLearnMiscISVMachine_Type;
+bool init_BobLearnMiscISVMachine(PyObject* module);
+int PyBobLearnMiscISVMachine_Check(PyObject* o);
+
+
+
 #endif // BOB_LEARN_EM_MAIN_H
diff --git a/bob/learn/misc/test_jfa.py b/bob/learn/misc/test_jfa.py
index 1434a3e443a33810cae6062261a1b8a6f59eff7f..594ada436ef1d66614302066b4e20d275b9b2a31 100644
--- a/bob/learn/misc/test_jfa.py
+++ b/bob/learn/misc/test_jfa.py
@@ -63,23 +63,28 @@ def test_JFABase():
   U = numpy.array([[1, 2], [3, 4], [5, 6], [7, 8], [9, 10], [11, 12]], 'float64')
   V = numpy.array([[6, 5], [4, 3], [2, 1], [1, 2], [3, 4], [5, 6]], 'float64')
   d = numpy.array([0, 1, 0, 1, 0, 1], 'float64')
-  m = JFABase(ubm)
-  assert m.dim_ru == 1
-  assert m.dim_rv == 1
+  m = JFABase(ubm, ru=1, rv=1)
+  
+  _,_,ru,rv = m.shape 
+  assert ru == 1
+  assert rv == 1
 
   # Checks for correctness
   m.resize(2,2)
   m.u = U
   m.v = V
-  m.d = d
+  m.d = d  
+  n_gaussians,dim,ru,rv = m.shape
+  supervector_length    = m.supervector_length
+  
   assert (m.u == U).all()
   assert (m.v == V).all()
-  assert (m.d == d).all()
-  assert m.dim_c == 2
-  assert m.dim_d == 3
-  assert m.dim_cd == 6
-  assert m.dim_ru == 2
-  assert m.dim_rv == 2
+  assert (m.d == d).all()  
+  assert n_gaussians        == 2
+  assert dim                == 3
+  assert supervector_length == 6
+  assert ru                 == 2
+  assert rv                 == 2
 
   # Saves and loads
   filename = str(tempfile.mkstemp(".hdf5")[1])
@@ -95,21 +100,21 @@ def test_JFABase():
   assert m == mc
 
   # Variant
-  mv = JFABase()
+  #mv = JFABase()
   # Checks for correctness
-  mv.ubm = ubm
-  mv.resize(2,2)
-  mv.u = U
-  mv.v = V
-  mv.d = d
-  assert (m.u == U).all()
-  assert (m.v == V).all()
-  assert (m.d == d).all()
-  assert m.dim_c == 2
-  assert m.dim_d == 3
-  assert m.dim_cd == 6
-  assert m.dim_ru == 2
-  assert m.dim_rv == 2
+  #mv.ubm = ubm
+  #mv.resize(2,2)
+  #mv.u = U
+  #mv.v = V
+  #mv.d = d
+  #assert (m.u == U).all()
+  #assert (m.v == V).all()
+  #assert (m.d == d).all()
+  #assert m.dim_c == 2
+  #assert m.dim_d == 3
+  #assert m.dim_cd == 6
+  #assert m.dim_ru == 2
+  #assert m.dim_rv == 2
 
   # Clean-up
   os.unlink(filename)
@@ -120,27 +125,30 @@ def test_ISVBase():
   weights = numpy.array([0.4, 0.6], 'float64')
   means = numpy.array([[1, 6, 2], [4, 3, 2]], 'float64')
   variances = numpy.array([[1, 2, 1], [2, 1, 2]], 'float64')
-  ubm = GMMMachine(2,3)
-  ubm.weights = weights
-  ubm.means = means
+  ubm           = GMMMachine(2,3)
+  ubm.weights   = weights
+  ubm.means     = means
   ubm.variances = variances
 
   # Creates a ISVBase
   U = numpy.array([[1, 2], [3, 4], [5, 6], [7, 8], [9, 10], [11, 12]], 'float64')
   d = numpy.array([0, 1, 0, 1, 0, 1], 'float64')
-  m = ISVBase(ubm)
-  assert m.dim_ru == 1
+  m = ISVBase(ubm, ru=1)
+  _,_,ru = m.shape
+  assert ru == 1
 
   # Checks for correctness
   m.resize(2)
   m.u = U
   m.d = d
+  n_gaussians,dim,ru = m.shape
+  supervector_length = m.supervector_length
   assert (m.u == U).all()
-  assert (m.d == d).all()
-  assert m.dim_c == 2
-  assert m.dim_d == 3
-  assert m.dim_cd == 6
-  assert m.dim_ru == 2
+  assert (m.d == d).all()  
+  assert n_gaussians        == 2
+  assert dim                == 3
+  assert supervector_length == 6
+  assert ru                 == 2
 
   # Saves and loads
   filename = str(tempfile.mkstemp(".hdf5")[1])
@@ -156,18 +164,18 @@ def test_ISVBase():
   assert m == mc
 
   # Variant
-  mv = ISVBase()
+  #mv = ISVBase()
   # Checks for correctness
-  mv.ubm = ubm
-  mv.resize(2)
-  mv.u = U
-  mv.d = d
-  assert (m.u == U).all()
-  assert (m.d == d).all()
-  assert m.dim_c == 2
-  assert m.dim_d == 3
-  assert m.dim_cd == 6
-  assert m.dim_ru == 2
+  #mv.ubm = ubm
+  #mv.resize(2)
+  #mv.u = U
+  #mv.d = d
+  #assert (m.u == U).all()
+  #assert (m.d == d).all()
+  #ssert m.dim_c == 2
+  #assert m.dim_d == 3
+  #assert m.dim_cd == 6
+  #assert m.dim_ru == 2
 
   # Clean-up
   os.unlink(filename)
@@ -175,12 +183,12 @@ def test_ISVBase():
 def test_JFAMachine():
 
   # Creates a UBM
-  weights = numpy.array([0.4, 0.6], 'float64')
-  means = numpy.array([[1, 6, 2], [4, 3, 2]], 'float64')
+  weights   = numpy.array([0.4, 0.6], 'float64')
+  means     = numpy.array([[1, 6, 2], [4, 3, 2]], 'float64')
   variances = numpy.array([[1, 2, 1], [2, 1, 2]], 'float64')
-  ubm = GMMMachine(2,3)
-  ubm.weights = weights
-  ubm.means = means
+  ubm           = GMMMachine(2,3)
+  ubm.weights   = weights
+  ubm.means     = means
   ubm.variances = variances
 
   # Creates a JFABase
@@ -198,11 +206,14 @@ def test_JFAMachine():
   m = JFAMachine(base)
   m.y = y
   m.z = z
-  assert m.dim_c == 2
-  assert m.dim_d == 3
-  assert m.dim_cd == 6
-  assert m.dim_ru == 2
-  assert m.dim_rv == 2
+  n_gaussians,dim,ru,rv = m.shape
+  supervector_length    = m.supervector_length  
+  
+  assert n_gaussians        == 2
+  assert dim                == 3
+  assert supervector_length == 6
+  assert ru                 == 2
+  assert rv                 == 2
   assert (m.y == y).all()
   assert (m.z == z).all()
 
@@ -220,18 +231,18 @@ def test_JFAMachine():
   assert m == mc
 
   # Variant
-  mv = JFAMachine()
+  #mv = JFAMachine()
   # Checks for correctness
-  mv.jfa_base = base
-  m.y = y
-  m.z = z
-  assert m.dim_c == 2
-  assert m.dim_d == 3
-  assert m.dim_cd == 6
-  assert m.dim_ru == 2
-  assert m.dim_rv == 2
-  assert (m.y == y).all()
-  assert (m.z == z).all()
+  #mv.jfa_base = base
+  #m.y = y
+  #m.z = z
+  #assert m.dim_c == 2
+  #assert m.dim_d == 3
+  #assert m.dim_cd == 6
+  #assert m.dim_ru == 2
+  #assert m.dim_rv == 2
+  #assert (m.y == y).all()
+  #assert (m.z == z).all()
 
   # Defines GMMStats
   gs = GMMStats(2,3)
@@ -250,23 +261,26 @@ def test_JFAMachine():
   eps = 1e-10
   x_ref = numpy.array([0.291042849767692, 0.310273618998444], 'float64')
   score_ref = -2.111577181208289
-  score = m.forward(gs)
-  assert numpy.allclose(m.__x__, x_ref, eps)
+  score = m(gs)
+  assert numpy.allclose(m.x, x_ref, eps)
   assert abs(score_ref-score) < eps
 
   # x and Ux
   x = numpy.ndarray((2,), numpy.float64)
   m.estimate_x(gs, x)
-  x_py = estimate_x(m.dim_c, m.dim_d, ubm.mean_supervector, ubm.variance_supervector, U, n, sumpx)
+  n_gaussians, dim,_,_ = m.shape
+  x_py = estimate_x(n_gaussians, dim, ubm.mean_supervector, ubm.variance_supervector, U, n, sumpx)
   assert numpy.allclose(x, x_py, eps)
 
   ux = numpy.ndarray((6,), numpy.float64)
   m.estimate_ux(gs, ux)
-  ux_py = estimate_ux(m.dim_c, m.dim_d, ubm.mean_supervector, ubm.variance_supervector, U, n, sumpx)
+  n_gaussians, dim,_,_ = m.shape  
+  ux_py = estimate_ux(n_gaussians, dim, ubm.mean_supervector, ubm.variance_supervector, U, n, sumpx)
   assert numpy.allclose(ux, ux_py, eps)
-  assert numpy.allclose(m.__x__, x, eps)
+  assert numpy.allclose(m.x, x, eps)
 
   score = m.forward_ux(gs, ux)
+
   assert abs(score_ref-score) < eps
 
   # Clean-up
@@ -283,23 +297,26 @@ def test_ISVMachine():
   ubm.means = means
   ubm.variances = variances
 
-  # Creates a JFABaseMachine
+  # Creates a ISVBaseMachine
   U = numpy.array([[1, 2], [3, 4], [5, 6], [7, 8], [9, 10], [11, 12]], 'float64')
-  V = numpy.array([[0], [0], [0], [0], [0], [0]], 'float64')
+  #V = numpy.array([[0], [0], [0], [0], [0], [0]], 'float64')
   d = numpy.array([0, 1, 0, 1, 0, 1], 'float64')
   base = ISVBase(ubm,2)
   base.u = U
-  base.v = V
+  #base.v = V
   base.d = d
 
   # Creates a JFAMachine
   z = numpy.array([3,4,1,2,0,1], 'float64')
   m = ISVMachine(base)
   m.z = z
-  assert m.dim_c == 2
-  assert m.dim_d == 3
-  assert m.dim_cd == 6
-  assert m.dim_ru == 2
+  
+  n_gaussians,dim,ru    = m.shape
+  supervector_length    = m.supervector_length  
+  assert n_gaussians          == 2
+  assert dim                  == 3
+  assert supervector_length   == 6
+  assert ru                   == 2
   assert (m.z == z).all()
 
   # Saves and loads
@@ -316,14 +333,17 @@ def test_ISVMachine():
   assert m == mc
 
   # Variant
-  mv = ISVMachine()
+  mv = ISVMachine(base)
   # Checks for correctness
-  mv.isv_base = base
+  #mv.isv_base = base
   m.z = z
-  assert m.dim_c == 2
-  assert m.dim_d == 3
-  assert m.dim_cd == 6
-  assert m.dim_ru == 2
+
+  n_gaussians,dim,ru    = m.shape
+  supervector_length    = m.supervector_length  
+  assert n_gaussians        == 2
+  assert dim                == 3
+  assert supervector_length == 6
+  assert ru                 == 2
   assert (m.z == z).all()
 
   # Defines GMMStats
@@ -344,12 +364,13 @@ def test_ISVMachine():
   x_ref = numpy.array([0.291042849767692, 0.310273618998444], 'float64')
   score_ref = -3.280498193082100
 
-  score = m.forward(gs)
-  assert numpy.allclose(m.__x__, x_ref, eps)
+  score = m(gs)
+  assert numpy.allclose(m.x, x_ref, eps)  
   assert abs(score_ref-score) < eps
 
   # Check using alternate forward() method
-  Ux = numpy.ndarray(shape=(m.dim_cd,), dtype=numpy.float64)
+  supervector_length = m.supervector_length
+  Ux = numpy.ndarray(shape=(supervector_length,), dtype=numpy.float64)
   m.estimate_ux(gs, Ux)
   score = m.forward_ux(gs, Ux)
   assert abs(score_ref-score) < eps
@@ -357,14 +378,16 @@ def test_ISVMachine():
   # x and Ux
   x = numpy.ndarray((2,), numpy.float64)
   m.estimate_x(gs, x)
-  x_py = estimate_x(m.dim_c, m.dim_d, ubm.mean_supervector, ubm.variance_supervector, U, n, sumpx)
+  n_gaussians,dim,_    = m.shape  
+  x_py = estimate_x(n_gaussians, dim, ubm.mean_supervector, ubm.variance_supervector, U, n, sumpx)
   assert numpy.allclose(x, x_py, eps)
 
   ux = numpy.ndarray((6,), numpy.float64)
   m.estimate_ux(gs, ux)
-  ux_py = estimate_ux(m.dim_c, m.dim_d, ubm.mean_supervector, ubm.variance_supervector, U, n, sumpx)
+  n_gaussians,dim,_    = m.shape  
+  ux_py = estimate_ux(n_gaussians, dim, ubm.mean_supervector, ubm.variance_supervector, U, n, sumpx)
   assert numpy.allclose(ux, ux_py, eps)
-  assert numpy.allclose(m.__x__, x, eps)
+  assert numpy.allclose(m.x, x, eps)
 
   score = m.forward_ux(gs, ux)
   assert abs(score_ref-score) < eps
diff --git a/doc/index.rst b/doc/index.rst
index 019f09966e76afdd71ac1a2215ba1f80101e060a..04bf7cd0bdb654f06e62096b97f38f52a99c6cde 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -28,6 +28,8 @@ References
 -----------
 
 .. [Reynolds2000] *Reynolds, Douglas A., Thomas F. Quatieri, and Robert B. Dunn*. **Speaker Verification Using Adapted Gaussian Mixture Models**, Digital signal processing 10.1 (2000): 19-41.
+..   [Vogt2008]   *R. Vogt, S. Sridharan*. **'Explicit Modelling of Session Variability for Speaker Verification'**, Computer Speech & Language, 2008, vol. 22, no. 1, pp. 17-38
+..   [McCool2013] *C. McCool, R. Wallace, M. McLaren, L. El Shafey, S. Marcel*. **'Session Variability Modelling for Face Authentication'**, IET Biometrics, 2013
 
 
 Indices and tables
diff --git a/setup.py b/setup.py
index 7939d6a45931ced56b616eedc1d88d54a7ca26e1..187d11ccf38b0d4026ab960b0075de5d0d3f0664 100644
--- a/setup.py
+++ b/setup.py
@@ -58,12 +58,17 @@ setup(
           "bob/learn/misc/cpp/GMMMachine.cpp",
           "bob/learn/misc/cpp/GMMStats.cpp",
           #"bob/learn/misc/cpp/IVectorMachine.cpp",
-          #"bob/learn/misc/cpp/JFAMachine.cpp",
           "bob/learn/misc/cpp/KMeansMachine.cpp",
-          #"bob/learn/misc/cpp/LinearScoring.cpp",
+          "bob/learn/misc/cpp/LinearScoring.cpp",
           #"bob/learn/misc/cpp/PLDAMachine.cpp",
           #"bob/learn/misc/cpp/ZTNorm.cpp",
 
+          "bob/learn/misc/cpp/FABase.cpp",
+          "bob/learn/misc/cpp/JFABase.cpp",
+          "bob/learn/misc/cpp/ISVBase.cpp",
+          "bob/learn/misc/cpp/JFAMachine.cpp",
+          "bob/learn/misc/cpp/ISVMachine.cpp",
+
           #"bob/learn/misc/cpp/EMPCATrainer.cpp",
           "bob/learn/misc/cpp/GMMBaseTrainer.cpp",
           #"bob/learn/misc/cpp/IVectorTrainer.cpp",
@@ -108,6 +113,10 @@ setup(
           "bob/learn/misc/gmm_base_trainer.cpp",
           "bob/learn/misc/ML_gmm_trainer.cpp",
           "bob/learn/misc/MAP_gmm_trainer.cpp",
+          "bob/learn/misc/jfa_base.cpp",
+          "bob/learn/misc/isv_base.cpp",
+          "bob/learn/misc/jfa_machine.cpp",
+          "bob/learn/misc/isv_machine.cpp",
 
           "bob/learn/misc/main.cpp",
         ],