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

name = "SuperLU_MT"
version = v"3.1"

# Collection of sources required to build Elemental
sources = [
    ArchiveSource("https://portal.nersc.gov/project/sparse/superlu/superlu_mt_3.1.tar.gz",
                  "407b544b9a92b2ed536b1e713e80f986824cf3016657a4bfc2f3e7d2a76ecab6"),
    DirectorySource("bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd SuperLU_MT_*

cp MAKE_INC/make.linux.openmp make.inc

atomic_patch -p1 ../patches/01-fix-makefile.patch

# Build object files suitable for a shared library, and link the libraries we need.
echo "CFLAGS += -fPIC -fopenmp \$(BLASLIB)" >> make.inc
echo "NOOPTS += -fPIC -fopenmp \$(BLASLIB)" >> make.inc

# We need to add a suffix to BLAS symbols.
SYMBOLS=()
for fun in isamax sasum saxpy scopy strsv idamax dasum daxpy dcopy dtrsv ctrsv ztrsv; do
  SYMBOLS+=("-D$fun=${fun}_")
done
echo "CFLAGS += ${SYMBOLS[@]}" >> make.inc

# Don't use 64-bit integers on non-64-bit systems.
if [[ "$nbits" != 64 ]]; then
  sed -i "s/^CFLAGS.*+= -D_LONGINT$//" make.inc
fi

# Clang doesn't play nicely with OpenMP.
if [[ "$target" == *-freebsd* || "$target" == *-apple-* ]]; then
  CC=gcc
fi

# Weird sed delimiters because some variables contain slashes.
sed -i "s~^BLASLIB.*~BLASLIB = -L$libdir -lopenblas~" make.inc
sed -i "s~^CC.*~CC = $CC~" make.inc
sed -i "s~^FORTRAN.*~FORTRAN = $FC~" make.inc

make superlulib "-j$nproc"
cp lib/single "$libdir/libsuperlumts.$dlext"
cp lib/double "$libdir/libsuperlumtd.$dlext"
cp lib/complex "$libdir/libsuperlumtc.$dlext"
cp lib/complex16 "$libdir/libsuperlumtz.$dlext"
cp SRC/*.h "$prefix/include"
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libsuperlumts", :libsuperlumts),
    LibraryProduct("libsuperlumtd", :libsuperlumtd),
    LibraryProduct("libsuperlumtc", :libsuperlumtc),
    LibraryProduct("libsuperlumtz", :libsuperlumtz),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("OpenBLAS32_jll"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#__precompile__(true)
## Code to make writing of d3 javascript d3-like.

module DThree

import Base.getindex, Base.push!

using JSON
using Mustache
using Blink

export asis, @asis_str, D3, browse
export js, loadurl, loadfile

type AsIs
    x
end
asis(x) = AsIs(x)
macro asis_str(x)
    asis(x)
end

## D3 instances can chain commands with interface nearly same as d3.js
## d3.selectAll("p").style("color","#00f").render()
## render passes JavaScript to browser. If you want to assign to variable name, pass as argument.
## If method not defined, then can do something like d3._(:selectAll, "p"). ...
## by default this uses JavaScript d3 object as receiver. This can be changed with d3.receiver("chart"). ...
## arguments are converted via to_json, except:
## * functions are treated as callbacks into julia. These are asynchronous.
## * use asis"x" or asis("x") to treat object "as is". This is needed to quote JavaScript functions
## TODO: PyCall this baby so members can be added from a list of symbols
type D3
    cmd
    _var
    var
    receiver
    eval_js
    render
    _
    __
    select
    selectAll
    data
    datum
    enter
    append
    style
    color
    attr
    text
    domain
    range
    rangeRoundBounds
    scale_linear
    scale_ordinal
    scale_category10
    time_scale
    svg_axis
    svg_line
    svg_area
    scale
    ticks
    tickFormat
    format
    axisLabel
    orient
    call
    max
    transition
    duration
    function D3()
        this = new("d3", nothing)

        tmp = tempname() * ".html"

        ## for setting var nm = ...
        this.var      = (value::Union{Void, AbstractString}) -> begin this._var = value; this end

        ## for setting reciever.meth1(...).meth2. ... where default is d3.
        this.receiver = (value::AbstractString) -> begin this.cmd = value; this end

        ## return JavaScript as string, reset cmd
        this.render = () -> begin
            out = this.cmd
            if !isa(this._var, Void)
                out = "var $(this._var) = " * out * ";"
            end
            this.receiver("d3")
            this.var(nothing)
            return(out)
        end
        this._ = (meth, args...) ->  begin
            cmd = this.cmd
            meth = replace(string(meth), r"_", ".")
            args = map(u -> prep(u), args)
            args = join(args, ", ")
            this.cmd = "$(this.cmd).$meth($args)"
            this
        end
        this.__ = (prop) -> begin # add a property
            this.cmd = "$(this.cmd).$prop"
            this
        end
        ## generate in @eval map meth by replacing _ with symbol(replace("scale_linear", "_", "."))
        this.select    = (args...) -> this._(:select, args...)
        this.selectAll = (args...) -> this._(:selectAll, args...)
        this.data      = (args...) -> this._(:data, args...)
        this.datum     = (args...) -> this._(:datum, args...)
        this.enter     = (args...) -> this._(:enter, args...)
        this.append    = (args...) -> this._(:append, args...)
        this.style     = (args...) -> this._(:style, args...)
        this.color     = (args...) -> this._(:color, args...)
        this.attr      = (args...) -> this._(:attr, args...)
        this.text      = (args...) -> this._(:text, args...)
        this.domain    = (args...) -> this._(:domain, args...)
        this.range     = (args...) -> this._(:range, args...)
        this.rangeRoundBounds     = (args...) -> this._(:rangeRoundBands, args...)
        this.scale_linear = (args...) -> this._("scale.linear", args...)
        this.scale_ordinal = (args...) -> this._("scale.ordinal", args...)
        this.scale_category10 = (args...) -> this._("scale.category10", args...)
        this.time_scale= (args...) -> this._("time.scale", args...)
        this.svg_axis  = (args...) -> this._("svg.axis", args...)
        this.svg_line  = (args...) -> this._("svg.line", args...)
        this.svg_area  = (args...) -> this._("svg.area", args...)
        this.scale     = (args...) -> this._(:scale, args...)
        this.ticks     = (args...) -> this._(:ticks, args...)
        this.tickFormat= (args...) -> this._(:tickFormat, args...)
        this.format    = (args...) -> this._(:format, args...)
        this.axisLabel= (args...) -> this._(:axisLabel, args...)
        this.orient    = (args...) -> this._(:orient, args...)
        this.call      = (args...) -> this._(:call, args...)
        this.max       = (args...) -> this._(:max, args...)
        this.transition= (args...) -> this._(:transition, args...)
        this.duration  = (args...) -> this._(:duration, args...)

        this
    end
end

## used for easier? access D3()[:one](2)[:three]("three") ...
getindex(x::D3, i::Symbol) = (args...) -> x._(i, args...)

prep(x::Any) = json(x)
prep(x::AsIs) = string(x.x)
## NEED prep(d3, x::Function) to callback into `julia`
prep(x::D3) = prep(asis(x.render()))

include("colors.jl")
include("display.jl")
include("blink.jl")
#include("plots.jl")
include("plotly.jl")

end ## module

ext: jl


text:
@doc raw"""
`levenshteinAlign` returns the Levenshtein-optimal pairwise alignment
between two sequences (abstract string or vector) as a two-column
matrix
"""
function levenshteinAlign(s1::T, s2::T) where {T<:Union{AbstractString,Vector}}
    n, m = length.([s1, s2])
    dp = Matrix{Int}(undef, n + 1, m + 1)
    pointer = Matrix{Int}(undef, n + 1, m + 1)
    dp[:, 1] .= 0:n
    dp[1, :] .= 0:m
    pointer[:, 1] .= 1
    pointer[1, :] .= 2
    pointer[1, 1] = -1
    for i = 1:m
        for j = 1:n
            insrt = dp[j, i+1] + 1
            dlt = dp[j+1, i] + 1
            mtch = dp[j, i] + (s1[j] != s2[i])
            states = [mtch, insrt, dlt]
            dp[j+1, i+1] = minimum(states)
            pointer[j+1, i+1] = argmin(states) - 1
        end
    end
    a1 = Union{eltype(s1), Missing}[]
    a2 = Union{eltype(s2), Missing}[]
    i, j = m + 1, n + 1
    while pointer[j, i] >= 0
        p = pointer[j, i]
        if p == 0
            pushfirst!(a1, s1[j-1])
            pushfirst!(a2, s2[i-1])
            i -= 1
            j -= 1
        elseif p == 1
            pushfirst!(a1, s1[j-1])
            pushfirst!(a2, missing)
            j -= 1
        else
            pushfirst!(a1, missing)
            pushfirst!(a2, s2[i-1])
            i -= 1
        end
    end
    hcat(a1, a2)
end

@doc raw"""
`levenshteinDistance` computes the edit distance between two iterables
(abstract string or vector)
"""
function levenshteinDistance(
    s1::T,
    s2::T,
) where {T<:Union{AbstractString,Vector}}
    n, m = length.([s1, s2])
    dp = Matrix{Int}(undef, n + 1, m + 1)
    dp[:, 1] .= 0:n
    dp[1, :] .= 0:m
    for i = 1:m
        for j = 1:n
            insrt = dp[j, i+1] + 1
            dlt = dp[j+1, i] + 1
            mtch = dp[j, i] + (s1[j] != s2[i])
            states = [mtch, insrt, dlt]
            dp[j+1, i+1] = minimum(states)
        end
    end
    dp[end, end]
end

@doc raw"""
`ldn` computes the normalized edit distance between two iterables
(abstract string or vector), i.e.
$L(x,y) /max(|x|, |y|)$
"""
function ldn(s1::T, s2::T) where {T<:Union{AbstractString,Vector}}
    levenshteinDistance(s1, s2) / maximum(length.([s1, s2]))
end


text:
# Copyright 2007, Timothy A. Davis, Univ. of Florida
# Copyright 2017, Andrey Prokopenko, ORNL

# Adapted from Tim Davis' meshnd code
function meshnd(m, n, k = 1)
    """MESHND creation and nested dissection of a regular 2D or 3D mesh.

    G = meshnd(m,n)   constructs a m-by-n 2D mesh G
    G = meshnd(m,n,k) constructs a m-by-n-by-k 3D mesh G

    Example:
    G = meshnd(4,5)

    returns
    G =
    1     2     3     4     5
    6     7     8     9    10
    11    12    13    14    15
    16    17    18    19    20

    With no inputs, a few example meshes are generated and plotted.

    See also nested, numgrid.
    """

    if k == 1
        # create the m-by-n-by-k mesh in "natural" (row-major) order.  This is how
        # a typical 2D mesh is ordered.  A column-major order would be better, since
        # in that case G(:) would equal 1:(m*n) ... but let's stick with tradition.

        # FIXME: this does not work as expected
        # G = reshape(1:(m*n*k), n, m, k)'
        G = zeros(Int64, m, n, k)
        for ii = 1:m
            for jj = 1:n
                G[ii, jj, 1] = (ii-1)*m + jj
            end
        end
    else
        # create the m-by-n-by-k mesh in column-major order. The first m-by-n-by-1
        # slice is in column-major order, followed by all the other slices 2 to k.
        G = reshape(1:(m*n*k), m, n, k)
    end

    return G

end

max_stars_repo_name: aprokop/ca-preconditioners

max_issues_repo_name: aprokop/ca-preconditioners

max_forks_repo_name: aprokop/ca-preconditioners


text:
"""
$(TYPEDEF)

Stores information on requested and measured [`Configuration`](@ref)s.

$(TYPEDFIELDS)
"""
const Measurement = NamedTuple{
    (
        :measured,
        :cost,
        :technique,
        :date_requested,
        :date_measured
    ),
    Tuple{
        Bool,
        Union{Float64, Missing},
        Union{String, Missing},
        Union{DateTime, Missing},
        Union{DateTime, Missing}
    }
}

Measurement() = (
    measured = false,
    cost = missing::Union{Float64, Missing},
    technique = missing::Union{String, Missing},
    date_requested = now()::Union{DateTime, Missing},
    date_measured = missing::Union{DateTime, Missing}
)

"""
$(TYPEDSIGNATURES)

Pushes a new [`Configuration`](@ref) to a `DataFrame`, appending an
empty [`Measurement`](@ref).
"""
function push!(target_table::DataFrame, c::Configuration)
    push!(
        target_table,
        merge(
            Measurement(),
            NamedTuple{keys(c)}(
                (p.current for p in c.parameters)
            )
        )
    )
end

"""
$(TYPEDSIGNATURES)

Creates   an    empty   `DataFrame`   with    column   types   from    a   given
[`Configuration`](@ref) and from an empty [`Measurement`](@ref).
"""
function configuration_table(c::Configuration)
    push!(
        DataFrame(
            NamedTuple{
                (keys(c)..., fieldnames(Measurement)...)
            }(
                Array{T, 1}() for T in ((typeof(p.current) for p in values(c))...,
                                        fieldtypes(Measurement)...)
            )
        ),
        c
    )
end


text:
using Documenter
using GitForge

makedocs(;
    modules=[GitForge],
    format=Documenter.HTML(),
    pages=["Home" => "index.md"],
    repo="https://github.com/christopher-dG/GitForge.jl/blob/{commit}{path}#L{line}",
    sitename="GitForge.jl",
    authors="Chris de Graaf",
)

deploydocs(; repo="github.com/christopher-dG/GitForge.jl")


text:
@valueequals struct Rectangle
    x::Union{Int64, Float64}
    y::Union{Int64, Float64}
    w::Union{Int64, Float64}
    h::Union{Int64, Float64}

    Rectangle(x::Integer, y::Integer, w::Integer, h::Integer) = new(Int64(x), Int64(y), Int64(w), Int64(h))
    Rectangle(x, y, w, h) = new(Float64(x), Float64(y), Float64(w), Float64(h))
end

# AABB check
function checkCollision(box1::Rectangle, box2::Rectangle)
    return (
        box1.x < box2.x + box2.w &&
        box1.x + box1.w > box2.x &&
        box1.y < box2.y + box2.h &&
        box1.y + box1.h > box2.y
    )
end

function intersection(box1::Rectangle, box2::Rectangle)
    x = max(box1.x, box2.x)
    y = max(box1.y, box2.y)
    w = min(box1.x + box1.w, box2.x + box2.w) - x
    h = min(box1.y + box1.h, box2.y + box2.h) - y
    return w > 0 && h > 0 ? Rectangle(x, y, w, h) : nothing
end

# Find top left corner and bottom right corner
function bounds(rects::Array{Rectangle, 1})
    tlx = tly = typemax(Int)
    brx = bry = typemin(Int)

    for rect in rects
        tlx = min(tlx, rect.x)
        tly = min(tly, rect.y)

        brx = max(brx, rect.x + rect.w)
        bry = max(bry, rect.y + rect.h)
    end

    return tlx, tly, brx, bry
end

# Rectangle that covers all input rectangles
function coverRectangles(rects::Array{Rectangle, 1})
    tlx, tly, brx, bry = bounds(rects)

    return Rectangle(tlx, tly, brx - tlx, bry - tly)
end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
import GenericInstruments: query

function connect!(address; mode=GI.VI_NO_LOCK, timeout=UInt32(10000))
		# resource manager
		rmg = GI.viOpenDefaultRM()
		#Pointer for the instrument handle
		vi = GI.ViPSession(0)
		GI.check_status(GI.viOpen(rmg, address, mode, timeout, vi))
		handle = vi.x
end

function disconnect!(handle)
	GI.viClose(handle)
end

function write(handle, cmd::String)
	GI.viWrite(handle, cmd*"\n")
	return nothing
end

function read(handle, bufSize::UInt32=0x00000100)
	strip(GI.viRead(handle; bufSize=bufSize), ['\r', '\n'])
end

function query(handle, cmd::String, delay::Real=0.0)
    #flush read buffer
    #GenericInstruments.viFlush(psu.handle, GenericInstruments.VI_READ_BUF_DISCARD)
	write(handle, cmd)
	sleep(delay)
	x = read(handle)
	#GenericInstruments.viFlush(psu.handle, GenericInstruments.VI_READ_BUF_DISCARD)
	return x
end

# aliases
function visaWrite(handle, cmd)
	write(handle, cmd)
end

function visaWrite(address::String, cmd)
	handle = connect!(address)
	write(handle, cmd)
	disconnect!(handle)
end

function visaRead(handle, bufSize::UInt32=0x00000100)
	read(handle, bufSize)
end

function visaRead(address::String, bufSize::UInt32=0x00000100)
	handle = connect!(address)
	ro = read(handle, bufSize)
	disconnect!(handle)
	return ro
end

function query(address::String, cmd::String)
	handle = connect!(address)
	ro = ""
	try
		ro = query(handle, cmd)
	catch
		ro = query(handle, cmd)
	end
	disconnect!(handle)
	return ro
end

# Helper functions to find instruments
function find_resources(expr::AbstractString="?*::INSTR")
	resmgr = GenericInstruments.viOpenDefaultRM()
	addresses = GenericInstruments.viFindRsrc(resmgr, expr)
	for address in addresses
		try
			instr_name = query(address::String, "*IDN?")
			println("Found $instr_name on address: $address")
		catch
			nothing
		end
	end
end

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
interpolate_from_dict(ex::Expr, dict) =
    Expr(ex.head, interpolate_from_dict.(ex.args, Ref(dict))...)
interpolate_from_dict(ex::Symbol, dict::Dict) = get(dict, ex, ex)
interpolate_from_dict(ex::Any, dict) = ex

layout_simbols = (
    Symbol(/),
    Symbol(*),
    :vstack,
    :lvstack,
    :leftalign,
    :hstack,
    :center,
    :rightalign,
    :lvstack,
    :rvstack,
    :cvstack,
    :pad,
    :pad!,
    :vertical_pad,
    :vertical_pad!,
)

function collect_elements(exp::Expr)
    if exp.args[1] ∉ layout_simbols && length(exp.args) > 2
        n, w, h = exp.args
        return :($n($w, $h))
    elseif exp.args[1] ∉ layout_simbols
        return nothing
    else
        symbols = map(x -> x isa Symbol ? x : collect_elements(x), exp.args)
        symbols = filter(s -> s ∉ layout_simbols && !isnothing(s), symbols)
        return reduce(vcat, symbols)
    end
end

function clean_layout_symbol(s::Symbol)
    s[1] ∉ layout_simbols ? s[1] : s
end

function clean_layout_expr(exp::Expr)
    if exp.args[1] ∉ layout_simbols
        return exp.args[1]
    else
        exp.args = map(a -> a isa Expr ? clean_layout_expr(a) : a, exp.args)
    end
    return exp
end

placeholder(s, w, h, c) = PlaceHolder(
    w,
    h;
    style = c,
    text = "{bold underline bright_blue}$s{/bold underline bright_blue} {white}($w × $h){/white}",
)

function extract_renderable_from_kwargs(s, w, h; kwargs...)
    ren = get(kwargs, s, nothing)
    if !isnothing(ren)
        ren = ren isa AbstractRenderable ? ren : RenderableText(ren)

        # check renderable has the right size
        msg = "While creating a Compository layout, the layout element"
        @assert ren.measure.w == w highlight(
            msg * " :$s has width $w but the renderable passed has width $(ren.measure.w)",
        )
        @assert ren.measure.h == h highlight(
            msg *
            " :$s has height $h but the renderable passed has height $(ren.measure.h)",
        )
    end
    return ren
end


text:
using Transformers
using Transformers.Basic
using Test

import Flux
using Flux: gradient

const tests = [
    "transformer",
    "nntopo",
    "embed",
    "basic",
    "gpt",
    "bert",
    "pretrain",
]

if v"1.0.0" <= VERSION  < v"1.4.0"
    @test isempty(detect_ambiguities(Transformers, Transformers.GenerativePreTrain, Transformers.Basic, Transformers.Datasets))
end

@testset "Transformers" begin
  if Flux.use_cuda[]
      @info "Test CUDA"
    include("test_cuda.jl")
  end

  for t in tests
    fp = joinpath(dirname(@__FILE__), "test_$t.jl")
    @info "Test $(uppercase(t))"
    include(fp)
  end
end


text:
type Hmm
    e_to_i::Dict{String, Integer}
    i_to_e::Dict{Integer, String}
    s_to_i::Dict{String, Integer}
    i_to_s::Dict{Integer, String}
    pi_::Vector
    A::Matrix
    B::Matrix
    data::Matrix
end

function normalize(x)
    return x./sum(x)
end

# forward-backward algorithm for scoring and conditional probability

# Smoothing.
# Input:  The HMM (state and observation maps, and probabilities)
#         A list of T observations: E(0), E(1), ..., E(T-1)
#
# Ouptut: The posterior probability distribution over each state given all
#         of the observations: P(X(k)|E(0), ..., E(T-1) for 0 <= k <= T-1.
#
#         These distributions should be returned as a list of lists.

function scoring(model::Hmm)
    n_sample = size(model.data,1)
    n_state = size(model.A, 1)
    alpha = zeros(n_state, n_sample)
    alpha[:, 1] = vec(model.pi_) .* model.B[:,model.e_to_i[model.data[1,2]]]
    for i = 2:size(data,1)
        alpha[:, i] = vec(alpha[:, i-1]' * model.A) .* model.B[:,model.e_to_i[model.data[i, 2]]]
    end
    score_alpha = sum(alpha[:, end])
    beta = zeros(n_state, n_sample)
    beta[:, n_sample] =  ones(3)
    for i = (n_sample-1):-1:1
        beta[:, i] = model.A * beta[:, i+1] .* model.B[:, model.e_to_i[model.data[i+1, 2]]]
    end
    score_beta = sum(beta[:, 1] .* model.pi_ .* model.B[:, model.e_to_i[model.data[1,2]]])

    phi = alpha .* beta
    for i = 1:size(phi,2)
        phi[:,i] = normalize(phi[:, i])
    end
    return phi
end

# vitorbi

function matching(mode::Hmm)
    n_sample = size(model.data,1)
    n_state = size(model.A, 1)
    phi = zeros(n_state, n_sample)
    phi[:, 1] = normalize(model.pi_ .* model.B[:, model.e_to_i[model.data[1,2]]])
    eta = zeros(n_state, n_sample)
    eta[:, 1] = zeros(3)
    for i = 2:n_sample
        for j = 1:n_state
            temp = phi[:, i-1] .* model.A[:, j]
            eta[j,i] = indmax(temp)
            phi[j,i] = maximum(temp)
        end
        phi[:, i] = phi[:, i] .* model.B[:, model.e_to_i[model.data[i,2]]]
        phi[:, i] = normalize(phi[:, i])
    end
    @show eta[:, 1:10]
    state_optim = zeros(n_sample)
    state_optim[n_sample] = indmax(phi[:, n_sample])
    for i = (n_sample-1):-1:1
        state_optim[i] = eta[convert(Int,state_optim[i+1]), i+1]
    end
    state_optim = map(x -> model.i_to_s[convert(Int,x)], state_optim)
    return state_optim
end

function learning(model::Hmm)
    #learning its parameters
    #pi_, A, B
    #given the obs
end

function load_data(filename)
    num_ = readline(filename)
    num = parse(Int64, num_)
    data = []
    count = 1
    open(filename, "r") do f
        for line in readlines(f)
            if count == 1
                count += 1
                continue
            end
            str = split(chomp(line), ",")
            push!(data,str')
        end
        data = reduce(vcat,data)
    end
    return data
end

function test_Hmm()
    X_train, X_test, y_train, y_test = make_cla()
    model = Hmm()
    train!(model,X_train, y_train)
    predictions = predict(model,X_test)
    print("classification accuracy", accuracy(y_test, predictions))
end

####################################
## there is some problem to be fixed
####################################

function test_HMM()
    # state map
    weatherStateMap = Dict([("sunny", 1), ("rainy", 2), ("foggy", 3)])
    weatherStateIndex = Dict([(1, "sunny"), (2, "rainy"), (3, "foggy")])

    # observation map
    weatherObsMap = Dict([("no", 1), ("yes", 2)])
    weatherObsIndex = Dict([(1, "no"), (2, "yes")])

    # prior prob
    weatherprob = [0.5,0.25,0.25]

    # trasition probabilities

    weather_trans = [0.8 0.05 0.15;
                     0.2 0.6 0.2;
                     0.2 0.3 0.5]

    # obs

    weather_obs = [0.9 0.1 ;0.2 0.8 ;0.7 0.3]

    data = load_data("data/weather-test1-1000.txt")
    model = Hmm(weatherObsMap, weatherObsIndex, weatherStateMap,
                weatherStateIndex, weatherprob,
                weather_trans, weather_obs, data)

    phi = scoring(model)

    matching_state = matching(model)
end


text:
using Documenter, Triangulate, Literate, Pluto

function rendernotebook(name)
    input=joinpath(@__DIR__,"..","examples",name*".jl")
    output=joinpath(@__DIR__,"src",name*".html")
    session = Pluto.ServerSession();
    notebook = Pluto.SessionActions.open(session, input; run_async=false)
    html_contents = Pluto.generate_html(notebook)
    write(output, html_contents)
end

function make_all()
    rendernotebook("pluto-examples")

    makedocs(
        sitename="Triangulate.jl",
        modules = [Triangulate],
        clean = false,
        doctest = false,
        authors = "Juergen Fuhrmann, Francesco Furiani, Konrad Simon",
        repo="https://github.com/JuliaGeometry/Triangulate.jl",
        pages=[
            "Home"=>"index.md",
            "triangle-h.md",
            "changes.md",
            "allindex.md",
            "examples.md"
        ]
    )

    if !isinteractive()
        deploydocs(repo = "github.com/JuliaGeometry/Triangulate.jl.git")
    end
end

make_all()


text:
using ZStd
using Base.Test

@testset "helpers" begin
    @test typeof(ZStd.MAX_COMPRESSION) == Int
    @test ZStd.MAX_COMPRESSION > 0

    @test ZStd.check_zstd_error(UInt64(0)) == UInt64(0)
    @test_throws ZStd.ZStdError ZStd.check_zstd_error(typemax(UInt64))

    let
        err = try
            ZStd.check_zstd_error(typemax(UInt64))
        catch ex
            ex
        end
        err_txt = sprint(showerror, err)
        @test startswith(err_txt, "ZStd: Error (generic)")
    end

    @test typeof(ZStd.maxcompressedsize(1)) == Int
    @test ZStd.maxcompressedsize(1) > 0
    @test typeof(ZStd.maxcompressedsize(UInt(1))) == UInt
    @test ZStd.maxcompressedsize(UInt(1)) > UInt(0)

    @test ZStd.ZSTD_VERSION == v"1.3.3"
end


text:
using ReinforcementLearning

function random_mdp()
	g = Dict()
	s = State(1)

	mdp = MDP(20, 20)
	printmdp(mdp)

	println("\nPolicy Iteration")
	policy, V = policy_iteration(mdp; Ɣ=0.8, verbose=false)

	ss = getAllStates(mdp)
	for s in ss
		println("State: $(s), Value: $(V[s]), Action: $(policy.mapping[s][1][1])")
	end

	println("\nValue Iteration")
	policy, V = synchronous_value_iteration(mdp; Ɣ=0.8, verbose=true)
	for s in ss
		println("State: $(s), Value: $(V[s]), Action: $(policy.mapping[s][1][1])")
	end

	println("\nGauss-Seidel Value Iteration")
	policy, V = gauss_seidel_value_iteration(mdp; Ɣ=0.8, verbose=true)
	for s in ss
		println("State: $(s), Value: $(V[s]), Action: $(policy.mapping[s][1][1])")
	end

	V = iterative_policy_evaluation(mdp, policy, Ɣ=0.8; verbose=false)

	println("\nPolicy Evaluation")
	for s in ss
		println("State: $(s), Value: $(V[s])")
	end

end

function twostates()
	ss = [State(1), State(2)]
	as = [Action(1), Action(2), Action(3)]

	graph = Dict()
	graph[ss[1]] = Dict()
	graph[ss[2]] = Dict()
	graph[ss[1]][as[1]] = [(ss[1], 5.0, 0.5), (ss[2], 5.0, 0.5)]
	graph[ss[1]][as[2]] = [(ss[2], 10.0, 1.0)]
	graph[ss[1]][as[3]] = [(ss[1], 0.0, 1.0)]
	graph[ss[2]][as[1]] = [(ss[2], -1.0, 1.0)]

	mdp = MDP(2, 3, graph)
	printmdp(mdp)

	mapping = Dict()
	mapping[ss[1]] = [(as[1], 1.0)]
	mapping[ss[2]] = [(as[1], 1.0)]
	policy = Policy(mapping)
	V = iterative_policy_evaluation(mdp, policy, Ɣ=0.8; verbose=false)

	println("Policy Evaluation")
	for s in ss
		println("State: $(s), Value: $(V[s])")
	end

	println("Policy Iteration")
	policy, V = policy_iteration(mdp; Ɣ=0.8, verbose=false)
	for s in ss
		println("State: $(s), Value: $(V[s]), Action: $(policy.mapping[s][1][1])")
	end

	println("Value Iteration")
	policy, V = synchronous_value_iteration(mdp; Ɣ=0.8, verbose=false)
	for s in ss
		println("State: $(s), Value: $(V[s]), Action: $(policy.mapping[s][1][1])")
	end

end

function main()
	println("*** RANDOM MDP ***")
	random_mdp()

	#=
	println("*** Example MDP ***")
	twostates()
	=#
end

main()

max_stars_repo_name: ozanarkancan/ReinforcementLearning

max_issues_repo_name: ozanarkancan/ReinforcementLearning

max_forks_repo_name: ozanarkancan/ReinforcementLearning


text:
# ------------------------------------------------------------
# Truncated cosine and its Fourier transform

# [-1/2,1/2]
function tcos(x)
	if 0.0 <= x <= 0.5
		return cos(2*pi*x)
	else
		return 0.0
	end
end

function ftcos(xi)
	if abs(xi) == 1.0
		return 0.25 + 0.0*im
	else
		return im*xi*( 1 + exp(-pi*xi*im) ) / (2*pi*(1 - xi^2))
	end
end

# Overload tcos and ftcos to vector input
@vectorize_1arg Float64 tcos
@vectorize_1arg Float64 ftcos

# ------------------------------------------------------------
# Reconstruct truncated cosine in Haar basis

using GeneralizedSampling

J = 5
M = 2^(J+2)
# Both GeneralizedSampling and Winston (below) have a grid function
xi = GeneralizedSampling.grid(M, 0.5)
f = ftcos(xi)

T = Freq2Wave(xi, "haar", J)
wcoef = T \ f

# ------------------------------------------------------------
# Plot reconstruction

using IntervalWavelets
using Winston

x, yw = weval( real(wcoef), "haar", 10 )
# Waiting for IntervalWavelets to return x on [-0.5, 0.5] instead of [-1,1]
x -= 0.5

plot(x, yw, x, tcos(x))


text:
# This file was generated, do not modify it. # hide
using Gnuplot,Random
Random.seed!(124)
n = 30
x, y, z = randn(n), randn(n), randn(n)
@gp x y 5z "with points notit lt 1 pt 4 ps var "
save(term="pngcairo size 600,400", output="./code/burbuja1.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
export DiscordChannel

"""
A Discord channel.
More details [here](https://discordapp.com/developers/docs/resources/channel#channel-object).

Note: The name `Channel` is already used, hence the prefix.
"""
struct DiscordChannel <: DiscordObject
    id::Snowflake
    type::Optional{Int}
    guild_id::Optional{Snowflake}
    position::Optional{Int}
    permission_overwrites::Optional{Vector{Overwrite}}
    name::Optional{String}
    topic::OptionalNullable{String}
    nsfw::Optional{Bool}
    last_message_id::OptionalNullable{Snowflake}
    bitrate::Optional{Int}
    user_limit::Optional{Int}
    rate_limit_per_user::Optional{Int}
    recipients::Optional{Vector{User}}
    icon::OptionalNullable{String}
    owner_id::Optional{Snowflake}
    application_id::Optional{Snowflake}
    parent_id::OptionalNullable{Snowflake}
    last_pin_timestamp::OptionalNullable{DateTime}  # Not supposed to be nullable.
end
@boilerplate DiscordChannel :constructors :docs :lower :merge :mock


text:
using Manopt, Manifolds, Test, Random
using LinearAlgebra: I, eigvecs, tr
Random.seed!(42)

@testset "Riemannian quasi-Newton Methods" begin
    @testset "Mean of 3 Matrices" begin
        # Mean of 3 matrices
        A = [18.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0]
        B = [0.0 0.0 0.0 0.009; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0]
        C = [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; -5.0 0.0 0.0 0.0]
        ABC = [A, B, C]
        x_solution = mean(ABC)
        F(::Euclidean, x) = 0.5 * norm(A - x)^2 + 0.5 * norm(B - x)^2 + 0.5 * norm(C - x)^2
        gradF(::Euclidean, x) = -A - B - C + 3 * x
        M = Euclidean(4, 4)
        x = zeros(Float64, 4, 4)
        x_lrbfgs = quasi_Newton(
            M, F, gradF, x; stopping_criterion=StopWhenGradientNormLess(10^(-6))
        )
        @test norm(x_lrbfgs - x_solution) ≈ 0 atol = 10.0^(-14)
        # with Options
        lrbfgs_o = quasi_Newton(
            M,
            F,
            gradF,
            x;
            stopping_criterion=StopWhenGradientNormLess(10^(-6)),
            return_options=true,
        )
        @test lrbfgs_o.x == x_lrbfgs
        # with Cached Basis
        x_lrbfgs_cached = quasi_Newton(
            M,
            F,
            gradF,
            x;
            stopping_criterion=StopWhenGradientNormLess(10^(-6)),
            basis=get_basis(M, x, DefaultOrthonormalBasis()),
        )
        @test x_lrbfgs_cached == x_lrbfgs

        x_lrbfgs_cached_2 = quasi_Newton(
            M,
            F,
            gradF,
            x;
            stopping_criterion=StopWhenGradientNormLess(10^(-6)),
            basis=get_basis(M, x, DefaultOrthonormalBasis()),
            memory_size=-1,
        )
        @test isapprox(M, x_lrbfgs_cached_2, x_lrbfgs; atol=1e-5)

        x_clrbfgs = quasi_Newton(
            M,
            F,
            gradF,
            x;
            cautious_update=true,
            stopping_criterion=StopWhenGradientNormLess(10^(-6)),
        )
        @test norm(x_clrbfgs - x_solution) ≈ 0 atol = 10.0^(-14)

        x_rbfgs_Huang = quasi_Newton(
            M,
            F,
            gradF,
            x;
            memory_size=-1,
            stepsize=WolfePowellBinaryLinesearch(
                ExponentialRetraction(), ParallelTransport()
            ),
            stopping_criterion=StopWhenGradientNormLess(10^(-6)),
        )
        @test norm(x_rbfgs_Huang - x_solution) ≈ 0 atol = 10.0^(-14)

        for T in [InverseBFGS(), BFGS(), InverseDFP(), DFP(), InverseSR1(), SR1()]
            for c in [true, false]
                x_direction = quasi_Newton(
                    M,
                    F,
                    gradF,
                    x;
                    direction_update=T,
                    cautious_update=c,
                    memory_size=-1,
                    stopping_criterion=StopWhenGradientNormLess(10^(-12)),
                )
                @test norm(x_direction - x_solution) ≈ 0 atol = 10.0^(-14)
            end
        end
    end
    @testset "Rayleigh Quotient Minimzation" begin
        n = 4
        rayleigh_atol = 1e-8
        A = [2.0 1.0 0.0 3.0; 1.0 3.0 4.0 5.0; 0.0 4.0 3.0 2.0; 3.0 5.0 2.0 6.0]
        A = (A + A') / 2
        M = Sphere(n - 1)
        F(::Sphere, X) = X' * A * X
        gradF(::Sphere, X) = 2 * (A * X - X * (X' * A * X))
        x_solution = abs.(eigvecs(A)[:, 1])

        x = Matrix{Float64}(I, n, n)[n, :]
        x_lrbfgs = quasi_Newton(
            M,
            F,
            gradF,
            x;
            basis=get_basis(M, x, DefaultOrthonormalBasis()),
            memory_size=-1,
            stopping_criterion=StopWhenGradientNormLess(1e-9),
        )
        @test norm(abs.(x_lrbfgs) - x_solution) ≈ 0 atol = rayleigh_atol

        x_clrbfgs = quasi_Newton(
            M,
            F,
            gradF,
            x;
            cautious_update=true,
            stopping_criterion=StopWhenGradientNormLess(1e-9),
        )

        x_cached_lrbfgs = quasi_Newton(
            M,
            F,
            gradF,
            x;
            basis=get_basis(M, x, DefaultOrthonormalBasis()),
            memory_size=-1,
            stopping_criterion=StopWhenGradientNormLess(1e-9),
        )
        @test norm(abs.(x_cached_lrbfgs) - x_solution) ≈ 0 atol = rayleigh_atol

        for T in [
                InverseDFP(),
                DFP(),
                Broyden(0.5),
                InverseBroyden(0.5),
                Broyden(0.5, :Davidon),
                Broyden(0.5, :InverseDavidon),
                InverseBFGS(),
                BFGS(),
            ],
            c in [true, false]

            x_direction = quasi_Newton(
                M,
                F,
                gradF,
                x;
                direction_update=T,
                cautious_update=c,
                memory_size=-1,
                stopping_criterion=StopWhenGradientNormLess(5 * 1e-8),
            )
            @test norm(abs.(x_direction) - x_solution) ≈ 0 atol = rayleigh_atol
        end
    end
    @testset "Brockett" begin
        struct GradF
            A::Matrix{Float64}
            N::Diagonal{Float64,Vector{Float64}}
        end
        function (gradF::GradF)(::Stiefel, X::Array{Float64,2})
            AX = gradF.A * X
            XpAX = X' * AX
            return 2 .* AX * gradF.N .- X * XpAX * gradF.N .- X * gradF.N * XpAX
        end

        n = 4
        k = 2
        M = Stiefel(n, k)
        A = [2.0 1.0 0.0 3.0; 1.0 3.0 4.0 5.0; 0.0 4.0 3.0 2.0; 3.0 5.0 2.0 6.0]
        F(::Stiefel, X) = tr((X' * A * X) * Diagonal(k:-1:1))
        gradF = GradF(A, Diagonal(Float64.(collect(k:-1:1))))

        x = Matrix{Float64}(I, n, n)[:, 2:(k + 1)]
        x_inverseBFGSCautious = quasi_Newton(
            M,
            F,
            gradF,
            x;
            memory_size=8,
            vector_transport_method=ProjectionTransport(),
            retraction_method=QRRetraction(),
            cautious_update=true,
            stopping_criterion=StopWhenGradientNormLess(1e-6),
        )

        x_inverseBFGSHuang = quasi_Newton(
            M,
            F,
            gradF,
            x;
            memory_size=8,
            stepsize=WolfePowellBinaryLinesearch(QRRetraction(), ProjectionTransport()),
            vector_transport_method=ProjectionTransport(),
            retraction_method=QRRetraction(),
            cautious_update=true,
            stopping_criterion=StopWhenGradientNormLess(1e-6),
        )
        @test isapprox(M, x_inverseBFGSCautious, x_inverseBFGSHuang; atol=2e-5)
    end
end


text:
function create_noloss_dict(mn_data; existing_override=Dict())
   noloss_override = deepcopy(existing_override)
   try
      for (conv_idx, converter) in mn_data["converter"]
         if !("converter" in keys(noloss_override))
            noloss_override["converter"] = Dict()
         end
         for (bus, c1) in converter["c1"]
            if !(conv_idx in keys(noloss_override["converter"]))
               noloss_override["converter"][conv_idx] = Dict()
            end
            if !("c1" in keys(noloss_override["converter"][conv_idx]))
               noloss_override["converter"][conv_idx]["c1"] = Dict()
            end
            noloss_override["converter"][conv_idx]["c1"][bus] = 0.0
         end
         for (bus, c2) in converter["c2"]
            if !(conv_idx in keys(noloss_override["converter"]))
               noloss_override["converter"][conv_idx] = Dict()
            end
            if !("c2" in keys(noloss_override["converter"][conv_idx]))
               noloss_override["converter"][conv_idx]["c2"] = Dict()
            end
            noloss_override["converter"][conv_idx]["c2"][bus] = 0.0
         end
         for (bus, c3) in converter["c3"]
            if !(conv_idx in keys(noloss_override["converter"]))
               noloss_override["converter"][conv_idx] = Dict()
            end
            if !("c3" in keys(noloss_override["converter"][conv_idx]))
               noloss_override["converter"][conv_idx]["c3"] = Dict()
            end
            noloss_override["converter"][conv_idx]["c3"][bus] = 0.0
         end
         for (bus, sw1) in converter["sw1"]
            if !(conv_idx in keys(noloss_override["converter"]))
               noloss_override["converter"][conv_idx] = Dict()
            end
            if !("sw1" in keys(noloss_override["converter"][conv_idx]))
               noloss_override["converter"][conv_idx]["sw1"] = Dict()
            end
            noloss_override["converter"][conv_idx]["sw1"][bus] = 0.0
         end
         for (bus, sw2) in converter["sw2"]
            if !(conv_idx in keys(noloss_override["converter"]))
               noloss_override["converter"][conv_idx] = Dict()
            end
            if !("sw2" in keys(noloss_override["converter"][conv_idx]))
               noloss_override["converter"][conv_idx]["sw2"] = Dict()
            end
            noloss_override["converter"][conv_idx]["sw2"][bus] = 0.0
         end
         for (bus, sw3) in converter["sw3"]
            if !(conv_idx in keys(noloss_override["converter"]))
               noloss_override["converter"][conv_idx] = Dict()
            end
            if !("sw3" in keys(noloss_override["converter"][conv_idx]))
               noloss_override["converter"][conv_idx]["sw3"] = Dict()
            end
            noloss_override["converter"][conv_idx]["sw3"][bus] = 0.0
         end
      end
   catch KeyError
      println("No converters in the network: no loss parameters to change.")
   end
   for (subnet_idx, subnet) in mn_data["sn"]
      for (bus_idx,bus) in subnet["bus"]
         if "converter_c1" in keys(bus)
            if !("sn" in keys(noloss_override))
               noloss_override["sn"] = Dict()
            end
            if !(subnet_idx in keys(noloss_override["sn"]))
               noloss_override["sn"][subnet_idx] = Dict()
            end
            if !("bus" in keys(noloss_override["sn"][subnet_idx]))
               noloss_override["sn"][subnet_idx]["bus"] = Dict()
            end
            if !(bus_idx in keys(noloss_override["sn"][subnet_idx]["bus"]))
               noloss_override["sn"][subnet_idx]["bus"][bus_idx] = Dict()
            end
            for (conv_idx, conv_c1) in bus["converter_c1"]
               if !("converter_c1" in keys(noloss_override["sn"][subnet_idx]["bus"][bus_idx]))
                  noloss_override["sn"][subnet_idx]["bus"][bus_idx]["converter_c1"] = Dict()
               end
               noloss_override["sn"][subnet_idx]["bus"][bus_idx]["converter_c1"][conv_idx] = 0.0
            end
         end
         if "converter_c2" in keys(bus)
            if !("sn" in keys(noloss_override))
               noloss_override["sn"] = Dict()
            end
            if !(subnet_idx in keys(noloss_override["sn"]))
               noloss_override["sn"][subnet_idx] = Dict()
            end
            if !("bus" in keys(noloss_override["sn"][subnet_idx]))
               noloss_override["sn"][subnet_idx]["bus"] = Dict()
            end
            if !(bus_idx in keys(noloss_override["sn"][subnet_idx]["bus"]))
               noloss_override["sn"][subnet_idx]["bus"][bus_idx] = Dict()
            end
            for (conv_idx, conv_c2) in bus["converter_c2"]
               if !("converter_c2" in keys(noloss_override["sn"][subnet_idx]["bus"][bus_idx]))
                  noloss_override["sn"][subnet_idx]["bus"][bus_idx]["converter_c2"] = Dict()
               end
               noloss_override["sn"][subnet_idx]["bus"][bus_idx]["converter_c2"][conv_idx] = 0.0
            end
         end
         if "converter_c3" in keys(bus)
            if !("sn" in keys(noloss_override))
               noloss_override["sn"] = Dict()
            end
            if !(subnet_idx in keys(noloss_override["sn"]))
               noloss_override["sn"][subnet_idx] = Dict()
            end
            if !("bus" in keys(noloss_override["sn"][subnet_idx]))
               noloss_override["sn"][subnet_idx]["bus"] = Dict()
            end
            if !(bus_idx in keys(noloss_override["sn"][subnet_idx]["bus"]))
               noloss_override["sn"][subnet_idx]["bus"][bus_idx] = Dict()
            end
            for (conv_idx, conv_c3) in bus["converter_c3"]
               if !("converter_c3" in keys(noloss_override["sn"][subnet_idx]["bus"][bus_idx]))
                  noloss_override["sn"][subnet_idx]["bus"][bus_idx]["converter_c3"] = Dict()
               end
               noloss_override["sn"][subnet_idx]["bus"][bus_idx]["converter_c3"][conv_idx] = 0.0
            end
         end
         if "converter_sw1" in keys(bus)
            if !("sn" in keys(noloss_override))
               noloss_override["sn"] = Dict()
            end
            if !(subnet_idx in keys(noloss_override["sn"]))
               noloss_override["sn"][subnet_idx] = Dict()
            end
            if !("bus" in keys(noloss_override["sn"][subnet_idx]))
               noloss_override["sn"][subnet_idx]["bus"] = Dict()
            end
            if !(bus_idx in keys(noloss_override["sn"][subnet_idx]["bus"]))
               noloss_override["sn"][subnet_idx]["bus"][bus_idx] = Dict()
            end
            for (conv_idx, conv_sw1) in bus["converter_sw1"]
               if !("converter_sw1" in keys(noloss_override["sn"][subnet_idx]["bus"][bus_idx]))
                  noloss_override["sn"][subnet_idx]["bus"][bus_idx]["converter_sw1"] = Dict()
               end
               noloss_override["sn"][subnet_idx]["bus"][bus_idx]["converter_sw1"][conv_idx] = 0.0
            end
         end
         if "converter_sw2" in keys(bus)
            if !("sn" in keys(noloss_override))
               noloss_override["sn"] = Dict()
            end
            if !(subnet_idx in keys(noloss_override["sn"]))
               noloss_override["sn"][subnet_idx] = Dict()
            end
            if !("bus" in keys(noloss_override["sn"][subnet_idx]))
               noloss_override["sn"][subnet_idx]["bus"] = Dict()
            end
            if !(bus_idx in keys(noloss_override["sn"][subnet_idx]["bus"]))
               noloss_override["sn"][subnet_idx]["bus"][bus_idx] = Dict()
            end
            for (conv_idx, conv_sw2) in bus["converter_sw2"]
               if !("converter_sw2" in keys(noloss_override["sn"][subnet_idx]["bus"][bus_idx]))
                  noloss_override["sn"][subnet_idx]["bus"][bus_idx]["converter_sw2"] = Dict()
               end
               noloss_override["sn"][subnet_idx]["bus"][bus_idx]["converter_sw2"][conv_idx] = 0.0
            end
         end
         if "converter_sw3" in keys(bus)
            if !("sn" in keys(noloss_override))
               noloss_override["sn"] = Dict()
            end
            if !(subnet_idx in keys(noloss_override["sn"]))
               noloss_override["sn"][subnet_idx] = Dict()
            end
            if !("bus" in keys(noloss_override["sn"][subnet_idx]))
               noloss_override["sn"][subnet_idx]["bus"] = Dict()
            end
            if !(bus_idx in keys(noloss_override["sn"][subnet_idx]["bus"]))
               noloss_override["sn"][subnet_idx]["bus"][bus_idx] = Dict()
            end
            for (conv_idx, conv_sw3) in bus["converter_sw3"]
               if !("converter_sw3" in keys(noloss_override["sn"][subnet_idx]["bus"][bus_idx]))
                  noloss_override["sn"][subnet_idx]["bus"][bus_idx]["converter_sw3"] = Dict()
               end
               noloss_override["sn"][subnet_idx]["bus"][bus_idx]["converter_sw3"][conv_idx] = 0.0
            end
         end
      end
   end
   return noloss_override
end


text:
module Scalar

export scalar

"""
    scalar(x)

Convert `x` to a ScalarType.
"""
scalar

mutable struct NumberScalar{T<:Number} <: Number
    v::T
end
scalar(v::Number) = NumberScalar(v)
@inline Base.getindex(s::NumberScalar) = s.v
@inline Base.setindex!(s::NumberScalar{T}, v::T) where {T<:Number} = s.v = v

mutable struct RealScalar{T<:Real} <: Real
    v::T
end
scalar(v::Real) = RealScalar(v)
@inline Base.getindex(s::RealScalar) = s.v
@inline Base.setindex!(s::RealScalar{T}, v::T) where {T<:Real} = s.v = v

"""
    ScalarType{T}

ScalarType is a type of number used to represent a scalar value in a model.
Its value is accessible via `getindex` and can be updated via `setindex!`.
"""
const ScalarType{T} = Union{NumberScalar{T}, RealScalar{T}}
const ST = ScalarType

Base.convert(::Type{T}, s::ST{T}) where {T<:Number} = s[]
Base.convert(::Type{T}, s::ST) where {T<:Number} = convert(T, s[])
Base.convert(::Type{T}, s::ST) where {P, T<:ST{P}} = T(convert(P, s[]))

Base.promote_rule(::Type{<:ST{S}}, ::Type{T}) where {S<:Number, T<:Number} =
    promote_type(S, T)
Base.promote_rule(::Type{<:ST{S}}, ::Type{<:ST{T}}) where {S<:Number, T<:Number} =
    promote_type(S, T)

# getindex
@inline function Base.getindex(x::ST, i::Integer)
    @boundscheck isone(i) || throw(BoundsError(x, i))
    return x[]
end
@inline function Base.getindex(x::ST, I::Integer...)
    @boundscheck all(isone, I) || throw(BoundsError(x, I))
    return x[]
end
@inline Base.getindex(x::ST, ::CartesianIndex{0}) = x[] # to fix ambiguity
Base.@propagate_inbounds Base.getindex(x::ST, I::CartesianIndex) = x[I.I...]

# setindex
@inline Base.setindex!(x::ST{T}, v) where {T} = setindex!(x, convert(T, v)) # type conversion
@inline function Base.setindex!(x::ST, v, i::Integer)
    @boundscheck isone(i) || throw(BoundsError(x, i))
    return setindex!(x, v)
end
@inline function Base.setindex!(x::ST, v, I::Integer...)
    @boundscheck all(isone, I) || throw(BoundsError(x, I))
    return setindex!(x, v)
end
@inline Base.setindex!(x::ST, v, ::CartesianIndex{0}) = setindex!(x, v) # to fix ambiguity
Base.@propagate_inbounds Base.setindex!(x::ST, v, I::CartesianIndex) = setindex!(x, v, I.I...)

# mathematical operations
for op in (:+, :-, :conj, :real, :imag, :float, :exp, :log)
    @eval @inline Base.$op(x::ST) = $op(x[])
end

for op in (:+, :-, :*, :/, :\, :^, :(==))
    @eval @inline Base.$op(x::ST{T}, y::ST{T}) where {T} = $op(x[], y[])
end

# show
Base.show(io::IO, x::ST) = print(io, "scalar(", x[], ")")
Base.show(io::IO, ::MIME"text/plain", x::ST) = show(io, x[])

end # module Scalar


text:
# coding: utf-8

# In[1]:

#global list is the key
#otherwise we will lose the list throughout iterations
global factors=[]

# In[2]:

function factorization(num)

    #negative and float should be excluded
    if (num<0) || !(isinteger(num))

        error("negative or float is not allowed")

    end

    #if n is smaller than 4
    #prime number it is
    if num>4

        #exclude 1 and itself
        #the largest factor of n can not exceed the half of n
        #because 2 is the smallest factor
        #the range of factors we are gonna try starts from 2 to fld(num,2)+1
        #int function to solve the odd number problem
        for i in 2:(fld(num,2)+1)

            #the factorization process
            if num%i==0

                push!(factors,i)

                #return is crucial
                #if the number is not a prime number
                #it will stop function from appending non-prime factors
                #the next few lines will be ignored
                return factorization(Int32(num/i))

            end

        end

    end

    #append the last factor
    #it could be n itself if n is a prime number
    #in that case there is only one element in the list
    #or it could be the last indivisible factor of n which is also a prime number
    push!(factors,num)

    if length(factors)==1

        println(num," is a prime number")
        empty!(factors)

    end

end

# In[3]:

factorization(71392643)

# In[4]:

print(factors)


text:
using AGH
using Octavian, Test

@testset "Octavian validity" begin
    a = rand(2, 3)
    b = rand(3, 2)
    c = matmul(a, b)
    @test c == a * b
end


text:
using ForeignCallbacks
using Test
using InteractiveUtils

struct Message
    id::Int
    data::Ptr{Cvoid}
end

struct Message2
    id::Int
    data::Ref{Int}
end

@testset "Constructor" begin
    @test_throws AssertionError ForeignCallbacks.SingleConsumerStack{Ref{Int}}()
    @test_throws AssertionError ForeignCallbacks.SingleConsumerStack{Base.RefValue{Int}}()
    @test_throws AssertionError ForeignCallbacks.SingleConsumerStack{Message2}()
    @test_throws AssertionError ForeignCallbacks.SingleConsumerStack{Array}()
    @test_throws AssertionError ForeignCallbacks.SingleConsumerStack{Array{Int, 1}}()
end

@testset "SingleConsumerStack" begin
    lfq = ForeignCallbacks.SingleConsumerStack{Int}()

    @test ForeignCallbacks.popall!(lfq) == []
    push!(lfq, 1)
    @test ForeignCallbacks.popall!(lfq) == [1]
    @test ForeignCallbacks.popall!(lfq) == []

    GC.@preserve lfq begin
        ptr = Base.pointer_from_objref(lfq)
        ForeignCallbacks.unsafe_push!(ptr, 2)
    end
    # TODO: Test no load from TLS in `unsafe_push!`
    @test ForeignCallbacks.popall!(lfq) == [2]
    @test ForeignCallbacks.popall!(lfq) == []
end

@testset "callback" begin
    ch = Channel{Int}(Inf)
    callback = ForeignCallbacks.ForeignCallback{Int}() do val
        put!(ch, val)
        return
    end

    GC.@preserve callback begin
        token = ForeignCallbacks.ForeignToken(callback)
        ForeignCallbacks.notify!(token, 1)
        @test fetch(ch) === 1
    end
end

@testset "IR" begin
    let llvm = sprint(io->code_llvm(io, push!, Tuple{ForeignCallbacks.SingleConsumerStack{Int}, Int}))
        @test !contains(llvm, "%thread_ptr")
        @test !contains(llvm, "%pgcstack")
        @test !contains(llvm, "%gcframe")
    end
    let llvm = sprint(io->code_llvm(io, ForeignCallbacks.unsafe_push!, Tuple{Ptr{Cvoid}, Int}))
        @test !contains(llvm, "%thread_ptr")
        @test !contains(llvm, "%pgcstack")
        @test !contains(llvm, "%gcframe")
    end
    let llvm = sprint(io->code_llvm(io, ForeignCallbacks.notify!, Tuple{ForeignCallbacks.ForeignToken, Int}))
        @test !contains(llvm, "%thread_ptr")
        @test !contains(llvm, "%pgcstack")
        @test !contains(llvm, "%gcframe")
    end
    let llvm = sprint(io->code_llvm(io, ForeignCallbacks.unsafe_push!, Tuple{Ptr{Cvoid}, Message}))
        @test !contains(llvm, "%thread_ptr")
        @test !contains(llvm, "%pgcstack")
        @test !contains(llvm, "%gcframe")
    end
    let llvm = sprint(io->code_llvm(io, ForeignCallbacks.notify!, Tuple{ForeignCallbacks.ForeignToken, Message}))
        @test !contains(llvm, "%thread_ptr")
        @test !contains(llvm, "%pgcstack")
        @test !contains(llvm, "%gcframe")
    end
end

if Threads.nthreads() == 1 && Sys.CPU_THREADS > 1
    @info "relaunching with" threads = Sys.CPU_THREADS
    cmd = `$(Base.julia_cmd()) --threads=$(Sys.CPU_THREADS) $(@__FILE__)`
    @test success(pipeline(cmd, stdout=stdout, stderr=stderr))
    exit()
end

function producer!(lfq)
    for i in 1:100
        push!(lfq, i)
        yield()
    end
end

function unsafe_producer!(lfq)
    for i in 1:100
        GC.@preserve lfq begin
            ptr = Base.pointer_from_objref(lfq)
            ForeignCallbacks.unsafe_push!(ptr, i)
        end
        yield()
    end
end

function consumer!(lfq)
    acc = 0

    done = false
    while !done
        for x in ForeignCallbacks.popall!(lfq)
            acc += x
        end
        done = acc == sum(1:100)*2*Threads.nthreads()
        yield()
    end
end

@testset "Queue threads" begin
    @test Threads.nthreads() == Sys.CPU_THREADS
    let lfq = ForeignCallbacks.SingleConsumerStack{Int}()
        @sync begin
            for n in 1:2*Threads.nthreads()
                Threads.@spawn producer!(lfq)
            end
            Threads.@spawn consumer!(lfq)
        end
        @test true
    end

    let lfq = ForeignCallbacks.SingleConsumerStack{Int}()
        @sync begin
            for n in 1:2*Threads.nthreads()
                Threads.@spawn unsafe_producer!(lfq)
            end
            Threads.@spawn consumer!(lfq)
        end
        @test true
    end
end

@testset "Callback threads" begin
    ch = Channel{Int}(Inf)
    callback = ForeignCallbacks.ForeignCallback{Int}() do val
        put!(ch, val)
        return
    end

    GC.@preserve callback begin
        token = ForeignCallbacks.ForeignToken(callback)
        for n in 1:2*Threads.nthreads()
            Threads.@spawn ForeignCallbacks.notify!(token, 1)
        end
    end

    acc = 0
    while acc < 2*Threads.nthreads()
        acc += fetch(ch)
    end
    @test acc == 2*Threads.nthreads()
end


text:
using RecipesBase

@recipe function f(b::Curve)
    grid --> false
    axis --> false
    ticks --> false
    legend --> false
    aspect_ratio --> :equal
    linewidth --> 2.0
    @series begin
        xy(u) = b(u)
        lim = b.limits
        Δ = diam(lim)/100
        θ = lim.lo:Δ:lim.hi
        Tuple.(xy.(θ))
    end
end

@recipe function f(b::Curve, lim::Interval)
    grid --> false
    axis --> false
    ticks --> false
    legend --> false
    aspect_ratio --> :equal
    linewidth --> 2.0
    @series begin
        xy(u) = b(u)
        Δ = diam(lim)/100
        θ = lim.lo:Δ:lim.hi
        Tuple.(xy.(θ))
    end
end

@recipe function f(a::Tuple{SVector{D, T}, SVector{D, T}}) where {D, T}
    grid --> false
    axis --> false
    ticks --> false
    legend --> false
    aspect_ratio --> :equal
    linewidth --> 1.0
    linecolor --> :red
    markercolor --> :red
    markershape --> :circle
    markersize --> 5.0
    markerstrokecolor --> :red
    @series begin
        Tuple.([a...])
    end
end


text:
module DeformableRegistration

# helping functions
include("helpers/matrixFreeHelpers.jl")
include("helpers/registrationOptions.jl");
#include("helpers/objectiveFunctionCreation.jl")

# submodules
include("ImageProcessing.jl")
include("Transformation.jl")
include("Interpolation.jl")

include("Distance.jl")
include("Regularizer.jl")
include("Optimization.jl")

include("Examples.jl")

try
    include("Visualization.jl")
catch err
    println("caught: ", err)
    println("This is 'normal' during automated testing: ")
end

end


text:
using SymbolicUtils
using SymbolicUtils.Rewriters
using Symbolics
# using SymbolicNumericIntegration

using PyCall
sympy = pyimport("sympy")

@syms x 𝐞 a b c d e p t m n z

axion_rules = [
    @rule ^(𝐞, ~k) => exp(~k)
    @rule 𝐞 => MathConstants.e
]

function convert_axiom(name::AbstractString)
    fd = open(name, "r")
    re = r"\[.*\]"
    L = []

    D = Dict{Any,Int}()

    for (lineno,line) in enumerate(readlines(fd))
        line = strip(line)
        if length(line) < 3 || line[1:2] == "--" continue end
        ma = match(re, line)

        if ma != nothing
            line = ma.match
            println(line)
            line = replace(line, "%e" => 𝐞)
            line  = replace(line, "%i" => im)

            try
                expr = Meta.parse(line)
                P = eval(expr)

                print(lineno, ": ", P[1], " => ")
                P[1] = Prewalk(PassThrough(Chain(axion_rules)))(P[1])
                P[4] = Prewalk(PassThrough(Chain(axion_rules)))(P[4])

                for ν in Symbolics.get_variables(P[1])
                    if !isequal(ν, x) && !haskey(D, ν)
                        D[ν] = length(D) + 2
                    end
                end

                P[1] = substitute(P[1], D)
                P[4] = substitute(P[4], D)

                println(P[1])
                push!(L, P)
            catch err
                printstyled(lineno, ": ", err, '\n'; color=:red)
                printstyled('\t', line, '\n'; color=:red)
            end
        end
    end

    close(fd)
    return L
end

function convert_axiom(sym)
    name =
        if sym == :Apostle || sym == 1
            "Apostol Problems.input"
        elseif sym == :Bondarenko || sym == 2
            "Bondarenko Problems.input"
        elseif sym == :Bronstein || sym == 3
            "Bronstein Problems.input"
        elseif sym == :Charlwood || sym == 4
            "Charlwood Problems.input"
        elseif sym == :Hearn || sym == 5
            "Hearn Problems.input"
        elseif sym == :Hebisch || sym == 6
            "Hebisch Problems.input"
        elseif sym == :Jeffrey || sym == 7
            "Jeffrey Problems.input"
        elseif sym == :Moses || sym == 8
            "Moses Problems.input"
        elseif sym == :Stewart || sym == 9
            "Stewart Problems.input"
        elseif sym == :Timofeev || sym == 10
            "Timofeev Problems.input"
        elseif sym == :Welz || sym == 11
            "Welz Problems.input"
        elseif sym == :Wester || sym == 12
            "Wester Problems.input"
        end

        name = joinpath("test/0", name)
        println(name)
        return convert_axiom(name)
end

function test_point(complex_plane, radius)
    if complex_plane
        return radius * sqrt(rand()) * cis(2π*rand())
    else
        return Complex(radius * (2*rand() - 1))
    end
end

function accept_integral(sol, ans, x; radius=1.0, abstol=1e-3, n=5)
    try
        ϵ = zeros(n)
        for i = 1:n
            x₀ = test_point(true, radius)
            ϵ[i] = abs(substitute(sol - ans, Dict(x => x₀)))
        end
        return maximum(ϵ) - minimum(ϵ) < abstol
    catch e
    end
    return false
end

function test_axiom(L, try_sympy=true; kwargs...)
    n_ok = 0
    n_fail = 0
    n_diff = 0
    n_catch = 0

    for (i,p) in enumerate(L)
        try
            printstyled(i, ": "; color=:yellow)
            eq = p[1]
            x = p[2]
            ans = p[4]

            printstyled(eq, '\n'; color=:green)
            sol = integrate(eq, x; kwargs...)[1]

            if isequal(sol, 0)
                printstyled("\tFailure\n"; color=:red)
                n_fail += 1
            elseif accept_integral(sol, ans, x)
                printstyled("\tSuccess:\t", sol, '\n'; color=:cyan)
                n_ok += 1
            else
                printstyled("\tDiscrepancy:\t", sol, '\n'; color=:yellow)
                n_diff += 1
            end

            if try_sympy
                s = pythonize(eq)
                py = sympy.integrate(s, sympy.Symbol(string(x)))
                printstyled("\tSymPy      :\t", string(py)[10:end], '\n'; color=:magenta)
            end

            printstyled("\tAnswer: \t", ans, '\n'; color=:blue)
        catch e
            printstyled(i, ": ", e, '\n'; color=:red)
            n_catch += 1
        end
    end

    println()
    printstyled("Success     = ", n_ok, '\n'; color=:green)
    printstyled("Failure     = ", n_fail, '\n'; color=:red)
    printstyled("Discrepancy = ", n_diff, '\n'; color=:yellow)
    printstyled("Exception   = ", n_catch, '\n'; color=:cyan)
end

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

pythonize(eq::SymbolicUtils.Add) = "(" * join(pythonize.(arguments(eq)), ")+(") * ")"
pythonize(eq::SymbolicUtils.Mul) = "(" * join(pythonize.(arguments(eq)), ")*(") * ")"

function pythonize(eq::SymbolicUtils.Pow)
    a = pythonize.(arguments(eq))
    "(" * a[1] * ")**(" * a[2] * ")"
end

pythonize(eq::SymbolicUtils.Term) = string(operation(eq)) * '(' * pythonize(arguments(eq)[1]) * ')'

function pythonize(eq)
    s = string(eq)
    s = replace(s, "π" => "pi")
    s = replace(s, "im" => "j")
    s
end


text:
using Documenter, CmdStan

DOC_ROOT = rel_path_cmdstan("..", "docs")

makedocs( root = DOC_ROOT,
  modules = [CmdStan],
  sitename = "StanJulia/CmdStan.jl",
  authors = "Rob J Goedman",
  pages = Any[
      "Home" => "INTRO.md",
      "Installation" => "INSTALLATION.md",
      "Walkthrough" => "WALKTHROUGH.md",
      "Versions" => "VERSIONS.md",
      "Index" => "index.md"
  ]
)

deploydocs(
  root = DOC_ROOT,
  repo = "github.com/StanJulia/CmdStan.jl.git",
 )


text:
import Logging
import TerminalLoggers
Logging.global_logger(TerminalLoggers.TerminalLogger())

import UnPack
import JSON
import ArgParse
import TurbulenceConvection

import CloudMicrophysics
const CM = CloudMicrophysics
const CM0 = CloudMicrophysics.Microphysics_0M
const CM1 = CloudMicrophysics.Microphysics_1M

import ClimaCore
const CC = ClimaCore
import SciMLBase

import OrdinaryDiffEq
const ODE = OrdinaryDiffEq
import StaticArrays: SVector

const tc_dir = dirname(dirname(pathof(TurbulenceConvection)))

include(joinpath(tc_dir, "driver", "initial_conditions.jl"))
include(joinpath(tc_dir, "diagnostics", "compute_diagnostics.jl"))
include(joinpath(tc_dir, "driver", "parameter_set.jl"))
include(joinpath(tc_dir, "driver", "Cases.jl"))
include(joinpath(tc_dir, "driver", "dycore.jl"))
import .Cases

struct Simulation1d{IONT, G, S, GM, C, EDMF, D, TIMESTEPPING, STATS, PS}
    io_nt::IONT
    grid::G
    state::S
    gm::GM
    Case::C
    edmf::EDMF
    diagnostics::D
    TS::TIMESTEPPING
    Stats::STATS
    param_set::PS
    skip_io::Bool
    adapt_dt::Bool
    cfl_limit::Float64
    dt_min::Float64
end

function Simulation1d(namelist)
    TC = TurbulenceConvection
    param_set = create_parameter_set(namelist)

    FT = Float64
    skip_io = namelist["stats_io"]["skip"]
    adapt_dt = namelist["time_stepping"]["adapt_dt"]
    cfl_limit = namelist["time_stepping"]["cfl_limit"]
    dt_min = namelist["time_stepping"]["dt_min"]

    grid = TC.Grid(FT(namelist["grid"]["dz"]), namelist["grid"]["nz"])
    Stats = skip_io ? nothing : TC.NetCDFIO_Stats(namelist, grid)
    case = Cases.get_case(namelist)
    ref_params = Cases.reference_params(case, grid, param_set, namelist)

    gm = TC.GridMeanVariables(namelist, grid, param_set)
    Sur = TC.SurfaceBase(Cases.get_surface_type(case); namelist, ref_params)
    Fo = TC.ForcingBase{Cases.get_forcing_type(case)}()
    Rad = TC.RadiationBase{Cases.get_radiation_type(case)}()

    Case = Cases.CasesBase(case, namelist, grid, param_set, Sur, Fo, Rad)
    edmf = TC.EDMF_PrognosticTKE(namelist, grid, param_set)
    TS = TC.TimeStepping(namelist)

    n_updrafts = edmf.n_updrafts

    cspace = TC.center_space(grid)
    fspace = TC.face_space(grid)

    cent_prog_fields() = TC.FieldFromNamedTuple(cspace, cent_prognostic_vars(FT, n_updrafts))
    face_prog_fields() = TC.FieldFromNamedTuple(fspace, face_prognostic_vars(FT, n_updrafts))
    aux_cent_fields = TC.FieldFromNamedTuple(cspace, cent_aux_vars(FT, n_updrafts))
    aux_face_fields = TC.FieldFromNamedTuple(fspace, face_aux_vars(FT, n_updrafts))
    diagnostic_cent_fields = TC.FieldFromNamedTuple(cspace, cent_diagnostic_vars(FT, n_updrafts))
    diagnostic_face_fields = TC.FieldFromNamedTuple(fspace, face_diagnostic_vars(FT, n_updrafts))

    prog = CC.Fields.FieldVector(cent = cent_prog_fields(), face = face_prog_fields())
    aux = CC.Fields.FieldVector(cent = aux_cent_fields, face = aux_face_fields)
    diagnostics = CC.Fields.FieldVector(cent = diagnostic_cent_fields, face = diagnostic_face_fields)

    # `nothing` goes into State because OrdinaryDiffEq.jl owns tendencies.
    state = TC.State(prog, aux, nothing)

    TC.compute_ref_state!(state, grid, param_set; ref_params...)

    io_nt = (;
        ref_state = TC.io_dictionary_ref_state(),
        aux = TC.io_dictionary_aux(),
        diagnostics = io_dictionary_diagnostics(),
    )

    return Simulation1d(
        io_nt,
        grid,
        state,
        gm,
        Case,
        edmf,
        diagnostics,
        TS,
        Stats,
        param_set,
        skip_io,
        adapt_dt,
        cfl_limit,
        dt_min,
    )
end

function TurbulenceConvection.initialize(sim::Simulation1d, namelist)
    TC = TurbulenceConvection
    state = sim.state
    Cases.initialize_profiles(sim.Case, sim.grid, sim.gm, state)
    satadjust(sim.gm, sim.grid, sim.state)

    Cases.initialize_surface(sim.Case, sim.grid, state, sim.param_set)
    Cases.initialize_forcing(sim.Case, sim.grid, state, sim.gm, sim.param_set)
    Cases.initialize_radiation(sim.Case, sim.grid, state, sim.gm, sim.param_set)

    initialize_edmf(sim.edmf, sim.grid, state, sim.Case, sim.gm, sim.TS)

    sim.skip_io && return nothing
    TC.initialize_io(sim.io_nt.ref_state, sim.Stats)
    TC.io(sim.io_nt.ref_state, sim.Stats, state) # since the reference prog is static

    TC.initialize_io(sim.io_nt.aux, sim.Stats)
    TC.initialize_io(sim.io_nt.diagnostics, sim.Stats)

    # TODO: depricate
    TC.initialize_io(sim.gm, sim.Stats)
    TC.initialize_io(sim.Case, sim.Stats)
    TC.initialize_io(sim.edmf, sim.Stats)

    TC.open_files(sim.Stats)
    TC.write_simulation_time(sim.Stats, sim.TS.t)

    TC.io(sim.io_nt.aux, sim.Stats, state)
    TC.io(sim.io_nt.diagnostics, sim.Stats, sim.diagnostics)

    # TODO: depricate
    TC.io(sim.gm, sim.grid, state, sim.Stats)
    TC.io(sim.Case, sim.grid, state, sim.Stats)
    TC.io(sim.edmf, sim.grid, state, sim.Stats, sim.TS, sim.param_set)
    TC.close_files(sim.Stats)

    return
end

include("callbacks.jl")

function run(sim::Simulation1d; time_run = true)
    TC = TurbulenceConvection
    iter = 0
    grid = sim.grid
    state = sim.state
    prog = state.prog
    aux = state.aux
    TS = sim.TS
    diagnostics = sim.diagnostics
    sim.skip_io || TC.open_files(sim.Stats) # #removeVarsHack

    t_span = (0.0, sim.TS.t_max)
    params = (;
        edmf = sim.edmf,
        grid = grid,
        gm = sim.gm,
        aux = aux,
        io_nt = sim.io_nt,
        case = sim.Case,
        diagnostics = diagnostics,
        TS = sim.TS,
        Stats = sim.Stats,
        skip_io = sim.skip_io,
        adapt_dt = sim.adapt_dt,
        cfl_limit = sim.cfl_limit,
        dt_min = sim.dt_min,
    )

    callback_io = ODE.DiscreteCallback(condition_io, affect_io!; save_positions = (false, false))
    callback_cfl = ODE.DiscreteCallback(condition_every_iter, monitor_cfl!; save_positions = (false, false))
    callback_cfl = sim.edmf.Precip.precipitation_model == "clima_1m" ? (callback_cfl,) : ()
    callback_dtmax = ODE.DiscreteCallback(condition_every_iter, dt_max!; save_positions = (false, false))
    callback_filters = ODE.DiscreteCallback(condition_every_iter, affect_filter!; save_positions = (false, false))
    callback_adapt_dt = ODE.DiscreteCallback(condition_every_iter, adaptive_dt!; save_positions = (false, false))
    callback_adapt_dt = sim.adapt_dt ? (callback_adapt_dt,) : ()

    callbacks = ODE.CallbackSet(callback_adapt_dt..., callback_dtmax, callback_cfl..., callback_filters, callback_io)

    prob = ODE.ODEProblem(TC.∑tendencies!, state.prog, t_span, params; dt = sim.TS.dt)

    # TODO: LES_driven_SCM is currently unstable w.r.t. higher order moments (HOM).
    # So, we tell OrdinaryDiffEq.jl to not perform NaNs check on the solution
    # so that it doesn't abort early (as the HOM prognostic variables are 1-way coupled)
    unstable_check_kwarg(::Cases.LES_driven_SCM) = (; unstable_check = (dt, u, p, t) -> false)
    unstable_check_kwarg(case) = ()

    kwargs = (;
        progress_steps = 100,
        save_start = false,
        saveat = last(t_span),
        callback = callbacks,
        progress = true,
        unstable_check_kwarg(sim.Case.case)...,
        progress_message = (dt, u, p, t) -> t,
    )
    alg = ODE.Euler()

    if time_run
        sol = @timev ODE.solve(prob, alg; kwargs...)
    else
        sol = ODE.solve(prob, alg; kwargs...)
    end

    sim.skip_io || TC.close_files(sim.Stats) # #removeVarsHack
    if first(sol.t) == sim.TS.t_max
        return :success
    else
        return :simulation_aborted
    end
end

main(namelist; kwargs...) = @timev main1d(namelist; kwargs...)

nc_results_file(stats::TC.NetCDFIO_Stats) = stats.path_plus_file
nc_results_file(::Nothing) = @info "The simulation was run without IO, so no nc files were exported"

function main1d(namelist; time_run = true)
    # TODO: generalize convesion of arrays from namelist to `SVector`s.
    _p = namelist["turbulence"]["EDMF_PrognosticTKE"]["general_ent_params"]
    namelist["turbulence"]["EDMF_PrognosticTKE"]["general_ent_params"] = SVector{length(_p)}(_p)
    sim = Simulation1d(namelist)
    TurbulenceConvection.initialize(sim, namelist)
    if time_run
        return_code = @timev run(sim; time_run)
    else
        return_code = run(sim)
    end
    return_code == :success && println("The simulation has completed.")
    return nc_results_file(sim.Stats), return_code
end

function parse_commandline()
    s = ArgParse.ArgParseSettings(; description = "Run case input")

    ArgParse.@add_arg_table! s begin
        "case_name"
        help = "The case name"
        arg_type = String
        required = true
    end

    return ArgParse.parse_args(s)
end

if abspath(PROGRAM_FILE) == @__FILE__

    args = parse_commandline()
    case_name = args["case_name"]

    namelist = open("namelist_" * "$case_name.in", "r") do io
        JSON.parse(io; dicttype = Dict, inttype = Int64)
    end
    main(namelist)
end


text:
using DataInterpolations, FileIO, CSV, DataFrames, Dates, Statistics
using StatsPlots, JLD2, Flux, ForwardDiff, GalacticOptim, LinearAlgebra, LogExpFunctions
using Parameters, Distributions, OrdinaryDiffEq, DiffEqCallbacks, Sundials, Tullio
using RecursiveArrayTools, Interpolations, MCMCChains, Suppressor, Interpolations
using RCall
import KenyaCoVaccines

##

@load("data/p_ID.jld2")
@load("data/N_kenya.jld2")
model = @suppress_err load("modelfits/Nairobi_model.jld2")["model"]
age_symptoms = model.p_symp
veff1 = 1 - model.VE_severe_disease_risk[2]
veff2 = 1 - model.VE_severe_disease_risk[3]
veff_death1 = 0.2
veff_death2 = 0.1

include("generation_methods.jl");
include("inference_methods.jl");

@load("data/fitted_mort_by_county_age_variant.jld2")
@load("data/fitted_rel_risk_ICU.jld2")
@load("data/fitted_rel_risk_hosp.jld2")

## Get vaccine scenarios

t_0 = (Date(2021, 9, 27) - Date(2020, 2, 20)).value #Vaccines start Sept 2021
t_end = (Date(2021, 9, 27) + Month(18) - Date(2020, 2, 20)).value #Option where vaccines finish in 18 months
t_end_rapid = (Date(2021, 9, 27) + Month(6) - Date(2020, 2, 20)).value #Option where vaccines finish in 6 months
ts = [0; t_0; t_end; t_end + 100] .+ 14 # Time points for 18 month deployment lagged 14 days to account for delay between jab and maximum effectiveness
ts_rapid = [0; t_0; t_end_rapid; t_end_rapid + 100] .+ 14 # Time points for 6 month deployment

# Number of people DOUBLE jabbed by time points under three coverages
fully_vaccinated_30perc = [0; 0; 15.9e6 / 2; 15.9e6 / 2]
fully_vaccinated_50perc = [0; 0; 25e6 / 2; 25e6 / 2]
fully_vaccinated_70perc = [0; 0; 35e6 / 2; 35e6 / 2]

## Create vaccine callbacks

include("create_vac_callbacks.jl");

## Get basic models for each county

fitfiles = readdir("modelfits", join=true)
model_vect = [@suppress_err load(filename)["model"] for filename in fitfiles]

##Test model

model = model_vect[30]
perc_over_18 = sum(N_kenya[5:end, model.areaname]) / sum(N_kenya[5:end, :])
KenyaCoVaccines.change_prob_immune_escape_with_waning!(model, perc_over_18; startdate=Date(2020, 12, 1), enddate=Date(2021, 12, 1))
VE_acq_base = model.VE_acquisition
VE_inf_base = model.VE_infectiousness
VE_sev_base = model.VE_severe_disease_risk

##

θs = model |> flatten_chain
# function immune_escape!(integrator)
#     println(size(integrator.u))
#     R = @view integrator.u[:, 7, :]
#     #Reduce generation time by 30%
#     gen_time_scale = 0.7
#     integrator.p[1] = 1 * (1 / gen_time_scale) * integrator.p[1]# β₀
#     integrator.p[6] = (1 / gen_time_scale) * integrator.p[6] # α
#     integrator.p[8:10] = (1 / gen_time_scale) * integrator.p[8:10]# αP, γA, γM

#     #Reduce vaccine protection against acquiring and transmitting infection by 50%
#     integrator.p[19:24] = 0.5 * integrator.p[19:24] #ve_ac1, ve_ac2, ve_ac3, ve_inf1, ve_inf2, ve_inf3
#     #Reduce protection from reinfection by 50%
#     reinf_red = 0.5
#     integrator.p[28] = (1.0 * reinf_red + 0.16 * (1 - reinf_red)) / 0.16
#     #50% of R -> W at emergence of immune escape variant but set loss of complete immunity to 0
#     integrator.u[:, 7, :] .-= 0.5 .* R
#     integrator.u[:, 8, :] .+= 0.5 .* R
#     integrator.p[11] = 0.0 # ω

# end

fs = falses(6, 11, 5)
fs[:, 6, :] .= true
idxs_R_cvode = fs[:]
fs = falses(6, 11, 5)
fs[:, 7, :] .= true
idxs_W1_cvode = fs[:]
fs = falses(6, 11, 5)
fs[:, 8, :] .= true
idxs_W2_cvode = fs[:]

function immune_escape_CVODE!(integrator)

    #Reduce generation time by 30%
    gen_time_scale = 0.7
    integrator.p[1] = 1 * (1 / gen_time_scale) * integrator.p[1]# β₀
    integrator.p[6] = (1 / gen_time_scale) * integrator.p[6] # α
    integrator.p[8:10] = (1 / gen_time_scale) * integrator.p[8:10]# αP, γA, γM

    #Reduce vaccine protection against acquiring and transmitting infection by 50%
    integrator.p[19:24] = 0.5 * integrator.p[19:24] #ve_ac1, ve_ac2, ve_ac3, ve_inf1, ve_inf2, ve_inf3
    #Reduce protection from reinfection in W1 by 0%
    reinf_red = 0.0
    integrator.p[28] = (1.0 * reinf_red + 0.16 * (1 - reinf_red)) / 0.16
    #50% of R -> W2 and W1 -> W2 at emergence of immune escape variant but set loss of complete immunity to 0
    integrator.u[idxs_W2_cvode] .+= 0.5 .* integrator.u[idxs_R_cvode]
    integrator.u[idxs_W2_cvode] .+= 0.5 .* integrator.u[idxs_W1_cvode]
    integrator.u[idxs_R_cvode] .-= 0.5 .* integrator.u[idxs_R_cvode]
    integrator.u[idxs_W1_cvode] .-= 0.5 .* integrator.u[idxs_W1_cvode]

    # integrator.p[11] = 0.0 # ω

end

cb_imm_esc_cvode = PresetTimeCallback([349.0], immune_escape_CVODE!, save_positions=(false, false))

# cb_imm_esc = CallbackSet()
solver = CVODE_BDF(linear_solver=:GMRES)
# solver = Rosenbrock23(autodiff = false)
##Generate random vaccine effectivenesses
VE_acqs = [[0.0, rand(Uniform(0.5, 0.6)), rand(Uniform(0.6, 0.85))] for i = 1:2000]
VE_infs = [[0.0, 0.0, rand(Uniform(0.0, 0.69))] for i = 1:2000]
VE_sevs = [[0.0, rand(Uniform(0.8, 0.9)), rand(Uniform(0.95, 0.99))] for i = 1:2000]
VE_deaths = [[0.0, rand(Uniform(0.9, 0.95)), rand(Uniform(0.95, 0.99))] for i = 1:2000]

@save("vac_effs.jld2", VE_acqs, VE_infs, VE_sevs, VE_deaths)
#Matrix versions
VE_waning = 0.0
VE_sevs_mat = 1 .- [vecvec_to_mat(VE_sevs) vecvec_to_mat(VE_sevs)[:, end] VE_waning .* vecvec_to_mat(VE_sevs)[:, end]]
VE_deaths_mat = 1 .- [vecvec_to_mat(VE_deaths) vecvec_to_mat(VE_deaths)[:, end] VE_waning .* vecvec_to_mat(VE_deaths)[:, end]]

##
rep = 1062
@time sol_vw = solve_model_waning_vac(θs[rep], model, Date(2023, 6, 1), cb_imm_esc_cvode, vac_cbs[7], solver, VE_acqs[rep], VE_infs[rep], VE_sevs[rep])
# ι_all = clamp!(KenyaCoVaccines.get_incidence_time_array_by_age_dose(sol_vw, 10), 0.0, Inf) #All infections
# ι_dis_weighted = clamp!(KenyaCoVaccines.get_incidence_time_array_by_age_dose(sol_vw, 11), 0.0, Inf) #Weighted by infection risk
E = [sum(u[:, 2, :]) for u in sol_vw.u]
A = [sum(u[:, 3, :]) for u in sol_vw.u]
R = [sum(u[:, 6, :]) for u in sol_vw.u]
W = [sum(u[:, 7, :]) for u in sol_vw.u]
W2 = [sum(u[:, 8, :]) for u in sol_vw.u]
F = [sum(u[:, 9, :]) for u in sol_vw.u]
C = [sum(u[:, 10, :]) for u in sol_vw.u]
C2 = [sum(u[:, 11, :]) for u in sol_vw.u]

V1 = [sum(u[:, 1:8, 2]) for u in sol_vw.u]
V2 = [sum(u[:, 1:8, 3]) for u in sol_vw.u]
WV1 = [sum(u[:, 1:8, 4]) for u in sol_vw.u]
WV2 = [sum(u[:, 1:8, 5]) for u in sol_vw.u]
##
plot(R ./ sum(model.N))
plot(W)
plot((R .+ W .+ W2) / sum(model.N), lab="")
plot(E)
plot(A, lab="")
plot(W2 ./ sum(model.N))

plot(diff(C))
plot(V1, lab="")
plot(V2, lab="")
plot(WV1, lab="")
plot(WV2, lab="")
plot((V1 .+ V2 .+ WV1 .+ WV2) ./ sum(model.N[2:6]))
plot((V2 .+ V1) ./ sum(model.N[2:6]))

##
# findall(sol_ie.u[end] .< 0)
##
# plot(diff(C), lab="All infections - no vaccine")
plot(C, lab="All infections - 50% rapid vaccine")
# plot!(diff(C2), lab="Dis risk weighted infections")
# plot!(diff(F), lab="First time infections")
plot!(C, lab="All infections - 70% rapid vaccine")
plot!(xlims=(30, 800),
    ylabel="Daily infections",
    xlabel="Days after 1/12/20",
    title="Nairobi infections with waning immunity",
    legend=:bottomright)
plot!(sum(ι_all, dims=[2, 3])[:])

plot!(sum(ι_dis_weighted, dims=[2, 3])[:])
plot!(sum(ι_all2, dims=[2, 3])[:])

vline!([349], lab="Variant introduction")
vline!([427], lab="")
vline!([230], lab="First jabs start")
##
@time infs = project_daily_infections_by_age_and_vaccine_waning_vac(θs[2], model, Date(2023, 6, 1), cb_imm_esc_cvode, vac_cbs[7], solver)
ι_all2 = sum(infs[1], dims=[2, 3])[:]
ι_dis_weighted2 = sum(infs[2], dims=[2, 3])[:]
plot!(ι_all2)
plot!(ι_dis_weighted2)

##
θs = model |> flatten_chain
infections_by_age_dose = map((θ, VE_acq, VE_inf, VE_sev) -> project_daily_infections_by_age_and_vaccine_waning_vac(θ, model, Date(2023, 6, 1), cb_imm_esc_cvode, vac_cbs[7], solver, VE_acq, VE_inf, VE_sev),
    θs, VE_acqs, VE_infs, VE_sevs)
# infections_by_age_dose = θs .|> θ -> project_daily_infections_by_age_and_vaccine_waning_vac(θ, model, Date(2023, 6, 1), cb_imm_esc_cvode, vac_cbs[7], solver)

##

infs_gathered = gather_infection_type_for_econ_analysis_waning_vac(model, Date(2023, 6, 1), cb_imm_esc_cvode, vac_cbs[1], solver, VE_acqs, VE_infs, VE_sevs, VE_deaths;
    mort_scale=mort_scales[30], mort_age=mort_age, mort_var=mort_var,
    H_ICU_to_death_fit=H_ICU_to_death_fit, H_hosp=H_hosp, H_var_ab=H_var_ab, H_var_delta=H_var_delta)

##
A_pred = get_credible_intervals(sum(infs_gathered.asymptomatic_infs, dims=2)[:, 1, :])
M_pred = get_credible_intervals(sum(infs_gathered.mild_infections, dims=2)[:, 1, :])
sev_pred = get_credible_intervals(sum(infs_gathered.severe_infs, dims=2)[:, 1, :])
crit_pred = get_credible_intervals(sum(infs_gathered.critical_infs, dims=2)[:, 1, :])
deaths_pred = get_credible_intervals(sum(infs_gathered.deadly_infs, dims=2)[:, 1, :])

plot!(A_pred.pred, ribbon=(A_pred.lb, A_pred.ub))
plot(deaths_pred.pred, ribbon=(deaths_pred.lb, deaths_pred.ub))

##
scenarios_names = ["no_vac", "30_perc", "50_perc", "70_perc", "30_perc_rapid", "50_perc_rapid", "70_perc_rapid"]

for scen_num = 1:7
    println("Beginning waning vac. eff. scenario $(scen_num)")
    for k = 1:length(model_vect)
        perc_over_18 = sum(N_kenya[5:end, model.areaname]) / sum(N_kenya[5:end, :])
        KenyaCoVaccines.change_prob_immune_escape_with_waning!(model_vect[k], perc_over_18; startdate=Date(2020, 12, 1), enddate=Date(2023, 12, 1))
    end

    infs_types = mapreduce((model, mort_scale) -> gather_infection_type_for_econ_analysis_waning_vac(model, Date(2023, 6, 1), cb_imm_esc_cvode, vac_cbs[scen_num], solver, VE_acqs, VE_infs, VE_sevs, VE_deaths;
            mort_scale=mort_scale, mort_age=mort_age, mort_var=mort_var,
            H_ICU_to_death_fit=H_ICU_to_death_fit, H_hosp=H_hosp, H_var_ab=H_var_ab, H_var_delta=H_var_delta),
        add_infs, model_vect, mort_scales)

    kenya_inc_A = infs_types.asymptomatic_infs
    kenya_inc_M = infs_types.mild_infections
    kenya_inc_sev = infs_types.severe_infs
    kenya_inc_crit = infs_types.critical_infs
    kenya_inc_deaths = infs_types.deadly_infs

    savefilename_inc_A = "generation_of_data_for_econ_analysis/waning_imm_scenario_$(scen_num)/kenya_inc_A_" * scenarios_names[scen_num] * "_w_vac.rda"
    savefilename_inc_M = "generation_of_data_for_econ_analysis/waning_imm_scenario_$(scen_num)/kenya_inc_M_" * scenarios_names[scen_num] * "_w_vac.rda"
    savefilename_inc_sev = "generation_of_data_for_econ_analysis/waning_imm_scenario_$(scen_num)/kenya_inc_sev_" * scenarios_names[scen_num] * "_w_vac.rda"
    savefilename_inc_crit = "generation_of_data_for_econ_analysis/waning_imm_scenario_$(scen_num)/kenya_inc_crit_" * scenarios_names[scen_num] * "_w_vac.rda"
    savefilename_inc_deaths = "generation_of_data_for_econ_analysis/waning_imm_scenario_$(scen_num)/kenya_inc_deaths_" * scenarios_names[scen_num] * "_w_vac.rda"

    @rput kenya_inc_A kenya_inc_M kenya_inc_sev kenya_inc_crit kenya_inc_deaths
    @rput savefilename_inc_A savefilename_inc_M savefilename_inc_sev savefilename_inc_crit savefilename_inc_deaths

    R"""
        save(kenya_inc_A,file = savefilename_inc_A)
        save(kenya_inc_M,file = savefilename_inc_M)
        save(kenya_inc_sev,file = savefilename_inc_sev)
        save(kenya_inc_crit,file = savefilename_inc_crit)
        save(kenya_inc_deaths,file = savefilename_inc_deaths)
    """
end

##Test randomized draws

R"""
load(file = "generation_of_data_for_econ_analysis/waning_imm_scenario_6/kenya_inc_A_50_perc_rapid_w_vac.rda")
"""
kenya_inc_A_sc_6 = @rget kenya_inc_A
kenya_inc_A_sc_6 = sum(kenya_inc_A_sc_6, dims=2)[:, 1, :]
kenya_inc_A = 0
R"""
kenya_inc_A = 0 #Free memory
"""
R"""
load(file = "generation_of_data_for_econ_analysis/waning_imm_scenario_7/kenya_inc_A_70_perc_rapid_w_vac.rda")
"""
kenya_inc_A_sc_7 = @rget kenya_inc_A
kenya_inc_A_sc_7 = sum(kenya_inc_A_sc_7, dims=2)[:, 1, :]
kenya_inc_A = 0
R"""
kenya_inc_A = 0 #Free memory
"""

diff_A = (sum(kenya_inc_A_sc_6, dims=1).-sum(kenya_inc_A_sc_7, dims=1))[:]
boxplot(diff_A)

pred_A_sc_6 = get_credible_intervals(kenya_inc_A_sc_6)
pred_A_sc_7 = get_credible_intervals(kenya_inc_A_sc_7)
pred_A_diff_cum = get_credible_intervals(cumsum(kenya_inc_A_sc_6 .- kenya_inc_A_sc_7, dims=1))
pred_A_diff = get_credible_intervals(kenya_inc_A_sc_6 .- kenya_inc_A_sc_7)

plot(pred_A_sc_6.pred, ribbon=(pred_A_sc_6.lb, pred_A_sc_6.ub),
    lab="50% rapid", title="Asymptomatic infections", fillalpha=0.1)
plot!(pred_A_sc_7.pred, ribbon=(pred_A_sc_7.lb, pred_A_sc_7.ub),
    lab="70% rapid", fillalpha=0.1)

plot(pred_A_diff.pred, ribbon=(pred_A_diff.lb, pred_A_diff.ub),
    lab="", title="Daily Difference between 50% and 70%", fillalpha=0.2, lw=3)

## No waning

R"""
load(file = "generation_of_data_for_econ_analysis/immune_escape_scenario_6/kenya_inc_A_50_perc_rapid_imm_esc.rda")
"""
kenya_inc_A_sc_6_nw = @rget kenya_inc_A
kenya_inc_A_sc_6_nw = sum(kenya_inc_A_sc_6_nw, dims=2)[:, 1, :]
kenya_inc_A = 0
R"""
kenya_inc_A = 0 #Free memory
"""
R"""
load(file = "generation_of_data_for_econ_analysis/immune_escape_scenario_7/kenya_inc_A_70_perc_rapid_imm_esc.rda")
"""
kenya_inc_A_sc_7_nw = @rget kenya_inc_A
kenya_inc_A_sc_7_nw = sum(kenya_inc_A_sc_7_nw, dims=2)[:, 1, :]
kenya_inc_A = 0
R"""
kenya_inc_A = 0 #Free memory
"""

pred_A_sc_6_nw = get_credible_intervals(kenya_inc_A_sc_6_nw)
pred_A_sc_7_nw = get_credible_intervals(kenya_inc_A_sc_7_nw)
pred_A_diff_cum_nw = get_credible_intervals(cumsum(kenya_inc_A_sc_6_nw .- kenya_inc_A_sc_7_nw, dims=1))
pred_A_diff_nw = get_credible_intervals(kenya_inc_A_sc_6_nw .- kenya_inc_A_sc_7_nw)

plot(pred_A_sc_6_nw.pred, ribbon=(pred_A_sc_6_nw.lb, pred_A_sc_6_nw.ub),
    lab="50% rapid", title="Asymptomatic infections (no vac waning)", fillalpha=0.1)
plot!(pred_A_sc_7_nw.pred, ribbon=(pred_A_sc_7_nw.lb, pred_A_sc_7_nw.ub),
    lab="70% rapid", fillalpha=0.1)

plot(pred_A_diff_nw.pred, ribbon=(pred_A_diff_nw.lb, pred_A_diff_nw.ub),
    lab="", title="Daily Difference between 50% and 70% (no vac waning)", fillalpha=0.2, lw=3)

plot(pred_A_diff_cum_nw.pred, ribbon=(pred_A_diff_cum_nw.lb, pred_A_diff_cum_nw.ub),
    lab="", title="Cum. Difference between 50% and 70% (no vac waning)", fillalpha=0.2, lw=3)


text:
###############################################################################
#
#   FreeSymbol/FreeGroupElem/FreeGroup definition
#
###############################################################################

struct FreeSymbol <: GSymbol
   id::Symbol
   pow::Int
end

FreeGroupElem = GroupWord{FreeSymbol}

mutable struct FreeGroup <: AbstractFPGroup
   gens::Vector{FreeSymbol}

   function FreeGroup(gens::Vector{T}) where {T<:GSymbol}
      G = new(gens)
      G.gens = gens
      return G
   end
end

export FreeGroupElem, FreeGroup

###############################################################################
#
#   Type and parent object methods
#
###############################################################################

elem_type(::Type{FreeGroup}) = FreeGroupElem

parent_type(::Type{FreeGroupElem}) = FreeGroup

###############################################################################
#
#   FreeSymbol constructors
#
###############################################################################

FreeSymbol(s::Symbol) = FreeSymbol(s,1)
FreeSymbol(s::String) = FreeSymbol(Symbol(s))

FreeGroup(n::Int, symbol::String="f") = FreeGroup([Symbol(symbol,i) for i in 1:n])

FreeGroup(a::Vector) = FreeGroup(FreeSymbol.(a))

###############################################################################
#
#   Parent object call overloads
#
###############################################################################

function (G::FreeGroup)()
   id = FreeGroupElem(FreeSymbol[])
   id.parent = G
   return id
end

function (G::FreeGroup)(w::GroupWord{FreeSymbol})
   if length(w) > 0
      for s in w.symbols
         i = findfirst(g -> g.id == s.id, G.gens)
         i == 0 && throw(DomainError(
            "Symbol $s does not belong to $G."))
         s.pow % G.gens[i].pow == 0 || throw(DomainError(
            "Symbol $s doesn't belong to $G."))
      end
   end
   w.parent = G
   return w
end

(G::FreeGroup)(s::FreeSymbol) = G(FreeGroupElem(s))

###############################################################################
#
#   Basic manipulation
#
###############################################################################

hash(s::FreeSymbol, h::UInt) = hash(s.id, hash(s.pow, hash(FreeSymbol, h)))

change_pow(s::FreeSymbol, n::Int) = FreeSymbol(s.id, n)

length(s::FreeSymbol) = abs(s.pow)

###############################################################################
#
#   String I/O
#
###############################################################################

function show(io::IO, G::FreeGroup)
   print(io, "Free group on $(length(G.gens)) generators: ")
   join(io, G.gens, ", ")
end

###############################################################################
#
#   Comparison
#
###############################################################################

###############################################################################
#
#   Inversion
#
###############################################################################

inv(s::FreeSymbol) = change_pow(s, -s.pow)


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

using Markdown
using InteractiveUtils

# ╔═╡ 2dd190ae-3888-11eb-2e68-090f885617d6
using DataStructures

# ╔═╡ aa862a18-3bc3-11eb-17a4-a15284e47f56
using LightGraphs

# ╔═╡ 47c147e6-3888-11eb-1574-cd2d6aea7b1c
file = open(f->read(f, String), "day_11.txt")

# ╔═╡ 019df412-388a-11eb-17c1-2bfe6376975b
#first match get command, second get value of comand
regex = r"([.#L]+)"

# ╔═╡ d29eb8e8-3bc4-11eb-19db-213471925ca2
input =[x.captures[1] for x in eachmatch(regex,file)]

# ╔═╡ 428e6964-3bc7-11eb-2bde-9b9d607bec8b
input_size = (size(input)[1], length(input[1]))

# ╔═╡ 9ec584ca-3bc4-11eb-3d79-6bbd386e9715
begin
	seats = DefaultDict('.')
	for (index_y,line) in enumerate(input)
		for (index_x,item) in enumerate(line)
			seats[(index_y,index_x)] = item
		end
	end
end

# ╔═╡ 8c2fbd94-3bc4-11eb-1c00-136e0add7fcb
seats

# ╔═╡ 6382c796-3bc5-11eb-087a-cde7d6b75abb
begin
	function addjacenty(pos::Tuple{Int64,Int64},dict::DefaultDict)
		neig_pos = [(-1,-1),(-1,0),(-1,1),(1,-1),(1,0),(1,1),(0,-1),(0,1)]
		return [dict[(pos[1]+x,pos[2]+y)] for (x,y) in neig_pos]
	end

	is_floor(pos::Tuple{Int64,Int64},dict::DefaultDict) = '.'

	function is_empty(pos::Tuple{Int64,Int64},dict::DefaultDict)
		if all([(x!= '#') for x in addjacenty(pos,dict)])
			return '#'
		else
			return 'L'
		end
	end

	function is_occupy(pos::Tuple{Int64,Int64},dict::DefaultDict)
		if sum([(x== '#') for x in addjacenty(pos,dict)]) ≥ 4
			return 'L'
		else
			return '#'
		end
	end

	one_step_seat_func = Dict([('#',is_occupy),('L',is_empty),('.',is_floor)])
	function one_step_seat(pos::Tuple{Int64,Int64},dict::DefaultDict)
		old_value = dict[pos]

		new_value = one_step_seat_func[old_value](pos,dict)
		return new_value
	end
end

# ╔═╡ 33fcb568-3bc5-11eb-16b6-97ee834bb5cf
function one_step_plane(dict)

	change = false
	new_seats = DefaultDict('.')

	for y in 1:input_size[1]
		for x in 1:input_size[2]
			pos = (y,x)
			new_seats[pos] = one_step_seat(pos,dict)

			if(new_seats[pos] != dict[pos])
				change = true
			end
		end
	end

	return new_seats,change
end

# ╔═╡ 0f1dc8d6-3bca-11eb-048c-b3dd68dee015

begin
	seats1 = copy(seats)
	change = true
	while(change == true)
		(seats1,change) = one_step_plane(seats1)
	end

	sum([(value == '#') for (_, value) in seats1])

end

#29 + 22 + 49

# ╔═╡ 7902a486-3bcc-11eb-0691-b9eedcc0b2c1
md"# Question 2"

# ╔═╡ f4b9ed26-3bd8-11eb-39ab-7d627dfa0712
push!([[5]],[3])

# ╔═╡ 7e524662-3bcc-11eb-28c9-45b61768ed7d
begin
	function addjacenty2(pos::Tuple{Int64,Int64},dict::DefaultDict)
		neig_vert = [(y,pos[2]) for y in 1:input_size[1]]
		neig_hor = [(pos[1],x) for x in 1:input_size[2]]
		max_offset = min(input_size[1],input_size[2])
		neig_dir1 = [(pos[1]+q,pos[2]+q) for q in -max_offset:max_offset]
		neig_dir2 = [(pos[1]-q,pos[2]-q) for q in -max_offset:max_offset]

		resp = []
		push!(resp,neig_vert)
		push!(resp,neig_hor)
		push!(resp,neig_dir1)
		push!(resp,neig_dir2)

		unique!(resp)
		filter!(e->e !=pos,resp)
		return [dict[(pos[1]+x,pos[2]+y)] for (x,y) in neig_pos]
	end

	function is_empty2(pos::Tuple{Int64,Int64},dict::DefaultDict)
		if all([(x!= '#') for x in addjacenty2(pos,dict)])
			return '#'
		else
			return 'L'
		end
	end

	is_floor2 = is_floor

	function is_occupy2(pos::Tuple{Int64,Int64},dict::DefaultDict)

		neig_vert = [(y,pos[2]) for y in 1:input_size[1]]
		neig_hor = [(pos[1],x) for x in 1:input_size[2]]
		max_offset = min(input_size[1],input_size[2])
		neig_dir1 = [(pos[1]+q,pos[2]+q) for q in -max_offset:max_offset]
		neig_dir2 = [(pos[1]-q,pos[2]-q) for q in -max_offset:max_offset]

		resp = []
		append!(resp,neig_vert)
		append!(resp,neig_hor)
		append!(resp,neig_dir1)
		append!(resp,neig_dir2)

		unique!(resp)
		filter!(e->e !=pos,resp)

		if sum([(x== '#') for x in addjacenty2(pos,dict)]) ≥ 5
			return 'L'
		else
			return '#'
		end
	end

	one_step_seat_func2 = Dict([('#',is_occupy2),('L',is_empty2),('.',is_floor2)])
	function one_step_seat2(pos::Tuple{Int64,Int64},dict::DefaultDict)
		old_value = dict[pos]

		new_value = one_step_seat_func2[old_value](pos,dict)
		return new_value
	end
end

# ╔═╡ f56f80ce-3bcd-11eb-3cac-4d18887d9b36
function one_step_plane2(dict)

	change = false
	new_seats = DefaultDict('.')

	for y in 1:input_size[1]
		for x in 1:input_size[2]
			pos = (y,x)
			new_seats[pos] = one_step_seat2(pos,dict)

			if(new_seats[pos] != dict[pos])
				change = true
			end
		end
	end

	return new_seats,change
end

# ╔═╡ ea5878ba-3bcf-11eb-11fc-6f7fa62e6e88
begin
	seats2 = copy(seats)
	change2 = true
	while(change2 == true)
		(seats2,change2) = one_step_plane2(seats2)
	end

	sum([(value == '#') for (_, value) in seats2])
end

# ╔═╡ Cell order:
# ╠═2dd190ae-3888-11eb-2e68-090f885617d6
# ╠═aa862a18-3bc3-11eb-17a4-a15284e47f56
# ╠═47c147e6-3888-11eb-1574-cd2d6aea7b1c
# ╠═019df412-388a-11eb-17c1-2bfe6376975b
# ╠═d29eb8e8-3bc4-11eb-19db-213471925ca2
# ╠═428e6964-3bc7-11eb-2bde-9b9d607bec8b
# ╠═9ec584ca-3bc4-11eb-3d79-6bbd386e9715
# ╠═8c2fbd94-3bc4-11eb-1c00-136e0add7fcb
# ╠═6382c796-3bc5-11eb-087a-cde7d6b75abb
# ╠═33fcb568-3bc5-11eb-16b6-97ee834bb5cf
# ╠═0f1dc8d6-3bca-11eb-048c-b3dd68dee015
# ╠═7902a486-3bcc-11eb-0691-b9eedcc0b2c1
# ╠═f4b9ed26-3bd8-11eb-39ab-7d627dfa0712
# ╠═7e524662-3bcc-11eb-28c9-45b61768ed7d
# ╠═f56f80ce-3bcd-11eb-3cac-4d18887d9b36
# ╠═ea5878ba-3bcf-11eb-11fc-6f7fa62e6e88


text:
function foo(x::Int64, y) 42 end
function foo(x, y::Int64) 1729 end
println(foo(1,2))

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
@testset "plotting" begin
    x = rand(10)
    f = GP(SqExponentialKernel())
    gp = f(x, 0.1)

    plt1 = sampleplot(gp, 10)
    @test plt1.n == 10

    # Check recipe dispatches for `FiniteGP`s
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), gp)
    @test length(rec) == 1 && length(rec[1].args) == 2 # one series with two arguments
    @test rec[1].args[1] == x
    @test rec[1].args[2] == gp
    @test isempty(rec[1].plotattributes) # no default attributes

    z = 1 .+ x
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), z, gp)
    @test length(rec) == 1 && length(rec[1].args) == 2 # one series with two arguments
    @test rec[1].args[1] == z
    @test rec[1].args[2] == zero(x)
    # 3 default attributes
    attributes = rec[1].plotattributes
    @test sort!(collect(keys(attributes))) == [:fillalpha, :linewidth, :ribbon]
    @test attributes[:fillalpha] == 0.3
    @test attributes[:linewidth] == 2
    @test attributes[:ribbon] == sqrt.(cov_diag(gp))

    # Check recipe dispatches for `AbstractGP`s
    # with `AbstractVector` and `AbstractRange`:
    for x in (rand(10), 0:0.01:1)
        rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), f, x)
        @test length(rec) == 1 && length(rec[1].args) == 1 # one series with one argument
        @test rec[1].args[1] isa AbstractGPs.FiniteGP
        @test rec[1].args[1].x == x
        @test rec[1].args[1].f == f
        @test isempty(rec[1].plotattributes) # no default attributes

        z = 1 .+ x
        rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), z, f, x)
        @test length(rec) == 1 && length(rec[1].args) == 2 # one series with two arguments
        @test rec[1].args[1] == z
        @test rec[1].args[2] isa AbstractGPs.FiniteGP
        @test rec[1].args[2].x == x
        @test rec[1].args[2].f == f
        @test isempty(rec[1].plotattributes) # no default attributes
    end

    # with minimum and maximum:
    xmin = rand()
    xmax = 4 + rand()
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), f, xmin, xmax)
    @test length(rec) == 1 && length(rec[1].args) == 1 # one series with one argument
    @test rec[1].args[1] isa AbstractGPs.FiniteGP
    @test rec[1].args[1].x == range(xmin, xmax; length=1_000)
    @test rec[1].args[1].f == f
    @test isempty(rec[1].plotattributes) # no default attributes

    z = range(0, 1; length=1_000)
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), z, f, xmin, xmax)
    @test length(rec) == 1 && length(rec[1].args) == 2 # one series with two arguments
    @test rec[1].args[1] == z
    @test rec[1].args[2] isa AbstractGPs.FiniteGP
    @test rec[1].args[2].x == range(xmin, xmax; length=1_000)
    @test rec[1].args[2].f == f
    @test isempty(rec[1].plotattributes) # no default attributes

    # Check dimensions
    @test_throws DimensionMismatch plot(rand(5), gp)
end


text:
export CLT, learn_chow_liu_tree, parent_vector
using Graphs: SimpleGraph, SimpleDiGraph, complete_graph, kruskal_mst,
    bfs_tree, center, connected_components, induced_subgraph, src, dst
using SimpleWeightedGraphs: SimpleWeightedGraph
using MetaGraphs: MetaDiGraph, set_prop!, props, nv, ne, edges, add_edge!, vertices

#####################
# Learn a Chow-Liu tree from (weighted) data
#####################

"""
Chow-Liu Tree
"""
const CLT = MetaDiGraph

"""
learn a Chow-Liu tree from training set `train_x`, with Laplace smoothing factor `α`, specifying the tree root by `clt_root`
return a `CLT`
"""
function learn_chow_liu_tree(train_x; α = 1.0, clt_root="graph_center",
        weight=ones(Float64, num_examples(train_x)))::CLT
    features_num = num_features(train_x)

    # calculate mutual information
    (dis_cache, MI) = mutual_information(train_x, weight; α = α)

    # maximum spanning tree/ forest
    g = SimpleWeightedGraph(complete_graph(features_num))
    mst_edges = kruskal_mst(g,- MI)
    tree = SimpleGraph(features_num)
    map(mst_edges) do edge
        add_edge!(tree, src(edge), dst(edge))
    end

    # Build rooted tree / forest
    if clt_root == "graph_center"
        clt = SimpleDiGraph(features_num)
        for c in filter(c -> (length(c) > 1), connected_components(tree))
            sg, vmap = induced_subgraph(tree, c)
            sub_root = vmap[center(sg)[1]]
            clt = union(clt, bfs_tree(tree, sub_root))
        end
    elseif clt_root == "rand"
        roots = [rand(c) for c in connected_components(tree)]
        clt = SimpleDiGraph(features_num)
        for root in roots clt = union(clt, bfs_tree(tree, root)) end
    else
        error("Cannot learn CLT with root $clt_root")
    end

    clt = MetaDiGraph(clt)
    parent = parent_vector(clt)
    for (c, p) in enumerate(parent)
        set_prop!(clt, c, :parent, p)
    end

    for v in vertices(clt)
        p = parent[v]
        cpt_matrix = get_cpt(p, v, dis_cache)
        set_prop!(clt, v, :cpt, cpt_matrix)
    end

    return clt
end

"""
Calculate CPT of `child` conditioned on `parent` from `dis_cache`
"""
function get_cpt(parent, child, dis_cache)
    if parent == 0
        p = dis_cache.marginal[child, :]
        return Dict(0=>p[1], 1=>p[2])
    else
        p = dis_cache.pairwise[child, parent, :] ./ [dis_cache.marginal[parent, :]; dis_cache.marginal[parent, :]]
        @. p[isnan(p)] = 0; @. p[p==Inf] = 0; @. p[p == -Inf] = 0
        return Dict((0,0)=>p[1], (1,0)=>p[3], (0,1)=>p[2], (1,1)=>p[4]) #p(child|parent)
    end
end

"Get parent vector of a tree"
function parent_vector(tree::CLT)::Vector{Int64}
    v = zeros(Int64, nv(tree)) # parent of roots is 0
    foreach(e->v[dst(e)] = src(e), edges(tree))
    return v
end

import LogicCircuits: print_tree
"Print edges and vertices of a ChowLiu tree"
function print_tree(clt::CLT)
    for e in edges(clt) print(e); print(" ");end
    if clt isa MetaDiGraph
        for v in vertices(clt) print(v); print(" "); println(props(clt, v)) end
    end
end