/**********************************************************************/
/*                                                                    */
/*           FILE: anovar.c                                           */
/*           WRITTEN BY: Jonathan G. Fiscus                           */
/*           DATE: April 14 1989                                      */
/*                  NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY    */
/*                  SPEECH RECOGNITION GROUP                          */
/*                                                                    */
/*           USAGE: This uses the rank structure to perform           */
/*                  a Friedman two-way analysis of variance           */
/*                  by ranks and generates a report.                  */
/*                                                                    */
/*           SOURCE:Applied Non Parametric Statistics by Daniel       */
/*                                                                    */
/**********************************************************************/

#include "sctk.h"

#define FRI_TITLE	"Friedman Two-way Analysis of Variance by Ranks"
#define RANK_TITLEA	"COMPARISON MATRIX: Comparing All Systems"
#define RANK_TITLEB	"COMPARISON MATRIX: Comparing All Speakers"
#define RANK_TITLE_1	"Using a Multiple Comparison Test"
#define ADJUST_THRESH	0.50000
#define MULTI_COMP_ALPHA	0.05


static int analyze_rank_sum(RANK *, int, int, double *, int, int, FILE *);
static void print_result_of_analyze_rank_sum(int, double, double, char *,FILE *);
static double dsum_sqr(double *, int);
static double calc_adjust_for_trt(RANK *);
static double calc_adjust_for_blks(RANK *);
static void do_multi_comp_for_anovar(int, int, char **, double *, double *, char *, int *, int **, int, FILE *);
static void calc_stat_ranks(int, int, int **, double*);
static double calc_comp_value(int, int);
static void print_ANOVAR_comp_matrix(int, int **, char **, char *, int *, FILE *);

/****************************************************************/
/*   main procedure to perform a two-way analysis of variance by*/
/*   ranks.  first the sum of the treatment and block ranks     */
/*   are summed then if there is a difference between at least  */
/*   one of the blocks or treatments, then a multiple comparison*/
/*   is performed                                               */
/****************************************************************/
void compute_anovar(RANK *rank, int verbose, int report, int ***out_sys_winner, char *outroot, int feedback, double ***out_conf)
{
    double *sum_trt_rank, *sum_blk_rank;
    char pad[FULL_SCREEN];
    int i,j, **trt_winner, **blk_winner;
    double **conf;
    FILE *fp = stdout;
    
    if (report || verbose){
	char *f = rsprintf("%s.anovar",outroot);
	if ((fp=(strcmp(outroot,"-") == 0) ? stdout : fopen(f,"w")) ==
	    (FILE *)0){
	    fprintf(stderr,"Warning: Open of %s for write failed.  "
		           "Using stdout instead.\n",f);
	    fp = stdout;
	} else
	    if (feedback >= 1) printf("        Output written to '%s'\n",f);
    }

    alloc_2dimZ(conf,rank->n_trt,rank->n_trt,double,0.0);
    alloc_2dimZ(trt_winner,rank->n_trt,rank->n_trt,int,NO_DIFF);
    alloc_2dimZ(blk_winner,rank->n_blk,rank->n_blk,int,NO_DIFF);
    *out_sys_winner = trt_winner;
    *out_conf = conf;

    alloc_singZ(sum_trt_rank,rnk_trt(rank),double,0.0);
    alloc_singZ(sum_blk_rank,rnk_blks(rank),double,0.0);

    /*  sum the ranks */
    for (i=0;i<rnk_blks(rank);i++)
        for (j=0;j<rnk_trt(rank);j++){
            sum_trt_rank[j] += Vrnk_t_rank(rank,i,j);
            sum_blk_rank[i] += Vrnk_b_rank(rank,i,j);
	}

    if (verbose){
        set_pad(pad,FRI_TITLE, FULL_SCREEN);
        fprintf(fp,"\n\n\n%s%s\n",pad,FRI_TITLE);
        fprintf(fp,"%s----------------------------------------------\n\n",pad);
        set_pad_cent_n(pad,SCREEN_WIDTH, FULL_SCREEN);
    fprintf(fp,"%s\tHo: Testing the hypothesis that all recognizers are the same",
                 pad);
    }
    /* if at least one difference, do multiple comparison */
    if (analyze_rank_sum(rank,rnk_trt(rank),rnk_blks(rank),
			 sum_trt_rank,
                         FOR_TREATMENTS,verbose,fp) == TEST_DIFF)
	clear_2dim(trt_winner,rank->n_trt,rank->n_trt,NO_DIFF);
        do_multi_comp_for_anovar(rnk_trt(rank),rnk_blks(rank),
                                 rnk_t_name(rank),sum_trt_rank,
                                 ovr_t_rank(rank),RANK_TITLEA,
                                 srt_t_rank(rank),trt_winner,verbose,fp);
    if (verbose) form_feed(fp);

    if (verbose){
        set_pad(pad,FRI_TITLE, FULL_SCREEN);
        fprintf(fp,"\n\n\n%s%s\n",pad,FRI_TITLE);
        fprintf(fp,"%s----------------------------------------------\n\n",pad);
        set_pad_cent_n(pad,SCREEN_WIDTH, FULL_SCREEN);
       fprintf(fp,"%s\tHo: Testing the hypothesis that all speakers are the same",
                 pad);
    }
    /* if at least one difference, do multiple comparison */
    if (analyze_rank_sum(rank,rnk_blks(rank),rnk_trt(rank),sum_blk_rank,
                         FOR_BLOCKS,verbose,fp) == TEST_DIFF)
	clear_2dim(blk_winner,rank->n_blk,rank->n_blk,NO_DIFF);
        do_multi_comp_for_anovar(rnk_blks(rank),rnk_trt(rank),
                                 rnk_b_name(rank),sum_blk_rank,
                                 ovr_b_rank(rank),RANK_TITLEB,
                                 srt_b_rank(rank),blk_winner, verbose,fp);
    if (verbose && (fp == stdout)) form_feed(fp);
    free_2dimarr(blk_winner,rank->n_blk,int); 
    free_singarr(sum_trt_rank,double);
    free_singarr(sum_blk_rank,double);
    if (fp != stdout) fclose(fp);
}

