/** File gg_featext.c
 *  Transcripted from Motif to GTK+ by P. Vincent to
 *  Replace featext.c
 */


#include <config.h>
#include <cmath.h>
#include <gtk/gtk.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "globals.h"
#include "utils.h"
#include "graphs.h"
#include "noxprotos.h"

#include "gw_list.h"

#include "gg_gutil.h"
#include "gg_gtkinc.h"
#include "gg_protos.h"
#include "fourier.h"


static gint popup_width = 350  ,list_height = 250;

static void gg_featext_AAC_CB 	(GtkWidget *w ,gint reponse);
void gg_fext_routine( int gto, int feature, int abs_src, int abs_set, int abs_graph);
int gg_getmedian (int grno, int setno, int sorton, double *median);
int gg_get_zero_crossing (int setl, double *xv, double *yv, double *crossing);
int gg_get_rise_time (int setl, double *xv, double *yv, double min, double max,
		      double *width);
int gg_get_fall_time( int setl, double *xv, double *yv, double min, double max,
		      double *width );
int gg_mute_linear_regression(int n, double *x, double *y, double *slope, 
			      double *intercept);
int gg_get_half_max_width(int len, double *x, double *y, double *width);
int gg_get_barycenter( int n, double *x, double *y, double *barycenter );
void gg_get_max_pos( double *a, double *b, int n, double max, double *d );
double gg_linear_interp( double x1, double y1, double x2, double y2, double x );
int gg_dbl_comp( const void *a, const void *b );


typedef struct _gg_Featext_ui {
  GtkWidget *top;
  GtkWidget *tograph;
  GraphMenu3Struct *graph_menu3;
  GtkWidget *feature;
  GtkWidget *xval;
  GtkWidget *absic_graph;
  GraphMenu3Struct *absic_graph_menu3;
  GtkWidget *absic_set;
  SetsMenu3Struct  *absic_sets_menu3;
  GtkWidget *legload_rc;
} gg_Featext_ui;

static gg_Featext_ui feui;

/**
 *  Replace   do_fext_toggle   featext.c   85
 */
void gg_fext_toggle_CB (GtkWidget *w ,gpointer p)
{
  int n = gg_get_int (w);
  
  if (n ==2 || n == 3) {
    gtk_widget_set_sensitive (feui.absic_graph ,TRUE);
    gtk_widget_set_sensitive (feui.absic_set   ,TRUE);
      } else {
    gtk_widget_set_sensitive (feui.absic_graph ,FALSE);
    gtk_widget_set_sensitive (feui.absic_set   ,FALSE);
  }
}

void gg_update_absic_set (void)
{
  printf ("gg_update_absic_set\n");
  gg_update_all_sets_lists ();

}

/**
 *  Replace   create_featext_frame      featext.c 95
 */
