/* QRegressionApp.java (January 2004)
 * This program was written by J. Hasbun on December 2003 at the
 * State Univ. of West Georgia in order to demonstrate a way to do
 * linear regresssion. It is possible to linearize
 * certain non-linear fits, and we include that here. The program is written
 * as part of our endeavors to expand the OSP applications.
 * J. Hasbun would like to acknowledge the student Maxwell Perkins who helped
 * to write this code.
 * This version was further enhanced by Max Perkins with Dr. J. Hasbun
 * to include quadratic fitting, table editing, and copy, paste, and delete
 * functions. Max Perkins was sponsored by the NASA Space
 * Grant Consortium through Dr. Ben de Mayo of the Physics Department
 * at the State University of West Georgia. Spring 2004.
 */

package org.opensourcephysics.jahasbun.QE;

import org.opensourcephysics.display.*;
import org.opensourcephysics.jahasbun.QE.display.*;
import org.opensourcephysics.controls.*;
import org.opensourcephysics.davidson.applets.*;
import java.text.NumberFormat;
import java.awt.Color;
import javax.swing.Timer;
import org.opensourcephysics.jahasbun.QE.applets.*;
import org.opensourcephysics.jahasbun.QE.controls.*;

public class QERegressionApp implements Calculation {
   private Control myControl;
   private EAppletOSPControl myAppletControl;
//-- plotting declaration
   PlottingPanel plotPanel    = new PlottingPanel ("X axis ", "Y axis", "Line of Regression");
   DrawingFrame  drawingFrame = new DrawingFrame ("Plot", plotPanel);
//-- tables declaration
   EDataTable dataTableb = new EDataTable();
   EDataTable dataTablea = new EDataTable();
   EDataTableFrame frameb= new EDataTableFrame("Line Fit Calculation", dataTableb);
   EDataTableFrame framea= new EDataTableFrame("Input Data", dataTablea);
//-- data arrays declaration
   EDatasetManager datasetsa = new EDatasetManager();
   EDatasetManager datasetsb = new EDatasetManager();
   EDataset dataseta = new EDataset ();
   EDataset datasetb = new EDataset ();
//-- other variables
   String message;
   int n, Npts, M, InitialRows;
   String Model;
   NumberFormat nf = NumberFormat.getInstance();// used in doBestFitt()
   double x, y, m, b, r, a, c, syx;
   double eps=1.e-5;
   double [] xx,yy,y2;

   public QERegressionApp(){
//-- plotting properties
     drawingFrame.setSize(310, 310);
     drawingFrame.setLocation(1, 1);
     plotPanel.enableInspector(true);
//-- dataseta,b properties
     dataseta.setSorted(true);//sorts the data
     dataseta.setConnected(false);
     dataseta.setMarkerShape(EDataset.CIRCLE);
     dataseta.setMarkerColor(Color.red, Color.red, Color.red);
     datasetb.setConnected(true);
     datasetb.setMarkerShape(EDataset.NO_MARKER);
//-- datasetsb and it table properties
     datasetsb.setXYColumnNames(0, "Xc", "Yc");
     dataTableb.setRowNumberVisible(true);
     dataTableb.setRowSelectionAllowed(true);
     dataTableb.setMaximumFractionDigits(5);
     frameb.setSize(200, 250); //200, 293
     frameb.setLocation(1, 320);
//-- datasetsa and it table properties
     datasetsa.setSorted(true);
     datasetsa.setXYColumnNames(0, "Xd", "Yd");
     dataTablea.setRowNumberVisible(true);
     dataTablea.setRowSelectionAllowed(true);
     dataTablea.setMaximumFractionDigits(5);
     dataTablea.add(datasetsa);
     dataTablea.refreshTable();
     framea.setSize(200, 250);
     framea.setLocation(200, 320);
// -- doBestFitt() information
     nf.setMaximumFractionDigits(2);
     frameb.show();
     drawingFrame.show();
     framea.show();


     if(!org.opensourcephysics.display.OSPFrame.appletMode) {
        framea.setDefaultCloseOperation(EDataTableFrame.DO_NOTHING_ON_CLOSE);
        frameb.setDefaultCloseOperation(EDataTableFrame.DO_NOTHING_ON_CLOSE);
        drawingFrame.setDefaultCloseOperation(DrawingFrame.DO_NOTHING_ON_CLOSE);
     }
   }

