/*      Example of the use of the class Regression demonstrating
*       the use of the non-linear regression method, Regression.simplexPlot
*       in fitting data to the function, y = a + b.exp(-c.x1) + f.exp(-g.x2)
*       with a fixed and b, c, f and g unknown
*       The second derivatives used in calculating the Covariance matrix
*       are calculated analytically 
*
*       Michael Thomas Flanagan
*       http://www.ee.ucl.ac.uk/~mflanaga/java/Regression.html
*       January 2012
*
*       
*   Copyright (c) 2012 - 2014

*
*   PERMISSION TO COPY:
*   Permission to use, copy and modify this software and its documentation for 
*   NON-COMMERCIAL purposes is granted, without fee, provided that an acknowledgement 
*   to the author, Dr Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, 
*   appears in all copies and associated documentation or publications. 
*   
*   Public listing of the source codes on the internet is not permitted. 
*   
*   Redistribution of the source codes or of the flanagan.jar file is not permitted. 
*   
*   Redistribution in binary form of all or parts of these classes is not permitted. 
*   
*   Dr Michael Thomas Flanagan makes no representations about the suitability 
*   or fitness of the software for any or for a particular purpose. 
*   Dr Michael Thomas Flanagan shall not be liable for any damages suffered 
*   as a result of using, modifying or distributing this software or its derivatives.
*
***************************************************************************************/

// Class to evaluate the function y = a + b.exp(-c.x1) + f.exp(-g.x2)
// where a is fixed and best estimates are required for b, c, f and g.
import flanagan.analysis.*;

class Funct implements RegressionFunction{

        public double a=0.0D;

        public double function(double[] p, double[] x){
            double y = a + p[0]*Math.exp(-p[1]*x[0]) + p[2]*Math.exp(-p[3]*x[1]);
            return y;
        }
}

// Class to evaluate dy/db, dy/dc, dy/df, dy/dg, d2y/db2, d2y/dbdc, d2y/dbdf,
//  d2y/dbdg, d2y/dc2, d2y/dcdf, d2y/dcdg, d2y/df2, d2y/dfdg and d2y/dg2
// for the function y = a + b.exp(-c.x1) + f.exp(-g.x2)
// where a is fixed and best estimates are required for b, c, d and f.
class FunctDeriv implements RegressionDerivativeFunction{

        public double[] function(double[] p, double[] x, int ii, int jj){
            double[] derivs = new double[3];
            
            switch(ii){
                case 0: derivs[0] = Math.exp(-p[1]*x[0]);
                        switch(jj){
                            case 0: derivs[1] = Math.exp(-p[1]*x[0]);
                                    derivs[2] = 0.0;
                                    break;
                            case 1: derivs[1] = -p[0]*x[0]*Math.exp(-p[1]*x[0]);   
                                    derivs[2] = x[0]*Math.exp(-p[1]*x[0]);
                                    break;
                            case 2: derivs[1] = Math.exp(-p[3]*x[1]);  
                                    derivs[2] = 0.0;
                                    break;
                            case 3: derivs[1] = -p[2]*x[1]*Math.exp(-p[3]*x[1]);
                                    derivs[2] = 0.0;
                                    break;
                        }
                        break;
                case 1: derivs[0] = -p[0]*x[0]*Math.exp(-p[1]*x[0]);
                        switch(jj){
                            case 0: derivs[1] = Math.exp(-p[1]*x[0]);
                                    derivs[2] = x[0]*Math.exp(-p[1]*x[0]);
                                    break;
                            case 1: derivs[1] = -p[0]*x[0]*Math.exp(-p[1]*x[0]);   
                                    derivs[2] = p[0]*x[0]*x[0]*Math.exp(-p[1]*x[0]);
                                    break;
                            case 2: derivs[1] = Math.exp(-p[3]*x[1]);  
                                    derivs[2] = 0.0;
                                    break;
                            case 3: derivs[1] = -p[2]*x[1]*Math.exp(-p[3]*x[1]);
                                    derivs[2] = 0.0;
                                    break;
                        }
                        break;
                case 2: derivs[0] = Math.exp(-p[3]*x[1]); 
                        switch(jj){
                            case 0: derivs[1] = Math.exp(-p[1]*x[0]);
                                    derivs[2] = 0.0;
                                    break;
                            case 1: derivs[1] = -p[0]*x[0]*Math.exp(-p[1]*x[0]);   
                                    derivs[2] = 0.0;
                                    break;
                            case 2: derivs[1] = Math.exp(-p[3]*x[1]);  
                                    derivs[2] = 0.0;
                                    break;
                            case 3: derivs[1] = -p[2]*x[1]*Math.exp(-p[3]*x[1]);
                                    derivs[2] = -x[1]*Math.exp(-p[3]*x[1]);
                                    break;
                        }
                        break;
                case 3: derivs[0] = -p[2]*x[1]*Math.exp(-p[3]*x[1]);
                        switch(jj){
                            case 0: derivs[1] = Math.exp(-p[1]*x[0]);
                                    derivs[2] = 0.0;
                                    break;
                            case 1: derivs[1] = -p[0]*x[0]*Math.exp(-p[1]*x[0]);   
                                    derivs[2] = 0.0;
                                    break;
                            case 2: derivs[1] = Math.exp(-p[3]*x[1]);  
                                    derivs[2] = -x[1]*Math.exp(-p[3]*x[1]);
                                    break;
                            case 3: derivs[1] = -p[2]*x[1]*Math.exp(-p[3]*x[1]);
                                    derivs[2] = p[2]*x[1]*x[1]*Math.exp(-p[3]*x[1]);
                                    break;
                        }
                        break;
            }
            return derivs;
        }
}

