/* merge.c

   Compile with:  gcc -g merge.c -lm -o N.data/merge
                  gcc -g merge.c -DANCHOR -lm -o A.data/merge 

   Example run:   merge intv_file1 intv_file2 > intv_file3
                  (see bat_merge for a more contextual example)

   Merges two interval files of anchor values, to produce an interval file
   with a finer interleaving. For instance, merging:
   ::::: tmp1 :::::
   0.001  0.003  15 20 35
   0.004  0.006  12 10 22
   with:
   ::::: tmp2 :::::
   0.001  0.002  10 10 20
   0.003  0.006  20 20 40
   produces:
            0.001  0.002  25   30      55
            0.003  0.003  35   40      75
            0.004  0.006  32   30      62
   NOTE: used in combined region assessments.
*/

#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <ctype.h>
#include <sys/types.h>
#include <sys/param.h>

#define  EPSILON          0.001
#ifndef ANCHOR
#define  LOWER_BOUND        0.00
#define  UPPER_BOUND      100.00
#else 
#define  LOWER_BOUND     -100.00
#define  UPPER_BOUND      100.00
#endif

#define Usage "%s <file1> <file2> [L=] [U=]"

#define max(x,y) (((x)>=(y)) ? (x) : (y))
#define min(x,y) (((x)<=(y)) ? (x) : (y))

void fatal(char *message);
void fatalf(char *msg, char *val);
FILE *ckopen(char *name, char *mode);

int  read_file(char *,double *,double *,int *,int *);
int  merge(double *,double *,double *,double *,double *,double *,
           int *,int *,int *,int *,int *,int *,int,int);
void print_results(double *,double *,int *,int *,int,int,int);
int  comparable(double,double,double,double);


int main(int argc, char *argv[])
{
    int    n1, n2, low, up, top;
    double From1[1000], From2[1000], To1[1000], To2[1000];
    double From[2000], To[2000];
    int    FP1[1000], FP2[1000], FN1[1000], FN2[1000];
    int    FP[2000], FN[2000];
    
    low = LOWER_BOUND;
    up = UPPER_BOUND;
    if ((argc<3) || (argc>5)) fatalf(Usage, argv[0]);
    while (argc>3) {
        if (!strncmp(argv[argc-1],"L=",2))
            low = atoi(argv[argc-1]+2);
        else if (!strncmp(argv[argc-1],"U=",2))        
            up = atoi(argv[argc-1]+2);
        argc--;
    } 
    n1 = read_file(argv[1], From1, To1, FP1, FN1);
    n2 = read_file(argv[2], From2, To2, FP2, FN2);
    top = merge(From1, From2, From, To1, To2, To, FP1, FP2, FP, FN1, FN2, FN, n1, n2); 

    print_results(From, To, FP, FN, top, low, up);

    return 0;
}

void print_results(double *from, double *to, int *fp, int *fn, int top, int low, int up)
{
     int i, save_i;
     char tmp[50];

     i = 0;
     while ((i<top) && (to[i]<low))  i++;
     if (i==top) { 
         sprintf(tmp, "No intervals in the range %d %d.", low, up);
         fatal(tmp);
     }
     printf("%-3.3f  %-3.3f  %-3d  %-3d    %3d\n", 
            max(from[i],low), min(to[i],up), fp[i], fn[i], fp[i]+fn[i]);
     save_i = i;
     ++i;
     while ((i<top) && (to[i]<=up))  {
         printf("%-3.3f  %-3.3f  %-3d  %-3d    %3d\n",
                from[i], to[i], fp[i], fn[i], fp[i]+fn[i]);
         i++;
     }
     
     if ((i<top) && (save_i!=i))
     printf("%-3.3f  %-3.3f  %-3d  %-3d    %3d\n",
            from[i], min(to[i],up), fp[i], fn[i], fp[i]+fn[i]);

     return;
}
int read_file(char *filename, double *from, double *to, int *fp, int *fn)
{ 
    FILE *file;
    int   count;
    char  buffer[200];

    file = ckopen(filename,"r");

    count = 0;
    while (fgets(buffer,200,file)!=NULL) {
       if (buffer[0]=='#') continue;
       sscanf(buffer,"%lf %lf %d %d", 
              &from[count], &to[count], &fp[count], &fn[count]);
       count++;
    }
    fclose(file);

    return count;
}
   
/* may alter the arrays, temporarily */ 
int merge(double *from1, double *from2, double *from, 
           double *to1, double *to2, double *to,
           int *fp1, int *fp2, int *fp, 
           int *fn1, int *fn2, int *fn, int M, int N)
{
     int    i, j, k;
     int    poscode;
 
     i = j = k = 0;
     while ((i<M) && (to1[i]<from2[j])) i++;
     if (i==M) fatal("No common intervals.");
     while ((j<N) && (to2[j]<from1[i])) j++;
     if (j==N) fatal("No common intervals.");

     while ((i<M) && (j<N) &&
            (poscode=comparable(from1[i], to1[i], from2[j], to2[j]))!=0) {
         switch (poscode) {
             case  1212: /* overlap 1, 2 */
                   from[k] = from2[j]; 
                   to[k] = to1[i];
                   fp[k] = fp1[i]+fp2[j];
                   fn[k] = fn1[i]+fn2[j];
                   k++; i++;
                   break;
             case  1221: /* 1 contains 2 */
                   from[k] = from2[j];
                   to[k] = to2[j];
                   fp[k] = fp1[i]+fp2[j];
                   fn[k] = fn1[i]+fn2[j];
                   k++; 
                   if (to1[i]!=to2[j]) j++; 
                   else { i++; j++; }
                   break;
             case  2121: /* overlap 2, 1 */
                   from[k] = from1[i];
                   to[k] = to2[j];
                   fp[k] = fp1[i]+fp2[j]; 
                   fn[k] = fn1[i]+fn2[j];
                   k++; j++;
                   break;
             case  2112: /* 2 contains 1 */
                   from[k] = from1[i]; 
                   to[k] = to1[i]; 
                   fp[k] = fp1[i]+fp2[j];
                   fn[k] = fn1[i]+fn2[j];
                   k++;
                   if (to1[i]!=to2[j]) i++;
                   else { i++; j++; }
                   break; 
            default: fatal("Unrecognized code.");  
         } 
     }
     return k;      
}

int comparable(double from1, double to1, double from2, double to2) 
{
    if ((from1<=from2) && (to2<=to1)) return 1221;
    if ((from2<=from1) && (to1<=to2)) return 2112; 
    if ((from2<=to1) && (from1<=from2)) return 1212;
    if ((from1<=to2) && (from2<=from1)) return 2121;

    return -1;
}


void fatal(char *msg)
{                
        (void)fprintf(stderr, "%s\n", msg);
        exit(1); 
}    
     
void fatalf(char *msg, char *val)
{    
        (void)fprintf(stderr, msg, val);  
        (void)putc('\n', stderr);
        exit(1);
}            
          
FILE *ckopen(char *name, char *mode)
{    
        FILE *fp;
     
        if ((fp = fopen(name, mode)) == NULL)
                fatalf("Cannot open %s.", name);
        return fp;
}