/* 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);
}