creating a function that changes equations at certain slope, usable in curve_fit

Question

I am currently working with fitting decline curves to real-world production data. I have had good luck with creating a hyperbolic and using curve_fit from scipy.optimize. The current function I use:

def hyp_func(x,qi,b,di):
    return qi*(1.0-b*di*x)**(-1.0/b)

What I would like to do now, is at a certain rate of decline, transition to an exponential function. How would i go about this and still be able to use in curve_fit (I think below works)? I am trying the code below, is this the way to do it? or is there a better way?

def hyp_func2(x,qi,b,di):
    dlim = -0.003
    hy = qi*(1.0-b*di*x)**(-1.0/b)
    hdy = di/(1.0-b*di*x)
    ex = x[hdy>dlim]
    qlim = qi*(dlim/di)**(1/b)
    xlim = ((qi/qlim)**b-1)/(b*-di)
    ey = qlim*np.exp(dlim*(ex-xlim))
    y = np.concatenate((hy[hdy<dlim],ey))
    return y

hy is the hyperbolic equation
hdy is the hy derivative
ex is the part of x after derivative hits dlim
ey is the exponential equation

I am still working out the equations, I am not getting a continuous function.

edit: data here, and updated equations

The approach [described here](https://stackoverflow.com/q/29382903/832621) can give you some insight — Saullo G. P. Castro, Sep 22 '17 at 15:55
Thanks, i will take a look at that, and learn more about `numpy.piecewise()` — Ditchbuster, Sep 22 '17 at 17:42

score 1 · Answer 1 · answered Sep 22 '17 at 04:55

Sorry to be the bearer of bad news, but if I understand what you are trying to do, I think it is very difficult to have scipy.optimize.curve_fit, or any of the other methods from scipy.optimize do what you are hoping to do.

Most fitting algorithms are designed to work with continuous variables, and usually (and curve_fit for sure) start off by making very small changes in parameter values to find the right direction and step size to take to improve the result.

But what you're looking for is a discrete variable as the breakpoint between one functional form (roughly, "power law") to another ("exponential") The algorithm won't normally make a large enough change in your di parameter to make a difference for which value is used as that breakpoint, and may decide that di does not affect the fit (your model used di in other ways too, so you might get lucky and di might have an affect on the fit.

I guess I didn't explain it well, I was trying to keep out of the weeds. The function should be continuous and basically, it should limit the decline based on `dlim`. Once the first part is declining smaller than `dlim`, change to an exponential equation which has a constant decline. The second part is all based on the first part so changes to the variables should change the overall function. — Ditchbuster, Sep 22 '17 at 17:59

score 1 · Answer 2 · answered Sep 26 '17 at 11:17

Assuming that qi>0 the slope is actually positive, so I do not get the choice of -0.003. Moreover I think the derivative is wrong. You can calculate exactly the value where the lope reaches a critical value. Now, from my experience you have two choices. If you define a piecewise function yourself, you usually run into trouble with function calls using numpy arrays. I typically use scipy.optimize.leastsq with a self-defined residual function. A second option is a continuous transition between the two functions. You can make that as sharp as you want, as value and slope already fit, by definition. The two solutions look as follows

import matplotlib.pyplot as plt import numpy as np

def hy(x,b,qi,di):
    return qi*(1.0-b*di*x)**(-1.0/b)


def abshy(x,b,qi,di):#same as hy but defined for all x
    return qi*abs(1.0-b*di*x)**(-1.0/b)


def dhy(x,b,qi,di):#derivative of hy
    return qi*di*(1.0-b*di*x)**(-(b+1.0)/b)

def get_x_from_slope(s,b,qi,di):#self explaining
    return (1.0-(s/(qi*di))**(-b/(b+1.0)))/(b*di)


def exh(x,xlim,qlim,dlim):#exponential part (actually no free parameters)
    return qlim*np.exp(dlim*(x-xlim))

def trans(x,b,qi,di, s0):#piecewise function
    x0=get_x_from_slope(s0,b,qi,di)
    if x<x0:
        out= hy(x,b,qi,di)
    else:
        H0=hy(x0,b,qi,di)
        out=exh(x,x0,H0,s0/H0)
    return out


def no_if_trans(x,b,qi,di, s0,sharpness=10):#continuous transition between the two functions
    x0=get_x_from_slope(s0,b,qi,di)
    H0=hy(x0,b,qi,di)
    weight=0.5*(1+np.tanh(sharpness*(x-x0)))
    return weight*exh(x,x0,H0,s0/H0)+(1.0-weight)*abshy(x,b,qi,di)


xList=np.linspace(0,5.5,90)
hyList=np.fromiter(( hy(x,2.2,1.2,.1) for x in xList ) ,np.float)
t1List=np.fromiter(( trans(x,2.2,1.2,.1,3.59) for x in xList ) ,np.float)
nt1List=np.fromiter(( no_if_trans(x,2.2,1.2,.1,3.59) for x in xList ) ,np.float)


fig1=plt.figure(1)
ax=fig1.add_subplot(1,1,1)
ax.plot(xList,hyList)
ax.plot(xList,t1List,linestyle='--')
ax.plot(xList,nt1List,linestyle=':')

ax.set_ylim([1,10])
ax.set_yscale('log')
plt.show()

There is almost no differences in the two solutions, but your options for using scipy fitting functions are slightly different. The second solution should easily work with curve_fit

creating a function that changes equations at certain slope, usable in curve_fit

2 Answers2