diff --git a/src/mlpack/methods/ann/connections/full_connection.hpp b/src/mlpack/methods/ann/connections/full_connection.hpp
index da95755953..dcccb9bb66 100644
--- a/src/mlpack/methods/ann/connections/full_connection.hpp
+++ b/src/mlpack/methods/ann/connections/full_connection.hpp
@@ -83,6 +83,16 @@ class FullConnection
     delta = (weights.t() * error);
   }
 
+  /*
+   * Calculate the gradient using the output delta and the input activation.
+   *
+   * @param gradient The calculated gradient.
+   */
+  void Gradient(MatType& gradient)
+  {
+    gradient = outputLayer.Delta() * inputLayer.InputActivation().t();
+  }
+
   //! Get the weights.
   MatType& Weights() const { return weights; }
   //! Modify the weights.
diff --git a/src/mlpack/methods/ann/connections/fullself_connection.hpp b/src/mlpack/methods/ann/connections/fullself_connection.hpp
index 37b6dc325d..470be4ee11 100644
--- a/src/mlpack/methods/ann/connections/fullself_connection.hpp
+++ b/src/mlpack/methods/ann/connections/fullself_connection.hpp
@@ -85,6 +85,16 @@ class FullselfConnection
     delta = (weights.t() * error);
   }
 
+  /*
+   * Calculate the gradient using the output delta and the input activation.
+   *
+   * @param gradient The calculated gradient.
+   */
+  void Gradient(MatType& gradient)
+  {
+    gradient = outputLayer.Delta() * inputLayer.InputActivation().t();
+  }
+
   //! Get the weights.
   MatType& Weights() const { return weights; }
   //! Modify the weights.
diff --git a/src/mlpack/methods/ann/connections/self_connection.hpp b/src/mlpack/methods/ann/connections/self_connection.hpp
index 6401c42bda..3b67a8a376 100644
--- a/src/mlpack/methods/ann/connections/self_connection.hpp
+++ b/src/mlpack/methods/ann/connections/self_connection.hpp
@@ -18,7 +18,7 @@ namespace ann /** Artificial Neural Network. */ {
 /**
  * Implementation of the self connection class. The self connection connects
  * every neuron from the input layer with the output layer in a multiplicative
- * way.
+ * way, except the elements on the main diagonal.
  *
  * @tparam InputLayerType Type of the connected input layer.
  * @tparam OutputLayerType Type of the connected output layer.
@@ -39,7 +39,7 @@ class SelfConnection
 {
  public:
   /**
-   * Create the FullConnection object using the specified input layer, output
+   * Create the SelfConnection object using the specified input layer, output
    * layer, optimizer and weight initialize rule.
    *
    * @param InputLayerType The input layer which is connected with the output
@@ -54,20 +54,26 @@ class SelfConnection
                  OutputLayerType& outputLayer,
                  OptimizerType& optimizer,
                  WeightInitRule weightInitRule = WeightInitRule()) :
-      inputLayer(inputLayer), outputLayer(outputLayer), optimizer(optimizer)
+      inputLayer(inputLayer),
+      outputLayer(outputLayer),
+      optimizer(optimizer),
+      connection(1 - arma::eye<MatType>(inputLayer.OutputSize(),
+          inputLayer.OutputSize()))
   {
-    weightInitRule.Initialize(weights, outputLayer.OutputSize(), 1);
+    weightInitRule.Initialize(weights, outputLayer.InputSize(),
+        inputLayer.OutputSize());
   }
 
   /**
    * Ordinary feed forward pass of a neural network, evaluating the function
    * f(x) by propagating the activity forward through f.
    *
-   * @param input Input data used for evaluating the specified activity function.
+   * @param input Input data used for evaluating the specified activity
+   * function.
    */
   void FeedForward(const VecType& input)
   {
-    outputLayer.InputActivation() += (weights % input);
+    outputLayer.InputActivation() += (weights % connection) * input;
   }
 
   /**
@@ -79,9 +85,17 @@ class SelfConnection
    */
   void FeedBackward(const VecType& error)
   {
-    // Calculating the delta using the partial derivative of the error with
-    // respect to a weight.
-    delta = (weights.t() * error);
+    delta = (weights % connection).t() * error;
+  }
+
+  /*
+   * Calculate the gradient using the output delta and the input activation.
+   *
+   * @param gradient The calculated gradient.
+   */
+  void Gradient(MatType& gradient)
+  {
+    gradient = outputLayer.Delta() * inputLayer.InputActivation().t();
   }
 
   //! Get the weights.
@@ -124,6 +138,9 @@ class SelfConnection
 
   //! Locally-stored detla object that holds the calculated delta.
   VecType delta;
+
+  //! Locally-stored connection multiplication type.
+  MatType connection;
 }; // class SelfConnection
 
 //! Connection traits for the self connection.
diff --git a/src/mlpack/methods/ann/ffnn.hpp b/src/mlpack/methods/ann/ffnn.hpp
index 4b5e951fe8..8cfb3479b6 100644
--- a/src/mlpack/methods/ann/ffnn.hpp
+++ b/src/mlpack/methods/ann/ffnn.hpp
@@ -44,7 +44,7 @@ class FFNN
      * @param outputLayer The outputlayer used to evaluate the network.
      */
     FFNN(const ConnectionTypes& network, OutputLayerType& outputLayer)
-        : network(network), outputLayer(outputLayer), err(0)
+        : network(network), outputLayer(outputLayer), trainError(0), seqNum(0)
     {
       // Nothing to do here.
     }
@@ -64,13 +64,8 @@ class FFNN
                      const VecType& target,
                      VecType& error)
     {
-      ResetActivations(network);
       seqNum++;
-
-      std::get<0>(std::get<0>(network)).InputLayer().InputActivation() = input;
-
-      FeedForward(network);
-      OutputError(network, target, error);
+      trainError += Evaluate(input, target, error);
     }
 
     /**
@@ -102,7 +97,7 @@ class FFNN
       ApplyGradients(network);
 
       // Reset the overall error.
-      err = 0;
+      trainError = 0;
       seqNum = 0;
     }
 
@@ -125,6 +120,26 @@ class FFNN
       OutputPrediction(network, output);
     }
 
+    /**
+     * Evaluate the trained network using the given input and compare the output
+     * with the given target vector.
+     *
+     * @param input Input data used to evaluate the trained network.
+     * @param target Target data used to calculate the network error.
+     * @param error The calulated error of the output layer.
+     * @tparam VecType Type of data (arma::colvec, arma::mat or arma::sp_mat).
+     */
+    template <typename VecType>
+    double Evaluate(const VecType& input, const VecType& target, VecType& error)
+    {
+      ResetActivations(network);
+
+      std::get<0>(std::get<0>(network)).InputLayer().InputActivation() = input;
+
+      FeedForward(network);
+      return OutputError(network, target, error);
+    }
+
     //! Get the error of the network.
     double Error() const { return trainError; }
 
@@ -218,9 +233,9 @@ class FFNN
      * Calculate the output error and update the overall error.
      */
     template<typename VecType, typename... Tp>
