multilayer_perceptron.cpp

/****************************************************************************************************************/
/*                                                                                                              */
/*   OpenNN: Open Neural MultilayerPerceptrons Library                                                          */
/*   www.artelnics.com/opennn                                                                                   */
/*                                                                                                              */
/*   M U L T I L A Y E R   P E R C E P T R O N   C L A S S                                                      */
/*                                                                                                              */
/*   Roberto Lopez                                                                                              */
/*   Artelnics - Making intelligent use of data                                                                 */
/*   [email protected]                                                                                 */
/*                                                                                                              */
/****************************************************************************************************************/

// OpenNN includes

#include "multilayer_perceptron.h"

namespace OpenNN
{

// DEFAULT CONSTRUCTOR


MultilayerPerceptron::MultilayerPerceptron(void)
{
   set();
}


// LAYERS CONSTRUCTOR


MultilayerPerceptron::MultilayerPerceptron(const Vector<PerceptronLayer>& new_layers)
{
   set(new_layers);
}


// NETWORK ARCHITECTURE CONSTRUCTOR


MultilayerPerceptron::MultilayerPerceptron(const Vector<size_t>& new_architecture)
{
   set(new_architecture);
}


// ONE LAYER CONSTRUCTOR


MultilayerPerceptron::MultilayerPerceptron(const size_t& new_inputs_number, const size_t& new_neurons_number)
{
   set(new_inputs_number, new_neurons_number);
}


// TWO LAYERS CONSTRUCTOR


MultilayerPerceptron::MultilayerPerceptron(const size_t& new_inputs_number, const size_t& new_hidden_neurons_number, const size_t& new_outputs_number)
{
   set(new_inputs_number, new_hidden_neurons_number, new_outputs_number);

   set_default();
}


// COPY CONSTRUCTOR


MultilayerPerceptron::MultilayerPerceptron(const MultilayerPerceptron& other_multilayer_perceptron)
{
   set(other_multilayer_perceptron);
}


// DESTRUCTOR


MultilayerPerceptron::~MultilayerPerceptron(void)
{
}


// ASSIGNMENT OPERATOR


MultilayerPerceptron& MultilayerPerceptron::operator = (const MultilayerPerceptron& other_multilayer_perceptron)
{
   if(this != &other_multilayer_perceptron) 
   {
      layers = other_multilayer_perceptron.layers;
   
      display = other_multilayer_perceptron.display;
   }

   return(*this);
}


// EQUAL TO OPERATOR

// bool operator == (const MultilayerPerceptron&) const method


bool MultilayerPerceptron::operator == (const MultilayerPerceptron& other_multilayer_perceptron) const
{
   if(layers == other_multilayer_perceptron.layers   
   && display == other_multilayer_perceptron.display)
   {
      return(true);   
   }
   else
   {
      return(false);
   }
}


// METHODS

// const Vector<PerceptronLayer>& get_layers(void) const method


const Vector<PerceptronLayer>& MultilayerPerceptron::get_layers(void) const 
{
   return(layers);
}


// const PerceptronLayer& get_layer(const size_t&) const method


const PerceptronLayer& MultilayerPerceptron::get_layer(const size_t& i) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t layers_number = get_layers_number();

   if(i >= layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "const PerceptronLayer get_layer(const size_t&) const method.\n"
             << "Index of layer must be less than number of layers.\n";

      throw std::logic_error(buffer.str());
   }

   #endif
 
   return(layers[i]);
}


// PerceptronLayer* get_layer_pointer(const size_t&) const method


PerceptronLayer* MultilayerPerceptron::get_layer_pointer(const size_t& i)
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t layers_number = get_layers_number();

   if(i >= layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "PerceptronLayer* get_layer_pointer(const size_t&) const method.\n"
             << "Index of layer must be less than number of layers.\n";

      throw std::logic_error(buffer.str());
   }

   #endif
 
   return(&layers[i]);
}


// size_t count_perceptrons_number(void) const method


size_t MultilayerPerceptron::count_perceptrons_number(void) const
{
   const Vector<size_t> layers_perceptrons_number = arrange_layers_perceptrons_numbers();

   return(layers_perceptrons_number.calculate_sum());
}


// Vector<size_t> count_cumulative_perceptrons_number(void) const method


Vector<size_t> MultilayerPerceptron::count_cumulative_perceptrons_number(void) const
{
   const size_t layers_number = get_layers_number();

   Vector<size_t> cumulative_neurons_number(layers_number);

   if(layers_number != 0)
   {
      const Vector<size_t> layers_size = arrange_layers_perceptrons_numbers();

      cumulative_neurons_number[0] = layers_size[0];

      for(size_t i = 1; i < layers_number; i++)
      {
         cumulative_neurons_number[i] = cumulative_neurons_number[i-1] + layers_size[i];   
      }
   }  

   return(cumulative_neurons_number);
}


// Vector<size_t> arrange_layers_parameters_number(void) const method


Vector<size_t> MultilayerPerceptron::arrange_layers_parameters_number(void) const
{
   const size_t layers_number = get_layers_number();

   Vector<size_t> layers_parameters_number(layers_number);

   for(size_t i = 0; i < layers_number; i++)
   {
      layers_parameters_number[i] = layers[i].count_parameters_number();
   }

   return(layers_parameters_number);
}


// Vector<size_t> arrange_layers_cumulative_parameters_number(void) const method


Vector<size_t> MultilayerPerceptron::arrange_layers_cumulative_parameters_number(void) const
{
   const size_t layers_number = get_layers_number();

   Vector<size_t> layers_cumulative_parameters_number(layers_number);

   layers_cumulative_parameters_number[0] = layers[0].count_parameters_number();

   for(size_t i = 1; i < layers_number; i++)
   {
      layers_cumulative_parameters_number[i] = layers_cumulative_parameters_number[i-1] + layers[i].count_parameters_number();   
   }

   return(layers_cumulative_parameters_number);
}


// Vector< Vector<double> > arrange_layers_biases(void) const method


Vector< Vector<double> > MultilayerPerceptron::arrange_layers_biases(void) const
{
   const size_t layers_number = get_layers_number();

   Vector< Vector<double> > layers_biases(layers_number);

   for(size_t i = 0; i < layers_number; i++)
   {
      layers_biases[i] = layers[i].arrange_biases();
   }

   return(layers_biases);
}


// Vector< Matrix<double> > arrange_layers_synaptic_weights(void) const method


Vector< Matrix<double> > MultilayerPerceptron::arrange_layers_synaptic_weights(void) const
{
   const size_t layers_number = get_layers_number();

   Vector< Matrix<double> > layers_synaptic_weights(layers_number);

   for(size_t i = 0; i < layers_number; i++)
   {
      layers_synaptic_weights[i] = layers[i].arrange_synaptic_weights();   
   }

   return(layers_synaptic_weights);
}


// Vector< Matrix<double> > get_layers_parameters(void) const method


Vector< Vector<double> > MultilayerPerceptron::get_layers_parameters(void) const
{
   const size_t layers_number = get_layers_number();

   Vector< Vector<double> > layers_parameters(layers_number);

   for(size_t i = 0; i < layers_number; i++)
   {
      layers_parameters[i] = layers[i].arrange_parameters();   
   }

   return(layers_parameters);
}


// size_t count_parameters_number(void) const method


size_t MultilayerPerceptron::count_parameters_number(void) const
{
   const size_t layers_number = get_layers_number();

   size_t parameters_number = 0;

   for(size_t i = 0; i < layers_number; i++)
   {
      parameters_number += layers[i].count_parameters_number();
   }

   return(parameters_number);
}


// Vector<double> arrange_parameters(void) const method


Vector<double> MultilayerPerceptron::arrange_parameters(void) const
{
   const size_t layers_number = get_layers_number();

   const size_t parameters_number = count_parameters_number();

   Vector<double> parameters(parameters_number);

   size_t position = 0;

   for(size_t i = 0; i < layers_number; i++)
   {
      const Vector<double> layer_parameters = layers[i].arrange_parameters();
      const size_t layer_parameters_number = layer_parameters.size();

      parameters.tuck_in(position, layer_parameters);
      position += layer_parameters_number; 
   }

   return(parameters);
}


// Vector<size_t> count_layers_parameters_numbers(void) const method


Vector<size_t> MultilayerPerceptron::count_layers_parameters_numbers(void) const
{
    const size_t layers_number = get_layers_number();

    Vector<size_t> layers_parameters_numbers(layers_number);

    for(size_t i = 0; i < layers_number; i++)
    {
        layers_parameters_numbers[i] = layers[i].count_parameters_number();
    }

    return(layers_parameters_numbers);
}


// size_t get_layer_index(const size_t&) const method


size_t MultilayerPerceptron::get_layer_index(const size_t& neuron_index) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t neurons_number = count_perceptrons_number();

   if(neuron_index >= neurons_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "int get_layer_index(const size_t&) const method.\n"
             << "Index of neuron must be less than number of neurons.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const Vector<size_t> cumulative_neurons_number = count_cumulative_perceptrons_number();

   const size_t layer_index = cumulative_neurons_number.calculate_cumulative_index(neuron_index);

   return(layer_index);
}


// size_t get_perceptron_index(const size_t&, const size_t&) const method


size_t MultilayerPerceptron::get_perceptron_index(const size_t& layer_index, const size_t& perceptron_position) const
{
   if(layer_index == 0)
   {
      return(perceptron_position);
   }
   else
   {
      const Vector<size_t> cumulative_neurons_number = count_cumulative_perceptrons_number();

      return(cumulative_neurons_number[layer_index-1] + perceptron_position);
   }
}


// size_t get_layer_bias_index(const size_t&, const size_t&) const method


