text:
@testset "test insert_parent!" begin
    mtg = read_mtg("files/simple_plant.mtg")
    template = MultiScaleTreeGraph.MutableNodeMTG("/", "Shoot", 0, 1)
    max_node_id = [max_id(mtg)]
    length_before = length(mtg)

    insert_parent!(mtg[1][1], template, node -> typeof(node.attributes)(), max_node_id) # providing max_id

    @test length(mtg) == length_before + 1
    @test mtg[1][1].MTG.link == template.link
    @test mtg[1][1].MTG.symbol == template.symbol
    @test mtg[1][1].MTG.index == template.index
    @test mtg[1][1].MTG.scale == template.scale
    @test mtg[1][1].id == max_node_id[1]  # max_id is already incremented, but we added one more node

    mtg_2 = read_mtg("files/simple_plant.mtg")
    insert_parent!(mtg_2[1][1], template) # not providing max_id

    @test mtg_2 == mtg
    @test mtg_2[1][1] == mtg[1][1]

    mtg_3 = read_mtg("files/simple_plant.mtg")
    # Provide the template as a function (should take the previous template from its child):
    insert_parent!(
        mtg_3[1][1],
        x -> (
            link = x[1].MTG.link,
            symbol = x[1].MTG.symbol,
            index = x[1].MTG.index,
            scale = x[1].MTG.scale
        )
    ) # not providing max_id

    @test mtg_3 == mtg
    @test mtg_3[1][1].MTG == mtg_3[1][1][1][1].MTG

    # Test inserting a new root:
    mtg_4 = read_mtg("files/simple_plant.mtg")
    insert_parent!(mtg_4, template)

    @test get_root(mtg_4).MTG == template

    # Test inserting with attributes:

    mtg_5 = read_mtg("files/simple_plant.mtg")

    # Just a constant value:
    insert_parent!(mtg_5, template, x -> Dict{Symbol,Any}(:Length => 1,))

    # Using descendants values
    insert_parent!(
        get_root(mtg_5),
        template,
        node -> Dict{Symbol,Any}(:Total_Length => sum(descendants(node, :Length, ignore_nothing = true)),)
    )

    mtg_5 = get_root(mtg_5)

    @test mtg_5.MTG == template
    @test mtg_5[:Total_Length] == 32.0
    @test mtg_5[1][:Length] == 1
end

@testset "test insert_parents!" begin
    mtg = read_mtg("files/simple_plant.mtg")
    template = MutableNodeMTG("/", "Shoot", 0, 1)
    length_before = length(mtg)
    insert_parents!(
        mtg,
        template,
        node -> Dict{Symbol,Any}(:Total_Length => sum(descendants(node, :Length, ignore_nothing = true)),),
        scale = 2
    )

    @test length(mtg) == length_before + 1
    @test mtg[1][1].MTG == template
    @test mtg[1][1].name == "node_8"
    @test mtg[1][1].attributes == Dict{Symbol,Any}(:Total_Length => 32.0)
end

@testset "test insert_child!" begin
    mtg = read_mtg("files/simple_plant.mtg")
    template = MultiScaleTreeGraph.MutableNodeMTG("/", "Shoot", 0, 1)
    max_node_id = [max_id(mtg)]
    length_before = length(mtg)

    mtg_orig = deepcopy(mtg)

    insert_child!(mtg[1], template, node -> typeof(node.attributes)(), max_node_id) # providing max_id
    insert_parent!(mtg_orig[1][1], template)

    @test length(mtg) == length(mtg_orig)
    @test mtg_orig[1][1] == mtg_orig[1][1]
end

@testset "test insert_children!" begin
    mtg = read_mtg("files/simple_plant.mtg")
    template = MutableNodeMTG("/", "Rachis", 0, 4)
    length_before = length(mtg)
    insert_children!(mtg, template, symbol = "Leaf")

    @test length(mtg) == length_before + 2

    @test get_node(mtg, 8).MTG == template
    @test get_node(mtg, 9).MTG == template
end

@testset "test insert_generation!" begin
    mtg = read_mtg("files/simple_plant.mtg")
    template = MultiScaleTreeGraph.MutableNodeMTG("/", "Shoot", 0, 1)
    max_node_id = [max_id(mtg)]
    length_before = length(mtg)

    mtg_orig = deepcopy(mtg)

    insert_child!(mtg[1], template, node -> typeof(node.attributes)(), max_node_id) # providing max_id
    insert_generation!(mtg_orig[1], template)

    @test length(mtg) == length(mtg_orig)
    @test mtg_orig[1][1] == mtg_orig[1][1]
end

@testset "test insert_generations!" begin
    mtg = read_mtg("files/simple_plant.mtg")
    template = MutableNodeMTG("/", "Rachis", 0, 4)
    length_before = length(mtg)
    insert_generations!(mtg, template, symbol = "Internode")

    @test length(mtg) == length_before + 2

    @test get_node(mtg, 8).MTG == template
    @test get_node(mtg, 9).MTG == template
end

@testset "test insert_sibling!" begin
    mtg = read_mtg("files/simple_plant.mtg")
    template = MultiScaleTreeGraph.MutableNodeMTG("/", "Individual", 0, 1)
    max_node_id = [max_id(mtg)]
    length_before = length(mtg)

    mtg_orig = deepcopy(mtg)

    insert_child!(mtg, template, node -> typeof(node.attributes)(), max_node_id) # providing max_id
    insert_sibling!(mtg_orig[1], template)

    @test length(mtg) == length(mtg_orig)
    @test mtg_orig[2] == mtg_orig[2]
end

@testset "test insert_siblings!" begin
    mtg = read_mtg("files/simple_plant.mtg")
    template = MutableNodeMTG("+", "Leaf", 0, 3)
    length_before = length(mtg)
    # Add a new leaf at each leaf node:
    insert_siblings!(mtg, template, symbol = "Leaf")

    @test length(mtg) == length_before + 2

    @test get_node(mtg, 8).MTG == template
    @test get_node(mtg, 9).MTG == template
end


text:
export read_file

"""
read_file(s)

Read entire SEGY file from stream 's'.
"""
function read_file(s::IO, warn_user::Bool; start_byte::Int = 3600,
                                           end_byte::Int = position(seekend(s)))

    # Read File Header
    fh = read_fileheader(s)

    # Move to start of block
    seek(s, start_byte)

    # Check datatype of file
    datatype = Float32
    if fh.bfh.DataSampleFormat == 1
        datatype = IBMFloat32
    elseif fh.bfh.DataSampleFormat != 5
        @error "Data type not supported ($(fh.bfh.DataSampleFormat))"
    end

    # Check fixed length trace flag
    (fh.bfh.FixedLengthTraceFlag!=1 & warn_user) && @warn "Fixed length trace flag set in stream: $s"

    ## Check for extended text header

    # Read traces
    ntraces = Int((end_byte - start_byte)/(240 + fh.bfh.ns*4))

    # Preallocate memory
    headers = Array{BinaryTraceHeader, 1}(undef, ntraces)
    data = Array{datatype, 2}(undef, fh.bfh.ns, ntraces)
    th_b2s = th_byte2sample()

    # Read each trace
    for trace in 1:ntraces

        read_trace!(s, fh.bfh, datatype, headers, data, trace, th_b2s)

    end

    return SeisBlock(fh, headers, data)
end

"""
read_file(s, keys)

Read entire SEGY file from stream 's', only reading the header values in 'keys'.
"""
function read_file(s::IO, keys::Array{String, 1}, warn_user::Bool;
                                                  start_byte::Int = 3600,
                                                  end_byte::Int = position(seekend(s)))

    # Read File Header
    fh = read_fileheader(s)

    # Move to start of block
    seek(s, start_byte)

    # Check datatype of file
    datatype = Float32
    if fh.bfh.DataSampleFormat == 1
        datatype = IBMFloat32
    else
        @error "Data type not supported ($(fh.bfh.DataSampleFormat))"
    end

    # Check fixed length trace flag
    (fh.bfh.FixedLengthTraceFlag!=1 & warn_user) && @warn "Fixed length trace flag set in stream: $s"

    ## Check for extended text header

    # Read traces
    ntraces = Int((end_byte - start_byte)/(240 + fh.bfh.ns*4))

    # Preallocate memory
    headers = Array{BinaryTraceHeader, 1}(undef, ntraces)
    data = Array{datatype, 2}(undef, fh.bfh.ns, ntraces)
    th_b2s = th_byte2sample()

    # Read each trace
    for trace in 1:ntraces

        read_trace!(s, fh.bfh, datatype, headers, data, trace, keys, th_b2s)

    end

    return SeisBlock(fh, headers, data)
end

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg, Base.BinaryPlatforms

name = "SoapyLMS7"
version = v"20.10.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/myriadrf/LimeSuite.git", "1480bfeaf4de211c40813fc5ca161b1b644778ec")
]

dependencies = [
    Dependency("soapysdr_jll"; compat="0.8.0"),
    Dependency("libusb_jll"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/LimeSuite
mkdir _build
cd _build

CMAKE_OPTIONS=(
    # Disable a bunch of things we don't care about
    -DENABLE_QUICKTEST=OFF
    -DENABLE_OCTAVE=OFF
    -DDOWNLOAD_IMAGES=FALSE
    -DENABLE_GUI=OFF
    -DENABLE_NOVENARF7=OFF

    # Disable SIMD for now; enable later if we need it
    # by expanding our microarchitectures
    -DDEFAULT_SIMD_FLAGS=cross

    # Enable the one thing we do care about
    -DENABLE_SOAPY_LMS7=ON
)

cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DCMAKE_BUILD_TYPE=Release \
      "${CMAKE_OPTIONS[@]}" \
      ..
make -j${nproc}
make install
if [[ "${target}" == *-apple-* ]]; then
    mv ${libdir}/SoapySDR/modules0.8/libLMS7Support.so  ${libdir}/SoapySDR/modules0.8/libLMS7Support.dylib
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> os(p) != "windows", supported_platforms(;experimental=true))
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = Product[
    # TODO: These should eventually be a separate output from this builder
    ExecutableProduct("LimeUtil", :LimeUtil),
    LibraryProduct("libLimeSuite", :libLimeSuite),
    LibraryProduct("libLMS7Support", :libLMS7Support, ["lib/SoapySDR/modules0.8/"]),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@testset "Validate vortex and symmetric_vortex" begin
    atol = 1000*eps()

    Nv = 20
    Nx = 3*Nv
    Ne = 5

    sensors = -11.0:0.5:10.0
    U = complex(0.0)

    config = let Nv = Nv,
                 U = U,
                 ss = sensors, Δt = 5e-3, δ = 1e-1,
                 ϵX = 1e-3, ϵΓ = 1e-3,
                 β = 1.0,
                 ϵY = 1e-2
        VortexConfig(Nv, U, ss, Δt, δ, ϵX, ϵΓ, β, ϵY)
    end
    # Test vortex routine
    x = rand(Nx)
    X = repeat(x, 1, Ne)
    X0 = deepcopy(X)
    t0 = 1.8
    sys = state_to_lagrange(x, config)
    sys0 = deepcopy(sys)
    sys₊ = deepcopy(sys)
    vels = self_induce_velocity(sys, t0)

    advect!(sys₊, sys0, vels, config.Δt)

    cachevels = allocate_velocity(state_to_lagrange(zeros(Nx), config))
    X1, t1 = vortex(deepcopy(X0), t0, 0, Nx, cachevels, config)

    for i=1:Ne
        @test isapprox(X1[:,i], lagrange_to_state(sys₊, config), atol = atol)
    end

    # Test that the symmetry is correclty enforced
    X = rand(Nx, Ne)
    X0 = deepcopy(X)
    t0 = 1.5

    cachevels = allocate_velocity(state_to_lagrange(zeros(Nx), config))

    X1, t1 = vortex(deepcopy(X0), t0, 0, Nx, cachevels, config)

    X2, t2 = symmetric_vortex(deepcopy(X0), t0, 0, Nx, cachevels, config)

    @test isapprox(X1, X2, atol = atol)
    @test isapprox(t1, t0 + config.Δt, atol = atol)
    @test isapprox(t2, t0 + config.Δt, atol = atol)
end

@testset "Validate vortex and symmetric_vortex with freestream" begin
    atol = 1000*eps()

    Nv = 20
    Nx = 3*Nv
    Ne = 5

    sensors = -11.0:0.5:10.0
    U = rand(ComplexF64)
    freestream = Freestream(U)

    config = let Nv = Nv,
                 U = U,
                 ss = sensors, Δt = 5e-3, δ = 1e-1,
                 ϵX = 1e-3, ϵΓ = 1e-3,
                 β = 1.0,
                 ϵY = 1e-2
        VortexConfig(Nv, U, ss, Δt, δ, ϵX, ϵΓ, β, ϵY)
    end
    # Test vortex routine
    x = rand(Nx)
    X = repeat(x, 1, Ne)
    X0 = deepcopy(X)
    t0 = 1.8
    sys = state_to_lagrange(x, config)
    sys0 = deepcopy(sys)
    sys₊ = deepcopy(sys)
    vels = self_induce_velocity(sys, t0)
    induce_velocity!(vels, sys, freestream, t0)

    advect!(sys₊, sys0, vels, config.Δt)

    cachevels = allocate_velocity(state_to_lagrange(zeros(Nx), config))
    X1, t1 = vortex(deepcopy(X0), t0, 0, Nx, cachevels, config, withfreestream = true)
    X2, t1 = vortex(deepcopy(X0), t0, 0, Nx, cachevels, config, withfreestream = false)

    X2[1:3:end,:] .+= real(config.U)*config.Δt
    X2[2:3:end,:] .+= imag(config.U)*config.Δt

    for i=1:Ne
        @test isapprox(X1[:,i], lagrange_to_state(sys₊, config), atol = atol)
        @test isapprox(X2[:,i], lagrange_to_state(sys₊, config), atol = atol)
    end

    # Test that the symmetry is correclty enforced
    X = rand(Nx, Ne)
    X0 = deepcopy(X)
    t0 = 1.5

    cachevels = allocate_velocity(state_to_lagrange(zeros(Nx), config))

    X1, t1 = vortex(deepcopy(X0), t0, 0, Nx, cachevels, config; withfreestream = true)

    X2, t2 = symmetric_vortex(deepcopy(X0), t0, 0, Nx, cachevels, config; withfreestream = true)

    @test isapprox(X1, X2, atol = atol)
    @test isapprox(t1, t0 + config.Δt, atol = atol)
    @test isapprox(t2, t0 + config.Δt, atol = atol)
end


text:
using JLD
using ProgressMeter

open("81.jl.out.vector") do f
    global nvocab, h

    nvocab, h = parse.(Int, split(readline(f)))
    global X = zeros(h, nvocab) # to be transposed later
    global vocab = Vector{String}(nvocab)

    @showprogress 1 "ahahaha" for (idx,line) in enumerate(readlines(f))
        s = split(line)
        word = s[1]
        fv = parse.(Float32, s[2:end])
        @assert length(fv) == h

        vocab[idx] = word
        X[:,idx] = fv
    end

    X = X'
    @assert size(X) == (nvocab, h)

    global vocab2idx = Dict()
    global idx2vocab = Dict()
    for (idx,word) in enumerate(vocab)
        vocab2idx[word] = idx
        idx2vocab[idx] = word
    end
end

if true
    jldopen("90.jld", "w") do file
        write(file, "X", X)
        write(file, "vocab", vocab)
        write(file, "vocab2idx", vocab2idx)
        write(file, "idx2vocab", idx2vocab)
    end
end


text:
using PyPlot, Revise, transD_GP, Random
## frequencies
nFreqsPerDecade     = 7
freqLowLimit        = 1e-1
freqHighLimit       = 1e5
freqs = 10 .^(log10(freqLowLimit):1/nFreqsPerDecade:log10(freqHighLimit))
## model
zfixed   = [-1e5,   0,    ]
rho      = [1e12,   100,  ]
nmax = 200
##  geometry
rRx = 100.
zRx = -0.011
zTx = -0.01
nkᵣeval = 200
## VMD
modelprimary = true
Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      nmax   = nmax,
                      zTx    = zTx,
                      rRx    = rRx,
                      freqs  = freqs,
                      zRx    = zRx,
                      nkᵣeval = nkᵣeval,
                      modelprimary = modelprimary)
transD_GP.AEM_VMD_HMD.getfieldFD!(Fvmd, zfixed, rho)
## now use a tiny loop radius - VMD approximation is worse as radius gets larger
rTx = 0.001
Floop = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      nmax   = nmax,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      freqs  = freqs,
                      zRx    = zRx,
                      nkᵣeval = nkᵣeval,
                      modelprimary = modelprimary)
transD_GP.AEM_VMD_HMD.getfieldFD!(Floop, zfixed, rho)
## Compare H vertical
f, ax = plt.subplots(2,1, sharex=true)
ax[1].loglog(freqs,abs.(imag(Fvmd.HFD_z)), label="imaginary, VMD")
ax[1].loglog(freqs,abs.(real(Fvmd.HFD_z)), label="real, VMD")
ax[1].loglog(freqs,abs.(imag(Floop.HFD_z)), "*", label="imaginary, loop")
ax[1].loglog(freqs,abs.(real(Floop.HFD_z)), "*", label="real, loop")
ax[1].legend()
ax[1].set_title("Compare with W&H Fig 4.2")
## Compare H radial
getradialH = true
Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      nmax   = nmax,
                      zTx    = 0,
                      rRx    = rRx,
                      freqs  = freqs,
                      zRx    = 0,
                      nkᵣeval = nkᵣeval,
                      modelprimary = false, # says W&H for within plane of observation !!
                      getradialH = getradialH)
transD_GP.AEM_VMD_HMD.getfieldFD!(Fvmd, zfixed, rho)
ax[2].loglog(freqs,abs.(imag(Fvmd.HFD_r)), label="imaginary, VMD")
ax[2].loglog(freqs,abs.(real(Fvmd.HFD_r)), label="real, VMD")
ax[2].set_title("Compare with W&H Fig 4.3")
ax[2].set_xlim(extrema(freqs))
ax[2].set_ylim(1e-12,1e-6)
ax[1].grid()
ax[2].grid()
ax[2].set_xlabel("Hz")
legend()
## now test loop field, receiver at center
zTx = -0.01
zRx = -0.02
rRx = 0.001
rTx = 50.0
Floopin = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      nmax   = nmax,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      freqs  = freqs,
                      zRx    = zRx,
                      nkᵣeval = nkᵣeval,
                      modelprimary = modelprimary)
transD_GP.AEM_VMD_HMD.getfieldFD!(Floopin, zfixed, rho)
##
figure()
loglog(freqs,pi*rTx^2*abs.(imag(Floopin.HFD_z)), label="imaginary, rx in = $(Floopin.rxwithinloop)")
loglog(freqs,pi*rTx^2*abs.(real(Floopin.HFD_z)), label="real, rx in = $(Floopin.rxwithinloop)")
xlim(extrema(freqs))
xlabel("Hz")
legend()
title("Compare with W&H Fig 4.7")
grid()


text:
using SymPy

wmotor1 = Sym("wmotor1");
wmotor2 = Sym("wmotor2");
wmotor3 = Sym("wmotor3");
wmotor4 = Sym("wmotor4");

## Propeller Forces
keng = 6.11e-3;
f1_body = [0;0;-keng*wmotor1];
f2_body = [0;0;-keng*wmotor2];
f3_body = [0;0;-keng*wmotor3];
f4_body = [0;0;-keng*wmotor4];

## Propeller Moments
kmeng = 1.5e-4;
m1_body = [0;0;-kmeng*wmotor1];
m2_body = [0;0; kmeng*wmotor2];
m3_body = [0;0;-kmeng*wmotor3];
m4_body = [0;0; kmeng*wmotor4];

## Moment arm from CG to Propeller Hub;
Larm = 0.15;
r1_arm_body = [ Larm*cosd(45); -Larm*sind(45); 0.0];
r2_arm_body = [ Larm*cosd(45);  Larm*sind(45); 0.0];
r3_arm_body = [-Larm*cosd(45);  Larm*sind(45); 0.0];
r4_arm_body = [-Larm*cosd(45); -Larm*sind(45); 0.0];

## Total Moment due to Motors
m1_body = m1_body + cross(r1_arm_body, f1_body);
m2_body = m2_body + cross(r2_arm_body, f2_body);
m3_body = m3_body + cross(r3_arm_body, f3_body);
m4_body = m4_body + cross(r4_arm_body, f4_body);

mtot = m1_body + m2_body + m3_body + m4_body;
ftot = f1_body + f2_body + f3_body + f4_body;

AA = jacobian([mtot;ftot[3]], [wmotor1, wmotor2, wmotor3, wmotor4]);


text:
using MBSinBLG, Quantica, Arpack, Parameters
using ElectronDisplay

function test(; α = 0, Ey = 4, kw...)
    p = Params(Ls = 0, scale = 40, λ = 5, α = α,
        μN = 0.6, Δ = 1, d = 0, τ = 1, EZ = SA[0, Ey, 0]);
    lnlist = 50:25:500
    x, y = cprvsln(lnlist, p, nev = 8; kw...)
    return x, y
end

function cprvsln(list, p; nev = 16, kw...)
    elist = zeros(Float64, nev, length(list))
    for i in 1:length(list)
        ph = MBSinBLG.rectangle_squid(reconstruct(p, Ln = list[i]); kw...)
        s = spectrum(ph(ϕ = π), method = ArpackPackage(sigma = 1e-7, nev = nev))
        elist[:,i] = s.energies
    end
    return list, elist
end

x_b, y_b = test(α = 0, Ey = 4, W = 1440)
x_d, y_d = test(α = 1, Ey = 4, W = 1440)

x_b2, y_b2 = test(α = 0, Ey = 4, W = 644)
x_d2, y_d2 = test(α = 1, Ey = 4, W = 644)

function plottest(x, y,  ylims = (-0.1, 0.1))
    fig = Figure(resolution = (500, 300), font = "Times New Roman")
    axa = Axis(fig[1,1])

    lines!(axa, x, y[1,:])
    for i in 2:size(y, 1)
        lines!(axa, x, y[i,:])
    end
    axa.xlabel = " Ln [nm]"
    axa.ylabel = "E [meV]"
    ylims!(axa, ylims)
    fig
end

max_stars_repo_name: fernandopenaranda/MBSinBLG

max_issues_repo_name: fernandopenaranda/MBSinBLG

max_forks_repo_name: fernandopenaranda/MBSinBLG


text:
using Markdown
using InteractiveUtils

using DrWatson

begin
  @quickactivate "SR2TuringPluto"
  using Turing
  using StatisticalRethinking
end

md"## Clip-05-01-27t.jl"

begin
  df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame)
  df.D = zscore(df.Divorce)
  df.M = zscore(df.Marriage)
  df.A = zscore(df.MedianAgeMarriage)
end

std(df.MedianAgeMarriage)

@model function m5_1_A(A, D)
    a ~ Normal(0, 0.2)
    bA ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, A, bA)
    D ~ MvNormal(μ, σ)
end

begin
    m5_1_At = m5_1_A(df.A, df.D)
  quap5_1_At = quap(m5_1_At)
  post5_1_At = DataFrame(rand(quap5_1_At.distr, 1000)', quap5_1_At.params)

  A_seq = range(-3, 3.2, length = 30)
  mu5_1_At = lin(post5_1_At.a', A_seq, post5_1_At.bA') |> meanlowerupper

  scatter(df.A, df.D, alpha = 0.4, legend = false)
  plot!(A_seq, mu5_1_At.mean, ribbon =
        (mu5_1_At.mean .- mu5_1_At.lower, mu5_1_At.upper .- mu5_1_At.mean))
end

@model function m5_1_M(M, D)
    a ~ Normal(0, 0.2)
    bM ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, M, bM)
    D ~ MvNormal(μ, σ)
end

begin
    m5_1_Mt = m5_1_M(df.M, df.D)
    quap5_1_Mt = quap(m5_1_Mt)
    post5_1_Mt = DataFrame(rand(quap5_1_Mt.distr, 1000)', quap5_1_Mt.params)

    M_seq = range(-3, 3.2, length = 30)
    mu5_1_Mt = lin(post5_1_Mt.a', M_seq, post5_1_Mt.bM') |> meanlowerupper

    scatter(df.M, df.D, alpha = 0.4, legend = false)
    plot!(M_seq, mu5_1_Mt.mean, ribbon =
        (mu5_1_Mt.mean .- mu5_1_Mt.lower, mu5_1_Mt.upper .- mu5_1_Mt.mean))
end

@model function m5_1_A_M(A, M, D)
    a ~ Normal(0, 0.2)
    bM ~ Normal(0, 0.5)
    bA ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, A, bA, M, bM)
    D ~ MvNormal(μ, σ)
end

