text:
include(joinpath("Thesis model", "glu_emp_selfrep_model.jl"))
import .Glu_EMP_Selfrep_Model

res_emp_selfrep = []
for var in [0.0001:0.0005:0.005;]
    push!(res_emp_selfrep, Glu_EMP_Selfrep_Model.glu_emp_selfrep_model(var))
end

################################################################

using Colors, ColorSchemes
using CairoMakie, FileIO
using CSV, DataFrames
using Statistics

imgpath = joinpath("Thesis model", "Figures", "EMP_Model")
resultspath = joinpath("Thesis model", "Results", "EMP_Model")

#CSV.write(joinpath(resultspath, "Results_EMP_Selfrep.csv"), res_emp[1])

x = [results_emp_selfrep["glc_ext"] for results_emp_selfrep in res_emp_selfrep]
y = [results_emp_selfrep["mu"] for results_emp_selfrep in res_emp_selfrep]

fig = Figure()
ax = Axis(fig[1, 1])
scatter_glc = scatter!(ax, x, y)
line_glc = lines!(ax, x, y)

ax.xlabel = "External glucose concetration [mM]"
ax.ylabel = "Biomass [1/h]"

fig

FileIO.save(joinpath(imgpath, "Glu_EMP_Selfrep_Growth.pdf"), fig)

## ##

x = [results_emp_selfrep["mu"] for results_emp_selfrep in res_emp_selfrep]
a = [results_emp_selfrep["g6p"] for results_emp_selfrep in res_emp_selfrep]
b = [results_emp_selfrep["glu"] for results_emp_selfrep in res_emp_selfrep]

fig = Figure()
ax = Axis(fig[1, 1])
scatter1 = scatter!(ax, x, a)
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b)
line2 = lines!(ax, x, b)

ax.xlabel = "Biomass [1/h]"
ax.ylabel = "Metabolite concentration"

fig[1, 1] = Legend(
    fig,
    [[scatter1, line1], [scatter2, line2]],
    ["Glucose-6-Phosphate", "Glutamate"],
    tellheight = false,
    tellwidth = false,
    halign = :left,
    valign = :bottom,
    labelsize = 14,
)
fig

FileIO.save(joinpath(imgpath, "Glu_EMP_G6P+GLU_Selfrep.pdf"), fig)

## Enzyme concentration band plot ##

x = [results_emp_selfrep["mu"] for results_emp_selfrep in res_emp_selfrep]
c = [results_emp_selfrep["pts"] for results_emp_selfrep in res_emp_selfrep]
d = [results_emp_selfrep["emp"] for results_emp_selfrep in res_emp_selfrep]
e = [results_emp_selfrep["pyk"] for results_emp_selfrep in res_emp_selfrep]
f = [results_emp_selfrep["ldh"] for results_emp_selfrep in res_emp_selfrep]
g = [results_emp_selfrep["burn"] for results_emp_selfrep in res_emp_selfrep]
h = [results_emp_selfrep["lp"] for results_emp_selfrep in res_emp_selfrep]

fracs = [c, d, e, f, g, h]
fraclabels = [
    "PTS",
    "EMP",
    "Pyruvate kinase",
    "Lactate dehydrogenase",
    "ATP burning enzyme",
    "Lactate permease",
]
function plot_proteome(ax, x, fracs, fraclabels)
    d = cumsum(hcat(fracs...)'; dims = 1)
    d = [zeros(size(d, 2))'; d ./ d[end, :]']
    for i = 1:length(fraclabels)
        band!(ax, x, d[i, :], d[i+1, :], label = fraclabels[i])
    end
end

f = Figure();
ax = Axis(f[1, 1]);
plot_proteome(ax, 1:10, fracs, fraclabels)
ax.xlabel = "Biomass [1/h]"
ax.ylabel = "Proteome fraction"
f[1, 2] = Legend(f, ax, "Protein", unique = true, framevisible = false)
f

FileIO.save(joinpath(imgpath, "Glu_EMP_Selfrep_Enzymeband.pdf"), f)

## Rates ##

x = [results_emp_selfrep["mu"] for results_em_selfrep in res_emp_selfrep]
i = [results_emp_selfrep["v_pts"] for results_em_selfrep in res_emp_selfrep]
j = [results_emp_selfrep["v_pyk"] for results_em_selfrep in res_emp_selfrep]
k = [results_emp_selfrep["v_ldh"] for results_em_selfrep in res_emp_selfrep]
l = [results_emp_selfrep["v_burn"] for results_em_selfrep in res_emp_selfrep]
m = [results_emp_selfrep["v_lp"] for results_em_selfrep in res_emp_selfrep]

fig = Figure()
ax = Axis(fig[1, 1])
scatter1 = scatter!(ax, x, i)
line1 = lines!(ax, x, i)
scatter2 = scatter!(ax, x, j)
line2 = lines!(ax, x, j)
scatter3 = scatter!(ax, x, k)
line3 = lines!(ax, x, k)
scatter4 = scatter!(ax, x, l)
line4 = lines!(ax, x, l)
scatter5 = scatter!(ax, x, m)
line5 = lines!(ax, x, m)

ax.xlabel = "Biomass [1/h]"
ax.ylabel = "Fluxes [mmol/gDW/h]"

fig[1, 1] = Legend(
    fig,
    [
        [scatter1, line1],
        [scatter2, line2],
        [scatter3, line3],
        [scatter4, line4],
        [scatter5, line5],
    ],
    ["v_pts", "v_pyk", "v_ldh", "v_burn", "v_lp"],
    tellheight = false,
    tellwidth = false,
    halign = :left,
    valign = :bottom,
    labelsize = 14,
)
fig

FileIO.save(joinpath(imgpath, "Glu_EMP_Selfrep_Rates.pdf"), fig)


text:
#reload("CvGeom")
using CvGeom
#reload("CvGeom")
using Base.Test

# write your own tests here

# some test data for homogeneous points
# some finite points first
a = [1, 2];       a1 = [1, 2, 1];       a2 = a1*2.0;         a3 = [.5, 1, .5]
b = [-3, 2];      b1 = [-3, 2, 1];      b2 = b1*-0.5;        b3 = [1, -2/3, -1/3]
c = [1/3, 2/3];   c1 = [1/3, 2/3, 1];   c2 = [1,2,3];        c3 = c1

# a point at infinite
d1 = [-3, 2, 0];      d2 = [1.5, -1, 0];   d3 = [1, -2/3, 0]

ABC = [a b c]
ABC1 = [a1 b1 c1]
ABC2 = [a2 b2 c2]
ABC3 = [a3 b3 c3]
D = [d1 d2 d3]

@test hom(a) == a1
@test hom(b) == b1
@test hom(c) == c1
@test hom(ABC) == ABC1

aa = Set({a1, a2, a3})
bb = Set({b1, b2, b3})
cc = Set({c1, c2, c3})
dd = Set({d1, d2, d3})

function test_ihom(u::Vector, uu::Set)
  for ui = uu
    @test ihom(ui) == u
  end
end

test_ihom(a, aa)
test_ihom(b, bb)
test_ihom(c, cc)

@test ihom(ABC1) == ABC
@test ihom(ABC2) == ABC
@test ihom(ABC3) == ABC

function test_homNorm(u1::Vector, uu::Set)
  for ui = uu
    @test homNorm(ui) == u1
  end
end

test_homNorm(a1, aa)
test_homNorm(b1, bb)
test_homNorm(c1, cc)

function test_homNormRobust(u3::Vector, uu::Set)
  for ui = uu
    @test homNormRobust(ui) == u3
  end
end

test_homNormRobust(a3, aa)
test_homNormRobust(b3, bb)
test_homNormRobust(c3, cc)

function isHomEq2(a::Vector, b::Vector)
  if size(a) != size(b)
    return false
  end

  return homNormRobust(a) == homNormRobust(b)

end

function test_isHomEq(uu::Set, vv::Set)
  for ui = uu
    for uj = uu
      @test isHomEq2(ui, uj)
    end
    for vj = vv
      @test !isHomEq2(ui, vj)
    end
  end
end

test_isHomEq(aa, bb)
test_isHomEq(bb, cc)
test_isHomEq(cc, dd)
test_isHomEq(dd, aa)

p1 = [1,7.5]
p2 = [13,102]
p3 = [1,1]
p4 = [-3,24]
p5 = [11,-7]

points = [p1 p2 p3 p4 p5]
c,A,d = CvGeom.fitConic(points)

A*c

function conic(c::Vector)
  @assert length(c) == 6
  a,b,c,d,e,f = c
  C = [ a  b/2 d/2
       b/2  c  e/2
       d/2 e/2  f ]
  return C
end

C = conic(c)
hom(p1)' * C * hom(p1)

#@test CvGeom.isPointOnLine

ext: jl


text:
slope = open("./inputs/day3.txt") do f
	readlines(f) |> l -> split.(l, "")
end

function count_trees(Δx, Δy)
	x = 1
	y = 1

	trees = 0

	while y <= length(slope)
		if slope[y][x] == "#"
			trees += 1
		end

		x += Δx
		y += Δy

		if x > length(slope[1])
			x %= length(slope[1])
		end
	end

	return trees
end

part1() = count_trees(3, 1)

part2() = ([(1, 1), (3, 1), (5, 1), (7, 1), (1, 2)] .|> args -> count_trees(args...)) |> prod

println("day 3:\n\tpart 1: ", part1(), "\n\tpart 2: ", part2())


text:
const BlockTriPertToeplitz{T} = BlockMatrix{T,Tridiagonal{Matrix{T},Vcat{Matrix{T},1,Tuple{Vector{Matrix{T}},Fill{Matrix{T},1,Tuple{OneToInf{Int}}}}}},
                                        BlockSizes{2,NTuple{2,Vcat{Int,1,Tuple{Int,Vcat{Int,1,Tuple{Vector{Int},InfStepRange{Int,Int}}}}}}}}

function BlockTridiagonal(adjA::Adjoint{T,BlockTriPertToeplitz{T}}) where T
    A = parent(adjA)
    BlockTridiagonal(Matrix.(adjoint.(A.blocks.du)),
                     Matrix.(adjoint.(A.blocks.d)),
                     Matrix.(adjoint.(A.blocks.dl)))
end

for op in (:-, :+)
    @eval begin
        function $op(A::BlockTriPertToeplitz{T}, λ::UniformScaling) where T
            TV = promote_type(T,eltype(λ))
            BlockTridiagonal(Vcat(convert.(AbstractVector{Matrix{TV}}, A.blocks.dl.args)...),
                             Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.d, Ref(λ)).args)...),
                             Vcat(convert.(AbstractVector{Matrix{TV}}, A.blocks.du.args)...))
        end
        function $op(λ::UniformScaling, A::BlockTriPertToeplitz{V}) where V
            TV = promote_type(eltype(λ),V)
            BlockTridiagonal(Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.dl.args))...),
                             Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, Ref(λ), A.blocks.d).args)...),
                             Vcat(convert.(AbstractVector{Matrix{TV}}, broadcast($op, A.blocks.du.args))...))
        end
        $op(adjA::Adjoint{T,BlockTriPertToeplitz{T}}, λ::UniformScaling) where T = $op(BlockTridiagonal(adjA), λ)
        $op(λ::UniformScaling, adjA::Adjoint{T,BlockTriPertToeplitz{T}}) where T = $op(λ, BlockTridiagonal(adjA))
    end
end

*(a::AbstractVector, b::AbstractFill{<:Any,2,Tuple{OneTo{Int},OneToInf{Int}}}) = ApplyArray(*,a,b)

sizes_from_blocks(A::AbstractVector, ::Tuple{OneToInf{Int}}) = BlockSizes((Vcat(1, 1 .+ cumsum(length.(A))),))
function sizes_from_blocks(A::AbstractMatrix, ::Tuple{OneTo{Int}, OneToInf{Int}})
    @assert size(A,1) == 1
    R,C = vec(size.(A,1)[:,1]), vec(size.(A,2))
    BlockSizes((vcat(1, 1 .+ cumsum(R)), Vcat(1, 1 .+ cumsum(C))))
end

function sizes_from_blocks(A::Diagonal, ::NTuple{2,OneToInf{Int}})
    sz = size.(A.diag, 1), size.(A.diag,2)
    BlockSizes(Vcat.(1,(c -> 1 .+ c).(cumsum.(sz))))
end

function sizes_from_blocks(A::Tridiagonal, ::NTuple{2,OneToInf{Int}})
    sz = size.(A.d, 1), size.(A.d,2)
    BlockSizes(Vcat.(1,(c -> 1 .+ c).(cumsum.(sz))))
end

_find_block(cs::Number, i::Integer) = i ≤ cs ? 1 : 0
function _find_block(cs::Vcat, i::Integer)
    n = 0
    for a in cs.args
        i < first(a) && return n
        if i ≤ last(a)
            return _find_block(a, i) + n
        end
        n += length(a)
    end
    return 0
end

print_matrix_row(io::IO,
        X::AbstractBlockVecOrMat, A::Vector,
        i::Integer, cols::AbstractVector{<:Infinity}, sep::AbstractString) = nothing

print_matrix_row(io::IO,
        X::Union{AbstractTriangular{<:Any,<:AbstractBlockMatrix},
                 Symmetric{<:Any,<:AbstractBlockMatrix},
                 Hermitian{<:Any,<:AbstractBlockMatrix}}, A::Vector,
        i::Integer, cols::AbstractVector{<:Infinity}, sep::AbstractString) = nothing

