forked from pa-m/sklearn
-
Notifications
You must be signed in to change notification settings - Fork 0
/
interpolate.go
92 lines (84 loc) · 2.46 KB
/
interpolate.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
package interpolate
import (
"math"
"sort"
)
type xy struct{ x, y []float64 }
func (s *xy) Len() int { return len(s.x) }
func (s *xy) Less(i, j int) bool { return s.x[i] < s.x[j] }
func (s *xy) Swap(i, j int) {
s.x[i], s.x[j] = s.x[j], s.x[i]
s.y[i], s.y[j] = s.y[j], s.y[i]
}
func (s *xy) XY(i int) (float64, float64) { return s.x[i], s.y[i] }
func interpolate2points(x0, y0, x1, y1 float64) func(float64) float64 {
return func(x float64) float64 {
if x1 == x0 {
return (y0 + y1) / 2.
}
return y0 + (x-x0)/(x1-x0)*(y1-y0)
}
}
// Interp1d return linear interpolation from x,y points
// mimics partly scipy.interpolate.interp1d
func Interp1d(x, y []float64) func(x float64) float64 {
var both *xy
if sort.Float64sAreSorted(x) {
both = &xy{x, y}
} else {
both = &xy{make([]float64, len(x)), make([]float64, len(x))}
copy(both.x, x)
copy(both.y, y)
sort.Sort(both)
}
return func(x float64) float64 {
for ix := range both.x[:len(both.x)-1] {
if x < both.x[ix+1] || ix == len(both.x)-2 {
return interpolate2points(both.x[ix], both.y[ix], both.x[ix+1], both.y[ix+1])(x)
}
}
return math.NaN()
}
}
type xyz struct{ x, y, z []float64 }
func (s *xyz) Len() int { return len(s.x) }
func (s *xyz) Less(i, j int) bool { return s.x[i] < s.x[j] || (s.x[i] == s.x[j] && s.y[i] < s.y[j]) }
func (s *xyz) Swap(i, j int) {
s.x[i], s.x[j] = s.x[j], s.x[i]
s.y[i], s.y[j] = s.y[j], s.y[i]
s.z[i], s.z[j] = s.z[j], s.z[i]
}
// Interp2d computes bilinear interpolation from x,y to z
func Interp2d(x, y, z []float64) func(x, y float64) float64 {
xsorted, ysorted, zsorted := make([]float64, len(x)), make([]float64, len(x)), make([]float64, len(x))
copy(xsorted, x)
copy(ysorted, y)
copy(zsorted, z)
all := &xyz{x: xsorted, y: ysorted, z: zsorted}
sort.Sort(all)
return func(x, y float64) float64 {
ix0 := 0
l := len(xsorted)
ix2 := l - 1
for i := range xsorted {
if xsorted[i] > xsorted[ix0] && xsorted[i] <= x {
ix0 = i
}
if xsorted[l-1-i] <= xsorted[ix2] && xsorted[l-1-i] >= x {
ix2 = l - 1 - i
}
}
ix1 := ix0
for ix1 < l && xsorted[ix1] == xsorted[ix0] {
ix1++
}
ix3 := ix2
for ix3 < l && xsorted[ix3] == xsorted[ix2] {
ix3++
}
zxlow := Interp1d(ysorted[ix0:ix1], zsorted[ix0:ix1])(y)
zxhi := Interp1d(ysorted[ix2:ix3], zsorted[ix2:ix3])(y)
//fmt.Println("xlow", xsorted[ix0], zxlow, "xhi", xsorted[ix2], zxhi)
return Interp1d([]float64{xsorted[ix0], xsorted[ix2]}, []float64{zxlow, zxhi})(x)
}
}