/* nag_tsa_multi_inp_model_estim (g13bec) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */

#include <nag.h>
#include <stdio.h>
#include <string.h>
#include <nag_string.h>
#include <nag_stdlib.h>
#include <nagg13.h>

#define XXY(I, J) xxy[(I) *tdxxy + J]

int main(void)
{
  Integer exit_status = 0;
  Integer i, inser, j, npara, nseries, nxxy, tdxxy;
  Nag_ArimaOrder arimav;
  Nag_G13_Opt options;
  Nag_TransfOrder transfv;
  double df, objf, *para = 0, rss, *sd = 0, *xxy = 0;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_tsa_multi_inp_model_estim (g13bec) Example Program Results\n");
  scanf(" %*[^\n]"); /* Skip heading in data file */

#define CM(I, J) options.cm[(J)+(I) *options.tdcm]
#define ZT(I, J) options.zt[(J)+(I) *options.tdzt]

  /*
   * Initialize the option structure.
   */
  /* nag_tsa_options_init (g13bxc).
   * Initialization function for option setting
   */
  nag_tsa_options_init(&options);
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "", &nxxy, &nseries,
        &options.max_iter);
  if (nxxy > 0 && nseries > 0) {
    /*
     * Set some specific option variables to the desired values.
     */

    options.criteria = Nag_Marginal;
    options.print_level = Nag_Soln_Iter_Full;
    /*
     * Allocate memory to the arrays in structure transfv containing
     * the transfer function model orders of the input series.
     */
    /* nag_tsa_transf_orders (g13byc), see above. */
    nag_tsa_transf_orders(nseries, &transfv, &fail);

    /*
     * Read the orders vector of the ARIMA model for the output noise
     * component into structure arimav.
     */
    scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT ""
          "%" NAG_IFMT "%" NAG_IFMT "", &arimav.p, &arimav.d, &arimav.q,
          &arimav.bigp, &arimav.bigd, &arimav.bigq, &arimav.s);
    /*
     * Read the transfer function model orders of the input series into
     * structure transfv.
     */
    inser = nseries - 1;

    for (j = 0; j < inser; ++j)
      scanf("%" NAG_IFMT "", &transfv.b[j]);
    for (j = 0; j < inser; ++j)
      scanf("%" NAG_IFMT "", &transfv.q[j]);
    for (j = 0; j < inser; ++j)
      scanf("%" NAG_IFMT "", &transfv.p[j]);
    for (j = 0; j < inser; ++j)
      scanf("%" NAG_IFMT "", &transfv.r[j]);

    npara = 0;
    for (i = 0; i < inser; ++i)
      npara = npara + transfv.q[i] + transfv.p[i];
    npara = npara + arimav.p + arimav.q + arimav.bigp + arimav.bigq + nseries;
    if (npara >= 1) {
      if (!(para = NAG_ALLOC(npara, double)) ||
          !(sd = NAG_ALLOC(npara, double)) ||
          !(xxy = NAG_ALLOC(nxxy * nseries, double)))
      {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
      tdxxy = nseries;
      for (i = 0; i < npara; ++i)
        scanf("%lf", &para[i]);
      for (i = 0; i < nxxy; ++i)
        for (j = 0; j < nseries; ++j)
          scanf("%lf", &XXY(i, j));

      /* nag_tsa_multi_inp_model_estim (g13bec), see above. */
      fflush(stdout);
      nag_tsa_multi_inp_model_estim(&arimav, nseries, &transfv, para,
                                    npara, nxxy, xxy, tdxxy, sd, &rss,
                                    &objf, &df, &options, &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_tsa_multi_inp_model_estim (g13bec)"
               ".\n%s\n", fail.message);
        exit_status = 1;
        goto END;
      }

      printf("\nThe correlation matrix is \n\n");
      for (i = 0; i < npara; ++i)
        for (j = 0; j < npara; ++j)
          printf("%10.4f%c", CM(i, j), (j % 5 == 4) ? '\n' : ' ');
      printf("\nThe residuals and the z and n values are\n\n");
      printf("   i         res[i]          z(t)        noise(t)\n\n");
      for (i = 0; i < nxxy; ++i) {
        if (i + 1 <= options.lenres) {
          printf("%4" NAG_IFMT "%15.3f", i + 1, options.res[i]);
          for (j = 0; j < nseries - 1; ++j)
            printf("%15.3f ", ZT(i, j));
          printf("%15.3f\n", options.noise[i]);
        }
      }
    }
    else {
      printf("npara is out of range: npara = %-3" NAG_IFMT "\n", npara);
      /* nag_tsa_free (g13xzc).
       * Freeing function for use with g13 option setting
       */
      nag_tsa_free(&options);
      /* nag_tsa_trans_free (g13bzc), see above. */
      nag_tsa_trans_free(&transfv);
      exit_status = 1;
      goto END;
    }
  }
  else {
    printf("One or both of nxxy and nseries are out of range:"
           " nxxy = %-3" NAG_IFMT "  while  nseries = %-3" NAG_IFMT "\n",
           nxxy, nseries);
    exit_status = 1;
    goto END;
  }
  /* nag_tsa_trans_free (g13bzc), see above. */
  nag_tsa_trans_free(&transfv);
  /* nag_tsa_free (g13xzc), see above. */
  nag_tsa_free(&options);

END:
  NAG_FREE(para);
  NAG_FREE(sd);
  NAG_FREE(xxy);

  return exit_status;
}