size_t MultilayerPerceptron::get_layer_bias_index(const size_t& layer_index, const size_t& perceptron_index) const
{  
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t layers_number = get_layers_number();

   if(layer_index >= layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "size_t get_layer_bias_index(const size_t&, const size_t&) const method.\n"
             << "Index of layer must be less than number of layers.\n";

      throw std::logic_error(buffer.str());
   }

   const size_t layer_perceptrons_number = layers[layer_index].get_perceptrons_number();

   if(perceptron_index >= layer_perceptrons_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "size_t get_layer_bias_index(const size_t&, const size_t&) const method.\n"
             << "Index of perceptron must be less than number of perceptrons in that layer.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   size_t layer_bias_index = 0;

   // Previous layers

   for(size_t i = 0; i < layer_index; i++)
   {
      layer_bias_index += layers[i].count_parameters_number();
   }

   // Previous layer neurons

   for(size_t j = 0; j < perceptron_index; j++)
   {
      layer_bias_index += layers[layer_index].get_perceptron(j).count_parameters_number();
   }

   return(layer_bias_index);
}


// size_t get_layer_synaptic_weight_index(const size_t&, const size_t&, const size_t&) const method


size_t MultilayerPerceptron::get_layer_synaptic_weight_index
(const size_t& layer_index, const size_t& perceptron_index, const size_t& input_index) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t layers_number = get_layers_number();

   if(layer_index >= layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "size_t get_layer_synaptic_weight_index(const size_t&, const size_t&, const size_t&) method.\n"
             << "Index of layer must be less than number of layers.\n";

      throw std::logic_error(buffer.str());
   }

   const size_t layer_perceptrons_number = layers[layer_index].get_perceptrons_number();

   if(perceptron_index >= layer_perceptrons_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "size_t get_layer_synaptic_weight_index(const size_t&, const size_t&, const size_t&) method.\n"
             << "Index of perceptron must be less than number of perceptrons in layer.\n";

      throw std::logic_error(buffer.str());
   }

   const size_t layer_inputs_number = layers[layer_index].get_inputs_number();

   if(input_index >= layer_inputs_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "size_t get_layer_synaptic_weight_index(const size_t&, const size_t&, const size_t&) method.\n"
             << "Index of inputs must be less than number of inputs in perceptron.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   size_t layer_synaptic_weight_index = 0;

   // Previous layers

   if(layer_index > 0)
   {
      for(size_t i = 0; i < layer_index-1; i++)
      {
         layer_synaptic_weight_index += layers[layer_index].count_parameters_number();
      }
   }

   // Previous layer neurons

   if(perceptron_index > 0)
   {
      for(size_t i = 0; i < perceptron_index-1; i++)
      {
         layer_synaptic_weight_index += layers[layer_index].get_perceptron(i).count_parameters_number();
      }
   }

   // Hidden neuron bias

   layer_synaptic_weight_index += 1;

   // Synaptic weight index

   layer_synaptic_weight_index += input_index;

   return(layer_synaptic_weight_index);
}


// Vector<size_t> arrange_parameter_indices(const size_t&) const method


Vector<size_t> MultilayerPerceptron::arrange_parameter_indices(const size_t& parameter_index) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t parameters_number = count_parameters_number();

   if(parameter_index >= parameters_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector<int> arrange_parameter_indices(const size_t&) const method.\n"
             << "Index of neural parameter must be less than number of multilayer perceptron parameters.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   return(arrange_parameters_indices().arrange_row(parameter_index));
}


// Matrix<size_t> arrange_parameters_indices(void) const method


Matrix<size_t> MultilayerPerceptron::arrange_parameters_indices(void) const
{
   size_t perceptron_parameters_number;

   const size_t layers_number = get_layers_number();
   const Vector<size_t> layers_size = arrange_layers_perceptrons_numbers();

   const size_t parameters_number = count_parameters_number();

   Matrix<size_t> parameters_indices(parameters_number, 3);
   
   size_t parameter_index = 0;

   for(size_t i = 0; i < layers_number; i++)
   {
      for(size_t j = 0; j < layers_size[i]; j++)
      {
         perceptron_parameters_number = layers[i].get_perceptron(j).count_parameters_number();

         for(size_t k = 0; k < perceptron_parameters_number; k++)
         {
            parameters_indices(parameter_index,0) = i;
            parameters_indices(parameter_index,1) = j;
            parameters_indices(parameter_index,2) = k;
            parameter_index++;
         }
      }
   }

   return(parameters_indices);
}


// Vector<Perceptron::ActivationFunction> get_layers_activation_function(void) const method


Vector<Perceptron::ActivationFunction> MultilayerPerceptron::get_layers_activation_function(void) const
{
   const size_t layers_number = get_layers_number();

   Vector<Perceptron::ActivationFunction> layers_activation_function(layers_number);

   for(size_t i = 0; i < layers_number; i++)
   {
      layers_activation_function[i] = layers[i].get_activation_function();
   }

   return(layers_activation_function);
}
   

// Vector<std::string> write_layers_activation_function(void) const method


Vector<std::string> MultilayerPerceptron::write_layers_activation_function(void) const
{
   const size_t layers_number = get_layers_number();

   Vector<std::string> layers_activation_function(layers_number);

   for(size_t i = 0; i < layers_number; i++)
   {
      layers_activation_function[i] = layers[i].write_activation_function();
   }

   return(layers_activation_function);
}


// const bool& get_display(void) const method


const bool& MultilayerPerceptron::get_display(void) const
{
   return(display);
}


// void set_default(void) method


void MultilayerPerceptron::set_default(void)
{
   // Multilayer perceptron architecture

   const size_t layers_number = get_layers_number();

   if(layers_number > 0)
   { 
      for(size_t i = 0; i < layers_number-1; i++)
      {
         layers[i].set_activation_function(Perceptron::HyperbolicTangent);         
      }

      layers[layers_number-1].set_activation_function(Perceptron::Linear);
   }
 
   // Display messages
   
   set_display(true);
}


// void set(void) method


void MultilayerPerceptron::set(void)
{
   layers.set();    
}


// void set(const Vector<PerceptronLayer>&) method


void MultilayerPerceptron::set(const Vector<PerceptronLayer>& new_layers)
{
   layers = new_layers;   
}


// void set(const Vector<size_t>&) method


void MultilayerPerceptron::set(const Vector<size_t>& new_architecture)
{
   std::ostringstream buffer;

   const size_t new_architecture_size = new_architecture.size();

   if(new_architecture_size == 0)
   {
      set();
   }
   else if(new_architecture_size == 1)
   {
      buffer << "OpenNN Exception: MultilayerPerceptron class.\n" 
             << "void set_architecture(const Vector<size_t>&) method.\n"
             << "Size of architecture cannot be one.\n";

      throw std::logic_error(buffer.str());   
   }
   else
   {
      for(size_t i = 0; i < new_architecture_size; i++)
      {
         if(new_architecture[i] == 0)
         {
            buffer << "OpenNN Exception: MultilayerPerceptron class.\n" 
                   << "void set_architecture(const Vector<size_t>&) method.\n"
                   << "Size " << i << " must be greater than zero.\n";

            throw std::logic_error(buffer.str());
         }
      }

      const size_t new_layers_number = new_architecture_size-1;
      layers.set(new_layers_number);
      
     // First layer

      for(size_t i = 0; i < new_layers_number; i++)
      {
         layers[i].set(new_architecture[i], new_architecture[i+1]);
      }
   
      // Activation

      for(size_t i = 0; i < new_layers_number-1; i++)
      {
         layers[i].set_activation_function(Perceptron::HyperbolicTangent);
      }

      layers[new_layers_number-1].set_activation_function(Perceptron::Linear);
   }
}   


// void set(const size_t&, const size_t&) method


void MultilayerPerceptron::set(const size_t& new_inputs_number, const size_t& new_perceptrons_number)
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   std::ostringstream buffer;

   if(new_inputs_number == 0)
   {  
      buffer << "OpenNN Exception: MultilayerPerceptron class.\n" 
             << "void set(const size_t&, const size_t&) method.\n"
             << "Number of inputs cannot be zero.\n";

      throw std::logic_error(buffer.str());
   }   
   else if(new_perceptrons_number == 0)
   {  
      buffer << "OpenNN Exception: MultilayerPerceptron class.\n" 
             << "void set_architecture(const size_t&, const size_t&) method.\n"
             << "Number of perceptrons cannot be zero.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   layers.set(1);

   layers[0].set(new_inputs_number, new_perceptrons_number);
}


// void set(const size_t&, const size_t&, const size_t&) method


void MultilayerPerceptron::set(const size_t& new_inputs_number, const size_t& new_hidden_neurons_number, const size_t& new_outputs_number)
{
  // Control sentence (if debug)

   #ifndef NDEBUG 

   std::ostringstream buffer;

   if(new_inputs_number == 0)
   {
      buffer << "OpenNN Exception: MultilayerPerceptron class.\n" 
             << "void set(const size_t&, const size_t&, const size_t&) method.\n"
             << "Number of inputs must be greater than zero.\n";

      throw std::logic_error(buffer.str());
   }
   else if(new_hidden_neurons_number == 0)
   {
      buffer << "OpenNN Exception: MultilayerPerceptron class.\n" 
             << "void set(const size_t&, const size_t&, const size_t&) method.\n"
             << "Number of hidden neurons must be greater than zero.\n";

      throw std::logic_error(buffer.str());
   }
   else if(new_outputs_number == 0)
   {
      buffer << "OpenNN Exception: MultilayerPerceptron class.\n" 
             << "void set(const size_t&, const size_t&, const size_t&) method.\n"
             << "Number of outputs must be greater than zero.\n";
   
      throw std::logic_error(buffer.str());
   }
   
   #endif

   layers.set(2);

   layers[0].set(new_inputs_number, new_hidden_neurons_number);
   layers[0].set_activation_function(Perceptron::HyperbolicTangent);

   layers[1].set(new_hidden_neurons_number, new_outputs_number);
   layers[1].set_activation_function(Perceptron::Linear);
}


// void set(const MultilayerPerceptron&) method


void MultilayerPerceptron::set(const MultilayerPerceptron& other_multilayer_perceptron)
{
   layers = other_multilayer_perceptron.layers;

   display = other_multilayer_perceptron.display;
}


// void set_inputs_number(const size_t&) method


void MultilayerPerceptron::set_inputs_number(const size_t& new_inputs_number)
{
   if(is_empty())
   {
      std::ostringstream buffer;
   
      buffer << "OpenNN Exception: MultilayerPerceptron class.\n" 
             << "void set_inputs_number(const size_t&) method.\n"
             << "Multilayer perceptron is empty.\n";

      throw std::logic_error(buffer.str());
   } 

   layers[0].set_inputs_number(new_inputs_number);
}


// void set_layers_perceptrons_number(const Vector<size_t>&) method


void MultilayerPerceptron::set_layers_perceptrons_number(const Vector<size_t>& new_layers_size)
{
   const Vector<size_t> inputs_number(1, get_inputs_number());

   set(inputs_number.assemble(new_layers_size));
}


// void set_layer_perceptrons_number(const size_t&, const size_t&) method


void MultilayerPerceptron::set_layer_perceptrons_number(const size_t& layer_index, const size_t& new_layer_perceptrons_number)
{
   const size_t layer_inputs_number = layers[layer_index].get_inputs_number();

   layers[layer_index].set_perceptrons_number(new_layer_perceptrons_number);
   
   layers[layer_index].set_inputs_number(layer_inputs_number);

   const size_t layers_number = get_layers_number();

   if(layer_index < layers_number-1)
   {
      layers[layer_index+1].set_inputs_number(new_layer_perceptrons_number);   
   }
}


// void set_layers(const Vector<PerceptronLayer>&) method


void MultilayerPerceptron::set_layers(const Vector<PerceptronLayer>& new_layers)
{
   layers = new_layers;
}


// void set_layers_biases(const Vector< Vector<double> >&) method


void MultilayerPerceptron::set_layers_biases(const Vector< Vector<double> >& new_layers_biases)
{
   const size_t layers_number = get_layers_number();

   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t size = new_layers_biases.size();

   if(size != layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "void set_layers_biases(const Vector< Vector<double> >&) method.\n"
             << "Size (" << size << ") must be equal to number of layers (" << layers_number << ").\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const Vector<size_t> layers_size = arrange_layers_perceptrons_numbers();

   for(size_t i = 0; i < layers_number; i++)
   {
      layers[i].set_biases(new_layers_biases[i]); 
   }
}


// void set_layers_synaptic_weights(const Vector< Matrix<double> >&) method


void MultilayerPerceptron::set_layers_synaptic_weights(const Vector< Matrix<double> >& new_layers_synaptic_weights)
{
   const size_t layers_number = get_layers_number();

   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t size = new_layers_synaptic_weights.size();

   if(size != layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "void set_layers_synaptic_weights(const Vector< Matrix<double> >&) method.\n"
             << "Size must be equal to number of layers.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   for(size_t i = 0; i < layers_number; i++)
   {
      layers[i].set_synaptic_weights(new_layers_synaptic_weights[i]);
   }

}


// void set_layers_parameters(const Vector< Vector<double> >&) method


void MultilayerPerceptron::set_layers_parameters(const Vector< Vector<double> >& new_layers_parameters)
{
   const size_t layers_number = get_layers_number();

   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t new_layers_parameters_size = new_layers_parameters.size();

   if(new_layers_parameters_size != layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "void set_layers_parameters(const Vector< Vector<double> >&) method.\n"
             << "Size of layer parameters must be equal to number of layers.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   for(size_t i = 0; i < layers_number; i++)
   {
      layers[i].set_parameters(new_layers_parameters[i]);
   }
}


// void set_parameters(const Vector<double>&) method


void MultilayerPerceptron::set_parameters(const Vector<double>& new_parameters)
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t size = new_parameters.size();

   const size_t parameters_number = count_parameters_number();

   if(size != parameters_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "void set_parameters(const Vector<double>&) method.\n"
             << "Size (" << size << ") must be equal to number of biases and synaptic weights (" << parameters_number << ").\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const size_t layers_number = get_layers_number();

   size_t layer_parameters_number;
   Vector<double> layer_parameters;

   size_t position = 0;

   for(size_t i = 0; i < layers_number; i++)
   {
      layer_parameters_number = layers[i].count_parameters_number();
      layer_parameters = layers[i].arrange_parameters();

      layer_parameters = new_parameters.take_out(position, layer_parameters_number);

      layers[i].set_parameters(layer_parameters);
      position += layer_parameters_number;
   }
}


// void set_layers_activation_function(const Vector<Perceptron::ActivationFunction>&) method


void MultilayerPerceptron::set_layers_activation_function
(const Vector<Perceptron::ActivationFunction>& new_layers_activation_function)
{ 
   const size_t layers_number = get_layers_number();

   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t new_layers_activation_function_size = new_layers_activation_function.size();

   if(new_layers_activation_function_size != layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n" 
             << "void set_layers_activation_function(const Vector<Perceptron::ActivationFunction>&) method.\n"
             << "Size of activation function of layers must be equal to number of layers.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   for(size_t i = 0; i < layers_number; i++)
   {
      layers[i].set_activation_function(new_layers_activation_function[i]);
   }
}


// void set_layer_activation_function(const size_t&, const Perceptron::ActivationFunction&) method


void MultilayerPerceptron::set_layer_activation_function
(const size_t& i, const Perceptron::ActivationFunction& new_layer_activation_function)
{
   // Control sentence (if debug)

   #ifndef NDEBUG

    const size_t layers_number = get_layers_number();

   if(i >= layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "void set_layer_activation_function(const size_t&, const Perceptron::ActivationFunction&) method.\n"
             << "Index of layer is equal or greater than number of layers.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

    layers[i].set_activation_function(new_layer_activation_function);
}


// void set_layers_activation_function(const Vector<std::string>&) method


void MultilayerPerceptron::set_layers_activation_function(const Vector<std::string>& new_layers_activation_function)
{
   const size_t layers_number = get_layers_number();

   // Control sentence (if debug)

   #ifndef NDEBUG

   const size_t size = new_layers_activation_function.size();

   if(size != layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "void set_layer_activation_function(const Vector<std::string>&) method.\n"
             << "Size of layers activation function is not equal to number of layers.\n";

      throw std::logic_error(buffer.str());
   }

   #endif


   for(size_t i = 0; i < layers_number; i++)
   {
      layers[i].set_activation_function(new_layers_activation_function[i]);
   }
}


// void set_display(const bool&) method


void MultilayerPerceptron::set_display(const bool& new_display)
{
   display = new_display;
}


// bool is_empty(void) const method


bool MultilayerPerceptron::is_empty(void) const
{
   if(layers.empty())
   {
      return(true);
   }
   else
   {
      return(false);
   }
}


// void grow_input(void) method


void MultilayerPerceptron::grow_input(void)
{
    const size_t layers_number = get_layers_number();

    if(layers_number > 0)
    {
//        layers[0].grow_input();
    }
}


// void grow_layer(const size_t&, const size_t&) method


void MultilayerPerceptron::grow_layer(const size_t& layer_index, const size_t& perceptrons_number)
{
    #ifndef NDEBUG

    if(layer_index >= (layers.size()-1))
    {
       std::ostringstream buffer;

       buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
              << "void grow_layer(const size_t&, const size_t&) method.\n"
              << "Index of layer is equal or greater than number of layers-1.\n";

       throw std::logic_error(buffer.str());
    }

    if(perceptrons_number == 0)
    {
       std::ostringstream buffer;

       buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
              << "void grow_layer(const size_t&, const size_t&) method.\n"
              << "Number of perceptrons must be greater than 0.\n";

       throw std::logic_error(buffer.str());
    }

    #endif

    for (size_t i = 0; i < perceptrons_number; i++)
    {
        layers[layer_index].grow_perceptrons(1);
        layers[layer_index+1].grow_inputs(1);
    }
}


// void prune_input(const size_t&) method


void MultilayerPerceptron::prune_input(const size_t& index)
{
    // Control sentence (if debug)

    #ifndef NDEBUG

    const size_t inputs_number = get_inputs_number();

    if(index >= inputs_number)
    {
       std::ostringstream buffer;

       buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
              << "void prune_input(const size_t&) method.\n"
              << "Index of input is equal or greater than number of inputs.\n";

       throw std::logic_error(buffer.str());
    }

    #endif

    const size_t layers_number = get_layers_number();

    if(layers_number > 0)
    {
        layers[0].prune_input(index);
    }
}


// void prune_output(const size_t&) method


void MultilayerPerceptron::prune_output(const size_t& index)
{
    // Control sentence (if debug)

    #ifndef NDEBUG

    const size_t outputs_number = get_outputs_number();

    if(index >= outputs_number)
    {
       std::ostringstream buffer;

       buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
              << "void prune_output(const size_t&) method.\n"
              << "Index of output is equal or greater than number of outputs.\n";

       throw std::logic_error(buffer.str());
    }

    #endif

    const size_t layers_number = get_layers_number();

    if(layers_number > 0)
    {
        layers[layers_number-1].prune_perceptron(index);
    }
}


// void prune_layer(const size_t&, const size_t&) method


void MultilayerPerceptron::prune_layer(const size_t& layer_index, const size_t& perceptron_index)
{
    #ifndef NDEBUG

    if(layer_index >= (layers.size()-1))
    {
       std::ostringstream buffer;

       buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
              << "void prune_layer(const size_t&, const size_t&) method.\n"
              << "Index of layer is equal or greater than number of layers-1.\n";

       throw std::logic_error(buffer.str());
    }

    if(perceptron_index >= layers[layer_index].get_perceptrons_number())
    {
       std::ostringstream buffer;

       buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
              << "void prune_layer(const size_t&, const size_t&) method.\n"
              << "Index of perceptron is equal or greater than number of perceptrons in the layer.\n";

       throw std::logic_error(buffer.str());
    }

    #endif

    layers[layer_index].prune_perceptron(perceptron_index);
    layers[layer_index+1].prune_input(perceptron_index);
}


// void initialize_random(void) method


void MultilayerPerceptron::initialize_random(void)
{
   const size_t architecture_size = rand()%10 + 2;

   Vector<size_t> architecture(architecture_size);

   for(size_t i = 0; i < architecture_size; i++)
   {
      architecture[i]  = rand()%10 + 1;
   }

   set(architecture);

   const size_t layers_number = get_layers_number();

   // Layers activation function

   for(size_t i = 0; i < layers_number; i++)
   {
      switch(rand()%5)
      {
         case 0:
         {
            layers[i].set_activation_function(Perceptron::Logistic);
         }
         break;

         case 1:
         {
            layers[i].set_activation_function(Perceptron::HyperbolicTangent);
         }
         break;

         case 2:
         {
            layers[i].set_activation_function(Perceptron::Threshold);
         }
         break;

         case 3:
         {
            layers[i].set_activation_function(Perceptron::SymmetricThreshold);
         }
         break;

         case 4:
         {
            layers[i].set_activation_function(Perceptron::Linear);
         }
         break;

         default:
         {
            std::ostringstream buffer;

            buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
                   << "void initialize_random(void) method.\n"
                   << "Unknown layer activation function.\n";
 
            throw std::logic_error(buffer.str());
         }
         break;
      }        
   }

   // Display messages
   
   set_display(true);
}


// void initialize_biases(const double&) method


void MultilayerPerceptron::initialize_biases(const double& value)
{
   const size_t layers_number = get_layers_number();

   for(size_t i = 0; i < layers_number; i++)
   {
      layers[i].initialize_biases(value);
   }
}


// void initialize_synaptic_weights(const double&) const method


void MultilayerPerceptron::initialize_synaptic_weights(const double& value) 
{
   const size_t layers_number = get_layers_number();
   const Vector<size_t> layers_size = arrange_layers_perceptrons_numbers();

   for(size_t i = 0; i < layers_number; i++)
   {
      layers[i].initialize_synaptic_weights(value);
   }
}


// void initialize_parameters(const double&) method


void MultilayerPerceptron::initialize_parameters(const double& value)
{
   const size_t parameters_number = count_parameters_number();

   const Vector<double> parameters(parameters_number, value);

   set_parameters(parameters);
}


// void randomize_parameters_uniform(void) method


void MultilayerPerceptron::randomize_parameters_uniform(void)
{
   const size_t parameters_number = count_parameters_number();

   Vector<double> parameters(parameters_number);

   parameters.randomize_uniform();

   set_parameters(parameters);
}


// void randomize_parameters_uniform(const double&, const double&) method


void MultilayerPerceptron::randomize_parameters_uniform(const double& minimum, const double& maximum)
{
   const size_t parameters_number = count_parameters_number();

   Vector<double> parameters(parameters_number);

   parameters.randomize_uniform(minimum, maximum);

   set_parameters(parameters);
}


// void randomize_parameters_uniform(const Vector<double>&, const Vector<double>&) method


void MultilayerPerceptron::randomize_parameters_uniform(const Vector<double>& minimum, const Vector<double>& maximum)
{
   const size_t parameters_number = count_parameters_number();

   Vector<double> parameters(parameters_number);

   parameters.randomize_uniform(minimum, maximum);

   set_parameters(parameters);
}


// void initializeparameters_uniform(const Vector< Vector<double> >&) method


void MultilayerPerceptron::randomize_parameters_uniform(const Vector< Vector<double> >& minimum_maximum)
{
   const size_t parameters_number = count_parameters_number();

   Vector<double> parameters(parameters_number);

   parameters.randomize_uniform(minimum_maximum[0], minimum_maximum[1]);

   set_parameters(parameters);
}


// void randomize_parameters_normal(void) method


void MultilayerPerceptron::randomize_parameters_normal(void)
{
   const size_t parameters_number = count_parameters_number();

   Vector<double> parameters(parameters_number);

   parameters.randomize_normal();

   set_parameters(parameters);
}


// void randomize_parameters_normal(const double&, const double&) method


void MultilayerPerceptron::randomize_parameters_normal(const double& mean, const double& standard_deviation)
{
   const size_t parameters_number = count_parameters_number();

   Vector<double> parameters(parameters_number);

   parameters.randomize_normal(mean, standard_deviation);

   set_parameters(parameters);
}


// void randomize_parameters_normal(const Vector<double>&, const Vector<double>&) method


void MultilayerPerceptron::randomize_parameters_normal(const Vector<double>& mean, const Vector<double>& standard_deviation)
{
   const size_t parameters_number = count_parameters_number();

   Vector<double> parameters(parameters_number);

   parameters.randomize_normal(mean, standard_deviation);

   set_parameters(parameters);
}


// void randomize_parameters_normal(const Vector< Vector<double> >&) method


void MultilayerPerceptron::randomize_parameters_normal(const Vector< Vector<double> >& mean_standard_deviation)
{
   const size_t parameters_number = count_parameters_number();

   Vector<double> parameters(parameters_number);

   parameters.randomize_normal(mean_standard_deviation[0], mean_standard_deviation[1]);

   set_parameters(parameters);
}


// void initialize_parameters(void) method


void MultilayerPerceptron::initialize_parameters(void)
{
   randomize_parameters_normal();
}


// void perturbate_parameters(const double&) method


void MultilayerPerceptron::perturbate_parameters(const double& perturbation)
{
    #ifndef NDEBUG

    if(perturbation < 0)
    {
       std::ostringstream buffer;

       buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
              << "void perturbate_parameters(const double&) method.\n"
              << "Perturbation must be equal or greater than 0.\n";

       throw std::logic_error(buffer.str());
    }

    #endif

    Vector<double> parameters = arrange_parameters();

    Vector<double> parameters_perturbation(parameters);

    parameters_perturbation.randomize_uniform(-perturbation,perturbation);

    parameters = parameters + parameters_perturbation;

    set_parameters(parameters);
}


// double calculate_parameters_norm(void) const const method


double MultilayerPerceptron::calculate_parameters_norm(void) const 
{
   const Vector<double> parameters = arrange_parameters();

   const double parameters_norm = parameters.calculate_norm();

   return(parameters_norm);
}


// Vector<double> calculate_outputs(const Vector<double>&) const method


Vector<double> MultilayerPerceptron::calculate_outputs(const Vector<double>& inputs) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t size = inputs.size();

   const size_t inputs_number = get_inputs_number();

   if(size != inputs_number) 
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector<double> calculate_outputs(const Vector<double>&) const method.\n"
             << "Size of inputs (" << size <<") must be equal to number of inputs (" << inputs_number << ").\n";

      throw std::logic_error(buffer.str());
   }   
   
   #endif

   const size_t layers_number = get_layers_number();

   Vector<double> outputs;

   if(layers_number == 0)
   {
      return(outputs);   
   }
   else
   {       
      outputs = layers[0].calculate_outputs(inputs);

      for(size_t i = 1; i < layers_number; i++)
      {
         outputs = layers[i].calculate_outputs(outputs);
      }
   }

   return(outputs);
}


// Matrix<double> calculate_Jacobian(const Vector<double>&) const method


Matrix<double> MultilayerPerceptron::calculate_Jacobian(const Vector<double>& inputs) const
{
   #ifndef NDEBUG 

   const size_t size = inputs.size();

   const size_t inputs_number = get_inputs_number();

   if(size != inputs_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Matrix<double> calculate_Jacobian(const Vector<double>&) const method.\n"
             << "Size must be equal to number of inputs.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const size_t layers_number = get_layers_number();

   if(layers_number == 0)
   {
      Matrix<double> Jacobian;

      return(Jacobian);
   }
   else
   {   
      const Vector< Matrix<double> > layers_Jacobian = calculate_layers_Jacobian(inputs);

      Matrix<double> Jacobian = layers_Jacobian[layers_number-1];

      for(int i = (int)layers_number-2; i > -1; i--)
      {
         Jacobian = Jacobian.dot(layers_Jacobian[i]);   
      }

      return(Jacobian);
   }
}


// Vector< Matrix<double> > calculate_Hessian_form(const Vector<double>&) const


Vector< Matrix<double> > MultilayerPerceptron::calculate_Hessian_form(const Vector<double>& inputs) const
{   
   #ifndef NDEBUG 

   const size_t size = inputs.size();

   const size_t inputs_number = get_inputs_number();

   if(size != inputs_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector< Matrix<double> > calculate_Hessian_form(const Vector<double>&) const method.\n"
             << "Size must be equal to number of inputs.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const size_t layers_number = get_layers_number();
   const Vector<size_t> layers_size = arrange_layers_perceptrons_numbers();
   const size_t outputs_number = get_outputs_number();

   // Calculate forward propagation of outputs, Jacobian and Hessian form

   const Vector< Matrix<double> > layers_Jacobian = calculate_layers_Jacobian(inputs);

   const Vector< Vector< Matrix<double> > > layers_Hessian_form = calculate_layers_Hessian_form(inputs);

   // Calculate multilayer_perceptron_pointer Hessian form

   Vector< Matrix<double> > Hessian_form = layers_Hessian_form[layers_number-1];
   Matrix<double> Jacobian = layers_Jacobian[layers_number-1];

   for(size_t layer_index = layers_number-2; layer_index != 0; layer_index--)
   {
      for(size_t output_index = 0; output_index < outputs_number; output_index++)
      {
         Hessian_form[output_index] = (layers_Jacobian[layer_index].calculate_transpose()).dot(Hessian_form[output_index]).dot(layers_Jacobian[layer_index]);

         for(size_t neuron_index = 0; neuron_index < layers_size[layer_index]; neuron_index++)
         {
            Hessian_form[output_index] += layers_Hessian_form[layer_index][neuron_index]*Jacobian(output_index,neuron_index);
         }
      }

     Jacobian = Jacobian.dot(layers_Jacobian[layer_index]);   
   }

//   Matrix<double> hidden_layer_Jacobian = layers_Jacobian[0];
//   Matrix<double> output_layer_Jacobian = layers_Jacobian[1];

//   const Vector< Matrix<double> > hidden_layer_Hessian_form = layers_Hessian_form[0];
//   const Vector< Matrix<double> > output_layer_Hessian_form = layers_Hessian_form[1];

//   const Vector< Matrix<double> > Hessian_form(outputs_number);

//   for(size_t output_index = 0; output_index < outputs_number; output_index++)
//   {
//      Hessian_form[output_index] = (hidden_layer_Jacobian.calculate_transpose()).dot(output_layer_Hessian_form[output_index]).dot(hidden_layer_Jacobian);

//      for(size_t neuron_index = 0; neuron_index < layers_size[0]; neuron_index++)
//      {
//         Hessian_form[output_index] += hidden_layer_Hessian_form[neuron_index]*output_layer_Jacobian[output_index][neuron_index]; 
//    }
//   }

   return(Hessian_form);
}


// Vector<double> calculate_layer_combination_combination(const size_t&, const Vector<double>&) const method


Vector<double> MultilayerPerceptron::calculate_layer_combination_combination(const size_t& layer_index, const Vector<double>& previous_layer_combination) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t layers_number = get_layers_number();

   if(layer_index == 0)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Matrix<double> calculate_layer_combination_combination(const size_t&, const Vector<double>&) const.\n"
             << "Index of layer must be greater than zero.\n";

      throw std::logic_error(buffer.str());
   }
   else if(layer_index >= layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Matrix<double> calculate_layer_combination_combination(const size_t&, const Vector<double>&) const.\n"
             << "Index of layer must be less than numnber of layers.\n";

      throw std::logic_error(buffer.str());
   }

   const size_t size = previous_layer_combination.size();
   const size_t previous_layer_perceptrons_number = layers[layer_index-1].get_perceptrons_number();

   if(size != previous_layer_perceptrons_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Matrix<double> calculate_layer_combination_combination(const size_t&, const Vector<double>&) const.\n"
             << "Size must be equal to size of previous layer.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   // Calculate combination to layer

   const Vector<double> previous_layer_activation = layers[layer_index-1].calculate_activations(previous_layer_combination);

   const Vector<double> layer_combination_combination = layers[layer_index].calculate_combinations(previous_layer_activation);

   return(layer_combination_combination);
}


// Matrix<double> calculate_layer_combination_combination_Jacobian(const size_t&, const Vector<double>&) const method


Matrix<double> MultilayerPerceptron::calculate_layer_combination_combination_Jacobian(const size_t& layer_index, const Vector<double>& previous_layer_activation_derivative) const
{
   const size_t previous_layer_perceptrons_number = layers[layer_index-1].get_perceptrons_number();

   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t layers_number = get_layers_number();

   if(layer_index == 0)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Matrix<double> calculate_layer_combination_combination_Jacobian(const size_t&, const Vector<double>&) const.\n"
             << "Index of layer must be greater than zero.\n";

      throw std::logic_error(buffer.str());
   }
   else if(layer_index >= layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Matrix<double> calculate_layer_combination_combination_Jacobian(const size_t&, const Vector<double>&) const.\n"
             << "Index of layer must be less than numnber of layers.\n";

      throw std::logic_error(buffer.str());
   }

   const size_t size = previous_layer_activation_derivative.size();

   if(size != previous_layer_perceptrons_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Matrix<double> calculate_layer_combination_combination_Jacobian(const size_t&, const Vector<double>&).\n"
             << "Size of activation derivative must be equal to size of previous layer.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const Matrix<double> layer_synaptic_weights = layers[layer_index].arrange_synaptic_weights();

   Matrix<double> previous_layer_activation_Jacobian(previous_layer_perceptrons_number, previous_layer_perceptrons_number, 0.0);
   previous_layer_activation_Jacobian.set_diagonal(previous_layer_activation_derivative);

   return(layer_synaptic_weights.dot(previous_layer_activation_Jacobian));
}


// Vector<double> calculate_interlayer_combination_combination(const size_t&, const size_t&, const Vector<double>&) const method


Vector<double> MultilayerPerceptron::calculate_interlayer_combination_combination(const size_t& domain_layer_index, const size_t& image_layer_index, const Vector<double>& domain_layer_combination) const
{
   if(domain_layer_index > image_layer_index)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector<double> calculate_interlayer_combination_combination(const size_t&, const size_t&, const Vector<double>&) const method.\n"
             << "Index of domain layer must be less or equal than index of image layer.\n";

      throw std::logic_error(buffer.str());
   }

   if(domain_layer_index == image_layer_index)
   {
      return(domain_layer_combination);
   }
   else
   {
      Vector<double> interlayer_combination_combination(domain_layer_combination);     

      for(size_t i = domain_layer_index+1; i <= image_layer_index; i++)
      {
         interlayer_combination_combination = calculate_layer_combination_combination(i, interlayer_combination_combination);
      }

      return(interlayer_combination_combination); 
   }
}


// Matrix<double> calculate_interlayer_combination_combination_Jacobian(const size_t&, const size_t&, const Vector<double>&) const method


Matrix<double> MultilayerPerceptron::calculate_interlayer_combination_combination_Jacobian(const size_t& domain_layer_index, const size_t& image_layer_index, const Vector<double>& domain_layer_combination) const
{
   Matrix<double> interlayer_combination_combination_Jacobian;

   if(domain_layer_index < image_layer_index)
   {
      const size_t size = image_layer_index-domain_layer_index;

      Vector< Vector<double> > layers_combination_combination(size);
      Vector< Matrix<double> > layers_combination_combination_Jacobian(size);

      layers_combination_combination[0] = calculate_layer_combination_combination(domain_layer_index+1, domain_layer_combination);      
      layers_combination_combination_Jacobian[0] = calculate_layer_combination_combination_Jacobian(domain_layer_index+1, layers[domain_layer_index].calculate_activations_derivatives(domain_layer_combination));

      for(size_t i = 1; i < size; i++)
      {
         layers_combination_combination[i] = calculate_layer_combination_combination(domain_layer_index+i+1, layers_combination_combination[i-1]);

         layers_combination_combination_Jacobian[i] = calculate_layer_combination_combination_Jacobian(domain_layer_index+i+1, layers_combination_combination[i-1]); 
      }

      interlayer_combination_combination_Jacobian = layers_combination_combination_Jacobian[size-1];

      for(int i = (int)size-2; i > -1; i--)
      {
         interlayer_combination_combination_Jacobian = interlayer_combination_combination_Jacobian.dot(layers_combination_combination_Jacobian[i]);   
      }
   }
   else if(domain_layer_index == image_layer_index)
   {
      const size_t image_layer_perceptrons_number = layers[image_layer_index].get_perceptrons_number();

      interlayer_combination_combination_Jacobian.set_identity(image_layer_perceptrons_number);
   }
   else
   {
      const size_t image_layer_perceptrons_number = layers[image_layer_index].get_perceptrons_number();
      const size_t domain_layer_perceptrons_number = layers[domain_layer_index].get_perceptrons_number();

      interlayer_combination_combination_Jacobian.set(image_layer_perceptrons_number, domain_layer_perceptrons_number, 0.0);      
   }

   return(interlayer_combination_combination_Jacobian); 
}


// Vector<double> calculate_output_layer_combination(const size_t&, const Vector<double>&) const method


Vector<double> MultilayerPerceptron::calculate_output_layer_combination(const size_t& layer_index, const Vector<double>& layer_combinations) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t size = layer_combinations.size();

   const size_t layer_perceptrons_number = layers[layer_index].get_perceptrons_number();

   if(size != layer_perceptrons_number) 
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector<double> calculate_output_layer_combination(const size_t&, const Vector<double>&) const method.\n"
             << "Size must be equal to layer size.\n";

      throw std::logic_error(buffer.str());
   }   
   
   #endif

   const size_t layers_number = get_layers_number();

   Vector<double> outputs = layers[layer_index].calculate_activations(layer_combinations);

   if(layer_index < layers_number)
   {
      for(size_t i = layer_index+1; i < layers_number; i++)
      {
         outputs = layers[i].calculate_outputs(outputs);
      }
   }

   return(outputs);
}


// Vector< Matrix<double> > calculate_output_layers_delta(const Vector< Vector<double> >&, const Vector<double>&) method


Vector< Matrix<double> > MultilayerPerceptron::calculate_output_layers_delta(const Vector< Vector<double> >& layers_activation_derivative) const
{
   // Control sentence (if debug)

   const size_t layers_number = get_layers_number();
   const Vector<size_t> layers_size = arrange_layers_perceptrons_numbers();

   // Control sentence (if debug)

   #ifndef NDEBUG 

   // Forward propagation activation derivative size

   const size_t layers_activation_derivative_size = layers_activation_derivative.size();

   if(layers_activation_derivative_size != layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector< Vector<double> > calculate_output_layers_delta(const Vector< Vector<double> >&) method.\n"
             << "Size of forward propagation activation derivative vector must be equal to number of layers.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   Vector< Matrix<double> > output_layers_delta(layers_number);

   if(layers_number > 0)
   {
      // Output layer

      output_layers_delta[layers_number-1] = layers[layers_number-1].arrange_activations_Jacobian(layers_activation_derivative[layers_number-1]);

      // Rest of layers

      for(int i = (int)layers_number-2; i >= 0; i--)
      {   
         output_layers_delta[i] = output_layers_delta[i+1].dot(layers[i+1].arrange_synaptic_weights()).dot(layers[i].arrange_activations_Jacobian(layers_activation_derivative[i]));
      }
   }

   return(output_layers_delta);
}


// Matrix< Vector< Matrix<double> > > calculate_output_interlayers_Delta(void) method 


Matrix< Vector< Matrix<double> > > MultilayerPerceptron::calculate_output_interlayers_Delta
(const Vector< Vector< Vector<double> > >& second_order_forward_propagation, 
 const Matrix< Matrix<double> >& interlayers_combination_combination_Jacobian,
 const Vector< Matrix<double> >& output_layers_delta) const
{
   const size_t layers_number = get_layers_number();
   const Vector<size_t> layers_size = arrange_layers_perceptrons_numbers();
   const size_t outputs_number = get_outputs_number();

   // Control sentence (if debug)

   #ifndef NDEBUG 

   std::ostringstream buffer;

   // Second order forward propagation

   const size_t second_order_forward_propagation_size = second_order_forward_propagation.size();

   if(second_order_forward_propagation_size != 3)
   {
      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Matrix< Vector< Matrix<double> > > calculate_output_interlayers_delta() method.\n"
             << "Size of second order forward propagation must be three.\n";

      throw std::logic_error(buffer.str());
   }

   // Interlayers combination-combination Jacobian

   const size_t interlayers_combination_combination_Jacobian_rows_number = interlayers_combination_combination_Jacobian.get_rows_number();

   if(interlayers_combination_combination_Jacobian_rows_number != layers_number)
   {
      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Matrix< Vector< Matrix<double> > > calculate_output_interlayers_delta() method.\n"
             << "Number of rows of interlayers combination-combination Jacobian must be equal to number of layers.\n";

      throw std::logic_error(buffer.str());
   }

   // Multilayer perceptron outputs layers delta

   const size_t output_layers_delta_size = output_layers_delta.size();

   if(output_layers_delta_size != layers_number)
   {
      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Matrix< Vector< Matrix<double> > > calculate_output_interlayers_delta() method.\n"
             << "Size of multilayer perceptron outputs layers delta must be equal to number of layers.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const Vector< Matrix<double> > layers_synaptic_weights = arrange_layers_synaptic_weights();

   const Vector< Vector<double> >& layers_activation_derivative = second_order_forward_propagation[1];
   const Vector< Vector<double> >& layers_activation_second_derivative = second_order_forward_propagation[2];

   // The Delta form consists of a square matrix of size the number of layers

   Matrix< Vector< Matrix<double> > > output_interlayers_Delta(layers_number, layers_number);

   // Each Delta element is a vector of size the number of outputs

   for(size_t i = 0; i < layers_number; i++)
   {
      for(size_t j = 0; j < layers_number; j++)
      {
         output_interlayers_Delta(i,j).set(outputs_number);
      }
   }

   // The subelements are matrices with rows and columns numbers the sizes of the first and the second layers

   for(size_t i = 0; i < layers_number; i++)
   {
      for(size_t j = 0; j < layers_number; j++)
      {
         for(size_t k = 0; k < outputs_number; k++)
         {
            output_interlayers_Delta(i,j)[k].set(layers_size[i], layers_size[j]);
         }
      }
   }


   if(layers_number > 0)
   {
      // Output-outputs layer (OK)

      output_interlayers_Delta(layers_number-1,layers_number-1) = layers[layers_number-1].arrange_activations_Hessian_form(layers_activation_second_derivative[layers_number-1]);

      // Rest of hidden layers

      double sum_1;
      double sum_2;

      for(size_t L = layers_number-1; L != 0; L--)
      {   
         for(size_t M = layers_number-1; M >= L; M--)
         {   
            if(!(L == layers_number-1 && M == layers_number-1))
            {
              for(size_t i = 0; i < outputs_number; i++)
               {
                  for(size_t j = 0; j < layers_size[L]; j++)
                  {  
                     sum_1 = 0.0;
         
                     for(size_t l = 0; l < layers_size[L+1]; l++)
                     {
                        sum_1 += output_layers_delta[L+1](i,l)*layers_synaptic_weights[L+1](l,j);
                     }

                     for(size_t k = 0; k < layers_size[M]; k++)
                     {
                        sum_2 = 0.0;

                        for(size_t l = 0; l < layers_size[L+1]; l++)
                        {
                           sum_2 += output_interlayers_Delta(L+1,M)[i](l,k)*layers_synaptic_weights[L+1](l,j);
                        }

                        output_interlayers_Delta(L,M)[i](j,k)
                        = layers_activation_second_derivative[L][j]
                        *interlayers_combination_combination_Jacobian(L,M)(j,k)*sum_1
                        + layers_activation_derivative[L][j]*sum_2;

                        output_interlayers_Delta(M,L)[i](k,j) = output_interlayers_Delta(L,M)[i](j,k);
                     }       
                  }
               }
            }
         }       
      }
   }

   return(output_interlayers_Delta);
}


// Vector<double> calculate_parametrs_outputs(const Vector<double>&, const Vector<double>&) const method


Vector<double> MultilayerPerceptron::calculate_outputs(const Vector<double>& inputs, const Vector<double>& parameters) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG

   const size_t size = inputs.size();

   const size_t inputs_number = get_inputs_number();

   if(size != inputs_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector<double> calculate_outputs(const Vector<double>&, const Vector<double>&) const method.\n"
             << "Size of inputs (" << size <<") must be equal to number of inputs (" << inputs_number << ").\n";

      throw std::logic_error(buffer.str());
   }

   const size_t parameters_size = parameters.size();

   const size_t parameters_number = count_parameters_number();

   if(parameters_size != parameters_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector<double> calculate_outputs(const Vector<double>&, const Vector<double>&) const method.\n"
             << "Size of parameters (" << parameters_size <<") must be equal to number of parameters (" << parameters_number << ").\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const size_t layers_number = get_layers_number();

   Vector<double> outputs;

   if(layers_number == 0)
   {
      return(outputs);
   }
   else
   {
      const Vector<size_t> layers_parameters_numbers = count_layers_parameters_numbers();

      const Vector<size_t> layers_cumulative_parameters_number = arrange_layers_cumulative_parameters_number();

      Vector<double> layer_parameters;

      layer_parameters = parameters.take_out(0, layers_parameters_numbers[0]);

      outputs = layers[0].calculate_outputs(inputs, layer_parameters);

      for(size_t i = 1; i < layers_number; i++)
      {
          layer_parameters = parameters.take_out(layers_cumulative_parameters_number[i-1], layers_parameters_numbers[i]);

         outputs = layers[i].calculate_outputs(outputs, layer_parameters);
      }
   }

   return(outputs);
}


// Matrix<double> calculate_parameters_Jacobian(const Vector<double>&, const Vector<double>&) const method


Matrix<double> MultilayerPerceptron::calculate_Jacobian(const Vector<double>& inputs, const Vector<double>&) const
{   
   #ifndef NDEBUG 

   const size_t inputs_number = get_inputs_number();
   const size_t size = inputs.size();

   if(size != inputs_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "void calculate_Jacobian(Vector<double>&, const Vector<double>&) const method.\n"
             << "Size must be equal to number of inputs.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   // Neural network stuff

   const size_t layers_number = get_layers_number();
   const Vector<size_t> layers_size = arrange_layers_perceptrons_numbers();

   // Calculate required quantities

   const Vector< Vector< Vector<double> > > first_order_forward_propagation = calculate_first_order_forward_propagation(inputs);
   const Vector< Vector<double> >& layers_activation = first_order_forward_propagation[0];
   const Vector< Vector<double> >& layers_activation_derivative = first_order_forward_propagation[1];

   const Vector< Vector<double> > layers_inputs = arrange_layers_input(inputs, layers_activation);

   const Vector< Matrix<double> > layers_combination_parameters_Jacobian = calculate_layers_combination_parameters_Jacobian(layers_inputs);

   const Vector< Matrix<double> > output_layers_delta = calculate_output_layers_delta(layers_activation_derivative);

   Matrix<double> parameters_Jacobian = output_layers_delta[0].dot(layers_combination_parameters_Jacobian[0]);

   for(size_t i = 1; i < layers_number; i++)
   {
      parameters_Jacobian = parameters_Jacobian.assemble_columns(output_layers_delta[i].dot(layers_combination_parameters_Jacobian[i]));
   }

   return(parameters_Jacobian);
}


// Vector< Matrix<double> > calculate_Hessian_form(const Vector<double>&, const Vector<double>&) const


Vector< Matrix<double> > MultilayerPerceptron::calculate_Hessian_form(const Vector<double>& inputs, const Vector<double>&) const
{
   // Neural network stuff

   const size_t layers_number = get_layers_number();
   const Vector<size_t> layers_size = arrange_layers_perceptrons_numbers();
   const size_t outputs_number = get_outputs_number();

   const size_t parameters_number = count_parameters_number();

   const Vector< Matrix<double> > layers_synaptic_weights = arrange_layers_synaptic_weights();

   // Calculate required quantities

   const Vector< Vector< Vector<double> > > second_order_forward_propagation = calculate_second_order_forward_propagation(inputs);

   const Vector< Vector<double> >& layers_activation = second_order_forward_propagation[0];
   const Vector< Vector<double> >& layers_activation_derivative = second_order_forward_propagation[1];
//   const Vector< Vector<double> >& layers_activation_second_derivative = second_order_forward_propagation[2];

   const Vector< Vector<double> > layers_inputs = arrange_layers_input(inputs, layers_activation);

   const Vector< Vector< Vector<double> > > perceptrons_combination_parameters_gradient = calculate_perceptrons_combination_parameters_gradient(layers_inputs);

   const Matrix< Matrix<double> > interlayers_combination_combination_Jacobian = calculate_interlayers_combination_combination_Jacobian(inputs);

   const Vector< Matrix<double> > output_layers_delta = calculate_output_layers_delta(layers_activation_derivative);

   const Matrix< Vector< Matrix<double> > > output_interlayers_Delta
   = calculate_output_interlayers_Delta(second_order_forward_propagation, interlayers_combination_combination_Jacobian, output_layers_delta);

   // Size multilayer_perceptron_pointer parameters Hessian form 

   Vector< Matrix<double> > parameters_Hessian_form(outputs_number);

   for(size_t i = 0; i < outputs_number; i++)
   {
      parameters_Hessian_form[i].set(parameters_number, parameters_number);
   }

   // Calculate Hessian form elements

   if(layers_number > 0)
   {
      const Matrix<size_t> parameters_indices = arrange_parameters_indices();

      size_t layer_j;
      size_t neuron_j;
      size_t parameter_j;

      size_t layer_k;
      size_t neuron_k;
      size_t parameter_k;

      Vector<double> neuron_j_combination_parameters_gradient;
      Vector<double> neuron_k_combination_parameters_gradient;

      for(size_t i = 0; i < outputs_number; i++)
      {
         for(size_t j = 0; j < parameters_number; j++)
         {
            layer_j = parameters_indices(j,0);
            neuron_j = parameters_indices(j,1);
            parameter_j = parameters_indices(j,2);

            for(size_t k = j; k < parameters_number; k++)
            {
               layer_k = parameters_indices(k,0);
               neuron_k = parameters_indices(k,1);
               parameter_k = parameters_indices(k,2);
                       
               parameters_Hessian_form[i](j,k) =
               output_interlayers_Delta(layer_j,layer_k)[i](neuron_j,neuron_k)
               *perceptrons_combination_parameters_gradient[layer_j][neuron_j][parameter_j]
               *perceptrons_combination_parameters_gradient[layer_k][neuron_k][parameter_k]
               ;

               if(layer_j != 0 && parameter_j != 0) // Neither the first layer nor parameter is a bias
               {
                  //parameters_Hessian_form[i](j,k) +=
                  //output_layers_delta[layer_j][i][neuron_j]
                  //*layers_activation_derivative[layer_j-1][parameter_j-1]
                  //*interlayers_combination_combination_Jacobian[layer_j-1][layer_k][parameter_j-1][neuron_k]
                  //*perceptrons_combination_parameters_gradient[layer_k][neuron_k][parameter_k]
                  ;               
               }

               if(layer_k != 0 && parameter_k != 0)
               {
                 parameters_Hessian_form[i](j,k) +=
                 output_layers_delta[layer_k](i,neuron_k)
                 *layers_activation_derivative[layer_k-1][parameter_k-1]
                 *interlayers_combination_combination_Jacobian(layer_k-1,layer_j)(parameter_k-1,neuron_j)
                 *perceptrons_combination_parameters_gradient[layer_j][neuron_j][parameter_j]
                 ;         
               }  

               parameters_Hessian_form[i](k,j) = parameters_Hessian_form[i](j,k);
            }
         }
      }
   }

   return(parameters_Hessian_form);
}


// Vector< Vector<double> > arrange_layers_input(const Vector<double>&, const Vector< Vector<double> >&) const method


Vector< Vector<double> > MultilayerPerceptron::arrange_layers_input(const Vector<double>& inputs, const Vector< Vector<double> >& layers_activation) const
{
   // Neural network stuff

   const size_t layers_number = get_layers_number();

   #ifndef NDEBUG 

   const size_t inputs_number = get_inputs_number();
   const size_t inputs_size = inputs.size();

   if(inputs_size != inputs_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector< Vector<double> > arrange_layers_input(const Vector<double>&, const Vector< Vector<double> >&) const method.\n"
             << "Size must be equal to number of inputs.\n";

      throw std::logic_error(buffer.str());
   }

   const size_t layers_activation_size = layers_activation.size();

   if(layers_activation_size != layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector< Vector<double> > arrange_layers_input(const Vector<double>&, const Vector< Vector<double> >&) const method.\n"
             << "Size must be equal to number of inputs.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   // Arrange layers inputs

   Vector< Vector<double> > layers_inputs(layers_number);

   if(layers_number != 0)
   {
      layers_inputs[0] = inputs;

      for(size_t j = 1; j < layers_number; j++)
      {
         layers_inputs[j] = layers_activation[j-1];
      }
   }

   return(layers_inputs);
}


// Vector< Vector<double> > arrange_layers_input(const Vector<double>&, const Vector< Vector< Vector<double> > >&) const method


Vector< Vector<double> > MultilayerPerceptron::arrange_layers_input(const Vector<double>& input, const Vector< Vector< Vector<double> > >& first_order_forward_propagation) const
{
    const size_t layers_number = get_layers_number();

    Vector< Vector<double> > layers_inputs(layers_number);

    const Vector< Vector<double> >& layers_activation = first_order_forward_propagation[0];

    layers_inputs[0] = input;

    for(size_t j = 1; j < layers_number; j++)
    {
       layers_inputs[j] = layers_activation[j-1];
    }

    return(layers_inputs);
}


// Vector< Vector<double> > calculate_layers_input(const Vector<double>&) const method


Vector< Vector<double> > MultilayerPerceptron::calculate_layers_input(const Vector<double>& inputs) const
{
   const size_t layers_number = get_layers_number();

   Vector< Vector<double> > layers_inputs(layers_number);

   layers_inputs[0] = inputs;

   for(size_t i = 1; i < layers_number; i++)
   {
      layers_inputs[i] = layers[i-1].calculate_outputs(layers_inputs[i-1]);
   }

   return(layers_inputs);
}


// Vector< Vector<double> > calculate_layers_combination(const Vector<double>&) const method


Vector< Vector<double> > MultilayerPerceptron::calculate_layers_combination(const Vector<double>& inputs) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t inputs_size = inputs.size();

   const size_t inputs_number = get_inputs_number();

   if(inputs_size != inputs_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector< Vector<double> > calculate_layers_combination(const Vector<double>&) const method.\n"
             << "Size must be equal to number of inputs.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const size_t layers_number = get_layers_number();

   Vector< Vector<double> > layers_combination(layers_number);

   if(layers_number > 0)
   {
      Vector< Vector<double> > layers_activation(layers_number);

      layers_combination[0] = layers[0].calculate_combinations(inputs);

      layers_activation[0] = layers[0].calculate_activations(layers_combination[0]);

      for(size_t i = 1; i < layers_number; i++)
      {
         layers_combination[i] = layers[i].calculate_combinations(layers_activation[i-1]);

         layers_activation[i] = layers[i].calculate_activations(layers_combination[i]);
      }
   }

   return(layers_combination);
}


// Vector< Matrix<double> > calculate_layers_combination_Jacobian(const Vector<double>&) const method


Vector< Matrix<double> > MultilayerPerceptron::calculate_layers_combination_Jacobian(const Vector<double>& inputs) const
{
   const size_t layers_number = get_layers_number();

   Vector< Matrix<double> > layers_combination_Jacobian(layers_number-2);

   Vector< Vector<double> > layers_output(layers_number);

   Vector< Matrix<double> > layers_Jacobian(layers_number);

   // Hidden and outputs layers Jacobian

   layers_output[0] = layers[0].calculate_outputs(inputs);   

   layers_combination_Jacobian[0] = layers[0].calculate_combinations_Jacobian(inputs);

   for(size_t i = 1; i < layers_number; i++)
   {
      layers_output[i] = layers[i].calculate_outputs(layers_output[i-1]);         

      layers_Jacobian[i] = layers[i].calculate_Jacobian(layers_output[i-1]);
   }

   Matrix<double> output_layer_Jacobian = layers[layers_number-1].calculate_Jacobian(layers_output[layers_number-1]);

   // Calculate forward propagation Jacobian

   layers_Jacobian[layers_number] = output_layer_Jacobian;

   for(int i = (int)layers_number-1; i > -1; i--)
   {
      layers_Jacobian[layers_number] = layers_Jacobian[layers_number].dot(layers_Jacobian[i]);   
   }

   for(int i = (int)layers_number-1; i > -1; i--)
   {
      layers_Jacobian[i] = layers_Jacobian[i];

      for(int j = (int)i-1; j > -1; j--)
      {
         layers_Jacobian[i] = layers_Jacobian[i].dot(layers_Jacobian[j]);     
      }      
   }

   return(layers_combination_Jacobian);
}


// Vector< Matrix<double> > calculate_layers_combination_parameters_Jacobian(const Vector< Vector<double> >&) const method


Vector< Matrix<double> > MultilayerPerceptron::calculate_layers_combination_parameters_Jacobian(const Vector< Vector<double> >& layers_inputs) const
{
   const size_t layers_number = get_layers_number();

   Vector< Matrix<double> > layers_combination_parameters_Jacobian(layers_number);

   const Vector<double> dummy;

   for(size_t i = 0; i < layers_number; i++)
   {
      layers_combination_parameters_Jacobian[i] = layers[i].calculate_combinations_Jacobian(layers_inputs[i], dummy);
   }

   return(layers_combination_parameters_Jacobian);
}


// Vector< Vector< Vector<double> > > calculate_perceptrons_combination_parameters_gradient(const Vector< Vector<double> >&) const method


Vector< Vector< Vector<double> > > MultilayerPerceptron::calculate_perceptrons_combination_parameters_gradient(const Vector< Vector<double> >& layers_inputs) const
{
   const size_t layers_number = get_layers_number();

   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t layers_input_size = layers_inputs.size();

   if(layers_input_size != layers_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector< Vector< Vector<double> > > calculate_perceptrons_combination_parameters_gradient(const Vector< Vector<double> >&) const method.\n"
             << "Size must be equal to number of layers.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const Vector<size_t> layers_inputs_number = get_layers_inputs_number();

   #ifndef NDEBUG 

   for(size_t i = 0; i < layers_number; i++)
   {
      if(layers_inputs[i].size() != layers_inputs_number[i])
      {
         std::ostringstream buffer;

         buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
                << "Vector< Vector< Vector<double> > > calculate_perceptrons_combination_parameters_gradient(const Vector< Vector<double> >&) const method.\n"
                << "Size of inputs to layer " << i << " must be equal to size of that layer.\n";

         throw std::logic_error(buffer.str());
      }
   }

   #endif

   const Vector<size_t> layers_size = arrange_layers_perceptrons_numbers();

   Vector<double> dummy_layer_parameters;

   Vector< Vector< Vector<double> > > perceptrons_combination_gradient(layers_number);

   for(size_t i = 0; i < layers_number; i++)
   {
      perceptrons_combination_gradient[i].set(layers_size[i]);

      for(size_t j = 0; j < layers_size[i]; j++)
      {
          const Perceptron& perceptron = layers[i].get_perceptron(j);

          perceptrons_combination_gradient[i][j] = perceptron.calculate_combination_gradient(layers_inputs[i], dummy_layer_parameters);
      }
   }

   return(perceptrons_combination_gradient);
}


// Vector< Vector<double> > calculate_layers_activation(const Vector<double>&) const method


Vector< Vector<double> > MultilayerPerceptron::calculate_layers_activation(const Vector<double>& inputs) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t inputs_size = inputs.size();

   const size_t inputs_number = get_inputs_number();

   if(inputs_size != inputs_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector< Vector<double> > calculate_layers_activation(const Vector<double>&) const method.\n"
             << "Size must be equal to number of inputs.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const size_t layers_number = get_layers_number();

   Vector< Vector<double> > layers_combination(layers_number);
   Vector< Vector<double> > layers_activation(layers_number);

   layers_combination[0] = layers[0].calculate_combinations(inputs);
   layers_activation[0] = layers[0].calculate_activations(layers_combination[0]);

   for(size_t i = 1; i < layers_number; i++)
   {
      layers_combination[i] = layers[i].calculate_combinations(layers_activation[i-1]);
      layers_activation[i] = layers[i].calculate_activations(layers_combination[i]);
   }

   return(layers_activation);
}


// Vector< Vector<double> > calculate_layers_activation_derivative(const Vector<double>&) const method


Vector< Vector<double> > MultilayerPerceptron::calculate_layers_activation_derivative(const Vector<double>& inputs) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t inputs_size = inputs.size();
   const size_t inputs_number = get_inputs_number();

   if(inputs_size != inputs_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector< Vector<double> > calculate_layers_activation_derivative(const Vector<double>&) const method.\n"
             << "Size must be equal to number of inputs.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const size_t layers_number = get_layers_number();

   Vector< Vector<double> > layers_combination(layers_number);
   Vector< Vector<double> > layers_activation(layers_number);
   Vector< Vector<double> > layers_activation_derivatives(layers_number);

   if(layers_number != 0)
   {   
      layers_combination[0] = layers[0].calculate_combinations(inputs);
      layers_activation[0] = layers[0].calculate_activations(layers_combination[0]);
      layers_activation_derivatives[0] = layers[0].calculate_activations_derivatives(layers_combination[0]);
  
      for(size_t i = 1; i < layers_number; i++)
      {
         layers_combination[i] = layers[i].calculate_combinations(layers_activation[i-1]);
         layers_activation[i] = layers[i].calculate_activations(layers_combination[i]);
         layers_activation_derivatives[i] = layers[i].calculate_activations_derivatives(layers_combination[i]);
      }
   }

   return(layers_activation_derivatives);
}


// Vector< Vector<double> > calculate_layers_activation_second_derivative(const Vector<double>&) const method


Vector< Vector<double> > MultilayerPerceptron::calculate_layers_activation_second_derivative(const Vector<double>& inputs) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t inputs_size = inputs.size();

   const size_t inputs_number = get_inputs_number();

   if(inputs_size != inputs_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector< Vector<double> > calculate_layers_activation_second_derivative(const Vector<double>&) const method.\n"
             << "Size must be equal to number of inputs.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const size_t layers_number = get_layers_number();

   Vector< Vector<double> > layers_combination(layers_number);
   Vector< Vector<double> > layers_activation(layers_number);
   Vector< Vector<double> > layers_activation_second_derivatives(layers_number);

   layers_combination[0] = layers[0].calculate_combinations(inputs);
   layers_activation[0] = layers[0].calculate_activations(layers_combination[0]);
   layers_activation_second_derivatives[0] = layers[0].calculate_activations_second_derivatives(layers_combination[0]);
  
   for(size_t i = 1; i < layers_number; i++)
   {
      layers_combination[i] = layers[i].calculate_combinations(layers_activation[i-1]);
   
      layers_activation[i] = layers[i].calculate_activations(layers_combination[i]);

      layers_activation_second_derivatives[i] = layers[i].calculate_activations_second_derivatives(layers_combination[i]);
   }

   return(layers_activation_second_derivatives);
}


// Vector< Matrix<double> > calculate_layers_Jacobian(const Vector<double>&) const method


Vector< Matrix<double> > MultilayerPerceptron::calculate_layers_Jacobian(const Vector<double>& inputs) const
{
   const size_t layers_number = get_layers_number();

   Vector<Vector<double> > layers_output(layers_number);
   Vector< Matrix<double> > layers_Jacobian(layers_number);

   layers_output[0] = layers[0].calculate_outputs(inputs);   
   layers_Jacobian[0] = layers[0].calculate_Jacobian(inputs);

   for(size_t i = 1; i < layers_number; i++)
   {
      layers_output[i] = layers[i].calculate_outputs(layers_output[i-1]);         

      layers_Jacobian[i] = layers[i].calculate_Jacobian(layers_output[i-1]);
   }

   return(layers_Jacobian);
}


// Vector< Matrix<double> > calculate_layers_Jacobian(const Vector<double>&) const method


Vector< Vector< Matrix<double> > > MultilayerPerceptron::calculate_layers_Hessian_form(const Vector<double>& inputs) const
{
   const size_t layers_number = get_layers_number();

   Vector<Vector<double> > layers_output(layers_number);
   Vector< Vector< Matrix<double> > > layers_Hessian_form(layers_number);

   layers_output[0] = layers[0].calculate_outputs(inputs);   
   layers_Hessian_form[0] = layers[0].calculate_Hessian_form(inputs);

   for(size_t i = 1; i < layers_number; i++)
   {
      layers_output[i] = layers[i].calculate_outputs(layers_output[i-1]);         

      layers_Hessian_form[i] = layers[i].calculate_Hessian_form(layers_output[i-1]);
   }

   return(layers_Hessian_form);
}


// Matrix< Matrix<double> > calculate_interlayers_combination_combination_Jacobian(const Vector<double>&) const method


Matrix< Matrix<double> > MultilayerPerceptron::calculate_interlayers_combination_combination_Jacobian(const Vector<double>& inputs) const  
{
   #ifndef NDEBUG 

   const size_t size = inputs.size();

   const size_t inputs_number = get_inputs_number();

   if(size != inputs_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Matrix< Matrix<double> > calculate_interlayers_combination_combination_Jacobian(const Vector<double>&) const method.\n"
             << "Size of inpouts must be equal to number of inputs.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const Vector< Vector<double> > layers_combination = calculate_layers_combination(inputs);

   return(calculate_interlayers_combination_combination_Jacobian(layers_combination));
}


// Matrix <Matrix<double> > calculate_interlayers_combination_combination_Jacobian(const Vector< Vector<double> >&) const method


Matrix < Matrix<double> > MultilayerPerceptron::calculate_interlayers_combination_combination_Jacobian(const Vector< Vector<double> >& layers_combination) const   
{
   const size_t layers_number = get_layers_number();

   #ifndef NDEBUG 

   size_t size;

   const Vector<size_t> layers_size = arrange_layers_perceptrons_numbers();

   for(size_t i = 0; i < layers_number; i++)
   {
      size = layers_combination[i].size();

      if(size != layers_size[i])
      {
         std::ostringstream buffer;

         buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
                << "Matrix< Matrix<double> > calculate_interlayers_combination_combination_Jacobian(const Vector< Vector<double> >&) const method.\n"
                << "Size must be equal to size of layer.\n";

         throw std::logic_error(buffer.str());
      }
   }

   #endif

   // Calculate interlayers combination combination Jacobian

   Matrix< Matrix<double> > interlayers_combination_combination_Jacobian(layers_number, layers_number);

   for(size_t image_index = 0; image_index < layers_number; image_index++)
   {
      for(size_t domain_index = 0; domain_index < layers_number; domain_index++)
      {
         interlayers_combination_combination_Jacobian(image_index,domain_index) =
         calculate_interlayer_combination_combination_Jacobian(domain_index, image_index, layers_combination[domain_index]);
      }
   }

   return(interlayers_combination_combination_Jacobian);
}


// Vector< Vector< Vector<double> > > calculate_first_order_forward_propagation(const Vector<double>&) const method


Vector< Vector< Vector<double> > > MultilayerPerceptron::calculate_first_order_forward_propagation(const Vector<double>& inputs) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t inputs_size = inputs.size();

   const size_t inputs_number = get_inputs_number();

   if(inputs_size != inputs_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector< Vector< Vector<double> > > calculate_first_order_forward_propagation(const Vector<double>&) const method.\n"
             << "Size must be equal to number of inputs.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const size_t layers_number = get_layers_number();

   Vector< Vector<double> > layers_combination(layers_number);

   Vector< Vector< Vector<double> > > first_order_forward_propagation(2);

   first_order_forward_propagation[0].set(layers_number);
   first_order_forward_propagation[1].set(layers_number);

   layers_combination[0] = layers[0].calculate_combinations(inputs);

   first_order_forward_propagation[0][0] = layers[0].calculate_activations(layers_combination[0]);

   first_order_forward_propagation[1][0] = layers[0].calculate_activations_derivatives(layers_combination[0]);

   for(size_t i = 1; i < layers_number; i++)
   {
      layers_combination[i] = layers[i].calculate_combinations(first_order_forward_propagation[0][i-1]);

      first_order_forward_propagation[0][i] = layers[i].calculate_activations(layers_combination[i]);

      first_order_forward_propagation[1][i] = layers[i].calculate_activations_derivatives(layers_combination[i]);
   }

   return(first_order_forward_propagation);
}


// Vector< Vector< Vector<double> > > calculate_second_order_forward_propagation(const Vector<double>&) const method


Vector< Vector< Vector<double> > > MultilayerPerceptron::calculate_second_order_forward_propagation(const Vector<double>& inputs) const
{
   // Control sentence (if debug)

   #ifndef NDEBUG 

   const size_t inputs_size = inputs.size();

   const size_t inputs_number = get_inputs_number();

   if(inputs_size != inputs_number)
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: MultilayerPerceptron class.\n"
             << "Vector< Vector< Vector<double> > > calculate_second_order_forward_propagation(const Vector<double>&) const method.\n"
             << "Size of multilayer perceptron inputs must be equal to number of inputs.\n";

      throw std::logic_error(buffer.str());
   }

   #endif

   const size_t layers_number = get_layers_number();

   Vector< Vector<double> > layers_combination(layers_number);

   Vector< Vector< Vector<double> > > second_order_forward_propagation(3);

   second_order_forward_propagation[0].set(layers_number);
   second_order_forward_propagation[1].set(layers_number);
   second_order_forward_propagation[2].set(layers_number);

   layers_combination[0] = layers[0].calculate_combinations(inputs);

   second_order_forward_propagation[0][0] = layers[0].calculate_activations(layers_combination[0]);
  
   second_order_forward_propagation[1][0] = layers[0].calculate_activations_derivatives(layers_combination[0]);

   second_order_forward_propagation[2][0] = layers[0].calculate_activations_second_derivatives(layers_combination[0]);

   for(size_t i = 1; i < layers_number; i++)
   {
      layers_combination[i] = layers[i].calculate_combinations(second_order_forward_propagation[0][i-1]);

      second_order_forward_propagation[0][i] = layers[i].calculate_activations(layers_combination[i]);

      second_order_forward_propagation[1][i] = layers[i].calculate_activations_derivatives(layers_combination[i]);

      second_order_forward_propagation[2][i] = layers[i].calculate_activations_second_derivatives(layers_combination[i]);
   }

   return(second_order_forward_propagation);
}


// std::string to_string(void) const method


std::string MultilayerPerceptron::to_string(void) const
{
   std::ostringstream buffer;

   buffer << "MultilayerPerceptron\n"
          << "Architecture: " << arrange_architecture() << "\n"
          << "Layers activation function: " << write_layers_activation_function() << "\n"
          << "Parameters: " << arrange_parameters() << "\n"
          << "Display: " << display << "\n";

   return(buffer.str());
}


// tinyxml2::XMLDocument* to_XML(void) const method


tinyxml2::XMLDocument* MultilayerPerceptron::to_XML(void) const
{
   std::ostringstream buffer;

   tinyxml2::XMLDocument* document = new tinyxml2::XMLDocument; 

   tinyxml2::XMLElement* multilayer_perceptron_element = document->NewElement("MultilayerPerceptron");

   document->InsertFirstChild(multilayer_perceptron_element);

   // Architecture
   {
      tinyxml2::XMLElement* architecture_element = document->NewElement("Architecture");
      multilayer_perceptron_element->LinkEndChild(architecture_element);

      std::string architecture_string = arrange_architecture().to_string();

      tinyxml2::XMLText* architecture_text = document->NewText(architecture_string.c_str());
      architecture_element->LinkEndChild(architecture_text);
   }

   // Layers activation function
   {
      tinyxml2::XMLElement* layers_activation_function_element = document->NewElement("LayersActivationFunction");
      multilayer_perceptron_element->LinkEndChild(layers_activation_function_element);

      std::string layers_activation_function_string = write_layers_activation_function().to_string();

      tinyxml2::XMLText* layers_activation_function_text = document->NewText(layers_activation_function_string.c_str());
      layers_activation_function_element->LinkEndChild(layers_activation_function_text);
   }

   // Parameters 
   {
      tinyxml2::XMLElement* parameters_element = document->NewElement("Parameters");
      multilayer_perceptron_element->LinkEndChild(parameters_element);

      const std::string parameters_string = arrange_parameters().to_string();

      tinyxml2::XMLText* parameters_text = document->NewText(parameters_string.c_str());
      parameters_element->LinkEndChild(parameters_text);
   }

   // Display
   {
      tinyxml2::XMLElement* display_element = document->NewElement("Display");
      multilayer_perceptron_element->LinkEndChild(display_element);

      buffer.str("");
      buffer << display;

      tinyxml2::XMLText* display_text = document->NewText(buffer.str().c_str());
      display_element->LinkEndChild(display_text);
   }

   return(document);
}


// void from_XML(const tinyxml2::XMLDocument&) method


void MultilayerPerceptron::from_XML(const tinyxml2::XMLDocument& document)
{
   const tinyxml2::XMLElement* root_element = document.FirstChildElement("MultilayerPerceptron");

   if(!root_element)
   {
      return;
   }

      // Architecture
      {
         const tinyxml2::XMLElement* architecture_element = root_element->FirstChildElement("Architecture");

         if(architecture_element)
         {
            const char* architecture_text = architecture_element->GetText();

            if(architecture_text)
            {
               Vector<size_t> new_architecture;
               new_architecture.parse(architecture_text);

               try
               {
                  set(new_architecture);
               }
               catch(const std::logic_error& e)
               {
                  std::cout << e.what() << std::endl;        
               }
            }
         }
      }

      // Layers activation function
      {
         const tinyxml2::XMLElement* layers_activation_function_element = root_element->FirstChildElement("LayersActivationFunction");

         if(layers_activation_function_element)
         {      
            const char* layers_activation_function_text = layers_activation_function_element->GetText();
   
            if(layers_activation_function_text)
            {
               Vector<std::string> new_layers_activation_function;                  
               new_layers_activation_function.parse(layers_activation_function_text);

               try
               {
                  set_layers_activation_function(new_layers_activation_function);
               }
               catch(const std::logic_error& e)
               {
                  std::cout << e.what() << std::endl;        
               }
            }
         }
      }

      // Parameters 
      {
         const tinyxml2::XMLElement* parameters_element = root_element->FirstChildElement("Parameters");

         if(parameters_element)
         {
            const char* parameters_text = parameters_element->GetText();

            if(parameters_text)
            {
               Vector<double> new_parameters;                  
               new_parameters.parse(parameters_text);

               try
               {
                  set_parameters(new_parameters);
               }
               catch(const std::logic_error& e)
               {
                  std::cout << e.what() << std::endl;        
               }
            }
         }
      }

      // Display 
      {         
         const tinyxml2::XMLElement* display_element = root_element->FirstChildElement("Display");

         if(display_element)
         {
            std::string new_display_string = display_element->GetText();     

            try
            {
               set_display(new_display_string != "0");
            }
            catch(const std::logic_error& e)
            {
               std::cout << e.what() << std::endl;       
            }
         }
      }  
}


// Matrix<std::string> write_information(void) const method


Matrix<std::string> MultilayerPerceptron::write_information(void) const
{
    std::ostringstream buffer;

    const size_t layers_number = get_layers_number();

    if(layers_number == 0)
    {
        Matrix<std::string> information;

        return(information);

    }
    else
    {
        Matrix<std::string> information(layers_number, 3);

        for(size_t i = 0; i < layers_number; i++)
        {
            // Inputs number

            buffer.str("");
            buffer << layers[i].get_inputs_number();

            information(i,0) = buffer.str();

            // Perceptrons number

            buffer.str("");
            buffer << layers[i].get_perceptrons_number();

            information(i,1) = buffer.str();

            // Activation function

            information(i,2) = layers[i].write_activation_function();
        }

        return(information);

    }
}


// std::string write_expression(const Vector<std::string>&, const Vector<std::string>&) const method


std::string MultilayerPerceptron::write_expression(const Vector<std::string>& inputs_name, const Vector<std::string>& outputs_name) const
{
   const size_t layers_number = get_layers_number();
   const Vector<size_t> layers_perceptrons_number = arrange_layers_perceptrons_numbers();

   std::ostringstream buffer;

   if(layers_number == 0)
   {
   }
   else if(layers_number == 1)
   {
      buffer << layers[0].write_expression(inputs_name, outputs_name) << "\n";        
   }
   else
   {
      Vector< Vector<std::string> > layers_outputs_name(layers_number);

      for(size_t i = 0; i < layers_number; i++)
      {
         layers_outputs_name[i].set(layers_perceptrons_number[i]);

         for(size_t j = 0; j < layers_perceptrons_number[i]; j++)
         {
            std::ostringstream new_buffer;
            new_buffer << "y_" << i+1 << "_" << j+1;
            layers_outputs_name[i][j] = new_buffer.str();
         }
      }

      buffer << layers[0].write_expression(inputs_name, layers_outputs_name[0]);        
   
      for(size_t i = 1; i < layers_number-1; i++)
      {
         buffer << layers[i].write_expression(layers_outputs_name[i-1], layers_outputs_name[i]);      
      }

      buffer << layers[layers_number-1].write_expression(layers_outputs_name[layers_number-2], outputs_name); 
   }

   return(buffer.str());
}

}

// OpenNN: Open Neural MultilayerPerceptrons Library.
// Copyright (c) 2005-2015 Roberto Lopez.
//
// This library 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 2.1 of the License, or any later version.
//
// This library 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 this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA