/* nag_anova_random}(g04bbc) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 *
 * Mark 6 revised, 2000.
 * Mark 8 revised, 2004.
 */

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

#define C(I, J)     c[(I) *tdc + J]
#define TABLE(I, J) table[(I) *tdtable + J]
int main(int argc, char *argv[])
{
  FILE       *fpin, *fpout;

  Integer    exit_status = 0, i, irdf, j, n, nblock, nt, tdc, tdtable;
  Integer    *irep = 0, *it = 0;
  NagError   fail;
  Nag_Blocks blocks;
  double     gmean, tol;
  double     *bmean = 0, *c = 0, *ef = 0, *r = 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_random (g04bbc) Example Program Results\n");

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

  if ((nblock >= 2 && !(n%nblock)) || (nblock == 0))
    {
      if (nblock != 0)
        {
          if (!(bmean = NAG_ALLOC(nblock, double)))
            {
              fprintf(fpout, "Allocation failure\n");
              exit_status = -1;
              goto END;
            }
          blocks = Nag_SerialBlocks;
        }
      else
        {
          blocks = Nag_NoBlocks;
        }

      if (!(c = NAG_ALLOC((nt)*(nt), double)) ||
          !(ef = NAG_ALLOC(nt, double)) ||
          !(r = NAG_ALLOC(n, double)) ||
          !(table = NAG_ALLOC((4)*(5), double)) ||
          !(tmean = NAG_ALLOC(nt, double)) ||
          !(y = NAG_ALLOC(n, double)) ||
          !(irep = NAG_ALLOC(nt, Integer)) ||
          !(it = NAG_ALLOC(n, Integer)))
        {
          fprintf(fpout, "Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tdc = nt;
      tdtable = 5;
    }
  else
    {
      fprintf(fpout, "Invalid nblock or n.\n");
      exit_status = 1;
      return exit_status;
    }
  for (i = 0; i < n; ++i)
    fscanf(fpin, "%lf", &y[i]);
  fscanf(fpin, "%*[^\n]");
  for (i = 0; i < n; ++i)
    fscanf(fpin, "%ld", &it[i]);
  fscanf(fpin, "%*[^\n]");

  tol = 5e-6;
  irdf = 0;

  /* nag_anova_random (g04bbc).
   * General block design or completely randomized design
   */
  nag_anova_random(n, y, blocks, nblock, nt, it, &gmean, bmean, tmean,
                   table, c, tdc, irep, r, ef, tol, irdf,
                   &fail);
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_anova_random (g04bbc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  fprintf(fpout, "\nANOVA table\n\n");
  fprintf(fpout, "   Source        df         SS          MS          F"
          "          Prob\n\n");
  fprintf(fpout, " Blocks      ");
  for (j = 0; j < 5; ++j)
    fprintf(fpout, "%10.4f  ", TABLE(0, j));
  fprintf(fpout, "\n");

  fprintf(fpout, " Treatments  ");
  for (j = 0; j < 5; ++j)
    fprintf(fpout, "%10.4f  ", TABLE(1, j));
  fprintf(fpout, "\n");

  fprintf(fpout, " Residual    ");
  for (j = 0; j < 3; ++j)
    fprintf(fpout, "%10.4f  ", TABLE(2, j));
  fprintf(fpout, "\n");

  fprintf(fpout, " Total       ");
  for (j = 1; j <= 2; ++j)
    fprintf(fpout, "%10.4f  ", TABLE(3, j-1));
  fprintf(fpout, "\n");

  fprintf(fpout, "\nEfficiency Factors\n\n");
  for (i = 0; i < nt; ++i)
    fprintf(fpout, "%10.5f", ef[i]);
  fprintf(fpout, "\n");

  fprintf(fpout, "\n%s%10.5f\n", "  Grand Mean", gmean);

  fprintf(fpout, "\nTreatment Means\n\n");
  for (i = 1; i <= nt; ++i)
    fprintf(fpout, "%10.5f", tmean[i-1]);
  fprintf(fpout, "\n");

  fprintf(fpout, "\nStandard errors of differences between means\n\n");
  for (i = 1; i < nt; ++i)
    {
      for (j = 0; j < i; ++j)
        fprintf(fpout, "%10.5f", C(i, j));
      fprintf(fpout, "\n");
    }
 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  if (bmean) NAG_FREE(bmean);
  if (c) NAG_FREE(c);
  if (ef) NAG_FREE(ef);
  if (r) NAG_FREE(r);
  if (table) NAG_FREE(table);
  if (tmean) NAG_FREE(tmean);
  if (y) NAG_FREE(y);
  if (irep) NAG_FREE(irep);
  if (it) NAG_FREE(it);
  return exit_status;
}