/* nag_anova_factorial (g04cac) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 *
 * Mark 6 revised, 2000.
 * Mark 7b revised, 2004.
 * Mark 8 revised, 2004.
 */

#include <nag.h>
#include <nagx04.h>
#include <nag_stdlib.h>
#include <stdio.h>
#include <nagg04.h>

#define MTERM 6

#define LFAC(I)     lfac[(I) -1]
#define IWK(I)      iwk[(I) -1]
#define IMEAN(I)    imean[(I) -1]
#define Y(I)        y[(I) -1]
#define TMEAN(I)    tmean[(I) -1]
#define SEMEAN(I)   semean[(I) -1]
#define R(I)        r[(I) -1]
#define E(I)        e[(I) -1]
#define BMEAN(I)    bmean[(I) -1]
#define TABLE(I, J) table[((I) -1) * (5) + ((J) -1)]
int main(int argc, char *argv[])
{
  FILE     *fpin, *fpout;

  Integer  c__27 = 27, exit_status = 0, i, *imean = 0, inter, irdf, j, k, l;
  Integer  *lfac = 0;
  Integer  mterm = MTERM, n, nblock, nfac, ntreat, num;
  NagError fail;
  double   *bmean = 0, *e = 0, *r = 0, *semean = 0, *table = 0, *tmean = 0;
  double   *y = 0;
  INIT_FAIL(fail);

  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);

  fprintf(fpout, "nag_anova_factorial (g04cac) Example Program Results\n\n");

  /* Skip heading in data file */
  fscanf(fpin, "%*[^\n]");

  fscanf(fpin, "%ld%ld%ld%ld%*[^\n]", &n, &nblock,
         &nfac, &inter);
  if (n >= 4 && nfac >= 1 && nblock >= 1 && !(n%nblock))
    {
      if (!(r = NAG_ALLOC(n, double)) ||
          !(y = NAG_ALLOC(n, double)) ||
          !(lfac = NAG_ALLOC(nfac, Integer)))
        {
          fprintf(fpout, "Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  else
    {
      fprintf(fpout, "Invalid n or nfac or nblock.\n");
      exit_status = 1;
      return exit_status;
    }
  for (j = 0; j < nfac; ++j)
    fscanf(fpin, "%ld", &lfac[j]);
  fscanf(fpin, "%*[^\n]");

  for (i = 0; i < n; ++i)
    fscanf(fpin, "%lf", &y[i]);
  fscanf(fpin, "%*[^\n]");

  irdf = 0;
  /* nag_anova_factorial (g04cac).
   * Complete factorial design
   */
  nag_anova_factorial(n, y, nfac, lfac, nblock, inter, irdf, &mterm, &table,
                      &tmean, &c__27, &e, &imean, &semean, &bmean, r,
                      &fail);
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_anova_factorial (g04cac).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  fprintf(fpout, "\n  ANOVA table\n\n");
  fprintf(fpout, "   Source      df         SS          MS           F"
          "          Prob\n\n");
  k = 0;
  if (nblock > 1)
    {
      ++k;
      fprintf(fpout, "%s   ", " Blocks  ");
      fprintf(fpout, "%4ld  ", (Integer) TABLE(1, 1));
      for (j = 2; j <= 5; ++j)
        fprintf(fpout, "%12.2f", TABLE(1, j));
      fprintf(fpout, "\n");
    }
  ntreat = mterm - 2 - k;
  for (i = 1; i <= ntreat; ++i)
    {
      fprintf(fpout, "%s%2ld ", " Effect  ", i);
      fprintf(fpout, "%4ld  ", (Integer) TABLE(k+i, 1));
      for (j = 2; j <= 5; ++j)
        fprintf(fpout, "%12.2f", TABLE(k+i, j));
      fprintf(fpout, "\n");
    }
  fprintf(fpout, "%s   ", " Residual");
  fprintf(fpout, "%4ld  ", (Integer) TABLE(mterm-1, 1));
  for (j = 2; j <= 3; ++j)
    fprintf(fpout, "%12.2f", TABLE(mterm-1, j));
  fprintf(fpout, "\n");

  fprintf(fpout, "%s   ", " Total   ");
  fprintf(fpout, "%4ld  ", (Integer) TABLE(mterm, 1));
  fprintf(fpout, "%12.2f\n\n", TABLE(mterm, 2));

  fprintf(fpout, "  Treatment Means and Standard Errors \n\n");
  k = 0;
  for (i = 0; i < ntreat; ++i)
    {
      l = imean[i];
      fprintf(fpout, "%s%2ld\n\n", " Effect ", i+1);

      num = 1;
      for (j = k; j < l; ++j)
        {
          fprintf(fpout, "%10.2f%s", tmean[j], num%8?"":"\n");
          num++;
        }

      fprintf(fpout, "\n\n%s%10.2f\n\n", " SE of difference in means  = ",
              semean[i]);
      k = l;
    }
  /* nag_anova_factorial_free (g04czc).
   * Memory freeing function for nag_anova_factorial (g04cac)
   */
  nag_anova_factorial_free(&table, &tmean, &e, &imean, &semean, &bmean);
 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  if (r) NAG_FREE(r);
  if (y) NAG_FREE(y);
  if (lfac) NAG_FREE(lfac);
  return exit_status;
}