begin
    m5_1_A_Mt = m5_1_A_M(df.A, df.M, df.D)
    quap5_1_A_Mt = quap(m5_1_A_Mt)
    post5_1_A_Mt = DataFrame(rand(quap5_1_A_Mt.distr, 1000)', quap5_1_A_Mt.params)

    A_M_seq = range(-3, 3.2, length = 30)
    mu5_1_A_Mt = lin(post5_1_A_Mt.a', A_M_seq, post5_1_A_Mt.bM') |> meanlowerupper

    scatter(df.M, df.D, alpha = 0.4, legend = false)
    plot!(M_seq, mu5_1_A_Mt.mean, ribbon =
        (mu5_1_A_Mt.mean .- mu5_1_A_Mt.lower, mu5_1_A_Mt.upper .- mu5_1_A_Mt.mean))
end

function plotcoef1(
  models::Vector{T},
  mnames::Vector{String},
  pars::Vector{Symbol};
  fig="", title="", func=nothing,
  sampler=NUTS(0.65), nsamples=2000, nchains=4) where {T <: DynamicPPL.Model}

  # mnames = [nameof(m) for m in models]
  levels = length(models) * (length(pars) + 1)
  colors = [:blue, :red, :green, :darkred, :black]

  s = Vector{NamedTuple}(undef, length(models))
  for (mindx, mdl) in enumerate(models)
    if isnothing(func)
        chns = mapreduce(c -> sample(mdl, sampler, nsamples),
            chainscat, 1:nchains)
        df = DataFrame(Array(chns), names(chns, [:parameters]))
        m, l, u = estimparam(df)
        d = Dict{Symbol, NamedTuple}()
        for (indx, par) in enumerate(names(chns, [:parameters]))
            d[par] = (mean=m[indx], lower=l[indx], upper=u[indx])
        end
        s[mindx] =   (; d...)
    else
        quap_mdl = quap(mdl)
        post = rand(quap_mdl.distr, 10_000)
        df = DataFrame(post', [keys(quap_mdl.coef)...])
        m, l, u = estimparam(df)
        d = Dict{Symbol, NamedTuple}()
        for (indx, par) in enumerate([keys(quap_mdl.coef)...])
            d[par] = (mean=m[indx], lower=l[indx], upper=u[indx])
        end
        s[mindx] =   (; d...)
    end
  end

  xmin = 0; xmax = 0.0
  for i in 1:length(s)
    for par in pars
      syms = Symbol.(keys(s[i]))
      if Symbol(par) in syms
        mp = s[i][par].mean
        xmin = min(xmin, s[i][par].lower)
        xmax = max(xmax, s[i][par].upper)
      end
    end
  end

  ylabs = String[]
  for j in 1:length(models)
    for i in 1:length(pars)
      l = length(String(pars[i]))
      str = repeat(" ", 9-l) * String(pars[i])
      append!(ylabs, [str])
    end
    l = length(mnames[j])
    str = mnames[j] * repeat(" ", 9-l)
    append!(ylabs, [str])
  end

  ys = [string(ylabs[i]) for i = 1:length(ylabs)]
  p = plot(xlims=(xmin, xmax), leg=false, framestyle=:grid)
  title!(title)
  yran = range(1, stop=length(ylabs), length=length(ys))
  yticks!(yran, ys)

  line = 0
  for (mindx, model) in enumerate(models)
    line += 1
    hline!([line] .+ length(pars), color=:darkgrey, line=(2, :dash))
    for (pindx, par) in enumerate(pars)
      line += 1
      syms = Symbol.(keys(s[mindx]))
      if Symbol(par) in syms
        ypos = (line - 1)
        mp = s[mindx][Symbol(par)].mean
        lower = s[mindx][Symbol(par)].lower
        upper = s[mindx][Symbol(par)].upper
        plot!([lower, upper], [ypos, ypos], leg=false, color=colors[pindx])
        scatter!([mp], [ypos], color=colors[pindx])
      end
    end
  end
  if length(fig) > 0
    savefig(p, fig)
  end
  (s, p)
end

s1, p1 = plotcoef(
    [m5_1_At, m5_1_Mt, m5_1_A_Mt],
    [:a, :bM, :bA, :σ]
)

s1 |> display

plot(p1)

gui()

s2, p2 = plotcoef(
    [m5_1_At, m5_1_Mt, m5_1_A_Mt],
    [:a, :bM, :bA, :σ];
    func=quap
)

s2 |> display

plot(p2)

gui()

s3, p3 = plotcoef(
    [m5_1_At],
    [:a, :bM, :bA, :σ]
)

s3 |> display

plot(p3)

gui()

s4, p4 = plotcoef(
    m5_1_At,
    [:a, :bM, :bA, :σ];
    func=quap
)

s4 |> display

plot(p4)

gui()


text:
# Test MathProgBase stuff

using MathProgBase
using Compat.LinearAlgebra
using Compat.SparseArrays

solver = IpoptSolver()
sol = linprog([-1,0],[2 1],'<',1.5,solver)
@test sol.status == :Optimal
@test abs(sol.objval - -0.75) <= 1e-6
@test norm(sol.sol - [0.75,0.0]) <= 1e-6

sol = linprog([-1,0],sparse([2 1]),'<',1.5,solver)
@test sol.status == :Optimal
@test abs(sol.objval - -0.75) <= 1e-6
@test norm(sol.sol - [0.75,0.0]) <= 1e-6

# test infeasible problem:
# min x
# s.t. 2x+y <= -1
# x,y >= 0
sol = linprog([1,0],[2 1],'<',-1,solver)
#    @test sol.status == :Infeasible

# test unbounded problem:
# min -x-y
# s.t. -x+2y <= 0
# x,y >= 0
sol = linprog([-1,-1],[-1 2],'<',[0],solver)
#    @test sol.status == :Unbounded

function mathprogbase_file(file::String)
    if VERSION >= v"0.7-"
        return joinpath(dirname(dirname(pathof(MathProgBase))), "test", file)
    else
        return joinpath(Pkg.dir("MathProgBase"), "test", file)
    end
end

include(mathprogbase_file("nlp.jl"))
nlptest(IpoptSolver())
nlptest_nohessian(IpoptSolver())
convexnlptest(IpoptSolver())
rosenbrocktest(IpoptSolver())

include(mathprogbase_file("quadprog.jl"))
quadprogtest(IpoptSolver())
qpdualtest(IpoptSolver())

# Test retoration only options
#
# Warm start with infeasible solution, force restoration on initial iteration.
# But limit to 0 iterations. Forces :UserLimit exit.
m = MathProgBase.NonlinearModel(IpoptSolver(start_with_resto="yes", resto_max_iter=0))
l = [1,1,1,1]
u = [5,5,5,5]
lb = [25, 40]
ub = [Inf, 40]
MathProgBase.loadproblem!(m, 4, 2, l, u, lb, ub, :Min, HS071())
MathProgBase.setwarmstart!(m,[0,15,15,11])
MathProgBase.optimize!(m)
@test MathProgBase.status(m) == :UserLimit


text:
# --------------------------------------------------------------------
# Requires
# --------------------------------------------------------------------
import .GLM
using StatsModels
using StatsBase
using Tables: columntable, istable

import StatsModels: TableRegressionModel, RegressionModel
import StatsBase: modelmatrix

const INNERMOD = Union{GLM.GeneralizedLinearModel, GLM.LinearModel}
# --------------------------------------------------------------------
# GLM Methods
# --------------------------------------------------------------------

# TODO: Find a good name for it
# invinvpseudohessian
# ....
invpseudohessian(m::TableRegressionModel{T}) where T<:INNERMOD = GLM.invchol(m.model.pp).*dispersion(m.model.rr)
invpseudohessian(m::T) where T<:INNERMOD = GLM.invchol(m.pp).*dispersion(m.rr)

chol(m::TableRegressionModel{T}) where T<:INNERMOD = chol(m.model)
chol(m::T) where T<:INNERMOD = m.pp.chol

modmatrix(m::TableRegressionModel{T}) where T<:INNERMOD = modmatrix(m.model)
modmatrix(m::T) where T<:GLM.GeneralizedLinearModel = sqrt.(m.rr.wrkwt).*modelmatrix(m)
function modmatrix(m::T) where T<:GLM.LinearModel
    X = modelmatrix(m)
    if !isempty(m.rr.wts)
        sqrt.(m.rr.wts).*X
    else
        copy(X)
    end
end

numobs(m::TableRegressionModel) = length(m.model.rr.y)
numobs(m::INNERMOD) = length(m.rr.y)
dof_resid(m::TableRegressionModel) = numobs(m) - length(coef(m))
dof_resid(m::INNERMOD) = numobs(m) - length(coef(m))

StatsModels.hasintercept(m::TableRegressionModel) = "(Intercept)" ∈ coefnames(m)
interceptindex(m::INNERMOD) = findfirst(map(x->allequal(x), eachcol(modmatrix(m))))
function StatsModels.hasintercept(m::INNERMOD)
    hasint = findfirst(map(x->allequal(x), eachcol(modmatrix(m))))
    hasint === nothing ? false : true
end

dispersion(m::TableRegressionModel{T}) where T<:GLM.GeneralizedLinearModel = dispersion(m.model.rr)
dispersion(m::TableRegressionModel{T}) where T<:GLM.LinearModel = 1
dispersion(m::GLM.GeneralizedLinearModel) = dispersion(m.rr)
dispersion(m::GLM.LinearModel) = 1
dispersion(rr::GLM.GlmResp{T1, T2, T3}) where {T1, T2, T3} = 1
dispersion(rr::GLM.LmResp) = 1
function dispersion(rr::GLM.GlmResp{T1, T2, T3}) where {T1, T2<:Union{GLM.Gamma, GLM.Bernoulli, GLM.InverseGaussian}, T3}
    sum(abs2, rr.wrkwt.*rr.wrkresid)/sum(rr.wrkwt)
end

resid(m::TableRegressionModel{T}) where T<:INNERMOD = resid(m.model)
function resid(m::T) where T<:GLM.LinearModel
    if !isempty(m.rr.wts)
        sqrt.(m.rr.wts).*residuals(m)
    else
        copy(residuals(m))
    end
end
function resid(m::T) where T<:GLM.GeneralizedLinearModel
    sqrt.(m.rr.wrkwt).*m.rr.wrkresid
end

momentmatrix(m::TableRegressionModel{T}) where T<:INNERMOD = momentmatrix(m.model)
momentmatrix(m::INNERMOD) = (modmatrix(m).*resid(m))./dispersion(m)

# TODO: move to the interface file
# hasresiduals(m::INNERMOD) = true
# hasmodelmatrix(m::TableRegressionModel{T}) where T<:INNERMOD = true

# --------------------------------------------------------------------
# HAC GLM Methods
# --------------------------------------------------------------------
function StatsBase.vcov(k::HAC, m; prewhite=false, dof_adjustment::Bool=true, scale::Real=1)
    return _vcov(k, m, prewhite, dof_adjustment, scale)
end

function _vcov(k::HAC, m, prewhite::Bool, dof_adjustment::Bool, scale::Real)
    setupkernelweights!(k, m)
    B  = invpseudohessian(m)
    mm = momentmatrix(m)
    A = lrvar(k, mm; prewhite=prewhite, demean=false)
    scale *= (dof_adjustment ? numobs(m)^2/dof_resid(m) : numobs(m))
    V = Symmetric((B*A*B).*scale)
    return V
end

vcovmatrix(
    k::HAC,
    m::RegressionModel,
    factorization=Cholesky;
    prewhite=false,
    dof_adjustment::Bool=true,
    scale::Real=1,
) = _vcovmatrix(k, m, prewhite, dof_adjustment, scale, factorization)

function _vcovmatrix(
    k::HAC,
    m::RegressionModel,
    prewhite::Bool,
    dof_adjustment::Bool,
    scale::Real,
    ::Type{Cholesky},
)
    V = _vcov(k, m, prewhite, dof_adjustment, scale)
    return CovarianceMatrix(cholesky(V, check=true), k, V)
end

function _vcovmatrix(
    k::HAC,
    m::RegressionModel,
    prewhite::Bool,
    dof_adjustment::Bool,
    scale::Real,
    ::Type{SVD},
)
    V = _vcov(k, m, prewhite, dof_adjustment, scale)
    return CovarianceMatrix(svd(V.data), k, V)
end

# --------------------------------------------------------------------
# HC GLM Methods
# --------------------------------------------------------------------
hatmatrix(m::TableRegressionModel{T}, x) where T<:INNERMOD = hatmatrix(m.model, x)

function hatmatrix(m::T, x) where T<:INNERMOD
    cf = m.pp.chol.UL::UpperTriangular
    rdiv!(x, cf)
    return sum(x.^2, dims = 2)
 end

adjfactor(k::HC0, m::RegressionModel, x) = one(first(x))
adjfactor(k::HC1, m::RegressionModel, x) = numobs(m)./dof_resid(m)
adjfactor(k::HC2, m::RegressionModel, x) = one(first(x)) ./(one(first(x)).-hatmatrix(m, x))
adjfactor(k::HC3, m::RegressionModel, x) = one(first(x))./(one(first(x)).-hatmatrix(m, x)).^2

function adjfactor(k::HC4, m::RegressionModel, x)
    n, p = size(x)
    tone = one(eltype(x))
    h = hatmatrix(m, x)
    @inbounds for j in eachindex(h)
        delta = min(4, n*h[j]/p)
        h[j] = tone/(tone-h[j])^delta
    end
    return h
end

function adjfactor(k::HC4m, m::RegressionModel, x)
    n, p = size(x)
    tone = one(eltype(x))
    h = hatmatrix(m, x)
    @inbounds for j in eachindex(h)
        delta = min(tone, n*h[j]/p) + min(1.5, n*h[j]/p)
        h[j] = tone/(tone-h[j])^delta
    end
    return h
end

function adjfactor(k::HC5, m::RegressionModel, x)
    n, p = size(x)
    tone = one(eltype(x))
    h = hatmatrix(m, x)
    mx = max(n*0.7*maximum(h)/p, 4)
    @inbounds for j in eachindex(h)
        alpha = min(n*h[j]/p, mx)
        h[j] = tone/sqrt((tone-h[j])^alpha)
    end
    return h
end

adjust!(m, adj::AbstractFloat) = m
adjust!(m, adj::AbstractMatrix) = m.*sqrt.(adj)

StatsBase.vcov(k::HC, m::RegressionModel; scale::Real=1) = _vcov(k, m, scale)

function _vcov(k::HC, m::RegressionModel, scale)
    B  = invpseudohessian(m)
    mm = momentmatrix(m)
    adj = adjfactor(k, m, modmatrix(m))
    mm = adjust!(mm, adj)
    scale *= length(adj) > 1 ? one(first(adj)) : adj
    A = mm'*mm
    return Symmetric((B*A*B).*scale)
end

vcovmatrix(k::HC, m::RegressionModel, factorization=Cholesky; scale::Real=1) =
    _vcovmatrix(k, m, scale, Val{:factorization})

function _vcovmatrix(k::HC, m::RegressionModel, scale::Real, ::Type{Cholesky})
    V = _vcov(k, m, scale)
    CovarianceMatrix(cholesky(V, check=true), k, V)
end

function _vcovmatrix(k::HC, m::RegressionModel, scale::Real, ::Type{SVD})
    V = _vcov(k, m, scale)
    CovarianceMatrix(svd(V), k, V)
end

# --------------------------------------------------------------------
# CRHC GLM Methods
# --------------------------------------------------------------------

function installcache(k::CRHC, m::RegressionModel)
    X = modmatrix(m)
    res = resid(m)
    f = categorize(k.cl)
    (X, res), sf = bysort((X, res), f)
    ci = clustersintervals(sf)
    p = size(X, 2)
    cf = chol(m)
    Shat = Matrix{eltype(res)}(undef,p,p)
    return CRHCCache(similar(X), X, res, similar(X, (0,0)), cf, Shat, ci, sf)
end

function StatsBase.vcov(k::CRHC, m::RegressionModel; scale::Real=1)
    knew = recast(k, m)
    length(knew.cl) == numobs(m) || throw(ArgumentError(k, "the length of the cluster variable must be $(numobs(m))"))
    cache = installcache(knew, m)
    return _vcov(knew, m, cache, scale)
end

function vcovmatrix(k::CRHC, m::RegressionModel, factorization = Cholesky; scale::Real=1)
    knew = recast(k, m)
    cache = installcache(knew, m)
    df = dofadjustment(knew, cache)
    return _vcovmatrix(knew, m, cache, scale, factorization)
end

function _vcov(k::CRHC, m::RegressionModel, cache::CRHCCache, scale::Real)
    B = invpseudohessian(m)
    res = adjustresid!(k, cache)
    cache.momentmatrix .= cache.modelmatrix.*res
    df = dofadjustment(k, cache)
    scale *= df
    Shat = clusterize!(cache)
    return Symmetric((B*Shat*B).*scale)
end

function _vcovmatrix(
    k::HC,
    m::RegressionModel,
    cache::CRHCCache,
    scale::Real,
    ::Type{Cholesky}
)
    V = _vcov(k, m, cache, scale)
    CovarianceMatrix(cholesky(V, check=true), k, V)
end

function _vcovmatrix(
    k::HC,
    m::RegressionModel,
    cache::CRHCCache,
    scale::Real,
    ::Type{SVD}
)
    V = _vcov(k, m, cache, cache)
    CovarianceMatrix(svd(V.data), k, V)
end

# -----------------------------------------------------------------------------
# CRHC GLM - Trick to use vcov(CRHC1(:cluster, df), ::GLM)
# -----------------------------------------------------------------------------
StatsBase.stderror(k::RobustVariance, m::RegressionModel; kwargs...) = sqrt.(diag(vcov(k, m; kwargs...)))
StatsBase.stderror(v::CovarianceMatrix) = sqrt.(diag(v.V))

# -----------------------------------------------------------------------------
# CRHC GLM - Trick to use vcov(CRHC1(:cluster, df), ::GLM)
# -----------------------------------------------------------------------------
recast(k::CRHC{T,D}, m::INNERMOD) where {T<:AbstractVector, D<:Nothing} = k
recast(k::CRHC{T,D}, m::TableRegressionModel) where {T<:AbstractVector, D<:Nothing} = k
reterm(k::CRHC{T,D}, m::TableRegressionModel) where {T, D} = tuple(k.cl...)

function recast(k::CRHC{T,D}, m::TableRegressionModel) where {T<:Symbol, D}
    @assert istable(k.df) "`df` must be a DataFrames"
    t = k.cl
    if length(k.df[!, t]) == length(m.mf.data[1])
        ## The dimension fit
        id = compress(categorical(k.df[:,t]))
        return renew(k, id)
    else
        f = m.mf.f
        frm = f.lhs ~ tuple(f.rhs.terms..., Term(t))
        nt = NamedTuple{tuple(StatsModels.termvars(frm)...)}(columntable(k.df))
        idx = StatsModels.missing_omit(nt)[2]
        id = compress(categorical(k.df[idx,t]))
        return renew(k, id)
    end
    # ct = columntable(clus)
    # length_unique = map(x->length(unique(x)), ct)
    # fg = 1:prod(length_unique)
    # #cl = map(x->compress(categorical(x)), eachcol(x))
    # clus[!, :clusid] .= size(clus, 2) > 1 ? zero(Int) : clus[!, tterms[1]]
    # if length(tterms) > 1
    #     for (i,j) in enumerate(groupby(clus, [tterms...]))
    #         j[:, :clusid] .= fg[i]
    #     end
    # end
    # id = compress(categorical(clus[!, :clusid]))
end

# -----------------------------------------------------------------------------
# optimalbandwidth method
# -----------------------------------------------------------------------------
function optimalbandwidth(
    k::HAC,
    m::TableRegressionModel{F};
    kwargs...
) where F<:INNERMOD
    setupkernelweights!(k, m)
    optimalbandwidth(k, m.model; kwargs...)
end

function optimalbandwidth(k::HAC, m::F; prewhite=false) where F<:INNERMOD
    mm = momentmatrix(m)
    setupkernelweights!(k, m)
    mmm, D = prewhiter(mm, Val{prewhite})
    return optimalbandwidth(k, mmm; prewhite=prewhite)
end

function setupkernelweights!(k, m::TableRegressionModel{T}) where T<:INNERMOD
    β = coef(m)
    resize!(k.weights, length(β))
    "(Intercept)" ∈ coefnames(m) ? (k.weights .= 1.0; k.weights[1] = 0.0) : k.weights .= 1.0
    return nothing
end

function setupkernelweights!(k, m::T) where T<:INNERMOD
    cf = coef(m)
    resize!(k.weights, length(cf))
    fill!(k.weights, 1)
    i = interceptindex(m)
    i !== nothing && (k.weights[i] = 0)
    return nothing
end


text:
include("../ds.jl")

struct compact_mnm_hyper <: distribution_hyper_params
    α::AbstractArray{Float64,1}
end

mutable struct compact_mnm_sufficient_statistics <: sufficient_statistics
    N::Float64
    points_sum::AbstractArray{Int64,1}
end

function calc_posterior(prior:: compact_mnm_hyper, suff_statistics::compact_mnm_sufficient_statistics)
    if suff_statistics.N == 0
        return prior
    end
    return compact_mnm_hyper(prior.α + suff_statistics.points_sum)
end

function sample_distribution(hyperparams::compact_mnm_hyper)
    cat_dist = rand(Dirichlet(hyperparams.α))
    return compact_mnm_dist(log.(cat_dist))
end

function create_sufficient_statistics(hyper::compact_mnm_hyper,posterior::compact_mnm_hyper,points::AbstractArray{Float64,2}, pts_to_group = 0)
    if length(points) == 0
        return compact_mnm_sufficient_statistics(size(points,2),zeros(Int64,length(hyper.α)))
    end
    pt_count = counts(Int.(points),length(hyper.α))
    return compact_mnm_sufficient_statistics(size(points,2),pt_count)
end

function log_multivariate_gamma(x::Number, D::Number)
    res::Float64 = D*(D-1)/4*log(pi)
    for d = 1:D
        res += logabsgamma(x+(1-d)/2)[1]
    end
    return res
end

function log_marginal_likelihood(hyper::compact_mnm_hyper, posterior_hyper::compact_mnm_hyper, suff_stats::compact_mnm_sufficient_statistics)
    D = length(suff_stats.points_sum)
    logpi = log(pi)
    val = logabsgamma(sum(hyper.α))[1] -logabsgamma(sum(posterior_hyper.α))[1] + sum((x-> logabsgamma(x)[1]).(posterior_hyper.α) - (x-> logabsgamma(x)[1]).(hyper.α))
    return val
end

function aggregate_suff_stats(suff_l::compact_mnm_sufficient_statistics, suff_r::compact_mnm_sufficient_statistics)
    return compact_mnm_sufficient_statistics(suff_l.N+suff_r.N, suff_l.points_sum + suff_r.points_sum)
end


text:
using Profile
using Random
using HTTP
using Cascadia
using Gumbo

# This is a sample crawler which crawls news data from TasnimNews Agency.

function fetch_urls(tag, date)
    base_url = "https://www.tasnimnews.com"
    url_params = "$base_url/fa/archive?service=$tag&sub=-1&date=$date"
    res = HTTP.request("GET", url_params)
    data = parsehtml(String(res.body))
    selector = Selector(".archive-result .content a")
    news_links = eachmatch(selector, data.root)
    links = []
    if size(news_links)[1] != 0 # Avoiding no news loop/errors with the given params
        for item in news_links
            push!(links, base_url * attrs(item)["href"])
        end
    end
    return links
end

function fetch_data(url)
    res = HTTP.request("GET", url)
    data = parsehtml(String(res.body))
    title_slc = Selector("h1.title")
    lead_slc = Selector("h3.lead")
    story_slc = Selector(".story")
    title = nodeText(eachmatch(title_slc, data.root)[1])
    spam = findfirst("صدای تسنیم", title)
    if spam != nothing
        return
    end
    lead = nodeText(eachmatch(lead_slc, data.root)[1])
    story = eachmatch(story_slc, data.root)
    story_part_slc = Selector("p")
    story_parts = map(item -> eachmatch(Selector("p"), item), story)
    full_news = ""
    println("fetching $title")
    for item in story_parts
        for part in item
            full_news *= " " * nodeText(eachmatch(Selector("p"), part)[1])
        end
    end
    title*",", lead*",", full_news

end

function fetch_news(t, d)
    # fectched_news = Dict()
    fectched_news = []
    y, m, d = d
    tag = t
    @time for i in y
        for j in m
            for k in d
                push!(fectched_news, map(url -> fetch_data(url), fetch_urls(tag, "$i%2F$j%2F$k")))
            end
        end
    end
    return fectched_news
end

function save_news(t, d, f_name)
    news_data = fetch_news(t, d)
    f = "$f_name.txt"
    count = 0
    open(f, "w") do file
        for day_news in news_data
            for news in day_news
                count = count+1
                if news != nothing
                    map(j -> write(file, j), news) # Write title, lead and story
                    write(file, "\n")
                end
            end
        end
        println("Finished! Wrote $count news to $f_name.txt")
    end
end

# Tags guide: politics("1"), sports("3"), culture("4"), economics("7")
tags = ["1", "3", "4", "7"]
date = (1399, 3, 1)
file_name = ["politics", "sports", "culture", "economics"]
save_news(tags[1], date, file_name[1])


text:
"""
Intervals (Rvls) have two boundaries, each one has *value* and *reach*
*value* is a Float represented quantity.
*reach* is Closed or Open.
"""

abstract AbstractInterval       <: Real
abstract Boundary               <: AbstractInterval
abstract BoundingValue             <: Boundary       # reserved, not used
abstract Grasp                     <: Boundary       # configuration of bounding extrema
                                                     # e.g. ClOp for lobound Closed, hibound Open)

type Self  <: Grasp end

type ClCl  <: Grasp end
type ClOp  <: Grasp end
type OpCl  <: Grasp end
type OpOp  <: Grasp end

# Rvl --> Intrvl, Irvl --> Rvl (ravel)

immutable Rvl{G<:Grasp, R<:Real} <: Real
    lo::R
    hi::R
end

lowerbound{G<:Grasp, R<:Real}(x::Rvl{G,R}) = x.lo
upperbound{G<:Grasp, R<:Real}(x::Rvl{G,R}) = x.hi

#=
aClCl = Rvl{ClCl,Float64}(1.0,2.0); aClCl == ClCl(1.0,2.0)
aOpCl = Rvl{OpCl,Float64}(1.0,2.0); aOpCl == OpCl(1.0,2.0)
aClOp = Rvl{ClOp,Float64}(1.0,2.0); aClOp == ClOp(1.0,2.0)
aOpOp = Rvl{OpOp,Float64}(1.0,2.0); aOpOp == OpOp(1.0,2.0)
=#

for B in (:Self, :ClCl, :ClOp, :OpCl, :OpOp)
  @eval begin
     function convert{R<:Real}(::Type{Rvl{($B),R}}, lo::R, hi::R)
        lo, hi = minmax(lo, hi)
        Rvl{$B,R}(lo,hi)
     end
     convert{R<:Real}(::Type{Rvl{($B),R}}, x::R) = convert(Rvl($B,R), x,x)
     function ($B){R<:Real}(lo::R, hi::R)
        lo, hi = minmax(lo, hi)
        Rvl{$B,R}(lo,hi)
     end
     function ($B){R1<:Real,R2<:Real}(lo::R1, hi::R2)
        ($B)(promote(lo,hi)...)
     end
     function ($B){R<:Real}(x::R)
        Rvl{$B,R}(x,x)
     end
  end
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the Lattice-Boltzmann equations for the D2Q9 scheme

equations = LatticeBoltzmannEquations2D(Ma=0.1, Re=1000)

"""
    initial_condition_lid_driven_cavity(x, t, equations::LatticeBoltzmannEquations2D)

Initial state for a lid-driven cavity flow setup. To be used in combination with
[`boundary_condition_lid_driven_cavity`](@ref) and [`boundary_condition_noslip_wall`](@ref).
"""
function initial_condition_lid_driven_cavity(x, t, equations::LatticeBoltzmannEquations2D)
  @unpack L, u0, nu = equations

  rho = 1
  v1 = 0
  v2 = 0

  return equilibrium_distribution(rho, v1, v2, equations)
end
initial_condition = initial_condition_lid_driven_cavity

"""
    boundary_condition_lid_driven_cavity(u_inner, orientation, direction, x, t,
                                         surface_flux_function,
                                         equations::LatticeBoltzmannEquations2D)

Boundary condition for a lid-driven cavity flow setup, where the top lid (+y boundary) is a moving
no-slip wall. To be used in combination with [`initial_condition_lid_driven_cavity`](@ref).
"""
function boundary_condition_lid_driven_cavity(u_inner, orientation, direction, x, t,
                                              surface_flux_function,
                                              equations::LatticeBoltzmannEquations2D)
  return boundary_condition_moving_wall_ypos(u_inner, orientation, direction, x, t,
                                             surface_flux_function, equations)
end

function boundary_condition_moving_wall_ypos(u_inner, orientation, direction, x, t,
                                             surface_flux_function,
                                             equations::LatticeBoltzmannEquations2D)
  @assert direction == 4 "moving wall assumed in +y direction"

  @unpack rho0, u0, weights, c_s = equations
  cs_squared = c_s^2

  pdf1 = u_inner[3] + 2 * weights[1] * rho0 * u0 / cs_squared
  pdf2 = u_inner[2] # outgoing
  pdf3 = u_inner[1] + 2 * weights[3] * rho0 * (-u0) / cs_squared
  pdf4 = u_inner[2]
  pdf5 = u_inner[5] # outgoing
  pdf6 = u_inner[6] # outgoing
  pdf7 = u_inner[5] + 2 * weights[7] * rho0 * (-u0) / cs_squared
  pdf8 = u_inner[6] + 2 * weights[8] * rho0 * u0 / cs_squared
  pdf9 = u_inner[9]

  u_boundary = SVector(pdf1, pdf2, pdf3, pdf4, pdf5, pdf6, pdf7, pdf8, pdf9)

  # Calculate boundary flux (u_inner is "left" of boundary, u_boundary is "right" of boundary)
  return surface_flux_function(u_inner, u_boundary, orientation, equations)
end
boundary_conditions = (
                       x_neg=boundary_condition_noslip_wall,
                       x_pos=boundary_condition_noslip_wall,
                       y_neg=boundary_condition_noslip_wall,
                       y_pos=boundary_condition_lid_driven_cavity,
                      )

solver = DGSEM(polydeg=5, surface_flux=flux_godunov)

coordinates_min = (0.0, 0.0)
coordinates_max = (1.0, 1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=4,
                periodicity=false,
                n_cells_max=10_000,)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 400.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=1000,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2macroscopic)

stepsize_callback = StepsizeCallback(cfl=1.0)

collision_callback = LBMCollisionCallback()

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback,
                        collision_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
## Load Packages
using Omega
using StatsBase
using Random
using Plots
using Distributions

Random.seed!(123)

## Define functions
function sample_weights(rng, hazards)
    result = sample(collect(keys(hazards)), Weights(collect(values(hazards))))
    return result
end

function initial_simulation(rng, t, transitions)

    """
    Simulates first run of Gillespie. Needs to be different because
    it takes in the intial prey, pred and rates.
    """

    ecology = Dict("prey" => prey_init(rng), "pred" => pred_init(rng))
    # theta = Dict("spawn_prey" => spawn_prey(rng),
    #              "prey2pred" => prey2pred(rng),
    #              "pred_dies" => pred_dies(rng))

    hazards = get_hazards_fix(ecology, theta)

    # Does this make sense?
    hazard_sample = ciid(sample_weights, hazards)
    transition = transitions[hazard_sample(rng)]
    t_new = t + sum(values(hazards))

    ecology_new_prey = ecology["prey"] + transition[1]
    ecology_new_pred = ecology["pred"] + transition[2]
    # Enforce only positive integers
    ecology_new_prey = max(1, ecology_new_prey)
    ecology_new_pred = max(1, ecology_new_pred)

    return Dict("prey" => ecology_new_prey,
                "pred" => ecology_new_pred,
                "t" => t_new)
end

function one_simulation(rng, n, theta, transitions)

    """
    Simulates all steps of gillespie except intial step. Takes into account
    the results of the previous step "n".
    """

    ecology = functions_list[n](rng)
    ecology_temp = copy(ecology)
    t = ecology["t"]
    delete!(ecology_temp, "t")
    hazards = get_hazards_fix(ecology_temp, theta)
    # Does this make sense?
    #hazard_sample = ciid(sample_weights, hazards)
    transition = transitions[sample(collect(keys(hazards)), Weights(collect(values(hazards))))]
    t_new = t + sum(values(hazards))

    ecology_new_prey = ecology_temp["prey"] + transition[1]
    ecology_new_pred = ecology_temp["pred"] + transition[2]
    # Enforce only positive integers
    ecology_new_prey = max(1, ecology_new_prey)
    ecology_new_pred = max(1, ecology_new_pred)

    return Dict("prey" => ecology_new_prey,
                "pred" => ecology_new_pred,
                "t" => t_new)
end

function get_hazards(ecology, theta)
    """
    Compute the hazard function given the current states.  Note that in this
        model the function depends only on state, not time, though this is not
        the case in general. "spawn_prey" represents the event of a prey being born,
        "prey2pred" represents a predator consuming a new prey and consequently spawning
        a new predator, "pred_dies" represents the death of a predator.

    args:
        ecology(dict): a dictionary containing species counts.
        theta(dict): A dictionary where keys represent events (reactions) and
            values are the location hyperparameters for rate parameter constants
            corresponding to those events.
    """

    return Dict(
        "spawn_prey" => theta["spawn_prey"] * ecology["prey"],
        "prey2pred" => theta["prey2pred"] * ecology["prey"] * ecology["pred"],
        "pred_dies" => theta["pred_dies"] * ecology["pred"]
        )
end

function get_hazards_fix(ecology, theta)
    """
    Same as above, for some reason rates become lists in subsequent runs so
    this fixes that. Juliaaaa issues 😫
    """
    return Dict(
        "spawn_prey" => theta["spawn_prey"][1] * ecology["prey"],
        "prey2pred" => theta["prey2pred"][1] * ecology["prey"] * ecology["pred"],
        "pred_dies" => theta["pred_dies"][1] * ecology["pred"]
        )
end

## Run simulation
# Define gillespie parameters
Pre = [[1, 0], [1, 1], [0, 1]]
Post = [[2, 0], [0, 2], [0, 0]]
transition_mat = Post - Pre
transitions = Dict("spawn_prey" => transition_mat[1,],
                    "prey2pred" => transition_mat[2,],
                    "pred_dies" => transition_mat[3,])
t0=0.0
N = 10000

# Random variables for starting values
prey_init = normal(100, 10)
pred_init = normal(50, 10)

# Random variables for rates
spawn_prey = normal(1, .1)
prey2pred = normal(.01, .001)
pred_dies = normal(.5, .05)

theta = Dict("spawn_prey" => rand(spawn_prey,1),
             "prey2pred" => rand(prey2pred,1),
             "pred_dies" => rand(pred_dies,1))

# Random variables for each step of simulation
temp = ciid(initial_simulation, t0, transitions) # initial step
functions_list = Any[]
push!(functions_list, temp)

for f in 2:N
    last = f - 1
    temp = ciid(one_simulation, last, theta, transitions) # individual step
    push!(functions_list, temp)
end

# Add initialized random vars to list
push!(functions_list,prey_init)
push!(functions_list, pred_init)
push!(functions_list, spawn_prey)
push!(functions_list, prey2pred)
push!(functions_list, pred_dies)

# sample
samples = rand(Tuple(x for x in functions_list), functions_list[5000] >= 120,
                5, alg = RejectionSample)

# extract run results and plot
prey_vals = []
pred_vals = []
for x in 1:N
    push!(prey_vals,samples[3][x]["prey"])
    push!(pred_vals,samples[3][x]["pred"])
end

plot(hcat(prey_vals,pred_vals))

# sample
samples = rand(Tuple(x for x in functions_list),prey_init > 120
                , alg = RejectionSample)

samples = rand(Tuple(x for x in functions_list),functions_list[5000] > 120
                , alg = RejectionSample)

# extract run results and plot
prey_vals = []
pred_vals = []
for x in 1:N
    push!(prey_vals,samples[x]["prey"])
    push!(pred_vals,samples[x]["pred"])
end

plot(hcat(prey_vals, pred_vals))


text:
let
    n = 50
    p = 3
    t = (1:n)/n
    Ut,Vt = egrss.generators(t,p);
    K = egrss.full(Ut,Vt)

    Wt = egrss.potrf(Ut,Vt)
    Lref = cholesky(K).L
    L = egrss.full_tril(Ut,Wt)
    @test( isapprox(L,Lref) )

    alpha = 1e-6
    Wt,c = egrss.potrf(Ut,Vt,alpha)
    Lref = cholesky(K + Diagonal(alpha*ones(n))).L
    L = egrss.full_tril(Ut,Wt,c)
    @test( isapprox(L,Lref) )

    d = collect((1:n)/n) .+ 1e-3
    Wt,c = egrss.potrf(Ut,Vt,d)
    Lref = cholesky(K + Diagonal(d)).L
    L = egrss.full_tril(Ut,Wt,c)
    @test( isapprox(L,Lref) )

end

max_stars_repo_name: martinandersen/EGRSS

max_issues_repo_name: martinandersen/EGRSS

max_forks_repo_name: martinandersen/EGRSS


text:
############################################################################################
"""
$(TYPEDEF)

SMC Diagnostics container, including diagnostics of kernels used in jittering step.

# Fields
$(TYPEDFIELDS)
"""
struct SMCDiagnostics{P,J,T<:AbstractFloat,G} <: AbstractDiagnostics
    "Diagnostics used for all Baytes kernels"
    base::BaytesCore.BaseDiagnostics{Vector{P}}
    "Weighted average of incremental log weights - can be used for marginal likelihood computation."
    ℓincrement::Float64
    "Log weights for resampling steps."
    ℓweights::Vector{Float64}
    "Normalized Log weights."
    ℓweightsₙ::Vector{Float64}
    "Diagnostics of jitter steps. If not resampled this iteration (i.e., resampled == false), contains jitterdiagnostics from previous step."
    jitterdiagnostics::Vector{J}
    "Number of jittering steps"
    jittersteps::Int64
    "Correlation from rejuvented continuous parameter"
    ρ::Vector{T}
    "ESS and accepted steps for SMC kernel"
    ESS::Float64
    "Boolean if step has been resampled."
    resampled::Bool
    "Generated quantities specified for objective"
    generated::G
    function SMCDiagnostics(
        base::BaytesCore.BaseDiagnostics{Vector{P}},
        ℓincrement::Float64,
        ℓweights::Vector{Float64},
        ℓweightsₙ::Vector{Float64},
        jitterdiagnostics::Vector{J},
        jittersteps::Int64,
        ρ::Vector{T},
        ESS::Float64,
        resampled::Bool,
        generated::G
    ) where {P,J,T<:AbstractFloat,G}
        return new{P,J,T,G}(
            base,
            ℓincrement,
            ℓweights,
            ℓweightsₙ,
            jitterdiagnostics,
            jittersteps,
            ρ,
            ESS,
            resampled,
            generated
        )
    end
end

############################################################################################
function generate_showvalues(diagnostics::D) where {D<:SMCDiagnostics}
    return function showvalues()
        return (:smc, "diagnostics"),
        (:iter, diagnostics.base.iter),
        (:Avgℓobjective, mean(diagnostics.ℓweights)),
        (:Temperature, diagnostics.base.temperature),
        (:ESS, diagnostics.ESS),
        (:resampled, diagnostics.resampled),
        (:AvgJitterCorrelation, mean(diagnostics.ρ))
    end
end

############################################################################################
#export
export SMCDiagnostics, generate_showvalues


text:
using Blink
import GeometryTypes

notinstalled = !AtomShell.isinstalled()
notinstalled && AtomShell.install()

window = Window()
vis = Visualizer()

if !haskey(ENV, "CI")
    open(vis)
end

if !(Sys.iswindows() && haskey(ENV, "CI"))
    # this gets stuck on windows CI, but I don't know why
    open(vis, window)
end

if !haskey(ENV, "CI")
    wait(vis)
end

# A custom geometry type to test that we can render arbitrary primitives
# by decomposing them into simple meshes. This replaces the previous test
# which did the same thing using a Polyhedron from Polyhedra.jl. The Polyhedra
# test was removed because it required too many external dependencies just to
# verify a simple interface.
struct CustomGeometry <: GeometryPrimitive{3, Float64}
end

GeometryTypes.isdecomposable(::Type{<:Face}, ::CustomGeometry) = true
function GeometryTypes.decompose(::Type{F}, c::CustomGeometry) where {F <: Face}
    [convert(F, Face(1,2,3))]
end
GeometryTypes.isdecomposable(::Type{<:Point}, ::CustomGeometry) = true
function GeometryTypes.decompose(::Type{P}, c::CustomGeometry) where {P <: Point}
    convert(Vector{P}, [Point(0., 0, 0), Point(0., 1, 0), Point(0., 0, 1)])
end

@testset "self-contained visualizer" begin
    cat_mesh_path = joinpath(dirname(pathof(GeometryTypes)), "..", "test", "data", "cat.obj")

    @testset "shapes" begin
        v = vis[:shapes]
        delete!(v)
        settransform!(v, Translation(1., 0, 0))
        @testset "box" begin
            setobject!(v[:box], HyperRectangle(Vec(0., 0, 0), Vec(0.1, 0.2, 0.3)))
            settransform!(v[:box], Translation(-0.05, -0.1, 0))
        end

        @testset "cylinder" begin
            setobject!(v[:cylinder], Mesh(
               Cylinder(Point(0., 0, 0), Point(0, 0, 0.2), 0.1),
               MeshLambertMaterial(color=colorant"lime")))
            settransform!(v[:cylinder], Translation(0, 0.5, 0.0))
        end

        @testset "sphere" begin
            setobject!(v[:sphere],
                HyperSphere(Point(0., 0, 0), 0.15),
                MeshLambertMaterial(color=colorant"maroon"))
            settransform!(v[:sphere], Translation(0, 1, 0.15))
        end

        @testset "ellipsoid" begin
            setobject!(v[:ellipsoid], HyperEllipsoid(Point(0., 1.5, 0), Vec(0.3, 0.1, 0.1)))
            settransform!(v[:ellipsoid], Translation(0, 0, 0.1))
        end

        @testset "cube" begin
            setobject!(v[:cube], HyperCube(Vec(-0.1, -0.1, 0), 0.2), MeshBasicMaterial())
            settransform!(v[:cube], Translation(0, 2.0, 0))
        end

        @testset "more complicated cylinder" begin
            setobject!(v[:cylinder2], Cylinder(Point(0, 2.5, 0), Point(0.1, 2.6, 0), 0.05))
            settransform!(v[:cylinder2], Translation(0, 0, 0.05))
        end

        @testset "triad" begin
            setobject!(v[:triad], Triad(0.2))
            settransform!(v[:triad], Translation(0, 3, 0.2))
        end

        @testset "cone" begin
            setobject!(v[:cone],
                Cone(Point(1., 1., 1.), Point(1., 1., 1.2), 0.1),
                MeshLambertMaterial(color=colorant"indianred"))
            settransform!(v[:cone], Translation(-1, 2.5, -1))
        end
    end

    @testset "meshes" begin
        v = vis[:meshes]
        @testset "cat" begin
            mesh = load(cat_mesh_path)
            setobject!(v[:cat], mesh)
            settransform!(v[:cat], Translation(0, -1, 0) ∘ LinearMap(RotZ(π)) ∘ LinearMap(RotX(π/2)))
        end

        @testset "cat_color" begin
            mesh = load(cat_mesh_path)
            color = RGBA{Float32}(0.5, 0.5, 0.5, 0.5)
            mesh_color = HomogenousMesh(vertices=mesh.vertices, faces=mesh.faces, color=color)
            object = Object(mesh_color)
            @test MeshCat.material(object).color == color
            mesh_color = setobject!(v[:cat_color], mesh_color)
            settransform!(v[:cat_color], Translation(0, -2.0, 0) ∘ LinearMap(RotZ(π)) ∘ LinearMap(RotX(π/2)))
        end

        @testset "mesh file geometries" begin
            settransform!(v[:mesh_file_geometries], Translation(0, -3.0, 0))
            base_path = joinpath(@__DIR__, "data", "meshes", "mesh_0_convex_piece_0")
            @testset "obj" begin
                path = base_path * ".obj"
                setobject!(v[:mesh_file_geometries, :obj], MeshFileGeometry(path))
            end
            @testset "dae" begin
                path = base_path * ".dae"
                setobject!(v[:mesh_file_geometries, :dae], MeshFileGeometry(path))
                settransform!(v[:mesh_file_geometries, :dae], Translation(0, -0.5, 0))
            end
            @testset "stl_ascii" begin
                path = base_path * ".ascii.stl"
                setobject!(v[:mesh_file_geometries, :stl_ascii], MeshFileGeometry(path))
                settransform!(v[:mesh_file_geometries, :stl_ascii], Translation(0, -1.0, 0))
            end
            @testset "stl_binary" begin
                path = base_path * ".binary.stl"
                setobject!(v[:mesh_file_geometries, :stl_binary], MeshFileGeometry(path))
                settransform!(v[:mesh_file_geometries, :stl_binary], Translation(0, -1.5, 0))
            end
        end

        @testset "mesh file objects" begin
            let v = v[:mesh_file_objects]
                settransform!(v, Translation(0, 3.0, -1.5))

                base_path = joinpath(@__DIR__, "data", "meshes")
                @testset "obj" begin
                    path = joinpath(base_path, "cube.obj")
                    setobject!(v[:obj], MeshFileObject(path))
                end

                @testset "dae" begin
                    path = joinpath(base_path, "cube.dae")
                    setobject!(v[:dae], MeshFileObject(path))
                    settransform!(v[:dae], Translation(0, -1.5, 0))
                end
            end
        end

        @testset "textured valkyrie" begin
            head = Mesh(
                load(joinpath(MeshCat.VIEWER_ROOT, "..", "data", "head_multisense.obj"), GLUVMesh),
                MeshLambertMaterial(
                    map=Texture(
                        image=PngImage(joinpath(MeshCat.VIEWER_ROOT, "..", "data", "HeadTextureMultisense.png"))
                    )
                ))
            setobject!(v[:valkyrie, :head], head)
            settransform!(v[:valkyrie, :head], Translation(0, 0.5, 0.5))
        end

    end

    @testset "points" begin
        v = vis[:points]
        settransform!(v, Translation(-1, 0, 0))
        @testset "random points" begin
            verts = rand(Point3f0, 100000);
            colors = reinterpret(RGB{Float32}, verts);
            setobject!(v[:random], PointCloud(verts, colors))
            settransform!(v[:random], Translation(-0.5, -0.5, 0))
        end
    end

    @testset "points with material (Issue #58)" begin
        v = vis[:points_with_material]
        settransform!(v, Translation(-1.5, -2.5, 0))
        material = PointsMaterial(color=RGBA(0,0,1,0.5))
        cloud = PointCloud(rand(Point3f0, 5000))
        obj = Object(cloud, material)
        @test MeshCat.threejs_type(obj) == "Points"
        setobject!(v, cloud, material)
    end

    @testset "lines" begin
        v = vis[:lines]
        settransform!(v, Translation(-1, -1, 0))
        @testset "LineSegments" begin
            θ = range(0, stop=2π, length=10)
            setobject!(v[:line_segments], LineSegments(Point.(0.5 .* sin.(θ), 0, 0.5 .* cos.(θ))))
        end
        @testset "Line" begin
            θ = range(0, stop=2π, length=10)
            setobject!(v[:line], Line(Point.(0.5 .* sin.(θ), 0, 0.5 .* cos.(θ))))
            settransform!(v[:line], Translation(0, 0.1, 0))
        end
        @testset "LineLoop" begin
            θ = range(0, stop=π, length=10)
            setobject!(v[:line_loop], LineLoop(Point.(0.5 .* sin.(θ), 0, 0.5 .* cos.(θ))))
            settransform!(v[:line_loop], Translation(0, 0.2, 0))
        end
    end

    @testset "Custom geometry primitives" begin
        primitive = CustomGeometry()
        setobject!(vis[:custom], primitive)
        settransform!(vis[:custom],  Translation(-0.5, 1.0, 0) ∘ LinearMap(UniformScaling(0.5)))
    end

    @testset "ArrowVisualizer" begin
        arrow_vis_1 = ArrowVisualizer(vis[:arrow1])
        show(devnull, arrow_vis_1)
        setobject!(arrow_vis_1)
        for l in [range(0, stop=1e-2, length=1000); 1.1 * eps(Float64)]
            settransform!(arrow_vis_1, Point(0, 1, 0), Vec(0, 0, l))
        end
        settransform!(arrow_vis_1, Point(0, 1, 0), Vec(0, 0, 1.1 * eps(Float64)))
        settransform!(arrow_vis_1, Point(0, 1, 0), Vec(1, 1, 1))
        setobject!(vis[:arrow1_base], HyperSphere(Point(0., 1., 0.), 0.015))
        setobject!(vis[:arrow1_tip], HyperSphere(Point(1., 2., 1.), 0.015))

        arrow_vis_2 = ArrowVisualizer(vis[:arrow2])
        setobject!(arrow_vis_2;
            shaft_material=MeshLambertMaterial(color=colorant"lime"),
            head_material=MeshLambertMaterial(color=colorant"indianred"))
        settransform!(arrow_vis_2, Point(0, 1, 0), Point(1, 1, 1))
    end

    @testset "Animation (old style)" begin
        anim = Animation()
        atframe(anim, vis[:shapes], 0) do frame_vis
            settransform!(frame_vis[:box], Translation(0., 0, 0))
        end
        atframe(anim, vis[:shapes], 30) do frame_vis
            settransform!(frame_vis[:box], Translation(2., 0, 0) ∘ LinearMap(RotZ(π/2)))
        end
        atframe(anim, vis, 0) do framevis
            setprop!(framevis["/Cameras/default/rotated/<object>"], "zoom", 1)
        end
        atframe(anim, vis, 30) do framevis
            setprop!(framevis["/Cameras/default/rotated/<object>"], "zoom", 0.5)
        end
        setanimation!(vis, anim)
    end

    @testset "Animation (new style)" begin
        anim1 = Animation()
        atframe(anim1, 0) do
            settransform!(vis[:shapes][:box], Translation(0., 0, 0))
        end
        atframe(anim1, 30) do
            settransform!(vis[:shapes][:box], Translation(2., 0, 0) ∘ LinearMap(RotZ(π/2)))
        end
        setanimation!(vis, anim1)
        anim2 = Animation()
        atframe(anim2, 0) do
            setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 1)
        end
        atframe(anim2, 30) do
            setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 0.5)
        end
        setanimation!(vis, anim2)
        anim_combined = merge(anim1, anim2)
        @test Set(Iterators.flatten((keys(anim1.clips), keys(anim2.clips)))) == Set(keys(anim_combined.clips))
        setanimation!(vis, anim_combined)
    end
end

@testset "AnimationTrack" begin
    track1 = MeshCat.AnimationTrack{Float64}("foo", "bar")
    @test track1.name == "foo"
    @test track1.jstype == "bar"
    insert!(track1, 0, 32.0)
    @test track1.events == [0 => 32.0]
    insert!(track1, 2, 64.0)
    @test track1.events == [0 => 32.0, 2 => 64.0]
    insert!(track1, 2, 65.0)
    @test track1.events == [0 => 32.0, 2 => 65.0]

    track2 = MeshCat.AnimationTrack{Float64}("foo", "bar")
    insert!(track2, 1, 17.0)
    insert!(track2, 2, 66.0)
    insert!(track2, 5, 1.0)

    merge!(track1, track2)
    @test track1.events == [0 => 32.0, 1 => 17.0, 2 => 66.0, 5 => 1.0]
end

@testset "setvisible!" begin
    v = vis[:box_to_hide]
    setobject!(v, HyperRectangle(Vec(0., 0, 0), Vec(0.1, 0.2, 0.3)))
    sleep(1)
    setvisible!(v, false)
    sleep(1)
    setvisible!(v, true)
end

sleep(5)

if !(Sys.iswindows() && haskey(ENV, "CI"))
    # this also fails on appveyor, and again I have no way to debug it
    notinstalled && AtomShell.uninstall()
