Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change internals to Diffeq #191

Merged
merged 10 commits into from
Jan 12, 2018
Merged

Change internals to Diffeq #191

merged 10 commits into from
Jan 12, 2018

Conversation

ChrisRackauckas
Copy link
Contributor

This sets up the internals to use DiffEq. All of the tests pass. There is a breaking change in that the steady state parts now return the solution array (with saveat matching the giving timespan, so it's the same but includes the initial point). Other than that, nothing changes. This requires DiffEqCallbacks v0.7.0 which should be released shortly.

As for speed, benchmarks show that as setup it's roughly the same on very non-stiff problems. The previous integrator's reltol=1e-7 takes as much steps as this one's reltol=1e-6, and with these the DiffEq one is ever so slightly faster. But true benchmarks should take into account error as well. The QO benchmarks should be run again like in #92 . Note that the current DiffEq infrusturcture takes a hit on very cheap problems due to Julia's keyword argument handling and lack of constant propogation, but those are both fixed in v0.7 so we'll get a free speedup when the next version of Julia is released.

But, the DiffEq DP5 uses PI adaptive time stepping which will be more stable and by far more efficient on semi-stiff equations. And that's mostly the point. I wouldn't expect anything major on simple non-stiff equations like in the tests, but if you increase the stiffness of the PDE, then DiffEq allows a lot more options to handle stiff equations (whereas explicit Runge-Kutta methods simply fail). However, the option to pass integration algorithms to override the internal default is not part of this PR (but it would be very simple). Also, by changing ODEProblem to things like DDEProblem or SDEProblem this (with different algorithm choices) will solve other types of differential equations as proposed in #164 .

@ChrisRackauckas
Copy link
Contributor Author

Re-running tests because the release that this needed is now in METADATA. Also, I just noticed that this doesn't do anything special for the single-trajectory Monte Carlo stuff, so I am not sure why tests pass there. Those jumps need to be DiscreteCallbacks.

@coveralls
Copy link

Coverage Status

Coverage decreased (-0.3%) to 99.233% when pulling 8d306d3 on ChrisRackauckas:diffeq into 74420ff on qojulia:master.

@Datseris
Copy link
Collaborator

Yeah maybe a good idea to not test on "latest" :P

@david-pl
Copy link
Member

david-pl commented Jan 11, 2018

I ran the benchmarks. They look very promising.

bench-master

And even though, as you said, this is not doing anything special for Monte Carlo approaches, it is also not slowing things down.

bench-mcwf

In the above, the orange line is your PR.

There is a problem though, that I noticed while benchmarking. In the benchmarks, we use an fout.
When using this code (similar to the benchmarking code):

b = FockBasis(9)
a = destroy(b)
at = create(b)
n = number(b)

tspan = [0:1.0:10;]

H = a + at
J = [a]
rates = [1.0]

Ψ₀ = coherentstate(b, 0)

tout, ρt = timeevolution.master(tspan, Ψ₀, H, J; rates=rates)
exp_n = real.(expect(n, ρt))

the amount of points stored correspond to length(tspan). However, using an fout like so

exp_n = Float64[]
fout(t, ρ) = push!(exp_n, real(expect(n, ρ)))
timeevolution.master(tspan, Ψ₀, H, J; rates=rates, fout=fout)

results in an exp_n of length length(tspan)+1. As it turns out, the first save is pushed twice to exp_n, which shouldn't be there of course. Also, the call to master returns a non-void object even if the fout is defined, whereas previously it returned nothing.

The build checks keep failing, because we still have support for Julia v0.5, for which your latest Callbacks is not available. We should probably drop support for v0.5 anyway.

Edit:
The call to fout when specifying the type for SavedValues in this line

out = DiffEqCallbacks.SavedValues(Float64,typeof(fout(tspan[1], state)))

in combination with the in-place push! used in the above example is the cause of the problem I think.

@ChrisRackauckas
Copy link
Contributor Author

The build checks keep failing, because we still have support for Julia v0.5, for which your latest Callbacks is not available. We should probably drop support for v0.5 anyway.

Support for v0.5 is locked in METADATA. No packages can release new versions with v0.5. I'll go ahead and upgrade the build scripts.

@ChrisRackauckas
Copy link
Contributor Author

in combination with the in-place push! used in the above example is the cause of the problem I think.

Oh, I forgot to mention that in the new version there is a breaking change that you no longer have to build arrays to push! into, you just return the values in fout and they will be the solution values given by the function. I could change this but it seemed much more intuitive IMO. See the way the tests are setup.

@david-pl
Copy link
Member

Oh, I could have caught that...
Hm, I agree that, e.g. like in one of the test,

function fout(t, rho)
    expect(op1, rho)
end

is more intuitive than using push!.

Properly inferring the type of fout (as you say in the TODO comment) would avoid this breaking though right?

@ChrisRackauckas
Copy link
Contributor Author

Properly inferring the type of fout (as you say in the TODO comment) would avoid this breaking though right?

Yes, it would be the same to the user, but I would just use Base.core.inference calls to know the output type instead of actually computing it once. Let me go ahead and add that right now. I didn't want to do it before I got everything working because it's a little bit tedious, but now that everything is together it makes sense to handle it.

@codecov
Copy link

codecov bot commented Jan 11, 2018

Codecov Report

Merging #191 into master will decrease coverage by 0.49%.
The diff coverage is 100%.

Impacted file tree graph

@@           Coverage Diff           @@
##           master   #191     +/-   ##
=======================================
- Coverage    99.5%    99%   -0.5%     
=======================================
  Files          32     32             
  Lines        2222   2215      -7     
=======================================
- Hits         2211   2193     -18     
- Misses         11     22     +11
Impacted Files Coverage Δ
src/steadystate.jl 93.33% <100%> (-6.67%) ⬇️
src/timeevolution_base.jl 100% <100%> (ø) ⬆️
src/mcwf.jl 100% <100%> (ø) ⬆️
src/timecorrelations.jl 100% <100%> (ø) ⬆️
src/ode_dopri.jl 94% <0%> (-4.67%) ⬇️
src/subspace.jl 88.23% <0%> (-2.95%) ⬇️
src/bases.jl 97.22% <0%> (-2.78%) ⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 74420ff...b5c73c8. Read the comment docs.

@coveralls
Copy link

Coverage Status

Coverage decreased (-0.3%) to 99.233% when pulling ea937ce on ChrisRackauckas:diffeq into 74420ff on qojulia:master.

@coveralls
Copy link

Coverage Status

Coverage decreased (-0.3%) to 99.233% when pulling ea937ce on ChrisRackauckas:diffeq into 74420ff on qojulia:master.

@ChrisRackauckas
Copy link
Contributor Author

I edited the inference to work without the fout call. Also users can pass an algorithm choice. It's an optional arg for type-stability, but it can be made a kwarg. That won't be stable of course until Julia v0.7, but if it's breaking it might as well just happen now?

@david-pl
Copy link
Member

I don't think the additional argument can break much. On the other hand, I don't know how much this instability will even affect things. So I think either way is fine.

@coveralls
Copy link

Coverage Status

Coverage decreased (-0.3%) to 99.234% when pulling 8c4c3e9 on ChrisRackauckas:diffeq into 74420ff on qojulia:master.

@ChrisRackauckas
Copy link
Contributor Author

MCWF is fixed to now use the DiffEq callbacks. alg is a new keyword argument for algorithms. Because of #189, the scope it's in is already type-unstable and probably won't be fixed before v0.7 so it might as well be in its "final form"

@ChrisRackauckas
Copy link
Contributor Author

I think that's everything. DDEs/SDEs/etc. can be a separate PR. But this should be everything switched over to DiffEq internals, with the breaking changes that fout now doesn't require user push!ing and instead saves the return, along with the SS saving at tstops which in the tests includes the first point.

@coveralls
Copy link

Coverage Status

Coverage decreased (-0.5%) to 99.007% when pulling b5c73c8 on ChrisRackauckas:diffeq into 74420ff on qojulia:master.

@coveralls
Copy link

coveralls commented Jan 11, 2018

Coverage Status

Coverage decreased (-0.3%) to 99.248% when pulling 8c4c3e9 on ChrisRackauckas:diffeq into 74420ff on qojulia:master.

@coveralls
Copy link

Coverage Status

Coverage decreased (-0.5%) to 99.007% when pulling b5c73c8 on ChrisRackauckas:diffeq into 74420ff on qojulia:master.

@david-pl
Copy link
Member

I ran the benchmarks on MCWF again. Looks like your final changes will give us quite a performance boost.

bench-mcwf

The green one here is your latest commit, the others are the same as before.

So let's merge this beauty!

@david-pl david-pl merged commit ff4af53 into qojulia:master Jan 12, 2018
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants