/* ======================================================= *
 * Copyright 1998-2005 Stephen C. Grubb                    *
 * http://ploticus.sourceforge.net                         *
 * Covered by GPL; see the file ./Copyright for details.   *
 * ======================================================= */

/* PROC CURVEFIT - Fit a curve to the data.
   Result becomes the new data that plotting procedures access.  */

#include "pl.h"

#define MAXPTS 1000	/* default max # input points for movingavg curve */
#define MAXORDER 21	/* max bspline order value (no limit for movingavg order) */
#define MAXBSP_IN 100  /* max # input points for bspline curve */

#define MOVINGAVG	'm'
#define REGRESSION	'r'
#define BSPLINE		'b'
#define SIMPLEAVG	'a'
#define INTERPOLATED	'i'

/* static double in[MAXPTS][2]; */

/* the PLV vector is used for curve points */
static int bspline(), mavg(), plainavg(), lregress();

static int dblcompare(const void *a, const void *b);  
/* static int dblcompare(double *f, double *g);   */

int
PLP_curvefit()
{
int i;
char attr[NAMEMAXLEN], val[256];
char *line, *lineval;
int nt, lvp;
int first;

int stat;

int order;
int xfield;
int yfield;
int npts; /* number of input points */
int nresultpoints; /* number of result points */
int irow;
double resolution;
int showresults;
char linedetails[256];
double linestart, linestop;
double calcstart, calcstop;
int calcrangegiven;
char curvetype[40];
char legendlabel[256]; /* raised (can contain urls for clickmap) scg 4/22/04 */
int linerangegiven;
int statsonly;
char selectex[256];
int selectresult;
char numstr[100];
int xsort;
int maxinpoints;
double *inpoints, *inp;
double drawx, drawy, prevdrawx, prevdrawy;
int doclip;

TDH_errprog( "pl proc curvefit" );

/* initialize */
order = 4;
xfield = -1;
yfield = -1;
nresultpoints = -1;
resolution = 5.0;
showresults = 0;
strcpy( linedetails, "" );
linestart = EDXlo;
linestop = EDXhi;
calcrangegiven = 0;
linerangegiven = 0;
strcpy( curvetype, "movingavg" );
statsonly = 0;
strcpy( selectex, "" );
strcpy( legendlabel, "" ); /* added scg 7/28/03 */
xsort = 0;
maxinpoints = MAXPTS;
doclip = 0;


/* get attributes.. */
first = 1;
while( 1 ) {
	line = getnextattr( first, attr, val, &lvp, &nt );
	if( line == NULL ) break;
	first = 0;
	lineval = &line[lvp];

	if( stricmp( attr, "xfield" )==0 ) xfield = fref( val ) - 1;
	else if( stricmp( attr, "yfield" )==0 ) yfield = fref( val ) - 1;
	else if( stricmp( attr, "order" )==0 ) order = atoi( val );
	else if( stricmp( attr, "resolution" )==0 ) resolution = atof( val );
	else if( stricmp( attr, "xsort" )==0 ) {
		if( strnicmp( val, YESANS, 1 )==0 ) xsort = 1;
		}
	else if( stricmp( attr, "linedetails" )==0 ) strcpy( linedetails, lineval );
	else if( stricmp( attr, "legendlabel" )==0 ) strcpy( legendlabel, lineval );
	else if( stricmp( attr, "select" )==0 ) strcpy( selectex, lineval );
	else if( stricmp( attr, "linerange" )==0 ) {
		if( lineval[0] != '\0' ) linerangegiven = 1;
		getrange( lineval, &linestart, &linestop, 'x', EDXlo, EDXhi );
		}
	else if( stricmp( attr, "calcrange" )==0 ) {
		if( lineval[0] != '\0' ) {
			calcrangegiven = 1;
			getrange( lineval, &calcstart, &calcstop, 'x', EDXlo, EDXhi );
			}
		else calcrangegiven = 0;
		}
	else if( stricmp( attr, "curvetype" )==0 ) strcpy( curvetype, val );
	else if( stricmp( attr, "maxinpoints" )==0 ) maxinpoints = atoi( val );
	else if( stricmp( attr, "showresults" )==0 ) {
		if( strnicmp( val, YESANS, 1 )==0 ) showresults = 1;
		else showresults = 0;
		}
	else if( stricmp( attr, "clip" )==0 ) {
		if( strnicmp( val, YESANS, 1 )==0 ) doclip = 1;
		else doclip = 0;
		}
	else if( stricmp( attr, "statsonly" )==0 ) {
		if( strnicmp( val, YESANS, 1 )==0 ) statsonly = 1;
		else statsonly = 0;
		}
	else Eerr( 1, "curvefit attribute not recognized", attr );
	}


/* overrides and degenerate cases */
if( Nrecords < 1 ) 
 	return( Eerr( 17, "No data has been read yet w/ proc getdata", "" ) );

if( yfield < 0 || yfield >= Nfields ) return( Eerr( 601, "yfield not specified or out of range", "" ) );
if( xfield < 0 || xfield >= Nfields ) return( Eerr( 601, "xfield not specified or out of range", "" ) );

if( !scalebeenset() )
         return( Eerr( 51, "No scaled plotting area has been defined yet w/ proc areadef", "" ) );
 
if( strnicmp( legendlabel, "#usefname", 9 )==0 ) getfname( yfield+1, legendlabel );


/* now do the computation work.. */
/* -------------------------- */

/* allocate memory for the input points list.. */
inpoints = (double *) malloc( maxinpoints * sizeof( double ) * 2 );
inp = inpoints;

/* process input data.. */
npts = 0;
for( irow = 0; irow < Nrecords; irow++ ) {

	if( selectex[0] != '\0' ) { /* process against selection condition if any.. */
                stat = do_select( selectex, irow, &selectresult );
                if( stat != 0 ) { Eerr( stat, "Select error", selectex ); continue; }
                if( selectresult == 0 ) continue; /* reject */
                }

	/* in[npts][1] = fda( irow, yfield, Y ); */
	inp += 1;  /* because we're getting Y */
	*inp = fda( irow, yfield, Y );
	if( Econv_error() ) { conv_msg( irow, yfield, "yfield" ); inp -= 1; continue; }   /* bug - inp-=1 added scg 2/3/05 */

	inp -= 1;  /* now back up one to get X */
	if( xfield < 0 ) *inp = (double)irow;   /* in[npts][0] = (int)irow; */
	else 	{
		/* in[npts][0] = fda( irow, xfield, X ); */
		*inp = fda( irow, xfield, X );
		if( Econv_error() ) { conv_msg( irow, xfield, "xfield" ); continue; }
		}

	/* compute curve only for points within calcrange */
	if( calcrangegiven ) {
		/* if( in[npts][0] < calcstart ) continue; */
		/* else if( in[npts][0] > calcstop ) break; */
		if( *inp < calcstart ) continue;
		else if( *inp > calcstop ) break;
		}

	if( curvetype[0] == BSPLINE && npts >= MAXBSP_IN ) {
		Eerr( 2599, "max of 100 input points allowed for bspline exceeded; curve truncated", "" );
		break;
		}
	if( npts >= maxinpoints ) {
		Eerr( 2599, "maxinpoints exceeded, curve truncated (see maxinpoints attribute)", "" );
		break;
		}
	npts++;
	inp+=2; /* to next point slot */
	}


/* sort points on x - added scg 11/2000 */
if( curvetype[0] != INTERPOLATED || xsort ) 
	qsort( inpoints, npts, sizeof(double)*2, dblcompare );


if( curvetype[0] == MOVINGAVG ) { 
	mavg( inpoints, npts, PLV, order );
	nresultpoints = npts;
	}


else if( curvetype[0] == REGRESSION ) { 
	double rlinestart, rlinestop;

	if( linerangegiven ) {
		rlinestart = linestart;
		rlinestop = linestop;
		}
	else	{
		/* rlinestart = in[0][0]; rlinestop = in[npts-1][0]; */
		inp = inpoints; rlinestart = *inp; inp += (npts-1)*2; rlinestop = *inp;
		/* fprintf( stderr, "[rlinestart=%g   rlinestop=%g ]\n", rlinestart, rlinestop ); */
		}
	lregress( inpoints, npts, PLV, rlinestart, rlinestop );
	nresultpoints = 2;

	/* vertical line (degenerate case) suppress.. */
	if( dat2d( 0, 0 ) == dat2d( 1, 0 ) ) nresultpoints = 0;  
	
	/* clip to plotting area.. */
	stat = Elineclip( &dat2d(0,0), &dat2d(0,1), &dat2d(1,0), &dat2d(1,1), EDXlo, EDYlo, EDXhi, EDYhi );
	if( stat != 0 ) nresultpoints = 0;
	}

else if( curvetype[0] == BSPLINE ) {  
	nresultpoints = npts * resolution;
	if( nresultpoints >= PLVhalfsize ) nresultpoints = PLVhalfsize-1;

	if( order >= MAXORDER ) {
		Eerr( 2158, "Using max bspline order of 20", "" );
		order = 20;
		}
	if( npts < order ) {
		if( inpoints != NULL ) free( inpoints );
		return( Eerr( 4892, "Must have at least 'order' data points", "" ) );
		}

	/* do the computation.. */
	bspline( inpoints, npts, PLV, nresultpoints, order );
	}

else if( curvetype[0] == INTERPOLATED ) { 
	stat = PL_smoothfit( inpoints, npts, PLV, &nresultpoints );
	}

else 	{ 	/* average curve (basic) */
	plainavg( inpoints, npts, PLV, order );
	nresultpoints = npts;
	}

if( inpoints != NULL ) free( inpoints );


/* curve has been generated.. now draw the line.. */
/* ---------------------------------------------- */
linedet( "linedetails", linedetails, 1.0 );
Emovu( dat2d(0,0), dat2d(0,1) );
if( showresults ) fprintf( PLS.diagfp, "// generated curve points follow:\n%g %g\n", dat2d(0,0), dat2d(0,1) );


for( i = 1; i < nresultpoints; i++ ) {

	drawx = dat2d(i,0);
	drawy = dat2d(i,1);

	/* draw only within linerange.. */
	if( i > 0 && drawx > linestart && (dat2d(i-1,0)) < linestart ) Emovu( drawx, drawy );

	else if( drawx < linestart ) continue;
	else if( drawx > linestop ) break;

	if( doclip && !statsonly && i > 0 ) {
		prevdrawx = dat2d(i-1,0);
		prevdrawy = dat2d(i-1,1);
		stat = Elineclip( &prevdrawx, &prevdrawy, &drawx, &drawy, EDXlo, EDYlo, EDXhi, EDYhi );
		if( stat ) goto BOTTOM;
		Emovu( prevdrawx, prevdrawy );
		}
	if( !statsonly ) Elinu( drawx, drawy );  
	BOTTOM:
	if( showresults ) fprintf( PLS.diagfp, "%g %g\n", drawx, drawy );
	}


/* set YFINAL and Xfinal */
i--;
Euprint( numstr, X, dat2d(i,0), "" );
setcharvar( "XFINAL", numstr );
Euprint( numstr, Y, dat2d(i,1), "" );
setcharvar( "YFINAL", numstr );




if( legendlabel[0] != '\0' ) PL_add_legent( LEGEND_LINE, legendlabel, "", linedetails, "", "" );

return( 0 );
}



/* =========================== */
static int mavg( in, npts, out, order )
double in[][2];  /* input points */
int npts;	 /* number of input points */
double out[][2]; /* output points (same n as input points array) */
int order;	 /* # of points to average */
{
int i, j;
int avgstart;
double accum;

for( i = 0; i < npts; i++ ) {
	avgstart = i - (order - 1);
	if( avgstart < 0 ) avgstart = 0;
	accum = 0.0;
	for( j = avgstart; j <= i; j++ ) accum += in[j][1];	
	out[i][0] = in[i][0];
	out[i][1] = accum / (double)(( i - avgstart ) + 1);
	}
return( 0 );
}

/* =========================== */
/* same as movingavg but averages in points to right as well.. */
static int plainavg( in, npts, out, order )
double in[][2];  /* input points */
int npts;	 /* number of input points */
double out[][2]; /* output points (same n as input points array) */
int order;	 /* # of points to average */
{
int i, j;
int avgstart, avgstop;
double accum;


for( i = 0; i < npts; i++ ) {
	avgstart = i - (order - 1);
	avgstop = i + (order - 1);
	if( avgstart < 0 ) avgstart = 0;
	if( avgstop >= npts ) avgstop = npts-1;
	accum = 0.0;
	for( j = avgstart; j <= avgstop; j++ ) accum += in[j][1];	
	out[i][0] = in[i][0];
	out[i][1] = accum / (double)(( avgstop - avgstart ) + 1  );
	}
return( 0 );
}

