/*---------------------------------------------------------------------------
   
   File name 	:	lwjit_class.c
   Author 		:	trogon
   Created on	:	May 2000
   Description	:	lwjitter classification

 *--------------------------------------------------------------------------*/


/*---------------------------------------------------------------------------
   								Includes
 ---------------------------------------------------------------------------*/

#include "eclipse.h"
#include "isaacp_lib.h"
#include "lwjit_class.h"
#include "lwjit_bb.h"


/*---------------------------------------------------------------------------
   								Defines
 ---------------------------------------------------------------------------*/

#define NEGLIG_OFF_DIFF 	0.1
#define SQR(x) ((x)*(x))

/*---------------------------------------------------------------------------
  							Function codes
 ---------------------------------------------------------------------------*/

static int get_offsets(char *, double *, double *) ;
static int abba_classification(lwjit_bb *,framelist	*) ;

/*---------------------------------------------------------------------------*/
/** 
  	@name lwjit_classification()
	@param bb Black board
	@return	0 if OK -1 if not
	@memo Read the header of the infiles and classify them
	@doc 

	This function filters out half-cycle frames from the input, and
	separates the remaining frames into two batches (chop_a and chop_b).
	The frames are tested for conformance with an abba sequence.
*/
/*---------------------------------------------------------------------------*/
int lwjit_classification(lwjit_bb *bb)
{
	framelist	*	fl ;

	/* Load frames */
	if ((fl = framelist_load(bb->frames_filelist)) == NULL) {
		e_error("loading frame list [%s]", bb->frames_filelist);
		bb->status_classification = ALGO_FAILED ;
		return -1 ;
	}

	/* Remove half-cycle frames */
	isaac_lw_filter_halfcycle(&fl) ;

	/* Odd number of frames ? */
	if (fl->n%2) e_warning("odd number of frames in input [%d]", fl->n) ;

	/* ABBA classification */
	if (abba_classification(bb, fl) == -1) {
		bb->status_classification = ALGO_FAILED ;
		framelist_del(fl);
		return -1 ;
	}
	bb->status_classification = ALGO_OK ;
	
	/* Free and return  */
	framelist_del(fl);
	return 0 ;
}

/*---------------------------------------------------------------------------*/
/** 
  	@name 	get_offsets
	@param	filename 	input fits file 
	@param	x 			pointer to x offset
	@param	y 			pointer to y offset
	@return	0 if OK -1 if not
	@memo	Get cum offsets
	@doc
		Uses isaac_get_cumoffsetx(), isaac_get_cumoffsety()
	@see	isaac_get_cumoffsetx(), isaac_get_cumoffsety()
*/
/*---------------------------------------------------------------------------*/
static int get_offsets(char * filename, double * x, double * y)
{
	char	*	s ;
    
	if ((s = isaac_get_cumoffsetx(filename)) == NULL) return -1 ;
	*x = (double)atof(s) ;
	if ((s = isaac_get_cumoffsety(filename)) == NULL) return -1 ;
	*y = (double)atof(s) ;
	return 0;
}


/*---------------------------------------------------------------------------*/
/** 
  	@name 	abba_classification
	@param 	bb 	Black board
	@param	fl	framelist
	@return	0 if OK -1 if not
	@memo 	Classify frames into ABBA sequences
	@doc 

	This function classifies frames from the input
	into two batches (chop_a and chop_b).
	The algorithm assumes that frames in the input are
	repetitions of ABBA, with possible single frame dropouts
	i.e. ABBAABA. The first frame must be an A frame.
	The algorithm attempts to discriminate the A and B frames
	by finding thw characteristic distance (chop throw) as the 
	most frequent distance between the offsets.
*/
/*---------------------------------------------------------------------------*/
static int abba_classification(
		lwjit_bb 	*	bb,
		framelist	*	fl)
{
	int				na, nb ;
	double 		*	offx,
				*	offy ;
	double			sel_dist ;
	double		*	dist ;
	int			*	dcount ;
	int				i_max,
					j_max,
					max_count ;
	double			throw_x,
					throw_y ;
	int				done ;
	int				change_seq ;

	int				i,
					j ;

	/* Get offsets */	
	offx = calloc(fl->n, sizeof(double)) ;
	offy = calloc(fl->n, sizeof(double)) ;
	for (i=0 ; i<fl->n ; i++) {
		if ((get_offsets(fl->name[i], &(offx[i]), &(offy[i]))) == -1) {
			free(offx);
			free(offx);
			return -1 ;
		}
	}

	/* Find all distances between offsets and count them */	
	dist	= calloc(fl->n * fl->n, sizeof(double)) ;
	dcount  = calloc(fl->n * fl->n, sizeof(int)) ;
	for (i=0 ; i<fl->n ; i++)
		for (j=0 ; j<fl->n ; j++) 
			dist[i+j*fl->n] = sqrt(SQR(offx[i]-offx[j])+SQR(offy[i]-offy[j])) ;

	for (i=0 ; i<fl->n*fl->n ; i++) {
		sel_dist = dist[i] ;
		for (j=0 ; j<fl->n*fl->n ; j++)
			if (fabs(sel_dist-dist[j]) <= NEGLIG_OFF_DIFF) dcount[i] ++ ;
	}

	/* Find Chop offsets as the nonzero distance with maximal occurrence */	
	max_count = 0 ;
	i_max     = 0 ; 
	j_max     = 0 ;
	for (i=0 ; i<fl->n ; i++) 
		for (j=0 ; j<fl->n ; j++) 
			if ((dcount[i+j*fl->n] > max_count) && (dist[i+j*fl->n] > 0.5)) {
				max_count = dcount[i+j*fl->n] ;
				i_max = i ;
				j_max = j ;
			}
	throw_x = offx[j_max] - offx[i_max] ;
	throw_y = offy[j_max] - offy[i_max] ;

	/* Test the throw */
	if ((throw_x == 0) && (throw_y == 0)) {
		e_error("Throw is equal to 0 - cannot classify") ;
		free(offx) ;
		free(offy) ;
		return -1 ;
	}
	
	/* Allocate maximal possible nb of chop frames (worst case: abab)*/
	bb->chop_a = calloc(fl->n-1, sizeof(char*)) ;
	bb->chop_b = calloc(fl->n-1, sizeof(char*)) ;
	
	/* Initialize	 */
	na         = 0 ;
	nb         = 0 ;
	i          = 0 ;
	done       = 0 ;
	change_seq = 0 ;

	while (i+1 < fl->n) {
		/* Consider the successive pairs of frames */
		/* AB ?  */
		if ((fabs(offx[i]-offx[i+1]+throw_x) < NEGLIG_OFF_DIFF) &&
				(fabs(offy[i]-offy[i+1]+throw_y) < NEGLIG_OFF_DIFF)) {
			bb->chop_a[na++] = strdup(fl->name[i]) ;
			bb->chop_b[nb++] = strdup(fl->name[i+1]) ;
			i += 2 ;
		/* BA ? */
		} else if ((fabs(offx[i]-offx[i+1]-throw_x) < NEGLIG_OFF_DIFF) &&
				(fabs(offy[i]-offy[i+1]-throw_y) < NEGLIG_OFF_DIFF)) {
			bb->chop_b[nb++] = strdup(fl->name[i]) ;
			bb->chop_a[na++] = strdup(fl->name[i+1]) ;
			i += 2 ;	
		/* AA or BB ? */
		} else i++ ;
	}
	bb->nnod = bb->nchop = na ;
	
	/* Free and return  */
	free(offx) ;
	free(offy) ;
	free(dist) ;
	free(dcount) ;
	
	if (na != nb) return -1 ;
	return 0 ;
}