/**---------------------------------------------------------------------------
 ** 
 ** sort.c -- Sorting data in vectors 
 ** 
 ** Author          : Pierre Jaccard
 ** Created On      : 1999/07/15 19:45:50
 ** Last Modified By: Pierre Jaccard
 ** Last Modified On: 1999/07/16 20:12:15
 ** Update Count    : 2
 ** Directory       : /home/pego/pcd1/codas3c/gfi/src/libs/math/
 ** Version         : 0.0
 ** Status          : Unknown
 ** ---------------------------------------------------------------------- ** 
 ** DESCRIPTION: 
 ** 
 **    Functions for sorting data and retrieving ranges in vectors. Some of
 **    these are still under construction. 
 ** 
 ** ---------------------------------------------------------------------- ** 
 ** REVISIONS: 
 ** ---------------------------------------------------------------------- ** 
 ** CHANGES: 
 **------------------------------------------------------------------------**/

#include "sort.h"

/* ------------------------------------------------------------------
	 The function `sort_dir_index()' does not work for negative FLOATs,
	 because of the sign byte. A solution to turn around this problem is
	 proposed here. An offset is applied to the floats so that they are all
   positive.

   RETURNS: non-zero on error
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int sort_float(/* Number of elements in vector */
               int n, 
               /* Place for sort indexes */
               int *indx, 
               /* Pointer to data */
               FLOAT *vec)
#else
int ()
#endif
{

	FLOAT min, offset;
	int i;

	/* Get the minimal value */
	min = MAXFLOAT;
	for(i=0; i<n; i++)
		min = (*(vec+i) < min) ? *(vec+i) : min;
	
	/* Offset all data so that they are all positive */
	offset = (FLOAT) -1.0;
	if(min < 0.0){
		offset = -min + MINFLOAT;
		for(i=0; i<n; i++){
			if(*(vec+i) < ADJ_BADFLOAT)
				*(vec+i) = *(vec+i) + offset;
		}
	}

	/* Sort the data */
	sort_dir_index(n, indx, (CHAR *) vec, (int) sizeof(FLOAT));

	/* Put back the correct values */
	if(offset > 0.0){
		for(i=0; i<n; i++){
			if(*(vec+i) < ADJ_BADFLOAT)
				*(vec+i) = *(vec+i) - offset;
		}
	}
	
	return(0);

}

/* ------------------------------------------------------------------
	 This function is similat to function sort_float(), except that it works fro 
   DOUBLE values.

   RETURNS: non-zero on error
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int sort_double(/* Number of data */
                int n, 
                /* Place for sort index values */
                int *indx, 
                /* Pointer to data */
                DOUBLE *vec)
#else
int ()
#endif
{

	FLOAT *fdata;
	int i, ierr=0;
	
	/* Allocate space for float vector */
	fdata = calloc((size_t) n, sizeof(FLOAT));
	ierr = (!fdata) ? INSUFFICIENT_MEMORY : 0;
	if(DBERROR(&ierr, "calloc(): Memory for FLOAT"))
		return(ierr);

	/* Convert double to float */
	for(i=0; i<n; i++)
		fdata[i] = (FLOAT) vec[i];
	
	/* Sort data */
	sort_float(n, indx, fdata);

	/* Release memory */
	free(fdata);

	return(0);

}

/* ------------------------------------------------------------------
   Merge the depths values from two different vectors, sort them and remove
   duplicate values. 

   RETURNS: number of elements in new vector, or negative on error.
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int merge_vectors_d(/* Number of elements in first vector */
                    int n0, 
                    /* Number of elements in second vector */
                    int n1, 
                    /* Pointer to data in first vector */
                    double *z0, 
                    /* Pointer to data in second vector */
                    double *z1, 
                    /* On output, contains sorted and unique occurences of
                       values appearing in both vectors. On input, memory must 
                       be allocated to hold data from both vectors */
                    double *z)
#else
int ()
#endif
{

	double *tmp=NULL;
	int *indx=NULL;
	
	double t;
	int ierr=0;
	int i, j, n;

	/* Copy all depths */
	j = 0;
	for(i=0; i<n0; i++){
		if(z0[i] < ADJ_BADFLOAT)
			z[j++] = z0[i];
	}
	for(i=0; i<n1; i++){
		if(z1[i] < ADJ_BADFLOAT)
			z[j++] = z1[i];
	}
	n = j;
	
	/* Allocate memory for index and sorted values */
	indx = (int *) calloc((size_t) n, sizeof(int));
	tmp  = (double *) calloc((size_t) n, sizeof(double));
	ierr = ((!tmp) || (!indx)) ? INSUFFICIENT_MEMORY : 0;
	if(DBERROR(&ierr, "calloc()"))
		goto on_error;

	/* Sort depths */
	ierr = sort_double(n, indx, z);
	if(DBERROR(&ierr, "sort_double()"))
		goto on_error;
	for(i=0; i<n; i++)
		tmp[i] = z[indx[i]];

	/* Remove duplicates */
	j = 0;
	t = BADDOUBLE;
	for(i=0; i<n; i++){
		if(abs_val(tmp[i] - t) > MINFLOAT)
			z[j++] = tmp[i];
		t = tmp[i];
	}
	n = j;

on_error:

	if(indx) free(indx);
	if(tmp)  free(tmp);

	if(ierr)
		n = -1;

	return(n);

}

