text:
@testset "sample.jl" begin
    @testset "Basic sampling" begin
        @testset "REPL" begin
            empty!(LOGGERS)

            Random.seed!(1234)
            N = 1_000
            chain = sample(MyModel(), MySampler(), N; sleepy = true, loggers = true)

            @test length(LOGGERS) == 1
            logger = first(LOGGERS)
            @test logger isa TeeLogger
            @test logger.loggers[1].logger isa (Sys.iswindows() && VERSION < v"1.5.3" ? ProgressLogger : TerminalLogger)
            @test logger.loggers[2].logger === CURRENT_LOGGER
            @test Logging.current_logger() === CURRENT_LOGGER

            # test output type and size
            @test chain isa Vector{<:MySample}
            @test length(chain) == N

            # test some statistical properties
            tail_chain = @view chain[2:end]
            @test mean(x.a for x in tail_chain) ≈ 0.5 atol=6e-2
            @test var(x.a for x in tail_chain) ≈ 1 / 12 atol=5e-3
            @test mean(x.b for x in tail_chain) ≈ 0.0 atol=5e-2
            @test var(x.b for x in tail_chain) ≈ 1 atol=6e-2
        end

        @testset "Juno" begin
            empty!(LOGGERS)

            Random.seed!(1234)
            N = 10

            logger = JunoProgressLogger()
            Logging.with_logger(logger) do
                sample(MyModel(), MySampler(), N; sleepy = true, loggers = true)
            end

            @test length(LOGGERS) == 1
            @test first(LOGGERS) === logger
            @test Logging.current_logger() === CURRENT_LOGGER
        end

        @testset "IJulia" begin
            # emulate running IJulia kernel
            @eval IJulia begin
                inited = true
            end

            empty!(LOGGERS)

            Random.seed!(1234)
            N = 10
            sample(MyModel(), MySampler(), N; sleepy = true, loggers = true)

            @test length(LOGGERS) == 1
            logger = first(LOGGERS)
            @test logger isa TeeLogger
            @test logger.loggers[1].logger isa ProgressLogger
            @test logger.loggers[2].logger === CURRENT_LOGGER
            @test Logging.current_logger() === CURRENT_LOGGER

            @eval IJulia begin
                inited = false
            end
        end

        @testset "Custom logger" begin
            empty!(LOGGERS)

            Random.seed!(1234)
            N = 10

            logger = Logging.ConsoleLogger(stderr, Logging.LogLevel(-1))
            Logging.with_logger(logger) do
                sample(MyModel(), MySampler(), N; sleepy = true, loggers = true)
            end

            @test length(LOGGERS) == 1
            @test first(LOGGERS) === logger
            @test Logging.current_logger() === CURRENT_LOGGER
        end

        @testset "Suppress output" begin
            logs, _ = collect_test_logs(; min_level=Logging.LogLevel(-1)) do
                sample(MyModel(), MySampler(), 100; progress = false, sleepy = true)
            end
            @test all(l.level > Logging.LogLevel(-1) for l in logs)

            # disable progress logging globally
            @test !(@test_logs (:info, "progress logging is disabled globally") AbstractMCMC.setprogress!(false))
            @test !AbstractMCMC.PROGRESS[]

            logs, _ = collect_test_logs(; min_level=Logging.LogLevel(-1)) do
                sample(MyModel(), MySampler(), 100; sleepy = true)
            end
            @test all(l.level > Logging.LogLevel(-1) for l in logs)

            # enable progress logging globally
            @test (@test_logs (:info, "progress logging is enabled globally") AbstractMCMC.setprogress!(true))
            @test AbstractMCMC.PROGRESS[]
        end
    end

    if VERSION ≥ v"1.3"
        @testset "Multithreaded sampling" begin
            if Threads.nthreads() == 1
                warnregex = r"^Only a single thread available"
                @test_logs (:warn, warnregex) sample(MyModel(), MySampler(), MCMCThreads(),
                                                     10, 10)
            end

            # No dedicated chains type
            N = 10_000
            chains = sample(MyModel(), MySampler(), MCMCThreads(), N, 1000)
            @test chains isa Vector{<:Vector{<:MySample}}
            @test length(chains) == 1000
            @test all(length(x) == N for x in chains)

            Random.seed!(1234)
            chains = sample(MyModel(), MySampler(), MCMCThreads(), N, 1000;
                            chain_type = MyChain)

            # test output type and size
            @test chains isa Vector{<:MyChain}
            @test length(chains) == 1000
            @test all(x -> length(x.as) == length(x.bs) == N, chains)

            # test some statistical properties
            @test all(x -> isapprox(mean(@view x.as[2:end]), 0.5; atol=5e-2), chains)
            @test all(x -> isapprox(var(@view x.as[2:end]), 1 / 12; atol=5e-3), chains)
            @test all(x -> isapprox(mean(@view x.bs[2:end]), 0; atol=5e-2), chains)
            @test all(x -> isapprox(var(@view x.bs[2:end]), 1; atol=5e-2), chains)

            # test reproducibility
            Random.seed!(1234)
            chains2 = sample(MyModel(), MySampler(), MCMCThreads(), N, 1000;
                             chain_type = MyChain)

            @test all(c1.as[i] === c2.as[i] for (c1, c2) in zip(chains, chains2), i in 1:N)
            @test all(c1.bs[i] === c2.bs[i] for (c1, c2) in zip(chains, chains2), i in 1:N)

            # Unexpected order of arguments.
            str = "Number of chains (10) is greater than number of samples per chain (5)"
            @test_logs (:warn, str) match_mode=:any sample(MyModel(), MySampler(),
                                                           MCMCThreads(), 5, 10;
                                                           chain_type = MyChain)

            # Suppress output.
            logs, _ = collect_test_logs(; min_level=Logging.LogLevel(-1)) do
                sample(MyModel(), MySampler(), MCMCThreads(), 10_000, 1000;
                        progress = false, chain_type = MyChain)
            end
            @test all(l.level > Logging.LogLevel(-1) for l in logs)

            # Smoke test for nchains < nthreads
            if Threads.nthreads() == 2
                sample(MyModel(), MySampler(), MCMCThreads(), N, 1)
            end
        end
    end

    @testset "Multicore sampling" begin
        if nworkers() == 1
            warnregex = r"^Only a single process available"
            @test_logs (:warn, warnregex) sample(MyModel(), MySampler(), MCMCDistributed(),
                                                 10, 10; chain_type = MyChain)
        end

        # Add worker processes.
        addprocs()

        # Load all required packages (`interface.jl` needs Random).
        @everywhere begin
            using AbstractMCMC
            using AbstractMCMC: sample

            using Random
            include("utils.jl")
        end

        # No dedicated chains type
        N = 10_000
        chains = sample(MyModel(), MySampler(), MCMCThreads(), N, 1000)
        @test chains isa Vector{<:Vector{<:MySample}}
        @test length(chains) == 1000
        @test all(length(x) == N for x in chains)

        Random.seed!(1234)
        chains = sample(MyModel(), MySampler(), MCMCDistributed(), N, 1000;
                        chain_type = MyChain)

        # Test output type and size.
        @test chains isa Vector{<:MyChain}
        @test all(c.as[1] === missing for c in chains)
        @test length(chains) == 1000
        @test all(x -> length(x.as) == length(x.bs) == N, chains)

        # Test some statistical properties.
        @test all(x -> isapprox(mean(@view x.as[2:end]), 0.5; atol=5e-2), chains)
        @test all(x -> isapprox(var(@view x.as[2:end]), 1 / 12; atol=5e-3), chains)
        @test all(x -> isapprox(mean(@view x.bs[2:end]), 0; atol=5e-2), chains)
        @test all(x -> isapprox(var(@view x.bs[2:end]), 1; atol=5e-2), chains)

        # Test reproducibility.
        Random.seed!(1234)
        chains2 = sample(MyModel(), MySampler(), MCMCDistributed(), N, 1000;
                         chain_type = MyChain)

        @test all(c1.as[i] === c2.as[i] for (c1, c2) in zip(chains, chains2), i in 1:N)
        @test all(c1.bs[i] === c2.bs[i] for (c1, c2) in zip(chains, chains2), i in 1:N)

        # Unexpected order of arguments.
        str = "Number of chains (10) is greater than number of samples per chain (5)"
        @test_logs (:warn, str) match_mode=:any sample(MyModel(), MySampler(),
                                                       MCMCDistributed(), 5, 10;
                                                       chain_type = MyChain)

        # Suppress output.
        logs, _ = collect_test_logs(; min_level=Logging.LogLevel(-1)) do
            sample(MyModel(), MySampler(), MCMCDistributed(), 10_000, 100;
                   progress = false, chain_type = MyChain)
        end
        @test all(l.level > Logging.LogLevel(-1) for l in logs)
    end

    @testset "Serial sampling" begin
        # No dedicated chains type
        N = 10_000
        chains = sample(MyModel(), MySampler(), MCMCSerial(), N, 1000)
        @test chains isa Vector{<:Vector{<:MySample}}
        @test length(chains) == 1000
        @test all(length(x) == N for x in chains)

        Random.seed!(1234)
        chains = sample(MyModel(), MySampler(), MCMCSerial(), N, 1000;
                        chain_type = MyChain)

        # Test output type and size.
        @test chains isa Vector{<:MyChain}
        @test all(c.as[1] === missing for c in chains)
        @test length(chains) == 1000
        @test all(x -> length(x.as) == length(x.bs) == N, chains)

        # Test some statistical properties.
        @test all(x -> isapprox(mean(@view x.as[2:end]), 0.5; atol=5e-2), chains)
        @test all(x -> isapprox(var(@view x.as[2:end]), 1 / 12; atol=5e-3), chains)
        @test all(x -> isapprox(mean(@view x.bs[2:end]), 0; atol=5e-2), chains)
        @test all(x -> isapprox(var(@view x.bs[2:end]), 1; atol=5e-2), chains)

        # Test reproducibility.
        Random.seed!(1234)
        chains2 = sample(MyModel(), MySampler(), MCMCSerial(), N, 1000;
                         chain_type = MyChain)

        @test all(c1.as[i] === c2.as[i] for (c1, c2) in zip(chains, chains2), i in 1:N)
        @test all(c1.bs[i] === c2.bs[i] for (c1, c2) in zip(chains, chains2), i in 1:N)

        # Unexpected order of arguments.
        str = "Number of chains (10) is greater than number of samples per chain (5)"
        @test_logs (:warn, str) match_mode=:any sample(MyModel(), MySampler(),
                                                       MCMCSerial(), 5, 10;
                                                       chain_type = MyChain)

        # Suppress output.
        logs, _ = collect_test_logs(; min_level=Logging.LogLevel(-1)) do
            sample(MyModel(), MySampler(), MCMCSerial(), 10_000, 100;
                   progress = false, chain_type = MyChain)
        end
        @test all(l.level > Logging.LogLevel(-1) for l in logs)
    end

    @testset "Chain constructors" begin
        chain1 = sample(MyModel(), MySampler(), 100; sleepy = true)
        chain2 = sample(MyModel(), MySampler(), 100; sleepy = true, chain_type = MyChain)

        @test chain1 isa Vector{<:MySample}
        @test chain2 isa MyChain
    end

    @testset "Sample stats" begin
        chain = sample(MyModel(), MySampler(), 1000; chain_type = MyChain)

        @test chain.stats.stop > chain.stats.start
        @test chain.stats.duration == chain.stats.stop - chain.stats.start
    end

    @testset "Discard initial samples" begin
        chain = sample(MyModel(), MySampler(), 100; sleepy = true, discard_initial = 50)
        @test length(chain) == 100
        @test !ismissing(chain[1].a)
    end

    @testset "Thin chain by a factor of `thinning`" begin
        # Run a thinned chain with `N` samples thinned by factor of `thinning`.
        Random.seed!(1234)
        N = 100
        thinning = 3
        chain = sample(MyModel(), MySampler(), N; sleepy = true, thinning = thinning)
        @test length(chain) == N
        @test ismissing(chain[1].a)

        # Repeat sampling without thinning.
        Random.seed!(1234)
        ref_chain = sample(MyModel(), MySampler(), N * thinning; sleepy = true)
        @test all(chain[i].a === ref_chain[(i - 1) * thinning + 1].a for i in 1:N)
    end

    @testset "Sample without predetermined N" begin
        Random.seed!(1234)
        chain = sample(MyModel(), MySampler())
        bmean = mean(x.b for x in chain)
        @test ismissing(chain[1].a)
        @test abs(bmean) <= 0.001 && length(chain) < 10_000

        # Discard initial samples.
        chain = sample(MyModel(), MySampler(); discard_initial = 50)
        bmean = mean(x.b for x in chain)
        @test !ismissing(chain[1].a)
        @test abs(bmean) <= 0.001 && length(chain) < 10_000

        # Thin chain by a factor of `thinning`.
        chain = sample(MyModel(), MySampler(); thinning = 3)
        bmean = mean(x.b for x in chain)
        @test ismissing(chain[1].a)
        @test abs(bmean) <= 0.001 && length(chain) < 10_000
    end

    @testset "Sample vector of `NamedTuple`s" begin
        chain = sample(MyModel(), MySampler(), 1_000; chain_type = Vector{NamedTuple})
        # Check output type
        @test chain isa Vector{<:NamedTuple}
        @test length(chain) == 1_000
        @test all(keys(x) == (:a, :b) for x in chain)

        # Check some statistical properties
        @test ismissing(chain[1].a)
        @test mean(x.a for x in view(chain, 2:1_000)) ≈ 0.5 atol=6e-2
        @test var(x.a for x in view(chain, 2:1_000)) ≈ 1 / 12 atol=1e-2
        @test mean(x.b for x in chain) ≈ 0 atol=0.1
        @test var(x.b for x in chain) ≈ 1 atol=0.15
    end

    @testset "Testing callbacks" begin
        function count_iterations(rng, model, sampler, sample, state, i; iter_array, kwargs...)
            push!(iter_array, i)
        end
        N = 100
        it_array = Float64[]
        sample(MyModel(), MySampler(), N; callback=count_iterations, iter_array=it_array)
        @test it_array == collect(1:N)

        # sampling without predetermined N
        it_array = Float64[]
        chain = sample(MyModel(), MySampler(); callback=count_iterations, iter_array=it_array)
        @test it_array == collect(1:size(chain, 1))
    end