-    void OutputError(std::tuple<Tp...>& t,
-                     const VecType& target,
-                     VecType& error)
+    double OutputError(std::tuple<Tp...>& t,
+                      const VecType& target,
+                      VecType& error)
     {
        // Calculate and store the output error.
       outputLayer.calculateError(std::get<0>(
@@ -229,12 +244,9 @@ class FFNN
 
       // Masures the network's performance with the specified performance
       // function.
-      err += PerformanceFunction::error(std::get<0>(
+      return PerformanceFunction::error(std::get<0>(
           std::get<sizeof...(Tp) - 1>(t)).OutputLayer().InputActivation(),
           target);
-
-      // Update the final training error.
-      trainError = err;
     }
 
     /*
@@ -353,8 +365,9 @@ class FFNN
     typename std::enable_if<I < sizeof...(Tp), void>::type
     Gradients(std::tuple<Tp...>& t)
     {
-      gradients[gradientNum++] += std::get<I>(t).OutputLayer().Delta() *
-          std::get<I>(t).InputLayer().InputActivation().t();
+      MatType gradient;
+      std::get<I>(t).Gradient(gradient);
+      gradients[gradientNum++] += gradient;
 
       Gradients<I + 1, Tp...>(t);
     }
@@ -395,11 +408,12 @@ class FFNN
     typename std::enable_if<I < sizeof...(Tp), void>::type
     Apply(std::tuple<Tp...>& t)
     {
+      // Take a mean gradient step over the number of inputs.
       if (seqNum > 1)
         gradients[gradientNum] /= seqNum;
 
       std::get<I>(t).Optimzer().UpdateWeights(std::get<I>(t).Weights(),
-          gradients[gradientNum], err);
+          gradients[gradientNum], trainError);
 
       // Reset the gradient storage.
       gradients[gradientNum++].zeros();
@@ -456,9 +470,6 @@ class FFNN
     //! The outputlayer used to evaluate the network
     OutputLayerType& outputLayer;
 
-    //! The current error of the network.
-    double err;
-
     //! The current training error of the network.
     double trainError;
 
diff --git a/src/mlpack/methods/ann/layer/binary_classification_layer.hpp b/src/mlpack/methods/ann/layer/binary_classification_layer.hpp
index c3e38f1d6b..c928f7dc61 100644
--- a/src/mlpack/methods/ann/layer/binary_classification_layer.hpp
+++ b/src/mlpack/methods/ann/layer/binary_classification_layer.hpp
@@ -81,5 +81,4 @@ class LayerTraits<BinaryClassificationLayer<MatType, VecType> >
 }; // namespace ann
 }; // namespace mlpack
 
-
 #endif
diff --git a/src/mlpack/methods/ann/layer/multiclass_classification_layer.hpp b/src/mlpack/methods/ann/layer/multiclass_classification_layer.hpp
index aa9a57eb42..40ca62be42 100644
--- a/src/mlpack/methods/ann/layer/multiclass_classification_layer.hpp
+++ b/src/mlpack/methods/ann/layer/multiclass_classification_layer.hpp
@@ -94,5 +94,4 @@ using ClassificationLayer = MulticlassClassificationLayer<MatType, VecType>;
 }; // namespace ann
 }; // namespace mlpack
 
-
 #endif
diff --git a/src/mlpack/methods/ann/optimizer/steepest_descent.hpp b/src/mlpack/methods/ann/optimizer/steepest_descent.hpp
index 581eec315b..8b43897788 100644
--- a/src/mlpack/methods/ann/optimizer/steepest_descent.hpp
+++ b/src/mlpack/methods/ann/optimizer/steepest_descent.hpp
@@ -3,7 +3,7 @@
  * @author Marcus Edel
  *
  * Intialization rule for the neural networks. This simple initialization is
- * performed by assigning a random matrix to the weight matrix. 
+ * performed by assigning a random matrix to the weight matrix.
  */
 #ifndef __MLPACK_METHOS_ANN_OPTIMIZER_STEEPEST_DESCENT_HPP
 #define __MLPACK_METHOS_ANN_OPTIMIZER_STEEPEST_DESCENT_HPP
@@ -31,8 +31,8 @@ class SteepestDescent
    */
   SteepestDescent(const size_t cols,
                   const size_t rows,
-                  const double lr = 1, 
-                  const double mom = 0.1) : 
+                  const double lr = 1,
+                  const double mom = 0.1) :
       lr(lr), mom(mom)
   {
     if (mom > 0)
@@ -69,5 +69,3 @@ class SteepestDescent
 }; // namespace mlpack
 
 #endif
-
-
diff --git a/src/mlpack/methods/ann/rnn.hpp b/src/mlpack/methods/ann/rnn.hpp
index c616aeb6e9..a4a085b583 100644
--- a/src/mlpack/methods/ann/rnn.hpp
+++ b/src/mlpack/methods/ann/rnn.hpp
@@ -48,7 +48,8 @@ class RNN
      * @param outputLayer The outputlayer used to evaluate the network.
      */
     RNN(const ConnectionTypes& network, OutputLayerType& outputLayer) :
-        network(network), outputLayer(outputLayer)
+        network(network), err(0),  trainError(0), seqNum(0),
+        outputLayer(outputLayer)
     {
       // Nothing to do here.
     }
@@ -71,31 +72,66 @@ class RNN
     {
       // Initialize the activation storage only once.
       if (!activations.size())
-        InitLayer(network, input, target);
+      {
+        InitLayer(network, input);
+      }
+      else
+      {
+        // Expand the activation storage to handle sequences of
+        // different length.
+        if (activations[0].n_cols < input.n_elem)
+        {
+          for (size_t i = 0; i < activations.size(); i++)
+          {
+            activations[i].insert_cols(activations[i].n_cols,
+                arma::zeros<MatType>(activations[i].n_rows,
+                input.n_elem - activations[i].n_cols));
+          }
+        }
+      }
 
-      // Reset the overall error.
-      err = 0;
-      error = MatType(target.n_elem, input.n_rows);
+      seqLen = input.n_rows / inputSize;
+      seqOutput = outputSize < target.n_elem ? true : false;
+      error = MatType(outputSize, outputSize < target.n_elem ? seqLen : 1);
 
       // Iterate through the input sequence and perform the feed forward pass.
-      for (seqNum = 0; seqNum < input.n_rows; seqNum++)
+      for (seqNum = 0; seqNum < seqLen; seqNum++)
       {
-        // Reset the network by zeroing the layer activations and set the input
-        // activation.
+        // Reset the network by zeroing the layer activations.
         ResetActivations(network);
+
+        // Set the current input activation.
         std::get<0>(std::get<0>(
-            network)).InputLayer().InputActivation() = input(seqNum);
+            network)).InputLayer().InputActivation() = input.submat(
+            seqNum * inputSize, 0, (seqNum + 1) * inputSize - 1, 0);
 
-        arma::colvec seqError = error.unsafe_col(seqNum);
-        FeedForward(network, target, seqError);
+        // Perform the forward pass and calculate the output error.
+        FeedForward(network);
+        if (seqOutput)
+        {
+          arma::colvec seqError = error.unsafe_col(seqNum);
+          arma::colvec seqTarget = target.subvec(seqNum * outputSize,
+              (seqNum + 1) * outputSize - 1);
 
-        // Save the network activation for the backward pass.
-        if (seqNum < (input.n_rows - 1))
+          OutputError(network, seqTarget, seqError);
+        }
+
+        // Save the network activation for the backward/forward pass and update
+        // the recurrent connections.
+        if (seqNum < (input.n_rows / inputSize - 1))
         {
           layerNum = 0;
           SaveActivations(network);
         }
       }
+
+      // Calculate the error only once for a non-sequence input.
+      if (!seqOutput)
+      {
+        seqNum = 0;
+        arma::colvec seqError = error.unsafe_col(seqNum);
+        OutputError(network, target, seqError);
+      }
     }
 
     /**
@@ -106,19 +142,20 @@ class RNN
      */
     void FeedBackward(const MatType& error)
     {
-      // Reset the network gradients by zeroing the storage.
-      for (size_t i = 0; i < gradients.size(); ++i)
-        gradients[i].zeros();
-
       // Reset the network deltas by zeroing the storage.
-      for (size_t i = 0; i < delta.size(); ++i)
-        delta[i].zeros();
+      for (size_t i = 0; i < delta.size(); i++)
+          delta[i].zeros();
 
       // Iterate through the input sequence and perform the feed backward pass.
-      for (seqNum = error.n_cols - 1; seqNum >= 0; seqNum--)
+      for (seqNum = seqLen - 1; seqNum >= 0; seqNum--)
       {
         gradientNum = 0;
-        FeedBackward(network, error.unsafe_col(seqNum));
+        deltaNum = 0;
+
+        // Perform the backward pass and update the gradient storage.
+        arma::colvec seqError = error.unsafe_col(seqOutput ? seqNum : 0);
+        FeedBackward(network, seqError);
+        UpdateGradients(network);
 
         // Load the network activation for the upcoming backward pass.
         if (seqNum > 0)
@@ -126,6 +163,7 @@ class RNN
           layerNum = 0;
           LoadActivations(network);
         }
+        else if (seqNum == 0) break;
       }
     }
 
@@ -137,10 +175,79 @@ class RNN
     {
       gradientNum = 0;
       ApplyGradients(network);
+
+      // Reset the overall error.
+      err = 0;
+      trainError = 0;
+      seqNum = 0;
+    }
+
+    /**
+     * Evaluate the network using the given input. The output activation is
+     * stored into the output parameter.
+     *
+     * @param input Input data used to evaluate the network.
+     * @param output Output data used to store the output activation
+     * @tparam VecType Type of data (arma::colvec, arma::mat or arma::sp_mat).
+     */
+    template <typename VecType>
+    void Predict(const VecType& input, VecType& output)
+    {
+      seqLen = input.n_rows / inputSize;
+
+      // Iterate through the input sequence and perform the feed forward pass.
+      for (seqNum = 0; seqNum < seqLen; seqNum++)
+      {
+        // Reset the network by zeroing the layer activations.
+        ResetActivations(network);
+
+        // Set the current input activation.
+        std::get<0>(std::get<0>(
+            network)).InputLayer().InputActivation() = input.submat(
+            seqNum * inputSize, 0, (seqNum + 1) * inputSize - 1, 0);
+
+        // Perform the forward pass and calculate the output error.
+        FeedForward(network);
+        if (seqOutput)
+        {
+          arma::colvec targetCol;
+          OutputPrediction(network, targetCol);
+          output = arma::join_cols(output, targetCol);
+        }
+
+        // Save the network activation for the backward/forward pass and update
+        // the recurrent connections.
+        if (seqNum < (input.n_rows / inputSize - 1))
+        {
+          layerNum = 0;
+          SaveActivations(network);
+        }
+      }
+
+      if (!seqOutput)
+        OutputPrediction(network, output);
+    }
+
+    /**
+     * Evaluate the trained network using the given input and compare the output
+     * with the given target vector.
+     *
+     * @param input Input data used to evaluate the trained network.
+     * @param target Target data used to calculate the network error.
+     * @param error The calulated error of the output layer.
+     * @tparam VecType Type of data (arma::colvec, arma::mat or arma::sp_mat).
+     */
+    template <typename VecType>
+    double Evaluate(const MatType& input,
+                     const VecType& target,
+                     MatType& error)
+    {
+      FeedForward(input, target, error);
+      return err;
     }
 
     //! Get the error of the network.