     public void doBestFitt2() {
/*
 * Written by J. Hasbun and Max Perkins for OSP purposes, Fall 2003
 * Regression line calculating program: y=m*x+b type
 * some non-linear cases are linearized, for example, the power fit
 * y=c*x^a is equivalent to ln(y)=ln(c)+a*ln(x)
 */
     int i;
     n = dataseta.getIndex();
     if( n < 2 ){ return; } // no calculations for less than 2 data points
     double[] allX  = dataseta.getXPoints();
     double[] allY  = dataseta.getYPoints();
     double sumOfX  = 0.0;
     double sumOfY  = 0.0;
     double sumOfXY = 0.0;
     double sumOfX2 = 0.0;
     double sumOfY2 = 0.0;
     double sumOfYYP2 = 0.0;
     double t;
     if( Npts < n ){
       Npts = n;
       myControl.setValue("Npts", Npts); ; //in case there is a lot of data
      }
     xx = new double [Npts];
     yy = new double [Npts];
     y2 = new double [n];
     double minX = dataseta.getXMin();
     double maxX = dataseta.getXMax();
     if(M==1){//power fit case y=c*x^a
       for(i = 0; i < n; i++) {
       if(Math.abs(allX[i])<eps){allX[i]=eps;}//to prevent function problems
       if(Math.abs(allY[i])<eps){allY[i]=eps;}//to prevent function problems
       allX[i]=Math.log(Math.abs(allX[i])); //the new x's
       allY[i]=Math.log(Math.abs(allY[i])); //the new y's
       }
     }
     if(M==2){//exponential fit case y=c*exp(a*x)
       for(i = 0; i < n; i++) {
       if(Math.abs(allX[i])<eps){allX[i]=eps;}//to prevent function problems
       if(Math.abs(allY[i])<eps){allY[i]=eps;}//to prevent function problems
       allX[i]=allX[i];                     //the new x's
       allY[i]=Math.log(Math.abs(allY[i])); //the new y's
       }
     }
     for(i = 0; i < n; i++) {
       if(Math.abs(allX[i])<eps){allX[i]=eps;}//to prevent function problems
       if(Math.abs(allY[i])<eps){allY[i]=eps;}//to prevent function problems
       sumOfX += allX[i];
       sumOfY += allY[i];
       sumOfXY += (allX[i] * allY[i]);
       sumOfX2 += (allX[i] * allX[i]);
       sumOfY2 += (allY[i] * allY[i]);
     }
//apply method of linear least squares regression analysis
     m = ((n * (sumOfXY)) - ((sumOfX) * (sumOfY))) / ((n * (sumOfX2)) - ((sumOfX) * (sumOfX)));
     b = ((sumOfY) - (m * (sumOfX))) / n;
     r=Math.pow(sumOfXY-(sumOfX * sumOfY/n),2);
     r=r/(sumOfX2-sumOfX*sumOfX/n);
     r=r/(sumOfY2-sumOfY*sumOfY/n);//goodness of fit
     myControl.setValue("Goodness of Fit", r);
     datasetb.clear();
     datasetsb.clear();
//=============== linear fit y=a*x+b =============
     if(M==0){
       a=m;
       myControl.setValue("Model a",a);
       myControl.setValue("Model b",b);
       myControl.setValue("Model c",0.0);
       double xstep = (maxX - minX) / (Npts - 1);
        for (i = 0; i < Npts; i++) {
          t = minX + i * xstep;
          xx[i] = t;
          yy[i] = a * t + b;
          //System.out.println("i="+i+" xx="+xx[i]+" yy="+yy[i]); // if needed
        }
       for(i = 0; i < n; i++) {
           t = allX[i];
           y2[i] = m * t + b;
           sumOfYYP2 += ((allY[i] - y2[i]) * (allY[i] - y2[i]));
       }
       syx = Math.sqrt(sumOfYYP2 / (n - 1)); // cannot use n-2 (div by zero)
       myControl.setValue("Standard Error", syx);
       message = "y = " + nf.format(a) + " * x + " + nf.format(b);
     }
//=============== power fit y=c*x^a =============
     if(M==1){
       c=Math.exp(b);
       a=m;
       myControl.setValue("Model a",a);
       myControl.setValue("Model b",0.0);
       myControl.setValue("Model c",c);
       double xstep = (maxX - minX) / (Npts - 1);
        for (i = 0; i < Npts; i++) {
        t = minX + i * xstep;
        if(Math.abs(t)<eps){t=eps;}//to prevent function problems
        xx[i] = t;
        yy[i] = c * Math.pow(Math.abs(t), a);
        //System.out.println("i="+i+" xx="+xx[i]+" yy="+yy[i]); // if needed
        }
       for(i = 0; i < n; i++) {
           t = allX[i];
           y2[i] = c * Math.pow(Math.abs(t), a);
           sumOfYYP2 += ((allY[i] - y2[i]) * (allY[i] - y2[i]));
       }
       syx = Math.sqrt(sumOfYYP2 / (n - 1)); // cannot use n - 2 (div by zero)
       myControl.setValue("Standard Error", syx);
       message = "y = " + nf.format(c) + " * x ^ " + nf.format(a);
     }
 //=============== power fit y=c*exp(a*x) ==========
      if(M==2){
        c=Math.exp(b);
        a=m;
        myControl.setValue("Model a",a);
        myControl.setValue("Model b",0.0);
        myControl.setValue("Model c",c);
        double xstep = (maxX - minX) / (Npts - 1);
         for (i = 0; i < Npts; i++) {
         t = minX + i * xstep;
         if(Math.abs(t)<eps){t=eps;}//to prevent function problems
         xx[i] = t;
         yy[i] = c * Math.exp( a * t );
         //System.out.println("i="+i+" xx="+xx[i]+" yy="+yy[i]); // if needed
         }
        for(i = 0; i < n; i++) {
           t = allX[i];
           y2[i] = c * Math.exp( a * t );
           sumOfYYP2 += ((allY[i] - y2[i]) * (allY[i] - y2[i]));
        }
        syx = Math.sqrt(sumOfYYP2 / (n - 1));
        myControl.setValue("Standard Error", syx);
        message = "y = " + nf.format(c) + " * exp( " + nf.format(a) + " * x )";
      }
      datasetb.append(xx, yy);
      datasetsb.append(0,xx,yy);
   }

  /* This method was written by Max Perkins with Dr. J. Hasbun in the Physics Department
   * at the State University of West Georgia. Max Perkins was sponsored by the NASA Space
   * Grant Consortium through Dr. Ben de Mayo of the Physics Department. Spring 2004.
   */
   public void doQuadratic() {
        double[] allX  = dataseta.getXPoints();
        double[] allY  = dataseta.getYPoints();
        double minX = dataseta.getXMin();
        double maxX = dataseta.getXMax();
        double e1 = 0.0 , e2 = 0.0, e3 = 0.0, e4 = 0.0;
        double f0 = 0.0 , f1 = 0.0, f2 = 0.0, t = 0.0;
        double D = 0.0;
        double sumOfX  = 0.0;
        double sumOfY  = 0.0;
        double sumOfXY = 0.0;
        double sumOfX2 = 0.0;
        double sumOfY2 = 0.0;
        double sumOfYYP2 = 0.0;
        int i;
        n= dataseta.getIndex();
        if( Npts < n ){
          Npts = n;
          myControl.setValue("Npts", Npts); ; //in case there is a lot of data
         }
        if( Npts < 3 || M != 3) { return; } // no calculations for less than 3 data points
                                         // and if model number isn't 3
        for(i = 0; i < n; i++) {
            e1 += allX[i];
            e2 += (allX[i] * allX[i]);
            e3 += (allX[i] * allX[i] * allX[i]);
            e4 += (allX[i] * allX[i] * allX[i] * allX[i]);

            f0 += allY[i];
            f1 += (allY[i] * allX[i]);
            f2 += (allY[i] * (allX[i] * allX[i]));
        }

        D = (n * e2 * e4) + (2 * e1 * e2 * e3) - (e2 * e2 * e2) - (n * (e3 * e3 )) - ((e1 * e1) * e4);
        c = ( (((e2 * e4) - (e3 * e3)) * f0) + (((e2 * e3) - (e1 * e4)) * f1) + (((e1 * e3) - (e2 * e2)) * f2) ) / D;
        b = ( (((e2 * e3) - (e1 * e4)) * f0) + (((n * e4) - (e2 * e2)) * f1) + (((e1 * e2) - (n * e3)) * f2) ) / D;
        a = ( (((e1 * e3) - (e2 * e2)) * f0) + (((e1 * e2) - (n * e3)) * f1) + (((n * e2) - (e1 * e1)) * f2) ) / D;

       myControl.setValue("Model a", a);
       myControl.setValue("Model b", b);
       myControl.setValue("Model c", c);

       xx = new double[Npts];
       yy = new double[Npts];
       y2 = new double [n];
       datasetb.clear();
       datasetsb.clear();
       double xstep = (maxX - minX) / (Npts - 1);
       for (i = 0; i < Npts; i++) {
          t = minX + i * xstep;
          xx[i] = t;
          yy[i] = a * t * t + b * t + c;
       }
       datasetb.append(xx, yy);
       datasetsb.append(0,xx,yy);
       for(i = 0; i < n; i++) {
           t = allX[i];
           y2[i] = (a * (t * t)) + (b * t) + c;
           sumOfYYP2 += ((allY[i] - y2[i]) * (allY[i] - y2[i]));
       }
       syx = Math.sqrt(sumOfYYP2 / (n - 2));
       myControl.setValue("Standard Error", syx);
       myControl.setValue("Goodness of Fit"," ");
       message = "y = " + nf.format(a) + "*x^2 + " + nf.format(b) + "*x + " + nf.format(c);
   }

