/* nag_tsa_multi_inp_update (g13bgc) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 8, 2004.
 */

#include <stdio.h>
#include <string.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagg13.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 CM(I, J)   options.cm[(I - 1) * options.tdcm + (J - 1)]
#define XXY(I, J)  xxy[(I - 1) * tdxxy + J - 1]
#define XXYN(I, J) xxyn[(I - 1) * tdxxyn + J - 1]
#define ZT(I, J)   options.zt[(I - 1) * options.tdzt + (J - 1)]

  INIT_FAIL(fail);

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

  /* Skip heading in data file */
  scanf("%*[^\n]");
  /* Initialise the option structure */
  /* nag_tsa_options_init (g13bxc).
   * Initialization function for option setting
   */
  nag_tsa_options_init(&options);
  scanf("%ld%ld%ld%ld", &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("%ld%ld%ld%ld%ld"
             "%ld%ld", &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("%ld", &transfv.b[j-1]);
        }
      for (j = 1; j <= inser; ++j)
        {
          scanf("%ld", &transfv.q[j-1]);
        }
      for (j = 1; j <= inser; ++j)
        {
          scanf("%ld", &transfv.p[j-1]);
        }
      for (j = 1; j <= inser; ++j)
        {
          scanf("%ld", &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_inp_model_estim (g13bec).
       * Estimation for time series models
       */
      fflush(stdout);
      nag_tsa_multi_inp_model_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_inp_model_estim "
                  "(g13bec).\n%s\n.", fail.message);
          exit_status = 1;
          goto END;
        }

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

      if (fail.code != NE_NOERROR)
        {
          printf("\nError from "
                  "nag_tsa_multi_inp_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("%2ld%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("%2ld", 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 = %-3ld while nser = %-3ld\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;
}