using System;
using System.Diagnostics;
using Deltares.Mathematics;

// =======================================================================================================================
// Class name: TDLevenbergMarquardtExtreme
//
// Description:
//
// Copyright (c) 2010 Deltares
//
// Date        ID   Modification
// 2010-02-09  Sel  Created
// 2010-02-10  Sel  Improved LevenbergMarquardt algorithm code
// 2010-02-10  Sel  Introduced MDVariationalCalculus
// 2010-02-15  Sel  Allocated Bounds in LM; Replaced AllocateMemory by SetVariableCount.
// 2010-05-20  Sel  Finetuned use of GA in general

namespace Deltares.Stability.Calculation2.Levenberg
{
    // GeoLib units
    public class LevenbergMarquardtExtreme : LevenbergMarquardt
    {
        // Private fields
        private IExtremeCalculation levenbergMarquardtCalculation = null;

        // Public properties
        public IExtremeCalculation LevenbergMarquardtCalculation
        {
            get
            {
                return levenbergMarquardtCalculation;
            }
            set
            {
                levenbergMarquardtCalculation = value;
            }
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2010-02-09  Sel  Created
        // =======================================================================================================================
        public override void SetVariableCount(int count)
        {
            base.SetVariableCount(count);
            functionValues = new double[count];
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2010-02-09  Sel  Created
        // =======================================================================================================================
        public override void DetermineDrift(double[] parameters)
        {
            // Determine derivative with Drift and store result in vector
            driftAfter = 0;
            DetermineDerivative(parameters);
            for (int variableIndex = 0; variableIndex < VariableCount; variableIndex ++)
            {
                driftAfter = driftAfter + Math.Pow(Convert.ToDouble(functionValues[variableIndex]), 2);
            }
            driftAfter = Math.Sqrt(driftAfter/VariableCount);

            // The drift never will vanish. If so, ComputeExtremeResult failed and the result is improper and must be discarded.
            if (driftAfter == 0)
            {
                driftAfter = 2*driftPrior;
            }
        }

        // =======================================================================================================================
        // Date        ID   Modification
        // 2010-01-19  Sel  Created
        // =======================================================================================================================
        public override void FillJacobian()
        {
            // This extreme determination is square
            Debug.Assert(ParameterCount == VariableCount);

            // Determine drift and copy derivative
            DetermineDrift(parameters);
            for (int variableIndex = 0; variableIndex < VariableCount; variableIndex ++)
            {
                jacobian.Vector[variableIndex] = functionValues[variableIndex];
            }

            // Determine variation of derivative and store result in matrix
            for (int parameterIndex = 0; parameterIndex < ParameterCount; parameterIndex ++)
            {
                // Apply a small variation on the parameter concerned
                double parameterValue = parameters[parameterIndex];
                parameters[parameterIndex] = parameterValue + FShift;

                // Determine derivative and store variation in matrix
                DetermineDerivative(parameters);

                for (int variableIndex = 0; variableIndex < VariableCount; variableIndex ++)
                {
                    jacobian.Matrix[variableIndex][parameterIndex] = (functionValues[variableIndex] - jacobian.Vector[variableIndex])/FShift;
                }

                // Restore parameter
                parameters[parameterIndex] = parameterValue;
            }

            // The numerical derivative is not precise. For a guess of the drift it is all right, but for the determination
            // of the extreme a correction of half the shift times the second derivative will improve the outcome.
            for (int variableIndex = 0; variableIndex < VariableCount; variableIndex ++)
            {
                jacobian.Vector[variableIndex] = jacobian.Vector[variableIndex] - FShift/2*jacobian.Matrix[variableIndex][variableIndex];
            }
        }

        /// <summary>
        /// Determines the derivative
        /// </summary>
        /// <param name="parameters">Parameters</param>
        /// <returns>Indication whether the underlying calculations were successful</returns>
        private void DetermineDerivative(double[] parameters)
        {
            Debug.Assert((levenbergMarquardtCalculation != null));

            // Determine function values and store the numeric derivative in vector
            double functionValue = levenbergMarquardtCalculation.ComputeExtremeResult(parameters);

            if (!IsValidResult(functionValue))
            {
                valid = false;
            }

            for (int parameterIndex = 0; parameterIndex < ParameterCount; parameterIndex ++)
            {
                // Apply a small variation, determine derivative and restore parameter
                double parameterValue = parameters[parameterIndex];
                parameters[parameterIndex] = parameterValue + FShift;
                double shiftedResult = levenbergMarquardtCalculation.ComputeExtremeResult(parameters);

                if (!IsValidResult(shiftedResult))
                {
                    valid = false;
                }

                functionValues[parameterIndex] = (levenbergMarquardtCalculation.ComputeExtremeResult(parameters) - functionValue)/FShift;
                parameters[parameterIndex] = parameterValue;
            }
        }
    }
}