/* ------------------------------------------------------------------
   Create an index vector relating input vector values to another vector.
   On output, the ith element in input vector can be found at the ith index
   value in the other vector. This function does not establish a relation
   for bad values. Index for bad values or values that are not present in
   the other vector are set to EOF.

   RETURNS: number of good indexes or negative on eror
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int index_vector_d(/* Number of elements in first vector */
                   int n0, 
                   /* Number of elements in second vector */
                   int n1, 
                   /* Pointer to data in first vector */
                   double *z0, 
                   /* Pointer to data in second vector */
                   double *z1, 
                   /* Vector of index values. Must be large enough to hold all 
                      data from first vector */
                   int *indx)
#else
int ()
#endif
{

	int i, j, n;

	/* Loop through each input value */
	for(i=0; i<n0; i++){
		indx[i] = EOF;
		if(z0[i] < ADJ_BADFLOAT){
			/* Loop through other vector until the same value is found */
			j=0;
			while((j < n1) && (abs_val(z1[j] - z0[i]) > MINFLOAT)) j++;
			if(j < n1)
				indx[i] = j;
		}
	}
	
	/* Count the number of good values */
	n = 0;
	for(i=0; i<n0; i++){
		if(indx[i] >= 0)
			n++;
	}
	
	return(n);
}

/* ------------------------------------------------------------------
   Build a mask for data inside defined ranges. The range vector contains
   pairs of min and max values definig a range. A mask of the same size as the 
   data vector is created, where elements set to 1 mean that the associated
   element in the data vector is within one of the ranges. 

   RETURNS: number of marks in the mask or negative on error.
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int overlap_vectors_d(/* Number of elements in range vector */
                      int n0, 
                      /* Number of elements in data vector */
                      int n1, 
                      /* Pointer to range data vector */
                      double *z0, 
                      /* Pointer to data vector */
                      double *z1, 
                      /* Pointer to mask data. It must be at least as large as 
                         the data vector. Values of 1 mean that the
                         corresponding data in the data vector is within
                         arange defined in the range vector */
                      int *m)
#else
int ()
#endif
{

	int *indx=NULL;

	int i, j, n;
	int imin, imax;
	int ierr=0;
	
	/* Create a table of index from values in z0 to values in z1 */
	indx = (int *) calloc((size_t) n0, sizeof(int));
	ierr = (!indx) ? INSUFFICIENT_MEMORY : 0;
	if(DBERROR(&ierr, "calloc(): index"))
		goto on_error;
	ierr = (index_vector_d(n0, n1, z0, z1, indx) <= 0) ? INVALID_NUMBER : 0;
	if(DBERROR(&ierr, "index_vector_d"))
		goto on_error;
	
	/* Initialize mask values */
	for(i=0; i<n1; i++)
		m[i] = 0;

	/* Go through the indexes and flag mask elements where values of z1 are
		 within good data ranges in z0 */
	i = 0;
	while(i < n0){
		/* Got to start of index range */
		imin = -1;
		while((i < n0) && (indx[i] < 0)) i++;
		if(i < n0) imin = i;
		/* Go to end of index range */
		while((i < n0) && (indx[i] >= 0)) i++;
		imax = i;
		/* Flag data within index range */
		if(imin >= 0){
			for(j=imin; j<imax; j++)
				m[indx[j]] = 1;
		}
	}

	/* Count the number of mark in mask */
	n = 0;
	for(i=0; i<n1; i++){
		if(m[i] > 0)
			n++;
	}

on_error:

	/* Release memory */
	if(indx) free(indx);

	/* Check for errors */
	if(ierr)
		n = -1;

	return(n);

}


