Can pMax and pMin be used as hard constraints?

[liuyxpp]

I have a problem that p is strictly constrained in a range, say [a, b], (which can be computed). When p is out of this range, the evaluation of the function will fail. In my case, a negative will be given to the log function. Currently, when I set pMin=a+eps(), and pMax=b-eps(), the parameter prediction step will still predict a value outside the range [a, b], which leads to the failing of the whole continuation process.

Sample output is

      From worker 2:    Continuation Step 339
      From worker 2:    Step size = 2.0000e-03
      From worker 2:    Parameter Aα = 8.2503e-01 ⟶  8.2741e-01 [guess]
      From worker 2:    --> Step Converged in 2 Nonlinear Iteration(s)
      From worker 2:    Parameter Aα = 8.2503e-01 ⟶  8.2741e-01
      From worker 2:    Predictor: BifurcationKit.SecantPred()
      From worker 2:    --> Computed 3 eigenvalues in 1 iterations, #unstable = 3
      From worker 2:    ──────────────────────────────────────────────────────────────────────
      From worker 2:    Continuation Step 340
      From worker 2:    Step size = 2.0000e-03
      From worker 2:    Parameter Aα = 8.2741e-01 ⟶  8.2979e-01 [guess]

For this problem, my b is 0.829285. As can be seen, the prediction of p in Step 340 is 0.82979 (>0.829285), which leads to the breaking of my continuation process.

I wonder if we can make the [pMin, pMax] constraint hard, so that p never goes outside of it. Or even better, BK can detect the failing of evaluation of the function when the prediction of p is outside the range, and reduce the step size accordingly to force the prediciton of p lies in the range.

[rveltz]

Strange, I would say that it works this way already, see here. Can you please confirm that you use the last version?

[liuyxpp]

Below is my Pluto environment

      From worker 2:          Status `/private/var/folders/49/vv1vdsyj4_9g74lfystrkb2w0000gn/T/jl_0Jxr0I/Project.toml`
      From worker 2:      [0f109fa4] BifurcationKit v0.1.7
      From worker 2:      [f6369f11] ForwardDiff v0.10.21
      From worker 2:      [d1acc4aa] IntervalArithmetic v0.20.0
      From worker 2:      [d2bf35a9] IntervalRootFinding v0.5.10
      From worker 2:      [b964fa9f] LaTeXStrings v1.2.1
      From worker 2:      [2774e3e8] NLsolve v4.5.1
      From worker 2:      [d96e819e] Parameters v0.12.3
      From worker 2:      [91a5bcdd] Plots v1.22.6
      From worker 2:      [7f904dfe] PlutoUI v0.7.16
      From worker 2:      [f2b01f46] Roots v1.3.5
      From worker 2:      [efcf1570] Setfield v0.8.0
      From worker 2:      [90137ffa] StaticArrays v1.2.13
      From worker 2:      [1f5e811d] TernaryPlots v0.1.0
      From worker 2:      [37e2e46d] LinearAlgebra
      From worker 2:      [44cfe95a] Pkg

FYI, here is my code related to BK

function continuation_binodal(model, initial, ϕAα,
                             ϕAα_min=0.0, ϕAα_max=1.0,
                             opts=BK.ContinuationPar())
    F(x, p) = collect(binodal_condition(p.Aα, x..., p.model))
    J(x, p) = ForwardDiff.jacobian(z->F(z,p), x)
    par = (Aα=ϕAα, model=model)
    x0 = zeros(3)
    x0 .= initial # initial can be Vector or Tuple, but not NamedTuple

    opts = BK.ContinuationPar(opts,
                            dsmax=0.02, dsmin=1e-5, ds=0.001, maxSteps=2000,
                            pMin=ϕAα_min,
                            pMax=ϕAα_max,
                            saveSolEveryStep=1,
                            maxBisectionSteps=50,
                            detectBifurcation=3)

    br, = BK.continuation(F, J, x0, par, (@lens _.Aα), opts,
            recordFromSolution = (x,p) -> (ϕAα=p, ϕBα=x[1], ϕAβ=x[2], ϕBβ=x[3]),
            bothside=false, verbosity=3)
    return br
end

[rveltz]

OK, I need to see this more in detail. In the mean time, you can perhaps clamp the value of p inside F before it goes into binodal_condition

[rveltz]

This minimum working example shows how it works. The F function is not allows to step outside [0,1]. You can see that the following does not error.

function continuation_binodal2(opts=BK.ContinuationPar())
    F(x, p) = (@assert 0 <= p <= 1 "Not allowed!"; return (-(p + 0.1) .+ x.^2))
    J(x, p) = ForwardDiff.jacobian(z->F(z,p), x)
    par = 0.5
    x0 = [sqrt(par)]

    opts = BK.ContinuationPar(opts,
                            dsmax=0.02, dsmin=1e-5, ds=-0.001, maxSteps=2000,
                            pMin=0.,
                            pMax=1.,
                            saveSolEveryStep=1,
                            maxBisectionSteps=50,
                            detectBifurcation=3)

    br, = BK.continuation(F, J, x0, par, (@lens _), opts,
            recordFromSolution = (x,p) -> (x[1]),
            bothside=false, verbosity=3, plot = true)
    return br
end

Perhaps you can provide a working example with a bug ?

[rveltz]

should we close this?