/* G02DA_T1W_F C++ Header Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 * Mark 28.3, 2022.
 */

#include <dco.hpp>
#include <iostream>
#include <math.h>
#include <nag.h>
#include <nagad.h>
#include <nagx04.h>
#include <stdio.h>
#include <string>

typedef double                   DCO_BASE_TYPE;
typedef dco::gt1s<DCO_BASE_TYPE> DCO_MODE;
typedef DCO_MODE::type           DCO_TYPE;

using namespace std;

int main()
{
  int               exit_status = 0;
  nag::ad::handle_t ad_handle;
  DCO_TYPE          tol;
  Integer           ifail = 0;

  cout << "G02DA_T1W_F C++ Header Example Program Results\n\n";

  tol = 0.000001;

  // Skip heading in data file
  string mystr;
  getline(cin, mystr);

  // Read problem sizes
  Integer n, m;
  cin >> n;
  cin >> m;

  // Allocate arrays depending on m and n.
  DCO_TYPE *x = 0, *y = 0, *y_in = 0;
  Integer * isx = 0;
  double *  dy  = 0;
  x             = new DCO_TYPE[m * n];
  y             = new DCO_TYPE[n];
  y_in          = new DCO_TYPE[n];
  dy            = new double[n];
  isx           = new Integer[m];

  // Create AD configuration data object
  ifail = 0;

  // Read model data
  double dd;
  for (int i = 0; i < n; ++i)
  {
    for (int j = 0; j < m; ++j)
    {
      cin >> dd;
      int k = i + j * n;
      x[k]  = dd;
    }
    cin >> dd;
    y_in[i] = dd;
    y[i]    = y_in[i];
  }

  // Calculate ip
  Integer ip = 0;
  for (int j = 0; j < m; ++j)
  {
    cin >> isx[j];
    if (isx[j] > 0)
    {
      ip++;
    }
  }
  // Mean = 'M'
  ip++;

  // Allocate arrays depending on ip
  DCO_TYPE *b = 0, *cov = 0, *h = 0, *p = 0, *q = 0, *res = 0, *se = 0, *wk = 0;
  double *  dbdy;
  Integer   lcov = (ip * ip + ip) / 2, lp = ip * (ip + 2), lq = n * (ip + 1);
  Integer   lwk = ip * ip + 5 * (ip - 1);
  b             = new DCO_TYPE[ip];
  cov           = new DCO_TYPE[lcov];
  h             = new DCO_TYPE[n];
  p             = new DCO_TYPE[lp];
  q             = new DCO_TYPE[lq];
  res           = new DCO_TYPE[n];
  se            = new DCO_TYPE[ip];
  wk            = new DCO_TYPE[lwk];
  dbdy          = new double[n * ip];

  // Perform Regression
  Integer  idf, irank;
  DCO_TYPE rss;
  logical  svd;
  ifail = 0;
  nag::ad::g02da(ad_handle, "M", "U", n, x, n, m, isx, ip, y, NULL, rss, idf, b,
                 se, cov, res, h, q, n, svd, irank, p, tol, wk, ifail);

  // Display results
  if (svd)
  {
    cout << "Model is not of full rank, rank = " << irank << endl;
  }
  cout << "Residual sum of squares         = " << dco::value(rss);
  cout << "\nDegrees of freedom              = " << idf << endl;
  cout << "\nVariable   Parameter estimate   Standard error\n\n";
  cout.setf(ios::scientific, ios::floatfield);
  cout.precision(2);
  for (int i = 0; i < ip; ++i)
  {
    cout.width(5);
    cout << i << "       ";
    cout.width(12);
    cout << dco::value(b[i]) << "       ";
    cout.width(12);
    cout << dco::value(se[i]) << endl;
  }

  cout << "\n\n Derivatives calculated: First order tangents\n";
  cout << " Computational mode    : algorithmic\n";
  cout << "\n Derivatives:\n\n";

  // Obtain derivatives
  for (int i = 0; i < n; i++)
  {
    y[i] = y_in[i];
  }
  for (int i = 0; i < n; i++)
  {
    dco::derivative(y[i]) = 1.0;
    nag::ad::g02da(ad_handle, "M", "U", n, x, n, m, isx, ip, y, NULL, rss, idf,
                   b, se, cov, res, h, q, n, svd, irank, p, tol, wk, ifail);
    dy[i] = dco::derivative(rss);

    for (int j = 0; j < ip; j++)
    {
      int k   = i + j * n;
      dbdy[k] = dco::derivative(b[j]);
    }
    for (int j = 0; j < n; ++j)
    {
      y[j] = y_in[j];
    }
  }
  cout << "    i        d(rss)/dy(i)\n";
  cout.precision(4);
  for (int i = 0; i < n; ++i)
  {
    cout.width(5);
    cout << i << "       ";
    cout.width(12);
    cout << dy[i] << endl;
  }

  // Print matrix routine
  cout << endl;

  NagError fail;
  INIT_FAIL(fail);
  x04cac(Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, n, ip, dbdy, n,
         "db/dy", 0, &fail);

  ifail = 0;

  delete[] b;
  delete[] cov;
  delete[] h;
  delete[] p;
  delete[] q;
  delete[] res;
  delete[] se;
  delete[] wk;
  delete[] dbdy;
  delete[] x;
  delete[] y;
  delete[] y_in;
  delete[] dy;
  delete[] isx;
  return exit_status;
}