-    double Error() const { return err; }
+    double Error() const { return trainError; }
 
   private:
     /**
@@ -179,6 +286,17 @@ class RNN
     {
       std::get<I>(t).OutputLayer().InputActivation().zeros(
           std::get<I>(t).OutputLayer().InputSize());
+
+      // Reset the recurrent connection only at the beginning of a new sequence.
+      if (seqNum == 0 && (ConnectionTraits<typename std::remove_reference<
+          decltype(std::get<I>(t))>::type>::IsSelfConnection ||
+          ConnectionTraits<typename std::remove_reference<decltype(
+          std::get<I>(t))>::type>::IsFullselfConnection))
+      {
+        std::get<I>(t).InputLayer().InputActivation().zeros(
+          std::get<I>(t).InputLayer().InputSize());
+      }
+
       Reset<I + 1, Tp...>(t);
     }
 
@@ -192,38 +310,13 @@ class RNN
      * connections, and one for the general case which peels off the first type
      * and recurses, as usual with variadic function templates.
      */
-    template<size_t I = 0,
-             typename TargetVecType,
-             typename ErrorVecType,
-             typename... Tp>
+    template<size_t I = 0, typename... Tp>
     typename std::enable_if<I == sizeof...(Tp), void>::type
-    FeedForward(std::tuple<Tp...>& t,
-                TargetVecType& target,
-                ErrorVecType& error)
-    {
-      // Calculate and store the output error.
-      outputLayer.calculateError(std::get<0>(
-          std::get<I - 1>(t)).OutputLayer().InputActivation(), target,
-          error);
+    FeedForward(std::tuple<Tp...>& /* unused */) { }
 
-      // Save the output activation for the upcoming feed backward pass.
-      activations.back().unsafe_col(seqNum) = std::get<0>(
-          std::get<I - 1>(t)).OutputLayer().InputActivation();
-
-      // Masures the network's performance with the specified performance
-      // function.
-      err = PerformanceFunction::error(std::get<0>(
-          std::get<I - 1>(t)).OutputLayer().InputActivation(), target);
-    }
-
-    template<size_t I = 0,
-            typename TargetVecType,
-            typename ErrorVecType,
-            typename... Tp>
+    template<size_t I = 0, typename... Tp>
     typename std::enable_if<I < sizeof...(Tp), void>::type
-    FeedForward(std::tuple<Tp...>& t,
-                TargetVecType& target,
-                ErrorVecType& error)
+    FeedForward(std::tuple<Tp...>& t)
     {
       Forward(std::get<I>(t));
 
@@ -232,7 +325,7 @@ class RNN
           std::get<0>(std::get<I>(t)).OutputLayer().InputActivation(),
           std::get<0>(std::get<I>(t)).OutputLayer().InputActivation());
 
-      FeedForward<I + 1, TargetVecType, ErrorVecType, Tp...>(t, target, error);
+      FeedForward<I + 1, Tp...>(t);
     }
 
     /**
@@ -254,6 +347,45 @@ class RNN
       Forward<I + 1, Tp...>(t);
     }
 
+    /*
+     * Calculate the output error and update the overall error.
+     */
+    template<typename VecType, typename... Tp>
+    void OutputError(std::tuple<Tp...>& t,
+                     const VecType& target,
+                     VecType& error)
+    {
+       // Calculate and store the output error.
+      outputLayer.calculateError(std::get<0>(
+          std::get<sizeof...(Tp) - 1>(t)).OutputLayer().InputActivation(),
+          target, error);
+
+      // Save the output activation for the upcoming feed backward pass.
+      activations.back().unsafe_col(seqNum) = std::get<0>(
+          std::get<sizeof...(Tp) - 1>(t)).OutputLayer().InputActivation();
+
+      // Masures the network's performance with the specified performance
+      // function.
+      err = PerformanceFunction::error(std::get<0>(
+          std::get<sizeof...(Tp) - 1>(t)).OutputLayer().InputActivation(),
+          target);
+
+      // Update the overall training error.
+      trainError += err;
+    }
+
+    /*
+     * Calculate and store the output activation.
+     */
+    template<typename VecType, typename... Tp>
+    void OutputPrediction(std::tuple<Tp...>& t, VecType& output)
+    {
+       // Calculate and store the output prediction.
+      outputLayer.outputClass(std::get<0>(
+          std::get<sizeof...(Tp) - 1>(t)).OutputLayer().InputActivation(),
+          output);
+    }
+
     /**
      * Run a single iteration of the feed backward algorithm, using the given
      * error of the output layer. Note that we iterate backward through the
@@ -265,11 +397,11 @@ class RNN
      * and recurses, as usual with variadic function templates.
      */
     template<size_t I = 0, typename VecType, typename... Tp>
-    typename std::enable_if<I == sizeof...(Tp), void>::type
+    typename std::enable_if<I == sizeof...(Tp) + 1, void>::type
     FeedBackward(std::tuple<Tp...>& /* unused */, VecType& /* unused */) { }
 
     template<size_t I = 1, typename VecType, typename... Tp>
-    typename std::enable_if<I < sizeof...(Tp), void>::type
+    typename std::enable_if<I < sizeof...(Tp) + 1, void>::type
     FeedBackward(std::tuple<Tp...>& t, VecType& error)
     {
       // Distinguish between the output layer and the other layer. In case of
@@ -282,20 +414,10 @@ class RNN
         std::get<0>(std::get<sizeof...(Tp) - I>(t)).OutputLayer().FeedBackward(
             activations.back().unsafe_col(seqNum), error,
             std::get<0>(std::get<sizeof...(Tp) - I>(t)).OutputLayer().Delta());
-
-        // Save the delta for the upcoming feed backward pass.
-        delta.back() += std::get<0>(
-            std::get<sizeof...(Tp) - I>(t)).OutputLayer().Delta();
-
-        // Save the gradient to update the weights at the end.
-        gradients.back() += std::get<0>(
-            std::get<sizeof...(Tp) - I>(t)).OutputLayer().Delta() *
-            std::get<0>(
-            std::get<sizeof...(Tp) - I>(t)).InputLayer().InputActivation().t();
       }
 
-      Backward(std::get<sizeof...(Tp) - I>(t), delta[delta.size() - I]);
-      UpdateGradients(std::get<sizeof...(Tp) - I - 1>(t));
+      Backward(std::get<sizeof...(Tp) - I>(t), std::get<0>(std::get<
+          sizeof...(Tp) - I>(t)).OutputLayer().Delta(), I, sizeof...(Tp));
 
       FeedBackward<I + 1, VecType, Tp...>(t, error);
     }
@@ -310,28 +432,71 @@ class RNN
      */
     template<size_t I = 0, typename VecType, typename... Tp>
     typename std::enable_if<I == sizeof...(Tp), void>::type
-    Backward(std::tuple<Tp...>& /* unused */, VecType& /* unused */) { }
+    Backward(std::tuple<Tp...>& /* unused */,
+             VecType& /* unused */,
+             const size_t /* unused */,
+             const size_t /* unused */) { }
 
     template<size_t I = 0, typename VecType, typename... Tp>
     typename std::enable_if<I < sizeof...(Tp), void>::type
-    Backward(std::tuple<Tp...>& t, VecType& error)
+    Backward(std::tuple<Tp...>& t,
+             VecType& error,
+             const size_t layer,
+             const size_t layerNum)
     {
       std::get<I>(t).FeedBackward(error);
 
+      // Update the recurrent delta.
+      if (ConnectionTraits<typename std::remove_reference<decltype(
+          std::get<I>(t))>::type>::IsSelfConnection)
+      {
+        std::get<I>(t).FeedBackward(delta[deltaNum]);
+        delta[deltaNum++] = std::get<I>(t).Delta();
+      }
+
       // We calculate the delta only for non bias layer and self connections.
       if (!(ConnectionTraits<typename std::remove_reference<decltype(
             std::get<I>(t))>::type>::IsSelfConnection ||
         LayerTraits<typename std::remove_reference<decltype(
             std::get<I>(t).InputLayer())>::type>::IsBiasLayer ||
         ConnectionTraits<typename std::remove_reference<decltype(
-            std::get<I>(t))>::type>::IsFullselfConnection))
+            std::get<I>(t))>::type>::IsFullselfConnection) && layer < layerNum)
       {
+        // Sum up the stored delta for recurrent connections.
+        if (recurrentLayer[layer])
+          std::get<I>(t).Delta() += delta[deltaNum];
+
+        // Perform the backward pass.
         std::get<I>(t).InputLayer().FeedBackward(
             std::get<I>(t).InputLayer().InputActivation(),
             std::get<I>(t).Delta(), std::get<I>(t).InputLayer().Delta());
+
+        // Update the delta storage for the next backward pass.
+        if (recurrentLayer[layer])
+          delta[deltaNum] = std::get<I>(t).InputLayer().Delta();
       }
 