/*****************************************************************/
/*        X2_r formula:                                          */
/*                       k                                       */
/*             12        --.  2                                  */
/*          -------  *   \   R   - 3b(t+1)                       */
/*         bt(t+1)       /    j                                  */
/*                       --'                                     */
/*                      j = 1                                    */
/*                                                               */
/*   b = num_blk                                                 */
/*   t = num_trt                                                 */
/*   R = Sum of the ranks for that column                        */
/*                                                               */
/*  this is the initial test to see if the treatment or ranks    */
/*  are different                                                */
/*****************************************************************/
static int analyze_rank_sum(RANK *rank, int trt_num, int blk_num, double *sum_arr, int orient, int verbose,FILE *fp)
{
    int df;
    double X2_r, adjust;
    char *subject_blk = "speaker";
    char *subject_trt = "recognition system";
  
    /* calculate the test  statistic */
    df = trt_num - 2;     /* #trt is one based so subtract 2 */
    /* make sure df isn't negative */;
    if (df < 0) df = 0;
    /* make sure df isn't greater than MAX_DF */
    if (df > MAX_DF) df = MAX_DF;

    X2_r = (12.00000 /
                (double)(blk_num*trt_num*(trt_num+1))) *
           (dsum_sqr(sum_arr,trt_num)) -
           3.00000 * (double)blk_num * (double)(trt_num+1);
    if (orient == FOR_TREATMENTS)
        adjust = calc_adjust_for_trt(rank);
    else
        adjust = calc_adjust_for_blks(rank);

    /* if two identical systems are entered, adjust will be 0.0 and core*/
    /* dump.  this avoids this */
    if (adjust == 0.0)
       X2_r = 0.0;
    else
       X2_r /= adjust;

    if (verbose) 
        if (orient == FOR_TREATMENTS)
            print_result_of_analyze_rank_sum(df, X2_r, adjust, subject_trt,fp);        else
            print_result_of_analyze_rank_sum(df, X2_r, adjust, subject_blk,fp);

    /* return the result */
    if (X2_r > X2.df[df].level[GEN_X2_PER])
        return(TEST_DIFF);
    else
        return(NO_DIFF);
}

