/* nag_tsa_multi_inputmod_update (g13bgc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.3, 2023.
 */

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

int main(void) {
  double df, objf, rss;
  Integer exit_status = 0, i, inser, j, kzef, nnv, npara, nser, nxxy, tdxxy,
          tdxxyn;
  double *para = 0, *sd = 0, *xxy = 0, *xxyn = 0;
  /* Nag types */
  Nag_ArimaOrder arimav;
  Nag_TransfOrder transfv;
  Nag_G13_Opt options;
  NagError fail;

#define XXY(I, J) xxy[(I - 1) * tdxxy + J - 1]
#define XXYN(I, J) xxyn[(I - 1) * tdxxyn + J - 1]

  INIT_FAIL(fail);

  printf("nag_tsa_multi_inputmod_update (g13bgc) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n]");
  /* 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 "%" NAG_IFMT "", &nxxy, &nser,
        &options.max_iter, &nnv);

  if (nxxy > 0 && nser > 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).
     * Allocates memory to transfer function model orders
     */
    nag_tsa_transf_orders(nser, &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 = nser - 1;
    for (j = 1; j <= inser; ++j) {
      scanf("%" NAG_IFMT "", &transfv.b[j - 1]);
    }
    for (j = 1; j <= inser; ++j) {
      scanf("%" NAG_IFMT "", &transfv.q[j - 1]);
    }
    for (j = 1; j <= inser; ++j) {
      scanf("%" NAG_IFMT "", &transfv.p[j - 1]);
    }
    for (j = 1; j <= inser; ++j) {
      scanf("%" NAG_IFMT "", &transfv.r[j - 1]);
    }
    scanf("%*[^\n]");

    npara = 0;
    for (i = 1; i <= inser; ++i) {
      npara += transfv.q[i - 1] + transfv.p[i - 1];
    }
    npara += arimav.p + arimav.q + arimav.bigp + arimav.bigq + nser;

    tdxxy = nser;
    tdxxyn = nser;
    /* Memory allocation */
    if (!(para = NAG_ALLOC(npara, double)) ||
        !(sd = NAG_ALLOC(npara, double)) ||
        !(xxy = NAG_ALLOC(nxxy * tdxxy, double)) ||
        !(xxyn = NAG_ALLOC(nnv * tdxxyn, double))) {
      printf("Memory allocation failure\n");
      exit_status = -1;
      goto END;
    }

    for (i = 1; i <= npara; ++i) {
      scanf("%lf", &para[i - 1]);
    }
    scanf("%*[^\n]");

    for (i = 1; i <= nxxy; ++i) {
      for (j = 1; j <= nser; ++j) {
        scanf("%lf", &XXY(i, j));
      }
    }
    scanf("%*[^\n]");

    for (i = 1; i <= nnv; ++i) {
      for (j = 1; j <= nser; ++j) {
        scanf("%lf", &XXYN(i, j));
      }
    }
    scanf("%*[^\n]");

    options.print_level = Nag_NoPrint;
    /* nag_tsa_multi_inputmod_estim (g13bec).
     * Estimation for time series models
     */
    fflush(stdout);
    nag_tsa_multi_inputmod_estim(&arimav, nser, &transfv, para, npara, nxxy,
                                 xxy, tdxxy, sd, &rss, &objf, &df, &options,
                                 &fail);

    if (fail.code != NE_NOERROR) {
      printf("\nError from nag_tsa_multi_inputmod_estim "
             "(g13bec).\n%s\n.",
             fail.message);
      exit_status = 1;
      goto END;
    }

    /* Calculate update */
    kzef = 1;
    /* nag_tsa_multi_inputmod_update (g13bgc).
     * Multivariate time series, update state set for
     * forecasting from multi-input model
     */
    nag_tsa_multi_inputmod_update(&arimav, nser, &transfv, para, npara, nnv,
                                  xxyn, tdxxyn, kzef, &options, &fail);

    if (fail.code != NE_NOERROR) {
      printf("\nError from "
             "nag_tsa_multi_inputmod_update (g13bgc).\n%s\n.",
             fail.message);
      exit_status = 1;
      goto END;
    }

    printf("\n\nThe residuals (after differencing)\n");
    for (i = 1; i <= nnv; i++) {
      printf("%2" NAG_IFMT "%12.4f\n", i, options.res[i - 1]);
    }
    printf("\n\n");

    printf("\nThe values of z(t) and n(t)\n");
    for (i = 1; i <= nnv; i++) {
      printf("%2" NAG_IFMT "", i);
      for (j = 1; j <= nser; j++) {
        printf(" %9.4f", XXYN(i, j));
      }
      printf("\n");
    }
    printf("\n");
  } else {
    if (nxxy <= 0 || nser <= 0) {
      printf("One or both of nxxy and nser are out of range: "
             "nxxy = %-3" NAG_IFMT " while nser = %-3" NAG_IFMT "\n",
             nxxy, nser);
      exit_status = -1;
      goto END;
    }
  }

END:
  /* nag_tsa_trans_free (g13bzc).
   * Freeing function for the structure holding the transfer
   * function model orders
   */
  nag_tsa_trans_free(&transfv);
  /* nag_tsa_free (g13xzc).
   * Freeing function for use with g13 option setting
   */
  nag_tsa_free(&options);
  NAG_FREE(para);
  NAG_FREE(sd);
  NAG_FREE(xxy);
  NAG_FREE(xxyn);

  return exit_status;
}