// Class to fit data to the function y = a + b.exp(-c.x1) + f.exp(-g.x2)
// where a is fixed and best estimates are required for b, c, f and g.

public class RegressionExampleThree{

   public static void main(String[] arg){
        // x1 data array
        double[] x1 = {0.0D, 0.5D, 1.0D, 1.5D, 2.0D, 2.5D, 3.0D, 3.5D};
        // x2 data array
        double[] x2 = {0.0D, 3.0D, 6.0D, 9.0D, 12.0D, 15.0D, 18.0D, 21.0D};
        // observed y data array
        double[] yArray = {13.9D, 12.0D, 9.8D, 10.01D, 8.0D, 7.91D, 6.91D, 6.8D};
        // estimates of the standard  deviations of y
        double[] sdArray = {1.0D, 0.9D, 1.1D, 0.95D, 0.92D, 1.04D, 1.1D, 0.99D};

        // combine x1 and x2 data into common array
        int arrayLen = x1.length;
        double[][] xArray = new double[2][arrayLen];
        for(int i=0; i<arrayLen; i++){
            xArray[0][i]=x1[i];
            xArray[1][i]=x2[i];
        }

        // Create instances of the class holding the function, y = a + b.exp(-c.x1) + f.exp(-g.x2), evaluation method
        Funct funct = new Funct();

        // Assign values to constant a in the function y = a + b.exp(-c.x1) + f.exp(-g.x2)
        funct.a = 3.0D;
        
         // Create instances of the class holding the derivatives of the function, y = a + b.exp(-c.x1) + f.exp(-g.x2), evaluation method
        FunctDeriv functd = new FunctDeriv();

        //Initial estimates of b, c d and f in y = a + b.exp(-c.x1) + f.exp(-g.x2)
        double[] start = new double[4];
        start[0] = 8.0D;    // initial estimate of b
        start[1] = 0.5D;    // initial estimate of c
        start[2] = 3.0D;    // initial estimate of d
        start[3] = 0.1D;    // initial estimate of f

        //Initial step sizes for b, c d and f in y = a + b.exp(-c.x1) + f.exp(-g.x2)
        double[] step = new double[4];
        step[0] = 0.4D;    // step size for b
        step[1] = 0.2D;    // step size for c
        step[2] = 0.4D;    // step size for d
        step[3] = 0.2D;    // step size for f


        // create an instance of regression
        Regression reg = new Regression(xArray, yArray, sdArray);
        
        // call non-linear regression using default tolerance and maximum iterations
        reg.simplexPlot(funct, functd, start, step);

    }

}