python nodal/piecewise-linear power-law generator

Question

I need efficient python code that returns (not fits) a piecewise-linear (or, actually, piecewise power-law) continuous function for an arbitrary number of nodes (/poles/control points) defined by their positions plus (preferably) slopes rather than amplitudes. For example, for three pieces (four nodes) I have:

def powerLawFuncTriple(x,C,alpha,beta,gamma,xmin,xmax,x0,x1):

    """
    Extension of the two-power-law version described at
       http://en.wikipedia.org/wiki/Power_law#Broken_power_law
    """

    if x <= xmin or x > xmax:
        return 0.0
    elif xmin < x <= x0:
        n = C * (x ** alpha)
    elif x0 < x <= x1:
        n = C * x0**(alpha-beta) * (x ** beta)
    elif x1 < x <= xmax:
        n = C  * x0**(alpha-beta) * x1**(beta-gamma) * (x ** gamma)

    return n

Do any helpful functions already exist, otherwise what's an efficient way to code up code to generate these functions? Maybe this amounts to evaluating rather than fitting one of the scipy built-ins.

Somewhat related: Fitting piecewise function in Python

One possible answer may be:

def piecewise(x,n0,posns=[1.0,2.0,3.0],alpha=[-2.0,-1.5]):

    if x <= posns[0] or x > posns[-1]: return 0.0

    n=n0*x**alpha[0]
    np=len(posns)
    for ip in range(np):
        if posns[ip] < x <= posns[ip+1]: return n
        n *= (posns[ip+1]/float(x))**(alpha[ip]-alpha[ip+1])
    return n

But this must get slower as x increases. Would list comprehension or anything else speed up the loop?

Thanks!

score 0 · Accepted Answer · answered Jan 14 '15 at 05:58

In the end I've gone with the excellent CosmoloPy module:

http://pydoc.net/Python/CosmoloPy/0.1.103/cosmolopy.utils

class PiecewisePowerlaw()

"""
You can specify the intervals and power indices, and this class
will figure out the coefficients needed to make the function
continuous and normalized to unit integral.
"""

Seems to work swimmingly.

python nodal/piecewise-linear power-law generator

1 Answers1