text:
using Plots, Dates
a=["
23.908918080036315" "17:27:52.574";
"14.463959533978853" "17:27:52.861";
"18.202118982566667" "17:27:52.863";
"18.39031429901463" "17:27:52.864";
"16.325316060988722" "17:27:52.866";
"14.953975221442477" "17:27:52.867";
"14.359751560875704" "17:27:52.869";
"14.367776750909542" "17:27:52.87";
"13.954831110639821" "17:27:52.872";
"12.467447763806167" "17:27:52.873";
"11.95376710676135" "17:27:52.874";
"12.53474000766882" "17:27:52.876";
"11.51819201606951" "17:27:52.877";
"10.956046841355052" "17:27:52.878";
"11.279970412249579" "17:27:52.88";
"10.66473902526741" "17:27:52.881";
"10.152806325097346" "17:27:52.888";
"9.98390970086297" "17:27:52.89";
"9.844062663026032" "17:27:52.891";
"9.475675910733608" "17:27:52.893";
"55.450256" "17:30:10";
"63.70554" "17:30:10.428";
"76.23911" "17:30:10.74";
"90.1125" "17:30:11.148";
"89.81146" "17:30:11.431";
"100.01215" "17:30:13.905";
"105.650116" "17:30:15.737";
"108.89857" "17:30:16.086";
"110.386505" "17:30:16.431";
"110.86751" "17:30:16.736";
"110.638626" "17:30:17.03";
"109.861" "17:30:17.348";
"109.583786" "17:30:17.677";
"108.80472" "17:30:17.997";
"107.59737" "17:30:18.307";
"105.99181" "17:30:18.654";
"103.979004" "17:30:18.948";
"101.57057" "17:30:19.273";
"100.29155" "17:30:19.531";
"99.91356" "17:30:19.853";
"102.89426906310462" "17:30:19.859";
"72.4867903525124" "17:30:19.86";
"46.692183055353894" "17:30:19.862";
"26.958299441974365" "17:30:19.863";
"20.988848793598592" "17:30:19.865";
"27.5191529345425" "17:30:19.867";
"31.77945783687218" "17:30:19.868";
"32.19261654120688" "17:30:19.87";
"29.85377965190008" "17:30:19.871";
"25.631602189470126" "17:30:19.873";
"20.44294972252166" "17:30:19.874";
"15.961912012753626" "17:30:19.876";
"14.349491723574411" "17:30:19.877";
"15.701582686038629" "17:30:19.879";
"17.455227408685786" "17:30:19.88";
"17.947099309632282" "17:30:19.882";
"16.955630325454738" "17:30:19.883";
"15.055405773273886" "17:30:19.885";
"13.173866812931864" "17:30:19.887";
"12.066501027763492" "17:30:19.888";
"63.531586" "17:30:20.793";
"70.3219" "17:30:21.2";
"74.03856" "17:30:21.511";
"75.356346" "17:30:21.878";
"73.90744" "17:30:22.313";
"71.580154" "17:30:22.698";
"73.09373" "17:30:23.123";
"70.93234" "17:30:23.453";
"66.97326" "17:30:23.823";
"59.184372" "17:30:24.124";
"64.92076" "17:30:24.52";
"66.09159" "17:30:24.987";
"59.282707" "17:30:25.432";
"56.359097" "17:30:25.818";
"57.77132" "17:30:26.273";
"45.81396" "17:30:26.752";
"45.47974" "17:30:27.085";
"48.293953" "17:30:27.461";
"44.112705" "17:30:27.859";
"38.33211" "17:30:28.219";
"29.1106144112681" "17:30:28.221";
"16.863129246135824" "17:30:28.223";
"21.469221296822415" "17:30:28.224";
"24.451504471126327" "17:30:28.226";
"22.641558033904307" "17:30:28.228";
"18.358338817176094" "17:30:28.231";
"13.97222560636477" "17:30:28.232";
"11.742127581307642" "17:30:28.234";
"11.829261595741386" "17:30:28.236";
"12.144612967121358" "17:30:28.238";
"11.722685275885324" "17:30:28.24";
"10.795510934289448" "17:30:28.241";
"9.505458690477697" "17:30:28.243";
"8.268657781952445" "17:30:28.245";
"7.809709344835771" "17:30:28.247";
"8.409884375558383" "17:30:28.249";
"8.705589685159495" "17:30:28.251";
"8.094973905310225" "17:30:28.253";
"7.138194661308531" "17:30:28.254";
"6.451720566500355" "17:30:28.256";
"44.72667" "17:30:29.937";
"54.03601" "17:30:30.365";
"58.72908" "17:30:30.916";
"58.7531" "17:30:31.341";
"55.293983" "17:30:31.807";
"54.283268" "17:30:32.213";
"50.8247" "17:30:32.528";
"52.30549" "17:30:32.93";
"51.832645" "17:30:33.321";
"48.467236" "17:30:33.721";
"48.711605" "17:30:34.088";
"50.039333" "17:30:34.561";
"48.463642" "17:30:34.909";
"46.111095" "17:30:35.334";
"47.391342" "17:30:35.733";
"45.76744" "17:30:36.185";
"45.54458" "17:30:36.55";
"42.452587" "17:30:36.934";
"43.378273" "17:30:37.314";
"40.86766" "17:30:37.801";
"34.224546541904104" "17:30:37.802";
"20.690081533542763" "17:30:37.804";
"16.87089184783634" "17:30:37.808";
"19.41019725665371" "17:30:37.81";
"21.010021773896128" "17:30:37.811";
"18.786579205736082" "17:30:37.812";
"15.202535337026715" "17:30:37.814";
"12.334665054757062" "17:30:37.815";
"10.869808034602398" "17:30:37.817";
"10.213267972654487" "17:30:37.818";
"10.925783317270849" "17:30:37.82";
"10.569897780388688" "17:30:37.822";
"9.461448664958153" "17:30:37.823";
"9.088373175851764" "17:30:37.825";
"9.245948206079122" "17:30:37.826";
"8.53288385631793" "17:30:37.828";
"7.558470341474014" "17:30:37.829";
"6.508280726062224" "17:30:37.831";
"6.481304605910951" "17:30:37.833";
"7.2666132284831955" "17:30:37.835";
"41.850376" "17:30:38.674";
"66.26185" "17:30:39.372";
"54.60947" "17:30:39.963";
"32.510094" "17:30:40.606";
"51.074642" "17:30:41.063";
"28.664633" "17:30:41.755";
"30.306587" "17:30:42.22";
"15.147215" "17:30:43.039";
"23.782793" "17:30:43.606";
"19.002993" "17:30:44.461";
"8.787821" "17:30:45.269";
"13.101144" "17:30:45.776";
"13.728193" "17:30:46.577";
"6.2133594" "17:30:47.421";
"4.916186" "17:30:48.184";
"8.662348" "17:30:48.887";
"5.077702" "17:30:49.728";
"5.6022034" "17:30:50.426";
"8.474883" "17:30:51.172";
"4.9741273" "17:30:51.958";
"18.77114034968427" "17:30:51.959";
"14.046343415577093" "17:30:51.961";
"12.306814220721247" "17:30:51.962";
"10.005431410399392" "17:30:51.963";
"8.525267166193675" "17:30:51.964";
"10.869814181719663" "17:30:51.966";
"9.574841259740191" "17:30:51.967";
"6.1687697525020795" "17:30:51.968";
"7.728312099892402" "17:30:51.969";
"8.25230175374919" "17:30:51.971";
"6.541301837158007" "17:30:51.975";
"6.658201376578983" "17:30:51.977";
"6.546135258265819" "17:30:51.978";
"5.868650712480528" "17:30:51.979";
"6.318388838599809" "17:30:51.981";
"5.616184872491648" "17:30:51.982";
"4.971399821761165" "17:30:51.984";
"5.9276363567125205" "17:30:51.985";
"5.423998062420753" "17:30:51.987";
"4.738778674564512" "17:30:51.988";
"47.77442" "17:30:52.872";
"76.51193" "17:30:53.759";
"56.237415" "17:30:54.691";
"51.050514" "17:30:55.468";
"36.433197" "17:30:56.326";
"32.724205" "17:30:57.01";
"32.84369" "17:30:57.671";
"26.49535" "17:30:58.473";
"20.50503" "17:30:59.181";
"14.654626" "17:30:59.812";
"14.120778" "17:31:00.392";
"13.615862" "17:31:01.029";
"11.383756" "17:31:01.69";
"11.854716" "17:31:02.449";
"7.492182" "17:31:03.31";
"8.005376" "17:31:04.139";
"9.380755" "17:31:05.101";
"8.147045" "17:31:06.047";
"7.7330875" "17:31:06.722";
"6.198367" "17:31:07.602"
]

function parse_dates(x)
    t=split(x,(':','.'))
    w=map(x->parse(Float64,x),t)
    Dates.Time(w...)
end
function norm_time(x)
    abs_times = x
    start_time = abs_times[1]
    norm_times = Vector(undef,length(abs_times))
    for i in 1:length(abs_times)
        norm_times[i] = abs_times[i]-start_time
    end
    norm_times
end

loss_values = map(x->parse(Float64,x),a[:,1])
time_values = map(x->parse_dates(x),a[:,2])
#norm_time_vals = norm_time(time_values)

print(time_values[200]-time_values[1])
plot(time_values[1:20],loss_values[1:20],linecolor="blue")
plot!(time_values[21:40],loss_values[21:40],linecolor="red")
plot!(time_values[41:60],loss_values[41:60],linecolor="blue")
plot!(time_values[61:80],loss_values[61:80],linecolor="red")
plot!(time_values[81:100],loss_values[81:100],linecolor="blue")
plot!(time_values[101:120],loss_values[101:120],linecolor="red")
plot!(time_values[121:140],loss_values[121:140],linecolor="blue")
plot!(time_values[141:160],loss_values[141:160],linecolor="red")
plot!(time_values[161:180],loss_values[161:180],linecolor="blue")
plot!(time_values[181:200],loss_values[181:200],linecolor="red")
#Int64(norm_time_vals[1])

print(loss_values[200])


text:
import Winston.*

m=1e-3
u=1e-6
n=1e-9
p=1e-12
f=1e-15
k=1e3
K=1e3
Meg=1e6
M=1e6
MEG=1e6
meg=1e6
G=1e-9

vdiv(z1, z2)= (z2 + 0im) ./ (z1 + z2)
para(z1, z2)= (z1 + 0im) .* vdiva(z1, z2)
magn(x)=hypot(real(x), imag(x))
db(x)= 20 * log10(magn(x))
logsweep(x1, x2, npts)=logspace(log10(x1), log10(x2), npts)
zc(c,freq)= -(1im / (2pi * c)) ./ freq
zl(l,freq)= 1im * 2pi * l .* freq
zr(r,freq)= ones(Integer(length(freq))) * r

fstart = 1
fstop = 10meg
npts = 2810
freq = logsweep(fstart, fstop, npts)

# filt = db(magn(vdiv(zr(1k,freq), zc(.01u,freq))))
filt = db(vdiv(zr(100,freq)+zl(1m,freq), zc(.01u,freq)))

p = FramedPlot()
setattr(p, "xlog", true)
setattr(p, "xrange", (fstart, fstop))
setattr(p.frame1, "draw_grid", true)
setattr(p.x1, "label", "Frequency, Hz")
setattr(p.y1, "label", "LPF out, dB")

a = Curve(freq, filt, "color", "blue")
add(p, a)

ext: jl

max_stars_repo_name: tomacorp/PyPost

max_issues_repo_name: tomacorp/PyPost

max_forks_repo_name: tomacorp/PyPost


text:
"""
Formatting of incoming Redis Replies
"""
function getline(s::TCPSocket)
    l = chomp(readline(s))
    length(l) > 1 || throw(ProtocolException("Invalid response received: $l"))
    return l
end

function read_reply(conn::RedisConnectionBase)
    l = getline(conn.socket)
    reply = parseline(l, conn.socket)
    return reply
end

parse_error(l::AbstractString) = throw(ServerException(l))

function parse_bulk_string(s::TCPSocket, slen::Int)
    b = read(s, UInt8, slen+2) # add crlf
    if length(b) != slen + 2
        throw(ProtocolException(
            "Bulk string read error: expected $len bytes; received $(length(b))"
        ))
    else
        return bytestring(b[1:end-2])
    end
end

function parse_array(s::TCPSocket, slen::Int)
    a = Array{Any, 1}(slen)
    for i = 1:slen
        l = getline(s)
        r = parseline(l, s)
        a[i] = r
    end
    return a
end

function parseline(l::AbstractString, s::TCPSocket)
    reply_type = l[1]
    reply_token = l[2:end]
    if reply_type == '+'
        return reply_token
    elseif reply_type == ':'
        parse(Int, reply_token)
    elseif reply_type == '$'
        slen = parse(Int, reply_token)
        if slen == -1
            Nullable{AbstractString}()
        else
            parse_bulk_string(s, slen)
        end
    elseif reply_type == '*'
        slen = parse(Int, reply_token)
        if slen == -1
            Nullable{AbstractString}()
        else
            parse_array(s, slen)
        end
    elseif reply_type == '-'
        parse_error(reply_token)
    end