function BlockSkylineSizes(A::BlockTriPertToeplitz, (l,u)::NTuple{2,Int})
    N = max(length(A.blocks.du.args[1])+1,length(A.blocks.d.args[1]),length(A.blocks.dl.args[1]))
    block_sizes = Vector{Int}(undef, N) # assume square
    block_starts = BandedMatrix{Int}(undef, (N+l,N),  (l,u))
    block_strides = Vector{Int}(undef, N)
    for J=1:N
        block_starts[max(1,J-u),J] = J == 1 ? 1 :
                            block_starts[max(1,J-1-u),J-1]+block_sizes[J-1]*block_strides[J-1]

        for K=max(1,J-u)+1:J+l
            block_starts[K,J] = block_starts[K-1,J]+size(A[Block(K-1,J)],1)
        end
        block_strides[J] = block_starts[J+l,J] + size(A[Block(J+l,J)],1) - block_starts[max(1,J-u),J]
        block_sizes[J] = size(A[Block(J,J)],2)
    end

    block_stride∞ = 0
    for K=max(1,N+1-u):N+1+l
        block_stride∞ += size(A[Block(K,N+1)],1)
    end
    block_size∞ = size(A[Block(N+1,N+1)],2)

    bs∞ = fill(block_starts[max(1,N-u),N]+block_strides[N]*size(A[Block(N,N)],2):block_stride∞*block_size∞:∞, l+u+1)
    for k=2:l+u+1
        bs∞[k] = bs∞[k-1] .+ size(A[Block(N+1-u+k-1,N+1)],1)
    end

    BlockSkylineSizes(blocksizes(A),
                        _BandedMatrix(Hcat(block_starts.data, Vcat(adjoint.(bs∞)...)), ∞, l, u),
                        Vcat(block_strides, Fill(block_stride∞,∞)),
                        Fill(l,∞),Fill(u,∞))
end

function BlockBandedMatrix(A::BlockTriPertToeplitz{T}, (l,u)::NTuple{2,Int}) where T
    data = T[]
    append!(data,vec(A[Block.(1:1+l),Block(1)]))
    N = max(length(A.blocks.du.args[1])+1,length(A.blocks.d.args[1]),length(A.blocks.dl.args[1]))
    for J=2:N
        append!(data, vec(A[Block.(max(1,J-u):J+l),Block(J)]))
    end
    tl = mortar(Fill(vec(A[Block.(max(1,N+1-u):N+1+l),Block(N+1)]),∞))

    B = _BlockSkylineMatrix(Vcat(data,tl), BlockSkylineSizes(A, (l,u)))
end

BlockBandedMatrix(A::BlockTriPertToeplitz) = BlockBandedMatrix(A, blockbandwidths(A))


text:
#md # # Plot recipes examples
#nb # # BAT.jl Plot recipes examples
#jl # # BAT.jl Plot recipes examples

using BAT

# ## Generate samples
# We generate some generic multimodal samples to be plotted afterwards
using Distributions, IntervalSets

##Model definition to generate samples
struct GaussianShellDensity<:AbstractDensity
    r::Vector{Float64}
    sigma::Vector{Float64}
end

#define likelihood
function BAT.density_logval(target::GaussianShellDensity, params::Any)
    loglikelihood::Float64 = 0.
    for i in 1:length(params)

        result = -0.5 * (params[i][1]-target.r[i])^2/target.sigma[i]^2
        l1 = result - log(√2π * target.sigma[i])

        result2 = -0.5 * (params[i][1]+target.r[i])^2/target.sigma[i]^2
        l2 = result2 - log(√2π * target.sigma[i])^2

        loglikelihood += log(exp(l1) + 2*exp(l2))
    end

    return loglikelihood
end

likelihood = GaussianShellDensity([5.0, 5., 3.], [2, 2.4, 1.5])

prior = NamedTupleDist(
    λ1 = Normal(-3, 4.5),
    λ2 = -30.0..30.0,
    λ3 = Normal(3, 3.5)
)

posterior = PosteriorDensity(likelihood, prior);

#generate samples
samples, stats, chains = bat_sample(posterior, (10^5, 4), MetropolisHastings())

# ## Set up plotting
#  We set up plotting by using the ["Plots"](https://github.com/JuliaPlots/Plots.jl) package with the pyplot backend.
using Plots

# The plotting recipes have been tested with `pyplot`, `gr` and `plotly` backends.
# All available plotting options work correctly with `pyplot`.
# Only the colored 2D contour plots are currently not correctly supported with `gr` or `plotly` backends.

pyplot()

# ## Examples for 1D plots
# Below, all available seriestypes and plotting features for 1D representations of samples or prior are shown
# ### default 1D plot style for samples:
plot(samples, 2)
# The default seriestype for samples is `:smallest_intervals`.
# By default, the local mode of the histogram is indicated as a dotted black line.

# ### default 1D plot style for prior:
plot(prior, 2)
# The default seriestype for sprior is `:stephist`.

# Samples can either be plotted by their index (as shown above) or by using the parameter names given in NamedTupleDist. This can be done by passing either the posterior
plot(posterior, samples, :λ2)
# or the value shape(s) of the prior:
using ValueShapes
parshapes = valshape(prior)
plot(parshapes, samples, :λ2)

# Prior can also be plotted by their index or by using the parameter names given in NamedTupleDist:
plot(prior, 1)
# or
plot(prior, :λ1)

# Plot prior + posterior
plot(posterior, samples, :λ1)
plot!(prior, :λ1)

# ### Further seriestypes available:
# #### central intervals
plot(samples, 2, seriestype = :central_intervals)

# ### Customizing smallest/central interval plots:
# The probability intervals to be highlighted can be specified using the `intervals` keyword. Their colors (in same order) need to be specified using the `colors` keyword argument.
plot(samples, 2, seriestype=:smallest_intervals, intervals=[0.5, 0.1, 0.3, 0.99], colors=[:grey, :red, :blue, :orange])

# #### histogram
plot(samples, 2, seriestype = :histogram)
# alias: `:hist == :histogram`

# #### step histogram
plot(samples, 2, seriestype = :stephist)

# ### Keyword arguments for [attributes supported by *Plots.jl*](https://docs.juliaplots.org/latest/attributes/#magic-arguments) can be passed:
plot(samples, 2, seriestype = :stephist, nbins=50, linecolor = :red, linewidth = 5, linealpha=0.3, xlim=(0,10))

# ### Plotting point estimators in 1D plots:
plot(samples, 2, globalmode=true, localmode=true, mean=true, std_dev=true)

# #### It is possible to customize the style of estimators by passing a dict with the respective attributes:
# By passing a boolean, the point estimators are plotted using their default styles shown above.
# The style of the point estimators can be modified by passing a dictionary specifying `linestyle`, `linecolor`, `linewidth` and `alpha` for *mean*, *globalmode* and *localmode*.
# The style of the standard deviation can be modified by specifying `fillcolor` and `fillalpha`.
plot(samples, 1, localmode=false, mean=Dict("linestyle" => :dot, "linecolor"=> :red, "linewidth"=>2, "alpha" => 0.7), std_dev=Dict("fillcolor" => :red, "fillalpha"=>0.2))

# ## Examples for 2D plots of the samples
# Below, all available seriestypes and plotting features for 2D representations of samples are shown
# ### default 2D plot style:
plot(samples, (1,2), nbins=200)
# The default 2D plotting style is a 3-color heatmap showing the smallest intervals containing 68.3, 95.5 and 99.7 perecent of the posterior probability.

# Samples can either be plotted by their index (as shown above) or by using the parameter names given in NamedTupleDist. This can be done by passing either the posterior
plot(posterior, samples, (:λ1, :λ2))
# or the value shape(s) of the prior:
parshapes = valshape(prior)
plot(parshapes, samples, (:λ1, :λ2))

# Prior can also be plotted by their index or by using the parameter names given in NamedTupleDist.
plot(prior, (1, 2))
# or
plot(prior, (:λ1, :λ2))

# Plot prior + posterior
plot(posterior, samples, (:λ1, :λ3))
plot!(prior, (:λ1, :λ3))

# ### Further seriestypes available:
# #### central intervals
plot(samples, (1,2), seriestype = :histogram, nbins=200)
# alias: `:histogram2d == :histogram`

# #### smallest intervals as colored contour lines:
# (currently only correctly supported with `pyplot()` backend)
plot(samples, (1,2), seriestype=:smallest_intervals_contour, bins=50)

# #### smallest intervals as filled contours:
# (currently only correctly supported with `pyplot()` backend)
plot(samples, (1,2), seriestype=:smallest_intervals_contourf, bins=50)

# ### Customizing smallest interval plots:
# The probability intervals to be highlighted can be specified using the `intervals` keyword. Their colors (in same order) need to be specified using the `colors` keyword argument.
plot(samples, (1,2), seriestype=:smallest_intervals, nbins=200, intervals=[0.7, 0.2], colors=[:blue, :red])

# #### scatter
# for large number of samples, this plotting style needs a lot of time and resources
#plot(samples, (1,2), seriestype = :scatter)

# ### Plotting point estimators in 2D plots:
plot(samples, (1,2), seriestype=:smallest_intervals, nbins=200, mean=true, std_dev=true, localmode=true, globalmode=true)

# #### It is possible to customize the style of estimators by passing a dict with the respective attributes:
# By passing a boolean, the point estimators are plotted using their default styles shown above.
# The style of the point estimators can be modified by passing a dictionary specifying `markershape`, `markercolor`, `markersize`, `markeralpha`, `markerstrokecolor`, `markerstrokestyle`, `markerstrokewidth` and `markerstrokealpha` for *mean*, *globalmode* and *localmode*.
# If `std_dev==true`, the standard deviation of the mean value will be displayed as x- and y-errorbars.
plot(samples, (1,2), seriestype=:smallest_intervals, nbins=200, localmode=Dict("markershape"=> :diamond, "markeralpha"=>1, "markercolor"=>:red, "markersize"=>5))

# ### Marginal plots
plot(samples, (1,2), seriestype = :marginal)

# ### Customizing marginal plots:
# The marginal plots can be modified by passing dictionaries to the keyword arguments `upper`, `right` and `diagonal`.
# The dictionaries for `upper` and `right` can contain the 1D seriestypes and plot options shown above.
# The dictionary for `diagonal` can use the 2D seriestypes and plot options shown above.
# Nested dictonaries are possible (e.g. for modifying point estimators)
plot(samples, (1,2), seriestype = :marginal, diagonal = Dict("seriestype"=>:histogram, "mean"=>Dict("markercolor"=>:green)), upper=Dict("mean" => true, "seriestype" => :smallest_intervals, "colors"=>[:blue, :grey, :orange]), right=Dict("seriestype" => :stephist, "mean"=>true))

# ### Overview plot
plot(samples)
#
plot(prior)
#
plot(samples)
plot!(prior)

# By default the 1D and 2D plots for the first 5 parameters are shown.
# The keyword argument `params` allows to specify the considered parameters.
plot(samples, params=[1, 3])

# ### Customizing overview plots:
# The overview plots can be modified by passing dictionaries to the keyword arguments `upper`, `lower` and `diagonal`.
# The dictionaries for `upper` and `lower` can contain the 2D seriestypes and plot options shown above.
# The dictionary for `diagonal` can use the 1D seriestypes and plot options shown above.
# Nested dictonaries are possible (e.g. for modifying point estimators)

plot(samples, mean=true, globalmode=true, legend=true, diagonal=Dict("seriestype"=>:stephist, "mean"=>Dict("linecolor" => :green, "linewidth" => 8)), lower = Dict("mean" => false, "colors"=>[:orange, :green, :grey]))

# ## Plots for MCMC diagnostics
# Plots histograms of the samples, the trace, a kernel density estimate and the autocorrelation function for each parameter per chain.
diagnostics = BAT.MCMCDiagnostics(samples, chains)
plot(diagnostics, params=[1])

# ### Customizing diagnostics plots:
plot(diagnostics,
    params=[1, 2],
    chains=[1, 2],
    diagnostics = [:histogram, :kde, :trace, :acf],
    histogram = Dict("seriestype" => :smallest_intervals, "legend" => :false),
    trace = Dict("linecolor" => :red),
    acf = Dict("lags" => collect(1:20), "title"=> "Autocorrelation"),
    description = true)

# ### available keyword arguments:
# * `params` - list of parameters to be plotted
# * `chains` - list of chains to be plotted
# * `diagnostics` - list of MCMC diagnostics to be plotted
# * `:histogram` - 1D histograms of samples
# * `:kde` - Kernel density estimate (using [*KernelDensity.jl*](https://github.com/JuliaStats/KernelDensity.jl))
# * `:trace` - Trace plot
# * `:acf` - Autocorrelation function (using [*StatsBase.autocor*](http://juliastats.github.io/StatsBase.jl/stable/signalcorr/#StatsBase.autocor))
# * `description::Bool = true` - show description (current chain, parameter, number of samples) as first column of plots
# * `histogram::Dict` - options for histogram plots (supports all arguments for 1D plots for samples)
# * `kde::Dict` - options for kde plots
# * `trace::Dict` - options for trace plots
# * `acf::Dict` - options for acf plots

# ### special options arguments for `:kde` (see [*KernelDensity.jl*](https://github.com/JuliaStats/KernelDensity.jl))
# * `npoints`: number of interpolation points to use (default: npoints = 2048)
# * `boundary`: lower and upper limits of the kde as a tuple
# * `kernel`: the distributional family from [*Distributions.jl*](https://github.com/JuliaStats/Distributions.jl) to use as the kernel (default = Distributions.Normal)
# * `bandwidth`: bandwidth of the kernel

# ### special keyword arguments for `:acf` (see [*StatsBase.autocor*](http://juliastats.github.io/StatsBase.jl/stable/signalcorr/#StatsBase.autocor))
# * `lags` - list of lags to be considered for ACF plots
# * `demean` - denotes whether the mean should be subtracted before computing the ACF


text:
using Printf, DataFrames, CSV, LinearAlgebra, Statistics, Plots, Plots.PlotMeasures, Dates
# pyplot()

# abbreviate the name of a state to 2 capital letters
function abbreviate_state(state)
    d = Dict(
        "Alabama" => "AL",
        "Alaska" => "AK",
        "Arizona" => "AZ",
        "Arkansas" => "AR",
        "California" => "CA",
        "Colorado" => "CO",
        "Connecticut" => "CT",
        "Delaware" => "DE",
        "District of Columbia" => "DC",
        "Florida" => "FL",
        "Georgia" => "GA",
        "Hawaii" => "HI",
        "Idaho" => "ID",
        "Illinois" => "IL",
        "Indiana" => "IN",
        "Iowa" => "IA",
        "Kansas" => "KS",
        "Kentucky" => "KY",
        "Louisiana" => "LA",
        "Maine" => "ME",
        "Maryland" => "MD",
        "Massachusetts" => "MA",
        "Michigan" => "MI",
        "Minnesota" => "MN",
        "Mississippi" => "MS",
        "Missouri" => "MO",
        "Montana" => "MT",
        "Nebraska" => "NE",
        "Nevada" => "NV",
        "New Hampshire" => "NH",
        "New Jersey" => "NJ",
        "New Mexico" => "NM",
        "New York" => "NY",
        "North Carolina" => "NC",
        "North Dakota" => "ND",
        "Ohio" => "OH",
        "Oklahoma" => "OK",
        "Oregon" => "OR",
        "Pennsylvania" => "PA",
        "Puerto Rico" => "PR",
        "Rhode Island" => "RI",
        "South Carolina" => "SC",
        "South Dakota" => "SD",
        "Tennessee" => "TN",
        "Texas" => "TX",
        "Utah" => "UT",
        "Vermont" => "VT",
        "Virginia" => "VA",
        "Washington" => "WA",
        "West Virginia" => "WV",
        "Wisconsin" => "WI",
        "Wyoming" => "WY"
    )
    get(d, state, state)
end

# load population data for a county
function load_pop_county(county::String, state::String)
    df = CSV.read("county_pop.csv", DataFrame)
    if county == "New York City"
        pop = df[(df.county.=="New York").&(df.state.=="New York"), :population][1] +
              df[(df.county.=="Kings").&(df.state.=="New York"), :population][1] +
              df[(df.county.=="Queens").&(df.state.=="New York"), :population][1] +
              df[(df.county.=="Bronx").&(df.state.=="New York"), :population][1] +
              df[(df.county.=="Richmond").&(df.state.=="New York"), :population][1]
    else
        pop = df[(df.county.==county).&(df.state.==state), :population][1]
    end
    pop
end

# load population data for a state
function load_pop_state(state::String)
    df = CSV.read("state_pop.csv", DataFrame)
    df[df.state.==state, :population][1]
end

# load population data for the nation
function load_pop_nation()
    df = CSV.read("state_pop.csv", DataFrame)
    sum(df[:, :population])
end

# load case and death data for a single county from a dataframe
function load_cd_county(df::DataFrame, county::String, state::String)
    Matrix(df[(df.county.==county).&(df.state.==state), [:cases, :deaths]])
end

# load case and death data for all states
function load_cd_states()
    df = CSV.read("covid-19-data/us-states.csv", DataFrame)
    date = df.date[end]
    states = unique(df.state)
    n_states = length(states) # not 50, includes territories, etc
    cd_data = Array{Any,1}(undef, n_states)
    for i = 1:n_states
        cd_data[i] = Matrix(df[df.state.==states[i], [:cases, :deaths]])
    end
    states, cd_data, date
end

# load case and death data for the nation
function load_cd_nation()
    df = CSV.read("covid-19-data/us.csv", DataFrame)
    cd_data = Matrix(df[:, [:cases, :deaths]])
end

# apply bidirectional first-order IIR to smooth data
function smooth_iir_1(y, alpha = 0.3)
    n = length(y)
    # first order IIR forward
    z0 = copy(y)
    z = y[1]
    for i = 1:n
        z = alpha * y[i] + (1 - alpha) * z
        z0[i] = z
    end
    # first order IIR backward
    z1 = copy(y)
    z = y[end]
    for i = n:-1:1
        z = alpha * y[i] + (1 - alpha) * z
        z1[i] = z
    end
    (z0 + z1) * 0.5
end

# estimate the differences of a time series and fill in gaps
function est_diff(y)
    n = length(y)
    d = diff(y)
    # remove negative dips
    d[d.<0.0] .= 0.0
    # redistribute spikes after zeros
    i0 = 1
    i = 2
    while i < n
        # X0 = hcat(ones(2 * m + 1), collect(-m:m))
        # if zero is found
        if d[i] == 0
            # find first non-zero
            while d[i] == 0.0 && i < n - 1
                i += 1
            end
            if d[i] == 0.0 && i == n - 1
                break
            end
            # redistribute d[i] to make a line from d[i0] to d[i]
            m = i - i0
            @assert(m > 0)
            t = (m * (m + 1.0)) / 2.0
            a = (d[i] - m * d[i0]) / t
            if a * m + d[i0] >= 0
                X = hcat(ones(m), collect(1:m))
                d[i0+1:i] = X * [d[i0], a]
            end
        end
        i0 = i
        i += 1
    end
    # remove huge spikes
    b = [-3, -2, -1, 1, 2, 3]
    for i = 4:n-4
        m = mean(d[b.+i])
        if d[i] > 5 * max(1, m)
            d[i] = m
        end
    end
    d
end

# sum of derivatives of sigmoids (generalized logistic function)
function sigmoid_deriv_sum(a, t)
    @assert(length(a) % 4 == 0)
    b = reshape(a, 4, :)
    m = size(b, 2)
    y = zeros(length(t))
    for i = 1:m
        e = exp.((-5.0 * b[2, i]) * (t .- b[3, i]))
        y .+= (4.0 * b[1, i]) * e ./ ((1.0 .+ e) .^ (1.0 + b[4, i]))
    end
    y
end

# Jacobian of the weighted residual in a sum of sigmoids fit
function sigmoid_deriv_sum_jac(a, t, w)
    @assert(length(a) % 4 == 0)
    b = reshape(a, 4, :)
    m = size(b, 2)
    n = length(t)
    jac = zeros(n, 4, m)
    u = (4.0 / sqrt(n)) * w
    for i = 1:m
        e = exp.(-5.0 * b[2, i] * (t .- b[3, i]))
        d = (1 .+ e)
        g = b[1, i] * u .* e
        h = g .* d .^ (-(2.0 + b[4, i])) .* (1.0 .- b[4, i] * e)
        jac[:, 1, i] = u .* e ./ (d .^ (1.0 + b[4, i]))
        jac[:, 2, i] = h .* (-5.0 * (t .- b[3, i]))
        jac[:, 3, i] = h * (5.0 * b[2, i])
        jac[:, 4, i] = g .* d .^ (-(1.0 + b[4, i])) .* (-log.(d))
    end
    reshape(jac, n, :)
end

# test for Jacobian of the residual in a sum of sigmoids fit
function test_sigmoid_deriv_sum_jac()
    n = 100
    t = collect(range(0, stop = 1, length = n))
    w = ones(n)
    a = [1.5, 0.6, 0.8, 1, 2, 0.3, 1, 1.1]
    y0 = rand(n)
    # direct computation of Jacobian of residual
    jac0 = sigmoid_deriv_sum_jac(a, t, w)
    # finite difference approximation
    jac1 = zeros(n, 8)
    eps = 1e-6
    for i = 1:8
        a1 = copy(a)
        a1[i] -= eps
        res0 = w .* (sigmoid_deriv_sum(a1, t) - y0) / sqrt(n)
        a1[i] += 2 * eps
        res1 = w .* (sigmoid_deriv_sum(a1, t) - y0) / sqrt(n)
        jac1[:, i] = (res1 - res0) / (2 * eps)
    end
    println(maximum(abs.(jac0 - jac1)))
end

# fit a sigmoid sum using the Gauss-Newton method to optimize weighted least squares
function fit_sigmoid_deriv_sum_wls_gn(x_init, t, w, y0_raw, iter_max, debug::Bool = false)
    # initialize return value
    x_opt = copy(x_init)

    # normalize data to fit
    n = length(x_init)
    y0_max = maximum(y0_raw)
    y0 = y0_raw / y0_max

    # Levenberg-Marquardt parameters
    eps_min = 1e-6
    eps_max = 1e+2
    eps_inc = 3.0
    eps_dec = 0.5

    # set up initial values
    eps = 1.0
    eps_it = 0
    x1 = copy(x_init)
    res = w .* (sigmoid_deriv_sum(x1, t) - y0) / sqrt(n)
    jac = sigmoid_deriv_sum_jac(x1, t, w)
    err1 = norm(res)
    err_opt = err1

    # regularization coefficients
    x_init[1:4:end] .= 0.0
    c_reg = 3e-2 * ones(n)
    c_reg[3:4:end] .= 0.0

    # Gauss-Newton iterations
    for it = 1:iter_max

        # update parameters
        x0 = copy(x1)
        x1 -= (jac' * jac + eps * I) \ (jac' * res) * 0.5 - c_reg .* (x1 - x_init)
        x1[1:4:end] = max.(min.(x1[1:4:end], 100.0), 0.01)
        x1[2:4:end] = max.(min.(x1[2:4:end], 10.0), 0.1)
        x1[3:4:end] = max.(min.(x1[3:4:end], 1.1), -0.1)
        x1[4:4:end] = max.(min.(x1[4:4:end], 10.0), 0.1)

        # update weighted residual and its Jacobian
        err0 = err1
        res = w .* (sigmoid_deriv_sum(x1, t) - y0) / sqrt(n)
        jac = sigmoid_deriv_sum_jac(x1, t, w)
        err1 = norm(res)

        # lock in parameters if the error goes down
        if err1 < err_opt
            err_opt = err1
            x_opt = copy(x1)
        end

        # update epsilon based on the error going up or down
        if err1 <= err0 + 1e-2
            eps *= eps_dec
            if eps < eps_min
                eps = eps_min
            end
            if eps == eps_min
                if eps_it == 3
                    break
                else
                    eps_it += 1
                end
            end
        else # backtrack
            if eps == eps_max
                if eps_it == 3
                    break
                else
                    eps_it += 1
                end
            else
                eps_it = 0
                err1 = err0
                x1 = copy(x0)
                eps *= eps_inc
                if eps > eps_max
                    eps = eps_max
                end
            end
        end

        # possibly output debug information
        if debug
            @printf("Iter: %2d,    Eps: %.6f,    Err: %.6f,    Best Err: %.6f\n",
                it, eps, err1, err_opt)
        end
    end
    if debug
        println()
    end

    x_opt[1:4:end] *= y0_max
    x_opt, err_opt / length(y0)
end

# fit sum of sigmoid functions (generalized logistics) using iteratively re-weighted least squares
function fit_sigmoid_deriv_sum(y0, n_ext)
    n = length(y0)
    t0 = collect(range(0, stop = 1, length = n))
    w0 = ones(n)
    w0[end-5:end] .= [0.8, 0.6, 0.4, 0.2, 0, 0]
    opt_a1 = []
    opt_b1 = Inf
    for n_sig = 1:10
        a0 = ones(4, n_sig)
        a0[1, :] .= 0.5
        for i_sig = 1:n_sig
            a0[3, i_sig] = (i_sig - 1.0) / n_sig
        end
        a1, err = fit_sigmoid_deriv_sum_wls_gn(a0[:], t0, w0, y0, 300)
        b1 = n_sig / (n / 60) + log(err) # bayesian information criterion
        if b1 < opt_b1
            opt_b1 = b1
            opt_a1 = a1
        end
    end
    t1 = [t0; 1.0 .+ (t0[2] - t0[1]) * collect(1:n_ext)]
    y1 = sigmoid_deriv_sum(opt_a1, t1)
    y1_max = sigmoid_deriv_sum(opt_a1, 10.0)
    y1, y1_max
end

function filter_first_of_month(t0)
    # find unique first days of each month
    t1 = unique(Dates.firstdayofmonth.(t0))
    if t1[1] < t0[1]
        t1 = t1[2:end]
    end
    t1
end

# hack to deal with twinx() undocumented/buggy behavior for dates: convert to int64
function dates_int64_xticks(t0)
    x = collect(0:length(t0)-1)
    x_offset = Dates.value(t0[1])
    tt = filter_first_of_month(t0)
    xt = zeros(length(tt))
    j = 1
    for i = 1:length(t0)
        if tt[j] == t0[i]
            xt[j] = x[i]
            j += 1
            if j > length(tt)
                break
            end
        end
    end
    xt, x, x_offset
end