/* ------------------------------------------------------------------
   Find common range of good data within two depths (or whatever)
   vectors. This function is not finished and is not made available at
   present. 
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int common_range_d(int n0, int n1, double *z0, double *z1, int *m0, int *m1)
#else
int ()
#endif
{

	double *z=NULL;
	int    *zm0=NULL, *zm1=NULL;
	int    *i0=NULL, *i1=NULL;
	int    *zm;

	int i, n, j0, j1;
	int ierr=0;

	/* Initialize overlapping masks */
	for(i=0; i<n0; i++)
		m0[i] = 0;
	for(i=0; i<n1; i++)
		m1[i] = 0;
	
	/* Create a vector with all different depths available in the two vectors */
	z  = (double *) calloc((size_t) (n0+n1), sizeof(double));
	ierr = (!z) ? INSUFFICIENT_MEMORY : 0;
	if(DBERROR(&ierr, "calloc(): merged vectors"))
		goto on_error;
	n = merge_vectors_d(n0, n1, z0, z1, z);
	ierr = (n <= 0) ? INVALID_NUMBER : 0;
	if(DBERROR(&ierr, "merge_d()"))
		goto on_error;

	/* Create a mask with defined ranges */
	zm0 = (int *) calloc((size_t) n, sizeof(int));
	zm1 = (int *) calloc((size_t) n, sizeof(int));
	ierr = (!zm0 || !zm1) ? INSUFFICIENT_MEMORY : 0;
	if(DBERROR(&ierr, "calloc(): masks"))
		goto on_error;
	ierr = (overlap_vectors_d(n0, n, z0, z, zm0) <= 0) ? INVALID_NUMBER : 0;
	if(DBERROR(&ierr, "overlap_vectors_d(): z0 on z"))
		goto on_error;
	ierr = (overlap_vectors_d(n1, n, z1, z, zm1) <= 0) ? INVALID_NUMBER : 0;
	if(DBERROR(&ierr, "overlap_vectors_d(): z1 on z"))
		goto on_error;
	zm = zm0;
	for(i=0; i<n; i++)
		zm[i] = (zm0[i] && zm1[i]);
	
	/* Create index array relating merged vector values with input vectors */
	i0 = (int *) calloc((size_t) n, sizeof(int));
	i1 = (int *) calloc((size_t) n, sizeof(int));
	ierr = ((!i0) || (!i1)) ? INSUFFICIENT_MEMORY : 0;
	if(DBERROR(&ierr, "calloc(): index vectors"))
		goto on_error;
	ierr = ( index_vector_d(n, n0, z, z0, i0) > 0 ) ? 0 : INVALID_NUMBER;
	if(DBERROR(&ierr, "index_vector_d(): z to z0"))
		goto on_error;
	ierr = ( index_vector_d(n, n1, z, z1, i1) > 0 ) ? 0 : INVALID_NUMBER;
	if(DBERROR(&ierr, "index_vector_d(): z to z1"))
		goto on_error;

	/* Update overlapping masks */
	for(i=0; i<n ; i++){
		if((zm[i] > 0) && (i0[i] >= 0))
			m0[i0[i]] = 1;
		if((zm[i] > 0) && (i1[i] >= 0))
			m1[i1[i]] = 1;
	}

	/* Calculate the number of common ranges */
	i = j0 = 0;
	while(i < n0){
		while((i < n0) && (m0[i] <= 0)) i++;
		if(i < n0) j0++;
		while((i < n0) && (m0[i] > 0)) i++;
	}
	i = j1 = 0;
	while(i < n1){
		while((i < n1) && (m1[i] <= 0)) i++;
		if(i < n1) j1++;
		while((i < n1) && (m1[i] > 0)) i++;
	}
	ierr = (!j0 || !j1 || (j0 != j1)) ? INVALID_NUMBER : 0;

	/* This will have to be updated to the correct value */
	n = -1;

on_error:

	/* Release memory */
	if(z)   free(z);
	if(zm0) free(zm0);
	if(zm1) free(zm1);
	if(i0)  free(i0);
	if(i1)  free(i1);
	
	if(ierr)
		n = -1;

	return(n);

}

/* ------------------------------------------------------------------
   Build index to start and end range values for vectors. Given a mask with
   data set to 1 for those indexes that should be considered, this function
   fills a start and end vector, with the value of the indexes corresponding
   to the start and the stop of a range of good data.

   RETURNS: number of elements in the start and end vectors, or negative on
            error. 
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int index_range_d(/* number of elements in vectors */
                  int n, 
                  /* Mask with marks set to 1 if data should be considered */
                  int *m, 
                  /* Index to start of ranges */
                  int *start, 
                  /* Index to end of ranges */
                  int *end)
#else
int ()
#endif
{

	int i, j, imin;
	
	/* Initialize start and end arrays */
	for(i=0; i<n; i++)
		start[i] = end[i] = -1;
	
	/* Loop through the mask */
	i = j = 0;
	while(i < n){
		/* Goto start of range */
		imin = -1;
		while((i < n) && (m[i] <= 0)) i++;
		if(i < n) imin = i;
		/* Goto end of range */
		while((i < n) && (m[i] > 0)) i++;
		/* Store start and end indexes */
		if(imin >= 0){
			start[j] = imin;
			end[j]   = i-1;
			j++;
		}
	}

	return(j);

}






