text:
module TestMultiReturningFunction6

using Modia
using Modia.StaticArrays
@usingModiaPlot

mutable struct Mbs
    phi1
    w1
    phi2
    w2
    derw1
    derw2

    Mbs() = new(0,0,0,0,0,0)
end

function setStates(mbs, phi1, w1, phi2, w2)
    mbs.phi1 = phi1
    mbs.w1   = w1
    mbs.phi2 = phi2
    mbs.w2   = w2
    return mbs
end

function setAccelerations(mbs::Mbs,derw1,derw2)
    mbs.derw1 = derw1
    mbs.derw2 = derw2
    return mbs
end

function getForces(mbs::Mbs)
    # m*L^2*derw + m*g*l*sin(phi) ) = u
    #   m=L=g=1
    # derw + sin(phi) = 0
    u1 = mbs.derw1 + sin(mbs.phi1)
    u2 = mbs.derw2 + sin(mbs.phi2)
    return SVector(u1,u2)
end

implicitDependency(x...) = nothing

Pendulum = Model(
    phi1 = Var(init=pi/2),
    w1   = Var(init=0.0),
    phi2 = Var(init=pi/4),
    w2   = Var(init=0.0),
    mbs  = Mbs(),
    equations = :[
        w1   = der(phi1)
        w2   = der(phi2)
        mbs1 = setStates(mbs,phi1,w1,phi2,w2)
        mbs2 = setAccelerations(mbs1,der(w1),der(w2))
        (0,0) = implicitDependency(getForces(mbs2), mbs1,der(w1),der(w2))
    ]
)

pendulum = @instantiateModel(Pendulum , unitless=true, log=false, logDetails=false, logCode=false, logStateSelection=false)

simulate!(pendulum, stopTime = 10.0, log=true,requiredFinalStates=[-1.0810242374611811, -0.9468374084033309, 0.13969589909383123, -0.7714994595239327])

plot(pendulum, ("phi1", "w1", "phi2", "w2"))

end


text:
using Match

function rw(ex)
  @match ex begin
    Expr(:comparison, [l,op,r], _) => "Comparison, lhs = $(l)"
    _ => "No match"
  end
end

rw( :(a == 1) )

ext: jl


text:
using ReinforcementLearning
using OpenSpiel

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:TabularCFR},
    ::Val{:OpenSpiel},
    game;
    n_iter = 300,
    seed = 123,
)
    env = OpenSpielEnv(game)
    rng = StableRNG(seed)
    π = TabularCFRPolicy(; rng = rng)

    description = "# Play `$game` in OpenSpiel with TabularCFRPolicy"
    Experiment(π, env, StopAfterStep(300, is_show_progress=!haskey(ENV, "CI")), EmptyHook(), description)
end

ex = E`JuliaRL_TabularCFR_OpenSpiel(kuhn_poker)`
run(ex)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
@testset "▶ types/drawing               " begin
    # SCATTER2D
    f = Figure()
    plot(rand(Float64, 3, 2))
    s = f.axes[1].drawings[1]
    @test s isa GPlot.Drawing2D
    @test s isa GPlot.Drawing
    @test s isa GPlot.Scatter2D
    @test isnothing(s.linestyles[1].lwidth)
    @test isnothing(s.markerstyles[1].color)
    @test isempty(s.labels)
    s.linestyles[1].lwidth = 0.5
    @test s.linestyles[1].lwidth == 0.5

    # FILL2D
    clf()
    fill_between!(1:5, rand(Float16, 5), 0.5)
    f = gcf().axes[1].drawings[1]
    @test f isa GPlot.Drawing2D
    @test f isa GPlot.Fill2D
    @test isnothing(f.xmin)
    @test isnothing(f.xmax)
    @test f.fillstyle.fill == colorant"cornflowerblue"

    # HIST2D
    clf()
    hist(rand(Int16, 5))
    h = gcf().axes[1].drawings[1]
    @test h isa GPlot.Drawing2D
    @test h isa GPlot.Hist2D
    @test isnothing(h.barstyle.color)
    @test h.horiz == false
    @test isnothing(h.bins)
    @test h.scaling == "none"

    # GROUPEDBAR2D
    f = Figure()
    bar(rand(Float32, 5, 3))
    g = f.axes[1].drawings[1]
    @test g isa GPlot.Drawing2D
    @test g isa GPlot.Bar2D
    @test g.stacked == false
    @test g.horiz == false
end

@testset "▶ /drawing                    " begin
    # DATA AGGREGATORS
    x = [5.0, missing, missing, 3.0, 2.0]
    z,h,n = G.plotdata(x)
    @test h == true
    @test n == 1
    @test checkzip(z, hcat(1:length(x), x))

    y = hcat([2.0, missing, 3.0, Inf, 2.0], [1.0, 2.0, NaN, 0, 1.0])
    z,h,n = G.plotdata(x, y)
    @test h == true
    @test n == 2
    @test checkzip(z, hcat(x, y))

    x = 1:1:10
    y1 = @. exp(x)
    y2 = @. sin(x)
    d = G.filldata(x,y1,y2)
    @test checkzip(d, hcat(x,y1,y2))

    x = [1, 2, missing, 5]
    d,h,n,r = G.histdata(x)
    @test checkzip(d, x)
    @test h == true
    @test n == 3
    @test r == (1.0, 5.0)

    # SCATTER2D
    # -- basic plot
    x = randn(5)
    y = randn(5)
    plot(x, y)
    el = gca().drawings[1]
    @test el isa GPlot.Scatter2D
    @test checkzip(el.data, hcat(x, y))
    plot!(2x, 2y)
    el2 = gca().drawings[2]
    plot!(1:3, exp.(1:3), lw=2, ls="--") # "--" is 9 in GLE
    el3 = gca().drawings[3]
    @test checkzip(el3.data, hcat(1:3, exp.(1:3)))
    @test el3.linestyles[1].lstyle == 9
    @test el3.linestyles[1].lwidth == 2.0
    plot!(exp.(1:3), marker="o")
    el4 = gca().drawings[4]
    @test isnothing(el4.linestyles[1].lstyle)
    @test el4.markerstyles[1].marker == "circle"
    @test isnothing(el.linestyles[1].smooth) # only 5 points

    plot(1:5, hcat(exp.(1:5), sin.(1:5)))
    el2 = gca().drawings[1]
    @test checkzip(el2.data, hcat(1:5, exp.(1:5), sin.(1:5)))
    xy = rand(Float32, 5, 3)
    plot(xy)
    @test checkzip(gca().drawings[1].data, hcat(1:size(xy, 1), xy))

    # -- "implicit"
    cla()
    plot(sin, 0, 1)
    x = range(0, stop=1, length=100)
    y = @. sin(x)
    @test checkzip(gca().drawings[1].data, hcat(x, y))

    # -- line
    cla()
    line([0,0],[1,1])
    @test checkzip(gca().drawings[1].data, [0 0; 1 1])

    erase!(gcf())
    # -- multiplot
    plot(x, y, 2y, 3y)
    el = gca().drawings[1]
    @test checkzip(el.data, hcat(x, y, 2y, 3y))

    erase!(gcf())
    # -- scatterplot
    scatter(x, y)
    scatter!(2x, 2y, 3y)
    el1 = gca().drawings[1]
    el2 = gca().drawings[2]
    @test el1.linestyles[1].lstyle == -1
    @test el1.markerstyles[1].marker == "fcircle"
    @test el2.linestyles[1].lstyle == -1
    @test el2.linestyles[2].lstyle == -1
    @test el2.markerstyles[1].marker == "fcircle"
    @test el2.markerstyles[2].marker == "fcircle"

    # FILL2D
    erase!(gcf())
    gcf().transparency = true
    x = 1:5
    y1 = exp.(x)
    y2 = sin.(x)
    fill_between(x, y1, y2, alpha=0.5)
    fill_between!(x, 0, y2, color="red", alpha=0.5)
    fill_between!(x, y1, 0)
    el1 = gca().drawings[1]
    el2 = gca().drawings[2]
    el3 = gca().drawings[3]
    @test checkzip(el1.data, hcat(x, y1, y2))
    @test el2.fillstyle.fill == RGBA{Colors.N0f8}(1.0,0.0,0.0,0.502)
    @test checkzip(el2.data, hcat(x, 0*y2, y2))
    @test checkzip(el3.data, hcat(x, y1, 0*y1))

    cla()
    fill_between(x, 1, 2)
    el = gca().drawings[1]
    @test checkzip(el.data, hcat(x, zero(x).+1, zero(x).+2))

    # HIST2D
    erase!(gcf())
    x = rand(Float32, 10)
    hist(x)
    @test checkzip(gca().drawings[1].data, x)

    # BAR
    cla()
    bar(x)
    y = rand(Float32, 10, 3)
    @test checkzip(gca().drawings[1].data, hcat(1:length(x), x))
    bar(x, x, x)
    @test checkzip(gca().drawings[1].data, hcat(x, x, x))
    bar(x, y)
    @test checkzip(gca().drawings[1].data, hcat(x, y))
    bar!(x, x, x)
    @test checkzip(gca().drawings[2].data, hcat(x, x, x))
end

@testset "▶ set_prop/drawing            " begin
    f = Figure()
    x, y = 1:2, exp.(1:2)
    plot(x, y, label="blah")
    @test gca().drawings[1].labels == ["blah"]
    legend(pos="top-left")
    @test gca().legend.position == "tl"
    @test_throws GPlot.OptionValueError legend(pos="blah")

    hist(y, bins=2, norm="pdf")
    @test gca().drawings[1].bins == 2
    @test gca().drawings[1].scaling == "pdf"
    @test_throws GPlot.OptionValueError hist(y, norm="blah")
    hist(y, horiz=true)
    @test gca().drawings[1].horiz

    x1 = 1:5
    y1 = exp.(x1)
    y2 = sin.(x1)
    fill_between(x1, y1, y2, min=0, max=6)
    @test gca().drawings[1].xmin == 0.0
    @test gca().drawings[1].xmax == 6.0

    bar(x1, y1, y1+y2, stacked=true, labels=["blah", "blih"])
    @test gca().drawings[1].stacked
    @test gca().drawings[1].labels == ["blah", "blih"]

    cla()
    fill_between(1:1:5, 1, 2, label="fill")
    @test gca().drawings[1].label == "fill"
end

@testset "▶ apply_gle/drawing           " begin
    g = G.GLE()
    f = G.Figure("blah", reset=true);
    G.add_axes2d!()

    # Scatter2D (see also apply_linestyle)
    plot([1, 2], exp.([1, 2]), lwidth=2.0, label="line")
    legend()
    G.apply_axes!(g, f.axes[1], f.id, 1); s = String(take!(g))
    isin(s, "data")
    isin(s, "d1=c1,c2")
    isin(s, "d1 line lwidth 2.0")
    isin(s, "begin key")
    isin(s, "compact")
    isin(s, "text \"line\" line lwidth 2.0")
    isin(s, "end key")

    clf()
    scatter([1, 2], exp.([1, 2]), color="blue")
    G.apply_axes!(g, f.axes[1], f.id, 1); s = String(take!(g))
    isin(s, "d1=c1,c2")
    isin(s, "d1 marker fcircle color rgba(0.0,0.0,1.0,1.0)")

    # FILL2D
    clf()
    fill_between([1, 2], 1, 2)
    G.apply_axes!(g, f.axes[1], f.id, 1); s = String(take!(g))
    isin(s, "d1=c1,c2 d2=c1,c3")
    isin(s, "fill d1,d2 color rgba(0.392,0.584,0.929,1.0)")

    # HIST2D
    clf()
    x = [0.1, 0.1, 0.2, 0.1, 0.3, 0.5, 0.7, 0.2, 0.1]
    hist(x, norm="pdf", col="blue", ecol="red", nbins=9)
    s = G.assemble_figure(gcf(), debug=true)
    isin(s, "let d2 = hist d1 from $(minimum(x)) to $(maximum(x)) bins 9")
    isin(s, "let d2 = d2*1.666")
    isin(s, "bar d2 width 0.066")
    isin(s, "color rgba(1.0,0.0,0.0,1.0) fill rgba(0.0,0.0,1.0,1.0)")

    clf()
    # color is not defined
    hist(x)
    s = G.assemble_figure(gcf(), debug=true)
    @test gca().drawings[1].barstyle.color == G.GP_ENV["PALETTE"][1]
    clf()
    hist(x, fill="red")
    s = G.assemble_figure(gcf(), debug=true)
    @test gca().drawings[1].barstyle.color == c"white"

    # BAR2D
    clf()
    x  = [1, 2, 3, 4, 5]
    y  = [3, 3, 4, 4, 3]
    y2 = [5, 5, 6, 6, 5]
    bar(x, y, y2, stacked=true)
    s = G.assemble_figure(gcf(), debug=true)
    isin(s, "bar d1")
    isin(s, "bar d2 from d1")

    # non-stacked
    clf()
    x  = [1, 2, 3, 4, 5]
    y  = [3, 3, 4, 4, 3]
    y2 = [5, 5, 6, 6, 5]
    bar(x, y, y2, stacked=false)
    s = G.assemble_figure(gcf(), debug=true)
    isin(s, "bar d1,d2 color $(G.col2str(G.GP_ENV["PALETTE"][1])),$(G.col2str(G.GP_ENV["PALETTE"][2]))")
    isin(s, "fill rgba(1.0,1.0,1.0,1.0),rgba(1.0,1.0,1.0,1.0)")
end


text:
#SignalProcessing: Discrete-time-domain dataset generators
#-------------------------------------------------------------------------------

