1
|
|
|
<?php |
2
|
|
|
|
3
|
|
|
/** |
4
|
|
|
* JPGraph v4.0.3 |
5
|
|
|
*/ |
6
|
|
|
|
7
|
|
|
namespace Amenadiel\JpGraph\Util; |
8
|
|
|
|
9
|
|
|
/** |
10
|
|
|
* File: JPGRAPH_REGSTAT.PHP |
11
|
|
|
* // Description: Regression and statistical analysis helper classes |
12
|
|
|
* // Created: 2002-12-01 |
13
|
|
|
* // Ver: $Id: jpgraph_regstat.php 1131 2009-03-11 20:08:24Z ljp $ |
14
|
|
|
* // |
15
|
|
|
* // Copyright (c) Asial Corporation. All rights reserved. |
16
|
|
|
*/ |
17
|
|
|
|
18
|
|
|
/** |
19
|
|
|
* @class Spline |
20
|
|
|
* // Create a new data array from an existing data array but with more points. |
21
|
|
|
* // The new points are interpolated using a cubic spline algorithm |
22
|
|
|
*/ |
23
|
|
|
class Spline |
24
|
|
|
{ |
25
|
|
|
// 3:rd degree polynom approximation |
26
|
|
|
|
27
|
|
|
private $xdata; |
28
|
|
|
private $ydata; // Data vectors |
29
|
|
|
private $y2; // 2:nd derivate of ydata |
30
|
|
|
private $n = 0; |
31
|
|
|
|
32
|
|
|
public function __construct($xdata, $ydata) |
33
|
|
|
{ |
34
|
|
|
$this->y2 = []; |
35
|
|
|
$this->xdata = $xdata; |
36
|
|
|
$this->ydata = $ydata; |
37
|
|
|
|
38
|
|
|
$n = safe_count($ydata); |
39
|
|
|
$this->n = $n; |
40
|
|
|
if ($this->n !== safe_count($xdata)) { |
41
|
|
|
JpGraphError::RaiseL(19001); |
42
|
|
|
//('Spline: Number of X and Y coordinates must be the same'); |
43
|
|
|
} |
44
|
|
|
|
45
|
|
|
// Natural spline 2:derivate == 0 at endpoints |
46
|
|
|
$this->y2[0] = 0.0; |
47
|
|
|
$this->y2[$n - 1] = 0.0; |
48
|
|
|
$delta[0] = 0.0; |
|
|
|
|
49
|
|
|
|
50
|
|
|
// Calculate 2:nd derivate |
51
|
|
|
for ($i = 1; $i < $n - 1; ++$i) { |
52
|
|
|
$d = ($xdata[$i + 1] - $xdata[$i - 1]); |
53
|
|
|
if ($d == 0) { |
54
|
|
|
JpGraphError::RaiseL(19002); |
55
|
|
|
//('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.'); |
56
|
|
|
} |
57
|
|
|
$s = ($xdata[$i] - $xdata[$i - 1]) / $d; |
58
|
|
|
$p = $s * $this->y2[$i - 1] + 2.0; |
59
|
|
|
$this->y2[$i] = ($s - 1.0) / $p; |
60
|
|
|
$delta[$i] = ($ydata[$i + 1] - $ydata[$i]) / ($xdata[$i + 1] - $xdata[$i]) - |
61
|
|
|
($ydata[$i] - $ydata[$i - 1]) / ($xdata[$i] - $xdata[$i - 1]); |
62
|
|
|
$delta[$i] = (6.0 * $delta[$i] / ($xdata[$i + 1] - $xdata[$i - 1]) - $s * $delta[$i - 1]) / $p; |
63
|
|
|
} |
64
|
|
|
|
65
|
|
|
// Backward substitution |
66
|
|
|
for ($j = $n - 2; $j >= 0; --$j) { |
67
|
|
|
$this->y2[$j] = $this->y2[$j] * $this->y2[$j + 1] + $delta[$j]; |
68
|
|
|
} |
69
|
|
|
} |
70
|
|
|
|
71
|
|
|
// Return the two new data vectors |
72
|
|
|
public function Get($num = 50) |
73
|
|
|
{ |
74
|
|
|
$n = $this->n; |
75
|
|
|
$step = ($this->xdata[$n - 1] - $this->xdata[0]) / ($num - 1); |
76
|
|
|
$xnew = []; |
77
|
|
|
$ynew = []; |
78
|
|
|
$xnew[0] = $this->xdata[0]; |
79
|
|
|
$ynew[0] = $this->ydata[0]; |
80
|
|
|
for ($j = 1; $j < $num; ++$j) { |
81
|
|
|
$xnew[$j] = $xnew[0] + $j * $step; |
82
|
|
|
$ynew[$j] = $this->Interpolate($xnew[$j]); |
83
|
|
|
} |
84
|
|
|
|
85
|
|
|
return [$xnew, $ynew]; |
86
|
|
|
} |
87
|
|
|
|
88
|
|
|
// Return a single interpolated Y-value from an x value |
89
|
|
|
public function Interpolate($xpoint) |
90
|
|
|
{ |
91
|
|
|
$max = $this->n - 1; |
92
|
|
|
$min = 0; |
93
|
|
|
|
94
|
|
|
// Binary search to find interval |
95
|
|
|
while ($max - $min > 1) { |
96
|
|
|
$k = ($max + $min) / 2; |
97
|
|
|
if ($this->xdata[$k] > $xpoint) { |
98
|
|
|
$max = $k; |
99
|
|
|
} else { |
100
|
|
|
$min = $k; |
101
|
|
|
} |
102
|
|
|
} |
103
|
|
|
|
104
|
|
|
// Each interval is interpolated by a 3:degree polynom function |
105
|
|
|
$h = $this->xdata[$max] - $this->xdata[$min]; |
106
|
|
|
|
107
|
|
|
if ($h == 0) { |
108
|
|
|
JpGraphError::RaiseL(19002); |
109
|
|
|
//('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.'); |
110
|
|
|
} |
111
|
|
|
|
112
|
|
|
$a = ($this->xdata[$max] - $xpoint) / $h; |
113
|
|
|
$b = ($xpoint - $this->xdata[$min]) / $h; |
114
|
|
|
|
115
|
|
|
return $a * $this->ydata[$min] + $b * $this->ydata[$max] + |
116
|
|
|
(($a * $a * $a - $a) * $this->y2[$min] + ($b * $b * $b - $b) * $this->y2[$max]) * ($h * $h) / 6.0; |
117
|
|
|
} |
118
|
|
|
} |
119
|
|
|
|
120
|
|
|
// EOF |
121
|
|
|
|