 public void calculate() {
// does the main calculations for line fitting
     //x = myControl.getDouble("x");
     //y = myControl.getDouble("y");
     /*
     plotPanel.clear();
     dataTablea.clear();
     dataTableb.clear();
     dataseta.append(x,y);
     datasetsa.append(0,x,y);//just one pair of columns
     if( M <= 2)
         doBestFitt2 ();
     else if( M == 3)
         doQuadratic ();
//--- plot data and fit
     plotPanel.addDrawable(dataseta);
     plotPanel.addDrawable(datasetb);
     plotPanel.setMessage(message);
     plotPanel.render();
     drawingFrame.show();
//--- calculation table
     datasetsb.clear();//clear previous calculated table
     datasetsb.append(0,xx,yy);//just one pair of columns
     dataTableb.add(datasetsb);
     dataTableb.refreshTable();
     frameb.show();
//--- entered data table
     dataTablea.add(datasetsa);
     //dataTablea.add(dataseta);//old line before EDatasetManager
     dataTablea.refreshTable();
     framea.show();*/

     switch(M) {
         case 0:
         case 1:
         case 2:
             doBestFitt2(); break;
         case 3:
             doQuadratic(); break;
     }
   }

   // Used in refreshDataset
   private boolean isRelevant(char theChar) {
     if(theChar == '1' || theChar == '2' || theChar == '3' ||
        theChar == '4' || theChar == '5' || theChar == '6' ||
        theChar == '7' || theChar == '8' || theChar == '9' ||
        theChar == '0' || theChar == '-' || theChar == '.' ||
        theChar == '+')
            return true;
     return false;
   }

   // Used in replot(), this gets the modified cell values
   // and takes out the spaces (if both columns in the row
   // are blank)
   private void refreshDataset() {
     double [] tx = datasetsa.getXPoints(0); // zero for 1st dataset (the only one)
     double [] ty = datasetsa.getYPoints(0);

     String xVal = null;
     String yVal = null;

     dataseta.clear();
     datasetsa.clear();

     for(int j = 0; j < tx.length; j++) {
        xVal = Double.toString(tx[j]);
        yVal = Double.toString(ty[j]);

        // if both x and y are blank, skip (this clears that part of the table,
        // putting the ones ahead in its place)
        if(!isRelevant(xVal.charAt(0)) && !isRelevant(yVal.charAt(0))){
            continue;
        }
        // else append
        else {
             dataseta.append(tx[j], ty[j]);
             datasetsa.append(0, tx[j], ty[j]);
        }
     }
   }

   public void replot() {
     refreshDataset();
     datasetsa.initDataTable(InitialRows);
     dataTablea.clear();
     dataTablea.add(datasetsa);
     dataTablea.refreshTable();

     calculate();

     dataTableb.clear();
     dataTableb.add(datasetsb);
     dataTableb.refreshTable();

     plotPanel.clear();
     plotPanel.addDrawable(dataseta);
     plotPanel.addDrawable(datasetb);
     plotPanel.setMessage(message);
     plotPanel.render();
     plotPanel.enableInspector(true);

     drawingFrame.show();
   }