void gg_create_featext_popup (void)
{
  GtkWidget *vb ,*sep;

  if (feui.top == NULL) {
    feui.top     = gg_CreateAACDialog ("Feature Extraction" ,gg_featext_AAC_CB ,0 ,0 ,&vb);
    feui.tograph = gg_CreateGraphChoice ("Results to graph:" ,GTK_SELECTION_BROWSE
					 ,0 ,0 ,NULL
					 ,&(feui.graph_menu3));
    gtk_box_pack_start_defaults (GTK_BOX (vb) ,feui.tograph);
    feui.feature = gg_combo_new (vb ,"Feature:" ,25
					 , "Y minimum"               /*  0 */
					 , "Y maximum"       	     /*  1 */
					 , "Y average"       	     /*  2 */
					 , "Y std. dev."     	     /*  3 */
					 , "Y median"        	     /*  4 */
					 , "X minimum"       	     /*  5 */
					 , "X maximum"       	     /*  6 */
					 , "X average"       	     /*  7 */
					 , "X std. dev."     	     /*  8 */
					 , "X median"        	     /*  9 */
					 , "Frequency"       	     /* 10 */
					 , "Period"          	     /* 11 */
					 , "Zero crossing"   	     /* 12 */
					 , "Rise time"       	     /* 13 */
					 , "Fall time"       	     /* 14 */
					 , "Slope"           	     /* 15 */
					 , "Y intercept"     	     /* 16 */
					 , "Set length"      	     /* 17 */
					 , "Half maximal width"	     /* 18 */
					 , "Barycenter X"    	     /* 19 */
					 , "Barycenter Y"    	     /* 20 */
					 , "X(Y max)"        	     /* 21 */
					 , "Y(X max)"        	     /* 22 */
					 , "Integral");              /* 23 */

    feui.xval    = gg_combo_new (vb ,"X values from:" ,5
					 ,"Index"
					 ,"Legends"
					 ,"X from Set"
					 ,"Y from set");
					
    sep = gtk_hseparator_new ();
    gtk_box_pack_start (GTK_BOX (vb) ,sep ,FALSE ,TRUE ,1);	

    feui.absic_graph = gg_CreateGraphChoice ("Abscissa from graph:" ,GTK_SELECTION_BROWSE
					     ,0 ,0 ,NULL
					     ,&(feui.absic_graph_menu3));
    gtk_box_pack_start_defaults (GTK_BOX (vb) ,feui.absic_graph);

    feui.absic_sets_menu3 = NULL;
    feui.absic_set   = gg_CreateSetChoice   ("set:" ,GTK_SELECTION_BROWSE
					     ,popup_width ,list_height
					     ,NULL
					     ,&(feui.absic_sets_menu3)
					     ,NULL);
    gtk_box_pack_start_defaults (GTK_BOX (vb) ,feui.absic_set);

    gg_fext_toggle_CB (feui.xval ,NULL);
    g_signal_connect (feui.xval        ,"changed"       ,G_CALLBACK (gg_fext_toggle_CB) ,NULL);
    g_signal_connect (feui.absic_graph ,"gw_list_event" ,G_CALLBACK (gg_update_absic_set) ,NULL);
  }
  gtk_widget_show_all (feui.top);

}

/**
 *  Replace  do_fext_proc    	    featext.c 180
 */
static void gg_featext_AAC_CB 	(GtkWidget *w ,gint reponse)
{
  int gto ,feature ,abs_graph = -1 ,abs_set = -1 ,abs_src;

  if (reponse == GTK_RESPONSE_ACCEPT || reponse == GTK_RESPONSE_APPLY ) {  
    feature   = gg_get_int (feui.feature);
    if (gw_list_count_rows_selected (GW_LIST(feui.tograph)) < 1) {
       gto = get_cg();
    } else {
      gto       = gg_get_int (feui.tograph);
    }
    abs_src   = gg_get_int (feui.xval);
    if (abs_src == 2 || abs_src == 3) {
      abs_set   = gg_get_int (feui.absic_set);
      abs_graph = gg_get_int (feui.absic_graph);
    }
    gg_fext_routine (gto ,feature ,abs_src ,abs_set ,abs_graph); 
    gg_update_all ();
    gg_drawgraph ();
  }
  if (reponse == GTK_RESPONSE_ACCEPT || reponse == GTK_RESPONSE_CLOSE) {
    gtk_widget_hide (w);
  }
}

/**
 * Replace  fext_routine    	    featext.c 201
 */
