# adaptive step
# with noise
@everywhere using Distributions
@everywhere using StatsBase
@everywhere using DifferentialEquations
@everywhere using Interpolations
include("../src/SMCABC_RKHS_RJ.jl")
include("Models_rj.jl")
cd(@__DIR__)




# data generation for test---------
@everywhere function generate_dataset()
    timepoints = linspace(0,5, 500)
    μ = 0
    σ = 10
    kernel = Normal(μ, σ)
    x=rand(kernel, length(timepoints));

    itp = interpolate(x, BSpline(Linear()), OnGrid());
    sitp = Interpolations.scale(itp, timepoints);
    model = model_generator(0, sitp)

    y0 = [0.0,0.0]
    tspan = (timepoints[1], timepoints[end])
    params = [1., 0.05, 50., 1e+3]

    alg = CVODE_BDF(max_hnil_warns = 0, max_error_test_failures = 0)
    prob = ODEProblem(model, y0, tspan, params)
    res = solve(prob, alg, saveat=timepoints, maxiters=1e7, abstol = 1e-8, reltol = 1e-8)
    return [hcat(res.u...)']
end
@everywhere dataset = generate_dataset()
#--------------------------------
srand(9357478685)

@everywhere N_model = 2
N_particle = 5000
ϵ₀ = Inf
Δϵ_tol = 1e-8
α = 0.98
target = 1-α
N_sample = 100
Δϵ_window = 20
save_interval = 50
parallel = true
output = joinpath(@__DIR__, "result_RKHS_rj2")
@everywhere μ = 0
@everywhere σ = 10
@everywhere timepoints = linspace(0,5, 500)
@everywhere lbs_set = [[0.1, 0.005, 5., 10^2],[0.1, 0.005, 5., 10^2, 10^4]]
@everywhere ubs_set = [[10, 0.5, 500., 10^4],[10, 0.5, 500., 10^4, 10^6]]
@everywhere kernel = Normal(μ, σ)
@everywhere x=rand(kernel, length(timepoints));
@everywhere itp = interpolate(x, BSpline(Linear()), OnGrid());
@everywhere sitp = Interpolations.scale(itp, timepoints);
@everywhere models = [model_generator(i-1, sitp) for i in 1:N_model]
@everywhere obj_funcs = [f(param) = obj_func(dataset, models[i], param, lbs_set[i], ubs_set[i]) for i in 1:N_model]

run_rjsmc(N_particle, lbs_set, ubs_set, obj_funcs, ϵ₀, Δϵ_tol, α, Δϵ_window, output;
          save_interval=save_interval, target=target, N_sample=N_sample, parallel=parallel)

# plot -----------------------------------------------
# est_params = readcsv("result_rj/est_param_M2_655.csv")
#
# nrow=5
# ncol=1
# figsize=(5,10)
# fig, axes = subplots(nrow, ncol, figsize=figsize);
#
# for i in 1:5
#     axes[i][:hist](est_params[:,i], bins=100, normed=true, color="b");
#     axes[i][:axvline](params[i], color="r");
# end
# fig[:savefig](joinpath(output, "param_dist.png"), dpi=300);
# PyPlot.plt[:close]()
#
#
# y0 = [0.0,0.0]
# tspan = (timepoints[1], timepoints[end])
# params = [1., 0.05, 50., 1e+3, 1e+5]
#
# trajectories = zeros(size(est_params)[1], length(timepoints), 2)
# for i in 1:size(est_params)[1]
#     param = est_params[i,:]
#     prob = ODEProblem(model, y0, tspan, param)
#     res_temp = solve(prob, alg, saveat=timepoints, maxiters=1e7,
#                 abstol = 1e-8, reltol = 1e-8)
#     trajectories[i,:,:] = hcat(res_temp.u...)'
# end
#
# ave = mean(trajectories, Dims(1))
# sd = std(trajectories, Dims(1))
# y1=ave.-2sd
# y2=ave.+2sd
#
# prob = ODEProblem(model, y0, tspan, params)
# res = solve(prob, alg, saveat=timepoints, maxiters=1e7,
#             abstol = 1e-8, reltol = 1e-8)
#
# plot(timepoints, vec(ave[:,:,1]), color="red", ls="--", lw=1);
# fill_between(timepoints, vec(y1[:,:,1]), vec(y2[:,:,1]));
#
# plot(timepoints, res[1, :], lw=1, color="k", ls="--");
#
# savefig(joinpath(output, "predict.png"), dpi=300);
# PyPlot.plt[:close]()
