#!/usr/bin/env python3
# -*- python -*-

# Copyright © 2014  Mattias Andrée (maandree@member.fsf.org)
# 
# 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()