Subversion Repositories Applications.papyrus

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2150 mathias 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
?>