text:
mutable struct Trace{Tspl<:AbstractSampler, Tvi<:AbstractVarInfo, Tmodel<:Model}
    model::Tmodel
    spl::Tspl
    vi::Tvi
    ctask::CTask

    function Trace{SampleFromPrior}(model::Model, spl::AbstractSampler, vi::AbstractVarInfo)
        return new{SampleFromPrior,typeof(vi),typeof(model)}(model, SampleFromPrior(), vi)
    end
    function Trace{S}(model::Model, spl::S, vi::AbstractVarInfo) where S<:Sampler
        return new{S,typeof(vi),typeof(model)}(model, spl, vi)
    end
end

function Base.copy(trace::Trace)
    vi = deepcopy(trace.vi)
    res = Trace{typeof(trace.spl)}(trace.model, trace.spl, vi)
    res.ctask = copy(trace.ctask)
    return res
end

# NOTE: this function is called by `forkr`
function Trace(f, m::Model, spl::AbstractSampler, vi::AbstractVarInfo)
    res = Trace{typeof(spl)}(m, spl, deepcopy(vi))
    ctask = CTask() do
        res = f()
        produce(nothing)
        return res
    end
    task = ctask.task
    if task.storage === nothing
        task.storage = IdDict()
    end
    task.storage[:turing_trace] = res # create a backward reference in task_local_storage
    res.ctask = ctask
    return res
end

function Trace(m::Model, spl::AbstractSampler, vi::AbstractVarInfo)
    res = Trace{typeof(spl)}(m, spl, deepcopy(vi))
    reset_num_produce!(res.vi)
    ctask = CTask() do
        res = m(vi, spl)
        produce(nothing)
        return res
    end
    task = ctask.task
    if task.storage === nothing
        task.storage = IdDict()
    end
    task.storage[:turing_trace] = res # create a backward reference in task_local_storage
    res.ctask = ctask
    return res
end

# step to the next observe statement, return log likelihood
Libtask.consume(t::Trace) = (increment_num_produce!(t.vi); consume(t.ctask))

# Task copying version of fork for Trace.
function fork(trace :: Trace, is_ref :: Bool = false)
    newtrace = copy(trace)
    is_ref && set_retained_vns_del_by_spl!(newtrace.vi, newtrace.spl)
    newtrace.ctask.task.storage[:turing_trace] = newtrace
    return newtrace
end

# PG requires keeping all randomness for the reference particle
# Create new task and copy randomness
function forkr(trace::Trace)
    newtrace = Trace(trace.ctask.task.code, trace.model, trace.spl, deepcopy(trace.vi))
    newtrace.spl = trace.spl
    reset_num_produce!(newtrace.vi)
    return newtrace
end

current_trace() = current_task().storage[:turing_trace]

const Particle = Trace

"""
Data structure for particle filters
- effectiveSampleSize(pc :: ParticleContainer)
- normalise!(pc::ParticleContainer)
- consume(pc::ParticleContainer): return incremental likelihood
"""
mutable struct ParticleContainer{T<:Particle}
    "Particles."
    vals::Vector{T}
    "Unnormalized logarithmic weights."
    logWs::Vector{Float64}
end

function ParticleContainer(particles::Vector{<:Particle})
    return ParticleContainer(particles, zeros(length(particles)))
end

Base.collect(pc::ParticleContainer) = pc.vals
Base.length(pc::ParticleContainer) = length(pc.vals)
Base.@propagate_inbounds Base.getindex(pc::ParticleContainer, i::Int) = pc.vals[i]

# registers a new x-particle in the container
function Base.push!(pc::ParticleContainer, p::Particle)
    push!(pc.vals, p)
    push!(pc.logWs, 0.0)
    pc
end

# clones a theta-particle
function Base.copy(pc::ParticleContainer)
    # fork particles
    vals = eltype(pc.vals)[fork(p) for p in pc.vals]

    # copy weights
    logWs = copy(pc.logWs)

    ParticleContainer(vals, logWs)
end

"""
    reset_logweights!(pc::ParticleContainer)

Reset all unnormalized logarithmic weights to zero.
"""
function reset_logweights!(pc::ParticleContainer)
    fill!(pc.logWs, 0.0)
    return pc
end

"""
    increase_logweight!(pc::ParticleContainer, i::Int, x)

Increase the unnormalized logarithmic weight of the `i`th particle with `x`.
"""
function increase_logweight!(pc::ParticleContainer, i, logw)
    pc.logWs[i] += logw
    return pc
end

"""
    getweights(pc::ParticleContainer)

Compute the normalized weights of the particles.
"""
getweights(pc::ParticleContainer) = softmax(pc.logWs)

"""
    getweight(pc::ParticleContainer, i)

Compute the normalized weight of the `i`th particle.
"""
getweight(pc::ParticleContainer, i) = exp(pc.logWs[i] - logZ(pc))

"""
    logZ(pc::ParticleContainer)

Return the logarithm of the normalizing constant of the unnormalized logarithmic weights.
"""
logZ(pc::ParticleContainer) = logsumexp(pc.logWs)

"""
    effectiveSampleSize(pc::ParticleContainer)

Compute the effective sample size ``1 / ∑ wᵢ²``, where ``wᵢ```are the normalized weights.
"""
function effectiveSampleSize(pc::ParticleContainer)
    Ws = getweights(pc)
    return inv(sum(abs2, Ws))
end

"""
    resample_propagate!(pc::ParticleContainer[, randcat = resample_systematic, ref = nothing;
                        weights = getweights(pc)])

Resample and propagate the particles in `pc`.

Function `randcat` is used for sampling ancestor indices from the categorical distribution
of the particle `weights`. For Particle Gibbs sampling, one can provide a reference particle
`ref` that is ensured to survive the resampling step.
"""
function resample_propagate!(
    pc::ParticleContainer,
    randcat = Turing.Inference.resample_systematic,
    ref::Union{Particle, Nothing} = nothing;
    weights = getweights(pc)
)
    # check that weights are not NaN
    @assert !any(isnan, weights)

    # sample ancestor indices
    n = length(pc)
    nresamples = ref === nothing ? n : n - 1
    indx = randcat(weights, nresamples)

    # count number of children for each particle
    num_children = zeros(Int, n)
    @inbounds for i in indx
        num_children[i] += 1
    end

    # fork particles
    particles = collect(pc)
    children = similar(particles)
    j = 0
    @inbounds for i in 1:n
        ni = num_children[i]

        if ni > 0
            # fork first child
            pi = particles[i]
            isref = pi === ref
            p = isref ? fork(pi, isref) : pi
            children[j += 1] = p

            # fork additional children
            for _ in 2:ni
                children[j += 1] = fork(p, isref)
            end
        end
    end

    if ref !== nothing
        # Insert the retained particle. This is based on the replaying trick for efficiency
        # reasons. If we implement PG using task copying, we need to store Nx * T particles!
        @inbounds children[n] = ref
    end

    # replace particles and log weights in the container with new particles and weights
    pc.vals = children
    reset_logweights!(pc)

    pc
end

"""
    reweight!(pc::ParticleContainer)

Check if the final time step is reached, and otherwise reweight the particles by
considering the next observation.
"""
function reweight!(pc::ParticleContainer)
    n = length(pc)

    particles = collect(pc)
    numdone = 0
    for i in 1:n
        p = particles[i]

        # Obtain ``\\log p(yₜ | y₁, …, yₜ₋₁, x₁, …, xₜ, θ₁, …, θₜ)``, or `nothing` if the
        # the execution of the model is finished.
        # Here ``yᵢ`` are observations, ``xᵢ`` variables of the particle filter, and
        # ``θᵢ`` are variables of other samplers.
        score = Libtask.consume(p)

        if score === nothing
            numdone += 1
        else
            # Increase the unnormalized logarithmic weights, accounting for the variables
            # of other samplers.
            increase_logweight!(pc, i, score + getlogp(p.vi))

            # Reset the accumulator of the log probability in the model so that we can
            # accumulate log probabilities of variables of other samplers until the next
            # observation.
            resetlogp!(p.vi)
        end
    end

    # Check if all particles are propagated to the final time point.
    numdone == n && return true

    # The posterior for models with random number of observations is not well-defined.
    if numdone != 0
        error("mis-aligned execution traces: # particles = ", n,
              " # completed trajectories = ", numdone,
              ". Please make sure the number of observations is NOT random.")
    end

    return false
end

"""
    sweep!(pc::ParticleContainer, resampler)

Perform a particle sweep and return an unbiased estimate of the log evidence.

The resampling steps use the given `resampler`.

# Reference

Del Moral, P., Doucet, A., & Jasra, A. (2006). Sequential monte carlo samplers.
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(3), 411-436.
"""
function sweep!(pc::ParticleContainer, resampler)
    # Initial step:

    # Resample and propagate particles.
    resample_propagate!(pc, resampler)

    # Compute the current normalizing constant ``Z₀`` of the unnormalized logarithmic
    # weights.
    # Usually it is equal to the number of particles in the beginning but this
    # implementation covers also the unlikely case of a particle container that is
    # initialized with non-zero logarithmic weights.
    logZ0 = logZ(pc)

    # Reweight the particles by including the first observation ``y₁``.
    isdone = reweight!(pc)

    # Compute the normalizing constant ``Z₁`` after reweighting.
    logZ1 = logZ(pc)

    # Compute the estimate of the log evidence ``\\log p(y₁)``.
    logevidence = logZ1 - logZ0

    # For observations ``y₂, …, yₜ``:
    while !isdone
        # Resample and propagate particles.
        resample_propagate!(pc, resampler)

        # Compute the current normalizing constant ``Z₀`` of the unnormalized logarithmic
        # weights.
        logZ0 = logZ(pc)

        # Reweight the particles by including the next observation ``yₜ``.
        isdone = reweight!(pc)

        # Compute the normalizing constant ``Z₁`` after reweighting.
        logZ1 = logZ(pc)

        # Compute the estimate of the log evidence ``\\log p(y₁, …, yₜ)``.
        logevidence += logZ1 - logZ0
    end

    return logevidence
end

struct ResampleWithESSThreshold{R, T<:Real}
    resampler::R
    threshold::T
end

function ResampleWithESSThreshold(resampler = Turing.Inference.resample_systematic)
    ResampleWithESSThreshold(resampler, 0.5)
end

function resample_propagate!(
    pc::ParticleContainer,
    resampler::ResampleWithESSThreshold,
    ref::Union{Particle,Nothing} = nothing;
    weights = getweights(pc)
)
    # Compute the effective sample size ``1 / ∑ wᵢ²`` with normalized weights ``wᵢ``
    ess = inv(sum(abs2, weights))

    if ess ≤ resampler.threshold * length(pc)
        resample_propagate!(pc, resampler.resampler, ref; weights = weights)
    end

    pc
end


text:
using PyPlot
if isinstalled("seaborn")
    import Seaborn
end

function get_indiv(model::Model, paramset::Int)::Vector{Float64}
    best_generation::Int64 = readdlm(
        joinpath(
            model.path,
            "fitparam",
            "$paramset",
            "generation.dat"
        )
    )[1,1]
    best_indiv::Vector{Float64} = readdlm(
        joinpath(
            model.path,
            "fitparam",
            "$paramset",
            "fit_param$best_generation.dat"
        ),
    )[:,1]
    return best_indiv
end

function load_param(
        model::Model,
        paramset::Int)::Tuple{Array{Float64,1},Array{Float64,1}}
    best_indiv::Vector{Float64} = get_indiv(model, paramset)
    (p, u0) = model.update_param(best_indiv)
    return p, u0
end

function get_executable(model::Model)::Vector{Int}
    n_file::Vector{Int} = []
    fitparam_files::Vector{String} = readdir(
        joinpath(
            model.path,
            "fitparam"
        )
    )
    for file in fitparam_files
        if occursin(r"\d", file)
            push!(n_file, parse(Int64, file))
        end
    end
    empty_folder::Vector{Int} = []
    for (i, nth_param_set) in enumerate(n_file)
        if !isfile(
            joinpath(
                model.path,
                "fitparam",
                "$nth_param_set",
                "generation.dat"
            )
        )
            push!(empty_folder, i)
        end
    end
    for i in sort(empty_folder, rev=true)
        deleteat!(n_file, i)
    end
    return n_file
end

function validate!(model::Model, nth_param_set::Int64)
    (p, u0) = load_param(model, nth_param_set)
    if model.sim.simulate!(p, u0) isa Nothing
        return model, true
    else
        print("Simulation failed. #$nth_param_set\n")
        return model, false
    end
    return model
end

function get_norm_max(
        i::Int, j::Int, obs_name::String, simulations_all::Array{Float64,4})::Float64
    if length(model.sim.normalization) > 0
        norm_max::Float64 = (
            model.sim.normalization[obs_name]["timepoint"] !== nothing ? maximum(
                simulations_all[
                    i,
                    j,
                    model.sim.normalization[obs_name]["timepoint"],
                    [model.cond2idx(c) for c in model.sim.normalization[obs_name]["condition"]]
                ]
            ) : maximum(
                simulations_all[
                    i,
                    j,
                    :,
                    [model.cond2idx(c) for c in model.sim.normalization[obs_name]["condition"]]
                ]
            )
        )
        return norm_max
    else
        return 0.0
    end
end

function plot_timecourse(
        model::Model,
        n_file::Vector{Int},
        viz_type::String,
        show_all::Bool,
        stdev::Bool,
        simulations_all::Array{Float64,4},
        save_format::String)
    if !isdir(
        joinpath(
            model.path,
            "figure",
            "simulation",
            "$viz_type"
        )
    )
        mkpath(
            joinpath(
            model.path,
            "figure",
            "simulation",
            "$viz_type"
        )
        )
    end

    cmap = [
        "#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
        "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf"
    ]
    shape = [
        "o", "v", "^", "<", ">", "8", "s", "p", "*", "h", "H", "D", "d", "P", "X"
    ]

    # rcParams
    rc("figure", figsize=(4, 3))
    rc("font", size=18)
    rc("axes", linewidth=1.5)
    rc("xtick.major", width=1.5)
    rc("ytick.major", width=1.5)
    rc("lines", linewidth=1.8)
    rc("lines", markersize=12)

    for (i, obs_name) in enumerate(model.observables)
        gca().spines["right"].set_visible(false)
        gca().spines["top"].set_visible(false)
        gca().yaxis.set_ticks_position("left")
        gca().xaxis.set_ticks_position("bottom")

        if viz_type != "experiment"
            if show_all
                for j in eachindex(n_file)
                    if length(model.sim.normalization) > 0
                        norm_max = get_norm_max(i, j, obs_name, simulations_all)
                    end
                    for (l, condition) in enumerate(model.sim.conditions)
                        plot(
                            model.sim.t,
                            simulations_all[i,j,:,l] ./ ifelse(
                                length(model.sim.normalization) == 0 || maximum(simulations_all[i,j,:,l]) == 0.0,
                                1.0,
                                norm_max
                            ),
                            color=cmap[l],
                            lw=0.5,alpha=0.35
                        )
                    end
                end
            end
            if viz_type == "average"
                normalized = Array{Float64,4}(
                    undef,
                    length(model.observables),length(n_file),length(model.sim.t),length(model.sim.conditions)
                )
                @inbounds for j in eachindex(n_file)
                    if length(model.sim.normalization) > 0
                        norm_max = get_norm_max(i, j, obs_name, simulations_all)
                    end
                    @simd for l in eachindex(model.sim.conditions)
                        normalized[i,j,:,l] = (
                            simulations_all[i,j,:,l] ./ ifelse(
                                length(model.sim.normalization) == 0 || maximum(simulations_all[i,j,:,l]) == 0.0,
                                1.0,
                                norm_max
                            )
                        )
                    end
                end
                if length(model.sim.normalization) > 0 && model.sim.normalization[obs_name]["timepoint"] === nothing
                    mean_norm_max::Float64 = maximum(
                        vcat(
                            [
                                [
                                    mean(
                                        filter(
                                            !isnan,normalized[i,:,k,l]
                                        )
                                    ) for k in eachindex(model.sim.t)
                                ] for l in eachindex(model.sim.normalization[obs_name]["condition"])
                            ]...
                        )
                    )
                    for j in eachindex(n_file)
                        for k in eachindex(model.sim.t)
                            for l in eachindex(model.sim.conditions)
                                if !isnan(mean_norm_max) && mean_norm_max != 0.0
                                    @inbounds normalized[i,j,k,l] /= mean_norm_max
                                end
                            end
                        end
                    end
                end
                for (l, condition) in enumerate(model.sim.conditions)
                    plot(
                        model.sim.t,[
                            mean(
                                filter(
                                    !isnan,normalized[i,:,k,l]
                                )
                            ) for k in eachindex(model.sim.t)
                        ],
                        color=cmap[l],
                        label=condition
                    )
                end
                if stdev
                    for (l, condition) in enumerate(model.sim.conditions)
                        y_mean = [
                            mean(
                                filter(
                                    !isnan,normalized[i,:,k,l]
                                )
                            ) for k in eachindex(model.sim.t)
                        ]
                        y_std = [
                            std(
                                filter(
                                    !isnan,normalized[i,:,k,l]
                                )
                            ) for k in eachindex(model.sim.t)
                        ]
                        fill_between(
                            model.sim.t,
                            y_mean - y_std, y_mean + y_std,
                            color=cmap[l],
                            lw=0,alpha=0.1
                        )
                    end
                end
            else
                norm_max = length(model.sim.normalization) > 0 ? (
                    model.sim.normalization[obs_name]["timepoint"] !== nothing ? maximum(
                        model.sim.simulations[
                            i,
                            model.sim.normalization[obs_name]["timepoint"],
                            [model.cond2idx(c) for c in model.sim.normalization[obs_name]["condition"]]
                        ]
                    ) : maximum(
                        model.sim.simulations[
                            i,
                            :,
                            [model.cond2idx(c) for c in model.sim.normalization[obs_name]["condition"]]
                        ]
                    )
                ) : 1.0
                for (l, condition) in enumerate(model.sim.conditions)
                    plot(
                        model.sim.t,
                        model.sim.simulations[i,:,l] ./ ifelse(
                            length(model.sim.normalization) == 0 || maximum(model.sim.simulations[i,:,l]) == 0.0,
                            1.0,
                            norm_max
                        ),
                        color=cmap[l],
                        label=condition
                    )
                end
            end
        end

        if isassigned(model.exp.experiments, i)
            exp_t = model.exp.get_timepoint(obs_name)
            if isassigned(model.exp.error_bars, i)
                for (l, condition) in enumerate(model.sim.conditions)
                    if condition in keys(model.exp.experiments[i])
                        exp_data = errorbar(
                            exp_t,
                            model.exp.experiments[i][condition],
                            yerr=model.exp.error_bars[i][condition],
                            lw=1,markerfacecolor="None",
                            color=cmap[l],
                            markeredgecolor=cmap[l],
                            ecolor=cmap[l],
                            fmt=shape[l],
                            capsize=8,
                            clip_on=false
                        )
                        for capline in exp_data[2]
                            capline.set_clip_on(false)
                        end
                        for barlinecol in exp_data[3]
                            barlinecol.set_clip_on(false)
                        end
                    end
                end
            else
                for (l, condition) in enumerate(model.sim.conditions)
                    if condition in keys(model.exp.experiments[i])
                        plot(
                            exp_t,
                            model.exp.experiments[i][condition],
                            shape[l],
                            color=cmap[l],
                            markerfacecolor="None",
                            markeredgecolor=cmap[l],
                            clip_on=false
                        )
                    end
                end
            end
        end
        xlabel("Time")
        ylabel(replace(obs_name, "_" => " "))
        savefig(
            joinpath(
                model.path,
                "figure",
                "simulation",
                "$viz_type",
                "$obs_name." * "$save_format"
            ),
            dpi=save_format == "pdf" ? nothing : 600,
            bbox_inches="tight"
        )
        close()
    end
end

function save_param_range(model::Model, n_file::Vector{Int}, save_format::String)
    search_idx::Tuple{Array{Int64,1},Array{Int64,1}} = model.search_idx()
    popt::Matrix{Float64} = zeros(length(n_file), length(search_idx[1]))
    @inbounds for (i, nth_param_set) in enumerate(n_file)
        best_indiv = get_indiv(model, nth_param_set)
        popt[i, :] = best_indiv[1:length(search_idx[1])]
    end

    # fig, ax = subplots(figsize=(4, length(search_idx[1]) / 5))
    # fig, ax = subplots(figsize=(8, 75))
    # rcParams
    rc("figure", figsize=(4, length(search_idx[1]) / 5))
    rc("font", size=6)

    # sns.despine()
    gca().spines["right"].set_visible(false)
    gca().spines["top"].set_visible(false)
    gca().yaxis.set_ticks_position("left")
    gca().xaxis.set_ticks_position("bottom")

    ax = Seaborn.boxplot(
        data=popt, orient="h", linewidth=0.5, fliersize=1, palette="Set2"
    )

    ax.set_xlabel("Parameter value")
    ax.set_ylabel("")
    ax.set_yticklabels([model.parameters.NAMES[i] for i in search_idx[1]])
    ax.set_xscale("log")

    savefig(
        joinpath(
            model.path,
            "figure",
            "param_range." * "$save_format"
        ),
        dpi=save_format == "pdf" ? nothing : 600,
        bbox_inches="tight"
    )
    close()
end

function run_simulation(
        model::Model;
        viz_type::String,
        show_all::Bool=false,
        stdev::Bool=false,
        save_format::String="pdf")
    if !isdir(
        joinpath(
            model.path,
            "figure"
        )
    )
        mkdir(
            joinpath(
                model.path,
                "figure"
            )
        )
    end

    if !(viz_type in ["best","average","original","experiment"])
        try
            parse(Int64, viz_type)
        catch
            error(
                "Avairable viz_type are: 'best','average','original','experiment','n(=1,2,...)'"
            )
        end
    end

    n_file::Vector{Int} = viz_type in ["original", "experiment"] ? [] : get_executable(model)

    simulaitons_all::Array{Float64,4} = fill(
        NaN,
        (
            length(model.observables),
            length(n_file),
            length(model.sim.t),
            length(model.sim.conditions)
        )
    )
    if viz_type != "experiment"
        if length(n_file) > 0
            if length(n_file) == 1 && viz_type == "average"
                error("viz_type should be best, not $viz_type")
            end
            for (i, nth_param_set) in enumerate(n_file)
                (model, is_successful) = validate!(model, nth_param_set)
                if is_successful
                    for j in eachindex(model.observables)
                        @inbounds simulaitons_all[j,i,:,:] = model.sim.simulations[j,:,:]
                    end
                end
            end
            best_fitness_all::Vector{Float64} = fill(Inf, length(n_file))
            for (i, nth_param_set) in enumerate(n_file)
                if isfile(
                    joinpath(
                        model.path,
                        "fitparam",
                        "$nth_param_set",
                        "best_fitness.dat"
                    )
                )
                    best_fitness_all[i] = readdlm(
                        joinpath(
                            model.path,
                            "fitparam",
                            "$nth_param_set",
                            "best_fitness.dat"
                        )
                    )[1,1]
                end
            end
            best_param_set::Int = n_file[argmin(best_fitness_all)]
            if viz_type == "best"
                model, _ = validate!(model, best_param_set)
            elseif viz_type != "average" && parse(Int64, viz_type) <= length(n_file)
                model, _ = validate!(model, parse(Int64, viz_type))
            elseif viz_type != "average" && parse(Int64, viz_type) > length(n_file)
                error(
                    @sprintf(
                        "n (%d) must be smaller than n_fitparam (%d)",
                        parse(Int64, viz_type), length(n_file)
                    )
                )
            end

            if isinstalled("seaborn") && length(n_file) > 1
                save_param_range(model, n_file, save_format)
            end

        else
            p::Vector{Float64} = param_values()
            u0::Vector{Float64} = initial_values()
            if model.sim.simulate!(p, u0) !== nothing
                error(
                    "Simulation failed."
                )
            end
        end
    end
    plot_timecourse(
        model,n_file,viz_type,show_all,stdev,simulaitons_all, save_format
    )
end

ext: jl


text:
plt = Observable{Any}()
ui = dom"div"(plt)
# Ploting recon
map!(t-> begin
        @manipulate for slice = 1:size(image,3)
            aux = abs.(image) * prod(size(image)[1:2])
            plot_image(angle.(image[:,:,slice]),zmin=-π,zmax=π;darkmode,title="Reconstruction ($slice/$(size(image,3)))")
        end
    end
    , plt, img_obs)
content!(w, "div#content", ui)


text:
println(y)
y = 1

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using TinyNERF, Flux, CUDA, JSON, Images, StatsBase

transforms = JSON.parsefile("local/transforms_train.json")

W, H = 800, 800

const focal = W / (2*tan(transforms["camera_angle_x"]/2))
ṙ = tan(transforms["camera_angle_x"]/2) * 2 / W * 2/√12

poses = []
images = []