-      Backward<I + 1, VecType, Tp...>(t, error);
+      Backward<I + 1, VecType, Tp...>(t, error, layer, layerNum);
+    }
+
+    /**
+     * Helper function to update the gradient storage.
+     *
+     * enable_if (SFINAE) is used to select between two template overloads of
+     * the get function - one for when I is equal the size of the tuple of
+     * connections, and one for the general case which peels off the first type
+     * and recurses, as usual with variadic function templates.
+     */
+    template<size_t I = 0, typename... Tp>
+    typename std::enable_if<I == sizeof...(Tp), void>::type
+    UpdateGradients(std::tuple<Tp...>& /* unused */) { }
+
+    template<size_t I = 0, typename... Tp>
+    typename std::enable_if<I < sizeof...(Tp), void>::type
+    UpdateGradients(std::tuple<Tp...>& t)
+    {
+      Gradients(std::get<I>(t));
+      UpdateGradients<I + 1, Tp...>(t);
     }
 
     /**
@@ -343,16 +508,17 @@ class RNN
      */
     template<size_t I = 0, typename... Tp>
     typename std::enable_if<I == sizeof...(Tp), void>::type
-    UpdateGradients(std::tuple<Tp...>& /* unused */) { }
+    Gradients(std::tuple<Tp...>& /* unused */) { }
 
     template<size_t I = 0, typename... Tp>
     typename std::enable_if<I < sizeof...(Tp), void>::type
-    UpdateGradients(std::tuple<Tp...>& t)
+    Gradients(std::tuple<Tp...>& t)
     {
-      gradients[gradientNum++] += std::get<I>(t).OutputLayer().Delta() *
-          std::get<I>(t).InputLayer().InputActivation().t();
+      MatType gradient;
+      std::get<I>(t).Gradient(gradient);
+      gradients[gradientNum++] += gradient;
 
-      UpdateGradients<I + 1, Tp...>(t);
+      Gradients<I + 1, Tp...>(t);
     }
 
     /**
@@ -365,14 +531,14 @@ class RNN
      * and recurses, as usual with variadic function templates.
      */
     template<size_t I = 0, typename... Tp>
-    typename std::enable_if<I == sizeof...(Tp) - 1, void>::type
+    typename std::enable_if<I == sizeof...(Tp), void>::type
     ApplyGradients(std::tuple<Tp...>& /* unused */) { }
 
     template<size_t I = 0, typename... Tp>
-    typename std::enable_if<I < sizeof...(Tp) - 1, void>::type
+    typename std::enable_if<I < sizeof...(Tp), void>::type
     ApplyGradients(std::tuple<Tp...>& t)
     {
-      Gradients(std::get<I>(t));
+      Apply(std::get<I>(t));
       ApplyGradients<I + 1, Tp...>(t);
     }
 
@@ -387,16 +553,19 @@ class RNN
      */
     template<size_t I = 0, typename... Tp>
     typename std::enable_if<I == sizeof...(Tp), void>::type
-    Gradients(std::tuple<Tp...>& /* unused */) { }
+    Apply(std::tuple<Tp...>& /* unused */) { }
 
     template<size_t I = 0, typename... Tp>
     typename std::enable_if<I < sizeof...(Tp), void>::type
-    Gradients(std::tuple<Tp...>& t)
+    Apply(std::tuple<Tp...>& t)
     {
       std::get<I>(t).Optimzer().UpdateWeights(std::get<I>(t).Weights(),
-          gradients[gradientNum++], err);
+          gradients[gradientNum], trainError);
 
-      Gradients<I + 1, Tp...>(t);
+      // // Reset the gradient storage.
+      gradients[gradientNum++].zeros();
+
+      Apply<I + 1, Tp...>(t);
     }
 
     /**
@@ -408,21 +577,49 @@ class RNN
      * connections, and one for the general case which peels off the first type
      * and recurses, as usual with variadic function templates.
      */
-    template<size_t I = 0, typename VecType, typename... Tp>
+    template<size_t I = 0, typename... Tp>
     typename std::enable_if<I == sizeof...(Tp), void>::type
-    InitLayer(std::tuple<Tp...>& /* unused */,
-              const MatType& input,
-              const VecType& target)
+    InitLayer(std::tuple<Tp...>& t, const MatType& input)
     {
-      activations.push_back(new MatType(target.n_elem, input.n_elem));
+      recurrentLayer.push_back(false);
+      outputSize = std::get<0>(std::get<I - 1>(t)).OutputLayer().OutputSize();
+      activations.push_back(new MatType(outputSize, input.n_elem));
     }
 
-    template<size_t I = 0, typename VecType, typename... Tp>
+    template<size_t I = 0, typename... Tp>
     typename std::enable_if<I < sizeof...(Tp), void>::type
-    InitLayer(std::tuple<Tp...>& t, const MatType& input, const VecType& target)
+    InitLayer(std::tuple<Tp...>& t, const MatType& input)
     {
-      Layer(std::get<I>(t), input);
-      InitLayer<I + 1, VecType, Tp...>(t, input, target);
+      if (I == 0)
+        inputSize = std::get<0>(std::get<I>(t)).InputLayer().InputSize();
+
+      recurrentLayer.push_back(false);
+      Recurrent(std::get<sizeof...(Tp) - I - 1>(t));
+
+      Layer(std::get<I>(t), input, I);
+      InitLayer<I + 1, Tp...>(t, input);
+    }
+
+    template<size_t I = 0, typename... Tp>
+    typename std::enable_if<I == sizeof...(Tp), void>::type
+    Recurrent(std::tuple<Tp...>& /* unusded */) { }
+
+    template<size_t I = 0, typename... Tp>
+    typename std::enable_if<I < sizeof...(Tp), void>::type
+    Recurrent(std::tuple<Tp...>& t)
+    {
+      if (ConnectionTraits<typename std::remove_reference<decltype(
+              std::get<I>(t))>::type>::IsSelfConnection ||
+          ConnectionTraits<typename std::remove_reference<decltype(
+            std::get<I>(t))>::type>::IsFullselfConnection)
+      {
+        recurrentLayer.back() = true;
+        delta.push_back(new VecTypeDelta(std::get<I>(t).Weights().n_rows));
+      }
+      else
+      {
+        Recurrent<I + 1, Tp...>(t);
+      }
     }
 
     /**
@@ -436,30 +633,21 @@ class RNN
      */
     template<size_t I = 0, typename VecType, typename... Tp>
     typename std::enable_if<I == sizeof...(Tp), void>::type
-    Layer(std::tuple<Tp...>& /* unused */, const VecType& /* unused */) { }
+    Layer(std::tuple<Tp...>& /* unusded */,
+          const VecType& /* unused */,
+          const size_t /* unsued */) { }
 
     template<size_t I = 0, typename VecType, typename... Tp>
     typename std::enable_if<I < sizeof...(Tp), void>::type
-    Layer(std::tuple<Tp...>& t, const VecType& input)
+    Layer(std::tuple<Tp...>& t, const VecType& input, const size_t layer)
     {
       activations.push_back(new MatType(
         std::get<I>(t).InputLayer().OutputSize(), input.n_elem));
 
       gradients.push_back(new MatType(std::get<I>(t).Weights().n_rows,
-          std::get<I>(t).Weights().n_cols));
+          std::get<I>(t).Weights().n_cols, arma::fill::zeros));
 
