# (C) Copyright IBM Corp. 2019, 2020, 2021, 2022.
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
from typing import Optional
import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline as ius
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve
"""differentiation.py
This module contains derivative approximation algorithms
commonly used in numerical analysis
"""
# This is like an identifier
[docs]class Derivative(object):
"""Parent class for the derivative operations
It is used for identification purposes
"""
def __init__(self):
pass
# This is like an identifier
[docs]class TimeDerivative(Derivative):
"""Numerical time derivative based on centered differences.
"""
def __init__(self):
super().__init__()
def __call__(self, u:np.ndarray=None, delta:float=1) -> np.ndarray:
return np.gradient(u, delta)
# This is like an identifier
[docs]class SpaceDerivative(Derivative):
def __init__(self):
super().__init__()
[docs]class LeleDerivative:
def __init__(self, N:int=1, h:float=1) -> None:
"""10th-order derivative algorithm purposed by Lele.
:param N: number of mesh nodes
:type param: int
:param h: mesh size
:type h: float
"""
self.RightVector = np.array(N)
# Constructing the left operator
LeftMatrix = lil_matrix((N + 4, N + 4))
RightMatrix = lil_matrix((N + 6, N + 6))
for j in range(2, N + 2):
i = j - 2
delta_1 = self.delta_1(i, N)
delta_2 = self.delta_2(i, N)
LeftMatrix[j, j - 1] = delta_1
LeftMatrix[j, j + 1] = delta_1
LeftMatrix[j, j - 2] = delta_2
LeftMatrix[j, j + 2] = delta_2
LeftMatrix[j, j] = 1.
self.LeftMatrix = LeftMatrix[2:-2, 2:-2]
for j in range(3, N + 3):
i = j - 3
delta_3 = self.delta_3(i, N, h)
delta_4 = self.delta_4(i, N, h)
delta_5 = self.delta_5(i, N, h)
delta_0 = self.delta_0(i, N, h)
RightMatrix[j, j - 1] = -delta_5
RightMatrix[j, j + 1] = delta_5
RightMatrix[j, j - 2] = -delta_4
RightMatrix[j, j + 2] = delta_4
RightMatrix[j, j - 3] = -delta_3
RightMatrix[j, j + 3] = delta_3
RightMatrix[j, j] = delta_0
self.RightMatrix = RightMatrix[3:-3, 3:-3]
[docs] def solve(self, f:np.ndarray) -> np.ndarray:
"""Performing Lele differentiation
:param f: variable values to be differentiated
:type f: np.ndarray
:return: the derivatives of f
:rtype: np.ndarray
"""
b = self.RightMatrix @ f
print("Performing Lele Derivation.")
return spsolve(self.LeftMatrix, b)
[docs] def delta_1(self, j:int, N:int) -> float:
"""term delta_1 from the Lele's expression
:param j: node index
:type j: int
:param N: total number of nodes
:type N: int
:return: the evaluation of delta_1
:rtype: float
"""
if j == 0 or j == N - 1:
return 2
elif j == 1 or j == N - 2:
return 1. / 4
elif j == 2 or j == N - 3:
return 4.7435 / 10.67175
elif j == 3 or j == N - 4:
return 4.63271875 / 9.38146875
else:
return 1. / 2
[docs] def delta_2(self, j:int, N:int) -> float:
"""term delta_2 from the Lele's expression
:param j: node index
:type j: int
:param N: total number of nodes
:type N: int
:return: the evaluation of delta_2
:rtype: float
"""
if j == 0 or j == N - 1:
return 0
elif j == 1 or j == N - 2:
return 0
elif j == 2 or j == N - 3:
return 0.2964375 / 10.67175
elif j == 3 or j == N - 4:
return 0.451390625 / 9.38146875
else:
return 1. / 20
[docs] def delta_3(self, j:int, N:int, h:float) -> float:
"""term delta_3 from the Lele's expression
:param j: node index
:type j: int
:param N: total number of nodes
:type N: int
:param h: mesh discretization size
:type h: float
:return: the evaluation of delta_3
:rtype: float
"""
if j == 0 or j == N - 1:
return 0
elif j == 1 or j == N - 2:
return 0
elif j == 2 or j == N - 3:
return 0
elif j == 3 or j == N - 4:
return 6 * (0.015 / 9.38146875) / (6 * h)
else:
return 1. / 100 * (1 / (6 * h))
[docs] def delta_4(self, j:int, N:int, h:float) -> float:
"""term delta_4 from the Lele's expression
:param j: node index
:type j: int
:param N: total number of nodes
:type N: int
:param h: mesh discretization size
:type h: float
:return: the evaluation of delta_4
:rtype: float
"""
if j == 0:
return 1. / (2 * h)
if j == N - 1:
return -1. / (2 * h)
elif j == 1 or j == N - 2:
return 0
elif j == 2 or j == N - 3:
return 4 * (1.23515625 / 10.67175) / (4 * h)
elif j == 3 or j == N - 4:
return 4 * (1.53 / 9.38146875) / (4 * h)
else:
return (101. / 105) * (1. / (4 * h))
[docs] def delta_5(self, j:int, N:int, h:float) -> float:
"""term delta_5 from the Lele's expression
:param j: node index
:type j: int
:param N: total number of nodes
:type N: int
:param h: mesh discretization size
:type h: float
:return: the evaluation of delta_5
:rtype: float
"""
if j == 0:
return 2. / h
elif j == N - 1:
return -2. / h
elif j == 1 or j == N - 2:
return 3. / 2 * (1. / (2 * h))
elif j == 2 or j == N - 3:
return 2 * (7.905 / 10.67175) / (2 * h)
elif j == 3 or j == N - 4:
return 2 * (6.66984375 / 9.38146875) / (2 * h)
else:
return 17. / 12 * (1. / (2 * h))
[docs] def delta_0(self, j:int, N:int, h:float) -> float:
"""term delta_0 from the Lele's expression
:param j: node index
:type j: int
:param N: total number of nodes
:type N: int
:param h: mesh discretization size
:type h: float
:return: the evaluation of delta_0
:rtype: float
"""
if j == 0:
return -5 / (2 * h)
elif j == N - 1:
return 5 / (2 * h)
else:
return 0
[docs]class CenteredDerivative:
def __init__(self, config:dict=None) -> None:
"""Performing second-order centered differentiation
:param config: Configuration dictionary
:type config: dict
:return: nothing
"""
self.step = config['step']
assert self.step, "A value for the differentiation step must be provided."
[docs] def solve(self, data:np.ndarray=None, axis:int=0) -> np.ndarray:
"""It solves the centered derivative
:param data: the values to be differentiated
:type data: np.ndarray
:param axis: axis to perform the differentiation
:type axis: int
:return: the derivatives of data
:rtype: np.ndarray
"""
print("Performing Second-Order Centered Derivation.")
return np.gradient(data, self.step, axis=axis)
def __call__(self, data:np.ndarray=None) -> np.ndarray:
"""__call__ wrapper for executing self.solve
:param data: the values to be differentiated
:type data: np.ndarray
:return: the derivatives of data
:rtype: np.ndarray
"""
return self.solve(data=data)
# It interpolates the derivatives using splines and derive it
[docs]class CollocationDerivative:
def __init__(self, config:dict=None) -> None:
"""Derivative using splines to interpolate the function to be differentiated
:param config: a dictionary containing some parameters to be passed to ScipY engine
:type config: dict
:return: nothing
"""
self.step = None
self.t = None
if 'step' in config.keys():
self.step = config['step']
self.original_shape = None
def _guaratee_correct_shape(self, data:np.ndarray) -> np.ndarray:
"""This method ensures the input as one-dimensional arrays
since the scipy.interpolate.InterpolatedUnivariateSpline class
requires it.
:param data: The input dataset to be differentiated
:type data: np.ndarray
:return: the input data properly reshaped (if necessary)
:rtype: np.ndarray
"""
if len(data.shape) <= 2:
return data
else:
self.original_shape = data.shape
immutable_dim = data.shape[0]
collapsible_dims = data.shape[1:]
collapsed_shape = np.prod(collapsible_dims)
return data.reshape((immutable_dim, collapsed_shape))
[docs] def solve(self, data:np.ndarray=None, x: Optional[np.ndarray]=None) -> np.ndarray:
"""Performing collocation differentiation
:param data: values to be differentiated
:type data: np.ndarray
:param x: the axis to be used for executing the differentiation (Optional)
:return: the derivatives of the input data along the chosen axis
:rtype: np.ndarray
"""
print("Performing Collocation Derivation.")
intp_list = list()
N = data.shape[0]
if self.step:
self.t = [self.step * ti for ti in range(N)]
elif x is not None:
self.t = x
else:
raise Exception("There is no way for evaluating derivatives.")
data = self._guaratee_correct_shape(data)
for i in range(data.shape[1]):
intp = ius(self.t, data[:, i])
dintp = intp.derivative()
intp_list.append(dintp(self.t)[:, None])
intp_array = np.hstack(intp_list)
return intp_array.reshape(self.original_shape)
[docs] def interpolate_and_solve(self, data:np.ndarray=None,
x_grid:np.ndarray=None, x:np.ndarray=None) -> (np.ndarray, np.ndarray):
"""Performing collocation differentiation
:param data: values to be differentiated
:type data: np.ndarray
:param x_grid: the grid in which the input is defined
:type x_grid: np.ndarray
:param x: the grid in which to interpolate the input
:return: a pair (interpolated, derivated)
:rtype: (np.ndarray, np.ndarray)
"""
assert data.shape[0] <= x.shape[0], "In order to perform interpolation, it is necessary dim(x) > dim(data)."
print("Performing Interpolation and Collocation Derivation.")
intp_list = list()
diff_intp_list = list()
self.t = x_grid
data = self._guaratee_correct_shape(data)
for i in range(data.shape[1]):
intp = ius(self.t, data[:, i])
dintp = intp.derivative()
intp_list.append(intp(x)[:, None])
diff_intp_list.append(dintp(x)[:, None])
intp_array = np.hstack(intp_list)
diff_intp_array = np.hstack(diff_intp_list)
if self.original_shape is not None:
return intp_array.reshape((-1,) + self.original_shape[1:]), \
diff_intp_array.reshape((-1,) + self.original_shape[1:])
else:
return intp_array.reshape(self.original_shape), diff_intp_array.reshape(self.original_shape)
def __call__(self, data:np.ndarray=None) -> np.ndarray:
"""__call__ wrapper for executing self.solve
:param data: the values to be differentiated
:type data: np.ndarray
:return: the derivatives of data
:rtype: np.ndarray
"""
derivative_data = self.solve(data=data)
return derivative_data