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