-      // We calculate the delta only for non bias layer and self connections.
-      if (!(ConnectionTraits<typename std::remove_reference<decltype(
-              std::get<I>(t))>::type>::IsSelfConnection ||
-          LayerTraits<typename std::remove_reference<decltype(
-              std::get<I>(t).InputLayer())>::type>::IsBiasLayer ||
-          ConnectionTraits<typename std::remove_reference<decltype(
-            std::get<I>(t))>::type>::IsFullselfConnection))
-      {
-        delta.push_back(new VecTypeDelta(std::get<I>(t).Weights().n_rows));
-      }
-
-      Layer<I + 1, VecType, Tp...>(t, input);
+      Layer<I + 1, VecType, Tp...>(t, input, layer);
     }
 
     /**
@@ -554,20 +742,23 @@ class RNN
       Save<I + 1, Tp...>(t);
     }
 
+    //! The layer we are using to build the network.
+    ConnectionTypes network;
+
     //! The current error of the network.
     double err;
 
+    //! The current training error of the network.
+    double trainError;
+
     //! The activation storage we are using to perform the feed backward pass.
     boost::ptr_vector<MatType> activations;
 
     //! The gradient storage we are using to perform the feed backward pass.
     boost::ptr_vector<MatType> gradients;
 
-    //! The detla storage we are using to perform the feed backward pass.
-    boost::ptr_vector<VecTypeDelta> delta;
-
     //! The index of the current sequence number.
-    long int seqNum;
+    size_t seqNum;
 
     //! The index of the currently activate layer.
     size_t layerNum;
@@ -575,11 +766,29 @@ class RNN
     //! The index of the currently activate gradient.
     size_t gradientNum;
 
-    //! The layer we are using to build the network.
-    ConnectionTypes network;
+    //! The index of the currently activate delta.
+    size_t deltaNum;
+
+    //! Locally stored network output size.
+    size_t outputSize;
+
+    //! Locally stored network input size.
+    size_t inputSize;
+
+    //! Locally stored parameter that indicates if the input is a sequence.
+    bool seqOutput;
 
     //! The outputlayer used to evaluate the network
     OutputLayerType& outputLayer;
+
+    //! Locally stored number of samples in one input sequence.
+    size_t seqLen;
+
+    //! The recurrentLayer storage we are using to perform the backward pass.
+    std::vector<bool> recurrentLayer;
+
+    //! The detla storage we are using to perform the feed backward pass.
+    boost::ptr_vector<VecTypeDelta> delta;
 }; // class RNN
 
 //! Network traits for the FFNN network.
@@ -599,4 +808,3 @@ class NetworkTraits<RNN<ConnectionTypes, OutputLayerType, PerformanceFunction> >
 }; // namespace mlpack
 
 #endif
-
diff --git a/src/mlpack/methods/ann/trainer/trainer.hpp b/src/mlpack/methods/ann/trainer/trainer.hpp
index f5b42a384a..01c1301623 100644
--- a/src/mlpack/methods/ann/trainer/trainer.hpp
+++ b/src/mlpack/methods/ann/trainer/trainer.hpp
@@ -87,8 +87,6 @@ class Trainer
 
       while(true)
       {
-
-        // Randomly shuffle the index sequence if not in batch mode.
         if (shuffle)
           index = arma::shuffle(index);
 
@@ -98,7 +96,7 @@ class Trainer
         if (validationError <= tolerance)
           break;
 
-        if (maxEpochs > 0 && ++epoch > maxEpochs)
+        if (maxEpochs > 0 && ++epoch >= maxEpochs)
           break;
       }
     }
@@ -145,8 +143,8 @@ class Trainer
       {
         net.FeedForward(data.unsafe_col(index(i)),
             target.unsafe_col(index(i)), error);
-        trainingError += net.Error();
 
+        trainingError += net.Error();
         net.FeedBackward(error);
 
         if (((i + 1) % batchSize) == 0)
@@ -172,8 +170,8 @@ class Trainer
 
       for (size_t i = 0; i < data.n_cols; i++)
       {
-        net.FeedForward(data.unsafe_col(i), target.unsafe_col(i), error);
-        validationError += net.Error();
+         validationError += net.Evaluate(data.unsafe_col(i),
+            target.unsafe_col(i), error);
       }
 
       validationError /= data.n_cols;
diff --git a/src/mlpack/tests/feedforward_network_test.cpp b/src/mlpack/tests/feedforward_network_test.cpp
index 6f8e338b51..5f67fbc01d 100644
--- a/src/mlpack/tests/feedforward_network_test.cpp
+++ b/src/mlpack/tests/feedforward_network_test.cpp
@@ -36,30 +36,50 @@ using namespace mlpack::ann;
 
 BOOST_AUTO_TEST_SUITE(FeedForwardNetworkTest);
 
-
 /**
  * Train and evaluate a vanilla network with the specified structure.
  */
-template<typename WeightInitRule,
-         typename PerformanceFunction,
-         typename OptimizerType,
-         typename OutputLayerType,
-         typename PerformanceFunctionType,
-         typename MatType = arma::mat,
-         typename VecType = arma::colvec
+template<
+    typename WeightInitRule,
+    typename PerformanceFunction,
+    typename OptimizerType,
+    typename OutputLayerType,
+    typename PerformanceFunctionType,
+    typename MatType = arma::mat,
+    typename VecType = arma::colvec
 >
 void BuildVanillaNetwork(MatType& trainData,
                          MatType& trainLabels,
                          MatType& testData,
                          MatType& testLabels,
-                         size_t hiddenLayerSize,
-                         size_t maxEpochs,
-                         double classificationErrorThreshold,
-                         double ValidationErrorThreshold,
+                         const size_t hiddenLayerSize,
+                         const size_t maxEpochs,
+                         const double classificationErrorThreshold,
+                         const double ValidationErrorThreshold,
                          WeightInitRule weightInitRule = WeightInitRule())
 {
+  /*
+   * Construct a feed forward network with trainData.n_rows input nodes,
+   * hiddenLayerSize hidden nodes and trainLabels.n_rows output nodes. The
+   * network structure looks like:
+   *
+   *  Input         Hidden        Output
+   *  Layer         Layer         Layer
+   * +-----+       +-----+       +-----+
+   * |     |       |     |       |     |
+   * |     +------>|     +------>|     |
+   * |     |     +>|     |       |     |
+   * +-----+     | +--+--+       +-----+
+   *             |
+   *  Bias       |
+   *  Layer      |
+   * +-----+     |
+   * |     |     |
+   * |     +-----+
+   * |     |
+   * +-----+
+   */
   BiasLayer<> biasLayer0(1);
-  BiasLayer<> biasLayer1(1);
 
   NeuronLayer<PerformanceFunction> inputLayer(trainData.n_rows);
   NeuronLayer<PerformanceFunction> hiddenLayer0(hiddenLayerSize);
@@ -169,7 +189,7 @@ BOOST_AUTO_TEST_CASE(VanillaNetworkTest)
                       SteepestDescent<>,
                       BinaryClassificationLayer<>,
                       MeanSquaredErrorFunction<> >
-      (dataset, labels, dataset, labels, 100, 450, 0.6, 90, randInitB);
+      (dataset, labels, dataset, labels, 100, 100, 0.6, 10, randInitB);
 
   // Vanilla neural net with tanh activation function.
   BuildVanillaNetwork<RandomInitialization<>,
@@ -177,7 +197,7 @@ BOOST_AUTO_TEST_CASE(VanillaNetworkTest)
                     SteepestDescent<>,
                     BinaryClassificationLayer<>,
                     MeanSquaredErrorFunction<> >
-    (dataset, labels, dataset, labels, 10, 450, 0.6, 90, randInitB);
+    (dataset, labels, dataset, labels, 10, 200, 0.6, 20, randInitB);
 }
 
 /**
@@ -237,13 +257,13 @@ BOOST_AUTO_TEST_CASE(VanillaNetworkConvergenceTest)
  * Train a vanilla network with the specified structure step by step and
  * evaluate the network.
  */
