Re: [julia-users] { } vector syntax is discontinued?
Thanks for the info; it still isn't clear to me how the (generic) call to
this function should be changed:
function cvode{f,r,T}(::Type{f},::Type{r},d::Array{Int64},p::Vector{T}, y0
::Vector{Float64}, t::Vector{Float64}; reltol::Float64=1e-4, abstol::Float64
=1e-6)
(from https://github.com/sdwfrost/PDMP.jl/blob/master/src/cvode.jl)
On Monday, August 22, 2016 at 10:19:14 AM UTC+1, Mauro wrote:
>
> It used to mean Any[], but not anymore (so `Any[]` is the new syntax).
> It was still usable in 0.4 with a deprecation warning, now in 0.5 an
> error is thrown. It will get a new meaning in the 0.6 release, but it
> is not decided for what yet. See e.g.:
> https://github.com/JuliaLang/julia/issues/8470
>
> On Mon, 2016-08-22 at 11:13, Simon Frost >
> wrote:
> > Dear All,
> >
> > Apologies if this is mentioned somewhere, but I couldn't find the answer
> in
> > various searches. What does '{ } vector syntax is discontinued' mean,
> and
> > what is the new syntax?
> >
> > Best
> > Simon
>
[julia-users] { } vector syntax is discontinued?
Dear All,
Apologies if this is mentioned somewhere, but I couldn't find the answer in
various searches. What does '{ } vector syntax is discontinued' mean, and
what is the new syntax?
Best
Simon
[julia-users] Re: Performance issues with stochastic simulation
Dear Chris, I made a mistake in the old benchmarks - updated ones are now better. Knowing how to make function in-place would be great though! Best Simon On Thursday, July 28, 2016 at 11:20:58 AM UTC+1, Simon Frost wrote: > > Dear Chris, > > I've managed to get a little better performance by adding a loop and > avoiding some type conversions. > > https://gist.github.com/sdwfrost/8a0e926a5e16d7d104bd2bc1a5f9ed0b > > I'm still getting about a 7-fold slowdown from using the function call (as > an immutable object). How would one make my function F in-place, still > using a similar generic API? > > Best > Simon > > On Thursday, July 21, 2016 at 4:35:56 PM UTC+1, Chris Rackauckas wrote: >> >> You can change line 70 to be in place with a loop: >> >> for i in 1:length(x) >> x[i] = x[i] + deltax[i] >> end >> >> I don't think you can do >> >> x[:] =x .+deltax >> >> as fancy syntax here since the x is part of the statement though (you can >> check). This should cut out an allocation here and bring down the time. >> >> Do you need to use a WeightVec? If you do (for future things), keep the >> WeightVec separate from the Vector so that the types aren't changing. let >> wpf always be the WeightVec you make from pf. Otherwise pf isn't type >> stable. It would be best if you could make F in-place as well since this is >> where your bottleneck is. >> >> On Thursday, July 21, 2016 at 7:56:51 AM UTC-7, Simon Frost wrote: >>> >>> Dear All, >>> >>> I'm having some issues with code speed for some Gillespie type >>> simulations. The toy model is described here: >>> >>> >>> http://phylodynamics.blogspot.co.uk/2013/06/comparing-performance-of-r-and-rcpp-for.html >>> http://phylodynamics.blogspot.co.uk/2013/06/an-sir-model-in-julia.html >>> >>> I get good performance with my vanilla Julia code, but a more generic >>> implementation is slower: >>> >>> http://github.com/sdwfrost/Gillespie.jl >>> >>> The gist is here: >>> >>> https://gist.github.com/sdwfrost/1b4bce19faf2d7b8624cac048a36f32d >>> >>> Lines 57 and 70 appear to be the culprit: >>> >>> https://github.com/sdwfrost/Gillespie.jl/blob/master/src/SSA.jl >>> >>> I've tried some devectorisation, but in my hackery, I appear to get side >>> effects, where the argument x0 passed to the ssa function is modified. Any >>> tips? >>> >>> Best >>> Simon >>> >>
[julia-users] Re: Performance issues with stochastic simulation
Dear Chris, I've managed to get a little better performance by adding a loop and avoiding some type conversions. https://gist.github.com/sdwfrost/8a0e926a5e16d7d104bd2bc1a5f9ed0b I'm still getting about a 7-fold slowdown from using the function call (as an immutable object). How would one make my function F in-place, still using a similar generic API? Best Simon On Thursday, July 21, 2016 at 4:35:56 PM UTC+1, Chris Rackauckas wrote: > > You can change line 70 to be in place with a loop: > > for i in 1:length(x) > x[i] = x[i] + deltax[i] > end > > I don't think you can do > > x[:] =x .+deltax > > as fancy syntax here since the x is part of the statement though (you can > check). This should cut out an allocation here and bring down the time. > > Do you need to use a WeightVec? If you do (for future things), keep the > WeightVec separate from the Vector so that the types aren't changing. let > wpf always be the WeightVec you make from pf. Otherwise pf isn't type > stable. It would be best if you could make F in-place as well since this is > where your bottleneck is. > > On Thursday, July 21, 2016 at 7:56:51 AM UTC-7, Simon Frost wrote: >> >> Dear All, >> >> I'm having some issues with code speed for some Gillespie type >> simulations. The toy model is described here: >> >> >> http://phylodynamics.blogspot.co.uk/2013/06/comparing-performance-of-r-and-rcpp-for.html >> http://phylodynamics.blogspot.co.uk/2013/06/an-sir-model-in-julia.html >> >> I get good performance with my vanilla Julia code, but a more generic >> implementation is slower: >> >> http://github.com/sdwfrost/Gillespie.jl >> >> The gist is here: >> >> https://gist.github.com/sdwfrost/1b4bce19faf2d7b8624cac048a36f32d >> >> Lines 57 and 70 appear to be the culprit: >> >> https://github.com/sdwfrost/Gillespie.jl/blob/master/src/SSA.jl >> >> I've tried some devectorisation, but in my hackery, I appear to get side >> effects, where the argument x0 passed to the ssa function is modified. Any >> tips? >> >> Best >> Simon >> >
[julia-users] Re: Coveralls and coverage issues
Odd; refresh on Chrome wasn't refreshing properly. Fancy looking at my code speed question ;) ? On Thursday, July 21, 2016 at 3:49:01 PM UTC+1, Chris Rackauckas wrote: > > Click refresh when you're on the repo readme? It updated on my screen, > refresh to make sure you're not displaying the site from cache. > > On Thursday, July 21, 2016 at 7:41:47 AM UTC-7, Simon Frost wrote: >> >> Dear Chris, >> >> Yes, I am an idiot ;) >> >> Any idea why the badge isn't updating? >> >> Best >> Simon >> >> On Thursday, July 21, 2016 at 9:06:51 AM UTC+1, Chris Rackauckas wrote: >>> >>> Look at the files it's trying to cover... it's DataFrames.jl :) >>> >>> I sent you a pull request to fix your travis.yml to be for your package. >>> >>> On Thursday, July 21, 2016 at 12:16:35 AM UTC-7, Simon Frost wrote: >>>> >>>> Dear All, >>>> >>>> I'm trying to get code coverage working, but despite having some tests >>>> - at the moment, just running examples - I get 0% coverage >>>> >>>> http://github.com/sdwfrost/Gillespie.jl >>>> >>>> Is this because I'm just using 'include' in runtests.jl? >>>> >>>> Best >>>> Simon >>>> >>>
[julia-users] Performance issues with stochastic simulation
Dear All, I'm having some issues with code speed for some Gillespie type simulations. The toy model is described here: http://phylodynamics.blogspot.co.uk/2013/06/comparing-performance-of-r-and-rcpp-for.html http://phylodynamics.blogspot.co.uk/2013/06/an-sir-model-in-julia.html I get good performance with my vanilla Julia code, but a more generic implementation is slower: http://github.com/sdwfrost/Gillespie.jl The gist is here: https://gist.github.com/sdwfrost/1b4bce19faf2d7b8624cac048a36f32d Lines 57 and 70 appear to be the culprit: https://github.com/sdwfrost/Gillespie.jl/blob/master/src/SSA.jl I've tried some devectorisation, but in my hackery, I appear to get side effects, where the argument x0 passed to the ssa function is modified. Any tips? Best Simon
[julia-users] Re: Coveralls and coverage issues
Dear Chris, Yes, I am an idiot ;) Any idea why the badge isn't updating? Best Simon On Thursday, July 21, 2016 at 9:06:51 AM UTC+1, Chris Rackauckas wrote: > > Look at the files it's trying to cover... it's DataFrames.jl :) > > I sent you a pull request to fix your travis.yml to be for your package. > > On Thursday, July 21, 2016 at 12:16:35 AM UTC-7, Simon Frost wrote: >> >> Dear All, >> >> I'm trying to get code coverage working, but despite having some tests - >> at the moment, just running examples - I get 0% coverage >> >> http://github.com/sdwfrost/Gillespie.jl >> >> Is this because I'm just using 'include' in runtests.jl? >> >> Best >> Simon >> >
[julia-users] Coveralls and coverage issues
Dear All, I'm trying to get code coverage working, but despite having some tests - at the moment, just running examples - I get 0% coverage http://github.com/sdwfrost/Gillespie.jl Is this because I'm just using 'include' in runtests.jl? Best Simon
[julia-users] Re: Using callable types or FastAnonymous with Sundials
Dear Dan,
Changing the function to (say) g doesn't help. I want cvode to use J as a
function (as I'm overloading call); this is really just doing what
FastAnonymous does under the hood. I think it may be because Sundials
passes the derivatives as an argument to the function, and modifies them in
place.
Best
Simon
**
using Sundials
function cvode{f}(::Type{f}, y0::Vector{Float64}, t::Vector{Float64};
reltol::Float64=1e-4, abstol::Float64=1e-6)
neq = length(y0)
mem = Sundials.CVodeCreate(Sundials.CV_BDF, Sundials.CV_NEWTON)
flag = Sundials.CVodeInit(mem, cfunction(Sundials.cvodefun, Int32,
(Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Ref{Function})),
t[1], Sundials.nvector(y0))
flag = Sundials.CVodeSetUserData(mem, f)
flag = Sundials.CVodeSStolerances(mem, reltol, abstol)
flag = Sundials.CVDense(mem, neq)
yres = zeros(length(t), length(y0))
yres[1,:] = y0
y = copy(y0)
tout = [0.0]
for k in 2:length(t)
flag = Sundials.CVode(mem, t[k], y, tout, Sundials.CV_NORMAL)
yres[k,:] = y
end
Sundials.CVodeFree([mem])
return yres
end
function g(t, y, ydot)
ydot[1] = 0.1*(-72-y[1])+0.1*1.4*exp((y[1]+48)/1.4)+10
ydot[3] = 0.
ydot[2] = 0.
end
immutable J; end
call(::Type{J},t, y, ydot) = g(t, y, ydot)
t = [0.1:0.0001:1]
res = Sundials.cvode(g, [-60.0, 0.0, 0.0, t); # this works, passing a
Function type
res = cvode(J, [-60.0, 0.0, 0.0], t); # this gives ReadOnlyMemoryError
On Saturday, December 19, 2015 at 12:52:23 AM UTC, Dan wrote:
>
> The parameter for the `cvode` function is `f` and so is the function you
> want to use. These get confused, and it tries to use the "function" `J`
> instead. Changing the parameter name to something other than `{f}` should
> work
>
> On Wednesday, December 16, 2015 at 4:26:41 PM UTC+2, Simon Frost wrote:
>>
>> Dear Julia Users,
>>
>> I'm trying to speed up some code that employs passing functions as
>> arguments. One part of the code solves an ODE; if I use CVODE from
>> Sundials, and rewrite the function to accept callable types, I get a
>> ReadOnlyMemoryError - as I don't know the Sundials API, can someone help me
>> with where I'm going wrong? Code below.
>>
>> Best
>> Simon
>>
>> **
>>
>> using Sundials
>>
>> function cvode{f}(::Type{f}, y0::Vector{Float64}, t::Vector{Float64};
>> reltol::Float64=1e-4, abstol::Float64=1e-6)
>> neq = length(y0)
>> mem = Sundials.CVodeCreate(Sundials.CV_BDF, Sundials.CV_NEWTON)
>> flag = Sundials.CVodeInit(mem, cfunction(Sundials.cvodefun, Int32,
>> (Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Ref{Function})),
>> t[1], Sundials.nvector(y0))
>> flag = Sundials.CVodeSetUserData(mem, f)
>> flag = Sundials.CVodeSStolerances(mem, reltol, abstol)
>> flag = Sundials.CVDense(mem, neq)
>> yres = zeros(length(t), length(y0))
>> yres[1,:] = y0
>> y = copy(y0)
>> tout = [0.0]
>> for k in 2:length(t)
>> flag = Sundials.CVode(mem, t[k], y, tout, Sundials.CV_NORMAL)
>> yres[k,:] = y
>> end
>> Sundials.CVodeFree([mem])
>> return yres
>> end
>>
>> function f(t, y, ydot)
>> ydot[1] = 0.1*(-72-y[1])+0.1*1.4*exp((y[1]+48)/1.4)+10
>> ydot[3] = 0.
>> ydot[2] = 0.
>> end
>>
>> immutable J; end
>> call(::Type{J},t, y, ydot) = f(t, y, ydot)
>>
>> t = [0.1:0.0001:1]
>> res = cvode(J, [-60.0, 0.0, 0.0], t);
>>
>
[julia-users] Using callable types or FastAnonymous with Sundials
Dear Julia Users,
I'm trying to speed up some code that employs passing functions as
arguments. One part of the code solves an ODE; if I use CVODE from
Sundials, and rewrite the function to accept callable types, I get a
ReadOnlyMemoryError - as I don't know the Sundials API, can someone help me
with where I'm going wrong? Code below.
Best
Simon
**
using Sundials
function cvode{f}(::Type{f}, y0::Vector{Float64}, t::Vector{Float64};
reltol::Float64=1e-4, abstol::Float64=1e-6)
neq = length(y0)
mem = Sundials.CVodeCreate(Sundials.CV_BDF, Sundials.CV_NEWTON)
flag = Sundials.CVodeInit(mem, cfunction(Sundials.cvodefun, Int32,
(Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Ref{Function})),
t[1], Sundials.nvector(y0))
flag = Sundials.CVodeSetUserData(mem, f)
flag = Sundials.CVodeSStolerances(mem, reltol, abstol)
flag = Sundials.CVDense(mem, neq)
yres = zeros(length(t), length(y0))
yres[1,:] = y0
y = copy(y0)
tout = [0.0]
for k in 2:length(t)
flag = Sundials.CVode(mem, t[k], y, tout, Sundials.CV_NORMAL)
yres[k,:] = y
end
Sundials.CVodeFree([mem])
return yres
end
function f(t, y, ydot)
ydot[1] = 0.1*(-72-y[1])+0.1*1.4*exp((y[1]+48)/1.4)+10
ydot[3] = 0.
ydot[2] = 0.
end
immutable J; end
call(::Type{J},t, y, ydot) = f(t, y, ydot)
t = [0.1:0.0001:1]
res = cvode(J, [-60.0, 0.0, 0.0], t);
Re: [julia-users] Does anyone have a fork/branch of Sundials.jl that works on 0.4/0.5?
Dear Alex, Thanks for letting me know about the PR. As for the minor API differences, e.g. cvode in Sundials.jl returns a matrix, whereas (e.g.) ode23s returns a tuple with time steps and the dependent variable y separately. Not a big deal really, but it would be nice to avoid having to deal with two different types of output. Best, Simon On Wednesday, October 28, 2015 at 12:30:09 PM UTC, Alex wrote: > > Hi Simon, > > Christian Haargaard just created a PR ( > https://github.com/JuliaLang/Sundials.jl/pull/56), which means that there > is now (or should be) a working branch. > > May I ask what you mean by "minor API differences"? > > Best, > > Alex. > > On Saturday, 17 October 2015 10:58:10 UTC+2, Simon Frost wrote: > >> Dear Mauro, >> >> I think the fixes, at least for cvode, have already been done in one of >> the branches, but there isn't a PR yet. Apart from some minor API >> differences, which hopefully will be ironed out in the future, cvode in >> Sundials.jl is a lot faster than ode23s in ODE.jl. >> >> Best >> Simon >> >
Re: [julia-users] Does anyone have a fork/branch of Sundials.jl that works on 0.4/0.5?
Dear Mauro, I think the fixes, at least for cvode, have already been done in one of the branches, but there isn't a PR yet. Apart from some minor API differences, which hopefully will be ironed out in the future, cvode in Sundials.jl is a lot faster than ode23s in ODE.jl. Best Simon
[julia-users] Does anyone have a fork/branch of Sundials.jl that works on 0.4/0.5?
Dear All, Sundials.jl master isn't compatible with 0.4/0.5, I presume due to the changes in Ref/Ptr between 0.3 and 0.4. There's quite a few forks and branches, which look as though at least in some cases, this has been fixed - I care mostly about cvode. Rather than wade through installation and testing of a whole bunch of branches, it would be great if someone could point me in the right direction of a repo to clone, as I don't know when there's going to be a PR. Best wishes, Simon
[julia-users] Why is |> deprecated?
Dear All, I tried (and failed) to search for this, but I wanted to know why the forward pipe operator |> has been deprecated in favour of pipe(...,...,)? Best wishes, Simon
Re: [julia-users] Latest on wrapping C structs for use in Julia
Dear All, I've updated my gist of wrapping BEAGLE here: https://gist.github.com/sdwfrost/5c574857bd91648fb7ee The script runs fine, and the transition matrix qmat is being calculated correctly, but my log likelihood is still -Inf (rather than c. -84); is there anything obviously wrong with my Julia code? The C++ code I'm trying to mirror is here: https://github.com/sdwfrost/beagle-lib/blob/master/examples/standalone/hellobeagle/src/hello.cpp and the .h https://github.com/sdwfrost/beagle-lib/blob/master/libhmsbeagle/beagle.h Best Simon
Re: [julia-users] Latest on wrapping C structs for use in Julia
Dear Isaiah, Thanks - I noted before that changing types to immutable in some cases fixed things, but I didn't think that this would apply to BeagleOperation. I'll read over the new docs carefully. As for the second point, the signature in beagle.h is /** * @brief Set a state frequency buffer * * This function copies a state frequency array into an instance buffer. * * @param instance Instance number (input) * @param stateFrequenciesIndex Index of state frequencies buffer (input) * @param inStateFrequenciesState frequencies array (stateCount) (input) * * @return error code */ BEAGLE_DLLEXPORT int beagleSetStateFrequencies(int instance, int stateFrequenciesIndex, const double* inStateFrequencies); I'll follow up with the BEAGLE devs to see why your fix makes the script run, as the output (logL) isn't correct. Best Simon
Re: [julia-users] Latest on wrapping C structs for use in Julia
Dear Isaiah,
Thanks! I've put my code here:
https://gist.github.com/sdwfrost/5c574857bd91648fb7ee
The code currently dies here:
operations = [
BeagleOperation(3, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 0, 0, 1, 1),
BeagleOperation(4, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 2, 2, 3, 3)
]
ccall((:beagleUpdatePartials, "libhmsbeagle"),
Int,
(Cint,Ptr{BeagleOperation},Cint,Cint),
instance,operations,2,BEAGLE_OP_NONE)
I've tried various things to pass an array of structs as a pointer, but
nothing avoids a segfault.
Best
Simon
On Saturday, April 18, 2015 at 2:29:14 PM UTC+1, Isaiah wrote:
>
> Hi Simon,
>
> Is the use of an array here ... a way round converting things to
>> pointers?
>
>
> Yes, see:
>
> http://docs.julialang.org/en/release-0.3/manual/calling-c-and-fortran-code/#passing-pointers-for-modifying-inputs
>
> (you could do `retInfo = BeagleInstanceDetails(...)` and then pass
> `&retInfo` in ccall, but then you would not see the modifications)
>
> Also, is there any reason that resourceNumber::Int should be rather than
>> resourceNumber::Cint?
>
>
> Oversight on my part. When interfacing with C it is a good idea to use the
> C* aliases because they will take care of some cross-platform peculiarities.
>
> kills my IJulia kernel...are there any good ways to debug the external C
>> calls?
>
>
> The best way is to use a debugger, gdb or lldb (though admittedly a bit of
> a learning curve). Try running your test under the basic Julia shell; that
> should give you a backtrace at least. If you want to post your code
> somewhere (e.g. gist.github.com) I could have a look.
>
> Isaiah
>
>
>
> On Sat, Apr 18, 2015 at 7:12 AM, Simon Frost > wrote:
>
>> Dear Isaiah,
>>
>> Thanks, this has all really helped. Is the use of an array here:
>>
>> retInfo = [BeagleInstanceDetails(0,0,0,0,0)]
>>
>> a way round converting things to pointers? I had initial problems as I
>> defined
>>
>> retInfo = BeagleInstanceDetails(0,0,0,0,0)
>>
>> and didn't know how to convert it.
>>
>> Also, is there any reason that resourceNumber::Int should be rather than
>> resourceNumber::Cint?
>>
>> I've wrapped the complete example, but the last function call
>>
>> ccall((:beagleCalculateRootLogLikelihoods, "libhmsbeagle"),
>> Int,
>> (Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Cint,Ptr{Cdouble}),
>>
>> instance,rootIndex,categoryWeightIndex,stateFrequencyIndex,cumulativeScaleIndex,1,logL)
>>
>> kills my IJulia kernel...are there any good ways to debug the external C
>> calls?
>>
>> Best
>> Simon
>>
>
>
Re: [julia-users] Latest on wrapping C structs for use in Julia
Dear Isaiah,
Thanks, this has all really helped. Is the use of an array here:
retInfo = [BeagleInstanceDetails(0,0,0,0,0)]
a way round converting things to pointers? I had initial problems as I
defined
retInfo = BeagleInstanceDetails(0,0,0,0,0)
and didn't know how to convert it.
Also, is there any reason that resourceNumber::Int should be rather than
resourceNumber::Cint?
I've wrapped the complete example, but the last function call
ccall((:beagleCalculateRootLogLikelihoods, "libhmsbeagle"),
Int,
(Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Cint,Ptr{Cdouble}),
instance,rootIndex,categoryWeightIndex,stateFrequencyIndex,cumulativeScaleIndex,1,logL)
kills my IJulia kernel...are there any good ways to debug the external C
calls?
Best
Simon
[julia-users] Latest on wrapping C structs for use in Julia
Dear All, I'm trying to wrap the BEAGLE library for use in Julia: https://github.com/beagle-dev/beagle-lib I've used SWIG in the past, which works fine for Python: https://github.com/beagle-dev/beagle-lib/tree/master/examples/swig_python and I'm currently hacking my way through R: https://github.com/sdwfrost/beagle-lib/tree/master/examples/swig_r The API isn't big at all, but does require me to initialise some pointers to structs, which are filled in by reference. A toy example in C is here: https://github.com/sdwfrost/beagle-lib/blob/master/examples/standalone/hellobeagle/src/hello.cpp I've looked through the docs in Julia, as well as Clang.jl (I've made some wrappers previously using an old build, although I can't get Clang.jl working at the moment), and StrPack.jl, and I was wondering whether anyone had any pointers (no pun intended) to the current best practice for wrapping a library like this? I work best from examples, so any suggestions for a library with a similar API would be much appreciated. Best Simon
[julia-users] Re: Implementing Gillespie's Stochastic Simulation Algorithm
I just added a draft implementation along the lines of ODE.jl: https://github.com/sdwfrost/Gillespie.jl
[julia-users] Re: Implementing Gillespie's Stochastic Simulation Algorithm
Dear David, It probably would be easier to pass a function that, given the current states and the parameters, returned a Float64 vector of rates. The main issue that I have is that eval can't see all the variables defined in the function scope, so my code runs fine if run globally, but not as a function. If I go along that route, it would probably be best to define an appropriate model type, rather like Tim Vaughan's MASTER program: https://github.com/CompEvol/MASTER/wiki/Trajectory-specification I'll tweak my code to work with functions, and push to github when I get a chance. Best Simon
[julia-users] Implementing Gillespie's Stochastic Simulation Algorithm
Dear All,
I'm implementing Gillespie's direct method for stochastic simulation:
http://en.wikipedia.org/wiki/Gillespie_algorithm
I'm loosely following the API for the now-orphaned R package GillespieSSA:
http://www.jstatsoft.org/v25/i12
http://artax.karlin.mff.cuni.cz/r-help/library/GillespieSSA/html/ssa.html
In Julia, the analogous function call is something like the following.
function
ssa(x0::Vector{Int64},a::Vector{Expr},nu::Matrix{Float64,2},parms::Expr,tf::Float64,ignore_negative_state::Bool,console_interval::Int64,census_interval::Int64,verbose::Bool,max_wall_time::Float64)
However, as I've never done metaprogramming in Julia before, I was
wondering whether anyone had any input to make the model definition as
clean as possible (i.e. without cluttering the syntax with symbols, while
not polluting the environment)?
Best
Simon
[julia-users] Re: ODEs: ode23 works but ode45 doesn't?
Dear Alasdair,
There is all sorts of oddness with the ODE API, which appears to be in a
bit of flux at the moment, so I haven't filed any issues. For example,
ode23 really should only accept a span, but will accept a vector of any
length. It would make sense to use types to enforce this; perhaps using
(F::Function, tfinal::{T}, y0::AbstractVector{T})
for a final time
(F::Function, tspan::(T,T), y0::AbstractVector{T})
for a span
(F::Function, tvec::AbstractVector{T}, y0::AbstractVector{T})
for pre-specified times, and so on.
Also, ode23 will accept an Array{Float64,1} for initial values, but ode4,
which I'm currently using as it allows a fixed set of times to be specified
out of the box (needed as I'm fitting an ODE model to data) accepts a
column vector Array{Float64,2}.
Best
Simon
On Tuesday, January 7, 2014 11:04:18 AM UTC, Alasdair McAndrew wrote:
>
> Actually, I've looked at
> https://github.com/JuliaLang/ODE.jl/blob/master/test/test_ode.jl and
> attempted:
>
> t,x =
> ode.ode45((t,x)->[-beta*x[1]*x[2],beta*x[1]*x[2]-gamma*x[2],gamma*x[2]],[0,14],[760.,3.,0.]);
>
> where x[1], x[2] and x[3] are S, I, R respectively, and the parameters
> beta and gamma were defined earlier. This seems to work well. I suppose
> there's a way of using S, I, R instead of x[1], x[2], x[3] here...
>
> Sundials may well be the perfect solution, but I couldn't get it to work,
> and the documentation is daunting, to say the least!
>
> On Tuesday, January 7, 2014 9:58:54 PM UTC+11, Ivar Nesje wrote:
>>
>> This seems likely to be caused by changes in Julia. The fix will probably
>> involve a bug report to https://github.com/JuliaLang/ODE.jl and some
>> code reading and testing. If you see that the commit count on that
>> repository is 11 it is quite possible that it does not receive the care ODE
>> solvers deserve in a language like Julia.
>>
>> I am not sure, but Sundials.jl
>> might
>> be an option, but when I tried it last time it was no automatic script for
>> compiling and installing the dependencies.
>>
>>
>>
>> kl. 11:23:33 UTC+1 tirsdag 7. januar 2014 skrev Alasdair McAndrew
>> følgende:
>>>
>>> I have copied the SIR model here:
>>> http://phylodynamics.blogspot.co.uk/2013/07/differential-equation-modeling-with.htmlwhich
>>> works fine with ode23. But ode45, which would seem to have the same
>>> syntax, throws an error:
>>>
>>> julia> sol = ode.ode45((t,x)->SIR(t,x,p),t,inits);
>>> ERROR: InexactError()
>>> in setindex! at array.jl:471
>>> in oderkf at /opt/julia/usr/share/julia/site/v0.3/ODE/src/ODE.jl:217
>>> in ode45_dp at /opt/julia/usr/share/julia/site/v0.3/ODE/src/ODE.jl:277
>>>
>>> Does anybody know what's going on here, and how I can overcome it? In
>>> other words, how do I use ode45? Thanks!
>>>
>>>
>>>
