/***********************************************************************
 * Program to compute three-factor analysis of variance
 * with equal replication.
 *
 * Author: M.V.Wickerhauser
 * Date: 28 March 2005
 *
 */
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>


/* ********** STATIC (LOCAL) FUNCTIONS ********** */

/* ************************************************************
 * Make and return a 1-dimensional array X[a]
 */
static float *
array1d(int a) {
  float *x;

  x = (float *)malloc(a*sizeof(float)); assert(x);
  return x;
}

/* ************************************************************
 * Make and return a 2-dimensional array X[a][b]
 */
static float **
array2d(int a, int b) {
  float **x;
  int i;

  x = (float **)malloc(a*sizeof(float *)); assert(x);
  for(i=0; i<a; i++)  x[i] = array1d(b);
  return x;
}

/* ************************************************************
 * Make and return a 3-dimensional array X[a][b][c]
 */
static float ***
array3d(int a, int b, int c) {
  float ***x;
  int i;

  x = (float ***)malloc(a*sizeof(float **)); assert(x);
  for(i=0; i<a; i++) x[i] = array2d(b,c);
  return x;
}

/* ************************************************************
 * Make and return a 4-dimensional array X[a][b][c][n]
 */
static float ****
array4d(int a, int b, int c, int n) {
  float ****x;
  int i;

  x = (float ****)malloc(a*sizeof(float ***)); assert(x);
  for(i=0; i<a; i++)  x[i] = array3d(b,c,n);
  return x;
}

/* ************************************************************
 * Sum a 1-dimensional array X[n]
 */
static float
sum1d(float *x, int n) {
  float sum;
  int l;
  sum=0;
  for(l=0; l<n; l++) sum += x[l];
  return sum;
}

/* ************************************************************
 * Sum a 4-dimensional array X[a][b][c][n]
 */
static float
sum4d(float ****x, int a, int b, int c, int n) {
  float sum;
  int i, j, k;
  sum=0;
  for(i=0; i<a; i++) 
    for(j=0; j<b; j++) 
      for(k=0; k<c; k++) sum += sum1d(x[i][j][k], n);
  return sum;
}

/* ************************************************************
 * Return the square of the input
 */
static float
square(float x) {
  return x*x;
}


static void 
usage(char *name) 
{
  fprintf(stderr, "Usage: %s a b c n < datafile\n", name);
  fprintf(stderr, " where    a  = number of Factor A levels (index i)\n");
  fprintf(stderr, "          b  = number of Factor B levels (index j)\n");
  fprintf(stderr, "          c  = number of Factor C levels (index k)\n");
  fprintf(stderr, "          n  = number of replications per cell (index l)\n");
  fprintf(stderr, "    datafile = text file of observations indexed ijkl in the order below:\n");
  fprintf(stderr, "              A1                    A2          ...           Aa         \n");
  fprintf(stderr, "     -------------------   -------------------       ------------------- \n");
  fprintf(stderr, "      B1    B2  ...  Bb     B1    B2  ...  Bb   ...   B1    B2  ...  Bb  \n");
  fprintf(stderr, "     ----  ----     ----   ----  ----     ----       ----  ----     ---- \n");
  fprintf(stderr, " C1| 1111  1211 ... 1b11   2111  2211 ... 2b11       a111  a211 ... ab11 \n");
  fprintf(stderr, "   | 1112  1212 ... 1b12   2112  2212 ... 2b12       a112  a212 ... ab12 \n");
  fprintf(stderr, "   | ...                   ...                       ...                 \n");
  fprintf(stderr, "   | 111n  121n ... 1b1n   211n  221n ... 2b1n       a11n  a21n ... ab1n \n");
  fprintf(stderr, " _______________________________________________________________________ \n");
  fprintf(stderr, " C2| 1121  1221 ... 1b21   2121  2221 ... 2b21       a121  a221 ... ab21 \n");
  fprintf(stderr, "   | 1122  1222 ... 1b22   2122  2222 ... 2b22       a122  a222 ... ab22 \n");
  fprintf(stderr, "   | ...                   ...                       ...                 \n");
  fprintf(stderr, "   | 112n  122n ... 1b2n   212n  222n ... 2b2n       a12n  a22n ... ab2n \n");
  fprintf(stderr, " ...                                            ...                      \n");
  fprintf(stderr, " _______________________________________________________________________ \n");
  fprintf(stderr, " Cc| 11c1  12c1 ... 1bc1   21c1  22c1 ... 2bc1       a1c1  a2c1 ... abc1 \n");
  fprintf(stderr, "   | 11c2  12c2 ... 1bc2   21c2  22c2 ... 2bc2       a1c2  a2c2 ... abc2 \n");
  fprintf(stderr, "   | ...                   ...                       ...                 \n");
  fprintf(stderr, "   | 11cn  12cn ... 1bcn   21cn  22cn ... 2bcn       a1cn  a2cn ... abcn \n");
  exit(1);
}