-template<typename WeightInitRule,
-         typename PerformanceFunction,
-         typename OptimizerType,
-         typename OutputLayerType,
-         typename PerformanceFunctionType,
-         typename MatType = arma::mat,
-         typename VecType = arma::colvec
+template<
+    typename WeightInitRule,
+    typename PerformanceFunction,
+    typename OptimizerType,
+    typename OutputLayerType,
+    typename PerformanceFunctionType,
+    typename MatType = arma::mat
 >
 void BuildNetworkOptimzer(MatType& trainData,
                           MatType& trainLabels,
@@ -253,8 +273,28 @@ void BuildNetworkOptimzer(MatType& trainData,
                           size_t epochs,
                           WeightInitRule weightInitRule = WeightInitRule())
 {
+  /*
+   * Construct a feed forward network with trainData.n_rows input nodes,
+   * hiddenLayerSize hidden nodes and trainLabels.n_rows output nodes. The
+   * network structure looks like:
+   *
+   *  Input         Hidden        Output
+   *  Layer         Layer         Layer
+   * +-----+       +-----+       +-----+
+   * |     |       |     |       |     |
+   * |     +------>|     +------>|     |
+   * |     |     +>|     |       |     |
+   * +-----+     | +--+--+       +-----+
+   *             |
+   *  Bias       |
+   *  Layer      |
+   * +-----+     |
+   * |     |     |
+   * |     +-----+
+   * |     |
+   * +-----+
+   */
   BiasLayer<> biasLayer0(1);
-  BiasLayer<> biasLayer1(1);
 
   NeuronLayer<PerformanceFunction> inputLayer(trainData.n_rows);
   NeuronLayer<PerformanceFunction> hiddenLayer0(hiddenLayerSize);
diff --git a/src/mlpack/tests/performance_functions_test.cpp b/src/mlpack/tests/performance_functions_test.cpp
index df425c3fc0..82372b0ab9 100644
--- a/src/mlpack/tests/performance_functions_test.cpp
+++ b/src/mlpack/tests/performance_functions_test.cpp
@@ -47,4 +47,3 @@ BOOST_AUTO_TEST_CASE(SumSquaredErrorTest)
 }
 
 BOOST_AUTO_TEST_SUITE_END();
-
diff --git a/src/mlpack/tests/recurrent_network_test.cpp b/src/mlpack/tests/recurrent_network_test.cpp
new file mode 100644
index 0000000000..24f44fac4c
--- /dev/null
+++ b/src/mlpack/tests/recurrent_network_test.cpp
@@ -0,0 +1,646 @@
+/**
+ * @file feedforward_network_test.cpp
+ * @author Marcus Edel
+ *
+ * Tests the feed forward network.
+ */
+#include <mlpack/core.hpp>
+
+#include <mlpack/methods/ann/activation_functions/logistic_function.hpp>
+#include <mlpack/methods/ann/activation_functions/identity_function.hpp>
+#include <mlpack/methods/ann/activation_functions/softsign_function.hpp>
+#include <mlpack/methods/ann/activation_functions/tanh_function.hpp>
+#include <mlpack/methods/ann/activation_functions/rectifier_function.hpp>
+
+#include <mlpack/methods/ann/init_rules/random_init.hpp>
+#include <mlpack/methods/ann/init_rules/orthogonal_init.hpp>
+#include <mlpack/methods/ann/init_rules/oivs_init.hpp>
+#include <mlpack/methods/ann/init_rules/kathirvalavakumar_subavathi_init.hpp>
+#include <mlpack/methods/ann/init_rules/nguyen_widrow_init.hpp>
+
+#include <mlpack/methods/ann/layer/neuron_layer.hpp>
+#include <mlpack/methods/ann/layer/bias_layer.hpp>
+#include <mlpack/methods/ann/layer/binary_classification_layer.hpp>
+#include <mlpack/methods/ann/layer/multiclass_classification_layer.hpp>
+
+#include <mlpack/methods/ann/connections/full_connection.hpp>
+#include <mlpack/methods/ann/connections/self_connection.hpp>
+#include <mlpack/methods/ann/connections/fullself_connection.hpp>
+#include <mlpack/methods/ann/connections/connection_traits.hpp>
+
+#include <mlpack/methods/ann/trainer/trainer.hpp>
+
+#include <mlpack/methods/ann/ffnn.hpp>
+#include <mlpack/methods/ann/rnn.hpp>
+
+
+#include <mlpack/methods/ann/performance_functions/mse_function.hpp>
+#include <mlpack/methods/ann/performance_functions/sse_function.hpp>
+#include <mlpack/methods/ann/performance_functions/cee_function.hpp>
+
+#include <mlpack/methods/ann/optimizer/steepest_descent.hpp>
+#include <mlpack/methods/ann/optimizer/rpropp.hpp>
+
+#include <boost/test/unit_test.hpp>
+#include "old_boost_test_definitions.hpp"
+
+using namespace mlpack;
+using namespace mlpack::ann;
+
+
+BOOST_AUTO_TEST_SUITE(RecurrentNetworkTest);
+
+// Be careful!  When writing new tests, always get the boolean value and store
+// it in a temporary, because the Boost unit test macros do weird things and
+// will cause bizarre problems.
+
+/**
+ * Construct a 2-class dataset out of noisy sines.
+ *
+ * @param data Input data used to store the noisy sines.
+ * @param labels Labels used to store the target class of the noisy sines.
+ * @param points Number of points/features in a single sequence.
+ * @param sequences Number of sequences for each class.
+ * @param noise The noise factor that influences the sines.
+ */
+void GenerateNoisySines(arma::mat& data,
+                        arma::mat& labels,
+                        const size_t points,
+                        const size_t sequences,
+                        const double noise = 0.3)
+{
+  arma::colvec x =  arma::linspace<arma::Col<double> >(0,
+      points - 1, points) / points * 20.0;
+  arma::colvec y1 = arma::sin(x + arma::as_scalar(arma::randu(1)) * 3.0);
+  arma::colvec y2 = arma::sin(x / 2.0 + arma::as_scalar(arma::randu(1)) * 3.0);
+
+  data = arma::zeros(points, sequences * 2);
+  labels = arma::zeros(2, sequences * 2);
+
+  for (size_t seq = 0; seq < sequences; seq++)
+  {
+    data.col(seq) = arma::randu(points) * noise + y1 +
+        arma::as_scalar(arma::randu(1) - 0.5) * noise;
+    labels(0, seq) = 1;
+
+    data.col(sequences + seq) = arma::randu(points) * noise + y2 +
+        arma::as_scalar(arma::randu(1) - 0.5) * noise;
+    labels(1, sequences + seq) = 1;
+  }
+}
+
+/**
+ * Train the vanilla network on a larger dataset.
+ */
+BOOST_AUTO_TEST_CASE(SequenceClassificationTest)
+{
+  // Generate 12 (2 * 6) noisy sines. A single sine contains 10 points/features.
+  arma::mat input, labels;
+  GenerateNoisySines(input, labels, 10, 6);
+
+  /*
+   * Construct a network with 1 input unit, 4 hidden units and 2 output units.
+   * The hidden layer is connected to itself. The network structure looks like:
+   *
+   *  Input         Hidden        Output
+   * Layer(1)      Layer(4)      Layer(2)
+   * +-----+       +-----+       +-----+
+   * |     |       |     |       |     |
+   * |     +------>|     +------>|     |
+   * |     |    ..>|     |       |     |
+   * +-----+    .  +--+--+       +-----+
+   *            .     .
+   *            .     .
+   *            .......
+   */
+  NeuronLayer<LogisticFunction> inputLayer(1);
+  NeuronLayer<LogisticFunction> hiddenLayer0(4);
+  NeuronLayer<LogisticFunction> recurrentLayer0(hiddenLayer0.InputSize());
+  NeuronLayer<LogisticFunction> hiddenLayer1(2);
+  BinaryClassificationLayer<> outputLayer;
+
+  SteepestDescent< > conOptimizer0(inputLayer.InputSize(),
+      hiddenLayer0.InputSize(), 1, 0);
+  SteepestDescent< > conOptimizer2(hiddenLayer0.InputSize(),
+      hiddenLayer0.InputSize(), 1, 0);
+  SteepestDescent< > conOptimizer3(hiddenLayer0.InputSize(),
+      hiddenLayer1.OutputSize(), 1, 0);
+
+  NguyenWidrowInitialization<> randInit;
+
+  FullConnection<
+      decltype(inputLayer),
+      decltype(hiddenLayer0),
+      decltype(conOptimizer0),
+      decltype(randInit)>
+      layerCon0(inputLayer, hiddenLayer0, conOptimizer0, randInit);
+
+  SelfConnection<
+    decltype(recurrentLayer0),
+    decltype(hiddenLayer0),
+    decltype(conOptimizer2),
+    decltype(randInit)>
+    layerCon2(recurrentLayer0, hiddenLayer0, conOptimizer2, randInit);
+
+  FullConnection<
+      decltype(hiddenLayer0),
+      decltype(hiddenLayer1),
+      decltype(conOptimizer3),
+      decltype(randInit)>
+      layerCon4(hiddenLayer0, hiddenLayer1, conOptimizer3, randInit);
+
+  auto module0 = std::tie(layerCon0, layerCon2);
+  auto module1 = std::tie(layerCon4);
+  auto modules = std::tie(module0, module1);
+
+  RNN<decltype(modules),
+      decltype(outputLayer),
+      MeanSquaredErrorFunction<> > net(modules, outputLayer);
+
+  // Train the network for 1000 epochs.
+  Trainer<decltype(net)> trainer(net, 1000);
+  trainer.Train(input, labels, input, labels);
+
+  // Ask the network to classify the trained input data.
+  arma::colvec output;
+  for (size_t i = 0; i < input.n_cols; i++)
+  {
+    net.Predict(input.unsafe_col(i), output);
+
+    bool b = arma::all((output == labels.unsafe_col(i)) == 1);
+    BOOST_REQUIRE_EQUAL(b, 1);
+  }
+}
+
+/**
+ * Train and evaluate a vanilla feed forward network and a recurrent network
+ * with the specified structure and compare the two networks output and overall
+ * error.
+ */
+template<
+    typename WeightInitRule,
+    typename PerformanceFunction,
+    typename OptimizerType,
+    typename OutputLayerType,
+    typename PerformanceFunctionType,
+    typename MatType = arma::mat
+>
+void CompareVanillaNetworks(MatType& trainData,
+                            MatType& trainLabels,
+                            MatType& testData,
+                            MatType& testLabels,
+                            const size_t hiddenLayerSize,
+                            const size_t maxEpochs,
+                            WeightInitRule weightInitRule = WeightInitRule())
+{
+  BiasLayer<> biasLayer0(1);
+
+  NeuronLayer<PerformanceFunction> inputLayer(trainData.n_rows);
+  NeuronLayer<PerformanceFunction> hiddenLayer0(hiddenLayerSize);
+  NeuronLayer<PerformanceFunction> hiddenLayer1(trainLabels.n_rows);
+
+  OutputLayerType outputLayer;
+
+  OptimizerType ffnConOptimizer0(trainData.n_rows, hiddenLayerSize);
+  OptimizerType ffnConOptimizer1(1, hiddenLayerSize);
+  OptimizerType ffnConOptimizer2(hiddenLayerSize, trainLabels.n_rows);
+
+  OptimizerType rnnConOptimizer0(trainData.n_rows, hiddenLayerSize);
+  OptimizerType rnnConOptimizer1(1, hiddenLayerSize);
+  OptimizerType rnnConOptimizer2(hiddenLayerSize, trainLabels.n_rows);
+
+  FullConnection<
+    decltype(inputLayer),
+    decltype(hiddenLayer0),
+    decltype(ffnConOptimizer0),
+    decltype(weightInitRule)>
+    ffnLayerCon0(inputLayer, hiddenLayer0, ffnConOptimizer0, weightInitRule);
+
+  FullConnection<
+    decltype(inputLayer),
+    decltype(hiddenLayer0),
+    decltype(rnnConOptimizer0),
+    decltype(weightInitRule)>
+    rnnLayerCon0(inputLayer, hiddenLayer0, rnnConOptimizer0, weightInitRule);
+
+  FullConnection<
+    decltype(biasLayer0),
+    decltype(hiddenLayer0),
+    decltype(ffnConOptimizer1),
+    decltype(weightInitRule)>
+    ffnLayerCon1(biasLayer0, hiddenLayer0, ffnConOptimizer1, weightInitRule);
+
+  FullConnection<
+    decltype(biasLayer0),
+    decltype(hiddenLayer0),
+    decltype(rnnConOptimizer1),
+    decltype(weightInitRule)>
+    rnnLayerCon1(biasLayer0, hiddenLayer0, rnnConOptimizer1, weightInitRule);
+
+  FullConnection<
+      decltype(hiddenLayer0),
+      decltype(hiddenLayer1),
+      decltype(ffnConOptimizer2),
+      decltype(weightInitRule)>
+      ffnLayerCon2(hiddenLayer0, hiddenLayer1, ffnConOptimizer2, weightInitRule);
+
+  FullConnection<
+      decltype(hiddenLayer0),
+      decltype(hiddenLayer1),
+      decltype(rnnConOptimizer2),
+      decltype(weightInitRule)>
+      rnnLayerCon2(hiddenLayer0, hiddenLayer1, rnnConOptimizer2, weightInitRule);
+
+  auto ffnModule0 = std::tie(ffnLayerCon0, ffnLayerCon1);
+  auto ffnModule1 = std::tie(ffnLayerCon2);
+  auto ffnModules = std::tie(ffnModule0, ffnModule1);
+
+  auto rnnModule0 = std::tie(rnnLayerCon0, rnnLayerCon1);
+  auto rnnModule1 = std::tie(rnnLayerCon2);
+  auto rnnModules = std::tie(rnnModule0, rnnModule1);
+
+  /*
+   * Construct a feed forward network with trainData.n_rows input units,
+   * hiddenLayerSize hidden units and trainLabels.n_rows output units. The
+   * network structure looks like:
+   *
+   *  Input         Hidden        Output
+   *  Layer         Layer         Layer
+   * +-----+       +-----+       +-----+
+   * |     |       |     |       |     |
+   * |     +------>|     +------>|     |
+   * |     |       |     |       |     |
+   * +-----+       +--+--+       +-----+
+   */
+  FFNN<decltype(ffnModules), decltype(outputLayer), PerformanceFunctionType>
+      ffn(ffnModules, outputLayer);
+
+  /*
+   * Construct a recurrent network with trainData.n_rows input units,
+   * hiddenLayerSize hidden units and trainLabels.n_rows output units. The
+   * hidden layer is connected to itself. The network structure looks like:
+   *
+   *  Input         Hidden        Output
+   *  Layer         Layer         Layer
+   * +-----+       +-----+       +-----+
+   * |     |       |     |       |     |
+   * |     +------>|     +------>|     |
+   * |     |    ..>|     |       |     |
+   * +-----+    .  +--+--+       +-----+
+   *            .     .
+   *            .     .
+   *            .......
+   */
+  RNN<decltype(rnnModules), decltype(outputLayer), PerformanceFunctionType>
+      rnn(rnnModules, outputLayer);
+
+  // Train the network for maxEpochs epochs or until we reach a validation error
+  // of less then 0.001.
+  Trainer<decltype(ffn)> ffnTrainer(ffn, maxEpochs, 1, 0.001, false);
+  Trainer<decltype(rnn)> rnnTrainer(rnn, maxEpochs, 1, 0.001, false);
+
+  for (size_t i = 0; i < 5; i++)
+  {
+    rnnTrainer.Train(trainData, trainLabels, testData, testLabels);
+    ffnTrainer.Train(trainData, trainLabels, testData, testLabels);
+
+    if (!arma::is_finite(ffnTrainer.ValidationError()))
+      continue;
+
+    BOOST_REQUIRE_CLOSE(ffnTrainer.ValidationError(),
+        rnnTrainer.ValidationError(), 1e-3);
+  }
+}
+
+/**
+ * Train a vanilla feed forward and recurrent network on a sequence with len
+ * one. Ideally the recurrent network should produce the same output as the
+ * recurrent network. The self connection shouldn't affect the output when using
+ * a sequence with a length of one.
+ */
+BOOST_AUTO_TEST_CASE(FeedForwardRecurrentNetworkTest)
+{
+  arma::mat input;
+  arma::mat labels;
+
+  RandomInitialization<> randInit(1, 1);
+
+  // Test on a non-linearly separable dataset (XOR).
+  input << 0 << 1 << 1 << 0 << arma::endr
+        << 1 << 0 << 1 << 0 << arma::endr;
+  labels << 0 << 0 << 1 << 1;
+
+  // Vanilla neural net with logistic activation function.
+  CompareVanillaNetworks<RandomInitialization<>,
+                      LogisticFunction,
+                      SteepestDescent<>,
+                      BinaryClassificationLayer<>,
+                      MeanSquaredErrorFunction<> >
+      (input, labels, input, labels, 10, 10, randInit);
+
+  // Vanilla neural net with identity activation function.
+  CompareVanillaNetworks<RandomInitialization<>,
+                      IdentityFunction,
+                      SteepestDescent<>,
+                      BinaryClassificationLayer<>,
+                      MeanSquaredErrorFunction<> >
+      (input, labels, input, labels, 1, 1, randInit);
+
+  // Vanilla neural net with rectifier activation function.
+  CompareVanillaNetworks<RandomInitialization<>,
+                    RectifierFunction,
+                    SteepestDescent<>,
+                    BinaryClassificationLayer<>,
+                    MeanSquaredErrorFunction<> >
+    (input, labels, input, labels, 10, 10, randInit);
+
+  // Vanilla neural net with softsign activation function.
+  CompareVanillaNetworks<RandomInitialization<>,
+                    SoftsignFunction,
+                    SteepestDescent<>,
+                    BinaryClassificationLayer<>,
+                    MeanSquaredErrorFunction<> >
+    (input, labels, input, labels, 10, 10, randInit);
+
+  // Vanilla neural net with tanh activation function.
+  CompareVanillaNetworks<RandomInitialization<>,
+                    TanhFunction,
+                    SteepestDescent<>,
+                    BinaryClassificationLayer<>,
+                    MeanSquaredErrorFunction<> >
+    (input, labels, input, labels, 10, 10, randInit);
+}
+
+/**
+ * Generate a random Reber grammar.
+ *
+ * For more information, see the following thesis.
+ *
+ * @code
+ * @misc{Gers2001,
+ *   author = {Felix Gers},
+ *   title = {Long Short-Term Memory in Recurrent Neural Networks},
+ *   year = {2001}
+ * }
+ * @endcode
+ *
+ * @param transitions Reber grammar transition matrix.
+ * @param reber The generated Reber grammar string.
+ */
+void GenerateReber(const arma::Mat<char>& transitions, std::string& reber)
+{
+  size_t idx = 0;
+  reber = "B";
+
+  do
+  {
+    const int grammerIdx = rand() % 2;
+    reber += arma::as_scalar(transitions.submat(idx, grammerIdx, idx,
+        grammerIdx));
+
+    idx = arma::as_scalar(transitions.submat(idx, grammerIdx + 2, idx,
+        grammerIdx + 2)) - '0';
+  } while (idx != 0);
+
+  reber =  "BPTVVE";
+}
+
+/**
+ * Generate a random embedded Reber grammar.
+ *
+ * @param transitions Embedded Reber grammar transition matrix.
+ * @param reber The generated embedded Reber grammar string.
+ */
+void GenerateEmbeddedReber(const arma::Mat<char>& transitions,
+                           std::string& reber)
+{
+  GenerateReber(transitions, reber);
+  const char c = (rand() % 2) == 1 ? 'P' : 'T';
+  reber = c + reber + c;
+  reber = "B" + reber + "E";
+}
+
+/**
+ * Convert a Reber symbol to a unit vector.
+ *
+ * @param symbol Reber symbol to be converted.
+ * @param translation The converted symbol stored as unit vector.
+ */
+void ReberTranslation(const char symbol, arma::colvec& translation)
+{
+  arma::Col<char> symbols;
+  symbols << 'B' << 'T' << 'S' << 'X' << 'P' << 'V' << 'E' << arma::endr;
+  const int idx = arma::as_scalar(arma::find(symbols == symbol, 1, "first"));
+
+  translation = arma::zeros<arma::colvec>(7);
+  translation(idx) = 1;
+}
+
+/**
+ * Convert a unit vector to a Reber symbol.
+ *
+ * @param translation The unit vector to be converted.
+ * @param symbol The converted unit vector stored as Reber symbol.
+ */
+void ReberReverseTranslation(const arma::colvec& translation, char& symbol)
+{
+  arma::Col<char> symbols;
+  symbols << 'B' << 'T' << 'S' << 'X' << 'P' << 'V' << 'E' << arma::endr;
+  const int idx = arma::as_scalar(arma::find(translation == 1, 1, "first"));
+
+  symbol = symbols(idx);
+}
+
+/**
+ * Given a Reber string, return a Reber string with all reachable next symbols.
+ *
+ * @param translation The unit vector to be converted.
+ * @param symbol The converted unit vector stored as Reber symbol.
+ */
+void GenerateNextReber(const arma::Mat<char>& transitions,
+                       const std::string& reber, std::string& nextReber)
+{
+  size_t idx = 0;
+
+  for (size_t grammer = 1; grammer < reber.length(); grammer++)
+  {
+    const int grammerIdx = arma::as_scalar(arma::find(
+        transitions.row(idx) == reber[grammer], 1, "first"));
+
+    idx = arma::as_scalar(transitions.submat(idx, grammerIdx + 2, idx,
+        grammerIdx + 2)) - '0';
+  }
+
+  nextReber = arma::as_scalar(transitions.submat(idx, 0, idx, 0));
+  nextReber += arma::as_scalar(transitions.submat(idx, 1, idx, 1));
+}
+
+BOOST_AUTO_TEST_CASE(ReberGrammarTest)
+{
+  // Reber state transition matrix. (The last two columns are the indices to the
+  // next path).
+  arma::Mat<char> transitions;
+  transitions << 'T' << 'P' << '1' << '2' << arma::endr
+              << 'X' << 'S' << '3' << '1' << arma::endr
+              << 'V' << 'T' << '4' << '2' << arma::endr
+              << 'X' << 'S' << '2' << '5' << arma::endr
+              << 'P' << 'V' << '3' << '5' << arma::endr
+              << 'E' << 'E' << '0' << '0' << arma::endr;
+
+  const size_t trainReberGrammarCount = 1000;
+  const size_t testReberGrammarCount = 10;
+
+  std::string trainReber, testReber;
+  arma::field<arma::mat> trainInput(1, trainReberGrammarCount);
+  arma::field<arma::mat> trainLabels(1, trainReberGrammarCount);
+  arma::field<arma::mat> testInput(1, testReberGrammarCount);
+  arma::field<arma::mat> testLabels(1, testReberGrammarCount);
+  arma::colvec translation;
+
+  // Generate the training data.
+  for (size_t i = 0; i < trainReberGrammarCount; i++)
+  {
+    GenerateReber(transitions, trainReber);
+
+    for (size_t j = 0; j < trainReber.length() - 1; j++)
+    {
+      ReberTranslation(trainReber[j], translation);
+      trainInput(0, i) = arma::join_cols(trainInput(0, i), translation);
+
+      ReberTranslation(trainReber[j + 1], translation);
+      trainLabels(0, i) = arma::join_cols(trainLabels(0, i), translation);
+    }
+  }
+
+  // Generate the test data.
+  for (size_t i = 0; i < testReberGrammarCount; i++)
+  {
+    GenerateReber(transitions, testReber);
+
+    for (size_t j = 0; j < testReber.length() - 1; j++)
+    {
+      ReberTranslation(testReber[j], translation);
+      testInput(0, i) = arma::join_cols(testInput(0, i), translation);
+
+      ReberTranslation(testReber[j + 1], translation);
+      testLabels(0, i) = arma::join_cols(testLabels(0, i), translation);
+    }
+  }
+
+  /*
+   * Construct a network with 7 input units, 5 hidden units and 7 output units.
+   * The hidden layer is connected to itself. The network structure looks like:
+   *
+   *  Input         Hidden        Output
+   * Layer(7)      Layer(5)      Layer(7)
+   * +-----+       +-----+       +-----+
+   * |     |       |     |       |     |
+   * |     +------>|     +------>|     |
+   * |     |    ..>|     |       |     |
+   * +-----+    .  +--+--+       +-----+
+   *            .     .
+   *            .     .
+   *            .......
+   */
+  NeuronLayer<LogisticFunction> inputLayer(7);
+  NeuronLayer<LogisticFunction> hiddenLayer0(5);
+  NeuronLayer<LogisticFunction> recurrentLayer0(hiddenLayer0.InputSize());
+  NeuronLayer<LogisticFunction> hiddenLayer1(7);
+  BinaryClassificationLayer<> outputLayer;
+
+  SteepestDescent< > conOptimizer0(inputLayer.InputSize(),
+      hiddenLayer0.InputSize());
+  SteepestDescent< > conOptimizer2(hiddenLayer0.InputSize(),
+      hiddenLayer0.InputSize());
+  SteepestDescent< > conOptimizer3(hiddenLayer0.InputSize(),
+      hiddenLayer1.OutputSize());
+
+  NguyenWidrowInitialization<> randInit;
+
+  FullConnection<
+      decltype(inputLayer),
+      decltype(hiddenLayer0),
+      decltype(conOptimizer0),
+      decltype(randInit)>
+      layerCon0(inputLayer, hiddenLayer0, conOptimizer0, randInit);
+
+  SelfConnection<
+    decltype(recurrentLayer0),
+    decltype(hiddenLayer0),
+    decltype(conOptimizer2),
+    decltype(randInit)>
+    layerCon2(recurrentLayer0, hiddenLayer0, conOptimizer2, randInit);
+
+  FullConnection<
+      decltype(hiddenLayer0),
+      decltype(hiddenLayer1),
+      decltype(conOptimizer3),
+      decltype(randInit)>
+      layerCon4(hiddenLayer0, hiddenLayer1, conOptimizer3, randInit);
+
+  auto module0 = std::tie(layerCon0, layerCon2);
+  auto module1 = std::tie(layerCon4);
+  auto modules = std::tie(module0, module1);
+
+  RNN<decltype(modules),
+      decltype(outputLayer),
+      MeanSquaredErrorFunction<> > net(modules, outputLayer);
+
+  // Train the network for (500 * trainReberGrammarCount) epochs.
+  Trainer<decltype(net)> trainer(net, 1, 1, 0, false);
+
+  arma::mat inputTemp, labelsTemp;
+  for (size_t i = 0; i < 500; i++)
+  {
+    for (size_t j = 0; j < trainReberGrammarCount; j++)
+    {
+      inputTemp = trainInput.at(0, j);
+      labelsTemp = trainLabels.at(0, j);
+      trainer.Train(inputTemp, labelsTemp, inputTemp, labelsTemp);
+    }
+  }
+
+  double error = 0;
+
+  // Ask the network to predict the next Reber grammar in the given sequence.
+  for (size_t i = 0; i < testReberGrammarCount; i++)
+  {
+    arma::colvec output;
+    arma::colvec input = testInput.at(0, i);
+
+    net.Predict(input, output);
+
+    const size_t reberGrammerSize = 7;
+    std::string inputReber = "";
+
+    size_t reberError = 0;
+    for (size_t j = 0; j < (output.n_elem / reberGrammerSize); j++)
+    {
+      if (arma::sum(output.subvec(j * reberGrammerSize, (j + 1) *
+          reberGrammerSize - 1)) != 1) break;
+
+      char predictedSymbol, inputSymbol;
+      std::string reberChoices;
+
+      ReberReverseTranslation(output.subvec(j * reberGrammerSize, (j + 1) *
+          reberGrammerSize - 1), predictedSymbol);
+      ReberReverseTranslation(input.subvec(j * reberGrammerSize, (j + 1) *
+          reberGrammerSize - 1), inputSymbol);
+      inputReber += inputSymbol;
+
+      GenerateNextReber(transitions, inputReber, reberChoices);
+
+      if (reberChoices.find(predictedSymbol) != std::string::npos)
+        reberError++;
+    }
+
+    if (reberError != (output.n_elem / reberGrammerSize))
+      error += 1;
+  }
+
+  error /= testReberGrammarCount;
+
+  BOOST_REQUIRE_LE(error, 0.2);
+}
+
+BOOST_AUTO_TEST_SUITE_END();