Separate into different functions the calculations of time derivative…

…s and observed variables. Modify inputs in-place. Assume that forcing are generated using the Grid.jl package.

src/codegeneration.jl

@@ -119,49 +119,49 @@ end
 ####################################################################################
 ####################################################################################
 
-# Create a return line when the output is a tuple of 2 arrays
-function create_return_line_julia(states, observed)
- return_line = "(["
- for i in states
- if i != states[end]
- return_line = return_line * i * ","
- else
- return_line = return_line * i
- end
- end
- return_line = return_line * "], ["
- for i in observed
- if i != observed[end]
- return_line = return_line * i * ", "
- else
- return_line = return_line * i
- end
+# Create the function in Julia on the Equations section of the model Dict
+function create_derivatives_julia(model::OdeSorted, states, name)
+ code = ""
+ names_derivatives = ["d_"*i*"_dt" for i in states]
+ print(names_derivatives)
+ for level in 1:length(model.SortedEquations)
+ for (lhs, rhs) in model.SortedEquations[level]
+ if level == 1
+ code *= string(rhs.Expr) * "\n"
+ elseif in(lhs, names_derivatives)
+ c = findin(names_derivatives, [lhs])
+ code *= "ydot[$c] = " * string(rhs.Expr) * "\n"
+ else
+ code *= lhs * " = " * string(rhs.Expr) * "\n"
+ end
  end
- return_line = return_line * "])"
+ end
+ # Construct the function
+ paste("\n","@inbounds function derivatives_$name(time::Float64, states::Vector{Float64},
+ params::Vector{Float64}, forcs::Vector{Float64}, ydot::Vector{Float64})\n",
+ code,"nothing \n end")
 end
 
 # Create the function in Julia on the Equations section of the model Dict
-function create_function_julia(model::OdeSorted, observed, name)
+function create_observed_julia(model::OdeSorted, observed, name)
  code = ""
+ c = 1
  for level in 1:length(model.SortedEquations)
  for (lhs, rhs) in model.SortedEquations[level]
  if level == 1
- code *= string(rhs.Expr) * "\n"
+ code *= string(rhs.Expr) * "\n"
+ elseif in(lhs, observed)
+ code *= "$lhs = " * string(rhs.Expr) * "\n"
+ code *= "obs[$c] = $lhs\n"
+ c += 1
  else
-  code *= lhs * " = " * string(rhs.Expr) * "\n"
+ code *= lhs * " = " * string(rhs.Expr) * "\n"
  end
  end
  end
- # Determine what the time derivatives are
- time_derivatives = String[]
- for i in collect(keys(model.States))
- push!(time_derivatives, "d_"*i*"_dt")
- end
- # Return line
- return_line = create_return_line_julia(time_derivatives,observed)
- # Return the output
- code = replace(code, "1.0 *", "")
- return paste("\n","@inbounds function $name(time::Float64, states::Array{Float64,1}, params::Array{Float64,1}, forcs::Array{Float64,1})\n", code, return_line,"end")
+ paste("\n","@inbounds function observed_$name(time::Float64, states::Vector{Float64},
+ params::Vector{Float64}, forcs::Vector{Float64}, obs::Vector{Float64})\n",
+ code, "nothing \n end")
 end
 
 function generate_code_Julia!(ode_model::OdeSource; unit_analysis = true, name = "autogenerated_model", file = "autogenerated_model", jacobian = false, sensitivities = false)