void gg_fext_routine (int gto, int feature, int abs_src, int abs_set, int abs_graph )
{
	int i, cs, ns, fts, ncurves, extract_err;
	double datum, dummy, *absy;
	double y1, y2;
	int iy1, iy2;
	char tbuf[1024];
	float *abscissa;
        double xmin, xmax, ymin, ymax;
        int cg = get_cg();
        double *x;

	abscissa = xmalloc(number_of_sets(cg)*SIZEOF_FLOAT);
	
	if( !is_graph_active( gto )	){
		gg_errwin("Graph for results must be active");
	    return;
	}
	if( (ns=nextset( gto ) )== -1 ) {
		gg_errwin("gg_fext_routine: can't create new set");
	    return;
	}
	ncurves = nactive(get_cg());
	switch( abs_src ) {
		case 0:		/* use index */
			for( i=0; i<ncurves; i++ )
				abscissa[i] = i;
			break;	
		case 1:		/* use legend label */
			cs = 0;
			for( i=0; i<ncurves; i++ ){
				while( !is_set_active( get_cg(), cs ) )
					cs++;
				if(!sscanf( get_legend_string(get_cg(), cs), "%f", &abscissa[i]))
					break;
				cs++;
			}
			if( i != ncurves ) {
				gg_errwin("Bad legend label");
				return;
			}
			break;
		case 2:		/* use X from set */
			if( !is_set_active( abs_graph, abs_set ) ){
	    		gg_errwin("Abscissa set not active");
	    		return;
			}
			if( getsetlength( abs_graph, abs_set ) < ncurves ) {
				gg_errwin("Not enough points in set");
				return;
			}
			absy = getx( abs_graph, abs_set );
			for( i=0; i<ncurves; i++ )
				abscissa[i] = absy[i];
			break;			
		case 3:										/* use Y from set */
			if( !is_set_active( abs_graph, abs_set ) ){
	    		gg_errwin("Abscissa set not active");
	    		return;
			}
			if( getsetlength( abs_graph, abs_set ) < ncurves ) {
				gg_errwin("Not enough points in set");
				return;
			}
			absy = gety( abs_graph, abs_set );
			for( i=0; i<ncurves; i++ )
				abscissa[i] = absy[i];
			break;
	}

	cs = 0;
	tbuf[0] = '\0';
	for( i=0; i<ncurves; i++ ) {
		while( !is_set_active( get_cg(), cs ) )
			cs++;
		extract_err = 0;
			
		getsetminmax(get_cg(), cs, &xmin, &xmax, &ymin, &ymax);
                switch( feature ) {
			case 0:			/* Y minimum */
				datum = ymin;		
				break;
			case 1: 		/* Y maximum */
				datum = ymax;		
				break;
			case 2: 		/* Y mean    */
				stasum(gety(cg, cs), getsetlength(cg, cs), &datum, &dummy);
				break;
			case 3:			/* Y std dev */
				stasum(gety(cg, cs), getsetlength(cg, cs), &dummy, &datum);
				break;
			case 4: 		/* Y median  */
				gg_getmedian( cg, cs, DATA_Y, &datum );
				break;
			case 5:			/* X minimum */
				datum = xmin;		
				break;
			case 6: 		/* X maximum */
				datum = xmax;		
				break;
			case 7: 		/* X mean    */
				stasum(getx(cg, cs), getsetlength(cg, cs), &datum, &dummy);
				break;
			case 8:			/* X std dev */
				stasum(getx(cg, cs), getsetlength(cg, cs), &dummy, &datum);
				break;
			case 9:			/* X median  */
				gg_getmedian( cg, cs, DATA_X, &datum );
				break;
			case 10: 		/* frequency and period */
			case 11:
			  do_fourier_featext (cg ,cs ,&fts);
			  minmax (gety (cg ,fts) ,getsetlength(cg ,fts),&y1,&y2,&iy1,&iy2);
			  x = getx(cg ,fts);
			  if( feature == 10 )
			    datum = x[iy2];
			  else
			    datum = 1./x[iy2];
			  killset( cg ,fts );				/* get rid of Fourier set */
			  break;
			case 12:		/* first zero crossing */
				if( gg_get_zero_crossing( getsetlength( cg, cs ), 
							  getx( cg, cs ),gety( cg, cs ), &datum ) ){
					sprintf( tbuf+strlen(tbuf), 
						 "Unable to find zero crossing of set %d\n", cs );
					gg_errwin( tbuf );
					extract_err = 1;
				}
				break;
			case 13:		/* rise time   */
				getsetminmax(cg, cs, &xmin, &xmax, &ymin, &ymax);
                                if( gg_get_rise_time( getsetlength(cg,cs), getx(cg,cs), 
					gety(cg,cs), ymin, ymax, &datum ) ){
					sprintf( tbuf+strlen(tbuf), 
							"Unable to find rise time of set %d\n", cs );
					gg_errwin( tbuf );
					extract_err = 1;
				}
				break;
			case 14: 		/* fall time   */
				getsetminmax(cg, cs, &xmin, &xmax, &ymin, &ymax);
                                if (gg_get_fall_time( getsetlength(cg,cs), getx(cg,cs), 
					gety(cg,cs), ymin, ymax, &datum ) ){
					sprintf( tbuf+strlen(tbuf), 
									"Unable to find fall time of set %d\n", cs );
					extract_err = 1;
					gg_errwin( tbuf );
				}
				break;
			case 15:		/* slope       */
				if (gg_mute_linear_regression( getsetlength( cg, cs ), 
					getx( cg, cs ),gety( cg, cs ), &datum, &dummy ) ) {
					sprintf( tbuf+strlen(tbuf), 
										"Unable to find slope of set %d\n", cs );
					gg_errwin( tbuf );
					extract_err = 1;
				}
				break;
			case 16:		/* Y intercept */
				if (gg_mute_linear_regression( getsetlength( cg, cs ), 
						getx( cg, cs ), gety( cg, cs ), &dummy, &datum ) ) {
					sprintf( tbuf+strlen(tbuf), 
						"Unable to find y-intercept of set %d\n", cs );
					gg_errwin( tbuf );
					extract_err = 1;
				}
				break;
			case 17:		/* set length  */
				datum = getsetlength( cg, cs );
				break;
			case 18:		/* half maximal widths */
                                if (gg_get_half_max_width(getsetlength(cg, cs), getx(cg,cs), 
					   gety(cg,cs), &datum) != RETURN_SUCCESS) {
					sprintf( tbuf+strlen(tbuf), 
						"Unable to find half maximal width of set %d\n", cs );
					extract_err = 1;
					gg_errwin( tbuf );
				}
				break;
			case 19:		/* Barycenter X */
				gg_get_barycenter( getsetlength( cg, cs ), gety(cg,cs), 
									getx(cg,cs), &datum );
				break;
			case 20:		/* Barycenter Y */
				gg_get_barycenter( getsetlength( cg, cs ), getx(cg,cs), 
									gety(cg,cs), &datum );
				break;
			case 21:		/* X of Maximum Y */
				getsetminmax(cg, cs, &xmin, &xmax, &ymin, &ymax);
                                gg_get_max_pos( gety(cg, cs), getx( cg, cs ),
							getsetlength( cg, cs ), ymax, &datum ); 
				break;
			case 22:		/* Y of Maximum X */
				getsetminmax(cg, cs, &xmin, &xmax, &ymin, &ymax);
                                gg_get_max_pos( getx(cg, cs), gety( cg, cs ),
							getsetlength( cg, cs ), xmax, &datum ); 
				break;
			case 23:		/* cumulative sum */
				datum = do_int(cg, cs, 1);
				break;
		}
		if( !extract_err )
			add_point(gto, ns, abscissa[i], datum);
		cs++;
	}

	/* set comment */	
	switch( feature ) {
		case 0:			/* Y minimum */
			sprintf(tbuf,"Y minima of graph %d",cg); 
			break;
		case 1: 		/* Y maximum */
			sprintf(tbuf,"Y maxima of graph %d",cg);
			break;
		case 2: 		/* Y mean    */
			sprintf(tbuf,"Y means of graph %d",cg);
			break;
		case 3:			/* Y std dev */
			sprintf(tbuf,"Y std. dev.'s of graph %d",cg);
			break;
		case 4:			/* Y median  */
			sprintf(tbuf,"Y medians of graph %d",cg);
			break;
		case 5:			/* X minimum */
			sprintf(tbuf,"X minima of graph %d",cg); 
			break;
		case 6: 		/* X maximum */
			sprintf(tbuf,"X maxima of graph %d",cg);
			break;
		case 7: 		/* X mean    */
			sprintf(tbuf,"X means of graph %d",cg);
			break;
		case 8:			/* X std dev */
			sprintf(tbuf,"X std. dev.'s of graph %d",cg);
			break;
		case 9:			/* X median  */
			sprintf(tbuf,"X medians of graph %d",cg);
			break;
		case 10: 		/* frequency and period */
			sprintf(tbuf,"frequencies of graph %d",cg);
			break;
		case 11:
			sprintf(tbuf,"periods of graph %d",cg);
			break;
		case 12:		/* first zero crossing */
			sprintf(tbuf,"zero crossings of graph %d",cg);
			break;
		case 13:		/* rise time */
			sprintf(tbuf,"rise times of graph %d",cg);
			break;
		case 14: 		/* fall time */
			sprintf(tbuf,"fall times of graph %d",cg);
			break;
		case 15: 		/* slopes     */
			sprintf(tbuf,"slopes of graph %d",cg);
			break;
		case 16: 		/* Y intercepts */
			sprintf(tbuf,"Y intercepts of graph %d",cg);
			break;
		case 17: 		/* set lengths */
			sprintf(tbuf,"set lengths of graph %d",cg);
			break;
		case 18: 		/* 1/2 maximal widths */
			sprintf(tbuf,"half maximal widths of graph %d",cg);
			break;
		case 19: 		/* barycenter X */
			sprintf(tbuf,"X barycenters of graph %d",cg);
			break;
		case 20: 		/* barycenter Y */
			sprintf(tbuf,"Y barycenters of graph %d",cg);
			break;
		case 21:		/* X of maximum Y */
			sprintf(tbuf,"X positions of maximum Y's of graph %d",cg);
			break;
		case 22:		/* Y of maximum X */
			sprintf(tbuf,"Y positions of maximum X's of graph %d",cg);
			break;
		case 23:		/* integral */
			sprintf(tbuf,"integrals of sets of graph %d",cg);
			break;
	}
	set_set_hidden(gto, ns, FALSE);
        setcomment( gto, ns, tbuf );
	gg_view_comment (TRUE);
	xfree( abscissa );
}



