What to do with integer inputs?

[ivarne]

A common dilemma in Julia, is what to do when the user provides integers to a function where they probably want do do the calculation in floating point. Many users do not care about the difference between 13. and 13, and does not want the "wrong" result when they forget the ..

Suggestions

1. Try to solve the problem in integer space, and fail with InexactError 99.9 % of the time
2. Convert to float internally, but convert back to Int before returning
3. Error
4. Silently convert T <: Integer to float, and return floating point results.

It seems to be the general opinion that ODE solvers are performance sensitive, and we want them to be type stable, both internally and externally, to improve performance. We then have to decide what we want to do when the user supplies integers to tspan and y0.

[tomasaschan]

It's good that we make a consious decision on this, but I can't think of a single use case when the user wouldn't want and expect option 4.

[ivarne]

You are pretty leet hacker if you want anything else than 4. I just wanted to show the alternatives. 2 is just plain stupid, 3 is correct and 1 is the generic solution where we just do what we're are asked to do.

[stevengj]

4 — The common convention in Julia is, in a function (like ODE integration) where only floating-point numeric types make sense, other numeric arguments are promoted to floating-point as needed.

(See e.g. any of the transcendental math functions like sin, the quadgk integration function, etc.)

This isn't any less "correct" than 3. It's not a question of correctness, it's a question of defining the function in the most useful way.

[stevengj]

4 also applies to rational inputs, for the same reasons.

[stevengj]

See e.g. how quadgk does the promotion. You do the implementation assuming everything is the same floating-point type (or complex version thereof), and then write a separate generic method that promotes all arguments to the same type as needed.

[tomasaschan]

