Skip to content
Snippets Groups Projects
  • Bentriou Mahmoud's avatar
    5d886fc4
    The package becomes more meta to reach higher performance. · 5d886fc4
    Bentriou Mahmoud authored
    This commit groups the change operated to the creation of models and
    simulate function of a ContinuousTimeModel.
    The general idea is to create a concrete type and a simulate function
     per model creation by metaprogramming.
    - Now, ContinuousTimeModel is an abstract type. Each creation of a model
    defines a concrete type T <: ContinuousTimeModel by meta programming.
    - f! and isabsorbing ContinuousTimeModel fields are Symbols.
    - simulate(::ContinuousTimeModel) is run by multiple dispatch, according
    to the type of the model.
    
    Can't run the whole tests for now but unit/simulate_available_models.jl
    runs properly (i've updated the list of models in this commit), and I've
    manually checked in the repl that simulations run correctly (distributed
    / plots).
    5d886fc4
    History
    The package becomes more meta to reach higher performance.
    Bentriou Mahmoud authored
    This commit groups the change operated to the creation of models and
    simulate function of a ContinuousTimeModel.
    The general idea is to create a concrete type and a simulate function
     per model creation by metaprogramming.
    - Now, ContinuousTimeModel is an abstract type. Each creation of a model
    defines a concrete type T <: ContinuousTimeModel by meta programming.
    - f! and isabsorbing ContinuousTimeModel fields are Symbols.
    - simulate(::ContinuousTimeModel) is run by multiple dispatch, according
    to the type of the model.
    
    Can't run the whole tests for now but unit/simulate_available_models.jl
    runs properly (i've updated the list of models in this commit), and I've
    manually checked in the repl that simulations run correctly (distributed
    / plots).
square_wave_oscillator.jl 2.45 KiB

d=3
k=7
dict_var_square = Dict(:A => 1, :HIGH => 2, :LOW => 3)
dict_params_square = Dict(:p1 => 1, :w1 => 2, :p2 => 3, :w2 => 4, :p3 => 5, :w3 => 6, :d => 7)
l_tr_square = [:Tu, :t1, :t2, :t3]
p_square = [1, 0, 0, 0, 0, 0, 5.0]
x0_square = [1, 0, 1]
t0_square = 0.0
function SquareWave_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},
                        xn::Vector{Int}, tn::Float64, p::Vector{Float64})
    if p[1] == 0 && p[3] == 0 && p[5] == 0
        copyto!(xnplus1, xn)
        return nothing
    end
    # If in state LOW
    if xn[3] == 1
        xnplus1[1] = 10
        xnplus1[2] = 1
        xnplus1[3] = 0
        l_t[1] += p[7]
        l_tr[1] = :Tu
    # If in state HIGH
    else
        possible_transitions = (:t1, :t2, :t3)
        l_priority = (p[1], p[3], p[5])
        l_weights = (p[2], p[4], p[6])
        max_priority = maximum(l_priority)
        l_idx = zeros(Int, 0)
        for i = eachindex(l_priority)
            if l_priority[i] >= max_priority
                push!(l_idx, i)
            end
        end
        if length(l_idx) == 1
            transition = l_idx[1]
        else
            wsum = sum(l_weights[l_idx])
            b_inf = 0.0
            b_sup = l_weights[l_idx[1]]
            transition = 0
            u = rand()
            for i in l_idx
                if b_inf < wsum*u < b_sup
                    transition = i
                    break
                end
                @inbounds b_inf += l_weights[i]
                @inbounds b_sup += l_weights[i]
            end
        end
        xnplus1[1] = 1
        xnplus1[2] = 0
        xnplus1[3] = 1
        l_t[1] += 2^(transition-1) * p[7]
        l_tr[1] = possible_transitions[transition]
    end
end
isabsorbing_SquareWave(p::Vector{Float64}, xn::Vector{Int}) = (p[1] == 0 && p[3] == 0 && p[5] == 0)
g_square_wave = [:A, :HIGH, :LOW]

@everywhere @eval $(MarkovProcesses.generate_code_model_type_def(:SquareWaveModel))
@everywhere @eval $(MarkovProcesses.generate_code_model_type_constructor(:SquareWaveModel))
@everywhere @eval $(MarkovProcesses.generate_code_simulation(:SquareWaveModel, :SquareWave_f!, :isabsorbing_SquareWave))

square_wave_oscillator = SquareWaveModel(d, k, dict_var_square, dict_params_square, l_tr_square, 
                                         p_square, x0_square, t0_square,
                                         :SquareWave_f!, :isabsorbing_SquareWave; 
                                         g=g_square_wave, time_bound = 105.0)