/* nag_wfilt (c09aac) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */
/* Pre-processor includes */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagc09.h>

int main(void)
{
  /* Constants */
  Integer licomm = 100;
  /*Integer scalar and array declarations */
  Integer exit_status = 0;
  Integer i, n, nf, nnz, nwc, nwlmax, ny;
  Integer *dwtlev = 0, *icomm = 0;
  NagError fail;
  Nag_Wavelet wavnamenum;
  Nag_WaveletMode modenum;
  /*Double scalar and array declarations */
  double *c = 0, *x = 0, *y = 0;
  /*Character scalar and array declarations */
  char mode[24], wavnam[20];

  INIT_FAIL(fail);

  printf("nag_wfilt (c09aac) Example Program Results\n\n");
  fflush(stdout);

  /*     Skip heading in data file */
  scanf("%*[^\n] ");
  /*     Read n - length of input data sequence */
  scanf("%" NAG_IFMT "%*[^\n] ", &n);
  if (!(x = NAG_ALLOC(n, double)) ||
      !(y = NAG_ALLOC(n, double)) || !(icomm = NAG_ALLOC(licomm, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  /*     Read Wavelet name (wavnam) and end mode (mode) */
  scanf("%19s%23s%*[^\n] ", wavnam, mode);
  /*
   * nag_enum_name_to_value (x04nac).
   * Converts NAG enum member name to value
   */
  wavnamenum = (Nag_Wavelet) nag_enum_name_to_value(wavnam);
  modenum = (Nag_WaveletMode) nag_enum_name_to_value(mode);
  if (n >= 2) {
    printf("   Parameters read from file :: \n");
    printf("     Wavelet  :%15s\n", wavnam);
    printf("     End mode :%15s\n", mode);
    printf("     N        :%15" NAG_IFMT "\n\n", n);
    /*        Read data array and write it out */
    printf("%s\n", "   Input Data           X :");
    for (i = 0; i < n; i++) {
      scanf("%lf ", &x[i]);
      printf("%8.3f%s", x[i], (i + 1) % 8 ? " " : "\n");
    }
    scanf("%*[^\n] ");
    printf("\n");
    /*
     * nag_wfilt (c09aac)
     * Wavelet filter query
     */
    nag_wfilt(wavnamenum, Nag_MultiLevel, modenum, n, &nwlmax, &nf, &nwc,
              icomm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_wfilt (c09aac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    if (!(c = NAG_ALLOC(nwc, double)) ||
        !(dwtlev = NAG_ALLOC(nwlmax + 1, Integer)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    /*        Perform Discrete Wavelet transform */
    /*
     * nag_mldwt (c09ccc)
     * one-dimensional multi-level discrete wavelet transform (mldwt)
     */
    nag_mldwt(n, x, nwc, c, nwlmax, dwtlev, icomm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_mldwt (c09ccc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    printf("   Length of wavelet filter :            %10" NAG_IFMT "\n", nf);
    printf("   Number of Levels :                    %10" NAG_IFMT "\n",
           nwlmax);
    printf("   Number of coefficients in each level : \n");
    for (i = 0; i < nwlmax + 1; i++)
      printf("    %8" NAG_IFMT "%s", dwtlev[i], (i + 1) % 8 ? " " : "\n");
    printf("\n");
    printf("   Total number of wavelet coefficients :%10" NAG_IFMT "\n", nwc);
    nnz = 0;
    for (i = 0; i < nwlmax + 1; i++)
      nnz = nnz + dwtlev[i];
    printf("\n");
    printf("   Wavelet coefficients C : \n");
    for (i = 0; i < nnz; i++)
      printf("%8.3f%s", c[i], (i + 1) % 8 ? " " : "\n");
    printf("\n");
    /*           Reconstruct original data */
    ny = n;
    /*
     * nag_imldwt (c09cdc)
     * one-dimensional inverse multi-level discrete wavelet transform
     * (imldwt)
     */
    nag_imldwt(nwlmax, nwc, c, n, y, icomm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_imldwt (c09cdc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    printf("\n");
    printf("   Reconstruction       Y : \n");
    for (i = 0; i < ny; i++)
      printf("%8.3f%s", y[i], (i + 1) % 8 ? " " : "\n");
  }

END:
  NAG_FREE(c);
  NAG_FREE(dwtlev);
  NAG_FREE(x);
  NAG_FREE(y);
  NAG_FREE(icomm);

  return exit_status;
}