/* file: ibs.c              Albert Yang, MD                25 August 2004

-------------------------------------------------------------------------------
ibs: calculates information-based similarity index (IBS) between two
data sets

This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.  See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 59 Temple
Place - Suite 330, Boston, MA 02111-1307, USA.

You may contact the author by e-mail (ccyang@physionet.org/).  For updates to
this software, please visit PhysioNet (http://www.physionet.org/).
-----------------------------------------------------------------------------

Compile this program by

    gcc -o ibs -O ibs.c -lm

Input usage

ibs m file1 file2

m: m-tuple binary sequence
file1 and file2: input time series

The program performs two main calculations:
1. Determine the word rank order frequency (WROF) list of a given time
series.  
2. Calculate the information-based similarity between two WROF lists.

Output file:
the information-based similarity index using m-tuple binary sequences.

*/


#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAXDAT 150000

int m, nw;

static int rr1[MAXDAT];
static int rr2[MAXDAT];


void swap (int array[2][nw], int i, int j) {
	int k, temp;
	for(k=0 ; k<2 ; k++) { 
           temp = array[k][i];
	   array[k][i] = array[k][j];
	   array[k][j] = temp;
	}   
}


int partition(int array[2][nw], int left, int right, int ic) {
	int pivot = array[ic][left];
	int i = left-1; 
	int j = right+1;

	while (1) {
		do {j--;
		} while (array[ic][j] > pivot);
		
		do {i++;
		} while (array[ic][i] < pivot);
		
		if (i < j) 
			swap(array, i, j); // swap ith and jth elements of array
		else
			return j;
	}
}

void quicksort(int array[2][nw], int left, int right, int ic) {
	int pivot_index;
	if (left < right) {
		pivot_index = partition(array, left, right, ic);
		quicksort(array, left, pivot_index, ic);
		quicksort(array, pivot_index+1, right, ic);
	}
}

float distance(int rr1[], int rr2[], int nlin1, int nlin2)
{
  float d;
  int i, j, mtuple, ns1 = nlin1 - m + 1, ns2 = nlin2 - m + 1;
  float f1, f2, dr, p, pi, pj;
  
  int w1[2][nw];
  int w2[2][nw];
  int temp[nw];
  
  for (i=0 ; i<nw; i++){
     w1[0][i]=i;
     w1[1][i]=0;
     w2[0][i]=i;
     w2[1][i]=0;
  }   
  
  for (i=0 ; i<ns1; i++){
     mtuple = 0;
     for (j=1 ; j<=m ; j++){mtuple = mtuple + (1 << (m-j))*rr1[i+j-1];}
     w1[1][mtuple] = w1[1][mtuple] + 1;
  }

  for (i=0 ; i<ns2; i++){
     mtuple = 0;
     for (j=1 ; j<=m ; j++){mtuple = mtuple + (1 << (m-j))*rr2[i+j-1];}
     w2[1][mtuple] = w2[1][mtuple] + 1;
  }
  
  for (i=0 ; i<nw; i++){ temp[i] = w1[1][i]; }
  quicksort( w1, 0, nw - 1, 1 );
  for (i=0 ; i<nw; i++){ w1[1][i] = nw - i - 1; }
  quicksort( w1, 0, nw - 1, 0 );
  for (i=0 ; i<nw; i++){ w1[0][i] = temp[i]; }
  
  
  for (i=0 ; i<nw; i++){ temp[i] = w2[1][i]; }
  quicksort( w2, 0, nw - 1, 1 );
  for (i=0 ; i<nw; i++){ w2[1][i] = nw - i - 1; }
  quicksort( w2, 0, nw - 1, 0 );
  for (i=0 ; i<nw; i++){ w2[0][i] = temp[i]; }

 
  /*Calculate Total Occurrence*/
  f1=0;
  f2=0;
  for (i=0 ; i<nw ; i++)
  {
    f1 = f1 + w1[0][i];
    f2 = f2 + w2[0][i];
  }
  
  /*Calculate Distance(W1,W2)*/
  dr = 0;
  p = 0;
  for (i=0 ; i<nw ; i++)
  {
    pi = w1[0][i]/f1;
    pj = w2[0][i]/f2;
    if (pi > 0) pi = -pi * log( pi ); else pi = 0;
    if (pj > 0) pj = -pj * log( pj ); else pj = 0;
    dr = dr + abs( w1[1][i] - w2[1][i] ) * (pi + pj);
    p = p + ( pi + pj );

  }

  if (p > 0) d = dr / p / (nw-1); else d = -1;  
   
  return (d);
}


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

     int i, nlin1, nlin2;
     float rrn,rrn1;
     
     m=atol(argv[1]);
     nw=1 << m;

     /*
     for (i = 1; i < argc; ) {
          printf("%s", argv[i]);
          if (++i < argc)
               putchar(' ');
     putchar('\n');
     }
     */
      
     FILE *RF1;   /* File variable for RR interval time series 1 */
     FILE *RF2;   /* File variable for RR interval time series 2 */
     
     RF1 = fopen(argv[2], "r");
     RF2 = fopen(argv[3], "r");

     /* Read RR interval time series 1 */
     i = 0;
     fscanf(RF1, "%f", &rrn);
     while(!feof(RF1)) 
     { 
       nlin1 = i;  
       fscanf(RF1, "%f", &rrn1);
       if (rrn1 > rrn) rr1[i] = 1; else rr1[i] = 0; 
       rrn = rrn1; 
       i++;
     }

     /* Read RR interval time series 2*/ 
     i = 0;      
     fscanf(RF2, "%f", &rrn);
     while(!feof(RF2)) 
     { 
       nlin2 = i;  
       fscanf(RF2, "%f", &rrn1);
       if (rrn1 > rrn) rr2[i] = 1;  else rr2[i] = 0; 
       rrn = rrn1; 
       i++;
     }
     
     
    int fclose( FILE *stream );
    
    /*
    for(i=0 ; i<nlin2 ; i++) { 
      printf("%d\n", rr2[i]);
    } 
    */
     
    
    printf("%f\n",distance(rr1, rr2, nlin1, nlin2));  

    exit(0);

}