/* ============================== */
/* LREGRESS - linear regression curve */
static int lregress( in, npts, out, start, stop )
double in[][2];
int npts;	 /* number of input points */
double out[][2]; /* output points (n=2 since it's a straight line) */
double start, stop; /* X values - result line start and stop */
{
int i;
double sumx, sumxx, sumy, sumyy, sumxy;
double b, m;
double numer, denom, denomleft, denomright;
char buf[128], tok[128];
double sqrt();
char *GL_autoroundf();

sumx = sumxx = sumy = sumyy = sumxy = 0.0;
for( i = 0; i < npts; i++ ) {
	sumx += in[i][0];
	sumxx += (in[i][0] * in[i][0]);
	sumy += in[i][1];
	sumyy += (in[i][1] * in[i][1]);
	sumxy += (in[i][0] * in[i][1]);
	}
/* compute x intercept (b) */
numer = (sumy * sumxx) - (sumx * sumxy);
denom = ( (double)npts * sumxx ) - (sumx * sumx);
b = numer / denom; 

/* compute slope (m) */
numer = ((double)npts * sumxy) - (sumx * sumy);
denom = ((double)npts * sumxx) - (sumx * sumx);
if( denom == 0.0 ) m = 0.0;  /* ? */
else m = numer / denom;

out[0][0] = start;
out[0][1] = (m * start) + b;
out[1][0] = stop;
out[1][1] = (m * stop) + b;

strcpy( tok, GL_autoroundf( m, 0 ) );
sprintf( buf, "Y = %s + %sX", GL_autoroundf(b,0), tok );
TDH_setvar( "REGRESSION_LINE", buf );

/* compute r (pearson correlation coeff) */
/* numer = ((double) nvalues * sumxy) - (sumx * sumy ); */
denomleft = ((double) npts * sumxx)  - (sumx * sumx);
denomright = ((double) npts * sumyy)  - (sumy * sumy);
denom = sqrt( denomleft * denomright );

if( denom == 0.0 ) strcpy( buf, "(none)" );
else sprintf( buf, "%s", GL_autoroundf( (numer/denom), 0 ) );
TDH_setvar( "CORRELATION", buf );

return( 0 );
}

