diff --git a/bob/fusion/base/algorithm/Algorithm.py b/bob/fusion/base/algorithm/Algorithm.py
index d83cb9d08ef5b188f1b2a0dddf92a39585c23dbf..f08ab93265b30545ecf0c7742c6e05d2b43d88e4 100644
--- a/bob/fusion/base/algorithm/Algorithm.py
+++ b/bob/fusion/base/algorithm/Algorithm.py
@@ -38,60 +38,136 @@ class Algorithm(object):
     self.preprocessors = preprocessors
     self._kwargs = kwargs
     self._kwargs['preprocessors'] = preprocessors
+    if classifier is not self:
+      self._kwargs['classifier'] = classifier
 
   def train_preprocessors(self, X):
+    """Train preprocessors in order.
+    X: numpy.ndarray with the shape of (n_samples, n_systems)."""
     if self.preprocessors is not None:
       for preprocessor in self.preprocessors:
         X = preprocessor.fit_transform(X)
 
   def preprocess(self, scores):
+    """
+    scores: numpy.ndarray with the shape of (n_samples, n_systems).
+    returns the transformed scores."""
     if self.preprocessors is not None:
       for preprocessor in self.preprocessors:
         scores = preprocessor.transform(scores)
     return scores
 
-  def train(self, train, devel=None):
+  def train(self, train_neg, train_pos, devel_neg=None, devel_pos=None):
     """If you use development data for training you need to override this
     method.
-    train: A :py:meth:`tuple` of length 2 containing
-      the negatives and positives. negatives and positives should be
-      numpy.ndarray with the shape of (n_samples, n_systems).
-    devel: same as train but used for development (validation)
+
+    train_neg: numpy.ndarray
+      Negatives training data should be numpy.ndarray with the shape of
+      (n_samples, n_systems).
+    train_pos: numpy.ndarray
+      Positives training data should be numpy.ndarray with the shape of
+      (n_samples, n_systems).
+    devel_neg, devel_pos: numpy.ndarray
+      Same as ``train`` but used for development (validation).
     """
-    (negatives, positives) = train
-    train_scores = np.vstack((negatives, positives))
-    neg_len = negatives.shape[0]
+    train_scores = np.vstack((train_neg, train_pos))
+    neg_len = train_neg.shape[0]
     y = np.zeros((train_scores.shape[0],), dtype='bool')
     y[neg_len:] = True
     self.classifier.fit(train_scores, y)
 
   def fuse(self, scores):
     """
-    scores: A numpy.ndarray with the shape of (n_samples, n_systems).
+    scores: numpy.ndarray
+      A numpy.ndarray with the shape of (n_samples, n_systems).
+
+    **Returns:**
+
+    fused_score: numpy.ndarray
+      The fused scores in shape of (n_samples,).
     """
     return self.classifier.decision_function(scores)
 
   def __str__(self):
-    """__str__() -> info
+    """Return all parameters of this class (and its derived class) in string.
 
-    This function returns all parameters of this class (and its derived class).
 
     **Returns:**
 
