# Source code for PySpice.Math.Calculus

####################################################################################################
#
# PySpice - A Spice Package for Python
# Copyright (C) 2014 Fabrice Salvaire
#
# 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/>.
#
####################################################################################################

"""This module provides algorithms to compute the derivative of a function sampled on an uniform
grid.
"""

####################################################################################################

import fractions

import numpy as np

####################################################################################################

from PySpice.Math import odd

####################################################################################################

[docs]def compute_exact_finite_difference_coefficients(derivative_order, grid, x0=0):

"""This function compute the finite difference coefficients for the given derivative order and
grid.  The parameter *x* specifies where is computed the derivative on the grid.  The grid is
given as a list of integer offsets.

The algorithm is derived from the article: Generation of Finite Difference Formulas on Arbitrary Space
Grids, Bengt Fornberg, Mathematics of computation, volume 51, number 184, october 1988
"""

N = len(grid)

# d[m,n,v]
d = [[[0
for v in range(N)]
for n in range(N)]
for m in range(derivative_order +1)]

d[0][0][0] = fractions.Fraction(1,1)
c1 = 1
for n in range(1, N):
c2 = 1
for v in range(n):
c3 = grid[n] - grid[v]
c2 *= c3
if n <= derivative_order:
d[n][n-1][v] = 0
for m in range(min(n, derivative_order) +1):
d[m][n][v] = ( (grid[n] - x0)*d[m][n-1][v] - m*d[m-1][n-1][v] ) / c3
for m in range(min(n, derivative_order) +1):
d[m][n][n] = fractions.Fraction(c1,c2)*( m*d[m-1][n-1][n-1] - (grid[n-1] - x0)*d[m][n-1][n-1] )
c1 = c2

return d[-1][-1]

####################################################################################################

def compute_finite_difference_coefficients(derivative_order, grid):
return [float(x) for x in compute_exact_finite_difference_coefficients(derivative_order, grid)]

####################################################################################################

_coefficient_cache = dict(centred={}, forward={}, backward={})

def get_finite_difference_coefficients(derivative_order, accuracy_order, grid_type):

if derivative_order < 1:
raise ValueError("Wrong derivative order")

if odd(accuracy_order) or accuracy_order < 2:
raise ValueError("Wrong accuracy order")

if grid_type == 'centred':
window_size = accuracy_order // 2
grid = list(range(-window_size, window_size +1))
elif grid_type == 'forward':
grid = list(range(derivative_order + accuracy_order))
elif grid_type == 'backward':
grid = list(range(-(derivative_order + accuracy_order) +1, 1))
grid = list(reversed(grid)) # Fixme: why ?
else:
raise ValueError("Wrong grid type")

key = '{}-{}'.format(derivative_order, accuracy_order)
coefficients = _coefficient_cache[grid_type].get(key, None)
if coefficients is None:
coefficients = compute_finite_difference_coefficients(derivative_order, grid)
_coefficient_cache[grid_type][key] = coefficients

return grid, coefficients

####################################################################################################

[docs]def simple_derivative(x, values):
""" Compute the derivative as a simple slope. """
return x[:-1], np.diff(values)/np.diff(x)

####################################################################################################

[docs]def derivative(x, values, derivative_order=1, accuracy_order=4):

"""Compute the derivative at the given derivative order and accuracy order. The precision of the
Taylor expansion is :math:\mathcal{O}(dx^{accuracy}).
"""

dx = np.diff(x)
# if not np.all(dx == dx[0]):
#     raise ValueError("Sampling is not uniform")
dx = dx[0]

values_size, = values.shape

derivative = np.zeros(values_size, dtype=values.dtype)

grid, coefficients = get_finite_difference_coefficients(derivative_order, accuracy_order, 'centred')
window_size = grid[-1]
# print grid, coefficients
vector_size = values_size - 2*window_size
if not vector_size:
raise ValueError("The size of the value's array is not sufficient for the given accuracy order")
lower_index = window_size
upper_index = values_size - window_size
derivative_view = derivative[window_size:-window_size]
for offset, coefficient in zip(grid, coefficients):
if coefficient:
# print offset, lower_index + offset, upper_index + offset
derivative_view += values[lower_index + offset:upper_index + offset] * coefficient

grid, coefficients = get_finite_difference_coefficients(derivative_order, accuracy_order, 'forward')
# print grid, coefficients
grid_size = len(grid)
upper_index = window_size
derivative_view = derivative[:window_size]
for offset, coefficient in zip(grid, coefficients):
# print offset, offset, window_size+offset
derivative_view += values[offset:upper_index + offset] * coefficient

grid, coefficients = get_finite_difference_coefficients(derivative_order, accuracy_order, 'backward')
# print grid, coefficients
grid_size = len(grid)
lower_index = values_size - window_size
upper_index = values_size
derivative_view = derivative[-window_size:]
for offset, coefficient in zip(grid, coefficients):
# print offset, lower_index + offset, upper_index + offset
derivative_view += values[lower_index + offset:upper_index + offset] * coefficient

return derivative / dx**derivative_order