#* **************************************************************************
#*
#* CDTK, Chemical Dynamics Toolkit
#* A modular system for chemical dynamics applications and more
#*
#* Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016
#* Oriol Vendrell, DESY, <oriol.vendrell@desy.de>
#*
#* Copyright (C) 2017, 2018, 2019
#* Ralph Welsch, DESY, <ralph.welsch@desy.de>
#*
#* Copyright (C) 2020, 2021, 2022, 2023
#* Ludger Inhester, DESY, ludger.inhester@cfel.de
#*
#* This file is part of CDTK.
#*
#* CDTK 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/>.
#*
#* **************************************************************************
from numpy import array, zeros, identity, dot, arange
from numpy.linalg import solve
from .Mathematics import richardson
"""
Routines to take numerical first and second derivatives of a given function The
Richardson method is available to produce the best possible extrapolation of the finite
difference formulas. In case the richardson method is not necessary the global variable
Fast can be set to True. In such a case the richardson routine evaluates the finite
difference only once using the initial value of the step and returns.
"""
[docs]
def gradient(f,x,ms=None,step=0.1,fast=True):
"""Return an array with the gradient in x along the directions in modes
f -- function with respect the gradient is taken
x -- coordinates to take the gradient
ms -- list of arrays with orthonormal directions
step -- step to take in the numerical differentiation
"""
if type(ms) == type(None):
ms = identity(len(x))
g=[]
for m in ms:
g.append(_directionalDerivative(f,x,m,step,fast=fast))
return array(g)
[docs]
def hessian(f,x,ms=None,step=0.1,fast=True):
"""Return an array with the hessian in x along the directions in modes
f -- function with respect the hessian taken
x -- coordinates to take the hessian
ms -- list of arrays with orthonormal directions
step -- step to take in the numerical differentiation
"""
if type(ms) == type(None):
ms = identity(len(x))
h=zeros((len(ms),len(ms)),float)
for i in range(len(ms)):
for j in range(0,i+1):
if i == j:
h[i][j] = _secondDerivative(f,x,ms[i],step,fast=fast)
else:
h[i][j] = _crossedDerivative(f,x,ms[i],ms[j],step,fast=fast)
h[j][i] = h[i][j]
return h
def _directionalDerivative(func,x,m,step=0.1,fast=True):
"""Here m is just an array representing the direction to take
the derivative
"""
def findiffd1(hh):
return (func(x+m*hh)-func(x-m*hh))/(2.0*hh)
d,err,hf = richardson(findiffd1,step,fast=fast)
return d
def _secondDerivative(func,x,m,step=0.1,fast=True):
def findiffd2(hh):
return (func(x+m*hh) + func(x-m*hh) - 2*func(x))/(hh**2)
d,err,hf = richardson(findiffd2,step,fast=fast)
return d
def _crossedDerivative(func,x,m1,m2,step=0.1,fast=True):
def findiffd2d2(hh):
x1i2i = x - hh*m1 - hh*m2
x1i2f = x - hh*m1 + hh*m2
x1f2i = x + hh*m1 - hh*m2
x1f2f = x + hh*m1 + hh*m2
e1i2i = func(x1i2i)
e1i2f = func(x1i2f)
e1f2i = func(x1f2i)
e1f2f = func(x1f2f)
# e1i2i = _energy(func,x1i2i)
# e1i2f = _energy(func,x1i2f)
# e1f2i = _energy(func,x1f2i)
# e1f2f = _energy(func,x1f2f)
return ((e1f2f - e1f2i) - (e1i2f - e1i2i))/(4.0*hh**2)
d,err,hf = richardson(findiffd2d2,step)
return d
def _gradientN(f,x,n,ms=None,step=0.3,fast=True):
"""Return an array with the Nth derivative in x along the
directions in modes
f -- function with respect the Nth der. is taken
x -- coordinates to take the gradient
n -- order of the derivative
ms -- list of arrays with orthonormal directions
step -- step to take in the numerical differentiation
"""
if type(ms) == type(None):
ms = identity(len(x))
g=[]
for m in ms:
g.append(_nthDerivR(f,x,n,m,step,fast=fast))
return array(g)
def _nthDerivR(f,x,n,m,step=0.3,fast=True):
h = step
def findiff(hh):
return _nthDeriv(f,x,n,m,hh)
d,err,hf = richardson(findiff,step,fast=True)
return d
def _nthDeriv(f,x,n,m,step=0.1):
hh = step
if n > 1:
return (_nthDeriv(f,x+m*hh,n-1,m,hh)
- _nthDeriv(f,x-m*hh,n-1,m,hh))/(2.0*hh)
if n == 1:
return (f(x+m*hh)-f(x-m*hh))/(2.0*hh)
if n == 0:
return f(x)