-    info : str
+    info: str
       A string containing the full information of all parameters of this
         (and the derived) class.
     """
     return "%s(%s)" % (str(self.__class__), ", ".join(
-      ["%s=%s" % (key, value) for key, value in
-       self._kwargs.items() if value is not None]))
+        ["%s=%s" % (key, value) for key, value in
+         self._kwargs.items() if value is not None]))
 
   def save(self, model_file):
-    """If your class cannot be pickled, you need to override this method."""
-    with open(model_file, "wb") as f:
-      pickle.dump(self, f)
+    """Save the instance of the algorithm.
+
+    model_file: str
+      A path to save the file. Please note that file objects
+      are not accepted. The filename MUST end with ".pkl".
+      Also, an algorithm may save itself in multiple files with different
+      extensions such as model_file and model_file[:-3]+'hdf5'.
+    """
+    # support for bob machines
+    if hasattr(self, "machine"):
+      self.save_bob(model_file)
+    else:
+      with open(model_file, "wb") as f:
+        pickle.dump(self, f)
 
   def load(self, model_file):
-    """If your class cannot be pickled, you need to override this method."""
+    """Load the algorithm the same way it was saved.
+    A new instance will be returned.
+
+    **Returns:**
+
+    loaded_algorithm: Algorithm
+      A new instance of the loaded algorithm.
+    """
+    with open(model_file, "rb") as f:
+      temp = pickle.load(f)
+    if isinstance(temp, Algorithm):
+      return temp
+    else:
+      return self.load_bob(model_file)
+
+  def _get_hdf5_file(self, model_file):
+    return model_file[:-3] + 'hdf5'
+
+  def save_bob(self, model_file):
+    # dump preprocessors in a pickle file because
+    # we don't know how they look like
+    # saves the class to create it later.
+    with open(model_file, 'wb') as f:
+      pickle.dump(self.preprocessors, f)
+      pickle.dump(type(self), f)
+      # just for consistent string representation
+      pickle.dump(self._kwargs, f)
+
+    d5 = bob.io.base.HDF5File(self._get_hdf5_file(model_file), "w")
+    try:
+      self.machine.save(d5)
+    finally:
+      d5.close()
+
+  def load_bob(self, model_file):
+    # load preprocessors and the class
     with open(model_file, "rb") as f:
-      return pickle.load(f)
+      preprocessors = pickle.load(f)
+      myclass = pickle.load(f)
+      _kwargs = pickle.load(f)
+
+    myinstance = myclass(preprocessors=preprocessors)
+    # just for consistent string representation
+    myinstance._kwargs.update(_kwargs)
+
+    d5 = bob.io.base.HDF5File(self._get_hdf5_file(model_file))
+    try:
+      myinstance.machine.load(d5)
+    finally:
+      d5.close()
+
+    return myinstance
diff --git a/bob/fusion/base/algorithm/LLR.py b/bob/fusion/base/algorithm/LLR.py
new file mode 100644
index 0000000000000000000000000000000000000000..e6373449d52fa1e529b344246f2244bf0b4dd943
--- /dev/null
+++ b/bob/fusion/base/algorithm/LLR.py
@@ -0,0 +1,37 @@
+#!/usr/bin/env python
+
+from __future__ import division
+from __future__ import absolute_import
+
+import bob.learn.linear
+
+from .Algorithm import Algorithm
+
+import logging
+logger = logging.getLogger("bob.fusion.base")
+
+
+class LLR(Algorithm):
+  """LLR Score fusion using Bob"""
+
+  def __init__(self,
+               trainer=None,
+               *args, **kwargs):
+    self.trainer = trainer if trainer else \
+        bob.learn.linear.CGLogRegTrainer()
+    # this is needed to be able to load the machine
+    self.machine = bob.learn.linear.Machine()
+    super(LLR, self).__init__(
+        classifier=self,
+        trainer=str(type(self.trainer)),
+        *args, **kwargs)
+
+  def train(self, train_neg, train_pos, devel_neg=None, devel_pos=None):
+    # Trainning the LLR machine
+    self.machine = self.trainer.train(train_neg, train_pos)
+
+  def decision_function(self, scores):
+    scores = self.machine(scores)
+    if scores.ndim == 2 and scores.shape[1] == 1:
+      scores = scores.ravel()
+    return scores
diff --git a/bob/fusion/base/algorithm/MLP.py b/bob/fusion/base/algorithm/MLP.py
index 8de2c08cfba9233a052555439dfaeed935e44662..f707761d84b5bd9553a11f5a1d8a80d31a41076f 100644
--- a/bob/fusion/base/algorithm/MLP.py
+++ b/bob/fusion/base/algorithm/MLP.py
@@ -50,13 +50,13 @@ class MLP(Algorithm):
     self.trainer = self.trainer if self.trainer and not force else \
         bob.learn.mlp.RProp(1, bob.learn.mlp.SquareError(
             self.machine.output_activation), machine=self.machine,
-          train_biases=False)
-    self._kwargs = {
-      'seed': self.seed,
-      'mlp_shape': self.mlp_shape,
-      'machine': self.machine,
-      'train': self.train,
-    }
+            train_biases=False)
+    self._kwargs.update({
+        'seed': self.seed,
+        'mlp_shape': self.mlp_shape,
+        'machine': str(self.machine),
+        'trainer': str(type(self.trainer))})
+    self._kwargs.update(self._my_kwargs)
 
   def prepare_train(self, train, devel):
     (negatives, positives) = train
@@ -64,7 +64,8 @@ class MLP(Algorithm):
     if n_systems != self.mlp_shape[0]:
       logger.warn(
         'Reinitializing the MLP machine with the shape of {} to {} to match th'
-        'e input size.'.format(self.mlp_shape, [n_systems]+self.mlp_shape[1:]))
+        'e input size.'.format(self.mlp_shape,
+                               [n_systems] + self.mlp_shape[1:]))
       self.mlp_shape = [n_systems] + self.mlp_shape[1:]
       self.n_systems = n_systems
       self.hidden_layers = self.mlp_shape[1:-1]
@@ -77,10 +78,12 @@ class MLP(Algorithm):
         trainer=self.trainer,
         **self._my_kwargs)
 
-  def train(self, train, devel=None):
-    if devel is None:
-      devel = train
-    self.prepare_train(train, devel)
+  def train(self, train_neg, train_pos, devel_neg=None, devel_pos=None):
+    if devel_neg is None:
+      devel_neg = train_neg
+    if devel_pos is None:
+      devel_pos = train_pos
+    self.prepare_train((train_neg, train_pos), (devel_neg, devel_pos))
     self.machine, self.analyzer = self.train_helper()
 
   def decision_function(self, scores):
@@ -88,31 +91,3 @@ class MLP(Algorithm):
     if scores.ndim == 2 and scores.shape[1] == 1:
       scores = scores.ravel()
     return scores
-
-  def _get_hdf5_file(self, model_file):
-    return model_file[:-3] + 'hdf5'
-
-  def save(self, model_file):
-    d5 = bob.io.base.HDF5File(self._get_hdf5_file(model_file), "w")
-    try:
-      self.machine.save(d5)
-    finally:
-      d5.close()
-
-    # dump preprocessors in a pickle file because
-    # we don't know how they look like
-    with open(model_file, 'wb') as f:
-      pickle.dump(self.preprocessors, f)
-
-  def load(self, model_file):
-    d5 = bob.io.base.HDF5File(self._get_hdf5_file(model_file))
-    try:
-      self.machine.load(d5)
-    finally:
-      d5.close()
-
-    # load preprocessors
-    with open(model_file, "rb") as f:
-      self.preprocessors = pickle.load(f)
-
-    return self
diff --git a/bob/fusion/base/algorithm/Weighted_Sum.py b/bob/fusion/base/algorithm/Weighted_Sum.py
index 0a70e0e4daf422641d302ff71c6f2f6362fa6d62..7a5fcc70386f9130b45c1b58c725963dc96bc486 100644
--- a/bob/fusion/base/algorithm/Weighted_Sum.py
+++ b/bob/fusion/base/algorithm/Weighted_Sum.py
@@ -29,12 +29,3 @@ class Weighted_Sum(Algorithm):
       return numpy.mean(scores, axis=1)
     else:
       return numpy.sum(scores * self.weights, axis=1)
-
-  def closed_form(self, x1, y):
-    if self.weights is None:
-      return 2*y - x1
-    else:
-      w1 = self.weights[0]
-      w2 = self.weights[1]
-      x2 = (y - x1*w1)/w2
-      return x2
diff --git a/bob/fusion/base/algorithm/__init__.py b/bob/fusion/base/algorithm/__init__.py
index 39cd044f1ebed29a61df2ea8a757f92351221b26..25f3fc12a4df60c9fcd68cdf3622472347ef47e4 100644
--- a/bob/fusion/base/algorithm/__init__.py
+++ b/bob/fusion/base/algorithm/__init__.py
@@ -1,6 +1,7 @@
 from .Algorithm import Algorithm
 from .Weighted_Sum import Weighted_Sum
 from .MLP import MLP
+from .LLR import LLR
 
 # gets sphinx autodoc done right - don't remove it
 __all__ = [_ for _ in dir() if not _.startswith('_')]
diff --git a/bob/fusion/base/algorithm/mlp_train_helper.py b/bob/fusion/base/algorithm/mlp_train_helper.py
index 8ae2f8b43628e902273c435f6c01010cac4219dd..6a3fb8829a54a25513d3c4676917936f4fcce0e5 100644
--- a/bob/fusion/base/algorithm/mlp_train_helper.py
+++ b/bob/fusion/base/algorithm/mlp_train_helper.py
@@ -192,13 +192,8 @@ class MLPTrainer(object):
     self.max_iter = max_iter
     self.no_improvements = no_improvements
     self.valley_condition = valley_condition
-    self.machine = machine if machine else \
-        bob.learn.mlp.Machine(self.mlp_shape)
-    self.machine.randomize()
-    self.trainer = trainer if trainer else \
-        bob.learn.mlp.RProp(batch_size, bob.learn.mlp.SquareError(
-          self.machine.output_activation), machine=self.machine,
-          train_biases=False)
+    self.machine = machine
+    self.trainer = trainer
 
   def __call__(self):
     return self.make_mlp()
@@ -232,8 +227,6 @@ class MLPTrainer(object):
     # machine = bob.learn.mlp.Machine(self.shape)
     # machine.activation = bob.learn.activation.HyperbolicTangent() #the
     # defaults are anyway Hyperbolic Tangent for hidden and output layer
-    # machine.randomize()
-    # import ipdb; ipdb.set_trace()
     self.machine.input_subtract, self.machine.input_divide = \
         shuffler.stdnorm()
 
@@ -277,7 +270,8 @@ class MLPTrainer(object):
           min_devel_rmse = avg_devel_rmse
           logger.info("%d: New valley stop threshold set to %.4e",
                       iteration, avg_devel_rmse / VALLEY_CONDITION)
-        if stop_condition(min_devel_rmse, avg_devel_rmse, last_devel_rmse):
+        if stop_condition(min_devel_rmse, avg_devel_rmse, last_devel_rmse) \
+           or numpy.allclose(avg_devel_rmse / VALLEY_CONDITION, 0):
           logger.info("%d: Stopping on devel valley condition", iteration)
           logger.info("%d: Best machine happened on iteration %d with average "
                       "devel. RMSE of %.4e", iteration, best_machine_iteration,
@@ -288,7 +282,6 @@ class MLPTrainer(object):
 
         # train for 'epoch' times w/o stopping for tests
         for i in range(self.epoch):
-          # import ipdb; ipdb.set_trace()
           shuffler(data=shuffled_input, target=shuffled_target)
           self.trainer.batch_size = len(shuffled_input)
           self.trainer.train(
diff --git a/bob/fusion/base/config/algorithm/llr.py b/bob/fusion/base/config/algorithm/llr.py
index 9a4c61726809a71050926bc66ed9693e557fb8c6..5b1676598b800258c434b75d56de73c145f81ffa 100644
--- a/bob/fusion/base/config/algorithm/llr.py
+++ b/bob/fusion/base/config/algorithm/llr.py
@@ -2,8 +2,5 @@
 
 import bob.fusion.base
 from sklearn.preprocessing import StandardScaler
-from sklearn.linear_model import LogisticRegression
 
-algorithm = bob.fusion.base.algorithm.Algorithm(
-  preprocessors=[StandardScaler()],
-  classifier=LogisticRegression())
+algorithm = bob.fusion.base.algorithm.LLR(preprocessors=[StandardScaler()])
diff --git a/bob/fusion/base/config/algorithm/llr_skl.py b/bob/fusion/base/config/algorithm/llr_skl.py
new file mode 100644
index 0000000000000000000000000000000000000000..9a4c61726809a71050926bc66ed9693e557fb8c6
--- /dev/null
+++ b/bob/fusion/base/config/algorithm/llr_skl.py
@@ -0,0 +1,9 @@
+#!/usr/bin/env python
+
+import bob.fusion.base
+from sklearn.preprocessing import StandardScaler
+from sklearn.linear_model import LogisticRegression
+
+algorithm = bob.fusion.base.algorithm.Algorithm(
+  preprocessors=[StandardScaler()],
+  classifier=LogisticRegression())
diff --git a/bob/fusion/base/script/__init__.py b/bob/fusion/base/script/__init__.py
index b22516fe54dc4528352cf701281c2bc4790b0a3c..8c298d25dcd43aa5a6f83a775ddf852ea6a713df 100644
--- a/bob/fusion/base/script/__init__.py
+++ b/bob/fusion/base/script/__init__.py
@@ -1,2 +1,2 @@
-from . import fuse
+from . import bob_fuse
 from . import plot_fusion_decision_boundary
diff --git a/bob/fusion/base/script/fuse.py b/bob/fusion/base/script/bob_fuse.py
similarity index 95%
rename from bob/fusion/base/script/fuse.py
rename to bob/fusion/base/script/bob_fuse.py
index 556d6f2ffcb25ff0361dea04f6680c747dc73eac..f14dc75eb423965d264acb37f47c6d9d7f34af94 100755
--- a/bob/fusion/base/script/fuse.py
+++ b/bob/fusion/base/script/bob_fuse.py
@@ -28,7 +28,7 @@ def fuse(args, command_line_parameters):
   score_lines_list_dev = [load_score(path, ncolumns=args.score_type)
                           for path in args.dev_files]
   scores_dev = get_all_scores(score_lines_list_dev)
-  trainer_scores = get_negatives_positives_all(score_lines_list_dev)
+  (neg, pos) = get_negatives_positives_all(score_lines_list_dev)
   if args.eval_files:
     score_lines_list_eval = [load_score(path, ncolumns=args.score_type)
                              for path in args.eval_files]
@@ -55,9 +55,7 @@ def fuse(args, command_line_parameters):
   # preprocess data
   scores_dev = algorithm.preprocess(scores_dev)
   scores_eval = algorithm.preprocess(scores_eval)
-  neg, pos = trainer_scores
   neg, pos = algorithm.preprocess(neg), algorithm.preprocess(pos)
-  trainer_scores = (neg, pos)
 
   # train the model
   if utils.check_file(args.model_file, args.force, 1000):
@@ -65,7 +63,7 @@ def fuse(args, command_line_parameters):
       "- Fusion: model '%s' already exists.", args.model_file)
     algorithm = algorithm.load(args.model_file)
   else:
-    algorithm.train(trainer_scores)
+    algorithm.train(neg, pos)
     algorithm.save(args.model_file)
 
   # fuse the scores (dev)
diff --git a/bob/fusion/base/script/plot_fusion_decision_boundary.py b/bob/fusion/base/script/plot_fusion_decision_boundary.py
index f5817fcea1b57dd53999b0d803221067ef739482..9452d91510db54e489346d692800907e511e4a53 100644
--- a/bob/fusion/base/script/plot_fusion_decision_boundary.py
+++ b/bob/fusion/base/script/plot_fusion_decision_boundary.py
@@ -34,7 +34,7 @@ import numpy
 import bob.fusion.base
 import bob.core
 from bob.measure.load import load_score, get_negatives_positives,\
-  get_all_scores
+    get_all_scores
 from bob.measure import eer_threshold
 from ..tools import grouping
 
@@ -68,8 +68,8 @@ def plot_boundary_decision(algorithm, scores, score_labels, threshold,
 
   if scores.shape[1] > 2:
     raise NotImplementedError(
-      "Currently plotting the decision boundary for more than two systems "
-      "is not supported.")
+        "Currently plotting the decision boundary for more than two systems "
+        "is not supported.")
 
   import matplotlib.pyplot as plt
   plt.gca()  # this is necessary for subplots to work.
@@ -79,8 +79,8 @@ def plot_boundary_decision(algorithm, scores, score_labels, threshold,
   x_min, x_max = X[:, i1].min() - x_pad, X[:, i1].max() + x_pad
   y_min, y_max = X[:, i2].min() - y_pad, X[:, i2].max() + y_pad
   xx, yy = numpy.meshgrid(
-    numpy.linspace(x_min, x_max, resolution),
-    numpy.linspace(y_min, y_max, resolution))
+      numpy.linspace(x_min, x_max, resolution),
+      numpy.linspace(y_min, y_max, resolution))
   temp = numpy.c_[xx.ravel(), yy.ravel()]
   temp = algorithm.preprocess(temp)
   Z = (algorithm.fuse(temp) > threshold).reshape(xx.shape)
@@ -92,16 +92,16 @@ def plot_boundary_decision(algorithm, scores, score_labels, threshold,
     negatives, positives = grouping(negatives, positives, **kwargs)
     X = numpy.concatenate((negatives, positives), axis=0)
     Y = numpy.concatenate(
-      (numpy.zeros(negatives.shape[0], dtype=numpy.bool8),
-       numpy.ones(positives.shape[0], dtype=numpy.bool8)),
-      axis=0)
+        (numpy.zeros(negatives.shape[0], dtype=numpy.bool8),
+         numpy.ones(positives.shape[0], dtype=numpy.bool8)),
+        axis=0)
 
   negatives, positives = X[numpy.logical_not(Y)], X[Y]
   colors = plt.cm.viridis(numpy.linspace(0, 1, 2))
   for i, X in enumerate((negatives, positives)):
     plt.scatter(
-      X[:, 0], X[:, 1], marker=markers[i], alpha=alpha,
-      c=colors[i], label=legends[i])
+        X[:, 0], X[:, 1], marker=markers[i], alpha=alpha,
+        c=colors[i], label=legends[i])
   plt.legend()
 
   if thres_system1 is not None:
@@ -126,7 +126,8 @@ def main(command_line_parameters=None):
 
   # load the scores
   score_lines_list = [
-    load_score(path, int(args['--score-type'])) for path in args['SCORE_FILE']]
+      load_score(path, int(args['--score-type'])) for
+      path in args['SCORE_FILE']]
   scores = get_all_scores(score_lines_list)
   score_lines = numpy.array(score_lines_list[0])
   score_lines['score'] = algorithm.fuse(algorithm.preprocess(scores))
@@ -135,11 +136,11 @@ def main(command_line_parameters=None):
 
   # plot the decision boundary
   plot_boundary_decision(
-    algorithm, scores, score_labels, threshold,
-    do_grouping=True,
-    npoints=int(args['--group']),
-    seed=0,
-    gformat=args['--grouping']
+      algorithm, scores, score_labels, threshold,
+      do_grouping=True,
+      npoints=int(args['--group']),
+      seed=0,
+      gformat=args['--grouping']
   )
   plt.savefig(args['--output'])
   plt.close()
diff --git a/bob/fusion/base/test/test_algorithm.py b/bob/fusion/base/test/test_algorithm.py
new file mode 100644
index 0000000000000000000000000000000000000000..359e7336c957faa9e1763ed4c564c51daec08c56
--- /dev/null
+++ b/bob/fusion/base/test/test_algorithm.py
@@ -0,0 +1,188 @@
+#!/usr/bin/env python
+import warnings
+import numpy
+from numpy import array
+from tempfile import NamedTemporaryFile
+import bob.fusion.base
+import bob.learn.linear
+from sklearn.preprocessing import StandardScaler
+from sklearn.linear_model import LogisticRegression
+
+
+NEG = array([[-1.23594765, -2.59984279],
+             [-2.02126202, -0.7591068],
+             [-1.13244201, -3.97727788],
+             [-2.04991158, -3.15135721],
+             [-3.10321885, -2.5894015],
+             [-2.85595643, -1.54572649],
+             [-2.23896227, -2.87832498],
+             [-2.55613677, -2.66632567],
+             [-1.50592093, -3.20515826],
+             [-2.6869323, -3.85409574]])
+POS = array([[0.44701018, 3.6536186],
+             [3.8644362, 2.25783498],
+             [5.26975462, 1.54563433],
+             [3.04575852, 2.81281615],
+             [4.53277921, 4.46935877],
+             [3.15494743, 3.37816252],
+             [2.11221425, 1.01920353],
+             [2.65208785, 3.15634897],
+             [4.23029068, 4.20237985],
+             [2.61267318, 2.69769725]])
+X = numpy.vstack((NEG, POS))
+TEST = array([[-1.04855297, -1.42001794],
+              [-1.70627019, 1.9507754],
+              [-0.50965218, -0.4380743],
+              [-1.25279536, 0.77749036],
+              [-1.61389785, -0.21274028],
+              [-0.89546656, 0.3869025],
+              [-0.51080514, -1.18063218],
+              [-0.02818223, 0.42833187],
+              [0.06651722, 0.3024719],
+              [-0.63432209, -0.36274117]])
+
+
+def run_steps(algorithm):
+  algorithm.train_preprocessors(X)
+  neg = algorithm.preprocess(NEG)
+  pos = algorithm.preprocess(POS)
+  algorithm.train(neg, pos)
+  fused = algorithm.fuse(TEST)
+  with NamedTemporaryFile(suffix='.pkl') as f:
+    algorithm.save(f.name)
+    loaded_algorithm = algorithm.load(f.name)
+
+  try:
+    assert str(algorithm) == str(loaded_algorithm)
+  except Exception:
+    warnings.warn("String comparison of algorithms do not match which is OK.")
+    print(str(algorithm))
+    print(str(loaded_algorithm))
+  if algorithm.preprocessors:
+    assert len(algorithm.preprocessors) == len(loaded_algorithm.preprocessors)
+  assert fused.ndim == 1
+
+  return neg, pos, fused, loaded_algorithm
+
+
+def test_algorithm_llr_sklearn():
+  algorithm = bob.fusion.base.algorithm.Algorithm(
+      preprocessors=[
+        StandardScaler(**{'copy': True, 'with_mean': True, 'with_std': True})],
+      classifier=LogisticRegression(**{'C': 1.0,
+                                       'class_weight': None,
+                                       'dual': False,
+                                       'fit_intercept': True,
+                                       'intercept_scaling': 1,
+                                       'max_iter': 100,
+                                       'multi_class': 'ovr',
+                                       'n_jobs': 1,
+                                       'penalty': 'l2',
+                                       'random_state': None,
+                                       'solver': 'liblinear',
+                                       'tol': 0.0001,
+                                       'verbose': 0,
+                                       'warm_start': False}))
+  neg, pos, fused, loaded_algorithm = run_steps(algorithm)
+  assert str(algorithm) == "<class 'bob.fusion.base.algorithm.Algorithm.Algorithm'>(preprocessors=[StandardScaler(copy=True, with_mean=True, with_std=True)], classifier=LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,\n          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,\n          penalty='l2', random_state=None, solver='liblinear', tol=0.0001,\n          verbose=0, warm_start=False))"
+
+  assert numpy.allclose(algorithm.preprocessors[0].mean_, array([0.52676307, 0.09832188]))
+  assert numpy.allclose(algorithm.preprocessors[0].scale_, array([2.857145, 2.98815147]))
+
+  assert numpy.allclose(neg, array([[-0.61694829, -0.90295445],
+                                    [-0.89180811, -0.28694284],
+                                    [-0.58072134, -1.36392007],
+                                    [-0.90183545, -1.08752154],
+                                    [-1.27049272, -0.89946022],
+                                    [-1.18395093, -0.5501891],
+                                    [-0.96800314, -0.99614993],
+                                    [-1.07901413, -0.92520328],
+                                    [-0.71143886, -1.10552634],
+                                    [-1.12479253, -1.32269654]]))
+  assert numpy.allclose(pos, array([[-0.02791349, 1.18979803],
+                                    [1.16818472, 0.72269198],
+                                    [1.6600458, 0.48435043],
+                                    [0.88164775, 0.90841923],
+                                    [1.4021046, 1.4627896],
+                                    [0.91986384, 1.09761526],
+                                    [0.5549075, 0.3081777],
+                                    [0.74386312, 1.02338423],
+                                    [1.29623369, 1.37344375],
+                                    [0.73006799, 0.86989411]]))
+
+  assert numpy.allclose(algorithm.classifier.intercept_, array([0.04577333]))
+  assert numpy.allclose(algorithm.classifier.classes_, array([False, True], dtype=bool))
+  assert numpy.allclose(algorithm.classifier.coef_, array([[1.33489128, 1.38092354]]))
+
+  assert numpy.allclose(fused, array([-3.31486708, 0.4619598, -1.23950404, -0.55291754, -2.40238289,
+                                      -0.61529441, -2.26645877, 0.59964668, 0.55225715, -1.30189552]))
+
+  assert numpy.allclose(algorithm.preprocessors[0].mean_, loaded_algorithm.preprocessors[0].mean_)
+  assert numpy.allclose(algorithm.preprocessors[0].scale_, loaded_algorithm.preprocessors[0].scale_)
+  assert numpy.allclose(algorithm.classifier.intercept_, loaded_algorithm.classifier.intercept_)
+  assert numpy.allclose(algorithm.classifier.classes_, loaded_algorithm.classifier.classes_)
+  assert numpy.allclose(algorithm.classifier.coef_, loaded_algorithm.classifier.coef_)
+
+
+def test_algorithm_llr_bob():
+  algorithm = bob.fusion.base.algorithm.LLR(
+    preprocessors=[StandardScaler(**{'copy': True,
+                                     'with_mean': True,
+                                     'with_std': True})],
+    trainer=None)
+  neg, pos, fused, loaded_algorithm = run_steps(algorithm)
+
+  algorithm = bob.fusion.base.algorithm.LLR(
+    preprocessors=[StandardScaler(**{'copy': True,
+                                     'with_mean': True,
+                                     'with_std': True})],
+    trainer=bob.learn.linear.CGLogRegTrainer(
+      prior=0.5, convergence_threshold=1e-5,
+      max_iterations=10000, reg=1., mean_std_norm=False))
+  neg, pos, fused, loaded_algorithm = run_steps(algorithm)
+  assert str(algorithm) == "<class 'bob.fusion.base.algorithm.LLR.LLR'>(trainer=<type 'bob.learn.linear.CGLogRegTrainer'>, preprocessors=[StandardScaler(copy=True, with_mean=True, with_std=True)])"
+
+  assert numpy.allclose(algorithm.machine.biases, array([0.04577333]), atol=0.05)
+  assert numpy.allclose(algorithm.machine.weights, array([[1.33489128, 1.38092354]]), atol=0.05)
+
+  assert numpy.allclose(fused, array([-3.31486708, 0.4619598, -1.23950404, -0.55291754, -2.40238289,
+                                      -0.61529441, -2.26645877, 0.59964668, 0.55225715, -1.30189552]), atol=0.05)
+
+  assert numpy.allclose(algorithm.machine.biases, loaded_algorithm.machine.biases)
+  assert numpy.allclose(algorithm.machine.weights, loaded_algorithm.machine.weights)
+
+
+def test_weighted_sum_1():
+  algorithm = bob.fusion.base.algorithm.Weighted_Sum()
+  neg, pos, fused, loaded_algorithm = run_steps(algorithm)
+  assert str(algorithm) == "<class 'bob.fusion.base.algorithm.Weighted_Sum.Weighted_Sum'>()"
+  assert numpy.allclose(fused, numpy.mean(TEST, axis=1))
+  assert algorithm.weights == loaded_algorithm.weights
+
+
+def test_weighted_sum_2():
+  weights = [0.3, 0.7]
+  algorithm = bob.fusion.base.algorithm.Weighted_Sum(weights=weights)
+  neg, pos, fused, loaded_algorithm = run_steps(algorithm)
+  assert str(
+    algorithm) == "<class 'bob.fusion.base.algorithm.Weighted_Sum.Weighted_Sum'>(weights=[0.3, 0.7])"
+  assert numpy.allclose(fused, numpy.sum(TEST * weights, axis=1))
+  assert algorithm.weights == loaded_algorithm.weights
+
+
+def test_mlp():
+  algorithm = bob.fusion.base.algorithm.MLP(
+      n_systems=2, hidden_layers=[3], seed=0, batch_size=10, epoch=1,
+      max_iter=1000, no_improvements=0, valley_condition=0.9,
+      preprocessors=[StandardScaler(**{'copy': True,
+                                       'with_mean': True,
+                                       'with_std': True})])
+  assert numpy.allclose(algorithm.machine.weights[0], array([[0.0097627, 0.01856892, 0.04303787],
+                                                             [0.06885315, 0.02055267, 0.07158912]]))
+  assert numpy.allclose(algorithm.machine.weights[1], array([[0.02471274],
+                                                             [0.02917882],
+                                                             [-0.02312366]]))
+  assert numpy.allclose(algorithm.machine.biases[0], array([0.00897664, 0.06945035, -0.01526904]))
+  assert numpy.allclose(algorithm.machine.biases[1], array([-0.01248256]))
+  _, _, fused, loaded_algorithm = run_steps(algorithm)
+  assert numpy.allclose(fused, [-1, 1, -1, -1, -1, -1, -1, 1, 1, -1], atol=0.001)
diff --git a/bob/fusion/base/test/test_scripts.py b/bob/fusion/base/test/test_scripts.py
new file mode 100644
index 0000000000000000000000000000000000000000..4265cc3e6c16c09774190fa55d609cd9fe0808e4
--- /dev/null
+++ b/bob/fusion/base/test/test_scripts.py
@@ -0,0 +1 @@
+#!/usr/bin/env python
diff --git a/buildout.cfg b/buildout.cfg
index 4184edc32de25d3f054a905fcb058ac689a015fe..20f0513df5c63c49a249e8403a919f768e578310 100644
--- a/buildout.cfg
+++ b/buildout.cfg
@@ -7,6 +7,7 @@ parts = scripts
 eggs = bob.fusion.base
        gridtk
        ipdb
+       pytest
 
 extensions = bob.buildout
              mr.developer
diff --git a/requirements.txt b/requirements.txt
index eb97d3f8cfdaa2f41157f6b5b65bc5410372653b..410dd7631d2af68b97ca9d8461a6866b1aea03a8 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -3,6 +3,7 @@ numpy
 bob.core
 bob.extension
 bob.measure
+bob.learn.linear
 bob.learn.activation
 bob.learn.mlp
 bob.bio.base
diff --git a/setup.py b/setup.py
index 362aaab9d706087b1cb437375bc8292158a32ccd..77f706d0e927a7d8ed810dd1da0a52bbcb62f0e9 100644
--- a/setup.py
+++ b/setup.py
@@ -101,13 +101,14 @@ setup(
 
       # scripts should be declared using this entry:
       'console_scripts': [
-        'fuse.py     = bob.fusion.base.script.fuse:main',
+        'bob_fuse.py     = bob.fusion.base.script.bob_fuse:main',
         'plot_fusion_decision_boundary.py = bob.fusion.base.script.plot_fusion_decision_boundary:main',
       ],
 
       'bob.fusion.algorithm': [
         'mean        = bob.fusion.base.config.algorithm.mean:algorithm',
         'llr         = bob.fusion.base.config.algorithm.llr:algorithm',
+        'llr-skl     = bob.fusion.base.config.algorithm.llr_skl:algorithm',
         'plr-2       = bob.fusion.base.config.algorithm.plr_2:algorithm',
         'mlp         = bob.fusion.base.config.algorithm.mlp:algorithm',
       ],