diff --git a/doc/conf.py b/doc/conf.py
index 1460c68d5d8ae132a5fa0f947ecb6088723bc7e6..befc2e60c007b7fca49fac952f3ae3ae93903309 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -201,7 +201,7 @@ latex_font_size = '10pt'
 # Grouping the document tree into LaTeX files. List of tuples
 # (source start file, target name, title, author, documentclass [howto/manual]).
 latex_documents = [
-  ('index', 'xbob_learn_linear.tex', u'Bob Machines',
+  ('index', 'xbob_learn_linear.tex', u'Bob Linear Machines and Trainers',
    u'Biometrics Group, Idiap Research Institute', 'manual'),
 ]
 
@@ -241,7 +241,7 @@ rst_epilog = """
 # One entry per manual page. List of tuples
 # (source start file, name, description, authors, manual section).
 man_pages = [
-    ('index', 'xbob_learn_linear', u'Bob Machines Documentation', [u'Idiap Research Institute'], 1)
+    ('index', 'xbob_learn_linear', u'Bob Linear Machines and Trainers', [u'Idiap Research Institute'], 1)
 ]
 
 # Default processing flags for sphinx
diff --git a/doc/index.rst b/doc/index.rst
index 10bf9ce98812feaa19f1831ec71f64b39b65239f..8402a551eeaf1b74e2bfddead63d23576d564c93 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -2,17 +2,25 @@
 .. Andre Anjos <andre.anjos@idiap.ch>
 .. Fri 13 Dec 2013 12:50:06 CET
 ..
-.. Copyright (C) 2011-2013 Idiap Research Institute, Martigny, Switzerland
+.. Copyright (C) 2011-2014 Idiap Research Institute, Martigny, Switzerland
 
-==============
- Bob Machines
-==============
+==================================
+ Bob Linear Machines and Trainers
+==================================
 
 .. todolist::
 
-This module contains base functionality from Bob bound to Python, available in
-the C++ counter-part ``bob::machine``. It includes machines from our Machine
-Learning core.
+This package includes the definition of a linear machine, which is capable of
+either projecting the input data into maximally spread representations,
+linearly, or providing linear separation planes for multi-class data samples.
+The package includes the machine definition *per se* and a selection of
+different trainers for specialized purposes:
+
+ * Principal Component Analysis
+ * Fisher's Linear Discriminant Analysis
+ * (Conjugate Gradient) Logistic Regression
+ * Whitening filter
+ * Within-class covariance normalization (WCCN)
 
 Documentation
 -------------
diff --git a/setup.py b/setup.py
index 24249a19ab762a3191f8d133e0862d675951314a..c0f5a21f90e9cec7005d2323e954cd78a615cbc8 100644
--- a/setup.py
+++ b/setup.py
@@ -19,7 +19,7 @@ include_dirs = [
     xbob.learn.activation.get_include()
     ]
 
-packages = ['bob-machine >= 1.2.2', 'bob-trainer >= 1.2.2']
+packages = ['bob-machine >= 2.0.0a2', 'bob-io >= 2.0.0a2']
 version = '2.0.0a0'
 
 setup(
@@ -64,6 +64,12 @@ setup(
           "xbob/learn/linear/pca.cpp",
           "xbob/learn/linear/lda.cpp",
           "xbob/learn/linear/main.cpp",
+          "xbob/learn/linear/cpp/machine.cpp",
+          "xbob/learn/linear/cpp/pca.cpp",
+          "xbob/learn/linear/cpp/lda.cpp",
+          "xbob/learn/linear/cpp/logreg.cpp",
+          "xbob/learn/linear/cpp/wccn.cpp",
+          "xbob/learn/linear/cpp/whitening.cpp",
           ],
         packages = packages,
         include_dirs = include_dirs,
diff --git a/xbob/learn/linear/cpp/lda.cpp b/xbob/learn/linear/cpp/lda.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e9614ebef56af55ffb6eb0a582b0a68f9f0f276a
--- /dev/null
+++ b/xbob/learn/linear/cpp/lda.cpp
@@ -0,0 +1,199 @@
+/**
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @date Sat Jun 4 21:38:59 2011 +0200
+ *
+ * @brief Implements a multi-class Fisher/LDA linear machine Training using
+ * Singular Value Decomposition (SVD). For more information on Linear Machines
+ * and associated methods, please consult Bishop, Machine Learning and Pattern
+ * Recognition chapter 4.
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#include <boost/format.hpp>
+#include <bob/math/pinv.h>
+#include <bob/math/eig.h>
+#include <bob/math/linear.h>
+#include <bob/math/stats.h>
+
+#include <xbob.learn.linear/lda.h>
+
+namespace bob { namespace learn { namespace linear {
+
+  FisherLDATrainer::FisherLDATrainer
+    (bool use_pinv, bool strip_to_rank)
+    : m_use_pinv(use_pinv),
+      m_strip_to_rank(strip_to_rank)
+  {
+  }
+
+  FisherLDATrainer::FisherLDATrainer
+    (const FisherLDATrainer& other)
+    : m_use_pinv(other.m_use_pinv),
+      m_strip_to_rank(other.m_strip_to_rank)
+  {
+  }
+
+  FisherLDATrainer::~FisherLDATrainer()
+  {
+  }
+
+  FisherLDATrainer& FisherLDATrainer::operator=
+    (const FisherLDATrainer& other)
+    {
+      if (this != &other) // avoid auto assignment
+      {
+        m_use_pinv = other.m_use_pinv;
+        m_strip_to_rank = other.m_strip_to_rank;
+      }
+      return *this;
+    }
+
+  bool FisherLDATrainer::operator==
+    (const FisherLDATrainer& other) const
+    {
+      return m_use_pinv == other.m_use_pinv && \
+                         m_strip_to_rank == other.m_strip_to_rank;
+    }
+
+  bool FisherLDATrainer::operator!=
+    (const FisherLDATrainer& other) const
+    {
+      return !(this->operator==(other));
+    }
+
+  bool FisherLDATrainer::is_similar_to
+    (const FisherLDATrainer& other, const double r_epsilon,
+     const double a_epsilon) const
+    {
+      return this->operator==(other);
+    }
+
+  /**
+   * Returns the indexes for sorting a given blitz::Array<double,1>
+   */
+  struct compare_1d_blitz {
+    const blitz::Array<double,1>& v_;
+    compare_1d_blitz(const blitz::Array<double,1>& v): v_(v) { }
+    bool operator() (size_t i, size_t j) { return v_(i) < v_(j); }
+  };
+
+  static std::vector<size_t> sort_indexes(const blitz::Array<double,1>& v) {
+
+    // initialize original index locations
+    std::vector<size_t> idx(v.size());
+    for (size_t i = 0; i != idx.size(); ++i) idx[i] = i;
+
+    // sort indexes based on comparing values in v
+    std::sort(idx.begin(), idx.end(), compare_1d_blitz(v));
+
+    return idx;
+  }
+
+  void FisherLDATrainer::train
+    (Machine& machine, blitz::Array<double,1>& eigen_values,
+     const std::vector<blitz::Array<double, 2> >& data) const
+    {
+      // if #classes < 2, then throw
+      if (data.size() < 2) {
+        boost::format m("The number of arrays in the input data == %d whereas for LDA you should provide at least 2");
+        m % data.size();
+        throw std::runtime_error(m.str());
+      }
+
+      // checks for arrayset data type and shape once
+      int n_features = data[0].extent(1);
+
+      for (size_t cl=0; cl<data.size(); ++cl) {
+        if (data[cl].extent(1) != n_features) {
+          boost::format m("The number of features/columns (%d) in array at position %d of your input differs from that of array at position 0 (%d)");
+          m % data[cl].extent(1) % n_features;
+          throw std::runtime_error(m.str());
+        }
+      }
+
+      int osize = output_size(data);
+
+      // Checks that the dimensions are matching
+      if (machine.inputSize() != (size_t)data[0].extent(1)) {
+        boost::format m("Number of features at input data set (%d columns) does not match machine input size (%d)");
+        m % data[0].extent(1) % machine.inputSize();
+        throw std::runtime_error(m.str());
+      }
+      if (machine.outputSize() != (size_t)osize) {
+        boost::format m("Number of outputs of the given machine (%d) does not match the expected number of outputs calculated by this trainer = %d");
+        m % machine.outputSize() % osize;
+        throw std::runtime_error(m.str());
+      }
+      if (eigen_values.extent(0) != osize) {
+        boost::format m("Number of eigenvalues on the given 1D array (%d) does not match the expected number of outputs calculated by this trainer = %d");
+        m % eigen_values.extent(0) % osize;
+        throw std::runtime_error(m.str());
+      }
+
+      blitz::Array<double,1> preMean(n_features);
+      blitz::Array<double,2> Sw(n_features, n_features);
+      blitz::Array<double,2> Sb(n_features, n_features);
+      bob::math::scatters_(data, Sw, Sb, preMean);
+
+      // computes the generalized eigenvalue decomposition
+      // so to find the eigen vectors/values of Sw^(-1) * Sb
+      blitz::Array<double,2> V(Sw.shape());
+      blitz::Array<double,1> eigen_values_(n_features);
+
+      if (m_use_pinv) {
+
+        //note: misuse V and Sw as temporary place holders for data
+        bob::math::pinv_(Sw, V); //V now contains Sw^-1
+        bob::math::prod_(V, Sb, Sw); //Sw now contains Sw^-1*Sb
+        blitz::Array<std::complex<double>,1> Dtemp(eigen_values_.shape());
+        blitz::Array<std::complex<double>,2> Vtemp(V.shape());
+        bob::math::eig_(Sw, Vtemp, Dtemp); //V now contains eigen-vectors
+
+        //sorting: we know this problem on has real eigen-values
+        blitz::Range a = blitz::Range::all();
+        blitz::Array<double,1> Dunordered(blitz::real(Dtemp));
+        std::vector<size_t> order = sort_indexes(Dunordered);
+        for (int i=0; i<n_features; ++i) {
+          eigen_values_(i) = Dunordered(order[i]);
+          V(a,i) = blitz::real(Vtemp(a,order[i]));
+        }
+      }
+      else {
+        bob::math::eigSym_(Sb, Sw, V, eigen_values_);
+      }
+
+      // Convert ascending order to descending order
+      eigen_values_.reverseSelf(0);
+      V.reverseSelf(1);
+
+      // limit the dimensions of the resulting projection matrix and eigen values
+      eigen_values = eigen_values_(blitz::Range(0,osize-1));
+      V.resizeAndPreserve(V.extent(0), osize);
+
+      // normalizes the eigen vectors so they have unit length
+      blitz::Range a = blitz::Range::all();
+      for (int column=0; column<V.extent(1); ++column) {
+        math::normalizeSelf(V(a,column));
+      }
+
+      // updates the machine
+      machine.setWeights(V);
+      machine.setInputSubtraction(preMean);
+
+      // also set input_div and biases to neutral values...
+      machine.setInputDivision(1.0);
+      machine.setBiases(0.0);
+    }
+
+  void FisherLDATrainer::train(Machine& machine,
+      const std::vector<blitz::Array<double,2> >& data) const {
+    blitz::Array<double,1> throw_away(output_size(data));
+    train(machine, throw_away, data);
+  }
+
+  size_t FisherLDATrainer::output_size(const std::vector<blitz::Array<double,2> >& data) const {
+    return m_strip_to_rank ? std::min(data.size()-1, (size_t)data[0].extent(1)) : data[0].extent(1);
+  }
+
+}}}
diff --git a/xbob/learn/linear/cpp/logreg.cpp b/xbob/learn/linear/cpp/logreg.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..230ac487a4e48a422b09a51df1bff71d08944319
--- /dev/null
+++ b/xbob/learn/linear/cpp/logreg.cpp
@@ -0,0 +1,227 @@
+/**
+ * @date Sat Sep 1 19:26:00 2012 +0100
+ * @author Laurent El Shafey <laurent.el-shafey@idiap.ch>
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#include <bob/core/logging.h>
+#include <bob/math/linear.h>
+#include <limits>
+
+#include <xbob.learn.linear/logreg.h>
+
+namespace bob { namespace learn { namespace linear {
+
+  CGLogRegTrainer::CGLogRegTrainer(const double prior,
+      const double convergence_threshold, const size_t max_iterations,
+      const double lambda, const bool mean_std_norm):
+    m_prior(prior),
+    m_convergence_threshold(convergence_threshold),
+    m_max_iterations(max_iterations),
+    m_lambda(lambda),
+    m_mean_std_norm(mean_std_norm)
+  {
+    if(prior<=0. || prior>=1.)
+    {
+      boost::format m("Prior (%f) not in the range ]0,1[.");
+      m % prior;
+      throw std::runtime_error(m.str());
+    }
+  }
+
+  CGLogRegTrainer::CGLogRegTrainer(const CGLogRegTrainer& other):
+    m_prior(other.m_prior),
+    m_convergence_threshold(other.m_convergence_threshold),
+    m_max_iterations(other.m_max_iterations),
+    m_lambda(other.m_lambda),
+    m_mean_std_norm(other.m_mean_std_norm)
+  {
+  }
+
+  CGLogRegTrainer::~CGLogRegTrainer() {}
+
+  CGLogRegTrainer& CGLogRegTrainer::operator=
+    (const CGLogRegTrainer& other)
+    {
+      if(this != &other)
+      {
+        m_prior = other.m_prior;
+        m_convergence_threshold = other.m_convergence_threshold;
+        m_max_iterations = other.m_max_iterations;
+        m_lambda = other.m_lambda;
+        m_mean_std_norm = other.m_mean_std_norm;
+      }
+      return *this;
+    }
+
+  bool CGLogRegTrainer::operator==(const CGLogRegTrainer& b) const {
+    return (this->m_prior == b.m_prior &&
+        this->m_convergence_threshold == b.m_convergence_threshold &&
+        this->m_max_iterations == b.m_max_iterations &&
+        this->m_lambda == b.m_lambda &&
+        this->m_mean_std_norm == b.m_mean_std_norm);
+  }
+
+  bool CGLogRegTrainer::operator!=(const CGLogRegTrainer& b) const {
+    return !(this->operator==(b));
+  }
+
+  void CGLogRegTrainer::train(Machine& machine, const blitz::Array<double,2>& negatives, const blitz::Array<double,2>& positives) const {
+
+    // Checks for arraysets data type and shape once
+    bob::core::array::assertSameDimensionLength(negatives.extent(1), positives.extent(1));
+
+    // Data is checked now and conforms, just proceed w/o any further checks.
+    size_t n_samples1 = positives.extent(0);
+    size_t n_samples2 = negatives.extent(0);
+    size_t n_samples = n_samples1 + n_samples2;
+    size_t n_features = positives.extent(1);
+
+    // Defines useful ranges
+    blitz::Range rall = blitz::Range::all();
+    blitz::Range rd = blitz::Range(0,n_features-1);
+    blitz::Range r1 = blitz::Range(0,n_samples1-1);
+    blitz::Range r2 = blitz::Range(n_samples1,n_samples-1);
+
+    // indices for blitz iterations
+    blitz::firstIndex i;
+    blitz::secondIndex j;
+
+    blitz::Array<double,1> mean(n_features);
+    blitz::Array<double,1> std_dev(n_features);
+    // mean and variance of the training data
+    if (m_mean_std_norm){
+      // collect all samples in one matrix
+      blitz::Array<double,2> all_data(n_samples,n_features);
+      for(size_t i=0; i<n_samples1; ++i)
+        all_data(i,rall) = positives(i,rall);
+      for(size_t i=0; i<n_samples2; ++i)
+        all_data(i+n_samples1, rall) = negatives(i,rall);
+
+      // compute mean and std-dev from samples
+      mean = blitz::mean(all_data(j,i), j);
+      blitz::Array<double,2> squared(blitz::pow2(all_data));
+      std_dev = blitz::sqrt((blitz::sum(squared(j,i), j) - n_samples*blitz::pow2(mean)) / n_samples);
+    } else {
+      mean = 0.;
+      std_dev = 1.;
+    }
+
+    // Creates a large blitz::Array containing the samples
+    // x = |positives - negatives|, of size (n_features+1,n_samples1+n_samples2)
+    //     |1.  -1. |
+    blitz::Array<double,2> x(n_features+1, n_samples);
+    x(n_features,r1) = 1.;
+    x(n_features,r2) = -1.;
+    for(size_t i=0; i<n_samples1; ++i)
+      x(rd,i) = (positives(i,rall) - mean(rall)) / std_dev(rall);
+    for(size_t i=0; i<n_samples2; ++i)
+      x(rd,i+n_samples1) = -(negatives(i,rall) - mean(rall)) / std_dev(rall);
+
+    // Ratio between the two classes and weights vector
+    double prop = (double)n_samples1 / (double)n_samples;
+    blitz::Array<double,1> weights(n_samples);
+    weights(r1) = m_prior / prop;
+    weights(r2) = (1.-m_prior) / (1.-prop);
+
+    // Initializes offset vector
+    blitz::Array<double,1> offset(n_samples);
+    const double logit = log(m_prior/(1.-m_prior));
+    offset(r1) = logit;
+    offset(r2) = -logit;
+
+    // Initializes gradient and w vectors
+    blitz::Array<double,1> g_old(n_features+1);
+    blitz::Array<double,1> w_old(n_features+1);
+    blitz::Array<double,1> g(n_features+1);
+    blitz::Array<double,1> w(n_features+1);
+    g_old = 0.;
+    w_old = 0.;
+    g = 0.;
+    w = 0.;
+
+    // Initialize working arrays
+    blitz::Array<double,1> s1(n_samples);
+    blitz::Array<double,1> u(n_features+1);
+    blitz::Array<double,1> tmp_n(n_samples);
+    blitz::Array<double,1> tmp_d(n_features+1);
+
+    // Iterates...
+    static const double ten_epsilon = 10*std::numeric_limits<double>::epsilon();
+    for(size_t iter=0; ; ++iter)
+    {
+      // 1. Computes the non-weighted version of the likelihood
+      // s1 = sum_{i=1}^{n}(1./(1.+exp(-y_i (w^T x_i + logit))
+      //   where - the x blitz::Array contains -y_i x_i values
+      //         - the offset blitz::Array contains -y_i logit values
+      s1 = 1. / (1. + blitz::exp(blitz::sum(w(j)*x(j,i), j) + offset));
+      // 2. Likelihood weighted by the prior/proportion
+      tmp_n = s1 * weights;
+      // 3. Gradient g of this weighted likelihood wrt. the weight vector w
+      bob::math::prod(x, tmp_n, g);
+      g -= m_lambda * w; // Regularization
+
+      // 4. Conjugate gradient step
+      if(iter == 0)
+        u = g;
+      else
+      {
+        tmp_d = (g-g_old);
+        double den = blitz::sum(u * tmp_d);
+        if(den == 0)
+          u = 0.;
+        else
+        {
+          // Hestenes-Stiefel formula: Heuristic to set the scale factor beta
+          //   (chosen as it works well in practice)
+          // beta = g^t(g-g_old) / (u_old^T (g - g_old))
+          double beta = blitz::sum(tmp_d * g) / den;
+          u = g - beta * u;
+        }
+      }
+
+      // 5. Line search along the direction u
+      // a. Compute ux
+      bob::math::prod(u,x,tmp_n);
+      // b. Compute u^T H u
+      //      = sum_{i} weights(i) sigmoid(w^T x_i) [1-sigmoid(w^T x_i)] (u^T x_i) + lambda u^T u
+      double uhu = blitz::sum(blitz::pow2(tmp_n) * weights * s1 * (1.-s1)) + m_lambda*blitz::sum(blitz::pow2(u));
+      // Terminates if uhu is close to zero
+      if(fabs(uhu) < ten_epsilon)
+      {
+        bob::core::info << "# CGLogReg Training terminated: convergence after " << iter << " iterations (u^T H u == 0)." << std::endl;
+        break;
+      }
+      // c. Compute w = w_old - (g^T u)/(u^T H u) u
+      w = w + blitz::sum(u*g) / uhu * u;
+
+      // Terminates if convergence has been reached
+      if(blitz::max(blitz::fabs(w-w_old)) <= m_convergence_threshold)
+      {
+        bob::core::info << "# CGLogReg Training terminated: convergence after " << iter << " iterations." << std::endl;
+        break;
+      }
+      // Terminates if maximum number of iterations has been reached
+      if(m_max_iterations > 0 && iter+1 >= m_max_iterations)
+      {
+        bob::core::info << "# CGLogReg terminated: maximum number of iterations (" << m_max_iterations << ") reached." << std::endl;
+        break;
+      }
+
+      // Backup previous values
+      g_old = g;
+      w_old = w;
+    }
+
+    // Updates the LinearMachine
+    machine.resize(n_features, 1);
+    machine.setInputSubtraction(mean);
+    machine.setInputDivision(std_dev);
+
+    blitz::Array<double,2>& w_ = machine.updateWeights();
+    w_(rall,0) = w(rd); // Weights: first D values
+    machine.setBiases(w(n_features)); // Bias: D+1 value
+  }
+
+}}}
diff --git a/xbob/learn/linear/cpp/machine.cpp b/xbob/learn/linear/cpp/machine.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..86be6859361c48a652e90dbb27a0f74da85d8849
--- /dev/null
+++ b/xbob/learn/linear/cpp/machine.cpp
@@ -0,0 +1,233 @@
+/**
+ * @date Tue May 31 09:22:33 2011 +0200
+ * @author Andre Anjos <andre.anjos@idiap.ch>
+ *
+ * @brief Implements a Machine
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#include <cmath>
+#include <boost/make_shared.hpp>
+#include <boost/format.hpp>
+
+#include <bob/core/array_copy.h>
+#include <bob/math/linear.h>
+
+#include <xbob.learn.linear/machine.h>
+
+namespace bob { namespace learn { namespace linear {
+
+  Machine::Machine(const blitz::Array<double,2>& weight)
+    : m_input_sub(weight.extent(0)),
+    m_input_div(weight.extent(0)),
+    m_bias(weight.extent(1)),
+    m_activation(boost::make_shared<bob::machine::IdentityActivation>()),
+    m_buffer(weight.extent(0))
+  {
+    m_input_sub = 0.0;
+    m_input_div = 1.0;
+    m_bias = 0.0;
+    m_weight.reference(bob::core::array::ccopy(weight));
+  }
+
+  Machine::Machine():
+    m_input_sub(0),
+    m_input_div(0),
+    m_weight(0, 0),
+    m_bias(0),
+    m_activation(boost::make_shared<bob::machine::IdentityActivation>()),
+    m_buffer(0)
+  {
+  }
+
+  Machine::Machine(size_t n_input, size_t n_output):
+    m_input_sub(n_input),
+    m_input_div(n_input),
+    m_weight(n_input, n_output),
+    m_bias(n_output),
+    m_activation(boost::make_shared<bob::machine::IdentityActivation>()),
+    m_buffer(n_input)
+  {
+    m_input_sub = 0.0;
+    m_input_div = 1.0;
+    m_weight = 0.0;
+    m_bias = 0.0;
+  }
+
+  Machine::Machine(const Machine& other):
+    m_input_sub(bob::core::array::ccopy(other.m_input_sub)),
+    m_input_div(bob::core::array::ccopy(other.m_input_div)),
+    m_weight(bob::core::array::ccopy(other.m_weight)),
+    m_bias(bob::core::array::ccopy(other.m_bias)),
+    m_activation(other.m_activation),
+    m_buffer(m_input_sub.shape())
+  {
+  }
+
+  Machine::Machine (bob::io::HDF5File& config) {
+    load(config);
+  }
+
+  Machine::~Machine() {}
+
+  Machine& Machine::operator=
+    (const Machine& other) {
+      if(this != &other)
+      {
+        m_input_sub.reference(bob::core::array::ccopy(other.m_input_sub));
+        m_input_div.reference(bob::core::array::ccopy(other.m_input_div));
+        m_weight.reference(bob::core::array::ccopy(other.m_weight));
+        m_bias.reference(bob::core::array::ccopy(other.m_bias));
+        m_activation = other.m_activation;
+        m_buffer.resize(m_input_sub.shape());
+      }
+      return *this;
+    }
+
+  bool Machine::operator==(const Machine& b) const {
+    return (bob::core::array::isEqual(m_input_sub, b.m_input_sub) &&
+        bob::core::array::isEqual(m_input_div, b.m_input_div) &&
+        bob::core::array::isEqual(m_weight, b.m_weight) &&
+        bob::core::array::isEqual(m_bias, b.m_bias) &&
+        m_activation->str() == b.m_activation->str());
+  }
+
+  bool Machine::operator!=(const Machine& b) const {
+    return !(this->operator==(b));
+  }
+
+  bool Machine::is_similar_to(const Machine& b, const double r_epsilon, const double a_epsilon) const {
+
+    return (bob::core::array::isClose(m_input_sub, b.m_input_sub, r_epsilon, a_epsilon) &&
+        bob::core::array::isClose(m_input_div, b.m_input_div, r_epsilon, a_epsilon) &&
+        bob::core::array::isClose(m_weight, b.m_weight, r_epsilon, a_epsilon) &&
+        bob::core::array::isClose(m_bias, b.m_bias, r_epsilon, a_epsilon) &&
+        m_activation->str() == b.m_activation->str());
+  }
+
+  void Machine::load (bob::io::HDF5File& config) {
+
+    //reads all data directly into the member variables
+    m_input_sub.reference(config.readArray<double,1>("input_sub"));
+    m_input_div.reference(config.readArray<double,1>("input_div"));
+    m_weight.reference(config.readArray<double,2>("weights"));
+    m_bias.reference(config.readArray<double,1>("biases"));
+    m_buffer.resize(m_input_sub.extent(0));
+
+    //switch between different versions - support for version 1
+    if (config.hasAttribute(".", "version")) { //new version
+      config.cd("activation");
+      m_activation = bob::machine::load_activation(config);
+      config.cd("..");
+    }
+    else { //old version
+      uint32_t act = config.read<uint32_t>("activation");
+      m_activation = bob::machine::make_deprecated_activation(act);
+    }
+
+  }
+
+  void Machine::resize (size_t input, size_t output) {
+
+    m_input_sub.resizeAndPreserve(input);
+    m_input_div.resizeAndPreserve(input);
+    m_buffer.resizeAndPreserve(input);
+    m_weight.resizeAndPreserve(input, output);
+    m_bias.resizeAndPreserve(output);
+
+  }
+
+  void Machine::save (bob::io::HDF5File& config) const {
+
+    config.setAttribute(".", "version", 1);
+    config.setArray("input_sub", m_input_sub);
+    config.setArray("input_div", m_input_div);
+    config.setArray("weights", m_weight);
+    config.setArray("biases", m_bias);
+    config.createGroup("activation");
+    config.cd("activation");
+    m_activation->save(config);
+    config.cd("..");
+
+  }
+
+  void Machine::forward_ (const blitz::Array<double,1>& input, blitz::Array<double,1>& output) const {
+
+    m_buffer = (input - m_input_sub) / m_input_div;
+    bob::math::prod_(m_buffer, m_weight, output);
+    for (int i=0; i<m_weight.extent(1); ++i)
+      output(i) = m_activation->f(output(i) + m_bias(i));
+
+  }
+
+  void Machine::forward (const blitz::Array<double,1>& input, blitz::Array<double,1>& output) const {
+
+    if (m_weight.extent(0) != input.extent(0)) { //checks input dimension
+      boost::format m("mismatch on the input dimension: expected a vector of size %d, but you input one with size = %d instead");
+      m % m_weight.extent(0) % input.extent(0);
+      throw std::runtime_error(m.str());
+    }
+    if (m_weight.extent(1) != output.extent(0)) { //checks output dimension
+      boost::format m("mismatch on the output dimension: expected a vector of size %d, but you input one with size = %d instead");
+      m % m_weight.extent(1) % output.extent(0);
+      throw std::runtime_error(m.str());
+    }
+    forward_(input, output);
+
+  }
+
+  void Machine::setWeights (const blitz::Array<double,2>& weight) {
+
+    if (weight.extent(0) != m_input_sub.extent(0)) { //checks 1st dimension
+      boost::format m("mismatch on the weight shape (number of rows): expected a weight matrix with %d row(s), but you input one with %d row(s) instead");
+      m % m_input_sub.extent(0) % weight.extent(0);
+      throw std::runtime_error(m.str());
+    }
+    if (weight.extent(1) != m_bias.extent(0)) { //checks 2nd dimension
+      boost::format m("mismatch on the weight shape (number of columns): expected a weight matrix with %d column(s), but you input one with %d column(s) instead");
+      m % m_bias.extent(0) % weight.extent(1);
+      throw std::runtime_error(m.str());
+    }
+    m_weight.reference(bob::core::array::ccopy(weight));
+
+  }
+
+  void Machine::setBiases (const blitz::Array<double,1>& bias) {
+
+    if (m_weight.extent(1) != bias.extent(0)) {
+      boost::format m("mismatch on the bias shape: expected a vector of size %d, but you input one with size = %d instead");
+      m % m_weight.extent(1) % bias.extent(0);
+      throw std::runtime_error(m.str());
+    }
+    m_bias.reference(bob::core::array::ccopy(bias));
+
+  }
+
+  void Machine::setInputSubtraction (const blitz::Array<double,1>& v) {
+
+    if (m_weight.extent(0) != v.extent(0)) {
+      boost::format m("mismatch on the input subtraction shape: expected a vector of size %d, but you input one with size = %d instead");
+      m % m_weight.extent(0) % v.extent(0);
+      throw std::runtime_error(m.str());
+    }
+    m_input_sub.reference(bob::core::array::ccopy(v));
+
+  }
+
+  void Machine::setInputDivision (const blitz::Array<double,1>& v) {
+
+    if (m_weight.extent(0) != v.extent(0)) {
+      boost::format m("mismatch on the input division shape: expected a vector of size %d, but you input one with size = %d instead");
+      m % m_weight.extent(0) % v.extent(0);
+      throw std::runtime_error(m.str());
+    }
+    m_input_div.reference(bob::core::array::ccopy(v));
+
+  }
+
+  void Machine::setActivation (boost::shared_ptr<bob::machine::Activation> a) {
+    m_activation = a;
+  }
+
+}}}
diff --git a/xbob/learn/linear/cpp/pca.cpp b/xbob/learn/linear/cpp/pca.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d113ca83d02e0df418061a3be6b65719938ada3b
--- /dev/null
+++ b/xbob/learn/linear/cpp/pca.cpp
@@ -0,0 +1,188 @@
+/**
+ * @date Tue Jan 18 17:07:26 2011 +0100
+ * @author André Anjos <andre.anjos@idiap.ch>
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ *
+ * @brief Principal Component Analysis implemented with Singular Value
+ * Decomposition or using the Covariance Method. Both are implemented using
+ * LAPACK. Implementation.
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#include <algorithm>
+#include <blitz/array.h>
+#include <boost/format.hpp>
+#include <bob/math/stats.h>
+#include <bob/math/svd.h>
+#include <bob/math/eig.h>
+
+#include <xbob.learn.linear/pca.h>
+
+namespace bob { namespace learn { namespace linear {
+
+  PCATrainer::PCATrainer(bool use_svd)
+    : m_use_svd(use_svd), m_safe_svd(false)
+  {
+  }
+
+  PCATrainer::PCATrainer(const PCATrainer& other)
+    : m_use_svd(other.m_use_svd), m_safe_svd(other.m_safe_svd)
+  {
+  }
+
+  PCATrainer::~PCATrainer() {}
+
+  PCATrainer& PCATrainer::operator= (const PCATrainer& other) {
+    if (this != &other) {
+      m_use_svd = other.m_use_svd;
+      m_safe_svd = other.m_safe_svd;
+    }
+    return *this;
+  }
+
+  bool PCATrainer::operator== (const PCATrainer& other) const {
+
+    return m_use_svd == other.m_use_svd &&
+      m_safe_svd == other.m_safe_svd;
+
+  }
+
+  bool PCATrainer::operator!= (const PCATrainer& other) const {
+
+    return !(this->operator==(other));
+
+  }
+
+  bool PCATrainer::is_similar_to (const PCATrainer& other, const double r_epsilon,
+      const double a_epsilon) const {
+    return this->operator==(other);
+  }
+
+  /**
+   * Sets up the machine calculating the PC's via the Covariance Matrix
+   */
+  static void pca_via_covmat(Machine& machine,
+      blitz::Array<double,1>& eigen_values, const blitz::Array<double,2>& X,
+      int rank) {
+
+    /**
+     * computes the covariance matrix (X-mu)(X-mu)^T / (len(X)-1) and then solves
+     * the generalized eigen-value problem taking into consideration the
+     * covariance matrix is symmetric (and, by extension, hermitian).
+     */
+    blitz::Array<double,1> mean(X.extent(1));
+    blitz::Array<double,2> Sigma(X.extent(1), X.extent(1));
+    bob::math::scatter_(X, Sigma, mean);
+    Sigma /= (X.extent(0)-1); //unbiased variance estimator
+
+    blitz::Array<double,2> U(X.extent(1), X.extent(1));
+    blitz::Array<double,1> e(X.extent(1));
+    bob::math::eigSym_(Sigma, U, e);
+    e.reverseSelf(0);
+    U.reverseSelf(1);
+
+    /**
+     * sets the linear machine with the results:
+     */
+    machine.setInputSubtraction(mean);
+    machine.setInputDivision(1.0);
+    machine.setBiases(0.0);
+    if (e.size() == eigen_values.size()) {
+      eigen_values = e;
+      machine.setWeights(U);
+    }
+    else {
+      eigen_values = e(blitz::Range(0,rank-1));
+      machine.setWeights(U(blitz::Range::all(), blitz::Range(0,rank-1)));
+    }
+
+  }
+
+  /**
+   * Sets up the machine calculating the PC's via SVD
+   */
+  static void pca_via_svd(Machine& machine, blitz::Array<double,1>& eigen_values,
+      const blitz::Array<double,2>& X, int rank, bool safe_svd) {
+
+    // removes the empirical mean from the training data
+    blitz::Array<double,2> data(X.extent(1), X.extent(0));
+    blitz::Range a = blitz::Range::all();
+    for (int i=0; i<X.extent(0); ++i) data(a,i) = X(i,a);
+
+    // computes the mean of the training data
+    blitz::secondIndex j;
+    blitz::Array<double,1> mean(X.extent(1));
+    mean = blitz::mean(data, j);
+
+    // applies the training data mean
+    for (int i=0; i<X.extent(0); ++i) data(a,i) -= mean;
+
+    /**
+     * computes the singular value decomposition using lapack
+     *
+     * note: Lapack already organizes the U,Sigma,V**T matrixes so that the
+     * singular values in Sigma are organized by decreasing order of magnitude.
+     * You **don't** need sorting after this.
+     */
+    const int rank_1 = (rank == (int)X.extent(1))? X.extent(1) : X.extent(0);
+    blitz::Array<double,2> U(X.extent(1), rank_1);
+    blitz::Array<double,1> sigma(rank_1);
+    bob::math::svd_(data, U, sigma, safe_svd);
+
+    /**
+     * sets the linear machine with the results:
+     *
+     * note: eigen values are sigma^2/X.extent(0) diagonal
+     *       eigen vectors are the rows of U
+     */
+    machine.setInputSubtraction(mean);
+    machine.setInputDivision(1.0);
+    machine.setBiases(0.0);
+    blitz::Range up_to_rank(0, rank-1);
+    machine.setWeights(U(a,up_to_rank));
+
+    //weight normalization (if necessary):
+    //norm_factor = blitz::sum(blitz::pow2(V(all,i)))
+
+    // finally, we set also the eigen values in this version
+    eigen_values = (blitz::pow2(sigma)/(X.extent(0)-1))(up_to_rank);
+  }
+
+  void PCATrainer::train(Machine& machine, blitz::Array<double,1>& eigen_values,
+      const blitz::Array<double,2>& X) const {
+
+    // data is checked now and conforms, just proceed w/o any further checks.
+    const int rank = output_size(X);
+
+    // Checks that the dimensions are matching
+    if (machine.inputSize() != (size_t)X.extent(1)) {
+      boost::format m("Number of features at input data set (%d columns) does not match machine input size (%d)");
+      m % X.extent(1) % machine.inputSize();
+      throw std::runtime_error(m.str());
+    }
+    if (machine.outputSize() != (size_t)rank) {
+      boost::format m("Number of outputs of the given machine (%d) does not match the maximum covariance rank, i.e., min(#samples-1,#features) = min(%d, %d) = %d");
+      m % machine.outputSize() % (X.extent(0)-1) % X.extent(1) % rank;
+      throw std::runtime_error(m.str());
+    }
+    if (eigen_values.extent(0) != rank) {
+      boost::format m("Number of eigenvalues on the given 1D array (%d) does not match the maximum covariance rank, i.e., min(#samples-1,#features) = min(%d,%d) = %d");
+      m % eigen_values.extent(0) % (X.extent(0)-1) % X.extent(1) % rank;
+      throw std::runtime_error(m.str());
+    }
+
+    if (m_use_svd) pca_via_svd(machine, eigen_values, X, rank, m_safe_svd);
+    else pca_via_covmat(machine, eigen_values, X, rank);
+  }
+
+  void PCATrainer::train(Machine& machine, const blitz::Array<double,2>& X) const {
+    blitz::Array<double,1> throw_away_eigen_values(output_size(X));
+    train(machine, throw_away_eigen_values, X);
+  }
+
+  size_t PCATrainer::output_size (const blitz::Array<double,2>& X) const {
+    return (size_t)std::min(X.extent(0)-1,X.extent(1));
+  }
+
+}}}
diff --git a/xbob/learn/linear/cpp/wccn.cpp b/xbob/learn/linear/cpp/wccn.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..80ca519549c4b488dd19adb0c5cb673f58a21149
--- /dev/null
+++ b/xbob/learn/linear/cpp/wccn.cpp
@@ -0,0 +1,104 @@
+/**
+ * @author Elie Khoury <Elie.Khoury@idiap.ch>
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @date Tue Apr 2 21:08:00 2013 +0200
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#include <boost/make_shared.hpp>
+#include <bob/math/inv.h>
+#include <bob/math/lu.h>
+#include <bob/math/stats.h>
+
+#include <xbob.learn.linear/wccn.h>
+
+namespace bob { namespace learn { namespace linear {
+
+  WCCNTrainer::WCCNTrainer() {
+  }
+
+  WCCNTrainer::WCCNTrainer(const WCCNTrainer& other) {
+  }
+
+  WCCNTrainer::~WCCNTrainer() {}
+
+  WCCNTrainer& WCCNTrainer::operator= (const WCCNTrainer& other) {
+    return *this;
+  }
+
+  bool WCCNTrainer::operator== (const WCCNTrainer& other) const {
+    return true;
+  }
+
+  bool WCCNTrainer::operator!= (const WCCNTrainer& other) const {
+    return !(this->operator==(other));
+  }
+
+  bool WCCNTrainer::is_similar_to (const WCCNTrainer& other,
+      const double r_epsilon, const double a_epsilon) const {
+    return true;
+  }
+
+
+  void WCCNTrainer::train(Machine& machine,
+      const std::vector<blitz::Array<double, 2> >& data) const {
+
+    const size_t n_classes = data.size();
+    // if #classes < 2, then throw
+    if (n_classes < 2) {
+      boost::format m("number of classes should be >= 2, but you passed %u");
+      m % n_classes;
+      throw std::runtime_error(m.str());
+    }
+
+    // checks for data type and shape once
+    const int n_features = data[0].extent(1);
+
+    for (size_t cl=0; cl<n_classes; ++cl) {
+      if (data[cl].extent(1) != n_features) {
+        boost::format m("number of features (columns) of array for class %u (%d) does not match that of array for class 0 (%d)");
+        m % cl % data[cl].extent(1) % n_features;
+        throw std::runtime_error(m.str());
+      }
+    }
+
+    // machine dimensions
+    const size_t n_inputs = machine.inputSize();
+    const size_t n_outputs = machine.outputSize();
+
+    // Checks that the dimensions are matching
+    if ((int)n_inputs != n_features) {
+      boost::format m("machine input size (%u) does not match the number of columns in input array (%d)");
+      m % n_inputs % n_features;
+      throw std::runtime_error(m.str());
+    }
+    if ((int)n_outputs != n_features) {
+      boost::format m("machine output size (%u) does not match the number of columns in output array (%d)");
+      m % n_outputs % n_features;
+      throw std::runtime_error(m.str());
+    }
+
+    // 1. Computes the mean vector and the Scatter matrix Sw and Sb
+    blitz::Array<double,1> mean(n_features);
+    blitz::Array<double,2> buf1(n_features, n_features); // Sw
+    blitz::Array<double,2> buf2(n_features, n_features); // Sb
+    bob::math::scatters(data, buf1, buf2, mean); // buf1 = Sw; buf2 = Sb
+
+    // 2. Computes the inverse of (1/N * Sw), Sw is the within-class covariance matrix
+    buf1 /= n_classes;
+    bob::math::inv(buf1, buf2); // buf2 = (1/N * Sw)^{-1}
+
+  // 3. Computes the Cholesky decomposition of the inverse covariance matrix
+  bob::math::chol(buf2, buf1); //  buf1 = cholesky(buf2)
+
+  // 4. Updates the linear machine
+  machine.setInputSubtraction(0); // we do not substract the mean
+  machine.setInputDivision(1.);
+  machine.setWeights(buf1);
+  machine.setBiases(0);
+  machine.setActivation(boost::make_shared<bob::machine::IdentityActivation>());
+
+  }
+
+}}}
diff --git a/xbob/learn/linear/cpp/whitening.cpp b/xbob/learn/linear/cpp/whitening.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1bda8e811df3bd7f6ecbc9f2eeacef5dfef31989
--- /dev/null
+++ b/xbob/learn/linear/cpp/whitening.cpp
@@ -0,0 +1,90 @@
+/**
+ * @date Tue Apr 2 21:08:00 2013 +0200
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#include <boost/make_shared.hpp>
+#include <bob/math/inv.h>
+#include <bob/math/lu.h>
+#include <bob/math/stats.h>
+
+#include <xbob.learn.linear/whitening.h>
+
+namespace bob { namespace learn { namespace linear {
+
+  WhiteningTrainer::WhiteningTrainer()
+  {
+  }
+
+  WhiteningTrainer::WhiteningTrainer(const WhiteningTrainer& other)
+  {
+  }
+
+  WhiteningTrainer::~WhiteningTrainer() {}
+
+  WhiteningTrainer& WhiteningTrainer::operator= (const WhiteningTrainer& other)
+  {
+    return *this;
+  }
+
+  bool WhiteningTrainer::operator== (const WhiteningTrainer& other) const
+  {
+    return true;
+  }
+
+  bool WhiteningTrainer::operator!= (const WhiteningTrainer& other) const
+  {
+    return !(this->operator==(other));
+  }
+
+  bool WhiteningTrainer::is_similar_to (const WhiteningTrainer& other,
+      const double r_epsilon, const double a_epsilon) const
+  {
+    return true;
+  }
+
+  void WhiteningTrainer::train(Machine& machine, const blitz::Array<double,2>& ar) const {
+    // training data dimensions
+    const size_t n_samples = ar.extent(0);
+    const size_t n_features = ar.extent(1);
+    // machine dimensions
+    const size_t n_inputs = machine.inputSize();
+    const size_t n_outputs = machine.outputSize();
+
+    // Checks that the dimensions are matching
+    if (n_inputs != n_features) {
+      boost::format m("machine input size (%u) does not match the number of columns in input array (%d)");
+      m % n_inputs % n_features;
+      throw std::runtime_error(m.str());
+    }
+    if (n_outputs != n_features) {
+      boost::format m("machine output size (%u) does not match the number of columns in output array (%d)");
+      m % n_outputs % n_features;
+      throw std::runtime_error(m.str());
+    }
+
+    // 1. Computes the mean vector and the covariance matrix of the training set
+    blitz::Array<double,1> mean(n_features);
+    blitz::Array<double,2> cov(n_features,n_features);
+    bob::math::scatter(ar, cov, mean);
+    cov /= (double)(n_samples-1);
+
+    // 2. Computes the inverse of the covariance matrix
+    blitz::Array<double,2> icov(n_features,n_features);
+    bob::math::inv(cov, icov);
+
+    // 3. Computes the Cholesky decomposition of the inverse covariance matrix
+    blitz::Array<double,2> whiten(n_features,n_features);
+    bob::math::chol(icov, whiten);
+
+    // 4. Updates the linear machine
+    machine.setInputSubtraction(mean);
+    machine.setInputDivision(1.);
+    machine.setWeights(whiten);
+    machine.setBiases(0);
+    machine.setActivation(boost::make_shared<bob::machine::IdentityActivation>());
+  }
+
+}}}
diff --git a/xbob/learn/linear/include/xbob.learn.linear/api.h b/xbob/learn/linear/include/xbob.learn.linear/api.h
index 4942935194f60784bc30e43bf08fce05c45db5d6..6edb7b6aa1927cd5e367cb4ce1b3874a6910f70f 100644
--- a/xbob/learn/linear/include/xbob.learn.linear/api.h
+++ b/xbob/learn/linear/include/xbob.learn.linear/api.h
@@ -10,9 +10,9 @@
 
 #include <Python.h>
 #include <xbob.learn.linear/config.h>
-#include <bob/machine/LinearMachine.h>
-#include <bob/trainer/PCATrainer.h>
-#include <bob/trainer/FisherLDATrainer.h>
+#include <xbob.learn.linear/machine.h>
+#include <xbob.learn.linear/pca.h>
+#include <xbob.learn.linear/lda.h>
 
 #define XBOB_LEARN_LINEAR_MODULE_PREFIX xbob.learn.linear
 #define XBOB_LEARN_LINEAR_MODULE_NAME _library
@@ -50,7 +50,7 @@ enum _PyBobLearnLinear_ENUM{
 
 typedef struct {
   PyObject_HEAD
-  bob::machine::LinearMachine* cxx;
+  bob::learn::linear::Machine* cxx;
 } PyBobLearnLinearMachineObject;
 
 #define PyBobLearnLinearMachine_Type_TYPE PyTypeObject
@@ -67,7 +67,7 @@ typedef struct {
 
 typedef struct {
   PyObject_HEAD
-  bob::trainer::PCATrainer* cxx;
+  bob::learn::linear::PCATrainer* cxx;
 } PyBobLearnLinearPCATrainerObject;
 
 #define PyBobLearnLinearPCATrainer_Type_TYPE PyTypeObject
@@ -81,7 +81,7 @@ typedef struct {
 
 typedef struct {
   PyObject_HEAD
-  bob::trainer::FisherLDATrainer* cxx;
+  bob::learn::linear::FisherLDATrainer* cxx;
 } PyBobLearnLinearFisherLDATrainerObject;
 
 #define PyBobLearnLinearFisherLDATrainer_Type_TYPE PyTypeObject
diff --git a/xbob/learn/linear/include/xbob.learn.linear/lda.h b/xbob/learn/linear/include/xbob.learn.linear/lda.h
new file mode 100644
index 0000000000000000000000000000000000000000..0b23c8fb67367565f3f4bce7073f6f259dbd7b55
--- /dev/null
+++ b/xbob/learn/linear/include/xbob.learn.linear/lda.h
@@ -0,0 +1,149 @@
+/**
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @author Andre Anjos <andre.anjos@idiap.ch>
+ * @date Sat Jun 4 21:38:59 2011 +0200
+ *
+ * @brief Implements a multi-class Fisher/LDA linear machine Training using
+ * Singular Value Decomposition (SVD). For more information on Linear Machines
+ * and associated methods, please consult Bishop, Machine Learning and Pattern
+ * Recognition chapter 4.
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#ifndef BOB_LEARN_LINEAR_FISHER_H
+#define BOB_LEARN_LINEAR_FISHER_H
+
+#include <vector>
+#include "machine.h"
+
+namespace bob { namespace learn { namespace linear {
+
+  /**
+   * @brief Sets a linear machine to perform the Fisher/LDA decomposition.\n
+   *
+   * References:\n
+   * 1. Bishop, Machine Learning and Pattern Recognition chapter 4.\n
+   * 2. http://en.wikipedia.org/wiki/Linear_discriminant_analysis
+   */
+  class FisherLDATrainer {
+    public: //api
+
+      /**
+       * @brief Initializes a new Fisher/LDA trainer. The training stage will
+       * place the resulting fisher components in the linear machine and set it
+       * up to extract the variable means automatically.
+       *
+       * @param strip_to_rank Specifies what is the final size of the prepared
+       * LinearMachine. The default setting (<code>true</code>), makes the
+       * trainer return only the K-1 eigen-values/vectors limiting the output
+       * to the rank of Sw^1 Sb. If you set this value to <code>false</code>,
+       * the it returns all eigen-values/vectors of Sw^1 Sb, including the ones
+       * that are supposed to be zero.
+       *
+       * @param use_pinv If set (to <code>true</code>), then use the
+       * pseudo-inverse to compute Sw^-1 and then a generalized eigen-value
+       * decomposition instead of using the default (more stable) version of
+       * the eigen-value decomposition, starting from Sb and Sw.
+       */
+      FisherLDATrainer(bool use_pinv = false, bool strip_to_rank = true);
+
+      /**
+       * @brief Destructor
+       */
+      virtual ~FisherLDATrainer();
+
+      /**
+       * @brief Copy constructor.
+       */
+      FisherLDATrainer(const FisherLDATrainer& other);
+
+      /**
+       * @brief Assignment operator
+       */
+      FisherLDATrainer& operator=(const FisherLDATrainer& other);
+
+      /**
+       * @brief Equal to
+       */
+      bool operator==(const FisherLDATrainer& other) const;
+
+      /**
+       * @brief Not equal to
+       */
+      bool operator!=(const FisherLDATrainer& other) const;
+
+      /**
+       * @brief Similar to
+       */
+      bool is_similar_to(const FisherLDATrainer& other,
+          const double r_epsilon=1e-5, const double a_epsilon=1e-8) const;
+
+      /**
+       * @brief Gets the pseudo-inverse flag
+       */
+      bool getUsePseudoInverse () const { return m_use_pinv; }
+
+      /**
+       * @brief Sets the pseudo-inverse flag
+       */
+      void setUsePseudoInverse (bool v) { m_use_pinv = v; }
+
+      /**
+       * @brief Gets the strip-to-rank flag
+       */
+      bool getStripToRank () const { return m_strip_to_rank; }
+
+      /**
+       * @brief Sets the strip-to-rank flag
+       */
+      void setStripToRank (bool v) { m_strip_to_rank = v; }
+
+      /**
+       * @brief Trains the LinearMachine to perform Fisher/LDA discrimination.
+       * The resulting machine will have the eigen-vectors of the
+       * Sigma-1 * Sigma_b product, arranged by decreasing energy.
+       *
+       * Each input arrayset represents data from a given input class.
+       *
+       * Note we set only the N-1 eigen vectors in the linear machine since the
+       * last eigen value should be zero anyway. You can compress the machine
+       * output further using resize() if necessary.
+       */
+      void train(Machine& machine,
+          const std::vector<blitz::Array<double,2> >& X) const;
+
+      /**
+       * @brief Trains the LinearMachine to perform Fisher/LDA discrimination.
+       * The resulting machine will have the eigen-vectors of the
+       * Sigma-1 * Sigma_b product, arranged by decreasing energy. You don't
+       * need to sort the results. Also returns the eigen values of the
+       * covariance matrix product so you can use that to choose which
+       * components to keep.
+       *
+       * Each input arrayset represents data from a given input class.
+       *
+       * Note we set only the N-1 eigen vectors in the linear machine since the
+       * last eigen value should be zero anyway. You can compress the machine
+       * output further using resize() if necessary.
+       */
+      void train(Machine& machine, blitz::Array<double,1>& eigen_values,
+          const std::vector<blitz::Array<double,2> >& X) const;
+
+      /**
+       * @brief Returns the expected size of the output given the data.
+       *
+       * This number could be either K-1 (K = number of classes) or the number
+       * of columns (features) in X, depending on the seetting of
+       * <code>strip_to_rank</code>.
+       */
+      size_t output_size(const std::vector<blitz::Array<double,2> >& X) const;
+
+    private:
+      bool m_use_pinv; ///< use the 'pinv' method for LDA
+      bool m_strip_to_rank; ///< return rank or full matrix
+  };
+
+}}}
+
+#endif /* BOB_LEARN_LINEAR_FISHER_H */
diff --git a/xbob/learn/linear/include/xbob.learn.linear/logreg.h b/xbob/learn/linear/include/xbob.learn.linear/logreg.h
new file mode 100644
index 0000000000000000000000000000000000000000..ae4df85fb0f30fa7839ee60e6d9b2ddc8ac3672c
--- /dev/null
+++ b/xbob/learn/linear/include/xbob.learn.linear/logreg.h
@@ -0,0 +1,123 @@
+/**
+ * @author Laurent El Shafey <laurent.el-shafey@idiap.ch>
+ * @date Sat Sep 1 19:16:00 2012 +0100
+ *
+ * @brief Linear Logistic Regression trainer using a conjugate gradient
+ * approach.
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#ifndef BOB_LEARN_LINEAR_LOGREG_H
+#define BOB_LEARN_LINEAR_LOGREG_H
+
+#include <boost/format.hpp>
+#include "machine.h"
+
+namespace bob { namespace learn { namespace linear {
+
+  /**
+   * Trains a Linear Logistic Regression model using a conjugate gradient
+   * approach. The objective function is normalized with respect to the
+   * proportion of elements in class 1 to the ones in class 2, and
+   * then weighted with respect to a given synthetic prior, P, as this is
+   * done in the FoCal toolkit.
+   * References:
+   *   1/ "A comparison of numerical optimizers for logistic regression",
+   *   T. Minka, Unpublished draft, 2003 (revision in 2007),
+   *   http://research.microsoft.com/en-us/um/people/minka/papers/logreg/
+   *   2/ FoCal, http://www.dsp.sun.ac.za/~nbrummer/focal/
+   */
+  class CGLogRegTrainer {
+
+    public: //api
+
+      /**
+       * Default constructor.
+       * @param prior The synthetic prior. It should be in the range ]0.,1.[
+       * @param convergence_threshold The threshold to detect the convergence
+       *           of the iterative conjugate gradient algorithm
+       * @param max_iterations The maximum number of iterations of the
+       *           iterative conjugate gradient algorithm (0 <-> infinity)
+       * @param lambda The regularization factor
+       * @param mean_std_norm Compute mean and standard deviation in training
+       *           data and set the input_subtract and input_divide parameters
+       *           of the resulting machine
+       */
+      CGLogRegTrainer(const double prior=0.5,
+        const double convergence_threshold=1e-5,
+        const size_t max_iterations=10000,
+        const double lambda=0.,
+        const bool mean_std_norm=false);
+
+      /**
+       * Copy constructor
+       */
+      CGLogRegTrainer(const CGLogRegTrainer& other);
+
+      /**
+       * Destructor
+       */
+      virtual ~CGLogRegTrainer();
+
+      /**
+       * Assignment operator
+       */
+      CGLogRegTrainer& operator=(const CGLogRegTrainer& other);
+
+      /**
+       * @brief Equal to
+       */
+      bool operator==(const CGLogRegTrainer& b) const;
+      /**
+       * @brief Not equal to
+       */
+      bool operator!=(const CGLogRegTrainer& b) const;
+
+      /**
+       * Getters
+       */
+      double getPrior() const { return m_prior; }
+      double getConvergenceThreshold() const { return m_convergence_threshold; }
+      size_t getMaxIterations() const { return m_max_iterations; }
+      double getLambda() const { return m_lambda; }
+      bool getNorm() const { return m_mean_std_norm; }
+
+      /**
+       * Setters
+       */
+      void setPrior(const double prior)
+      { if(prior<=0. || prior>=1.)
+        {
+          boost::format m("Prior (%f) not in the range ]0,1[.");
+          m % prior;
+          throw std::runtime_error(m.str());
+        }
+        m_prior = prior; }
+      void setConvergenceThreshold(const double convergence_threshold)
+      { m_convergence_threshold = convergence_threshold; }
+      void setMaxIterations(const size_t max_iterations)
+      { m_max_iterations = max_iterations; }
+      void setLambda(const double lambda)
+      { m_lambda = lambda; }
+      void setNorm(const bool mean_std_norm) { m_mean_std_norm = mean_std_norm; }
+
+      /**
+       * Trains the LinearMachine to perform Linear Logistic Regression
+       */
+      virtual void train(Machine& machine,
+          const blitz::Array<double,2>& negatives,
+          const blitz::Array<double,2>& positives) const;
+
+    private:
+      // Attributes
+      double m_prior;
+      double m_convergence_threshold;
+      size_t m_max_iterations;
+      double m_lambda;
+      bool m_mean_std_norm;
+  };
+
+}}}
+
+#endif /* BOB_LEARN_LINEAR_LOGREG_H */
diff --git a/xbob/learn/linear/include/xbob.learn.linear/machine.h b/xbob/learn/linear/include/xbob.learn.linear/machine.h
new file mode 100644
index 0000000000000000000000000000000000000000..d6addcd0c36abfdc969c2b7ec6fc72fa9071b8b6
--- /dev/null
+++ b/xbob/learn/linear/include/xbob.learn.linear/machine.h
@@ -0,0 +1,273 @@
+/**
+ * @author Andre Anjos <andre.anjos@idiap.ch>
+ * @date Tue May 31 09:22:33 2011 +0200
+ *
+ * A machine that implements the liner projection of input to the output using
+ * weights, biases and sums:
+ * output = sum(inputs * weights) + bias
+ * It is possible to setup the machine to previously normalize the input taking
+ * into consideration some input bias and division factor. It is also possible
+ * to set it up to have an activation function.
+ * A linear classifier. See C. M. Bishop, "Pattern Recognition and Machine
+ * Learning", chapter 4
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#ifndef BOB_LEARN_LINEAR_MACHINE_H
+#define BOB_LEARN_LINEAR_MACHINE_H
+
+#include <blitz/array.h>
+#include <bob/io/HDF5File.h>
+#include <bob/machine/Activation.h>
+
+namespace bob { namespace learn { namespace linear {
+
+  /**
+   * A linear classifier. See C. M. Bishop, "Pattern Recognition and Machine
+   * Learning", chapter 4 for more details.
+   */
+  class Machine {
+
+    public: //api
+
+      /**
+       * Default constructor. Builds an otherwise invalid 0 x 0 linear machine.
+       * This is equivalent to construct a Machine with two size_t
+       * parameters set to 0, as in Machine(0, 0).
+       */
+      Machine ();
+
+      /**
+       * Constructor, builds a new Linear machine. Weights and biases are
+       * not initialized.
+       *
+       * @param input Size of input vector
+       * @param output Size of output vector
+       */
+      Machine (size_t input, size_t output);
+
+      /**
+       * Builds a new machine with a set of weights. Each column of the weight
+       * matrix should represent a direction in which the input is projected
+       * to.
+       */
+      Machine(const blitz::Array<double,2>& weight);
+
+      /**
+       * Copies another machine
+       */
+      Machine (const Machine& other);
+
+      /**
+       * Starts a new Machine from an existing Configuration object.
+       */
+      Machine (bob::io::HDF5File& config);
+
+      /**
+       * Just to virtualise the destructor
+       */
+      virtual ~Machine();
+
+      /**
+       * Assigns from a different machine
+       */
+      Machine& operator= (const Machine& other);
+
+      /**
+        * @brief Equal to
+        */
+      bool operator==(const Machine& b) const;
+      /**
+        * @brief Not equal to
+        */
+      bool operator!=(const Machine& b) const;
+      /**
+       * @brief Similar to.
+       */
+      bool is_similar_to(const Machine& b, const double r_epsilon=1e-5,
+        const double a_epsilon=1e-8) const;
+
+      /**
+       * Loads data from an existing configuration object. Resets the current
+       * state.
+       */
+      void load (bob::io::HDF5File& config);
+
+      /**
+       * Saves an existing machine to a Configuration object.
+       */
+      void save (bob::io::HDF5File& config) const;
+
+      /**
+       * Forwards data through the network, outputs the values of each linear
+       * component the input signal is decomposed at.
+       *
+       * The input and output are NOT checked for compatibility each time. It
+       * is your responsibility to do it.
+       */
+      void forward_ (const blitz::Array<double,1>& input,
+          blitz::Array<double,1>& output) const;
+
+      /**
+       * Forwards data through the network, outputs the values of each linear
+       * component the input signal is decomposed at.
+       *
+       * The input and output are checked for compatibility each time the
+       * forward method is applied.
+       */
+      void forward (const blitz::Array<double,1>& input,
+          blitz::Array<double,1>& output) const;
+
+      /**
+       * Resizes the machine. If either the input or output increases in size,
+       * the weights and other factors should be considered uninitialized. If
+       * the size is preserved or reduced, already initialized values will not
+       * be changed.
+       *
+       * Tip: Use this method to force data compression. All will work out
+       * given most relevant factors to be preserved are organized on the top
+       * of the weight matrix. In this way, reducing the system size will
+       * supress less relevant projections.
+       */
+      void resize(size_t n_input, size_t n_output);
+
+      /**
+       * Returns the number of inputs expected by this machine
+       */
+      inline size_t inputSize () const { return m_weight.extent(0); }
+
+      /**
+       * Returns the number of outputs generated by this machine
+       */
+      inline size_t outputSize () const { return m_weight.extent(1); }
+
+      /**
+       * Returns the input subtraction factor
+       */
+      inline const blitz::Array<double, 1>& getInputSubtraction() const
+      { return m_input_sub; }
+
+      /**
+       * Sets the current input subtraction factor. We will check that the
+       * number of inputs (first dimension of weights) matches the number of
+       * values currently set and will raise an exception if that is not the
+       * case.
+       */
+      void setInputSubtraction(const blitz::Array<double,1>& v);
+
+      /**
+       * Returns the current input subtraction factor in order to be updated.
+       * @warning Use with care. Only trainers should use this function for
+       * efficiency reasons.
+       */
+      inline blitz::Array<double, 1>& updateInputSubtraction()
+      { return m_input_sub; }
+
+      /**
+       * Sets all input subtraction values to a specific value.
+       */
+      inline void setInputSubtraction(double v) { m_input_sub = v; }
+
+      /**
+       * Returns the input division factor
+       */
+      inline const blitz::Array<double, 1>& getInputDivision() const
+      { return m_input_div; }
+
+      /**
+       * Sets the current input division factor. We will check that the number
+       * of inputs (first dimension of weights) matches the number of values
+       * currently set and will raise an exception if that is not the case.
+       */
+      void setInputDivision(const blitz::Array<double,1>& v);
+
+      /**
+       * Returns the current input division factor in order to be updated.
+       * @warning Use with care. Only trainers should use this function for
+       * efficiency reasons.
+       */
+      inline blitz::Array<double, 1>& updateInputDivision()
+      { return m_input_div; }
+
+
+      /**
+       * Sets all input division values to a specific value.
+       */
+      inline void setInputDivision(double v) { m_input_div = v; }
+
+      /**
+       * Returns the current weight representation. Each column should be
+       * considered as a vector from which each of the output values is derived
+       * by projecting the input onto such a vector.
+       */
+      inline const blitz::Array<double, 2>& getWeights() const
+      { return m_weight; }
+
+      /**
+       * Sets the current weights. We will check that the number of outputs and
+       * inputs matches the expected values inside this machine and will raise
+       * an exception if that is not the case.
+       */
+      void setWeights(const blitz::Array<double,2>& weight);
+
+      /**
+       * Returns the current weight representation in order to be updated.
+       * Each column should be considered as a vector from which each of the
+       * output values is derived by projecting the input onto such a vector.
+       * @warning Use with care. Only trainers should use this function for
+       * efficiency reasons.
+       */
+      inline blitz::Array<double, 2>& updateWeights()
+      { return m_weight; }
+
+      /**
+       * Sets all weights to a single specific value.
+       */
+      inline void setWeights(double v) { m_weight = v; }
+
+      /**
+       * Returns the biases of this classifier.
+       */
+      inline const blitz::Array<double, 1>& getBiases() const
+      { return m_bias; }
+
+      /**
+       * Sets the current biases. We will check that the number of biases
+       * matches the number of weights (first dimension) currently set and
+       * will raise an exception if that is not the case.
+       */
+      void setBiases(const blitz::Array<double,1>& bias);
+
+      /**
+       * Sets all output bias values to a specific value.
+       */
+      inline void setBiases(double v) { m_bias = v; }
+
+      /**
+       * Returns the currently set activation function
+       */
+      inline boost::shared_ptr<bob::machine::Activation> getActivation() const { return m_activation; }
+
+      /**
+       * Sets the activation function for each of the outputs.
+       */
+      void setActivation(boost::shared_ptr<bob::machine::Activation> a);
+
+    private: //representation
+
+      typedef double (*actfun_t)(double); ///< activation function type
+
+      blitz::Array<double, 1> m_input_sub; ///< input subtraction
+      blitz::Array<double, 1> m_input_div; ///< input division
+      blitz::Array<double, 2> m_weight; ///< weights
+      blitz::Array<double, 1> m_bias; ///< biases for the output
+      boost::shared_ptr<bob::machine::Activation> m_activation; ///< currently set activation type
+
+      mutable blitz::Array<double, 1> m_buffer; ///< a buffer for speed
+
+  };
+
+}}}
+
+#endif /* BOB_LEARN_LINEAR_MACHINE_H */
diff --git a/xbob/learn/linear/include/xbob.learn.linear/pca.h b/xbob/learn/linear/include/xbob.learn.linear/pca.h
new file mode 100644
index 0000000000000000000000000000000000000000..65e17546c88aade89a6caae7e72315c9dbefc34c
--- /dev/null
+++ b/xbob/learn/linear/include/xbob.learn.linear/pca.h
@@ -0,0 +1,145 @@
+/**
+ * @author André Anjos <andre.anjos@idiap.ch>
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @date Tue Jan 18 17:07:26 2011 +0100
+ *
+ * @brief Principal Component Analysis implemented with Singular Value
+ * Decomposition or using the Covariance Method. Both are implemented using
+ * LAPACK.
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#ifndef BOB_LEARN_LINEAR_PCA_H
+#define BOB_LEARN_LINEAR_PCA_H
+
+#include "machine.h"
+
+namespace bob { namespace learn { namespace linear {
+
+  /**
+   * @brief Sets a linear machine to perform the Karhunen-Loève Transform
+   * (KLT) on a given dataset using either Singular Value Decomposition (SVD),
+   * the default, or the Covariance Method.
+   *
+   * References:
+   * 1. Eigenfaces for Recognition, Turk & Pentland, Journal of Cognitive
+   *    Neuroscience (1991) Volume: 3, Issue: 1, Publisher: MIT Press,
+   *    Pages: 71-86
+   * 2. http://en.wikipedia.org/wiki/Singular_value_decomposition
+   * 3. http://en.wikipedia.org/wiki/Principal_component_analysis
+   */
+  class PCATrainer {
+
+    public: //api
+
+      /**
+       * @brief Initializes a new PCA trainer using SVD or, as an option, the
+       * Covariance Matrix method. The training stage will place the resulting
+       * principal components in the linear machine and set it up to extract
+       * the variable means automatically. As an option, you may preset the
+       * trainer so that the normalization performed by the resulting linear
+       * machine also divides the variables by the standard deviation of each
+       * variable ensemble.
+       */
+      PCATrainer(bool use_svd=true);
+
+      /**
+       * @brief Copy constructor
+       */
+      PCATrainer(const PCATrainer& other);
+
+      /**
+       * @brief Destructor
+       */
+      virtual ~PCATrainer();
+
+      /**
+       * @brief Assignment operator
+       */
+      PCATrainer& operator=(const PCATrainer& other);
+
+      /**
+       * @brief Equal to
+       */
+      bool operator==(const PCATrainer& other) const;
+
+      /**
+       * @brief Not equal to
+       */
+      bool operator!=(const PCATrainer& other) const;
+
+      /**
+       * @brief Gets the SVD/CovMat flag. <code>true</code> means use SVD
+       * method.
+       */
+      bool getUseSVD () const { return m_use_svd; }
+
+      /**
+       * @brief Sets the SVD/CovMat flag. <code>true</code> means use SVD
+       * method.
+       */
+      void setUseSVD (bool value) { m_use_svd = value; }
+
+      /**
+       * @brief Gets the SVD LAPACK function used. <code>true</code> means
+       * use the LAPACK DGESVD method and false the DGESDD method.
+       * To be useful, the SVD method should be used by enabling the UseSVD
+       * flag.
+       */
+      bool getSafeSVD () const { return m_safe_svd; }
+
+      /**
+       * @brief Sets the SVD LAPACK function used. <code>true</code> means
+       * use the LAPACK DGESVD method and false the DGESDD method..
+       * To be useful, the SVD method should be used by enabling the UseSVD
+       * flag.
+       */
+      void setSafeSVD (bool value) { m_safe_svd = value; }
+
+      /**
+       * @brief Similar to
+       */
+      bool is_similar_to(const PCATrainer& other,
+          const double r_epsilon=1e-5, const double a_epsilon=1e-8) const;
+
+      /**
+       * @brief Trains the LinearMachine to perform the KLT. The resulting
+       * machine will have the eigen-vectors of the covariance matrix arranged
+       * by decreasing energy automatically. You don't need to sort the results.
+       */
+      virtual void train(Machine& machine,
+          const blitz::Array<double,2>& X) const;
+
+      /**
+       * @brief Trains the LinearMachine to perform the KLT. The resulting
+       * machien will have the eigen-vectors of the covariance matrix arranged
+       * by decreasing energy automatically. You don't need to sort the results.
+       * Also returns the eigen values of the covariance matrix so you can use
+       * that to choose which components to keep.
+       */
+      virtual void train(Machine& machine,
+          blitz::Array<double,1>& eigen_values,
+          const blitz::Array<double,2>& X) const;
+
+      /**
+       * @brief Calculates the maximum possible rank for the covariance matrix
+       * of X, given X.
+       *
+       * This determines what is the maximum number of non-zero eigen values
+       * that can be generated by this trainer. It should be used to setup
+       * Machines and input vectors prior to feeding them into this trainer.
+       */
+      size_t output_size(const blitz::Array<double,2>& X) const;
+
+    private: //representation
+
+      bool m_use_svd; ///< if this trainer should be using SVD or Covariance
+      bool m_safe_svd; ///< if svd is set, tells which LAPACK function to use
+                       ///  among dgesdd (false) and dgesvd (true)
+
+  };
+
+}}}
+
+#endif /* BOB_LEARN_LINEAR_PCA_H */
diff --git a/xbob/learn/linear/include/xbob.learn.linear/wccn.h b/xbob/learn/linear/include/xbob.learn.linear/wccn.h
new file mode 100644
index 0000000000000000000000000000000000000000..e9baba963e0bfc841148f6d5407ff07352e5d064
--- /dev/null
+++ b/xbob/learn/linear/include/xbob.learn.linear/wccn.h
@@ -0,0 +1,79 @@
+/**
+ * @author Elie Khoury <Elie.Khoury@idiap.ch>
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @date Tue Apr 9 22:10:00 2013 +0200
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#ifndef BOB_LEARN_LINEAR_WCCN_H
+#define BOB_LEARN_LINEAR_WCCN_H
+
+#include "machine.h"
+
+namespace bob { namespace learn { namespace linear {
+
+  /**
+   * @brief Sets a linear machine to perform a WCCN transform\n
+   *
+   * Reference 1:\n
+   * "Within-class covariance normalization for SVM-based speaker recognition",
+   *   Andrew O. Hatch, Sachin Kajarekar, and Andreas Stolcke,
+   *   ICSLP, 2006\n
+   * Reference 2:\n
+   * "N. Dehah, P. Kenny, R. Dehak, P. Dumouchel, P. Ouellet",
+   *   Front-end factor analysis for speaker verification,
+   *   IEEE TASLP, 2011\n
+   *
+   * Given a training set X, this will compute the W matrix such that:\n
+   *   \f$W = cholesky(inv(cov(X_{n},X_{n}^{T})))\f$, where \f$X_{n}\f$
+   *   corresponds to the center data
+   */
+  class WCCNTrainer {
+
+    public: //api
+
+      /**
+       * @brief Initializes a new WCCN trainer.
+       */
+      WCCNTrainer();
+
+      /**
+       * @brief Copy constructor
+       */
+      WCCNTrainer(const WCCNTrainer& other);
+
+      /**
+       * @brief Destructor
+       */
+      virtual ~WCCNTrainer();
+
+      /**
+       * @brief Assignment operator
+       */
+      WCCNTrainer& operator=(const WCCNTrainer& other);
+
+      /**
+       * @brief Equal to
+       */
+      bool operator==(const WCCNTrainer& other) const;
+      /**
+       * @brief Not equal to
+       */
+      bool operator!=(const WCCNTrainer& other) const;
+      /**
+       * @brief Similar to
+       */
+      bool is_similar_to(const WCCNTrainer& other, const double r_epsilon=1e-5,
+          const double a_epsilon=1e-8) const;
+
+      /**
+       * @brief Trains the LinearMachine to perform the WCCN
+       */
+      virtual void train(Machine& machine, const std::vector<blitz::Array<double, 2>>& data) const;
+
+  };
+
+}}}
+
+#endif /* BOB_LEARN_LINEAR_WCCN_H */
diff --git a/xbob/learn/linear/include/xbob.learn.linear/whitening.h b/xbob/learn/linear/include/xbob.learn.linear/whitening.h
new file mode 100644
index 0000000000000000000000000000000000000000..93103f8c2914620aad51eb211151f1f40b296792
--- /dev/null
+++ b/xbob/learn/linear/include/xbob.learn.linear/whitening.h
@@ -0,0 +1,73 @@
+/**
+ * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
+ * @date Tue Apr 2 21:04:00 2013 +0200
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#ifndef BOB_LEARN_LINEAR_WHITENING_H
+#define BOB_LEARN_LINEAR_WHITENING_H
+
+#include "machine.h"
+
+namespace bob { namespace learn { namespace linear {
+
+  /**
+   * @brief Sets a linear machine to perform a Whitening transform\n
+   *
+   * Reference:\n
+   * "Independent Component Analysis: Algorithms and Applications",
+   *   Aapo Hyvärinen, Erkki Oja,
+   *   Neural Networks, 2000, vol. 13, p. 411--430\n
+   *
+   * Given a training set X, this will compute the W matrix such that:\n
+   *   \f$W = cholesky(inv(cov(X_{n},X_{n}^{T})))\f$, where \f$X_{n}\f$
+   *   corresponds to the center data
+   */
+  class WhiteningTrainer {
+
+    public: //api
+
+      /**
+       * @brief Initializes a new Whitening trainer.
+       */
+      WhiteningTrainer();
+
+      /**
+       * @brief Copy constructor
+       */
+      WhiteningTrainer(const WhiteningTrainer& other);
+
+      /**
+       * @brief Destructor
+       */
+      virtual ~WhiteningTrainer();
+
+      /**
+       * @brief Assignment operator
+       */
+      WhiteningTrainer& operator=(const WhiteningTrainer& other);
+
+      /**
+       * @brief Equal to
+       */
+      bool operator==(const WhiteningTrainer& other) const;
+      /**
+       * @brief Not equal to
+       */
+      bool operator!=(const WhiteningTrainer& other) const;
+      /**
+       * @brief Similar to
+       */
+      bool is_similar_to(const WhiteningTrainer& other, const double r_epsilon=1e-5,
+          const double a_epsilon=1e-8) const;
+
+      /**
+       * @brief Trains the LinearMachine to perform the Whitening
+       */
+      virtual void train(Machine& machine, const blitz::Array<double,2>& data) const;
+  };
+
+}}}
+
+#endif /* BOB_LEARN_LINEAR_WHITENING_H */
diff --git a/xbob/learn/linear/lda.cpp b/xbob/learn/linear/lda.cpp
index 320c27adda42543477fd7307b66c8f0f9266216d..c37c0b110395aca034cf9422dd7613067146e896 100644
--- a/xbob/learn/linear/lda.cpp
+++ b/xbob/learn/linear/lda.cpp
@@ -10,7 +10,6 @@
 #define XBOB_LEARN_LINEAR_MODULE
 #include <xbob.blitz/cppapi.h>
 #include <xbob.blitz/cleanup.h>
-#include <bob/trainer/FisherLDATrainer.h>
 #include <xbob.learn.linear/api.h>
 #include <structmember.h>
 
@@ -146,7 +145,7 @@ static int PyBobLearnLinearFisherLDATrainer_init_bools
   if (strip_to_rank_ == -1) return -1;
 
   try {
-    self->cxx = new bob::trainer::FisherLDATrainer(use_pinv_?true:false,
+    self->cxx = new bob::learn::linear::FisherLDATrainer(use_pinv_?true:false,
         strip_to_rank_?true:false);
   }
   catch (std::exception& ex) {
@@ -177,7 +176,7 @@ static int PyBobLearnLinearFisherLDATrainer_init_copy
   auto copy = reinterpret_cast<PyBobLearnLinearFisherLDATrainerObject*>(other);
 
   try {
-    self->cxx = new bob::trainer::FisherLDATrainer(*(copy->cxx));
+    self->cxx = new bob::learn::linear::FisherLDATrainer(*(copy->cxx));
   }
   catch (std::exception& ex) {
     PyErr_SetString(PyExc_RuntimeError, ex.what());
@@ -366,7 +365,7 @@ static PyObject* PyBobLearnLinearFisherLDATrainer_Train
   if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|O!", kwlist,
         &X, &PyBobLearnLinearMachine_Type, &machine)) return 0;
 
-  /** 
+  /**
   // Note: strangely, if you pass dict.values(), this check does not work
   if (!PyIter_Check(X)) {
     PyErr_Format(PyExc_TypeError, "`%s' requires an iterable for parameter `X', but you passed `%s' which does not implement the iterator protocol", Py_TYPE(self)->tp_name, Py_TYPE(X)->tp_name);
@@ -384,7 +383,7 @@ static PyObject* PyBobLearnLinearFisherLDATrainer_Train
 
   while (PyObject* item = PyIter_Next(iterator)) {
     auto item_ = make_safe(item);
-    
+
     PyBlitzArrayObject* bz = 0;
 
     if (!PyBlitzArray_Converter(item, &bz)) {
diff --git a/xbob/learn/linear/machine.cpp b/xbob/learn/linear/machine.cpp
index fbd85d0d48fdfd80208083f74192ef91ac875c74..79c7cd12caa366279f9fb3ecba344460ed4383b4 100644
--- a/xbob/learn/linear/machine.cpp
+++ b/xbob/learn/linear/machine.cpp
@@ -67,7 +67,7 @@ static int PyBobLearnLinearMachine_init_sizes
         &input_size, &output_size)) return -1;
 
   try {
-    self->cxx = new bob::machine::LinearMachine(input_size, output_size);
+    self->cxx = new bob::learn::linear::Machine(input_size, output_size);
   }
   catch (std::exception& ex) {
     PyErr_SetString(PyExc_RuntimeError, ex.what());
@@ -101,7 +101,7 @@ static int PyBobLearnLinearMachine_init_weights(PyBobLearnLinearMachineObject* s
   }
 
   try {
-    self->cxx = new bob::machine::LinearMachine
+    self->cxx = new bob::learn::linear::Machine
       (*PyBlitzArrayCxx_AsBlitz<double,2>(weights));
   }
   catch (std::exception& ex) {
@@ -132,7 +132,7 @@ static int PyBobLearnLinearMachine_init_hdf5(PyBobLearnLinearMachineObject* self
   auto h5f = reinterpret_cast<PyBobIoHDF5FileObject*>(config);
 
   try {
-    self->cxx = new bob::machine::LinearMachine(*(h5f->f));
+    self->cxx = new bob::learn::linear::Machine(*(h5f->f));
   }
   catch (std::exception& ex) {
     PyErr_SetString(PyExc_RuntimeError, ex.what());
@@ -162,7 +162,7 @@ static int PyBobLearnLinearMachine_init_copy
   auto copy = reinterpret_cast<PyBobLearnLinearMachineObject*>(other);
 
   try {
-    self->cxx = new bob::machine::LinearMachine(*(copy->cxx));
+    self->cxx = new bob::learn::linear::Machine(*(copy->cxx));
   }
   catch (std::exception& ex) {
     PyErr_SetString(PyExc_RuntimeError, ex.what());
@@ -1012,7 +1012,7 @@ PyObject* PyBobLearnLinearMachine_NewFromSize
 
   PyBobLearnLinearMachineObject* retval = (PyBobLearnLinearMachineObject*)PyBobLearnLinearMachine_new(&PyBobLearnLinearMachine_Type, 0, 0);
 
-  retval->cxx = new bob::machine::LinearMachine(input, output);
+  retval->cxx = new bob::learn::linear::Machine(input, output);
 
   return reinterpret_cast<PyObject*>(retval);
 
diff --git a/xbob/learn/linear/main.cpp b/xbob/learn/linear/main.cpp
index e9d864c2257df3a561783f02e88458855c00f0ab..673ccec995f2d5048e2a21ef56e6c2b232134e71 100644
--- a/xbob/learn/linear/main.cpp
+++ b/xbob/learn/linear/main.cpp
@@ -2,7 +2,7 @@
  * @author Andre Anjos <andre.anjos@idiap.ch>
  * @date Fri 13 Dec 2013 12:35:59 CET
  *
- * @brief Bindings to bob::machine
+ * @brief Bindings to bob::learn::linear
  */
 
 #define XBOB_LEARN_LINEAR_MODULE
@@ -20,7 +20,7 @@ static PyMethodDef module_methods[] = {
     {0}  /* Sentinel */
 };
 
-PyDoc_STRVAR(module_docstr, "bob::machine's linear machine and trainers");
+PyDoc_STRVAR(module_docstr, "Bob's Linear machine and trainers");
 
 int PyXbobLearnLinear_APIVersion = XBOB_LEARN_LINEAR_API_VERSION;
 
diff --git a/xbob/learn/linear/pca.cpp b/xbob/learn/linear/pca.cpp
index 0edff39d5b0a19cc49caf3a695d7efba6be505a2..eaeafb3511e791385d7285d8118dec84e1005c56 100644
--- a/xbob/learn/linear/pca.cpp
+++ b/xbob/learn/linear/pca.cpp
@@ -11,7 +11,6 @@
 #include <xbob.blitz/cppapi.h>
 #include <xbob.blitz/cleanup.h>
 #include <bob/config.h>
-#include <bob/trainer/PCATrainer.h>
 #include <xbob.learn.linear/api.h>
 #include <structmember.h>
 
@@ -157,7 +156,7 @@ static int PyBobLearnLinearPCATrainer_init_bool
   if (use_svd_ == -1) return -1;
 
   try {
-    self->cxx = new bob::trainer::PCATrainer(use_svd_?true:false);
+    self->cxx = new bob::learn::linear::PCATrainer(use_svd_?true:false);
   }
   catch (std::exception& ex) {
     PyErr_SetString(PyExc_RuntimeError, ex.what());
@@ -187,7 +186,7 @@ static int PyBobLearnLinearPCATrainer_init_copy
   auto copy = reinterpret_cast<PyBobLearnLinearPCATrainerObject*>(other);
 
   try {
-    self->cxx = new bob::trainer::PCATrainer(*(copy->cxx));
+    self->cxx = new bob::learn::linear::PCATrainer(*(copy->cxx));
   }
   catch (std::exception& ex) {
     PyErr_SetString(PyExc_RuntimeError, ex.what());