end


text:
abstract type AbstractPrimalHeuristic end
abstract type RoundingHeuristic <: AbstractPrimalHeuristic end
abstract type AllBlockHeuristic <: AbstractPrimalHeuristic end

function add!(htype::Type{T}, LD::AbstractLagrangeDual) where T <: AbstractPrimalHeuristic
    push!(LD.heuristics, htype)
end

"""
    run!

This run a rounding heuristic.
"""
function run!(::Type{RoundingHeuristic}, LD::AbstractLagrangeDual, val, ub, lb)
    #set weights of each model if none exists
    num_all_blocks = parallel.sum(num_blocks(LD))
    all_block_ids = Set([ks[1] for (ks, i) in LD.var_to_index])
    if num_all_blocks != LD.block_model.combined_weights.count
        LD.block_model.combined_weights = Dict(block_id => 1/num_all_blocks for block_id in all_block_ids)
    end

    # get the mean value of the coupling variables at the root process and broadcast to other processes
    new_primal_solution = Dict() #maps coupling_id to value
    if parallel.is_root()
        for (ks, i) in LD.var_to_index
            coupling_id = ks[2]
            block_id = ks[1]
            new_primal_solution[coupling_id] = haskey(new_primal_solution, coupling_id) ? new_primal_solution[coupling_id] + LD.block_model.combined_weights[block_id] * val[block_id][i] : LD.block_model.combined_weights[block_id] * val[block_id][i]
        end
        parallel.bcast(new_primal_solution)
    else
        new_primal_solution = parallel.bcast(nothing)
    end

    #fix variables with new primal solution and enforce integrality through rounding
    for variables in LD.block_model.coupling_variables
        if JuMP.is_integer(variables.ref) || JuMP.is_binary(variables.ref)
            JuMP.fix(variables.ref, round(new_primal_solution[variables.key.coupling_id]), force=true)
            new_primal_solution[variables.key.coupling_id] = round(new_primal_solution[variables.key.coupling_id])
        else
            JuMP.fix(variables.ref, new_primal_solution[variables.key.coupling_id], force=true)
            new_primal_solution[variables.key.coupling_id] = new_primal_solution[variables.key.coupling_id]
        end
    end

    #obtain primal bound by solving the subproblems in parallel
    cur_primal_bound = 0.0
    for (id,m) in block_model(LD)
        JuMP.optimize!(m)
        if ! (JuMP.termination_status(m)  in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED])
            cur_primal_bound = + Inf
            break
        else
            cur_primal_bound += JuMP.objective_value(m)
        end
    end
    cur_primal_bound_sum = parallel.sum(cur_primal_bound)

    #update primal bound and solution
    if cur_primal_bound_sum < LD.block_model.primal_bound
        LD.block_model.primal_bound = cur_primal_bound_sum
        LD.block_model.primal_solution = copy(new_primal_solution)
        if parallel.is_root()
            println("  Found new best primal bound: $(LD.block_model.primal_bound)")
        end
    end

    #unfix variables and recover their original bounds
    for var in coupling_variables(LD)
        JuMP.unfix(var.ref)
        JuMP.set_lower_bound(var.ref, lb[var.key.block_id][index_of_λ(LD, var)])
        JuMP.set_upper_bound(var.ref, ub[var.key.block_id][index_of_λ(LD, var)])
    end
