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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
|
#!/usr/bin/env python3
# -*- python -*-
# Copyright © 2014, 2015, 2016, 2017 Mattias Andrée (maandree@kth.se)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Test of monotone cubic interpolation
# using the Fritsch–Carlson method.
# Does not overshoot, but regular
# cubic interpolation with linear
# replacement for overshots is better.
# Load matplotlib.pyplot,
# it can take some time so
# print information about it.
print('Loading matplotlib.pyplot...')
import matplotlib.pyplot as plot
print('Done loading matplotlib.pyplot')
# Modules used for input data generation
from math import *
from random import *
def main():
# Create a page with graphs
fig = plot.figure()
# Add graphs
add_graph(fig, 221, [i / 15 for i in range(16)])
add_graph(fig, 222, [sin(6 * i / 15) for i in range(16)])
add_graph(fig, 223, [(i / 15) ** 0.5 for i in range(16)])
#add_graph(fig, 221, [random() for i in range(16)])
#add_graph(fig, 222, [random() for i in range(16)])
#add_graph(fig, 223, [random() for i in range(16)])
#add_graph(fig, 224, [random() for i in range(16)])
#add_graph(fig, 224, [(-1) ** (i // 2) for i in range(16)])
add_graph(fig, 224, [i / 15 * (-1) ** (i // 2) for i in range(16)])
# Show graphs
plot.show()
def add_graph(fig, graph_pos, input_values):
'''
Add a graph
@param fig:Figure The page to which to add the graph
@param graph_pos:int Where to place the graph
@param input_values:list<float> The input values for each point
'''
# Interpolate data
output_values = interpolate(input_values)
# Number of input points
n = len(input_values)
# Number of output points
m = len(output_values)
# Create graph
graph = fig.add_subplot(graph_pos)
# Plot interpolated data
graph.plot([i / (m - 1) for i in range(m)], output_values, 'b-')
# Plot input data
graph.plot([i / (n - 1) for i in range(n)], input_values, 'ro')
def interpolate(small, tension = 0):
'''
Interpolate data
@param small:list<float> The input values for each point
@param tension:float A [0, 1] value of the interpolation tension
@return :list<float> The values for each point in a scaled up version
'''
large = [None] * len(small) ** 2
small_, large_ = len(small) - 1, len(large) - 1
# Basis functions
#h00 = lambda t : (1 + 2 * t) * (1 - t) ** 2
h10 = lambda t : t * (1 - t) ** 2
h01 = lambda t : t ** 2 * (3 - 2 * t)
h11 = lambda t : t ** 2 * (t - 1)
# Tension coefficent
c_ = 1 - tension
## Interpolant selection
# Compute the slopes of the secant
# lines between successive points
ds = [small[i + 1] - small[i] for i in range(small_)]
# Initialize the tangents at every
# data point as the average of the secants
ms = [ds[0]] + [(ds[i - 1] + ds[i]) / 2 for i in range(1, small_)] + [ds[small_ - 1]]
βlast = 0
for i in range(small_):
if ds[i] == 0:
# Two successive values are equal, ms[i],
# must be zero to preserve monotonicity,
# no idea to do further work on them.
ms[i], βlast = 0, -1
continue
# Look for local extremums
α, β = ms[i] / ds[i], ms[i + 1] / ds[i]
if (α < 0) or (βlast < 0):
# Local extremum found,
# ensure piecewise monotonicity
ms[i], β = 0, -1
elif α ** 2 + β ** 2 > 9:
# Otherwise, prevent overshoot and ensure
# monotonicity by restricting the (α, β)
# vector to a circle of radius 3.
τ = 3 / (α ** 2 + β ** 2) ** 0.5
ms[i], ms[i + 1] = τ * α * ds[i], τ * β * ds[i]
βlast = β
## Interpolation
for i in range(len(large)):
# Scaling
j = i * small_ / large_
# Floor, weight, ceiling
j, w, k = int(j), j % 1, min(int(j) + 1, small_)
# Points
pj, pk = small[j], small[k]
# Tangents
mj, mk = c_ * ms[j], c_ * ms[k]
# Interpolation
large[i] = pj + h10(w) * mj + h01(w) * (pk - pj) + h11(w) * mk
return large
# Plot interpolation
main()
|
