API for Parameters #2
Comments
|
This scares me, as it has the potential for adding an additional api idiom to Julia that is likely to be different from other similar api's (such as root solvers, optimization solvers, anything that takes a callback). I don't really see why not arguing for a closure isn't the best solution (outside of the mentioned speed regression, which seems like premature optimization, with a strong possibility of making requiring the parameter call to be the default, even if the speed regression is fixed. I think this since so many people will be used to the "old way" of doing this from C/Fortran etc). My feeling is that API choice is rarely a good solution, as it trades convenience for consistency. I do highly support pushing Parameters.jl though :) it is the best (I am dying for your Dict unpack to land). |
|
I think it's moreso because we have to have one, so it's about finding the right choice. It's not essential for projects like ODE.jl which are just solving ODEs, but the moment you start to go past that, you need the parameters. You need the parameters to be able to take parameter derivatives for sensitivity analysis (which is why Sundials does So it's more of, do we just put parameters as the fourth argument (like Sundials)? Is that our thing? Make it a type so that the fields can be anything? That seems like a reasonable choice, I just think we should make a good choice, know the downsides , and stick with it. Sticking with it leads to consistency. If Sundials and DifferentialEquations.jl does it, then ODE.jl should be made compatible, which since it's domain is simpler it can just use a closure to make |
|
It wouldn't be pushing Parameters.jl, but if the choice is to go by defining parameters as types (somehow, but how do you deal with namespace issues? Like if I want to define parameters in a macro and give that parameter |
|
I guess I am missing the issue. If I pass a callback to a bifurcation routine, then I definitely need to pass the parameters, but I don't if I am solving an ODE. I don't see why the later should be made consistent with the former, they have different API's. Similar for sensitivity analysis, that routine would expect a callback with a different signature. Having a universal API for all of these domains doesn't feel right to me. If I am doing a maximum likelihood problem I need to pass the data and the parameters, but that doesn't mean an optimization callback should always take a data argument. This feels like it would require ever increasing argument lists if every possible application area needs to be modelled. |
|
But in modern applications these things are intertwined. That's why Sundials gives you the sensitivity analysis results at the same time as solving the ODE when calling The functions for differential equations only go that far. Is there any application you can think of where more than the parameters is needed as arguments to the function? I can't think of any. There is more that can be associated with them, which is why I went the way of the If there was a PR to ODE.jl so that it first looks at the function for the number of arguments, and if it's 4 it closes the optional userdata into it (like Sundials), would that be accepted? Because I'm going to do something because I'm going to keep pushing DifferentialEquations.jl forward, and I'm going to start making bifurcation analysis tools and have them be compatible: I just want to make sure I head in the right direction with it and, while keeping compatibility with Sundials, see if ODE.jl will follow. |
|
It is hard for me to fully grasp the Sundials approach since it is written in C and therefore closures are not really an option (though like you are suggesting it could be that even if closures are possible then the parameters argument is still required for the general solution). Maybe this is the only way forward. I just am sad that soon the examples people share will be mixed up with parameter types being passed, or closures, or some kind of macro dsl ... I feel Julia does things a variety of ways far to often, since it can, but this makes the burden of learning the language and the libraries brutal, and makes it hard for me to recommend it to colleagues over more restrained languages like matlab or python. I also hate that this will be so different from using integration routines and optimization, but I admit you have greater knowledge here, so it might be required. |
|
The reason for this issue is that in some situations one needs additional information about the ODE (parameters, Jacobian, etc). This is similar to optimization problems. IMO the Julian way to solve this would be via "problem types", where by default additional information is discarded or obtained in some automatic way. The user can provide additional information by overwriting functions and we can make use of multiple dispatch. I think this is the approach taken by JuliaOpt. However, to make the libraries more accessible and the learning curve not so steep, we would still provide a "simplified interface" like it is now. In many cases this will also be sufficient. IMHO the possibility to do things different ways in Julia is actually its strength. In particular, you can adjust your code with respect to simplicity vs. performance, but you don't have to leave the language (like in Matlab or python). Of course macros et al make the code more complicated, but if you really need the performance you will embrace the possibilities. In this case it's not so much about performance, but rather about functionality, which we cannot offer with the current API. |
|
I think parameters are slightly different though. You can specify a Jacobian and stuff in a problem type. But to take parameter derivatives for sensitivity and things like that, you need them to be part of the function. I think parameters are the only thing that need to specifically be part of the function for their full use. So most parts of the ODE definition will be handled as fields in problem types (as I am doing already), but I think parameters need to be slightly different. What makes me uncomfortable though is that if it's passing a type as the parameter function testdefine()
type TempType end
tmp = TempType()
end
testdefine()That will error since you cannot define a type inside a local scope. Would that be necessary? |
|
Speaking of JuliaOpt / JuMP, one of the things they did was turn to macros for user definitions. This is to hide some of the complexity that's actually going on. I have some very new macros for doing that here: http://juliadiffeq.github.io/DifferentialEquations.jl/latest/man/function_definition_macros.html (they should generate functions compatible with all four major ODE packages). I plan to expand it a little bit like; f,p = @ode_define begin
dy₁ = -k₁*y₁+k₃*y₂*y₃
dy₂ = k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
dy₃ = k₂*y₂^2
end k₁=>0.04 k₂=>3e7 k₃=1e4where type TmpParameters
k₁
k₂
endwith |
|
I totally agree with @gabrielgellner, there is not enough arguments for not using closures. They may not be very intuitive at first for some new users (mostly coming from different programming languages) but these users don't come back to complain after we explain them how the closure works. There might be some performance drawbacks but from my understanding these drawbacks only show up in the global scope and are bugs in Julia, so in a sense, they are only temporary. Designing an API basing on comments by newcomers and to circumvent a bug in Julia would be a huge mistake. When I was implementing an ODE suite in C or Fortran I was forced to add parameters (no closures in these languages) and I ended up saving the parameters in types and/or passing the parameters to every subroutine that had to evaluate the equation; it was painful to look at and needlessly complicated. All in all the approach proposed by @ChrisRackauckas looks to me like something coming straight from C/Fortran. @acroy we have problem types in ODE.jl (the forked version) but as @ChrisRackauckas commented, it is a slightly different issue. Also we try to keep these types hidden from the user by defining simpler layer of "pedestrian" API.
The API we are discussing is a very simple one and I don't think we have to resort to macros at this point. On top of that, if new users find closures counterintuitive they will have trouble with macros as well. |
|
What about just having the front ends accept both |
|
How would you distinguish between EDIT: we could have |
|
Is it not possible to have a modular approach: the solver use closures and the tools which need direct access to the parameters, such as bifurcation and sensitivity analysis tools, can carry the parameters, and create a closure on the fly when passing computations off to solvers. Ideally a lower level part of a stack should not need to care about a higher level part. |
|
@pwl You can just read the number of parameters and swap out the function that's actually used using a closure (for safety, you can isolate this and then have it call your inner real function). I already do that, and I thought you were doing that on your branch with changing functions from not in-place to in-place. |
|
@pwl: The point is that you can't use a closure since you need the explicit value of the parameter(s) in the solver, not only in the rhs function. And exactly since you can't dispatch on function signatures, one would have to introduce a type. All of this doesn't mean that the current interface has to go away. |
"""
`numparameters(f)`
Returns the number of parameters of `f` for the method which has the most parameters.
"""
function numparameters(f)
if length(methods(f))>1
warn("Number of methods for f is greater than 1. Choosing based off of method with most parameters")
end
numparm = maximum([length(m.sig.parameters) for m in methods(f)]) #in v0.5, all are generic
return (numparm-1) #-1 in v0.5 since it adds f as the first parameter
endCurrently, when defining the problem type, I have this read the parameters and change a function which is not in place |
|
@acroy could you provide an example of how a solver you describe would work? It is hard for me to argue if I don't see the specific use case. |
|
The keyword is "sensitivity analysis" (there are plenty of examples out there), where you basically need the derivative of |
|
Well, you guys know much better the internals of ODE 2.0 etc, but I don't see why an ODE with parameters is not a special subtype of an ODE problem? |
|
For the ignorant ones, like myself, there is quite a nice, short description of sensitivity analysis in CVODES manual: http://computation.llnl.gov/sites/default/files/public/cvs_guide.pdf. Is my understanding correct that sensitivity analysis needs to be built into the solver and cannot be bolted-on on top, right? If so, and if we want to support it eventually in ODE.jl, then a way to pass parameters is indeed needed. One thing to note is that the ICs are also parameters, so it would probably be best to pass IC and other parameters in the same way. |
|
@mauro3 Why not have a separate type EDIT: we could then have a top level interface defined as |
|
While ICs mathematically are parameters (or functions), they are different since they don't need to be passed to the function each timestep. And sensitivity analysis is just one place. One application that is done a lot in uncertainty quantification is to put noise on parameters and solve the ODE a bunch of times (note that this is different than solving the SDE since it's measuring the numerical uncertainty) (this is related to sensitivity analysis). Another application is types of control problems where it's part of the solver to feedback onto a parameter depending on error estimates. |
|
It looks like we are mixing two things here: parameters for the sake of sensitivity analysis and parameters in the sense of constants coming from the user; let's call them sensitivity parameters and user parameters respectively. These two are completely independent use cases of parameters. The user parameters can be (and are) implemented with closures, there is no need for any special API in this case. Modulo some bugs and minor user complains this is the most elegant solution in sight. On the other hand, sensitivity parameters have to be part of the definition of an ODE problem as in One reason for this is because most users (myself included) don't deal with sensitivity analysis so there is no need to pass around a void argument. The other reason is that sensitivity parameters seem to be more strictly defined (I think they can only be scalars, so To summarize I vote for defaulting to Am I missing some use cases? @ChrisRackauckas You mentioned that one could add noise through |
|
Sundials already supports sensitivity analysis in The uncertainty quantification usage is adding noise to The parameters are mathematically the same thing in both of the cases you mention. The only difference is whether you use them explicitly (to be able to check/affect them), or use them as constants. And I think you could do things like pass an array as a parameter as long as it's not a parameter that's used for the sensitivity analysis. Using a type for parameters would let you choose the parameters for this analysis by field name, defaulting to I don't think there's a need to default to |
|
I am still not convinced that Sundials has To reiterate my fear, in hopefully a more clear manner, is that if all routines (ones that need them or not) allow a "parameter" argument, then that will be used by end users whether it is needed or not (giving a closure friendly version as well will not matter, both will be used, likely passing the parameter will be used more often as it is more like "older" languages). And then when teaching/using these packages you will need to convince people to use closures versus passing arguments, and all the subtleties between using one api versus the other. My experience is that this will just scare people way. I have found that having, what seems, like a more complex api is not always a problem for end users. But having many api's is. That is why I don't like having many ways to define the same RHS. If you need to learn to do the u[1] = blah style definition, the complexity needs to be learned once. If sometimes you use a macro to define the types and you can so foo = blah, but it gets translated to u[1] = blah behind the scenes, you are adding a tremendous level of complexity for teaching and reasoning about the package, in a sense you now require everyone to learn two api's versus one (even if it seems like everything gets reduced to the same base form) if you have any hope of reading other peoples code. My point about DSL with macros etc., is not that I don't like macros and Julia's flexibility, but rather fear of a Lisp like world were every library has its own semantics, or Perl where there is a million and one ways to do anything, so it is really hard to learn/read other peoples code. My feeling is that sometimes a less optimal solution gives a far better user experience because when looked at from a more high level it gives consistency across a language and makes everything feel more discoverable. Mathematica is an amazing example of this. Anyway I will just lurk now. I love all the packages in JuliaDiffEQ, so whatever is decided is fine. But I just want to say, as an end user, that API consistency (both in JuliaDiffEQ, and across Julia itself, hence the importance of closures in callbacks) is far more important to me than having an optimal API. |
|
@ihnorton had a really great idea of implementing the parameterized version via call overloading. For example, it would look like: type LotkaVoltera <: ParameterizedFunction
a::Float64
b::Float64
end
f = LotkaVoltera(0.0,0.0)
(p::LotkaVoltera)(t,u,du) = begin
du[1] = p.a * u[1] - p.b * u[1]*u[2]
du[2] = -3 * u[2] + u[1]*u[2]
endThen This then ties in really nicely with re-defining functions and the function definition macros. For example, a definition macro could look like f = @ode_def LotkaVoltera begin
dx = a*x - b*x*y
dy = -c*y + d*x*y
end a=>1.5 b=>1 c=3 d=1and then The end result is that all of this is a superset of the existing functionality which is compatible with the existing functionality. You can teach people to define equations for differential equations using the Does that hit all of your points @gabrielgellner? |
|
This is a very neat idea! I still don't think that macros are necessary here but the type structure is simple and elegant. Is there anything against having EDIT: we should still use closures by default, leaving |
|
That's a good idea. I disagree that the macros aren't necessary. It depends on how intense your models are. For simple models, sure you can define the ParameterizedFunction without the macro and you'll be fine. But the reason that I came up with the macros is because I actually had to implement the following for a recent publication. This isn't an exception: systems bio is filled with these kinds of differential equations. I'll be keeping them around because I am definitely going to use them to enhance readability in these situations. (A model I am working on now has about 50 reactants. If I use reactant names, I can easily double check the equations since I know the chemical reactions. If I use function f(t,y,dy)
# y(1) = snailt
# y(2) = SNAIL
# y(3) = miR34t
# y(4) = SR1 # abundance of SNAIL/miR34 complex
# y(5) = zebt
# y(6) = ZEB
# y(7) = miR200t
# y(8) = ZR1 # abundance of ZEB/miR200 complex with i copies of miR200 bound on the sequence of ZEB1
# y(9) = ZR2
# y(10) = ZR3
# y(11) = ZR4
# y(12) = ZR5
# y(13) = tgft
# y(14) = TGF
# y(15) = tgfR # abundance of TGF/miR200 complex
# y(16) = Ecad
# y(17) = Ncad
# y(18) = Ovol2
TGF0=.5(t>100)
#ODEs
dy[1]=k0_snail+k_snail*(((y[14]+TGF0)/J_snail0))^2/(1+(((y[14]+TGF0)/J_snail0))^2+(y[19]/J_SO)^nSO)/(1+y[2]/J_snail1)-kd_snail*(y[1]-y[4])-kd_SR1*y[4]
dy[2]=k_SNAIL*(y[1]-y[4])-kd_SNAIL*y[2]
dy[3]=k0_34+k_34/(1+((y[2]/J1_34))^2+((y[6]/J2_34))^2)-kd_34*(y[3]-y[4])-kd_SR1*y[4]+lamdas*kd_SR1*y[4]
dy[4]=Timescale*(Ks*(y[1]-y[4])*(y[3]-y[4])-y[4])
dy[5]=k0_zeb+k_zeb*((y[2]/J_zeb))^2/(1+((y[2]/J_zeb))^2+((y[19]/J_2z))^nO)-kd_zeb*(y[5]-(5*y[8]+10*y[9]+10*y[10]+5*y[11]+y[12]))-dk_ZR1*5*y[8]-dk_ZR2*10*y[9]-dk_ZR3*10*y[10]-dk_ZR4*5*y[11]-dk_ZR5*y[12]
dy[6]=k_ZEB*(y[5]-(5*y[8]+10*y[9]+10*y[10]+5*y[11]+y[12]))-kd_ZEB*y[6]
dy[7]=k0_200+k_200/(1+((y[2]/J1_200))^3+((y[6]/J2_200))^2)-kd_200*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])-dk_ZR1*5*y[8]-dk_ZR2*2*10*y[9]-dk_ZR3*3*10*y[10]-dk_ZR4*4*5*y[11]-dk_ZR5*5*y[12]+lamda1*dk_ZR1*5*y[8]+lamda2*dk_ZR2*2*10*y[9]+lamda3*dk_ZR3*3*10*y[10]+lamda4*dk_ZR4*4*5*y[11]+lamda5*dk_ZR5*5*y[12]-kd_tgfR*y[15]+lamdatgfR*kd_tgfR*y[15]
dy[8]=Timescale*(K1*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])*(y[5]-(5*y[8]+10*y[9]+10*y[10]+5*y[11]+y[12]))-y[8])
dy[9]=Timescale*(K2*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])*y[8]-y[9])
dy[10]=Timescale*(K3*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])*y[9]-y[10])
dy[11]=Timescale*(K4*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])*y[10]-y[11])
dy[12]=Timescale*(K5*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])*y[11]-y[12])
dy[13]=k_tgf-kd_tgf*(y[13]-y[15])-kd_tgfR*y[15]
dy[14]=k_OT+k_TGF*(y[13]-y[15])-kd_TGF*y[14]
dy[15]=Timescale*(TGF_flg+KTGF*(y[7]-(5*y[8]+2*10*y[9]+3*10*y[10]+4*5*y[11]+5*y[12])-y[15])*(y[13]-y[15])-y[15])
dy[16]=GE*(k_ecad0+k_ecad1/(((y[2]/J_ecad1))^2+1)+k_ecad2/(((y[6]/J_ecad2))^2+1)-kd_ecad*y[16])
dy[17]=k_ncad0+k_ncad1*(((y[2]/J_ncad1))^2)/(((y[2]/J_ncad1))^2+1)+k_ncad2*(((y[6]/J_ncad2))^2)/(((y[6]/J_ncad2)^2+1)*(1+y[19]/J_ncad3))-kd_ncad*y[17]
dy[18]=k0O+kO/(1+((y[6]/J_O))^nzo)-kdO*y[18]
dy[19]=kOp*y[18]-kd_Op*y[19]
end |
|
I really like the parameterized solution. I could see myself using that even for regular ODE problems with the nice macro syntax. What I really like about this is how regular it is. It doesn't feel like a special case, but rather a really cool way of using Julia's call overloading to get at a kind of functor like object. Again my issue with the macro solutions was how it felt like a new dialect, using the call overloading approach and the macro together feels like it could easily become the default way of defining such problems. Really really cool solution. I will need to play with it a bit. |
|
One interesting thing that comes out of this: an easy way of specifying the Jacobian. It's a little hacky: function Base.ctranspose(p::LotkaVoltera)(t,u,J)
J[1,1] = a-b
J[1,2] = -b
J[2,1] = 1
J[2,2] = -2
endand then f = @ode_def LotkaVoltera begin
dx = a*x - b*x*y
dy = -c*y + d*x*y
end a=>1.5 b=>1 c=3 d=1to return an overloaded type where This could be a way for the macro to automatically supply the Jacobian to a Rosenbrock solver. |
|
Oh man that looks so sweet. Though I worry that LotkaVoltera example is just to make me love it too much :) |
|
I made a prototype. Check out ParameterizedFunctions.jl. The macro does calculate the Jacobian automatically as well (using SymPy.jl). This needs to be tested with the ODE solvers, but the tests in the package shows that it makes the right functions and Jacobians. It also has a macro for the Finite Element Methods. |
|
I like the first example, but the Jacobian-example gets in the cute but too obscure territory, IMO. |
|
So the bottom line of this discussion is that we can completely avoid Because there are many ways to include parameters we shouldn't be too specific with the API until we have a pure Julia solver with actual sensitivity analysis or similar features. Until then, let's keep the API simply as Can we close this issue now and move the discussion on |
|
Yeah sounds good. |
|
@mauro3 I came up with that because I was wondering if I could silently calculate and return the Jacobian as part of the construction of |
I want to consolidate loosely the topics:
The idea is for a general API change from
(t,y,dy)to(t,y,dy,p)where p is some type that holds the parameters. It seems to be a common feature request among users. Indeed, it's not truly necessary because on can use a closure here, but it can be tricky to do correctly due to JuliaLang/julia#15276. Also, to be compatible with tools which require "jiggling the parameters", like dynamical systems tools (bifurcation diagrams), sensitivity analysis, etc. the functions have to explicitly have the parameterThus I think it's best for across the board compatibility for the ODE packages to use
(t,y,dy,p). This would make the same functions usable in not only ODE solvers but for other related purposes. The ODE solvers would could have a simple change to write a closure and continue as usual. However, it's not clear what that type should be for the parameters. I would suggest using types for holding the parameters and embracing https://github.com/mauro3/Parameters.jl as a standard way to define such types (use it in examples).@pwl @mauro3 @acroy
The text was updated successfully, but these errors were encountered: