/*
 *                    Copyright (c) 1995 Giovanni Flammia 
 *
 *                    Massachusetts Institute of Technology 
 *                    Laboratory for Computer Science
 *                    Spoken Language Systems Group
 *                    http://www.sls.lcs.mit.edu 
 *  
 *      	      Contact: gflammia AT alum.mit.edu
 *
 *                    File: kappa.c
 *
 *                    Compute generalized kappa coefficient from
 *                    data classified by multiple coders into mutually
 *                    exclusive and exhaustive set of categories.
 *  
 *                    For more information about the formulae used,
 *                    read the appendix of Giovanni Flammia PhD Thesis
 *                    "Discourse Segmentation of Spoken Dialogue: An Empirical Approach" 
 *                    PDF document located at: http://www.theredesign.com/Technology/Dialogue
 *
 *                    Directly implements the computations described in the paper:
 *                    Uebersax, J.S.: "A Generalized Kappa Coefficient"
 *                    Educational and Psychological Measurement
 *                    1982, Vol. 42, pp. 181--183.
 *
 * By using this software, you ("Licensee") agree to the following license:
 *
 * 1. Redistribution and use in source and binary forms, with or without
 *    modification, are permitted provided that the following conditions
 *    are met:
 *
 * 2. You are granted a non-exclusive, non-transferable limited license 
 *    to use this Software free of charge for 
 *    non-commercial purposes only.  
 *
 * 3. Redistributions of source or binary form of this code must 
 *    retain the above copyright notice, this list of conditions 
 *    and the following disclaimer.
 *
 * 4. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 
 *    OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 
 *    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN 
 *    NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, 
 *    DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR 
 *    OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE 
 *    USE OR OTHER DEALINGS IN THE SOFTWARE.  
 *                    
 */


  /* To compile: just use cc or gcc

     Usage: kappa <r_matrix_text_file>

     Required input file format for the r_matrix_text_file:

     -------------------------------------------------------------------------
     input_directory = /where-all-the-data-comes-from/
     num_categories = 3
     categories = [ Initiative Response Other ] (optional)
     annotated_files = [ (optional)
     aa_30 40 
     aa_5 58 
     aa_55 38 
     st1t1_51 43 
     st1t4_18 66 
     st2t1_32 119 
     ]
     num_coders     = 5
     coders         = [ (optional)
     christine
     jwc
     laura
     raylau
     rob
     ]
     num_points     = 364 (i.e. number of lines for the R matrix )
     R_matrix = [ 
     0 1 0 1 0 1 0 1 0 1 
     1 0 1 0 0 1 1 0 0 1 
     ]

     Format of the R matrix: 
     The R matrix tells how each line in the data has been coded by
     each one of the coders. It is a 3-dimensional matrix 
     i.e. num_points * num_coders * num_categories 

     Each line is divided into num_coders * num_categories columns:

     For example, if coder 1 has assigned cat 2 to line m and coder 2 
     has assigned cat 1 to line m, then the corresponding line in R_matrix
     should look like this:

              coder[1]                               coder[2]          
              cat[1] cat[2] ... cat[num_categ]       cat[1] cat[2] ... cat[num_categ] 
     line[m]  0      1      ... 0                    1      0          0              

     R[c][k][m] = 1 if coder c assigned class k to data point m, 
     else R[c][k][m] = 0 (also if coder c has not assigned any category
     to line m). 

      
     */




#include <stdio.h>
#include <time.h>
#include <sys/param.h>
#include <sys/types.h>
#include <string.h>
#include <ctype.h>
#include <math.h>

int     num_categories,num_coders,num_points;


double compute_term1(R,i,j,k)
  int ***R;
  int i;
  int j;
  int k;
{
  
  int ii,jj;         /* range over 0 .. num_coders-1 */
  int kk;            /* range over 0 .. num_categories-1 */
  int l,ll;          /* range over 0 .. num_points-1 */
  int denom_i=0;
  int denom_j=0;
  int num_i=0;
  int num_j=0;
  
  for(kk=0;kk<num_categories;kk++) {
    for(ll=0;ll<num_points;ll++) {
      denom_i += R[i][kk][ll];
      denom_j += R[j][kk][ll];
    }
  }
   
  for(l=0;l<num_points;l++) {
    num_i += R[i][k][l];
    num_j += R[j][k][l];
  }

  return(((double)num_i/(double)denom_i)*((double)num_j/(double)denom_j));
  
}



double compute_term2(R,i,j,l)
  int ***R;
  int i;
  int j;
  int l;
{
  int ii,jj;         /* range over 0 .. num_coders-1 */
  int k,kk;          /* range over 0 .. num_categories-1 */
  int ll;            /* range over 0 .. num_points-1 */
  int sum_i=0;
  int sum_j=0;

  for(k=0;k<num_categories;k++) {
    sum_i += R[i][k][l];
    sum_j += R[j][k][l];
  }

  return( (double)sum_i * (double)sum_j );
  
}



double compute_denom(R)
  int ***R;
{
 
  int     i;              /* range over 0 .. num_coders-1 */
  int     k,kk;           /* range over 0 .. num_categories-1 */
  int     l,ll;           /* range over 0 .. num_points-1 */
  int     *sum;
  double  result=0;

  sum = (int *)calloc(num_points,sizeof(int));
  for (l=0;l<num_points;l++) {
    sum[l]=0;
    for(i=0;i<num_coders;i++) {
      for(k=0;k<num_categories;k++) {
	sum[l] += R[i][k][l];
      }
    }
  }
  for (l=0;l<num_points;l++) {
    result += (double)sum[l] * (double)( sum[l] - 1);
  }

  return( result );

  
}



double compute_partial_pchance(R,i,j,denom)
  int     ***R;
  int     i;
  int     j;
  double  denom;
{
  int ii,jj;          /* range over 0 .. num_coders-1 */
  int k,kk;           /* range over 0 .. num_categories-1 */
  int l,ll;           /* range over 0 .. num_points-1 */
  double term1=0;
  double term2=0;

  for(k=0;k<num_categories;k++) {
    term1 += compute_term1(R,i,j,k);
  }

  for(l=0;l<num_points;l++) {
    term2 += compute_term2(R,i,j,l);
  }

  return( term1 * ( term2 / denom ) );

}
  


double compute_partial_pobs(R,k,l)
  int ***R;
  int k;
  int l;
{
  int i; /* range over 0 .. num_coders-1 */
  int sum=0;

  for (i=0;i<num_coders;i++) {
    sum += R[i][k][l];
  }

  return( (double)sum * (double)(sum - 1) );
}



void main(argc, argv)
	int	argc; 
	char	**argv;
{

  FILE    *file;
  int     nd,x,frequency;
  double  pobs,pchance,kappa,num,denom,partial_pchance,a_priori;
  char    info[4000],line[2048],*line_ptr;
  char    word1[2048],word2[2048];
  int     is_data;
  float   *average_frequency;

  /* The R matrix defined as in the Uebersax paper: */
  int     i,j,ii,jj;       /* range over 0 .. num_coders-1 */
  int     k,kk;            /* range over 0 .. num_categories-1 */
  int     l,ll,n,m,nn,mm;  /* range over 0 .. num_points-1 */
  int     ***R;            /* num_coders x num_categories x num_points */


  /* Check command line argument */
  if (argc<2) {
    fprintf(stderr,"Usage: %s <r_matrix_text_file>\n",argv[0]);
    exit();
  }

  if ((file=fopen(argv[1],"r"))==NULL) {
    fprintf(stderr,"%s: File %s not found.\n",argv[0],argv[1]);
    exit();
  }

  /* Read R matrix data and put it into info and R data structures 
     Input format:
     
     input_directory = /usr/users/flammia/eurospeech_95/data/by_coders/group_1
     num_categories = 2
     categories = [ Open_Segment Other ]
     annotated_files = [ 
     aa_30 40 
     aa_5 58 
     aa_55 38 
     st1t1_51 43 
     st1t4_18 66 
     st2t1_32 119 
     ]
     num_coders     = 5
     coders         = [
     christine
     jwc
     laura
     raylau
     rob
     ]
     num_points     = 364
     R_matrix = [ 
     0 1 0 1 0 1 0 1 0 1 
     1 0 1 0 0 1 1 0 0 1 
  
     */

  is_data=0;
  while(fgets(line, 2048, file) != NULL && strcmp(line,"\n")!=0) {
    x=0;
    if ((nd = sscanf(line,"%s %s %d",word1,word2,&x))==3 &&
	strcmp(word2,"=")==0 &&
	x>0) {
      /* read num_categories: */
      if (strcmp(word1,"num_categories")==0) {
	num_categories = x;
      } else if (strcmp(word1,"num_coders")==0) {
	/* read num_coders: */
	num_coders = x;
      } else if (strcmp(word1,"num_points")==0) {
	/* read num_points */
	num_points = x;
      } 
    }
    else if (nd>=2 && strcmp(word1,"R_matrix")==0) {
      is_data=1;
      R = (int ***)calloc(num_coders,sizeof(int **));
      for(i=0;i<num_coders;i++) {
	R[i] = (int **)calloc(num_categories,sizeof(int *));
	for (k=0;k<num_categories;k++) {
	  R[i][k] = (int *)calloc(num_points,sizeof(int));
	  for (m=0;m<num_points;m++) {
	    R[i][k][m]=0;
	  }
	}
      }
      m = -1;
    }
    if (is_data==0 && strlen(line)>=1 && strcmp(line,"\n")!=0) {
      strcat(info,"\n");
      strcat(info,line);
    } else if (is_data==1) {
      /* Read data: for each line:
	 coder[1]                      coder[2]                    ...
 	 cat[1]...cat[num_categ]       cat[1]...cat[num_categ]     ...
         0 0  ... 1 .....      0       0 1 ....           0  0     ...
	 R[c][k][m] = 1 if coder c assigned class k to data point m,
      */ 
      line_ptr=line;
      i=0;
      k=0;
      if (line_ptr && (nd=sscanf(line_ptr,"%d",&x))==1) {
	m++;
	/* fprintf(stderr,"R[][][%d] %s",m,line_ptr); */
      }
      while (m>=0 && m<num_points && line_ptr && (nd=sscanf(line_ptr,"%d",&x))==1) {
	R[i][k][m] = x;
	/* fprintf(stderr,"R[%d][%d][%d]=%d ",i,k,m,R[i][k][m]); */
	k++;
	if (k==num_categories) {
	  i++;
	  k=0;
	}
	nd=sscanf(line_ptr,"%s",word1);
	line_ptr=line_ptr+strlen(word1)+1;
      }
      /* fprintf(stderr,"\n"); */
    }
  }
  fprintf(stderr,"num_categories = %d num_coders = %d num_points = %d\n",
	  num_categories,num_coders,num_points);
  
  /* 1. Do kappa analysis for all coders: */

  /* 1.a. compute chance probability of agreement pchance for all coders */
  pchance = 0;
  denom = compute_denom(R);
  average_frequency = (float *)calloc(num_categories,sizeof(float));
  for (k=0;k<num_categories;k++) {
    average_frequency[k] = 0.0;
  }
  for (i=0;i<num_coders;i++) {
    for(j=0;j<num_coders;j++) {
      if (i!=j) {
	partial_pchance = compute_partial_pchance(R,i,j,denom);
	/* fprintf(stderr,"PChance[%d][%d] = %f\n",i,j,partial_pchance); */
	pchance += partial_pchance;
      }
    }
    for (k=0;k<num_categories;k++) {
      frequency=0;
      for(l=0;l<num_points;l++) {
	frequency += R[i][k][l];

      }
      a_priori = (double)frequency/(double)num_points;
      fprintf(stderr,"Frequency[%d][%d] = %f\n",i,k,a_priori);
    }
  }
  for (k=0;k<num_categories;k++) {
    for (l=0;l<num_points;l++) {
      for(i=0;i<num_coders;i++) {
	average_frequency[k] += (float)R[i][k][l];
      }
    }
  }
  for(k=0;k<num_categories;k++) {
    average_frequency[k] = average_frequency[k]/((float)num_points*(float)num_coders);
    fprintf(stderr,"Average_frequency[%d] = %f\n",k,average_frequency[k]);
  }
  fprintf(stderr,"PChance = %f\n",pchance);

  /* 1.b. compute observed probability of agreement pobs for all coders */
  num = 0;
  for (k=0;k<num_categories;k++) {
    for(l=0;l<num_points;l++) {
      num += compute_partial_pobs(R,k,l);
    }
  }
  pobs = num / denom;
  fprintf(stderr,"PObs = %f\n",pobs);

  kappa = ( pobs - pchance ) / ( 1.0 - pchance );
  fprintf(stderr,"Kappa = %f\n",kappa);

  /* Future extensions... */

  /* 1.c do statistical test for assessing whether pobs is different from */
  /* pchance: if there are only two coders compute the Z-test and if there
  /* are more than two coders do a Monte Carlo randomized experiment */
 
  /* 1.d from the R_matrix isolate the areas of agreement from the areas of
     disagreement in terms of rows of the R_matrix */

  /* 2. If there are more than two coders: Run the analysis 1.a to 1.d again
     once for each pair  of coders */

}