e2pdf / vendors / svggraph / BestFitCurve.php

BestFitCurve.php in E2Pdf – Export Pdf Tool for WordPress 1.32.26, at vendors/svggraph/BestFitCurve.php

246 lines 6.0 KB

Raw Download Zip

1	<?php
2	/**
3	* Copyright (C) 2023 Graham Breach
4	*
5	* This program is free software: you can redistribute it and/or modify
6	* it under the terms of the GNU Lesser General Public License as published by
7	* the Free Software Foundation, either version 3 of the License, or
8	* (at your option) any later version.
9	*
10	* This program is distributed in the hope that it will be useful,
11	* but WITHOUT ANY WARRANTY; without even the implied warranty of
12	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13	* GNU Lesser General Public License for more details.
14	*
15	* You should have received a copy of the GNU Lesser General Public License
16	* along with this program. If not, see <http://www.gnu.org/licenses/>.
17	*/
18	/**
19	* For more information, please contact <graham@goat1000.com>
20	*/
21
22	namespace Goat1000\SVGGraph;
23
24	/**
25	* Class for calculating a curved best-fit line
26	*/
27	class BestFitCurve {
28
29	protected $graph;
30	protected $points;
31	protected $line;
32	protected $projection;
33	protected $types;
34
35	public function __construct(&$graph, $points, $types = null)
36	{
37	$this->graph =& $graph;
38	$this->line = new PathData;
39	$this->projection = new PathData;
40	$this->points = $points;
41	$this->types = $types === null \|\| is_array($types) ? $types : [$types];
42	}
43
44	/**
45	* Calculates the line and projection
46	*/
47	public function calculate(BoundingBox $area, $limit_start, $limit_end,
48	$project_start, $project_end)
49	{
50	// can't draw a line through fewer than 2 points
51	$count = count($this->points);
52	if($count < 2)
53	return false;
54
55	$old_scale = bcscale(50);
56	$b = [];
57	$y = new Matrix($count, 1);
58	foreach($this->points as $p)
59	$b[] = $p->y;
60	$y->load($b);
61
62	$supported_types = ['straight', 'quadratic', 'cubic', 'quartic', 'quintic'];
63
64	// choose which functions to fit
65	if($this->types !== null) {
66	$types = [];
67	foreach($this->types as $type) {
68	if(!in_array($type, $supported_types))
69	throw new \Exception("Unknown curve type '{$type}'");
70	$types[] = $type;
71	}
72	} else {
73	$types = ['quintic'];
74	switch($count)
75	{
76	case 2 : $types = ['straight'];
77	break;
78	case 3 : $types = ['quadratic'];
79	break;
80	case 4 : $types = ['cubic'];
81	break;
82	case 5 : $types = ['quartic'];
83	}
84	}
85
86	// fit the functions, measure the error
87	$results = [];
88	$errors = [];
89	foreach($types as $t) {
90	$v = $this->vandermonde($t);
91	$result = $this->solve($v, $y);
92	if($result !== null) {
93	$errors[$t] = $this->error($v, $result);
94	$results[$t] = $result;
95	}
96	}
97
98	if(!empty($errors)) {
99
100	// sort by error, best first
101	uasort($errors, 'bccomp');
102	$best = null;
103	foreach($errors as $k => $v) {
104	$best = $k;
105	break;
106	}
107
108	$r = $results[$best];
109	$c = $r->asArray();
110
111	// plot the function
112	$fn = new Algebraic($best);
113	$fn->setCoefficients($c);
114	$this->buildPaths($fn, $area, $limit_start, $limit_end,
115	$project_start, $project_end);
116	}
117	bcscale($old_scale);
118	}
119
120	/**
121	* Calculates the error
122	*/
123	protected function error(Matrix $v, Matrix $r)
124	{
125	$old_scale = bcscale(50);
126	$tr = $r->transpose();
127	$vr = $v->multiply($tr);
128	$i = 0;
129	$err = "0";
130	foreach($this->points as $p) {
131	$diff = bcsub($vr($i, 0), $p->y);
132	if(bccomp($diff, "0") == -1)
133	$err = bcsub($err, $diff);
134	else
135	$err = bcadd($err, $diff);
136	++$i;
137	}
138	bcscale($old_scale);
139	return $err;
140	}
141
142	/**
143	* Solves the normal equation
144	*/
145	protected function solve(Matrix $v, Matrix $y)
146	{
147	$v_t = $v->transpose();
148	$v_t_v = $v_t->multiply($v);
149	$v_t_y = $v_t->multiply($y);
150	return $v_t_v->gaussian_solve($v_t_y);
151	}
152
153	/**
154	* Returns the Vandermonde matrix for the curve type
155	*/
156	protected function vandermonde($type)
157	{
158	$old_scale = bcscale(50);
159
160	// find size of matrix
161	$a = new Algebraic($type);
162	$test = $a->vandermonde(1);
163	$m = count($this->points);
164	$n = count($test) + 1;
165	$v = new Matrix($m, $n);
166
167	$i = 0;
168	foreach($this->points as $p)
169	{
170	$v($i, 0, 1);
171	$bcx = sprintf("%20.20F", $p->x);
172	$cols = $a->vandermonde($bcx);
173	foreach($cols as $k => $value)
174	$v($i, $k + 1, $value);
175	++$i;
176	}
177	bcscale($old_scale);
178	return $v;
179	}
180
181	/**
182	* Builds the line and projection paths.
183	* For vertical lines, $slope = null and $y_int = $x
184	*/
185	protected function buildPaths(&$fn, $area, $limit_start, $limit_end,
186	$project_start, $project_end)
187	{
188
189	// initialize min and max points of line
190	$x_min = $limit_start === null ? 0 : max($limit_start, 0);
191	$x_max = $limit_end === null ? $area->width() : min($limit_end, $area->width());
192	$y_min = 0;
193	$y_max = $area->height();
194	$line = new PathData;
195	$projection = new PathData;
196
197	$step = 1;
198	if($project_start)
199	$this->buildPath($projection, $fn, 0, $x_min, $y_min, $y_max, $step, $area);
200	$this->buildPath($line, $fn, $x_min, $x_max, $y_min, $y_max, $step, $area);
201	if($project_end)
202	$this->buildPath($projection, $fn, $x_max, $area->width(), $y_min, $y_max, $step, $area);
203
204	$this->projection = $projection;
205	$this->line = $line;
206	return true;
207	}
208
209	/**
210	* Builds a single path section between $x1 and $x2
211	*/
212	private function buildPath(&$path, &$fn, $x1, $x2, $y_min, $y_max, $step, $area)
213	{
214	$cmd = 'M';
215	for($x = $x1; $x <= $x2; $x += $step) {
216	$y = $fn($x);
217	if($y < $y_min \|\| $y > $y_max) {
218	$cmd = 'M';
219	continue;
220	}
221	$path->add($cmd, $area->x1 + $x, $area->y2 - $y);
222	switch($cmd) {
223	case 'M' : $cmd = 'L'; break;
224	case 'L' : $cmd = ''; break;
225	}
226	}
227	}
228
229	/**
230	* Returns the best-fit line as PathData
231	*/
232	public function getLine()
233	{
234	return $this->line;
235	}
236
237	/**
238	* Returns the projection line(s) as PathData
239	*/
240	public function getProjection()
241	{
242	return $this->projection;
243	}
244	}
245
246