FFNet.h

Go to the documentation of this file.

/*!
 * 
 *
 * \brief       Implements a Feef-Forward multilayer perceptron
 * 
 * 
 *
 * \author      O. Krause
 * \date        2010-2014
 *
 *
 * \par Copyright 1995-2017 Shark Development Team
 * 
 * <BR><HR>
 * This file is part of Shark.
 * <http://shark-ml.org/>
 * 
 * Shark is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published 
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * Shark is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with Shark.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#ifndef SHARK_MODELS_FFNET_H
#define SHARK_MODELS_FFNET_H

#include <shark/Models/AbstractModel.h>
#include <shark/Models/Neurons.h>
#include <boost/serialization/vector.hpp>

namespace shark{
    
struct FFNetStructures{
    enum ConnectionType{
        Normal, //< Layerwise connectivity without shortcuts
        InputOutputShortcut, //< Normal with additional shortcuts from input to output neuron
        Full //< Every layer is fully connected to all neurons in the lower layer 
    };
};

//! \brief Offers the functions to create and to work with a feed-forward network.
//!
//! A feed forward network consists of several layers. every layer consists of a linear 
//! function with optional bias whose response is modified by a (nonlinear) activation function.
//! starting from the input layer, the output of every layer is the input of the next.
//! The two template arguments govern the activation functions of the network. 
//! The activation functions are typically sigmoidal.
//! All hidden layers share one activation function, while the output layer can be chosen to use
//! a different one, for example to allow the last output to be unbounded, in which case a 
//! linear output function is used.
//! It is not possible to use arbitrary activation functions but Neurons following in the structure
//! in Models/Neurons.h Especially it holds that the derivative of the activation function
//! must have the form f'(x) = g(f(x)).
//!
//! This network class allows for several different topologies of structure. The layer-wise structure
//! outlined above is the default one, but the network also allows for shortcuts. most typically 
//! an input-output shortcut is used, that is a shortcut that connects the input neurons directly 
//! with the output using linear weights. But also a fully connected structure is possible, where
//! every layer is fed as input to every successive layer instead of only the next one.
template<class HiddenNeuron,class OutputNeuron>
class FFNet :public AbstractModel<RealVector,RealVector>
{
    struct InternalState: public State{
        //!  \brief Used to store the current results of the activation
        //!         function for all neurons for the last batch of patterns \f$x\f$.
        //!
        //!  There is one value for input+hidden+output units for every element of the batch.
        //!  For every value, the following holds:
        //!  Given a network with \f$M\f$ neurons, including
        //! \f$c\f$ input and \f$n\f$ output neurons the single
        //! values for \f$z\f$ are given as:
        //! <ul>
        //!     <li>\f$z_i = x_i,\ \mbox{for\ } 0 \leq i < c\f$</li>
        //!     <li>\f$z_i = g_{hidden}(x),\ \mbox{for\ } c \leq i < M - n\f$</li>
        //!     <li>\f$z_i = y_{i-M+n} = g_{output}(x),\ \mbox{for\ } M - n \leq
        //!                  i < M\f$</li>
        //! </ul>
        RealMatrix responses;
        
        void resize(std::size_t neurons, std::size_t patterns){
            responses.resize(neurons,patterns);
        }
    };
    

public:
    
    //! Creates an empty feed-forward network. After the constructor is called,
    //! one version of the #setStructure methods needs to be called
    //! to define the network topology.
    FFNet()
    :m_numberOfNeurons(0),m_inputNeurons(0),m_outputNeurons(0){
        m_features|=HAS_FIRST_PARAMETER_DERIVATIVE;
        m_features|=HAS_FIRST_INPUT_DERIVATIVE;
    }

    //! \brief From INameable: return the class name.
    std::string name() const
    { return "FFNet"; }

    //! \brief Number of input neurons.
    std::size_t inputSize()const{
        return m_inputNeurons;
    }
    //! \brief Number of output neurons.
    std::size_t outputSize()const{
        return m_outputNeurons;
    }
    //! \brief Total number of neurons, that is inputs+hidden+outputs.
    std::size_t numberOfNeurons()const{
        return m_numberOfNeurons;
    }
    //! \brief Total number of hidden neurons.
    std::size_t numberOfHiddenNeurons()const{
        return numberOfNeurons() - inputSize() -outputSize();
    }

    //! \brief Returns the matrices for every layer used by eval.
    std::vector<RealMatrix> const& layerMatrices()const{
        return m_layerMatrix;
    }
    
    //! \brief Returns the weight matrix of the i-th layer.
    RealMatrix const& layerMatrix(std::size_t layer)const{
        return m_layerMatrix[layer];
    }
    
    void setLayer(std::size_t layerNumber, RealMatrix const& m, RealVector const& bias){
        SIZE_CHECK(m.size1() == bias.size());
        SIZE_CHECK(m.size1() == m_layerMatrix[layerNumber].size1());
        SIZE_CHECK(m.size2() == m_layerMatrix[layerNumber].size2());
        m_layerMatrix[layerNumber] = m;
        std::size_t start = 0;
        for(std::size_t i = 0; i != layerNumber; ++i){
            start += m_layerMatrix[i].size1();
        }
        noalias(subrange(m_bias,start,start+bias.size())) = bias;
        //set backprop matrices
        setParameterVector(parameterVector());
    }

    //! \brief Returns the matrices for every layer used by backpropagation.
    std::vector<RealMatrix> const& backpropMatrices()const{
        return m_backpropMatrix;
    }
    
    //! \brief Returns the direct shortcuts between input and output neurons.
    //!
    //! This does not necessarily exist.
    RealMatrix const& inputOutputShortcut() const{
        return m_inputOutputShortcut;
    }
    
    /// \brief Returns the activation function of the hidden units.
    HiddenNeuron const& hiddenActivationFunction()const{
        return m_hiddenNeuron;
    }
    /// \brief Returns the activation function of the output units.
    OutputNeuron const& outputActivationFunction()const{
        return m_outputNeuron;
    }
    
    /// \brief Returns the activation function of the hidden units.
    HiddenNeuron& hiddenActivationFunction(){
        return m_hiddenNeuron;
    }
    /// \brief Returns the activation function of the output units.
    OutputNeuron& outputActivationFunction(){
        return m_outputNeuron;
    }

    //! \brief Returns the bias values for hidden and output units.
    //!
    //! This is either empty or a vector of size numberOfNeurons()-inputSize().
    //! the first entry is the value of the first hidden unit while the last outputSize() units
    //! are the values of the output units.
    const RealVector& bias()const{
        return m_bias;
    }
    
    ///\brief Returns the portion of the bias vector of the i-th layer.
    RealVector bias(std::size_t layer)const{
        std::size_t start = 0;
        for(std::size_t i = 0; i != layer; ++i){
            start +=layerMatrices()[i].size1();
        }
        return subrange(m_bias,start,start+layerMatrices()[layer].size1());
    }
    
    //! \brief Returns the total number of parameters of the network. 
    std::size_t numberOfParameters()const{
        std::size_t numParams = m_inputOutputShortcut.size1()*m_inputOutputShortcut.size2();
        numParams += bias().size();
        for(std::size_t i = 0; i != layerMatrices().size(); ++i){
            numParams += layerMatrices()[i].size1()*layerMatrices()[i].size2();
        }
        return numParams;
    }

    //! returns the vector of used parameters inside the weight matrix
    RealVector parameterVector() const{
        RealVector parameters(numberOfParameters());
        std::size_t pos = 0;
        for(auto const& mat: m_layerMatrix){
            auto vec = to_vector(mat);
            noalias(subrange(parameters,pos,pos+vec.size())) = vec;
            pos += vec.size();
        }
        noalias(subrange(parameters,pos,parameters.size())) = m_bias | to_vector(m_inputOutputShortcut);
        return parameters;
    }
    //! uses the values inside the parameter vector to set the used values inside the weight matrix
    void setParameterVector(RealVector const& newParameters){
        //set the forward propagation weights
        std::size_t pos = 0;
        for(auto& mat: m_layerMatrix){
            auto vec = to_vector(mat);
            
            noalias(vec) = subrange(newParameters,pos,pos+vec.size());
            pos += vec.size();
        }
        noalias(m_bias) = subrange(newParameters,pos, pos + m_bias.size());
        noalias(to_vector(m_inputOutputShortcut)) = subrange(newParameters,pos + m_bias.size(), newParameters.size());
        
        //we also have to update the backpropagation weights
        //this is more or less an inversion. for all connections of a neuron i with a neuron j, i->j
        //the backpropagation matrix has an entry j->i.
        
        // we start with all neurons in layer i, looking at all layers j > i and checking whether 
        // they are connected, in this case we transpose the part of the matrix which is connecting
        // layer j with layer i and copying it into the backprop matrix.
        // we assume here, that either all neurons in layer j are connected to all neurons in layer i
        // or that there are no connections at all between the layers.
        std::size_t layeriStart = 0;
        for(std::size_t layeri = 0; layeri != m_layerMatrix.size(); ++layeri){
            std::size_t columni = 0;
            std::size_t neuronsi = inputSize();
            if(layeri > 0)
                neuronsi = m_layerMatrix[layeri-1].size1();
            
            std::size_t layerjStart = layeriStart + neuronsi;
            for(std::size_t layerj = layeri; layerj != m_layerMatrix.size(); ++layerj){
                std::size_t neuronsj = m_layerMatrix[layerj].size1();
                //only process, if layer j has connections with layer i
                if(layerjStart-m_layerMatrix[layerj].size2() <= layeriStart){
                    
                    //Start of the weight columns to layer i in layer j.
                    //parantheses are important to protect against underflow
                    std::size_t weightStartj = layeriStart -(layerjStart - m_layerMatrix[layerj].size2());
                    noalias(columns(m_backpropMatrix[layeri],columni,columni+neuronsj)) 
                    = trans(columns(m_layerMatrix[layerj],weightStartj,weightStartj+neuronsi)); 
                }
                columni += neuronsj;
                layerjStart += neuronsj; 
            }
            layeriStart += neuronsi;
        }
    }

    //! \brief Returns the output of all neurons after the last call of eval
    //!
    //!     \param  state last result of eval
    //!     \return Output value of the neurons.
    RealMatrix const& neuronResponses(State const& state)const{
        InternalState const& s = state.toState<InternalState>();
        return s.responses;
    }
    
    boost::shared_ptr<State> createState()const{
        return boost::shared_ptr<State>(new InternalState());
    }

    ///\brief Returns the response of the i-th layer given the input of that layer.
    ///
    /// this is useful if only a portion of the network needs to be evaluated
    /// be aware that this only works without shortcuts in the network
    void evalLayer(std::size_t layer,RealMatrix const& patterns,RealMatrix& outputs)const{
        std::size_t numPatterns = patterns.size1();
        std::size_t numOutputs = m_layerMatrix[layer].size1();
        outputs.resize(numPatterns,numOutputs);
        outputs.clear();

        //calculate activation. first compute the linear part and the optional bias and then apply
        // the non-linearity
        noalias(outputs) = prod(patterns,trans(layerMatrix(layer)));
        if(!bias().empty()){
            noalias(outputs) += repeat(bias(layer),numPatterns);
        }
        // if this is the last layer, use output neuron response
        if(layer < m_layerMatrix.size()-1) {
            noalias(outputs) = m_hiddenNeuron(outputs);
        }
        else {
            noalias(outputs) = m_outputNeuron(outputs);
        }
    }
    
    ///\brief Returns the response of the i-th layer given the input of that layer.
    ///
    /// this is useful if only a portion of the network needs to be evaluated
    /// be aware that this only works without shortcuts in the network
    Data<RealVector> evalLayer(std::size_t layer, Data<RealVector> const& patterns)const{
        int batches = (int) patterns.numberOfBatches();
        Data<RealVector> result(batches);
        SHARK_PARALLEL_FOR(int i = 0; i < batches; ++i){
            evalLayer(layer,patterns.batch(i),result.batch(i));
        }
        return result;
    }
    
    void eval(RealMatrix const& patterns,RealMatrix& output, State& state)const{
        InternalState& s = state.toState<InternalState>();
        std::size_t numPatterns = patterns.size1();
        //initialize the input layer using the patterns.
        s.resize(numberOfNeurons(),numPatterns);
        s.responses.clear();
        noalias(rows(s.responses,0,m_inputNeurons)) = trans(patterns);
        std::size_t beginNeuron = m_inputNeurons;
        
        for(std::size_t layer = 0; layer != m_layerMatrix.size();++layer){
            RealMatrix const& weights = m_layerMatrix[layer];
            //number of rows of the layer is also the number of neurons
            std::size_t endNeuron = beginNeuron + weights.size1();
            //some subranges of vectors
            //inputs are the last n neurons, where n is the number of columns of the matrix
            auto const input = rows(s.responses,beginNeuron - weights.size2(),beginNeuron);
            //the neurons responses
            auto responses = rows(s.responses,beginNeuron,endNeuron);

            //calculate activation. first compute the linear part and the optional bias and then apply
            // the non-linearity
            noalias(responses) = prod(weights,input);
            if(!bias().empty()){
                //the bias of the layer is shifted as input units can not have bias.
                auto bias = subrange(m_bias,beginNeuron-inputSize(),endNeuron-inputSize());
                noalias(responses) += trans(repeat(bias,numPatterns));
            }
            SHARK_CRITICAL_REGION{//beware Dropout Neurons!
                // if this is the last layer, use output neuron response instead
                if(layer < m_layerMatrix.size()-1) {
                    noalias(responses) = m_hiddenNeuron(responses);
                }
                else {
                    //add shortcuts if necessary
                    if(m_inputOutputShortcut.size1() != 0){
                        noalias(responses) += prod(m_inputOutputShortcut,trans(patterns));
                    }
                    noalias(responses) = m_outputNeuron(responses);
                }
            }
            //go to the next layer
            beginNeuron = endNeuron;
        }
        //Sanity check
        SIZE_CHECK(beginNeuron == m_numberOfNeurons);

        //copy output layer into output
        output.resize(numPatterns,m_outputNeurons);
        noalias(output) = trans(rows(s.responses,m_numberOfNeurons-outputSize(),m_numberOfNeurons));
    }
    using AbstractModel<RealVector,RealVector>::eval;

    void weightedParameterDerivative(
        BatchInputType const& patterns, RealMatrix const& coefficients, State const& state, RealVector& gradient
    )const{
        SIZE_CHECK(coefficients.size2() == m_outputNeurons);
        SIZE_CHECK(coefficients.size1() == patterns.size1());
        std::size_t numPatterns=patterns.size1();
        
        //initialize delta using coefficients and clear the rest. also don't compute the delta for
        // the input neurons as they are not needed.
        RealMatrix delta(numberOfNeurons(),numPatterns,0.0);
        auto outputDelta = rows(delta,delta.size1()-outputSize(),delta.size1());
        noalias(outputDelta) = trans(coefficients);

        computeDelta(delta,state,false);
        computeParameterDerivative(delta,state,gradient);
        
    }
    
    void weightedInputDerivative(
        BatchInputType const& patterns, RealMatrix const& coefficients, State const& state, BatchInputType& inputDerivative
    )const{
        SIZE_CHECK(coefficients.size2() == m_outputNeurons);
        SIZE_CHECK(coefficients.size1() == patterns.size1());
        std::size_t numPatterns=patterns.size1();
        
        //initialize delta using coefficients and clear the rest
        //we compute the full set of delta values here. the delta values of the inputs are the inputDerivative
        RealMatrix delta(numberOfNeurons(),numPatterns,0.0);
        auto outputDelta = rows(delta,delta.size1()-outputSize(),delta.size1());
        noalias(outputDelta) = trans(coefficients);

        computeDelta(delta,state,true);
        inputDerivative.resize(numPatterns,inputSize());
        noalias(inputDerivative) = trans(rows(delta,0,inputSize()));
    }
    
    virtual void weightedDerivatives(
        BatchInputType const & patterns,
        BatchOutputType const & coefficients,
        State const& state,
        RealVector& parameterDerivative,
        BatchInputType& inputDerivative
    )const{
        SIZE_CHECK(coefficients.size2() == m_outputNeurons);
        SIZE_CHECK(coefficients.size1() == patterns.size1());
        std::size_t numPatterns = patterns.size1();
        
        
        //compute full delta and thus the input derivative
        RealMatrix delta(numberOfNeurons(),numPatterns,0.0);
        auto outputDelta = rows(delta,delta.size1()-outputSize(),delta.size1());
        noalias(outputDelta) = trans(coefficients);
        
        computeDelta(delta,state,true);
        inputDerivative.resize(numPatterns,inputSize());
        noalias(inputDerivative) = trans(rows(delta,0,inputSize()));
        
        //reuse delta to compute the parameter derivative
        computeParameterDerivative(delta,state,parameterDerivative);
    }
    
    //! \brief Calculates the derivative for the special case, when error terms for all neurons of the network exist.
    //!
    //! This is useful when the hidden neurons need to meet additional requirements.
    //! The value of delta is changed during computation and holds the results of the backpropagation steps.
    //! The format is such that the rows of delta are the neurons and the columns the patterns.
    void weightedParameterDerivativeFullDelta(
        RealMatrix const& patterns, RealMatrix& delta, State const& state, RealVector& gradient
    )const{
        InternalState const& s = state.toState<InternalState>();
        SIZE_CHECK(delta.size1() == m_numberOfNeurons);
        SIZE_CHECK(delta.size2() == patterns.size1());
        SIZE_CHECK(s.responses.size2() == patterns.size1());
        
        computeDelta(delta,state,false);
        //now compute the parameter derivative from the delta values
        computeParameterDerivative(delta,state,gradient);
    }

    //!  \brief Creates a connection matrix for a network.
    //!
    //! Automatically creates a network with several layers, with
    //! the numbers of neurons for each layer defined by \em layers.
    //! layers must be at least size 2, which will result in a network with no hidden layers.
    //! the first and last values correspond to the number of inputs and outputs respectively.
    //!
    //! The network supports three different tpes of connection models:
    //! FFNetStructures::Normal corresponds to a layerwise connection between consecutive
    //! layers. FFNetStructures::InputOutputShortcut additionally adds a shortcut between
    //! input and output neurons. FFNetStructures::Full connects every layer to every following
    //! layer, this includes also the shortcuts for input and output neurons. Additionally
    //! a bias term an be used.
    //!
    //! While Normal and Full only use the layer matrices, inputOutputShortcut also uses
    //! the corresponding matrix variable (be aware that in the case of only one hidden layer,
    //! the shortcut between input and output leads to the same network as the Full - in that case
    //! the Full topology is chosen for optimization reasons)
    //!
    //! \param  layers contains the numbers of neurons for each layer of the network.
    //! \param  connectivity type of connection used between layers
    //! \param  biasNeuron if set to \em true, connections from
    //!                    all neurons (except the input neurons)
    //!                    to the bias will be set.
    void setStructure(
        std::vector<size_t> const& layers,
        FFNetStructures::ConnectionType connectivity = FFNetStructures::Normal,
        bool biasNeuron = true
    ){
        SIZE_CHECK(layers.size() >= 2);
        m_layerMatrix.resize(layers.size()-1);//we don't model the input layer
        m_backpropMatrix.resize(layers.size()-1);//we don't model the output layer
        
        //small optimization for networks with only 3 layers
        //in this case, we don't need an explicit shortcut as we can integrate it into
        //the big matrices
        if(connectivity == FFNetStructures::InputOutputShortcut && layers.size() ==3)
            connectivity = FFNetStructures::Full;
        
        
        m_inputNeurons = layers.front();
        m_outputNeurons = layers.back();
        m_numberOfNeurons = 0;
        for(std::size_t i = 0; i != layers.size(); ++i){
            m_numberOfNeurons += layers[i];
        }
        if(biasNeuron){
            m_bias.resize(m_numberOfNeurons - m_inputNeurons);
        }
        
        if(connectivity == FFNetStructures::Full){
            //connect to all previous layers.
            std::size_t numNeurons = layers[0];
            for(std::size_t i = 0; i != m_layerMatrix.size(); ++i){
                m_layerMatrix[i].resize(layers[i+1],numNeurons);
                m_backpropMatrix[i].resize(layers[i],m_numberOfNeurons-numNeurons);
                numNeurons += layers[i+1];
                
            }
            m_inputOutputShortcut.resize(0,0);
        }else{
            //only connect with the previous layer
            for(std::size_t i = 0; i != m_layerMatrix.size(); ++i){
                m_layerMatrix[i].resize(layers[i+1],layers[i]);
                m_backpropMatrix[i].resize(layers[i],layers[i+1]);
            }
            
            //create a shortcut from input to output when desired
            if(connectivity == FFNetStructures::InputOutputShortcut){
                m_inputOutputShortcut.resize(m_outputNeurons,m_inputNeurons);
            }
        }
    }
    
    //!  \brief Creates a connection matrix for a network with a
    //!         single hidden layer
    //!
    //!  Automatically creates a network with
    //!  three different layers: An input layer with \em in input neurons,
    //!  an output layer with \em out output neurons and one hidden layer
    //!  with \em hidden neurons, respectively.
    //!
    //!  \param  in         number of input neurons.
    //!  \param  hidden    number of neurons of the second hidden layer.
    //!  \param  out        number of output neurons.
    //!  \param  connectivity  Type of connectivity between the layers
    //!  \param  bias       if set to \em true, connections from
    //!                     all neurons (except the input neurons)
    //!                     to the bias will be set.
    void setStructure(
        std::size_t in,
        std::size_t hidden,
        std::size_t out,
        FFNetStructures::ConnectionType connectivity = FFNetStructures::Normal,
        bool bias      = true
    ){
        std::vector<size_t> layer(3);
        layer[0] = in;
        layer[1] = hidden;
        layer[2] = out;
        setStructure(layer, connectivity, bias);
    }
    
    //!  \brief Creates a connection matrix for a network with two
    //!         hidden layers.
    //!
    //!  Automatically creates a network with
    //!  four different layers: An input layer with \em in input neurons,
    //!  an output layer with \em out output neurons and two hidden layers
    //!  with \em hidden1 and \em hidden2 hidden neurons, respectively.
    //!
    //!  \param  in         number of input neurons.
    //!  \param  hidden1    number of neurons of the first hidden layer.
    //!  \param  hidden2    number of neurons of the second hidden layer.
    //!  \param  out        number of output neurons.
    //!  \param  connectivity  Type of connectivity between the layers
    //!  \param  bias       if set to \em true, connections from
    //!                     all neurons (except the input neurons)
    //!                     to the bias will be set.
    void setStructure(
        std::size_t in,
        std::size_t hidden1,
        std::size_t hidden2,
        std::size_t out,
        FFNetStructures::ConnectionType connectivity = FFNetStructures::Normal,
        bool bias      = true
    ){
        std::vector<size_t> layer(4);
        layer[0] = in;
        layer[1] = hidden1;
        layer[2] = hidden2;
        layer[3] = out;
        setStructure(layer, connectivity, bias);
    }

    //! From ISerializable, reads a model from an archive
    void read( InArchive & archive ){
        archive>>m_inputNeurons;
        archive>>m_outputNeurons;
        archive>>m_numberOfNeurons;
        archive>>m_layerMatrix;
        archive>>m_backpropMatrix;
        archive>>m_inputOutputShortcut;
        archive>>m_bias;
    }

    //! From ISerializable, writes a model to an archive
    void write( OutArchive & archive ) const{
        archive<<m_inputNeurons;
        archive<<m_outputNeurons;
        archive<<m_numberOfNeurons;
        archive<<m_layerMatrix;
        archive<<m_backpropMatrix;
        archive<<m_inputOutputShortcut;
        archive<<m_bias;
    }


private:
    
    void computeDelta(
        RealMatrix& delta, State const& state, bool computeInputDelta
    )const{
        SIZE_CHECK(delta.size1() == numberOfNeurons());
        InternalState const& s = state.toState<InternalState>();

        //initialize output neurons using coefficients
        auto outputDelta = rows(delta,delta.size1()-outputSize(),delta.size1());
        auto outputResponse = rows(s.responses,delta.size1()-outputSize(),delta.size1());
        noalias(outputDelta) *= m_outputNeuron.derivative(outputResponse);

        //iterate backwards using the backprop matrix and propagate the errors to get the needed delta values
        //we stop until we have filled all delta values. Thus we might not necessarily compute all layers.
        
        //last neuron of the current layer that we need to compute
        //we don't need (or can not) compute the values of the output neurons as they are given from the outside
        std::size_t endNeuron = delta.size1()-outputSize();
        std::size_t layer = m_backpropMatrix.size()-1;
        std::size_t endIndex = computeInputDelta? 0: inputSize();
        while(endNeuron > endIndex){
            
            RealMatrix const& weights = m_backpropMatrix[layer];
            std::size_t beginNeuron = endNeuron - weights.size1();//first neuron of the current layer
            //get the delta and response values of this layer
            auto layerDelta = rows(delta,beginNeuron,endNeuron);
            auto layerDeltaInput = rows(delta,endNeuron,endNeuron+weights.size2());
            auto layerResponse = rows(s.responses,beginNeuron,endNeuron);

            noalias(layerDelta) += prod(weights,layerDeltaInput);//add the values to the maybe non-empty delta part
            if(layer != 0){
                noalias(layerDelta) *= m_hiddenNeuron.derivative(layerResponse);
            }
            //go a layer backwards
            endNeuron=beginNeuron;
            --layer;
        }
        
        //add the shortcut deltas if necessary
        if(inputOutputShortcut().size1() != 0)
            noalias(rows(delta,0,inputSize())) += prod(trans(inputOutputShortcut()),outputDelta);
    }
    
    void computeParameterDerivative(RealMatrix const& delta, State const& state, RealVector& gradient)const{
        SIZE_CHECK(delta.size1() == numberOfNeurons());
        InternalState const& s = state.toState<InternalState>();
        // calculate error gradient
        //todo: take network structure into account to prevent checking all possible weights...
        gradient.resize(numberOfParameters());
        std::size_t pos = 0;
        std::size_t layerStart = inputSize();
        for(std::size_t layer = 0; layer != layerMatrices().size(); ++layer){
            //obtain input, delta and gradients for the current layer
            std::size_t layerRows =  layerMatrices()[layer].size1();
            std::size_t layerColumns =  layerMatrices()[layer].size2();
            std::size_t params = layerRows*layerColumns;
            auto gradMatrix  = to_matrix(subrange(gradient,pos,pos+params),layerRows,layerColumns);
            auto deltaLayer = rows(delta,layerStart,layerStart+layerRows);
            auto inputLayer = rows(s.responses,layerStart-layerColumns,layerStart);
            noalias(gradMatrix) = prod(deltaLayer, trans(inputLayer));
            
            pos += params;
            layerStart += layerRows;
        }
        //check whether we need the bias derivative
        if(!bias().empty()){
            //calculate bias derivative
            for (std::size_t neuron = m_inputNeurons; neuron < m_numberOfNeurons; neuron++){
                gradient(pos) = sum(row(delta,neuron));
                pos++;
            }
        }
        //compute shortcut derivative
        if(inputOutputShortcut().size1() != 0){
            std::size_t params = inputSize()*outputSize();
            auto gradMatrix  = to_matrix(subrange(gradient,pos,pos+params),outputSize(),inputSize());
            auto deltaLayer = rows(delta,delta.size1()-outputSize(),delta.size1());
            auto inputLayer = rows(s.responses,0,inputSize());
            noalias(gradMatrix) = prod(deltaLayer, trans(inputLayer));
        }
        
    }
    

    //! \brief Number of all network neurons.
    //!
    //! This is the total number of neurons in the network, i.e.
    //! input, hidden and output neurons.
    std::size_t m_numberOfNeurons;
    std::size_t m_inputNeurons;
    std::size_t m_outputNeurons;

    //! \brief represents the connection matrix using a layered structure for forward propagation
    //!
    //! a layer is made of neurons with consecutive indices which are not
    //! connected with each other. In other words, if there exists a k i<k<j such
    //! that C(i,k) = 1 or C(k,j) = 1 or C(j,i) = 1 than the neurons i,j are not in the same layer.
    //! This is the forward view, meaning that the layers holds the weights which are used to calculate
    //! the activation of the neurons of the layer.
    std::vector<RealMatrix> m_layerMatrix;
    
    //! \brief optional matrix directly connecting input to output
    //!
    //! This is only filled when the network has an input-output shortcut but not a full layer connection.
    RealMatrix m_inputOutputShortcut;
    
    //!\brief represents the backwards view of the network as layered structure.
    //!
    //! This is the backward view of the Network which is used for the backpropagation step. So every
    //! Matrix contains the weights of the neurons which are activated by the layer.
    std::vector<RealMatrix> m_backpropMatrix;

    //! bias weights of the neurons
    RealVector m_bias;

    //!Type of hidden neuron. See Models/Neurons.h for a few choices
    HiddenNeuron m_hiddenNeuron;
    //! Type of output neuron. See Models/Neurons.h for a few choices
    OutputNeuron m_outputNeuron;
};


}
#endif