@stevengj I'm not sure quadgk's approach is feasible here; since we want to be able to duck-type y0, we can't use generics to restrict the actual implementation to calls with the same kinds of floats. Won't the built-in promotion rules be enough for us, especially if we call float on the input arguments the first thing that happens? (And I assume here that float(x) is a no-op if x is already a floating-point number, please correct me if I'm wrong...)

[ivarne]

I think we should explicitly check for nonsensical y0 types before doing anything but duck typing (that result in NoMethodErrors).

It looks like everybody agrees that if a user supplies machine integers Int as y0, it's a simple "mistake" that we can automatically fix by converting to Float64. Other Integer and Rational types are much more likely to be a deliberate choice the user made, and I think there is something to be said for not doing conversions we do not know the user wants.

if(y0 <: Int || eltype(y0) <: Int)
    y0 = float(y0)
end

[tomasaschan]

@ivarne Sure, but at the same time it's reasonable to assume that any user of an ODE library realizes that the solution space for any ODE problem is spanned by (a tensor power of) the real numbers. Trying to restrict the solution to an ODE to rational numbers just doesn't make sense - at least not with the algorithms we have implemented. Thus, if the input arguments do not support (and behave well) under float conversion, there's something fundamentally wrong with them.

[ivarne]

You guys probably have a better understanding about what kind of numeric types it makes sense for ODE solvers to run on. I just wanted to challenge your assumptions to see if there are fundamental reasons why we should always convert everything to T <: FloatingPoint or if that is just old habits. For Int there is convenience because the Julia parser assumes numbers without . to be integers. I am not sure the convenience argument is that strong for other types.

We could (in the future) create a keyword arg to disable the float() conversion if someone actually finds it problematic, but that would make our interface type unstable.

[stevengj]

@tlycken, a simple technique would be to first apply quadgk-like promotion rules to tspan to get a floating-point type T, and then compute y0 * one(T) to get the initial condition in the corresponding vector-space type (since we necessarily require the y0 to live in a normed vector space, multiplication by floating-point scalars must be a supported operation).

@ivarne, I don't understand why you think an integer input is a "mistake". Is it a mistake to type sin(3) or sin(4//3)? Integer or rational-valued initial conditions are perfectly sensible and produce a well-defined result. The fact that the results are in general transcendental, however, means that the output should be a floating-point type regardless of the input type. (Correspondingly, ODE quadrature algorithms are not sensible to apply to non-fp types.)

[tomasaschan]

@stevengj Yeah, that's a good approach to get a correctly typed y0 from the provided one. But what should the method signature of the "stricter" method (i.e. the one called after promotion floats) be? How do we define a type on y0 that restricts input of non-floats, but in no other way restricts the type of the argument?

[stevengj]

@tlycken, the whole point of duck typing is that you don't need to declare the type. When the odeXX routine is called, Julia will specialize it as needed for the type that was passed.

(Nor do I think that you need to call a separate method after promotion of y0. Just do y0T = one(T) * y0 and continue from there. The analogous thing in quadgk is the evalrule function: regardless of the type returned by the integration rule f, the first thing we do is to sum and multiply by the integration weights, which are in the desired precision, to get Ig and Ik. The resulting function is fast and type-stable.)

In other words, you define something like:

function odeXX{T<:FloatingPoint}(f::Function, y0, tspan::AbstractVector{T}; ....)
   y0T = one(T) * y0
   ....
end
function odeXX(f::Function, y0, tspan; kws...)
    s = start(tspan)
    done(tspan, s) && throw(ArgumentError("tspan must not be empty"))
    (v, s) = next(tspan, s)
    T = typeof(float(v))
    while !done(tspan, s)
         (v, s) = next(tspan, s)
         T = promote_type(T, typeof(v))
    end
    odeXX(f, y0, T[t for t in tspan]; kws...)
end
odeXX(f::Function, y0, t1, t2, ts...; kws) = odeXX(f, y0, [t1, t2, ts...])

where I have used start/done/next to support arbitrary iterable types for tspan.

[ivarne]

@stevengj I do think sin(3) is kind of a mistake on the users part, but just like with 1/2 it is obvious that we must convert to Float64 before calling the Float64 version of the function because the user probably intended to do sin(3.). For rationals the picture is more blurry, because nobody will use a rational by accident.

When it comes to ODE solvers are generic algorithms that does (to my knowledge) generic operations where it might make sense to allow exotic types to propagate. When the algorithms makes assumptions about the numeric properties of the underlying type, they can raise errors if the given type does not make sense.

Also I have not checked with SIUnits or other wrapping types, but as we do not have Multiple inheritance, they will often be unable to inherit from everything you want them to. I just feel that it is useless to restrict values to T <: FloatingPoint in the API.

[stevengj]

You need to be careful to avoid dispatch loops. Actually, even the snippet I gave above has a problem for Complex types. Perhaps a better solution would be to do, instead of odeXX{T <: FloatingPoint}:

function odeXX{Tn <: Number}(f::Function, y0, tspan::AbstractVector{Tn}; kws...)
    T = typeof(float(one(Tn)))
    tspanT = T[t for t in tspan]
    !isreal(tspanT) && throw(ArgumentError("only real tspan are supported"))
    y0T = one(T) * y0
    .... continue ODE integration using tspanT and y0T
end

which will work with any Number type as long as float and isreal is defined for it, and otherwise throws an error. I'm not completely happy with the prospect of making an unnecessary copy of tspan when the input is already a floating-point type, but this could be avoided with a bit more code.

ODE integration over (e.g.) complex tspan could potentially be supported as integration over straight-line segments in the complex plane. However, this requires some care to make it actually work, and is likely to require a separate implementation for performance reasons.

But I don't see any sensible use case for not doing some kind of floating-point promotion in an ODE solver, whether via float or some similar function. e.g. a slightly fancier thing would be to do T = typeof(float(one(Tn)) * tspan[1]), which will preserve the units for the SIQuantity type from SIUnits.

There is a perfectly sensible reason to use Rational values as initial conditions, and as constants in general: even if the computations are done in floating point, this allows you to represent a rational constant so that it gets promoted to the closest floating-point value in whatever precision is used for tspan. Basically, it can be used analogous to the MathConst type.

[stevengj]

(PS. I disagree that sin(3) is a mistake. Automatic promotion is a hugely useful feature, and the extensibility of promotion rules is one of the awesome things about Julia. Whether you agree with it or not, it is pervasive in Base, and it would be very un-Julian not to do automatic promotion here. The only possible alternative would be to throw an error.)

[ivarne]

(@stevengj Good to see that we actually agree on the important stuff. Promotion rules that extends to user defined types is definitely awesome, but the discussion is not about promotion in general. Julia has not defined sin(3) nor 1/2 == 0.5 to be a mistake and I do not think it should be. It is likely that unsuspecting users want the Float64 answer, and knowledgeable users understand that they have to use div to avoid conversion. The question is if Rational and all other non Floatingpoint numbers should behave like a MathConst.)

[stevengj]

Using them as (effectively) a MathConst (to use your phrase) is the only sensible meaning in this context. There is no reasonable way to use integer or rational values in a numerical ODE solver except by promoting them to floating point.

Given that fact, it is unhelpful to throw an error.

[stevengj]

On a separate note, an easier way to prevent dispatch loops here is to just call an odeXX_ function with a different name (e.g. with an underscore) at the bottom level.

[ChrisRackauckas]

There is no reasonable way to use integer or rational values in a numerical ODE solver except by promoting them to floating point.

Yes there is. With fixed timestep methods this is clearly possible. With adaptive methods, you just have to check your timestepping to change by rational values. One scheme for this is do reject/accept by doubling/halving the timestep.

Using integer and rational values for time can be good for ODE solvers which form the backend to DDE solvers since those are not very robust to small intervals that are due to propagation of discontinuities (this arises since 1/3 + 1/3 + 1/3 != 1 in floating point).