/*
 * Copyright 1988 Anant Agarwal
 *
 * Permission to use, copy, modify, and distribute this software and its
 * documentation for any purpose and without fee is hereby granted, provided
 * that the above copyright notice appear in all copies and that both that
 * copyright notice and this permission notice appear in supporting
 * documentation.  No representations about the suitability of this software
 * for any purpose.  It is provided "as is" without express or implied 
 * warranty.
 *
 * Author:  Anant Agarwal, MIT Laboratory for Computer Science
 */

#include <strings.h>	
#include <stdio.h>
#include <math.h>
#include "plotio.h"
#include "defs.h"
#include "structs.h"
#include "plot.h"



extern FILE *fout;
extern double getdegrees(), atan(), charheight(), charwidth();
extern double log10(), mylog(), exp(), log();



/*
 * Do least squares linear fit. Weights w[i] not implemented. 
 * If plotting on a log or log-log graph, the data is reduced prior to 
 * fitting and then exponetiated again prior to plotting. The idea is 
 * produce a good stright line fit through the plotted points as displayed.
 */

regression(c)
int c;/* current curve */
{
pointType *ptPtr; /*current point */
int numpts;/* current point */
double maxx, minx;
double sumx, sumy, avgx, avgy;
double num, den, xdiff;
double slope, intercept;
Boolean logx, logy;

	if (axis[xpt].scale.type == logstype) logx = TRUE;
	else logx = FALSE;
	if (axis[ypt].scale.type == logstype) logy = TRUE;
	else logy = FALSE;

	/* print the points and find average and domain */
	numpts = (curves[c])->lastpoint;
	maxx = -BIG; minx = BIG;
	sumx = 0.0; sumy = 0.0;
	ptPtr = (curves[c])->fstPtr;
	while (ptPtr != nil)	
	    {
		/* symbol part */
		InitObjSymbol(obj[curobj], c);
		(obj[curobj])->from[xpt] = ptPtr->pt[xpt];
		(obj[curobj])->from[ypt] = ptPtr->pt[ypt];
		GetNewObjcurobj();
		sumy = sumy + (logy ? log(ptPtr->pt[ypt]) : ptPtr->pt[ypt]);
		sumx = sumx + (logx ? log(ptPtr->pt[xpt]) : ptPtr->pt[xpt]);
		if( ptPtr->pt[xpt] > maxx ) maxx=ptPtr->pt[xpt];
		if( ptPtr->pt[xpt] < minx ) minx=ptPtr->pt[xpt];
		ptPtr = ptPtr->nxtPtr;
	    }
	if( numpts != -1 )
	    {
		avgx = sumx / (numpts + 1);
		avgy = sumy / (numpts + 1);
	    }


	/* We use the X - Xavg form since it is more numerically stable */
	/* It requires two passes over the data though */
	num = 0; /* numerator = sum(y[i](x[i]-xavg)w[i]) */
	den = 0; /* demonimator = sum(w[i](x[i]-xavg)**2) */
	ptPtr = (curves[c])->fstPtr;
	while (ptPtr != nil)	
	    {
		xdiff =	(logx ? log(ptPtr->pt[xpt]) : ptPtr->pt[xpt] ) - avgx;
		num = num + ( logy ?  log(ptPtr->pt[ypt]) : ptPtr->pt[ypt]  )  * xdiff;
		den = den + xdiff * xdiff;
		ptPtr = ptPtr->nxtPtr;
	    }
	if( (den != 0.0) && (numpts != 0) ) /* not(1 pt or vertical line) */
	    {
		/* non-vertical line */
		slope = num / den;
		intercept = avgy - slope * avgx;
		if( logy )
		    {
			LinSt( minx, exp(intercept+slope *
					( logx ? log(minx) : minx) ),
					 (curves[c])->gray ,
					(curves[c])->interp,
			      		(curves[c])->type,
			      		c);
			(obj[curobj])->from[xpt] = minx;
			(obj[curobj])->from[ypt] = exp( intercept + slope *
						( logx ? log(minx) : minx) );
			(obj[curobj])->un.line.to[xpt] = maxx;
			(obj[curobj])->un.line.to[ypt] = exp(intercept+slope *
						( logx ? log(maxx) : maxx) );
		    }
		else
		    {
			LinSt( minx, intercept+slope *
						( logx ? log(minx) : minx),
						(curves[c])->gray,	
						(curves[c])->interp,
			      			(curves[c])->type,
			      			c);
			(obj[curobj])->from[xpt] = minx;
			(obj[curobj])->from[ypt] = intercept + slope *
						( logx ? log(minx) : minx);
			(obj[curobj])->un.line.to[xpt] = maxx;
			(obj[curobj])->un.line.to[ypt] = intercept + slope *
						( logx ? log(maxx) : maxx);
		    }
		InitCurLineObj(obj[curobj], c);
		GetNewObjcurobj();
		LinEnd((curves[c])->gray,(curves[c])->interp);

		/* output least squares fit to postscript file */
		if(( !XView ) && (!IView) )
		    {
			if( logy )
				printf("%% Regression: log(y) = ");
			else
				printf("%% Regression: y = ");
			if( logx )
				printf("%10g * log(x) + %10g\n",slope,intercept);
			else
				printf("%10g * x + %10g\n",slope,intercept);
		    }
	    }
	 else
	    {
		if( numpts != 0 )
		    {
			fprintf(stderr,
				"Vertical line from least squares fit.\n");
			printf("%% Regression: Vertical line. Graph omitted.\n");
		    }
		else
		    {
			printf("%% Regression: One data point. Graph omitted.\n");
		    }
	    }
    }