/**
 *  Replace   (identical to) mute_linear_regression   featext.c 498
* linear regression without posting results to log
*/
int gg_mute_linear_regression(int n, double *x, double *y, double *slope, 
			double *intercept)
{
    double xbar, ybar;		/* sample means */
    double sdx, sdy;		/* sample standard deviations */
    double SXX, SYY, SXY;	/* sums of squares */
    int i;

    if (n <= 3) {
		return 1;
    }
    xbar = ybar = 0.0;
    SXX = SYY = SXY = 0.0;
    for (i = 0; i < n; i++) {
		xbar = xbar + x[i];
		ybar = ybar + y[i];
    }
    xbar = xbar / n;
    ybar = ybar / n;
    for (i = 0; i < n; i++) {
		SXX = SXX + (x[i] - xbar) * (x[i] - xbar);
		SYY = SYY + (y[i] - ybar) * (y[i] - ybar);
		SXY = SXY + (x[i] - xbar) * (y[i] - ybar);
    }
    sdx = sqrt(SXX / (n - 1));
    sdy = sqrt(SYY / (n - 1));
    if (sdx == 0.0) {
		return 2;
    }
    if (sdy == 0.0) {
		return 2;
    }
    *slope = SXY / SXX;
    *intercept = ybar - *slope * xbar;
    return 0;
}
/**
 *  Replace  (identical to) get_rise_time            featext.c 539
 * assume graph starts off at ymin and rises to ymax 
 * Determine time to go from 10% to 90% of rise
 */
