//      e04gd Example Program Text
//      C# version, NAG Copyright 2008
using System;
using NagLibrary;
using System.Globalization;
using System.IO;
namespace NagDotNetExamples
{
  public class E04GDE
  {
    static int m = 15;
    static int nt = 3;
    static int ldfjac = m;
    public static double[,] t = new double[m, nt];
    public static double[] y = new double[m];
    static string datafile = "ExampleData/e04gde.d";
    // as a command line argument. It defaults to the named file specified below otherwise.
    //
    static void Main(String[] args)
    {
      if (args.Length == 1)
      {
        datafile = args[0];
      }
      StartExample();
    }
    public static void StartExample()
    {
      E04.E04YA_LSQFUN lsqfunE04YA = new E04.E04YA_LSQFUN(lsqfun);
      E04.E04GD_LSQFUN lsqfunE04GD = new E04.E04GD_LSQFUN(lsqfun);
      E04.E04GD_LSQMON lsqmonE04GD = new E04.E04GD_LSQMON(lsqmon);
      try
      {
        DataReader sr = new DataReader(datafile);
      double eta,  fsumsq,  stepmx,  xtol;
      int i,  ifail,  iprint,  j,  maxcal,  nf,  niter;
      int n = 3;
      int ldv = n;
      double[,] fjac = new double[ldfjac, n];
      double[] fvec = new double[m];
      double[] g = new double[n];
      double[] s = new double[n];
      double[,] v = new double[ldv, n];
      double[] x = new double[n];
      int[] iw = new int[1];
      Console.WriteLine("e04gd Example Program Results");
      //      Skip heading in data file
      sr.Reset();
      //      Observations of tj (j = 1, 2, 3) are held in t[i-1, j-1]
      //      (i = 1, 2, . . . , 15)
      for (i = 1; i <= m; i++)
      {
        sr.Reset();
        y[i - 1] = double.Parse(sr.Next(), CultureInfo.InvariantCulture);
        for (j = 1; j <= nt; j++)
        {
          t[i - 1, j - 1] = double.Parse(sr.Next(), CultureInfo.InvariantCulture);
        }
      }
      //      Check lsqfunE04YA by calling e04ya at an arbitrary point. 
      x[0] = 0.190e0;
      x[1] = -1.340e0;
      x[2] = 0.880e0;
      E04.e04ya(m, n, lsqfunE04YA, x, fvec, fjac, out ifail);
      //
      if (ifail != 0)
      {
        Console.WriteLine("");
        Console.WriteLine(" ** e04ya returned with ifail = {0, 3}", ifail);
        goto L40;
      }
      //      Continue setting parameters for e04gd
      //      Set iprint to 1 to obtain output from lsqfunE04YA at each iteration *
      iprint = -1;
      maxcal = 50 * n;
      eta = 0.90e0;
      xtol = 10.00e0 * Math.Sqrt(X02.x02aj());
      //      We estimate that the minimum will be within 10 units of the
      //      starting point
      stepmx = 10.00e0;
      //      Set up the starting point
      x[0] = 0.50e0;
      x[1] = 1.00e0;
      x[2] = 1.50e0;
      //
      E04.e04gd(m, n, lsqfunE04GD, lsqmonE04GD, iprint, maxcal, eta, xtol, stepmx, x,
      out fsumsq, fvec, fjac, s, v, out niter, out nf, out ifail);
      //
      if (ifail != 0)
      {
        Console.WriteLine("");
        Console.WriteLine(" ** e04gd returned with ifail = {0, 3}", ifail);
        if ((ifail < 0) || (ifail == 1))
        {
          goto L40;
        }
      }
      //
      Console.WriteLine("");
      Console.WriteLine("  {0}{1,12:f4}", "On exit, the sum of squares is", fsumsq);
      Console.Write(" " + " {0}", "at the point");
      for (j = 1; j <= n; j++)
      {
        Console.Write(" " + " {0, 12:f4}", x[j - 1]);
      }
      Console.WriteLine(" ");
      lsqgrd(m, n, fvec, fjac, ldfjac, g);
      Console.Write(" " + " {0}", "The corresponding gradient is");
      for (j = 1; j <= n; j++)
      {
        Console.Write(" " + " {0, 12:e3}", g[j - 1]);
      }
      Console.WriteLine(" ");
      Console.WriteLine(" {0}", "                           (machine dependent)");
      Console.WriteLine(" {0}\n", "and the residuals are");
      for (i = 1; i <= m; i++)
      {
        Console.WriteLine(" " + " {0, 12:e1}", fvec[i - 1]);
      }
      Console.WriteLine(" ");
      //
      L40: ;
      //
    }
    catch (Exception e) {
    Console.WriteLine(e.Message);
    Console.WriteLine( "Exception Raised");
  }
}
//
public static void lsqfun(ref int iflag, int m, int n, double[] xc, double[] fvec,
double[,] fjac)
{
  //      Routine to evaluate the residuals and their 1st derivatives.
  //      An static variable could be updated here to count the
  //      number of calls of lsqfun with iflag set to 1 (since nf
  //      in lsqmon only counts calls with iflag set to 2)
  double denom,  dummy;
  int i = 0;
  for (i = 1; i <= m; i++)
  {
    denom = xc[1] * t[i - 1, 1] + xc[2] * t[i - 1, 2];
    if (iflag == 2)
    {
      fvec[i - 1] = xc[0] + t[i - 1, 0] / denom - y[i - 1];
    }
    fjac[i - 1, 0] = 1.00e0;
    dummy = -1.00e0 / (denom * denom);
    fjac[i - 1, 1] = t[i - 1, 0] * t[i - 1, 1] * dummy;
    fjac[i - 1, 2] = t[i - 1, 0] * t[i - 1, 2] * dummy;
  }
}
//
public static void lsqmon(int m, int n, double[] xc, double[] fvec, double[,] fjac,
double[] s, int igrade, int niter, int nf)
{
  //      Monitoring method
  double fsumsq,  gtg;
  int j = 0;
  double[] g = new double[3];
  int ifail;
  fsumsq = F06.f06ea(m, fvec, 1, fvec, 1, out ifail);
  lsqgrd(m, n, fvec, fjac, ldfjac, g);
  gtg = F06.f06ea(n, g, 1, g, 1, out ifail);
  //      A static variable giving the number of calls of
  //      lsqfun with iflag set to 1 could be printed here
  Console.WriteLine("");
  Console.WriteLine(" {0}", " Itns    F evals          SUMSQ             GTG        grade");
  Console.WriteLine("  {0,4}            {1,5}            {2,13:e5}            {3,9:e1}            {4,3}", niter, nf, fsumsq, gtg, igrade);
  Console.WriteLine("");
  Console.WriteLine(" {0}", "       x                    g           Singular values");
  for (j = 1; j <= n; j++)
  {
    Console.WriteLine("  {0,13:e5}                    {1,9:e1}                    {2,9:e1}", xc[j - 1], g[j - 1], s[j - 1]);
  }
  //
}
//
public static void lsqgrd(int m, int n, double[] fvec, double[,] fjac, int ldfjac,
double[] g)
{
  //      Routine to evaluate gradient of the sum of squares
  double sum = 0.0;
  int i,  j;
  for (j = 1; j <= n; j++)
  {
    sum = 0.00e0;
    for (i = 1; i <= m; i++)
    {
      sum = sum + fjac[i - 1, j - 1] * fvec[i - 1];
    }
    g[j - 1] = sum + sum;
  }
}
}
}