// $Header: /home/cvs/jakarta-jmeter/src/components/org/apache/jmeter/visualizers/Spline3.java,v 1.8 2004/02/13 01:48:46 sebb Exp $
/*
 * Copyright 2001-2004 The Apache Software Foundation.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 * 
*/

package org.apache.jmeter.visualizers;

import org.apache.jorphan.logging.LoggingManager;
import org.apache.log.Logger;


/*
 * TODO :
 *  - implement ImageProducer interface
 *  - suggestions ;-)
 */

/**
 *  This class implements the representation of an interpolated Spline curve.
 *  

 *  The curve described by such an object interpolates an arbitrary number
 *  of fixed points called nodes. The distance between two nodes
 *  should currently be constant. This is about to change in a later version
 *  but it can last a while as it's not really needed.
 *  Nevertheless, if you need the feature, just
 *  write me a note
 *  and I'll write it asap.
 *  

 *  The interpolated Spline curve can't be described by an polynomial analytic
 *  equation, the degree of which would be as high as the number of nodes,
 *  which would cause extreme oscillations of the curve on the edges.
 *  

 *  The solution is to split the curve accross a lot of little
 *  intervals : an interval starts at one node
 *  and ends at the next one. Then, the interpolation is done on each
 *  interval, according to the following conditions :
 *  

 *   
the interpolated curve is degree 3 : it's a cubic curve ;
 *   
the interpolated curve contains the two points delimiting
 *    the interval. This condition obviously implies the curve is continuous ;
 *   
the interpolated curve has a smooth slope : the curvature has
 *    to be the same on the left and the right sides of each node ;
 *   
the curvature of the global curve is 0 at both edges.
 *  
 *  Every part of the global curve is represented by a cubic (degree-3)
 *  polynomial, the coefficients of which have to be computed in order
 *  to meet the above conditions.
 *  

 *  This leads to a n-unknow n-equation system to resolve.
 *  One can resolve an equation system by several manners ;
 *  this class uses the Jacobi iterative method, particularly well
 *  adapted to this situation, as the diagonal of the system matrix
 *  is strong compared to the other elements. This implies the algorithm
 *  always converges ! This is not the case of the Gauss-Seidel algorithm,
 *  which is quite faster (it uses intermediate results of each iteration
 *  to speed up the convergence) but it doesn't converge in all the cases
 *  or it converges to a wrong value. This is not acceptable and
 *  that's why the Jacobi method is safer. Anyway, the gain of speed
 *  is about a factor of 3 but, for a 100x100 system, it means 10 ms
 *  instead of 30 ms, which is a pretty good reason not to explore the question
 *  any further :)
 *  

 *  Here is a little piece of code showing how to use this class :
 *  
 *  // ...
 *  float[] nodes = {3F, 2F, 4F, 1F, 2.5F, 5F, 3F};
 *  Spline3 curve = new Spline3(nodes);
 *  // ...
 *  public void paint(Graphics g) {
 *      int[] plot = curve.getPlots();
 *      for (int i = 1; i < n; i++) {
 *          g.drawLine(i - 1, plot[i - 1], i, plot[i]);
 *      }
 *  }
 *  // ...
 *  
 *  Have fun with it !

 *  Any comments, feedback, bug reports or suggestions will be
 *  appreciated.
 *
 *  @author Jean-Pierre Norguet
 *  @version $Revison$ updated $Date: 2004/02/13 01:48:46 $
 */
public class Spline3
{
    transient private static Logger log = LoggingManager.getLoggerForClass();

    protected float[][] _coefficients;
    protected float[][] _A;
    protected float[] _B;
    protected float[] _r;
    protected float[] _rS;
    protected int _m;  // number of nodes
    protected int _n;  // number of non extreme nodes (_m-2)
    final static protected float DEFAULT_PRECISION = (float) 1E-1;
    final static protected int DEFAULT_MAX_ITERATIONS = 100;
    protected float _minPrecision = DEFAULT_PRECISION;
    protected int _maxIterations = DEFAULT_MAX_ITERATIONS;

    /**
     *  Creates a new Spline curve by calculating the coefficients of
     *  each part of the curve, i.e. by resolving the equation system
     *  implied by the interpolation condition on every interval.
     *  @param r an array of float containing the vertical coordinates
     *   of the nodes to interpolate ; the vertical coordinates start at 0
     *   and are equidistant with a step of 1.
     */
    public Spline3(float[] r)
    {
        int n = r.length;

        // the number of nodes is defined by the length of r
        this._m = n;
        // grab the nodes
        this._r = new float[n];
        for (int i = 0; i < n; i++)
        {
            _r[i] = r[i];
        }
        // the number of non extreme nodes is the number of intervals
        // minus 1, i.e. the length of r minus 2
        this._n = n - 2;
        // computes interpolation coefficients
        try
        {
            long startTime = System.currentTimeMillis();

            this.interpolation();
            if (log.isDebugEnabled())
            {
                long endTime = System.currentTimeMillis();
                long elapsedTime = endTime - startTime;

                log.debug("New Spline curve interpolated in ");
                log.debug(elapsedTime + " ms");
            }
        }
        catch (Exception e)
        {
            log.error("Error when interpolating : ", e);
        }

    }

