/*********************************************************************
 *                                                                   *
 * COMPRESSSOURCES.C  -  Compresses sources                          *
 * (c) 2004 by Bernd Vollmer, CDS, Observatoire de Strasbourg        *
 *                                                                   *
 *********************************************************************/
#include <stdio.h>
#include "nrutil.h" 
#include <math.h>
#include "code.h"

local int fcompare(bodyptr,bodyptr);

void compresssource(bodyptr p)
{
  int i,j, k, pl, ps, pb, strc;
  real sum, error, sum1, error1, factor, freq;
  comppareptr q;
  bodyptr ppp[maxparents],sss[maxsons],bbb[maxbrothers], pq;
  short pflag[maxparents],bflag[maxbrothers],sflag[maxsons];

  for (i=0;i<Nbrothers(p);i++) {             /* remember brothers */
    bflag[i]=1;
    bbb[i]=(bodyptr) Brothers(p)[i];
  }
  for (i=0;i<Nparents(p);i++) {              /* remember parents */ 
    pflag[i]=1;
    ppp[i]=(bodyptr) Parents(p)[i];
  }
  for (i=0;i<Nsons(p);i++) {                 /* remember sons */
    sflag[i]=1;
    sss[i]=(bodyptr) Sons(p)[i];
  }
  pl=ps=pb=0;
  
  if (Nbrothers(p) > 0) {                    /* check brothers */ 
    freq=Frequency(p);
    for (i=0;i<Nbrothers(p);i++) {
      if (Frequency(Brothers(p)[i]) != freq) {  /* not at the same frequency */
	if (bflag[i] == 1) {
	  bflag[i]=2;
	  pq=(bodyptr) Brothers(p)[i];
	  if (Nbrothers(p) > 1) {
	    for (j=0;j<Nbrothers(p);j++) {
	      if (Frequency(Brothers(p)[j]) == Frequency(pq)) {
		if (Nparents(pq) > 0 ) {
		  for (k=0;k<Nparents(pq);k++) {
		    strc=strcmp(Name(Brothers(p)[j]),Name(Parents(pq)[k])); /* check for parent of brother = brother */
		    if (strc == 0 && bflag[j] == 1) bflag[j]=0;
		  }
		}
		if (Nbrothers(pq) > 0 ) {
		  for (k=0;k<Nbrothers(pq);k++) {
		    strc=strcmp(Name(Brothers(p)[j]),Name(Brothers(pq)[k])); /* check for brother of brother = brother */
		    if (strc == 0 && bflag[j] == 1) bflag[j]=0;
		  }
		}
		if (Nsons(pq) > 0 ) {
		  for (k=0;k<Nsons(pq);k++) {
		    strc=strcmp(Name(Brothers(p)[j]),Name(Sons(pq)[k])); /* check for son of brother = brother */
		    if (strc == 0 && bflag[j] == 1) bflag[j]=0;
		  }
		}
	      }
	    }
	  }
	  
	  if (Nparents(p) > 0) {              
	    for (j=0;j<Nparents(p);j++) {
	      if (Frequency(Parents(p)[j]) == Frequency(pq)) {
		if (Nparents(pq) > 0 ) {
		  for (k=0;k<Nparents(pq);k++) {
		    strc=strcmp(Name(Parents(p)[j]),Name(Parents(pq)[k])); /* check for parent of brother = parent */
		    if (strc == 0 && pflag[j] == 1) pflag[j]=0;
		  }
		}
		if (Nbrothers(pq) > 0 ) {
		  for (k=0;k<Nbrothers(pq);k++) {
		    strc=strcmp(Name(Parents(p)[j]),Name(Brothers(pq)[k])); /* check for brother of bother = parent */
		    if (strc == 0 && pflag[j] == 1) pflag[j]=0;
		  }
		}
		if (Nsons(pq) > 0 ) {
		  for (k=0;k<Nsons(pq);k++) {
		    strc=strcmp(Name(Parents(p)[j]),Name(Sons(pq)[k])); /* check for son of brother = parent */
		    if (strc == 0 && pflag[j] == 1) pflag[j]=0;
		  }
		}
	      }
	    }
	  }
	  
	  if (Nsons(p) > 0) {
	    for (j=0;j<Nsons(p);j++) {
	      if (Frequency(Sons(p)[j]) == Frequency(pq)) {
		if (Nparents(pq) > 0 ) {
		  for (k=0;k<Nparents(pq);k++) {
		    strc=strcmp(Name(Sons(p)[j]),Name(Parents(pq)[k])); /* check for parent of brother = son */
		    if (strc == 0 && sflag[j] == 1) sflag[j]=0;
		  }
		}
		if (Nbrothers(pq) > 0 ) {
		  for (k=0;k<Nbrothers(pq);k++) {
		    strc=strcmp(Name(Sons(p)[j]),Name(Brothers(pq)[k])); /* check for brother of brother =son */
		    if (strc == 0 && sflag[j] == 1) sflag[j]=0;
		  }
		}
		if (Nsons(pq) > 0 ) {
		  for (k=0;k<Nsons(pq);k++) {
		    strc=strcmp(Name(Sons(p)[j]),Name(Sons(pq)[k])); /* check for son of brother = son */
		    if (strc == 0 && sflag[j] == 1) sflag[j]=0;
		  }
		}
	      }
	    }
	  }
	}
      }
    }
  }


  if (Nparents(p) > 0) {                   /* check parents */
    freq=Frequency(p);
    for (i=0;i<Nparents(p);i++) {
      if (Frequency(Parents(p)[i]) != freq) {
	if (pflag[i] == 1) {
	  pflag[i]=2;	
	  pq=(bodyptr) Parents(p)[i];
	  if (Nbrothers(p) > 0) {
	    for (j=0;j<Nbrothers(p);j++) {
	      if (Frequency(Brothers(p)[j]) == Frequency(pq)) {
		if (Nparents(pq) > 0 ) {
		  for (k=0;k<Nparents(pq);k++) {
		    strc=strcmp(Name(Brothers(p)[j]),Name(Parents(pq)[k])); /* check for parent of parent = brother */
		    if (strc == 0 && bflag[j] == 1) bflag[j]=0;
		  }
		}
		if (Nbrothers(pq) > 0 ) {
		  for (k=0;k<Nbrothers(pq);k++) {
		    strc=strcmp(Name(Brothers(p)[j]),Name(Brothers(pq)[k])); /* check for brother of parent = brother */
		    if (strc == 0 && bflag[j] == 1) bflag[j]=0;
		  }
		}
		if (Nsons(pq) > 0 ) {
		  for (k=0;k<Nsons(pq);k++) {
		    strc=strcmp(Name(Brothers(p)[j]),Name(Sons(pq)[k])); /* check for son of parent = brother */
		    if (strc == 0 && bflag[j] == 1) bflag[j]=0;
		  }
		}
	      }
	    }
	  }
	  
	  if (Nparents(p) > 1) {
	    for (j=0;j<Nparents(p);j++) {
	      if (Frequency(Parents(p)[j]) == Frequency(pq)) {
		if (Nparents(pq) > 0 ) {
		  for (k=0;k<Nparents(pq);k++) {
		    strc=strcmp(Name(Parents(p)[j]),Name(Parents(pq)[k])); /* check for parent of parent = parent */
		    if (strc == 0 && pflag[j] == 1) pflag[j]=0;
		  }
		}
		if (Nbrothers(pq) > 0 ) {
		  for (k=0;k<Nbrothers(pq);k++) {
		    strc=strcmp(Name(Parents(p)[j]),Name(Brothers(pq)[k])); /* check for brother of parent = parent */
		    if (strc == 0 && pflag[j] == 1) pflag[j]=0;
		  }
		}
		if (Nsons(pq) > 0 ) {
		  for (k=0;k<Nsons(pq);k++) {
		    strc=strcmp(Name(Parents(p)[j]),Name(Sons(pq)[k])); /* check for son of parent = parent */
		    if (strc == 0 && pflag[j] == 1) pflag[j]=0;
		  }
		}
	      }
	    }
	  }
	  
	  if (Nsons(p) > 0) {
	    for (j=0;j<Nsons(p);j++) {
	      if (Frequency(Sons(p)[j]) == Frequency(pq)) {
		if (Nparents(pq) > 0 ) {
		  for (k=0;k<Nparents(pq);k++) {
		    strc=strcmp(Name(Sons(p)[j]),Name(Parents(pq)[k])); /* check for parent of parent = son */
		    if (strc == 0 && sflag[j] == 1) sflag[j]=0;
		  }
		}
		if (Nbrothers(pq) > 0 ) {
		  for (k=0;k<Nbrothers(pq);k++) {
		    strc=strcmp(Name(Sons(p)[j]),Name(Brothers(pq)[k])); /* check for brother of parent = son */
		    if (strc == 0 && sflag[j] == 1) sflag[j]=0;
		  }
		}
		if (Nsons(pq) > 0 ) {
		  for (k=0;k<Nsons(pq);k++) {
		    strc=strcmp(Name(Sons(p)[j]),Name(Sons(pq)[k])); /* check for son of parent = son */
		    if (strc == 0 && sflag[j] == 1) sflag[j]=0;
		  }
		}
	      }
	    }
	  }
	}
      }
    }
  }
  
  
  if (Nsons(p) > 0) {                    /* check sons */
    freq=Frequency(p);
    for (i=0;i<Nsons(p);i++) {
      if (Frequency(Sons(p)[i]) != freq) {
	if (sflag[i] == 1) {
	  sflag[i]=2;	
	  pq=(bodyptr) Sons(p)[i];
	  if (Nbrothers(p) > 0) {
	    for (j=0;j<Nbrothers(p);j++) {
	      if (Frequency(Brothers(p)[j]) == Frequency(pq)) {
		if (Nparents(pq) > 0 ) {
		  for (k=0;k<Nparents(pq);k++) {
		    strc=strcmp(Name(Brothers(p)[j]),Name(Parents(pq)[k])); /* check for parent of son = brother */
		    if (strc == 0 && bflag[j] == 1) bflag[j]=0;
		  }
		}
		if (Nbrothers(pq) > 0 ) {
		  for (k=0;k<Nbrothers(pq);k++) {
		    strc=strcmp(Name(Brothers(p)[j]),Name(Brothers(pq)[k])); /* check for brother of son = brother */
		    if (strc == 0 && bflag[j] == 1) bflag[j]=0;
		  }
		}
		if (Nsons(pq) > 0 ) {
		  for (k=0;k<Nsons(pq);k++) {
		    strc=strcmp(Name(Brothers(p)[j]),Name(Sons(pq)[k])); /* check for son of son = brother */
		    if (strc == 0 && bflag[j] == 1) bflag[j]=0;
		  }
		}
	      }
	    }
	  }
	  
	  if (Nparents(p) > 0) {
	    for (j=0;j<Nparents(p);j++) {
	      if (Frequency(Parents(p)[j]) == Frequency(pq)) {
		if (Nparents(pq) > 0 ) {
		  for (k=0;k<Nparents(pq);k++) {
		    strc=strcmp(Name(Parents(p)[j]),Name(Parents(pq)[k])); /* check for parent of son = parent */
		    if (strc == 0 && pflag[j] == 1) pflag[j]=0;
		  }
		}
		if (Nbrothers(pq) > 0 ) {
		  for (k=0;k<Nbrothers(pq);k++) {
		    strc=strcmp(Name(Parents(p)[j]),Name(Brothers(pq)[k])); /* check for brother of son = parent */
		    if (strc == 0 && pflag[j] == 1) pflag[j]=0;
		  }
		}
		if (Nsons(pq) > 0 ) {
		  for (k=0;k<Nsons(pq);k++) {
		    strc=strcmp(Name(Parents(p)[j]),Name(Sons(pq)[k])); /* check for son of son = parent */
		    if (strc == 0 && pflag[j] == 1) pflag[j]=0;
		  }
		}
	      }
	    }
	  }
	  
	  if (Nsons(p) > 1) {
	    for (j=0;j<Nsons(p);j++) {
	      if (Frequency(Sons(p)[j]) == Frequency(pq)) {
		if (Nparents(pq) > 0 ) {
		  for (k=0;k<Nparents(pq);k++) {
		    strc=strcmp(Name(Sons(p)[j]),Name(Parents(pq)[k])); /* check for parent of son = son */
		    if (strc == 0 && sflag[j] == 1) sflag[j]=0;
		  }
		}
		if (Nbrothers(pq) > 0 ) {
		  for (k=0;k<Nbrothers(pq);k++) {
		    strc=strcmp(Name(Sons(p)[j]),Name(Brothers(pq)[k])); /* check for brother of son = son */
		    if (strc == 0 && sflag[j] == 1) sflag[j]=0;
		  }
		}
		if (Nsons(pq) > 0 ) {
		  for (k=0;k<Nsons(pq);k++) {
		    strc=strcmp(Name(Sons(p)[j]),Name(Sons(pq)[k])); /* check for son of son = son */
		    if (strc == 0 && sflag[j] == 1) sflag[j]=0;
		  }
		}
	      }
	    }
	  }
	}
      }
    }
    }
  
  j=0;
  for (i=0;i<Nparents(p);i++) {              /* reconfigure parents */
    if (pflag[i] > 0) {
      Parents(p)[j]=(nodeptr) ppp[i];
      j++;
    } else ncompress++;
  }
  Nparents(p)=j;

  j=0;
  for (i=0;i<Nbrothers(p);i++) {             /* reconfigure brothers */  
    if (bflag[i] > 0) {
      Brothers(p)[j]=(nodeptr) bbb[i];
      j++;
    } else ncompress++;
  }
  Nbrothers(p)=j;

  j=0;
  for (i=0;i<Nsons(p);i++) {                 /* reconfigure sons */
    if (sflag[i] > 0) {
      Sons(p)[j]=(nodeptr) sss[i];
      j++;
    } else ncompress++;
  }
  Nsons(p)=j;




  pl=0;
  if (Nparents(p) > 1) {
    for (i=0;i<Nparents(p);i++) {
      pflag[i]=1;
      ppp[i]=(bodyptr) Parents(p)[i];
    }
    j=0;
    for (i=0;i<Nparents(p);i++) {            /* search for double parents */
      for (j=i+1;j<Nparents(p);j++) {
	strc=strcmp(Name(Parents(p)[i]),Name(Parents(p)[j]));
	if (strc == 0) {
	  pflag[i]=0;
	  pl++;
	}
      }
    }
    Nparents(p)=j;
  }

  if (Nbrothers(p) > 1) {
    for (i=0;i<Nbrothers(p);i++) {           /* search for double brothers */ 
      bflag[i]=1;
      bbb[i]=(bodyptr) Brothers(p)[i];
    }
    j=0;
    for (i=0;i<Nbrothers(p);i++) {
      for (j=i+1;j<Nbrothers(p);j++) {
	strc=strcmp(Name(Brothers(p)[i]),Name(Brothers(p)[j]));
	if (strc == 0) {
	  bflag[i]=0;
	  pl++;
	}
      }
    }
    Nbrothers(p)=j;
  }

  if (Nsons(p) > 1) {                        /* search for double sons */
    for (i=0;i<Nsons(p);i++) {
      sflag[i]=1;
      sss[i]=(bodyptr) Sons(p)[i];
    }
    j=0;
    for (i=0;i<Nsons(p);i++) {
      for (j=i+1;j<Nsons(p);j++) {
	strc=strcmp(Name(Sons(p)[i]),Name(Sons(p)[j]));
	if (strc == 0) {
	  sflag[i]=0;
	  pl++;
	}
      }
    }
    Nsons(p)=j;
  }

   j=0;
   for (i=0;i<Nparents(p);i++) {                 /* reconfigure parents */
    if (pflag[i] > 0) {
      Parents(p)[j]=(nodeptr) ppp[i];
      j++;
    } else ncompress++;
  }
  Nparents(p)=j;

  j=0;
  for (i=0;i<Nbrothers(p);i++) {                 /* reconfigure brothers */
    if (bflag[i] > 0) {
      Brothers(p)[j]=(nodeptr) bbb[i];
      j++;
    } else ncompress++;
  }
  Nbrothers(p)=j;

  j=0;
  for (i=0;i<Nsons(p);i++) {                     /* reconfigure sons */
    if (sflag[i] > 0) {
      Sons(p)[j]=(nodeptr) sss[i];
      j++;
    } else ncompress++;
  }
  Nsons(p)=j;
  
  if (pl > 0) {
    ndouble++;                                   /* number of doubles */ 
  }

}



/*
 *  quicksort subroutine
 */ 
	
	 
local int fcompare(bodyptr a,bodyptr b)
{
  if ((Compfrequency(a)-Compfrequency(b)) > 0.0) return(1);
  else if ((Compfrequency(a)-Compfrequency(b)) < 0.0) return(-1);
  else return(0);
}