Trying to plot freezes Julia?

[HelgavonLichtenstein]

Hello,

I am trying to create a plot using Plots.jl and have run some code but it has already been over 15 minutes with no response from Julia. The screen looks like:

Is this normal? It seems unusual to me. I am using Julia version 1.0 and have updated all packages.

My code (which has been checked by others in Using DifferentialEquations.jl for a simple food-chain ) is:

using DifferentialEquations

f = @ode_def_bare FoodChain begin
   du1 = u1*(1-u1)-(1/e)*y2*x2*u1^(1+q)/(u012^(1+q)+u1^(q+1))*u2
   du2 = y2*x2*u1^(1+q)/(u012^(1+q)+u1^(1+q))*u2-(1/e)*y3*x3*u2^(1+q)/(u023^(1+q)+u2^(1+q))*u3-x2*u2
   du3 = y3*x3*u2^(1+q)/(u023^(1+q)+u2^(1+q))*u3-x3*u3
end e y2 x2 q u012 y3 x3 u023
u0 = [0.01,0.01,0.01]
tspan = (0.0,10000.0)
p = (1,2,3,4,5,6,7,8)
prob = ODEProblem(f,u0,tspan,p)
sol=solve(prob,saveat=0.05)

using Plots
plot(sol,xlabel = "Time" ,ylabel = "Density", layout = 3,1)

I am trying to create a graph series like this, in Julia. This plot was created in R and took much less time.

It is very probable that I have written something wrong in the plotting instructions. I am new to this and want to learn by doing, but if Julia does not generate plots I cannot experiment

I was able to run some examples from http://docs.juliadiffeq.org/latest/tutorials/ode_example.html#Example-1-:-Solving-Scalar-Equations-1 but now Julia is not responding. What am I doing wrong in running Plot.jl?

[Sukera]

I think you’re missing a ‘)’:

sol=solve(prob,saveat=0.05

[HelgavonLichtenstein]

Good eye! It was a copy/paste error, I ran it with the ‘)’ in Julia.

[Sukera]

(Took some time to compile the DifferentialEquations package)

The layout should be a tuple:

plot(sol,xlabel = “Time” ,ylabel = “Density”, layout = (3,1))

I can only currently plot with UnicodePlots, but with that it works fine and displays a plot (whether the plot itself is correct I do not know).

[e3c6]

You could try the GR backend with rasterized output, Plots.gr(fmt=:png). This allows me to plot hundreds of thousands of points without performance issues.

[HelgavonLichtenstein]

Sorry, I don’t know what this means. I know GR is a package, and this suggests it works with Plots.jl?

[HelgavonLichtenstein]

Thank you, it produces a graph now, but alas, something else is amiss…

[stillyslalom]

Try the plotdensity kwarg: http://docs.juliadiffeq.org/latest/basics/plot.html#Density-1

plot(sol, plotdensity=1000)

[HelgavonLichtenstein]

Adding
plot(sol,xlabel = "Time" ,ylabel = "Density", plotdensity=1000, layout = (3,1))

Does not seem to make the graph change. I think the plot density should be 200 000 (tstart = 0, tmax = 10 000, increases in increments of 0.05) but this does not change the look of the graphs either.

[e3c6]

Plots.jl supports several backends that do the actual plotting. Instead of:

using Plots
plot(sol,xlabel = "Time" ,ylabel = "Density", layout = 3,1)

You can use the gr backend by doing this:

using Plots
Plots.gr(fmt=:png)
plot(sol,xlabel = "Time" ,ylabel = "Density", layout = 3,1)

In my experience gr handles large number of points very fast, but maybe this is not your issue.

[stillyslalom]

It doesn’t change the look of the graphs because the additional 199,000 points aren’t visible due to limited screen resolution. Your screen is only a few thousand pixels wide, and plot density beyond that is overkill.

The solver’s behaving correctly–is the equation itself free of typos? Do you have the original R code handy?

