Code Coverage
 
Lines
Functions and Methods
Classes and Traits
Total
100.00% covered (success)
100.00%
50 / 50
100.00% covered (success)
100.00%
2 / 2
 CRAP
100.00% covered (success)
100.00%
1 / 1
CubicSplineInterpolation
100.00% covered (success)
100.00%
50 / 50
100.00% covered (success)
100.00%
2 / 2
 10
100.00% covered (success)
100.00%
1 / 1
 __construct
100.00% covered (success)
100.00%
37 / 37
100.00% covered (success)
100.00%
1 / 1
 5
 interpolate
100.00% covered (success)
100.00%
13 / 13
100.00% covered (success)
100.00%
1 / 1
 5
1<?php
2/**
3 * Jingga
4 *
5 * PHP Version 8.1
6 *
7 * @package   phpOMS\Math\Numerics\Interpolation
8 * @copyright Dennis Eichhorn
9 * @license   OMS License 2.0
10 * @version   1.0.0
11 * @link      https://jingga.app
12 */
13declare(strict_types=1);
14
15namespace phpOMS\Math\Numerics\Interpolation;
16
17use phpOMS\Math\Matrix\Matrix;
18use phpOMS\Math\Matrix\Vector;
19
20/**
21 * Cubic spline interpolation.
22 *
23 * @package phpOMS\Math\Numerics\Interpolation
24 * @license OMS License 2.0
25 * @link    https://jingga.app
26 * @since   1.0.0
27 */
28final class CubicSplineInterpolation implements InterpolationInterface
29{
30    /**
31     * Points for spline interpolation
32     *
33     * @var array<array{x:int|float, y:int|float}>
34     * @since 1.0.0
35     */
36    private array $points = [];
37
38    /**
39     * Parameter a of cubic spline
40     *
41     * @var Vector
42     * @since 1.0.0
43     */
44    private Vector $solveA;
45
46    /**
47     * Parameter b of cubic spline
48     *
49     * @var Matrix
50     * @since 1.0.0
51     */
52    private Matrix $solveB;
53
54    /**
55     * Parameter c of cubic spline
56     *
57     * @var Vector
58     * @since 1.0.0
59     */
60    private Vector $solveC;
61
62    /**
63     * Constructor.
64     *
65     * @param array<array{x:int|float, y:int|float}> $points              Points to create the interpolation with
66     * @param float                                  $leftCurvature       Left point curvature
67     * @param int                                    $leftDerivativeType  Derivative type for the left point
68     * @param float                                  $rightCurvature      Right point curvature
69     * @param int                                    $rightDerivativeType Derivative type for the right point
70     *
71     * @since 1.0.0
72     */
73    public function __construct(
74        array $points,
75        float $leftCurvature = 0.0,
76        int $leftDerivativeType = DerivativeType::FIRST,
77        float $rightCurvature = 0.0,
78        int $rightDerivativeType = DerivativeType::FIRST
79    ) {
80        $this->points = $points;
81
82        $n      = \count($this->points);
83        $b      = [];
84        $matrix = new Matrix($n, $n);
85
86        for ($i = 1; $i < $n - 1; ++$i) {
87            $matrix->set($i, $i - 1, 1.0 / 3.0 * ($this->points[$i]['x'] - $this->points[$i - 1]['x']));
88            $matrix->set($i, $i, 2.0 / 3.0 * ($this->points[$i + 1]['x'] - $this->points[$i - 1]['x']));
89            $matrix->set($i, $i + 1, 1.0 / 3.0 * ($this->points[$i + 1]['x'] - $this->points[$i]['x']));
90
91            $b[$i] = ($this->points[$i + 1]['y'] - $this->points[$i]['y']) / ($this->points[$i + 1]['x'] - $this->points[$i]['x'])
92                - ($this->points[$i]['y'] - $this->points[$i - 1]['y']) / ($this->points[$i]['x'] - $this->points[$i - 1]['x']);
93        }
94
95        if ($leftDerivativeType === DerivativeType::FIRST) {
96            $matrix->set(0, 0, 2.0 * ($this->points[1]['x'] - $this->points[0]['x']));
97            $matrix->set(0, 1, 1.0 * ($this->points[1]['x'] - $this->points[0]['x']));
98
99            $b[0] = 3.0 * (($this->points[1]['y'] - $this->points[0]['y']) / ($this->points[1]['x'] - $this->points[0]['x']) - $rightCurvature);
100        } else {
101            $matrix->set(0, 0, 2.0);
102            $matrix->set(0, 1, 0.0);
103
104            $b[0] = $leftCurvature;
105        }
106
107        if ($rightDerivativeType === DerivativeType::FIRST) {
108            $matrix->set($n - 1, $n - 1, 2.0 * ($this->points[$n - 1]['x'] - $this->points[$n - 2]['x']));
109            $matrix->set($n - 1, $n - 2, 1.0 * ($this->points[$n - 1]['x'] - $this->points[$n - 2]['x']));
110
111            $b[$n - 1] = 3.0 * ($rightCurvature - ($this->points[$n - 1]['y'] - $this->points[$n - 2]['y']) / ($this->points[$n - 1]['x'] - $this->points[$n - 2]['x']));
112        } else {
113            $matrix->set($n - 1, $n - 1, 2.0);
114            $matrix->set($n - 1, $n - 2, 0.0);
115
116            $b[$n - 1] = $rightCurvature;
117        }
118
119        $bVector = new Vector($n);
120        $bVector->setMatrixV($b);
121
122        $this->solveB = $matrix->solve($bVector);
123        $this->solveA = new Vector($n);
124        $this->solveC = new Vector($n);
125
126        for ($i = 0; $i < $n - 1; ++$i) {
127            $this->solveA->setV($i, 1.0 / 3.0 * ($this->solveB->get($i + 1) - $this->solveB->get($i)) / ($this->points[$i + 1]['x'] - $this->points[$i]['x']));
128            $this->solveC->setV($i,
129                ($this->points[$i + 1]['y'] - $this->points[$i]['y']) / ($this->points[$i + 1]['x'] - $this->points[$i]['x'])
130                - 1.0 / 3.0 * (2 * $this->solveB->get($i) + $this->solveB->get($i + 1)) * ($this->points[$i + 1]['x'] - $this->points[$i]['x']));
131        }
132
133        $h = $this->points[$n - 1]['x'] - $this->points[$n - 2]['x'];
134
135        $this->solveA->setV($n - 1, 0.0);
136        $this->solveC->setV($n - 1, 3.0 * $this->solveA->getV($n - 2) * $h ** 2 + 2.0 * $this->solveB->get($n - 2) * $h + $this->solveC->getV($n - 2));
137
138        /**
139         * linear extrapolation at start and end point
140         * $this->solveB->setV($n - 1, 0.0)
141         */
142    }
143
144    /**
145     * {@inheritdoc}
146     */
147    public function interpolate(int | float $x) : float
148    {
149        $n    = \count($this->points);
150        $xPos = $n - 1;
151
152        foreach ($this->points as $key => $point) {
153            if ($x <= $point['x']) {
154                $xPos = $key;
155
156                break;
157            }
158        }
159
160        $xPos = \max($xPos - 1, 0);
161        $h    = $x - $this->points[$xPos]['x'];
162
163        if ($x < $this->points[0]['x']) {
164            return ($this->solveB->get(0) * $h + $this->solveC->getV(0)) * $h + $this->points[0]['y'];
165            /**
166             * linear extrapolation at start and end point
167             * ($this->solveC->getV(0)) * $h + $this->points[0]['y'];
168             */
169        } elseif ($x > $this->points[$n - 1]['x']) {
170            return ($this->solveB->get($n - 1) * $h + $this->solveC->getV($n - 1) * $h + $this->points[$n - 1]['y']);
171        }
172
173        return (($this->solveA->getV($xPos) * $h + $this->solveB->get($xPos)) * $h + $this->solveC->getV($xPos)) * $h + $this->points[$xPos]['y'];
174    }
175}