Python functions for iterative application of the Faa di Bruno formula

FaaDiBrunoFromBellPolynomials.py

## 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 )

FaaDiBrunoFromPartitions.py

## 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} )