end

"""
Formatting of outgoing commands to the Redis server
"""
function pack_command(command)
    packed_command = "*$(length(command))\r\n"
    for token in command
        ltoken = ifelse(typeof(token) <: Number, length(string(token)), length(token))
        packed_command = string(packed_command, "\$$(ltoken)\r\n", token, "\r\n")
    end
    packed_command
end

function execute_command_without_reply(conn::RedisConnectionBase, command)
    is_connected(conn) || throw(ConnectionException("Socket is disconnected"))
    send_command(conn, pack_command(command))
end

function execute_command(conn::RedisConnectionBase, command)
    execute_command_without_reply(conn, command)
    read_reply(conn)
end

baremodule SubscriptionMessageType
    const Message = 0
    const Pmessage = 1
    const Other = 2
end

immutable SubscriptionMessage
    message_type
    channel::AbstractString
    message::AbstractString

    function SubscriptionMessage(reply::AbstractArray)
        notification = reply
        message_type = notification[1]
        if message_type == "message"
            new(SubscriptionMessageType.Message, notification[2], notification[3])
        elseif message_type == "pmessage"
            new(SubscriptionMessageType.Pmessage, notification[2], notification[4])
        else
            new(SubscriptionMessageType.Other, "", "")
        end
    end
end


text:
@testset "Test time series data application" begin
    sys = PSB.build_system(PSB.PSITestSystems, "c_sys5")
    pmi_data = PMI.get_pm_data(sys)

    mn_data =
        PMI.apply_time_series(pmi_data, sys, last(PSY.get_forecast_initial_times(sys)), 3:5)

    @test mn_data["multinetwork"]
    @test length(mn_data["nw"]) == 3

    pmi_mn_data = PMI.get_pm_data(
        sys,
        start_time = last(PSY.get_forecast_initial_times(sys)),
        time_periods = 3:5,
    )

    @test mn_data == pmi_mn_data

    tp_data = PMI.get_pm_data(
        sys,
        start_time = first(PSY.get_forecast_initial_times(sys)),
        period = 5,
    )

    @test tp_data != pmi_data

    PMI.apply_time_period!(pmi_data, sys, first(PSY.get_forecast_initial_times(sys)), 5)

    @test tp_data == pmi_data

    #TODO: add tests to verify data is being applied correctly
end

max_stars_repo_name: NREL-SIIP/PowerModelsInterface

max_issues_repo_name: NREL-SIIP/PowerModelsInterface

max_forks_repo_name: NREL-SIIP/PowerModelsInterface


text:
module FastHistograms

import Base: eltype, @propagate_inbounds
using ComputedFieldTypes
using LoopVectorization
using StaticArrays

include("traits.jl")
include("bin_search.jl")
include("bin_update.jl")
include("real_histogram.jl")
include("text_histogram.jl")

export create_fast_histogram, bin_search, increment_bins!, counts, zero!

"""
    create_fast_histogram(
        ::BinType,
        ::BinSearchAlgorithm,
        ::HistogramParallelization,
        args...
    )

Creates a histogram with the given `BinType`, `BinSearchAlgorithm`, and `HistogramParallelization` traits.
Methods of this function will also require additional arguments (here `args...`) that depend on the combination of
traits selected.
"""
create_fast_histogram

"""
FastHistograms declares and implements a minimal histogram interface with a focus on speed.

```julia-repl
julia> using FastHistograms, Random

# Create a 2D histogram for 8-bit integer data.
julia> h = create_fast_histogram(
    # Use fixed-width bins with an optimized bin search algorithm (Arithmetic)
    #  for fixed-width bins.
    FastHistograms.FixedWidth(),
    FastHistograms.Arithmetic(),
    # Don't use any parallelization because our data are small.
    FastHistograms.NoParallelization(),
    [(0x00, 0xff, 4), (0x00, 0xff, 4)],
);

# Create two random images to compute the joint histogram for
julia> img1 = rand(0x00:0xff, 32, 32);

julia> img2 = rand(0x00:0xff, 32, 32);

# Compute the histogram bin counts
julia> increment_bins!(h, img1, img2)

# Get the bin counts
julia> counts(h)
4×4 Matrix{Int64}:
 61  64  67  64
 65  59  72  65
 61  66  71  61
 53  67  63  65
```
"""
FastHistograms

end