end

"""
    run!

This run a heuristic that simply fixes the coupling variable values.
"""
function run!(::Type{AllBlockHeuristic}, LD::AbstractLagrangeDual, val, ub, lb)
    all_block_ids = Set([ks[1] for (ks, i) in LD.var_to_index])
    all_coupling_ids = Set([ks[2] for (ks, i) in LD.var_to_index])

    # broadcast the val
    if parallel.is_root()
        parallel.bcast(val)
    else
        val = parallel.bcast(nothing)
    end

    #iterate over all blocks. Fix the coupling variables to the optimal solution of each block.
    for block_id in all_block_ids
        new_primal_solution = Dict() #maps coupling_id to value
        for coupling_id in all_coupling_ids
            new_primal_solution[coupling_id] = val[block_id][LD.var_to_index[block_id, coupling_id]]
        end

        #fix variables with new primal solution and enforce integrality
        for variables in LD.block_model.coupling_variables
            JuMP.fix(variables.ref, new_primal_solution[variables.key.coupling_id], force=true)
        end

        #obtain primal bound by solving the subproblems in parallel
        cur_primal_bound = 0.0
        for (id,m) in block_model(LD)
            JuMP.optimize!(m)
            if ! (JuMP.termination_status(m)  in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED])
                cur_primal_bound = + Inf
                break
            else
                cur_primal_bound += JuMP.objective_value(m)
            end
        end

        #update bounds and solution
        cur_primal_bound_sum = parallel.sum(cur_primal_bound)
        if cur_primal_bound_sum < LD.block_model.primal_bound
            LD.block_model.primal_bound = cur_primal_bound_sum
            LD.block_model.primal_solution = copy(new_primal_solution)

            if parallel.is_root()
                println("  Found new best primal bound: $(LD.block_model.primal_bound)")
            end
        end
    end

    #unfix variables
    for var in coupling_variables(LD)
        JuMP.unfix(var.ref)
        JuMP.set_lower_bound(var.ref, lb[var.key.block_id][index_of_λ(LD, var)])
        JuMP.set_upper_bound(var.ref, ub[var.key.block_id][index_of_λ(LD, var)])
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayRewriteRule <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    ruleSequence::Any # spec type: Union{ Nothing, Int32 } # spec name: ruleSequence
    conditions::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayRewriteRuleCondition} } # spec name: conditions
    actionSet::Any # spec type: Union{ Nothing, ApplicationGatewayRewriteRuleActionSet } # spec name: actionSet

    function ApplicationGatewayRewriteRule(;name=nothing, ruleSequence=nothing, conditions=nothing, actionSet=nothing)
        o = new()
        validate_property(ApplicationGatewayRewriteRule, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(ApplicationGatewayRewriteRule, Symbol("ruleSequence"), ruleSequence)
        setfield!(o, Symbol("ruleSequence"), ruleSequence)
        validate_property(ApplicationGatewayRewriteRule, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(ApplicationGatewayRewriteRule, Symbol("actionSet"), actionSet)
        setfield!(o, Symbol("actionSet"), actionSet)
        o
    end
end # type ApplicationGatewayRewriteRule

const _property_map_ApplicationGatewayRewriteRule = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("ruleSequence")=>Symbol("ruleSequence"), Symbol("conditions")=>Symbol("conditions"), Symbol("actionSet")=>Symbol("actionSet"))
const _property_types_ApplicationGatewayRewriteRule = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("ruleSequence")=>"Int32", Symbol("conditions")=>"Vector{ApplicationGatewayRewriteRuleCondition}", Symbol("actionSet")=>"ApplicationGatewayRewriteRuleActionSet")
Base.propertynames(::Type{ ApplicationGatewayRewriteRule }) = collect(keys(_property_map_ApplicationGatewayRewriteRule))
Swagger.property_type(::Type{ ApplicationGatewayRewriteRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayRewriteRule[name]))}
Swagger.field_name(::Type{ ApplicationGatewayRewriteRule }, property_name::Symbol) =  _property_map_ApplicationGatewayRewriteRule[property_name]

