# <<BEGIN-copyright>>
# <<END-copyright>>
"""
This module contains miscellaneous routines for supporting the endl2dmath class from module endl2dmathClasses.
Module's global variables::
endl2d_repr_xFormat : Default format string used for converting an x-element of the data of an endl2dmath object to a string.
endl2d_repr_yFormat : Default format string used for converting an y-element of the data of an endl2dmath object to a string.
"""
import math
import types
from fudge.core import fudgemisc
import endl2dmathClasses
import endlParameters
from fudge.core.math import fudgemath
endl2d_repr_xFormat = "%14.7e"
endl2d_repr_yFormat = "%12.5e"
[docs]def get2dmathData( object, callingRoutine, msg ) :
"""If object is a python list, then it is returned. Else, if object is a sub-class
of endl2dmath its data is return otherwise a raise is executed."""
if( type( object ) == type( [] ) ) :
data = object
else :
data = valid2dClassType( object, callingRoutine, msg ).data
return( data )
[docs]def valid2dClassType( object, callingRoutine, msg ) :
"""Returns the first argument, object, if it is a subclass of the endl2dmath class; else, triggers a raise."""
if( isinstance( object, endl2dmathClasses.endl2dmath ) ) : return( object )
raise Exception( "\nError in %s: invalid type = %s for %s" % ( callingRoutine, type( object ), msg ) )
[docs]def check2dPoint( p ) :
"""check2dPoint( p )\n Checks that p is the list [ number, number ]. A raise is generated if p is not."""
if( ( type( p ) != type( [] ) ) or ( len( p ) != 2 ) or not( fudgemath.isNumber( p[0] ) ) or not( fudgemath.isNumber( p[1] ) ) ) :
raise Exception( "\nError in check2dPoint: data point not list of [ number, number ]:", p )
[docs]def check2dData( data, allowNegativeX = True, allowSameX = False, allowZeroX = True, positiveY = True, printWarning = True, printErrors = True,
xCloseEps = None, formatCode = 12, maxAbsFloatValue = None ) :
"""Checks that data (or data.data if data is an endl2dmath instance) is a list containing
valid endl2dmath points. If positiveY is true (default) then a negative y value
generates a raise. If positiveY is false then the number of negative y values
is counted and the count is returned. The x values must be in increasing order
else a raise is generated. If allowZeroX is False, a zero first x value will generate a
raise. Also all x values must be positive unless allowNegativeX
is true. If printWarning is true than close x value warnings are printed. X values
x1 and x2 are considered close if x2 - x1 < xCloseEps * ( x1 + x2 ). If on input xCloseEps
is None, the it is set to endlParameters.endlEpsx * pow( 0.1, formatCode - 14 ) / 3.
If on input formatCode is 12 then it is set to 14.
This function returns an integer and an array. The integer is the number of negative
y values and the array is a list of all indices with close x-values. If printWarning
is 'True' than warnings are printed (warning will not execute a raise). If printErrors
is 'True' then errors are printed, and if at least one error exist, a raise will be
executed after all data has been checked."""
points = get2dmathData( data, "check2dData", "" )
l = len( points )
w = 0 # Number of close energy warnings.
ne = 0 # Number of points with negative y values.
badXIndicies = []
messages = []
if( xCloseEps is None ) :
if( formatCode == 12 ) : formatCode = 14
xCloseEps = endlParameters.endlEpsx / 3. * pow( 0.1, ( formatCode - 14 ) )
if( l > 0 ) :
pl = None
pn = points[0]
i = 0
if( pn[0] == 0. ) and ( not allowZeroX ) :
s = 'check2dData: zero x value at index = 0'
messages.append( s )
elif( ( pn[0] < 0. ) and ( not allowNegativeX ) ) :
s = 'check2dData: negative x value at index = %d, x = %e' % ( i, pn[0] )
messages.append( s )
if( printErrors ) : fudgemisc.printWarning( '\n'.join( fudgemisc.checkMessagesToString( s, indentation = ' ' ) ) )
while 1 :
check2dPoint( pn )
fudgemath.checkNumber( pn[0], "check2dData: x[%d]" % i, messages = messages, indentation = ' ', printErrors = printErrors,
maxAbsFloatValue = maxAbsFloatValue )
fudgemath.checkNumber( pn[1], "check2dData: y[%d]" % i, messages = messages, indentation = ' ', printErrors = printErrors,
maxAbsFloatValue = maxAbsFloatValue )
if( pl is not None ) :
checkCloseness = True
if( pl[0] >= pn[0] ) :
checkCloseness = False
if( not allowSameX ) :
s = 'check2dData: x value x[%d] >= x[%d] (%.8e >= %.8e)' % ( i - 1, i, pl[0], pn[0] )
messages.append( s )
if( printErrors ) : fudgemisc.printWarning( '\n'.join( fudgemisc.checkMessagesToString( s, indentation = ' ' ) ) )
if( checkCloseness and ( pn[0] - pl[0] ) < ( ( pn[0] + pl[0] ) * xCloseEps ) ) :
sn = endl2d_repr_xFormat % pn[0]
sl = endl2d_repr_xFormat % pl[0]
if( sn[-4:] == 'e-00' ) : sn = sn[:-4] + 'e+00'
if( sl[-4:] == 'e-00' ) : sl = sl[:-4] + 'e+00'
if( sn == sl ) :
badXIndicies.append( i )
s = 'check2dData: x values %.16e and %.16e very close' % ( pl[0], pn[0] )
messages.append( s )
if( ( w < 2 ) and printWarning ) : fudgemisc.printWarning( ' Warning in %s' % s )
w = w + 1
if( pn[1] < 0 ) :
if( positiveY ) :
s = 'check2dData: negative y value at index = %d (x = %e y = %e)' % ( i, pn[0], pn[1] )
messages.append( s )
if( printErrors ) : fudgemisc.printWarning( '\n'.join( fudgemisc.checkMessagesToString( s, indentation = ' ' ) ) )
ne = ne + 1
i = i + 1
if( i == l ) : break
pl = pn
pn = points[i]
if( printWarning and ( w > 1 ) ) : fudgemisc.printWarning( " Warning in check2dData: %d x values very close" % w )
if( printErrors and ( ( len( messages ) - w ) > 0 ) ) : raise Exception( ' bad 2d data' )
return ne, badXIndicies, messages
[docs]def interpolate2dPoints( interpolation, x, x1, y1, x2, y2 ) :
"""This function interpolates y at x where x1 <= x <= x2 using interpolation."""
y = y1
if( y1 != y2 ) :
if( x1 == x2 ) : # A why to handle a step function.
y = 0.5 * ( y1 + y2 )
elif( x == x1 ) :
pass
elif( x == x2 ) :
y = y2
else :
if( interpolation == 0 ) : # linear-linear interpolation.
y = ( ( x - x1 ) * y2 + ( x2 - x ) * y1 ) / ( x2 - x1 )
elif ( interpolation == 1 ) : # log-linear interpolation.
y = ( y2 - y1 ) * math.log( x / x1 ) / math.log( x2 / x1 ) + y1
elif ( interpolation == 2 ) : # linear-log interpolation.
y = y1 * math.exp( math.log( y2 / y1 ) * ( x - x1 ) / ( x2 - x1 ) )
elif ( interpolation == 3 ) : # log-log interpolation.
y = y1 * math.exp( math.log( y2 / y1 ) * math.log( x / x1 ) / math.log( x2 / x1 ) )
else :
raise Exception( "\nError in endl2dmathmisc.interpolate_points: unsupported interpolation = %s" % `interpolation` )
return( y )
[docs]def interpolate_X( data, x ) :
"""For internal use only."""
l = len( data )
if( l < 1 ) : return( 0. )
y = 0.
if( data[0][0] <= x <= data[-1][0] ) :
i = 0
while( i < l ) :
if( x <= data[i][0] ) : break
i += 1
if( x == data[i][0] ) :
y = data[i][1]
else :
d1 = data[i-1]
d2 = data[i]
x1 = d1[0];
x2 = d2[0];
y1 = d1[1];
y2 = d2[1];
if( x2 == x1 ) :
y = y1
else :
y = ( y1 * ( x2 - x ) + y2 * ( x - x1 ) ) / ( x2 - x1 )
return( y )
[docs]def dulledges2d( d ) :
"""Calls dullLowerEdge2d( d ) and dullUpperEdge2d( d )."""
dullLowerEdge2d( d )
dullUpperEdge2d( d )
return d
[docs]def dullLowerEdge2d( d ) :
"""Dulls the lower edge of d if the y value at that edge is not 0. d must be an endl2dmath instance."""
if ( len( d.data ) == 0 ) or ( d.data[0][1] == 0. ) : return d
x = d.data[0][0]
if ( x > endlParameters.endlSmallx ) :
if ( len( d.data ) > 1 ) and ( d.data[1][0] > ( 1. + 2. * endlParameters.endlEpsx ) * x ) :
d.data[0][1] = d[( 1. + endlParameters.endlEpsx ) * x]
d.data[0][0] = ( 1. + endlParameters.endlEpsx ) * x
d.data.insert( 0, [ ( 1 - endlParameters.endlEpsx ) * x, 0. ] )
return d
[docs]def dullUpperEdge2d( d ) :
"""Dulls the upper edge of d if the y value at that edge is not 0. d must be an endl2dmath instance."""
if ( len( d.data ) == 0 ) or ( d.data[-1][1] == 0. ) : return d
x = d.data[-1][0]
if ( x < endlParameters.endlLargex ) :
if ( len( d.data ) > 1 ) and ( d.data[-2][0] < ( 1. - 2. * endlParameters.endlEpsx ) * x ) :
d.data[-1][1] = d[( 1. - endlParameters.endlEpsx ) * x]
d.data[-1][0] = ( 1. - endlParameters.endlEpsx ) * x
d.data.append( [ ( 1 + endlParameters.endlEpsx ) * x, 0. ] )
return d
[docs]def gauss2d( xc, xw, f = 1e-3, xsigmax = 4, xMin = None, xMax = None ) :
"""Returns an endl2dmath instance the has a Gaussian shape centered at xc with
width xw (i.e., Exp( -( (x - xc) / xw )^2 / 2 ). The points are generated
such that linear interpolation is accurate to f. The x value's range is limited
by xMin, xMax and abs(x - xc) / xw <= xsigmax."""
def addpoints( xmin, xmax, x, a, f, DoPoints ) :
fs = 1.
if ( x > 1. ) : fs = -1.
y = ( 1. + fs * f ) * math.exp( -x * x / 2. )
xl = -x
yl = y
if ( DoPoints[0] ) :
if ( xl < xmin ) :
DoPoints[0] = 0
if ( len( a ) > 0 ) :
yl = ( ( xmin - a[0][0] ) * yl + ( xl - xmin ) * a[0][1] ) / float( xl - a[0][0] )
xl = xmin
if ( xl >= xmin ) : a.insert( 0, [ xl, yl ] )
if ( DoPoints[1] ) :
if ( x > xmax ) :
DoPoints[1] = 0
if ( len( a ) > 0 ) :
y = ( ( xmax - a[-1][0] ) * y + ( x - xmax ) * a[-1][1] ) / float( x - a[-1][0] )
x = xmax
if ( x <= xmax ) : a.append( [ x, y ] )
def nextpoint( xs_, xe_, f ) :
xs = xs_
xe = xe_
fs = -1
if ( xs < 1 ) : fs = 1
x1 = xs
y1 = math.exp( -x1 * x1 / 2. )
for i in range( 40 ) :
x = ( xs + xe ) / 2.
if ( xe - x < 1e-14 * x ) : break
s = ( math.exp( ( x1 - x ) * ( x1 + x ) / 2. ) - 1. ) / ( x - x1 )
z = 0.5 * ( - 1. / s + x1 )
xr = z - fs * math.sqrt( z * z - 1 )
r = 1. - ( s * ( xr - x1 ) + 1. ) * y1 * math.exp( xr * xr / 2. )
if( fs * r < 2. * f ) :
xs = x
else :
xe = x
return x
f = max( min( 1e-2, f ), 1e-8 )
if( xMin is None ) :
xmin = -xsigmax
else :
xmin = ( xMin - xc ) / float( xw )
if( xMax is None ) :
xmax = xsigmax
else :
xmax = ( xMax - xc ) / float( xw )
if( xmin < 0 ) and ( xmin < -xsigmax ) : xmin = -xsigmax
if( xsigmax < xmax ) : xmax = xsigmax
if ( xmin <= 0 ) and ( xmax >= 0 ) :
xamin = 0
else :
xamin = min( abs( xmin ), abs( xmax ) )
xamax = max( abs( xmin ), abs( xmax ) )
DoPoints = [ 1, 1 ]
a = []
if ( xamin != 0. ) : # The curve does not go through x = 0.
addpoints( xmin, xmax, xamin, a, f, DoPoints )
x1 = xamin
else :
a.append( [0., 1.] ) # The center and the next point.
s = math.sqrt( 2. * ( math.sqrt( 1. + 6. * f ) - 1. ) / 3. )
x0 = 2. * s / ( 1. - s * s )
s = ( math.exp( -x0 * x0 / 2. ) -1. ) / x0
c = s * ( s * x0 + 1. )
b = x0 + c
c = x0 * x0 / 2. + c * x0 + 0.5
e = ( -b + math.sqrt( b * b + 4. * f * c ) ) / ( 2. * c )
x1 = x0 + e
addpoints( xmin, xmax, x1, a, f, DoPoints )
e = pow( 3. * f / 2., 1. / 3. ) # Treat the points around x = 1 as special.
e = pow( 3. * f / 2. / ( 1. - 9. * e / 8. + 63. * e * e / 40. ), 2. / 3. )
xs = 0.
xe = 1.
for i in range( 40 ) :
x = ( xs + xe ) / 2.
if ( xe - x < 1e-14 * x ) : break
d = 1. - math.exp( ( x * x - 1. ) / 2. ) * ( 1. - ( 1. - e ) * ( x - 1. ) ) + f
if ( d < 0. ) :
xs = x
else :
xe = x
x = 1. - .9 * ( 1 - x )
x1m = x # Point < 1.
x1p = 1. + ( 1. - x ) # Point > 1.
x = x1
if ( x < x1m ) :
while ( x < 1. ) :
xs = x
x = nextpoint( xs, x1m, f )
if( abs( x1m - x ) < 1e-2 * ( x - xs ) ) : break
addpoints( xmin, xmax, x, a, f, DoPoints )
if ( x < x1m ) :
addpoints( xmin, xmax, x1m, a, f, DoPoints )
x = x1m
if ( x < x1p ) :
addpoints( xmin, xmax, x1p, a, f, DoPoints )
x = x1p
while ( x < xamax ) :
x = nextpoint( x, xamax + 1, f )
addpoints( xmin, xmax, x, a, f, DoPoints )
for e in a : e[0] = e[0] * xw + xc
return endl2dmathClasses.endl2dmath( a, checkDataType = 0 )
[docs]def gauss2dDull( xc, xw, f = 1e-3, xsigmax = 4, xMin = None, xMax = None ) :
"""Returns an endl2dmath instance of gauss2d( xc, xw, f, xsigmax, xMin, xMax ) that is dulled using dulledges2d( )."""
return dulledges2d( gauss2d( xc, xw, f, xsigmax, xMin, xMax ) )
[docs]def point2d( x ) :
"""Returns an endl2dmath instance of [ x, 1. ]."""
return endl2dmathClasses.endl2dmath( [ [ x, 1. ] ], checkDataType = 0 )
[docs]def point2dDull( x ) :
"""Returns an endl2dmath instance of point2d( x ) that is dulled using dulledges2d( )."""
return dulledges2d( point2d( x ) )
[docs]def flattop2d( xMin, xMax ) :
"""Returns an endl2dmath instance of a flattop from xMin to xMax with height 1. 0 < xMin < xMax."""
if( xMin >= xMax ) : raise Exception( "\nError in flattop2d: ( xMin < 0 ) or xMin >= xMax" )
return endl2dmathClasses.endl2dmath( [ [ xMin, 1. ], [ xMax, 1. ] ], checkDataType = 0 )
[docs]def flattop2dDull( xMin, xMax ) :
"""Returns the endl2dmath instances returned by flattop2d( xMin, xMax ) and that is dulled using dulledges2d( )."""
return dulledges2d( flattop2d( xMin, xMax ) )
[docs]def ramp2d( xMin, xMax, down = 0 ) :
"""Returns an endl2dmath instance of a ramp from xMin to xMax with height
0 at xMin and 1 at xMax if down = 0, and 1 at xMin and 0 at xMax otherwise."""
if( xMin >= xMax ) : raise Exception( "\nError in ramp2d: xMin >= xMax" )
if( down != 0 ) :
return endl2dmathClasses.endl2dmath( [ [ xMin, 1. ], [ xMax, 0. ] ], checkDataType = 0 )
else :
return endl2dmathClasses.endl2dmath( [ [ xMin, 0. ], [ xMax, 1. ] ], checkDataType = 0 )
[docs]def ramp2dDull( xMin, xMax, down = 0 ) :
"""Returns the endl2dmath instances returned by ramp2d( xMin, xMax, ramp ) and that is dulled
using dulledges2d( )."""
return dulledges2d( ramp2d( xMin, xMax, ramp ) )
[docs]def triangle2d( x0, x1, x2 ) :
"""Returns an endl2dmath instance of a triangle that starts at x0 with height 0,
ramps up to 1 at x1 and ramps down to 0 at x2."""
if( ( x0 > x1 ) or ( x1 > x2 ) ) :
raise Exception( "\nError in triangle2d: triangle points must have x0 <= x1 <= x2" )
return endl2dmathClasses.endl2dmath( [ [ x0, 0. ], [ x1, 1. ], [ x2, 0. ] ], checkDataType = 0 )
[docs]def convertFunctionToLinLin( f, xMin, xMax, accuracy ):
'''
Returns a new endl2dmath object that represents the function f(x) as a lin-lin table to
the specified absolute accuracy on the interval xMin to xMax.
'''
def convertFunctionToLinLin2( f, xMin, yMin, xMax, yMax, accuracy, data, level ) :
if( level > 16 ) : return
xMid, yMid = 0.5 * ( xMin + xMax ), 0.5 * ( yMin + yMax )
y = f( xMid )
if( abs( yMid - y ) > abs( y * accuracy ) ) :
convertFunctionToLinLin2( f, xMin, yMin, xMid, y, accuracy, data, level + 1 )
data.append( [ xMid, y ] )
convertFunctionToLinLin2( f, xMid, y, xMax, yMax, accuracy, data, level + 1 )
data = [ [ xMin, f( xMin ) ] ]
convertFunctionToLinLin2( f, xMin, data[0][1], xMax, f( xMax ), accuracy, data, 0 )
data.append( [ xMax, f( xMax ) ] )
return( endl2dmathClasses.endl2dmath( data = data ) )
[docs]def convertLogLogToLinLin( self, accuracy, diSectionMax = 3 ) :
"""Returns a new endl2dmath object of self converted from log-log to lin-lin to specified accuracy."""
if( accuracy < 1e-4 ) : accuracy = 1e-4
if( diSectionMax > 8 ) : diSectionMax = 8
def convertLogLogToLinLin2( self, x1, y1, x2, y2, level ) :
"""For internal use."""
level_ = level + 1
if( level_ > diSectionMax ) : return
u2 = x2 / x1
v2 = y2 / y1
logXs = math.log( u2 )
logYs = math.log( v2 )
logYXs = logYs / logXs
u = logYXs * ( u2 - v2 ) / ( ( 1 - logYXs ) * ( v2 - 1 ) )
vLog = math.pow( u, logYXs )
vLin = ( u2 - u + v2 * ( u - 1 ) ) / ( u2 - 1 )
if( abs( vLin / vLog - 1. ) > accuracy ) :
x = x1 * u
y = y1 * vLog
self.setValue( x, y )
convertLogLogToLinLin2( self, x1, y1, x, y, level_ )
convertLogLogToLinLin2( self, x, y, x2, y2, level_ )
new = self.copyData( )
x1 = None
for x2, y2 in self.data :
if( x1 is not None ) :
if( ( x1 <= 0 ) or ( y1 <= 0 ) or ( x2 <= 0 ) or ( y2 <= 0 ) ) :
raise Exception( "One of x1 = %e y1 = %e x2 = %d or y2 = %e is <= 0." % ( x1, y1, x2, y2 ) )
if( ( x1 != x2 ) and ( y1 != y2 ) ) : convertLogLogToLinLin2( new, x1, y1, x2, y2, 0 )
x1, y1 = x2, y2
return( new )