/* ======================== */
/* BSPLINE curve */
static int
bspline( in, npts, out, ncv, order )
double in[][2];  /* input points */
int npts;	 /* number of input points */
double out[][2]; /* output points */
int ncv;         /* number of output points to generate */
int order;	 /* order of the curve  (2..n_in) */
{
int i, j, k, n;
int nknot; /* size of knot vector */
double t; /* parameter */
double N[MAXBSP_IN][MAXORDER]; /* weighting function */ 
double knot[MAXBSP_IN+MAXORDER]; /* knot vector */
double tmp1, tmp2, tmp3, tmp4, tmp5, tmp6;
int n_out;



/* generate the knot vector */
/* ------------------------ */
for( i = 0; i < MAXBSP_IN; i++ ) knot[i] = 0.0; /* init */
n = npts - 1;
nknot = n + 1 + order;
for( i = 0; i < order; i++ ) knot[i] = 0;
j = 1;
for( ; i < nknot-order; i++ ) knot[i] = j++;
for( ; i < nknot; i++ ) knot[i] = j;
/* printf( "Knot= [ " ); 
 *  for( i = 0; i < nknot; i++ ) printf( "%g ", knot[i] ); 
 * printf( "]\n" ); 
 */


t = 0.0;
n_out = 0;
for( j = 0; j < ncv; j++ ) { /* for each point to be generated.. */

    /* do the N(?,1) set.. */
    for( i = 0; i <= n+1; i++ ) {
	if( knot[i] <= t && t < knot[i+1] ) N[i][1] = 1.0;
	else N[i][1] = 0.0;
	}

    /* do middle N.. */
    for( k = 2; k <= order; k++ ) {
	for( i = 0; i <= npts; i++ ) {
	    tmp1 = ((t - knot[i])*N[i][k-1]);
	    tmp2 = (knot[i+k-1] - knot[i]);
	    tmp3 = ((knot[i+k] - t)*N[i+1][k-1]);
	    tmp4 = (knot[i+k] - knot[i+1]);
	    if( tmp2 == 0.0 ) { tmp5 = 0.0; }
	    else { tmp5 = tmp1 / tmp2; }
	    if( tmp4 == 0.0 ) { tmp6 = 0.0; }
	    else { tmp6 = tmp3 / tmp4; }
	    N[i][k] = tmp5 + tmp6;
	    }
	}

    if( j == ncv-1 ) N[n][order] = 1.0; /* for last point */

    tmp1 = tmp2 = tmp3 = 0;
    for( i = 0; i < npts; i++ ) {
	tmp1 += (in[i][0]*N[i][order]);
	tmp2 += (in[i][1]*N[i][order]);
	tmp3 += 0.0; /* (in[i][2]*N[i][order]); */
	}

    /* put curve into D2 */
    out[n_out][0] = tmp1;
    out[n_out][1] = tmp2;
    /* out[n_out][2] = tmp3; */
    n_out++;
    t += ( knot[nknot-1] / (double)(ncv - 1) );
    }

return( 1 );
}

/* ============================= */
static int 
dblcompare( a, b )
const void *a, *b;

/* dblcompare( f, g )
 * double *f, *g;
 */  /* changed to eliminate gcc warnings  scg 5/18/06 */

{
double *f, *g;
f = (double *)a;
g = (double *)b;
if( *f > *g ) return( 1 );
if( *f < *g ) return( -1 );
return( 0 );
}


/* ======================================================= *
 * Copyright 1998-2005 Stephen C. Grubb                    *
 * http://ploticus.sourceforge.net                         *
 * Covered by GPL; see the file ./Copyright for details.   *
 * ======================================================= */