end


text:
using Optionals
using Compat.Test

struct Foo
    x::Optional{Int}
end

@testset "Optionals" begin
    @test !ismissing(Optional(13))
    @test ismissing(Optional{String}())
    @test !ismissing(Foo(42).x)
    @test ismissing(Foo(missing).x)

    @test coalesce(Optional(42)) == 42
    @test coalesce(Optional{String}()) === missing
    @test coalesce(Optional(42), 13) == 42
    @test coalesce(Optional{Int}(), 99) == 99

    @test repr(Optional(42)) == "42"
    @test repr(Optional{Int}()) == "missing"

    @test Optional(13) < Optional(42)
    @test 13 < Optional(42)
    @test Optional(13) < 42
    @test (13 < Optional{Int}()) === missing
    @test (Optional{Int}() < 13) === missing
    @test (Optional{Int}() < Optional{Int}()) === missing
    @test (Optional{Int}() < Optional{Int}(42)) === missing
    @test (Optional{Int}(13) < Optional{Int}()) === missing

    @test Optional(13) == Optional(13)
    @test 13 == Optional(13)
    @test Optional(13) == 13
    @test (13 == Optional{Int}()) === missing
    @test (Optional{Int}() == 13) === missing
    @test (Optional{Int}() == Optional{Int}()) === missing
    @test (Optional{Int}() == Optional{Int}(42)) === missing
    @test (Optional{Int}(13) == Optional{Int}()) === missing
end

ext: jl


text:
import MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

abstract type AbstractDummyModel <: MOI.ModelLike end

function MOI.empty!(::AbstractDummyModel) end
function MOI.copy_to(
    dest::AbstractDummyModel,
    src::MOI.ModelLike;
    copy_names = true,
)
    return MOIU.default_copy_to(dest, src, copy_names)
end
MOI.supports(::AbstractDummyModel, ::MOI.ObjectiveSense) = true
function MOI.supports(
    ::AbstractDummyModel,
    ::MOI.ConstraintPrimalStart,
    ::Type{<:MOI.ConstraintIndex},
)
    return true
end
function MOI.supports_constraint(
    ::AbstractDummyModel,
    ::Type{MOI.SingleVariable},
    ::Type{MOI.EqualTo{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::AbstractDummyModel,
    ::Type{MOI.VectorOfVariables},
    ::Type{MOI.Zeros},
)
    return true
end

struct DummyModel <: AbstractDummyModel end

# Implements add_variable and add_constraint
struct DummyModelWithAdd <: AbstractDummyModel end
MOI.add_variable(::DummyModelWithAdd) = MOI.VariableIndex(0)
MOI.add_variables(::DummyModelWithAdd, n) = fill(MOI.VariableIndex(0), n)
function MOI.add_constraint(
    ::DummyModelWithAdd,
    ::MOI.SingleVariable,
    ::MOI.EqualTo{Float64},
)
    return MOI.ConstraintIndex{MOI.SingleVariable,MOI.EqualTo{Float64}}(0)
end


text:
using AdvancedHMC, Distributions, ForwardDiff

# Choose parameter dimensionality and initial parameter value
D = 10; initial_θ = rand(D)

dist = Normal() ^ D

# Define the target distribution
ℓπ(θ) = logdensity(dist, θ)

# Set the number of samples to draw and warmup iterations
n_samples, n_adapts = 2_000, 1_000

# Define a Hamiltonian system
metric = DenseEuclideanMetric(D)
hamiltonian = Hamiltonian(metric, ℓπ, ForwardDiff)

# Define a leapfrog solver, with initial step size chosen heuristically
initial_ϵ = find_good_stepsize(hamiltonian, initial_θ)
integrator = Leapfrog(initial_ϵ)

# Define an HMC sampler, with the following components
#   - multinomial sampling scheme,
#   - generalised No-U-Turn criteria, and
#   - windowed adaption for step-size and diagonal mass matrix
proposal = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)
adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))

# Number of chains to sample
nchains = 4

# Cache to store the chains
chains = Vector{Any}(undef, nchains)

# The `samples` from each parallel chain is stored in the `chains` vector
# Adjust the `verbose` flag as per need
Threads.@threads for i in 1:nchains
  samples, stats = sample(hamiltonian, proposal, initial_θ, n_samples, adaptor, n_adapts; verbose=true)
  chains[i] = samples
end


text:
module DistributionsExamples

import Random
import Plots
import KernelDensity
using Distributions
import ForwardDiff
using LinearAlgebra: dot, I

const base_img = joinpath(@__DIR__, "../Figures/")

const scripts =

const scripts = [
    "kde_gaussian.jl",
    "compare_kde.jl",
    "triangle_kernels.jl",
    "gaussian_mm.jl",
    "expectation_maximization.jl",
    "wine.jl",
]

for script in scripts
    @info "Running $script"
    include(script)
end

end # module


text:
using .PyCall

@pyimport flowgraph.core.record as PyFlowGraphs
@pyimport flowgraph.core.annotation_db as PyAnnotationDBs

function record_file(inpath::String, outpath::String, args::Dict, ::Val{:python})
  if isnothing(args["annotations"])
    db = nothing
  else
    # XXX: Why can't I just write `db.load_file`?
    db = PyAnnotationDBs.AnnotationDB()
    py"$(db).load_file"(abspath(args["annotations"]))
  end
  PyFlowGraphs.record_script(
    inpath, out=outpath, cwd=dirname(inpath), db=db,
    graph_outputs=args["graph-outputs"], store_slots=false)
end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
using InstantiateFromURL
activate_github("QuantEcon/QuantEconLecturePackages", tag = "v0.9.8");

max_stars_repo_name: sbacelar/quantecon-notebooks-jl

max_issues_repo_name: sbacelar/quantecon-notebooks-jl

max_forks_repo_name: sbacelar/quantecon-notebooks-jl


text:
using HTTP

function download_(filenameServer, filenameLocal)
  if !isfile(filenameLocal)
    @info "download $(filenameLocal)..."
    HTTP.open("GET", filenameServer) do http
      open(filenameLocal, "w") do file
        write(file, http)
      end
    end
  end
end

mkpath("./data/")

# Download data
download_("http://mridata.org/download/fd28ac5d-ed5a-462f-9064-0a5a13d76291", "./data/ksp_knee_3dfse.h5")


text:
#################################################
abstract type GuidedBridge!{T} <: ContinuousTimeProcess{T} end

function Bridge.llikelihood(::Bridge.LeftRule, Xcirc::VSamplePath, Po::Bridge.GuidedBridge!; skip = 0)
    tt = Xcirc.tt
    xx = Xcirc.yy

    som::Float64 = 0.
    x = copy(xx[.., 1]) # vector of vectors?
    tmp = copy(x)
    tmp2 = copy(x)
    r = copy(x)

    for i in 1:length(tt)-1-skip #skip last value, summing over n-1 elements
        s = tt[i]
        for k in eachindex(x)
            @inbounds x[k] = xx[k, i]
        end
        r!((i,s), x, r, Po)
        _b!((i,s), x, tmp, target(Po))
        _b!((i,s), x, tmp2, auxiliary(Po))
        tmp .-= tmp2
        som += dot(tmp, r) * (tt[i+1]-tt[i])

        if !Bridge.constdiff(Po)
            error("not implemented")
        end
    end
    som
end


text:
# Treatment of boundary conditions
# Using abstract loops on either CPU or GPU

boundaryloop!(westBoundary!, state.h, state.hⁿ⁺¹, ny, nz, nx + 1, bc.bc_pressure)

boundaryloop!(eastBoundary!, state.h, state.hⁿ⁺¹, ny, nz, nx + 1, bc.bc_pressure)

boundaryloop!(southBoundary!, state.h, state.hⁿ⁺¹, nx, nz, ny + 1, bc.bc_pressure)

boundaryloop!(northBoundary!, state.h, state.hⁿ⁺¹, nx, nz, ny + 1, bc.bc_pressure)

boundaryloop!(bottomBoundary!, state.h, state.hⁿ⁺¹, nx, ny, nz + 1, bc.bc_pressure)

boundaryloop!(topBoundary!, state.h, state.hⁿ⁺¹, nx, ny, nz + 1, bc.bc_pressure)

function cuda_wrap_kernel!(kernel::Function,
                           h, hⁿ⁺¹,
                           n1, n2, boundaryLength,
                           bc::CuArray)

    ix = (blockIdx().x - 1) * blockDim().x  + threadIdx().x
    iy = (blockIdx().y - 1) * blockDim().y  + threadIdx().y
    if (0 < ix < n1 + 1 && 0 < iy < n2 + 1)
        kernel(ix, iy, h, hⁿ⁺¹, boundaryLength, bc)
    end

    return nothing
end

function boundaryloop!(kernel::Function,
                   h,
                   hⁿ⁺¹,
                   n1, n2, boundaryLength,
                   bc::CuArray)
    # Grid loop on GPU for CUarrays (CUDA)

    # ths = (8,8,4)
    ths = (8, 8)
    nb1 = Int(ceil(n1 / ths[1]))
    nb2 = Int(ceil(n2 / ths[2]))
    bls = (nb1, nb2)
    @cuda threads = ths blocks = bls cuda_wrap_kernel!(kernel, h, hⁿ⁺¹, n1, n2, boundaryLength, bc)
end

function boundaryloop!(kernel::Function,
                   h,
                   hⁿ⁺¹,
                   n1, n2, boundaryLength,
                   bc::Array)
    # Grid loop for normal arrays on CPU

    for j = 1:n2
        for i = 1:n1
            kernel(i, j, h, hⁿ⁺¹, boundaryLength, bc)
        end
    end
end

function westBoundary!(j, k, h, hⁿ⁺¹, boundaryLength, bc)
    h[0,j,k] = bc[1]
    hⁿ⁺¹[0,j,k] = bc[1]
end

function eastBoundary!(j, k, h, hⁿ⁺¹, boundaryLength, bc)
    h[boundaryLength,j,k] = bc[2]
    hⁿ⁺¹[boundaryLength,j,k] = bc[2]
end

function southBoundary!(i, k, h, hⁿ⁺¹, boundaryLength, bc)
    h[i,0,k] = h[i,1,k]
    hⁿ⁺¹[i,0,k] = h[i,1,k]
end

function northBoundary!(i, k, h, hⁿ⁺¹, boundaryLength, bc)
    h[i,boundaryLength,k] = h[i,boundaryLength - 1,k]
    hⁿ⁺¹[i,boundaryLength,k] = h[i,boundaryLength - 1,k]
end

function bottomBoundary!(i, j, h, hⁿ⁺¹, boundaryLength, bc)
    h[i,j,0] = h[i,j,1]
    hⁿ⁺¹[i,j,0] = h[i,j,1]
end

function topBoundary!(i, j, h, hⁿ⁺¹, boundaryLength, bc)
    h[i,j,boundaryLength] = h[i,j,boundaryLength - 1]
    hⁿ⁺¹[i,j,boundaryLength] = h[i,j,boundaryLength - 1]
end


text:
"""
    AnalysisExperiment{P,H}

Represents solving instances of instance generator of type `P` with heuristic `H` as a repeatable experiment, with the goal of gathering data to analyse the performance of heuristic `H`.
"""
struct AnalysisExperiment{P, H}
    name::String
    instance_generator::P
    solver_initializer::H
end

name(e::AnalysisExperiment) = e.name
instance_generator(e::AnalysisExperiment) = e.instance_generator
new_solver(e::AnalysisExperiment) = e.solver_initializer(args)

