<?php 
/*=======================================================================
// File:    JPGRAPH_REGSTAT.PHP
// Description: Regression and statistical analysis helper classes
// Created:     2002-12-01
// Author:    Johan Persson (johanp@aditus.nu)
// Ver:        $Id: jpgraph_regstat.php 464 2006-02-04 12:07:57Z ljp $
//
// Copyright (c) Aditus Consulting. All rights reserved.
//========================================================================
*/

//------------------------------------------------------------------------
// CLASS Spline
// Create a new data array from an existing data array but with more points.
// The new points are interpolated using a cubic spline algorithm
//------------------------------------------------------------------------
class Spline {
    // 3:rd degree polynom approximation

    private $xdata,$ydata;   // Data vectors
    private $y2;         // 2:nd derivate of ydata    
    private $n=0;

    function Spline($xdata,$ydata) {
    $this->y2 = array();
    $this->xdata = $xdata;
    $this->ydata = $ydata;

    $n = count($ydata);
    $this->n = $n;
    if( $this->n !== count($xdata) ) {
        JpGraphError::RaiseL(19001);
//('Spline: Number of X and Y coordinates must be the same');
    }

    // Natural spline 2:derivate == 0 at endpoints
    $this->y2[0]    = 0.0;
    $this->y2[$n-1] = 0.0;
    $delta[0] = 0.0;

    // Calculate 2:nd derivate
    for($i=1; $i < $n-1; ++$i) {
        $d = ($xdata[$i+1]-$xdata[$i-1]);
        if( $d == 0  ) {
        JpGraphError::RaiseL(19002);
//('Invalid input data for spline. Two or more consecutive input X-values are equal. Each input X-value must differ since from a mathematical point of view it must be a one-to-one mapping, i.e. each X-value must correspond to exactly one Y-value.');
        }
        $s = ($xdata[$i]-$xdata[$i-1])/$d;
        $p = $s*$this->y2[$i-1]+2.0;
        $this->y2[$i] = ($s-1.0)/$p;
        $delta[$i] = ($ydata[$i+1]-$ydata[$i])/($xdata[$i+1]-$xdata[$i]) - 
                 ($ydata[$i]-$ydata[$i-1])/($xdata[$i]-$xdata[$i-1]);
        $delta[$i] = (6.0*$delta[$i]/($xdata[$i+1]-$xdata[$i-1])-$s*$delta[$i-1])/$p;
    }

    // Backward substitution
    for( $j=$n-2; $j >= 0; --$j ) {
        $this->y2[$j] = $this->y2[$j]*$this->y2[$j+1] + $delta[$j];
    }
    }

    // Return the two new data vectors
    function Get($num=50) {
    $n = $this->n ;
    $step = ($this->xdata[$n-1]-$this->xdata[0]) / ($num-1);
    $xnew=array();
    $ynew=array();
    $xnew[0] = $this->xdata[0];
    $ynew[0] = $this->ydata[0];
    for( $j=1; $j < $num; ++$j ) {
        $xnew[$j] = $xnew[0]+$j*$step;
        $ynew[$j] = $this->Interpolate($xnew[$j]);
    }
    return array($xnew,$ynew);
    }

    // Return a single interpolated Y-value from an x value
    function Interpolate($xpoint) {

    $max = $this->n-1;
    $min = 0;

    // Binary search to find interval
    while( $max-$min > 1 ) {
        $k = ($max+$min) / 2;
        if( $this->xdata[$k] > $xpoint ) 
        $max=$k;
        else 
        $min=$k;
    }    

    // Each interval is interpolated by a 3:degree polynom function
    $h = $this->xdata[$max]-$this->xdata[$min];

    if( $h == 0  ) {
        JpGraphError::RaiseL(19002);
//('Invalid input data for spline. Two or more consecutive input X-values are equal. Each input X-value must differ since from a mathematical point of view it must be a one-to-one mapping, i.e. each X-value must correspond to exactly one Y-value.');
    }


    $a = ($this->xdata[$max]-$xpoint)/$h;
    $b = ($xpoint-$this->xdata[$min])/$h;
    return $a*$this->ydata[$min]+$b*$this->ydata[$max]+
         (($a*$a*$a-$a)*$this->y2[$min]+($b*$b*$b-$b)*$this->y2[$max])*($h*$h)/6.0;
    }
}

//------------------------------------------------------------------------
// CLASS Bezier
// Create a new data array from a number of control points
//------------------------------------------------------------------------
class Bezier {
/**
 * @author Thomas Despoix, openXtrem company
 * @license released under QPL
 * @abstract Bezier interoplated point generation,
 * computed from control points data sets, based on Paul Bourke algorithm :
 * http://astronomy.swin.edu.au/~pbourke/curves/bezier/
 */
    private $datax = array();
    private $datay = array();
    private $n=0;
 
    function Bezier($datax, $datay, $attraction_factor = 1) {
    // Adding control point multiple time will raise their attraction power over the curve
    $this->n = count($datax);
    if( $this->n !== count($datay) ) {
        JpGraphError::RaiseL(19003);
//('Bezier: Number of X and Y coordinates must be the same');
    }
    $idx=0;
    foreach($datax as $datumx) {
        for ($i = 0; $i < $attraction_factor; $i++) {
        $this->datax[$idx++] = $datumx;
        }
    }
       $idx=0;
    foreach($datay as $datumy) {
        for ($i = 0; $i < $attraction_factor; $i++) {
        $this->datay[$idx++] = $datumy;
        }
    }
    }

    function Get($steps) {
    $datax = array();
    $datay = array();
    for ($i = 0; $i < $steps; $i++) {
        list($datumx, $datumy) = $this->GetPoint((double) $i / (double) $steps);       
        $datax[] = $datumx;
        $datay[] = $datumy;
    }
   
    $datax[] = end($this->datax);
    $datay[] = end($this->datay);
   
    return array($datax, $datay);
    }
 
    function GetPoint($mu) {
    $n = $this->n - 1;
    $k = 0;
    $kn = 0;
    $nn = 0;
    $nkn = 0;
    $blend = 0.0;
    $newx = 0.0;
    $newy = 0.0;

    $muk = 1.0;
    $munk = (double) pow(1-$mu,(double) $n);

    for ($k = 0; $k <= $n; $k++) {
        $nn = $n;
        $kn = $k;
        $nkn = $n - $k;
        $blend = $muk * $munk;
        $muk *= $mu;
        $munk /= (1-$mu);
        while ($nn >= 1) {
        $blend *= $nn;
        $nn--;
        if ($kn > 1) {
            $blend /= (double) $kn;
            $kn--;
        }
        if ($nkn > 1) {
            $blend /= (double) $nkn;
            $nkn--;
        }
        }
        $newx += $this->datax[$k] * $blend;
        $newy += $this->datay[$k] * $blend;
    }

    return array($newx, $newy);
    }
}

// EOF
?>