/* G02DA_P0W_F C++ Header Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 * Mark 30.2, 2024.
 */

#include <iostream>
#include <dco.hpp>
#include <nag.h>
#include <nagad.h>
#include <stdio.h>
#include <string>
using namespace std;

int main()
{
  int    exit_status = 0;
  double tol         = 0.000001;

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

  // 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.
  double * x = 0, *y = 0;
  Integer *isx = 0;
  x            = new double[m * n];
  y            = new double[n];
  isx          = new Integer[m];

  // Read model data
  for (int i = 0; i < n; ++i)
  {
    for (int j = 0; j < m; ++j)
    {
      int k = i + j * n;
      cin >> x[k];
    }
    cin >> y[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
  double *b = 0, *cov = 0, *h = 0, *p = 0, *q = 0, *res = 0, *se = 0, *wk = 0;
  Integer lcov = (ip * ip + ip) / 2, lp = ip * (ip + 2), lq = n * (ip + 1);
  Integer lwk = ip * ip + 5 * (ip - 1);
  b           = new double[ip];
  cov         = new double[lcov];
  h           = new double[n];
  p           = new double[lp];
  q           = new double[lq];
  res         = new double[n];
  se          = new double[ip];
  wk          = new double[lwk];

  // Perform Regression
  Integer           idf, irank;
  double            rss, wt[1];
  logical           svd;
  nag::ad::handle_t ad_handle;
  Integer           ifail = 0;

  nag::ad::g02da(ad_handle, "M", "U", n, x, n, m, isx, ip, y, wt, 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         = " << 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 << b[i] << "       ";
    cout.width(12);
    cout << se[i] << endl;
  }

  delete[] x;
  delete[] y;
  delete[] isx;
  delete[] b;
  delete[] cov;
  delete[] h;
  delete[] p;
  delete[] q;
  delete[] res;
  delete[] se;
  delete[] wk;

  return exit_status;
}