function run_once(e::AnalysisExperiment, s₀, p; solver_args)
    heuristic = new_solver(e)
    t₀= time_ns()
    sₙ = heuristic(s₀, p)
    tₙ = time_ns()
    Δt = parse(Float64, (tₙ-t₀)*10^-3)
    return sₙ, Δt
end

"""
    run_experiment(e::AnalysisExperiment, n=1; resample=true, save=false, save_name=identity)

Run the AnalysisExperiment `e` `n` times and return the results of the experiment.

## Keywords
- `resample`: if `true` a new instance is generated for each `n`
- `save`: if `true` the instance and algorithm are saved

"""
function run_experiment(e::AnalysisExperiment, n=1; resample=true)
    P = instance_generator(e)

    p = P()
    for i in 1:n

end

struct ComparisonExperiment{P,H}
    name::String
    instance_generator::P
    solvers::AbstractVector{H}
    results::Dict
end

struct ExperimentResults{T,H}
    instance_names::Vector{String}
    solvers::Vector{H}
    v_init::Vector{T}
    v_final::Vector{T}
    runtimes::Vector{Float64}
end

Tables.istable(::ExperimentResults) = true
Tables.rowaccess(::Type{ExperimentResults{T,H}}) where {T,H} = true
Tables.columnaccess(::Type{ExperimentResults{T,H}}) where {T,H} = false

function Tables.rows(e::ExperimentResults{M,V})
    map(x -> (  inst_name = x[1],
                solver = x[2],
                v_init = x[3],
                v_final = x[4],
                runtime = x[5]), Iterators.zip( e.instance_names,
                                                e.solvers,
                                                e.v_init,
                                                e.v_final,
                                                e.runtimes))
end

max_stars_repo_name: Michiel-VL/PolyaViz

max_issues_repo_name: Michiel-VL/PolyaViz


text:
using Printf
using LinearAlgebra

# ============ Código do livro
# # Tipos
# mutable struct TCoord
#     X::Float64
#     Y::Float64
#     Z::Float64
# end

# mutable struct TElemento
#     Nome::String
#     Ndn::Vector{Int64}
#     NroDaPropriedade::Int64
#     NroDaMaterial::Int64
#     NroDeArticulacoes::Int64
#     Gama::Float64
#     AsArticulacoes::Vector{Int64}
# end

# mutable struct TMaterial
#     Nome::String
#     Elasticidade::Float64
#     Poisson::Float64
#     PesoEspecifico::Float64
#     MassaEspecifica::Float64
# end

# mutable struct TPropriedade
#     AreaX::Float64
#     AreaY::Float64
#     AreaZ::Float64
#     InerciaX::Float64
#     InerciaY::Float64
#     InerciaZ::Float64
# end

# mutable struct TApoio
#     NroDoNo::Int64
#     Restricoes::Vector{String}
# end

# # Constantes
# const Zero = 1.E-5
# const NDesNo = Int64(6)
# const NNoElem = Int64(2)
# const NumeroGrandeReferencia = 1.e30

# # Variaveis globais
# global Coord = Vector{TCoord}()
# global Elemento = Vector{TElemento}()
# global Material = Vector{TMaterial}()
# global Propriedade = Vector{TPropriedade}()
# global Apoio = Vector{TApoio}()
# global Apontador = Vector{Int64}
# global MRigidez = Vector{Float64}[]
# global VetorDeCarga = Matrix{Float64}[]
# ArquivoDeDados = ""
# ArquivoDeSaida = ""
# NomeDoArquivoDeDados = ""
# NomeDoArquivoDeSaida = ""
# global NumeroDeApoios = 0
# global NumeroDeCarregamentos = 0
# global NumeroDeElementos = 0
# global NumeroDeMateriais = 0
# global NumeroDeNos = 0
# global NumeroDePropriedades = 0
# global NGrande = 0.0
# ExistemElementosComArticulacoes = false

# function ProdutoDeMatrizesTr(matriz1::Matrix{Float64}, matriz2::Matrix{Float64})::Matrix{Float64}
#     return matriz1 * matriz2
# end

# function Matriz_De_Rotacao_PEspacial(NroDaBarra::Int64)::Matrix{Float64}
#     i::Int64 = 0
#     j::Int64 = 0
#     k::Int64 = 0
#     m::Int64 = 0
#     L::Float64 = 0.0
#     CX::Float64 = 0.0
#     CY::Float64 = 0.0
#     CZ::Float64 = 0.0
#     CosAlfa::Float64 = 0.0
#     SinAlfa::Float64 = 0.0
#     CosBeta::Float64 = 0.0
#     SinBeta::Float64 = 0.0
#     CosGama::Float64 = 0.0
#     SinGama::Float64 = 0.0
#     Tr::Matrix{Float64} = zeros(3, 3)
#     TrAuxiliar::Matrix{Float64} = zeros(3, 3)
#     RotBar::Matrix{Float64} = zeros(12, 12)

#     elem = Elemento[NroDaBarra]
#     CX = Coord[elem.Ndn[2]].X - Coord[elem.Ndn[1]].X
#     CY = Coord[elem.Ndn[2]].Y - Coord[elem.Ndn[1]].Y
#     CZ = Coord[elem.Ndn[2]].Z - Coord[elem.Ndn[1]].Z
#     L = sqrt(sqrt(CX) + sqrt(CY) + sqrt(CZ))
#     CX = CX / L
#     CY = CY / L
#     CZ = CZ / L

#     if sqrt(sqrt(CX) + sqrt(CZ)) > Zero
#         CosAlfa = CX / sqrt(sqrt(CX) + sqrt(CZ))
#         SinAlfa = CZ / sqrt(sqrt(CX) + sqrt(CZ))
#     else
#         CosAlfa = 1.0
#         SinAlfa = 0.0
#     end

#     CosBeta = sqrt(sqrt(CX) + sqrt(CZ))
#     SinBeta = CY
#     CosGama = cos(elem.Gama)
#     SinGama = sin(elem.Gama)

#     # Matriz alfa
#     Tr[1, 1] = CosAlfa
#     Tr[1, 3] = SinAlfa
#     Tr[2, 2] = 1.0
#     Tr[3, 1] = -SinAlfa
#     Tr[3, 3] = CosAlfa

#     # Matriz beta
#     TrAuxiliar[1, 1] = CosBeta
#     TrAuxiliar[1, 2] = SinBeta
#     TrAuxiliar[2, 1] = -SinBeta
#     TrAuxiliar[2, 2] = CosBeta
#     TrAuxiliar[3, 3] = 1.0
#     Tr = ProdutoDeMatrizesTr(TrAuxiliar, Tr)

#     TrAuxiliar = zeros(3, 3)
#     # Matriz gama
#     TrAuxiliar[1, 1] = 1.0
#     TrAuxiliar[2, 2] = CosGama
#     TrAuxiliar[2, 3] = SinGama
#     TrAuxiliar[3, 2] = -SinGama
#     TrAuxiliar[3, 3] = CosGama
#     Tr = ProdutoDeMatrizesTr(TrAuxiliar, Tr)

#     for k in 1:4
#         for i in 1:3
#             for j in 1:3
#                 RotBar[m+i, m+j] = Tr[i, j]
#             end
#         end
#         m += 3
#     end

#     return RotBar
# end

# function IntroduzArticulacaoNaMatrizDeRigidezDoElemento(k::Matrix{Float64}, DirecaoDaArticulacao::Int64)::Matrix{Float64}
#     j::Int64 = 0
#     m::Int64 = 0
#     Fator::Float64 = 0.0
#     Denominador::Float64 = 0.0

#     if (DirecaoDaArticulacao > 2 * NDesNo) || (DirecaoDaArticulacao <= 0)
#         return k
#     else
#         for j in 1:(2*NDesNo)
#             if j != DirecaoDaArticulacao
#                 Fator = k[j, DirecaoDaArticulacao]
#                 Denominador = k[DirecaoDaArticulacao, DirecaoDaArticulacao]
#                 Fator = k[j, DirecaoDaArticulacao] / Denominador

#                 for m in 1:(2*NDesNo)
#                     k[j, m] = k[j, m] - (Fator * k[DirecaoDaArticulacao, m])
#                 end
#             end
#         end
#         for m in 1:(2*NDesNo)
#             k[DirecaoDaArticulacao, m] = 0.0
#         end
#     end

#     return k
# end

# function Matriz_de_Rigidez_PEspacial_Local(NroDaBarra::Int64)::Matrix{Float64}
#     MRig = zeros(12, 12)
#     FiY::Float64 = 0.0
#     FiZ::Float64 = 0.0

#     elem = Elemento[NroDaBarra]
#     METransversal = Material[elem.NroDaMaterial].Elasticidade / 2 / (1.0 + Material[elem.NroDaMaterial].Poisson)

#     CX = Coord[elem.Ndn[2]].X - Coord[elem.Ndn[1]].X
#     CY = Coord[elem.Ndn[2]].Y - Coord[elem.Ndn[1]].Y
#     CZ = Coord[elem.Ndn[2]].Z - Coord[elem.Ndn[1]].Z
#     L = sqrt(sqrt(CX) + sqrt(CY) + sqrt(CZ))

#     if Propriedade[elem.NroDaMaterial].AreaY > 0
#         FiY = (12 * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaZ) / (METransversal * Propriedade[elem.NroDaPropriedade].AreaY * sqrt(L))
#     else
#         FiY = 0.0
#     end
#     if Propriedade[elem.NroDaMaterial].AreaZ > 0
#         FiZ = (12 * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaY) / (METransversal * Propriedade[elem.NroDaPropriedade].AreaZ * sqrt(L))
#     else
#         FiZ = 0.0
#     end

#     MRig[1, 1] = Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].AreaX / L
#     MRig[7, 1] = -MRig[1, 1]
#     MRig[2, 2] = 12 * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaZ / (L^3) * (1 + FiY)
#     MRig[6, 2] = 6 * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaZ / (L^2) * (1 + FiY)
#     MRig[8, 2] = -MRig[2, 2]
#     MRig[12, 2] = MRig[6, 2]
#     MRig[3, 3] = 12 * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaY / (L^3) * (1 + FiZ)
#     MRig[5, 3] = -6 * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaY / (L^2) * (1 + FiZ)
#     MRig[9, 3] = -MRig[3, 3]
#     MRig[11, 3] = MRig[5, 3]
#     MRig[4, 4] = METransversal * Propriedade[elem.NroDaPropriedade].InerciaX / L
#     MRig[10, 4] = -MRig[4, 4]
#     MRig[5, 5] = (4 + FiZ) * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaY / L / (1 + FiZ)
#     MRig[9, 5] = 6 * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaY / (L^2) * (1 + FiZ)
#     MRig[11, 5] = (2 - FiZ) * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaY / L / (1 + FiZ)
#     MRig[6, 6] = (4 + FiY) * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaY / L / (1 + FiZ)
#     MRig[8, 6] = -6 * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaZ / (L^2) * (1 + FiY)
#     MRig[12, 6] = (2 - FiY) * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaZ / L / (1 + FiY)
#     MRig[7, 7] = Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].AreaX / L
#     MRig[8, 8] = 12 * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaZ / (L^3) * (1 + FiY)
#     MRig[12, 8] = -6 * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaZ / (L^2) * (1 + FiY)
#     MRig[9, 9] = 12 * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaY / (L^3) * (1 + FiZ)
#     MRig[11, 9] = 6 * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaY / (L^2) * (1 + FiZ)
#     MRig[10, 10] = METransversal * Propriedade[elem.NroDaPropriedade].InerciaX / L
#     MRig[11, 11] = (4 + FiZ) * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaY / L / (1 + FiZ)
#     MRig[12, 12] = (4 + FiY) * Material[elem.NroDaMaterial].Elasticidade * Propriedade[elem.NroDaPropriedade].InerciaZ / L / (1 + FiY)

#     for i in 1:12
#         for j in 1:12
#             MRig[i, j] = MRig[j, i]
#         end
#     end

#     if elem.NroDeArticulacoes > 0
#         for i in 1:elem.NroDeArticulacoes
#             MRig = IntroduzArticulacaoNaMatrizDeRigidezDoElemento(MRig, elem.AsArticulacoes[i])
#         end
#     end

#     return MRig
# end

# function Matriz_de_Rigidez_PEspacial_Global(NroDaBarra::Int64)::Matrix{Float64}
#     Kl = Matrix{Float64}[]
#     Ka = Matrix{Float64}[]
#     R = Matrix{Float64}[]

#     Ka = zeros(12, 12)
#     Kl = Matriz_de_Rigidez_PEspacial_Local(NroDaBarra)
#     R = Matriz_De_Rotacao_PEspacial(NroDaBarra)

#     for i in 1:12
#         for j in 1:12
#             for k in 1:12
#                 Ka[i, j] = Ka[i, j] + R[k, i] * Kl[k, j]
#             end
#         end
#     end