@@ -179,22 +179,6 @@ function generate_code_Julia!(ode_model::OdeSource; unit_analysis = true, name =
  # Sort the equations
  sorted_model = sort_equations(ode_model);
 
- # Created compressed model (only if Jacobian or Sensitivities are required!)
- jacobian_function = "jacobian_$name() = nothing"
- sensitivity_function = "() -> ()"
- sensitivity_jacobian_function = "() -> ()"
- if jacobian || sensitivities
- compressed_model = compress_model(sorted_model, level = 2)
- jacobian && (jacobian_function = generate_jacobian_function(compressed_model, name))
- sensitivities && ((sensitivity_function, sensitivity_jacobian_function) = generate_extended_system(compressed_model, name))
- end
-
- # Check the units
- unit_analysis && check_units(sorted_model)
-
- # Generate the rhs function (compressed at level 2)
- model_function = create_function_julia(sorted_model,observed,name)
-
  # Create the default arguments
  named_states = OrderedDict{String, Any}()
  for (key,val) in sorted_model.States
@@ -210,8 +194,27 @@ function generate_code_Julia!(ode_model::OdeSource; unit_analysis = true, name =
  c += 1
  forcings[key] = (float(value.Time), float(value.Value)*value.Units.f)
  end
+
+ # Created compressed model (only if Jacobian or Sensitivities are required!)
+ jacobian_function = "jacobian_$name() = nothing"
+ sensitivity_function = "() -> ()"
+ sensitivity_jacobian_function = "() -> ()"
+ if jacobian || sensitivities
+ compressed_model = compress_model(sorted_model, level = 2)
+ jacobian && (jacobian_function = generate_jacobian_function(compressed_model, name))
+ sensitivities && ((sensitivity_function, sensitivity_jacobian_function) = generate_extended_system(compressed_model, name))
+ end
+
+ # Check the units
+ unit_analysis && check_units(sorted_model)
+
+ # Generate the rhs function
+ model_function = create_derivatives_julia(sorted_model,collect(keys(named_states)), name)
+ # Generate function with the observed values
+ observed_function = create_observed_julia(sorted_model,observed,name)
+
  write_model_Julia!(named_states,named_parameters, forcings, observed,
- model_function,jacobian_function, name, file)
+ model_function,observed_function, jacobian_function, name, file)
  nothing
 end
 
@@ -220,6 +223,7 @@ function write_model_Julia!(States::OrderedDict{String, Any},
  Forcings::OrderedDict{String, Any},
  Observed::Array{String, 1},
  Model::String,
+ Observed_model::String,
  Jacobian::String,
  name::String,
  file::String)
@@ -241,8 +245,9 @@ function write_model_Julia!(States::OrderedDict{String, Any},
  end
  println(f, "Observed = $Observed")
  println(f, "$Model")
+ println(f, "$(Observed_model)")
  println(f, "$Jacobian")
- println(f, "ODEDSL.DataTypes.OdeModel(States, Parameters, Forcings, Observed, $name, jacobian_$name)")
+ println(f, "ODEDSL.DataTypes.OdeModel(States, Parameters, Forcings, Observed, derivatives_$name, observed_$name, jacobian_$name)")
  println(f, "end")
  close(f)
  nothing

src/datatypes.jl

@@ -283,8 +283,9 @@ type OdeModel
  States::OrderedDict{String, Any}
  Parameters::OrderedDict{String, Any}
  Forcings::OrderedDict{String, Any}
- Observed::Array{String, 1}
+ Observed_names::Array{String, 1}
  Model::Function
+ Observed::Function
  Jacobian::Function
 end
 

src/simulation.jl

@@ -17,28 +17,13 @@ using Grid
  cols::Ptr{Ptr{Float64}}
 end
 
-@inbounds function linear_interpolate(x::Vector{Float64}, y::Vector{Float64}, xout::Float64)
- if xout >= x[end]
- return y[end]
- elseif xout <= x[1]
- return y[1]
- else
- index = searchsortedfirst(x, xout) - 1
- inty = (y[index + 1] - y[index])/(x[index + 1] - x[index])*(xout -x[index]) + y[index]
- return inty
- end
-end
-
 # We must convert ydot into an array and assign element-by-element
 # This is because we are modifying by-reference an NVector in Sundials...
 @inbounds function interface_ode(t::Sundials.realtype, y::Sundials.N_Vector, ydot::Sundials.N_Vector, user_data::Array{Any, 1})
  y = Sundials.asarray(y)
  ydot = Sundials.asarray(ydot)
  forcings = Float64[user_data[3][i][t] for i in 1:length(user_data[3])]
- derivatives = user_data[1](t, y, user_data[2], forcings)[1]
- for i in 1:length(derivatives)
- ydot[i] = derivatives[i]
- end
+ user_data[1](t, y, user_data[2], forcings, ydot)
  return int32(0)
 end
 
@@ -67,6 +52,7 @@ end
  forcings_data,
  settings::Dict{Any,Any},
  model::Function,
+ observer::Function,
  jacobian::Function)
  neq = length(states)
 
@@ -167,14 +153,12 @@ end
  # Make a first call to the model to check that everything is ok and retrieve the number of observed variables
  forcings = zeros(Float64, length(forcings_data))
  for i in 1:length(forcings) forcings[i] = forcings_data[i][times[1]] end
- first_call = model(times[1], states, parameters, forcings)
-
- # Fill up the output matrix with the values for the initial time
+ observed = zeros(Float64, settings["nobs"])
+ observer(times[1], states, parameters, forcings, observed)
 
- nder = length(first_call[2])
  # Use time as column -> Because of column-ordered
- output = zeros(Float64, (neq + nder + 1, length(times)))
- output[:,1] = [times[1], states, first_call[2]]
+ output = zeros(Float64, (neq + settings["nobs"] + 1, length(times)))
+ output[:,1] = [times[1], states, observed]
 
  # Main time loop. Each timestep call cvode. Handle exceptions and fill up output
  t = times[1]
@@ -217,11 +201,11 @@ end
  end
  end
  # If we have observed variables we call the model function again
- if nder > 0
+ if settings["nobs"] > 0
  for i in 1:length(times)
  for h in 1:length(forcings) forcings[h] = forcings_data[h][times[i]] end
- model_call = model(times[i], vec(output[2:(neq + 1), i]), parameters, forcings)
- for h in 1:nder output[neq + 1 + h, i] = model_call[2][h] end
+ observer(times[i], vec(output[2:(neq + 1), i]), parameters, forcings, observed)
+ for h in 1:settings["nobs"] output[neq + 1 + h, i] = observed[h] end
  end
  end
 
@@ -242,14 +226,15 @@ function simulate(model::OdeModel, times, settings)
  !haskey(settings, "maxnonlin") && (settings["maxnonlin"] = 12)
  !haskey(settings, "maxconvfail") && (settings["maxconvfail"] = 12)
  !haskey(settings, "jacobian") && (settings["jacobian"] = false)
+ settings["nobs"] = length(model.Observed_names)
  # The following conversions should be performed at storage time...
  forcs = [InterpIrregular(i[2][1], i[2][2], BCnearest, InterpLinear) for i in model.Forcings]
  states = Float64[i[2] for i in model.States]
  parameters = Float64[i[2] for i in model.Parameters]
  simulation = convert(DataFrame, simulate(times, states, parameters,
- forcs, settings, model.Model, model.Jacobian))
+ forcs, settings, model.Model, model.Observed, model.Jacobian))
  # Use names for indexing...
  names!(simulation, convert(Array{Symbol, 1}, [symbol("$i") for i in
- ["time", collect(keys(model.States)), model.Observed]]))
+ ["time", collect(keys(model.States)), model.Observed_names]]))
  return simulation
 end