int gg_get_rise_time( int setl, double *xv, double *yv, double min, double max,
		double *width )
{
	int x10=0, x90;
	double amp10, amp90;
	
	amp10 = min + (max-min)*0.1;
	amp90 = min + (max-min)*0.9;
	while( x10<setl && yv[x10]<amp10  )
		x10++;
	
	if( x10==setl || x10==0)
		return 1;
	
	x90 = x10+1;
	
	while( x90<setl && yv[x90]<amp90 )
		x90++;
	
	*width = gg_linear_interp( yv[x90-1], xv[x90-1], yv[x90], xv[x90], amp90 ) -
			 gg_linear_interp( yv[x10-1], xv[x10-1], yv[x10], xv[x10], amp10 );
	return 0;
}

/**
 *  Replace  (identical to) get_fall_time            featext.c 565
 *  assume graph starts off at ymax and drops to ymin 
 *  Determine time to go from 90% to 10% of drop	
 */
int gg_get_fall_time( int setl, double *xv, double *yv, double min, double max,
		double *width )
{
	int x10, x90=0;
	double amp10, amp90;
	
	amp10 = min + (max-min)*0.1;
	amp90 = min + (max-min)*0.9;
	while( x90<setl && yv[x90]>amp90 )
		x90++;
	
	if( x90==setl || x90==0)
		return 1;
	
	x10= x90+1;
	
	while( x10<setl && yv[x10]>amp10  )
		x10++;

	if( x10==setl )
		return 1;
	
	
	*width = gg_linear_interp( yv[x10-1], xv[x10-1], yv[x10], xv[x10], amp10 )-
	         gg_linear_interp( yv[x90-1], xv[x90-1], yv[x90], xv[x90], amp90 );
	return 0;
}