function check_required(o::ApplicationGatewayRewriteRule)
    true
end

function validate_property(::Type{ ApplicationGatewayRewriteRule }, name::Symbol, val)
end


text:
# joOuterProd

## helper module
module joOuterProd_etc

    using JOLI: LocalVecOrMat

    @inline frwrd(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Complex,VDT<:Complex,xDT<:Number} = U*(adjoint(V)*x)
    @inline frwrd(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Complex,VDT<:Real,xDT<:Number} = U*(transpose(V)*x)
    @inline frwrd(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Real,VDT<:Complex,xDT<:Number} = U*(adjoint(V)*x)
    @inline frwrd(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Real,VDT<:Real,xDT<:Number} = U*(transpose(V)*x)

    @inline trnsp(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Complex,VDT<:Complex,xDT<:Number} = conj(V)*(transpose(U)*x)
    @inline trnsp(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Complex,VDT<:Real,xDT<:Number} = V*(transpose(U)*x)
    @inline trnsp(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Real,VDT<:Complex,xDT<:Number} = conj(V)*(transpose(U)*x)
    @inline trnsp(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Real,VDT<:Real,xDT<:Number} = V*(transpose(U)*x)

    @inline adjnt(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Complex,VDT<:Complex,xDT<:Number} = V*(adjoint(U)*x)
    @inline adjnt(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Complex,VDT<:Real,xDT<:Number} = V*(adjoint(U)*x)
    @inline adjnt(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Real,VDT<:Complex,xDT<:Number} = V*(transpose(U)*x)
    @inline adjnt(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Real,VDT<:Real,xDT<:Number} = V*(transpose(U)*x)

    @inline cnjgt(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Complex,VDT<:Complex,xDT<:Number} = conj(U)*(transpose(V)*x)
    @inline cnjgt(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Complex,VDT<:Real,xDT<:Number} = conj(U)*(transpose(V)*x)
    @inline cnjgt(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Real,VDT<:Complex,xDT<:Number} = U*(transpose(V)*x)
    @inline cnjgt(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT},x::LocalVecOrMat{xDT}) where
        {UDT<:Real,VDT<:Real,xDT<:Number} = U*(transpose(V)*x)

end
using .joOuterProd_etc

export joOuterProd
"""
    julia> op = joOuterProd(U,V;[DDT=...,][RDT=...,][name=...])

Memory efficient implementation of operator A = U*V'

# Signature

    joOuterProd(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT};
        DDT::DataType=joFloat,RDT::DataType=promote_type(UDT,VDT,DDT),name::String="joOuterProd")
            where {UDT<:Number,VDT<:Number}

# Arguments

- `U`: left vector or matrix
- `V`: right vector or matrix
- keywords
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Notes

- 2nd dimension of U and V must match

# Examples

    A=joOuterProd(rand(4),rand(5))
    A=joOuterProd(rand(4,2),rand(5,2))

examples with DDT/RDT

    A=joOuterProd(rand(4),rand(5); DDT=Float32)
    A=joOuterProd(rand(4),rand(5); DDT=Float32,RDT=Float64)

"""
function joOuterProd(U::LocalVecOrMat{UDT},V::LocalVecOrMat{VDT};
    DDT::DataType=joFloat,RDT::DataType=promote_type(UDT,VDT,DDT),name::String="joOuterProd") where {UDT<:Number,VDT<:Number}

    size(U,2)==size(V,2) || throw(joLinearFunctionException("joOuterProd: 2nd dimension of U & V must match"))
    m = size(U,1); n=size(V,1);

    return joLinearFunctionFwd(m,n,
        v1->jo_convert(RDT,joOuterProd_etc.frwrd(U,V,v1),false),
        v2->jo_convert(DDT,joOuterProd_etc.trnsp(U,V,v2),false),
        v3->jo_convert(DDT,joOuterProd_etc.adjnt(U,V,v3),false),
        v4->jo_convert(RDT,joOuterProd_etc.cnjgt(U,V,v4),false),
        DDT,RDT;
        name=name,
        fMVok=true
        )

end


text:
@testset "Widget" begin
    x = randn(Float32, 10, 2)
    y = rand(Float32, 1, 2)

    model = Dense(10, 1) |> Solaris.track

    function loss(x, y)
        xs = Flux.unstack(x, 3)
        ys = Flux.unstack(y, 3)
        ŷs = model.(xs)
        l = 0.0f0
        for t = 1:length(ŷs)
            l += Flux.mse(ys[t], ŷs[t])
        end
        return l / length(ŷs)
    end
    ps = Flux.params(model)
    data = repeat([(x, y)], 10)
    opt = ADAM()
    cb = () -> Flux.reset!(model)
    Solaris.tracker_mode()
    Flux.train!(loss, ps, data, opt, cb = cb)

    ps1 = Flux.params(model)
    Flux.train!(loss, ps1, data, opt, cb = cb)

    model |> Solaris.untrack
    Solaris.zygote_mode()
end


text:
module AKCalc

using Dates
using Tables
using DataFrames
using JSON

export @ak_str

# accessors
export resource_types, factory_recipes, workshop_recipes, operator_bases
export enemy_bases, enemies, tiles, maps, stages, stage_reports
export operators, resources

# AK types
export ResourceItemType, ResourceClassType, ResourceType, Resource
export FactoryRecipe, WorkshopRecipe
export AtkType, DmgType, Resistances
export PromotionPhase, OperatorBase, OperatorLevel, OperatorSkills, Operator
export EnemyFlags, EnemyRank, EnemyGrade, EnemyBase, EnemyStats
export Tile, Map, DropRarity, Stage
export StageReport, RecruitReport

# loaded early so the loops below can gather type instances
include("ak_types.jl")

for s in instances(ResourceItemType)
  @eval export $(Symbol(s))
end
for s in instances(ResourceClassType)
  @eval export $(Symbol(s))
end

export GDMem, GDKengxxiao, GDAceShip
export PDMem, PDSQLite
export SDMem

include("gamedata.jl")
include("playerdata.jl")
include("statsdata.jl")
include("queries.jl")
include("config.jl")
include("repl.jl")

using .GameData
using .PlayerData
using .StatsData
using .REPLMode: @ak_str

const DEFAULT_IO = Ref{IO}()

function __init__()
  DEFAULT_IO[] = stderr
  if isdefined(Base, :active_repl)
    REPLMode.repl_init(Base.active_repl)
  else
    atreplinit() do repl
      if isinteractive() && repl isa REPL.LineEditREPL
        isdefined(repl, :interface) || (repl.interface = REPL.setup_interface(repl))
        REPLMode.repl_init(repl)
      end
    end
  end
  return nothing
end

end


text:
M = 6;
c = 1e-5;
λ = 0.99

# Function
f(x) = sin.(x) .* x

# Training
x = [1, 3, 5, 6, 7, 8, 9, 10];
y = f(x);

# Test
X = collect(range(0, 10, length = 301));
Y = f(X);

# Kernel
K = RBF(10, 1);

#m_rls  = RLS(y, x, λ);
m_krls = KRLS(y, x, K, M, λ, c, "B2P");
#m_gp   = GP(y, x, x, K, c);
m_rgp = RGP(y, x, K, M, c);

using Plots
plot(X, Y, label = "f(x)", color = "black", legend = :topleft)
scatter!(x, y, label = "Observations", color = "red")
plot!(x, m_krls.predictions, label = "KRLS", color = "green")
scatter!(x[m_krls.basis], y[m_krls.basis], label = "Dictionary", color = "blue")

ŷ = predict(m_krls.μ, m_krls.Q, reshape(X, :, 1), m_krls.xb, k_krls)
plot!(X, ŷ, label = "Prediction", color = "pink")