   public void example() {
     double [] xe = {0.1,1.5,3.2,4.3,5.56,6.34,7.98,9.12,10.1,
                     10.8,11.2,12.1,13.4,14.56};
     double [] ye = {0.01,2.25,10.24,18.49,40.3,73.5,79.7,107.23,
                     123.9,131.7,157.23,192.1,214.19,253.283};
     plotPanel.clear();
     dataTablea.clear();
     dataTableb.clear();
     dataseta.clear();
     datasetsa.clear();
     dataseta.append(xe,ye);
     datasetsa.append(0,xe,ye);//just one pair of columns
// the fitter program
     calculate();
//--- plot data and fit
     plotPanel.addDrawable(dataseta);
     plotPanel.addDrawable(datasetb);
     plotPanel.setMessage(message);
     plotPanel.render();
//--- calculation table
     dataTableb.add(datasetsb);
     dataTableb.refreshTable();
     frameb.show();
//--- entered data table
     InitialRows = xe.length;
     myControl.setValue("Initial Rows", InitialRows);
     datasetsa.initDataTable(InitialRows);
     dataTablea.add(datasetsa);
     //dataTablea.add(dataseta);//old line before EDatasetManager
     dataTablea.refreshTable();
     framea.show();
   }

   public void clear  () {
     plotPanel.clear();
     message="";
     plotPanel.setMessage(message);
     plotPanel.render();
     dataTablea.clear();
     dataTablea.add(datasetsa);
     dataTablea.refreshTable();
     dataTableb.clear();
     dataTableb.refreshTable();
     dataseta.clear();
     datasetb.clear();
     datasetsa.clear();
     datasetsb.clear();
     resetCalculation();
   }

   public void initialize () {
     Npts = myControl.getInt("Npts");
     InitialRows = myControl.getInt("Initial Rows");
     M = myControl.getInt("Model number");
     if(M < 0 || M > 3){
        M=0;
        myControl.setValue("Model number",M);
     }
     if(M==0){
        Model="y=a*x+b";
     }
     if(M==1){
        Model="y=c*x^a";
     }
     if(M==2){
        Model="y=c*exp(a*x)";
     }
     if(M==3) {
        Model="y=a*x^2+b*x+c";
     }
     myControl.setValue("Model a"," ");
     myControl.setValue("Model b"," ");
     myControl.setValue("Model c"," ");
     myControl.setValue("Goodness of Fit"," ");
     myControl.setValue("Standard Error"," ");
     myControl.setValue("Present Model",Model);
     datasetsa.initDataTable(InitialRows);
     dataTablea.refreshTable();
   }

   public void resetCalculation() {
     myControl.clearMessages();
     myControl.println ("Set 'Model number' to 0 for linear model y=a*x+b");
     myControl.println ("Set 'Model number' to 1 for nonlinear power fit model y=c*x^a");
     myControl.println ("Set 'Model number' to 2 for nonlinear exp fit model y=c*exp(a*x)");
     myControl.println ("Set 'Model number' to 3 for quadratic fit model y=a*x^2+b*x+c");
     myControl.println ("Goodness of fit applies only for 'Model number' between 0 and 2");
     myControl.println ("Change 'Initial Rows' for more input data table rows");
     myControl.println ("Change 'Npts' for more Line Fit points");
     myControl.println ("'Init' if the parameters are changed");
     myControl.println ("'Replot' if the model or input data are changed");
     myControl.println ("'Ex' does an example");
     myControl.println ("'Clr' clears everything");
     n=0;
     x=1.0;
     y=1.0;
     m=0.0;//slope in linear model
     b=0.0;//intercept in linear model
     a=0.0;//parameter used in linearized models
     b=0.0;
     c=0.0;//parameter used in linearized models
     r=0.0;//Goodness of Fit
     syx=0.0;//Standard Error
     xx=null;
     yy=null;
     myControl.setValue("Model a", a);
     myControl.setValue("Model b", b);
     myControl.setValue("Model c", c);
     myControl.setValue("Goodness of Fit",r);
     myControl.setValue("Standard Error",syx);
     InitialRows = 100;
     Npts=10;
     M=0;//model number
     Model="y=a*x+b";//present model
     myControl.setValue("Present Model",Model);
     myControl.setValue("Initial Rows",InitialRows);
     myControl.setValue("Npts",Npts);
     myControl.setValue("Model number",M);
     initialize ();
   }

