rachford_rice about chemicals HOT 14 CLOSED

Hi Caleb,

I'm pretty happy about the way flexsolve is working now. I implemented the switch from secant/aitken_secant to IQ_interpolation (instead of bisection). Now I have a parameter called "checkroot" which, if True, tells the solver to make sure both xtol and ytol are satisfied. I also have a new parameter called "checkiter" that tells the solver whether to raise an error if no solution could be found within given maximum number of iterations.

All new parameters may be a little confusing, but you can check out the documentation and it should be straight forward:
https://github.com/yoelcortes/flexsolve

I also made some classes that allow for easier testing of solvers. Just pip install flexsolve==0.3.24 and you can try them. I used your test cases for this. I may have gone overboard, but I think its pretty nice:
solver_tests.ipynb.txt

If you increase the number of iterations, flexsolve's secant/aitken_secant achieves 67 passed cases.

Thanks for all you help on this, I really appreciate it.
I'll move on with making progress on the chemicals package.

from chemicals.

BTW, I added pytests for all solvers, including jit compiled solvers:
https://github.com/yoelcortes/flexsolve/blob/master/tests/test_scalar_solvers.py

All of flexsolve's jitted solvers work just as well as their non-jitted versions (according to the tests).
Let me know if you have any comments or questions!

from chemicals.

Hi Yoel,
I think your new flexsolve tests are really great. I am thrilled I could help. A numerical solver isn't something I have really gotten excited about before, but I am really happy you enjoy working on them.

I think all the options you have added are good, and really cool. Some, like checkroot I have a hard time imagining myself using but some of your points like passing the last point into another solver could be useful in the right circumstances. I have done something a little similar - to try to detect discontinuities and move to trying a whole different equation set (two-phase EOS vs. one phase EOS).

If you really want to spend lots of time on the topic, you could think about publishing flexsolve, or seeing if there are others who would like to try to get the most optimized solver possible and get a publication and then get the algorithm all over the place :) I know the Julia community has been pretty interested in numerical solvers, there is : https://github.com/JuliaMath/Roots.jl and https://github.com/JuliaNLSolvers/NLsolve.jl ; and scipy also has toms748 https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.toms748.html#scipy.optimize.toms748 which I have confirmed to be better in some cases in terms of number of iterations (but the implementation is actually quite slow - minimum solve time is ~1 ms).

I haven't found time to get any code out lately, but I hope I do soon!

from chemicals.

https://github.com/CalebBell/chemicals/blob/master/tests/test_rachford_rice.py#L357 is the analytical failure of rachford rice. By the way, sometimes you edit your comments and add a question but I don't get an email notification so I don't see them. Please feel free to comment multiple times in a row to avoid this.

from chemicals.

Hi Caleb, excellent work on the Rachford Rice solutions! I never thought about having a minimum or maximum vapor fraction so that it can solve. I've actually had trouble with this before. I'm not too sure if the solution I've found to this problem is even correct. I'll look into implementing what you've made.

I've put some work on flexsolve recently and I think its just about ready for other users. I would like to ask you to check it out and give me your comments. There are considerable differences to scipy's solvers and some aspects that I consider even better... Perhaps you may find some inspiration, or tryout flexsolve eventually. Here are the mayor differences:

• Gives you the option to pass the upper and lower bound of the objective functions, saving two objective function calls. This is especially relevant when the objective function call has an inner loop. For example, solving for things like temperature at a given vapor fraction and pressure when you know the dew point and bubble point. I think this would be useful in Rashford Rice (although not as significant)
• Termination is based on absolute tolerance of the objective function (i.e. f(x) < ytol). I prefer this because usually when we solve a root, its easier to think of the absolute tolerance we would like. But also having an rtol and using numpy's allclose is something I'm planning on adding in the near future.
• Optionally raise an error if solver could not converge (if checkroot=True). It is useful to set checkroot=False when you would like to pass the last guess in another solver, or if you are willing to accept a solution under a smaller tolerance than originally provided.

Some aspects that I think you will like:

• No subclassed exceptions
• Compatible with numba
• Fast and conventional solvers (Aitken, Wegstein, *Inverse quadratic interpolation)
• Accelerated fixed-point methods have exception handling to try normal fixed point iteration when accelerated guess doesn't work (due to infeasible values) and then jumps back to the accelerated algorithm.