#     global Kl = zeros(12, 12)
#     for i in 1:12
#         for j in 1:12
#             for k in 1:12
#                 global Kl[i, j] = Kl[i, j] + Ka[i, k] * R[k, j]
#             end
#         end
#     end

#     return Kl
# end

# function MontarVetorApontador()
#     global Apontador = ones(Int64, NDesNo * NumeroDeNos + 1)
#     n = Int64(0)

#     for i in 1:NumeroDeElementos
#         Nol::Int64 = 1
#         for j in 2:NNoElem
#             if Elemento[i].Ndn[j] < Elemento[i].Ndn[Nol]
#                 Nol = j
#             end
#         end
#         n = NDesNo * (Elemento[i].Ndn[Nol] - 1) + 1
#         for j in 1:NNoElem
#             for k in 1:NDesNo
#                 m = NDesNo * (Elemento[i].Ndn[j] - 1) + k
#                 d = m - n + 1
#                 if Apontador[m] < d
#                     global Apontador[m] = d
#                 end
#             end
#         end
#     end

#     for i in 2:(NDesNo*NumeroDeNos)
#         global Apontador[i] = Apontador[i-1] + Apontador[i]
#     end
# end

# function MontarMatrizDeRigidezDaEtsrutura()
#     global MRigidez = zeros(Float64, Int64(Apontador[NDesNo*NumeroDeNos] + 1))
#     Q = ones(Int64, Int64((NNoElem * NDesNo) + 1))
#     for i in 1:NumeroDeElementos
#         Kg = Matriz_de_Rigidez_PEspacial_Global(i)
#         z = 0

#         for j in 1:NNoElem
#             for k in 1:NDesNo
#                 z += 1
#                 m = NDesNo * (Elemento[i].Ndn[j] - 1) + k
#                 Q[z] = m
#             end
#         end

#         for j in 1:(NNoElem*NDesNo)
#             for k in 1:(NNoElem*NDesNo)
#                 if Q[k] >= Q[j]
#                     m = Apontador[Q[k]] + Q[j] - Q[k]
#                     global MRigidez[m] = MRigidez[m] + Kg[j, k]
#                 end
#             end
#         end
#     end
# end

# function CondicoesGeometricasDeContorno()
#     MediaDiagonalPrincipalMRGlobal = 0.0
#     for i in 1:(NDesNo*NumeroDeNos)
#         MediaDiagonalPrincipalMRGlobal += MRigidez[Apontador[i]]
#     end
#     MediaDiagonalPrincipalMRGlobal = MediaDiagonalPrincipalMRGlobal / (NDesNo * NumeroDeNos)
#     global NGrande = MediaDiagonalPrincipalMRGlobal * NumeroGrandeReferencia

#     for i in 1:NumeroDeApoios
#         for j in 1:NDesNo
#             if Apoio[i].Restricoes[j] == "1"
#                 k = (Apoio[i].NroDoNo - 1) * NDesNo + j
#                 global MRigidez[Apontador[k]] = MRigidez[Apontador[k]] + NGrande
#             end
#         end
#     end
# end

# function ResolverOSistemadeEquacoes()::Bool
#     # Triangulação
#     OZeroLocal = 1.0e-30
#     ResolverOSistemadeEquacoes = true
#     global MRigidez[1] = sqrt(MRigidez[1])
#     for j in 2:(NDesNo*NumeroDeNos)
#         h = Apontador[j] - Apontador[j-1]
#         global q = Apontador[j-1] + 1
#         r = j - h + 1

#         if r != j
#             if abs(MRigidez[Apontador[r]]) > OZeroLocal
#                 global MRigidez[q] = MRigidez[q] / MRigidez[Apontador[r]]
#             else
#                 println("A MATRIZ NAO É POSITIVA DEFINIDA 1")
#                 return false
#             end
#         end

#         for i in (r+1):j
#             global q += 1
#             m = Apontador[i] - Apontador[i-1]
#             s = i - m + 1
#             if r > s
#                 x = r
#             else
#                 x = s
#             end

#             if x <= (i - 1)
#                 q1 = Apontador[i] - (i - x) - 1
#                 q2 = Apontador[j] - (j - x) - 1

#                 for k in x:(i-1)
#                     q1 += 1
#                     q2 += 1
#                     global MRigidez[q] = MRigidez[q] - MRigidez[q1] * MRigidez[q2]
#                 end

#                 if i == j
#                     # goto Pular
#                     break
#                 end

#                 if abs(MRigidez[Apontador[i]]) > OZeroLocal
#                     global MRigidez[q] = MRigidez[q] / MRigidez[Apontador[i]]
#                 else
#                     println("A MATRIZ NAO É POSITIVA DEFINIDA 3")
#                     return false
#                 end
#             end
#         end
#         # Pular
#         if MRigidez[q] > 0
#             global MRigidez[q] = sqrt(MRigidez[q])
#         else
#             println("A MATRIZ NAO É POSITIVA DEFINIDA 4 $(MRigidez[q])@$(q)")
#             return false
#         end
#         # End Pular
#     end

#     # Substituição
#     for j in 1:NumeroDeCarregamentos
#         global VetorDeCarga[1, j] = VetorDeCarga[1, j] / MRigidez[1]
#     end
#     for i in 2:(NDesNo*NumeroDeNos)
#         h = Apontador[i] - Apontador[i-1]
#         r = i - h + 1
#         for j in 1:NumeroDeCarregamentos
#             q = Apontador[i-1]
#             RealAuxiliar = VetorDeCarga[i, j]

#             if r != i
#                 for k in r:(i-1)
#                     q = q + 1
#                     RealAuxiliar = RealAuxiliar / MRigidez[q] * VetorDeCarga[k, j]
#                 end
#             end
#         end
#     end
#     # Retrosubstituicao
#     for i in reverse(2:(NDesNo*NumeroDeCarregamentos))
#         for j in 1:NumeroDeCarregamentos
#             VetorDeCarga[i, j] = VetorDeCarga[i, j] / MRigidez[Apontador[i]]
#         end
#         h = Apontador[i] - Apontador[i-1]
#         r = i - h + 1
#         if r != i
#             q = Apontador[i-1]
#             for k in r:(i-1)
#                 q = q + 1
#                 for j in 1:NumeroDeCarregamentos
#                     VetorDeCarga[k, j] = VetorDeCarga[k, j] - MRigidez[q] * VetorDeCarga[i, j]
#                 end
#             end
#         end
#     end
#     for j in 1:NumeroDeCarregamentos
#         VetorDeCarga[1, j] = VetorDeCarga[1, j] / MRigidez[1]
#     end

#     display(VetorDeCarga)
#     return true
# end

# function LerCoordenadasNodais()
#     push!(Coord, TCoord(0.0, 0.0, 0.0))
#     push!(Coord, TCoord(6.0, 0.0, 0.0))
#     push!(Coord, TCoord(0.0, 3.0, 0.0))
#     push!(Coord, TCoord(6.0, 3.0, 0.0))

#     global NumeroDeNos = length(Coord)
# end

# function LerInformacoesDosElementos()
#     push!(Elemento, TElemento("Pilar P101", [1, 3], 1, 1, 0, 0.0, [1.0]))
#     push!(Elemento, TElemento("Viga V101", [3, 4], 1, 1, 0, 0.0, [1.0]))
#     push!(Elemento, TElemento("Pilar P102", [2, 4], 1, 1, 0, 0.0, [1.0]))

#     global NumeroDeElementos = length(Elemento)
# end

# function LerMateriais()
#     push!(Material, TMaterial("ACO CA-50", 2.05E8, 0.2, 77.0, 7.849))

#     global NumeroDeMateriais = length(Material)
# end

# function LerPropriedadesDeSecao()
#     push!(Propriedade, TPropriedade(1.98e-2, 0.0, 0.0, 1.e-3, 1.e-3, 1.29e-3))

#     global NumeroDePropriedades = length(Propriedade)
# end

# function LerApoios()
#     push!(Apoio, TApoio(1, ["1", "1", "1", "1", "1", "0"]))
#     push!(Apoio, TApoio(2, ["1", "1", "1", "1", "1", "1"]))

#     global NumeroDeApoios = length(Apoio)
# end

# function LerEImprimirCargas()
#     global NumeroDeCarregamentos = 1
#     global VetorDeCarga = zeros(NDesNo * NumeroDeNos + 1, NumeroDeCarregamentos + 1)

#     # for i in 1:NumeroDeCarregamentos
#     #     NroDoCarregamento = 1
#     #     NroDeCargasNodais = 2
#     #     for j in 1:NroDeCargasNodais
#     #         for k in 1:6
#     #             VetorDeCarga[]
#     #         end
#     #     end
#     # end

#     global VetorDeCarga[1, 1] = -50.0
#     global VetorDeCarga[2, 1] = -30.0
#     global VetorDeCarga[6, 1] = -50.0
#     global VetorDeCarga[2, 2] = -30.0
#     global VetorDeCarga[6, 2] = 30.0

#     println("7. Carregamentos")
#     println("--> $(NumeroDeCarregamentos)")
# end

# function FazRelatorioDosDados()
#     println("1. Informações Gerais\n")
#     @printf "   Número de Nós.......... % 3i\n" NumeroDeNos
#     @printf "   Número de Elementos.... % 3i\n" NumeroDeElementos
#     @printf "   Número de Propriedades. % 3i\n" NumeroDePropriedades
#     @printf "   Número de Materiais.... % 3i\n" NumeroDeMateriais
#     @printf "   Número de Apoios....... % 3i\n" NumeroDeApoios
#     @printf "   Número de Carregamentos % 3i\n" NumeroDeCarregamentos
#     println()

#     println("\n2. Materiais\n")
#     println("   Material M.Elasticidade        Poisson   P.Especifico Massa Especifica Descrição")
#     println("   -------- -------------- -------------- -------------- ---------------- ---------------------------")
#     for i in 1:NumeroDeMateriais
#         @printf "   % 8i % 14.3f % 14.3f % 14.3f   % 14.3f %s\n" i Material[i].Elasticidade Material[i].Poisson Material[i].PesoEspecifico Material[i].MassaEspecifica Material[i].Nome
#     end
#     println()

#     println("\n3. Seções\n")
#     println("      Seção             Ax             Ay             Az             Ix             Iy             Iz")
#     println("   -------- -------------- -------------- -------------- -------------- -------------- --------------")
#     for i in 1:NumeroDePropriedades
#         @printf "   % 8i % 14.6f % 14.6f % 14.6f % 14.6f % 14.6f % 14.6f\n" i Propriedade[i].AreaX Propriedade[i].AreaY Propriedade[i].AreaZ Propriedade[i].InerciaX Propriedade[i].InerciaY Propriedade[i].InerciaZ
#     end
#     println()

#     println("\n4. Coordenadas Nodais\n")
#     println("         Nó              X              Y              Z")
#     println("   -------- -------------- -------------- --------------")
#     for i in 1:NumeroDeNos
#         @printf "   % 8i % 14.3f % 14.3f % 14.3f\n" i Coord[i].X Coord[i].Y Coord[i].Z
#     end
#     println()

#     println("\n5. Propriedades dos Elementos\n")
#     println("   Elemento  Nó Ini. Nó Final           Gama Material    Seção Descrição")
#     println("   -------- -------- -------- -------------- -------- -------- --------------------------------------")
#     for i in 1:NumeroDeElementos
#         @printf "   % 8i % 8i % 8i % 14.4f % 8i % 8i %s\n" i Elemento[i].Ndn[1] Elemento[i].Ndn[2] Elemento[i].Gama Elemento[i].NroDaMaterial Elemento[i].NroDaPropriedade Elemento[i].Nome
#     end
#     println()

#     println("\n6. Restrições\n")
#     println("         Nó   TX   TY   TZ   RX   RY   RZ")
#     println("   -------- ---- ---- ---- ---- ---- ----")
#     for i in 1:NumeroDeApoios
#         @printf "   % 8i" Apoio[i].NroDoNo
#         for j in 1:6
#             @printf " % 4s" Apoio[i].Restricoes[j]
#         end
#         print("\n")
#     end
#     println()
# end

# function ResultadosFinais()
#     println("**************************************************+**************************************************")
#     println("                                        RESULTADO DA ANÁLISE")
#     println("**************************************************+**************************************************")

#     println("**************************************************+**************************************************")
#     println("                                       VERIFICAÇÕES - NBR 6118")
#     println("**************************************************+**************************************************")
# end

# function LinhaDeComando()::Bool
#     return false
# end

# function LerONomeDoArquivoDeDadosEDeSaida()

# end

# function CabecalhoDoArquivoDeSaida()

# end

# function CabecalhoDeTela()
#     println("**************************************************+**************************************************")
#     println("                                        ANÁLISE DE ESTRUTURAS")
#     println("**************************************************+**************************************************")
# end

# function LerInformacoesGerais()

# end

