/// <summary>
///   Calculates partial derivatives for all weights of the network.
/// </summary>
///
/// <param name="input">The input vector.</param>
/// <param name="desiredOutput">Desired output vector.</param>
/// <param name="outputIndex">The current output location (index) in the desired output vector.</param>
///
/// <returns>Returns summary squared error of the last layer.</returns>
///
private double CalculateDerivatives(double[] input, double[] desiredOutput, int outputIndex)
{
    // assume, that all neurons of the network have the same activation function
    IActivationFunction function = network[0][0].ActivationFunction;

    double[] previousLayerOutput;

    // Start by the output layer first
    int outputLayerIndex = network.LayersCount - 1;
    ActivationLayer outputLayer = network[outputLayerIndex];

    // If we have only one single layer, the previous layer outputs is given by the input layer
    previousLayerOutput = (outputLayerIndex == 0) ? input : network[outputLayerIndex - 1].Output;

    // Assume single output neuron
    ActivationNeuron outputNeuron = outputLayer[outputIndex];
    double[] neuronWeightDerivatives = weightDerivatives[outputLayerIndex][outputIndex];

    double output = outputNeuron.Output;
    double e = desiredOutput[outputIndex] - output;
    double derivative = function.Derivative2(output);

    // Set derivative for each weight in the neuron
    for (int i = 0; i < previousLayerOutput.Length; i++)
        neuronWeightDerivatives[i] = derivative * previousLayerOutput[i];

    // Set derivative for the current threshold (bias) term
    thresholdsDerivatives[outputLayerIndex][outputIndex] = derivative;


    // Now, proceed to the hidden layers
    for (int layerIndex = network.LayersCount - 2; layerIndex >= 0; layerIndex--)
    {
        int nextLayerIndex = layerIndex + 1;

        ActivationLayer layer = network[layerIndex];
        ActivationLayer nextLayer = network[nextLayerIndex];

        // If we are in the first layer, the previous layer is just the input layer
        previousLayerOutput = (layerIndex == 0) ? input : network[layerIndex - 1].Output;

        // Now, we will compute the derivatives for the current layer applying the chain
        //  rule. To apply the chain-rule, we will make use of the previous derivatives
        //  computed for the inner layers (forming a calculation chain, hence the name).

        // So, for each neuron in the current layer:
        for (int neuronIndex = 0; neuronIndex < layer.NeuronsCount; neuronIndex++)
        {
            ActivationNeuron neuron = layer[neuronIndex];

            neuronWeightDerivatives = weightDerivatives[layerIndex][neuronIndex];

            double[] layerDerivatives = thresholdsDerivatives[layerIndex];
            double[] nextLayerDerivatives = thresholdsDerivatives[layerIndex + 1];

            double sum = 0;

            // The chain-rule can be stated as (f(w*g(x))' = f'(w*g(x)) * w*g'(x)
            //
            // We will start computing the second part of the product. Since the g'
            //  derivatives have already been computed in the previous computation,
            //  we will be summing all previous function derivatives and weighting
            //  them using their connection weight (sinapses).
            //
            // So, for each neuron in the next layer:
            for (int j = 0; j < nextLayerDerivatives.Length; j++)
            {
                // retrieve the weight connecting the output of the current
                //   neuron and the activation function of the next neuron.
                double weight = nextLayer[j][neuronIndex];

                // accumulate the sinapse weight * next layer derivative
                sum += weight * nextLayerDerivatives[j];
            }

            // Continue forming the chain-rule statement
            derivative = sum * function.Derivative2(neuron.Output);

            // Set derivative for each weight in the neuron
            for (int i = 0; i < previousLayerOutput.Length; i++)
                neuronWeightDerivatives[i] = derivative * previousLayerOutput[i];

            // Set derivative for the current threshold
            layerDerivatives[neuronIndex] = derivative;

            // The threshold derivatives also gather the derivatives for
            // the layer, and thus can be re-used in next calculations.
        }
    }

    // return error
    return e;
}