Flexsolve's inverse quadratic interpolation solver (IQ_interpolation), is different from scipy's brentq in the following ways:

• Interpolation is always used as long as the guess is within the bounds. Otherwise, false position is used.
• When the guess is getting stuck, the next guess is (x_lb + x_guess + x_ub) / 3. In brentq, the next guess is bisection (i.e. (x_lb + x_ub) / 2).
• When false position cannot be used due to local min or max, bisection is used (this is the same in brentq).
• If checkbounds=False when the root does not lie within the given bracket, it tries to solve it anyway and returns the value that minimizes the objective function error out of all attempted guesses (including the brackets). I don't really like this behavior in general, but it is useful when you don't get the bounds right due to small numerical errors and the solution lies slightly outside the bounds. Also, sometimes the solution was skipped due to the numerical accuracy of inner solvers and nonlinear conditions (from activity coefficients)... choosing the best value that has been evaluated is an easy fix I came up with.
• You can pass an initial guess

Additional notes:

• I've found that flexsolve's IQ_interpolation is consistently faster that brentq and brenth for VLE in thermosteam by requiring less iterations.
• Since the last time you've had a look at flexsolve, I've deprecated some creative solvers I had fun making (e.g. a secant with Wegstein acceleration) as well as a yval parameter for root offset. I removed this to make the package as conventional as possible and easier to manage.

from chemicals.

Hi Yoel,

That's really cool! Diversity in using numerical solvers is great. While there are lots of works about it, I am sure the best algorithms are not yet known.
I could add some comments about my own biases for which solvers are best and worst. But in the end it is the results that matter.

Here is a Jupyter notebook you may find interesting, I developed it when I was playing around with a few things. I am sure I have not made a fair comparison against flexsolve. I hope you find it useful.

I will just add that while I find using an absolute tolerance to be very useful and do it quite a bit myself, setting a default value is tough - maybe someone is playing with numbers on the order of 1e-40, and then the solver will think it is done right away.

Univariate root finding test functions - flexsolve.ipynb.txt

I will try to give some more detailed feedback when I try it out more. It would also be cool to try and get a standard python "benchmark" for numerical solvers. For example SciPy has one for linear programming problems.

from chemicals.

Dope! Thanks a bunch for looking into this and helping me out with test functions. I also hope to have more tests (both pytest and doctests) in the near future.

I just updated my post to include this latest update for bounded solvers:
"If checkbounds=False when the root does not lie within the given bracket, it tries to solve it anyway and returns the value that minimizes the objective function error out of all attempted guesses (including the brackets). I don't really like this behavior in general, but it is useful when you don't get the bounds right due to small numerical errors and the solution lies slightly outside the bounds. Also, sometimes the solution was skipped due to the numerical accuracy of inner solvers and nonlinear conditions (from activity coefficients)... choosing the best value that has been evaluated is an easy fix I came up with."

Looking forward to your comments!

from chemicals.

I'll also go ahead and add an rtol value so that I can get these tests passing.
EDIT: Just realized that this would not affect the tests, so I probably won't add this until later.

from chemicals.

I made xtol an optional argument and fixed a few details in flexsolve. Also, both Wegstein and Aitken solve functions of the form f(x) = x, not f(x) = 0. I corrected the form in the file:

[removed old file and image]
When Wegstein or Aitken do not solve, sometimes its because of a floating point error or division by zero... I should find a way to prevent this. I'll check out the other cases individually as well and see if there is anything I can do.

Thanks so much for the tests, It helped make flexsolve better already.

from chemicals.

I introduced a bug on the last update to the accelerated fixed point solvers, here is the updated version (with the bugfix):

Univariate_root_finding_methods_flexsolve_tests.ipynb.txt

from chemicals.

Hi Yoel,
I didn't realize the wegstein and .aitken solvers worked that way. Those are some awesome numbers!

I think the only other trick you might not be aware of is to add in an optional bisection to the solvers - let the solver do what it wants, but if it has already tried two points that bracket the root and later it tries to jump outside of that range, do a bisection step instead, because the derivative approximation is probably misbehaving.

I really do agree that using a tolerance on something concrete, like enthalpy or temperature, is much better when you have a specific application in mind. When you don't, the xtol type check can be very useful for getting about as tight of an answer as you can. But as some of these tests show, some functions are pathological and will never show convergence in the xtol direction.

I have argument require_xtol in the secant in fluids. When it is False, the solver returns when only the y tolerance is satisfied. This prevents those pathological cases from being an issue. Using those two techniques and the following parameters {'xtol': 1e-14, 'bisection': True, 'require_xtol': False, 'ytol': 1e-12}, fluids.numerics.secant scores 71 passes. I think all the other cases fail because the derivative approximation is zero.

For me, the solver is secondary - I am concerned about the public APIs, and have no interest in trying to expose numerical methods to anyone. Keep up the good work!

On another topic, a surprising amount of my code uses numerical solvers, and it seems numba is not presently able to cache those. There is also a bug, when I turn on caching numba crashes using fluids.compressible.isothermal_gas with caching turned on. I am not sure how to proceed here. Have you had any success?

from chemicals.

Hui Caleb, thanks for the quick reply.

Wegstein and Aitken-Steffensen can work either way (f(x) = 0 or f(x) = x), I just used the forms that were originally published, but I believe Scipy uses the f(x)=0 form for Aitken.

I just updated flexsolve to use xtol=0. as its default value; I'll make a checkxtol argument at later time and just make sure xtol is not 0. when checkxtol is True (thanks for your suggestion!). I'll try out adding switching from secant/aitken-secant to IQ_interpolation to make the open solvers more robust (based on your idea to use bisection).

I believe I have a solution to your cache problem. The problem is that numba tries to compile code that it never gets to with a value of None. I haven't actually run anything from fluids; but I can tell by looking at the code for secant:

def secant(func, x0, args=(), maxiter=100, low=None, high=None, damping=1.0,
           xtol=1.48e-8, ytol=None, x1=None, require_eval=False, 
           f0=None, f1=None, bisection=False, same_tol=1.0, kwargs={},
           require_xtol=True):
    ...
    # Lines like this are problematic; ytol is compiled on first run as None
    # Then it tries to compile lt(float, none); even thought it didn't actually need to run it.
    if (ytol is not None and abs(q0) < ytol and not require_xtol) or q0 == 0.0: 
        return p0

I would suggest the following change (which is the solution used in flexsolve):

def secant(func, x0, args=(), maxiter=100, low=None, high=None, damping=1.0,
           xtol=1.48e-8, ytol=0., x1=None, require_eval=False, 
           f0=None, f1=None, bisection=False, same_tol=1.0, kwargs={},
           require_xtol=True):
    ...
    # Note that now ytol is compiled as a number
    if abs(q0) < ytol and not require_xtol or q0 == 0.0: 
        return p0

You'll have to be care of these None cases with other functions too. Here is a concrete test case that pinpoints this problem:

from numba import njit
import flexsolve as flx

@njit(cache=True)
def f(x, c):
    return 0.4*x**3 + 0.3*x + c

@njit(cache=True)
def solve(c=None):
    if not c: c = 2. # This tells numba its a number
    if c > 2: c = 2. # Now numba knows to do gt(float, float)
    return flx.aitken_secant(f, 0.1, args=(c,))

val = solve() # This works

@njit(cache=True)
def solve_fail(c=None):
    # Numba cannot compile c as a number on first run
    if not c and c > 2: c = 2. # tries to compile gt(none, float)
    return flx.aitken_secant(f, 0.1, args=(c,))

val = solve_fail() # This fails

from chemicals.

Thanks for the nice comments! I may look into publishing flexsolve, and possibly making a julia version of some of the solvers and making pull request with Julia's Roots package. As you've seen, numerical methods is pretty fun for me. But I'm not sure when I will get around to this, my priorities are really BioSTEAM related, so first I'll get to finishing up the chemicals package, publishing thermosteam, and then possibly publishing flexsolve.

BTW, I think I may have confused you with the old and new definitions for the checkroots parameter. Now checkroots=True checks for both tolerances (xtol and ytol) for termination, and if checkroots=False either tolerance will be used for termination. The new checkiter parameter decides whether or not to raise an runtime error if the solver could not terminate before the maximum number of iterations (this parameter is basically the old checkroots). Confusing changes, but now it should be stable jaja.

Thanks again, I'll get back to contributing to chemicals over the weekend too!

from chemicals.

Hi Yoel,
I am going to close this issue. I would appreciate it if you could see how it does vs. your own implementation with IQ_interpolate at some point too!

Sincerely,
Caleb

from chemicals.