Source code for yade.linterpolation

# encoding: utf-8
#
# © 2009 Václav Šmilauer <eudoxos@arcig.cz>
#

"""
Module for rudimentary support of manipulation with piecewise-linear
functions (which are usually interpolations of higher-order functions,
whence the module name). Interpolation is always given as two lists
of the same length, where the x-list must be increasing.

Periodicity is supported by supposing that the interpolation
can wrap from the last x-value to the first x-value (which
should be 0 for meaningful results).

Non-periodic interpolation can be converted to periodic one
by padding the interpolation with constant head and tail using
the sanitizeInterpolation function.

There is a c++ template function for interpolating on such sequences in
pkg/common/Engine/PartialEngine/LinearInterpolate.hpp (stateful, therefore
fast for sequential reads).

TODO: Interpolating from within python is not (yet) supported.
"""
from __future__ import print_function

from builtins import range
[docs]def revIntegrateLinear(I,x0,y0,x1,y1): """Helper function, returns value of integral variable x for linear function f passing through (x0,y0),(x1,y1) such that 1. x∈[x0,x1] 2. ∫_x0^x f dx=I and raise exception if such number doesn't exist or the solution is not unique (possible?) """ from math import sqrt dx,dy=x1-x0,y1-y0 if dy==0: # special case, degenerates to linear equation return (x0*y0+I)/y0 a=dy/dx; b=2*(y0-x0*dy/dx); c=x0**2*dy/dx-2*x0*y0-2*I det=b**2-4*a*c; assert(det>0) p,q=(-b+sqrt(det))/(2*a),(-b-sqrt(det))/(2*a) pOK,qOK=x0<=p<=x1,x0<=q<=x1 if pOK and qOK: raise ValueError("Both radices within interval!?") if not pOK and not qOK: raise ValueError("No radix in given interval!") return p if pOK else q
[docs]def integral(x,y): """Return integral of piecewise-linear function given by points x0,x1,… and y0,y1,… """ assert(len(x)==len(y)) sum=0 for i in range(1,len(x)): sum+=(x[i]-x[i-1])*.5*(y[i]+y[i-1]) return sum
[docs]def xFractionalFromIntegral(integral,x,y): """Return x within range x0…xn such that ∫_x0^x f dx==integral. Raises error if the integral value is not reached within the x-range. """ sum=0 for i in range(1,len(x)): diff=(x[i]-x[i-1])*.5*(y[i]+y[i-1]) if sum+diff>integral: return revIntegrateLinear(integral-sum,x[i-1],y[i-1],x[i],y[i]) else: sum+=diff raise ValueError("Integral not reached within the interpolation range!")
[docs]def xFromIntegral(integralValue,x,y): """Return x such that ∫_x0^x f dx==integral. x wraps around at xn. For meaningful results, therefore, x0 should == 0 """ from math import floor period=x[-1]-x[0] periodIntegral=integral(x,y) numPeriods=floor(integralValue/periodIntegral) xFrac=xFractionalFromIntegral(integralValue-numPeriods*periodIntegral,x,y) #print '### wanted _%g; period=%g; periodIntegral=_%g (numPeriods=%g); rests _%g (xFrac=%g)'%(integralValue,period,periodIntegral,numPeriods,integralValue-numPeriods*periodIntegral,xFrac) #print '### returning %g*%g+%g=%g'%(period,numPeriods,xFrac,period*numPeriods+xFrac) return period*numPeriods+xFrac
[docs]def sanitizeInterpolation(x,y,x0,x1): """Extends piecewise-linear function in such way that it spans at least the x0…x1 interval, by adding constant padding at the beginning (using y0) and/or at the end (using y1) or not at all.""" xx,yy=[],[] if x0<x[0]: xx+=[x0]; yy+=[y[0]] xx+=x; yy+=y if x1>x[-1]: xx+=[x1]; yy+=[y[-1]] return xx,yy
if __name__=="main": xx,yy=sanitizeInterpolation([1,2,3],[1,1,2],0,4) print(xx,yy) print(integral(xx,yy)) # 5.5 print(revIntegrateLinear(.625,1,1,2,2)) # 1.5 print(xFractionalFromIntegral(1.625,xx,yy)) # 1.625 print(xFractionalFromIntegral(2.625,xx,yy)) # 2.5