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

[WIP] Adding pooled callbacks to VectorContinousCallback #523

Draft
wants to merge 7 commits into
base: master
Choose a base branch
from

Conversation

freemin7
Copy link

I need help with this. I am currently running into:
ERROR (update): MethodError: Cannot convert an object of type Array{Int64,1} to an object of type Int64
at /.julia/packages/OrdinaryDiffEq/JrtsK/src/integrators/integrator_utils.jl:251
Looks like integrator.vector_event_last_time expects to be an Integer.

@freemin7
Copy link
Author

I am not convinced that it handles non root-finding callbacks correctly. There are currently no test on callbacks at all.
This WIP PI is related #519 .

@kanav99
Copy link
Contributor

kanav99 commented May 31, 2020

https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/src/integrators/integrator_utils.jl#L251 this is where you are failing. You have to make it's type Union{Int,Array{Int}}

@ChrisRackauckas
Copy link
Member

How come this isn't caught by the tests?

@kanav99
Copy link
Contributor

kanav99 commented May 31, 2020

Because they are behind the pool_events and we never reach them, default value of pool_events is false

@ChrisRackauckas
Copy link
Member

Ahh okay, so yeah it definitely needs a downstream test.

@freemin7
Copy link
Author

Well test/callbacks.jl is quiet sparse in general. I am not sure if normal callbacks get tested appropriately by it.

@devmotion
Copy link
Member

@@ -209,14 +213,13 @@ function VectorContinuousCallback(condition,affect!,len;
affect_neg! = affect!,
interp_points=10,
dtrelax=1,
abstol=10eps(),reltol=0)
abstol=10eps(),reltol=0, pooltol=missing, pool_events=false)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

IMO usually in DiffEq default values are set to nothing (as indicated by the docstring) or a sensible default, maybe one could use the same as abstol?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is very prototypie but thanks for your input, the doc strings aren't really updated ...

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The convention for the pooltol is also still in discussion, so i won't merge that proposal but keep in mind that doc strings would need to be added

Comment on lines +760 to +769
if callback.pool_events
tmp = get_condition(integrator, callback, integrator.dt + new_t)
if callback.pooltol isa Missing
# This is still dubious
pool_tol = eps(integrator.t + new_t) + eps(typeof(tmp[end]))
else
pool_tol = callback.pooltol
end
min_event_idx = findall(x-> abs(x) < pool_tol, tmp)
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm a bit surprised that one has to redo all calculations here. Is it not possible to change lines such as

if integrator.tdir * Θ < integrator.tdir * min_t
min_event_idx = idx
min_t = Θ
end
such that min_t and min_event_idx are not replaced but that one gets a list of min_t for every index, which then allows to select all indices within some tolerance later? And to change
min_event_idx = event_idx[1]
and
min_event_idx = event_idx[1]
to not just select the first index?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes we discussed that, it is the initial PR, so lets get that perfect

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The code is not optimized since it does not work yet.
I have decided to keep the pooled callback logic out of branches for different types of time resolving, so there is just a single place to change it.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this is necessarily the right direction. You're not looking for if the other values are near zero at the same time, you're looking for if there is another zero crossing nearby. One decent approximation would be to find the earliest event, and the move forward in time by some tolerance, and check the condition. You'd then take all of the events that were triggered in that time, and you'd know whether they are up or down crossings. Multiple crossings are extremely unlikely if this is on the size of a*eps(t). This would be insensitive to the relative size of the values.

This should show up in my repo only
Updated my branch to keep up
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