
@everywhere module SingleParticle
using StatsBase
using KernelDensity
using Distributions

export Particle, init_particle, init_particle2, bounds_check, bounds_check2, judge, MCstep, MCstep, MCstep_parallel, MCstep_parallel2

mutable struct Particle
    x_present::Array{Float64, 1}
    obj_func::Function # return relative weight
    obj_func_val::Float64
    weight::Float64
    cov::Array{Float64, 2}
end

function init_particle(lbs, ubs, obj_func, weight)
    """
    obj_func: must return -1 as a flag if has an error.
    """
    @assert length(lbs) == length(ubs)
    num_param = length(lbs)
    obj_func_val = -1
    kernel = Uniform(0, 10)
    x_present = zeros(length(lbs))   # generate parameter
    while obj_func_val < 0
        for i in 1:num_param
            x_present[i] = rand(kernel)
        end
        obj_func_val = obj_func(x_present)
    end
    return Particle(x_present, obj_func, obj_func_val, weight, eye(num_param))
end

function bounds_check(x_next)
    flag = true
    if x_next < 0 || x_next > 10
        flag = false
    end
    return flag
end

function judge(particle, proposed, F_obj_new, prior, ϵ)
    flag = false
    prob_present = prod([pdf(p, particle.x_present[i]) for (i, p) in enumerate(prior)])
    prob_proposed = prod([pdf(p, proposed[i]) for (i, p) in enumerate(prior)])
    MH_ratio = (F_obj_new < ϵ) * prob_proposed / ((particle.obj_func_val < ϵ) * prob_present)
    if rand() < MH_ratio
        flag = true
    else
        nothing
    end
    return flag
end

function MCstep(particle, prior, ϵ)
    """
    particle transition based on MH algorithm.
    """
    proposed = zeros(particle.x_present)
    kernel = MultivariateNormal(particle.x_present, particle.cov)
    F_obj_new = -1
    flag_MH = 0

    while F_obj_new < 0
        flag_bound = false
        while flag_bound == false
            proposed = rand(kernel)
            # proposed = np.random[:multivariate_normal](particle.x_present, np.array(particle.cov))
            for (i, x) in enumerate(proposed)
                flag_bound = bounds_check(x)
                if flag_bound == false
                    break
                end
            end
        end
        F_obj_new = particle.obj_func(proposed)
    end

    if judge(particle, proposed, F_obj_new, prior, ϵ) == true
        particle.x_present = proposed
        particle.obj_func_val = F_obj_new
        flag_MH = 1
    end

    return particle, flag_MH
end

function MCstep_parallel(particle, prior, ϵ)
    """
    particle transition based on MH algorithm.
    """
    paramset = particle.x_present
    obj_func_val = particle.obj_func_val
    proposed = zeros(particle.x_present)
    kernel = MultivariateNormal(particle.x_present, particle.cov)
    F_obj_new = -1
    flag_MH = 0

    while F_obj_new < 0
        flag_bound = false
        while flag_bound == false
            proposed = rand(kernel)
            # proposed = np.random[:multivariate_normal](particle.x_present, np.array(particle.cov))
            for (i, x) in enumerate(proposed)
                flag_bound = bounds_check(x)
                if flag_bound == false
                    break
                end
            end
        end
        F_obj_new = particle.obj_func(proposed)
    end

    if judge(particle, proposed, F_obj_new, prior, ϵ) == true
        paramset = proposed
        obj_func_val = F_obj_new
        flag_MH = 1
    end

    return paramset, obj_func_val, flag_MH
end


end
