add AutoPiecewiseLinFit class

[trueParadise]

Hi Charles
First Thanks for your contributions to this project!
I write a class to demonstrate my thoughts to approach how to determine the number of segments. Basically, I use polynomial to fit the graph first and then find the roots, then use those roots with your custom deviations to define bounds for speeding up convergence. However, I think convergence can be further speeded up by restricting breakpoints to int number. But I didn't finish this part so far.
Looking for your suggestion about my thoughts. you can email me ([email protected])

King regards,
Yucheng

[cjekel]

Is the purpose of lb and ub to give some 'wiggle room' or 'slack' where break points are specified?

I think you need to make lb ==ub and then default this to be some percentage of the total distance of x.

[trueParadise]

yes, it is 'wiggle room' or 'slack' where breakpoints are specified. I agree with you to set lb==ub = some percentage of the total distance of x

[cjekel]

Hi Yucheng,

Thank you for this ambitious software contribution! I think there are many users that will find this really helpful. This is a functionality that has been requested dozens of times.

This would close #17 and #88

The principle of restricting breakpoints to their own much smaller bounds will definitely make the optimization cheaper. This reduces the design space of possibilities by a lot.

I think this PR is going to take a bit of work on both of us to get this in. Are you in a position where you can put more time into this? If so I think this will be an excellent contribution that will really help a lot of people.

Questions/comments about your method

• Do you take all of the imaginary roots to seed your breakpoints? or just a subset?
• I made a comment about the root +- lb. I think this would be better as a percentage of x. Can you just comment on the purpose of this?
• What happens if you use a smaller degree polynomial?
• I think you may want to perform a search to find some optimal degree polynomial
• What are other strategies you have explored to select 'bounds' on breakpoints? Did this method work the best?
• Since you have reduced the design space a lot, we might want to just run a single local optimizer with your new bounds. This should be faster than the .fit() call that uses differential evolution.

Overall design of pwlf automatic methods

It's going to seem that I'm very picky here. Please feel free to disagree.

Automatic pwlf methods should feel like a natural extension to pwlf. They should be simple and a user should be able to get the result of an automatic fit with one function call.

I think we should use auto.py as the placeholder for automatic methods. I think we should prepare that there may be more than one automatic method in the future.

In terms of what goes into auto.py. I think we really need to preserve a lot of the functionality that is already in PiecewiseLinFit that many users are familiar with. We have a couple choices here, as I can see a future where we have more than one automatic method. I want your feedback here, and feedback from users.

• We could do class inheritance to preserve all functionality within PiecewiseLinFit for automatic methods.
• We could just make single methods auto.py that will return an instance of PiecewiseLinFit that is automatically fit.
• We could add automatic fitting methods into the PiecewiseLinFit as PiecewiseLinFit.auto(method="polynomial_breakpoint_reduction")

On working with integers as breakpoints

I've seen a number of people use these models with integers only. It turns out that integer optimization is not necessarily faster than floating point optimization.

The ways I see to solve integer optimization problem:

• Compute the total Combinatorial possibilities. If it is possible to solve all of these possibilities, do so and you will have your solution.
• If the above is not possible, try some discrete optimization algorithm (like a genetic algorithm) for a cost that you can afford. This will give you the best approximate solution.

Thanks again for this! I'm sure people will find it useful in it's current form. However, with a few modifications I think it can become even more powerful!

Thanks,
CJ

[trueParadise]

Hi Charles,
Thank you for your acknowledgment and confirmation and I am glad that it could help.
Yes, I agree with you that the current form is not a perfect one. There is still some room for improvement, so I am happy to put more time into it and make it better.

Here are points you've mention

• Question about imaginary roots. Keep imaginary seem unreasonable. There are two reasons.

  1. a = np.array([1+100000j, 3+0j, 5+6j])

r = list(filter(lambda x : x < 5 and x > 1, a))

r=[(1+100000j), (3+0j)]

     I am not sure how the filter works for comparing a real number and a complex number. I think they are incomparable directly. However, I am not sure what imaginary roots mean for polynomial. Can we show those roots on the Cartesian coordinate system? Do you know about that?

  2. if roots include imaginary, the results seem to get worst and include some unreasonable breakpoints.
     with imaginary roots(lb=0, ub=1) | without imaginary roots (lb=0, ub=1)
     :-------------------------:|:-------------------------:
     |

• Question about lb and ub
  Yes, lb and ub give some 'wiggle room' or 'slack' where breakpoints are specified. I agree with you. LB and UB should be the same. LB == UB == some ratio of x However, I have no idea how to find good a ratio, do you have any idea?
• Polynomial degree
  When I use a low Polynomial degree, the ability of the Polynomial model to fit will decline. For my dataset, I think the degree should be greater than 14 or 16 maybe, but I did not tweak this number. I think there are two ways to do this
  1. Find pieces you are not happy with and use Polynomial to do sub-fitting. and find breakpoints in this interval. I suggest that the line should contain enough data(like 10% of all your dataset or else). if the line contains not have enough data, it may lead your model too complicated.
     
     
  2. Increase the polynomial degree and find the elbow point of the error curve. I think it should work. Do you have any other method to find the optimal Polynomial degree?
• Speed
  I think the speed of differential evolution with bound is ok, but not very fast. One good thing is that differential evolution can run parallel. Maybe there is another optimizer that is more suitable for this situation, but I am not an expert in Optimization. Maybe you have a good idea! But I think parallelism is very important, new methods should support parallel.
• Integer breakpoints
  I think solving all integer breakpoints combinations is expensive, but we can parallel this process. I am not familiar with the Genetic algorithm, could you explain a little be why you call it a discrete optimization algorithm?