    /**
     *  Computes the coefficients of the Spline interpolated curve,
     *  on each interval. The matrix system to resolve is AX=B
     */
    protected void interpolation()
    {
        // creation of the interpolation structure
        _rS = new float[_m];
        _B = new float[_n];
        _A = new float[_n][_n];
        _coefficients = new float[_n + 1][4];
        // local variables
        int i = 0, j = 0;

        // initialize system structures (just to be safe)
        for (i = 0; i < _n; i++)
        {
            _B[i] = 0;
            for (j = 0; j < _n; j++)
            {
                _A[i][j] = 0;
            }
            for (j = 0; j < 4; j++)
            {
                _coefficients[i][j] = 0;
            } 
        }
        for (i = 0; i < _n; i++)
        {
            _rS[i] = 0;
        }
        // initialize the diagonal of the system matrix (A) to 4
        for (i = 0; i < _n; i++)
        {
            _A[i][i] = 4;
        }
        // initialize the two minor diagonals of A to 1
        for (i = 1; i < _n; i++)
        {
            _A[i][i - 1] = 1;
            _A[i - 1][i] = 1;
        }
        // initialize B
        for (i = 0; i < _n; i++)
        {
            _B[i] = 6 * (_r[i + 2] - 2 * _r[i + 1] + _r[i]);
        }
        // Jacobi system resolving
        this.jacobi();  // results are stored in _rS
        // computes the coefficients (di, ci, bi, ai) from the results
        for (i = 0; i < _n + 1; i++)
        {
            // di (degree 0)
            _coefficients[i][0] = _r[i];
            // ci (degree 1)
            _coefficients[i][1] = _r[i + 1] - _r[i]
                    - (_rS[i + 1] + 2 * _rS[i]) / 6;
            // bi (degree 2)
            _coefficients[i][2] = _rS[i] / 2;
            // ai (degree 3)
            _coefficients[i][3] = (_rS[i + 1] - _rS[i]) / 6;
        }
    }

    /**
     *  Resolves the equation system by a Jacobi algorithm.
     *  The use of the slower Jacobi algorithm instead of Gauss-Seidel is
     *  choosen here because Jacobi is assured of to be convergent
     *  for this particular equation system, as the system matrix
     *  has a strong diagonal.
     */
    protected void jacobi()
    {
        // local variables
        int i = 0, j = 0, iterations = 0;
        // intermediate arrays
        float[] newX = new float[_n];
        float[] oldX = new float[_n];

        // Jacobi convergence test
        if (!converge())
        {
            if (log.isDebugEnabled())
            {
                log.debug(
                    "Warning : equation system resolving is unstable");
            }
        }
        // init newX and oldX arrays to 0
        for (i = 0; i < _n; i++)
        {
            newX[i] = 0;
            oldX[i] = 0;
        }
        // main iteration
        while ((this.precision(oldX, newX) > this._minPrecision)
                && (iterations < this._maxIterations))
        {
            for (i = 0; i < _n; i++)
            {
                oldX[i] = newX[i];
            }
            for (i = 0; i < _n; i++)
            {
                newX[i] = _B[i];
                for (j = 0; j < i; j++)
                {
                    newX[i] = newX[i] - (_A[i][j] * oldX[j]);
                }
                for (j = i + 1; j < _n; j++)
                {
                    newX[i] = newX[i] - (_A[i][j] * oldX[j]);
                }
                newX[i] = newX[i] / _A[i][i];
            }
            iterations++;
        }
        if (this.precision(oldX, newX) < this._minPrecision)
        {
            if (log.isDebugEnabled())
            {
                log.debug("Minimal precision (");
                log.debug(this._minPrecision + ") reached after ");
                log.debug(iterations + " iterations");
            }
        }
        else if (iterations > this._maxIterations)
        {
            if (log.isDebugEnabled())
            {
                log.debug("Maximal number of iterations (");
                log.debug(this._maxIterations + ") reached");
                log.debug("Warning : precision is only ");
                log.debug("" + this.precision(oldX, newX));
                log.debug(", divergence is possible");
            }
        }
        for (i = 0; i < _n; i++)
        {
            _rS[i + 1] = newX[i];
        }
    }

