text:
using Plots
gr
x = 1:10
y = rand(10)
plot(x,y)

y = rand(10,2)
plot(x,y)

z = rand(10)
plot!(x,z)

p = plot(x,y)
z = rand(10)
plot!(p, x, z)

# NG
# plot!(p, linewidth = 10, label = ["L1", "L2"], title="Two Lines")

# -----------
# Attribute
# -----------
p = plot(x, y, linewidth = 5, label = ["L1" "L2"], title="Two Lines")
xlabel!(p, "X")
#attribute!(p, linewidth = 1)

# plotattr()
# plotattr(:Axis)
# plotattr(:Subplot)
# plotattr(:Series)
# plotattr(:Plot)
# plotattr("linewidth")
# plotattr("seriestype")
plot(x,y, seriestype=:scatter)
scatter(x,y)

# --------
# Layout
# --------
y = rand(10,4)
plot(x,y, layout=(4,1))

p1 = plot(x,y)
p2 = scatter(x,y)
p3 = plot(x,y, xlabel = "Labelled", lw = 3, title = "Subtitle")
p4 = histogram(x, y)
plot(p1, p2, p3, p4, layout = (2,2), legend = false)

# -----------------------------------
# Plot Recipes and Recipe Libraries
# http://docs.juliaplots.org/latest/tutorial/#Plot-Recipes-and-Recipe-Libraries-1
# -----------------------------------
using StatsPlots
using DataFrames
df = DataFrame(a = 1:10, b = 10*rand(10), c = 10*rand(10))
@df df plot(:a, [:b :c])
@df df scatter(:a, :b, title = "DataFrame Scatter Plot!")

using Distributions
plot(Normal(3,5), lw = 3)

using RDatasets
iris = dataset("datasets", "iris")
@df iris marginalhist(:PetalLength, :PetalWidth)


text:
using CairoMakie, Random
Random.seed!(113)
n = 100
x = randn(n)
y = cumsum(x) .- minimum(cumsum(x))
leftPoints =  Point2f.(tuple.(1, LinRange(minimum(y), maximum(y),n)))
rightPoints = Point2f.(tuple.(n, LinRange(minimum(y), maximum(y),n)))
fig = Figure(resolution = (600, 400))
ax = Axis(fig[1,1])
band!(ax, leftPoints, rightPoints, color = 1:length(leftPoints),
    colormap = :devon)
band!(ax, 1:n, y, fill(maximum(y), n); color = :white)
lines!(ax, y; color = :grey20)
foreach(ele -> translate!(ax.elements[ele], 0, 0, 9),
    [:ygridlines, :xgridlines])
save("bandHorizontal.png", fig, px_per_unit = 2)
fig

ext: jl


text:
SUITE["dot"] = BenchmarkGroup()
SUITE["dcontract"] = BenchmarkGroup()
SUITE["otimes"] = BenchmarkGroup()
SUITE["other"] = BenchmarkGroup()
SUITE["promotion"] = BenchmarkGroup()

for dim in (1,2,3)
    for T in (Float64, Float32)
        v1 = tensor_dict[(dim, 1, T)]
        V2 = tensor_dict[(dim, 2, T)]
        V4 = tensor_dict[(dim, 4, T)]
        V2sym = symtensor_dict[(dim, 2, T)]
        V4sym = symtensor_dict[(dim, 4, T)]

        # dot
        SUITE["dot"]["dim $dim - order 1 1 - $T"]             = @benchmarkable dot($v1, $v1)
        SUITE["dot"]["dim $dim - order 2 1 - $T"]             = @benchmarkable dot($V2, $v1)
        SUITE["dot"]["dim $dim - order 2sym 1 - $T"]          = @benchmarkable dot($V2sym, $v1)
        SUITE["dot"]["dim $dim - order 2 2 - $T"]             = @benchmarkable dot($V2, $V2)
        SUITE["dot"]["dim $dim - order 2sym 2sym - $T"]       = @benchmarkable dot($V2sym, $V2sym)
        SUITE["dot"]["dim $dim - order 2sym 2 - $T"]          = @benchmarkable dot($V2sym, $V2)

        # dcontract
        SUITE["dcontract"]["dim $dim - order 2 2 - $T"]       = @benchmarkable dcontract($V2, $V2)
        SUITE["dcontract"]["dim $dim - order 2sym 2sym - $T"] = @benchmarkable dcontract($V2sym, $V2sym)
        SUITE["dcontract"]["dim $dim - order 2sym 2 - $T"]    = @benchmarkable dcontract($V2sym, $V2)

        SUITE["dcontract"]["dim $dim - order 4 2 - $T"]       = @benchmarkable dcontract($V4, $V2)
        SUITE["dcontract"]["dim $dim - order 4sym 2 - $T"]    = @benchmarkable dcontract($V4sym, $V2)
        SUITE["dcontract"]["dim $dim - order 2 4sym - $T"]    = @benchmarkable dcontract($V2, $V4sym)
        SUITE["dcontract"]["dim $dim - order 2 4 - $T"]       = @benchmarkable dcontract($V2, $V4)
        SUITE["dcontract"]["dim $dim - order 4 2sym - $T"]    = @benchmarkable dcontract($V4, $V2sym)
        SUITE["dcontract"]["dim $dim - order 4sym 2sym - $T"] = @benchmarkable dcontract($V4, $V2sym)

        SUITE["dcontract"]["dim $dim - order 4 4 - $T"]       = @benchmarkable dcontract($V4, $V4)
        SUITE["dcontract"]["dim $dim - order 4sym 4 - $T"]    = @benchmarkable dcontract($V4sym, $V4)
        SUITE["dcontract"]["dim $dim - order 4sym 4sym - $T"] = @benchmarkable dcontract($V4sym, $V4sym)

        # otimes
        SUITE["otimes"]["dim $dim - order 1 1 - $T"]          = @benchmarkable otimes($v1, $v1)
        SUITE["otimes"]["dim $dim - order 2 2 - $T"]          = @benchmarkable otimes($V2, $V2)
        SUITE["otimes"]["dim $dim - order 2sym 2sym - $T"]    = @benchmarkable otimes($V2sym, $V2sym)
        SUITE["otimes"]["dim $dim - order 2sym 2 - $T"]       = @benchmarkable otimes($V2sym, $V2)

        # other
        for (i, V2t) in enumerate((V2, V2sym))
            symstr = i == 2 ? "sym" : ""
            for f in (norm, trace, vol, det, inv, transpose, symmetric, skew, eig, mean, dev)
                i == 1 && f == eig && continue
                SUITE["other"]["$f - dim $dim - order 2$(symstr) - $T"] = @benchmarkable $f($V2t)
            end
        end

        for (i, V4t) in enumerate((V4, V4sym))
            symstr = i == 2 ? "sym" : ""
            for f in (norm, symmetric)
                SUITE["other"]["$f - dim $dim - order 4$(symstr) - $T"] = @benchmarkable $f($V4t)
            end
        end

        SUITE["promotion"]["dim $dim - order 2"] = @benchmarkable promote($V2, $V2sym)
        SUITE["promotion"]["dim $dim - order 4"] = @benchmarkable promote($V4, $V4sym)
    end
end


text:
module Camera

export SimpleCamera, FovCamera, getray

using CthulhuVision.Math
using CthulhuVision.Random
using CthulhuVision.Light

struct SimpleCamera
    origin     :: Vector3
    lowerleft  :: Vector3
    horizontal :: Vector3
    vertical   :: Vector3
end

@inline function getray(cam::SimpleCamera, u::Float32, v::Float32, rng::UniformRNG) :: Ray
    direction = cam.lowerleft + u*cam.horizontal + v*cam.vertical - cam.origin
    Ray(cam.origin, direction)
end

@inline unithost(v::Vector3) = 1.0f0 / sqrt(v.x * v.x + v.y * v.y + v.z * v.z) * v

@inline function randominunitdisk(rng::UniformRNG) :: Vector3
    p = Vector3(1.0f0, 1.0f0, 0.0f0)

    while dot(p, p) >= 1.0f0
        p = 2.0f0 * Vector3(next(rng), next(rng), 0.0f0) - Vector3(1.0f0, 1.0f0, 0.0f0)
    end

    p
end

struct FovCamera
    origin     :: Vector3
    lowerleft  :: Vector3
    horizontal :: Vector3
    vertical   :: Vector3
    u          :: Vector3
    v          :: Vector3
    w          :: Vector3
    lensradius :: Float32

    function FovCamera(lookfrom::Vector3, lookat::Vector3, vup::Vector3, vfov::Float32, aspect::Float32, aperture::Float32, focusdist::Float32)
        lensradius = aperture / 2.0f0
        theta = vfov * pi / 180.0f0
        halfheight = tan(theta / 2.0f0)
        halfwidth = aspect * halfheight

        origin = lookfrom
        w = unithost(lookfrom - lookat)
        u = unithost(cross(vup, w))
        v = cross(w, u)

        lowerleft  = origin - halfwidth * focusdist * u - halfheight * focusdist * v - focusdist * w
        horizontal = 2.0f0 * halfwidth * focusdist * u
        vertical   = 2.0f0 * halfheight * focusdist * v

        new(origin, lowerleft, horizontal, vertical, u, v, w, lensradius)
    end
end

@inline function getray(cam::FovCamera, s::Float32, t::Float32, rng::UniformRNG) :: Ray
    rd = cam.lensradius * randominunitdisk(rng)
    offset = rd.x * cam.u + rd.y * cam.v
    direction = cam.lowerleft + s*cam.horizontal + t*cam.vertical - cam.origin - offset
    Ray(cam.origin + offset, direction)
end

end


text:
using Plots
using Statistics
using LaTeXStrings
function std_calculater(arr,stdList)
    hight=[]
    for col in eachcol(arr)
        count=0
        for element in col
            if element!=0
                count+=1
            end
        end
        push!(hight, count)
    end
    std_num=std(hight)
    push!(stdList, std_num)
    return stdList
end
L=200
arr=zeros((L,L))
color=1
stdList=[]
t_interval=ceil.(Int, exp.(1:0.5:10))
for i in 1:18
    falling_n=t_interval[i+1]-t_interval[i]
    for each in 1:falling_n
        column=rand(1:200)
        row=1
        for j in 1:200
            if arr[row, column]== 0.0
                arr[row, column]= color
                break
            else
                row+=1
            end
        end
        if i%(10*200*color)==0
            color+=1
        end
    end
    stdList=std_calculater(arr,stdList)
end
BigStdList=[]
run_number=100
for num in 1:run_number
    L=200
    arr=zeros((L,L))
    color=1
    stdList=[]
    time=ceil.(Int, exp.(0:0.5:10))
    for i in 1:20
        falling_n=time[i+1]-time[i]
        for each in 1:falling_n
            column=rand(1:200)
            row=1
            for j in 1:200
                if arr[row, column]== 0.0
                    arr[row, column]= color
                    break
                else
                    row+=1
                end
            end
            if i%(10*200*color)==0
                color+=1
            end
        end
        stdList=std_calculater(arr,stdList)
    end
    push!(BigStdList, stdList)
end
BiggerMeanList=[]
BiggerStdList=[]
for i in 1:19
    tempvalue=[]
    for j in 1:100
        push!(tempvalue,BigStdList[j][i])
    end
    push!(BiggerStdList, std(tempvalue))
    push!(BiggerMeanList, mean(tempvalue))
end
function LineFit(Time, Mean)
    a=[hcat(log.(Time)) reshape(ones(19), 19, 1)]
    b=hcat(log.(Mean))
    line=(a\b)
    x = 1:10
    y = x .* line[1] .+ line[2]
    return x, y, line
end
X,Y,Line=LineFit(t_interval, BiggerMeanList)
plot(log.(t_interval), log.(BiggerMeanList),yerr=BiggerStdList, xlabel="time", ylabel="w")
plot!(X,Y,c= :black,label = L"y = %$(round(Line[1],digits= 2))x + %$(round(Line[2],digits= 2))")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter3\\Fig\\RandomBallisticDeposition3.png")


text:
module NoBang

using Base.Iterators: Pairs
using Base: ImmutableDict
using Requires
using Setfield: setproperties

if VERSION < v"1.1-"
    using Future: copy!
end

include("singletoncontainers.jl")
include("base.jl")
include("linearalgebra.jl")

function __init__()
    @require StaticArrays="90137ffa-7385-5640-81b9-e52037218182" begin
        include("staticarrays.jl")
    end
end

end  # module


text:
function clc()
    if Sys.islinux()
        Base.run(`clear`)
    end
    return
end
clc();

max_stars_repo_name: manuelmontesinos/Econometrics

max_issues_repo_name: manuelmontesinos/Econometrics

max_forks_repo_name: manuelmontesinos/Econometrics


text:
using BinaryBuilder

cuda_version = v"10.2.89"   # NOTE: could be less specific

sources_linux_x64 = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.2.39/10.2_20200626/cudnn-10.2-linux-x64-v8.0.2.39.tgz",
                  "c9cbe5c211360f3cfbc0fb104f0e9096b37e53f89392525679f049276b2f701f")
]
sources_linux_ppc64le = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.2.39/10.2_20200626/cudnn-10.2-linux-ppc64le-v8.0.2.39.tgz",
                  "c32325ff84a8123491f2e58b3694885a9a672005bc21764b38874688c0e43262")
]
sources_windows = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.2.39/10.2_20200626/cudnn-10.2-windows10-x64-v8.0.2.39.zip",
                  "ba7d544fba2aaa7e6445cc41d6890a890567ea14ca942c13115f7d5654ea451f")
]

include("../common.jl")

max_stars_repo_name: jpsamaroo/Yggdrasil

max_issues_repo_name: jpsamaroo/Yggdrasil

max_forks_repo_name: jpsamaroo/Yggdrasil


text:
"""
< DEFINITION of DATA TYPES for JuliaPBF >
"""
module StructPBF
import Base.+, Base.*, Base.-, Base.isless, Base.abs2 # To add generic functions
export
    AnalysisDataStruct,
    SimulationDataStruct

    include("BasicStruct.jl")
    include("AnalysisDataStruct.jl")
    include("SimulationDataStruct.jl")

    """시뮬레이션 데이터"""
    struct AnalysisDataStruct
        timeStep            ::TimeStep
        endTime             ::Float64
        gravity             ::Vec2
        analysisBox         ::Box2
        nMaxParticles       ::Int64
        maxParticleSize     ::Float64
        kernelRadius        ::Float64
        tensileInstability  ::TensileInstability
        a_initialBox        ::Vector{InitialBox}
        a_phases            ::Vector{PhaseData}
    end

    struct SimulationDataStruct
        particles           ::Vector{ParticleData} # For AoS
        grids               ::GridData
        phases              ::Vector{PhaseData}
        boundaries          ::Vector{BoundaryData}
        alivedParticles     ::Vector{Int64} # Fluids / Solids
    end

end

max_stars_repo_name: Hipgineer/JuliaPBF

max_issues_repo_name: Hipgineer/JuliaPBF

max_forks_repo_name: Hipgineer/JuliaPBF


text:
using BinaryBuilder

name = "METIS"
version = v"5.1.1" # <-- This is a lie, we're bumping to 5.1.1 to create a Julia v1.6+ release with experimental platforms

# Collection of sources required to build METIS
sources = [
    ArchiveSource("http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/metis-5.1.0.tar.gz",
                  "76faebe03f6c963127dbb73c13eab58c9a3faeae48779f049066a21c087c5db2"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
# Patches from https://github.com/msys2/MINGW-packages/tree/master/mingw-w64-metis
script = raw"""
cd $WORKSPACE/srcdir/metis-*
if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
    atomic_patch -p1 $WORKSPACE/srcdir/patches/0001-mingw-w64-does-not-have-sys-resource-h.patch
    atomic_patch -p1 $WORKSPACE/srcdir/patches/0002-mingw-w64-do-not-use-reserved-double-underscored-names.patch
    atomic_patch -p1 $WORKSPACE/srcdir/patches/0003-WIN32-Install-RUNTIME-to-bin.patch
    atomic_patch -p1 $WORKSPACE/srcdir/patches/0004-Fix-GKLIB_PATH-default-for-out-of-tree-builds.patch
fi
mkdir -p build
cd build/
cmake $WORKSPACE/srcdir/metis-5.1.0/ \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
    -DCMAKE_VERBOSE_MAKEFILE=1 \
    -DGKLIB_PATH=$WORKSPACE/srcdir/metis-5.1.0/GKlib \
    -DSHARED=1
make -j${nproc} install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libmetis", :libmetis)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[]

# Build the tarballs
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:
# sort clause
#include("misc.jl")

# for equalclause, they should be sorted.
#==
function renamevarinlit(vars, lit)
# dont complete
 for v in vars
  lit = replace(lit, string(v), "*")
 end
 return lit
end

function gelit(var1,L1,var2,L2)
 sl1 = renamevarinlit(var1,string(L1))
 sl2 = renamevarinlit(var2,string(L2))
 return(sl1<=sl2)
end

# sort literals for uniquness
function sortcls(vars, cls)
 nargs=sort(cls, lt=(x,y)->gelit(vars,x,vars,y))
 return vars, nargs
end

==#
#### readclause

#from a file
function readclausefromfile(fname)
 lines=readlines(fname)
 clss = []
 proc = []
 for line in lines
  if length(line)>0
    if line[1]=='['
      push!(clss,Meta.parse(line))
    elseif line[1] == '!'
      push!(proc,line[2:end])
    elseif line[1] == '<'
      iclss, iproc = readclausefromfile(strip(line[2:end],[' ','\t']))
      append!(clss, iclss)
      append!(proc, iproc)
    end
  end
 end
 return clss, proc
end

function readkpclausefromfile(fname)
 lines=readlines(fname)
 clss = []
 proc = []
 for line in lines
  if length(line)>0
    if line[1]=='['
      push!(clss,kpparse(line))
    elseif line[1] == '!'
      push!(proc,line[2:end])
    elseif line[1] == '<'
      iclss, iproc = readkpclausefromfile(strip(line[2:end],[' ','\t']))
      append!(clss, iclss)
      append!(proc, iproc)
    end
  end
 end
 return clss, proc
end

function printcnf(fname)
 cnfs, proc=readclausefromfile(fname)
 @show proc
 for p in proc
   println("!$p")
 end
println()
 for ix in 1:length(cnfs)
   cnf = cnfs[ix]
   println("$(cidof(ix)): $cnf")
 end
end

#### readcore
function readcore(fname)
 println("readcore fname=$fname")
 cls,proc = readclausefromfile(fname)
 createcore(fname, cls, proc)
end

function readcoredir(dirname)
 cores = []
 for (r,ds,fs) in walkdir(dirname)
   for f in fs
     push!(cores,readcore("$r/$f"))
   end
 end
 cores
end

####
iscap(x) = isuppercase(x[1])
isinvar(x) = '!' == String(x)[end]

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
## Generate samples from the prior distribution
function genSamplesPrior(model_def, bayinf_def, nsamps, mu=[],coo=[])

    sam = ["lognormal", "normal", "uniform", "bernoulli", "binomial", "neg_binomial", "poisson", "student_t", "cauchy",
    "double_exponential", "gumbel", "chi_square", "exponential", "gamma", "weibull", "frechet",  "rayleigh", "pareto", "beta", "von_mises", "dirichlet"];

    println("-------------------------------------- WARNING --------------------------------------")
    println("This function will only work if the contents from bayinf_def[","\"","Priors","\"","][","\"","pridis","\"","] are")
    println("distributions with numbers defined in the parameters. If the parameters given to the ")
    println("distributions are variables instead of numbers an error will be prompted since sampling will")
    println("not work. So no hierarchical priors. Sorry for the inconvenience. ")
    println(" ")
    println("The only exception is the multivariate normal case where you can introduce the mean vector")
    println("and the covariance matrix as 4th and 5th arguments in the function in this order.")
    println("However, truncations will not be possible and only 1 MultiVariate normal can be assessed. ")
    println("")
    println("Also, sampling from all types of stan distributions is not implemented yet. For now the ")
    println("available options are: ")
    println(join(sam, ", "))
    println(" ")
    println("-------------------------------------------------------------------------------------")
    println(" ")

    origsamp = Array{String,1}(undef, length(bayinf_def["Priors"]["pridis"]));

    for i in 1:length(bayinf_def["Priors"]["pridis"])
        origsamp[i] = replace(bayinf_def["Priors"]["pridis"][i], bayinf_def["Priors"]["pridis"][i][1:findfirst("~", bayinf_def["Priors"]["pridis"][i])[1]] => "")
        origsamp[i] = replace(origsamp[i], ";" => "")
        origsamp[i] = replace(origsamp[i], "\n" => "")
        if occursin("multi_normal", bayinf_def["Priors"]["pridis"][i])
            if mu == [] || coo ==[]
                println("Process Stopped, please give means and covariances for the multinormal");
                return
            else
                origsamp[i] = replace(bayinf_def["Priors"]["pridis"][i], bayinf_def["Priors"]["pridis"][i][findfirst("~", bayinf_def["Priors"]["pridis"][i])[1]:end] => string("= MvNormal(",mu,", ",coo,")"))
            end
        elseif occursin("lognormal", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "lognormal" => "LogNormal");
        elseif occursin("normal", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "normal" => "Normal");
        elseif occursin("uniform", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "uniform" => "Uniform");
        elseif occursin("bernoulli", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "bernoulli" => "Bernoulli");
        elseif occursin("binomial", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "binomial" => "Binomial");
        elseif occursin("neg_binomial", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "neg_binomial" => "NegativeBinomial");
        elseif occursin("poisson", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "poisson" => "Poisson");
        elseif occursin("student_t", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "student_t" => "TDist");
        elseif occursin("cauchy", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "cauchy" => "Cauchy");
        elseif occursin("double_exponential", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "double_exponential" => "Laplace");
    #     elseif occursin("logistic", bayinf_def["Priors"]["pridis"][i])
    #         origsamp[i] = replace(origsamp[i], "" => "");
        elseif occursin("gumbel", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "gumbel" => "Gumbel");
        elseif occursin("chi_square", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "chi_square" => "Chisq");
        elseif occursin("exponential", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "exponential" => "Exponential");
        elseif occursin("gamma", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "gamma" => "Gamma");
        elseif occursin("weibull", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "weibull" => "Weibull");
        elseif occursin("frechet", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "frechet" => "Frechet");
        elseif occursin("rayleigh", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "rayleigh" => "Rayleigh");
        elseif occursin("pareto", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "pareto" => "Pareto");
        elseif occursin("beta", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "beta" => "Beta");
        elseif occursin("von_mises", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "von_mises" => "VonMises");
        elseif occursin("dirichlet", bayinf_def["Priors"]["pridis"][i])
            origsamp[i] = replace(origsamp[i], "dirichlet" => "Dirichlet");
        end
    end

    for i in 1:length(bayinf_def["Priors"]["pars"])
        if !occursin("MvNormal", origsamp[i])
            if occursin("lower=", bayinf_def["Priors"]["pars"][i])
                ind1 = findfirst("lower=", bayinf_def["Priors"]["pars"][i])[end]+1;
                ind2 = "";
                try
                    ind2 = findfirst(",", bayinf_def["Priors"]["pars"][i])[1]-1;
                catch
                    ind2 = findfirst(">", bayinf_def["Priors"]["pars"][i])[1]-1;
                end
                lob = bayinf_def["Priors"]["pars"][i][ind1:ind2]
            else
                lob = [];
            end

            if occursin("upper=", bayinf_def["Priors"]["pars"][i])
                ind3 = findfirst("upper=", bayinf_def["Priors"]["pars"][i])[end]+1;
                ind4 = findfirst(">", bayinf_def["Priors"]["pars"][i])[1]-1;
                upb = bayinf_def["Priors"]["pars"][i][ind3:ind4]
            else
                upb = [];
            end

            if lob == [] && upb !=[]
                lob = "convert(Float64,*Inf)"
            elseif lob != [] && upb ==[]
                upb = "convert(Float64,Inf)"
            end

            if lob != [] && upb != []
                origsamp[i] = string("truncated(",origsamp[i],", ",lob,", ",upb,")")
            else
                origsamp[i] = "";
            end
        end
    end

    origsampTR = Array{String,1}(undef, length(bayinf_def["Priors"]["pars"]));

    for i in 1:length(bayinf_def["Priors"]["pars"])
        for j in 1:length(bayinf_def["Priors"]["transpars"])
            if occursin(model_def["parName"][i], bayinf_def["Priors"]["transpars"][j])
                origsampTR[i] = replace(bayinf_def["Priors"]["transpars"][j], ";"=>"")
                origsampTR[i] = replace(origsampTR[i], "\n"=>"")
                origsampTR[i] = replace(origsampTR[i], origsampTR[i][1:findfirst("=", origsampTR[i])[1]]=>"")
            end
        end
        if !isassigned(origsampTR,i)
            origsampTR[i] = ""
        end
    end

    pretransP = zeros(model_def["nPar"],nsamps);
    for i in 1:length(origsamp)
        tpp = Meta.parse(string("rand(", origsamp[i], ",", nsamps, ")"));
        valp = @eval $tpp;
        pretransP[i,:] = valp;
    end

    priorSamps = zeros(nsamps, model_def["nPar"]);
    for i in 1:model_def["nPar"]
        if origsampTR != ""
            tmpex = origsampTR[i];
            if occursin("*", origsampTR[i])
                tmpex = replace(tmpex, "*"=>".*")
            end
            if occursin("/", origsampTR[i])
                tmpex = replace(tmpex, "/"=>"./")
            end
            if occursin("+", origsampTR[i])
                tmpex = replace(tmpex, "+"=>".+")
            end
            if occursin("-", origsampTR[i])
                tmpex = replace(tmpex, "-"=>".-")
            end
            if occursin("^", origsampTR[i])
                tmpex = replace(tmpex, "^"=>".^")
            end

            for j in 1:model_def["nPar"]
                if occursin(model_def["parName"][j], tmpex)
                    tmpex = replace(tmpex, model_def["parName"][j] => string(pretransP[j,:]));
                end
            end
            tppT = Meta.parse(tmpex);
            valpT = @eval $tppT;
            priorSamps[:,i] = valpT;
        else
            priorSamps[:,i] = pretransP[i,:];
        end
    end

    return priorSamps
end

##
# --------------------------------- RELATIVE ENTROPY FUNCTION --------------------------------- #

# Calculation of the entropy of prior and posterior results for an inference result using the Taylor series approximation
# for Gaussian Mixtures proposed in "On Entropy Approximation for Gaussian Mixture Random Vectors". For all entropy calculations,
# x represents a random vector, MU the vectors of means for the Gaussian mixtures, E the covariance matrices for the Gaussian
# mixtures and w the vector of weights for the Gaussian Mixtures. As inputs the function takes a character string with the
# name of the experimental profile to be analysed.

## Calculation of upper bound for the posterior Entropy

function H_Upper(w,E)

    comp = length(w);
    dims = size(E[1])[1];
    Hu = zeros(1,comp);

    for i in 1:comp
        hu = w[i]*(-log(w[i])+0.5*log(det(Matrix(E[i]))*(2*pi*exp(1))^dims))
        Hu[1,i] = hu;
    end

    return(sum(Hu))
end

## PDF of a multivariate Gaussian distribution

function mvGauss(x, MU, E)

    gx = (1/sqrt(det(Matrix(E))*(2*pi)^2))*exp((-1/2)*(x-MU)'*inv(Matrix(E))*(x-MU))

    return(gx)

end

## Lower bound for the posterior entropy
function H_Lower(w, E, MU)

    comp = length(w);
    Hl = zeros(1,comp);

    for i in 1:comp
        inter = zeros(1,comp);
        for j in 1:comp
            insu = w[j]*BOMBS.mvGauss(MU[i,:], MU[j,:], E[i]+E[j]);
            inter[1,j] = insu;
        end
        Hl[1,i] = w[i]*log(sum(inter))
    end

    return(-sum(Hl))

end

## Function for a multivariate PDF of a Gaussian Mixture
function GaussMix(x, MU, E, w)

    comp = length(w);
    FX = zeros(1,comp);

    for i in 1:comp
        fx = w[i]*(1/sqrt(det(Matrix(E[i]))*(2*pi)^2))*exp((-1/2)*(x-MU[i,:])'*inv(Matrix(E[i]))*(x-MU[i,:]));
        FX[i] = fx;
    end

    return(sum(FX))

end

## Zero-order Taylor Series Expansion

function ZOTSE(MU, E, w)

    comp = length(w);
    ZO = zeros(1,comp);

    for i in 1:comp
        zo = w[i]*log(GaussMix(MU[i,:], MU, E, w));
        ZO[i] = zo;
    end

    return(-sum(ZO))
end

## Function to calculate the most computationaly expensive part of the Taylor Expansion component
function GaussMix2(x)

    comp = length(mvPro);
    FX = zeros(1,comp);

    for i in 1:comp
        fx = mvPro[i]*(1/sqrt(det(Matrix(mvCovar[i]))*(2*pi)^2))*exp((-1/2)*(x-mvMean[i,:])'*inv(Matrix(mvCovar[i]))*(x-mvMean[i,:]));
        FX[i] = fx;
    end

    return(sum(FX))

end

function FMix(x, MU, E, w)

    comp = length(w);
    dims = size(E[1])[1];

    F_x = Array{Any,1}(undef,comp);

    for j in 1:comp
        inter = (1/GaussMix(x, MU, E, w))*(x-MU[j,:])*(Calculus.gradient(GaussMix2, x))'+(x-MU[j,:])*(inv(Matrix(E[j]))*(x-MU[j,:]))'-Diagonal(ones(dims));
        fx = w[j]*inv(Matrix(E[j]))*(inter)*(mvGauss(x , MU[j,:], E[j]));
        F_x[j] = fx;
    end
    enn = sum(F_x)/GaussMix(x, MU, E, w)

    return(enn)
end

## Calculation of the Second-Order Taylor Series expanssion term

function SOTSE(MU, E, w)

    comp = length(w);
    dims = size(E[1])[1];

    H2 = zeros(1,comp);

    for i in 1:comp
        println(string("              Iteration ", i, " of ", comp))
        inter2 = (w[i]/2)*sum(FMix(MU[i,:],MU,E,w).*Matrix(E[i]))
        H2[i] = inter2;
    end

    return(sum(H2))

end

## Function to estimate the ENTROPY
# Function to compute entropy from distribution samples
function computeH(sampl, model_def, tag = "")

    if length(size(sampl)) == 1
        sampl = reshape(sampl,length(sampl), 1);
    end

    if size(sampl)[1]==model_def["nPar"]
        sampl = convert(Array, sampl');
    elseif size(sampl)[2]==model_def["nPar"]
        nothing
    else
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the samples that you have introduced. There is something odd in there. ")
        return
    end

    cudi = pwd(); #@__DIR__;
    if !isdir(string(cudi,"\\Results"))
        mkdir(string(cudi,"\\Results"))
    end
    if !isdir(string(cudi, "\\Results\\", model_def["NameF"],"_",today()))
        mkdir(string(cudi, "\\Results\\", model_def["NameF"],"_",today()))
    end

    # Estimate joint posterior as Gaussian Mixtures
    # Perform clustering in a range between 1 to 40 components and use LLK to asses the optimum value
    grid_search = ScikitLearnBase.fit!(GridSearchCV(GMM(n_components=1, kind=:full), Dict(:n_components=>collect(1:40))), sampl)

    # Get and store best results
    comps = length(grid_search.best_estimator_.w); # Number of components
    weis = grid_search.best_estimator_.w; # Weights for each component
    meaas = grid_search.best_estimator_.μ; # Means for each component
    covas = [inv(Matrix(grid_search.best_estimator_.Σ[i]))*inv(Matrix(grid_search.best_estimator_.Σ[i])') for i in 1:length(grid_search.best_estimator_.Σ)]; # Covariance Matrix for each component

    entro_res = Dict();
    entro_res["BestGC"] = Dict();
    entro_res["BestGC"]["c"] = comps;
    entro_res["BestGC"]["w"] = weis;
    entro_res["BestGC"]["mu"] = meaas;
    entro_res["BestGC"]["cov"] = covas;

    global mvMean = meaas;
    global mvCovar = covas;
    global mvPro = weis;

    # Upper bound for the posterior entropy
    H_U_Post = H_Upper(weis, covas);
    entro_res["HUpper"] = H_U_Post;

    # Lower bound for the posterior entropy
    H_L_Post = H_Lower(weis, covas, meaas);
    entro_res["HLower"] = H_L_Post;

    # Zero-order Taylor series result for the approximation of the posterior entropy
    ZO_Expan = ZOTSE(meaas, covas, weis);

    # Second-order Taylor series result term for the approximation of the posterior entropy
    SO_Expan = SOTSE(meaas, covas, weis);

    # Aproximation of the Entropy for the posterior
    H_posterior = ZO_Expan-SO_Expan;
    entro_res["HTheta"] = H_posterior;

    # Assess if entropy estimation is between bounds
    if H_posterior>=H_L_Post && H_posterior<=H_U_Post
        eebb = true;
    else
        eebb = false;
    end
    entro_res["BetweenBounds"] = eebb;

    entro_res["Cov"] = cov(sampl);

    JLD.save(string(cudi, "\\Results\\", model_def["NameF"],"_",today(), "\\EntropyResults_", tag, ".jld"), "entro_res", entro_res)

    println("")
    println("----------------------------------------- ENTROPY RESULTS -----------------------------------------")
    println("Entropy results are saved in the directory: ")
    println(string("    ",cudi, "\\Results\\", model_def["NameF"],"_",today(), ""))
    println("Under the name: ")
    println(string("            EntropyResults_", tag, ".jld"))
    println("--------------------------------------------------------------------------------------")
    println("")

    return(entro_res)

end

## Function to compute entropy gain between prior and posterior (based on samples)
# In this case bayinf_def can only have the entry Priors (nothing else will be checked)
function computeHgain(prior, posterior, model_def, tag = "")

    cudi = pwd(); #@__DIR__;
    if !isdir(string(cudi,"\\Results"))
        mkdir(string(cudi,"\\Results"))
    end
    if !isdir(string(cudi, "\\Results\\", model_def["NameF"],"_",today()))
        mkdir(string(cudi, "\\Results\\", model_def["NameF"],"_",today()))
    end

    if typeof(prior) <: Dict
        pris = genSamplesPrior(model_def, prior, 10000);
    else
        pris = prior;
    end;

    entro_prior = computeH(pris, model_def, "prior");
    entro_posterior = computeH(posterior, model_def, "posterior");

    Relative_H = entro_prior["HTheta"]-entro_posterior["HTheta"];

    Hdiff_res = Dict();
    Hdiff_res["HPrior"] = entro_prior;
    Hdiff_res["HPosterior"] = entro_posterior;
    Hdiff_res["HDiff"] = Relative_H;

    JLD.save(string(cudi, "\\Results\\", model_def["NameF"],"_",today(), "\\EntropyGainResults_", tag, ".jld"), "Hdiff_res", Hdiff_res)

    println("")
    println("----------------------------------------- ENTROPY RESULTS -----------------------------------------")
    println("Entropy gain results are saved in the directory: ")
    println(string("    ",cudi, "\\Results\\", model_def["NameF"],"_",today(), ""))
    println("Under the name: ")
    println(string("            EntropyGainResults_", tag, ".jld"))
    println("--------------------------------------------------------------------------------------")
    println("")

    return(Hdiff_res)

end

max_stars_repo_name: DavidGomezC/BOMBSbeta

max_issues_repo_name: DavidGomezC/BOMBSbeta

max_forks_repo_name: DavidGomezC/BOMBSbeta


text:
module TestDoctest

using UnionArrays
using Documenter: doctest
using Test

@testset begin
    doctest(UnionArrays)
end

end  # module


text:
@testset "Scheduler Tests" begin
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
# This file contains a list of the currently available covariance functions

abstract type Kernel end

"""
    KernelData

Data to be used with a kernel object to calculate a covariance matrix,
which is independent of kernel hyperparameters.

See also [`EmptyData`](@ref).
"""
abstract type KernelData end

"""
    EmptyData <: KernelData

Default empty `KernelData`.
"""
struct EmptyData <: KernelData end

KernelData(k::Kernel, X::MatF64) = EmptyData()
kernel_data_key(k::Kernel, X::MatF64) = "EmptyData"

"""
    cov(k::Kernel, X₁::Matrix{Float64}, X₂::Matrix{Float64})

Create covariance matrix from kernel `k` and matrices of observations `X₁` and `X₂`, where
each column is an observation.
"""
function Statistics.cov(k::Kernel, X₁::MatF64, X₂::MatF64)
    d(x1,x2) = cov(k, x1, x2)
    return map_column_pairs(d, X₁, X₂)
end

"""
    cov!(cK::Matrix{Float64}, k::Kernel, X₁::Matrix{Float64}, X₂::Matrix{Float64})

Like [`cov(k, X₁, X₂)`](@ref), but stores the result in `cK` rather than a new matrix.
"""
function cov!(cK::MatF64, k::Kernel, X₁::MatF64, X₂::MatF64)
    d(x1,x2) = cov(k, x1, x2)
    return map_column_pairs!(cK, d, X₁, X₂)
end

"""
    cov(k::Kernel, X::Matrix{Float64}[, data::KernelData = KernelData(k, X)])

Create covariance function from kernel `k`, matrix of observations `X`, where each column is
an observation, and kernel data `data` constructed from input observations.
"""
Statistics.cov(k::Kernel, X::MatF64, data::EmptyData) = cov(k, X)
cov!(k::Kernel, X::MatF64, data::EmptyData) = cov!(cK, k, X)

function cov!(cK::MatF64, k::Kernel, X::MatF64)
    dim, nobsv = size(X)
    @inbounds for j in 1:nobsv
        cK[j,j] = cov_ij(k, X, j, j, dim)
        for i in 1:j-1
            cK[i,j] = cov_ij(k, X, i, j, dim)
            cK[j,i] = cK[i,j]
        end
    end
    return cK
end
function Statistics.cov(k::Kernel, X::MatF64)
    dim, nobsv = size(X)
    cK = Array{Float64}(undef, nobsv, nobsv)
    cov!(cK, k, X)
end
function cov!(cK::MatF64, k::Kernel, X::MatF64, data::KernelData)
    dim, nobsv = size(X)
    @inbounds for j in 1:nobsv
        cK[j,j] = cov_ij(k, X, data, j, j, dim)
        for i in 1:j-1
            cK[i,j] = cov_ij(k, X, data, i, j, dim)
            cK[j,i] = cK[i,j]
        end
    end
    return cK
end
function Statistics.cov(k::Kernel, X::MatF64, data::KernelData)
    dim, nobsv = size(X)
    cK = Array{Float64}(undef, nobsv, nobsv)
    cov!(cK, k, X, data)
end

@inline cov_ij(k::Kernel, X::MatF64, i::Int, j::Int, dim::Int) = cov(k, @view(X[:,i]), @view(X[:,j]))
@inline cov_ij(k::Kernel, X::MatF64, data::EmptyData, i::Int, j::Int, dim::Int) = cov_ij(k, X, i, j, dim)

############################
##### Kernel Gradients #####
############################
@inline @inbounds function dKij_dθ!(dK::VecF64, kern::Kernel, X::MatF64,
                                    i::Int, j::Int, dim::Int, npars::Int)
    for p in 1:npars
        dK[p] = dKij_dθp(kern, X, i, j, p, dim)
    end
end
@inline @inbounds function dKij_dθ!(dK::VecF64, kern::Kernel, X::MatF64, data::KernelData,
                                    i::Int, j::Int, dim::Int, npars::Int)
    for iparam in 1:npars
        dK[iparam] = dKij_dθp(kern, X, data, i, j, iparam, dim)
    end
end

function grad_slice!(dK::MatF64, k::Kernel, X::MatF64, data::KernelData, p::Int)
    dim = size(X,1)
    @inbounds for j in 1:size(X, 2)
        @simd for i in 1:j
            dK[i,j] = dKij_dθp(k,X,data,i,j,p,dim)
            dK[j,i] = dK[i,j]
        end
    end
    return dK
end

function grad_slice!(dK::MatF64, k::Kernel, X::MatF64, data::EmptyData, p::Int)
    dim = size(X, 1)
    @inbounds for j in 1:size(X, 2)
        @simd for i in 1:j
            dK[i,j] = dKij_dθp(k,X,i,j,p,dim)
            dK[j,i] = dK[i,j]
        end
    end
    return dK
end

# Calculates the stack [dk / dθᵢ] of kernel matrix gradients
function grad_stack!(stack::AbstractArray, k::Kernel, X::MatF64, data::KernelData)
    @inbounds for p in 1:num_params(k)
        grad_slice!(view(stack, :, :, p), k, X, data, p)
    end
    stack
end

grad_stack!(stack::AbstractArray, k::Kernel, X::MatF64) =
    grad_stack!(stack, k, X, KernelData(k, X))

grad_stack(k::Kernel, X::MatF64) = grad_stack(k, X, KernelData(k, X))

function grad_stack(k::Kernel, X::MatF64, data::KernelData)
    nobs = size(X, 2)
    stack = Array{Float64}(undef, nobs, nobs, num_params(k))
    grad_stack!(stack, k, X, data)
end

include("stationary.jl")
include("distance.jl")
include("lin.jl")               # Linear covariance function
include("se.jl")                # Squared exponential covariance function
include("rq.jl")                # Rational quadratic covariance function
include("mat.jl")               # Matern covariance function
include("periodic.jl")          # Periodic covariance function
include("poly.jl")              # Polnomial covariance function
include("noise.jl")             # White noise covariance function
include("const.jl")             # Constant (bias) covariance function

# Composite kernels
include("composite_kernel.jl")
include("pair_kernel.jl")
include("sum_kernel.jl")
include("prod_kernel.jl")

# Wrapped kernels
include("masked_kernel.jl")     # Masked kernels (apply to subset of X dims)
include("fixed_kernel.jl")      # Fixed kernels (fix some hyperparameters)


text:
using Distributions, Klara

p = BasicContUnvParameter(:p, pdf=Normal())

model = likelihood_model([p], isindexed=false)

mcsampler = AM(1.)

tuner = VanillaMCTuner()

mcrange = BasicMCRange(nsteps=10000, burnin=1000)

v0 = Dict(:p=>3.11)

outopts = Dict(:monitor=>[:value, :logtarget], :diagnostics=>[:accept])

job = BasicMCJob(model, mcsampler, mcrange, v0, outopts=outopts)

@time run(job)

chain = output(job)

mean(chain)

acceptance(chain)


text:
# AdaptivePenalization penalty
# ------------------------------------------------------------
@with_kw mutable struct AdaptiveData{T <: AbstractFloat}
    r::T = 1.0
end

@with_kw mutable struct AdaptiveParameters{T <: AbstractFloat}
    ζ::T = 0.1
    γ₁::T = 1e-5
    γ₂::T = 0.01
    γ₃::T = 0.25
    σ::T = 1e-5
    α::T = 0.95
    θ::T = 1.1
    ν::T = 0.1
    β::T = 1.1
    η::T = 1.25
end

"""
    AdaptivePenalization

Functor object for using adaptive penalty in a progressive-hedging algorithm. Create by supplying an [`Adaptive`](@ref) object through `penalty` in the `ProgressiveHedgingSolver` factory function and then pass to a `StochasticPrograms.jl` model.

...
# Parameters
- `ζ::T = 0.1`: Used to calculate the initial penalty. Non-anticipativity in the initial decision is enforced more as ζ increase.
- `γ₁::T = 1e-5`: Tolerance for primal changes being significant
- `γ₂::T = 0.01`: Tolerance for primal changes dominating dual changes
- `γ₃::T = 0.25`: Tolerance for dual changes dominating primal changes
- `σ::T = 1e-5`: Tolerance for the quadratic penalty dominating the Lagrangian
- `α::T = 0.95`: Penalty decrease after primal changes dominating dual changes
- `θ::T = 1.1`: Penalty increase after dual changes dominating primal changes
- `ν::T = 0.1`: Tolerance for significant non-anticipativity violation
- `β::T = 1.1`: Penalty increase after increased non-anticipativity violation
- `η::T = 1.25`: Default penalty increase in the default case
...
"""
struct AdaptivePenalization{T <: AbstractFloat} <: AbstractPenalization
    data::AdaptiveData{T}
    parameters::AdaptiveParameters{T}

    function AdaptivePenalization(r::AbstractFloat; kw...)
        T = typeof(r)
        return new{T}(AdaptiveData{T}(; r = r), AdaptiveParameters{T}(;kw...))
    end
end
function penalty(::AbstractProgressiveHedgingSolver, penalty::AdaptivePenalization)
    return penalty.data.r
end
function initialize_penalty!(ph::AbstractProgressiveHedgingSolver, penalty::AdaptivePenalization)
    update_dual_gap!(ph)
    @unpack δ₂ = ph.data
    @unpack ζ = penalty.parameters
    penalty.data.r = max(1., 2*ζ*abs(calculate_objective_value(ph)))/max(1., δ₂)
end
function update_penalty!(ph::AbstractProgressiveHedgingSolver, penalty::AdaptivePenalization)
    @unpack δ₁, δ₂ = ph.data
    @unpack r = penalty.data
    @unpack γ₁, γ₂, γ₃, σ, α, θ, ν, β, η = penalty.parameters

    δ₂_prev = length(ph.dual_gaps) > 0 ? ph.dual_gaps[end] : Inf

    μ = if δ₁/norm(ph.ξ,2)^2 >= γ₁
        if (δ₁-δ₂)/(1e-10 + δ₂) > γ₂
            α
        elseif (δ₂-δ₁)/(1e-10 + δ₁) > γ₃
            θ
        else
            1.
        end
    elseif δ₂ > δ₂_prev
        if (δ₂-δ₂_prev)/δ₂_prev > ν
            β
        else
            1.
        end
    else
        η
    end
    penalty.data.r = μ*r
end

# API
# ------------------------------------------------------------
"""
    Adaptive

Factory object for [`AdaptivePenalization`](@ref). Pass to `penalty` in the `ProgressiveHedgingSolver` factory function. See ?AdaptivePenalization for parameter descriptions.

"""
struct Adaptive{T <: AbstractFloat} <: AbstractPenalizer
    r::T
    parameters::Dict{Symbol,Any}

    function Adaptive(r::AbstractFloat; kw...)
        T = typeof(r)
        return new{T}(r, Dict{Symbol,Any}(kw))
    end
end
Adaptive(; kw...) = Adaptive(1.0; kw...)

function (adaptive::Adaptive)()
    return AdaptivePenalization(adaptive.r; adaptive.parameters...)
end

function str(::Adaptive)
    return "adaptive penalty"
end


text:
module Client

using PyCall: PyObject

using DockerPy: Collection, docker, @pyinterface
using DockerPy.Images: Image, pull, push, search
using DockerPy.Containers: Container, start, stop
using DockerPy.Volumes: Volume

import DockerPy.Images
import DockerPy.Containers

export DockerClient,
    images, containers, volumes, events, datausage, info, login, ping, version

mutable struct DockerClient
    o::PyObject
end
DockerClient(;
    base_url = nothing,
    version = "1.35",
    timeout = 60,
    tls = false,
    user_agent = "docker-sdk-python/$version",
    credstore_env = nothing,
) = DockerClient(docker.client.DockerClient(
    base_url = base_url,
    version = version,
    timeout = timeout,
    tls = tls,
    user_agent = user_agent,
    credstore_env = credstore_env,
))

images(x::DockerClient) = collect(x.images)
containers(x::DockerClient) = collect(x.containers)
volumes(x::DockerClient) = collect(x.volumes)
events(x::DockerClient) = PyObject(x).events()
datausage(x::DockerClient) = Dict{String,Any}(PyObject(x).df())
info(x::DockerClient) = Dict{String,Any}(PyObject(x).info())
login(
    x::DockerClient,
    username;
    password = nothing,
    email = nothing,
    registry = nothing,
    reauth = false,
    dockercfg_path = nothing,
) = PyObject(x).login(
    username,
    password = password,
    email = email,
    registry = registry,
    dockercfg_path = dockercfg_path,
)
ping(x::DockerClient) = PyObject(x).ping()
version(x::DockerClient) = Dict{String,Any}(PyObject(x).version())
Base.show(io::IO, x::DockerClient) = print(io, "DockerClient(\"$(objectid(x))\")")

Images.pull(x::DockerClient, repository; kwargs...) = pull(x.images, repository; kwargs...)
Images.push(x::DockerClient, repository; kwargs...) = push(x.images, repository; kwargs...)
Images.search(x::DockerClient, term) = search(x.images, term)
Base.rm(x::DockerClient, image::Image; kwargs...) = rm(x.images, image; kwargs...)
Containers.Container(x::DockerClient, image::Image; kwargs...) =
    Container(x.containers, image; kwargs...)

function Base.getproperty(x::DockerClient, name::Symbol)
    if name == :images
        return Collection{Image}(PyObject(x).images)
    elseif name == :containers
        return Collection{Container}(PyObject(x).containers)
    elseif name == :volumes
        return Collection{Volume}(PyObject(x).volumes)
    else
        return getproperty(PyObject(x), name)
    end
end # function Base.getproperty

@pyinterface DockerClient

end # module Client


text:
__precompile__()

module Halton

using Distributions
using Primes

export HaltonSeq!, HaltonDraws!

"""
    `HaltonSeq!{T<:AbstractFloat}(H::Vector{T}, B::Int; skip::Int=500)`

Replaces `H` with entries from Halton low discrepancy sequence with base `B`.
Elements in `H` take values in the interval (0, 1).
Keyword argument `skip` is the number initial "burn-in" elements to drop.
"""
function HaltonSeq!(H::AbstractArray{<:AbstractFloat}, B::Integer; skip::Integer=500)
  isprime(B) || error("base number not prime")
  H!(H, B, skip)
end

"""
    `HaltonDraws!{T<:AbstractFloat}(H::Vector{T}, B::Int, [skip::Int=500, distr=Normal()])`

Replaces `H` with draws from a distribution `Distributions.dist()`.
Draws are generated by using Halton sequence with base `B` as the quantiles drawn.
Keyword argument `skip` is the number initial "burn-in" elements to drop.
"""
function HaltonDraws!(H::AbstractArray, B::Integer; skip::Integer=500, distr = Normal())
  HaltonSeq!(H, B, skip=skip)
  H .= Distributions.quantile.(distr, H)
end

####

## Algorithm for generating Halton sequences
function H!(H::AbstractArray{T}, b::IT, skip::Integer) where {T<:AbstractFloat, IT<:Integer}
  # Fill H with Halton Sequence based on b
  S = skip + length(H)
  # set D to get generated seq >= S
  D = ceil(IT, log(S) / log(b))
  # placeholders
  d = zeros(T, D+1)
  r = zeros(T, D+1)

  # based on algorithm found in https://www.researchgate.net/publication/229173824_Fast_portable_and_reliable_algorithm_for_the_calculation_of_Halton_numbers
  for nn in 1:S
    ii = 1
    while d[ii] == b-1
      d[ii] = zero(T)
      ii += 1
    end
    d[ii] += 1
    if ii>=2
      r[ii - 1] = (d[ii] + r[ii]) / b
    end
    if ii >= 3
      for jj in (ii-1) : -1 : 2
        r[jj-1] = r[jj] / b
      end
    end
    if nn>skip
      H[nn-skip] = (d[1] + r[1]) / b
    end
  end
  return H
end

end


text:
# mbe_demo_c_test.jl
#
# This script tests the Julia implementation against the C implementation.

# Include our custom models.
include("models.jl")

# We work inside a module because it's more flexible than the main workspace.
module MBEDemoCTest

# Use/import the modules we'll need.
using ..MBEModels  # The MBEModels module included above (and "above" this module)
using Random
import Libdl

# Set some constants (we could randomize these too...).
I  = 5.
κc = 0.5
μc = 1.
ρ  = 2.
α  = 10.

# We'll do this many random tests.
n = 100

# Connect to the library containing our C version of the controller.
c_lib = Libdl.dlopen("windows/lib-reuc.dll")
c_fcn = Libdl.dlsym(c_lib, :reuc)

try

    # Creat a place to store the differences.
    Δτ = zeros(2, n)

    # Seed the random number generator so that this test is repeatable.
    Random.seed!(1)

    for k = 1:n

        # Create a random state and target.
        q_TI   = normalize(randn(4))
        q_BI   = normalize(randn(4))
        ω_BI_B = [randn(1); randn(1); 0.]

        # Run the Julia implementation.
        τ_B = MBEModels.rate_controller(I, κc, μc, ρ, α, q_TI, q_BI, ω_BI_B)

        # Run the C implementation.
        τ_B_c = [0.; 0.; 0.]
        ccall(c_fcn, # Function to call
            Nothing, # No return value
            # Function interface:
            (Float64, Float64, Float64, Float64, Float64, # Parameters
            Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, # Quaternions and rate
            Ptr{Float64}), # Output torque
            # Actual values to pass:
            I, κc, μc, ρ, α,
            q_TI, q_BI, ω_BI_B,
            τ_B_c)

        # Store the differences.
        Δτ[:,k] = τ_B[1:2] - τ_B_c[1:2]

    end

    display(Δτ)
    println("Maximum difference in torque: ", maximum(Δτ[:]), " Nm")

finally
    Libdl.dlclose(c_lib) # Close out the library.
end

end # module


text:
@testset "Building a frame graph" begin
    frame = Frame(device)
    fg = FrameGraph(device, frame)

    add_pass!(identity, fg, :gbuffer, RenderPass((0,0,1920,1080)))
    add_pass!(identity, fg, :lighting, RenderPass((0,0,1920,1080)))
    add_pass!(identity, fg, :adapt_luminance, RenderPass((0,0,1920,1080)))
    add_pass!(identity, fg, :combine, RenderPass((0,0,1920,1080)))
    # can't add a pass more than once
    @test_throws ErrorException add_pass!(identity, fg, :combine, RenderPass((0,0,1920,1080)))

    add_resource!(fg, :vbuffer, BufferResourceInfo(1024))
    add_resource!(fg, :ibuffer, BufferResourceInfo(1024))
    add_resource!(fg, :average_luminance, ImageResourceInfo(Vk.FORMAT_R32G32B32A32_SFLOAT))
    add_resource!(fg, :emissive, AttachmentResourceInfo(Vk.FORMAT_R32G32B32A32_SFLOAT))
    add_resource!(fg, :albedo, AttachmentResourceInfo(Vk.FORMAT_R32G32B32A32_SFLOAT))
    add_resource!(fg, :normal, AttachmentResourceInfo(Vk.FORMAT_R32G32B32A32_SFLOAT))
    add_resource!(fg, :pbr, AttachmentResourceInfo(Vk.FORMAT_R32G32B32A32_SFLOAT))
    add_resource!(fg, :color, AttachmentResourceInfo(Vk.FORMAT_R32G32B32A32_SFLOAT))
    add_resource!(fg, :output, AttachmentResourceInfo(Vk.FORMAT_R32G32B32A32_SFLOAT))
    add_resource!(fg, :depth, AttachmentResourceInfo(Vk.FORMAT_D32_SFLOAT))
    # can't add a resource more than once
    @test_throws ErrorException add_resource!(fg, :depth, AttachmentResourceInfo(Vk.FORMAT_D32_SFLOAT))

    # imported
    add_resource!(fg, :shadow_main, AttachmentResourceInfo(Vk.FORMAT_D32_SFLOAT))
    add_resource!(fg, :shadow_near, AttachmentResourceInfo(Vk.FORMAT_D32_SFLOAT))
    add_resource!(fg, :bloom_downsample_3, AttachmentResourceInfo(Vk.FORMAT_R32G32B32A32_SFLOAT))

    usages = @resource_usages begin
        emissive::Color, albedo::Color, normal::Color, pbr::Color, depth::Depth = gbuffer(vbuffer::Buffer::Vertex, ibuffer::Buffer::Index)
        color::Color = lighting(emissive::Color, albedo::Color, normal::Color, pbr::Color, depth::Depth, shadow_main::Texture, shadow_near::Texture)
        average_luminance::Image::Storage = adapt_luminance(average_luminance::Image::Storage, bloom_downsample_3::Texture)
        output::Color = combine(color::Color, average_luminance::Texture)
    end

    add_resource_usage!(fg, usages)
    Lava.resolve_attributes!(fg)

    @test Lava.buffer_usage(fg, :vbuffer) == Vk.BUFFER_USAGE_VERTEX_BUFFER_BIT
    @test Lava.image_usage(fg, :depth) == Vk.IMAGE_USAGE_DEPTH_STENCIL_ATTACHMENT_BIT
    @test Lava.image_usage(fg, :shadow_main) == Vk.IMAGE_USAGE_SAMPLED_BIT

    for resource in keys(fg.resources)
        @test Int(Lava.resource_attribute(fg, resource, :usage)) ≠ 0
    end
end

function color_attachment(device)
    usage = Vk.IMAGE_USAGE_COLOR_ATTACHMENT_BIT | Vk.IMAGE_USAGE_TRANSFER_SRC_BIT | Vk.IMAGE_USAGE_TRANSFER_DST_BIT
    color_image = allocate!(ImageBlock(device, (1920, 1080), Vk.FORMAT_R16G16B16A16_SFLOAT, usage), MEMORY_DOMAIN_DEVICE)
    color_attachment = Attachment(View(color_image), WRITE)
end

function add_color_attachment(fg::FrameGraph)
    color_attachm = color_attachment(fg.device)
    register(fg.frame, :color, color_attachm)
    add_resource!(fg, :color, AttachmentResourceInfo(Lava.format(color_attachm)))
end

function save_test_render(filename, data, h::UInt; tmp = false)
    filename = render_file(filename; tmp)
    ispath(filename) && rm(filename)
    save(filename, data')
    @test stat(filename).size > 0
    @test hash(data) == h
end

function program_1(device, vdata)
    prog = Program(device, ShaderSpecification(shader_file("rectangle.vert"), GLSL), ShaderSpecification(shader_file("rectangle.frag"), GLSL))

    fg = FrameGraph(device)
    add_color_attachment(fg)
    add_pass!(fg, :main, RenderPass((0,0,1920,1080))) do rec
        set_program(rec, prog)
        ds = draw_state(rec)
        set_draw_state(rec, @set ds.program_state.primitive_topology = Vk.PRIMITIVE_TOPOLOGY_TRIANGLE_STRIP)
        draw(rec, TargetAttachments([:color]), vdata, collect(1:4); alignment = 8)
    end

    usage = @resource_usages begin
        color::Color = main()
    end
    add_resource_usage!(fg, usage)
    clear_attachments(fg, :main, [:color => (0.08, 0.05, 0.1, 1.)])
    fg
end

function program_2(device, vdata)
    prog = Program(device, ShaderSpecification(shader_file("texture.vert"), GLSL), ShaderSpecification(shader_file("texture.frag"), GLSL))

    fg = FrameGraph(device)
    add_color_attachment(fg)

    normal = load(texture_file("normal.png"))
    normal = convert(Matrix{RGBA{Float16}}, normal)
    normal_map = image(device, normal, Vk.FORMAT_R16G16B16A16_SFLOAT; usage = Vk.IMAGE_USAGE_SAMPLED_BIT)
    register(fg.frame, :normal_map, normal_map)
    add_resource!(fg, :normal_map, ImageResourceInfo(Lava.format(normal_map)))

    add_pass!(fg, :main, RenderPass((0,0,1920, 1080))) do rec
        set_program(rec, prog)
        ds = draw_state(rec)
        set_draw_state(rec, @set ds.program_state.primitive_topology = Vk.PRIMITIVE_TOPOLOGY_TRIANGLE_STRIP)
        set_material(rec,
            (0.1f0, 1f0), # uv scaling coefficients
            Texture(:normal_map, setproperties(DEFAULT_SAMPLING, (magnification = Vk.FILTER_LINEAR, minification = Vk.FILTER_LINEAR))),
        )
        draw(rec, TargetAttachments([:color]), vdata, collect(1:4))
    end

    usage = @resource_usages begin
        color::Color = main(normal_map::Texture)
    end
    add_resource_usage!(fg, usage)
    clear_attachments(fg, :main, [:color => (0.08, 0.05, 0.1, 1.)])
    fg
end

function program_3(device, positions, ppm)
    prog = Program(device, ShaderSpecification(shader_file("glyph.vert"), GLSL), ShaderSpecification(shader_file("glyph.frag"), GLSL))

    fg = FrameGraph(device)
    add_color_attachment(fg)

    font = OpenTypeFont(font_file("juliamono-regular.ttf"))
    curves = map(OpenType.curves(font['C'])) do curve
        map(curve) do points
            map(Float32, points)
        end
    end
    curve_buffer = buffer(device, curves)
    vdata = [(pos..., UInt32(0), UInt32(length(curves))) for pos in positions]
    register(fg.frame, :curve_buffer, curve_buffer)

    add_pass!(fg, :main, RenderPass((0,0,1920, 1080))) do rec
        set_program(rec, prog)
        ds = draw_state(rec)
        set_draw_state(rec, @set ds.program_state.primitive_topology = Vk.PRIMITIVE_TOPOLOGY_TRIANGLE_STRIP)
        set_material(rec,
            RGBA{Float32}(1., 1., 0., 1.), # text color
            ppm, # pixel per em
            device_address(curve_buffer);
            alignment = 8
        )
        draw(rec, TargetAttachments([:color]), vdata, collect(1:4); alignment = 8)
    end

    usage = @resource_usages begin
        color::Color = main()
    end
    add_resource_usage!(fg, usage)
    clear_attachments(fg, :main, [:color => (0.08, 0.05, 0.1, 1.)])
    fg
end

@testset "Rendering" begin
    vdata = [
        (-0.5f0, -0.5f0, RGB{Float32}(1., 0., 0.)),
        (0.5f0, -0.5f0, RGB{Float32}(1., 1., 1.)),
        (-0.5f0, 0.5f0, RGB{Float32}(0., 1., 0.)),
        (0.5f0, 0.5f0, RGB{Float32}(0., 0., 1.)),
    ]
    fg = program_1(device, vdata)
    snoop = Lava.SnoopCommandBuffer()
    render(fg; command_buffer = snoop, submit = false)
    @test length(snoop.records) == 8
    @test fg.frame.gd.index_list == [1, 2, 3, 4]
    ib = collect(UInt32, fg.frame.gd.index_buffer[])
    @test ib == UInt[0, 1, 2, 3] && sizeof(ib) == 16
    @test fg.frame.gd.allocator.last_offset == 24 * 4

    fg = program_1(device, vdata)
    @test wait(render(fg))
    data = collect(RGBA{Float16}, image(fg.frame.resources[:color].data), device)
    save_test_render("colored_rectangle.png", data, 0x9430efd8e0911300)

    vdata = [
        (-0.5f0, -0.5f0, 0f0, 0f0),
        (0.5f0, -0.5f0, 1f0, 0f0),
        (-0.5f0, 0.5f0, 0f0, 1f0),
        (0.5f0, 0.5f0, 1f0, 1f0),
    ]
    fg = program_2(device, vdata)
    @test wait(render(fg))
    data = collect(RGBA{Float16}, image(fg.frame.resources[:color].data), device)
    save_test_render("distorted_normal_map.png", data, 0x9eda4cb9b969b269)

    positions = [
        (-0.8f0, -0.8f0),
        (0.8f0, -0.8f0),
        (-0.8f0, 0.8f0),
        (0.8f0, 0.8f0),
    ]
    fg = program_3(device, positions, 12f0)
    @test wait(render(fg))
    # FIXME: debug why font shader fails
    # data = collect(RGBA{Float16}, image(fg.frame.resources[:color].data), device)
    # save_test_render("glyph_A.png", data, hash(data))
end


text:
using SynchrotronKernel
#using QuadGK
global const X_MAX = 70.0  # -> F(x) < 2e-6
global const X_MIN = 1e-20 # -> F(x) < 5e-6

"""
    min_synch_loud_momentum(ν0::Real, B::Real)

Computes the minimum synchrotron loud momentum as described in Donnert+16, MNRAS 462, 2014–2032 (2016), Eq. 22.
"""
min_synch_loud_momentum(ν0::Real, B::Real) = √( ν0 / (C_crit * B)  )

"""
    ν_over_ν_crit(p::Real, B::Real, ν::Real)

Computes the fraction ``x = \frac{ν}{ν_c}`` needed for the synchrotron Kernel.
See Donnert+16, MNRAS 462, 2014–2032 (2016), Eq. 19.
"""
ν_over_ν_crit(p::Real, B::Real, ν::Real, sinθ::Real=1.0) = ν / ( C_crit * B * sinθ * p^2 )

"""
    dsa_spectral_index_momentum(M::Real)

Compute the spectral index from DSA for a distribution in momentum space.
"""
function dsa_spectral_index_momentum(M::Real)
    xs = ( γ + 1 ) * M^2 / ( ( γ - 1 ) * M^2 + 2 )
    return 3.0 * xs / (xs - 1.0 )
end

"""
    spectrum(p::Real, p_min::Real, n0::Real, s::Real)

Number density `n` for a given momentum `p` interpolated from `n0` at minimum momentum `p_min`.
"""
spectrum(p::Real, p_min::Real, n0::Real, s::Real) = p < p_min ? 0.0 : n0 * (p/p_min)^(-s)

"""
    integrate_θ(x_in::Real, θ_steps::Integer=100)

Pitch angle integration in Donnert+16, Eq. 17.
"""
function integrate_θ(x_in::Real, θ_steps::Integer=64)

    dθ = 0.5π / θ_steps
    K = 0.0
    @inbounds for θ ∈ LinRange(0.0, 0.5π, θ_steps)
        sinθ = sin(θ)
        x  = x_in / sinθ
        K += sinθ^2 * synchrotron_kernel(x) * dθ
    end
    return K
end

"""
    integrate_θ(x_in::Real, θ_steps::Integer=100)

Pitch angle integration in Donnert+16, Eq. 17 using Simpson's rule.
"""
function integrate_θ_simpson(x_in::Real, θ_steps::Integer=50)

    dθ = 0.5π / θ_steps
    #θ = LinRange(0.0, 0.5π, θ_steps)

    F     = Vector{Float64}(undef, θ_steps)
    F_mid = Vector{Float64}(undef, θ_steps)

    # sin(0.0) = 0.0
    F[1]     = 0.0
    F_mid[1] = sin(0.5dθ)^2 * synchrotron_kernel(x_in/sin(0.5dθ))
    K = 0.0
    @inbounds for i = 2:θ_steps
        # boundary
        sinθ = sin((i-1)*dθ)
        x  = x_in / sinθ
        F[i] = sinθ^2 * synchrotron_kernel(x)

        # mid point
        sinθ = sin((i-0.5)*dθ)
        x    = x_in / sinθ
        F_mid[i] = sinθ^2 * synchrotron_kernel(x)

        # Simpson rule: https://en.wikipedia.org/wiki/Simpson%27s_rule
        K += dθ / 6.0 * ( F[i] + F[i-1] + 4F_mid[i] )
    end
    return K
end

"""
    get_F( p::Real, p_inj::Real, ϵ_cr0::Real, s::Real, B_cgs::Real, ν0::Real)

Helper function to get synchrotron emission at momentum `p`.
"""
function get_F( p::Real, p_inj::Real, ϵ_cr0::Real, s::Real, B_cgs::Real, ν0::Real, integrate_pitch_angle::Bool)

    # energy density at momentum p
    N_E = spectrum(p, p_inj, ϵ_cr0, s)
    # x from Eq. 19 (without sinθ -> for integration)
    x   = ν_over_ν_crit(p, B_cgs, ν0)

    if integrate_pitch_angle
        K = integrate_θ_simpson(x)
    else
        K = synchrotron_kernel(x)
    end

    # integral Donnert+16, Eq. 17 evaluated at momenum p
    return N_E * K

end

"""
    get_log_mid(p_low::Real, p_up::Real)

Helper function to compute the momentum in the middle of the bin.
"""
get_log_mid(q_low::Real, q_up::Real) = 10.0^( 0.5 * ( log10(q_low) + log10(q_up) ) )

"""
    spectral_synchrotron_emission(rho_cgs::Real, B_cgs::Real,
                                  T_K::Real, Mach::Real;
                                  xH::Real=0.76, dsa_model::Integer=1,
                                  ν0::Real=1.44e9,
                                  K_ep::Real=0.01,
                                  Emin::Real=5.0e4,
                                  Emax::Real=1.e10,
                                  p_inj::Real=0.1, # in [me*c]
                                  integrate_pitch_angle::Bool=true,
                                  convert_to_mJy::Bool=false,
                                  N_sample_bins::Integer=128)

Computes the synchrotron emission for a powerlaw spectrum as described in Donnert+16, MNRAS 462, 2014–2032 (2016).

# Arguments
- `rho_cgs::Array{<:Real}`: Density in ``g/cm^3``.
- `B_cgs::Array{<:Real}`:   Magnetic field in Gauss.
- `T_K::Array{<:Real}`:     Temperature in Kelvin.
- `Mach::Array{<:Real}`:    Mach number.

# Keyword Arguments
- `xH::Float64 = 0.76`:               Hydrogen fraction of the simulation, if run without chemical model.
- `dsa_model::Integer=1`:             Diffuse-Shock-Acceleration model. Takes values `0...4`, see next section.
- `ν0::Real=1.44e9`:                  Observation frequency in ``Hz``.
- `K_ep::Real=0.01`:                  Ratio of CR proton to electron energy density.
- `integrate_pitch_angle::Bool=true`: Explicitly integrates over the pitch angle. IF `false` assumes ``sin(θ) = 1``.
- `convert_to_mJy::Bool=false`:       Convert the result from ``[erg/cm^3/Hz/s]`` to ``mJy/cm``.

# DSA models
- `0`: [`KR07_acc`](@ref). Efficiency model from Kang, Ryu, Cen, Ostriker 2007, http://arxiv.org/abs/0704.1521v1.
- `1`: [`KR13_acc`](@ref). Efficiency model from Kang&Ryu 2013, doi:10.1088/0004-637X/764/1/95 .
- `2`: [`Ryu19_acc`](@ref). Efficiency model from Ryu et al. 2019, https://arxiv.org/abs/1905.04476 .
- `3`: [`CS14_acc`](@ref). Efficiency model from Caprioli&Spitkovsky 2015, doi: 10.1088/0004-637x/783/2/91 .
- `4`: [`P16_acc`](@ref). Constant efficiency as in Pfrommer+ 2016, doi: 10.1093/mnras/stw2941 .

"""
function spectral_synchrotron_emission(rho_cgs::Real, B_cgs::Real,
                                       T_K::Real, Mach::Real;
                                       xH::Real=0.76, dsa_model::Integer=1,
                                       ν0::Real=1.44e9,
                                       K_ep::Real=0.01,
                                       Emin::Real=5.0e4,
                                       Emax::Real=1.e10,
                                       p_inj::Real=0.1, # in [me*c]
                                       integrate_pitch_angle::Bool=true,
                                       convert_to_mJy::Bool=false,
                                       N_sample_bins::Integer=128)

    if ( B_cgs == 0 || Mach < 2 )
        return 0
    end

    if dsa_model == 0
        acc_function = KR07_acc
    elseif dsa_model == 1
        acc_function = KR13_acc
    elseif dsa_model == 2
        acc_function = Ryu19_acc
    elseif dsa_model == 3
        acc_function = CS14_acc
    elseif dsa_model == 4
        acc_function = P16_acc
    else
        error("Invalid DSA model selection!")
    end

    # in dimensionless momentum!
    p_inj = 0.1
    p_min = Emin * eV2cgs / ( m_e * c_light^2 )
    p_max = Emax * eV2cgs / ( m_e * c_light^2 )

    # prefactor to Eq. 17
    j_ν_prefac = q_e * √(3) / (m_e * c_light^2)

    # include magnetic field into this
    j_ν_prefac *= B_cgs

    # if run without pitch angle integration
    # sinθ = 1.0 -> integral factor π/2
    if !integrate_pitch_angle
        j_ν_prefac *= 0.5π
    end

    # get spectral index s ∈ ]-∞, 4.0[
    s  = dsa_spectral_index_momentum(Mach)

    # energy density of thermal gas [erg/cm^3]
    ϵ_th = EpsNtherm(rho_cgs, T_K, xH=xH)

    # # CR Energy density of whole powerlaw
    # ϵ_cr0 = K_ep * get_rel_energy_density(Mach, acc_function) * ϵ_th
    # # compensate for injection momentum
    # ϵ_cr0 *= (s - 2) / (p_inj^( 2 - s ) )
    # ! This gives the same result as in Smac2

    # this gives the same result as the analytic solution
    ϵ_cr0 = cre_spec_norm_particle(Mach, acc_function) * ϵ_th

    # width of momentum bins
    di    = log10(p_max/p_min) / (N_sample_bins-1)

    # momentum bins
    p     = @. p_min * 10.0^( di * 0:N_sample_bins-1)

    # width of momentum bins
    dp = zeros(N_sample_bins)
    dp[1] = p[1] - p_min * 10.0^(-1di)
    @inbounds for i = 2:N_sample_bins
        dp[i] = p[i] - p[i-1]
    end

    # storage array for integrand of Eq. 17
    F        = Vector{Float64}(undef, N_sample_bins)
    F_mid    = Vector{Float64}(undef, N_sample_bins)

    F[1]     = get_F(p[1], p_inj, ϵ_cr0, s, B_cgs, ν0, integrate_pitch_angle)
    # not needed
    F_mid[1] = 0.0

    # store total synchrotron emissivity
    jν = 0.0

    # Simpson rule
    @inbounds for i = 2:N_sample_bins

        # end of bin
        F[i]     = get_F(p[i], p_inj, ϵ_cr0, s, B_cgs, ν0, integrate_pitch_angle)

        # middle of bin
        p_mid    = get_log_mid(p[i-1], p[i])
        F_mid[i] = get_F(p_mid, p_inj, ϵ_cr0, s, B_cgs, ν0, integrate_pitch_angle)

        # Simpson rule: https://en.wikipedia.org/wiki/Simpson%27s_rule
        jν += dp[i] / 6.0 * ( F[i] + F[i-1] + 4F_mid[i] )

    end

    if convert_to_mJy
        j_ν_prefac *= mJy_factor
    end

    return j_ν_prefac * jν

end

"""
    spectral_synchrotron_emission( n_p::Vector{<:Real},
                                   p::Vector{<:Real},
                                   B_cgs::Real;
                                   ν0::Real=1.44e9
                                   integrate_pitch_angle::Bool=false,
                                   convert_to_mJy::Bool=false)

Computes the synchrotron emission (in ``[erg/cm^3/Hz/s]``) for a CR spectrum as described in Donnert+16, MNRAS 462, 2014–2032 (2016), Eq. 17.

# Arguments
- `n_p::Vector{<:Real}`: Number density in ``1/cm^3`` for momenta `p`.
- `p::Vector{<:Real}`:   Momenta `p` for number densities.
- `B_cgs::Real`:         Magnetic field strength (absolute value) in Gauss.

# Keyword Arguments
- `ν0::Real=1.44e9`:                  Observation frequency in ``Hz``.
- `integrate_pitch_angle::Bool=false`: Explicitly integrates over the pitch angle. If `false` assumes ``sin(θ) = 1``.
- `convert_to_mJy::Bool=false`:       Convert the result from ``[erg/cm^3/Hz/s]`` to ``mJy/cm``.

"""
function spectral_synchrotron_emission(n_p::Vector{<:Real},
                                       p::Vector{<:Real},
                                       B_cgs::Real;
                                       ν0::Real=1.44e9,
                                       integrate_pitch_angle::Bool=false,
                                       convert_to_mJy::Bool=false)

    if ( B_cgs == 0 ) || sum(n_p) == 0
        return 0
    end

    # avoid error for n_p = -0.0
    n_p[ n_p .== -0.0 ] .= 0.0

    # if there are Infs or NaNs present we get a wrong result, so skip this particle
    # ToDo: Check if there is a better way to fix this!
    if (length(findall( isnan.(n_p) )) > 0) ||
       (length(findall( isinf.(n_p) )) > 0) ||
       (length(findall(   n_p .< 0  )) > 0)
        return 0
    end

    # get the number of momentums for which the energy density is defined
    Nbins = length(n_p)

    # prefactor to Eq. 17
    j_ν_prefac = q_e * √(3) / (m_e * c_light^2)

    # include magnetic field into this
    j_ν_prefac *= B_cgs

    # if run without pitch angle integration
    # sinθ = 1.0 -> integral factor π/2
    if !integrate_pitch_angle
        j_ν_prefac *= 0.5π
    end

    # width of momentum bins
    dp = Vector{Float64}(undef, Nbins)
    @inbounds for i = 1:Nbins
        dp[i] = p[i+1] - p[i]
    end

    # storage array for integrand of Eq. 17
    F        = Vector{Float64}(undef, Nbins)
    F_mid    = Vector{Float64}(undef, Nbins-1)

    @inbounds for i = 1:Nbins-1

        # beginning of bin

        # x from Eq. 19
        x   = ν_over_ν_crit(p[i], B_cgs, ν0)
        # pitch-angle integral
        if integrate_pitch_angle
            K = integrate_θ_simpson(x)
        else
            K = synchrotron_kernel(x)
        end

        # energy density at momentum p * integrated synchrotron kernel
        F[i] = n_p[i] * K

        # middle of bin
        N_E_mid  = get_log_mid( n_p[i],  n_p[i+1] )

        p_mid = get_log_mid( p[i], p[i+1] )
        x     = ν_over_ν_crit(p_mid, B_cgs, ν0)

        if integrate_pitch_angle
            K_mid = integrate_θ_simpson(x)
        else
            K_mid = synchrotron_kernel(x)
        end

        F_mid[i] = N_E_mid * K_mid

    end

    # last bin seperate
    x   = ν_over_ν_crit(p[end], B_cgs, ν0)
    # pitch-angle integral
    if integrate_pitch_angle
        K = integrate_θ_simpson(x)
    else
        K = synchrotron_kernel(x)
    end

    # energy density at momentum p * integrated synchrotron kernel
    F[end] = n_p[end] * K

    # store total synchrotron emissivity
    jν = 0.0
    @inbounds for i = 1:Nbins-1
        # Simpson rule: https://en.wikipedia.org/wiki/Simpson%27s_rule
        jν += dp[i] / 6.0 * ( F[i] + F[i+1] + 4F_mid[i] )
    end

    if convert_to_mJy
        j_ν_prefac *= mJy_factor
    end

    return j_ν_prefac * jν

end


text:
function
directRP!(
    RP::Vector,
    y::Vector,
    m::Int,
    h::Real,
    N::Int,
    ε::Real,
    ulₘ::Vector,
    urₘ::Vector
)
    # Функция для внутреннего использования, поэтому входные массивы размерности N-1
    @assert length(y)   == N-1  "Вектор искомой функции размерности N-1"
    @assert length(RP)  == N-1  "Вектор функции правой части размерности N-1"

    @assert length(ulₘ) > m > 0
    @assert length(urₘ) > m > 0

    @assert N > 4
    @assert ε > 0

    RP[1] = ε * (y[2] - 2 * y[1] + ulₘ[m]) / h^2 +
        y[1] * (y[2] - ulₘ[m]) / (2 * h) - y[1];

    for n in 2:N-2
        RP[n] = ε * (y[n+1] - 2 * y[n] + y[n-1]) / h^2 +
            y[n] * (y[n+1] - y[n-1]) / (2 * h) - y[n];
    end

    RP[N-1] = ε * (urₘ[m] - 2 * y[N-1] + y[N-2]) / h^2 +
        y[N-1] * (urₘ[m] - y[N-2]) / (2 * h) - y[N-1];

    return RP;
end

function
DRP!(
    dl::Vector,
    dd::Vector,
    du::Vector,
    y::Vector,
    m::Int,
    h::Real,
    N::Int,
    ε::Real,
    ulₘ::Vector,
    urₘ::Vector
)
    # Функция для внутреннего использования, поэтому входные массивы размерности N-1
    @assert length(y)   == N-1  "Вектор искомой функции размерности N-1"

    @assert length(ulₘ) > m > 0
    @assert length(urₘ) > m > 0

    @assert m > 0
    @assert h > 0
    @assert N > 3
    @assert ε > 0

    @assert length(dl) == N-2   "Вектор поддиагональных элементов размерностью N-2"
    @assert length(dd) == N-1   "Вектор    диагональных элементов размерностью N-1"
    @assert length(du) == N-2   "Вектор поддиагональных элементов размерностью N-2"

    dd[1] = -2 * ε / h^2 + (y[2] - ulₘ[m]) / (2 * h) - 1;
    du[1] = ε / h^2 - y[1] / (2 * h);

    for n in 2:N-2
        dl[n - 1] = ε / h^2 - y[n] / (2 * h);
        dd[n] = - 2 * ε / h^2 + (y[n+1] - y[n-1])/ (2 * h) - 1;
        du[n] = ε / h^2 + y[n] / (2 * h);
    end

    dl[N-2] = ε / h^2 + y[N-1] / (2 * h);
    dd[N-1] = - 2 * ε / h^2 + (urₘ[m] - y[N-2])/ (2 * h) - 1;

    return dl, dd, du;
end

function
solve!(
    u::Matrix,
    u₀::Vector,
    N::Int,
    h::Real,
    Tₘ::Vector,
    M::Int,
    τ::Real,
    ε::Real,
    ulₘ::Vector,
    urₘ::Vector;
    α::Complex = complex(0.5, 0.5))

    if size(u) != (N+1, M+1)
        throw(ArgumentError("""size(u) == $(size(u)), N == $(N), M == $(M)
        Матрица `u` должена иметь размерность (N+1, M+1)."""))
    end

    if length(u₀) != N+1
        throw(ArgumentError("""length(u₀) == $(length(u₀)), N == $(N)
        Массив `u₀` должен иметь размерность N+1."""))
    end

    if length(Tₘ) != M+1
        throw(ArgumentError("""length(Tₘ) == $(length(Tₘ)), M == $(M)
        Массив `Tₘ` должен иметь размерность M+1"""))
    end
    if length(ulₘ) != M+1
        throw(ArgumentError("""length(ulₘ) == $(length(ulₘ)), M == $(M)
        Массив `ulₘ` должен иметь размерность M+1"""))
    end
    if length(urₘ) != M+1
        throw(ArgumentError("""length(y) == $(length(y)), M == $(M)
        Массив `ulₘ` должен иметь размерность M+1"""))
    end

    u[1  , :]   = ulₘ;                  # Запись левого  ГУ на каждом шаге
    u[N+1, :]   = urₘ;                  # Запись правого ГУ на каждом шаге
    u[:, 1] = u₀;                       # Запись начального условия, это перепишет ГУ.

    y = u₀[2:end-1];
    rp = zero(y);
    dl = zeros(N-2);
    dd = zeros(N-1);
    du = zeros(N-2);
    j = Tridiagonal(zeros(N-2), zeros(N-1), zeros(N-2));

    for m in 1:M

        directRP!(rp, y, m, h, N, ε, ulₘ, urₘ);
        DRP!(dl, dd, du, y, m, h, N, ε, ulₘ, urₘ);

        j = Tridiagonal(dl, dd, du);

        W = (I - α  * τ * j) \ rp;
        y = y .+ τ * real(W);

        u[2:N, m+1] = y;
    end

    return u;
end


text:
using CSJLab
using Test

using FileIO
using Images

@testset "CSJLab.jl" begin
    # Write your own tests here.
end


text:
export  excelNodeCompositionReport

function excelNodeCompositionReport( mpSim::MPsim, timeGrid::Vector{T},
    nodes::AbstractString...; filename::AbstractString = "compReport",
    overwrite::Bool = true, timeFactor::Real = 12.0 ) where T <: Real

    filename, overwrite = setupFile( filename, overwrite )
    tStart = now()
    compReport = nodeCompositionReport( mpSim, timeGrid, nodes... )
    reportGenerationTime = ( now() - tStart ).value / 1000.0

    if isempty( compReport )
        return
    end  # if isempty( compReport )

    writeNodeCompositionData( compReport, nodes, mpSim, timeFactor,
        reportGenerationTime, filename, overwrite )

end  # excelNodeCompositionReport( mpSim, timeGrid, nodes, filename, overwrite,
     #   timeFactor )

excelNodeCompositionReport( mpSim::MPsim, timeRes::Real,
    nodes::AbstractString...; filename::AbstractString = "compReport",
    overwrite::Bool = true, timeFactor::Real = 12.0 ) =
    excelNodeCompositionReport( mpSim, MP.generateTimeGrid( mpSim, timeRes ),
    nodes...; filename = filename, overwrite = overwrite,
    timeFactor = timeFactor )

include( "private/nodecomposition.jl" )


text:
# using PMDPs
# using Test

using StaticArrays
using Distributions
using PMDPs.CountingProcesses

@testset "PricingProblem.jl" begin
    pp = simple_pp()

    @test isa(pp, PMDPs.PMDPProblem)
    @test PMDPs.selling_period_end(pp)==8

    @test size(pp) == (3,2,5)
end


text:
"""
- cdfE calculated empirical cumulative distribution function F(y<=yi)
- To get the CCDF of Y (complementary cumulative distribution function):
	- ccdf(Y) = 1- cdf(Y);
Example:
    [cdfY,y] = cdfE(Y,nbins); % nbins: number of bins (i.e. number of y points)
    [cdfY,y] = cdfE(Y,y); % y is the mid-points of the edges of bins. % This is expected to be the fastest.
                                               % Noted that the 'EdgeScale' has no effect on output results since y is specified.
    [cdfY,y] = cdfE(Y);
"""
function cdfE(Y; nedges::Int=100)
	edges(l,r,pt) = collect(range(l,r,length = pt));
	y = edges(minimum(Y), maximum(Y), nedges);
	cdfY = cdfE_core(y, Y, nedges);
	(y, cdfY);
end

function cdfE(Y, y::Vector; edgescale=:linear)
	nedges = length(y);
	cdfY = cdfE_core(y, Y, nedges)
end

function cdfE_core(y, Y, nedges)
	Y = [Y...];
	cdfY = fill(NaN,nedges)
	invlenY = 1/length(Y);
	for i = 1:nedges
		cdfY[i] = sum(Y .<= y[i])*invlenY;
	end
	cdfY
end


text:
"""
    zigzag_ordering(array::Array{Int64,1})

`zigzag_ordering` reorders an unsorted array in the zigzag fashion. This procedure implies
that between every second step a break will be generated by switching the `flag-boolean`.

# Arguments
- `array::Array{Int64,1}`: Unsorted array

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr = [4, 3, 7, 8, 6, 2, 1, 10, 13, 3]
julia> ClassicAlgorithmsCollections.zigzag_ordering(arr)
[3, 7, 4, 8, 2, 6, 1, 13, 3, 10]
```
"""
function zigzag_ordering(array::Array{Int64,1})
    # Intialize the flag
    flag = true
    for i in 1:length(array)-1
        # If flag is true, then check if next array element is greater than previous
        if flag
            if array[i] > array[i+1]
                array[i], array[i+1] = array[i+1], array[i]
            end
            # If flag is true, then check if next array element is smaller than previous
        else
            if array[i] < array[i+1]
                array[i], array[i+1] = array[i+1], array[i]
            end
        end
        # Swapt the flag after every turn
        flag = !flag
    end
    return array
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
# VectorAffineDecisionFunction #
# ========================== #
mutable struct VectorAffineDecisionConstraintBridge{T, S} <: MOIB.Constraint.AbstractBridge
    constraint::CI{MOI.VectorAffineFunction{T}, S}
    decision_function::VectorAffineDecisionFunction{T}
    set::S
end

function MOIB.Constraint.bridge_constraint(::Type{VectorAffineDecisionConstraintBridge{T, S}},
                                           model,
                                           f::VectorAffineDecisionFunction{T},
                                           set::S) where {T, S}
    # All decisions have been mapped to either the decision part constant or the variable part terms
    # at this point.
    F = MOI.VectorAffineFunction{T}
    # Calculate total constant
    constants = f.variable_part.constants +
        f.known_part.constants
    # Add the bridged constraint
    constraint = MOI.add_constraint(model,
                                    MOI.VectorAffineFunction(f.variable_part.terms, constants),
                                    set)
    # Save the constraint index, the decision function, and the set, to allow modifications
    return VectorAffineDecisionConstraintBridge{T, S}(constraint, f, set)
end

function MOI.supports_constraint(::Type{<:VectorAffineDecisionConstraintBridge{T}},
                                 ::Type{<:VectorAffineDecisionFunction{T}},
                                 ::Type{<:MOI.AbstractVectorSet}) where {T}
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:VectorAffineDecisionConstraintBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:VectorAffineDecisionConstraintBridge{T, S}}) where {T, S}
    return [(MOI.VectorAffineFunction{T}, S)]
end
function MOIB.Constraint.concrete_bridge_type(::Type{<:VectorAffineDecisionConstraintBridge{T}},
                                              ::Type{<:VectorAffineDecisionFunction{T}},
                                              S::Type{<:MOI.AbstractVectorSet}) where T
    return VectorAffineDecisionConstraintBridge{T, S}
end

MOI.get(b::VectorAffineDecisionConstraintBridge{T, S}, ::MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, S}) where {T, S} = 1
MOI.get(b::VectorAffineDecisionConstraintBridge{T, S}, ::MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, S}) where {T, S} = [b.constraint]

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::VectorAffineDecisionConstraintBridge{T}) where T
    f = bridge.decision_function
    # Remove mapped variables to properly unbridge
    from_decision(v) = begin
        result = any(t -> t.variable_index == v, f.decision_part.terms)
    end
    g = VectorAffineDecisionFunction(
        MOIU.filter_variables(v -> !from_decision(v), f.variable_part),
        copy(f.decision_part),
        copy(f.known_part))
    return g
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintSet,
                 bridge::VectorAffineDecisionConstraintBridge{T}) where T
    return bridge.set
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal,
                 bridge::VectorAffineDecisionConstraintBridge{T}) where T
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::VectorAffineDecisionConstraintBridge{T}) where T
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.delete(model::MOI.ModelLike, bridge::VectorAffineDecisionConstraintBridge)
    MOI.delete(model, bridge.constraint)
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintFunction,
                 bridge::VectorAffineDecisionConstraintBridge{T,S},
                 f::VectorAffineDecisionFunction{T}) where {T,S}
    # Update bridge functions and function constant
    bridge.decision_function = f
    # Recalculate total constants and modify constraint function
    constants = f.variable_part.constants +
        f.known_part.constants
    MOI.set(model, MOI.ConstraintFunction(), bridge.constraint,
            MOI.VectorAffineFunction(f.variable_part.terms, constants))
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::VectorAffineDecisionConstraintBridge{T,S}, change::S) where {T,S}
    f = bridge.decision_function
    bridge.set = change
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint, change)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::VectorAffineDecisionConstraintBridge{T,S}, change::MOI.VectorConstantChange) where {T,S}
    f = bridge.decision_function
    # Modify variable part of decision function
    f.variable_part.constants .= change.new_constant
    # Recalculate total constants and modify constraint function
    constants = f.variable_part.constants +
        f.known_part.constants
    MOI.set(model, MOI.ConstraintFunction(), bridge.constraint,
            MOI.VectorAffineFunction(f.variable_part.terms, constants))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::VectorAffineDecisionConstraintBridge{T,S}, change::MOI.MultirowChange) where {T,S}
    f = bridge.decision_function
    # Modify variable part of decision function
    modify_coefficients!(f.variable_part.terms, change.decision, change.new_coefficients)
    # Modify decision part if variable originated from a decision
    if !isempty(findall(t -> t.scalar_term.variable_index == change.variable, f.decision_part.terms))
        modify_coefficients!(f.decision_part.terms, change.variable, change.new_coefficients)
    end
    # Modify the variable part of the constraint as usual
    MOI.modify(model, bridge.constraint, change)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::VectorAffineDecisionConstraintBridge{T,S}, change::DecisionMultirowChange) where {T,S}
    f = bridge.decision_function
    # Update the decision coefficient
    modify_coefficients!(f.decision_part.terms, change.decision, change.new_coefficients)
    # Update mapped variable through MultiRowChange
    MOI.modify(model, bridge, MOI.MultiRowChange(change.decision, change.new_coefficients))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::VectorAffineDecisionConstraintBridge{T,S}, change::KnownCoefficientChange) where {T,S}
    f = bridge.decision_function
    # Query known value
    known_value = MOI.get(model, MOI.VariablePrimal(), change.known)
    seen = Int[]
    for i in findall(t -> t.scalar_term.variable_index == change.decision, f.known_part.terms)
        # Update known part of constraint constant
        coefficient = f.known_part.terms[i].scalar_term.coefficient
        f.known_part.constants[i] +=
            (change.new_coefficient - coefficient) * known_value
        push!(seen, i)
    end
    for i in filter(i -> i ∉ seen, collect(1:MOI.output_dimension(f)))
        f.known_part.constants[i] += change.new_coefficient * known_value
    end
    # Update the known decision coefficient
    modify_coefficients!(f.known_part.terms, change.known, change.new_coefficients)
    # Recalculate total constants and modify constraint function
    constants = f.variable_part.constants +
        f.known_part.constants
    MOI.set(model, MOI.ConstraintFunction(), bridge.constraint,
            MOI.VectorAffineFunction(f.variable_part.terms, constants))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::VectorAffineDecisionConstraintBridge{T,S}, change::KnownValueChange) where {T,S}
    f = bridge.decision_function
    seen = Int[]
    for i in findall(t -> t.scalar_term.variable_index == change.known, f.known_part.terms)
        # Update known part of constraint constant
        coefficient = f.decision_part.terms[i].scalar_term.coefficient
        f.known_part.constants[i] +=
            coefficient * change.value_difference
    end
    # Recalculate total constant and shift constraint set
    constant = f.variable_part.constant +
        f.known_part.constant
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint,
            MOIU.shift_constant(bridge.set, -constant))
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::VectorAffineDecisionConstraintBridge{T,S}, change::KnownValuesChange) where {T,S}
    f = bridge.decision_function
    known_vals = zeros(T, MOI.output_dimension(f))
    for term in f.known_part.terms
        i = term.output_index
        scalar_term = term.scalar_term
        known_vals[i] += scalar_term.coefficient * MOI.get(model, MOI.VariablePrimal(), scalar_term.variable_index)
    end
    # Update known part of objective
    f.known_part.constants .= known_vals
    # Recalculate total constants and modify constraint function
    constants = f.variable_part.constants +
        f.known_part.constants
    MOI.set(model, MOI.ConstraintFunction(), bridge.constraint,
            MOI.VectorAffineFunction(f.variable_part.terms, constants))
    return nothing
end


text:
#! /usr/bin/env julia

using BoundingBoxes

@boundingbox Bounds2, "x", "y"

a = Bounds2(Float64)

# warmup rand
rand()

# dummy

#update with random values
loop1() =  for i = 1:10000
    a.x_min = rand()
    a.y_min = rand()
    a.x_max = rand()
    a.y_max = rand()
end

@time loop1()
@time loop1()

#now make immutable
immutable Bounds2I{T}
    x_min::T
    y_min::T
    x_max::T
    y_max::T
end

b = Bounds2I{Float64}(Inf,Inf,-Inf,-Inf)
#update with random values
loop2() = for i = 1:10000
    b = Bounds2I{Float64}(rand(),rand(),rand(),rand())
end

@time loop2()
@time loop2()

max_stars_repo_name: sjkelly/JuliaPlayground

max_issues_repo_name: sjkelly/JuliaPlayground

max_forks_repo_name: sjkelly/JuliaPlayground


text:
module TuringDash

include("render_funcs.jl")
include("utils.jl")
using Dash, DashBootstrapComponents
using DashHtmlComponents, DashCoreComponents
using PlotlyJS
using Turing
using Random
using CSV, DataFrames
using Base64

#Initialize channel to get data during optimization
chn = Base.Channel{Vector{Float64}}(Inf)
x_graph = []
y_graph = []
y_graph1 = []
itr = 0
eval_f = nothing
range_value = []
df = nothing

function make_app()
    global chn, x_graph, y_graph, itr, eval_f
    global range_value

    chn = Base.Channel{Vector{Float64}}(Inf)
    x_graph = []
	y_graph = []
    y_graph1 = []
    itr = 0
    df = nothing
    nperiteration = 1

    external_stylesheets = "https://cdnjs.cloudflare.com/ajax/libs/font-awesome/4.7.0/css/font-awesome.min.css"
    app = dash(external_stylesheets=[dbc_themes.BOOTSTRAP,external_stylesheets],
	    suppress_callback_exceptions=true)

    app.layout = html_div([
        render_all()
    ])
    callback!(app,
        Output("dtu-params-plot-loc", "children"),
        Input("dtu-btn-turing", "n_clicks"),
        Input("dtu-upload-datan", "contents"),
        State("dtu-def-data-drop", "value"),
        State("dtu-user-funcs", "value"),
        State("dtu-user-func-init", "value"),
        State("dtu-upload-datan", "filename")
    ) do n_clicks, contents, drop_val, user_funcs, funcs_init, filename
        global chn
        global chn, x_graph, y_graph, y_graph1, itr, df
        if (contents isa Nothing) && (drop_val == "0")
            throw(PreventUpdate())
        else
            if contents isa Nothing
                df = load_default(drop_val)
            else
                df = parse_contents(contents, filename)
            end
        end
        ctx = callback_context()
        if isempty(ctx.triggered)
            throw(PreventUpdate())
        end
        if !(ctx.triggered[1].prop_id == "dtu-btn-turing.n_clicks")
            throw(PreventUpdate())
        end

        if (user_funcs == "") || (funcs_init == "")
            throw(PreventUpdate())
        end

        x_graph = []
        y_graph = []
        y_graph1 = []
        itr = 0

        for colname in Symbol.(names(df))
            @eval $colname = df.$colname
        end
        # evealuate model
        eval(Meta.parse(user_funcs))

        #evaluate eval_definition
        model = eval(Meta.parse(funcs_init))

        alg = MH(
            [0.25 0.05;
             0.05 0.50]
        )
        rng = Random.GLOBAL_RNG

        nsamples = 500
        @async begin
            r = sample(rng, model, alg, nperiteration; chain_type=MCMCChains.Chains, save_state=true, progress=false)
            put!(chn, Array(r)[end,1:2])
            for _ in (nperiteration + 1):nperiteration:nsamples
                r = Turing.Inference.resume(r, nperiteration, save_state=true, progress=false)
                put!(chn, Array(r)[end,1:2])
            end
        end
        # @show nsamples
        return render_liveupdates()
    end

    callback!(app,
        Output("dtu-params-state1", "figure"),
        Output("dtu-params-state2", "figure"),
        Output("interval-component", "disabled"),
        Input("interval-component", "n_intervals"),
        State("dtu-params-state1", "figure"),
        State("dtu-params-state2", "figure"),
        prevent_initial_call=true) do n, state_fig1, state_fig2
        st = false
        global x_graph, y_graph, y_graph1, chn, itr
        if isready(chn)
            val = take!(chn)
            itr += nperiteration
            append!(x_graph, itr)
            append!(y_graph, val[1])
			append!(y_graph1, val[2])
        else
            st = true
        end
        scatter_plt1 = scatter(
            ;x=x_graph, y=y_graph,
            mode="lines", line_color="blue"
        )
        scatter_plt2 = scatter(
            ;x=x_graph, y=y_graph1,
            mode="lines", line_color="blue"
        )
        fig1 = PlotlyJS.plot([scatter_plt1], Layout(
            title="Trace Plot: Var 1",
            xaxis_title="Iteration",
            yaxis_title="Value",
        ))
        fig2 = PlotlyJS.plot([scatter_plt2], Layout(
            title="Trace Plot: Var 2",
            xaxis_title="Iteration",
            yaxis_title="Value",
        ))
        return fig1, fig2, st
    end
    return app
end

end # module


text:
"""
All outputs must inherit from Output.

Simulation outputs are decoupled from simulation behaviour and in
many cases can be used interchangeably.
"""
abstract type Output{T} <: AbstractVector{T} end

"Generic ouput constructor. Converts init array to vector of frames."
(::Type{F})(init::AbstractMatrix; kwargs...) where F <: Output =
    F(; frames=[deepcopy(init)], kwargs...)
(::Type{F})(init::NamedTuple; kwargs...) where F <: Output =
    F(; frames=[deepcopy(init)], kwargs...)

(::Type{F})(o::T; kwargs...) where F <: Output where T <: Output =
    F(; frames=frames(o), starttime=starttime(o), stoptime=stoptime(o), kwargs...)

# Forward base methods to the frames array
Base.parent(o::Output) = frames(o)
Base.length(o::Output) = length(frames(o))
Base.size(o::Output) = size(frames(o))
Base.firstindex(o::Output) = firstindex(frames(o))
Base.lastindex(o::Output) = lastindex(frames(o))
Base.@propagate_inbounds Base.getindex(o::Output, i) =
    getindex(frames(o), i)
Base.@propagate_inbounds Base.setindex!(o::Output, x, i) = setindex!(frames(o), x, i)
Base.push!(o::Output, x) = push!(frames(o), x)
Base.append!(o::Output, x) = append!(frames(o), x)

"""
Mixin of basic fields for all outputs
"""
@premix @default_kw struct Output{T<:AbstractVector}
    frames::T      | []
    running::Bool  | false
    starttime::Any | 1
    stoptime::Any  | 1
end

# Getters and setters
frames(o::Output) = o.frames
starttime(o::Output) = o.starttime
stoptime(o::Output) = o.stoptime
tspan(o::Output) = (stoptime(o), starttime(o))
isrunning(o::Output) = o.running
setrunning!(o::Output, val) = o.running = val
setstarttime!(output, x) = output.starttime = x
setstoptime!(output, x) = output.stoptime = x

# Placeholder methods for graphic functions that are
# ignored in simple outputs
settimestamp!(o::Output, f) = nothing
fps(o::Output) = nothing
setfps!(o::Output, x) = nothing
showfps(o::Output) = nothing
isasync(o::Output) = false
isstored(o::Output) = true
isshowable(o::Output, f) = false
finalize!(o::Output, args...) = nothing
delay(o::Output, f) = nothing
showframe(o::Output, args...) = nothing

# Frame strorage and updating
frameindex(o::Output, data::AbstractSimData) = frameindex(o, currentframe(data))
frameindex(o::Output, f) = isstored(o) ? f : oneunit(f)

zeroframes(init::AbstractArray, nframes) = [zero(init) for f in 1:nframes]
zeroframes(init::NamedTuple, nframes) =
    [map(layer -> zero(layer), init) for f in 1:nframes]

@inline blockdo!(data::SimData, frame::AbstractArray, index, f) =
    return @inbounds frame[index...] = data[index...]

storeframe!(output, data) = storeframe!(output, data, frameindex(output, data))
storeframe!(output, data::SimData, f) = begin
    checkbounds(output, f)
    blockrun!(data, output[f], f)
end
storeframe!(output, multidata::MultiSimData, f) = begin
    checkbounds(output, f)
    for key in keys(multidata)
        # TODO use blocks for MutiSimData?
        # blockrun!(data(multidata)[key], output[f][key], f)
        for i in CartesianIndices(output[f][key])
            output[f][key][i] = data(multidata)[key][i]
        end
    end
end
# Replicated frames
storeframe!(output, data::AbstractVector{<:SimData}, f) = begin
    for j in 1:size(output[1], 2), i in 1:size(output[1], 1)
        replicatesum = zero(eltype(output[1]))
        for d in data
            replicatesum += d[i, j]
        end
        output[f][i, j] = replicatesum / length(data)
    end
end

allocateframes!(o::Output, init, framerange) = begin
    append!(frames(o), [similar(init) for f in framerange])
    o
end

"""
Frames are preallocated and reused.
"""
initframes!(o::Output, init) = begin
    first(o) .= init
    for f = (firstindex(o) + 1):lastindex(o)
        @inbounds o[f] .= zero(eltype(init))
    end
    o
end
initframes!(o::Output, init::NamedTuple) = begin
    for key in keys(init)
        first(o)[key] .= init[key]
    end
    for f = (firstindex(o) + 1):lastindex(o)
        for key in keys(init)
            @inbounds o[f][key] .= zero(eltype(init[key]))
        end
    end
    o
end