[Sommaire] [Imprimer]

   1  <?php 
   2  /*=======================================================================
   3  // File:    JPGRAPH_REGSTAT.PHP
   4  // Description: Regression and statistical analysis helper classes
   5  // Created:     2002-12-01
   6  // Ver:        $Id: jpgraph_regstat.php 781 2006-10-08 08:07:47Z ljp $
   7  //
   8  // Copyright (c) Aditus Consulting. All rights reserved.
   9  //========================================================================
  10  */
  11  
  12  //------------------------------------------------------------------------
  13  // CLASS Spline
  14  // Create a new data array from an existing data array but with more points.
  15  // The new points are interpolated using a cubic spline algorithm
  16  //------------------------------------------------------------------------
  17  class Spline {
  18      // 3:rd degree polynom approximation
  19  
  20      private $xdata,$ydata;   // Data vectors
  21      private $y2;         // 2:nd derivate of ydata    
  22      private $n=0;
  23  
  24      function Spline($xdata,$ydata) {
  25      $this->y2 = array();
  26      $this->xdata = $xdata;
  27      $this->ydata = $ydata;
  28  
  29      $n = count($ydata);
  30      $this->n = $n;
  31      if( $this->n !== count($xdata) ) {
  32          JpGraphError::RaiseL(19001);
  33  //('Spline: Number of X and Y coordinates must be the same');
  34      }
  35  
  36      // Natural spline 2:derivate == 0 at endpoints
  37      $this->y2[0]    = 0.0;
  38      $this->y2[$n-1] = 0.0;
  39      $delta[0] = 0.0;
  40  
  41      // Calculate 2:nd derivate
  42      for($i=1; $i < $n-1; ++$i) {
  43          $d = ($xdata[$i+1]-$xdata[$i-1]);
  44          if( $d == 0  ) {
  45          JpGraphError::RaiseL(19002);
  46  //('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.');
  47          }
  48          $s = ($xdata[$i]-$xdata[$i-1])/$d;
  49          $p = $s*$this->y2[$i-1]+2.0;
  50          $this->y2[$i] = ($s-1.0)/$p;
  51          $delta[$i] = ($ydata[$i+1]-$ydata[$i])/($xdata[$i+1]-$xdata[$i]) - 
  52                   ($ydata[$i]-$ydata[$i-1])/($xdata[$i]-$xdata[$i-1]);
  53          $delta[$i] = (6.0*$delta[$i]/($xdata[$i+1]-$xdata[$i-1])-$s*$delta[$i-1])/$p;
  54      }
  55  
  56      // Backward substitution
  57      for( $j=$n-2; $j >= 0; --$j ) {
  58          $this->y2[$j] = $this->y2[$j]*$this->y2[$j+1] + $delta[$j];
  59      }
  60      }
  61  
  62      // Return the two new data vectors
  63      function Get($num=50) {
  64      $n = $this->n ;
  65      $step = ($this->xdata[$n-1]-$this->xdata[0]) / ($num-1);
  66      $xnew=array();
  67      $ynew=array();
  68      $xnew[0] = $this->xdata[0];
  69      $ynew[0] = $this->ydata[0];
  70      for( $j=1; $j < $num; ++$j ) {
  71          $xnew[$j] = $xnew[0]+$j*$step;
  72          $ynew[$j] = $this->Interpolate($xnew[$j]);
  73      }
  74      return array($xnew,$ynew);
  75      }
  76  
  77      // Return a single interpolated Y-value from an x value
  78      function Interpolate($xpoint) {
  79  
  80      $max = $this->n-1;
  81      $min = 0;
  82  
  83      // Binary search to find interval
  84      while( $max-$min > 1 ) {
  85          $k = ($max+$min) / 2;
  86          if( $this->xdata[$k] > $xpoint ) 
  87          $max=$k;
  88          else 
  89          $min=$k;
  90      }    
  91  
  92      // Each interval is interpolated by a 3:degree polynom function
  93      $h = $this->xdata[$max]-$this->xdata[$min];
  94  
  95      if( $h == 0  ) {
  96          JpGraphError::RaiseL(19002);
  97  //('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.');
  98      }
  99  
 100  
 101      $a = ($this->xdata[$max]-$xpoint)/$h;
 102      $b = ($xpoint-$this->xdata[$min])/$h;
 103      return $a*$this->ydata[$min]+$b*$this->ydata[$max]+
 104           (($a*$a*$a-$a)*$this->y2[$min]+($b*$b*$b-$b)*$this->y2[$max])*($h*$h)/6.0;
 105      }
 106  }
 107  
 108  //------------------------------------------------------------------------
 109  // CLASS Bezier
 110  // Create a new data array from a number of control points
 111  //------------------------------------------------------------------------
 112  class Bezier {
 113  /**
 114   * @author Thomas Despoix, openXtrem company
 115   * @license released under QPL
 116   * @abstract Bezier interoplated point generation,
 117   * computed from control points data sets, based on Paul Bourke algorithm :
 118   * http://astronomy.swin.edu.au/~pbourke/curves/bezier/
 119   */
 120      private $datax = array();
 121      private $datay = array();
 122      private $n=0;
 123   
 124      function Bezier($datax, $datay, $attraction_factor = 1) {
 125      // Adding control point multiple time will raise their attraction power over the curve
 126      $this->n = count($datax);
 127      if( $this->n !== count($datay) ) {
 128          JpGraphError::RaiseL(19003);
 129  //('Bezier: Number of X and Y coordinates must be the same');
 130      }
 131      $idx=0;
 132      foreach($datax as $datumx) {
 133          for ($i = 0; $i < $attraction_factor; $i++) {
 134          $this->datax[$idx++] = $datumx;
 135          }
 136      }
 137         $idx=0;
 138      foreach($datay as $datumy) {
 139          for ($i = 0; $i < $attraction_factor; $i++) {
 140          $this->datay[$idx++] = $datumy;
 141          }
 142      }
 143      $this->n *= $attraction_factor;
 144      }
 145  
 146      function Get($steps) {
 147      $datax = array();
 148      $datay = array();
 149      for ($i = 0; $i < $steps; $i++) {
 150          list($datumx, $datumy) = $this->GetPoint((double) $i / (double) $steps);       
 151          $datax[] = $datumx;
 152          $datay[] = $datumy;
 153      }
 154     
 155      $datax[] = end($this->datax);
 156      $datay[] = end($this->datay);
 157     
 158      return array($datax, $datay);
 159      }
 160   
 161      function GetPoint($mu) {
 162      $n = $this->n - 1;
 163      $k = 0;
 164      $kn = 0;
 165      $nn = 0;
 166      $nkn = 0;
 167      $blend = 0.0;
 168      $newx = 0.0;
 169      $newy = 0.0;
 170  
 171      $muk = 1.0;
 172      $munk = (double) pow(1-$mu,(double) $n);
 173  
 174      for ($k = 0; $k <= $n; $k++) {
 175          $nn = $n;
 176          $kn = $k;
 177          $nkn = $n - $k;
 178          $blend = $muk * $munk;
 179          $muk *= $mu;
 180          $munk /= (1-$mu);
 181          while ($nn >= 1) {
 182          $blend *= $nn;
 183          $nn--;
 184          if ($kn > 1) {
 185              $blend /= (double) $kn;
 186              $kn--;
 187          }
 188          if ($nkn > 1) {
 189              $blend /= (double) $nkn;
 190              $nkn--;
 191          }
 192          }
 193          $newx += $this->datax[$k] * $blend;
 194          $newy += $this->datay[$k] * $blend;
 195      }
 196  
 197      return array($newx, $newy);
 198      }
 199  }
 200  
 201  // EOF
 202  ?>