/usr/lib/python2.7/dist-packages/woo/linterpolation.py is in python-woo 1.0+dfsg1-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 | # 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.
"""
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
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
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 "Integral not reached within the interpolation range!"
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
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
|