Instantly share code, notes, and snippets.

Embed
What would you like to do?
Python functions for iterative application of the Faa di Bruno formula
## Original author: David Hoyle (david.hoyle@hoyleanalytics.org)
##
## Date: 2017-05-20
##
## Licence: CC-BY
from scipy.special import gammaln
import numpy as np
import math
import sympy
def calcFaaDiBrunoFromBellPolynomials( baseDerivatives, l, verbose=True ):
nMax = len( baseDerivatives )
# initialize derivatives at the fixed point of the base function
derivatives_atFixedPoint_tmp = np.zeros( (l, nMax) )
# set starting derivatives to the supplied base derivatives
derivatives_atFixedPoint_tmp[0, ] = baseDerivatives
# set number of function iterations
nIterations = l-1
# if we need to iterate then do so
if( nIterations > 0 ):
# generate and cache Bell polynomials
symbols_tmp = sympy.symbols('x:' + str(nMax+5) ) # We will generate some extra symbols for good measure
bellPolynomials = {}
for n in range(1, nMax+1):
for k in range(1, n+1):
bp_tmp = sympy.bell(n, k, symbols_tmp)
bellPolynomials[str(n) + '_' + str(k)] = bp_tmp
# Loop, iteratively applying Faa di Bruno formula
if( verbose ):
print( "Iterating Faa di Bruno formula" )
for iteration in range(nIterations):
if( verbose ):
print( "Evaluating derivatives for function iteration " + str(iteration+1) )
for n in range(1, nMax+1):
sum_tmp = 0.0
for k in range(1, n+1):
# retrieve kth derivative of base function at previous iteration
f_k_tmp = derivatives_atFixedPoint_tmp[0, k-1]
#evaluate arguments of Bell polynomials
bellPolynomials_key = str( n ) + '_' + str( k )
bp_tmp = bellPolynomials[bellPolynomials_key]
replacements = [( symbols_tmp[i], derivatives_atFixedPoint_tmp[iteration, i] ) for i in range(n-k+1) ]
sum_tmp = sum_tmp + ( f_k_tmp * bp_tmp.subs(replacements) )
derivatives_atFixedPoint_tmp[iteration+1, n-1] = sum_tmp
return( derivatives_atFixedPoint_tmp )
## Original author: David Hoyle (david.hoyle@hoyleanalytics.org)
##
## Date: 2017-05-20
##
## Licence: CC-BY
from sympy.utilities.iterables import partitions
from scipy.special import gammaln
import numpy as np
import math
# function to calculate the Hardy-Ramanujan asymptotic approximation to the number of partitions of n
# function returns the log of the approximate number of partitions
def calcHardyRamanujanApprox( n ):
log_nPartitions = np.pi * np.sqrt( 2.0 * float( n )/ 3.0 ) - np.log( 4.0 * float( n ) * np.sqrt( 3.0 ) )
return( log_nPartitions )
def calcFaaDiBrunoFromPartitions(baseDerivatives, l, thresholdForExp=30.0, verbose=True):
# get number of derivatives available
n = len( baseDerivatives )
# set up arrays for holding derivatives (on log scale)
# of the base function, the current iteration and next iteration
baseDerivativesLog = np.zeros( n )
baseDerivativesSign = np.ones( n )
currentDerivativesLog = np.zeros( n )
currentDerivativesSign = np.ones( n )
for k in range( n ):
if( np.abs( baseDerivatives[k] ) > 1.0e-12 ):
baseDerivativesLog[k] = np.log( np.abs(baseDerivatives[k]) )
baseDerivativesSign[k] = np.sign( baseDerivatives[k])
# get current derivatives
currentDerivativesLog[k] = np.log( np.abs( baseDerivatives[k] ) )
currentDerivativesSign[k] = np.sign( baseDerivatives[k] )
else:
baseDerivativesLog[k] = float( '-Inf' )
baseDerivativesSign[k] = 1.0
currentDerivativesLog[k] = float( '-Inf' )
currentDerivativesSign[k] = 1.0
nextDerivativesLog = np.zeros( n )
nextDerivativesSign = np.ones( n )
# initialize array for holding derivatives at each iteration
iteratedFunctionDerivativesLog = np.zeros( (l, n) )
iteratedFunctionDerivativesSign = np.ones( (l, n) )
# store base derivatives in first row of array
iteratedFunctionDerivativesLog[0, :] = baseDerivativesLog
iteratedFunctionDerivativesSign[0, :] = baseDerivativesSign
# set number of function iterations
nIterations = l-1
# if we need to iterate then do so
if( nIterations > 0 ):
# evaluate approximate number of paritions required
log_nPartitions = calcHardyRamanujanApprox( n )
if( verbose==True ):
print( "You have requested evaluation up to derivative " + str(n) )
if( log_nPartitions > np.log( 1.0e6 ) ):
print( "Warning: There are approximately " + str(int(np.round(np.exp(log_nPartitions)))) + " partitions of " + str(n) )
print( "Warning: Evaluation will be costly both in terms of memory and run-time" )
# store partitions
pStore = {}
for k in range( n ):
# get partition iterator
pIterator = partitions(k+1)
pStore[k] = [p.copy() for p in pIterator]
# loop over function iterations
for iteration in range( nIterations ):
if( verbose==True ):
print( "Evaluating derivatives for function iteration " + str(iteration+1) )
for k in range( n ):
faaSumLog = float( '-Inf' )
faaSumSign = 1
# get partitions
partitionsK = pStore[k]
for pidx in range( len(partitionsK) ):
p = partitionsK[pidx]
sumTmp = 0.0
sumMultiplicty = 0
parityTmp = 1
for i in p.keys():
value = float(i)
multiplicity = float( p[i] )
sumMultiplicty += p[i]
sumTmp += multiplicity * currentDerivativesLog[i-1]
sumTmp -= gammaln( multiplicity + 1.0 )
sumTmp -= multiplicity * gammaln( value + 1.0 )
parityTmp *= np.power( currentDerivativesSign[i-1], multiplicity )
sumTmp += baseDerivativesLog[sumMultiplicty-1]
parityTmp *= baseDerivativesSign[sumMultiplicty-1]
# now update faaSum on log scale
if( sumTmp > float( '-Inf' ) ):
if( faaSumLog > float( '-Inf' ) ):
diffLog = sumTmp - faaSumLog
if( np.abs(diffLog) <= thresholdForExp ):
if( diffLog >= 0.0 ):
faaSumLog = sumTmp
faaSumLog += np.log( 1.0 + (float(parityTmp*faaSumSign) * np.exp( -diffLog )) )
faaSumSign = parityTmp
else:
faaSumLog += np.log( 1.0 + (float(parityTmp*faaSumSign) * np.exp( diffLog )) )
else:
if( diffLog > thresholdForExp ):
faaSumLog = sumTmp
faaSumSign = parityTmp
else:
faaSumLog = sumTmp
faaSumSign = parityTmp
nextDerivativesLog[k] = faaSumLog + gammaln( float(k+2) )
nextDerivativesSign[k] = faaSumSign
# update accounting for proceeding to next function iteration
currentDerivativesLog[0:n] = nextDerivativesLog[0:n]
currentDerivativesSign[0:n] = nextDerivativesSign[0:n]
iteratedFunctionDerivativesLog[iteration+1, 0:n] = currentDerivativesLog[0:n]
iteratedFunctionDerivativesSign[iteration+1, 0:n] = currentDerivativesSign[0:n]
return( {'logDerivatives':iteratedFunctionDerivativesLog, 'signDerivatives':iteratedFunctionDerivativesSign} )
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment