From 8f7276823543630333d1aaba6caf54ced1cbc2cc Mon Sep 17 00:00:00 2001 From: AleMorales Date: Sat, 28 Feb 2015 14:03:22 +0100 Subject: [PATCH] Separate into different functions the calculations of time derivatives and observed variables. Modify inputs in-place. Assume that forcing are generated using the Grid.jl package. --- src/codegeneration.jl | 103 ++++++++++++++++++++++-------------------- src/datatypes.jl | 3 +- src/simulation.jl | 39 +++++----------- 3 files changed, 68 insertions(+), 77 deletions(-) diff --git a/src/codegeneration.jl b/src/codegeneration.jl index 2694685..be3ce51 100644 --- a/src/codegeneration.jl +++ b/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 diff --git a/src/datatypes.jl b/src/datatypes.jl index d821b45..df11a80 100644 --- a/src/datatypes.jl +++ b/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 diff --git a/src/simulation.jl b/src/simulation.jl index 8bc8ac1..c34ac15 100644 --- a/src/simulation.jl +++ b/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