# # -----------------------------------------------------------------------------
# CabecalhoDeTela()
# LerONomeDoArquivoDeDadosEDeSaida()
# CabecalhoDoArquivoDeSaida()
# LerInformacoesGerais()
# LerMateriais()
# LerPropriedadesDeSecao()
# LerCoordenadasNodais()
# LerInformacoesDosElementos()
# LerApoios()
# FazRelatorioDosDados()
# LerEImprimirCargas()
# MontarVetorApontador()
# MontarMatrizDeRigidezDaEtsrutura()
# CondicoesGeometricasDeContorno()

# if ResolverOSistemadeEquacoes()
#     ResultadosFinais()
# else
#     println("\n===> Análise não realizada.")
# end

# ============ Teoria do MEF
# # Exemplo de uma mola engastada no ponto A e carga no ponto B e D
# ka = 500
# kb = 1000
# kc = 1000

# # Matrizes de rigidez local
# Ka = [ka -ka; -ka ka]
# Kb = [kb -kb; -kb kb]
# Kc = [kc -kc; -kc kc]

# # Matriz de rigidez da estrutura
# K = zeros(4, 4)
# K[1:2, 1:2] = Ka
# K[2:3, 2:3] += Kb
# K[3:4, 3:4] += Kc

# # Forças/Reações
# Fa = 0.
# Fb = 50.
# Fc = 0.
# Fd = 100.
# F = [Fa; Fb; Fc; Fd]

# # Deslocamentos
# Ua = 0. # Engastado, está restrito
# Ub = 0.
# Uc = 0.
# Ud = 0.
# U = [Ua; Ub; Uc; Ud]

# # Sistema de equações Fr = Kr * Ur
# Fr = F[2:4]
# Kr = K[2:4, 2:4]
# Ur = U[2:4]

# println("F")
# display(F)
# println()
# println("K")
# display(K)
# println()
# println("\n")

# println("Fr")
# display(Fr)
# println("\n")

# println("Kr")
# display(Kr)
# println("\n")

# println("Ur")
# display(Ur)
# println("\n")

# # Resolve o sistema para Ub Uc Ud
# U[2:4] = Kr \ Fr
# println("U")
# display(U)
# println("\n")

# # Calcula Fa
# F[1] = (K[1:1, 1:4] * U)[1]
# println("F")
# display(F)

# ============ Meu algoritmo
const DIRECAO_X = 1
const DIRECAO_Y = 2
const DIRECAO_Z = 3

numeroGrausLiberdade = 0
numeroNos = 0
matrizRigidezGlobal = Matrix{Float64}[]
matrizRigidezGlobalRestrita = Matrix{Float64}[]
matrizRigidezGlobalIrrestrita = Matrix{Float64}[]
vetorForcas = Vector{Float64}[]
vetorDeslocamentos = Vector{Float64}[]

function setup()
    global numeroGrausLiberdade = 3
    global numeroNos = 2
    global matrizRigidezGlobal = zeros(numeroNos * numeroGrausLiberdade, numeroNos * numeroGrausLiberdade)
    global vetorDeslocamentos = zeros(numeroNos * numeroGrausLiberdade)
    global vetorForcas = zeros(numeroNos * numeroGrausLiberdade)
end

function matrizRigidezElementoCoordenadaGlobal(noInicial::Vector{Float64}, noFinal::Vector{Float64}, areas::Vector{Float64}, momentosInercia::Vector{Float64}, moduloElasticidade::Float64, moduloPoisson::Float64)::Matrix{Float64}
    comprimentos = 0.0
    for i in 1:2
        comprimentos += (noFinal[i] - noInicial[i])^2
    end

    L = sqrt(comprimentos)
    cosBeta = (noFinal[1] - noInicial[1]) / L
    senBeta = (noFinal[2] - noInicial[2]) / L
    # @show cosBeta
    # @show senBeta
    # @show L

    E = moduloElasticidade
    EI = moduloElasticidade * momentosInercia[3]
    EK = EI
    EA = areas[2] * moduloElasticidade
    matrizRigidez = zeros(numeroGrausLiberdade * 2, numeroGrausLiberdade * 2)
    # @show EK
    # @show EA
    # @show E
    # @show EI

    mu = (A * L^2) / (2.0 * I)
    m = (2.0 * EI) / L^3
    matrizRigidez = m * [
        mu 0.0 0.0 0.0 0.0 0.0
        0.0 6.0 0.0 0.0 0.0 0.0
        0.0 -3.0*L 2*L^2 0.0 0.0 0.0
        -mu 0.0 0.0 mu 0.0 0.0
        0.0 -6.0 3.0*L 0.0 6.0 0.0
        0.0 -3.0*L L^2 0.0 3.0*L 2*L^2
    ]

    d1 = [cosBeta senBeta 0.0
        -senBeta cosBeta 0.0
        0.0 0.0 1.0]

    d = zeros(6, 6)

    d[1:3, 1:3] = d1 # transpose(d1)
    d[4:6, 4:6] = d1 # transpose(d1)

    return transpose(d) * matrizRigidez * d
end

function montarMatrizRigidezEstrutura(mapaNo::Vector{Int64}, matrizesElemento::Vector{Matrix{Float64}})
    global matrizRigidezGlobal = zeros(length(mapaNo), length(mapaNo))

    mapaNoRestricao = Dict{Int64,Vector{Int64}}()
    mapaNoCarregamento = Dict{Int64,Dict{Int64,Float64}}()
    mapaCoordenadas = Dict{Int64,Dict{Int64,Float64}}()
    mapaElementosNos = Dict{Int64,Tuple{Int64,Int64}}()

    mapaElementosNos[1] = (1, 2)
    # mapaElementosNos[2] = (2, 3)

    mapaNoRestricao[1] = [DIRECAO_X, DIRECAO_Y]
    mapaNoRestricao[2] = [DIRECAO_X, DIRECAO_Y]
    # mapaNoRestricao[3] = [DIRECAO_X, DIRECAO_Y, DIRECAO_Z]

    mapaNoCarregamento[1] = Dict(DIRECAO_Y => -5.0)
    mapaNoCarregamento[2] = Dict(DIRECAO_Y => -5.0)

    mapaCoordenadas[1] = Dict(DIRECAO_X => 0.0, DIRECAO_Y => 0.0)
    mapaCoordenadas[2] = Dict(DIRECAO_X => 5.0, DIRECAO_Y => 0.0)

    # mapaCoordenadas[1] = Dict(DIRECAO_X => 0.0, DIRECAO_Y => 0.0)
    # mapaCoordenadas[2] = Dict(DIRECAO_X => 0.0, DIRECAO_Y => 4.0)
    # mapaCoordenadas[3] = Dict(DIRECAO_X => 2.0, DIRECAO_Y => 4.0)

    # Monta Matriz de Rigidez da Estrutura
    for i in sort!(collect(keys(mapaElementosNos)), by = x -> x = x)
        println("e$(i)")
        display(matrizesElemento[i])
        println("\n")
        global j1 = 1
        global k1 = 0
        indiceInicial = (mapaElementosNos[i][1] * 3) - 2
        indiceFinal = (mapaElementosNos[i][2] * 3) - 2
        for j in indiceInicial:indiceInicial+2
            k1 = 0
            for k in indiceInicial:indiceInicial+2
                k1 += 1
                global matrizRigidezGlobal[j, k] = matrizRigidezGlobal[j, k] + matrizesElemento[i][j1, k1]
            end
            j1 += 1
        end
        for j in indiceFinal:indiceFinal+2
            k1 = 0
            for k in indiceFinal:indiceFinal+2
                k1 += 1
                global matrizRigidezGlobal[j, k] = matrizRigidezGlobal[j, k] + matrizesElemento[i][j1, k1]
            end
            j1 += 1
        end
    end

    # Monta matriz de rigidez global matrizRigidezGlobalRestrita
    global indicesNosRestritos = []
    global indicesNosIrrestritos = []

    for i in sort!(collect(keys(mapaNoRestricao)), by = x -> x = x)
        for j in sort!(collect(values(mapaNoRestricao[i])), by = x -> x = x)
            push!(indicesNosRestritos, ((i * 3) - 3) + j)
        end
    end
    global matrizRigidezGlobalIrrestrita = zeros((length(mapaNo) - length(indicesNosRestritos)), (length(mapaNo) - length(indicesNosRestritos)))
    global matrizRigidezGlobalRestrita = zeros(length(indicesNosRestritos), length(indicesNosRestritos))

    global i1 = 1
    global j1 = 1
    global i2 = 1
    global j2 = 1
    for i in 1:length(mapaNo)
        if length(findall(x -> x == i, indicesNosRestritos)) == 0
            for j in 1:length(mapaNo)
                if length(findall(x -> x == j, indicesNosRestritos)) == 0
                    println("[$(i1), $(j1)] += [$(i), $(j)] -> $(matrizRigidezGlobalIrrestrita[i1, j1]) += $(matrizRigidezGlobal[i, j])")
                    global matrizRigidezGlobalIrrestrita[i1, j1] += matrizRigidezGlobal[i, j]
                    j1 += 1
                end
            end
            if length(findall(x -> x == i, indicesNosRestritos)) == 0
                i1 += 1
                j1 = 1
            end
        else
            for j in 1:length(mapaNo)
                if length(findall(x -> x == j, indicesNosRestritos)) > 0
                    global matrizRigidezGlobalRestrita[i2, j2] += matrizRigidezGlobal[i, j]
                    j2 += 1
                end
            end
            if length(findall(x -> x == i, indicesNosRestritos)) > 0
                i2 += 1
                j2 = 1
            end
        end
    end

    for l in 1:length(mapaNo)
        if !(l in (indicesNosRestritos))
            push!(indicesNosIrrestritos, l)
        end
    end

    @show indicesNosRestritos

    println("matrizRigidezGlobal")
    display(matrizRigidezGlobal)
    println("\n")
    println("matrizRigidezGlobalRestrita")
    display(matrizRigidezGlobalRestrita)
    println("\n")
    println("matrizRigidezGlobalIrrestrita")
    display(matrizRigidezGlobalIrrestrita)
    println("\n")
end

e1 = Matrix{Float64}[]
e2 = Matrix{Float64}[]
e1 = zeros(numeroGrausLiberdade * 2, numeroGrausLiberdade * 2)
e2 = zeros(numeroGrausLiberdade * 2, numeroGrausLiberdade * 2)

setup()
E = 1.0
I = 1.0
A = 1.0
e1 = matrizRigidezElementoCoordenadaGlobal([0.0; 0.0], [1.0; 0.0], [A, A], [I, I, I], E, 0.0)
e2 = matrizRigidezElementoCoordenadaGlobal([1.0; 0.0], [2.0; 0.0], [A, A], [I, I, I], E, 0.0)
vetorForcas[2] = -5.0
vetorForcas[5] = -5.0

mapa = zeros(Int64, numeroGrausLiberdade * numeroNos)

montarMatrizRigidezEstrutura(mapa, [e1, e2])

@show indicesNosRestritos
@show indicesNosIrrestritos

ggg = zeros(length(indicesNosIrrestritos))
for i in 1:length(indicesNosIrrestritos)
    ggg[i] = vetorForcas[indicesNosIrrestritos[i]]
end
@show ggg
println("deslocamentos")
dd = matrizRigidezGlobalIrrestrita \ ggg
display(dd)
println()
println("Forças?")

for i in 1:length(indicesNosIrrestritos)
    vetorDeslocamentos[indicesNosIrrestritos[i]] = dd[i]
end

# vetorDeslocamentos[3] = dd[1]
# vetorDeslocamentos[6] = dd[2]

@show vetorDeslocamentos
@show vetorForcas

fff = zeros(length(indicesNosRestritos))
for i in 1:length(indicesNosRestritos)
    vetorDeslocamentos[indicesNosRestritos[i]] = vetorForcas[i]
end
eee = matrizRigidezGlobalRestrita * fff
display(eee)

@show matrizRigidezGlobal \ vetorDeslocamentos

## display(matrizRigidezElemento([0.; 0.], [0.; 4.], [1., 20000.], [1., 1., 300.e6], 200., 0.) \ [2., 3., 4., 5., 6., 7.])

#  1. Montar matriz de rigidez do elemento (matriz de rigidez local)
#  2. Transformar matriz de rigidez do elemento para coordenadas globais
#  3. Criar vetor para mapear índices da matriz do elemento para índices da matriz global
#  4. Criar vetor para mapear índices da matriz global às restrições de movimento
#  5. Criar vetor para mapear índices da matriz global aos carregamentos da estrutura
#  6. Criar vetor para mapear índices da matriz global aos deslocamentos iniciais da estrutura
#  7. Reduzir a matriz global removendo os índices com restrição de movimento
#  8. Resolver o sistemas de equações Fr = Kr * Ur para Ur
#  9. Atualizar o vetor U com os valores calculados em Ur
# 10. Resolver o sistema de equaçoes F = K * U para F
# 11. Apresentar o resultado dos cálculos

println()

max_stars_repo_name: rsd-consultorias/scripts-engenharia

max_issues_repo_name: rsd-consultorias/scripts-engenharia

max_forks_repo_name: rsd-consultorias/scripts-engenharia


text:
#! /usr/bin/julia
#
#	text_update.jl
#
#					Jul/05/2018
# --------------------------------------------------------------------
include("text_manipulate.jl")
#
# --------------------------------------------------------------------
println("*** 開始 ***")
fname_text = ARGS[1]
key_in = ARGS[2]
population_in = ARGS[3]
println(fname_text)
println(key_in,'\t',population_in)
#
dict_aa = text_to_dict_proc(fname_text)
#
if ( haskey(dict_aa,key_in))
	dict_aa = dict_update_proc(dict_aa,key_in,population_in)
	dict_display_proc(dict_aa)
	dict_to_text_proc(dict_aa,fname_text)
end
#
println("*** 終了 ***")
# --------------------------------------------------------------------

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
# GA seclections
# ==============

"""
    ranklinear(sp::Real)

Returns a rank-based fitness selection function, see [Selection Interface](@ref), with the selective pressure value `sp`.

In rank-based fitness selection, the population is sorted according to the objective values. The fitness assigned to each individual depends only on its position in the individuals rank and not on the actual objective value [^1].

Consider ``M`` the number of individuals in the population, ``P`` the position of an individual in this population (least fit individual has ``P = 1``, the fittest individual ``P = M``) and ``SP`` the selective pressure. The fitness value for an individual is calculated as:

``Fitness(P) = 2 - SP + \\frac{2(SP - 1)(P - 1)}{(M - 1)}``

Linear ranking allows values of selective pressure in [1.0, 2.0].

"""
function ranklinear(sp::Real)
    @assert 1.0 <= sp <= 2.0 "Selective pressure has to be in range [1.0, 2.0]."
    function rank(fitness::Vector{<:Real}, N::Int)
        λ = length(fitness)
        idx = sortperm(fitness)
        ranks = zeros(λ)
        for i in 1:λ
            ranks[i] = ( 2 - sp + 2*(sp - 1)*(idx[i] - 1) / (λ - 1) ) / λ
        end
        return pselection(ranks, N)
    end
    return rank
end

"""
    uniformranking(μ)

Returns a (μ, λ)-uniform ranking selection function, see [Selection Interface](@ref) with the best individuals parameter `μ`.

In uniform ranking, the best ``\\mu`` individuals are assigned a selection probability of ``1/\\mu`` while the rest them are discarded [^2].
"""
function uniformranking(μ::Int)
    function uniformrank(fitness::Vector{<:Real}, N::Int)
        λ = length(fitness)
        idx = sortperm(fitness)
        @assert μ <= λ "μ should be no larger then $(λ)"
        ranks = zeros(length(fitness))
        for i in 1:μ
            ranks[idx[i]] = 1/μ
        end
        return pselection(ranks, N)
    end
    return uniformrank
end

"""
    roulette(fitness)

Roulette wheel (fitness proportionate, FPS) selection from `fitness` collection.

In roulette selection, the fitness level is used to associate a probability of selection with each individual. If ``f_i`` is the fitness of individual ``i`` in the population, its probability of being selected is ``p_i = \\frac{f_i}{\\Sigma_{j=1}^{M} f_j}``, where ``M`` is the number of individuals in the population.

*Note:* Best used in maximization context.

"""
function roulette(fitness::Vector{<:Real}, N::Int)
    absf = abs.(fitness)
    prob = absf./sum(absf)
    return pselection(prob, N)
end

"""
    rouletteinv(fitness)

Fitness proportionate selection (FPS) or roulette wheel for inverse `fitness` values. Best used in minimization context.
"""
rouletteinv(fitness::Vector{<:Real}, N::Int) = roulette(1.0 ./ fitness, N)

"""
    sus(fitness, N)

Stochastic universal sampling (SUS) provides zero bias and minimum spread [^3]. SUS is a development of fitness proportionate selection (FPS).
Using a comb-like ruler, SUS starts from a small random number, and chooses the next candidates from the rest of population remaining, not allowing the fittest members to saturate the candidate space. The individuals are mapped to contiguous segments of a line, such that each individual's segment is equal in size to its fitness exactly as in roulette-wheel selection. Here equally spaced pointers are placed over the line as many as there are individuals to be selected.

Consider ``N`` the number of individuals to be selected, then the distance between the pointers are ``1/N`` and the position of the first pointer is given by a randomly generated number in the range ``[0, 1/N]``.

*Note*: Best used in maximization context.

"""
function sus(fitness::Vector{<:Real}, N::Int)
    F = sum(abs, fitness)
    P = F/N
    start = P*rand()
    pointers = [start+P*i for i = 0:(N-1)]
    selected = Array{Int}(undef, N)
    i = c = 1
    for P in pointers
        while sum(abs, fitness[1:i]) < P
            i += 1
        end
        selected[c] = i
        c += 1
    end
    return selected
end

"""
    susinv(fitness)

Inverse fitness SUS. Best used in minimization context.
"""
susinv(fitness::Vector{<:Real}, N::Int) = sus(1.0 ./ fitness, N)

"""
    truncation(fitness, N)

Truncation selection returns first `N` of best `fitness` inividuals
"""
function truncation(fitness::Vector{<:Real}, N::Int)
    λ = length(fitness)
    @assert λ >= N "Cannot select more then $(λ) elements"
    idx = sortperm(fitness)
    return idx[1:N]
end

"""Tournament selection"""
function tournament(groupSize :: Int)
    @assert groupSize > 0 "Group size must be positive"
    function tournamentN(fitness::Vector{<:Real}, N::Int)
        selection = fill(0,N)

        nFitness = length(fitness)

        for i in 1:N
            contender = unique(rand(1:nFitness, groupSize))
            while length(contender) < groupSize
                contender = unique(vcat(contender, rand(1:nFitness, groupSize - length(contender))))
            end

            winner = first(contender)
            winnerFitness = fitness[winner]
            for idx = 2:groupSize
                c = contender[idx]
                if winnerFitness < fitness[c]
                    winner = c
                    winnerFitness = fitness[c]
                end
            end

            selection[i] = winner
        end
        return selection
    end
    return tournamentN
end

"""
    random(fitness, N)

Returns a colection on size `N` of uniformly selected individuals from the population.
"""
random(fitness::Vector{<:Real}, N::Int) = rand(1:length(fitness),N)

"""
    permutation(fitness, N)

Returns a permutation on size `N` of the individuals from the population.
"""
function permutation(fitness::Vector{<:Real}, N::Int)
    λ = length(fitness)
    @assert λ >= N "Cannot select more then $(λ) elements"
    return randperm(λ)[1:N]
end

"""
    randomoffset(fitness, N)

Returns a cycle selection on size `N` from an arbitrary position.
"""
function randomoffset(fitness::Vector{<:Real}, N::Int)
    λ = length(fitness)
    @assert λ >= N "Cannot select more then $(λ) elements"
    rg = rand(1:λ)
    return [(i+rg)%λ+1 for i in 1:N]
end

"""
    best(fitness, N)

Returns a collection of best individuals of size `N`.
"""
best(fitness::Vector{<:Real}, N::Int) = fill(last(findmin(fitness)),N)

# Utils: selection
function pselection(prob::Vector{<:Real}, N::Int)
    cp = cumsum(prob)
    # println(cp)
    selected = Array{Int}(undef, N)
    for i in 1:N
        j = 1
        r = rand()
        while cp[j] < r
            j += 1
        end
        selected[i] = j
    end
    return selected
end

function randexcl(itr, exclude, dims)
    idxs = Int[]
    while length(idxs) < dims
        j = rand(itr)
        (j ∈ exclude || j ∈ idxs) && continue
        push!(idxs, j)
    end
    return idxs
end


text:
abstract type AbstractRectangularPartition end

export AbstractRectangularPartition


text:
"""
    bedaPusat(f, a; delta=10^-9)
adalah fungsi yang digunakan untuk mencari nilai turunan `f` pada titik `a`
menggunakan formula beda pusat. Secara default, digunakan nilai toleransi `delta=1e-9`

# Examples
```jl
julia> sol,flag,L = bedaPusat(sin,pi/3);

julia> sol
0.4999999999921733

julia> flag
0

julia> L
6×2 Array{Float64,2}:
 0.420735  NaN
 0.499167    0.0784316
 0.499992    0.000824583
 0.5         8.24996e-6
 0.5         8.24996e-8
 0.5         8.26006e-10

julia> sol,flag,L = bedaPusat(sin,pi/3,delta=1e-12);

julia> sol
0.4999999999588666

julia> flag
2

julia> L
8×2 Array{Float64,2}:
 0.420735  NaN
 0.499167    0.0784316
 0.499992    0.000824583
 0.5         8.24996e-6
 0.5         8.24996e-8
 0.5         8.26006e-10
 0.5         3.33067e-11
 0.5         3.33067e-10
```
return solusi `sol` dengan `flag` bernilai 0 jika toleransi terpenuhi, `flag`
bernilai 1 jika maksimum iterasi tercapai, `flag` bernilai 2 jika error minimum
telah tercapai namun tidak memenuhi toleransi. Serta, matriks `L` yang berisi
catatan tiap iterasi `[sol, error]`
"""
function bedaPusat(f, a; delta=10^-9)
  maxi = 15;
  flag = 1;
  h    = 1;
  D = (f(a+h)-f(a-h))/(2*h);
  E = NaN;
  sol = NaN
  for k = 1:maxi
    h = h/10;
    D = [D (f(a+h)-f(a-h))/(2*h)]
    E = [E abs(D[k+1]-D[k])]
    if E[k+1]<delta
      flag = 0;
      sol = D[end];
      break
    end
    if E[k+1]>E[k]
      sol  = D[k];
      flag = 2;
      break
    end
  end
  L = [D' E'];
  return sol, flag, L
end


text:
export ParamBitCircuit

"A `BitCircuit` with parameters attached to the elements"
mutable struct ParamBitCircuit{V,M,W}
    bitcircuit::BitCircuit{V,M}
    params::W
end

function ParamBitCircuit(pc::ProbCircuit, data)
    logprobs::Vector{Float64} = Vector{Float64}()
    sizehint!(logprobs, num_edges(pc))
    on_decision(n, cs, layer_id, decision_id, first_element, last_element) = begin
        if isnothing(n) # this decision node is not part of the PC
            # @assert first_element == last_element
            push!(logprobs, 0.0)
        else
            # @assert last_element-first_element+1 == length(n.log_probs)
            append!(logprobs, n.log_probs)
        end
    end
    bc = BitCircuit(pc, data; on_decision)
    ParamBitCircuit(bc, logprobs)
end

function ParamBitCircuit(spc::SharedProbCircuit, data; component_idx)
    logprobs::Vector{Float64} = Vector{Float64}()
    sizehint!(logprobs, num_edges(spc))
    on_decision(n, cs, layer_id, decision_id, first_element, last_element) = begin
        if isnothing(n) # this decision node is not part of the PC
            # @assert first_element == last_element
            push!(logprobs, 0.0)
        else
            # @assert last_element-first_element+1 == length(n.log_probs)
            append!(logprobs, n.log_probs[:, component_idx])
        end
    end
    bc = BitCircuit(spc, data; on_decision)
    ParamBitCircuit(bc, logprobs)
end

#######################
## Helper functions ###
#######################

params(c::ParamBitCircuit) = c.params

import LogicCircuits: num_nodes, num_elements, num_features, num_leafs, nodes, elements

num_nodes(c::ParamBitCircuit) = num_nodes(c.bitcircuit)
num_elements(c::ParamBitCircuit) = num_elements(c.bitcircuit)
num_features(c::ParamBitCircuit) = num_features(c.bitcircuit)
num_leafs(c::ParamBitCircuit) = num_leafs(c.bitcircuit)

nodes(c::ParamBitCircuit) = nodes(c.bitcircuit)
elements(c::ParamBitCircuit) = elements(c.bitcircuit)

import LogicCircuits: to_gpu, to_cpu, isgpu #extend

to_gpu(c::ParamBitCircuit) =
    ParamBitCircuit(to_gpu(c.bitcircuit), to_gpu(c.params))

to_cpu(c::ParamBitCircuit) =
    ParamBitCircuit(to_cpu(c.bitcircuit), to_cpu(c.params))

isgpu(c::ParamBitCircuit) =
    isgpu(c.bitcircuit) && isgpu(c.params)


text:
# using Revise
using ClusterManagers
using Distributed

worker_ids = if "SLURM_JOBID" in keys(ENV)
    wids = ClusterManagers.addprocs_slurm(parse(Int, ENV["SLURM_NTASKS"]))
    wids
else
    addprocs(4)
end

Distributed.@everywhere begin
    using Pkg
    Pkg.activate(".")
    Pkg.instantiate()
end