/*****************************************************************/
/*   print to stdout the results of the rank_sum test            */
/*****************************************************************/
static void print_result_of_analyze_rank_sum(int df, double X2_r, double adjust, char *subject, FILE *fp)
{
    char pad[FULL_SCREEN];
    int i;

    set_pad_cent_n(pad,SCREEN_WIDTH, FULL_SCREEN);
    fprintf(fp,"\n\n\n%s%35sReject if\n",pad,"");
    fprintf(fp,"%s%26sX2_r > X2 of %s %s (%2.3f)\n", pad, "", 
                          X2_pct_str(GEN_X2_PER),
                          X2_df_str(df),
                          X2_score(df,GEN_X2_PER));
    fprintf(fp,"\n"); 
    fprintf(fp,"%s%26sadjustment = %2.3f\n",pad,"",adjust);
    fprintf(fp,"%s%30s  X2_r = %2.3f\n",pad,"",X2_r);
    if (adjust < ADJUST_THRESH){
        fprintf(fp,"\n\n%s\t\t*** Warning:  ties adjustment may have severely\n",
                                                                   pad);
        fprintf(fp,"%s\t\t              exagerated the X2_r value\n\n",pad);
    }
    fprintf(fp,"%sANALYSIS:\n%s--------\n",pad,pad);
    if (X2_r > X2_score(df,GEN_X2_PER)){
     fprintf(fp,"%s\tThe test statistic X2_r shows, with %s confidence, that at\n"
                                            ,pad,X2.neg_per_str[GEN_X2_PER]),
     fprintf(fp,"%s\tleast one %s is significantly different.\n",pad,subject);
        fprintf(fp,"\n");
        fprintf(fp,"%s\tFurther, the probablity of there being a difference is\n",
               pad);
        for (i=GEN_X2_PER;i<MAX_X2_PER;i++)
            if (X2_r < X2_score(df,i+1))
                break;
        if (i==MAX_X2_PER)
            fprintf(fp,"%s\tgreater that %s.\n",pad,X2_neg_pct_str(i));
        else
            fprintf(fp,"%s\tbetween %s to %s.\n",pad,X2_neg_pct_str(i),
                                             X2_neg_pct_str(i+1));
        fprintf(fp,"\n\n");
    }
    else{
       fprintf(fp,"%s\tThe test statistic X2_r shows that at the %s confidence\n"
                          ,pad,X2_neg_pct_str(GEN_X2_PER)),
       fprintf(fp,"%s\tinterval, the %ss are not significantly\n",pad,subject);
        fprintf(fp,"%s\tdifferent.\n\n",pad);
        fprintf(fp,"%s\tFurther, the probablity of there being a difference is\n"
                          ,pad);
        for (i=GEN_X2_PER;i>MIN_X2_PER;i--)
            if (X2_r > X2_score(df,i-1))
                break;
        if (i==MIN_X2_PER)
            fprintf(fp,"%s\tless than %s.\n",pad,X2_neg_pct_str(i));
        else
            fprintf(fp,"%s\tbetween %s to %s.\n",pad,X2_neg_pct_str(i-1),
                                         X2_neg_pct_str(i));
        fprintf(fp,"\n\n");
    }
}

/********************************************************************/
/*  return the doubleing point value of the sum of the squares       */
/********************************************************************/
static double dsum_sqr(double *arr, int len)
{
    int i;
    double tot=0.000000;

    for (i=0;i<len;i++)
        tot += (arr[i]*arr[i]);
    return(tot);
}

/********************************************************************/
/*     Ties adjustments                                             */
/*              b                                                   */
/*             --.                                                  */
/*    1   -    \   T                                                */
/*             /    i                                               */
/*             --'                                                  */
/*          -----------                                             */
/*                2                                                 */
/*            bt(t -1)                                              */
/*                                                                  */
/*               3                                                  */
/*     T =  Sum t   - Sum t        where t = num ties in ith block  */
/*      i        i         i                                        */
/*                                                                  */
/*   to adjust for ties, the test statistic must be taken down      */
/*   by a certain factor, this function returns it                  */
/********************************************************************/
static double calc_adjust_for_trt(RANK *rank)
{
    double *sort_arr;
    int sort_num, *sort_count, i, j, sum_T=0;

    alloc_singarr(sort_arr,rnk_trt(rank),double);
    alloc_singarr(sort_count,rnk_trt(rank),int);

    /* find the ties ans sum Ti */
    for (i=0;i<rnk_blks(rank);i++){
        for (j=0;j<rnk_trt(rank); j++){
            sort_arr[j]=0.0;
            sort_count[j]=0;
	}
        for (j=0;j<rnk_trt(rank); j++){
            sort_num = 0;
            while ((sort_arr[sort_num] != 0.0) && 
                   (Vrnk_t_rank(rank,i,j) != sort_arr[sort_num]))
                sort_num++;
            sort_arr[sort_num] = Vrnk_t_rank(rank,i,j);
            sort_count[sort_num]++;
	}
        for (j=0;j<rnk_trt(rank); j++)
            if (sort_count[j] > 0)
                sum_T += (sort_count[j] * sort_count[j] * sort_count[j])
                           - sort_count[j];
    }
    free_singarr(sort_arr,double);
    free_singarr(sort_count,int);

    /* return the adjustment number */
    return( 1.000000 - (sum_T/
                        (double)(rnk_blks(rank)*
                                rnk_trt(rank)*
                                ((rnk_trt(rank) * rnk_trt(rank)) - 1))));
}

