/* nag_tsa_multi_inputmod_forecast (g13bjc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.2, 2023.
 */

#include <nag.h>
#include <stdio.h>
#include <string.h>

#define PARX(I, J) parx[(I)*tdparx + J]
#define XXY(I, J) xxy[(I)*tdxxy + J]
#define MRX(I, J) mrx[(I)*tdmrx + J]

int main(void) {
  Integer exit_status = 0;
  Integer i, inser, j, ldparx, *mrx = 0, n, nev, nfv, npara;
  Integer nseries, tdmrx, tdparx, tdxxy;
  Nag_ArimaOrder arimav;
  Nag_G13_Opt options;
  Nag_TransfOrder transfv;
  double *fsd = 0, *fva = 0, *para = 0, *parx = 0, *rmsxy = 0;
  double *xxy = 0;
  NagError fail;

  INIT_FAIL(fail);

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

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

  /*
   * Initialize the option-setting function.
   */
  /* nag_tsa_options_init (g13bxc).
   * Initialization function for option setting
   */
  nag_tsa_options_init(&options);
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "", &nev, &nfv, &nseries);
  if (nseries > 0 && nev > 0 && nfv > 0) {
    /*
     * Set option variable to the desired value.
     */
    options.cfixed = Nag_TRUE;
    /*
     * 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);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_tsa_transf_orders (g13byc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    /*
     * 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;
    ldparx = 8;
    if (npara >= 1) {
      if (!(fsd = NAG_ALLOC(nfv, double)) || !(fva = NAG_ALLOC(nfv, double)) ||
          !(para = NAG_ALLOC(npara, double)) ||
          !(parx = NAG_ALLOC(ldparx * (nseries - 1), double)) ||
          !(rmsxy = NAG_ALLOC(nseries, double)) ||
          !(xxy = NAG_ALLOC((nev + nfv) * (nseries), double)) ||
          !(mrx = NAG_ALLOC(7 * (nseries - 1), Integer))) {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
      tdmrx = nseries - 1;
      tdparx = nseries - 1;
      tdxxy = nseries;

      for (i = 0; i < npara; ++i)
        scanf("%lf", &para[i]);
      n = nev + nfv;
      for (i = 0; i < n; ++i)
        for (j = 0; j < nseries; ++j)
          scanf("%lf", &XXY(i, j));
      for (i = 0; i < nseries; ++i)
        scanf("%lf", &rmsxy[i]);
      for (i = 0; i < 7; ++i)
        for (j = 0; j < inser; ++j)
          scanf("%" NAG_IFMT "", &MRX(i, j));
      for (i = 0; i < 5; ++i)
        for (j = 0; j < inser; ++j)
          scanf("%lf", &PARX(i, j));

      /* nag_tsa_multi_inputmod_forecast (g13bjc), see above. */
      fflush(stdout);
      nag_tsa_multi_inputmod_forecast(
          &arimav, nseries, &transfv, para, npara, nev, nfv, xxy, tdxxy, rmsxy,
          mrx, tdmrx, parx, ldparx, tdparx, fva, fsd, &options, &fail);

      if (fail.code == NE_NOERROR || fail.code == NE_SOLUTION_FAIL_CONV ||
          fail.code == NE_MAT_NOT_POS_DEF) {
        printf("%1" NAG_IFMT " sets of observations were processed.\n", nev);
        printf("\nThe residual mean square for the output ");
        printf("series is %10.4f\n\n", rmsxy[nseries - 1]);
        printf("The forecast values and their standard errors are\n\n");
        printf("\n   i      fva      fsd\n\n");
        for (i = 0; i < nfv; ++i)
          printf("%4" NAG_IFMT "%10.3f%10.4f\n", i + 1, fva[i], fsd[i]);
        printf("\nThe values of z(t) and noise(t) are\n\n");
        printf("   i       z1         z2         z3         z4"
               "         z5       noise\n\n");
        for (i = 0; i < n; ++i) {
          printf("%4" NAG_IFMT "", i + 1);
          for (j = 0; j < nseries - 1; ++j)
            printf("%10.3f ", ZT(i, j));
          printf("%10.3f\n", options.noise[i]);
        }
      } else {
        printf("Error from nag_tsa_multi_inputmod_forecast (g13bjc)."
               "\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }
    } 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 more of nseries, nev and nfv are out of range:"
           " nseries = %-3" NAG_IFMT ",  nv = %-3" NAG_IFMT " while  "
           "nfv = %-3" NAG_IFMT "q\n",
           nseries, nev, nfv);
    exit_status = 1;
    goto END;
  }
  /* nag_tsa_free (g13xzc), see above. */
  nag_tsa_free(&options);
  /* nag_tsa_trans_free (g13bzc), see above. */
  nag_tsa_trans_free(&transfv);

END:
  NAG_FREE(fsd);
  NAG_FREE(fva);
  NAG_FREE(para);
  NAG_FREE(parx);
  NAG_FREE(rmsxy);
  NAG_FREE(xxy);
  NAG_FREE(mrx);

  return exit_status;
}