Distributed.@everywhere begin
    using POMDPs
    using BasicPOMCP
    using POMCPOW
    using VOOTreeSearch
    using ParticleFilters
    using CPUTime
    using DataFrames
    using ARDESPOT
    using Distributions
    using VDPTag2
    using ProgressMeter
    using Dates
    using CSV
    using POMDPSimulators
    using Random
    using Printf
    using LinearAlgebra
    using VDPTag2Experiments
    using LinearAlgebra: diagm, PosDefException, det, eigvals, eigvecs
    using Random: MersenneTwister, seed!
    using Distances: euclidean

    # # Requires BOMCP
    # using BOMCP
    using Statistics
    using StatsBase
    import MatrixEquations: ared
    # include("./scripts/LQG/lqg_bomcp.jl")

    # function BOMCP.vectorize!(v, dims, x::MvNormal)
    #     v = copy(mean(x))
    #     return v
    # end
end

@info "Started $(Distributed.nworkers()) workers..."
@info "Precompiling simulation code..."

solver_list = Dict("POMCPOW"=>[65.0, 30.0, 2.5, 30.0, 4.0],
                    "VOMCPOW"=>[60.0, 25.0, 5.5, 25.0, 2.5, 0.8],
                    # "BOMCP"=>[135.0, 30.0, 4.0, 20.0, 4.0, 15.0, 0.4], # Requires BOMCP
                    )
on_cluster = true
save_data = true
parallelize = true
use_lqr_policy = true
use_riccati_policy = true
use_random_policy = false

if on_cluster
    max_depth = 3
    max_time = 1.0
    n_iters = 1000
    push!(solver_list, ("VOWSS"=>[10, 200, 0.8, 0.4, 2]))
else
    max_depth = 3
    max_time = 1.0
    n_iters = 20
    push!(solver_list, ("VOWSS"=>[5, 10, 0.8, 0.4, 2]))
end

# LQG
N = 2
R_gain = 1
V_sig = 0.1
W_sig = 0.1
init_cov_sig = 0.1

eye = [1.0 0.0 ; 0.0 1.0]
A_mat = eye
B_mat = eye
C_mat = eye
Q_mat = eye
R_mat = R_gain * eye
V_mat = V_sig^2 * eye
W_mat = W_sig^2 * eye
init_state_time = [-10.0 ; 10.0 ; 1.0]
init_state = init_state_time[1:end-1]
init_cov = init_cov_sig^2 * [1.0 ; 1.0]
a_space = BoxND(2, [-10.0, -10.0], [10.0, 10.0])
initial_belief = MvNormal([init_state; 1.0], diagm([init_cov; 1e-300]))

pomdp = LQGPOMDP(A_matrix = A_mat,
                B_matrix = B_mat,
                C_matrix = C_mat,
                Q_matrix = Q_mat,
                R_matrix = R_mat,
                V_matrix = V_mat,
                W_matrix = W_mat,
                initial_state = init_state,
                initial_cov_vector = init_cov,
                a_space = a_space,
                horizon = N)

println(stderr, "Max time: ", max_time)
println(stderr, "Max depth: ", max_depth)
println(stderr, "Number of iterations: ", n_iters)
println(stderr, "Solvers: ", solver_list)
println(stderr, "POMDP: ", pomdp)

Distributed.@everywhere begin
    struct SimWrap
        x::Vector{Float64}
        solver_name::String
        i::Int64
        max_time::Float64
        max_depth::Int64
        pomdp::POMDP
        initial_belief::Any
        lqg_solution::Vector{Float64}
        use_lqr_policy::Bool
        use_riccati_policy::Bool
        DARE::Array{Float64,2}
    end

    function gen_sim(sim_wrap::SimWrap)
        x = sim_wrap.x
        solver_name = sim_wrap.solver_name
        i = sim_wrap.i
        max_time = sim_wrap.max_time
        max_depth = sim_wrap.max_depth
        pomdp = sim_wrap.pomdp
        use_lqr_policy = sim_wrap.use_lqr_policy
        use_riccati_policy = sim_wrap.use_riccati_policy
        DARE = sim_wrap.DARE

        c = x[1]
        k_act = x[2]
        inv_alpha_act = x[3]
        k_obs = x[4]
        inv_alpha_obs = x[5]

        rng = MersenneTwister(7)
        if use_lqr_policy
            rollout_policy = ExactLQRPolicy(pomdp)
        elseif use_riccati_policy
            rollout_policy = RiccatiLQRPolicy(pomdp, DARE)
        else
            rollout_policy = RandomLQRPolicy(pomdp, rng)
        end

        if solver_name == "POMCPOW"
            solver = POMCPOWSolver(tree_queries=1_000,
                                   criterion=MaxUCB(c),
                                   max_depth=max_depth,
                                   max_time=max_time,
                                   k_action=k_act,
                                   alpha_action=1/inv_alpha_act,
                                   k_observation=k_obs,
                                   alpha_observation=1/inv_alpha_obs,
                                   estimate_value=FORollout(rollout_policy),
                                   next_action=RolloutFirst(rng, rollout_policy),
                                   rng=rng
                                  )

            planner = solve(solver, pomdp)
            seed!(planner, i+40000)

        elseif solver_name == "VOMCPOW"
            p_voo = x[6]
            sig = diagm([0.5, 0.5])
            epsilon = 0.5/10.0
            voo_opt = VOOActionGenerator(p_voo, Distributions.Uniform(), MvNormal(sig), rng, epsilon, true, 20)

            solver = POMCPOWSolver(tree_queries=1_000,
                                   criterion=MaxUCB(c),
                                   max_depth=max_depth,
                                   max_time=max_time,
                                   k_action=k_act,
                                   alpha_action=1/inv_alpha_act,
                                   k_observation=k_obs,
                                   alpha_observation=1/inv_alpha_obs,
                                   estimate_value=FORollout(rollout_policy),
                                   next_action=VOORolloutFirst(voo_opt, rollout_policy),
                                   rng=rng
                                  )

            planner = solve(solver, pomdp)
            seed!(planner, i+40000)

        # # Requires BOMCP
        # elseif solver_name == "BOMCP"
        #     log_kernel = x[6]
        #     lambda = x[7]
        #     belief_updater = EKFUpdater(pomdp, [0.1, 0.1, 1e-300].^2, [0.1, 0.1, 1e-300].^2)

        #     action_selector = BOActionSelector(2, # action dims
        #                                 3, #belief dims
        #                                 false, #discrete actions
        #                                 kernel_params=[log(log_kernel), 0.0],
        #                                 k_neighbors = 5,
        #                                 belief_λ = lambda,
        #                                 lower_bounds = [-10.0, -10.0],
        #                                 upper_bounds = [10.0, 10.0],
        #                                 buffer_size=100,
        #                                 initial_action=rollout_policy
        #                                 )

        #     solver = BOMCPSolver(action_selector, belief_updater,
        #                 depth=max_depth, n_iterations=100,
        #                 max_time = max_time,
        #                 exploration_constant=c,
        #                 k_belief = k_obs,
        #                 alpha_belief = 1/inv_alpha_obs,
        #                 k_action = k_act,
        #                 alpha_action = 1/inv_alpha_act,
        #                 estimate_value=BOMCP.RolloutEstimator(rollout_policy),
        #                 rng=rng
        #                 )

        #     planner = solve(solver, pomdp)
        #     seed!(planner, i+40000)

        elseif solver_name == "VOWSS"
            state_width = x[1]
            action_width = x[2]
            vowss_p = x[3]
            action_width_decay = x[4]
            cont_action_dim = x[5]

            vowss_sig = diagm([0.5, 0.5])
            vowss_epsilon = 0.5/10.0
            voo = VOOActionGenerator(vowss_p, Distributions.Uniform(), MvNormal(vowss_sig), rng, vowss_epsilon)
            vowss_opt = VOWSSOptions(max_depth, state_width, action_width, action_width_decay, true, voo)

            solver = VOOSparseSamplingSolver(vowss_opt, rng)
            planner = solve(solver, deepcopy(pomdp))
            seed!(planner.rng, i+40000)
        end

        best_action = POMDPs.action(planner, sim_wrap.initial_belief)
        output = (solver = solver_name, action_x = best_action[1], action_y = best_action[2], dist = euclidean(best_action, sim_wrap.lqg_solution))
        return output
    end
end

DARESolution, _, _ = ared(A_mat, B_mat, R_mat, Q_mat);
lqg_solution = solve_LQG_one_step(pomdp)
datestring = Dates.format(now(), "e_d_u_Y_HH_MM")

try
    ########## LQR Policy
    if use_lqr_policy
        sims = SimWrap[]
        results_csv = DataFrame(solver=String[], action_x=Float64[], action_y=Float64[], dist=Float64[])
        for solver_name in keys(solver_list)
            for j in 1:n_iters
                sim_wrap = SimWrap(solver_list[solver_name], solver_name, j, max_time, max_depth, deepcopy(pomdp), deepcopy(initial_belief), lqg_solution, true, false, DARESolution)
                push!(sims, sim_wrap)
            end
        end
        println(stderr, "Starting simulations...")
        if parallelize
            results = @showprogress pmap(gen_sim, sims)
        else
            results = @showprogress map(gen_sim, sims)
        end
        results = POMDPSimulators.create_dataframe(results)
        for row in eachrow(results)
            result_row = [row.solver, row.action_x, row.action_y, row.dist]
            push!(results_csv, result_row)
        end
        combined = by(results, :solver) do df
            DataFrame(mean_x=mean(df[:action_x]), mean_y=mean(df[:action_y]), mean_dist=mean(df[:dist]), se_dist=(std(df[:dist])/sqrt(n_iters)))
        end
        for row in eachrow(combined)
            println(stderr, @sprintf("- %s mean action: [%6.3f, %6.3f], dist to LQG: %6.3f ± %6.3f", row.solver, row.mean_x, row.mean_y, row.mean_dist, row.se_dist))
        end

        if save_data
            fname_csv = joinpath(@__DIR__, "../data", "LQG_action_compare_LQR_results"*"_"*datestring*".csv")
            CSV.write(fname_csv, results_csv)
            println(stderr,"Saving parameters to ", fname_csv)
        end
    end

    ########## Saturation Policy
    if use_riccati_policy
        sims = SimWrap[]
        results_csv = DataFrame(solver=String[], action_x=Float64[], action_y=Float64[], dist=Float64[])
        for solver_name in keys(solver_list)
            if use_lqr_policy && solver_name == "VOWSS"
                continue # Don't simulate VOWSS again if doing both types of simulation
            end
            for j in 1:n_iters
                sim_wrap = SimWrap(solver_list[solver_name], solver_name, j, max_time, max_depth, deepcopy(pomdp), deepcopy(initial_belief), lqg_solution, false, true, DARESolution)
                push!(sims, sim_wrap)
            end
        end
        println(stderr, "Starting simulations...")
        if parallelize
            results = @showprogress pmap(gen_sim, sims)
        else
            results = @showprogress map(gen_sim, sims)
        end
        results = POMDPSimulators.create_dataframe(results)
        for row in eachrow(results)
            result_row = [row.solver, row.action_x, row.action_y, row.dist]
            push!(results_csv, result_row)
        end
        combined = by(results, :solver) do df
            DataFrame(mean_x=mean(df[:action_x]), mean_y=mean(df[:action_y]), mean_dist=mean(df[:dist]), se_dist=(std(df[:dist])/sqrt(n_iters)))
        end
        for row in eachrow(combined)
            println(stderr, @sprintf("- %s mean action: [%6.3f, %6.3f], dist to LQG: %6.3f ± %6.3f", row.solver, row.mean_x, row.mean_y, row.mean_dist, row.se_dist))
        end

        if save_data
            fname_csv = joinpath(@__DIR__, "../data", "LQG_action_compare_Riccati_results"*"_"*datestring*".csv")
            CSV.write(fname_csv, results_csv)
            println(stderr,"Saving parameters to ", fname_csv)
        end
    end

    ########## Random Policy
    if use_random_policy
        sims = SimWrap[]
        results_csv = DataFrame(solver=String[], action_x=Float64[], action_y=Float64[], dist=Float64[])
        for solver_name in keys(solver_list)
            if solver_name == "VOWSS"
                continue # Don't simulate VOWSS again if doing both types of simulation
            end
            for j in 1:n_iters
                sim_wrap = SimWrap(solver_list[solver_name], solver_name, j, max_time, max_depth, deepcopy(pomdp), deepcopy(initial_belief), lqg_solution, false, false, DARESolution)
                push!(sims, sim_wrap)
            end
        end
        println(stderr, "Starting simulations...")
        if parallelize
            results = @showprogress pmap(gen_sim, sims)
        else
            results = @showprogress map(gen_sim, sims)
        end
        results = POMDPSimulators.create_dataframe(results)
        for row in eachrow(results)
            result_row = [row.solver, row.action_x, row.action_y, row.dist]
            push!(results_csv, result_row)
        end
        combined = by(results, :solver) do df
            DataFrame(mean_x=mean(df[:action_x]), mean_y=mean(df[:action_y]), mean_dist=mean(df[:dist]), se_dist=(std(df[:dist])/sqrt(n_iters)))
        end
        for row in eachrow(combined)
            println(stderr, @sprintf("- %s mean action: [%6.3f, %6.3f], dist to LQG: %6.3f ± %6.3f", row.solver, row.mean_x, row.mean_y, row.mean_dist, row.se_dist))
        end

        if save_data
            fname_csv = joinpath(@__DIR__, "../data", "LQG_action_compare_Random_results"*"_"*datestring*".csv")
            CSV.write(fname_csv, results_csv)
            println(stderr,"Saving parameters to ", fname_csv)
        end
    end