/********************************************************************/
/*  see calc_adjust_for_trt but substitute B for t                  */
/********************************************************************/
static double calc_adjust_for_blks(RANK *rank)
{
    double *sort_arr;
    int sort_num, *sort_count, i, j, sum_T=0;

    alloc_singarr(sort_arr,rnk_blks(rank),double);
    alloc_singarr(sort_count,rnk_blks(rank),int);

    /* find the ties ans sum Ti */
    for (i=0;i<rnk_trt(rank);i++){
        for (j=0;j<rnk_blks(rank); j++){
            sort_arr[j]=0.0;
            sort_count[j]=0;
	}
        for (j=0;j<rnk_blks(rank); j++){
            sort_num = 0;
            while ((sort_arr[sort_num] != 0.0) && 
                   (Vrnk_b_rank(rank,j,i) != sort_arr[sort_num]))
                sort_num++;
            sort_arr[sort_num] = Vrnk_b_rank(rank,j,i);
            sort_count[sort_num]++;
	}
        for (j=0;j<rnk_trt(rank); j++)
            if (sort_count[j] > 0)
                sum_T += (sort_count[j] * sort_count[j] * sort_count[j])
                           - sort_count[j];
    }
    free_singarr(sort_arr,double);
    free_singarr(sort_count,int);
    /* return the adjustment number */
    return( 1.000000 - (sum_T/
                        (double)(rnk_trt(rank)*
                                rnk_blks(rank)*
                                ((rnk_blks(rank) * rnk_blks(rank)) - 1))));
}

/********************************************************************/
/*  for the treatments, calculate and do the MULTIPLE comparison    */
/*  print everything out if requested                               */
/********************************************************************/
static void do_multi_comp_for_anovar(int trt_num, int blk_num, char **trt_names, double *sum_arr, double *ovr_rank_arr, char *title, int *srt_ptr, int **stat_sum, int verbose, FILE *fp)
{
    calc_stat_ranks(trt_num,blk_num,stat_sum,sum_arr);

    if (verbose)
        print_ANOVAR_comp_matrix(trt_num,stat_sum,trt_names,title,srt_ptr,fp);
}

/*********************************************************************/
/*      Multiple comparison test to rank the systems                 */
/*                             ________                              */
/*                            / bt(t+1)                              */
/*     | R  - R  |    >=   z / --------                              */
/*        j    i            V     6                                  */ 
/*                                                                   */
/*  R = the Sum of the nth  column ranks                             */
/*  z = Z score the alpha = 0.05                                     */
/*                          ----                                     */
/*                         t(t-1)                                    */
/*                                                                   */
/*  go through all possible comparisons and rank the treatments for  */
/*  according to the number of treatments beaten or beaten by        */
/*                                                                   */
/*********************************************************************/
static void calc_stat_ranks(int trt_num, int blk_num, int **stat_sum, double *sum_arr)
{
    int comp1, comp2, result;
    double Zcomp_value, compare_value;
    double sqr_root;

    /* calculate the comparison value for the confidence interval */
    Zcomp_value = calc_comp_value(trt_num,blk_num);
    sqr_root = sqrt(blk_num * trt_num * (trt_num - 1) / 6.000000);
    compare_value = (double)Zcomp_value * (double)sqr_root;

    /* go through the comparison and build an array to rank */
    for (comp1=0; comp1 <(trt_num-1); comp1++)
        for (comp2=comp1+1; comp2<(trt_num); comp2++){
            if (fabs(sum_arr[comp1] - sum_arr[comp2]) >= compare_value)
                if (sum_arr[comp1] > sum_arr[comp2])
                    result = 1;
                else
                    result = -1;
            else
                result = 0;
            stat_sum[comp1][comp2] += result;
            stat_sum[comp2][comp1] += result * (-1);
	}
}

