tanhevg / GpABC.jl

MIT License
54 stars 15 forks source link

3 genes example in the manual: SimulatedABCRejection() fails #72

Closed EvaJanouskova closed 2 years ago

EvaJanouskova commented 2 years ago

I tried to run the first In[1]--In[4] of the example.

If I run only In[1] and In[2] and then try to run simulator_function() it leads to an

ERROR: MethodError: Cannotconvertan object of type ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#ODE_3GeneReg#7", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, RK4, OrdinaryDiffEq.InterpolationData{ODEFunction{true, var"#ODE_3GeneReg#7", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.RK4Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats} to an object of type Float64 Closest candidates are: convert(::Type{T}, ::Static.StaticFloat64{N}) where {N, T<:AbstractFloat} at ~/.julia/packages/Static/pkxBE/src/float.jl:26 convert(::Type{T}, ::LLVM.GenericValue, ::LLVM.LLVMType) where T<:AbstractFloat at ~/.julia/packages/LLVM/YSJ2s/src/execution.jl:39 convert(::Type{T}, ::LLVM.ConstantFP) where T<:AbstractFloat at ~/.julia/packages/LLVM/YSJ2s/src/core/value/constant.jl:111 ...

I made the simulator_function() run only, when I changed the function GeneReg() to return Obs only (without the type).

Then I tried to run In[3] and In[4]. Nevertheless, the SimulatedABCRejection() fails.

Here is the exact code I used:

`

~/GIT/GIT_Schistoxpkg.jl_1.2.15/src/ABCexample.jl

#

ABC settings

# using GpABC using OrdinaryDiffEq using Distances using Distributions using Plots using StatsBase using Printf pyplot()

true_params = [2.0, 1.0, 15.0, 1.0, 1.0, 1.0, 100.0, 1.0, 1.0, 1.0] # nominal parameter values priors = [Uniform(0.2, 5.), Uniform(0.2, 5.), Uniform(10., 20.), Uniform(0.2, 2.), Uniform(0.2, 2.), Uniform(0.2, 2.), Uniform(75., 125.), Uniform(0.2, 2.), Uniform(0.2, 2.), Uniform(0.2, 2.)] param_indices = [1, 3, 9] #indices of the parameters we want to estimate priors = priors[param_indices]

#

ODE solver settings

# Tspan = (0.0, 10.0) x0 = [3.0, 2.0, 1.0] solver = RK4() saveat = 0.1

#

Returns the solution to the toy model as solved by OrdinaryDiffEq

# GeneReg = function(params::AbstractArray{Float64,1}, Tspan::Tuple{Float64,Float64}, x0::AbstractArray{Float64,1}, solver::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm, saveat::Float64)

if size(params,1) != 10 throw(ArgumentError("GeneReg needs 10 parameters, $(size(params,1)) were provided")) end

function ODE_3GeneReg(dx, x, par, t) dx[1] = par[1]/(1+par[7]x[3]) - par[4]x[1] dx[2] = par[2]par[8]x[1]/(1+par[8]x[1]) - par[5]x[2] dx[3] = par[3]par[9]x[1]par[10]x[2]./(1+par[9]x[1])./(1+par[10]x[2]) - par[6]*x[3] end

prob = ODEProblem(ODE_3GeneReg, x0 ,Tspan, params) Obs = solve(prob, solver, saveat=saveat)

return Obs ######## #instead of: return Array{Float64, 2}(Obs) end

#

A function that simulates the model

# function simulator_function(var_params) params = copy(true_params) params[param_indices] .= var_params GeneReg(params, Tspan, x0, solver, saveat) end

simulator_function([2.0, 15.0, 1.0])

#

Get reference data and plot it

# reference_data = simulator_function(true_params[param_indices]) plot(reference_data', xlabel="t", ylabel="C(t)", linewidth=2, labels=["u1(t)" "u2(t)" "u3(t)"])

#

Simulation

# n_particles = 2000 threshold = 1.0 sim_result = SimulatedABCRejection(reference_data, simulator_function, priors, threshold, n_particles; max_iter=convert(Int, 2e6), write_progress=false) plot(sim_result) `

From the SimulatedABCRejection() I get an ERROR: DimensionMismatch("parent has 101 elements, which is incompatible with size (303,)") Stacktrace: [1] _throw_dmrs(n::Int64, str::String, dims::Tuple{Int64}) @ Base ./reshapedarray.jl:181 [2] _reshape @ ./reshapedarray.jl:176 [inlined] [3] reshape @ ./reshapedarray.jl:112 [inlined] [4] reshape @ ./reshapedarray.jl:116 [inlined] [5] keep_all_summary_statistic(data::ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#ODE_3GeneReg#7", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, RK4, OrdinaryDiffEq.InterpolationData{ODEFunction{true, var"#ODE_3GeneReg#7", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.RK4Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats}) @ GpABC ~/.julia/packages/GpABC/o0EN1/src/abc/summary_stats.jl:107 [6] SimulatedABCRejection(reference_data::ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, var"#ODE_3GeneReg#7", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, RK4, OrdinaryDiffEq.InterpolationData{ODEFunction{true, var"#ODE_3GeneReg#7", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.RK4Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats}, simulator_function::Function, priors::Vector{Uniform{Float64}}, threshold::Float64, n_particles::Int64; summary_statistic::String, distance_function::Euclidean, max_iter::Int64, kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:write_progress,), Tuple{Bool}}}) @ GpABC ~/.julia/packages/GpABC/o0EN1/src/abc/simulation.jl:54 [7] top-level scope @ ~/GIT/GIT_Schistoxpkg.jl_1.2.15/src/ABCexample.jl:76

Maybe smt has been changed in julia or even in the SimulatedABCRejection() since the last example update (June 2020)?

tanhevg commented 2 years ago

Hello!

Always good to get user feedback.

Sorry about this. Turns out in the most recent release I have updated the dependency versions but have forgotten to update the examples. I have just pushed another update, so this should be fixed in master now. To be precise, the last line in GeneReg function should read:

return hcat(Obs.u...)

Could you please confirm that this is now resolved on your end and close the issue.

EvaJanouskova commented 2 years ago

I'm sorry for the late reply. The example works now on my side. Thanks!