# plot daily cases/deaths per day
function plot_dcd(cd, title_str, date, pop, set_ymax, n_days_ext, n_days_plot)
    # cutoff for fit
    cd = cd[1:end-3, :]
    n_days_plot = min(size(cd, 1), n_days_plot)
    cd = cd[end-(n_days_plot-1):end, :]
    @assert(size(cd, 1) == n_days_plot)
    # normalize by population (per 100,000)
    c0 = cd[:, 1] / pop * 1e5
    d0 = cd[:, 2] / pop * 1e5
    # compute new cases
    dc0 = est_diff(c0)
    dd0 = est_diff(d0)
    # remove trailing zeros
    n_days = length(dc0)
    n_days_zero = 0
    while dc0[end-n_days_zero] == 0.0
        n_days_zero += 1
        # exclude case when cases go smoothly to zero
        if 0.0 < dc0[end-n_days_zero] && dc0[end-n_days_zero] < 1.0
            n_days_zero = 0
            break
        end
    end
    n_days_plot -= n_days_zero
    dc0 = dc0[1:end-n_days_zero]
    dd0 = dd0[1:end-n_days_zero]
    t0 = collect(date-Day(n_days - 1):Day(1):date-Day(n_days_zero))
    t1 = collect(date-Day(n_days - 1):Day(1):date-Day(n_days_zero)+Day(n_days_ext))
    # smooth data
    dc0 = smooth_iir_1(dc0)
    dd0 = smooth_iir_1(dd0, 0.15)
    # fit sigmoid derivative to new cases
    dc1, dc1_max = fit_sigmoid_deriv_sum(dc0, n_days_ext)
    dd1, dd1_max = fit_sigmoid_deriv_sum(dd0, n_days_ext)
    # current cases
    curr_cases = round(dc0[end], digits = 1)
    title_str = title_str * "  [$curr_cases]"
    # plot cases
    if set_ymax
        c_ymax = 150.0
        d_ymax = 2.0
    else
        c_ymax = max(1.05 * max(maximum(dc0), maximum(dc1[1:n_days_plot])), 10.0)
        d_ymax = max(2.00 * max(maximum(dd0), maximum(dd1[1:n_days_plot])), 1.0)
    end
    xt, x, x_offset = dates_int64_xticks(t0)
    p = plot(x, dc0, w = 2, color = :black, xticks = false, ylabel = "Daily Cases (per 100k)", yguidefont = font(8),
        title = title_str, titlefontsize = 10, legend = false, ylims = (0, c_ymax))
    tx = twinx()
    plot!(tx, dd0, w = 2, color = :red, grid = :off, xticks = false, legend = false,
        ylabel = "Daily Deaths (per 100k)", yguidefont = font(8, :red), ylims = (0, d_ymax))
    xt, x, x_offset = dates_int64_xticks(t1)
    plot!(x, dc1, color = :grey, xticks = xt, xtickfontsize = 8, legend = false, ylims = (0, c_ymax))
    plot!(tx, dd1, color = :orange, grid = :off, xticks = false, legend = false, ylims = (0, d_ymax))
    plot!(xformatter = x -> Dates.monthabbr(Date(Dates.UTD(x + x_offset))))
    p
end

# plot cases/deaths over time and new cases versus total cases
function plot_county_list(df, county_list, plot_national = true)
    set_ymax = false
    n_days_ext = 5
    n_days_plot = 150
    n_cnty = length(county_list)
    n_plot = n_cnty + plot_national
    date = df[end, :date]
    p = Array{Any,1}(undef, n_plot)
    for i = 1:n_cnty
        county, state = county_list[i]
        cd = load_cd_county(df, county, state)
        pop = load_pop_county(county, state)
        state_abbr = abbreviate_state(state)
        if county == "New York City"
            title_str = "New York City, NY"
        else
            title_str = "$county County, $state_abbr"
        end
        p[i] = plot_dcd(cd, title_str, date, pop, set_ymax, n_days_ext, n_days_plot)
    end
    if plot_national
        cd = load_cd_nation()
        pop = load_pop_nation()
        p[end] = plot_dcd(cd, "United States", date, pop, set_ymax, n_days_ext, n_days_plot)
    end
    plot(p..., layout = (n_plot, 1), size = (600, 250 * n_plot), right_margin = 18mm, left_margin = 8mm)
end

# load data
print("Loading Covid-19 county cases/deaths data...")
t = @elapsed df = CSV.read("covid-19-data/us-counties.csv", DataFrame)
@printf("finished in %.2f sec.\n", t)

# plot data
# display(plot_county_list(df, [("Winnebago", "Illinois"),
#                               ("Boone", "Illinois"),
#                               ("DeKalb", "Illinois")]))

display(plot_county_list(df, [("Santa Clara", "California"),
    ("Santa Barbara", "California"),
    ("Los Angeles", "California"),
    ("New York City", "New York")]))


text:
# Autogenerated wrapper script for openbabel_jll for aarch64-apple-darwin
export libinchi, libopenbabel, obabel, obconformer, obdistgen, obenergy, obfit, obfitall, obgen, obgrep, obminimize, obmm, obprobe, obprop, obrms, obrotamer, obrotate, obspectrophore, obsym, obtautomer, obthermo, roundtrip