    /**
     *  Test if the Jacobi resolution of the equation system converges.
     *  It's OK if A has a strong diagonal.
     */
    protected boolean converge()
    {
        boolean converge = true;
        int i = 0, j = 0;
        float lineSum = 0F;

        for (i = 0; i < _n; i++)
        {
            if (converge)
            {
                lineSum = 0;
                for (j = 0; j < _n; j++)
                {
                    lineSum = lineSum + Math.abs(_A[i][j]);
                }
                lineSum = lineSum - Math.abs(_A[i][i]);
                if (lineSum > Math.abs(_A[i][i]))
                {
                    converge = false;
                }
            }
        }
        return converge;
    }

    /**
     *  Computes the current precision reached.
     */
    protected float precision(float[] oldX, float[] newX)
    {
        float N = 0F, D = 0F, erreur = 0F;
        int i = 0;

        for (i = 0; i < _n; i++)
        {
            N = N + Math.abs(newX[i] - oldX[i]);
            D = D + Math.abs(newX[i]);
        }
        if (D != 0F)
        {
            erreur = N / D;
        }
        else
        {
            erreur = Float.MAX_VALUE;
        }
        return erreur;
    }

    /**
     *  Computes a (vertical) Y-axis value of the global curve.
     *  @param t abscissa
     *  @return computed ordinate
     */
    public float value(float t)
    {
        int i = 0, splineNumber = 0;
        float abscissa = 0F, result = 0F;

        // verify t belongs to the curve (range [0, _m-1])
        if ((t < 0) || (t > (_m - 1)))
        {
            if (log.isDebugEnabled())
            {
                log.debug(
                    "Warning : abscissa "
                        + t
                        + " out of bounds [0, "
                        + (_m - 1)
                        + "]");
            }
            // silent error, consider the curve is constant outside its range
            if (t < 0)
            {
                t = 0;
            }
            else
            {
                t = _m - 1;
            }
        }
        // seek the good interval for t and get the piece of curve on it
        splineNumber = (int) Math.floor(t);
        if (t == (_m - 1))
        {
            // the upper limit of the curve range belongs by definition
            // to the last interval
            splineNumber--;
        }
        // computes the value of the curve at the pecified abscissa
        // and relative to the beginning of the right piece of Spline curve
        abscissa = t - splineNumber;
        // the polynomial calculation is done by the (fast) Euler method
        for (i = 0; i < 4; i++)
        {
            result = result * abscissa;
            result = result + _coefficients[splineNumber][3 - i];
        }
        return result;
    }

    /**
     *  Manual check of the curve at the interpolated points.
     */
    public void debugCheck()
    {
        int i = 0;

        for (i = 0; i < _m; i++)
        {
            log.info("Point " + i + " : ");
            log.info(_r[i] + " =? " + value(i));
        }
    }

    /**
     *  Computes drawable plots from the curve for a given draw space.
     *  The values returned are drawable vertically and from the
     *  bottom of a Panel.
     *  @param width width within the plots have to be computed
     *  @param height height within the plots are expected to be drawed
     *  @return drawable plots within the limits defined, in an array of
     *   int (as many int as the value of the width parameter)
     */
    public int[] getPlots(int width, int height)
    {
        int[] plot = new int[width];
        // computes auto-scaling and absolute plots
        float[] y = new float[width];
        float max = java.lang.Integer.MIN_VALUE;
        float min = java.lang.Integer.MAX_VALUE;

        for (int i = 0; i < width; i++)
        {
            y[i] = value(((float) i) * (_m - 1) / width);
            if (y[i] < min)
            {
                min = y[i];
            }
             
            if (y[i] > max)
            {
                max = y[i];
            } 
        }
        if (min < 0)
        {
            min = 0; // shouldn't draw negative values
        }
        // computes relative auto-scaled plots to fit in the specified area
        for (int i = 0; i < width; i++)
        {
            plot[i] =
                (int) Math.round(((y[i] - min) * (height - 1)) / (max - min));
        }
        return plot;
    }

    public void setPrecision(float precision)
    {
        this._minPrecision = precision;
    }

    public float getPrecision()
    {
        return this._minPrecision;
    }

    public void setToDefaultPrecision()
    {
        this._minPrecision = DEFAULT_PRECISION;
    }

    public float getDefaultPrecision()
    {
        return DEFAULT_PRECISION;
    }

    public void setMaxIterations(int iterations)
    {
        this._maxIterations = iterations;
    }

    public int getMaxIterations()
    {
        return this._maxIterations;
    }

    public void setToDefaultMaxIterations()
    {
        this._maxIterations = DEFAULT_MAX_ITERATIONS;
    }

    public int getDefaultMaxIterations()
    {
        return DEFAULT_MAX_ITERATIONS;
    }

}