/* nag_anova_factorial (g04cac) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.1, 2023.
 *
 */

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

#define MTERM 6

#define TABLE(I, J) table[((I)-1) * (5) + ((J)-1)]
int main(void) {

  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);

  printf("nag_anova_factorial (g04cac) Example Program Results\n\n");

  /* Skip heading in data file */
  scanf("%*[^\n]");
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\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))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  } else {
    printf("Invalid n or nfac or nblock.\n");
    exit_status = 1;
    return exit_status;
  }
  for (j = 0; j < nfac; ++j)
    scanf("%" NAG_IFMT "", &lfac[j]);
  scanf("%*[^\n]");

  for (i = 0; i < n; ++i)
    scanf("%lf", &y[i]);
  scanf("%*[^\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) {
    printf("Error from nag_anova_factorial (g04cac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

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

  printf("%s   ", " Total   ");
  printf("%4" NAG_IFMT "  ", (Integer)TABLE(mterm, 1));
  printf("%12.2f\n\n", TABLE(mterm, 2));

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

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

    printf("\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:
  NAG_FREE(r);
  NAG_FREE(y);
  NAG_FREE(lfac);
  return exit_status;
}