/* file: pnnlist.c	Joe E. Mietus	19 November 2002
			Last revised:   23 February 2003
-------------------------------------------------------------------------------
pnnlist: derive pNNx statistics from an annotation interval list
Copyright (C) 2003 Joe E. Mietus

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 (joe at physionet dot org) or postal mail
(BIDMC Room KB-28, 300 Brookline Ave., Boston, MA 02215 USA).  For updates to
this software, please visit PhysioNet (http://www.physionet.org/).
_______________________________________________________________________________

A detailed description of the application of this algorithm can be
found in: JE Mietus, C-K Peng, I Henry, R L Goldsmith and AL Goldberger.
"The pNNx files: re-examining a widely used heart rate variability measure."
Heart 2002; vol. 88, pp. 378-380. 
*/

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

#define ABS(A) ((A) < 0 ? -(A) : (A))
#define ROFFERR 1e-9

char *pname, *prog_name();

main(argc, argv)
int argc;
char *argv[];
{
    char ann, lastann[2];
    int pflag, sflag, dblcmp();
    double rr, lastrr, inc, nninc, mult, *pnn = NULL;
    long i, maxdat, sum, tot, totn, totp;
    void bsd_qsort(), help();

    pname = prog_name(argv[0]);

    inc = pflag = sflag = 0;
    mult = 1000;

    for (i = 0; ++i < argc && *argv[i] == '-'; ) {
        switch(argv[i][1]) {
	    case 'i':
		if (i < argc-1)
		    inc = atof(argv[++i]);
		else {
		    fprintf(stderr, "%s: increment must follow -i\n", pname);
		    exit(1);
		}
		break;
	    case 'p':
		pflag = 1;
		mult = 100;
		break;
	    case 's':
		sflag = 1;
		break;
            default:
		help();
		exit(1);
	}
    }
    if (inc < 0) {
        fprintf(stderr, "%s -i inc : increment must be greater than 0\n",
		pname);
	exit(1);
    }
    else if (inc > 0) {
      if (pflag)
	  inc /= 100;
      else
	  inc /= 1000;
    }

    if (scanf("%lf %c", &rr, &ann) != 2) {
        fprintf(stderr, "%s : improperly formatted data\n\n", pname);
	help();
        exit(2);
    }

    i = maxdat = 0;
    lastann[0] = 'X';
    lastann[1] = ann;

    while (scanf("%lf %c", &rr, &ann) == 2) {

        if (lastann[0] == 'N' && lastann[1] == 'N' && ann == 'N') {
	    if (i >= maxdat) {
		double *s;

		maxdat += 50000;	/* allow the pnn array to grow (the
					   increment is arbitrary) */
		if ((s = realloc(pnn, maxdat * sizeof(double))) == NULL) {
		    fprintf(stderr,
			    "%s: truncating input after %ld increments\n",
			    pname, i);
		    break;
		}
		pnn = s;
	    }
	    
	    if (pflag)
		pnn[i] = (rr-lastrr)/lastrr;
	    else
		pnn[i] = rr-lastrr;

	    if (!sflag)
		pnn[i] = ABS(pnn[i]);

	    i++;
	}

	lastrr = rr;
	lastann[0] = lastann[1];
	lastann[1] = ann;
    }

    tot = i;
    if (tot <= 0) {
	fprintf(stderr, "%s: no consecutive NN intervals in input\n", pname);
	exit(3);
    }

    bsd_qsort(pnn, tot, sizeof(double), dblcmp);

    if (sflag) { /* signed distributions */
        for (totn = totp = i = 0; i < tot; i++) {
	    if (pnn[i] <= 0)
	        totn++;
	    if (pnn[i] >= 0)
	        totp++;
	}
	if (inc == 0) {
            printf("%g\t0\n", mult*pnn[0]);
            for (sum = i = 1; i < tot; i++) {
	        if (pnn[i] < 0) {
		    if (pnn[i] - pnn[i-1] > ROFFERR)
	                printf("%g\t%g\n", mult*pnn[i], 100*(double)sum/totn);
		    sum++;
	        }
	        else if (pnn[i] == 0 && pnn[i-1] != 0) {
		    printf("%g\t%g\n", mult*pnn[i], 100*(double)sum/totn);
		    sum = 1;
	        }
	        else if (pnn[i] == 0 && pnn[i-1] == 0)
		    sum++;
	        else {
		    if (pnn[i] - pnn[i-1] > ROFFERR)
		        printf("%g\t%g\n", mult*pnn[i-1],
			       100*(1-(double)sum/totp));
		    sum++;
	        }
	    }
	    printf("%g\t%g\n", mult*pnn[i-1], 100*(1-(double)sum/totp));
        }
	else {  /* inc > 0 */
	    for (nninc = 0; nninc > pnn[0]; nninc -= inc)
	        ;
            printf("%g\t0\n", mult*nninc);
	    nninc += inc;
            for (sum = i = 0; i < tot; i++) {
	        if (pnn[i] < 0) {
		    if (pnn[i] - nninc > ROFFERR) {
	                printf("%g\t%g\n", mult*nninc, 100*(double)sum/totn);
			nninc += inc;
			while (pnn[i] - nninc > ROFFERR) {
	                    printf("%g\t%g\n", mult*nninc,
				   100*(double)sum/totn);
			    nninc += inc;
			}
		    }
		    sum++;
	        }
	        else if (pnn[i] == 0 && pnn[i-1] != 0) {
		    printf("%g\t%g\n", mult*nninc, 100*(double)sum/totn);
		    sum = 1;
	        }
	        else if (pnn[i] == 0 && pnn[i-1] == 0)
		    sum++;
	        else {
		    if (pnn[i] - nninc > ROFFERR) {
		        printf("%g\t%g\n", mult*nninc,
			       100*(1-(double)sum/totp));
		        nninc += inc;
		        while (pnn[i] - nninc > ROFFERR) {
		            printf("%g\t%g\n", mult*nninc,
				   100*(1-(double)sum/totp));
		            nninc += inc;
			}
		    }
		    sum++;
	        }
	    }
	    printf("%g\t%g\n", mult*nninc, 100*(1-(double)sum/totp));
	}
    }

    else {  /* unsigned distributions */
        if (inc == 0) {
            for (sum = i = 1; i < tot; i++) {
	        if (pnn[i] - pnn[i-1] > ROFFERR)
	            printf("%g\t%g\n", mult*pnn[i-1], 100*(1-(double)sum/tot));
	        sum++;
            }
            printf("%g\t%g\n", mult*pnn[i-1], 100*(1-(double)sum/tot));
        }
	else {  /* inc > 0 */
            for (nninc = 0, sum = i = 0; i < tot; i++) {
	        if (pnn[i] - nninc > ROFFERR) {
	            printf("%g\t%g\n", mult*nninc, 100*(1-(double)sum/tot));
		    nninc += inc;
		    while (pnn[i] - nninc > ROFFERR) {
	                printf("%g\t%g\n", mult*nninc,
			       100*(1-(double)sum/tot));
		        nninc += inc;
		    }
		}
	        sum++;
            }
            printf("%g\t%g\n", mult*nninc, 100*(1-(double)sum/tot));
	}
    }
    if (pnn) free(pnn);
    exit(0);
}


