using System;
using System.Text;
using System.Runtime.InteropServices;  /* Provides mappings between C# and native code */

namespace NumericalAlgorithmsGroup
{
    class MultipleLinearRegression
    {
      // Declaration of the routine G02DAF in the NAG DLL
      // Points to note:
      //   (i) the character string argument "mean" requires an extra "invisible"
      //       argument "length_mean" which is passed by value and must be set to the
      //       length of the string argument. The same goes for string argument "weight".
      //       These "invisible arguments" are inserted by the Fortran compiler
      //       and are not required when compiling a Fortran program, which is
      //       why they are not documented in the NAG Fortran Library documentation.
      //   (ii) all scalar arguments other than the string length arguments are
      //        passed by reference, as required by Fortran
      [DllImport("FLDLL214Z_nag.dll")]/* Import from DLL, the C# compiler provides a rudimentry check of the signature */
        public static extern void G02DAF(string mean, int length_mean,
                                         string weight, int length_weight,
                                         ref int n, double[,] x, ref int ldx,
                                         ref int m, int[] isx, ref int ip,
                                         double[] y, double[] wt, ref double rss,
                                         ref int idf, double[] b, double[] se,
                                         double[] cov, double[] res, double[] h,
                                         double[,] q, ref int ldq, ref int svd,
                                         ref int irank, double[] p, ref double tol,
                                         double[] work, ref int ifail);

      public static void Main()
        {
          // We solve the example problem presented in the documentation of
          // routine G02DAF in the NAG Fortran Library Manual.
          int n = 12;
          int m = 4;

          // Notice that array X is initialised in transposed format because
          // the NAG Fortran library has arrays in column-major order
      double [,] x = new double [,] {
	{ 1.0, 0.0, 0.0, 0.0 },
	{ 0.0, 0.0, 0.0, 1.0 },
	{ 0.0, 1.0, 0.0, 0.0 },
	{ 0.0, 0.0, 1.0, 0.0 },
	{ 0.0, 0.0, 0.0, 1.0 },
	{ 0.0, 1.0, 0.0, 0.0 },
	{ 0.0, 0.0, 0.0, 1.0 },
	{ 1.0, 0.0, 0.0, 0.0 },
	{ 0.0, 0.0, 1.0, 0.0 },
	{ 1.0, 0.0, 0.0, 0.0 },
	{ 0.0, 0.0, 1.0, 0.0 },
	{ 0.0, 1.0, 0.0, 0.0 }
      };
     
      // Notice that array X is transposed to XT transposed because
      // the NAG Fortran library has arrays in column-major order
      double [,] xt = new double [4,12];
     for (int i=0; i<n; i++)
       for (int j=0; j<m; j++)
	 xt[j,i]=x[i,j];

          double[] y = new double[] {33.63, 39.62, 38.18, 41.46, 38.02,
                                     35.83, 35.99, 36.58, 42.92, 37.80,
                                     40.43, 37.89};
          int[] sx = new int[m];
          for (int i = 0; i < m; i++)
            sx[i] = 1;

          int ip = m + 1;
          double[] b = new double[ip];
          double[] se = new double[ip];
          double[] cov = new double[ip * (ip + 1) / 2];
          double[] res = new double[n];
          double[] h = new double[n];
          double[,] q = new double[ip + 1, n];
          double[] p = new double[2 * ip + ip * ip];
          double[] com_ar = new double[5 * (ip - 1) + ip * ip];

          int svd = 0;
          int rank = 0;
          int ifail = 1;

          double[] wt = new double[n];
          for (int i = 0; i < n; i++)
            wt[i] = 1.0;

          double rss = 0.0, tol = 0.00001;
          int idf = 0;

          com_ar[0] = 33.0;
          q[0, 0] = 33.0;
          G02DAF("m", 1, "w", 1, ref n, xt, ref n, ref m, sx, ref ip, y, wt,
                 ref rss, ref idf, b, se, cov, res, h, q, ref n, ref svd,
                 ref rank, p, ref tol, com_ar, ref ifail);

          if (svd != 0)
            Console.WriteLine("Model not of full rank, rank = {0}", rank);
          Console.WriteLine("Residual sum of squares = {0,12:e4}", rss);
          Console.WriteLine("Degrees of freedom = {0}", idf);
          Console.WriteLine("  Variable    Parameter estimate     Standard error");
          for (int j = 0; j < ip; j++)
            Console.WriteLine("   {0, 4} {1, 20:e4} {2, 20:e4}", j + 1, b[j], se[j]);
          Console.WriteLine();
          Console.WriteLine("     Obs           Residuals                h");
          for (int i = 0; i < n; i++)
            Console.WriteLine("   {0, 4} {1, 20:e4} {2, 20:e4}", i + 1, res[i], h[i]);
        }
    }
}