   public void setControl(Control control) {
     if(org.opensourcephysics.display.OSPFrame.appletMode) {
         initAppletMode();
         myControl = myAppletControl;
     }
     else
         myControl = control;
     resetCalculation();
   }

   /*
   * Needed for the applet
   */
   void initAppletMode (){
       if(org.opensourcephysics.display.OSPFrame.appletMode) {
        myAppletControl = new EAppletOSPControl(this);
        myAppletControl.addButton("resetCalculation", "         Reset         ", "Resets the initial conditions");
        myAppletControl.addButton("initialize", "        Initialize        ", "Initializes parameters");
        myAppletControl.addButton("replot", "         Replot         ","Replots");
        myAppletControl.addButton("example", "                   Example                   ","Does a random example");
        myAppletControl.addButton("clear", "                    Clear                   ", "Clears data and plot");
        myAppletControl.setLocation(415, 20);//785
        myAppletControl.setSize(397, 550);
        myAppletControl.setButtonPanelSize(0, 68);
        myAppletControl.setVisible(true);
        myAppletControl.show();
       }
     /*
     if (!org.opensourcephysics.display.OSPFrame.appletMode)return;
     if (myControl instanceof OSPControl) {
     ( (OSPControl) myControl).addButton("resetCalculation", "Reset",
                                         "Resets the initial conditions");
     ( (OSPControl) myControl).addButton("initialize", "Init",
                                         "Initializes parameters");
     //( (OSPControl) myControl).addButton("calculate", "+",
     //                                    "adds data to be analyzed");
     ( (OSPControl) myControl).addButton("replot", "Replot",
                                         "Replots");
     ( (OSPControl) myControl).addButton("example", "Ex",
                                         "Does an example");
     ( (OSPControl) myControl).addButton("clear", "Clr",
                                         "Clears data and plot");
     ( (OSPControl)myControl).setLocation(440, 30);//785
     ( (OSPControl)myControl).setSize(380, 550);
     ( (OSPControl)myControl).setVisible(true);
     ( (OSPControl)myControl).setDefaultCloseOperation(OSPControl.DO_NOTHING_ON_CLOSE);
     ( (OSPControl)myControl).show();
     }*/
   }

   public static void main(String[] args) {
     Calculation model = new QERegressionApp();
     EOSPControl myControl;

     // if you're not in applet mode, use regular OSPControl (er, the edited one: EOSPControl)
     if (!org.opensourcephysics.display.OSPFrame.appletMode) {
        myControl = new EOSPControl(model);
        myControl.addButton("resetCalculation", "         Reset         ", "Resets the initial conditions");
        myControl.addButton("initialize", "        Initialize        ", "Initializes parameters");
        myControl.addButton("replot", "         Replot         ","Replots");
        myControl.addButton("example", "                   Example                   ","Does a random example");
        myControl.addButton("clear", "                    Clear                   ", "Clears data and plot");
        myControl.setLocation(400, 20);//785
        myControl.setSize(397, 550);
        myControl.setButtonPanelSize(0, 68);
        myControl.setVisible(true);
        myControl.show();
        model.setControl(myControl);
     }
     /*
     Calculation model = new QERegressionApp();
     OSPControl myControl = new OSPControl(model); // makes it possible to use Buttons in applets
     myControl.addButton("resetCalculation", "Reset", "Resets the initial conditions");
     myControl.addButton("initialize", "Init", "Initializes parameters");
     //myControl.addButton("calculate", "+", "Adds data to be analyzed");
     myControl.addButton("replot", "Replot","Replots");
     myControl.addButton("example", "Ex","Does a random example");
     myControl.addButton("clear", "Clr", "Clears data and plot");
     myControl.setLocation(440, 30);//785
     //myControl.setSize(380, 550);
     myControl.setSize(370, 550);
     myControl.setVisible(true);
     myControl.show();
     model.setControl(myControl);*/
   }
}