How to fit a piecewise (alternating linear and constant segments) function to a parabolic function?

I do have a function, for example , but this can be something else as well, like a quadratic or logarithmic function. I am only interested in the domain of . The parameters of the function (a and k in this case) are known as well.

My goal is to fit a piecewise function to this, which contains alternating segments of linear functions (without intercepts) and constants. The first and last segments are both linear. And the number of segments should be pre-selected between around 9-29 (that is 5-15 linear steps + 4-14 constant plateaus).

Formally

• The input function:

• The fitted piecewise function:

I am looking for the optimal resulting parameters (c,r,b) if the segment numbers (n) are specified beforehand. The resulting constants (c) and the breakpoints (r) should be whole natural numbers, and the slopes (b) round two decimal point values.

I have tried to do the fitting numerically using the pwlf package using a segmented constant models, and further processed the resulting constant model with some graphical intuition to "slice" the constant steps with the slopes. It works to some extent, but I am sure this is suboptimal from both fitting perspective and computational efficiency. It takes multiple minutes to generate a fitting with 8 slopes on the range of 1-50000. I am sure there must be a better way to do this.

import numpy as np
import matplotlib.pyplot as plt
import pwlf

# The input function
def input_func(x,k,a):
    return np.power(x,1/a)*k

x = np.arange(1,5e4)
y = input_func(x, 1.8, 1.3)

plt.plot(x,y);

def pw_fit(func, x_r, no_seg, *fparams):
    
    # working on the specified range
    x = np.arange(1,x_r)
    
    y_input = func(x, *fparams)
    
    my_pwlf = pwlf.PiecewiseLinFit(x, y_input, degree=0)
    res = my_pwlf.fit(no_seg)
    yHat = my_pwlf.predict(x)
    
    # Function values at the breakpoints
    y_isec = func(res, *fparams)
    
    # Slope values at the breakpoints
    slopes = np.round(y_isec / res, decimals=2)
    slopes = slopes[1:]
    # For the first slope value, I use the intersection of the first constant plateau and the input function
    slopes = np.insert(slopes,0,np.round(y_input[np.argwhere(np.diff(np.sign(y_input - yHat))).flatten()[0]] / np.argwhere(np.diff(np.sign(y_input - yHat))).flatten()[0], decimals=2))

    plateaus = np.unique(np.round(yHat))
    
    # If due to rounding slope values (to two decimals), there is no change in a subsequent step, I just remove those segments
    to_del = np.argwhere(np.diff(slopes) == 0).flatten()
    slopes = np.delete(slopes,to_del + 1)
    plateaus = np.delete(plateaus,to_del)
    
    breakpoints = [np.ceil(plateaus[0]/slopes[0])]
    for idx, j in enumerate(slopes[1:-1]):
        breakpoints.append(np.floor(plateaus[idx]/j))
        breakpoints.append(np.ceil(plateaus[idx+1]/j))
    breakpoints.append(np.floor(plateaus[-1]/slopes[-1]))

    return slopes, plateaus, breakpoints

slo, plat, breaks = pw_fit(input_func, 50000, 8, 1.8, 1.3)

# The piecewise function itself
def pw_calc(x, slopes, plateaus, breaks):
    x = x.astype('float')
    
    cond_list = [x < breaks[0]]
    
    for idx, j in enumerate(breaks[:-1]):
        cond_list.append((j <= x) & (x < breaks[idx+1]))
        
    cond_list.append(breaks[-1] <= x)
    
    func_list = [lambda x: x * slopes[0]]
    
    for idx, j in enumerate(slopes[1:]):
        func_list.append(plateaus[idx])
        func_list.append(lambda x, j=j: x * j)
          
    return np.piecewise(x, cond_list, func_list)

y_output = pw_calc(x, slo, plat, breaks)

plt.plot(x,y,y_output);

from How to fit a piecewise (alternating linear and constant segments) function to a parabolic function?