ssa

1
package pchip
2

3
import "github.com/ready-steady/linear/system"
4

5
// Cubic is a cubic-spline interpolant.
6
type Cubic struct {
7
	Nodes   []float64
8
	Weights []float64
9
}
10

11
// NewCubic constructs a cubic-spline interpolant for a function y = f(x) given
12
// as a series of points (x, y). The x coordinates should be a strictly
13
// increasing sequence with at least two elements. The corresponding y
14
// coordinates can be multidimensional.
15
func NewCubic(x, y []float64) *Cubic {
16
	nn := len(x)
17
	nd := len(y) / nn
18

19
	dx := make([]float64, nn-1)
20
	dydx := make([]float64, (nn-1)*nd)
21
	for i := 0; i < (nn - 1); i++ {
22
		dx[i] = x[i+1] - x[i]
23
		for j := 0; j < nd; j++ {
24
			dydx[i*nd+j] = (y[(i+1)*nd+j] - y[i*nd+j]) / dx[i]
25
		}
26
	}
27

28
	s := &Cubic{}
29

30
	switch nn {
31
	case 2:
32
		s.Nodes = []float64{x[0], x[1]}
33
		s.Weights = make([]float64, nd*4)
34
		for j := 0; j < nd; j++ {
35
			s.Weights[j*4+2] = dydx[j]
36
			s.Weights[j*4+3] = y[j]
37
		}
38
	case 3:
39
		s.Nodes = []float64{x[0], x[2]}
40
		s.Weights = make([]float64, nd*4)
41
		for j := 0; j < nd; j++ {
42
			c1 := (dydx[nd+j] - dydx[j]) / (x[2] - x[0])
43
			s.Weights[j*4+1] = c1
44
			s.Weights[j*4+2] = dydx[j] - c1*dx[0]
45
			s.Weights[j*4+3] = y[j]
46
		}
47
	default:
48
		xb := x[2] - x[0]
49
		xe := x[nn-1] - x[nn-3]
50

51
		a := make([]float64, nn)
52
		for i := 0; i < (nn - 2); i++ {
53
			a[i] = dx[i+1]
54
		}
55
		a[nn-2] = xe
56

57
		b := make([]float64, nn)
58
		b[0] = dx[1]
59
		for i := 1; i < (nn - 1); i++ {
60
			b[i] = 2 * (dx[i] + dx[i-1])
61
		}
62
		b[nn-1] = dx[nn-3]
63

64
		c := make([]float64, nn)
65
		c[1] = xb
66
		for i := 2; i < nn; i++ {
67
			c[i] = dx[i-2]
68
		}
69

70
		d := make([]float64, nd*nn)
71
		for j := 0; j < nd; j++ {
72
			d[j*nn] = ((dx[0]+2*xb)*dx[1]*dydx[j] + dx[0]*dx[0]*dydx[nd+j]) / xb
73
			for i := 1; i < (nn - 1); i++ {
74
				d[j*nn+i] = 3 * (dx[i]*dydx[(i-1)*nd+j] + dx[i-1]*dydx[i*nd+j])
75
			}
76
			d[j*nn+nn-1] = (dx[nn-2]*dx[nn-2]*dydx[(nn-3)*nd+j] +
77
				(2*xe+dx[nn-2])*dx[nn-3]*dydx[(nn-2)*nd+j]) / xe
78
		}
79

80
		slopes := system.ComputeTridiagonal(a, b, c, d)
81

82
		s.Nodes = make([]float64, nn)
83
		copy(s.Nodes, x)
84

85
		s.Weights = make([]float64, (nn-1)*nd*4)
86
		for i, k := 0, 0; i < (nn - 1); i++ {
87
			for j := 0; j < nd; j++ {
88
				α := (dydx[i*nd+j] - slopes[j*nn+i]) / dx[i]
89
				β := (slopes[j*nn+i+1] - dydx[i*nd+j]) / dx[i]
90
				s.Weights[k] = (β - α) / dx[i]
91
				s.Weights[k+1] = 2*α - β
92
				s.Weights[k+2] = slopes[j*nn+i]
93
				s.Weights[k+3] = y[i*nd+j]
94
				k += 4
95
			}
96
		}
97
	}
98

99
	return s
100
}
101

102
// Evaluate interpolates the function values y = f(x). The x coordinates should
103
// be an increasing sequence whose elements belong to the range of the points
104
// that the interpolant has been constructed with.
105
func (s *Cubic) Evaluate(x []float64) []float64 {
106
	nn, np := len(s.Nodes), len(x)
107
	nd := len(s.Weights) / (4 * (nn - 1))
108

109
	y := make([]float64, np*nd)
110

111
	for i, l := 0, 0; i < np; i++ {
112
		for x[i] > s.Nodes[l+1] && l < (nn-2) {
113
			l++
114
		}
115

116
		z := x[i] - s.Nodes[l]
117
		for j, k := 0, l*nd*4; j < nd; j++ {
118
			y[i*nd+j] = z*(z*(z*s.Weights[k]+s.Weights[k+1])+s.Weights[k+2]) + s.Weights[k+3]
119
			k += 4
120
		}
121
	}
122

123
	return y
124
}
125