using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using alglib;


namespace ConsoleApplication1
{
    public class _Demo
    {
        public static int Main(string[] args)
        {
            int m = 0;
            int n = 0;
            int k = 0;
            double[] y = new double[0];
            double[,] x = new double[0, 0];
            double[] c = new double[0];
            lsfit.lsfitreport rep = new lsfit.lsfitreport();
            lsfit.lsfitstate state = new lsfit.lsfitstate();
            int info = 0;
            double epsf = 0;
            double epsx = 0;
            int maxits = 0;
            int i = 0;
            System.Console.Write("Fitting f = d-ae(-(x-b)^2/c^2)");
            System.Console.WriteLine();

            //
            // * using Hessian
            // * using alpha=1 and beta=0 as initial values
            // * using 1000 uniformly distributed points to fit to
            //
            // Notes:
            // * N - number of points
            // * M - dimension of space where points reside
            // * K - number of parameters being fitted
            //
            n = 17;
            m = 1;
            k = 4;
            y = new double[] {-0.013622, 0.012748, 0.071649, 0.128783, 0.139753, 0.126173, 0.080030
                , 0.036288, 0.017658, 0.011651, 0.001684, 0.002987, -0.000251, -0.000626, -0.005298, -0.010551, -0.013150};

            //
            // Prepare task matrix
            //
            x = new double[n, m];
            c = new double[k];
            for (i = 0; i <= n - 1; i++)
            {
                x[i, 0] = i;
            }
            // Initial Guess 
            c[0] = -0.150688;
            c[1] = 4.13677;
            c[2] = 2.44339;
            c[3] = -0.0037984;

            epsf = 0.0000001;
            epsx = 1e-8;
            maxits = 0;

            //
            // Solve
            //
            lsfit.lsfitnonlinearfgh(ref x, ref y, ref c, n, m, k, ref state);
            lsfit.lsfitnonlinearsetcond(ref state, epsf, epsx, maxits);
            while (lsfit.lsfitnonlineariteration(ref state))
            {

                //
                // F(x)      = d - a*exp(-((x-b)/c)^2); 
                //
                state.f = state.c[3] - state.c[0] / Math.Exp(AP.Math.Sqr(state.x[0] - state.c[1]) /AP.Math.Sqr(state.c[2]));
                               //
                // F(x)      = d - a*exp(-((x-b)/c)^2);
                // dF/da = -1/(exp((x-b)^2/c^2);
                // dF/db  = -2*a(x-b)/(c^2*exp((x-b)^2/c^2);
                // dF/dc = -(2*a*(x-b)^2/(c^3*exp((x-b)^2/c^2));
                // dF/dd = 1;
                //
                if (state.needfg | state.needfgh)
                {
                    state.g[0] = -1 / Math.Exp(AP.Math.Sqr(state.x[0] - state.c[1])/AP.Math.Sqr(state.c[2]));
                    state.g[1] = 2 * state.c[0] * (state.c[1] -state.x[0]) / AP.Math.Sqr(state.c[2]) * Math.Exp(AP.Math.Sqr(state.c[1] - state.x[0]) / AP.Math.Sqr(state.c[2]));
                    state.g[2] = -2 * state.c[0] * AP.Math.Sqr(state.c[1] - state.x[0]) / Math.Pow(state.c[2], 3) * Math.Exp(AP.Math.Sqr(state.c[1] - state.x[0]) / AP.Math.Sqr(state.c[2]));
                    state.g[3] = 1.0;
                }
                if (state.needfgh)
                {
                    state.h[0, 0] = 0.0;
                    state.h[0, 1] = 2 * (state.c[1] - state.x[0]) / AP.Math.Sqr(state.c[2]) * Math.Exp(AP.Math.Sqr((state.c[1] - state.x[0]) / state.c[2]));
                    state.h[0, 2] = -2 * AP.Math.Sqr(state.c[1] - state.x[0]) / Math.Pow(state.c[2], 3) * Math.Exp(AP.Math.Sqr((state.c[1] - state.x[0]) / state.c[2]));
                    state.h[0, 3] = 0.0;

                    state.h[1, 0] = 2 * (state.c[1] - state.x[0]) / AP.Math.Sqr(state.c[2]) * Math.Exp(AP.Math.Sqr((state.c[1] - state.x[0]) / state.c[2]));
                    state.h[1, 1] = 2 * state.c[0] / AP.Math.Sqr(state.c[2]) * Math.Exp(AP.Math.Sqr((state.c[1] - state.x[0]) / state.c[2])) - 4 * state.c[0] * AP.Math.Sqr(state.c[1] - state.x[0]) / Math.Pow(state.c[2], 4) * Math.Exp(AP.Math.Sqr((state.c[1] - state.x[0]) / state.c[2]));
                    state.h[1, 2] = -4 * state.c[0] * (state.c[1] - state.x[0]) / Math.Pow(state.c[2], 3) * Math.Exp(AP.Math.Sqr((state.c[1] - state.x[0]) / state.c[2])) + 4 * state.c[0] * Math.Pow((state.c[1] - state.x[0]), 3) / Math.Pow(state.c[2], 5) * Math.Exp(AP.Math.Sqr((state.c[1] - state.x[0]) / state.c[2]));
                    state.h[1, 3] = 0.0;

                    state.h[2, 0] = -2 * AP.Math.Sqr(state.c[1] - state.x[0]) / Math.Pow(state.c[2], 3) * Math.Exp(AP.Math.Sqr((state.c[1] - state.x[0]) / state.c[2]));
                    state.h[2, 1] = -4 * state.c[0] * (state.c[1] - state.x[0]) / Math.Pow(state.c[2], 3) * Math.Exp(AP.Math.Sqr((state.c[1] - state.x[0]) / state.c[2])) + 4 * state.c[0] * Math.Pow((state.c[1] - state.x[0]), 3) / Math.Pow(state.c[2], 5) * Math.Exp(AP.Math.Sqr((state.c[1] - state.x[0]) / state.c[2]));
                    state.h[2, 2] = 6 * state.c[0] * AP.Math.Sqr(state.c[1] - state.x[0]) / Math.Pow(state.c[2], 4) * Math.Exp(AP.Math.Sqr((state.c[1] - state.x[0]) / state.c[2])) - 4 * state.c[0] * Math.Pow((state.c[1] - state.x[0]), 4) / Math.Pow(state.c[2], 6) * Math.Exp(AP.Math.Sqr((state.c[1] - state.x[0]) / state.c[2]));
                    state.h[2, 3] = 0.0;

                    state.h[3, 0] = 0.0;
                    state.h[3, 1] = 0.0;
                    state.h[3, 2] = 0.0;
                    state.h[3, 3] = 0.0;
                }
            }
            lsfit.lsfitnonlinearresults(ref state, ref info, ref c, ref rep);
            System.Console.Write("a:   ");
            System.Console.Write("{0,0:F5}", c[0]);
            System.Console.WriteLine();
            System.Console.Write("b:    ");
            System.Console.Write("{0,0:F5}", c[1]);
            System.Console.WriteLine();
            System.Console.Write("c:    ");
            System.Console.Write("{0,0:F5}", c[2]);
            System.Console.WriteLine();
            System.Console.Write("d:    ");
            System.Console.Write("{0,0:F5}", c[3]);
            System.Console.WriteLine();
            System.Console.Write("rms.err: ");
            System.Console.Write("{0,0:F3}", rep.rmserror);
            System.Console.WriteLine();
            System.Console.Write("max.err: ");
            System.Console.Write("{0,0:F3}", rep.maxerror);
            System.Console.WriteLine();
            System.Console.Write("Termination type: ");
            System.Console.Write("{0,0:d}", info);
            System.Console.WriteLine();
            System.Console.WriteLine();
            System.Console.WriteLine();
            return 0;
        }
    }
}