Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Repo info
Activity
  • Dec 28 2016 09:12
    coveralls commented #28
  • Dec 28 2016 09:12
    coveralls commented #28
  • Dec 28 2016 09:12
    coveralls commented #28
  • Dec 28 2016 09:12
    coveralls commented #28
  • Dec 28 2016 09:12
    codecov-io commented #28
  • Dec 28 2016 09:07
    codecov-io commented #28
  • Dec 28 2016 09:05
    coveralls commented #28
  • Dec 28 2016 08:59
    ChrisRackauckas commented #28
  • Dec 28 2016 08:58
    ChrisRackauckas synchronize #28
  • Dec 28 2016 08:55
    ChrisRackauckas synchronize #28
  • Dec 28 2016 08:54
    coveralls commented #28
  • Dec 28 2016 08:54
    coveralls commented #28
  • Dec 28 2016 08:54
    coveralls commented #28
  • Dec 28 2016 08:54
    coveralls commented #28
  • Dec 28 2016 08:54
    codecov-io commented #28
  • Dec 28 2016 08:43
    ChrisRackauckas synchronize #28
  • Dec 22 2016 05:11
    ChrisRackauckas commented #28
  • Dec 12 2016 06:35
    coveralls commented #28
  • Dec 12 2016 06:35
    coveralls commented #28
  • Dec 12 2016 06:35
    codecov-io commented #28
Paweł Biernat
@pwl
just a thought!
Mauro
@mauro3

Yes, you're right, changing the return type depending on a keyword argument is bad.

Concerning your example, yes, there can be corner cases. The beauty of building on iterators, is that one can always drop down to them and do what one wants. Alternatively, we could only return a matrix when isa(y0, Vector), that should cover 95% of the use-cases.

When the number of output points is not known, then make a vector, push to that and reshape in the end. (Note that the reshape does not copy but just re-interprets the memory as matrix).

I still don't want to call something a solution which is not a solution at all but can only be used to, maybe, generate a solution.

Eventually having a solution type would be good (but in another PR). And having an option to store all the info to do dense output on the solution ex-post would be cool. However, it should be optional, as for big problems this may well be more info than one wants to store.
Paweł Biernat
@pwl
so, in preparation to implement having a custom solution type, we should have functions extracting the raw data out of the return type of solve, which, for now, could be whatever we feel like.
Mauro
@mauro3
gett, gety, getdy?
Paweł Biernat
@pwl
makes sense
Paweł Biernat
@pwl
getty to return (t,y)?
Paweł Biernat
@pwl
I would argue that if we are going to support custom user defined types (see my Point example) we can't treat them as a corner case. I would rather have the usual output from solve and a special solve_matrix for anyone who want's the matrix output.
Paweł Biernat
@pwl
I'm playing around with minimal type requirements and I was considering a type without getindex. In principle, we don't need the index for some solvers to work, we only need a type resembling an element of a (mathematical) vector space
so no need for length and getindex
of course it will work only for explicit solvers
any thoughts on that?
Paweł Biernat
@pwl
looks like the error estimate from the RK method is the bottleneck here, the formulas from Hairer&Wanner require component-wise estimates.
they cannot be replicated in a purely vector space manner (i.e. without referring to the components)
Christopher Rackauckas
@ChrisRackauckas
FYI many algorithms run more efficiently when not doing the component-wise estimate in L2. L infinity is more sensitive to it.
At least that's what I've seen from ~10 ODEs and ~5 SDEs
Paweł Biernat
@pwl
can you rephrase your first comment?
Christopher Rackauckas
@ChrisRackauckas
So I purposely leave the length scaling off by default.
If your error estimate is in L2, for many problems using the standard L2 norm instead of the Hairer scaled-L2 norm ( sqrt(sum(squares)/length)) can be more efficient.
You can run a bunch of benchmarks yourself and see that it's not obviously better to scale it. Your tolerance will be slightly off for large (100x100) problems.
Paweł Biernat
@pwl
yeah, I don't see the obvious reason for the length scaling in Hairer
Christopher Rackauckas
@ChrisRackauckas
It's so that the tolerance lines up.
But not scaling works pretty well too since it's all local...
Paweł Biernat
@pwl
but this is not the issue I have. The problem is the L2 norm is not actually a norm because the weights depend on the components.
Christopher Rackauckas
@ChrisRackauckas
The L2 norm is a norm...?
ODE.jl's main thing is that your standard tolerances are really weird.
Paweł Biernat
@pwl
sorry, I mean that in Hairer they have a formula for the error that looks like an L2 norm but it is actually not a norm, because the weights in the sqrt(sum(squares/weights)) depend on the state of the solution.
Christopher Rackauckas
@ChrisRackauckas
Oh yeah
At each timepoint it's a norm though
Paweł Biernat
@pwl
no, because it is never linear iny
ok, in a sense you are correct
Christopher Rackauckas
@ChrisRackauckas
I don't see why it would matter though.
Paweł Biernat
@pwl
anyway, what did you mean by the default tolerances?
Christopher Rackauckas
@ChrisRackauckas
your reltol and abstol
Paweł Biernat
@pwl
what would you like them to be, and why?
Christopher Rackauckas
@ChrisRackauckas
usually they are like 10^-6 and 10^-3 or something like that.
I thought that's what MATLAB does, and I know DOPRI (/all of ODEInterface) does that.
Paweł Biernat
@pwl
well, I won't argue for our current choice, I always specify the tolerances by hand anyway
but what's the rationale behind a particular choice, apart from the fact that other solvers use them?
Christopher Rackauckas
@ChrisRackauckas
10^-3 means you're right in the first 3 digits.
Anything beyond that is usually hard to see, at least when plotted.
It's a nice default, and if people want higher percision, they probably should be using the default alogrithm (DP5) anyways.
Mauro
@mauro3
@pwl: sorry I have not had time to look at our work recently. And will not have much time at all until early September.
Paweł Biernat
@pwl
sure, no worries, I'm pretty pressed with other stuff as well and work on ODE.jl only occasionally.
Mauro
@mauro3
tnx
Christopher Rackauckas
@ChrisRackauckas
ODE.jl fails to load on v0.4.6 for my tests (on your dev branch). Is this expected?
LoadError: LoadError: LoadError: type Array has no field x
 in anonymous at /home/travis/.julia/v0.4/ODE/src/tests/test_cases.jl:68
 in call at /home/travis/.julia/v0.4/ODE/src/base.jl:84
 in include at ./boot.jl:261
 in include_from_node1 at ./loading.jl:320
 in include at ./boot.jl:261
 in include_from_node1 at ./loading.jl:320
 in require at ./loading.jl:259
 in eval at /home/travis/.julia/v0.4/DifferentialEquations/src/DifferentialEquations.jl:3
 in init_package at /home/travis/.julia/v0.4/DifferentialEquations/src/general/backends.jl:23
 in init_package at /home/travis/.julia/v0.4/DifferentialEquations/src/general/backends.jl:16
 in initialize_backend at /home/travis/.julia/v0.4/DifferentialEquations/src/general/backends.jl:5
 in solve at /home/travis/.julia/v0.4/DifferentialEquations/src/ode/ode_solve.jl:230
This started about 20 hours ago.
Paweł Biernat
@pwl
It's been a while since I run tests on 0.4, I can reproduce this.
The issue is the new test case involving a custom type
On the other hand we were thinking of dropping the support for 0.4, now that 0.5 is so close.