/**
 *  Replace  (identical to)  get_zero_crossing        featext.c 594
 */
int gg_get_zero_crossing( int setl, double *xv, double *yv, double *crossing )
{
	int i=0;
	
	while( i<setl && yv[i] != 0. && yv[i]*yv[i+1]>0. )
		i++;
	
	if( i==setl )
		return 1;
	
	if( yv[i] == 0 )
		*crossing = xv[i];
	else
		*crossing = gg_linear_interp( yv[i], xv[i], yv[i+1], xv[i+1], 0 );

	return 0;
}

/**
 *  Replace  (identical to) get_half_max_width       featext.c 614
 * Get FWHM of the highest peak
 */
int gg_get_half_max_width(int len, double *x, double *y, double *width)
{
    int i, imin, imax;
    double ymin, ymax, yavg;
    double x_u, x_d;

    minmax(y, len, &ymin, &ymax, &imin, &imax);
    yavg = (ymin + ymax)/2.0;
	
    i = imax;
    while (i >= 0 && y[i] > yavg) {
        i--;
    }
    if (i < 0) {
        return RETURN_FAILURE;
    }
    x_d = gg_linear_interp(y[i], x[i], y[i + 1], x[i + 1], yavg);

    i = imax;
    while (i < len && y[i] > yavg) {
        i++;
    }
    if (i == len) {
        return RETURN_FAILURE;
    }
    x_u = gg_linear_interp(y[i - 1], x[i - 1], y[i], x[i], yavg);

    *width = fabs(x_u - x_d);
    
    return RETURN_SUCCESS;
}

/**
 *  Replace (identical to)   linear_interp            featext.c 647
 * linear interpolate between two points, return a y value given an x 
 */
double gg_linear_interp( double x1, double y1, double x2, double y2, double x )
{
	return y1 + ( x-x1 )*(y2-y1)/(x2-x1);
}




/**
 *  Replace  getmedian                featext.c 653
 */
/* get the median of the X or Y portion of a set */
int gg_getmedian (int grno, int setno, int sorton, double *median)
{
	int setlen;
	double *setdata;
	
	setlen = getsetlength( get_cg(), setno );
	setdata = (double *)xmalloc( setlen*sizeof(double) );
	if( sorton == DATA_X )
		memcpy( setdata, getx( grno, setno ), setlen*sizeof(double) );
	else
		memcpy( setdata, gety( grno, setno ), setlen*sizeof(double) );
	
	qsort( setdata, setlen, sizeof(double), gg_dbl_comp );
	
	if( setlen%2 )		/* odd set */
		*median = setdata[(setlen+1)/2-1];
	else
		*median = ( setdata[setlen/2-1] + setdata[setlen/2] )/2.;

	xfree( setdata );
	return 0;
}



/**
 *  Replace  (identical to)  dbl_comp                 featext.c 676
 */
int gg_dbl_comp( const void *a, const void *b )
{
	if( *(double *)a < *(double *)b )
		return -1;
	else if( *(double *)a > *(double *)b )
		return 1;
	else
		return 0;
}


/**
 *  Replace  get_barycenter           featext.c 686
 */
int gg_get_barycenter( int n, double *x, double *y, double *barycenter )
{
	double wsum=0;
	
	*barycenter = 0;
	for( n--; n>=0; n-- ) {
		wsum += x[n]*y[n];
		*barycenter += x[n];
	}
	*barycenter = wsum/(*barycenter);
	return 0;
}

/**
 * Replace  (identical to)  get_max_pos              featext.c 700
 *  given maximum of set a, find the corresponding entry in set b 
 */
void gg_get_max_pos( double *a, double *b, int n, double max, double *d )
{
	int i=-1;
	
	while( ++i<n && a[i] != max  );
	
	if( i==n )
		return;
	else
		*d = b[i];
}