using XML2_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("openbabel")
JLLWrappers.@declare_library_product(libinchi, "@rpath/libinchi.0.dylib")
JLLWrappers.@declare_library_product(libopenbabel, "@rpath/libopenbabel.7.dylib")
JLLWrappers.@declare_executable_product(obabel)
JLLWrappers.@declare_executable_product(obconformer)
JLLWrappers.@declare_executable_product(obdistgen)
JLLWrappers.@declare_executable_product(obenergy)
JLLWrappers.@declare_executable_product(obfit)
JLLWrappers.@declare_executable_product(obfitall)
JLLWrappers.@declare_executable_product(obgen)
JLLWrappers.@declare_executable_product(obgrep)
JLLWrappers.@declare_executable_product(obminimize)
JLLWrappers.@declare_executable_product(obmm)
JLLWrappers.@declare_executable_product(obprobe)
JLLWrappers.@declare_executable_product(obprop)
JLLWrappers.@declare_executable_product(obrms)
JLLWrappers.@declare_executable_product(obrotamer)
JLLWrappers.@declare_executable_product(obrotate)
JLLWrappers.@declare_executable_product(obspectrophore)
JLLWrappers.@declare_executable_product(obsym)
JLLWrappers.@declare_executable_product(obtautomer)
JLLWrappers.@declare_executable_product(obthermo)
JLLWrappers.@declare_executable_product(roundtrip)
function __init__()
    JLLWrappers.@generate_init_header(XML2_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libinchi,
        "lib/libinchi.0.4.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopenbabel,
        "lib/libopenbabel.7.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        obabel,
        "bin/obabel",
    )

    JLLWrappers.@init_executable_product(
        obconformer,
        "bin/obconformer",
    )

    JLLWrappers.@init_executable_product(
        obdistgen,
        "bin/obdistgen",
    )

    JLLWrappers.@init_executable_product(
        obenergy,
        "bin/obenergy",
    )

    JLLWrappers.@init_executable_product(
        obfit,
        "bin/obfit",
    )

    JLLWrappers.@init_executable_product(
        obfitall,
        "bin/obfitall",
    )

    JLLWrappers.@init_executable_product(
        obgen,
        "bin/obgen",
    )

    JLLWrappers.@init_executable_product(
        obgrep,
        "bin/obgrep",
    )

    JLLWrappers.@init_executable_product(
        obminimize,
        "bin/obminimize",
    )

    JLLWrappers.@init_executable_product(
        obmm,
        "bin/obmm",
    )

    JLLWrappers.@init_executable_product(
        obprobe,
        "bin/obprobe",
    )

    JLLWrappers.@init_executable_product(
        obprop,
        "bin/obprop",
    )

    JLLWrappers.@init_executable_product(
        obrms,
        "bin/obrms",
    )

    JLLWrappers.@init_executable_product(
        obrotamer,
        "bin/obrotamer",
    )

    JLLWrappers.@init_executable_product(
        obrotate,
        "bin/obrotate",
    )

    JLLWrappers.@init_executable_product(
        obspectrophore,
        "bin/obspectrophore",
    )

    JLLWrappers.@init_executable_product(
        obsym,
        "bin/obsym",
    )

    JLLWrappers.@init_executable_product(
        obtautomer,
        "bin/obtautomer",
    )

    JLLWrappers.@init_executable_product(
        obthermo,
        "bin/obthermo",
    )

    JLLWrappers.@init_executable_product(
        roundtrip,
        "bin/roundtrip",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module NSOS
using StatsBase
using Colors
using Compose
# package code goes here
export timeheadway
include("overtaking.jl")

end # module


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iotfleethub
using AWS.Compat
using AWS.UUIDs

"""
    create_application(application_name, role_arn)
    create_application(application_name, role_arn, params::Dict{String,<:Any})

Creates a Fleet Hub for AWS IoT Device Management web application.  Fleet Hub for AWS IoT
Device Management is in public preview and is subject to change.

# Arguments
- `application_name`: The name of the web application.
- `role_arn`: The ARN of the role that the web application assumes when it interacts with
  AWS IoT Core.  The name of the role must be in the form AWSIotFleetHub_random_string .

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"applicationDescription"`: An optional description of the web application.
- `"clientToken"`: A unique case-sensitive identifier that you can provide to ensure the
  idempotency of the request. Don't reuse this client token if a new idempotent request is
  required.
- `"tags"`: A set of key/value pairs that you can use to manage the web application
  resource.
"""
create_application(applicationName, roleArn; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("POST", "/applications", Dict{String, Any}("applicationName"=>applicationName, "roleArn"=>roleArn, "clientToken"=>string(uuid4())); aws_config=aws_config)
create_application(applicationName, roleArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("POST", "/applications", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("applicationName"=>applicationName, "roleArn"=>roleArn, "clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    delete_application(application_id)
    delete_application(application_id, params::Dict{String,<:Any})

Deletes a Fleet Hub for AWS IoT Device Management web application.  Fleet Hub for AWS IoT
Device Management is in public preview and is subject to change.

# Arguments
- `application_id`: The unique Id of the web application.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A unique case-sensitive identifier that you can provide to ensure the
  idempotency of the request. Don't reuse this client token if a new idempotent request is
  required.
"""
delete_application(applicationId; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("DELETE", "/applications/$(applicationId)", Dict{String, Any}("clientToken"=>string(uuid4())); aws_config=aws_config)
delete_application(applicationId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("DELETE", "/applications/$(applicationId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    describe_application(application_id)
    describe_application(application_id, params::Dict{String,<:Any})

Gets information about a Fleet Hub for AWS IoT Device Management web application.  Fleet
Hub for AWS IoT Device Management is in public preview and is subject to change.

# Arguments
- `application_id`: The unique Id of the web application.

"""
describe_application(applicationId; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("GET", "/applications/$(applicationId)"; aws_config=aws_config)
describe_application(applicationId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("GET", "/applications/$(applicationId)", params; aws_config=aws_config)

"""
    list_applications()
    list_applications(params::Dict{String,<:Any})

Gets a list of Fleet Hub for AWS IoT Device Management web applications for the current
account.  Fleet Hub for AWS IoT Device Management is in public preview and is subject to
change.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"nextToken"`: A token used to get the next set of results.
"""
list_applications(; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("GET", "/applications"; aws_config=aws_config)
list_applications(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("GET", "/applications", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Lists the tags for the specified resource.  Fleet Hub for AWS IoT Device Management is in
public preview and is subject to change.

# Arguments
- `resource_arn`: The ARN of the resource.

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("GET", "/tags/$(resourceArn)", params; aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds to or modifies the tags of the specified resource. Tags are metadata which can be used
to manage a resource.  Fleet Hub for AWS IoT Device Management is in public preview and is
subject to change.

# Arguments
- `resource_arn`: The ARN of the resource.
- `tags`: The new or modified tags for the resource.

"""
tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes the specified tags (metadata) from the resource.  Fleet Hub for AWS IoT Device
Management is in public preview and is subject to change.

# Arguments
- `resource_arn`: The ARN of the resource.
- `tag_keys`: A list of the keys of the tags to be removed from the resource.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_application(application_id)
    update_application(application_id, params::Dict{String,<:Any})

Updates information about a Fleet Hub for a AWS IoT Device Management web application.
Fleet Hub for AWS IoT Device Management is in public preview and is subject to change.

# Arguments
- `application_id`: The unique Id of the web application.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"applicationDescription"`: An optional description of the web application.
- `"applicationName"`: The name of the web application.
- `"clientToken"`: A unique case-sensitive identifier that you can provide to ensure the
  idempotency of the request. Don't reuse this client token if a new idempotent request is
  required.
"""
update_application(applicationId; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("PATCH", "/applications/$(applicationId)", Dict{String, Any}("clientToken"=>string(uuid4())); aws_config=aws_config)
update_application(applicationId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iotfleethub("PATCH", "/applications/$(applicationId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>string(uuid4())), params)); aws_config=aws_config)


text:
using ExecutableSpecifications
using ExecutableSpecifications.Gherkin

@testset "Asserts" begin
    @testset "Assert failure; Assert is 1 == 2; Failure has human readable string 1 == 2" begin
        try
            @expect 1 == 2
        catch ex
            @test ex.assertion == "1 == 2"
        end
    end

    @testset "Assert failure in included file; Assert is 1 == 2; Failure has human readable string 1 == 2" begin
        matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
            using ExecutableSpecifications

            @given "some precondition" begin
                @expect 1 == 2
            end
        """)
        given = Gherkin.Given("some precondition")
        context = ExecutableSpecifications.StepDefinitionContext()

        stepdefinition = ExecutableSpecifications.findstepdefinition(matcher, given)

        stepfailed = stepdefinition.definition(context)

        @test stepfailed.assertion == "1 == 2"
    end

    @testset "Assert failure in included file; Assert is isempty([1]); Failure has human readable string isempty([1])" begin
        matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
            using ExecutableSpecifications

            @given "some precondition" begin
                @expect isempty([1])
            end
        """)
        given = Gherkin.Given("some precondition")
        context = ExecutableSpecifications.StepDefinitionContext()

        stepdefinition = ExecutableSpecifications.findstepdefinition(matcher, given)

        stepfailed = stepdefinition.definition(context)

        @test stepfailed.assertion == "isempty([1])"
    end
end


text:
@gen function in_place_drift_randomness(trace, idx::Int, pos_var, rot_conc, sample_contact_faces)
    offset ~ gaussianVMF(IDENTITY_POSE, pos_var, rot_conc)
    offset
end

function in_place_drift_involution(
    trace,
    fwd_choices,
    fwd_ret,
    proposal_args;
    check = false,
)
    constraints = choicemap()
    bwd_choices = choicemap()

    idx, _,_, sample_contact_faces = proposal_args
    offset = fwd_choices[:offset]
    bwd_choices[:offset] = inverse_pose(offset)

    structure = get_structure(trace)

    poses = get_poses(trace)
    object_boxes = get_object_boxes(trace)

    new_pose = get_c_relative_to_a(poses[idx], offset)

    if isFloating(structure, idx)
        constraints[floating_pose_addr(idx)] = new_pose
    else
        parent_idx = parent(structure, idx)
        @assert !isnothing(parent_idx)

        if sample_contact_faces
            (parent_face, child_face), _ =
                get_closest_planar_contact(poses[parent_idx], object_boxes[parent_idx], new_pose, object_boxes[idx])
            parent_face = S.BOX_SURFACE_IDS[parent_face]
            child_face = S.BOX_SURFACE_IDS[child_face]
        else
            parent_face = trace[contact_addr(idx, :parent_face)]
            child_face = trace[contact_addr(idx, :child_face)]
        end
        pc = S.closestApproximatingContact(
            (poses[parent_idx] * S.getContactPlane(object_boxes[parent_idx], parent_face)) \
            (new_pose * S.getContactPlane(object_boxes[idx], child_face)),
        )
        contact = S.ShapeContact(parent_face, (), child_face, (), pc)
        cm = contact_submap_from_shape_contact(contact)
        set_submap!(constraints, (contact_addr(idx)), cm)

        for child_idx in children(structure, idx)
            if sample_contact_faces
                (parent_face, child_face), _ =
                    get_closest_planar_contact(new_pose, object_boxes[idx], poses[child_idx], object_boxes[child_idx])
                parent_face = S.BOX_SURFACE_IDS[parent_face]
                child_face = S.BOX_SURFACE_IDS[child_face]
            else
                parent_face = trace[contact_addr(child_idx, :parent_face)]
                child_face = trace[contact_addr(child_idx, :child_face)]
            end

            pc = S.closestApproximatingContact(
                (new_pose * S.getContactPlane(object_boxes[idx], parent_face)) \
                (poses[child_idx] * S.getContactPlane(object_boxes[child_idx], child_face)),
            )
            contact = S.ShapeContact(parent_face, (), child_face, (), pc)
            cm = contact_submap_from_shape_contact(contact)
            set_submap!(constraints, contact_addr(child_idx), cm)
        end
    end

    args = get_args(trace)
    argdiffs = map((_) -> NoChange(), args)
    new_trace, weight, _, _ = Gen.update(trace, args, argdiffs, constraints)
    return (new_trace, bwd_choices, weight)
end

export in_place_drift_randomness, in_place_drift_involution


text:
mutable struct Machine
    acc::Int
end
Machine() = Machine(0)

reset!(m::Machine) = (m.acc = 0; m)

function execute!(mach::Machine, program)
    i, N = 1, length(program)
    seen = falses(N)
    while i ≤ N && !seen[i]
        seen[i] = true
        op, arg = program[i]
        if op == "acc"
            mach.acc += arg
            i += 1
        elseif op == "jmp"
            i += arg
        elseif op == "nop"
            i += 1
        end
    end
    i != length(program) + 1, mach.acc
end

function main()
    part1 = -1

    lines = map(readlines("input.txt")) do line
        op, arg = split(line)
        (op, parse(Int, arg))
    end

    machine = Machine()

    _, part1 = execute!(machine, lines)

    for i in findall(l -> first(l) == "jmp", lines)
        reset!(machine)

        lines[i] = ("nop", last(lines[i]))
        early, value = execute!(machine, lines)
        lines[i] = ("jmp", last(lines[i]))
        if !early
            return part1, value
        end
    end

    for i in findall(l -> first(l) == "nop", lines)
        reset!(machine)

        lines[i] = ("jmp", last(lines[i]))
        early, value = execute!(machine, lines)
        lines[i] = ("nop", last(lines[i]))
        if !early
            return part1, value
        end
    end

    error("failed to debug")
end

let (part1, part2) = main()
    @info "Solutions" part1 part2
end

max_stars_repo_name: dglmoore/advent-of-code

max_issues_repo_name: dglmoore/advent-of-code

max_forks_repo_name: dglmoore/advent-of-code


text:
#batched CuArray gemm by BatchedRoutines.jl
function batched_gemm!(transA::AbstractChar,
                       transB::AbstractChar,
                       alpha::Float32,
                       A::CuArrays.CuArray{Float32, 3},
                       B::CuArrays.CuArray{Float32, 3},
                       beta::Float32,
                       C::CuArrays.CuArray{Float32, 3})
    CuArrays.CUBLAS.gemm_strided_batched!(transA, transB, alpha, A, B, beta, C)
end

function batched_gemm!(transA::AbstractChar,
                       transB::AbstractChar,
                       alpha::Float64,
                       A::CuArrays.CuArray{Float64, 3},
                       B::CuArrays.CuArray{Float64, 3},
                       beta::Float64,
                       C::CuArrays.CuArray{Float64, 3})
    CuArrays.CUBLAS.gemm_strided_batched!(transA, transB, alpha, A, B, beta, C)
end


text:
using Test, Printf

using CLIMA.TurbulenceConvection.FiniteDifferenceGrids
using CLIMA.TurbulenceConvection.StateVecs

n_elems_real = 10 # number of elements

grid = Grid(0.0, 1.0, n_elems_real)
FT = eltype(grid)
vars = ( (:ρ_0, DomainSubSet(gm=true)),
         (:a,   DomainSubSet(gm=true,en=true,ud=true)),
         (:w,   DomainSubSet(gm=true,en=true,ud=true)) )
dd = DomainDecomp(gm=1,en=1,ud=2)
q = StateVec(vars, grid, dd)
q_compare = deepcopy(q)
idx = DomainIdx(q)
gm, en, ud, sd, al = allcombinations(idx)

@testset "State Vectors" begin
  for k in over_elems(grid)
    q[:ρ_0, k] = 2.0
  end
  @test all([q[:ρ_0, k] for k in over_elems(grid)] .== 2.0)

  for k in over_elems(grid)
    q[:w, k, gm] = 3.0
  end
  @test all([q[:w, k, gm] for k in over_elems(grid)] .== 3.0)

  @test_throws BoundsError q[:w, 1, 1000] = 3.0
  @test_throws BoundsError q[:ρ_0, 1, en] = 3.0

  for k in over_elems(grid)
    ρ_0_e = q[:ρ_0, k]
    for i in alldomains(idx)
      w_0_e_i = q[:w, k, i]
    end
  end
  sprint(show, q)

  assign!(q, grid, 0.0)
  k = 1
  apply_Dirichlet!(q, :ρ_0, grid, FT(2), Zmin())
  @test q[:ρ_0, k] ≈ 4
  apply_Neumann!(q, :ρ_0, grid, -2/grid.Δz, Zmin())
  @test q[:ρ_0, k] ≈ 2

  k = grid.n_elem
  apply_Dirichlet!(q, :ρ_0, grid, FT(2), Zmax())
  @test q[:ρ_0, k] ≈ 4
  apply_Neumann!(q, :ρ_0, grid,  2/grid.Δz, Zmax())
  @test q[:ρ_0, k] ≈ 2

  q[:ρ_0, 1] = 1.0
  q[:ρ_0, 2] = 2.0
  q[:ρ_0, 3] = 3.0
  q[:a, 1] = 1.0
  q[:a, 2] = 2.0
  q[:a, 3] = 3.0
  ρα_0_cut = q[:ρ_0, Cut(2)].*q[:a, Cut(2)]
  @test all(ρα_0_cut .== ρα_0_cut)
  ρα_0_dual = q[:ρ_0, Dual(2)].*q[:a, Dual(2)]
  @test all(ρα_0_dual .== ρα_0_dual)

  q[:a, 3] = 1
  @test !isnan(q)
  @test !isinf(q)
  q[:a, 3] = 1/0
  @test isinf(q)
  q[:a, 3] = NaN
  @test isnan(q)

  assign!(q, grid, 2.0)
  vals = [q[ϕ,k,i] for ϕ in var_names(q) for i in over_sub_domains(q,ϕ) for k in over_elems(grid)]
  @test all([q[ϕ,k,i] for ϕ in var_names(q) for i in over_sub_domains(q,ϕ) for k in over_elems(grid)] .≈ 2.0)

  assign!(q, grid, 0.0)

  assign!(q, (:w), grid, 2.0)
  @test all([q[:w,k,i] for i in over_sub_domains(q,:w) for k in over_elems(grid)] .≈ 2.0)

  assign!(q, grid, 0.0)

  assign_real!(q, :w, grid, [2.0 for x in over_elems_real(grid)], gm)
  @test all([q[:w,k,gm] for i in over_sub_domains(q,:w) for k in over_elems_real(grid)] .≈ 2.0)
  @test all([q[:w,k,gm] for i in over_sub_domains(q,:w) for k in over_elems_ghost(grid)] .≈ 0.0)

  assign!(q, grid, 2.0)
  assign_ghost!(q, :w, grid, 0.0, gm)
  @test all([q[:w,k,gm] for i in over_sub_domains(q,:w) for k in over_elems_real(grid)] .≈ 2.0)
  @test all([q[:w,k,gm] for i in over_sub_domains(q,:w) for k in over_elems_ghost(grid)] .≈ 0.0)

  assign!(q_compare, var_names(q_compare), grid, FT(2.0))
  assign!(q, var_names(q_compare), grid, FT(1.0))

  D = compare(q, q, grid, eps(Float32))
  for (k,v) in D
    @test all(v)
  end
  D = compare(q, q_compare, grid, eps(Float32))
  for (k,v) in D
    @test !any(v)
  end

  assign!(q, grid, 0.0)
  assign_real!(q, :w, grid, [2.0 for x in over_elems_real(grid)], gm)
  extrap_0th_order!(q, :w, grid, gm)
  @test all([q[:w,k,gm] for k in over_elems(grid)] .≈ 2.0)

  for k in over_elems(grid)
    q[:a, k] = 2
  end

  assign_ghost!(q, :a, grid, 0.0)
  @test all(q[:a, k] ≈ 0.0 for k in over_elems_ghost(grid))

  extrap!(q, :a, grid)
  k = 1+grid.n_ghost
  @test q[:a, k-1] ≈ 2*q[:a, k] - q[:a, k+1]
  k = grid.n_elem-grid.n_ghost
  @test q[:a, k+1] ≈ 2*q[:a, k] - q[:a, k-1]

  assign!(q, grid, 0.0)
  ki_1 = first_interior(grid, Zmin())
  ki_2 = first_interior(grid, Zmax())
  kg_1 = first_ghost(grid, Zmin())
  kg_2 = first_ghost(grid, Zmax())

  q[:w,ki_1,gm] = 1.0
  q[:w,ki_2,gm] = 1.0
  extrap!(q, :w, grid, gm)
  @test q[:w,kg_1,gm] ≈ 2.0
  @test q[:w,kg_2,gm] ≈ 2.0

  # Extrapolate field
  q[:w,ki_1,gm] = 2.0
  q[:w,kg_2,gm] = 2.0
  extrap!(q, :w, grid, Zmin(), gm)
  q[:w,ki_1,gm] ≈ 4.0
  q[:w,kg_2,gm] ≈ 2.0

  q[:w,kg_1,gm] = 2.0
  q[:w,ki_2,gm] = 2.0
  extrap!(q, :w, grid, Zmax(), gm)
  q[:w,kg_1,gm] ≈ 2.0
  q[:w,kg_2,gm] ≈ 4.0

  # Extrapolate field with specified boundary value
  q[:w,ki_1,gm] = 0.0
  q[:w,kg_2,gm] = 0.0
  extrap!(q, :w, grid, 1.0, Zmin(), gm)
  q[:w,ki_1,gm] ≈ 2.0
  q[:w,kg_2,gm] ≈ 0.0

  q[:w,kg_1,gm] = 0.0
  q[:w,ki_2,gm] = 0.0
  extrap!(q, :w, grid, 1.0, Zmax(), gm)
  q[:w,kg_1,gm] ≈ 0.0
  q[:w,kg_2,gm] ≈ 2.0

  @test nice_string(:θ) == "theta"
  @test nice_string(:ρ) == "rho"
  @test nice_string(:∇) == "grad"
  @test nice_string(:εδ) == "entr-detr"

  @test var_string(q, :a, gm) == "a_gm"
  @test var_string(q, :a, en) == "a_en"
  @test var_string(q, :a, ud[1]) == "a_ud_1"
  @test var_string(q, :a, ud[2]) == "a_ud_2"
  @test var_suffix(q, :a, gm) == "_gm"
  @test var_suffix(q, :a, en) == "_en"
  @test var_suffix(q, :a, ud[1]) == "_ud_1"
  @test var_suffix(q, :a, ud[2]) == "_ud_2"
  @test_throws BoundsError var_suffix(q, :a, 1000)
end

output_root = joinpath("output", "tests", "StateVecFuncs")

n_subdomains = 3 # number of sub-domains
n_elems_real = 10 # number of elements

dd = DomainDecomp(gm=1,en=1,ud=1)
grid = Grid(0.0, 1.0, n_elems_real)
vars = ( (:ρ_0, DomainSubSet(gm=true)),
         (:a,   DomainSubSet(gm=true,en=true,ud=true)),
         (:w,   DomainSubSet(gm=true,en=true,ud=true)),
         (:ϕ,   DomainSubSet(gm=true,en=true,ud=true)),
         (:ψ,   DomainSubSet(gm=true,en=true,ud=true))
         )
q = StateVec(vars, grid, dd)
vars = ((:cv_ϕ_ψ, DomainSubSet(gm=true,en=true,ud=true)),
        (:TCV_ϕ_ψ, DomainSubSet(gm=true)))
tmp = StateVec(vars, grid, dd)

idx = DomainIdx(q)
gm, en, ud, sd, al = allcombinations(idx)

@testset "IO" begin
  for k in over_elems(grid)
    q[:ρ_0, k] = 1.0*k
    for i in al
      q[:w, k, i] = 2.0*k
      q[:a, k, i] = 3.0*k
    end
  end
  q_2 = deepcopy(q)
  export_state(q, grid, output_root, "q.csv")
  import_state!(q_2, grid, output_root, "q.csv")
  D = compare(q, q_2, grid, eps(Float32))
  for ϕ in var_names(q)
    @test all(D[ϕ])
  end
  rm(joinpath(output_root,"q.csv"))
end

@testset "Domain average, distribute, diagnose environment, total covariance" begin
  q[:a, 1, en] = 0.25
  q[:a, 1, ud[1]] = 0.75
  q[:w, 1, en] = 2
  q[:w, 1, ud[1]] = 2
  domain_average!(q, q, :w, :a, grid)
  @test q[:w, 1, gm] ≈ 2

  q[:w, 1, gm] = 2
  distribute!(q, grid, :w)
  @test q[:w, 1, en] ≈ 2
  for i in ud
    @test q[:w, 1, i] ≈ 2
  end

  q[:a, 1, ud[1]] = 0.1
  q[:a, 1, en] = 0.9
  q[:a, 1, gm] = 1.0

  q[:w, 1, gm] = 2
  for i in ud
    q[:w, 1, i] = 2
  end

  diagnose_environment!(q, grid, :a, (:w))
  @test q[:w, 1, en] ≈ 2

  q[:a, 1, 1] = 0.1
  q[:a, 1, 2] = 0.2
  q[:ϕ, 1, 1] = 1
  q[:ϕ, 1, 2] = 2
  q[:ψ, 1, 1] = 2
  q[:ψ, 1, 2] = 3
  tmp[:cv_ϕ_ψ, 1, 1] = 1.0
  tmp[:cv_ϕ_ψ, 1, 2] = 1.0
  decompose_ϕ_ψ(tcv) = tcv == :TCV_ϕ_ψ ?  (:ϕ , :ψ) : error("Bad init")
  total_covariance!(tmp, q, tmp, :TCV_ϕ_ψ, :cv_ϕ_ψ, :a, grid, decompose_ϕ_ψ)
  @test tmp[:TCV_ϕ_ψ, 1] ≈ 0.32
end

@testset "Boundary conditions" begin

  n_subdomains = 3 # number of sub-domains
  n_elems_real = 10 # number of elements

  dd = DomainDecomp(gm=1,en=1,ud=1)
  grid = Grid(0.0, 1.0, n_elems_real)
  FT = eltype(grid)
  vars = ( (:a,   DomainSubSet(gm=true,en=true,ud=true)),
           (:w,   DomainSubSet(gm=true,en=true,ud=true)),
           )
  q = StateVec(vars, grid, dd)
  vars = (
    (:ρ_0, DomainSubSet(gm=true,en=true,ud=true)),
    (:K, DomainSubSet(gm=true,en=true,ud=true)),
    )
  tmp = StateVec(vars, grid, dd)

  idx = DomainIdx(q)
  gm, en, ud, sd, al = allcombinations(idx)
  assign!(tmp, :K, grid, FT(3))
  assign!(tmp, :ρ_0, grid, FT(2))
  assign!(q, :a, grid, FT(0.1))
  @test bc_source(q, grid, tmp, :w, :ρ_0, :a, :K,
                  Zmin(), FT(2), Dirichlet(), DiffusionAbsorbed()) ≈ FT(240)
  @test bc_source(q, grid, tmp, :w, :ρ_0, :a, :K,
                  Zmin(), FT(2), Neumann(), DiffusionAbsorbed()) ≈ FT(-4)

  # TODO: Validate bc_source for `AdvectionAbsorbed`
  # @test bc_source(q, grid, tmp, :w, :ρ_0, :a, :K,
  #                 Zmin(), FT(2), Dirichlet(), AdvectionAbsorbed()) ≈
  # @test bc_source(q, grid, tmp,
  #                 :w, :ρ_0, :a, :K,
  #                 Zmin(), FT(2), Neumann(), AdvectionAbsorbed()) ≈

end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using Random, LinearAlgebra, Statistics, Optim, DataFrames, CSV, HTTP, GLM, ForwardDiff, FreqTables, Distributions

# read in data
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS4-mixture/nlsw88t.csv"
df = CSV.read(HTTP.get(url).body)
X = [df.age df.white df.collgrad]
Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4,
         df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
y = df.occ_code

# estimate multionomial logit model
include("mlogit_with_Z.jl")
θ_start = [ .0403744; .2439942; -1.57132; .0433254; .1468556; -2.959103; .1020574; .7473086; -4.12005; .0375628; .6884899; -3.65577; .0204543; -.3584007; -4.376929; .1074636; -.5263738; -6.199197; .1168824; -.2870554; -5.322248; 1.307477]
td = TwiceDifferentiable(b -> mlogit_with_Z(b, X, Z, y), θ_start; autodiff = :forward);
θ̂_optim_ad = optimize(td, θ_start, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50));
θ̂_mle_ad = θ̂_optim_ad.minimizer
H  = Optim.hessian!(td, θ̂_mle_ad)
θ̂_mle_ad_se = sqrt.(diag(inv(H)))

# get baseline predicted probabilities
function plogit(θ, X, Z, J)
    α = θ[1:end-1]
    γ = θ[end]
    K = size(X,2)
    bigα = [reshape(α,K,J-1) zeros(K)]
    P = exp.(X*bigα .+ Z*γ) ./ sum.(eachrow(exp.(X*bigα .+ Z*γ)))
    return P
end
println(θ̂_mle_ad)
P = plogit(θ̂_mle_ad, X, Z, length(unique(y)))

# compare with data
modelfitdf = DataFrame(data_freq = 100*convert(Array, prop(freqtable(df,:occ_code))), P = 100*vec(mean(P'; dims=2)))
display(modelfitdf)

# now run a counterfactual:
# suppose people don't care about wages when choosing occupation
# this means setting θ̂_mle_ad[end] = 0
θ̂_cfl = deepcopy(θ̂_mle_ad)
θ̂_cfl[end] = 0
P_cfl = plogit(θ̂_cfl, X, Z, length(unique(y)))
modelfitdf.P_cfl = 100*vec(mean(P_cfl'; dims=2))
modelfitdf.diff = modelfitdf.P_cfl .- modelfitdf.P
modelfitdf.ew   = vec(mean(Z'; dims=2))
display(modelfitdf)

# now do parametric bootstrap
Random.seed!(1234)
invHfix = UpperTriangular(inv(H)) .- Diagonal(inv(H)) .+ UpperTriangular(inv(H))' # this code makes inv(H) truly symmetric (not just numerically symmetric)
d = MvNormal(θ̂_mle_ad,invHfix)
B = 1000
cfl_bs = zeros(8,B)
for b = 1:B
    θ_draw = rand(d)
    Pb = plogit(θ_draw, X, Z, length(unique(y)))
    Pbs = 100*vec(mean(Pb'; dims=2))
    θ_draw[end] = 0
    P_cflb = plogit(θ_draw, X, Z, length(unique(y)))
    Pcfl_bs = 100*vec(mean(P_cflb'; dims=2))
    cfl_bs[:,b] = Pcfl_bs .- Pbs
end
bsdiffL = deepcopy(cfl_bs[:,1])
bsdiffU = deepcopy(cfl_bs[:,1])
for i=1:size(cfl_bs,1)
    bsdiffL[i] = quantile(cfl_bs[i,:],.025)
    bsdiffU[i] = quantile(cfl_bs[i,:],.975)
end
modelfitdf.bs_diffL = bsdiffL
modelfitdf.bs_diffU = bsdiffU
display(modelfitdf)


text:
#############################################################################
# chambolle_pock.jl
# Uses the chambolle_pock algorithm to solve optimization problems consisting
# of two proxable parts f and g.
# For reference look at:
#
# function combi_split(problem, params = CombiParams())
#   @assert ...
#   combi_constraint = problem.combi_constraints[1]
#   z = zeros(size(problem.variable))
#   for t = 1:params.max_iters
#     # compute proximal operator of objective at z
#     reg_objective = problem.objective + 1/2*sumsquares(problem.variable - z)
#     conv_prob = minimize(reg_objective, problem.convex_constraints) # = prox_{objective}(z)
#     solve!(conv_prob)
#     x = evaluate(problem.variable)
#
#     # compute proximal operator of combi_constraint at x
#     z = prox!(combi_constraint, x)
#   end
#   fix!(problem.variable, z) # does this work?
#   problem.solution = z
#   return z
# end
#############################################################################

TOL = 1e-3

export chambolle_pock

function chambolle_pock(f::AbstractExpr, p::AssocPoly,
                        max_iters = 10000 :: Number)
  # Initialization
  x = get_cv(f)[1]
  x0 = zeros(x.size[1])
  u0 = zeros(x.size[1])
  para = Variable(x.size[1])

  # First iteration, in preparation of warm starts
  fix!(para, x0 - u0)
  convex_prob = Problem(:minimize, f + 0.5 * norm(x - para)^2)
  solve!(convex_prob, SCSSolver(verbose = false))
  x1 = (x.value)[:, 1]
  x0 = copy(x1)
  u0 = (u0 + 2 * x1 - x0) - prox(p, (u0 + 2 * x1 - x0), u0)

  # Step 2 iterations
  for iter = 1:max_iters
    fix!(para, x0 - u0)
    solve!(convex_prob, SCSSolver(verbose=false), warmstart = true)
    x1 = (x.value)[:, 1]
    u0 = (u0 + 2 * x1 - x0) - prox(p, (u0 + 2 * x1 - x0), u0)
    if sum(abs.(x1 - x0)) < TOL
      break
    end
    x0 = copy(x1)
  end
  x.value = x0
  return x0
end


text:
"""
    NexProgressBar(nexSerial, name; pid=pid, cid=cid)

A Nextion Progress Bar UI control.

On Nextion side, it have the following properties:
- bco
- pco
- val
"""
struct NexProgressBar <: AbstractNexObject
    _nid::NexID

    viewable::IViewable
    numericalvalued::INumericalValued
    colourable::IColourable
    touchable::ITouchable

    function NexProgressBar(nexSerial::T, name::Name; pid=PageID(), cid=ComponentID()) where {T <: AbstractNexSerial}
        nid = NexID(nexSerial, name, pid, cid)
        new(nid, IViewable(nid), INumericalValued(nid, Int64, RangeNumber{Int64, 0:100}), IColourable(nid), ITouchable(nid))
        #new(nid, IViewable(nid), INumericalValued(nid), IColourable(nid), ITouchable(nid))
    end
end

function setproperty!(obj::NexProgressBar, property::Symbol, new_val)
    _setcommonproperty!(obj, property, new_val)
end

function getproperty(obj::NexProgressBar, property::Symbol)
    _getcommonproperty(obj, property)
end


text:
function irr(annuity::Array{Float64})

    tvmnpv(i,annuity)=begin
         n=collect(1:length(annuity));
         sum(annuity./(1+i).^n)
        end

	if tvmnpv(0.0,annuity) <= 1.0
        return 0.0
    end

	if tvmnpv(0.99, annuity) >= 0
        return 99.9
    end

    f(x)=tvmnpv(x, annuity)

    return find_zeros(f, 0, 1)

end

function npv(annuities::Array{Float64}, rate::Float64)
    pv = 0.0

    for (ix,p) in enumerate(annuities)

        pv += p/(1+rate)^(ix-1)

    end

    return pv

end

function period_pay(P::Float64, rate::Float64, n::Int64)

 return (P*rate*(1.0+rate)^n)/((1.0+rate)^n -1)

end

function amortize(rate::Float64, n::Int64, P::Float64)
    balance = Array{Float64,1}()
    p = period_pay(P,rate,n)

    for k in 1:n
        btmp = P*(1.0+rate)^k - p*(1.0-(1.0+rate)^k)/(-rate)
        push!(balance,btmp)
    end

    intp = vcat([P]*rate, balance[1:end-1]*rate)
    balance = balance[1:end]
    principal_paid = p .- intp

    return  principal_paid, intp, balance, p

end

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
@diffusion_process JansenRit{T} begin
    :dimensions
    process --> 6
    wiener --> 1
    :parameters
    (A, a, B, b, C, α1, α2, νmax, v, r, μy, σy) --> T
    :additional
    constdiff --> true
end
# θ
#param_names(P::JansenRit) = (:A, :a, :B, :b, :C, :α1, :α2, :νmax, :v, :r, :μy, :σy)
#θ =[3.25, 100.0, 22.0, 50.0, 135.0, 0.8, 0.25, 5.0, 6.0, 0.56, 220.0, 2000.0]
#P = JansenRit(θ...)

sigm(x, P::JansenRit) = P.νmax / (1 + exp(P.r*(P.v - x)))
μy(t, P::JansenRit) = P.μy #constant
C1(P::JansenRit) = P.C
C2(P::JansenRit) = P.α1*P.C
C3(P::JansenRit) = P.α2*P.C
C4(P::JansenRit) = P.α2*P.C

function b(t, x, P::JansenRit)
    @SVector [
        x[4],
        x[5],
        x[6],
        P.A*P.a*(sigm(x[2] - x[3], P)) - 2P.a*x[4] - P.a*P.a*x[1],
        P.A*P.a*(μy(t, P) + C2(P)*sigm(C1(P)*x[1], P)) - 2P.a*x[5] - P.a*P.a*x[2],
        P.B*P.b*(C4(P)*sigm(C3(P)*x[1], P)) - 2P.b*x[6] - P.b*P.b*x[3]]
end

function σ(t, x, P::JansenRit)
    @SMatrix [
        0.0;
        0.0;
        0.0;
        0.0;
        P.σy;
        0.0
    ]
end


text:
# Working on figuring out how to prevent rounding error in transit
# times.

# Compute in BigFloat:

#t0 =  4000.0
t0 =  -200.0
h  = 0.07
tmax = 400.0

time = t0
time_big = big(t0)
# Next use compensated summation:
time_comp = t0
# Utilize the notation from Kahan (1965):
s2 = 0.0

times = [time]
times_big = [time_big]
times_comp = [time_comp]
hbig = big(h)

while time < t0+tmax
  # Float64 times:
  time += h/16
  push!(times,time)
  # BigFloat times:
  time_big += hbig/16
  push!(times_big,time_big)
  # Kahan (1965) compensated summation algorithm:
  s2 += h/16
  tmp = time_comp + s2
  s2 = (time_comp - tmp) + s2
  time_comp = tmp
  push!(times_comp,time_comp)
end

using PyPlot

plot(times,times-convert(Array{Float64,1},times_big))
diff_double = maximum(abs.(times-convert(Array{Float64,1},times_big)))
plot(times,times_comp-convert(Array{Float64,1},times_big))
diff_comp = maximum(abs.(times_comp-convert(Array{Float64,1},times_big)))
println("Maximum difference with Float64 precision: ",diff_double)
println("Maximum difference with compensated summation: ",diff_comp)

max_forks_repo_name: langfzac/nbg-papers


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 46f5b944-052e-11eb-0b37-d5300eed0dcf
begin
	using Pkg
	Pkg.add("PlutoUI")
	Pkg.add("Formatting")
	using PlutoUI
	using CSV
	using DataFrames
	using Formatting
end

# ╔═╡ ec894576-02f2-11eb-35a1-b54a763d5020
md"# 🧬 Length of the Unrolled Human DNA

There are approx. 30 trillion cells in the human body. Every cell contains the DNA in its nucleus ([except some](https://wtamu.edu/~cbaird/sq/2013/08/22/why-does-every-cell-in-our-body-contain-dna/)). The unrolled DNA is around 2 meters long 🤯. Now, that is one hell of compaction.

Nature is the ultimate computer scientist.

"

# ╔═╡ ea2d2a0e-02f2-11eb-3db2-53747e1c5fef
# We setup some data in this block
begin
	num_cells = 70e12
	dna_length_metres = 2
	dna_length_km = 2/1000
	human_total_dna_km = dna_length_km * num_cells
	world_population = 7.8e12
end

# ╔═╡ d9c30548-02f3-11eb-0a52-47a9743386fb
md"The average number of cells in the human body is $(format(num_cells, commas=true))"

# ╔═╡ f6f903ec-02f3-11eb-1369-adddbbc08448
md"Total Human DNA ~ DNA Length x Cells in human body:  $(format(human_total_dna_km, commas=true)) KM"

# ╔═╡ a794ceae-0563-11eb-216a-53c433d9efe5
md"\"Amazing isn't it\", I said and expected it to end.

But, share this fact with the kids and see if it ends there.

Almost immediately, the questions will start to roll in....
"

# ╔═╡ 22321d50-02f4-11eb-37e4-9919d5e2d8f2
md"#### Q: If we unrolled the DNA from every cell of the living human and joined together, would it reach the Andromeda galaxy?
"

# ╔═╡ 06b7da50-0552-11eb-174e-2d8e1bd5fe8e
md"[Andromeda](https://en.wikipedia.org/wiki/Andromeda_Galaxy) is about 2.537 million light years from our home planet (and is on the [collision path](https://en.wikipedia.org/wiki/Andromeda%E2%80%93Milky_Way_collision) with our galaxy, in some 4.5 billion years)."

# ╔═╡ 3aeb5032-053f-11eb-2946-a7b9007195d4
md"
![](https://apod.nasa.gov/apod/image/1407/m31_bers_1824.jpg)
Image: NASA
"

# ╔═╡ e142105e-0300-11eb-347a-3b7b00da92fe
md"Let's find out."

# ╔═╡ bb038ac8-02fe-11eb-26a4-2d9291eed60f
begin
	light_year_km=9.461e12
	andromeda_distance=2.537e6
	andromeda_distance_km = light_year_km * andromeda_distance
end	;

# ╔═╡ a2cedd56-054f-11eb-0fa9-ed4a9cd82a94
md"Is the distance between Earth and Andromeda greater than the length of DNA from a single human?"

# ╔═╡ 2b3b24e0-02ff-11eb-2acb-af4c5180971f
andromeda_distance_km > human_total_dna_km

# ╔═╡ a87b8686-02f3-11eb-2a4b-156f41effb1d
mankind_dna_km = world_population * human_total_dna_km;

# ╔═╡ 6a8b0d82-02f4-11eb-37ad-9d1034b9b61f
md"Is the total Human DNA from all living Humans long enough to reach the Andromeda galaxy ?

Answer: $(mankind_dna_km >= andromeda_distance_km)"

# ╔═╡ ea38596e-02f8-11eb-0be4-836ed35af793
md"In fact, it is **$(format(round(mankind_dna_km/andromeda_distance_km), commas=true))** times the distance of Andromeda."

# ╔═╡ 7218783a-0538-11eb-2450-d7a9cfd7d91e
md"---"

# ╔═╡ db135b2a-054e-11eb-34c5-9bdd1526d5f4
md"Alright, so we know that the sum total length of every living human DNA is quite large. Some might say its Bigly."

# ╔═╡ 85e33882-053f-11eb-3893-85e54bfb6008
md"## Milky Way

#### Q: Is it longer than the length of our Galaxy?
"

# ╔═╡ 288b1584-053e-11eb-3bb4-898304daa9e2
md"
![](https://asd.gsfc.nasa.gov/blueshift/wp-content/uploads/2015/07/eso0932a.jpg)
Image: NASA
"

# ╔═╡ 77f90bde-053f-11eb-33aa-1d0485c98c9e
md"## Observable Universe

#### Q: It surely is going to be a lot less than the Universe?
"

# ╔═╡ 642f67b0-053f-11eb-1c04-41eab0f26fff
md"
![](https://upload.wikimedia.org/wikipedia/commons/thumb/d/db/Observable_Universe_Logarithmic_Map_%28horizontal_layout_english_annotations%29.png/2250px-Observable_Universe_Logarithmic_Map_%28horizontal_layout_english_annotations%29.png)
Image: Wikipedia
"

# ╔═╡ 9025091a-053f-11eb-38c0-a1754bcc23bf
md"## Extras

At this point, the tiny inbestigators started asking all sorts of question such as what if you only take the population of Maldives or India. How far will it take us?
"

# ╔═╡ a5aefac2-0565-11eb-29a5-45c409441b18
md"So, I gave them this. Go Play."

# ╔═╡ c3d2c5fa-0534-11eb-1ce2-f11fc7384a7c
md"Total length of DNA from the population of"

# ╔═╡ 6ab56e86-05ee-11eb-3a30-772a9e9a28ae
md"*Note: The interactive version will only work in the Pluto notebook*"

# ╔═╡ eed77dfa-053d-11eb-3ee0-6dc3cb973fbe
md"---"

# ╔═╡ dd33494e-053d-11eb-0c78-9f9ce7c643fe
md"**Setup Section**"

# ╔═╡ 50a55d3c-052e-11eb-0e94-f9562da4db97
# Datasource: https://data.worldbank.org/indicator/SP.POP.TOTL
world_pop = DataFrame(CSV.File("world_population_2019.csv"));

# ╔═╡ 888e99a4-0531-11eb-3bc3-8d847053be75
@bind country Select(vcat(world_pop[:,"Country Name"], "1 Human"))

# ╔═╡ 20df1f62-0534-11eb-275b-4b15029abe44
function pop_to_dna(population)
	 return population * human_total_dna_km
end

# ╔═╡ 7a16cfac-0533-11eb-3868-c36178506bc3
planet_distances_km = DataFrame(CSV.File("earth_distance.csv"));

# ╔═╡ b0a7118c-0534-11eb-04cb-3f273d884de3
@bind planet Select(planet_distances_km[:, "Name"])

# ╔═╡ 1fcf05c4-0535-11eb-2a29-fdd822b3f685
begin
	if country == "1 Human"
		current_pop = 1
	else
		idx = findfirst(x -> x == country, world_pop[:, "Country Name"])
		current_pop = world_pop[idx,"2019"]
	end
	planet_idx = findfirst(x -> x == planet, planet_distances_km[:, "Name"])
	planet_distance = planet_distances_km[planet_idx, "Distance_KM"]
	distance_times = pop_to_dna(current_pop)/planet_distance
end;

# ╔═╡ d0c06a54-0534-11eb-2f55-134fe93a4e20
md" is $(format(distance_times, commas=true)) times the distance from Earth to "

# ╔═╡ a2caf1bc-0531-11eb-28ac-41ab5b064e33
md"Population of $(uppercase(country)) is $(current_pop)"

# ╔═╡ f90bb6ce-0533-11eb-02b5-95ddddb853d6
# Width of Milky Way
milkyway_width_km = 105700 * light_year_km

# ╔═╡ 042f4fc0-02f5-11eb-05b5-b5452bc3d81c
md"

The width of Milky Way is **$(format(milkyway_width_km, commas=true))** KM.

We could go from one end of the Milky way to another $(format(mankind_dna_km/milkyway_width_km, commas=true)) times, for that amount of DNA.
"

# ╔═╡ 49e46e4a-0536-11eb-1d2a-cfcb1dbbf9fc
observable_universe_km = 93e9 * light_year_km

# ╔═╡ 0744796a-0539-11eb-0f1c-6bbcd56c3e02
md"The observable Universe is $(observable_universe_km) KM. And it is approximately the same length as the total length of entire world population DNA."

# ╔═╡ e6d30c24-053c-11eb-0226-b94f320cdeb1
md"Approximately ~ $(format( mankind_dna_km/observable_universe_km, commas=true )) times"

# ╔═╡ Cell order:
# ╟─ec894576-02f2-11eb-35a1-b54a763d5020
# ╠═ea2d2a0e-02f2-11eb-3db2-53747e1c5fef
# ╟─d9c30548-02f3-11eb-0a52-47a9743386fb
# ╟─f6f903ec-02f3-11eb-1369-adddbbc08448
# ╟─a794ceae-0563-11eb-216a-53c433d9efe5
# ╟─22321d50-02f4-11eb-37e4-9919d5e2d8f2
# ╟─06b7da50-0552-11eb-174e-2d8e1bd5fe8e
# ╟─3aeb5032-053f-11eb-2946-a7b9007195d4
# ╟─e142105e-0300-11eb-347a-3b7b00da92fe
# ╟─bb038ac8-02fe-11eb-26a4-2d9291eed60f
# ╟─a2cedd56-054f-11eb-0fa9-ed4a9cd82a94
# ╠═2b3b24e0-02ff-11eb-2acb-af4c5180971f
# ╟─a87b8686-02f3-11eb-2a4b-156f41effb1d
# ╟─6a8b0d82-02f4-11eb-37ad-9d1034b9b61f
# ╟─ea38596e-02f8-11eb-0be4-836ed35af793
# ╟─7218783a-0538-11eb-2450-d7a9cfd7d91e
# ╟─db135b2a-054e-11eb-34c5-9bdd1526d5f4
# ╟─85e33882-053f-11eb-3893-85e54bfb6008
# ╟─288b1584-053e-11eb-3bb4-898304daa9e2
# ╟─042f4fc0-02f5-11eb-05b5-b5452bc3d81c
# ╟─77f90bde-053f-11eb-33aa-1d0485c98c9e
# ╟─642f67b0-053f-11eb-1c04-41eab0f26fff
# ╟─0744796a-0539-11eb-0f1c-6bbcd56c3e02
# ╟─e6d30c24-053c-11eb-0226-b94f320cdeb1
# ╟─9025091a-053f-11eb-38c0-a1754bcc23bf
# ╟─a5aefac2-0565-11eb-29a5-45c409441b18
# ╟─c3d2c5fa-0534-11eb-1ce2-f11fc7384a7c
# ╟─888e99a4-0531-11eb-3bc3-8d847053be75
# ╟─d0c06a54-0534-11eb-2f55-134fe93a4e20
# ╟─b0a7118c-0534-11eb-04cb-3f273d884de3
# ╟─a2caf1bc-0531-11eb-28ac-41ab5b064e33
# ╟─6ab56e86-05ee-11eb-3a30-772a9e9a28ae
# ╟─eed77dfa-053d-11eb-3ee0-6dc3cb973fbe
# ╟─dd33494e-053d-11eb-0c78-9f9ce7c643fe
# ╠═46f5b944-052e-11eb-0b37-d5300eed0dcf
# ╠═50a55d3c-052e-11eb-0e94-f9562da4db97
# ╠═1fcf05c4-0535-11eb-2a29-fdd822b3f685
# ╠═20df1f62-0534-11eb-275b-4b15029abe44
# ╠═7a16cfac-0533-11eb-3868-c36178506bc3
# ╠═f90bb6ce-0533-11eb-02b5-95ddddb853d6
# ╠═49e46e4a-0536-11eb-1d2a-cfcb1dbbf9fc

max_stars_repo_path: Cells vs DNA/cells_vs_dna.jl

max_stars_repo_name: katnoria/data-books

max_issues_repo_path: Cells vs DNA/cells_vs_dna.jl

max_issues_repo_name: katnoria/data-books

max_forks_repo_path: Cells vs DNA/cells_vs_dna.jl

max_forks_repo_name: katnoria/data-books


text:
abstract type SearchDirection end
struct SearchBackward <: SearchDirection end
struct SearchForward <: SearchDirection end
struct SearchBoth <: SearchDirection end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
function info_analysis(info::Dict)
    D = get(info, :tree, nothing)
    if D !== nothing
        show(stdout, MIME("text/plain"), D)
        println("Number of action node expanded:", D.ba)
        println("Number of belief node expanded:", D.b)
    end
    depth = info[:depth]
    println("Times of exploration: ", length(depth))
    println(@sprintf("Depth of exploration: mean±std = %5.2f±%4.2f", mean(depth), std(depth)))
    println("90% Confidence interval = ", quantile(depth, (0.05, 0.95)))
    return nothing
end

function hist_analysis(hist::H; display_mean_and_std::Bool = false, layout=(1,2), font_size=12, margin=40px, figure_size=(900,400)) where H<:AbstractSimHistory
    infos = ainfo_hist(hist)

    median_d = Float64[]
    lower_d = Float64[]
    upper_d = Float64[]
    mean_d = Float64[]
    std_d = Float64[]

    for info in infos
        depth = info[:depth]
        l_d, m_d, u_d = quantile(depth, (0.05, 0.5, 0.95))
        push!(median_d, m_d)
        push!(lower_d, m_d-l_d)
        push!(upper_d, u_d-m_d)
        push!(mean_d, mean(depth))
        push!(std_d, std(depth))
    end
    if display_mean_and_std
        p1 = plot(median_d, ribbon=(lower_d, upper_d), xaxis="Steps", yaxis="Depth of Exploration", label="quantile", xtickfontsize=font_size, ytickfontsize=font_size, xguidefontsize=font_size, yguidefontsize=font_size, legendfontsize=font_size)
        plot!(p1, mean_d, ribbon=std_d, label="mean", xtickfontsize=font_size, ytickfontsize=font_size, xguidefontsize=font_size, yguidefontsize=font_size, legendfontsize=font_size)
    else
        p1 = plot(median_d, ribbon=(lower_d, upper_d), xaxis="Steps", yaxis="Depth of Exploration", legend=false, xtickfontsize=font_size, ytickfontsize=font_size, xguidefontsize=font_size, yguidefontsize=font_size, legendfontsize=font_size)
    end

    D = get(first(infos), :tree, nothing)
    if D === nothing
        display(p1)
    else
        num_anode = Int[]
        for info in infos
            D = info[:tree]
            push!(num_anode, D.ba)
        end
        base = 10^floor(log(10, mean(num_anode)))
        p2 = plot(num_anode, legend=false, xaxis="Steps", yaxis="Action Nodes Expanded", yformatter=y->y/base, xtickfontsize=font_size, ytickfontsize=font_size, xguidefontsize=font_size, yguidefontsize=font_size, legendfontsize=font_size)
        annotate!(p2, [(0.9, maximum(num_anode) + (maximum(num_anode)-minimum(num_anode)) * 0.07, Plots.text(@sprintf("\$\\times10^{%d}\$", round(Int, log(10,base))), font_size, :black, :center, "courier"))])
        display(plot(p1, p2, layout=layout, size=figure_size, margin=margin))
    end
    return nothing
end


text:
if length(ARGS) == 1
  global experiment_name = ARGS[1]
else
  println("Please call with one argument (the experiment name)")
end

println("Processing: ", experiment_name)

isjulia(s) = endswith(s, ".jld2")
is_experiment_file(s) = startswith(s, experiment_name^2) && isjulia(s)

trim_name(s) = s[length(experiment_name)+1:end]

if !isdir(experiment_name)
  println("Making experiment directory: ", experiment_name)
  mkdir(experiment_name)
end

for file in readdir()
  if is_experiment_file(file)
    dest = string(experiment_name, "/", trim_name(file))

    println("Old: ", file, "\nNew: ", dest)

    mv(file, dest)
  end
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
struct PSAModel{T} <: AbstractSampler
    nbits::Int
    w::Vector{T}
end

function PSAModel(nbits::Int)
    w = rand(nbits) .- 0.5
    PSAModel(nbits, w)
end

function get_logp(model::PSAModel, x::AbstractArray)
    xhat = sigmoid.(model.w)
    ϵ = 1E-7
  sum(x.*log.(xhat.+ϵ) + (1 .-x).*log.(1.0 .- xhat .+ ϵ), dims=1)
end

function gen_samples(model::PSAModel, nbatch::Int)
    x = zeros(Int, model.nbits, nbatch)
    for i = 1:model.nbits
        xhat = sigmoid.(model.w)
        x[i, :] .= xhat[i,:] .> rand(nbatch)
    end
    x
end

model_parameters(model::PSAModel) = (model.w, )

function model_dispatch!(model::PSAModel, θ)
    model.w .= θ[1]
end


text:
include("../oracles/assignment_oracle.jl")

using Test

@testset "Assignment Oracle" begin
    C = [2.0 6; 3 5]
    c = reshape(C, 4)
    (z, w) = assignment_oracle(c)
    W = reshape(w, (2,2))
    @test z ≈ 7.0
    @test W ≈ [1 0; 0 1]

    c_ham = assignment_hamming_reduce(c)
    (z2, w2) = assignment_oracle(c_ham)
    W2 = reshape(w2, (2,2))
    @test z2 ≈ 0
    @test W2 ≈ W
end

max_stars_repo_name: paulgrigas/SmartPredictThenOptimize

max_issues_repo_name: paulgrigas/SmartPredictThenOptimize

max_forks_repo_name: paulgrigas/SmartPredictThenOptimize


text:
using Documenter, Ledger

makedocs(;
    modules=[Ledger],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/alecloudenback/Ledger.jl/blob/{commit}{path}#L{line}",
    sitename="Ledger.jl",
    authors="Alec Loudenback",
    assets=String[],
)

max_stars_repo_name: alecloudenback/Ledger

max_issues_repo_name: alecloudenback/Ledger

max_forks_repo_name: alecloudenback/Ledger


text:
# The concept and trivial implementation of graph visitors

abstract AbstractGraphVisitor

# trivial implementation

# invoked when a vertex v is encountered for the first time
# this function returns whether to continue search
discover_vertex!(vis::AbstractGraphVisitor, v) = true

# invoked when the algorithm is about to examine v's neighbors
open_vertex!(vis::AbstractGraphVisitor, v) = nothing

# invoked when a neighbor is discovered & examined
examine_neighbor!(vis::AbstractGraphVisitor, u, v, color::Int, ecolor::Int) = nothing

# invoked when an edge is discovered & examined
examine_edge!(vis::AbstractGraphVisitor, e, color::Int) = nothing

# invoked when all of v's neighbors have been examined
close_vertex!(vis::AbstractGraphVisitor, v) = nothing

type TrivialGraphVisitor <: AbstractGraphVisitor
end

# This is the common base for BreadthFirst and DepthFirst
abstract AbstractGraphVisitAlgorithm

###########################################################
#
#   General algorithms based on graph traversal
#
###########################################################

# List vertices by the order of being discovered

type VertexListVisitor{V} <: AbstractGraphVisitor
    vertices::Vector{V}

    function VertexListVisitor(n::Integer)
        vs = Array(V, 0)
        sizehint(vs, n)
        new(vs)
    end
end

function discover_vertex!{V}(visitor::VertexListVisitor{V}, v::V)
    push!(visitor.vertices, v)
    true
end

function visited_vertices{V,E}(
    graph::AbstractGraph{V,E},
    alg::AbstractGraphVisitAlgorithm,
    sources)

    visitor = VertexListVisitor{V}(num_vertices(graph))
    traverse_graph(graph, alg, sources, visitor)
    visitor.vertices::Vector{V}
end

# Print visit log

type LogGraphVisitor{S<:IO} <: AbstractGraphVisitor
    io::S
end

function discover_vertex!(vis::LogGraphVisitor, v)
    println(vis.io, "discover vertex: $v")
    true
end

open_vertex!(vis::LogGraphVisitor, v) = println(vis.io, "open vertex: $v")
close_vertex!(vis::LogGraphVisitor, v) = println(vis.io, "close vertex: $v")

function examine_neighbor!(vis::LogGraphVisitor, u, v, vcolor::Int, ecolor::Int)
    println(vis.io, "examine neighbor: $u -> $v (vertexcolor = $vcolor, edgecolor= $ecolor)")
end

function examine_edge!(vis::LogGraphVisitor, e, color::Int)
    println(vis.io, "examine edge: $e")
end

function traverse_graph_withlog(g::AbstractGraph, alg::AbstractGraphVisitAlgorithm, sources, io::IO)
    visitor = LogGraphVisitor(io)
    traverse_graph(g, alg, sources, visitor)
end

traverse_graph_withlog(g::AbstractGraph, alg::AbstractGraphVisitAlgorithm,
    sources) = traverse_graph_withlog(g, alg, sources, STDOUT)


text:
import Base: @__doc__

"""
# Summary

abstract type HELICS.CWrapper
"""
abstract type CWrapper end

Base.convert(T::Type{<:CWrapper}, p::Ptr{Nothing}) = T(p)
Base.unsafe_convert(T::Type{Ptr{Nothing}}, t::CWrapper) = t.ptr

"""
# Summary

abstract type HELICS.Federate <: HELICS.CWrapper

# Subtypes

- [`HELICS.CombinationFederate`](@ref)
- [`HELICS.MessageFederate`](@ref)
- [`HELICS.ValueFederate`](@ref)

# Supertype Hierarchy

HELICS.Federate <: HELICS.CWrapper <: Any

"""
abstract type Federate <: CWrapper end

macro define(subtype, supertype)
    docstring = """
# Summary

struct [`HELICS.$subtype`](@ref) <: [`HELICS.$supertype`](@ref)
"""
    esc(quote
        struct $subtype <: $supertype
            ptr::Ptr{Nothing}
            function $subtype(ptr::Ptr{Nothing})
                ptr == C_NULL && error("Failed to create $($subtype). Received null pointer from HELICS C interface.")
                new(ptr)
            end
        end
        @doc $docstring $subtype
    end)
end

@define Broker CWrapper
@define Core CWrapper
@define FederateInfo CWrapper
@define ValueFederate Federate
@define MessageFederate Federate
@define CombinationFederate Federate
@define Publication CWrapper
@define Subscription CWrapper
@define Endpoint CWrapper
@define Filter CWrapper
@define Query CWrapper

@define Message CWrapper

const Input = Subscription


text:
mutable struct TraitLoggerXMLElement <: MyXMLElement
    el::XMLOrNothing
    treexml::TreeModelXMLElement
    traitxml::TraitLikelihoodXMLElement

    function TraitLoggerXMLElement(treexml::TreeModelXMLElement,
                                traitxml::TraitLikelihoodXMLElement)
        return new(nothing, treexml, traitxml)
    end

end

function make_xml(xml::TraitLoggerXMLElement)
    tree_el = make_xml(xml.treexml)
    trait_el = make_xml(xml.traitxml)

    el = new_element(bn.TRAIT_LOGGER)
    trait_name = attribute(trait_el, bn.TRAIT_NAME)
    set_attributes(el, [bn.ID => bn.TRAIT_LOGGER,
                    bn.TRAIT_NAME => trait_name,
                    bn.TAXON_EXPLICIT => bn.TRUE,
                    bn.NODES => bn.EXTERNAL])

    add_ref_el(el, tree_el)
    add_ref_el(el, trait_el)
    xml.el = el
    return el
end

function get_loggables(xml::TraitLoggerXMLElement)
    make_xml(xml)
    return [xml.el]
end


text:
using Quon
using Test

@testset "yang_baxter.jl" begin
    include("utils/yang_baxter.jl")
end

@testset "quon_const.jl" begin
    include("quon_const.jl")
end

@testset "yang-baxter-singular.jl" begin
    include("yang-baxter-singular.jl")
end

@testset "graphs/tait.jl" begin
    include("graphs/tait.jl")
end

@testset "plots.jl" begin
    include("plots.jl")
end

@testset "rules.jl" begin
    include("rules.jl")
end


text:
export QFTBlock

struct QFTBlock{N} <: PrimitiveBlock{N} end
mat(::Type{T}, q::QFTBlock{N}) where {T, N} = T.(applymatrix(q))

apply!(reg::DefaultRegister{B}, ::QFTBlock) where B = (reg.state = ifft!(invorder_firstdim(reg |> state), 1)*sqrt(1<<nactive(reg)); reg)
apply!(reg::DefaultRegister{B}, ::Daggered{N, <:QFTBlock}) where {B,N} = (reg.state = invorder_firstdim(fft!(reg|>state, 1)/sqrt(1<<nactive(reg))); reg)

# traits
ishermitian(q::QFTBlock{N}) where N = N==1
isreflexive(q::QFTBlock{N}) where N = N==1
isunitary(q::QFTBlock{N}) where N = true

function print_block(io::IO, pb::QFTBlock{N}) where N
    printstyled(io, "QFT(1-$N)"; bold=true, color=:blue)
end

function print_block(io::IO, pb::Daggered{N,<:QFTBlock}) where {N, T}
    printstyled(io, "IQFT(1-$N)"; bold=true, color=:blue)
end