/* 
file: pnnlist.c

Joe E. Mietus  Feb 24 2003

Usage: pnnlist [options]

List pnn probabilities.
Takes as standard input an RR interval file of
interval in sec and beat annotation and outputs
on standard output each unique NN increment in msec
and the percentage of NN interval increments greater
than that value

[options]
  [-i inc] : output probabilities at each NN increment
  [-p] : use percentage of previous interval
  [-s] : output separate pos and neg signed distributions

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 MAXBUF 172800
#define ABS(A) ((A) < 0 ? -(A) : (A))
#define ROFFERR 1e-9


main(argc, argv)
int argc;
char *argv[];
{

    char ann, lastann[2];
    int i, pflag, sflag, tot, totn, totp, sum, dblcmp();
    double rr, lastrr, inc, nninc, mult;
    static double pnn[MAXBUF];
    void ssort (void  *base, size_t nel, size_t width,
		int (*comp)(const void *, const void *));

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

    for (i=0; ++i < argc && *argv[i] == '-'; ) {
        switch(argv[i][1]) {
 	    case 'i':  if (++i >= argc || (inc = atof(argv[i])) <= 0) {
	                   fprintf(stderr, "%s : -i must be followed ", argv[0]);
	                   fprintf(stderr, "by a number greater than 0\n");
			   exit(1);
	               }
	               break;
	    case 'p':  pflag = 1;
	               mult = 100;
	               break;
	    case 's':  sflag = 1;
	               break;
            default :  usage(argv[0]);
                       exit(1);
	}
    }
    if (inc < 0) {
        fprintf(stderr, "%s -i inc : increment must be greater than 0\n", argv[0]);
	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", argv[0]);
        exit(2);
    }

    i = 0;
    lastrr = rr;
    lastann[0] = 'X';
    lastann[1] = ann;

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

        if (lastann[0] == 'N' && lastann[1] == 'N' && ann == 'N') {
	    
	    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;
    }

    if (i == MAXBUF && scanf("%lf %c", &rr, &ann) == 2)
        fprintf(stderr, "%s : input truncated at %d increments\n",
		argv[0], MAXBUF);

    tot = i;

    ssort(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));
	}
    }
}


usage(prog)
char *prog;
{
        fprintf(stderr, "Usage: %s [options]\n", prog);
	fprintf(stderr, "List pnn probabilities.\n");

	fprintf(stderr, "Takes as standard input an RR interval file of\n");
	fprintf(stderr, "interval in sec and beat annotation and outputs\n");
	fprintf(stderr, "on standard output each unique NN increment in msec\n");
	fprintf(stderr, "and the percentage of NN interval increments greater\n");
	fprintf(stderr, "than that value.\n");
        fprintf(stderr, "  [-i inc] : output probabilities at each NN increment\n");
        fprintf(stderr, "  [-p] : use percentage of previous interval\n");
        fprintf(stderr, "  [-s] : output separate pos and neg signed distributions\n");
}


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

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

}



/* +++Date last modified: 05-Jul-1997 */

/*
**  ssort()  --  Fast, small, qsort()-compatible Shell sort
**
**  by Ray Gardner,  public domain   5/90
*/

/*
 * Manuel Novoa III       Dec 2000
 *
 * There were several problems with the qsort code contained in uClibc.
 * It assumed sizeof(int) was 2 and sizeof(long) was 4.  It then had three
 * seperate quiicksort routines based on the width of the data passed: 2, 4,
 * or anything else <= 128.  If the width was > 128, it returned -1 (although
 * qsort should not return a value) and did no sorting.  On i386 with
 * -Os -fomit-frame-pointer -ffunction-sections, the text segment of qsort.o
 * was 1358 bytes, with an additional 4 bytes in bss.
 *
 * I decided to completely replace the existing code with a small
 * implementation of a shell sort.  It is a drop-in replacement for the
 * standard qsort and, with the same gcc flags as above, the text segment
 * size on i386 is only 183 bytes.
 *
 * Grabbed original file rg_ssort.c from snippets.org.
 * Modified original code to avoid possible overflow in wgap calculation.
 * Modified wgap calculation in loop and eliminated variables gap and wnel.
 */


#include <assert.h>

void ssort (void  *base,
            size_t nel,
            size_t width,
            int (*comp)(const void *, const void *))
{
	size_t wgap, i, j, k;
	char *a, *b, tmp;

	/* Note: still conceivable that nel * width could overflow! */
	assert(width > 0);

	if (nel > 1) {
		for (wgap = 0; ++wgap < (nel-1)/3 ; wgap *= 3) {}
		wgap *= width;
		nel *= width;			/* convert nel to 'wnel' */
		do {
			for (i = wgap; i < nel; i += width) {
				for (j = i - wgap; ;j -= wgap) {
					a = j + ((char *)base);
					b = a + wgap;
					if ( (*comp)(a, b) <= 0 ) {
						break;
					}
					k = width;
					do {
						tmp = *a;
						*a++ = *b;
						*b++ = tmp;
					} while ( --k );
					if (j < wgap) {
						break;
					}
				}
			}
			wgap = (wgap - width)/3;
		} while (wgap);
	}
}