int dblcmp(y1, y2)
double *y1, *y2;
{

    if (*y1 < *y2)
        return(-1);
    else if (*y1 > *y2)
        return(1);
    else
        return(0);

}

char *prog_name(s)
char *s;
{
    char *p = s + strlen(s);

#ifdef MSDOS
    while (p >= s && *p != '\\' && *p != ':') {
	if (*p == '.')
	    *p = '\0';		/* strip off extension */
	if ('A' <= *p && *p <= 'Z')
	    *p += 'a' - 'A';	/* convert to lower case */
	p--;
    }
#else
    while (p >= s && *p != '/')
	p--;
#endif
    return (p+1);
}

static char *help_strings[] = {
    "usage: %s [OPTIONS ...]\n",
 " Takes as standard input an annotation interval list, containing intervals",
 " in seconds and the (beat and non-beat) annotations that terminate each",
 " interval; and outputs on standard output each unique NN increment (x) in",
 " milliseconds, and the percentage of NN interval increments (pNNx) greater",
 " than x.",
 " ",
 "OPTIONS may include:",
 " -h      print this usage summary",
 " -i INC  compute and output pNNx for x = 0, INC, 2*INC, ... milliseconds",
 " -p      compute and output increments as percentage of initial intervals",
 " -s      compute and output separate distributions of positive and negative",
 "          intervals",
 " ",
 "For further information, see http://www.physionet.org/physiotools/pNNx/ .",
NULL
};

void help()
{
    int i;

    (void)fprintf(stderr, help_strings[0], pname);
    for (i = 1; help_strings[i] != NULL; i++)
	(void)fprintf(stderr, "%s\n", help_strings[i]);
}