text:
@testset "DeepONet" begin
    @testset "proper construction" begin
        deeponet = DeepONet((32,64,72), (24,48,72), σ, tanh)
        # Branch net
        @test size(deeponet.branch_net.layers[end].weight) == (72,64)
        @test size(deeponet.branch_net.layers[end].bias) == (72,)
        # Trunk net
        @test size(deeponet.trunk_net.layers[end].weight) == (72,48)
        @test size(deeponet.trunk_net.layers[end].bias) == (72,)
    end

    # Accept only Int as architecture parameters
    @test_throws MethodError DeepONet((32.5,64,72), (24,48,72), σ, tanh)
    @test_throws MethodError DeepONet((32,64,72), (24.1,48,72))

    # Just the first 16 datapoints from the Burgers' equation dataset
    a = [0.83541104, 0.83479851, 0.83404712, 0.83315711, 0.83212979, 0.83096755,
         0.82967374, 0.82825263, 0.82670928, 0.82504949, 0.82327962, 0.82140651,
         0.81943734, 0.81737952, 0.8152405, 0.81302771]
    sensors = collect(range(0, 1, length=16)')
    model = DeepONet((16, 22, 30), (1, 16, 24, 30), σ, tanh; init_branch=Flux.glorot_normal, bias_trunk=false)
    y = model(a, sensors)
    @test size(y) == (1, 16)

    mgrad = Flux.Zygote.gradient(() -> sum(model(a, sensors)), Flux.params(model))
    @test length(mgrad.grads) == 7

    # Output size of branch and trunk subnets should be same
    branch = Chain(Dense(16, 22), Dense(22, 30))
    trunk = Chain(Dense(1, 16), Dense(16, 24), Dense(24, 32))
    m = DeepONet(branch, trunk)
    @test_throws AssertionError DeepONet((32,64,70), (24,48,72), σ, tanh)
    @test_throws DimensionMismatch m(a, sensors)
end


text:
#####
##### Resampling
#####

resample!(::Nothing, args...) = nothing

struct Resampler{D}
    only_failed_particles :: Bool
    acceptable_failure_fraction :: Float64
    resample_failure_fraction :: Float64
    distribution :: D
    verbose :: Bool
end

function Resampler(; only_failed_particles = true,
                     acceptable_failure_fraction = 0.8,
                     resample_failure_fraction = 0.2,
                     verbose = false,
                     distribution = FullEnsembleDistribution())

    return Resampler(only_failed_particles,
                     acceptable_failure_fraction,
                     resample_failure_fraction,
                     distribution,
                     verbose)
end

abstract type EnsembleDistribution end

function ensemble_normal_distribution(X)
    μ = [mean(X, dims=2)...]
    Σ = cov(X, dims=2)
    return MvNormal(μ, Σ)
end

struct FullEnsembleDistribution <: EnsembleDistribution end
(::FullEnsembleDistribution)(X, G, successes) = ensemble_normal_distribution(X)

struct SuccessfulEnsembleDistribution <: EnsembleDistribution end
(::SuccessfulEnsembleDistribution)(X, G, successes) = ensemble_normal_distribution(X[:, successes])

function failed_particle_str(θ, k, error=nothing)
    first = string(@sprintf(" particle % 3d: ", k), param_str.(values(θ[k]))...)
    error_str = isnothing(error) ? "" : @sprintf(" error = %.6e", error)
    return string(first, error_str, '\n')
end

"""
    resample!(resampler::Resampler, X, G, eki)

Resamples the parameters `X` of the `eki` process based on the number of failed particles.
"""
function resample!(resampler::Resampler, X, G, eki)
    # `Nensemble` vector of bits indicating whether an ensemble member has failed.
    particle_failure = eki.mark_failed_particles(G)
    failures = findall(particle_failure) # indices of failed particles
    Nfailures = length(failures)
    failed_fraction = Nfailures / size(X, 2)

    if failed_fraction > 0
        # Print a nice message
        particles = Nfailures == 1 ? "particle" : "particles"

        priors = eki.inverse_problem.free_parameters.priors
        θ = transform_to_constrained(priors, X)
        failed_parameters_message = string("               ",  param_str.(keys(priors))..., '\n',
                                           (failed_particle_str(θ, k) for k in failures)...)

        @warn("""
              The forward map for $Nfailures $particles ($(100failed_fraction)%) failed.
              The failed particles are:
              $failed_parameters_message
              """)
    end

    if failed_fraction > resampler.acceptable_failure_fraction
        error("The forward map for $Nfailures particles ($(100failed_fraction)%) failed. Consider \n" *
              "    1. Increasing `Resampler.acceptable_failure_fraction` for \n" *
              "         EnsembleKalmanInversion.resampler::Resampler \n" *
              "    2. Reducing the time-step for `InverseProblem.simulation`, \n" *
              "    3. Evolving `InverseProblem.simulation` for less time \n" *
              "    4. Narrowing `FreeParameters` priors.")

    elseif failed_fraction >= resampler.resample_failure_fraction || !(resampler.only_failed_particles)
        # We are resampling!

        if resampler.only_failed_particles
            Nsample = Nfailures
            replace_columns = failures

        else # resample everything
            Nsample = size(G, 2)
            replace_columns = Colon()
        end

        found_X, found_G = find_successful_particles(eki, X, G, Nsample)

        @info "Replacing columns $replace_columns..."
        view(X, :, replace_columns) .= found_X
        view(G, :, replace_columns) .= found_G

        # Sanity...
        if resampler.verbose && resampler.only_failed_particles # print a helpful message about the failure replacements
            Nobs, Nensemble = size(G)
            y = eki.mapped_observations
            errors = [mapreduce((x, y) -> (x - y)^2, +, y, view(G, :, k)) / Nobs for k in failures]

            priors = eki.inverse_problem.free_parameters.priors
            new_θ = transform_to_constrained(priors, X)

            particle_strings = [failed_particle_str(new_θ, k, errors[i]) for (i, k) in enumerate(failures)]
            failed_parameters_message = string("               ",  param_str.(keys(priors))..., '\n',
                                               particle_strings...)

            @info """
            The replacements for failed particles are
            $failed_parameters_message
            """
        end
    end

    return nothing
end

"""
     find_successful_particles(eki, X, G, Nsample)

Generate `Nsample` new particles sampled from a multivariate Normal distribution parameterized
by the ensemble mean and covariance computed based on the `Nθ` × `Nensemble` ensemble
array `θ`, under the condition that all `Nsample` particles produce successful forward map
outputs.

`G` (`size(G) =  Noutput × Nensemble`) is the forward map output produced by `θ`.

Returns `Nθ × Nsample` parameter `Array` and `Noutput × Nsample` forward map output `Array`.
"""
function find_successful_particles(eki, X, G, Nsample)
    Nθ, Nensemble = size(X)
    Noutput = size(G, 1)

    Nfound = 0
    found_X = zeros(Nθ, 0)
    found_G = zeros(Noutput, 0)

    mark_failed_particles = eki.mark_failed_particles
    particle_failure = mark_failed_particles(G)
    successful_particles = findall(.!particle_failure)
    existing_sample_distribution = eki.resampler.distribution(X, G, successful_particles)

    while Nfound < Nsample
        @info "Searching for successful particles (found $Nfound of $Nsample)..."

        # Generate `Nensemble` new samples in unconstrained space.
        # Note that eki.inverse_problem.simulation
        # must run `Nensemble` particles no matter what.
        X_sample = rand(existing_sample_distribution, Nensemble)

        G_sample = inverting_forward_map(eki.inverse_problem, X_sample)

        particle_failure = mark_failed_particles(G_sample)
        success_columns = findall(.!particle_failure)
        @info "    ... found $(length(success_columns)) successful particles."

        found_X = cat(found_X, X_sample[:, success_columns], dims=2)
        found_G = cat(found_G, G_sample[:, success_columns], dims=2)
        Nfound = size(found_X, 2)
    end

    # Restrict found particles to requested size
    return found_X[:, 1:Nsample], found_G[:, 1:Nsample]
end


text:
## Code for the LEXICAL analysis ###################################################################################
import ColorTypes: RGB
import Base.copy
WHITESPACE = " \t\n\r"
SYMBOLS = "()<>[],*"

# location of a single Char:
mutable struct SourceLocation
    file_name::String
    line_num::Int
    col_num::Int

    SourceLocation(file_name::String = "",
                   line_num::Int32 = Int32(0),
                   col_num::Int32 = Int32(0)) =
                   new(file_name, line_num, col_num)
end

# different types of Tokens:
mutable struct Stop
    loc::SourceLocation
end

mutable struct Identifier
    loc::SourceLocation
    s::Union{String, Char}
end

mutable struct LiteralString
    loc::SourceLocation
    s::String
end

mutable struct LiteralNumber
    loc::SourceLocation
    number::Float32
end

mutable struct Symbol
    loc::SourceLocation
    symbol::Char
end

@enum KeywordEnum begin
    NEW = 1
    MATERIAL = 2
    PLANE = 3
    SPHERE = 4
    DIFFUSE = 5
    SPECULAR = 6
    UNIFORM = 7
    CHECKERED = 8
    IMAGE = 9
    IDENTITY = 10
    TRANSLATION = 11
    ROTATION_X = 12
    ROTATION_Y = 13
    ROTATION_Z = 14
    SCALING = 15
    CAMERA = 16
    ORTHOGONAL = 17
    PERSPECTIVE = 18
    FLOAT = 19
    AABOX = 20
    INT = 21
    LIGHTPOINT = 22
end

mutable struct Keyword
    loc::SourceLocation
    keyword::KeywordEnum
end

# generic Token:
mutable struct Token
    loc::SourceLocation
    value::Union{LiteralNumber, LiteralString, Keyword, Identifier, Symbol, Stop}
end

"""
    InputStream(stream,
                location,
                saved_char,
                save_location,
                saved_token,
                tabulation)
It creates a Input Stream, that can be passed to the READ FUNCTION.

## Arguments:
- `stream` is the text read (IO!);
- `location` is the SourceLocation;
- `saved_char` is the currently read **Char** (functionault = '€');
- `save_location` is the location of the saved_char;
- `saved_token` is the Token associated to the saved_char (it can be a **Token** or the functionault value **Nothing**);
- `tabulation` is the number of spaces that form a tabulation (functionault = 8).
"""
mutable struct InputStream
    stream::IO
    location::SourceLocation
    saved_char::Char
    saved_location::SourceLocation
    saved_token::Union{Token, Nothing}
    tabulation::Int

    InputStream( stream::IO,
                 location::SourceLocation = SourceLocation("", Int32(1), Int32(1)),
                 saved_char = '€',
                 saved_location = location,
                 saved_token = nothing,
                 tabulation = 8) =
                 new(stream, location, saved_char, saved_location, saved_token, tabulation)
end

"""
    GrammarError(msg::String, loc::SourceLocation)

New error message to dislpay when something went wrong in the lexical review.
"""
struct GrammarError <: Exception
    msg::String

    GrammarError(msg::String,
                 loc::SourceLocation = SourceLocation() ) =
                 new("$msg Stacktrace: in $(loc.file_name) at line $(loc.line_num) : $(loc.col_num)")
end

# function to create a copy of a SourceLocation:
copy(location::SourceLocation) = SourceLocation(location.file_name, Int32(location.line_num), Int32(location.col_num))

## FUNCTIONS ####################################################################################################################

# increment the positional indexes:
function _update_pos(stream::InputStream, ch::Char)
    if ch == '€'
        return
    elseif ch == '\n'
        stream.location.line_num += 1
        stream.location.col_num = 1
    elseif ch == '\t'
        stream.location.col_num += stream.tabulations
    else
        stream.location.col_num += 1
    end
end

"""
    read_char(stream::InputStream)

It reads Char one by one, updating the positional indexes.
"""
function read_char(stream::InputStream)

    if stream.saved_char != '€'
        ch = stream.saved_char
        stream.saved_char = '€'

    elseif eof(stream.stream)
        ch = '€'
    else
        ch = read(stream.stream, Char)
    end

    stream.saved_location = copy(stream.location)
    _update_pos(stream, ch)

    return ch
end

"""
    unread_char(stream::InputStream, ch::Char)

After reading the Char, it has to be unread, saving the location.
"""
function unread_char(stream::InputStream, ch)
    while true
        stream.saved_char == '€' && break
    end
    stream.saved_char = ch
    stream.location = copy(stream.saved_location)
end

# ignore all the whitespaces and comments:
function skip_whitespaces_and_comments(stream::InputStream)
    ch = read_char(stream)
    while occursin(ch, WHITESPACE) || ch == '#'
        if ch == '#'
            while (read_char(stream) in ['\r', '\n', '€']) == false
                nothing
            end
        end

        ch =read_char(stream)

        if ch == '€'
            return
        end
    end
    unread_char(stream, ch)
end

# is it a STRING?
function _parse_string_token(stream, token_location::SourceLocation)
    token = ""
    while true
        ch =read_char(stream)
        if ch == '"'
            break
        end
        if ch == '€'
            throw(GrammarError("Unterminated string", token_location))
        end
        token *= ch
    end
    return Token(token_location, LiteralString(token_location, token)) ## it returns a TOKEN!
end

# is it a FLOAT?
function _parse_float_token(stream, first_char::Char, token_location::SourceLocation)
    token = first_char
    while true
        ch = read_char(stream)
        if (isdigit(ch) || ch == '.' || ch in ['e', 'E']) == false
            unread_char(stream, ch)
            break
        end
        token *= ch
    end

    try
        token = parse(Float32, string(token))
    catch e
        throw(GrammarError("$(token) is an invalid floating-point number", token_location))
    end
    return Token(token_location, LiteralNumber(token_location, token)) ## it returns a TOKEN!
end

# is it a KEYWORD of a IDENTIFIER?
function _parse_keyword_or_identifier_token(stream, first_char::Char, token_location::SourceLocation)
    token = first_char
    while true
        ch = read_char(stream)

        if (isletter(ch) || isdigit(ch) || ch == '_') == false
            unread_char(stream, ch)
            break
        end
        token *= ch
    end
    # If it is a keyword, it must be listed in the KEYWORDS enum:
    for i in 1:length(instances(KeywordEnum))
        token == string(KeywordEnum(i)) && return Token(token_location, Keyword(token_location, KeywordEnum(i)))
    end
    # If not it must be an IDENTIFIER:
    return Token(token_location, Identifier(token_location, token)) ## it returns a TOKEN!

end

"""
    read_token(stream::InputStream)

It reads all the stream, searching for the **TOKENS**!
"""
function read_token(stream::InputStream)
    if stream.saved_token != nothing
        result = stream.saved_token
        stream.saved_token = nothing
        return result
    end

    skip_whitespaces_and_comments(stream) # now there aren't spaces or comments!

    ch = read_char(stream)
    if ch == '€'
        return Stop(stream.location) # there's nothing in the file => Stop!
    end

    # Now it searches for any Token that star with the Char ch! (after saving its location)
    token_location = stream.location

    if occursin(ch, SYMBOLS)
        # It's a one-character symbol, like "(" or ","
        return Token(token_location, Symbol(token_location, ch))

    elseif ch == '"'
        # it's a literal string (used for FILE NAMES)
        return _parse_string_token(stream, token_location)

    elseif isdigit(ch) || ch in ['+', '-', '.']
        # It's a floating-point number
        return _parse_float_token(stream, ch, token_location)

    elseif isletter(ch) || ch == '_'
        # It can be both a KEYWORD or a IDENTIFIER, thus:
        return _parse_keyword_or_identifier_token( stream, ch, token_location )

    else
        # It's some weird character, like "@` or `&`
        throw(GrammarError("Invalid character $ch", stream.location))
    end

end

"""
    unread(stream::InputStream, token::Token)

Unread the token: *look-ahead* funtion (necessary in LL(1) grammars).
"""
function unread_token(stream::InputStream, token::Token)
    if stream.saved_token == nothing
        stream.saved_token = token
    end
end

## PARSER #########################################################################################

# auxiliary ASSERT functions:
assert_is_keyword(token::Token, keyword::KeywordEnum)= isa(token.value, Keyword) && token.value.keyword == keyword

assert_is_identifier(token::Token, identifier::Union{String, Char}) = isa(token.value, Identifier) && token.value.s == identifier

assert_is_symbol(token::Token, symbol::Union{String, Char}) = occursin(symbol, SYMBOLS) && isa(token.value, ProjectMoana.Symbol) && token.value.symbol == symbol

assert_is_number(token::Token, number) = isa(token.value, LiteralNumber) && token.value.number == number
#assert_is_number(token::Token, number::Float32) = isa(token.value, LiteralNumber) && token.value.number == number

assert_is_string(token::Token, string::Union{String, Char}) = isa(token.value, LiteralString) && token.value.s == string

"""
    Scene(materials, world, camera, float_variables, overridden_variables)

It creates a complete scene.

## Arguments:

- materials is a Dict{String, Material} that collects all the materials;
- world is World containing all the shapes;
- camera is a Union{Camera, Nothing} representing the observer;
- float_variables is a Dict{String, Float32} that collects all the floating point variables;
- overridden_variables is a Set{String}, for controlling animations.

"""
mutable struct Scene

    materials::Dict{String, Material}
    world::World
    camera::Union{Camera, Nothing}
    float_variables::Dict{String, Float32}
    overridden_variables::Set{String}

    Scene(materials::Dict{String, Material} = Dict{String, Material}(),
          world::World = World(),
          camera::Union{Camera, Nothing} = nothing,
          float_variables::Dict{String, Float32} = Dict{String, Float32}(),
          overridden_variables::Set{String} = Set{String}() ) =
          new(materials, world, camera,float_variables, overridden_variables)
end

## EXPECT functions

function expect_symbol(input_file::InputStream, symbol::Char)
    token = read_token(input_file)
    try
        assert_is_symbol(token, symbol)
    catch e
        throw(GrammarError("got $token instead of $symbol", token.loc))
    end
end

function expect_keywords(input_file::InputStream, keywords::Array{KeywordEnum})
    token = read_token(input_file)
    try
        isa(token.value, Keyword)
    catch e
        throw(GrammarError("expected a keyword instead of $token", token.loc))
    end
    try
        token.value.keyword in keywords
    catch
        throw(GrammarError("expected one of the keywords $keywords instead of $token", token.loc))
    end
    return token.value.keyword
end

function expect_number(input_file::InputStream, scene::Scene)
    token = read_token(input_file)
    if isa(token.value, LiteralNumber)
        return token.value.number
    elseif isa(token.value, Identifier)
        variable_name = token.value.s
        if haskey(scene.float_variables, variable_name) == false
            throw(GrammarError("unknown variable $token", token.loc))
        end
        return scene.float_variables[variable_name]
    end
    throw(GrammarError("got $token instead of a number", token.loc))
end

function expect_string(input_file::InputStream)

    token = read_token(input_file)

    try
        isa(token, LiteralString) == false
    catch e
        throw(GrammarError("got $token instead of a string", token.loc))
    end
    return token.value.s
end

function expect_identifier(input_file::InputStream)

    token = read_token(input_file)

    try
        isa(token, Identifier)
    catch e
        throw(GrammarError("got $token instead of an identifier", token.loc))
    end
    return token.value.s
end

## PARSING functions ########################################################################

function parse_vector(input_file::InputStream, scene::Scene)
    expect_symbol(input_file, '[')
    x = expect_number(input_file, scene)
    expect_symbol(input_file, ',')
    y = expect_number(input_file, scene)
    expect_symbol(input_file, ',')
    z = expect_number(input_file, scene)
    expect_symbol(input_file, ']')

    return Vec(x, y, z)
end

function parse_color(input_file::InputStream, scene::Scene)
    expect_symbol(input_file, '<')
    red = expect_number(input_file, scene)
    expect_symbol(input_file, ',')
    green = expect_number(input_file, scene)
    expect_symbol(input_file, ',')
    blue = expect_number(input_file, scene)
    expect_symbol(input_file, '>')

    return RGB(red, green, blue)
end

function parse_pigment(input_file::InputStream, scene::Scene)
    keyword = expect_keywords(input_file, [UNIFORM, CHECKERED, IMAGE])
    expect_symbol(input_file, '(')

    if keyword == UNIFORM
        color = parse_color(input_file, scene)
        result = UniformPigment(color)
    elseif keyword == CHECKERED
        color1 = parse_color(input_file, scene)
        expect_symbol(input_file, ',')
        color2 = parse_color(input_file, scene)
        expect_symbol(input_file, ',')
        num_of_steps = Int32(expect_number(input_file, scene))
        result = CheckeredPigment(color1, color2, num_of_steps)
    elseif keyword == IMAGE
        file_name = expect_string(input_file)
        image_file = open(file_name, "r")
        image = read_pfm_image(image_file)
        result = ImagePigment(image)
    else
        throw(GrammarError("This line should be unreachable"))
    end
    expect_symbol(input_file, ')')
    return result
end

function parse_brdf(input_file::InputStream, scene::Scene)
    brdf_keyword = expect_keywords(input_file, [DIFFUSE, SPECULAR])
    expect_symbol(input_file, '(')
    pigment = parse_pigment(input_file, scene)
    expect_symbol(input_file, ')')

    if brdf_keyword == DIFFUSE
        return DiffuseBRDF(pigment)
    elseif brdf_keyword == SPECULAR
        return SpecularBRDF(pigment)
    else
        throw(GrammarError("This line should be unreachable"))
    end
end

function parse_material(input_file::InputStream, scene::Scene)
    name = expect_identifier(input_file)
    expect_symbol(input_file, '(')
    brdf = parse_brdf(input_file, scene)
    expect_symbol(input_file, ',')
    emitted_radiance = parse_pigment(input_file, scene)
    next_token = read_token(input_file)
    if isa(next_token.value, Symbol) && next_token.value.symbol == ')'
        return name, Material(brdf, emitted_radiance, 1.f0)
    else
        unread_token(input_file, next_token)
        expect_symbol(input_file, ',')
        emitted_intensity = expect_number(input_file, scene)
        expect_symbol(input_file, ')')
        return name, Material(brdf, emitted_radiance, emitted_intensity)
    end
end

function parse_transformation(input_file::InputStream, scene::Scene)
    result = Transformation()

    while true
        transformation_kw = expect_keywords(input_file, [
            IDENTITY,
            TRANSLATION,
            ROTATION_X,
            ROTATION_Y,
            ROTATION_Z,
            SCALING,
        ])

        if transformation_kw == IDENTITY
            nothing  # Do nothing (this is a primitive form of optimization!)
        elseif transformation_kw == TRANSLATION
            expect_symbol(input_file, '(')
            result *= translation(parse_vector(input_file, scene))
            expect_symbol(input_file, ')')
        elseif transformation_kw == ROTATION_X
            expect_symbol(input_file, '(')
            result *= rotation_x(expect_number(input_file, scene) * pi / 180.)
            expect_symbol(input_file, ')')
        elseif transformation_kw == ROTATION_Y
            expect_symbol(input_file, '(')
            result *= rotation_y(expect_number(input_file, scene) * pi / 180.)
            expect_symbol(input_file, ')')
        elseif transformation_kw == ROTATION_Z
            expect_symbol(input_file, '(')
            result *= rotation_z(expect_number(input_file, scene) * pi / 180.)
            expect_symbol(input_file, ')')
        elseif transformation_kw == SCALING
            expect_symbol(input_file, '(')
            result *= scaling(parse_vector(input_file, scene))
            expect_symbol(input_file, ')')
        end

        # We must peek the next token to check if there is another transformation that is being
        # chained or if the sequence ends. Thus, this is a LL(1) parser.
        next_kw = read_token(input_file)
        if (isa(next_kw.value, Symbol) && (next_kw.value.symbol == '*')) == false
            # Pretend you never read this token and put it back!
            unread_token(input_file, next_kw)
            break
        end
    end
    return result
end

function parse_sphere(input_file::InputStream, scene::Scene)
    expect_symbol(input_file, '(')

    material_name = expect_identifier(input_file)
    if haskey(scene.materials, material_name) == false
        throw(GrammarError("unknown material $material_name", input_file.location))
    end
    expect_symbol(input_file, ',')
    transformation = parse_transformation(input_file, scene)
    expect_symbol(input_file, ')')

    return Sphere(transformation, scene.materials[material_name])
end

function parse_plane(input_file::InputStream, scene::Scene)
    expect_symbol(input_file, '(')

    material_name = expect_identifier(input_file)
    if haskey(scene.materials, material_name) == false
        throw(GrammarError("unknown material $material_name", input_file.location))
    end
    expect_symbol(input_file, ',')
    transformation = parse_transformation(input_file, scene)
    expect_symbol(input_file, ')')

    return Plane(transformation, scene.materials[material_name])
end

function parse_aab(input_file::InputStream, scene::Scene)
    expect_symbol(input_file, '(')

    material_name = expect_identifier(input_file)
    if haskey(scene.materials, material_name) == false
        # We raise the exception here because input_file is pointing to the end of the wrong identifier
        throw(GrammarError("unknown material $material_name", input_file.location))
    end
    expect_symbol(input_file, ',')
    transformation = parse_transformation(input_file, scene)
    expect_symbol(input_file, ')')

    return AAB(transformation, scene.materials[material_name])
end

function parse_lightpoint(input_file::InputStream, scene::Scene)
    expect_symbol(input_file, '(')

    material_name = expect_identifier(input_file)
    if haskey(scene.materials, material_name) == false
        # We raise the exception here because input_file is pointing to the end of the wrong identifier
        throw(GrammarError("unknown material $material_name", input_file.location))
    end
    expect_symbol(input_file, ',')
    vector = parse_vector(input_file, scene)
    expect_symbol(input_file, ')')

    return LightPoint(toPoint(vector), scene.materials[material_name])
end

function parse_camera(input_file::InputStream, scene::Scene, width, height)
    expect_symbol(input_file, '(')
    type_kw = expect_keywords(input_file, [PERSPECTIVE, ORTHOGONAL])
    expect_symbol(input_file, ',')
    if type_kw == PERSPECTIVE
        distance = expect_number(input_file, scene)
        expect_symbol(input_file, ',')
    end
    transformation = parse_transformation(input_file, scene)

    aspect_ratio = 0
    next_token = read_token(input_file)
    if isa(next_token.value, Symbol) && next_token.value.symbol == ','
        aspect_ratio = expect_number(input_file, scene)
    end
    if aspect_ratio == 0
        aspect_ratio = width / height
        unread_token(input_file, next_token)
    end
    expect_symbol(input_file, ')')

    if type_kw == PERSPECTIVE
        result = PerspectiveCamera(aspect_ratio, transformation, distance)
    elseif type_kw == ORTHOGONAL
        result = OrthogonalCamera(aspect_ratio, transformation)
    end
    return result
end

## PARSING the whole scene #################################################################################

"""
    parse_scene(input_file, variables)

It parse the whole scene, making it ready to be rendered.

- input_file is a txt file, containing the scene to be parsed;
- variables is a Dict{String, Float32} containing the floating point variables.

"""
function parse_scene(input_file::InputStream, variables::Dict{String, Float32} = Dict{String, Float32}() )

    scene = Scene()
    scene.float_variables = copy(variables)
    scene.overridden_variables = Set{String}()

    while true
        what = read_token(input_file)
        if isa(what, Stop)
            break
        end
        if isa(what.value, Keyword) == false
            throw(GrammarError("expected a keyword instead of $what", what.loc))
        end
        if what.value.keyword == FLOAT || what.value.keyword == INT
            variable_name = expect_identifier(input_file)

            # Save this for the error message
            variable_loc = input_file.location

            expect_symbol(input_file, '(')
            variable_value = expect_number(input_file, scene)
            expect_symbol(input_file, ')')

            if haskey(scene.float_variables, variable_name) == true && haskey(scene.overridden_variables, variable_name) == false
                throw(GrammarError("variable «$variable_name» cannot be refunctionined",variable_loc))
            end

            if (variable_name in scene.overridden_variables) == false
                scene.float_variables[variable_name] = variable_value
            end

        elseif what.value.keyword == SPHERE
            add_shape(scene.world,parse_sphere(input_file, scene))
        elseif what.value.keyword == PLANE
            add_shape(scene.world, parse_plane(input_file, scene))
        elseif what.value.keyword == AABOX
            add_shape(scene.world, parse_aab(input_file, scene))
        elseif what.value.keyword == LIGHTPOINT
            add_shape(scene.world, parse_lightpoint(input_file, scene))
        elseif what.value.keyword == CAMERA
            if scene.camera != nothing
                throw(GrammarError("You cannot functionine more than one camera", what.location))
            end
            width = haskey(scene.float_variables, "WIDTH") ? scene.float_variables["WIDTH"] : 640
            height = haskey(scene.float_variables, "HEIGHT") ? scene.float_variables["HEIGHT"] : 480
            scene.camera = parse_camera(input_file, scene, width, height)
        elseif what.value.keyword == MATERIAL
            name, material = parse_material(input_file, scene)
            scene.materials[name] = material
        end
    end

    return scene
end


text:
@testset "382.linked-list-random-node.jl" begin
    node = ListNode{Int}([1, 2, 2, 3])
    @test begin
        freq = Dict{Int,Int}(1 => 0, 2 => 0, 3 => 0)
        for _ in 1:10000
            freq[get_random(node)] += 1
        end
        num = freq[1] + freq[2] + freq[3]
        !(freq[1] == 2 * freq[2] == freq[3]) &&
            abs(freq[1] / num - freq[3] / num) + abs(freq[1] / num - freq[2] / num / 2) <= 0.1
    end
end


text:
module PValueAdjust

  export Bonferroni
  export BenjaminiHochberg
  export BenjaminiYekutieli
  export Hochberg
  export Holm
  export padjust

  include("types.jl")
  include("helper.jl")
  include("bonferroni.jl")
  include("benjaminiHochberg.jl")
  include("hochberg.jl")
  include("holm.jl")
  include("benjaminiYekutieli.jl")
  include("padjust.jl")

end # module


text:
suml2loss(yhat, y) = sum((yhat .- y).^2)
sumabsloss(yhat, y) = sum(abs.(yhat .- y))
l2loss(yhat, y) = mean((yhat .- y).^2)
absloss(yhat, y) = mean(abs.(yhat .- y))
mape(yhat, y, eps = 1.0) =
    100.0 * mean(abs.(yhat .- (y .+ eps)) ./ abs.(y .+ eps))


text:
using Test
using LogicCircuits
using LogicCircuits: Element # test some internals

@testset "Trimmed SDD test" begin

    num_vars = 7
    manager = SddMgr(num_vars, :balanced)

    @test num_variables(manager) == num_vars
    @test num_nodes(manager) == 2*num_vars-1
    @test num_edges(manager) == 2*num_vars-2
    @test manager isa SddMgr

    @test varsubset(left_most_descendent(manager), manager)
    @test varsubset(manager.left, manager)
    @test varsubset(manager.right, manager)
    @test varsubset_left(manager.left, manager)
    @test varsubset_left(manager.left.left, manager)
    @test varsubset_left(manager.left.right, manager)
    @test varsubset_right(manager.right, manager)
    @test varsubset_right(manager.right.right, manager)
    @test varsubset_right(manager.right.left, manager)

    @test !varsubset(manager, left_most_descendent(manager))
    @test !varsubset_left(manager.right, manager)
    @test !varsubset_left(manager.right.left, manager)
    @test !varsubset_left(manager.right.right, manager)
    @test !varsubset_left(manager, manager)
    @test !varsubset_left(manager, manager.left)
    @test !varsubset_left(manager, manager.right)
    @test !varsubset_right(manager.left, manager)
    @test !varsubset_right(manager.left.right, manager)
    @test !varsubset_right(manager.left.left, manager)
    @test !varsubset_right(manager, manager)
    @test !varsubset_right(manager, manager.left)
    @test !varsubset_right(manager, manager.right)

    x = Var(1)
    y = Var(2)

    x_c = compile(manager, var2lit(x))
    y_c = compile(manager, var2lit(y))

    @test x_c != y_c

    @test variable(x_c) == x
    @test literal(x_c) == var2lit(x)
    @test vtree(x_c) ∈ manager
    @test ispositive(x_c)
    @test x_c == compile(manager, var2lit(x))

    @test variable(y_c) == y
    @test literal(y_c) == var2lit(y)
    @test vtree(y_c) ∈ manager
    @test ispositive(y_c)
    @test y_c == compile(manager, var2lit(y))

    notx = -var2lit(x)

    notx_c = compile(manager,notx)

    @test sat_prob(x_c) == 1//2
    @test sat_prob(notx_c) == 1//2
    @test model_count(x_c,num_vars) == BigInt(2)^(num_vars-1)
    @test model_count(notx_c,num_vars) == BigInt(2)^(num_vars-1)
    @test !notx_c(true)
    @test notx_c(false)

    @test variable(notx_c) == x
    @test literal(notx_c) == notx
    @test vtree(notx_c) ∈ manager
    @test isnegative(notx_c)
    @test notx_c == compile(manager, notx)

    true_c = compile(manager,true)

    @test istrue(true_c)
    @test constant(true_c) == true

    false_c = compile(manager,false)

    @test isfalse(false_c)
    @test constant(false_c) == false

    @test !true_c == false_c
    @test !false_c == true_c
    @test !x_c == notx_c
    @test !notx_c == x_c

    @test model_count(true_c,num_vars) == BigInt(2)^(num_vars)
    @test model_count(false_c,num_vars) == BigInt(0)

    v1 = compile(manager, Lit(1))
    v2 = compile(manager, Lit(2))
    v3 = compile(manager, Lit(3))
    v4 = compile(manager, Lit(4))
    v5 = compile(manager, Lit(5))
    v6 = compile(manager, Lit(6))
    v7 = compile(manager, Lit(7))
    @test_throws Exception compile(manager, Lit(8))

    p1 = [Element(true_c,v3)]
    @test canonicalize(p1, manager.left.right) === v3
    p2 = [Element(v1,true_c), Element(!v1,false_c)]
    @test canonicalize(p2, manager.left) === v1

    p3 = [Element(v1,v4), Element(!v1,v7)]
    n1 = canonicalize(p3, manager)
    p4 = [Element(!v1,v7), Element(v1,v4)]
    n2 = canonicalize(p4,manager)
    @test n1.vtree.left === manager.left
    @test n1.vtree.right === manager.right
    @test n1 === n2
    @test isdeterministic(n1)
    @test n1(true, false, false, true, false, false, false)
    @test !n1(false, true, false, true, false, false, false)

end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
"""
    RMSE(obs,sim)

    Returns the Root Mean Squared Error between observations `obs` and simulations `sim`.

    The closer to 0 the better.
"""
function RMSE(obs, sim)
    return sqrt(sum((obs .- sim).^2) / length(obs))
end

"""
    NRMSE(obs,sim)

    Returns the Normalized Root Mean Squared Error between observations `obs` and simulations `sim`.
    Normalization is performed using division by observations range (max-min).

    Output: Float/Particles
"""
function NRMSE(obs, sim)
    return sqrt(sum((obs .- sim).^2) / length(obs)) / (findmax(obs)[1] - findmin(obs)[1])
end

"""
    EF(obs,sim)

    Returns the Efficiency Factor between observations `obs` and simulations `sim` using NSE (Nash-Sutcliffe efficiency) model.
    More information can be found at https://en.wikipedia.org/wiki/Nash%E2%80%93Sutcliffe_model_efficiency_coefficient.

    The closer to 1 the better.
"""
function EF(obs, sim)
    SSres = sum((obs - sim).^2)
    SStot = sum((obs .- mean(obs)).^2)
    return 1 - SSres / SStot
end

"""
    dr(obs,sim)

    Returns the Willmott’s refined index of agreement dᵣ.
    Willmot et al. 2011. A refined index of model performance. https://rmets.onlinelibrary.wiley.com/doi/10.1002/joc.2419

    The closer to 1 the better.

"""
function dr(obs, sim)
    a = sum(abs.(obs .- sim))
    b = 2 * sum(abs.(obs .- mean(obs)))
    return 0 + (1 - a / b) * (a <= b) + (b / a - 1) * (a > b)
end


text:
if !isdefined(:__EXPRESSION_HASHES__)
    __EXPRESSION_HASHES__ = Set{UInt64}()
end

macro runonce(expr)
    h = hash(expr)
    return esc(quote
        if !in($h, __EXPRESSION_HASHES__)
            push!(__EXPRESSION_HASHES__, $h)
            $expr
        end
    end)
end

"""Same as `get` function, but evaluates default_expr only if needed"""
macro get(dict, key, default_expr)
    return quote
        if haskey($dict, $key)
            $dict[$key]
        else
            $default_expr
        end
    end
end

"""
Same as `@get`, but creates new object from `default_expr` if
it didn't exist before
"""
macro get_or_create(dict, key, default_expr)
    return quote
        if !haskey($dict, $key)
            $dict[$key] = $default_expr
        end
        $dict[$key]
    end
end

"""
Same as `@get`, but immediately exits function and return `default_expr`
if key doesn't exist.
"""
macro get_or_return(dict, key, default_expr)
    return quote
        if haskey($dict, $key)
            $dict[$key]
        else
            return $default_expr
            nothing  # not reachable, but without it code won't compile
        end
    end
end

"""
Get array of size `sz` from a `dict` by `key`. If element doesn't exist or
its size is not equal to `sz`, create and return new array
using `default_expr`. If element exists, but is not an error,
throw ArgumentError.
"""
macro get_array(dict, key, sz, default_expr)
    return quote
        if (haskey($dict, $key) && !isa($dict[$key], Array))
            local k = $key
            throw(ArgumentError("Key `$k` exists, but is not an array"))
        end
        if (!haskey($dict, $key) || size($dict[$key]) != $sz)
            # ensure $default_expr results in an ordinary array
            $dict[$key] = convert(Array, $default_expr)
        end
        $dict[$key]
    end
end

function logistic(x)
    return 1 ./ (1 + exp(-x))
end

const KNOWN_OPTIONS =
    [:debug, :gradient, :update, :sampler, :scorer, :reporter,
     :batch_size, :n_epochs, :n_gibbs,
     :lr, :momentum, :weight_decay_kind, :weight_decay_rate,
     :sparsity_cost, :sparsity_target,
     :randomize,
     # deprecated options
     :n_iter]
const DEPRECATED_OPTIONS = Dict(:n_iter => :n_epochs)

function check_options(opts::Dict)
    deprecated_keys = keys(DEPRECATED_OPTIONS)
    debug = get(opts, :debug, true)

    for opt in keys(opts)
        if debug && !in(opt, KNOWN_OPTIONS)
            warn("Option '$opt' is unknownm ignoring")
        end
        if debug && in(opt, deprecated_keys)
            warn("Option '$opt' is deprecated, " *
                 "use '$(DEPRECATED_OPTIONS[opt])' instead")
        end
    end
end

function split_evenly(n, len)
    n_parts = Int(ceil(n / len))
    parts = Array(Tuple, n_parts)
    for i=1:n_parts
        start_idx = (i-1)*len + 1
        end_idx = min(i*len, n)
        parts[i] = (start_idx, end_idx)
    end
    return parts
end

"""
`tofinite!` takes an array and
1. turns all NaNs to zeros
2. turns all Infs and -Infs to the largets and
   smallest representable values accordingly.
3. turns all zeros to the smallest representable
   non-zero values, if `nozeros` is true
"""
function tofinite!(x::Array; nozeros=false)
    for i in eachindex(x)
        if isnan(x[i])
            x[i] = 0.0
        elseif isinf(x[i])
            if x[i] > 0.0
                x[i] = prevfloat(x[i])
            else
                x[i] = nextfloat(x[i])
            end
        end

        if x[i] == 0.0 && nozeros
            x[i] = nextfloat(x[i])
        end
    end
end

function ensure_type(newT::DataType, A::AbstractArray)
    if eltype(A) != newT
        map(newT, A)
    else
        A
    end
end

function add!{T}(X::Array{T}, inc::T)
    @simd for i=1:length(X)
        @inbounds X[i] += inc
    end
end


text:
# NOTE distance in potentials is ångströms when fitting.
# But bohrs when CP2K calculates them. So that is why there is
# 0.52917721090 -terms in few places.

"""
GeneralAngle <: AbstractClusterPotential

General potential with customizable powers for both radial polynomials and cosθ.

# Fields
- `constants::Array{Float64,2}` : constants
- `ppowers::Vector{Int}`        : powers for radial polynomials
- `cpowers::Vector{Int}`        : powers for cosθ

To create potential use:
```julia
GeneralAngle(ppowers::AbstractVector{Int}, cpowers::AbstractVector{Int})
```
"""
mutable struct GeneralAngle <: AbstractClusterPotential
    constants::Array{Float64,2}
    "Powers for radial polynomials"
    ppowers::Vector{Int}
    "Powers for cosθ"
    cpowers::Vector{Int}
    function GeneralAngle(ppowers::AbstractVector{Int}, cpowers::AbstractVector{Int})
        new(zeros(length(cpowers), length(ppowers)), ppowers, cpowers)
    end
end

function Base.show(io::IO, potential::GeneralAngle; energy_unit="cm^-1")
    print(io, "GeneralAngle constants: ",potential.constants)
end

function calculate_potential(potential::GeneralAngle,
                                     cluster::AbstractCluster, indices)
    r = distances(cluster, indices[1], indices[2])
    θ = cluster_angle(cluster, indices[1], indices[2], indices[3])
    return potential(r,θ)
end

function potential_variables(potential::GeneralAngle,
                            cluster::AbstractCluster, indices)
    # Convert from Å to bohr
    r = distances(cluster, indices[1], indices[2]) ./ 0.52917721090
    θ = cluster_angle(cluster, indices[1], indices[2], indices[3])
    rr = r.^potential.ppowers
    cθ = cos(θ).^potential.cpowers
    out = [  x .* cθ for x in rr ]
    return vcat(out...)'
end

function get_potential!(potential::GeneralAngle, constants...)
    @assert length(constants) == length(potential.constants)
    for i in eachindex(constants)
        potential.constants[i] = constants[i]
    end
end

function (p::GeneralAngle)(r::Number,α::Number)
    # input is Å calculation in bohr
    c = cos.(α)
    cp = c.^p.cpowers
    rp = (r./0.52917721090).^p.ppowers
    sum(p.constants .* (cp * rp'))
end


text:
"""
    box(data)

Creates an `EChart` as a boxplot, optionally with outliers.

## Methods
```julia
box(data::AbstractVector{<:Union{Missing, Real}})
box(data::AbstractVector{<:AbstractVector{<:Union{Missing, Real}}})
box(df::AbstractDataFrame, data::Symbol)
box(df::AbstractDataFrame, data::Symbol, group::Symbol)
```

## Arguments
* `names::Union{AbstractVector, Void} = nothing` : provide names for series
* `outliers::Bool` : show outliers?
* `legend::Bool` : display legend?
* `horizontal::Bool = false` : show horizontal? (default: vertical)
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

Reasonable defaults set for different methods of `box`, such as displaying a legend when two or more series present.

## Examples
```julia
data = [850, 740, 900, 1070, 930, 850, 950, 980, 980, 880, 1000, 980, 930, 650, 760, 810, 1000, 1000, 960, 960]
b = box(data)
```
"""
function box(data::AbstractVector{<:AbstractVector{<:Union{Missing, Real}}};
 names::Union{AbstractVector, Nothing} = nothing,
 outliers::Bool = true,
 legend::Bool = false,
 horizontal::Bool = false,
 kwargs...)

    names == nothing ? names = [string(x) for x in 1:length(data)] : names
    length(names) != length(data) ? error("Number of series names must equal number of data series.") : nothing

    dataprep = [boxplotstat(x) for x in data]
    seriesdata = [x.summary for x in dataprep]

    ec = newplot(kwargs, ec_charttype = "box")

    ec.xAxis = [Axis(_type = "category", data = names, boundaryGap = true)]
    ec.yAxis = [Axis(_type = "value", scale = true)]

    ec.series = [BoxPlotSeries(name = "boxplot", _type = "boxplot", data = seriesdata)]

    if outliers
        #Format outliers to ECharts nested array format
        outliers = []
        for (i, s_) in enumerate(dataprep)
            if length(s_.outliers) > 0
                for outlier in s_.outliers
                    push!(outliers, [names[i], outlier])
                end
            end
        end
        #hack to call this BoxPlotSeries, since the data fields are the same
        push!(ec.series, BoxPlotSeries(name = "outliers", _type = "scatter", data = outliers))
    end

    #Add legend if requested
    legend ? legend!(ec) : nothing

    horizontal ? flip!(ec) : nothing

    return ec

end

#Single series
box(data::AbstractVector{<:Union{Missing, Real}};
  names::String = "",
  outliers::Bool = true,
  legend::Bool = false,
  horizontal::Bool = false,
  kwargs...) =

             box([data], names = [names], outliers = outliers, legend = legend, horizontal = horizontal, kwargs...)

#df with group arguments
function box(df::AbstractDataFrame, data::Symbol, group::Symbol;
         names::Union{AbstractVector, Nothing} = nothing,
         outliers::Bool = true,
         legend::Bool = false,
         horizontal::Bool = false,
         kwargs...)

    #Create grouped df
    subdf = groupby(df, group)

    #Get number of groups
    numgroups = length(subdf)

    #Get names
    names == nothing ? names = [unique(subdf[x][!, group])[1] for x in 1:numgroups] : names = names

    #Convert to array of array for existing method
    arrayarray = [convert(Array, subdf[x][!, data]) for x in 1:numgroups]

    ec = box(arrayarray, names = names, outliers = outliers, legend = legend, horizontal = horizontal, kwargs...)

    xaxis!(ec, name = string(group))
    yaxis!(ec, name = string(data))

    return ec

end

#df
function box(df::AbstractDataFrame, data::Symbol;
         names::Union{AbstractVector, Nothing} = nothing,
         outliers::Bool = true,
         legend::Bool = false,
         horizontal::Bool = false,
         kwargs...)

         converteddf = convert(Array, df[data])
         ec = box([converteddf], names = [names], outliers = outliers, legend = legend, horizontal = horizontal, kwargs...)

         yaxis!(ec, name = string(data))

         return ec

end


text:
# Encoding and decoding

# Code128 encoding with subtype Code128A
function _encode_code128a(msg::AbstractString)
    @inbounds for i = 1:ncodeunits(msg)
        codeunit(msg, i) ≥ 0x7f && throw(
            ArgumentError(
                "The given `msg` contains characters outside the range 0 - 126 and " *
                "cannot be enncoded in Code128"
            ),
        )
        codeunit(msg, i) ≥ 0x60 && throw(
            ArgumentError(
                "The given ascii `msg` contains lowercase letters or characters outside " *
                "the range 0 - 95 and cannot be encoded in subtype Code128A"
            )
        )
    end

    code = ["START A"]
    @inbounds for i = 1:ncodeunits(msg)
        if 0x20 ≤ codeunit(msg, i) ≤ 0x5f
            push!(code, string(msg[i]))
        elseif codeunit(msg, i) ≤ 0x1f
            push!(code, CODE128[codeunit(msg, i) + 0x41, :code128a])
        else
            throw(
                ArgumentError(
                    "Char in position $i in `msg` cannot be encoded in subtype Code128A"
                )
            )
        end
    end
    append!(code, ["CHECKSUM", "STOP"])

    return code
end

# Code128 encoding with subtype Code128B
function _encode_code128b(msg::AbstractString)
    @inbounds for i = 1:ncodeunits(msg)
        codeunit(msg, i) ≥ 0x7f && throw(
            ArgumentError(
                "The given `msg` contains characters outside the range 0 - 126 and " *
                "cannot be enncoded in Code128"
            )
        )
        codeunit(msg, i) ≤ 0x1f && throw(
            ArgumentError(
                "The given ascii `msg` contains symbology characters outside the range " *
                "32 - 127 and cannot be fully encoded in subtype Code128B"
            )
        )
    end

    code = [
        "START B"
        string.(collect(msg))
        "CHECKSUM"
        "STOP"
    ]
    return code
end

# Code128 encoding with subtype Code128C
function _encode_code128c(msg::AbstractString)
    all(isdigit, msg) || throw(
        ArgumentError(
            "The given `msg` contains characters which are not digits and cannot " *
            "be encoded in subtype Code128C"
        )
    )
    iseven(length(msg)) || throw(
        ArgumentError(
            "The given `msg` contains an odd number of digits and cannot be fully " *
            "encoded in subtype Code128C"
        )
    )

    code = [
        "START C"
        [msg[j:j+1] for j = 1:2:length(msg)]
        "CHECKSUM"
        "STOP"
    ]
    return code
end

# Optimized Code128 (actually GS1-128) mixed-subtype encoding following the rules in
# "GS1 General Specifications, Version 13, Issue 1, Jan-2013, Section 5.4.7.7.
# Use of Start, Code Set, and Shift symbols to Minimize Symbol Length (Informative),
# pages 268 to 269."
function _encode_code128(msg)
    isascii(msg) || throw(
        ArgumentError(
            "The given `msg` contains non-ascii characters and cannot " *
            "be fully encoded in Code128"
        )
    )
    code = String[]

    len_msg = length(msg)
    are_digits = [isdigit(d) for d in msg]

    # Determine start character:
    if (length(are_digits) == 2 && are_digits == [1, 1]) || (
        length(are_digits) ≥ 4 && are_digits[1:4] == [1, 1, 1, 1]
    ) # Determine whether there are enough digits
        subtype = :code128c
    elseif UInt8(first(msg)) ≤ 0x1f # check whether it is a symbology element (NUL to US)
        subtype = :code128a
    else
        subtype = :code128b
    end

    push!(code, "START $(uppercase(string(subtype)[end]))")
    ind = 1
    nextsubtype = subtype

    while ind ≤ length(msg)
        # Determine whether it needs to change or shift subtype
        if nextsubtype == :code128c && (len_msg == ind || (len_msg > ind && are_digits[ind:ind+1] != [1, 1]))
            if codeunit(msg, ind) ≤ 0x1f ||
                    (codeunit(msg, ind) ≤ 0x5f && len_msg > ind &&
                        0 ≤ codeunit(msg, ind+1) ≤ 0x1f)
                push!(code, "CODE A")
                nextsubtype = subtype = :code128a
            else
                push!(code, "CODE B")
                nextsubtype = subtype = :code128b
            end
        elseif nextsubtype != :code128c && are_digits[ind] && (
                (
                    findnext(iszero, are_digits, ind) !== nothing &&
                    (0, 2) ≤ reverse(divrem(findnext(iszero, are_digits, ind) - ind, 2)) < (1, 0)
                    # iseven(findnext(iszero, are_digits, ind) - ind)
                ) ||

                (
                    findnext(iszero, are_digits, ind) === nothing &&
                    (0, 2) ≤ reverse(divrem(len_msg - ind + 1, 2)) < (1, 0)
                    # len_msg > ind && isodd(len_msg - ind)
                )
            )
            push!(code, "CODE C")
            nextsubtype = subtype = :code128c
        elseif nextsubtype != :code128a && codeunit(msg, ind) ≤ 0x1f # check for symbology
            if nextsubtype == :code128b &&
                    len_msg > ind && 0x60 ≤ codeunit(msg, ind+1) ≤ 0x7e
                push!(code, "SHIFT A")
                nextsubtype = :code128a
            else
                push!(code, "CODE A")
                subtype = nextsubtype = :code128a
            end
        elseif nextsubtype != :code128b && 96 ≤ Int(msg[ind]) ≤ 126 # check for lowercase
            if nextsubtype == :code128a && len_msg > ind && 0 ≤ Int(msg[ind+1]) ≤ 31
                push!(code, "SHIFT B")
                nextsubtype = :code128b
            else
                push!(code, "CODE B")
                subtype = nextsubtype = :code128b
            end
        end

        # encode msg chunk
        if nextsubtype == :code128c
            push!(code, string(msg[ind:ind+1]))
            ind += 2
        elseif nextsubtype == :code128a
            if codeunit(msg, ind) ≤ 0x1f
                push!(code, CODE128[codeunit(msg, ind) + 65, :code128a])
            else
                push!(code, string(msg[ind]))
            end
            ind += 1
        else
            if codeunit(msg, ind) == 0x7f
                push!(code, "DEL")
            else
                push!(code, string(msg[ind]))
            end
            ind += 1
        end
        nextsubtype = subtype
    end
    append!(code, ["CHECKSUM", "STOP"])
    return code
end

"""
    encode(msg::AbstractString, encoding_type::Symbol)

Return the encoded sequence from the given `msg`, following the specifications determined
by the `encoding_type`.

Currently, only Code128 specification is available.

If `encoding_type` is either `:code128a`, `:code128b`, or `:code128c`, it returns the
encoding following the corresponding subtype. If `encoding_type` is `:code128`, it will
return an optimized encoding, possibily mixing different subtypes. This strategy
follows the specifications in "GS1 General Specifications, Version 13, Issue 1, Jan-2013,
Section 5.4.7.7. Use of Start, Code Set, and Shift symbols to Minimize Symbol Length
(Informative), pages 268 to 269."

The `msg` needs to be a string of ascii characteres to be encoded, otherwise the method
throws an `ArgumentError`.

# Examples

```jldoctest
julia> code = encode("000132", :code128c)
6-element Vector{String}:
 "START C"
 "00"
 "01"
 "32"
 "CHECKSUM"
 "STOP"

julia> code = encode("ABC", :code128a)
6-element Vector{String}:
 "START A"
 "A"
 "B"
 "C"
 "CHECKSUM"
 "STOP"

julia> code = encode("AaBC\x02", :code128)
9-element Vector{String}:
 "START B"
 "A"
 "a"
 "B"
 "C"
 "CODE A"
 "\x02"
 "CHECKSUM"
 "STOP"
```
"""
function encode(msg::AbstractString, encoding_type::Symbol)
    if encoding_type == :code128
        return _encode_code128(msg)
    elseif encoding_type == :code128a
        return _encode_code128a(msg)
    elseif encoding_type == :code128b
        return _encode_code128b(msg)
    elseif encoding_type == :code128c
        return _encode_code128c(msg)
    else
        throw(
            ArgumentError(
                "Encoding type `$(Meta.quot(mode))` not implemented"
            )
        )
    end
end

function _decode_code128(code::Vector{String})
    m = match(r"^START (A|B|C)$", first(code))
    m !== nothing || throw(
        ArgumentError(
            "First element of `code` should be either `START A`, `START B` or `START C`"
        ),
    )
    subtype = Symbol("code128$(lowercase(m.captures[1]))")
    nextsubtype = subtype
    msg = ""

    for c in code
        if startswith(c, r"^CODE [A|B|C]$")
            nextsubtype = subtype = Symbol("code128$(lowercase(c[end]))")
        elseif startswith(c, r"^SHIFT [A|B|C]$")
            nextsubtype = Symbol("code128$(lowercase(c[end]))")
        elseif nextsubtype == :code128a && !startswith(c, "CHECKSUM")
            val = CODE128.value[CODE128.code128a .== c][1]
            if val ≤ 0x3f # check if ≤ 63
                msg *= Char(val + 0x20) # add 32
            elseif val ≤ 0x5f # check if ≤ 95
                msg *= Char(val - 0x40) # subtract 64
            end
        elseif nextsubtype == :code128b && !startswith(c, "CHECKSUM")
            if CODE128.value[CODE128.code128b .== c][1] ≤ 94
                msg *= c
            elseif c == CODE128.code128b[96]
                msg *= "\x7f" # DEL
            end
        elseif nextsubtype == :code128c && all(isdigit, c)
            msg *= c
        end
    end
    return msg
end

"""
    decode(code::Vector{String}, encoding_type::Symbol)

Decode the gived encoded sequence, following the specifications determined by the
`encoding_type`.

Currently, only Code128 specification is available.

# Examples

```jldoctest
ulia> code = encode("000132", :code128)
6-element Vector{String}:
 "START C"
 "00"
 "01"
 "32"
 "CHECKSUM"
 "STOP"

julia> msg = decode(code, :code128)
"000132"

julia> code = encode("\x02ABC\x03", :code128)
8-element Vector{String}:
 "START A"
 "STX"
 "A"
 "B"
 "C"
 "ETX"
 "CHECKSUM"
 "STOP"

julia> msg = decode(code, :code128)
"\x02ABC\x03"

julia> code = encode("\x02Abc\x03", :code128)
10-element Vector{String}:
 "START A"
 "STX"
 "A"
 "CODE B"
 "b"
 "c"
 "CODE A"
 "ETX"
 "CHECKSUM"
 "STOP"
"""
function decode(code::Vector{String}, encoding_type::Symbol)
    if encoding_type in (:code128, :code128a, :code128b, :code128c)
        decoded = _decode_code128(code::Vector{String})
    else
        throw(
            ArgumentError(
                "Decoding type `$(Meta.quot(mode))` not implemented"
            )
        )
    end
    return decoded
end


text:
module RocketProxyObservableTest

using Test
using Rocket

@testset "ProxyObservable" begin

    println("Testing: proxy")

    struct DummyType end

    @testset begin
        @test_throws ErrorException proxy(Int, from(1:5), DummyType())
        @test_throws ErrorException proxy(Int, DummyType(), void(Int))
        @test_throws ErrorException proxy(Int, DummyType(), DummyType())
    end

    struct MyActorProxy <: ActorProxy end

    struct MyActor{L, A} <: Actor{L}
        actor :: A
    end

    Rocket.actor_proxy!(::Type{L}, proxy::MyActorProxy, actor::A) where { L, A } = MyActor{L, A}(actor)

    Rocket.on_next!(actor::MyActor{L}, data::L) where L = next!(actor.actor, data + 1)
    Rocket.on_error!(actor::MyActor, err)               = error!(actor.actor, err)
    Rocket.on_complete!(actor::MyActor)                 = complete!(actor.actor)

    @testset begin
        source = proxy(Int, from(1:5), MyActorProxy())
        actor  = keep(Int)

        subscribe!(source, actor)

        @test actor.values == [ 2, 3, 4, 5, 6 ]
    end

    struct MySourceProxy <: SourceProxy end

    struct MySource{L, S} <: Subscribable{L}
        source :: S
    end

    Rocket.source_proxy!(::Type{L}, proxy::MySourceProxy, source::S) where { L, S } = MySource{L, S}(source)

    function Rocket.on_subscribe!(source::MySource, actor)
        next!(actor, 0)
        return subscribe!(source.source, actor)
    end

    @testset begin
        source = proxy(Int, from(1:5), MySourceProxy())
        actor  = keep(Int)

        subscribe!(source, actor)

        @test actor.values == [ 0, 1, 2, 3, 4, 5 ]
    end

    struct MyActorSourceProxy <: ActorSourceProxy end

    struct MyActor2{L, A} <: Actor{L}
        actor :: A
    end

    Rocket.actor_proxy!(::Type{L}, proxy::MyActorSourceProxy, actor::A) where { L, A } = MyActor2{L, A}(actor)

    Rocket.on_next!(actor::MyActor2{L}, data::L) where L = next!(actor.actor, data + 2)
    Rocket.on_error!(actor::MyActor2, err)               = error!(actor.actor, err)
    Rocket.on_complete!(actor::MyActor2)                 = complete!(actor.actor)

    struct MySource2{L, S} <: Subscribable{L}
        source :: S
    end

    Rocket.source_proxy!(::Type{L}, proxy::MyActorSourceProxy, source::S) where { L, S } = MySource2{L, S}(source)

    function Rocket.on_subscribe!(source::MySource2, actor)
        next!(actor, 1)
        return subscribe!(source.source, actor)
    end

    @testset begin
        source = proxy(Int, from(1:5), MyActorSourceProxy())
        actor  = keep(Int)

        subscribe!(source, actor)

        @test actor.values == [ 3, 3, 4, 5, 6, 7 ]
    end

    @testset begin
        @test_throws Exception actor_proxy!(Any, DummyType(), DummyType())
        @test_throws Exception source_proxy!(Any, DummyType(), DummyType())

        @test_throws Exception Rocket.call_actor_proxy!(Any, DummyType(), DummyType())
        @test_throws Exception Rocket.call_actor_proxy!(Any, DummyType(), void(Any))
        @test_throws Exception Rocket.call_actor_proxy!(Any, MyActorProxy(), DummyType())
        @test_throws Exception Rocket.call_actor_proxy!(Any, MySourceProxy(), DummyType())
        @test_throws Exception Rocket.call_actor_proxy!(Any, MyActorSourceProxy(), DummyType())

        @test_throws Exception Rocket.call_source_proxy!(Any, DummyType(), DummyType())
        @test_throws Exception Rocket.call_source_proxy!(Any, DummyType(), never())
        @test_throws Exception Rocket.call_source_proxy!(Any, MyActorProxy(), DummyType())
        @test_throws Exception Rocket.call_source_proxy!(Any, MySourceProxy(), DummyType())
        @test_throws Exception Rocket.call_source_proxy!(Any, MyActorSourceProxy(), DummyType())
    end

    @testset begin
        source = proxy(Int, from(1:5), MySourceProxy())
        io = IOBuffer()

        show(io, source)

        printed = String(take!(io))

        @test occursin("ProxyObservable", printed)
        @test occursin(string(eltype(source)), printed)
    end

end

end


text:
using Documenter, Mortgage

makedocs(
	doctest = true,
	clean = true,
)

deploydocs(
    deps  = Deps.pip("mkdocs", "python-markdown-math"),
    repo  = "github.com/holder66/Mortgage.jl.git",
    julia = "release",
)


text:
using RetroSignalModel
using RetroSignalModel: rtgM4
using Documenter

makedocs(;
    modules=[RetroSignalModel],
    authors="stevengogogo <stevengogogo4321@gmail.com> and contributors",
    repo="https://github.com/stevengogogo/RetroSignalModel.jl/blob/{commit}{path}#L{line}",
    sitename="RetroSignalModel.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://stevengogogo.github.io/RetroSignalModel.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    branch = "gh-pages",
    repo="github.com/stevengogogo/RetroSignalModel.jl.git",
)


text:
using DataFrames: DataFrame
using FileIO

"""
    readMDTable(file::String; header::Bool=true)::DataFrame

Reads into a DataFrame any MD table file.\n
It's assumed the file will have a header, but can be specified with header=false if there isn't one.
"""
function readMDTable(file::String; header::Bool=true)::DataFrame
    mdTable::NamedTuple = NamedTuple()
    skip::Bool = header
    parsedHeader::Bool = false
    for line in eachline(file)
        if header
            mdTable = parseHeader(line)
            header = false
            parsedHeader = true
        elseif skip
            skip = false
        else
            if !parsedHeader
                mdTable = parseHeader(line, isHeader=false)
                parsedHeader = true
            end
            parseLine!(line, mdTable)
        end
    end
    return DataFrame(mdTable)
end

function parseHeader(line::String; isHeader::Bool=true)
    headersDirty = split(line, "|", keepempty=false)
    values::Array = []
    headersClean::Array = []
    count = 1
    for value in headersDirty
        key = lstrip(rstrip(value))
        if isHeader && key != ""
            push!(values, [])
            push!(headersClean, key)
        else
            push!(values, [])
            push!(headersClean, "x" * string(count))
            count += 1
        end
    end
    return NamedTuple{(Symbol.(headersClean)...,)}(values)
end

function parseLine!(line::String, mdTable::NamedTuple)
    values = split(line, "|", keepempty=false)
    if length(values) != length(mdTable)
        error("The number of values in a line must be the same as the number of headers of the file")
    end

    for i in 1:length(values)
        value = lstrip(rstrip(values[i]))
        push!(getindex(mdTable, i), value)
    end
end

"""
    load(f::File{format"MD"}; header::Bool=true)

Reads into a DataFrame any MD table file.
It's assumed the file will have a header, but can be specified with header=false if there isn't one.
"""
function load(f::File{format"MD"}; header::Bool=true)
    return readMDTable(f.filename, header=header)
end


text:
function acceleration(obj::AbstractParametricAcceleration, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "acceleration", FieldVector3D, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function acceleration(obj::AbstractParametricAcceleration, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "acceleration", Vector3D, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function get_events_detectors(obj::AbstractParametricAcceleration)
    return jcall(obj, "getEventsDetectors", Stream, ())
end

function get_field_events_detectors(obj::AbstractParametricAcceleration, arg0::Field)
    return jcall(obj, "getFieldEventsDetectors", Stream, (Field,), arg0)
end


text:
using Distributions

module Econometrics

## list packages whos namespace is used
using TimeData
using DataFrames
using Docile
using Dates
using GLM
using EconDatasets
using Distributions
using Gadfly
## using Winston
using NLopt
using JuMP

## required for testing
## using MAT

export                                  # important functions
AbstrModel,
AbstrMultivarModel,
AbstrUnivarModel,
CIR,
GARCH_1_1,
GARCH_1_1_Fit,
NChiSq,
NormIID,
TLSDist,
TlsIID,
bsDs,
bsCall,
bsPut,
bsDeltaCall,
bsDeltaPut,
bsGamma,
bsVega,
bsThetaCall,
bsThetaPut,
bsRhoCall,
bsRhoPut,
cirOls,
cirNllh,
cirNllhx,
disc2log,
estimate,
fit,
getParams,
implVola,
implVolaCall,
implVolaPut,
imputePreviousObs!,
localAppl,
log2disc,
plotLocalProperties,
price2ret,
ranks,
ret2price
## ishighest,
## islowest,
## ranks,
## plot

## distributions
include("distributions/t_loc_scale.jl")
include("distributions/nchisq.jl")

## models
include("models/Model.jl")
include("models/normiid.jl")
include("models/tlsiid.jl")
include("models/garch_types.jl")
include("models/garch.jl")
include("models/cir.jl")

## miscellaneous
include("bsOptions.jl")
include("localProperties.jl")
include("copula.jl")
include("returns.jl")
include("utils.jl")
include("autocorr.jl")

end # module


text:
function compute_bulkevolved_t!(bulkevol_t::BulkEvolved,
                                bulkconstrain::BulkConstrained, gauge_t::Gauge,
                                bulkevol::BulkEvolved, boundary::Boundary,
                                gauge::Gauge, sys::System, ::EvolTest0)

    B1_t, B2_t, G_t, phi_t = unpack(bulkevol_t)
    # B1  , B2  , G  , phi   = unpack(bulkevol)

    fill!(B1_t,  0)
    fill!(B2_t,  0)
    fill!(G_t,   0)
    fill!(phi_t, 0)

    nothing
end

function compute_bulkevolved_t!(bulkevol_t::BulkEvolved,
                                bulkconstrain::BulkConstrained, gauge_t::Gauge,
                                bulkevol::BulkEvolved, boundary::Boundary,
                                gauge::Gauge, sys::System{Inner}, evoleq::AffineNull)
    uu  = sys.ucoord
    Du  = sys.Du
    Dx  = sys.Dx
    Dy  = sys.Dy

    phi0  = evoleq.phi0
    phi02 = phi0 * phi0
    phi03 = phi0 * phi02

    Nu, Nx, Ny = size(sys)

    B1_t, B2_t, G_t, phi_t = unpack(bulkevol_t)

    # u = 0
    @fastmath @inbounds for j in 1:Ny
        @inbounds @simd for i in 1:Nx
            xi     = gauge.xi[1,i,j]

            B1     = bulkevol.B1[1,i,j]
            B2     = bulkevol.B2[1,i,j]
            G      = bulkevol.G[1,i,j]

            B1_u   = Du(bulkevol.B1, 1,i,j)
            B2_u   = Du(bulkevol.B2, 1,i,j)
            G_u    = Du(bulkevol.G,  1,i,j)

            B1d_u  = Du(bulkconstrain.B1d, 1,i,j)
            B2d_u  = Du(bulkconstrain.B2d, 1,i,j)
            Gd_u   = Du(bulkconstrain.Gd,  1,i,j)

            B1_t[1,i,j]  = B1d_u  + 2.5 * B1_u  + 4 * B1  * xi
            B2_t[1,i,j]  = B2d_u  + 2.5 * B2_u  + 4 * B2  * xi
            G_t[1,i,j]   = Gd_u   + 2.5 * G_u   + 4 * G   * xi
        end
    end

    # remaining inner grid points
    @fastmath @inbounds for j in 1:Ny
        @inbounds for i in 1:Nx
            xi    = gauge.xi[1,i,j]
            xi_t  = gauge_t.xi[1,i,j]
            @inbounds @simd for a in 2:Nu
                u      = uu[a]
                u4     = u * u * u * u

                B1     = bulkevol.B1[a,i,j]
                B2     = bulkevol.B2[a,i,j]
                G      = bulkevol.G[a,i,j]

                B1d    = bulkconstrain.B1d[a,i,j]
                B2d    = bulkconstrain.B2d[a,i,j]
                Gd     = bulkconstrain.Gd[a,i,j]
                A      = bulkconstrain.A[a,i,j]

                B1_u   = Du(bulkevol.B1, a,i,j)
                B2_u   = Du(bulkevol.B2, a,i,j)
                G_u    = Du(bulkevol.G,  a,i,j)

		B1_t[a,i,j] = ((u * B1_u + 4 * B1) *
                               (-2 * u * u * xi_t + A * u4 +
                                (xi * u + 1) * (xi * u + 1)) +
                               2 * B1d) / (2 * u) -
                               phi02 * u * (4 * B1 + B1_u * u) / 3

		B2_t[a,i,j] = ((u * B2_u + 4 * B2) *
                               (-2 * u * u * xi_t + A * u4 +
                                (xi * u + 1) * (xi * u + 1)) +
                               2 * B2d) / (2 * u) -
                               phi02 * u * (4 * B2 + B2_u * u) / 3

		G_t[a,i,j] = ((u * G_u + 4 * G) *
                               (-2 * u * u * xi_t + A * u4 +
                                (xi * u + 1) * (xi * u + 1)) +
                               2 * Gd) / (2 * u) -
                               phi02 * u * (4 * G + G_u * u) / 3
            end
        end
    end

    # if phi0 = 0 set phi_t to zero and return
    if abs(phi0) < 1e-9
        fill!(phi_t, 0)
        return
    end

    # otherwise, compute phi_t

    # u = 0
    @fastmath @inbounds for j in 1:Ny
        @inbounds @simd for i in 1:Nx
            xi     = gauge.xi[1,i,j]
            xi3    = xi * xi * xi
            xi_t   = gauge_t.xi[1,i,j]

            phi    = bulkevol.phi[1,i,j]

            phi_u  = Du(bulkevol.phi,1,i,j)

            phid_u = Du(bulkconstrain.phid,1,i,j)

            phi_t[1,i,j] = -xi3 / phi02 + 3 * xi * phi +  2//3 * xi +
                2 * xi * xi_t / phi02 + phid_u + 2 * phi_u
        end
    end

    # remaining inner grid points
    @fastmath @inbounds for j in 1:Ny
        @inbounds for i in 1:Nx
            xi     = gauge.xi[1,i,j]
            xi2    = xi * xi
            xi3    = xi * xi * xi
            xi_t   = gauge_t.xi[1,i,j]
            @inbounds @simd for a in 2:Nu
                u      = uu[a]
                u2     = u * u
                u4     = u2 * u2

                phi    = bulkevol.phi[a,i,j]
                phid   = bulkconstrain.phid[a,i,j]
                A      = bulkconstrain.A[a,i,j]

                phi_u  = Du(bulkevol.phi,a,i,j)

                phi_t[a,i,j] = (
                    - 6 * xi3 * u
                    + 3 * xi2 * (-3 + phi02 * u2 * (3 * phi + phi_u * u))
		    + 3 * A * u2 * (1 - 2 * xi * u
				    + phi02 * u2 * (3 * phi + phi_u * u))
		    + 2 * xi * u * (6 * xi_t + 9 * phi * phi02
				    + phi02 * (2 + 3 * phi_u * u))
		    + phi02 * (-2 + 6 * phid - (3 * phi + phi_u * u) *
                               (-3 + 6 * xi_t * u2 + 2 * phi02 * u2))
                ) / (6 * phi02 * u)
            end
        end
    end

    nothing
end

function compute_bulkevolved_t!(bulkevol_t::BulkEvolved,
                                bulkconstrain::BulkConstrained, gauge_t::Gauge,
                                bulkevol::BulkEvolved, boundary::Boundary,
                                gauge::Gauge, sys::System{Outer}, evoleq::AffineNull)
    uu  = sys.ucoord
    Du  = sys.Du
    Dx  = sys.Dx
    Dy  = sys.Dy

    phi0  = evoleq.phi0
    phi02 = phi0 * phi0
    phi03 = phi0 * phi02

    Nu, Nx, Ny = size(sys)

    B1_t, B2_t, G_t, phi_t = unpack(bulkevol_t)

    @fastmath @inbounds for j in 1:Ny
        @inbounds for i in 1:Nx
            xi_t  = gauge_t.xi[1,i,j]
            @inbounds @simd for a in 1:Nu
                u      = uu[a]
                u2     = u * u

                B1d    = bulkconstrain.B1d[a,i,j]
                B2d    = bulkconstrain.B2d[a,i,j]
                Gd     = bulkconstrain.Gd[a,i,j]
                A      = bulkconstrain.A[a,i,j]

                B1_u   = Du(bulkevol.B1, a,i,j)
                B2_u   = Du(bulkevol.B2, a,i,j)
                G_u    = Du(bulkevol.G,  a,i,j)

		B1_t[a,i,j] = B1d + u2 * (A/2 - xi_t) * B1_u
		B2_t[a,i,j] = B2d + u2 * (A/2 - xi_t) * B2_u
		G_t[a,i,j]  = Gd  + u2 * (A/2 - xi_t) * G_u
            end
        end
    end

    # if phi0 = 0 set phi_t to zero and return
    if abs(phi0) < 1e-9
        fill!(phi_t, 0)
        return
    end

    # otherwise, compute phi_t

    @fastmath @inbounds for j in 1:Ny
        @inbounds for i in 1:Nx
            xi_t   = gauge_t.xi[1,i,j]
            @inbounds @simd for a in 1:Nu
                u      = uu[a]
                u2     = u * u

                phid   = bulkconstrain.phid[a,i,j]
                A      = bulkconstrain.A[a,i,j]

                phi_u  = Du(bulkevol.phi, a,i,j)

		phi_t[a,i,j] = phid + u2 * (A/2 - xi_t) * phi_u
            end
        end
    end

    nothing
end

function sync_bulkevolved!(bulkevols_t, bulkconstrains, gauge_t::Gauge,
                           systems::SystemPartition, evoleq::AffineNull)
    Nsys = length(systems)

    @inbounds for i in 1:Nsys-1
        sync_bulkevolved!(bulkevols_t[i], bulkevols_t[i+1], bulkconstrains[i+1],
                          gauge_t, systems[i], systems[i+1], evoleq)
    end

    nothing
end

function sync_bulkevolved!(bulkevol1_t::BulkEvolved, bulkevol2_t::BulkEvolved,
                           bulkconstrain2::BulkConstrained, gauge_t::Gauge,
                           sys1::System{Outer}, sys2::System{Outer}, evoleq::AffineNull)
    u0 = sys2.ucoord[1]

    _, Nx, Ny = size(sys2)

    B11_t, B21_t, G1_t, phi1_t = unpack(bulkevol1_t)
    B12_t, B22_t, G2_t, phi2_t = unpack(bulkevol2_t)

    @fastmath @inbounds for j in 1:Ny
        @inbounds for i in 1:Nx
            xi_t   = gauge_t.xi[1,i,j]
            A      = bulkconstrain2.A[1,i,j]

            # characteristic speed
            c = u0 * u0 * (A/2 - xi_t)

            # if c > 0, mode is entering grid1 from grid2;
            # if c < 0, mode is entering grid2 from grid1.
            # we assume here that grids merely touch at the interface
            if c > 0
                B11_t[end,i,j]  = B12_t[1,i,j]
                B21_t[end,i,j]  = B22_t[1,i,j]
                G1_t[end,i,j]   = G2_t[1,i,j]
                phi1_t[end,i,j] = phi2_t[1,i,j]
            elseif c < 0
                B12_t[1,i,j]  = B11_t[end,i,j]
                B22_t[1,i,j]  = B21_t[end,i,j]
                G2_t[1,i,j]   = G1_t[end,i,j]
                phi2_t[1,i,j] = phi1_t[end,i,j]
            end
        end
    end

    nothing
end

function sync_bulkevolved!(bulkevol1_t::BulkEvolved, bulkevol2_t::BulkEvolved,
                           bulkconstrain2::BulkConstrained, gauge_t::Gauge,
                           sys1::System{Inner}, sys2::System{Outer}, evoleq::AffineNull)
    u0  = sys2.ucoord[1]
    u02 = u0 * u0
    u03 = u0 * u02
    u04 = u02 * u02

    phi0  = evoleq.phi0
    phi02 = phi0 * phi0
    phi03 = phi0 * phi02

    _, Nx, Ny = size(sys2)

    B11_t, B21_t, G1_t, phi1_t = unpack(bulkevol1_t)
    B12_t, B22_t, G2_t, phi2_t = unpack(bulkevol2_t)

    @fastmath @inbounds for j in 1:Ny
        @inbounds for i in 1:Nx
            xi_t   = gauge_t.xi[1,i,j]
            A      = bulkconstrain2.A[1,i,j]

            # characteristic speed
            c = u0 * u0 * (A/2 - xi_t)

            # if c > 0, mode is entering grid1 from grid2;
            # if c < 0, mode is entering grid2 from grid1.
            # we assume here that grids merely touch at the interface
            if c > 0
                B11_t[end,i,j]  = B12_t[1,i,j] / u04
                B21_t[end,i,j]  = B22_t[1,i,j] / u04
                G1_t[end,i,j]   = G2_t[1,i,j] / u04
            elseif c < 0
                B12_t[1,i,j]  = u04 * B11_t[end,i,j]
                B22_t[1,i,j]  = u04 * B21_t[end,i,j]
                G2_t[1,i,j]   = u04 * G1_t[end,i,j]
            end
        end
    end

    # if phi0 = 0 return
    if abs(phi0) < 1e-9
        return
    end

    @fastmath @inbounds for j in 1:Ny
        @inbounds for i in 1:Nx
            xi_t   = gauge_t.xi[1,i,j]
            A      = bulkconstrain2.A[1,i,j]

            # characteristic speed
            c = u0 * u0 * (A/2 - xi_t)

            # if c > 0, mode is entering grid1 from grid2;
            # if c < 0, mode is entering grid2 from grid1.
            # we assume here that grids merely touch at the interface
            if c > 0
                phi1_t[end,i,j] = phi2_t[1,i,j] / (u03 * phi03) + xi_t / (u0 * phi02)
            elseif c < 0
                phi2_t[1,i,j] = -phi0 * u02 * xi_t + u03 * phi03 * phi1_t[end,i,j]
            end
        end
    end

    nothing
end


text:
function estimate(mgc::GeneralizedCepstrum, x::AbstractArray;
                  norm=false,
                  kargs...)
    order = param_order(mgc)
    γ = glog_gamma(mgc)
    data = SPTK.gcep(x, order, γ; norm=norm, kargs...)
    SpectralParamState(mgc, data, true, norm)
end

function gcep(x::AbstractArray, order=25, γ=0.0; kargs...)
    estimate(MelGeneralizedCepstrum(order, 0.0, γ), x; kargs...)
end


text:
# This file was generated, do not modify it. # hide
xgbm = machine(xgb, X, y)
r = range(xgb, :gamma, lower=0, upper=10)
curve = learning_curve!(xgbm, range=r, resolution=30,
                        measure=cross_entropy);


text:
include("beta.jl")
include("bernoulli.jl")
include("dirichlet.jl")
include("categorical.jl")
include("gamma.jl")
include("gaussian.jl")
include("multivariate_gaussian.jl")

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
module SecantMethod

export secant_method

function secant_method(f::Function, x0::Number, x1::Number, args::Tuple=();
                       tol::AbstractFloat=1e-5, maxiter::Integer=50)
    for _ in 1:maxiter
       y1 = f(x1, args...)
       y0 = f(x0, args...)
       x = x1 - y1* (x1 - x0)/(y1 - y0)
       if abs(x - x1) < tol
          return x
       end
       x0 = x1
       x1 = x
    end
    error("Max iteration exceeded")
end

end


text:
module Arithmetic
export evaluate, Constant, BinaryOp, ADD, SUB, MUL, DIV

abstract type Expression end

struct Constant{T} <: Expression
  value::T
end

@enum Op begin
  ADD = 1
  SUB = 2
  MUL = 3
  DIV = 4
end

struct BinaryOp{T1, T2} <: Expression
  left::T1
  op::Op
  right::T2
end

evaluate(expr::Constant) = expr.value

function evaluate(expr::BinaryOp)
  left_val = evaluate(expr.left)
  right_val = evaluate(expr.right)
  if expr.op == ADD
    return left_val + right_val
  elseif expr.op == SUB
    return left_val - right_val
  elseif expr.op == MUL
    return left_val * right_val
  elseif expr.op == DIV
    return left_val / right_val
  else
    throw(error("Unknown op while evaluating BinaryOp"))
  end
end

end # module


text:
function halfcheetahState(; x::T=0.0, z::T=0.0, θ::T=0.0) where T
    mechanism = get_mechanism(:halfcheetah)
    initialize!(mechanism, :halfcheetah, x=x, z=z, θ=θ)

    Nb = length(mechanism.bodies)
    x = zeros(13 * Nb)

    for (i, body) in enumerate(mechanism.bodies)
        x2 = body.state.x2
        v15 = zeros(3)
        q2 = body.state.q2
        ϕ15 = zeros(3)
        x[13 * (i-1) .+ (1:13)] = [x2;  v15; vector(q2); ϕ15]
    end
    return x
end


text:
###############################################################################
# Callbacks
###############################################################################

function curl_write_cb(buff::Ptr{Cchar}, s::Csize_t, n::Csize_t, p::Ptr{Cvoid})::Csize_t
  ctxt = unsafe_pointer_to_objref(p)
  nbytes = s * n
  write(ctxt.resp.body, unsafe_string(buff, nbytes))
  ctxt.bytes_recd = ctxt.bytes_recd + nbytes

  nbytes
end

function writeptr(dst::Ptr{Cchar}, rd::ReadData, n::Csize_t)::Csize_t
  src = read(rd.src, n)
  n = length(src)
  ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt), dst, src, n)
  n
end

function curl_read_cb(out::Ptr{Cchar}, s::Csize_t, n::Csize_t, p::Ptr{Cvoid})::Csize_t
  ctxt = unsafe_pointer_to_objref(p)
  writeptr(out, ctxt.rd, s * n)
end


text:
using Replay

instructions = """
using UnicodePlots
histogram(randn(1000) .* 0.1, nbins = 15, closed = :right, xscale=:log10)
heatmap(collect(0:30) * collect(0:30)', xfact=0.1, yfact=0.1, xoffset=-1.5, colormap=:inferno)
"""

replay(instructions, julia_project=@__DIR__)


text:
# I took this from the ConjugatePriors package

# Used "Conjugate Bayesian analysis of the Gaussian distribution" by Murphy as
# a reference.  Note that there were some typos in that document so the code
# here may not correspond exactly.

struct NormalInverseGamma{T<:Real} <: ContinuousUnivariateDistribution
    mu::T
    v0::T     # scales variance of Normal
    shape::T
    scale::T

    function NormalInverseGamma{T}(mu::T, v0::T, sh::T, r::T) where T<:Real
    	v0 > zero(v0) && sh > zero(sh) && r > zero(r) || error("Both shape and scale must be positive")
    	new{T}(T(mu), T(v0), T(sh), T(r))
    end
end

function NormalInverseGamma(mu::Real, v0::Real, sh::Real, r::Real)
    T = promote_type(typeof(mu), typeof(v0), typeof(sh), typeof(r))
    return NormalInverseGamma{T}(T(mu),T(v0),T(sh),T(r))
end

mu(d::NormalInverseGamma) = d.mu
v0(d::NormalInverseGamma) = d.v0
shape(d::NormalInverseGamma) = d.shape
scale(d::NormalInverseGamma) = d.scale
rate(d::NormalInverseGamma) = 1. / d.scale

insupport(::Type{NormalInverseGamma}, x::T, sig2::T) where T<:Real =
    isfinite(x) && zero(sig2) <= sig2 < Inf

# Probably should guard agains dividing by and taking the log of 0.

function pdf(d::NormalInverseGamma, x::T, sig2::T) where T<:Real
    Zinv = d.scale.^d.shape / gamma(d.shape) / sqrt(d.v0 * 2.0*pi)
    return Zinv * 1.0/(sqrt(sig2)*sig2.^(d.shape+1.0)) * exp(-d.scale/sig2 - 0.5/(sig2*d.v0)*(x-d.mu).^2)
end

function logpdf(d::NormalInverseGamma, x::T, sig2::T) where T<:Real
    lZinv = d.shape*log(d.scale) - lgamma(d.shape) - 0.5*(log(d.v0) + log(2pi))
    return lZinv - 0.5*log(sig2) - (d.shape+1.)*log(sig2) - d.scale/sig2 - 0.5/(sig2*d.v0)*(x-d.mu).^2
end

function rand(d::NormalInverseGamma)
    # Guard against invalid precisions
    sig2 = rand(InverseGamma(d.shape, d.scale))
    if sig2 <= zero(Float64)
        sig2 = eps(Float64)
    end
    mu = rand(Normal(d.mu, sqrt(sig2*d.v0)))
    return mu, sig2
end


text:
using TuringModels
using Test

@testset "TuringModels" begin

include("../scripts/00/basic_example.jl")
include("../scripts/02/m2.1t.jl")
include("../scripts/04/m4.2t.jl")
#include("../scripts/08/m8.1t.jl")
#include("../scripts/08/m8.2t.jl")
#include("../scripts/08/m8.3t.jl")
#include("../scripts/08/m8.4t.jl")
#include("../scripts/10/m10.4t.jl")
#include("../scripts/11/m11.5t.jl")
#include("../scripts/12/m12.4t.jl")
#include("../scripts/13/m13.6t.jl")

end