• How to implement
  I think this kind of form is good for me! very clear! PiecewiseLinFit.auto(method="polynomial_breakpoint_reduction"). But I think we need some tools like cal_slopes() (my version is a little different from yours) parse_model() etc. Basically, we need parameters for each line(slope, intercept) also with its interval [begin, end]. I think we can use {[begin, end]: function ...} data struct because that is easy to use.

Thanks for your work! hope we can make this function better.

King regards,
Yucheng

[cjekel]

Instead of my_pwlf.fit() try fit_guess(x_root, bounds=self.bounds). I suspect this will be good enough, and much faster than the even the parallel differential evolution!

[cjekel]

I don't have a lot of intuition on how the polynomial roots are being used, and I'm still a bit confused here. Can you do some debugging steps, and show me how the polynomials are being used. I want to see the logic step by step.

0. Show me the resulting polynomial fits to some sample data. In essence, plot the polynomial on the data. Also show your result auto fit on the same plot.

coeff = np.polyfit(self.x_data, self.y_data, degree)
f = np.poly1d(coeff)
x_hat = np.linspace(self.x_data.min(), self.x_data.max(), 1000)
f_hat = f(x_hat)
# please plot(x_hat, f_hat)

1. Show me the result of xroot = np.roots(f_d)
2. Show me the result after both filters are applied. So xroot after

xroot = list(filter(lambda x : x.imag == 0, xroot))
xroot = list(filter(lambda x : x >= self.x_data.min() and x <= self.x_data.max(), xroot))

I think this will help me give you hints onto how to select polynomial degree, and other parameters.

---

Can we show those roots on the Cartesian coordinate system? Do you know about that?

So consider the polynomial y=x**2 + 1. It will have two complex roots at [0.0 + i, 0.0 - i], which correspond to the function minimum at x=0.

Question about lb and ub
Yes, lb and ub give some 'wiggle room' or 'slack' where breakpoints are specified. I agree with you. LB and UB should be the same. LB == UB == some ratio of x However, I have no idea how to find good a ratio, do you have any idea?

I would try something like +- 5% of the total length. You can then make an option that the user can specify which percent to use. Try it out, if it doesn't work you can tweak the default a bit.

Speed
I think the speed of differential evolution with bound is ok, but not very fast. One good thing is that differential evolution can run parallel. Maybe there is another optimizer that is more suitable for this situation, but I am not an expert in Optimization. Maybe you have a good idea! But I think parallelism is very important, new methods should support parallel.

I put a comment into the code, but I think you can get away with using fit_guess() instead of fit(). Can you give this a try and report if it works? My rationale is that I think you constrain the problem enough with reduced breakpoints that you don't need the exhaustive global search that differential evolution performs.

Integer breakpoints
I think solving all integer breakpoints combinations is expensive, but we can parallel this process. I am not familiar with the Genetic algorithm, could you explain a little be why you call it a discrete optimization algorithm?

Let's hold off on this.

How pwlf works by default is a continuous optimization problem, so breakpoints can be any real number between x.min() and x.max(). The reason I said integer breakpoints is a discrete optimization problem is that you would restrict possible breakpoints to the set of integers from x.min() to x.max(). There would then be some discrete (finite) number of possible integer breakpoints combinations, while the continuous approach is infinite.

How to implement
I think this kind of form is good for me! very clear! PiecewiseLinFit.auto(method="polynomial_breakpoint_reduction"). But I think we need some tools like cal_slopes() (my version is a little different from yours) parse_model() etc. Basically, we need parameters for each line(slope, intercept) also with its interval [begin, end]. I think we can use {[begin, end]: function ...} data struct because that is easy to use.

I'm cool with this! Let's also hold off on this, and do some tweaks to the method first. I'm pretty sure once we can quickly refactor to get this to elegantly implement into pwlf.

[trueParadise]

Hi Charles,
I added some pictures to show how degrees of polynomial work.
Moreover, there are some comments in autopwlf.py file.
It turns out that small polynomial degrees can not find extreme points precisely.

[trueParadise]

Sorry, there was some problem with my last commit.
this is a new commit.

[trueParadise]

Hi Charles,
Do you have any progress on this section?

[cjekel]

I just wanted to reach out and say I have not forgot about this. Getting this in is something that I think will improve this library!

I have just been completely swamped to help out here.

I should have more free time in the next couple months (at least I hope so!).

[ffffxue]

just want to tag along this PR if it's not too late.
Loved the idea of using the solutions to the derivatives for a high-degree polynomial fit. The problem I have found with it tho, is that it does a good job when I have a segmented line that the sign of the slope flipped from one segment to another (e.g. at one point in time, the slope changes from -1 to 1). However, for example, the method seems to struggle when there is an increasing slope or a decresing slope from one to another segment (e.g. at one point in time, the slope changes from 1 to 2).

Wondering what are your thoughts?

[cjekel]

just want to tag along this PR if it's not too late. Loved the idea of using the solutions to the derivatives for a high-degree polynomial fit. The problem I have found with it tho, is that it does a good job when I have a segmented line that the sign of the slope flipped from one segment to another (e.g. at one point in time, the slope changes from -1 to 1). However, for example, the method seems to struggle when there is an increasing slope or a decresing slope from one to another segment (e.g. at one point in time, the slope changes from 1 to 2).

Wondering what are your thoughts?

This is good feedback. Thank you @ffffxue

I think if i had a go-to way of auto fitting these models I would do it based on the variance in your data. It would shoot for a particular signal to noise ratio, but it's just in my head and not implemented yet.