/* ******************** MAIN PROGRAM ******************** */
int
main(int argc, char **argv) {
  float ****x;			/* Observations x[i][j][k][l] go here */
  float mX, *mA, *mB, *mC;	/* Total, and factor A, B, C means */
  float **mAB, **mBC, **mAC;	/* AxB, BxC, AxC interaction means */
  float ***mABC;		/* Cell mean == AxBxC interaction mean */
  float ssT, ssE, ssc;		/* Sums of squares */
  float ssA, ssB, ssC;		/* Sums of squares */
  float ssAB, ssBC, ssAC, ssABC; /* Sums of squares */
  float msT, msE, msc;		/* Mean squares */
  float msA, msB, msC;		/* Mean squares */
  float msAB, msBC, msAC, msABC; /* Mean squares */
  int dfT, dfE, dfc;		/* Degrees of freedom */
  int dfA, dfB, dfC;		/* Degrees of freedom */
  int dfAB, dfBC, dfAC, dfABC;	/* Degrees of freedom */
  int a, b, c, n, i, j, k, l;

  /* Check command line arguments: */
  if(argc!=5) usage(argv[0]);
  a=atoi(argv[1]);  b=atoi(argv[2]); c=atoi(argv[3]);  n=atoi(argv[4]);

  /* Make space for the data */
  x = array4d(a,b,c,n);

  /* Read the data as if it were tabulated as follows:
   *
   *              A1                    A2          ...           Aa        
   *     -------------------   -------------------       -------------------
   *      B1    B2  ...  Bb     B1    B2  ...  Bb   ...   B1    B2  ...  Bb 
   *     ----  ----     ----   ----  ----     ----       ----  ----     ----
   * C1| 1111  1211 ... 1b11   2111  2211 ... 2b11       a111  a211 ... ab11
   *   | 1112  1212 ... 1b12   2112  2212 ... 2b12       a112  a212 ... ab12
   *   | ...                   ...                       ...                
   *   | 111n  121n ... 1b1n   211n  221n ... 2b1n       a11n  a21n ... ab1n
   * _______________________________________________________________________      
   * C2| 1121  1221 ... 1b21   2121  2221 ... 2b21       a121  a221 ... ab21
   *   | 1122  1222 ... 1b22   2122  2222 ... 2b22       a122  a222 ... ab22
   *   | ...                   ...                       ...                
   *   | 112n  122n ... 1b2n   212n  222n ... 2b2n       a12n  a22n ... ab2n
   * ...                                            ...                     
   * _______________________________________________________________________
   * Cc| 11c1  12c1 ... 1bc1   21c1  22c1 ... 2bc1       a1c1  a2c1 ... abc1
   *   | 11c2  12c2 ... 1bc2   21c2  22c2 ... 2bc2       a1c2  a2c2 ... abc2
   *   | ...                   ...                       ...                
   *   | 11cn  12cn ... 1bcn   21cn  22cn ... 2bcn       a1cn  a2cn ... abcn
   *
   */
  for(k=0; k<c; k++)		/* slowest-changing index: blocks of rows */
    for(l=0; l<n; l++)		/* next slowest change: rows */
      for(i=0; i<a; i++)	/* next fastest change: blocks of columns */
	for(j=0; j<b; j++)	/* fastest-changing: columns == row elements */
	  scanf("%f", &x[i][j][k][l]);

  /* ******************** MEANS ******************** */
  /* TOTAL MEAN mX */
  mX = sum4d(x, a,b,c,n) / (a*b*c*n);

  printf("Total mean: %10.3f\n", mX);


  /* CELL MEANS == AxBxC INTERACTION MEANS */
  mABC = array3d(a,b,c);	/* Make space for the AxBxC means */
  for(i=0; i<a; i++) 
    for(j=0; j<b; j++)
      for(k=0; k<c; k++) mABC[i][j][k] = sum1d(x[i][j][k], n) / n;

  /* FACTOR A MEANS */
  mA = array1d(a);		/* Make space for the factor A means */
  for(i=0; i<a; i++) {
    mA[i] = 0;
    for(j=0; j<b; j++)
      for(k=0; k<c; k++) mA[i] += mABC[i][j][k];
    mA[i] /= b*c;
  }

  /* FACTOR B MEANS */
  mB = array1d(b);		/* Make space for the factor B means */
  for(j=0; j<b; j++) {
    mB[j] = 0;
    for(i=0; i<a; i++)
      for(k=0; k<c; k++) mB[j] += mABC[i][j][k];
    mB[j] /= a*c;
  }

  /* FACTOR C MEANS */
  mC = array1d(c);		/* Make space for the factor C means */
  for(k=0; k<c; k++) {
    mC[k] = 0;
    for(i=0; i<a; i++)
      for(j=0; j<b; j++) mC[k] += mABC[i][j][k];
    mC[k] /= a*b;
  }


  /* AxB INTERACTION MEANS */
  mAB = array2d(a,b);		/* Make space for the AxB means */
  for(i=0; i<a; i++) 
    for(j=0; j<b; j++) {
      mAB[i][j] = 0;
      for(k=0; k<c; k++) mAB[i][j] += mABC[i][j][k];
      mAB[i][j] /= c;
  }
  /* AxC INTERACTION MEANS */
  mAC = array2d(a,c);		/* Make space for the AxC means */
  for(i=0; i<a; i++) 
    for(k=0; k<c; k++) {
      mAC[i][k] = 0;
      for(j=0; j<b; j++) mAC[i][k] += mABC[i][j][k];
      mAC[i][k] /= b;
  }
  /* BxC INTERACTION MEANS */
  mBC = array2d(b,c);		/* Make space for the BxC means */
  for(j=0; j<b; j++)
    for(k=0; k<c; k++) {
      mBC[j][k] = 0;
      for(i=0; i<a; i++) mBC[j][k] += mABC[i][j][k];
      mBC[j][k] /= a;
  }



  /* ******************** SUMS OF SQUARES ******************** */
  /* TOTAL sum of squares */
  ssT = 0;
  for(i=0; i<a; i++)
    for(j=0; j<b; j++)
      for(k=0; k<c; k++)
	for(l=0; l<n; l++)
	  ssT += square(x[i][j][k][l] - mX);

  /* CELLS sum of squares */
  ssc = 0;
  for(i=0; i<a; i++)
    for(j=0; j<b; j++)
      for(k=0; k<c; k++)
	ssc += square(mABC[i][j][k] - mX);
  ssc *= n;

  /* ERROR sum of squares */
  ssE = 0;
  for(i=0; i<a; i++)
    for(j=0; j<b; j++)
      for(k=0; k<c; k++)
	for(l=0; l<n; l++)
	  ssE += square(x[i][j][k][l] - mABC[i][j][k]);


  printf("Total SS - (Cell SS + Error SS) = %f-(%f+%f)=%f\n",
	 ssT, ssc, ssE, ssT-(ssc+ssE));


  /* FACTOR A sum of squares */
  ssA = 0;
  for(i=0; i<a; i++) ssA += square(mA[i] - mX);
  ssA *= b*c*n;

  /* FACTOR B sum of squares */
  ssB = 0;
  for(j=0; j<b; j++) ssB += square(mB[j] - mX);
  ssB *= a*c*n;

  /* FACTOR C sum of squares */
  ssC = 0;
  for(k=0; k<c; k++) ssC += square(mC[k] - mX);
  ssC *= a*b*n;


  /* AxB INTERACTION sum of squares */
  ssAB = 0;
  for(i=0; i<a; i++)
    for(j=0; j<b; j++)
      ssAB += square(mAB[i][j] - mA[i] - mB[j] + mX);
  ssAB *= c*n;

  /* AxC INTERACTION sum of squares */
  ssAC = 0;
  for(i=0; i<a; i++)
    for(k=0; k<c; k++)
      ssAC += square(mAC[i][k] - mA[i] - mC[k] + mX);
  ssAC *= b*n;

  /* BxC INTERACTION sum of squares */
  ssBC = 0;
  for(j=0; j<b; j++)
    for(k=0; k<c; k++)
      ssBC += square(mBC[j][k] - mB[j] - mC[k] + mX);
  ssBC *= a*n;

  /* AxBxC INTERACTION sum of squares */
  ssABC = 0;
  for(i=0; i<a; i++)
    for(j=0; j<b; j++)
      for(k=0; k<c; k++)
	ssABC += square(mABC[i][j][k]
			- mAB[i][j] - mAC[i][k] - mBC[j][k]
			+ mA[i] + mB[j] + mC[k]
			- mX);
  ssABC *= n;


  /* ******************** DEGREES OF FREEDOM ******************** */

  dfT = a*b*c*n-1;
  dfc = a*b*c-1;
  dfE = a*b*c*(n-1);
  dfA = a-1;
  dfB = b-1;
  dfC = c-1;
  dfAB = (a-1)*(b-1);
  dfAC = (a-1)*(c-1);
  dfBC = (b-1)*(c-1);
  dfABC = (a-1)*(b-1)*(c-1);

  /* ******************** MEAN SQUARE DEVIATION ******************** */

  msT = ssT/dfT;
  msc = ssc/dfc;
  msA = ssA/dfA;     msB = ssB/dfB;     msC = ssC/dfC;
  msAB = ssAB/dfAB;  msBC = ssBC/dfBC;  msAC = ssAC/dfAC;
  msABC = ssABC/dfABC;
  if(dfE>0)  msE = ssE/dfE;


  /* ******************** OUTPUT ******************** */

  printf("-------------------------------------------------------------\n");
  printf("Source of variation    Sum of squares     DF      Mean Square\n");
  printf("-------------------------------------------------------------\n");

  printf("Total                   %12.3f     %3d     %12.3f\n", ssT, dfT, msT );
  printf("Cells                   %12.3f     %3d     %12.3f\n", ssc, dfc, msc );
  printf("Factor A                %12.3f     %3d     %12.3f\n", ssA, dfA, msA );
  printf("Factor B                %12.3f     %3d     %12.3f\n", ssB, dfB, msB );
  printf("Factor C                %12.3f     %3d     %12.3f\n", ssC, dfC, msC );
  printf("A x B Interaction       %12.3f     %3d     %12.3f\n", ssAB, dfAB, msAB );
  printf("A x C Interaction       %12.3f     %3d     %12.3f\n", ssAC, dfAC, msAC );
  printf("B x C Interaction       %12.3f     %3d     %12.3f\n", ssBC, dfBC, msBC );
  printf("A x B x C Interaction   %12.3f     %3d     %12.3f\n", ssABC, dfABC, msABC );
  if(dfE>0) {
    printf("Within-cell errors      %12.3f     %3d     %12.3f\n", ssE, dfE, msE );
    printf("-------------------------------------------------------------\n");
    printf("Source of F value            num DF    den DF            F   \n");
    printf("-------------------------------------------------------------\n");
    printf("(Factor A MS) / (Error MS)     %3d       %3d     %12.3f\n", dfA, dfE, msA/msE );
    printf("(Factor B MS) / (Error MS)     %3d       %3d     %12.3f\n", dfB, dfE, msB/msE );
    printf("(Factor C MS) / (Error MS)     %3d       %3d     %12.3f\n", dfC, dfE, msC/msE );
    printf(" (A x B  MS) / (Error MS)      %3d       %3d     %12.3f\n", dfAB, dfE, msAB/msE );
    printf(" (A x C  MS) / (Error MS)      %3d       %3d     %12.3f\n", dfAC, dfE, msAC/msE );
    printf(" (B x C  MS) / (Error MS)      %3d       %3d     %12.3f\n", dfBC, dfE, msBC/msE );
    printf(" (A x B x C  MS) / (Error MS)  %3d       %3d     %12.3f\n", dfABC, dfE, msABC/msE );
  }
  printf("-------------------------------------------------------------\n");

  return 0;
}