[HelgavonLichtenstein]

The original R code is

library("deSolve")
library("EMD") 

model.rma <-function(t, x, parms)
{
	b1<-x[1]
	b2<-x[2]
  b3<-x[3]
					
	db1 <- b1*(1-b1)-(1/e)*y2*x2*b1^(1+q)/(b012^(1+q)+b1^(q+1))*b2
	db2 <- y2*x2*b1^(1+q)/(b012^(1+q)+b1^(1+q))*b2-(1/e)*y3*x3*b2^(1+q)/(b023^(1+q)+b2^(1+q))*b3-x2*b2
  db3 <- y3*x3*b2^(1+q)/(b023^(1+q)+b2^(1+q))*b3-x3*b3

  return(list(c(db1,db2,db3)))
}

tmax <- 10000
b1.0 <- 0.01    # initial values
b2.0 <- 0.01
b3.0 <- 0.01

e <- 1          # assimilation efficiency
y2 <- 2.009     # maximum ingestion rate (relative to metabolic rate) of species 2
x2 <- 0.4       # metabolic rate of species 2
y3 <- 5         # maximum ingestion rate (relative to metabolic rate) of species 3
x3 <- 0.08      # metabolic rate of species 3
q <-  0         # Hill-Exponent = 1 + q
b012 <- 0.16129 # half-saturation density
b023 <- 0.5     # half-saturation density

#parms <-c(e,y2,x2,q,b012,y3,x3,b023)
xstart <-c(b1=b1.0, b2=b2.0, b3=b3.0)

times <-seq(0, tmax, 0.05)

out1 <-as.data.frame(lsoda(xstart, times, model.rma, parms))

par(mfrow=c(3,1))
par(las=1)
par(pty="m")  # maximises plotting area (instead of squares)

# plot time series
plot(
  c(out1$time,out1$time,out1$time),
  c(out1$b1,out1$b2,out1$b3),
  xlab="time",
  ylab="top predator density",
  type="n"
  )
lines(out1$time,out1$b3,col="red")

plot(
  c(out1$time,out1$time,out1$time),
  c(out1$b1,out1$b2,out1$b3),
  xlab="time",
  ylab="consumer density",
  type="n"
  )
lines(out1$time,out1$b2,col="orange")

plot(
  c(out1$time,out1$time,out1$time),
  c(out1$b1,out1$b2,out1$b3),
  xlab="time",
  ylab="prey density",
  type="n"
  )
lines(out1$time,out1$b1,col="darkgreen")

[HelgavonLichtenstein]

Ah I see! This will be useful for me in the future once I start dealing with more data. Thank you.

[Sukera]

It’s definitely a resolution/selective plotting issue with the plot output - when I plot only the first solution (plot(sol[1,:],xlabel = "Time" ,ylabel = "Density")) I get the reverse problem with the line being up instead of down below, i.e. seemingly only plotting local maxima instead of local minima. This might be a restriction of me plotting in the console though. Can you try only plotting one solution?

[HelgavonLichtenstein]

It looks like

[stillyslalom]

You didn’t use the same ODE parameters, so no wonder the plots look different. No need to set the timestep so small, either: DiffEq.jl does a good job with adaptive timestepping.

p = (1, 2.009, 0.4, 0, 0.16129, 5, 0.08, 0.5)
prob = ODEProblem(f,u0,tspan,p)
sol=solve(prob)

[Sukera]

That ought to do it! I suspected the data too now. Here’s some stats from before:

julia> mean(sol[1,:])
0.999537010691792

julia> median(sol[1,:])
0.9999999999999999

and after your code:

julia> mean(sol[1,:])
0.6955485905767033

julia> median(sol[1,:])
0.9031611072330993

[HelgavonLichtenstein]

Oh my. Yes you are correct, somehow parameters were switched. Thank you and sorry for making a silly mistake.