using System;
using System.Diagnostics;

// =======================================================================================================================
// Class name: TDLevenbergMarquardtFitting
// 
// Description:
// 
// Copyright (c) 2010-2011 Deltares
// 
// Date        ID   Modification
// 2010-02-03  Sel  Created
// 2010-02-05  Sel  Moved Variables to BasicExtremeAndFittingCalculation
// 2010-02-10  Sel  Improved LevenbergMarquardt algorithm code
// 2010-02-10  Sel  Introduced MDVariationalCalculus
// =======================================================================================================================

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

        public IFittingCalculation LevenbergMarquardtCalculation
        {
            get
            {
                return levenbergMarquardtCalculation;
            }
            set
            {
                levenbergMarquardtCalculation = value;
            }
        }

        public override void DetermineDrift(double[] parameters)
        {
            // Determine function with Drift and store result in vector
            driftAfter = 0;
            // Undeclared identifier: 'AParameters'
            functionValues = levenbergMarquardtCalculation.ComputeFittingResult(parameters);
            for (int variableIndex = 0; variableIndex <= VariableCount - 1; variableIndex++)
            {
                driftAfter = driftAfter + Math.Pow(Convert.ToDouble(functionValues[variableIndex]), 2);
            }
            driftAfter = Math.Sqrt(driftAfter/VariableCount);
        }

        public override void FillJacobian()
        {
            TMatrix LJacobian = null;
            Debug.Assert((levenbergMarquardtCalculation != null));
            // Allocate Non-square Jacobian
            LJacobian = new TMatrix();
            LJacobian.SetDimensions(VariableCount, ParameterCount);

            // Determine drift and copy function value
            DetermineDrift(parameters);
            for (int variableIndex = 0; variableIndex <= VariableCount - 1; variableIndex++)
            {
                LJacobian.Vector[variableIndex] = functionValues[variableIndex];
            }

            for (int parameterIndex = 0; parameterIndex <= ParameterCount - 1; parameterIndex++)
            {
                // Determine variation of function value and store result in matrix
                // Apply a small variation on the parameter concerned
                double parameterValue = parameters[parameterIndex];
                parameters[parameterIndex] = parameterValue + FShift;
                // Determine derivative and store variation in matrix
                functionValues = levenbergMarquardtCalculation.ComputeFittingResult(parameters);
                for (int variableIndex = 0; variableIndex <= VariableCount - 1; variableIndex++)
                {
                    LJacobian.Matrix[variableIndex][parameterIndex] = (functionValues[variableIndex] - LJacobian.Vector[variableIndex])/FShift;
                }
                // Restore parameter
                parameters[parameterIndex] = parameterValue;
            }
            // Quadrature of LJacobian yields square FJacobian
            jacobian = LJacobian.Quadrature();
        }
    }
}