M3/example #3

mauro3 · 2016-11-15T20:31:51Z

A simple example on how to implement the interface. This is based on top of branch m3/fns (PR #2)

ChrisRackauckas · 2016-11-15T20:32:47Z

examples/euler-algs.jl

+sol2 = solve(p, BwdEulerAlg, tsteps=0:dt:1)
+
+using Plots
+plot(sol.t, sol.u)


why not use the plot recipe? This way of plotting doesn't scale to vectors of vectors and is more verbose.

ChrisRackauckas · 2016-11-15T20:33:09Z

examples/euler-algs.jl

+
+function solve{uType,tType,isinplace}(p::ODEProblem{uType,tType,isinplace},
+                                      Alg::Type{BwdEulerAlg};
+                                      tsteps=error("provide tsteps"),


I thought we agreed on tsteps instead of tstops? At least that's what I put in the docs and can change it.

ChrisRackauckas · 2016-11-15T20:33:32Z

examples/euler-algs.jl

+
+ana = t -> exp(u)
+
+function solve{uType,tType,isinplace}(p::ODEProblem{uType,tType,isinplace},


There's no reason to restrict solve to a concrete type. Have it dispatch on AbstractODEProblem{uType,tType,isinplace} so it'll work for whatever ODE problem it's given (including test problems).

ChrisRackauckas · 2016-11-15T20:33:50Z

examples/euler-algs.jl

+using DiffEqBase
+
+# make a algorithm type
+abstract EulerAlgs <: OrdinaryDiffEqAlgorithm


Also, the algorithm should subtype AbstractODEAlgorithm. See how each package has an abstract.

https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/algorithms/ode_algorithms.jl#L130

ChrisRackauckas · 2016-11-15T20:35:24Z

examples/euler-algs.jl

+ff(t,u) = u
+ff(::Val{:jac}, t, u) = 1
+p = ODEProblem(ff, 1.0, (0.0,1.0))
+sol2 = solve(p, BwdEulerAlg, tsteps=0:dt:1)


Also, the example should show ODETestProblem, plot(sol,plotanalytic=true), and test_convergence. For the last:

dts = 1./2.^(8:-1:4) sim = test_convergence(dts,prob,BwdEulerAlg) plot(sim) sim.𝒪est[:final]

on a test problem (linear ODE should be fine). Check out OrdinaryDiffEq's tests + DiffEqProblemLibrary if you need help. If this is difficult, let me know and I'll improve the devdocs (though most of it should already be there? Just at a high level without an example).

ChrisRackauckas · 2016-11-15T20:36:27Z

This PR also has the function interface specs: accident?

mauro3 · 2016-11-15T20:38:48Z

I did base it on top of m3/fns because I used has_jac. But I can change this.

ChrisRackauckas · 2016-11-15T20:40:32Z

I did base it on top of m3/fns because I used has_jac

Makes sense.

mauro3 · 2016-11-15T20:42:25Z

examples/euler-algs.jl

+
+
+using Plots
+plot(sol, plot_analytic=true)


this gives me an error:

julia> include("euler-algs.jl") ERROR: LoadError: type ODETestSolution has no field timeseries_analytic in macro expansion at /home/mauro/.julia/v0.5/DiffEqBase/src/solutions/solution_interface.jl:62 [inlined]

What do I need to do?

That's my bad. That's a bug in the plot recipe since I used to call the outputs timeseries and timseries_analytic (but liked just using u / y like PR49). It should work just by changing every timseries_analytic to u_analytic in the plot recipe. I can do that right now.

Should be fixed on master now. I wish there was a good way to include the plot recipes in tests without having to make Travis install a whole plotting package.

Now I get:

ERROR: LoadError: BoundsError in getindex(::Float64, ::Int64) at ./number.jl:21 in macro expansion at /home/mauro/.julia/v0.5/DiffEqBase/src/solutions/solution_interface.jl:84 [inlined] in apply_recipe(::Dict{Symbol,Any}, ::DiffEqBase.ODETestSolution{Array{Float64,1},Array{Float64,1},Float64,FloatRange{Float64},Array{Any,1},DiffEqBase.ODETestProblem{Float64,Float64,false,##26#28},FwdEulerAlg}) at /home/mauro/.julia/v0.5/RecipesBase/src/RecipesBase.jl:238 in _process_userrecipes(::Plots.Plot{Plots.PyPlotBackend}, ::Dict{Symbol,Any}, ::Tuple{DiffEqBase.ODETestSolution{Array{Float64,1},Array{Float64,1},Float64,FloatRange{Float64},Array{Any,1},DiffEqBase.ODETestProblem{Float64,Float64,false,##26#28},FwdEulerAlg}}) at /home/mauro/.julia/v0.5/Plots/src/pipeline.jl:73 in _plot!(::Plots.Plot{Plots.PyPlotBackend}, ::Dict{Symbol,Any}, ::Tuple{DiffEqBase.ODETestSolution{Array{Float64,1},Array{Float64,1},Float64,FloatRange{Float64},Array{Any,1},DiffEqBase.ODETestProblem{Float64,Float64,false,##26#28},FwdEulerAlg}}) at /home/mauro/.julia/v0.5/Plots/src/plot.jl:173 in (::Plots.#kw##plot)(::Array{Any,1}, ::Plots.#plot, ::DiffEqBase.ODETestSolution{Array{Float64,1},Array{Float64,1},Float64,FloatRange{Float64},Array{Any,1},DiffEqBase.ODETestProblem{Float64,Float64,false,##26#28},FwdEulerAlg}) at ./<missing>:0 in include_from_node1(::String) at ./loading.jl:488 while loading /home/mauro/.julia/v0.5/DiffEqBase/examples/euler-algs.jl, in expression starting on line 89

Try it now. I just tried it out and plot_analytic it should be all compatible with the interface with dense plotting but true and false.

mauro3 · 2016-11-15T21:02:48Z

Also when running at the end of the script:

using DiffEqDevTools
dts = 1./2.^(8:-1:4)
sim = test_convergence(dts,p2,BwdEulerAlg)
plot(sim)
sim.𝒪est[:final]

(using master of DiffEqDevTools) I get:

julia> sim = test_convergence(dts,p2,BwdEulerAlg)
ERROR: MethodError: Cannot `convert` an object of type Type{BwdEulerAlg} to an object of type Array{Float64,1}
This may have arisen from a call to the constructor Array{Float64,1}(...),
since type constructors fall back to convert methods.
 in #solve#51(::Float64, ::Bool, ::Int64, ::DiffEqBase.ExplicitRKTableau, ::Bool, ::Void, ::Symbol, ::Bool, ::Bool, ::Array{Float64,1}, ::Array{Float64,1}, ::Bool, ::Float64, ::Rational{Int64}, ::Rational{Int64}, ::Void, ::Void, ::Rational{Int64}, ::Bool, ::Void, ::Void, ::Int64, ::Float64, ::Float64, ::Bool, ::OrdinaryDiffEq.#ODE_DEFAULT_NORM, ::OrdinaryDiffEq.#ODE_DEFAULT_ISOUTOFDOMAIN, ::Bool, ::Int64, ::String, ::Void, ::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.ODETestProblem{Float64,Float64,true,#ff}, ::Type{DiffEqBase.Tsit5}, ::Type{T}, ::Array{Any,1}, ::Array{Any,1}) at /home/mauro/.julia/v0.5/OrdinaryDiffEq/src/solve.jl:155
 in (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.ODETestProblem{Float64,Float64,true,#ff}, ::Type{DiffEqBase.Tsit5}, ::DataType, ::Array{Any,1}, ::Array{Any,1}) at ./<missing>:0 (repeats 2 times)
 in #solve#99(::Bool, ::Array{Any,1}, ::Function, ::DiffEqBase.ODETestProblem{Float64,Float64,true,#ff}, ::Type{T}, ::Vararg{Type{T},N}) at /home/mauro/.julia/v0.5/DiffEqBase/src/DiffEqBase.jl:62
 in (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.ODETestProblem{Float64,Float64,true,#ff}, ::Type{T}) at ./<missing>:0
 in (::DiffEqDevTools.##20#21{Bool,Bool,Array{Any,1},Array{Float64,1},DiffEqBase.ODETestProblem{Float64,Float64,true,#ff},DataType})(::Int64) at ./<missing>:0
 in collect(::Base.Generator{UnitRange{Int64},DiffEqDevTools.##20#21{Bool,Bool,Array{Any,1},Array{Float64,1},DiffEqBase.ODETestProblem{Float64,Float64,true,#ff},DataType}}) at ./array.jl:307
 in #test_convergence#19(::Bool, ::Bool, ::Array{Any,1}, ::Function, ::Array{Float64,1}, ::DiffEqBase.ODETestProblem{Float64,Float64,true,#ff}, ::Type{T}) at /home/mauro/.julia/v0.5/DiffEqDevTools/src/convergence.jl:41
 in test_convergence(::Array{Float64,1}, ::DiffEqBase.ODETestProblem{Float64,Float64,true,#ff}, ::Type{T}) at /home/mauro/.julia/v0.5/DiffEqDevTools/src/convergence.jl:40

This seems to call into OrdinaryDiffEq, which probably shouldn't happen.

ChrisRackauckas · 2016-11-15T21:05:59Z

This seems to call into OrdinaryDiffEq, which probably shouldn't happen.

That's because you subtyped OrdinaryDiffEqAlgorithm, and your specific algorithm doesn't exist, so it called my solve command. You just ran into this: SciML/Roadmap#9

In fact, this is exactly the example I gave. The solution that I know of is you have to import this algorithm in DiffEqDevTools, otherwise it won't know about your solve when it calls solve. That's the problem that I don't know how to solve: all of the higher level analysis tools (parameter estimation, sensitivity analysis, dev tools) will have to depend on each solver it's compatible with... is that okay? (but if you take the same convergence analysis code and copy paste it in the REPL, it's fine)

That fact might solve itself with extern modules.

mauro3 · 2016-11-15T21:15:22Z

Ah, maybe now I'll be able to follow that issue! I'll see again. Edit: see #5

Edit: but this also happens now that I subtype DiffEqBase.AbstractODEAlgorithm

[ci skip]

mauro3 · 2016-11-16T12:58:30Z

examples/euler-algs.jl

+immutable BwdEulerAlg <: EulerAlgs
+end
+
+function DiffEqBase.solve{uType,tType,isinplace}(p::AbstractODEProblem{uType,tType,isinplace},


I didn't extend DiffEqBase.solve but created a new function, because solve was not exported (#4) it didn't warn me either. Thus, below in the call to test_convergence another solve was called.

oh, the whole problem there was that solve wasn't exported? For some reason I thought it was...

The not exported meant that there was no warning, it would have worked otherwise.

ChrisRackauckas

This looks fine but #2 goes first.

ChrisRackauckas · 2016-11-16T16:36:34Z

Actually, shouldn't the example just go in devdocs? Seems weird to just put a random example here. Also that would give more room for explanation, and a tutorial would really help make it clear.

mauro3 · 2016-11-16T19:20:19Z

DevDocs would be fine too. I just wanted to test-drive the interface and had to put it somewhere for discussion. So maybe leave it in this PR until finished discussing and then move to the docs?

ChrisRackauckas · 2016-11-16T19:22:37Z

Sounds good

ChrisRackauckas · 2017-05-18T18:00:00Z

Would you mind turning this code into a ExampleSolverDiffEq.jl? I think having a simple (unregistered) repo for the devdocs to point to that is a working example for how to do this is a good idea. I'd swap this code over myself but I'm not sure how to do that and still give you the credit.

mauro3 · 2017-08-08T19:36:35Z

Ha, this was still valid! Good.

ChrisRackauckas · 2017-08-08T19:37:48Z

Yeah, as much as it seems like code changes around here often, the core actually hasn't changed much at all in a year. I guess this was a good proof of that haha.

ChrisRackauckas requested changes Nov 15, 2016

View reviewed changes

mauro3 commented Nov 15, 2016

View reviewed changes

mauro3 added 4 commits November 16, 2016 13:49

Added example

8c2288e

some updates

884361c

tstops

0875f90

fixed test_convergence

457a1ad

[ci skip]

mauro3 force-pushed the m3/example branch from d76b526 to 457a1ad Compare November 16, 2016 12:55

mauro3 commented Nov 16, 2016

View reviewed changes

ChrisRackauckas approved these changes Nov 16, 2016

View reviewed changes

ChrisRackauckas merged commit 457a1ad into master Aug 8, 2017

mauro3 deleted the m3/example branch August 8, 2017 19:36

DaniGlez mentioned this pull request Jul 19, 2023

switch callback rootfinding algorithm from InternalFalsi to InternalITP to Improve robustness #917

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

M3/example #3

M3/example #3

mauro3 commented Nov 15, 2016

ChrisRackauckas Nov 15, 2016

ChrisRackauckas Nov 15, 2016

ChrisRackauckas Nov 15, 2016

ChrisRackauckas Nov 15, 2016

ChrisRackauckas Nov 15, 2016

ChrisRackauckas commented Nov 15, 2016

mauro3 commented Nov 15, 2016

ChrisRackauckas commented Nov 15, 2016

mauro3 Nov 15, 2016

ChrisRackauckas Nov 15, 2016

ChrisRackauckas Nov 15, 2016

mauro3 Nov 15, 2016

ChrisRackauckas Nov 15, 2016

mauro3 commented Nov 15, 2016

ChrisRackauckas commented Nov 15, 2016 •

edited

Loading

mauro3 commented Nov 15, 2016 •

edited

Loading

mauro3 Nov 16, 2016

ChrisRackauckas Nov 16, 2016

mauro3 Nov 16, 2016

ChrisRackauckas left a comment

ChrisRackauckas commented Nov 16, 2016 •

edited

Loading

mauro3 commented Nov 16, 2016

ChrisRackauckas commented Nov 16, 2016

ChrisRackauckas commented May 18, 2017

mauro3 commented Aug 8, 2017

ChrisRackauckas commented Aug 8, 2017


		ana = t -> exp(u)

		function solve{uType,tType,isinplace}(p::ODEProblem{uType,tType,isinplace},

M3/example #3

M3/example #3

Conversation

mauro3 commented Nov 15, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ChrisRackauckas commented Nov 15, 2016

mauro3 commented Nov 15, 2016

ChrisRackauckas commented Nov 15, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mauro3 commented Nov 15, 2016

ChrisRackauckas commented Nov 15, 2016 • edited Loading

mauro3 commented Nov 15, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ChrisRackauckas left a comment

Choose a reason for hiding this comment

ChrisRackauckas commented Nov 16, 2016 • edited Loading

mauro3 commented Nov 16, 2016

ChrisRackauckas commented Nov 16, 2016

ChrisRackauckas commented May 18, 2017

mauro3 commented Aug 8, 2017

ChrisRackauckas commented Aug 8, 2017

ChrisRackauckas commented Nov 15, 2016 •

edited

Loading

mauro3 commented Nov 15, 2016 •

edited

Loading

ChrisRackauckas commented Nov 16, 2016 •

edited

Loading