#Creates step function of a single-pole system, using time from ref.x:
#-------------------------------------------------------------------------------
function Base.step(ref::DataF1, p::Pole; tdel=0, amp=1)
	#H(s) = 1/(s*tau+1)
	npoints = length(ref) #WANTCONST
	e = Base.MathConstants.e
	x = ref.x
	one_tau = value(:rad, p) #1/tau
	y = zero(ref.y)

	#Find first point where x>=tdel:
	istart = npoints+1
	for i in 1:npoints
		if x[i] >= tdel; istart=i; break; end
	end
	for i in istart:npoints
		y[i]=amp*(1-e^(-(x[i]-tdel)*one_tau))
	end
	return DataF1(x, y)
end

#Creates pulse response of a single-pole system, using time from ref.x:
# tpw: pulse width
#-------------------------------------------------------------------------------
function pulse(ref::DataF1, p::Pole; tdel=0, amp=1, tpw::Number=0)
	ensure(tpw > 0, ArgumentError("tpw must be positive"))
	return (step(ref, p, tdel=tdel, amp=amp) -
		step(ref, p, tdel=tdel+tpw, amp=amp))
end

#Creates a signal pattern from a sequence of data points, and a pulse response:
# seq: a binary data sequence
# Π: Pulse response
# tbit: bit period
function pattern(seq::Vector{T}, Π::DataF1; tbit::Number=0) where {T<:Number}
	ensure(tbit > 0, ArgumentError("tbit must be positive"))
	result = 0
	for i in 1:length(seq)
		result += seq[i]*xshift(Π, (i-1)*tbit)
	end
	return result
end
pattern(seq::DataF1, Π::DataF1; tbit::Number=0) = pattern(seq.y, Π, tbit = tbit)

#Broadcast 1-arguent function with reference DataF1 @ 2nd arg.
pattern(d1, d2::DataMD, args...; kwargs...) =
	broadcastMD(CastType(DataF1,2), pattern, d1, d2, args...; kwargs...)

#==HACK: Find a way to register these functions differently.
Options:
	Move functions to MDDatsets
	Create broadcasting macro
	...
==#

#Last line


text:
using Test
using BCVUMPS
using BCVUMPS:MCMC,magofβ,eneofβ
using Random

@testset "MCMC" begin
    Random.seed!(100)
    for β in [0.6,0.8]
        mag22,ene22 = MCMC(Ising22(1),12,β,10000,100000)
        mag33,ene33 = MCMC(Ising33(1),12,β,10000,100000)
        @test isapprox(mag22, mag33, atol=1e-2)
        @test isapprox(ene22, ene33, atol=1e-2)
        @test isapprox(mag22, magofβ(Ising(),β), atol=1e-2)
        @test isapprox(ene22, eneofβ(Ising(),β), atol=1e-2)
        @test isapprox(mag33, magofβ(Ising(),β), atol=1e-2)
        @test isapprox(ene33, eneofβ(Ising(),β), atol=1e-2)
    end
end


text:
import Pkg;
Pkg.activate(joinpath(@__DIR__, "../../"))  # environment: CytofResearch

# Load libraries
println("Loading libraries ...")
import CSV
using DataFrames
include("../cb/PreProcess.jl")
println("Finished loading libraries.")

# Path to data.
data_dir = "../data"
path_to_data = "$(data_dir)/cb_transformed.csv"

function separatesamples(Y)
  sample_ids = sort(unique(Y.sample_id))
  return Matrix.(select(Y[Y.sample_id .== i, :], Not(:sample_id))
                 for i in sample_ids)
end

# Read transformed CB data into Y as a data frame.
# Replace `missing` with `NaN`.
function final_preprocess(Y; p=0.9, digits=3)
  markernames = filter(x -> x != :sample_id, names(Y))
  y = separatesamples(Y)
  isgoodmarker = PreProcess.preprocess!(y, maxNanOrNegProp=p,
                                        maxPosProp=p, rowThresh=-6.0)

  includedmarkers = Symbol[]
  excludedmarkers = Symbol[]

  for j in 1:length(isgoodmarker)
    marker = markernames[j]
    if isgoodmarker[j]
      append!(includedmarkers, [marker])
    else
      append!(excludedmarkers, [marker])
    end
  end

  # Concatenate y and convert to data frame
  Y_reduced = DataFrame(vcat(y...), includedmarkers)

  # Convert NaN to missing
  Y_reduced = ifelse.(isnan.(Y_reduced), missing, Y_reduced)

  # Append a column for sample_id
  Y_reduced[!, :sample_id] = vcat([fill(i, size(y[i], 1))
                                   for i in 1:length(y)]...)

  return round.(Y_reduced, digits=digits), excludedmarkers
end

# Read data
Y = coalesce.(CSV.read(path_to_data), NaN)

# Excluded markers path.
path_to_excluded_markers = "excluded_markers.txt"
rm(path_to_excluded_markers, force=true)

for p in (0.85, 0.9, 0.95, 0.99)
  Y_reduced, excludedmarkers = final_preprocess(Y, p=p)

  open(path_to_excluded_markers, "a") do io
    println(io, "Excluded markers for p=$(p): $(String.(excludedmarkers))")
  end

  # Path to reduced data
  path_to_reduced_data = "$(data_dir)/cb_transformed_reduced_p$(p).csv"

  # Write reduced data to disk.
  println("Writing reduced-and-transformed CB data to: ")
  println(path_to_reduced_data)

  CSV.write(path_to_reduced_data, Y_reduced)
end
println("Done!")

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
using DataFrames
using Plots
using CSV

# First Julia DataScience program for Digimons

file = string("~/github/exercism/julia/datasets/digimon/DigiDB_digimonlist.csv")

df = CSV.read(file)

max_stars_repo_name: Devlious/LearningJulia

max_issues_repo_name: Devlious/LearningJulia

max_forks_repo_name: Devlious/LearningJulia


text:
@printf("\nTesting  calculus non linear operators\n")
##testing Scale
m = 3
x = randn(m)
r = randn(m)
A = Sigmoid(Float64,(m,),2)
op = 30*A

y, grad = test_NLop(op,x,r,verb)

Y = 30*(A*x)
@test norm(Y-y) <1e-8

m = 3
x = randn(m)
r = randn(m)
A = Pow(Float64,(m,),2)
op = -A

y, grad = test_NLop(op,x,r,verb)

Y = -A*x
@test norm(Y-y) <1e-8

#testing DCAT
n,m = 4,3
x = ArrayPartition(randn(n),randn(m))
r = ArrayPartition(randn(n),randn(m))
A = randn(n,n)
B = Sigmoid(Float64,(m,),2)
op = DCAT(MatrixOp(A),B)

y, grad = test_NLop(op,x,r,verb)

Y = ArrayPartition(A*x.x[1],B*x.x[2])
@test norm(Y-y) <1e-8

#testing HCAT
n,m = 4,3
x = ArrayPartition(randn(n),randn(m))
r = randn(m)
A = randn(m,n)
B = Sigmoid(Float64,(m,),2)
op = HCAT(MatrixOp(A),B)

y, grad = test_NLop(op,x,r,verb)

Y = A*x.x[1]+B*x.x[2]
@test norm(Y-y) <1e-8

m,n = 3,5
x = ArrayPartition(randn(m),randn(n))
r = randn(m)
A= Sin(Float64,(m,))
M = randn(m,n)
B = MatrixOp(M)
op = HCAT(A,B)

y, grad = test_NLop(op,x,r,verb)

Y = A*x.x[1]+M*x.x[2]
@test norm(Y-y) <1e-8

p = [2,1]
opP = AbstractOperators.permute(op,p)
xp = ArrayPartition(x.x[p]...)
J = Jacobian(opP,xp)'
println(size(J,1))
y, grad = test_NLop(opP,xp,r,verb)

#testing VCAT
n,m = 4,3
x = randn(m)
r = ArrayPartition(randn(n),randn(m))
A = randn(n,m)
B = Sigmoid(Float64,(m,),2)
op = VCAT(MatrixOp(A),B)

y, grad = test_NLop(op,x,r,verb)

Y = ArrayPartition(A*x,B*x)
@test norm(Y-y) <1e-8

#testing HCAT of VCAT
n,m1,m2,m3 = 4,3,2,7
x1 = randn(m1)
x2 = randn(m2)
x3 = randn(m3)
x = ArrayPartition(x1,x2,x3)
r = ArrayPartition(randn(n),randn(m1))
A1 = randn(n,m1)
A2 = randn(n,m2)
A3 = randn(n,m3)
B1 = Sigmoid(Float64,(m1,),2)
B2 = randn(m1,m2)
B3 = randn(m1,m3)
op1 = VCAT(MatrixOp(A1),B1)
op2 = VCAT(MatrixOp(A2),MatrixOp(B2))
op3 = VCAT(MatrixOp(A3),MatrixOp(B3))
op = HCAT(op1,op2,op3)

y, grad = test_NLop(op,x,r,verb)

Y = ArrayPartition(A1*x1+A2*x2+A3*x3,B1*x1+B2*x2+B3*x3)
@test norm(Y-y) <1e-8

#testing VCAT of HCAT
m1,m2,m3,n1,n2 = 3,4,5,6,7
x1 = randn(m1)
x2 = randn(n1)
x3 = randn(m3)
x = ArrayPartition(x1,x2,x3)
r = ArrayPartition(randn(n1),randn(n2))
A1 = randn(n1,m1)
B1 = Sigmoid(Float64,(n1,),2)
C1 = randn(n1,m3)
A2 = randn(n2,m1)
B2 = randn(n2,n1)
C2 = randn(n2,m3)
x = ArrayPartition(x1,x2,x3)
op1 = HCAT(MatrixOp(A1),         B1 ,MatrixOp(C1))
op2 = HCAT(MatrixOp(A2),MatrixOp(B2),MatrixOp(C2))
op = VCAT(op1,op2)

y, grad = test_NLop(op,x,r,verb)

Y = ArrayPartition(A1*x1+B1*x2+C1*x3,A2*x1+B2*x2+C2*x3)
@test norm(Y-y) <1e-8

#testing Compose

l,n,m = 5,4,3
x = randn(m)
r = randn(l)
A = randn(l,n)
C = randn(n,m)
opA = MatrixOp(A)
opB = Sigmoid(Float64,(n,),2)
opC = MatrixOp(C)
op = Compose(opA,Compose(opB,opC))

y, grad = test_NLop(op,x,r,verb)

Y = A*(opB*(opC*x))
@test norm(Y-y) <1e-8

## NN
m,n,l = 4,7,5
b = randn(l)
opS1 = Sigmoid(Float64,(n,),2)
x = ArrayPartition(randn(n,l),randn(n))
r = randn(n)

A1 = HCAT(LMatrixOp(b,n) ,Eye(n))
op = Compose(opS1,A1)
y, grad = test_NLop(op,x,r,verb)

###testing Reshape
n = 4
x = randn(n)
r = randn(n)
opS = Sigmoid(Float64,(n,),2)
op = Reshape(opS,2,2)

y, grad = test_NLop(op,x,r,verb)

Y = reshape(opS*x,2,2)
@test norm(Y-y) <1e-8

###testing BroadCast
n,l = 4,7
x = randn(n)
r = randn(n,l)
opS = Sigmoid(Float64,(n,),2)
op = BroadCast(opS,(n,l))

y, grad = test_NLop(op,x,r,verb)

Y = (opS*x).*ones(n,l)
@test norm(Y-y) <1e-8

n,l = 1,7
x = randn(n)
r = randn(n,l)
opS = Sigmoid(Float64,(n,),2)
op = BroadCast(opS,(n,l))

y, grad = test_NLop(op,x,r,verb)

Y = (opS*x).*ones(n,l)
@test norm(Y-y) <1e-8

##testing Sum
m = 5
x = randn(m)
r = randn(m)
A = randn(m,m)
opA = MatrixOp(A)
opB = Sigmoid(Float64,(m,),2)
op = Sum(opA,opB)

y, grad = test_NLop(op,x,r,verb)

Y = A*x+opB*x
@test norm(Y-y) <1e-8

## AffineAdd and NonLinearOperator
n = 10
d = randn(n)
T = AffineAdd(Exp(n),d,false)

r = randn(n)
x = randn(size(T,2))
y, grad = test_NLop(T,x,r,verb)
@test norm( y - (exp.(x)-d) ) < 1e-8

## AffineAdd and Compose NonLinearOperator
n = 10
d1 = randn(n)
d2 = randn(n)
T = Compose(AffineAdd(Sin(n),d2),AffineAdd(Eye(n),d1))

r = randn(n)
x = randn(size(T,2))
y, grad = test_NLop(T,x,r,verb)
@test norm( y - (sin.(x+d1)+d2) ) < 1e-8

n = 10
d1 = randn(n)
d2 = randn(n)
d3 = pi
T = Compose( AffineAdd(Sin(n),d3), Compose( AffineAdd(Exp(n),d2,false),AffineAdd(Eye(n),d1 ) ) )

r = randn(n)
x = randn(size(T,2))
y, grad = test_NLop(T,x,r,verb)
@test norm( y -( sin.(exp.(x+d1)-d2).+d3 )) < 1e-8

#### Axt_mul_Bx
n = 10
A,B = Eye(n),Sin(n)
P = Axt_mul_Bx(A,B)

x = randn(n)
r = randn(1)
y, grad = test_NLop(P,x,r,verb)
@test norm([(A*x)'*(B*x)]-y) < 1e-8

n,m = 3,4
A,B = MatrixOp(randn(n,m)), MatrixOp(randn(n,m))
P = Axt_mul_Bx(A,B)

x = randn(m)
r = randn(1)
y, grad = test_NLop(P,x,r,verb)
@test norm([(A*x)'*(B*x)]-y) < 1e-8

n,m,l = 3,7,5
A,B = MatrixOp(randn(n,m),l), MatrixOp(randn(n,m),l)
P = Axt_mul_Bx(A,B)
x = randn(m,l)
r = randn(l,l)
y, grad = test_NLop(P,x,r,verb)
@test norm((A*x)'*(B*x)-y) < 1e-8

n,m = 3,7
A,B = Sin(n,m), Cos(n,m)
P = Axt_mul_Bx(A,B)
x = randn(n,m)
r = randn(m,m)
y, grad = test_NLop(P,x,r,verb)
@test norm((A*x)'*(B*x)-y) < 1e-8

# testing with HCAT
m,n = 3,5
x = ArrayPartition(randn(m),randn(n))
r = randn(1)
b = randn(m)
A = AffineAdd(Sin(Float64,(m,)),b)
B = MatrixOp(randn(m,n))
op1 = HCAT(A,B)
C = Cos(Float64,(m,))
D = MatrixOp(randn(m,n))
op2 = HCAT(C,D)
P = Axt_mul_Bx(op1,op2)
y, grad = test_NLop(P,x,r,verb)
@test norm([(op1*x)'*(op2*x)]-y) < 1e-8

#test remove_displacement
y2, grad = test_NLop(remove_displacement(P),x,r,verb)
@test norm([(op1*x-b)'*(op2*x)]-y2) < 1e-8

# test permute
p = [2,1]
Pp = AbstractOperators.permute(P,p)
xp = ArrayPartition(x.x[p])
y2, grad = test_NLop(Pp,xp,r,verb)
@test norm(y2-y) < 1e-8

@test_throws Exception Axt_mul_Bx(Eye(2,2), Eye(2,1))
@test_throws Exception Axt_mul_Bx(Eye(2,2,2), Eye(2,2,2))

## Ax_mul_Bxt
n = 10
A,B = Eye(n),Sin(n)
P = Ax_mul_Bxt(A,B)
x = randn(n)
r = randn(n,n)
y, grad = test_NLop(P,x,r,verb)
@test norm((A*x)*(B*x)'-y) < 1e-9

n,m = 3,4
A,B = MatrixOp(randn(n,m)), MatrixOp(randn(n,m))
P = Ax_mul_Bxt(A,B)
x = randn(m)
r = randn(n,n)
y, grad = test_NLop(P,x,r,verb)
@test norm((A*x)*(B*x)'-y) < 1e-8

n,m,l = 3,7,5
A,B = MatrixOp(randn(n,m),l), MatrixOp(randn(n,m),l)
P = Ax_mul_Bxt(A,B)
x = randn(m,l)
r = randn(n,n)
y, grad = test_NLop(P,x,r,verb)
@test norm((A*x)*(B*x)'-y) < 1e-8

n,m = 3,7
A,B = Sin(n,m), Cos(n,m)
P = Ax_mul_Bxt(A,B)
x = randn(n,m)
r = randn(n,n)
y, grad = test_NLop(P,x,r,verb)
@test norm((A*x)*(B*x)'-y) < 1e-8

# testing with HCAT
m,n = 3,5
x = ArrayPartition(randn(m),randn(n))
r = randn(m,m)
b = randn(m)
A = AffineAdd(Sin(Float64,(m,)),b)
B = MatrixOp(randn(m,n))
op1 = HCAT(A,B)
C = Cos(Float64,(m,))
D = MatrixOp(randn(m,n))
op2 = HCAT(C,D)
P = Ax_mul_Bxt(op1,op2)
y, grad = test_NLop(P,x,r,verb)
@test norm((op1*x)*(op2*x)'-y) < 1e-8

#test remove_displacement
y2, grad = test_NLop(remove_displacement(P),x,r,verb)
@test norm((op1*x-b)*(op2*x)'-y2) < 1e-8

# test permute
p = [2,1]
Pp = AbstractOperators.permute(P,p)
xp = ArrayPartition(x.x[p])
y2, grad = test_NLop(Pp,xp,r,verb)
@test norm(y2-y) < 1e-8

@test_throws Exception Ax_mul_Bxt(Eye(2,2), Eye(2,1))
@test_throws Exception Ax_mul_Bxt(Eye(2,2,2), Eye(2,2,2))

## Ax_mul_Bx

n = 3
A,B = Eye(n,n), Eye(n,n)
P = Ax_mul_Bx(A,B)
x = randn(n,n)
r = randn(n,n)
y, grad = test_NLop(P,x,r,verb)
@test norm(x*x-y) < 1e-9

n = 3
A,B = Sin(n,n), Cos(n,n)
P = Ax_mul_Bx(A,B)
x = randn(n,n)
r = randn(n,n)
y, grad = test_NLop(P,x,r,verb)
@test norm((A*x)*(B*x)-y) < 1e-9

n = 3
A,B,C = Sin(n,n), Cos(n,n), Atan(n,n)
P = Ax_mul_Bx(A,B)
P2 = Ax_mul_Bx(C,P)
x = randn(n,n)
r = randn(n,n)
y, grad = test_NLop(P2,x,r,verb)
@test norm((C*x)*(A*x)*(B*x)-y) < 1e-9

n,l = 2,3
A,B = MatrixOp(randn(l,n),l), MatrixOp(randn(l,n),l)
P = Ax_mul_Bx(A,B)
x = randn(n,l)
r = randn(l,l)
y, grad = test_NLop(P,x,r,verb)
@test norm((A*x)*(B*x)-y) < 1e-8

@test_throws Exception Ax_mul_Bx(Eye(2), Eye(2))
@test_throws Exception Ax_mul_Bx(Eye(2,2), Eye(2,1))
@test_throws Exception Ax_mul_Bx(Eye(2,2,2), Eye(2,2,2))

# testing with HCAT
m,n = 3,5
x = ArrayPartition(randn(n,n),randn(m,n))
r = randn(n,n)
b = randn(n,n)
A = AffineAdd(Sin(Float64,(n,n)),b)
B = MatrixOp(randn(n,m),n)
op1 = HCAT(A,B)
C = Sin(Float64,(n,n))
D = MatrixOp(randn(n,m),n)
op2 = HCAT(C,D)
P = Ax_mul_Bx(op1,op2)
y, grad = test_NLop(P,x,r,verb)
@test norm((op1*x)*(op2*x)-y) < 1e-8

#test remove_displacement
y2, grad = test_NLop(remove_displacement(P),x,r,verb)
@test norm((op1*x-b)*(op2*x)-y2) < 1e-8

# test permute
p = [2,1]
Pp = AbstractOperators.permute(P,p)
xp = ArrayPartition(x.x[p])
y2, grad = test_NLop(Pp,xp,r,verb)
@test norm(y2-y) < 1e-8

#### some combos of Ax_mul_Bx etc...
n,m,l = 3,7,5
A,B = MatrixOp(randn(n,m),l), MatrixOp(randn(n,m),l)
P = Ax_mul_Bxt(A,B)
P2 = Axt_mul_Bx(A,P)
x = randn(m,l)
r = randn(l,n)
y, grad = test_NLop(P2,x,r,verb)
@test norm((A*x)'*((A*x)*(B*x)')-y) < 1e-8

n,m,l,k = 3,7,5,9
A,B = MatrixOp(randn(n,m),l), MatrixOp(randn(n,m),l)
C = MatrixOp(randn(k,m),l)
P = Axt_mul_Bx(A,B)
P2 = Ax_mul_Bx(C,P)
x = randn(m,l)
r = randn(k,l)
y, grad = test_NLop(P2,x,r,verb)
@test norm((C*x)*((A*x)'*(B*x))-y) < 1e-8

n,m,l,k = 3,7,5,9
A,B = MatrixOp(randn(n,m),l), MatrixOp(randn(n,m),l)
C = MatrixOp(randn(k,m),l)
P = Axt_mul_Bx(A,B)
P2 = Ax_mul_Bxt(C,P)
x = randn(m,l)
r = randn(k,l)
y, grad = test_NLop(P2,x,r,verb)
@test norm((C*x)*((A*x)'*(B*x))'-y) < 1e-8

#### HadamardProd

n = 3
A,B = Eye(n,n), Eye(n,n)
P = HadamardProd(A,B)
x = randn(n,n)
r = randn(n,n)
y, grad = test_NLop(P,x,r,verb)
@test norm(x.*x-y) < 1e-9

n,l = 3,2
A,B = Sin(n,l), Cos(n,l)
P = HadamardProd(A,B)
x = randn(n,l)
r = randn(n,l)
y, grad = test_NLop(P,x,r,verb)
@test norm((A*x).*(B*x)-y) < 1e-9

# testing with HCAT
m,n = 3,5
x = ArrayPartition(randn(m),randn(n))
r = randn(m)
b = randn(m)
A = AffineAdd(Sin(Float64,(m,)),b)
B = MatrixOp(randn(m,n))
op1 = HCAT(A,B)
C = Cos(Float64,(m,))
D = MatrixOp(randn(m,n))
op2 = HCAT(C,D)
P = HadamardProd(op1,op2)
y, grad = test_NLop(P,x,r,verb)
@test norm((op1*x).*(op2*x)-y) < 1e-9

#test remove_displacement
y2, grad = test_NLop(remove_displacement(P),x,r,verb)
@test norm((op1*x-b).*(op2*x)-y2) < 1e-8

# test permute
p = [2,1]
Pp = AbstractOperators.permute(P,p)
xp = ArrayPartition(x.x[p])
y2, grad = test_NLop(Pp,xp,r,verb)
@test norm(y2-y) < 1e-8

@test_throws Exception HadamardProd(Eye(2,2,2), Eye(1,2,2))


text:
"""
Check if a solution x is in the core of G. Namely,

    - ∑ x < v(N)
    - ∑ₛ x < v(S) ∀ S

"""
function isincore(x::Vector{Float64}, G::Game)
    if sum(x) != G.v(G.N) return false end

    for S in subsets(G.N)
        if sum(x[S]) < G.v(S) return false end
    end

    return true

end


text:
@testset "utility tests: get_rounded_zeros_and_ones!" begin
    A1 = [0 -0.99999; 0.999999 2.000000001; 1.0000001 -1.000001; -0.000001 0.0000001]
    A2 = [0 -0.99999; 0.999999 2.000000001; 1.0000001 -1.000001; -0.000001 0.0000001]

    LOpt.get_rounded_zeros_and_ones!(A1, 1E-5)
    @test A1[1,2] == -1
    @test A1[2,1] ==  1
    @test A1[2,2] >  2
    @test A1[3,1] == 1
    @test A1[3,2] == -1
    @test A1[4,1] == 0
    @test A1[4,2] == 0

    LOpt.get_rounded_zeros_and_ones!(A2, 1E-6)
    @test A2[1,2] > -1

end

@testset "utility tests: optimal_graph_edges" begin

    z_val = [0.0 1 1; 1 0 0; 1 0 0]
    edges = LOpt.optimal_graph_edges(z_val)
    @test edges[1] == (1,2)
    @test edges[2] == (1,3)

end

@testset "utility tests: algebraic_connectivity and fiedler_vector" begin

    adj_mat_1 = [0 1 2; 1 0 3; 2 3 0.0]
    graph_data = LOpt.GraphData(adj_mat_1)
    @test isapprox(graph_data.ac, 4.267949192, atol=1E-6)
    v = [0.7886751345948128, -0.5773502691896261, -0.21132486540518727]
    @test isapprox(v, graph_data.fiedler, atol=1E-6)

    adj_mat_2 = [0 1 2; 1 0 0; 2 0 0.0]
    graph_data = LOpt.GraphData(adj_mat_2)
    @test isapprox(graph_data.ac, 1.26794919243, atol=1E-6)
    v = [ -0.2113248654051879, 0.7886751345948118, -0.5773502691896268]
    @test isapprox(v, graph_data.fiedler, atol=1E-6)

    # Disconnected graph
    adj_mat_3 = [0 1 1 0; 1 0 0 0; 1 0 0 0; 0 0 0 0.0]
    graph_data = LOpt.GraphData(adj_mat_3)
    @test isapprox(graph_data.ac, 0, atol=1E-6)
    v = [ 0.5773502691896258, 0.5773502691896273, 0.5773502691896241, 0.0]
    @test isapprox(v, graph_data.fiedler, atol=1E-6)

end


text:
#TOOLS FILE
#Usefull function for the whole program
#---------------------------------------------------------------------

import Pkg

#---------------------------------------------------------------------
"""
	\tCheckPkg(PackagesList::Array{String})
	\t\tReturn Boolean
Check PackagesList to see if each specified package is already installed or not.\n
If not, the function will install the concerned package.
"""
function CheckPkg(p::Array{String})
	for i = 1:length(p)
		#Ctrl if the needed pkgs are arleady installed or not. DL if not.
		if(!haskey(Pkg.installed(), p[i]))
			println("* Installing "*p[i])
			if(Pkg.add(p[i]))
			else
				return false
			end
		else
			println("- "*p[i]*" already Installed")
		end
	end
	return true
end
#---------------------------------------------------------------------
"""
	\tCheckFiles(FilesList::Array{String})\n
	\t\tReturn Boolean
Check FilesList to see if each one is into the current directory or not.\n
"""
function CheckFiles(fl::Array{String})
	OK_array = map(x->(isfile(x)||isdir(x)),fl)
	#Total True == Number of Files ?
	if(sum(OK_array) == length(fl))
		println(sum(OK_array),"/",length(fl),"  files in the directory")
	else
		println("Warning, missing files !")
		for i in 1:length(OK_array)
			if(OK_array[i] == false)
				println("\t-> "*fl[i])
			end
		end
		return false
	end
	return true
end

#---------------------------------------------------------------------
function secure(status::Bool)
	#Allow an error to raise without crashing the software.
	asw = ""
	if(!status)
		println("\nWould you like to solve the issue and (R)etry ? Or (E)xit the program ?")
		asw = chomp(readline())
		if((asw[1]=='R')||(asw[1]=='r'))
			return true
		else
			println("\n---  P R O G R A M    E N D I N G  ---")
			exit()
		end
	end
	return false
end


text:
### A Pluto.jl notebook ###
# v0.18.4

using Markdown
using InteractiveUtils

# ╔═╡ 4067e9ab-1693-4c58-b0ab-99983f90b727
begin
    import Pkg
    # careful: this is _not_ a reproducible environment
    # activate the global environment
    Pkg.activate()

    using HetaSimulator, StatsPlots, NaNMath
end

# ╔═╡ 7021a4a7-86cf-4866-9790-6222cf01e896
md"# Dynamic QSP reporting with Julia and Pluto.jl"

# ╔═╡ b2148ea3-4200-43a6-9449-92f3fad1b67a
md"""
_This is an example of creation of a dynamic QSP report in Julia and Pluto.jl build based on HetaSimulator.jl. The content of this file and modeling platform is published in GitHub repository <https://github.com/insysbio/insulin-signaling-t2d>_
"""

# ╔═╡ 63d5bfbb-f86c-47a9-9c3e-4480ace0cedb
md"""
## Preamble

The QSP model which is used as an example of QSP model was published in the article:

>Brannmark C, Nyman E, Fagerholm S, Bergenholm L, Ekstrand EM, Cedersund G, Stralfors P. Insulin Signaling in Type 2 Diabetes: Experimental and modeling analyses reveal mechanisms of insulin resistance in human adipocytes. Journal of biological chemistry.. 2013 288(14):9867–9880. DOI: 10.1074/jbc.M112.432062

The SBML code was downloaded from BioModels <https://www.ebi.ac.uk/biomodels/BIOMD0000000448> and used as the part of the [Heta-based](https://hetalang.github.io/#/) modeling platform.

The report includes the steps to reproduce simulations from the original article demonstration facilities of the approach and necessary setups.

All necessary files can also be found in [the repository](https://github.com/insysbio/insulin-signaling-t2d).
"""

# ╔═╡ 0f2fec9b-6594-4cbe-8155-f48a4dab3fb7
md"""
## Preparations

For easier creation of the Heta-based platform install [heta compiler](https://hetalang.github.io/#/heta-compiler/?id=installation).

In command line interface run the code below to create heta platform template

```sh
heta init
```
The minimal content will be created.

Download the SBML model [from the database](https://www.ebi.ac.uk/biomodels/BIOMD0000000448) and copy it into `src/BIOMD0000000448.xml`
"""

# ╔═╡ 3d9e61b0-a65b-45f7-a7d9-0b4ef108bbb7
md"""
Update the `src/index.heta` with the following content:

```heta
// load SBML model as a content of the platform
include BIOMD0000000448_url.xml type sbml

// make the records of a model observable
block {output: true} begin
  measuredIRp;
  measuredIRint;
  measuredIRS1p;
  measuredIRS1307;
  measuredPKB308p;
  measuredPKB473p;
  measuredAS160p;
  measuredmTORC1a;
  measuredS6Kp;
  glucoseuptake;
  measuredmTORC2a;
  measuredS6p;
end

// make insulin as an input of the model
insulin @Const = 10; // nM

// make IR (insuline resistance) specific parameters as an input of the model
IR_total @Const = 100;    // reduce to 55%
GLUT4_total @Const = 100; // GLUT4 reduce to 50%
diabetes @Const = 1;      // reduce to 15%

// recalculate initial values for IR and base conditions
IR .= IR_total * 99.87/100; // 99.87
IRp .= 0;
IRins .= 0;
IRip .= IR_total * 0.02/100;      // 0.02
IRi .= IR_total * 0.11/100;       // 0.11
//
GLUT4 .= GLUT4_total * 73.48/100;   // 73.48
GLUT4m .= GLUT4_total * 26.52/100;  // 26.52

// variable parameters
k1a @Const = 0.6331;
k1basal @Const = 0.03683;
k1c @Const = 0.8768;
k1d @Const = 31.01;
```
"""

# ╔═╡ 0c9ee559-bb96-45e9-bd23-49d70e4a82f2
md"""
Install [Julia](https://julialang.org/downloads/) and add all necessary packages in Julia's console:

```julia
] add HetaSimulator StatsPlots NaNMath Pluto
```

Run Pluto notebook engine

```julia
using Pluto
Pluto.run()
```
"""

# ╔═╡ a9bc4164-4d07-4ab6-9e70-cfbcd20fedac
md"""
## Loading platform
"""

# ╔═╡ f25f0e29-7d7d-4496-b4a7-5ad519e06bba
md"""
Pluto.jl has a specific way to load Julia's modules. The default way by `load_platform` doesn't work.

The only way to load is splitting the whole process to three steps as shown: (1) compilation, (2) file update, (3) loading.

Step #2 inputs `using NaNMath` in the beginning as follows:

```julia
#=
    This code was generated by heta-compiler 0.6.12
=#
using NaNMath

__platform__ = (function()
...
```
"""

# ╔═╡ 063c8aa5-4751-4bb0-9c21-915de8215282
begin
	heta_build("..", julia_only=true, dist_dir="."); # 1. compilation of heta platform
	run(`sed -i '1s/^/using NaNMath \n/' ../_julia/model.jl`); # 2. _julia/model.jl update
	p = load_jlplatform("../_julia/model.jl", rm_out = true)  # 3. loading
end

# ╔═╡ 4218cca2-34ac-4108-a9ee-386f4a0d5662
md"""
Than the platform can be loaded.
"""

# ╔═╡ a783c6e9-a745-4e37-b745-e8366101cf80
m = models(p)[:nameless]

# ╔═╡ 4d7d9426-0743-4a4c-a364-25540f3798ca
md"""
## Default simulation

To simulate the observables only one line of code is required. Here we are using "chain" with `|>` syntax.
"""

# ╔═╡ e1bbe112-09eb-4f7c-a2a9-2a0cba056793
Scenario(m, (0.,30.)) |> sim |> plot

# ╔═╡ 26d95896-75a8-4d85-b935-3f05361ef1bf
md"Based on the original publication the insuline resistance (IR) can be set by updating three parameters."

# ╔═╡ cd5de49c-7d1d-4649-b57e-0fe95aad1f77
Scenario(m, (0.,30.), parameters = [:IR_total=>55., :GLUT4_total=>50, :diabetes=>0.15]) |> sim |> plot

# ╔═╡ da84dbae-be88-4638-956e-055ef056abee
md"""
## Simulation scenarios

Two simulation scenarios above can be created in tabular format and loaded into platform `p`. Create CSV file with the following content:

id | parameters.IR_total | parameters.GLUT4_total | parameters.diabetes | tspan
---|---|---|---|---
base | | | | 30
ir | 55 | 50 | 0.15| 30
"""

# ╔═╡ c0324eb1-1260-47b9-b23e-b0f3f49929c8
scn_csv = read_scenarios("../julia/scenarios.csv") # set csv file location here

# ╔═╡ 22e1278f-56f5-4f15-9666-e7ec538b3bba
add_scenarios!(p, scn_csv)

# ╔═╡ 3ab59992-95c6-4c21-bfda-ff0592a66f76
p

# ╔═╡ 9b23a1a8-e4f4-4b7c-ad1c-ceb110d341ee
md"""
## Advanced visualization

If you need to reshape and visualize different simulation results the easier way is to convert simulation results into `DataFrame`.
"""

# ╔═╡ e86db272-cf1c-4e80-bf43-bbb3ff6ee20b
results_df = sim(p) |> DataFrame # simulate all scenarios, then convert them

# ╔═╡ 2a715514-7c2c-457a-8d69-1a8d8046af5e
md"The following simulations reproduce the figures from the original paper."

# ╔═╡ f236c9c5-b720-4571-a4ff-ec0e12a89795
@df results_df plot(:t, :measuredIRp, group = :scenario, title = "Fig5 1B")

# ╔═╡ e2cca6b5-e9a1-449a-963c-c2c840f14b2e
@df results_df plot(:t, :measuredIRint, group = :scenario, title = "Fig5 1C")

# ╔═╡ a6d9b83b-9fde-4f67-ae3f-a102e5e25b4f
@df results_df plot(:t, :measuredIRS1p, group = :scenario, title = "Fig5 2B")

# ╔═╡ a29e5bee-b745-4bb9-8506-0f8025f7cac0
@df results_df plot(:t, :measuredIRS1307, group = :scenario, title = "Fig5 2C")

# ╔═╡ 80bac958-c5c0-4293-b480-7b3bb095ed49
@df results_df plot(:t, :measuredPKB308p, group = :scenario, title = "Fig5 3B")

# ╔═╡ b75a019f-2b55-4a6f-a8a7-70435a15f0fb
@df results_df plot(:t, :measuredPKB473p, group = :scenario, title = "Fig5 3C")

# ╔═╡ 45344e77-b326-408e-81d7-11de72852eed
md"""
## Titration-like simulations

The simulations of another type (not time dependence) can be performed and visualized by applying different simulation scenarios. For example in the original article the titration-like experiment is simulated: intake of insulin and measurements of different observables after 10 minutes.

To reproduce them the following set of scenarion should be created:

id | parameters.insulin | parameters.IR_total | parameters.GLUT4_total | parameters.diabetes | tspan | parameters.is_titr
---|---|---|---|---|---|---
base_ins_3 | 1.00E-03 | | | | 10 | 1
base_ins_25 | 3.16E-03 | | | | 10 | 1
base_ins_2 | 1.00E-02 | | | | 10 | 1
base_ins_15 | 3.16E-02 | | | | 10 | 1
base_ins_1 | 1.00E-01 | | | | 10 | 1
base_ins_05 | 3.16E-01 | | | | 10 | 1
base_ins0 | 1.00E+00 | | | | 10 | 1
base_ins05 | 3.16E+00 | | | | 10 | 1
base_ins1 | 1.00E+01 | | | | 10 | 1
base_ins15 | 3.16E+01 | | | | 10 | 1
base_ins2 | 1.00E+02 | | | | 10 | 1
ir_ins_3 | 1.00E-03 | 55 | 50 | 0.15| 10 | 1
ir_ins_25 | 3.16E-03 | 55 | 50 | 0.15| 10 | 1
ir_ins_2 | 1.00E-02 | 55 | 50 | 0.15| 10 | 1
ir_ins_15 | 3.16E-02 | 55 | 50 | 0.15| 10 | 1
ir_ins_1 | 1.00E-01 | 55 | 50 | 0.15| 10 | 1
ir_ins_05 | 3.16E-01 | 55 | 50 | 0.15| 10 | 1
ir_ins0 | 1.00E+00 | 55 | 50 | 0.15| 10 | 1
ir_ins05 | 3.16E+00 | 55 | 50 | 0.15| 10 | 1
ir_ins1 | 1.00E+01 | 55 | 50 | 0.15| 10 | 1
ir_ins15 | 3.16E+01 | 55 | 50 | 0.15| 10 | 1
ir_ins2 | 1.00E+02 | 55 | 50 | 0.15| 10 | 1
"""

# ╔═╡ e5721cf5-e4ba-4d45-93b8-d228aab15da0
scn_titr_csv = read_scenarios("../julia/titration.csv")

# ╔═╡ 915c245d-9b03-4fb8-a27a-dc3b958b7b9a
add_scenarios!(p, scn_titr_csv)

# ╔═╡ b24ce74c-ffff-4d46-824a-2b853a006369
results_titr = sim(p, saveat=[10.]);

# ╔═╡ 71c93c6c-8390-4df7-ade7-e23eef29f5ec
results_titr_df = DataFrame(results_titr, add_parameters=true);

# ╔═╡ 47ba0050-6f67-4027-9003-2470cdaced03
results_titr_subset = subset(results_titr_df, :is_titr=>x->x.===1.)

# ╔═╡ 07f91e70-aee0-459e-ba82-1dd01dbbfbf4
@df results_titr_subset plot(
    :insulin, :measuredmTORC2a,
    title = "Fig5 1A",
    xlabel="insulin, nM", ylabel="mTORC2a measured",
    group=:diabetes,
    xscale = :log10,
    legend = :topleft
)

# ╔═╡ 10c479f3-8c8e-4255-9c7c-48fefa78ce7a
md"""
## Multiple simulations

HetaSimulator includes the mechanism to run Monte-Carlo simulations based on parameter variability. For the demonstration purposes we will generate a random set of parameters: `k1a, k1basal, k1c, k1d`. This simulations mimic the uncertainty in the selected parameters.
"""

# ╔═╡ 28cdb580-be40-40d7-b050-3648190a5d83
results_mc = mc(
    p,
    [
        :k1a => LogNormal(0.6331, 0.5),
        :k1basal => LogNormal(0.03683, 0.5),
        :k1c => LogNormal(0.8768, 0.5),
        :k1d => LogNormal(31.01, 0.5)
    ],
    100,                      # number of
    scenarios = [:base, :ir], # scenarios names to simulate
    saveat = 0:(0.1):30,      # time points
    abstol = 1e-7
)

# ╔═╡ bf064fc2-3d1d-4f62-8218-a1bf3451d1a9
plot(results_mc, vars=[:measuredIRp, :measuredIRint])

# ╔═╡ b30af5e7-35a8-4d74-8ca9-dbfa75fe4f28
md"The generated dataset can be summarized to plot the median and 90% bounds for uncertainty predictions."

# ╔═╡ 7bb1add5-3600-4902-832c-6903e3815b5a
summary_mc = EnsembleSummary(results_mc; quantiles=[0.05,0.95]);

# ╔═╡ 3fd1c38b-472f-4afd-9c41-65f364c078fc
plot(summary_mc, idxs=[1,2]) # to plot observables 1 - IRp, 2 - IRint

# ╔═╡ Cell order:
# ╟─7021a4a7-86cf-4866-9790-6222cf01e896
# ╟─b2148ea3-4200-43a6-9449-92f3fad1b67a
# ╟─63d5bfbb-f86c-47a9-9c3e-4480ace0cedb
# ╟─0f2fec9b-6594-4cbe-8155-f48a4dab3fb7
# ╟─3d9e61b0-a65b-45f7-a7d9-0b4ef108bbb7
# ╟─0c9ee559-bb96-45e9-bd23-49d70e4a82f2
# ╟─a9bc4164-4d07-4ab6-9e70-cfbcd20fedac
# ╠═4067e9ab-1693-4c58-b0ab-99983f90b727
# ╟─f25f0e29-7d7d-4496-b4a7-5ad519e06bba
# ╠═063c8aa5-4751-4bb0-9c21-915de8215282
# ╟─4218cca2-34ac-4108-a9ee-386f4a0d5662
# ╠═a783c6e9-a745-4e37-b745-e8366101cf80
# ╟─4d7d9426-0743-4a4c-a364-25540f3798ca
# ╠═e1bbe112-09eb-4f7c-a2a9-2a0cba056793
# ╟─26d95896-75a8-4d85-b935-3f05361ef1bf
# ╠═cd5de49c-7d1d-4649-b57e-0fe95aad1f77
# ╟─da84dbae-be88-4638-956e-055ef056abee
# ╠═c0324eb1-1260-47b9-b23e-b0f3f49929c8
# ╠═22e1278f-56f5-4f15-9666-e7ec538b3bba
# ╠═3ab59992-95c6-4c21-bfda-ff0592a66f76
# ╟─9b23a1a8-e4f4-4b7c-ad1c-ceb110d341ee
# ╠═e86db272-cf1c-4e80-bf43-bbb3ff6ee20b
# ╟─2a715514-7c2c-457a-8d69-1a8d8046af5e
# ╠═f236c9c5-b720-4571-a4ff-ec0e12a89795
# ╠═e2cca6b5-e9a1-449a-963c-c2c840f14b2e
# ╠═a6d9b83b-9fde-4f67-ae3f-a102e5e25b4f
# ╠═a29e5bee-b745-4bb9-8506-0f8025f7cac0
# ╠═80bac958-c5c0-4293-b480-7b3bb095ed49
# ╠═b75a019f-2b55-4a6f-a8a7-70435a15f0fb
# ╟─45344e77-b326-408e-81d7-11de72852eed
# ╠═e5721cf5-e4ba-4d45-93b8-d228aab15da0
# ╠═915c245d-9b03-4fb8-a27a-dc3b958b7b9a
# ╠═b24ce74c-ffff-4d46-824a-2b853a006369
# ╠═71c93c6c-8390-4df7-ade7-e23eef29f5ec
# ╠═47ba0050-6f67-4027-9003-2470cdaced03
# ╠═07f91e70-aee0-459e-ba82-1dd01dbbfbf4
# ╟─10c479f3-8c8e-4255-9c7c-48fefa78ce7a
# ╠═28cdb580-be40-40d7-b050-3648190a5d83
# ╠═bf064fc2-3d1d-4f62-8218-a1bf3451d1a9
# ╟─b30af5e7-35a8-4d74-8ca9-dbfa75fe4f28
# ╠═7bb1add5-3600-4902-832c-6903e3815b5a
# ╠═3fd1c38b-472f-4afd-9c41-65f364c078fc


text:
import JSON
import RecursiveArrayTools
global filename = "../json/1-move.json"
function buildapuzzle(filename)
    try
        file = JSON.parsefile(filename)
        _start = VectorOfArray(file["start"])
        global start = convert(Array, _start)
        println("Start: $start")
        _goal = VectorOfArray(file["goal"])
        global goal = convert(Array, _goal)
        println("Goal: $goal")
        global n = file["n"]
        println("n: $n")
    catch e
        println("ERROR: Faulty JSON, try again.")
        return
    end
end

max_stars_repo_name: brice-allen/npuz

max_issues_repo_name: brice-allen/npuz

max_forks_repo_name: brice-allen/npuz


text:
"""
    nearest_node(g::OSMGraph,
                 points::Vector{GeoLocation},
                 n_neighbours::Integer=1
                 )::Tuple{Vector{Vector{Integer}}, Vector{Vector{AbstractFloat}}}

Finds nearest nodes from a vector of GeoLocations using a `NearestNeighbors.jl` KDTree.

# Arguments
- `g::OSMGraph`: Graph container.
- `points::Vector{GeoLocation}`: Vector of query points.
- `n_neighbours::Integer`: Number of neighbours to query for each point.

# Return
- Tuple of neighbours and straight line euclidean distances from each point `([[neighbours]...], [[dists]...])`.
"""
function nearest_node(g::OSMGraph{U, T, W},
                      points::Vector{GeoLocation},
                      n_neighbours::Integer=1
                      )::Tuple{Vector{Vector{T}},Vector{Vector{Float64}}} where {U, T, W}
    cartesian_locations = LightOSM.to_cartesian(points)
    idxs, dists = knn(g.kdtree, cartesian_locations, n_neighbours, true)
    neighbours = [[g.index_to_node[j] for j in i] for i in idxs]
    dists = collect(dists)
    return  neighbours, dists
end

nearest_node(g::OSMGraph, point::GeoLocation, n_neighbours::Integer=1) = nearest_node(g, [point], n_neighbours)
nearest_node(g::OSMGraph, points::Vector{<:Vector{<:AbstractFloat}}, n_neighbours::Integer=1) = nearest_node(g, GeoLocation(points), n_neighbours)
nearest_node(g::OSMGraph, point::Vector{<:AbstractFloat}, n_neighbours::Integer=1) = nearest_node(g, GeoLocation(point), n_neighbours)

"""
Finds nearest nodes from a vector of node ids, the query node id is excluded from neighbours.
"""
function nearest_node(g::OSMGraph{U, T, W},
                      nodes::Vector{<:Integer},
                      n_neighbours::Integer=1
                      )::Tuple{Vector{Vector{T}},Vector{Vector{Float64}}} where {U, T, W}
    locations = [g.nodes[n].location for n in nodes]
    n_neighbours += 1 # Closest node is always the input node itself, exclude self from result
    neighbours, dists = nearest_node(g, locations, n_neighbours)
    return [n[2:end] for n in neighbours], [d[2:end] for d in dists]
end

nearest_node(g::OSMGraph, node::Integer, n_neighbours::Integer=1) = nearest_node(g, [node], n_neighbours)
nearest_node(g::OSMGraph, nodes::Vector{<:Node}, n_neighbours::Integer=1) = nearest_node(g, [n.id for n in nodes], n_neighbours)
nearest_node(g::OSMGraph, node::Node, n_neighbours::Integer=1) = nearest_node(g, node.id, n_neighbours)


text:
abstract type MeanFunction end

"""
    ZeroMean{T<:Real} <: MeanFunction

Returns `zero(T)` everywhere.
"""
struct ZeroMean{T<:Real} <: MeanFunction end

"""
This is an AbstractGPs-internal workaround for AD issues; ideally we would just extend Base.map
"""
_map_meanfunction(::ZeroMean{T}, x::AbstractVector) where {T} = Zeros{T}(length(x))

ZeroMean() = ZeroMean{Float64}()

"""
    ConstMean{T<:Real} <: MeanFunction

Returns `c` everywhere.
"""
struct ConstMean{T<:Real} <: MeanFunction
    c::T
end

_map_meanfunction(m::ConstMean, x::AbstractVector) = Fill(m.c, length(x))

"""
    CustomMean{Tf} <: MeanFunction

A wrapper around whatever unary function you fancy. Must be able to be mapped over an
`AbstractVector` of inputs.
"""
struct CustomMean{Tf} <: MeanFunction
    f::Tf
end

_map_meanfunction(m::CustomMean, x::AbstractVector) = map(m.f, x)

max_stars_repo_name: JuliaGaussianProcesses/AbstractGP

max_issues_repo_name: JuliaGaussianProcesses/AbstractGP

max_forks_repo_name: JuliaGaussianProcesses/AbstractGP


text:
using LinearAlgebra, Test, WaveFD

@testset "compressedio, 3D Array, T=$(T), C=$(C), style=$(style)" for (T,C) in ((Float32, UInt32), (Float64, UInt32), (Float32, Float32), (Float64, Float64), (Float64, Float32)), style in (:file, :buffer)
    nz = 512
    ny = 256
    nx = 128

    bz = 32
    by = 32
    bx = 32

    x = -1 .+ 2 .* rand(T,nz,ny,nx)
    y = -1 .+ 2 .* rand(T,nz,ny,nx)

    C,P = WaveFD.comptype(C, T)
    comp = WaveFD.Compressor(T, P, C, (nz,ny,nx), (bz,by,bx), 1e-2, 2, false)
    open(comp)
    io = style == :file ? open("test.bin","w") : IOBuffer()
    WaveFD.compressedwrite(io, comp, 1, x)
    WaveFD.compressedwrite(io, comp, 2, y)
    style == :file && close(io)

    xx = similar(x)
    yy = similar(y)

    io = style == :file ? open("test.bin") : io
    WaveFD.compressedread!(io, comp, 1, xx)
    WaveFD.compressedread!(io, comp, 2, yy)
    close(io)

    snr_x = 10*log10(norm(x)^2 / norm(x .- xx)^2)
    snr_y = 10*log10(norm(y)^2 / norm(y .- yy)^2)
    err_x = norm(x .- xx) / length(x)
    err_y = norm(y .- yy) / length(y)
    if C == UInt32
        @test snr_x > 40
        @test snr_y > 40
    else
        @test err_x < 1e-6
        @test err_y < 1e-6
    end
    close(comp)
end

@testset "compressedio, 3D SubArray, T=$(T), C=$(C), style=$(style)" for (T,C) in ((Float32, UInt32), (Float64, UInt32), (Float32, Float32), (Float64, Float64), (Float64, Float32)), style in (:file, :buffer)
    nz = 512
    ny = 256
    nx = 128

    bz = 32
    by = 32
    bx = 32

    x = -1 .+ 2 .* rand(T,nz,ny,nx)
    y = -1 .+ 2 .* rand(T,nz,ny,nx)

    rng = (10:nz-10,10:ny-10,10:nx-10)

    x_sub = @view x[rng...]
    y_sub = @view y[rng...]

    nz, ny, nx = size(x_sub)

    C,P = WaveFD.comptype(C, T)
    comp = WaveFD.Compressor(T, P, C, (nz,ny,nx), (bz,by,bx), 1e-2, 2, true)
    open(comp)

    io = style == :file ? open("test.bin", "w") : IOBuffer()
    WaveFD.compressedwrite(io, comp, 1, x, rng)
    WaveFD.compressedwrite(io, comp, 2, y, rng)
    style == :file && close(io)

    xx = copy(x)
    yy = copy(y)

    xx_sub = @view xx[rng...]
    yy_sub = @view yy[rng...]

    xx_sub .= 0.0
    yy_sub .= 0.0

    io = style == :file ? open("test.bin") : io
    WaveFD.compressedread!(io, comp, 1, xx, rng)
    WaveFD.compressedread!(io, comp, 2, yy, rng)
    close(io)

    snr_x = 10*log10(norm(x)^2 / norm(x .- xx)^2)
    snr_y = 10*log10(norm(y)^2 / norm(y .- yy)^2)
    err_x = norm(x .- xx) / length(x)
    err_y = norm(y .- yy) / length(y)
    if C == UInt32
        @test snr_x > 40
        @test snr_y > 40
    else
        @test err_x < 1e-6
        @test err_y < 1e-6
    end
    close(comp)
end

@testset "compressedio, 2D Array, T=$(T), C=$(C), style=$(style)" for (T,C) in ((Float32, UInt32), (Float64, UInt32), (Float32, Float32), (Float64, Float64), (Float64, Float32)), style in (:file, :buffer)
    nz = 512
    nx = 128

    bz = 32
    bx = 32

    x = -1 .+ 2 .* rand(T,nz,nx)
    y = -1 .+ 2 .* rand(T,nz,nx)

    C,P = WaveFD.comptype(C, T)
    comp = WaveFD.Compressor(T, P, C, (nz,nx), (bz,bx), 1e-2, 2, false)
    open(comp)

    io = style == :file ? open("test.bin", "w") : IOBuffer()
    WaveFD.compressedwrite(io, comp, 1, x)
    WaveFD.compressedwrite(io, comp, 2, y)
    style == :file && close(io)

    xx = similar(x)
    yy = similar(y)

    io = style == :file ? open("test.bin") : io
    WaveFD.compressedread!(io, comp, 1, xx)
    WaveFD.compressedread!(io, comp, 2, yy)
    close(io)

    snr_x = 10*log10(norm(x)^2 / norm(x .- xx)^2)
    snr_y = 10*log10(norm(y)^2 / norm(y .- yy)^2)
    err_x = norm(x .- xx) / length(x)
    err_y = norm(y .- yy) / length(y)
    if C == UInt32
        @test snr_x > 40
        @test snr_y > 40
    else
        @test err_x < 1e-6
        @test err_y < 1e-6
    end
    close(comp)
end

@testset "compressedio 2D SubArray, T=$(T), C=$(C), style=$(style)" for (T,C) in ((Float32, UInt32), (Float64, UInt32), (Float32, Float32), (Float64, Float64), (Float64, Float32)), style in (:file, :buffer)
    nz = 512
    nx = 128

    bz = 32
    bx = 32

    x = -1 .+ 2 .* rand(T,nz,nx)
    y = -1 .+ 2 .* rand(T,nz,nx)

    rng = (10:nz-10,10:nx-10)

    x_sub = @view x[rng...]
    y_sub = @view y[rng...]

    nz, nx = size(x_sub)

    C,P = WaveFD.comptype(C, T)
    comp = WaveFD.Compressor(T, P, C, (nz,nx), (bz,bx), 1e-2, 2, true)
    open(comp)

    io = style == :file ? open("test.bin", "w") : IOBuffer()
    WaveFD.compressedwrite(io, comp, 1, x, rng)
    WaveFD.compressedwrite(io, comp, 2, y, rng)
    style == :file && close(io)

    xx = copy(x)
    yy = copy(y)

    xx_sub = @view xx[rng...]
    yy_sub = @view yy[rng...]

    xx_sub[:] .= 0.0
    yy_sub[:] .= 0.0

    io = style == :file ? open("test.bin") : io
    WaveFD.compressedread!(io, comp, 1, xx, rng)
    WaveFD.compressedread!(io, comp, 2, yy, rng)
    close(io)

    snr_x = 10*log10(norm(x)^2 / norm(x .- xx)^2)
    snr_y = 10*log10(norm(y)^2 / norm(y .- yy)^2)
    err_x = norm(x .- xx) / length(x)
    err_y = norm(y .- yy) / length(y)
    if C == UInt32
        @test snr_x > 40
        @test snr_y > 40
    else
        @test err_x < 1e-6
        @test err_y < 1e-6
    end
    close(comp)
end

@testset "compressedio, exception" for file in ("test.bin",)
    io = open(file,"w")
    @test_throws ErrorException WaveFD.compressedwrite_exception(io,0,1000)
    rm(file)
end

nothing

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
#=
inference:
- Julia version:
- Author: marcin
- Date: 2018-06-02
=#

# based on https://www.cs.helsinki.fi/u/bmmalone/probabilistic-models-spring-2014/Huang1996.pdf

# TODO generify types


text:
using Documenter, MonteCarlo

makedocs(
    modules = [MonteCarlo],
    doctest = false,
    sitename = "MonteCarlo.jl",
    pages = [
        "Introduction" => "index.md",
        "Manual" => [
            "Getting started" => "manual/gettingstarted.md",
            "Showcase" => "manual/showcase.md",
        ],
        "Physical models" => [
            "Ising model" => "models/ising.md",
            "Attractive Hubbard model" => "models/hubbardattractive.md",
        ],
        "Monte Carlo flavors" => [
            "MC" => "flavors/mc.md",
            "DQMC" => "flavors/dqmc.md",
        ],
        "Lattices" => "lattices.md",
        "Customize" => "customize.md",
        "Interfaces" => [
            "MC" => "interfaces/MC.md",
            "DQMC" => "interfaces/DQMC.md",
        ],
        "General exports" => ["General" => "methods/general.md"],
    ],
    # assets = ["assets/custom.css", "assets/custom.js"]
)

deploydocs(
    repo = "github.com/crstnbr/MonteCarlo.jl.git",
    push_preview = true,
    # target = "site",
)


text:
include("tile_entity.jl")

# Use the Maple function
# Then set the data from the placements (such as setting fields IE winged for strawberries)
struct EntityPlacement
    func::Union{Function, Type}
    placement::String
    data::Dict{String, Any}
    finalizer::Union{Function, Nothing}

    function EntityPlacement(func::Union{Function, Type}, placement::String="point", data::Dict{String, Any}=Dict{String, Any}(), finalizer::Union{Function, Nothing}=nothing)
        return new(func, placement, data, finalizer)
    end
end

const PlacementDict = Dict{String, EntityPlacement}

entityNameLookup = Dict{String, String}()

# The point the entity will be drawn relative too
# Make sure it returns in Int32/64
function position(entity::Maple.Entity)
    return (
        floor(Int, get(entity.data, "x", 0)),
        floor(Int, get(entity.data, "y", 0))
    )
end

# Entity is considered to have failed rendering if it returns exactly `false`, `nothing` is considered a success.
function renderCall(func, args...)
    if applicable(func, args...)
        if func(args...) != false
            return true
        end
    end

    return false
end

function attemptEntityRender(ctx::Cairo.CairoContext, layer::Layer, entity::Maple.Entity, room::Maple.Room)
    Cairo.save(ctx)

    res = renderCall(renderAbs, ctx, entity) ||
        renderCall(renderAbs, ctx, entity, room)

    translate(ctx, position(entity)...)
    res |= renderCall(render, ctx, entity) ||
        renderCall(render, ctx, entity, room)

    Cairo.restore(ctx)

    return res
end

function render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Entity, room::Maple.Room)
    return false
end

function renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Entity, room::Maple.Room)
    return false
end

function renderEntity(ctx::Cairo.CairoContext, layer::Layer, entity::Maple.Entity, room::Maple.Room; alpha::Number=1)
    success = false

    # Set global alpha here, passing alpha to the entity renderer is not sane
    setGlobalAlpha!(alpha)

    success = attemptEntityRender(ctx, layer, entity, room)

    # Reset global alpha again
    setGlobalAlpha!(1)

    if !success
        debug.log("Couldn't render entity '$(entity.name)' in room '$(room.name)'", "DRAWING_ENTITY_MISSING")
    end
end

function renderSelected(ctx::Cairo.CairoContext, entity::Maple.Entity)
    return false
end

function renderSelectedAbs(ctx::Cairo.CairoContext, entity::Maple.Entity)
    return false
end

function attemptEntitySelectionRender(ctx::Cairo.CairoContext, layer::Layer, entity::Maple.Entity, room::Maple.Room)
    Cairo.save(ctx)

    res = renderCall(renderSelectedAbs, ctx, entity) ||
        renderCall(renderSelectedAbs, layer, entity) ||
        renderCall(renderSelectedAbs, ctx, entity, room) ||
        renderCall(renderSelectedAbs, layer, entity, room)

    translate(ctx, position(entity)...)
    res |=
        renderCall(renderSelected, ctx, entity) ||
        renderCall(renderSelected, layer, entity) ||
        renderCall(renderSelected, ctx, entity, room) ||
        renderCall(renderSelected, layer, entity, room)

    Cairo.restore(ctx)

    return res
end

function renderEntitySelection(ctx::Cairo.CairoContext, layer::Layer, entity::Maple.Entity, room::Maple.Room; alpha::Number=1)
    # Set global alpha here, passing alpha to the entity renderer is not sane
    setGlobalAlpha!(alpha)

    success = attemptEntitySelectionRender(ctx, layer, entity, room)

    # Reset global alpha again
    setGlobalAlpha!(1)
end

# Helper function to respect debug option
function minimumSizeWrapper(target::Union{Maple.Entity, Maple.Trigger})
    if get(debug.config, "IGNORE_MINIMUM_SIZE", false)
        return 0, 0
    end

    return minimumSize(target)
end

# Helper function to respect debug option
function canResizeWrapper(target::Union{Maple.Entity, Maple.Trigger})
    if get(debug.config, "IGNORE_CAN_RESIZE", false)
        return true, true
    end

    return resizable(target)
end

minimumSize(entity::Maple.Entity) = 8, 8
resizable(entity::Maple.Entity) = false, false

nodeLimits(entity::Maple.Entity) = 0, 0

editingOptions(entity::Maple.Entity) = Dict{String, Any}()

function registerPlacements!(placements::PlacementDict, loaded::Array{String, 1})
    empty!(placements)

    for modul in loaded
        if hasModuleField(modul, "placements")
            merge!(placements, getModuleField(modul, "placements"))
        end
    end
end

function callFinalizer(map::Maple.Map, room::Room, ep::EntityPlacement, target::Union{Maple.Entity, Maple.Trigger})
    argsList = Tuple[
        (target, map, room),
        (target, room),
        (target,)
    ]

    for args in argsList
        if applicable(ep.finalizer, args...)
            return ep.finalizer(args...)
        end
    end
end

function updateEntityPosition!(target::Union{Entity, Trigger}, ep::EntityPlacement, map::Maple.Map, room::Room, x::Integer, y::Integer, nx::Integer=x + 8, ny::Integer=y + 8)
    merge!(target.data, ep.data)

    if ep.placement == "point"
        target.x = x
        target.y = y

    elseif ep.placement == "rectangle"
        rect = selectionRectangle(x, y, nx, ny)

        resizeW, resizeH = canResizeWrapper(target)
        minW, minH = minimumSizeWrapper(target)

        w = resizeW ? max(minW, rect.w - 1) : minW
        h = resizeH ? max(minH, rect.h - 1) : minH

        target.x = rect.x
        target.y = rect.y

        target.width = w
        target.height = h

    elseif ep.placement == "line"
        target.x = x
        target.y = y

        target.nodes = Tuple{Integer, Integer}[(nx, ny)]
    end

    if ep.finalizer !== nothing
        callFinalizer(map, room, ep, target)
    end

    return target
end

function generateEntity(map::Maple.Map, room::Room, ep::EntityPlacement, x::Integer, y::Integer, nx::Integer=x + 8, ny::Integer=y + 8)
    target = ep.func(x, y)

    return updateEntityPosition!(target, ep, map, room, x, y, nx, ny)
end

function updateCachedEntityPosition!(cache::Dict{String, T}, placements::PlacementDict, map::Maple.Map, room::Room, name::String, x::Integer, y::Integer, nx::Integer=x + 8, ny::Integer=y + 8) where T
    target = getCachedPlacement!(cache, placements, name)
    ep = placements[name]

    return updateEntityPosition!(target, ep, map, room, x, y, nx, ny)
end

function entityConfigOptions(entity::Union{Maple.Entity, Maple.Trigger}, ignores::Array{String, 1}=String[])
    addedNodes = false
    res = Form.Option[]

    options = editingOptions(entity)
    horizontalAllowed, verticalAllowed = canResizeWrapper(entity)
    nodeRange = nodeLimits(entity)

    key = isa(entity, Maple.Entity) ? "entities" : "triggers"
    tooltips = get(langdata, ["placements", key, entity.name, "tooltips"])
    names = get(langdata, ["placements", key, entity.name, "names"])

    # Add nothing keys for all Dict editing options
    # Merge entity data over afterwards, this makes it possible to "store" nothing values for editing later
    data = merge(
        Dict{String, Any}(
            attr => nothing for (attr, value) in options if isa(value, Dict)
        ),
        entity.data
    )

    for (attr, value) in data
        if !horizontalAllowed && attr == "width" || !verticalAllowed && attr == "height"
            continue
        end

        # Always ignore these, regardless of ignore argument
        if attr == "originX" || attr == "originY"
            continue
        end

        # Special cased below
        if attr == "nodes"
            continue
        end

        if attr in ignores
            continue
        end

        if get(debug.config, "TOOLTIP_ENTITY_MISSING", false) && !haskey(tooltips, Symbol(attr))
            if !(attr in debug.defaultIgnoredTooltipAttrs)
                println("Missing tooltip for '$(entity.name)' - $attr")
            end
        end

        name = expandTooltipText(get(names, Symbol(attr), ""))
        displayName = isempty(name) ? humanizeVariableName(attr) : name
        tooltip = expandTooltipText(get(tooltips, Symbol(attr), ""))
        attrOptions = get(options, attr, nothing)

        push!(res, Form.suggestOption(displayName, value, dataName=attr, tooltip=tooltip, choices=attrOptions, editable=true))
    end

    if nodeRange[2] != 0 && !("nodes" in ignores)
        push!(res, Form.ListOption("Node X;Node Y", get(entity.data, "nodes", Tuple{Int, Int}[]), dataName="nodes", minRows=nodeRange[1], maxRows=nodeRange[2]))
    end

    return res
end

function selection(entity::Maple.Entity)
    return false
end

function entitySelection(entity::Maple.Entity, room::Maple.Room, node::Integer=0)
    argsList = Tuple[
        (entity, room),
        (entity,)
    ]

    for args in argsList
        if applicable(selection, args...)
            return selection(args...)
        end
    end
end

function getSimpleEntityRng(entity::Maple.Entity)
    x, y = position(entity)
    seed = abs(x) << ceil(Int, log(2, abs(y) + 1)) | abs(y)

    return MersenneTwister(seed)
end

function getCachedPlacement!(cache::Dict{String, T}, placements::PlacementDict, name::String, map::Union{Maple.Map, Nothing}=Ahorn.loadedState.map, room::Union{Maple.Room, Nothing}=Ahorn.loadedState.room) where T <: Union{Entity, Trigger}
    if map !== nothing && room !== nothing
        if !haskey(cache, name)
            if haskey(placements, name)
                cache[name] = generateEntity(map, room, placements[name], 0, 0, 0, 0)
            end
        end
    end

    return get(cache, name, nothing)
end

function fillPlacementCache!(cache::Dict{String, T}, placements::PlacementDict, map::Union{Maple.Map, Nothing}=Ahorn.loadedState.map, room::Union{Maple.Room, Nothing}=Ahorn.loadedState.room) where T <: Union{Entity, Trigger}
    empty!(cache)

    for (name, placement) in placements
        getCachedPlacement!(cache, placements, name, map, room)
    end
end

const loadedEntities = joinpath.(abs"entities", readdir(abs"entities"))
loadModule.(loadedEntities)
const entityPlacements = PlacementDict()
const entityPlacementsCache = Dict{String, Entity}()


text:
# Types declarations for Krypto.jl
# Copyright (C): Jaka Smrekar (vinctux) <vinctux@outlook.com>, 2016

# Primitives
abstract CryptoAlgorithm
abstract SymmetricCrypto <: CryptoAlgorithm
abstract AsymmetricCrypto <: CryptoAlgorithm

# Asymmetric
abstract RSA <: AsymmetricCrypto
abstract ECC <: AsymmetricCrypto
abstract RLWE <: AsymmetricCrypto
abstract NTRU <: AsymmetricCrypto

# Symmetric
abstract AES <: SymmetricCrypto
abstract Blowfish <: SymmetricCrypto
abstract Twofish <: SymmetricCrypto
abstract Threefish <: SymmetricCrypto
abstract Serpent <: SymmetricCrypto
abstract Salsa20 <: SymmetricCrypto
abstract RABBIT <: SymmetricCrypto


text:
#=
  Bayesian Data Analysis of Iris dataset

5×7 DataFrame
 Row │ variable     mean     min     median  max        nmissing  eltype
     │ Symbol       Union…   Any     Union…  Any        Int64     DataType
─────┼────────────────────────────────────────────────────────────────────────────────────────────
   1 │ SepalLength  5.84333  4.3     5.8     7.9               0  Float64
   2 │ SepalWidth   3.05733  2.0     3.0     4.4               0  Float64
   3 │ PetalLength  3.758    1.0     4.35    6.9               0  Float64
   4 │ PetalWidth   1.19933  0.1     1.3     2.5               0  Float64
   5 │ Species               setosa          virginica         0  CategoricalValue{String, UInt8}

Mean value for each species: ["virginica", "setosa", "versicolor"]
 ("virginica", [6.587999999999999 2.9739999999999998 5.552 2.026])
 ("setosa", [5.006 3.428 1.462 0.24600000000000002])
 ("versicolor", [5.935999999999999 2.77 4.26 1.3259999999999998])

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

           a    1.1509    0.2081     0.0066    0.0106   484.2115    1.0090       97.6036
        c[1]   -0.1040    0.0589     0.0019    0.0026   439.0116    0.9990       88.4926
        c[2]   -0.0396    0.0595     0.0019    0.0036   384.6719    1.0042       77.5392
        c[3]    0.2247    0.0559     0.0018    0.0026   391.6563    1.0045       78.9471
        c[4]    0.6119    0.0907     0.0029    0.0046   380.7559    1.0086       76.7498
       sigma    0.2206    0.0133     0.0004    0.0007   604.2689    0.9992      121.8038

parameter: credible_intervals
Credible interval for a with mass 0.95: (0.762098..1.562943)
Credible interval for c[1] with mass 0.95: (-0.228162..-0.000480)
Credible interval for c[2] with mass 0.95: (-0.155606..0.083093)
Credible interval for c[3] with mass 0.95: (0.112843..0.334041)
Credible interval for c[4] with mass 0.95: (0.442641..0.800271)
Credible interval for sigma with mass 0.95: (0.194433..0.247127)

Confusion matrix:
3×3 Matrix{Float64}:
 50.0   0.0   0.0
  0.0  48.0   2.0
  0.0   2.0  48.0

Bad classfied:
4-element Vector{Any}:
  71 => (2 => 3)
  84 => (2 => 3)
 120 => (3 => 2)
 134 => (3 => 2)
num bad: 4

bads: credible_intervals
b:71 => (2 => 3)
x[71]: [5.9, 3.2, 4.8, 1.8]  true cat (y[71]): 2.0  (pred mean: 2.5812725160974286
Credible interval for y[71] with mass 0.3: (2.481020..2.645030)
b:84 => (2 => 3)
x[84]: [6.0, 2.7, 5.1, 1.6]  true cat (y[84]): 2.0  (pred mean: 2.548502053149339
Credible interval for y[84] with mass 0.3: (2.419795..2.590571)
b:120 => (3 => 2)
x[120]: [6.0, 2.2, 5.0, 1.5]  true cat (y[120]): 3.0  (pred mean: 2.479493724261082
Credible interval for y[120] with mass 0.3: (2.419840..2.583439)
b:134 => (3 => 2)
x[134]: [6.3, 2.8, 5.1, 1.5]  true cat (y[134]): 3.0  (pred mean: 2.442839883396137
Credible interval for y[134] with mass 0.3: (2.289724..2.465010)

Mean value for each species: 1:3
3-element Vector{Tuple{Int64, Matrix{Float64}}}:
 (1, [5.006 3.428 1.462 0.24600000000000002])
 (2, [5.935999999999999 2.77 4.26 1.3259999999999998])
 (3, [6.587999999999999 2.9739999999999998 5.552 2.026])

=#
using Turing, StatsPlots, Distributions
using DataFrames,RDatasets
include("jl_utils.jl")

# Read dataset
iris = dataset("datasets", "iris")

# Convert Categorical to integer
conv = Dict(
    "setosa" => 1,
    "versicolor" => 2,
    "virginica" => 3,
)
display(describe(iris))
iris2 = Matrix(transform(iris,[:Species] => ByRow(i->conv[i]))[:,[1,2,3,4,6]])
x = iris2[:,1:4]
y = iris2[:,5]

println("Mean value for each species: ", keys(conv))
display(keys(conv) .|> s -> (s,mean(x[iris.Species.==s,:], dims=1)))
println("\n")

@model function iris_model(x,y=missing, ::Type{T} = Float64) where {T}
    n,p = size(x)

    if y === missing
      y = Vector{T}(undef, n)
    end

    a ~ Normal(0,1)
    # Parameters:
    # - sepalLength
    # - sepalWidth
    # - petalLength
    # - petalWidth
    params ~ filldist(Normal(0,5),p)
    sigma ~ Uniform(0,1)

    # This works, but is slower (and not so neat)
    # for i in 1:n
    #     mu = a + c[1]*x[i,1] + c[2]*x[i,2] + c[3]*x[i,3] + c[4]*x[i,4]
    #     y[i] ~ Normal(mu,sigma)
    # end

    # Better (and neater)
    mu = a .+ x*params
    y ~ MvNormal(mu, sigma)

    return y

end

model = iris_model(x,y)

# println("Prior:")
# chns = sample(model, Prior(), 10_000)
# display(chns)

println("\nModel:")
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
# chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)
# chns = sample(model, HMC(0.1,6), 1_000)
chns = sample(model, NUTS(), 1_000)

display(chns)

println("\nparameter: credible_intervals")
for v in ["a", "params[1]", "params[2]", "params[3]", "params[4]","sigma"]
  credible_interval(chns,v,0.95)
end

# display(plot(chns))
println("\nPredictions:")

# mpred = iris_model(x[50:51,:],missing)
mpred = iris_model(x,missing)
pred = predict(mpred,chns)
display(pred)

cats = Dict(1=>1, 2=>2, 3=>3)
confusion, bads = confusion_matrix(y,pred,cats,false)
println("Confusion matrix:")
display(confusion)
println("\nBad classfied:")
display(bads)
println("num bad: $(length(bads))")
println("\nbads: credible_intervals")
for b in bads
  println("b:",b)
  ix = b[1]
  println("x[$ix]: $(x[ix,:])  true cat (y[$ix]): $(y[ix])  (pred mean: ", mean(pred["y[$ix]"]))
  credible_interval(pred,"y[$(b[1])]",0.3)

end

println("\nMean value for each species: ", 1:3)
display(1:3 .|> s -> (s,mean(x[y.==s,:], dims=1)))
println("\n")

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using JSON
using TOML

Base.:/(a::String, b::String) = joinpath(a, b)

function concept_exercises()
    base_path = "exercises" / "concept"
    collect(ex for ex in readdir(base_path) if isdir(base_path / ex))
end

function create_meta_configs()
    base_path = "exercises" / "concept"
    exercises = concept_exercises()

    Threads.@threads for exercise in exercises
        meta_path = base_path / exercise / ".meta"
        human_cfg = TOML.parsefile(meta_path / "config.toml")

        cfg = Dict(
            "language_versions" => "≥1.0",
            "files" => Dict(
                "test" => ["runtests.jl"],
                "exemplar" => [".meta/exemplar.jl"],
                "solution" => [isnothing(get(human_cfg, "solution", nothing)) ? "$exercise.jl" : human_cfg["solution"]],
            ),
            "authors" => [
                Dict("github_username" => x["github"], "exercism_username" => x["exercism"]) for x in human_cfg["authors"]
            ],
            "contributors" => [
                Dict("github_username" => x["github"], "exercism_username" => x["exercism"]) for x in human_cfg["contributors"]
            ],
            "forked_from" => get(human_cfg, "forked_from", nothing),
        )
        open(meta_path / "config.json", "w") do f
            JSON.print(f, cfg, 2)
        end
    end
end

function main()
    create_meta_configs()
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

max_stars_repo_name: tripitakit/julia

max_issues_repo_name: tripitakit/julia

max_forks_repo_name: tripitakit/julia


text:
using PyPlot
using LinearAlgebra, Random

"""
Gradient ascent to train a 2-layers ReLU neural net for the soft-min loss
INPUT: X (training input), Y (training output), m (nb neurons), both: training both layers or just the output
OUTPUT: Ws (training trajectory)
"""
function twonet(X, Y, m, stepsize, niter; both=true)
    (n,d) = size(X)
    W_init = randn(m, d+1)
    if !both
        W_init[:,end] .= 0
    end

    W     = copy(W_init)
    Ws    = zeros(m, d+1, niter) # store optimization path
    loss  = zeros(niter)
    margins = zeros(niter)
    betas = zeros(niter)

    for iter = 1:niter
        Ws[:,:,iter] = W
        act  =  max.( W[:,1:end-1] * X', 0.0) # (size m × n)
        out  =  (1/m) * sum( W[:,end] .* act , dims=1) # (size 1 × n)
        perf = Y .* out[:]
        margin = minimum(perf)
        temp = exp.(margin .- perf) # stabilization
        gradR = temp .* Y ./ sum(temp)' # size n
        grad_w1 = (W[:,end] .* float.(act .> 0) * ( X .* gradR  ))  # (size m × d)
        grad_w2 = act * gradR  # size m

        if both
            grad = cat(grad_w1, grad_w2, dims=2) # size (m × d+1)
            betas[iter] = sum(W.^2)/m
            loss[iter] = margin - log(sum(exp.(margin .- perf))/n)
            margins[iter] = margin/betas[iter]
            W = W + stepsize * grad/(sqrt(iter+1))

         else
            grad = cat(zeros(m,d), grad_w2, dims=2) # size (m × d+1)
            betas[iter] = maximum([1,sqrt(sum(W[:,end].^2)/m)])
            loss[iter] = margin - log(sum(exp.(margin .- perf))/n)
            margins[iter] = margin/(sqrt(sum(W[:,end].^2))/m)
            W = W + betas[iter] * stepsize * grad /(sqrt(iter+1))
        end
    end
    Ws, loss, margins, betas
end

"Coordinates of the 2d cluster centers, p is k^2 the number of clusters"
function cluster_center(p,k)
    p1 = mod.(p .- 1,k) .+ 1
    p2 = div.(p .- 1,k) .+ 1
    Delta = 1/(3k-1)
    x1 =  Delta*(1 .+ 3*(p1 .- 1)) .- 1/2
    x2 =  Delta*(1 .+ 3*(p2 .- 1)) .- 1/2
    return x1,x2
end

"""
Plot the classifier for a test case, comparing training both or output layer
"""
function illustration(k, n, m; stepsize= 0.5, niter=100000, name="decision")
# data distribution
sd = 0 # number of spurious dimensions
Delta = 1/(3k-1) # interclass distance
A = ones(k^2) # cluster affectation
A[randperm(k^2)[1:div(k^2,2)]] .= -1

# sample from it
P = rand(1:k^2,n) # cluster label
T = 2π*rand(n)  # shift angle
R = Delta*rand(n) # shift magnitude
X = cat(ones(n), cluster_center(P,k)[1] .+ R .* cos.(T),cluster_center(P,k)[2] + R .* sin.(T), (rand(n,sd) .- 1/2), dims=2)
Y = A[P]

# train neural network
Ws1, loss1, margins1, betas1 = twonet(X, Y, m, stepsize, niter; both=true)
Ws2, loss2, margins2, betas2 = twonet(X, Y, m, stepsize, niter; both=false)

# plots
X1 = X[(Y .== 1),:]
X2 = X[(Y .== -1),:]

figure(figsize=[2.5,2.5])
    plot(X1[:,2],X1[:,3],"+r")
    plot(X2[:,2],X2[:,3],"_b")
    plot(cluster_center(1:k^2,k)[1],cluster_center(1:k^2,k)[2],"ok")
    axis("equal");#axis("off")
    xticks([], []); yticks([], [])
    savefig(name * "setting.pdf",bbox_inches="tight")

figure(figsize=[2.5,2.5])
    f1(x1,x2,t) = (1/m) * sum( Ws1[:,end,t] .* max.( Ws1[:,1:3,t] * [1;x1;x2], 0.0)) # (size 1 × n)
    xs = -0.8:0.01:0.8
    tab1 = [f1(xs[i],xs[j],size(Ws1,3)) for i=1:length(xs), j=1:length(xs)]
    pcolormesh(xs', xs, tanh.(1000*tab1'),cmap="coolwarm",shading="gouraud",vmin=-1.0,vmax=1.0,edgecolor="face")
    xs = -0.8:0.005:0.8
    tab1 = [f1(xs[i],xs[j],size(Ws1,3)) for i=1:length(xs), j=1:length(xs)]
    contour(xs', xs, tanh.(1000*tab1'),levels =0, colors="k",antialiased = true,linewidths=2)
    plot(X1[:,2],X1[:,3],"+k")
    plot(X2[:,2],X2[:,3],"_k")
    axis("equal");axis("off");
    savefig(name * "_both.pdf",bbox_inches="tight")

 figure(figsize=[2.5,2.5])
    f2(x1,x2,t) = (1/m) * sum( Ws2[:,end,t] .* max.( Ws2[:,1:3,t] * [1;x1;x2], 0.0)) # (size 1 × n)
    xs = -0.8:0.01:0.8
    tab2 = [f2(xs[i],xs[j],size(Ws2,3)) for i=1:length(xs), j=1:length(xs)]
    pcolormesh(xs', xs, tanh.(1000*tab2'),cmap="coolwarm",shading="gouraud",vmin=-1.0,vmax=1.0,edgecolor="face")
    xs = -0.8:0.005:0.8
    tab2 = [f2(xs[i],xs[j],size(Ws2,3)) for i=1:length(xs), j=1:length(xs)]
    contour(xs', xs, tanh.(1000*tab2'),levels =0, colors="k",antialiased = true,linewidths=2)
    plot(X1[:,2],X1[:,3],"+k")
    plot(X2[:,2],X2[:,3],"_k")
    axis("equal");axis("off");
    savefig(name *"_output.pdf",bbox_inches="tight")
end

"""
Compute margin and generalisation error.
Run the experiment with the given list of d (input dimension), n (input points), m (width)
On n_distrib random input distributions and n_repeat random initialization.
"""
function experiment_2NN(n_distrib, ds, ns, n_repeat, ms; k = 3, niter = 5*10^5, stepsize = 1, n_test=4000)

    L1s = zeros(n_distrib, length(ds), length(ns), n_repeat, length(ms))
    L2s = zeros(n_distrib, length(ds), length(ns), n_repeat, length(ms))
    m1s = zeros(n_distrib, length(ds), length(ns), n_repeat, length(ms))
    m2s = zeros(n_distrib, length(ds), length(ns), n_repeat, length(ms))

    p = Progress(length(L1s))
    for i=1:n_distrib
        # define the distribution
        Delta = 1/(3k-1) # lower bound on the interclass distance
        A = ones(k^2)
        A[randperm(k^2)[1:div(k^2,2)]] .= -1 # cluster affectation

        for ii = 1: length(ds)
            d = ds[ii]
            # sample from the distribution (we then take a subset of it)
            nmax = maximum(ns)
            P = rand(1:k^2,nmax) # cluster label
            T = 2π*rand(nmax)  # shift angle
            R = Delta*rand(nmax) # shift magnitude
            sd = d-3 # number of spurious dimensions
            X = cat(ones(nmax),
                      cluster_center(P,k)[1] .+ R .* cos.(T),
                      cluster_center(P,k)[2] + R .* sin.(T),
                      rand(nmax,sd) .- 1/2,
                      dims=2)
            Y = A[P]
            for iii = 1:length(ns)
                n = ns[iii]
                for iv=1:n_repeat # same eexperiment but with another random init
                    for v=1:length(ms)
                        m = ms[v]

                        # train neural network
                        Ws1, loss1, margins1, betas1 = twonet(X[1:n,:], Y[1:n], m, stepsize, niter; both=true)
                        Ws2, loss2, margins2, betas2 = twonet(X[1:n,:], Y[1:n], m, stepsize, niter; both=false)

                        # test set
                        T_test = 2π*rand(n_test)  # shift angle
                        R_test = Delta*rand(n_test) # shift magnitude
                        P_test = rand(1:k^2,n_test) # cluster label
                        X_test = cat(ones(n_test), cluster_center(P_test,k)[1] .+ R_test .* cos.(T_test),
                        cluster_center(P_test,k)[2] + R_test .* sin.(T_test), (rand(n_test,sd) .- 1/2), dims=2)
                        Y_test = A[P_test]

                        # test error
                        preds1 = (1/m) * sum(Ws1[:,end,end] .* max.( Ws1[:,1:end-1,end] * X_test', 0.0), dims = 1)
                        L1 = sum((Y_test .* preds1') .< 0)/n_test # probability of error
                        preds2 = (1/m) * sum(Ws2[:,end,end] .* max.( Ws2[:,1:end-1,end] * X_test', 0.0), dims = 1)
                        L2 = sum((Y_test .* preds2') .< 0)/n_test # probability of error
                        L1s[i,ii,iii,iv,v] = L1
                        L2s[i,ii,iii,iv,v] = L2
                        m1s[i,ii,iii,iv,v] = maximum(margins1)
                        m2s[i,ii,iii,iv,v] = maximum(margins2)
                        GC.gc() # garbage collection
                        next!(p)
                    end
                end
            end
        end
    end
    return L1s, L2s, m1s, m2s
end


text:
using Plots
using ImageMagick
using FileIO
using Images

using DataStructures

import Plots: default
default(show=false, reuse=true)

ENV["GKSwstype"]="nul"

scale = 2
fnt = font("Computer Modern", 10*scale)
default(titlefont=fnt, guidefont=fnt, tickfont=fnt, legendfont=fnt)

function convert_to_float(a::RGB)
    a = float(a)
    return (a.r + a.g + a.b)/3.0
end

function maze_to_canvas(filename)
    img = load(filename)
    return convert_to_float.(img)
end

function inbounds(canvas_size, loc)
    if minimum(Tuple(loc)) < 1
        return false
    elseif loc[2] > canvas_size[2]
        return false
    elseif loc[1] > canvas_size[1]
        return false
    else
        return true
    end
end

function create_circle!(canvas, loc::CartesianIndex, radius, threshold,
                        old_val)
    for i = 1:size(canvas)[1]
        for j = 1:size(canvas)[2]
            if ((radius-threshold)^2
                < (i-loc[1])^2+(j-loc[2])^2
                < (radius+threshold)^2)
                canvas[i,j] = old_val
            end
        end
    end
end

function color!(canvas, loc::CartesianIndex, old_val, new_val)
    # bounds check
    if !inbounds(size(canvas), loc)
        return
    end

    if (old_val == new_val)
        return
    elseif (canvas[loc] != old_val)
        return
    else
        canvas[loc] = new_val
    end
end

function find_neighbors(canvas, loc::CartesianIndex, old_val, new_val)

    neighbors =  []
    possible_neighbors = [loc + CartesianIndex(0, 1),
                          loc + CartesianIndex(1, 0),
                          loc + CartesianIndex(0, -1),
                          loc + CartesianIndex(-1, 0)]
    for possible_neighbor in possible_neighbors
        if inbounds(size(canvas), possible_neighbor) &&
           canvas[possible_neighbor] == old_val
            push!(neighbors, possible_neighbor)
        end
    end

    return neighbors
end

function stack_fill!(canvas, loc::CartesianIndex, old_val, new_val;
                     output=false)
    if new_val == old_val
        return
    end

    s = Stack{CartesianIndex}()
    #s = CartesianIndex[]
    push!(s, loc)

    id = 0

    while length(s) > 0
        temp_loc = pop!(s)
        if canvas[temp_loc] == old_val
            color!(canvas, temp_loc, old_val, new_val)
            if output
                plt = heatmap(canvas, color=:coolwarm, ratio=:equal,
                              size=(1000,1000))
                savefig("out"*lpad(string(id[1]),5,string(0))*".png")
                id += 1
                println("Outputting: ", id)
            end
            possible_neighbors = find_neighbors(canvas, temp_loc,
                                                old_val, new_val)
            for neighbor in possible_neighbors
                push!(s,neighbor)
            end
        end

    end
end

function queue_fill!(canvas, loc::CartesianIndex, old_val, new_val;
                     output=false)
    if new_val == old_val
        return
    end

    q = Queue{CartesianIndex}()
    enqueue!(q, loc)
    color!(canvas, loc, old_val, new_val)

    id = 0

    while length(q) > 0
        temp_loc = dequeue!(q)
        if output
            plt = heatmap(canvas, color=:coolwarm, ratio=:equal,
                          size=(500,500))
            savefig("out"*lpad(string(id[1]),5,string(0))*".png")
            id += 1
            println("Outputting: ", id)
        end

        possible_neighbors = find_neighbors(canvas, temp_loc, old_val, new_val)
        for neighbor in possible_neighbors
            color!(canvas, neighbor, old_val, new_val)
            enqueue!(q,neighbor)
        end

    end
end

function maze_fill!(canvas, loc::CartesianIndex, old_val, new_val;
                     output=false, end_loc = CartesianIndex(-1,-1))
    if new_val == old_val
        return
    end

    q = Queue{CartesianIndex}()
    enqueue!(q, loc)
    color!(canvas, loc, old_val, new_val)

    id = 0

    stop_solve = false
    while !stop_solve
        temp_loc = dequeue!(q)
        if output
            plt = heatmap(canvas, ratio=:equal,
                          size=(500,500), legend=:bottomright)
            scatter!(plt, Tuple(loc), label="Start")
            scatter!(plt, Tuple(end_loc), label="Finish")
            savefig("out"*lpad(string(id[1]),5,string(0))*".png")
            id += 1
            println("Outputting: ", id)
        end

        possible_neighbors = find_neighbors(canvas, temp_loc, old_val, new_val)
        for neighbor in possible_neighbors
            color!(canvas, neighbor, old_val, new_val)
            enqueue!(q,neighbor)
        end

        if (length(q) == 0 || temp_loc == end_loc)
            stop_solve = true
        end

    end
end

function recursive_fill!(canvas, loc::CartesianIndex, old_val, new_val;
                         output=false, id=id)

    if (old_val == new_val)
        return
    elseif (canvas[loc] != old_val)
        return
    else
        canvas[loc] = new_val
    end

    if output
        plt = heatmap(canvas, color=:coolwarm, ratio=:equal, size=(1000,1000))
        savefig("out"*lpad(string(id[1]),5,string(0))*".png")
        id[1] += 1
        println("Outputting: ", id[1])
    end

    possible_neighbors = find_neighbors(canvas, loc, old_val, new_val)
    for possible_neighbor in possible_neighbors
        recursive_fill!(canvas, possible_neighbor, old_val, new_val;
                        output=output,id=id)
    end
end

function init_canvas()
    canvas = zeros(30,30)
    loc = CartesianIndex(15,15)
    create_circle!(canvas, loc, 12, 2, 0.75)
    return canvas
end

function create_grid(dims)
    tuple_array = Array{Tuple{Int, Int}}(undef, dims[1]*dims[2])
    grid = Array{Float64}(undef, dims)
    for i = 1:dims[1]
        for j = 1:dims[2]
            tuple_array[(i-1)*dims[2]+j] = (i,j)
            grid[i,j] = i+j
        end
    end

    lines = Array{Tuple{Int, Int}}(undef, 10,2)
    for i = 1:5
        lines[i,:]=[tuple_array[(i-1)*5+1], tuple_array[i*5]]
        lines[i+5,:]=[tuple_array[i], tuple_array[20+i]]
    end

    println(lines)

    plt = heatmap(grid, color=:coolwarm, size=(600,500), aspect_ratio=:equal,
                 xlims=(0.5,5.5), ylims=(0.5,5.5))
    for i = 1:size(lines)[1]
        plot!(plt, lines[i,:], label="", color=:black, linewidth=3)
    end
    #plot!(plt, lines, label="")
    scatter!(plt, tuple_array, label = "", markersize=15)

    savefig("grid.png")

end

function main()

    create_grid((5,5))
#=
    canvas = init_canvas()
    loc2 = CartesianIndex(15,15)
    #recursive_fill!(canvas, loc2, 0.0, 0.5; output=true, id = [0])
    #queue_fill!(canvas, loc2, 0.0, 0.5; output=true)
    #stack_fill!(canvas, loc2, 0.0, 0.5; output=true)

    plt = heatmap(canvas, color=:coolwarm, size=(600,500), aspect_ratio=:equal,
                  xlims=(1,30), ylims=(1,30))
    savefig("simple_circle.png")

    # For maze
    maze = maze_to_canvas("maze2.png")
    maze_start = CartesianIndex(50,50)
    maze_end = CartesianIndex(96,96)
    maze_fill!(maze, maze_start, 1.0, 0.5; end_loc = maze_end, output=true)
    heatmap(maze)
=#

end

main()


text:
module NestedTuples

include("utils.jl")
include("with.jl")
include("lazymerge.jl")
include("rmap.jl")
include("fold.jl")
include("placeholder.jl")
include("typelevel.jl")
include("lenses.jl")
# include("nested.jl")
include("keysort.jl")
include("exprify.jl")
include("leaves.jl")
include("leaves2.jl")
# include("arrays.jl")
# include("summarize.jl")
end


text:
# Autogenerated wrapper script for Xorg_xcb_util_wm_jll for aarch64-linux-gnu
export libxcb_ewmh, libxcb_icccm

using Xorg_xcb_util_jll
JLLWrappers.@generate_wrapper_header("Xorg_xcb_util_wm")
JLLWrappers.@declare_library_product(libxcb_ewmh, "libxcb-ewmh.so.2")
JLLWrappers.@declare_library_product(libxcb_icccm, "libxcb-icccm.so.4")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_xcb_util_jll)
    JLLWrappers.@init_library_product(
        libxcb_ewmh,
        "lib/libxcb-ewmh.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libxcb_icccm,
        "lib/libxcb-icccm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Firefly

export JointModel

include("joint.jl")
include("distributions.jl")

end