finally
    Distributed.rmprocs(worker_ids)
end


text:
@testset "503.next-greater-element-ii.jl" begin
    @test next_greater_elements([
        2, 25, 19, 17, 9, 18, 12, 27, 25, 1, 12, 17, 9, 21, 18, 24, 16, 19, 8, 21
    ]) == [25, 27, 27, 18, 18, 27, 27, -1, 27, 12, 17, 21, 21, 24, 24, 25, 19, 21, 21, 25]
    @test next_greater_elements([1, 2, 1]) == [2, -1, 2]
end


text:
"""
    licenses(conn::Connection,
             pat::GitHubPersonalAccessToken,
             schema::AbstractString = "gh_2007_\$(Dates.year(floor(now(), Year) - Day(1)))",
             )::Nothing

Uploads the licenses table to the database.
It includes every OSI-approved license that is machine readable with Licensee.
"""
function licenses()
    # Obtain all licences used by the Ruby gem: `licensee`.
    licensee = graphql(string(strip(replace(String(read(joinpath(@__DIR__, "assets", "graphql", "01_licensee.graphql"))), r"[\n\s]+" => " "))),
                       # The repository is https://github.com/licensee/licensee
                       vars = Dict("id" => "MDEwOlJlcG9zaXRvcnkyMzAyMjM3Nw==",
                                   # Path is https://github.com/licensee/licensee/tree/master/vendor/choosealicense.com/_licenses
                                   "expression" => "master:vendor/choosealicense.com/_licenses")) |>
        (obj -> JSON3.read(obj.Data).data.node.object) |>
        (obj -> (sha1 = obj.oid, spdx = [ uppercase(elem.name[1:end - 4]) for elem in obj.entries ]))
    # Queries the SPDX licenses (ID and whether OSI approved)
    spdx = request("GET",
                   "https://raw.githubusercontent.com/spdx/license-list-data/master/json/licenses.json") |>
        (obj -> JSON3.read(obj.body)) |>
        (obj -> (version = obj.licenseListVersion,
                 release_date = obj.releaseDate,
                 spdx = [ (spdx = license.licenseId,
                           name = license.name)
                           # Keep only OSI-approved licences
                           for license in obj.licenses if license.isOsiApproved ]))
    # Keep only licenses that are machine detectable with Licensee
    filter!(license -> uppercase(license.spdx) ∈ licensee.spdx, spdx.spdx)
    @unpack conn, schema = PARALLELENABLER
    execute(conn, "TRUNCATE TABLE $schema.licenses CASCADE;")
    execute(conn, "BEGIN;")
    load!(spdx.spdx, conn, "INSERT INTO $schema.licenses VALUES ($(join(("\$$i" for i in 1:2), ','))) ON CONFLICT DO NOTHING;")
    execute(conn, "COMMIT;")
    # Add the source data metadata
    execute(
        conn,
        """
        COMMENT ON TABLE $schema.licenses IS
        'OSI-approved machine detectable licenses (i.e., Licensee).
        The official list of Open Source Initiative (OSI) approved licenses is hosted at their website.
        For programmatic access we first obtain all licences from the SPDX data files.
        We filter the licenses based on whether it is OSI-approved.
        GitHub uses the Ruby Gem Licensee for systematically detecting the license of repositories.

        References:
        - https://spdx.org
        - https://licensee.github.io/licensee
        - https://help.github.com/en/github/creating-cloning-and-archiving-repositories/licensing-a-repository#detecting-a-license

        SPDX data version: $(spdx.version)
        Licensee data version: $(licensee.sha1)';
        COMMENT ON COLUMN $schema.licenses.spdx IS 'Software Package Data Exchange License ID';
        COMMENT ON COLUMN $schema.licenses.name IS 'Name of the license';
        """)
    nothing
end


text:
"""
    euclidean_graph(points::Matrix, G; L=1., p=2., cutoff=Inf, bc=:periodic)

Given the `d×N` matrix `points` builds an Euclidean graph of `N` vertices
according to the following procedure.

Defining the `d`-dimensional vectors `x[i] = points[:,i]`, an edge between
vertices `i` and `j` is inserted if `norm(x[i]-x[j], p) < cutoff`.
In case of negative `cutoff` instead every edge is inserted.
For `p=2` we have the standard Euclidean distance.
Set `bc=:periodic` to impose periodic boundary conditions in the box ``[0,L]^d``.
Set `bc=:open` for open boundary condition. In this case the keyword argument `L`
will be ignored.

Returns a graph and Dict containing the distance on each edge.
"""
function euclidean_graph(points::Matrix, ::Type{G} = Graph;
            L=1., p=2., cutoff=Inf, bc=:open) where G<:AGraph
    d, N = size(points)
    g = G(N)
    weights = Dict{Edge,Float64}()
    bc ∉ [:periodic,:open] && error("Not a valid boundary condition.")
    if bc == :periodic
        maximum(points) > L &&  error("Some points are outside the box of size $L.")
    end
    for i=1:N
        for j=i+1:N
            if bc == :open
                Δ = points[:,i] .- points[:,j]
            elseif bc == :periodic
                Δ = abs.(points[:,i] .- points[:,j])
                Δ = min.(L .- Δ, Δ)
            end
            dist = norm(Δ, p)
            if dist < cutoff
                e = Edge(i,j)
                add_edge!(g, e)
                weights[e] = dist
            end
        end
    end
    return g, weights
end


text:
# Baumgartner-Weis-Schindler statistic and hypothesis testing.
# Copyright (c) Robert Feldt 2013 (robert.feldt@gmail.com)
#
# For license details see the LICENSE.md file in the top dir of the repository
# where you found this file.
#
# The algorithms in this file are based on the definitions of the BWS in the paper:
#
#  Markus Neuhäuser, "Exact tests based on the Baumgartner-WeiB-Schindler statistic - A survey",
#  Statistical Papers, vol. 46, pp. 1-30, Springer Verlag, 2005.
#

export bws_statistic, bws_test_sampled

#############################################################################
## Part I. Calculating the BWS.
#############################################################################

# Calculate the BWS statistic from two samples of values.
function bws_statistic(x, y)
  # Get the ranks. The first 1:n are ranks for x, (n+1):end for y.
  rs = ranks(vcat(x, y))

  FeldtLib.bws_statistic_from_ranks(rs, length(x), length(y))
end

# Test if according to the Baumgartner-Weis-Schindler (BWS) statistic we
# can reject the hypothesis that the two set of samples x and y are likely to have
# the same location (mean).
function bws_test_sampled(x, y, numsamples = 10000)
  # Get the ranks. The first 1:n are ranks for x, (n+1):(n+m) for y.
  rs = ranks(vcat(x, y))

  FeldtLib.bws_test_sampled_from_ranks(rs, length(x), length(y), numsamples)
end

# Calculate the BWS statistic from a set of ranks of n values from one sample
# and m values from another sample. The first 1:n rank values in ranks are for
# the first sample and the rest are from the other.
function bws_statistic_from_ranks(ranks, n, m)
  # Pre-calc some values needed in both b calcs.
  nm = n+m
  np1 = n + 1
  mp1 = m + 1

  # Now calc the two b values
  b_x = (1/n) * (n / (m*nm)) * calc_bws_sum(ranks,   1,  n, n, np1, nm, 0)
  b_y = (1/m) * (m / (n*nm)) * calc_bws_sum(ranks, np1, nm, m, mp1, nm, n)

  (b_x + b_y) / 2
end

function calc_bws_sum(ranks, mini, maxi, n, np1, nm, subi = 0)
  sum = 0.0
  for(i in mini:maxi)
    k = i - subi
    nom = ( ranks[i] - (nm/n*k) )^2
    denom = (k/np1) * (1 - (k/np1))
    sum += (nom / denom)
  end
  sum
end

function bws_test_sampled_from_ranks(ranks, n, m, numsamples = 10000)
  rs = copy(ranks) # Copy since we will shuffle it below.

  actual_b = FeldtLib.bws_statistic_from_ranks(rs, n, m)

  # Array to save the b statistics
  samples = Array(Float64, numsamples)
  numlarger = 0

  # Set up for inline sorts. Not sure if this really speeds things up...
  ns = Array(Int64, n)
  ms = Array(Int64, m)
  np1 = n+1

  for(i in 1:numsamples)
    # Now shuffle the ranks and then ensure they are sorted for each group.
    # Not sure this is actually faster than letting julia decide on the copying
    # and gc of intermediate vectors/arrays... ;)
    shuffle!(rs)
    ns[1:n] = rs[1:n]
    ms[1:m] = rs[np1:end]
    sort!(ns)
    sort!(ms)
    rs[1:n] = ns
    rs[np1:end] = ms

    # Now we can calc the bws again and note if it was larger than before.
    samples[i] = FeldtLib.bws_statistic_from_ranks(rs, n, m)
    if samples[i] >= actual_b
      numlarger += 1
    end
  end

  # A larger (or equal) BWS statistic value indicates samples from distributions
  # that are "more apart" and thus that there is a difference between the samples.
  pvalue = numlarger / numsamples

  (actual_b, pvalue, numlarger, numsamples) # , mean(samples), std(samples))
end

#############################################################################
## Part II. Adding the same interface as used by hypothesis tests in
##          HypothesisTests.jl
#############################################################################
using HypothesisTests

type BWSTest <: HypothesisTests.HypothesisTest
  # Ranks of all values where first 1:n are from one sample and (n+1):end are
  # from the other.
  ranks::Array{Int64,1}
  n::Integer
  m::Integer

  b::Float64 # The BWS statistic itself

  # We save the numlarger and numsamples taken instead of the p-value since
  # we can then add additional iterations later if needed.
  numlarger::Integer
  numsamples::Integer
end
function BWSTest{S <: Number}(x::Vector{S}, y::Vector{S}, numsamples = 10000)
  rs = ranks(vcat(x, y))
  n = length(x)
  m = length(y)
  b, pvalue, numlarger, numsamples = bws_test_sampled_from_ranks(rs, n, m)
  BWSTest(rs, n, m, b, numlarger, numsamples)
end

function pvalue(bwst::BWSTest)
  bwst.numlarger / bwst.numsamples
end


text:
using Bridge
struct EulerMaruyamaBounded <: Bridge.SDESolver end
const EulerMaruyama = EulerMaruyamaBounded

_DIR = "general"
include(joinpath(_DIR, "types.jl"))
include(joinpath(_DIR, "coordinate_access.jl"))
include(joinpath(_DIR, "readjustments.jl"))

_DIR = "stochastic_process"
include(joinpath(_DIR, "bounded_diffusion_domain.jl"))
include(joinpath(_DIR, "guid_prop_bridge.jl"))
include(joinpath(_DIR, "path_to_wiener.jl"))

_DIR = "solvers"
include(joinpath(_DIR, "vern7.jl"))
include(joinpath(_DIR, "tsit5.jl"))
include(joinpath(_DIR, "rk4.jl"))
include(joinpath(_DIR, "ralston3.jl"))
include(joinpath(_DIR, "euler_maruyama_dom_restr.jl"))

_DIR = "transition_kernels"
include(joinpath(_DIR, "random_walk.jl"))
include(joinpath(_DIR, "pseudo_conjugate.jl"))

_DIR = "mcmc_extras"
include(joinpath(_DIR, "adaptation.jl"))
include(joinpath(_DIR, "blocking.jl"))
include(joinpath(_DIR, "first_passage_times.jl"))
include(joinpath(_DIR, "starting_pt.jl"))

_DIR = "examples"
include(joinpath(_DIR, "fitzhugh_nagumo.jl"))
include(joinpath(_DIR, "radial_ornstein_uhlenbeck.jl"))
include(joinpath(_DIR, "lorenz_system.jl"))
include(joinpath(_DIR, "lorenz_system_const_vola.jl"))
include(joinpath(_DIR, "prokaryotic_autoregulatory_gene_network.jl"))
include(joinpath(_DIR, "lotka_volterra.jl"))

_DIR = "mcmc"
include(joinpath(_DIR, "priors.jl"))
include(joinpath(_DIR, "mcmc_components.jl"))
include(joinpath(_DIR, "setup.jl"))
include(joinpath(_DIR, "workspace.jl"))
include(joinpath(_DIR, "fusion.jl"))
include(joinpath(_DIR, "conjugate_updt.jl"))
include(joinpath(_DIR, "updates.jl"))
include(joinpath(_DIR, "mcmc.jl"))
include(joinpath(_DIR, "repeated.jl"))