for frame in transforms["frames"]
    push!(poses, hcat(frame["transform_matrix"]...)')
    push!(images, float32.(channelview(load("./local/$(frame["file_path"]).png"))[1:3, :, :]))
end

poses = cat(poses..., dims=3)
images = cat(images..., dims=4) # channels × H × W × images

testpose = poses[:, :, 100]
testimg = gpu(images[:, :, :, 100])

poses = poses[:, :, 1:99]
images = images[:, :, :, 1:99]

mlp = TinyNERF.MLP()|>gpu
opt = ADAM(5e-4)

TinyNERF.train!(mlp, opt, images, poses, focal, ṙ, testimg, testpose; n_samples=64, batch_size=256, i_plot=500, n_iters=2000)


text:
"""
Circular Restricted Three-body Model tests.
"""
module PropagationTests

using LinearAlgebra, Unitful, UnitfulAstro, GeneralAstrodynamics, Test

@testset verbose=false "R2BP Propagation" begin

    orbit = KeplerianOrbit(0.4, 10_000, 0, 0, 0, 0, Mars, 0.0) |> CartesianOrbit
    final = propagate(orbit; save_everystep=false)[end]

    @test_broken final.state ≈ orbit.state

    @test all(isapprox.(position_vector(final), position_vector(orbit); atol=1e-8u"km"))
    @test all(isapprox.(velocity_vector(final), velocity_vector(orbit); atol=1e-8u"km/s"))

end

@testset verbose=false "CR3BP Propagation" begin

    @test isperiodic(halo(SunEarth;   L=2, Az=0.0005)...)
    @test isperiodic(halo(SunMars;    L=1, Az=0.0005)...)
    @test isperiodic(halo(EarthMoon;  L=2, Az=0.0005)...)
    @test isperiodic(halo(SunJupiter; L=1, Az=0.0005)...)

    orbit, T = halo(SunJupiter; L=2, Az=0)
    @test isapprox.(
        monodromy(orbit, T), [
            1177.6158450389235 -43.31366732210663 0.0 247.08544503455505 227.41064838834146 0.0;
           -1085.995406851377 40.86470821207657 0.0 -227.41064838735218 -209.97647807144125 0.0;
            0.0 0.0 0.9966422601737439 0.0 0.0 -0.09182654535515093;
            3446.1566018075323 -126.52938312132042 0.0 722.7945482654276 666.0424507143235 0.0;
           -2146.972890519174 79.03452084349573 0.0 -450.8572227441975 -413.95658856183604 0.0;
            0.0 0.0 0.07300944628976104 0.0 0.0 0.99664226017845
        ]; atol=1e-2
    ) |> all

end

end # module


text:
using CorpusLoaders
using Base.Test
using Base.Iterators
using MultiResolutionIterators
using InternedStrings

@testset "basic use" begin
    wk_gen = load(WikiCorpus())
    docs = collect(take(wk_gen, 10));

    @test typeof(docs) ==
        Vector{Document{CorpusLoaders.@NestedVector(String, 4), String}}
#    @test all(isa.(words, AbstractString))
#    @test "a" ∈ words
#    @test "the" ∈ words
    words = collect(flatten_levels(docs, (!lvls)(WikiCorpus, :word)))
    @test typeof(words) == Vector{String}
    @test length(words) >= length(docs)

    vocab = Set(words)
    @test "the" ∈ vocab
    @test "a" ∈ vocab

    #Make sure we're not catching anything we shouldn't
    @test "ENDOFARTICLE" ∉ vocab
    @test "ENDOFARTICLE." ∉ vocab
    @test "<doc" ∉ map(x->x[1:min(end,4)], vocab)

    docs_of_words = full_consolidate(flatten_levels(docs, (!lvls)(WikiCorpus, :doc, :word)))
    @test typeof(docs_of_words) == Vector{Document{Vector{String}, String}}
    @test length(docs_of_words) == length(docs)
    @test sum(length.(docs_of_words)) == length(words)
end


text:
using Random, Distributions, StatsBase, Plots
Random.seed!(0)

mu1, sigma1 = 10, 5
mu2, sigma2 = 40, 12
dist1, dist2 = Normal(mu1,sigma1), Normal(mu2,sigma2)
p = 0.3
mixRv() = (rand() <= p) ? rand(dist1) : rand(dist2)
mixCDF(x) = p*cdf(dist1,x) + (1-p)*cdf(dist2,x)

n = [30, 100]
data1 = [mixRv() for _ in 1:n[1]]
data2 = [mixRv() for _ in 1:n[2]]

empiricalCDF1 = ecdf(data1)
empiricalCDF2 = ecdf(data2)

xGrid = -10:0.1:80
plot(xGrid,empiricalCDF1.(xGrid), c=:blue, label="ECDF with n = $(n[1])")
plot!(xGrid,empiricalCDF2.(xGrid), c=:red, label="ECDF with n = $(n[2])")
plot!(xGrid, mixCDF.(xGrid), c=:black, label="Underlying CDF",
    xlims=(-10,80), ylims=(0,1),
    xlabel="x", ylabel="Probability", legend=:topleft)

max_stars_repo_name: P-Harvey/StatsWithJuliaBook

max_issues_repo_name: P-Harvey/StatsWithJuliaBook

max_forks_repo_name: P-Harvey/StatsWithJuliaBook


text:
@testset "Bus Constructors" begin
    tBus = Bus(nothing)
    tLoadZone = LoadZone(nothing)

    bus = Bus(
        1,
        "test",
        BusTypes.SLACK,
        0.0,
        0.0,
        (min = 0.0, max = 0.0),
        nothing,
        nothing,
        nothing,
    )
    @test PowerSystems.get_bustype(bus) == BusTypes.REF

    @test_throws(
        PowerSystems.DataFormatError,
        Bus(
            1,
            "test",
            BusTypes.ISOLATED,
            0.0,
            0.0,
            (min = 0.0, max = 0.0),
            nothing,
            nothing,
            nothing,
        )
    )
end

@testset "Generation Constructors" begin
    tThreePartCost = ThreePartCost(nothing)
    @test tThreePartCost isa IS.InfrastructureSystemsType
    tTwoPartCost = TwoPartCost(nothing)
    @test tTwoPartCost isa IS.InfrastructureSystemsType
    tThermalGen = ThermalStandard(nothing)
    @test tThermalGen isa PowerSystems.Component
    tHydroDispatch = HydroDispatch(nothing)
    @test tHydroDispatch isa PowerSystems.Component
    tHydroEnergyReservoir = HydroEnergyReservoir(nothing)
    @test tHydroEnergyReservoir isa PowerSystems.Component
    tRenewableFix = RenewableFix(nothing)
    @test tRenewableFix isa PowerSystems.Component
    tRenewableDispatch = RenewableDispatch(nothing)
    @test tRenewableDispatch isa PowerSystems.Component
    tRenewableDispatch = RenewableDispatch(nothing)
    @test tRenewableDispatch isa PowerSystems.Component
end

@testset "Source Constructors" begin
    tSource = Source(nothing)
    @test tSource isa PowerSystems.Component
end

@testset "Storage Constructors" begin
    tStorage = GenericBattery(nothing)
    @test tStorage isa PowerSystems.Component
end

@testset "Load Constructors" begin
    tPowerLoad = PowerLoad(nothing)
    @test tPowerLoad isa PowerSystems.Component
    tPowerLoadPF = PowerLoadPF(nothing)
    @test tPowerLoadPF isa PowerSystems.Component
    tPowerLoad = PowerLoad("init", true, Bus(nothing), nothing, 0.0, 0.0, 100.0, 0.0, 0.0)
    @test tPowerLoad isa PowerSystems.Component
    tPowerLoadPF = PowerLoadPF("init", true, Bus(nothing), nothing, 0.0, 0.0, 1.0, 100.0)
    @test tPowerLoadPF isa PowerSystems.Component
    tLoad = InterruptibleLoad(nothing)
    @test tLoad isa PowerSystems.Component
end

@testset "Branch Constructors" begin
    tLine = Line(nothing)
    @test tLine isa PowerSystems.Component
    tMonitoredLine = MonitoredLine(nothing)
    @test tMonitoredLine isa PowerSystems.Component
    tHVDCLine = HVDCLine(nothing)
    @test tHVDCLine isa PowerSystems.Component
    tVSCDCLine = VSCDCLine(nothing)
    @test tVSCDCLine isa PowerSystems.Component
    tTransformer2W = Transformer2W(nothing)
    @test tTransformer2W isa PowerSystems.Component
    tTapTransformer = TapTransformer(nothing)
    @test tTapTransformer isa PowerSystems.Component
    tPhaseShiftingTransformer = PhaseShiftingTransformer(nothing)
    @test tPhaseShiftingTransformer isa PowerSystems.Component
end

@testset "Service Constructors" begin
    tStaticReserve = StaticReserve{ReserveUp}(nothing)
    @test tStaticReserve isa PowerSystems.Service
    tVariableReserve = VariableReserve{ReserveDown}(nothing)
    @test tVariableReserve isa PowerSystems.Service
end

@testset "TimeSeriesData Constructors" begin
    tg = RenewableFix(nothing)
    data = PowerSystems.TimeSeries.TimeArray(
        [DateTime("01-01-01"), DateTime("01-01-01") + Hour(1)],
        [1.0, 1.0],
    )
    #SingleTimeSeries Tests
    ts = SingleTimeSeries("scalingfactor", Hour(1), DateTime("01-01-01"), 24)
    @test ts isa PowerSystems.TimeSeriesData
    ts = SingleTimeSeries(name = "scalingfactor", data = data)
    @test ts isa PowerSystems.TimeSeriesData
    # TODO 1.0
    #Probabilistic Tests
    #ts = Probabilistic("scalingfactor", Hour(1), DateTime("01-01-01"), [0.5, 0.5], 24)
    #@test ts isa PowerSystems.TimeSeriesData
    #ts = Probabilistic(name = "scalingfactor", percentiles = [1.0], data = data)
    #@test ts isa PowerSystems.TimeSeriesData
    ##Scenario Tests
    #ts = Scenarios("scalingfactor", Hour(1), DateTime("01-01-01"), 2, 24)
    #@test ts isa PowerSystems.TimeSeriesData
    #ts = Scenarios("scalingfactor", data)
    #@test ts isa PowerSystems.TimeSeriesData
end

@testset "Regulation Device" begin
    original_device = ThermalStandard(nothing)
    regulation = RegulationDevice(original_device)
    @test get_rating(regulation) == 0.0
    set_rating!(regulation, 10.0)
    @test get_rating(regulation) == 10.0
    regulation = RegulationDevice(original_device, droop = 0.5)
    @test get_droop(regulation) == 0.5
    @test get_participation_factor(regulation) == (up = 0.0, dn = 0.0)
    @test get_reserve_limit_up(regulation) == 0.0
    @test get_reserve_limit_dn(regulation) == 0.0
    @test get_inertia(regulation) == 0.0
end


text:
using NamedTupleTools
using Test
using OrderedCollections: OrderedDict, LittleDict

namesofvalues  = (:instrument, :madeby)
matchingvalues = ("violin", "Stradivarius")
nt = namedtuple(namesofvalues, matchingvalues)
@test nt == (instrument = "violin", madeby = "Stradivarius")

namesofvalues  = [:instrument, :madeby]
matchingvalues = ["violin", "Stradivarius"]
nt = namedtuple(namesofvalues, matchingvalues)
@test nt == (instrument = "violin", madeby = "Stradivarius")

namesofvalues  = (:instrument, :madeby)
matchingvalues = ["violin", "Stradivarius"]
nt = namedtuple(namesofvalues, matchingvalues)
@test nt == (instrument = "violin", madeby = "Stradivarius")

namesofvalues  = [:instrument, :madeby]
matchingvalues = ("violin", "Stradivarius")
nt = namedtuple(namesofvalues, matchingvalues)
@test nt == (instrument = "violin", madeby = "Stradivarius")

nt = NamedTuple{(:a, :b)}(1.0, "two")

@test fieldnames(nt) == (:a, :b)
@test fieldtypes(nt) == (Float64, String)
@test valtype(nt) == Tuple{Float64, String}
@test fieldtypes(nt) == (Float64, String)

namedtuple(:a, :b) == NamedTuple{(:a, :b),T} where T<:Tuple

ntproto1 = namedtuple(:a, :b, :c, :d)
ntproto2 = namedtuple(:a, :b)

@test ntproto1 == NamedTuple{(:a, :b, :c, :d),T} where T<:Tuple

nt1 = ntproto1(1, 2, 3, 4)
nt2 = ntproto2("one", "two")

ntnt1 = (x = (a = 1, b = 2, c = 3, d = 4), y = (a = "one", b = "two", c = "three"))
ntnt2 = (; y = (a = 1.0, b = "two"))
ntnt3 = (; x = (; a = :one))

proto1 = prototype(nt1)
proto2 = prototype(nt2)

@test isprototype(ntproto1) === true
@test isprototype(nt1) === false
@test isprototype(UnionAll) === false
@test isprototype(proto1) === true
@test isprototype(proto2) === true

@test nt2 === (a = "one", b = "two")
@test ntproto2(nt1) === (a = 1, b = 2)

@test delete(ntproto1, :a) === NamedTuple{(:b, :c, :d),T} where T<:Tuple
@test delete(ntproto1, :a, :c) === NamedTuple{(:b, :d),T} where T<:Tuple
@test delete(ntproto1, (:b, :c)) === NamedTuple{(:a, :d),T} where T<:Tuple

@test delete(nt1, :a) == (b = 2, c = 3, d = 4)
@test delete(nt1, :a, :c) == (b = 2, d = 4)
@test delete(nt1, (:a, :b, :c)) === (d = 4,)

@test merge(nt1, nt2) === (a = "one", b  = "two", c = 3, d = 4)
@test merge(ntnt1, ntnt2) == (
    x = (a = 1, b = 2, c = 3, d = 4),
    y = (a = 1.0, b = "two")
)
@test merge(ntnt1, ntnt2, ntnt3) == (
    x = (a = :one,),
    y = (a = 1.0, b = "two")
)

@test merge_recursive(nt1) == nt1
@test merge_recursive(nt1, nt2) === merge(nt1, nt2)
@test merge_recursive(nt1, nt2, nt1) == merge(nt1, nt2, nt1)
@test merge_recursive(ntnt1) == ntnt1
@test merge_recursive(ntnt1, ntnt2) == (
    x = (a = 1, b = 2, c = 3, d = 4),
    y = (a = 1.0, b = "two", c = "three")
)
@test merge_recursive(ntnt1, ntnt2, ntnt3) == (
    x = (a = :one, b = 2, c = 3, d = 4),
    y = (a = 1.0, b = "two", c = "three")
)

@test select(nt1, :a) == nt1[:a]
@test select(nt1, nt2) == (a=1,b=2)
@test select(nt1, nt2) == select(nt1, keys(nt2))
@test select((a = 1, b = [1, 2]), (:b,)) == (b = [1, 2],)

@test split(nt1, :a)[1] == (a = 1,)
@test split(nt1, :a)[2] == (b = 2, c = 3, d = 4)
@test split(nt1, (:b, :c))[1] == (b = 2, c = 3)
@test split(nt1, (:b, :c))[2] == (a = 1, d = 4)
@test merge(split(nt1, (:a, :b))...) == nt1

struct MyStruct
    tally::Int
    team::String
end

mystruct = MyStruct(5, "hometeam")
mynamedtuple = ntfromstruct(mystruct)

@test mynamedtuple == (tally = 5, team = "hometeam")

nt = (tally=5, team="hometeam")
ntstruct = structfromnt(MyStruct, nt)

@test ntstruct == mystruct

v = [:a => 1, :b => 2]
@test namedtuple(v) == NamedTuple{(:a, :b)}(1, 2)
v = ["a" => 1, "b" => 2]
@test namedtuple(v) == NamedTuple{(:a, :b)}(1, 2)
v = [:b => 2, :a => 1]
@test namedtuple(v) == NamedTuple{(:b, :a)}(2, 1)
v = ["b" => 2, "a" => 1]
@test namedtuple(v) == NamedTuple{(:b, :a)}(2, 1)

for DictType in [Dict, OrderedDict, LittleDict]
    let DT=DictType
        dict1 = DT(:a=>1)
        nt1 = (a = 1,)
        dict2 = DT(:a=>1, :b=>2//11, :c=>"three")
        nt2 = (a = 1, b = 2//11, c = "three")

        nt1_to_dict = convert(DT, nt1)
        @test nt1_to_dict == dict1
        @test nt1_to_dict isa DT
        @test namedtuple(dict1) == nt1

        nt2_to_dict = convert(DT, nt2)
        @test nt2_to_dict == dict2
        @test nt2_to_dict isa DT
        @test namedtuple(dict2) == nt2
    end
end

# Note: checking the types below requires an ordered dictionary
dict1 = LittleDict(:a=>1, :b=>2//11, :c=>"three")
nt1 = namedtuple(dict1)
@test nt1 isa NamedTuple{(:a, :b, :c),Tuple{Int64,Rational{Int64},String}}

dict2 = LittleDict("a"=>1, "b"=>2//11, "c"=>"three")
nt2 = namedtuple(dict2)
@test nt2 isa NamedTuple{(:a, :b, :c),Tuple{Int64,Rational{Int64},String}}

nt = (a = 1, b = 2)
a = 1; b = 2;
nt_ab = @namedtuple(a, b)
@test nt_ab == nt

nt = (a = 1, b = 2, c = 3)
@test @namedtuple(a, b, c = 3) == nt

nt1 = (a = 1, b = 2)
c = 3
nt = (a = 1, b = 2, c = 3)
@test @namedtuple(nt1..., c) == nt


text:
export BatchedArray, BatchedMatrix, BatchedVector, BatchedScalar, Element,
    element_axes,
    element_size

using LinearAlgebra

struct Element{T, N, P} <: AbstractArray{T, N}
    parent::P
    offset::Int
end

Element{T, N}(A::P, offset::Int) where {T, N, P} = Element{T, N, P}(A, offset)

Base.parent(A::Element) = A.parent
Base.axes(A::Element{T, N}) where {T, N} = ntuple(i->axes(parent(A), i), N)
Base.size(A::Element{T, N}) where {T, N} = ntuple(i->size(parent(A), i), N)
Base.strides(A::Element{T, N}) where {T, N} = ntuple(i->stride(parent(A), i), N)

Base.getindex(A::Element, inds::Int...) = _getindex(A, inds)
Base.setindex!(A::Element, v, inds::Int...) = _setindex!(A, v, inds)

_setindex!(A::Element{T, N}, v, inds::NTuple{N, Int}) where {T, N} = setindex!(parent(A), v, _offset(A, inds...) + A.offset + 1)
_getindex(A::Element{T, N}, inds::NTuple{N, Int}) where {T, N} = A.parent[_offset(A, inds...) + A.offset + 1]

function _offset(A::AbstractArray{T, N}, inds::Int...) where {T, N}
    offset = 0
    for (k, idx) in enumerate(inds)
        offset += (idx - 1) * stride(A, k)
    end
    return offset
end

struct BatchedArray{T, N, B, P} <: AbstractArray{Element{T, N}, B}
    parent::P
end

BatchedArray(n::Int, A::AbstractArray) = BatchedArray{eltype(A), n, ndims(A) - n, typeof(A)}(A)
BatchedArray(A::AbstractArray{T, 3}) where T = BatchedArray{T, 2, 1, typeof(A)}(A)
BatchedArray(A::AbstractArray{T, 2}) where T = BatchedArray{T, 1, 1, typeof(A)}(A)
BatchedArray(A::AbstractArray{T, 1}) where T = BatchedArray{T, 0, 1, typeof(A)}(A)

const BatchedMatrix{T, B, P} = BatchedArray{T, 3, B, P}
const BatchedVector{T, B, P} = BatchedArray{T, 2, B, P}
const BatchedScalar{T, B, P} = BatchedArray{T, 1, B, P}

Base.parent(A::BatchedArray) = A.parent
Base.axes(A::BatchedArray{T, N, B}) where {T, N, B} = ntuple(i->axes(parent(A), i+N), B)
Base.size(A::BatchedArray{T, N, B}) where {T, N, B} = ntuple(i->size(parent(A), i+N), B)
Base.strides(A::BatchedArray{T, N, B}) where {T, N, B} = ntuple(i->stride(parent(A), i+N), B)
Base.stride(A::BatchedArray{T, N}, idx::Int) where {T, N} = stride(parent(A), idx+N)

Base.pointer(A::Element) = pointer(A.parent, 1)
Base.getindex(A::BatchedArray{T, N}, inds::Int...) where {T, N} = Element{T, N}(parent(A), _offset(A, inds...))

function Base.setindex!(A::BatchedArray, v, inds...)
    return A[inds...] .= v
end

# creates an Array by default
function Base.similar(::Type{<:BatchedArray{T}}, element_axes::NTuple{N, <:Base.OneTo}, batch_axes::NTuple{B, <:Base.OneTo}) where {T, N, B}
    P = similar(Array{T}, (element_axes..., batch_axes...))
    return BatchedArray{T, N, B, typeof(P)}(P)
end

element_axes(A::BatchedArray{T, N}) where {T, N} = ntuple(i->axes(parent(A), i), N)
element_size(A::BatchedArray{T, N}) where {T, N} = ntuple(i->size(parent(A), i), N)
element_axes(A::AbstractArray) = axes(A)
element_size(A::AbstractArray) = size(A)


text:
const SVNT{N,T} = Union{NTuple{N,T}, <:StaticVector{N,T}}
const AV{T} = AbstractArray{T}

@propagate_inbounds function tuple_vgather(
        v::V, simd, idx = SIMD.Vec(ntuple( i -> N*(i-1), simd ))
    ) where {N,T,V<:AV{<:SVNT{N,T}}}

    vr = reinterpret(T, v)
    vo = ntuple(i -> vr[idx + i], Val(N))
    return vo
end

@propagate_inbounds function tuple_vgather_lazy(
        v::V, simd
    ) where {N,T,V<:AV{<:SVNT{N,T}}}

    n = length(v)
    m = n ÷ simd
    @boundscheck if n - m * simd != 0
        throw(DimensionMismatch("length of input ($n) % simd ($simd) != 0"))
    end
    vr = v |>
        x -> reinterpret(T, v) |>
        x -> reshape(x, (N,simd,m)) |>
        x -> PermutedDimsArray(x, (2,1,3)) |>
        x -> reshape(x, (N*n,)) |>
        x -> reinterpret(SVector{N,SIMD.Vec{simd,T}}, x)
    return vr
end

@propagate_inbounds function tuple_vscatter!(
        vo::VO, vi::VI, idx::SIMD.Vec{simd,I}
    ) where {N,T,simd,VO<:AV{T},VI<:SVNT{N,SIMD.Vec{simd,T}},I<:Integer}

    if @generated
        return quote
            @nexprs( $N, i -> vo[idx + i] = vi[i] )
            return
        end
    else
        for i in 1:N
            vo[idx + i] = vi[i]
        end
    end
end

@propagate_inbounds function tuple_vscatter!(
        vo::VO, vi::VI
    ) where {N,T,simd,VO<:AV{T},VI<:SVNT{N,SIMD.Vec{simd,T}}}

    idx = SIMD.Vec{simd,Int}(ntuple( i -> N*(i-1), Val(simd) ))
    for j in 1:length(vi)
        tuple_vscatter!( vo, vi[j], idx + (j-1)*N*simd )
    end
end

function get_vector_width(::V) where {N,T,simd,V<:AV{<:SVNT{N,SIMD.Vec{simd,T}}}}
    simd
end


text:
include("../src/utils/readfiles.jl")

function get_difference_1(basedir1::String, basedir2::String, tag::String)
    nproc = 441
    nspec = 4480
    me_gll = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec, nproc)
    s362ani_gll = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec, nproc)
    dummy = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)

    for iproc in 0:nproc - 1
        sem_io_read_gll_file_1!(basedir1, iproc, tag, dummy)
        me_gll[:,:,:,:,iproc + 1] = dummy
        sem_io_read_gll_file_1!(basedir2, iproc, tag, dummy)
        s362ani_gll[:,:,:,:,iproc + 1] = dummy
    end
    return s362ani_gll - me_gll, s362ani_gll, me_gll
end

max_stars_repo_name: ziyixi/SeisScripts

max_issues_repo_name: ziyixi/SeisScripts

max_forks_repo_name: ziyixi/SeisScripts


text:
# Non functional system data. This code is just for testing.
# Do not copy paste this code.
@testset "Test Dynamic Source" begin
    sys = System(100)
    bus = Bus(nothing)
    set_bustype!(bus, BusTypes.SLACK)
    add_component!(sys, bus)
    source = Source(nothing)
    set_bus!(source, bus)
    add_component!(sys, source)
    pvs = PeriodicVariableSource(nothing)
    add_component!(sys, pvs, source)
    @test get_components(PeriodicVariableSource, sys) !== nothing
end

@testset "Test Dynamic Source" begin
    sys = System(100)
    bus = Bus(nothing)
    set_bustype!(bus, BusTypes.REF)
    add_component!(sys, bus)
    source = Source(nothing)
    set_bus!(source, bus)
    add_component!(sys, source)
    pvs = PeriodicVariableSource(nothing)
    add_component!(sys, pvs, source)
    @test get_components(PeriodicVariableSource, sys) !== nothing
    sys2, result = validate_serialization(sys)
    @test result
end


text:
function create_system(origin::Origin{T}, joints::Vector{<:JointConstraint}, bodies::Vector{<:Body}, contacts::Vector{<:ContactConstraint}) where T
    adjacency = adjacency_matrix(joints, bodies, contacts)
    dims = length.([joints; bodies; contacts])
    system = System{T}(adjacency, dims, dims)

    for joint in joints
        joint.parent_id == origin.id && (joint.parent_id = 0)
    end

    origin.id = 0

    return system
end

function adjacency_matrix(joints::Vector{<:JointConstraint}, bodies::Vector{<:Body}, contacts::Vector{<:ContactConstraint})
    nodes = [joints; bodies; contacts]
    n = length(nodes)
    A = zeros(Bool, n, n)

    for node1 in nodes
        for node2 in nodes
            if typeof(node1) <: Constraint
                node2.id == node1.child_id && (A[node1.id, node2.id] = 1)
            elseif typeof(node2) <: Constraint
                node1.id == node2.parent_id && (A[node1.id, node2.id] = 1)
            # TODO these entries linking two bodies should be removed,
            # not sure why this is breaking some of the tests
            elseif typeof(node1) <: Body && typeof(node2) <: Body
                for joint in joints
                    if node1.id == joint.parent_id && node2.id == joint.child_id
                        A[node1.id, node2.id] = 1
                    end
                    if node2.id == joint.parent_id && node1.id == joint.child_id
                        A[node2.id, node1.id] = 1
                    end
                end
            end
        end
    end
    A = convert(Matrix{Int64}, A .| A')
    return A
end

get_entry(system, id1, id2) = system.matrix_entries[id1, id2]
get_entry(system, id) = system.vector_entries[id]

function recursivedirectchildren!(system, id::Integer)
    dirs = copy(children(system, id))
    dirslocal = copy(dirs)
    for child_id in dirslocal
        append!(dirs, recursivedirectchildren!(system, child_id))
    end
    return dirs
end