/*********************************************************************/
/*   Given the alpha, find the closest approximate Z index to that   */
/*   percentage                                                      */
/*********************************************************************/
static double calc_comp_value(int trt_num, int blk_num)
{
    double Z_area;
    int i;

    Z_area = MULTI_COMP_ALPHA / (double)(trt_num * (trt_num-1));

    if (Z_area < Z_exter(Z1tail,MAX_Z_PER))
        return(Z_score(Z1tail,MAX_Z_PER));

    for (i=MAX_Z_PER+1; i< MIN_Z_PER ; i++)
        if (Z_area < Z_exter(Z1tail,i))
            return(Z_score(Z1tail,i));

    return(Z_score(Z1tail,MIN_Z_PER));
}


static void print_ANOVAR_comp_matrix(int trt_num, int **stat_sum, char **trt_names, char *title, int *srt_ptr, FILE *fp)
{
    char pad[FULL_SCREEN],name_format[50],*hyphens="--------------------";
    char *spaces="                    ";
    int t, t2;
    int max_trt_name_len=4, hy_l=20;

    set_pad(pad,title, FULL_SCREEN);
    fprintf(fp,"\n\n%s%s\n",pad,title);
    set_pad(pad,RANK_TITLE_1, FULL_SCREEN);
    fprintf(fp,"%s%s\n\n\n\n",pad,RANK_TITLE_1);

    for (t=0; t<trt_num; t++)
        if (max_trt_name_len < strlen(trt_names[t]))
            max_trt_name_len = strlen(trt_names[t]);
    set_pad_cent_n(pad,(trt_num+1) * (max_trt_name_len+2+1), FULL_SCREEN);
    sprintf(name_format,"| %%-%ds ",max_trt_name_len);

    /* first line */
    fprintf(fp,"%s|%s",pad,(hyphens+(hy_l-(max_trt_name_len+2))));
    for (t=0; t<trt_num; t++)
       fprintf(fp,"%s",(hyphens+(hy_l-(max_trt_name_len+3))));
    fprintf(fp,"|\n");

    /* systems on top */
    fprintf(fp,"%s|%s",pad,(spaces+(hy_l-(max_trt_name_len+2))));
    for (t=0; t<trt_num; t++)
       fprintf(fp,name_format,trt_names[srt_ptr[t]]);
    fprintf(fp,"|\n"); 

    /* separation line */
    fprintf(fp,"%s|%s+",pad,(hyphens+(hy_l-(max_trt_name_len+2))));
    for (t=0; t<trt_num-1; t++)
       fprintf(fp,"%s",(hyphens+(hy_l-(max_trt_name_len+3))));
    fprintf(fp,"%s|\n",(hyphens+(hy_l-(max_trt_name_len+2))));

    /* systems */
    for (t=0; t<trt_num; t++){
       fprintf(fp,"%s",pad);
       fprintf(fp,name_format,trt_names[srt_ptr[t]]);
       for (t2=0; t2<trt_num; t2++){
           if (t < t2)
               switch (stat_sum[srt_ptr[t]][srt_ptr[t2]]) {
                   case -1: fprintf(fp,name_format,trt_names[srt_ptr[t]]);
                            break;
                   case 1:  fprintf(fp,name_format,trt_names[srt_ptr[t2]]);
                            break;
                   default: fprintf(fp,name_format,"same");
                            break;
	       }
           else
               fprintf(fp,name_format,"");
       }
       fprintf(fp,"|\n");
    }

    /* last line */
    fprintf(fp,"%s|%s",pad,(hyphens+(hy_l-(max_trt_name_len+2))));
    for (t=0; t<trt_num; t++)
       fprintf(fp,"%s",(hyphens+(hy_l-(max_trt_name_len+3))));
    fprintf(fp,"|\n");
}