Calculate jacobian hessian matrix in sympy exponential decay

Calculate Jacobian and Hessian matrix in python sympy for exponential decay function

See also:

1. https://en.wikipedia.org/wiki/Propagation_of_uncertainty
2. http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant
3. http://en.wikipedia.org/wiki/Hessian_matrix
4. http://maxima-online.org/articles/hessian.html
5. http://certik.github.io/scipy-2013-tutorial/html/tutorial/basic_operations.html
6. http://scipy-lectures.github.io/advanced/sympy.html
7. http://docs.sympy.org/dev/gotchas.html
8. https://github.com/sympy/sympy/wiki/Faq

Sumpy python installation

Consider for example installing Enthought Canopy

Tutorial with function for weighted difference between function evaluation with fitted parameters and measured values.

Created by:

Troels Emtekær Linnet
PhD student
Copenhagen University
SBiNLab

run by: python sympy_test.py

Tutorial with quadratic chi2 function

Created by:

Troels Emtekær Linnet
PhD student
Copenhagen University
SBiNLab

run by: python sympy_test.py

# Tutorial from: http://scipy-lectures.github.io/advanced/sympy.html

from sympy import *

# In contrast to other Computer Algebra Systems, in SymPy you have to declare symbolic variables explicitly:
x = Symbol('x')
y = Symbol('y')

# You can differentiate any SymPy expression using diff(func, var). Examples:
print("First some tests")
f = sin(x)
print("Function is", f)
print("df/dx is:", diff(f, x))

f = x**2
print("Function is", f)
print("df/dx is:", diff(f, x))

print("")

# Define exponential decay.
i0 = Symbol('i0')
times = Symbol('times')
r2eff = Symbol('r2eff')

print("Function to calculate intensity decay at time[i].")
exp_dec = i0 * exp( -times * r2eff)
print("Function f(r2eff, i0) is:", simplify(exp_dec))
print("df/d_i0 should be: exp(-r2eff*times)", simplify(diff(exp_dec, i0)) )
print("df/d_r2eff should be: -i0*times*exp(-r2eff*times)", simplify(diff(exp_dec, r2eff)) )

print("")

print("Function to calculate the chi2.")
# Function to calculate the chi2.
# The standard deviation of the measured values.
errors = Symbol('errors')
# The measured values.
values = Symbol('values')
# back_calc, the function which return the values, according to parameters.
back_calc = exp_dec
chi2 = (1.0 / errors * (values - back_calc))**2
print("The chi2 function is:", chi2)

print("")

print("Now calculate the Jacobian. The partial derivative matrix.\n")
print("Jacobian is m rows with function derivatives and n columns of parameters.")

# If you want to test for symbolic equality, one way is to subtract one expression from the other and run it through functions like expand(), simplify(), and trigsimp() and see if the equation reduces to 0.
# https://github.com/sympy/sympy/wiki/Faq - Why does SymPy say that two equal expressions are unequal?

print("Derivative of chi2 function for exponential decay, with respect to r2eff")
d_chi2_d_r2eff = diff(chi2, r2eff)
print("d_chi2_d_r2eff=", d_chi2_d_r2eff)
print("Should be: (2 * i0 * times * (values - i0 / exp(r2eff * times) ) ) / ( exp(r2eff * times) * self.errors**2) ")
val_d_chi2_d_r2eff = ( 2 * i0 * times * (values - i0 / exp(r2eff * times) ) ) / ( exp(r2eff * times) * errors**2)
print("Are they equal?: ", simplify(d_chi2_d_r2eff) == simplify(val_d_chi2_d_r2eff) )

print("")

print("Derivative of chi2 function for exponential decay, with respect to i0")
d_chi2_d_i0 = diff(chi2, i0)
print("d_chi2_d_i0=", d_chi2_d_i0)
print("Should be: - ( 2. * ( values - i0 / exp(r2eff * times) ) ) / ( exp(r2eff * times) * errors**2) ")
val_d_chi2_d_i0 = - ( 2. * ( values - i0 / exp(r2eff * times) ) ) / ( exp(r2eff * times) * errors**2) 
print("Are they equal?: ", simplify(d_chi2_d_i0) == simplify(val_d_chi2_d_i0) )

print("\n")

print("""Form the Jacobian matrix by:
------------------------------------------------------------------------------
from numpy import array, transpose

d_chi2_d_r2eff = %s
d_chi2_d_i0 = %s
jacobian_matrix = transpose(array( [d_chi2_d_r2eff , d_chi2_d_i0] ) )
------------------------------------------------------------------------------
""" % (d_chi2_d_r2eff, d_chi2_d_i0) )

print("")

print("Now calculate the Hessian. The second-order partial derivatives.\n")
print("See for example: http://maxima-online.org/articles/hessian.html")
print("If all second partial derivatives of f exist, then the Hessian matrix of f is the matrix:")
print("Hf(x)i,j = d**2 f(x) / ( dx_i * dx_j )")
print("We need to do:")
print("diff(f(r2eff, i0), r2eff, r2eff)")
print("diff(f(r2eff, i0), r2eff, i0)")
print("diff(f(r2eff, i0), i0, r2eff)")
print("diff(f(r2eff, i0), i0, i0)")
d2_chi2_d_r2eff_d_r2eff = diff(chi2, r2eff, r2eff)
d2_chi2_d_r2eff_d_i0 = diff(chi2, r2eff, i0)
d2_chi2_d_i0_d_r2eff = diff(chi2, i0, r2eff)
d2_chi2_d_i0_d_i0 = diff(chi2, i0, i0)
print("d2_chi2_d_r2eff_d_r2eff=", d2_chi2_d_r2eff_d_r2eff)
print("d2_chi2_d_r2eff_d_i0=", d2_chi2_d_r2eff_d_i0)
print("d2_chi2_d_i0_d_r2eff=", d2_chi2_d_i0_d_r2eff)
print("d2_chi2_d_i0_d_i0=", d2_chi2_d_i0_d_i0)

print("\n")

print("""Form the Hessian matrix by:
------------------------------------------------------------------------------
from numpy import array, transpose

d2_chi2_d_r2eff_d_r2eff = %s
d2_chi2_d_r2eff_d_i0 = %s
d2_chi2_d_i0_d_r2eff = %s
d2_chi2_d_i0_d_i0 = %s
hessian_matrix = transpose(array( [d2_chi2_d_r2eff_d_r2eff, d2_chi2_d_r2eff_d_i0, d2_chi2_d_i0_d_r2eff, d2_chi2_d_i0_d_i0] ) )
------------------------------------------------------------------------------
""" % (d2_chi2_d_r2eff_d_r2eff,d2_chi2_d_r2eff_d_i0, d2_chi2_d_i0_d_r2eff, d2_chi2_d_i0_d_i0) )