text:
using Oceananigans.Operators

# Define aliases for some spacings and areas at all locations.
for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ)

    x_spacing_alias = Symbol(:Δx, LX, LY, LZ)
    y_spacing_alias = Symbol(:Δy, LX, LY, LZ)
    z_spacing_alias = Symbol(:Δz, LX, LY, LZ)

    z_area_alias = Symbol(:Az, LX, LY, LZ)

    x_spacing_function = Symbol(:Δx, LX, LY, :ᵃ)
    y_spacing_function = Symbol(:Δy, LX, LY, :ᵃ)
    z_spacing_function = Symbol(:Δz, :ᵃᵃ, LZ)

    z_area_function = Symbol(:Az, LX, LY, :ᵃ)

    @eval begin
        const $x_spacing_alias = $x_spacing_function
        const $y_spacing_alias = $y_spacing_function
        const $z_spacing_alias = $z_spacing_function
        const $z_area_alias = $z_area_function
    end
end

abstract type AbstractGridMetric end

struct XSpacingMetric <: AbstractGridMetric end
struct YSpacingMetric <: AbstractGridMetric end
struct ZSpacingMetric <: AbstractGridMetric end

metric_function_prefix(::XSpacingMetric) = :Δx
metric_function_prefix(::YSpacingMetric) = :Δy
metric_function_prefix(::ZSpacingMetric) = :Δz

struct XAreaMetric <: AbstractGridMetric end
struct YAreaMetric <: AbstractGridMetric end
struct ZAreaMetric <: AbstractGridMetric end

metric_function_prefix(::XAreaMetric) = :Ax
metric_function_prefix(::YAreaMetric) = :Ay
metric_function_prefix(::ZAreaMetric) = :Az

struct VolumeMetric <: AbstractGridMetric end

metric_function_prefix(::VolumeMetric) = :V

# Convenient instances for users
const Δx = XSpacingMetric()
const Δy = YSpacingMetric()

"""
    Δz = ZSpacingMetric()

Instance of `ZSpacingMetric` that generates `BinaryOperation`s
between `AbstractField`s and the vertical grid spacing evaluated
at the same location as the `AbstractField`.

`Δx` and `Δy` play a similar role for horizontal grid spacings.

Example
=======

```jldoctest
julia> using Oceananigans

julia> using Oceananigans.AbstractOperations: Δz

julia> grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 2, 3)); c = CenterField(CPU(), grid);

julia> c_dz = c * Δz # returns BinaryOperation between Field and GridMetricOperation
BinaryOperation at (Center, Center, Center)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=1, Ny=1, Nz=1)
│   └── domain: x ∈ [0.0, 1.0], y ∈ [0.0, 2.0], z ∈ [-3.0, 0.0]
└── tree:
    * at (Center, Center, Center)
    ├── Field located at (Center, Center, Center)
    └── Δzᵃᵃᶜ at (Center, Center, Center)

julia> c .= 1;

julia> c_dz[1, 1, 1]
3.0
```
"""
const Δz = ZSpacingMetric()

const Ax = XAreaMetric()
const Ay = YAreaMetric()
const Az = ZAreaMetric()

"""
    volume = VolumeMetric()

Instance of `VolumeMetric` that generates `BinaryOperation`s
between `AbstractField`s and their cell volumes. Summing
this `BinaryOperation` yields an integral of `AbstractField`
over the domain.

Example
=======

```jldoctest
julia> using Oceananigans

julia> using Oceananigans.AbstractOperations: volume

julia> grid = RegularRectilinearGrid(size=(2, 2, 2), extent=(1, 2, 3)); c = CenterField(CPU(), grid);

julia> c .= 1;

julia> c_dV = c * volume
BinaryOperation at (Center, Center, Center)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=2, Ny=2, Nz=2)
│   └── domain: x ∈ [0.0, 1.0], y ∈ [0.0, 2.0], z ∈ [-3.0, 0.0]
└── tree:
    * at (Center, Center, Center)
    ├── Field located at (Center, Center, Center)
    └── Vᶜᶜᶜ at (Center, Center, Center)

julia> c_dV[1, 1, 1]
0.75

julia> sum(c_dV)
6.0
```
"""
const volume = VolumeMetric()

"""
    metric_function(loc, metric::AbstractGridMetric)

Returns the function associated with `metric::AbstractGridMetric`
at `loc`ation.
"""
function metric_function(loc, metric::AbstractGridMetric)
    code = Tuple(interpolation_code(ℓ) for ℓ in loc)
    prefix = metric_function_prefix(metric)
    metric_function_symbol = Symbol(prefix, code...)
    return eval(metric_function_symbol)
end

struct GridMetricOperation{X, Y, Z, A, G, T, M} <: AbstractOperation{X, Y, Z, A, G, T}
          metric :: M
            grid :: G
    architecture :: A

    function GridMetricOperation{X, Y, Z}(metric::M, grid::G) where {X, Y, Z, M, G}
        arch = architecture(grid)
        A = typeof(arch)
        T = eltype(grid)
        return new{X, Y, Z, A, G, T, M}(metric, grid, arch)
    end
end

@inline Base.getindex(gm::GridMetricOperation, i, j, k) = gm.metric(i, j, k, gm.grid)

# Special constructor for BinaryOperation
GridMetricOperation(L, metric, grid) = GridMetricOperation{L[1], L[2], L[3]}(metric_function(L, metric), grid)


text:
const Vec2 = SVector{2,Float64}

abstract type AbstractLD2 <: POMDP{Vec2, Vec2, Vec2} end

"""
2-Dimensional light-dark problem

    min_noise::Float64
        minimum standard deviation of observation noise

    min_noise_loc::Float64
        x location of minimum observation noise

    Q::Matrix{Float64}
    R::Matrix{Float64}
        cost function matrices: reward = -s'Qs -a'Ra

    term_radius::Float64
        problem will terminate if the state is within this radius of zero

    init_dist::Any
        initial state distribution

    discount::Float64
        discount factor
"""
@with_kw type LightDark2D <: AbstractLD2
    min_noise::Float64      = 0.0
    min_noise_loc::Float64  = 5.0
    Q::Matrix{Float64}      = diagm([0.5, 0.5])
    R::Matrix{Float64}      = diagm([0.5, 0.5])
    term_radius::Float64    = 1e-5
    init_dist::Any          = SymmetricNormal2([2.0, 2.0], 0.5)
    discount::Float64       = 1.0
    count::Int              = 0
end

# POMDPs.jl API functions:
generate_s(p::AbstractLD2, s::Vec2, a::Vec2, rng::AbstractRNG) = generate_s(p, s, a)
generate_o(p::AbstractLD2, s::Vec2, a::Vec2, sp::Vec2, rng::AbstractRNG) = generate_o(p, sp, rng)
observation(p::AbstractLD2, a::Vec2, sp::Vec2) = observation(p, sp)
isterminal(p::AbstractLD2, s::Vec2) = norm(s) <= p.term_radius
initial_state_distribution(p::AbstractLD2) = p.init_dist
reward(p::AbstractLD2, s::Vec2, a::Vec2, sp::Vec2) = -(dot(s, p.Q*s) + dot(a, p.R*a))
discount(p::AbstractLD2) = p.discount

immutable SymmetricNormal2
    mean::Vec2
    std::Float64
end
rand(rng::AbstractRNG, d::SymmetricNormal2) = d.mean + d.std*Vec2(randn(rng, 2))
pdf(d::SymmetricNormal2, s::Vec2) = exp(-0.5*sum((s-d.mean).^2)/d.std^2)/(2*pi*d.std^2)
mean(d::SymmetricNormal2) = d.mean
mode(d::SymmetricNormal2) = d.mean
Base.eltype(::Type{SymmetricNormal2}) = Vec2

# chose this on 2/6/17 because I like the bowtie particle patterns it produces
# unclear which one was actually used in the paper
# Masters Thesis by Pas assumes the sqrt version
obs_std(p::AbstractLD2, x::Float64) = sqrt(0.5*(p.min_noise_loc-x)^2 + p.min_noise)
# obs_std(p::AbstractLD2, x::Float64) = 0.5*(p.min_noise_loc-x)^2 + p.min_noise

function generate_s(p::AbstractLD2, s::Vec2, a::Vec2)
    p.count += 1
    return s+a
end
observation(p::AbstractLD2, sp::Vec2) = SymmetricNormal2(sp, obs_std(p, sp[1]))
generate_o(p::AbstractLD2, sp::Vec2, rng::AbstractRNG) = rand(rng, observation(p, sp))

ext: jl


text:
# Autogenerated wrapper script for UDUNITS_jll for i686-w64-mingw32
export libudunits2, udunits2

using Expat_jll
JLLWrappers.@generate_wrapper_header("UDUNITS")
JLLWrappers.@declare_library_product(libudunits2, "libudunits2.dll")
JLLWrappers.@declare_executable_product(udunits2)
function __init__()
    JLLWrappers.@generate_init_header(Expat_jll)
    JLLWrappers.@init_library_product(
        libudunits2,
        "bin\\libudunits2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        udunits2,
        "bin\\udunits2.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ParetoSets
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
let
    x = [0.5, 0.5]
    y = [1.0, 2.0]
    z = [1.5, 1.5]

    @test dominates(x, y)
    @test !dominates(y, x)
    @test dominates(x, z)
    @test !dominates(z, x)
    @test !dominates(y, z)
    @test !dominates(z, y)
end

let
    x = [0.01, 0.02]
    y = [0.02, 0.01]
    v = Vector{Float64}[]
    push!(v, x)
    push!(v, y)
    for i = 1:8
        push!(v, rand(2) .+ 0.02)
    end
    ids = naive_pareto(v)
    @test ids == [1, 2]

    ps = ParetoSet()
    push!(ps, v[1])
    @test pareto_ids(ps) == [1]
    @test pareto_ys(ps) == [x]
    push!(ps, v[2])
    @test pareto_ids(ps) == [1,2]
    @test pareto_ys(ps) == [x,y]
    append!(ps, v[3:end])
    @test pareto_ids(ps) == [1,2]
    @test pareto_ys(ps) == [x,y]
    @test pareto_ids(ps) == naive_pareto(v)
    @test dominated_ids(ps) == collect(3:length(v))
    @test dominated_ys(ps) == v[3:end]
end


text:
using JLD
using Distances

if !isdefined(:loaded) || !loaded
    println("Loading...")
    # jldopen("85.jld") do file
    jldopen("90.jld") do file
        global X, vocab, vocab2idx, idx2vocab
        X = read(file, "X")
        vocab = read(file, "vocab")
        vocab2idx = read(file, "vocab2idx")
        idx2vocab = read(file, "idx2vocab")
        global loaded = true
    end
end

# note: depends some globals
function similar_words(word; K=10, should_lower=true)
    if should_lower
        word = lowercase(word)
    end
    v = X[vocab2idx[word],:]
    distances = zeros(length(vocab))
    for i in 1:length(vocab)
        distances[i] = cosine_dist(v, X[i,:])
    end

    p = sortperm(distances)
    words = []
    for i in 1:K
        push!(words, (idx2vocab[p[i]], distances[p[i]]))
    end

    words
end

words = similar_words("England", K=11)
# words = similar_words("Spain", K=11)

for (idx, (word, score)) in enumerate(words)
    println("$idx: ($word, $score)")
end

#=
julia> @time include("./88.jl")
1: (England, 0.0)
2: (Scotland, 0.3268883016363757)
3: (Wales, 0.3797738775963655)
4: (Australia, 0.38495262864152646)
5: (Italy, 0.4457559448475833)
6: (Ireland, 0.45570685913751796)
7: (France, 0.47124095457277393)
8: (Netherlands, 0.49150145114866073)
9: (Spain, 0.4923261574483986)
10: (match, 0.4996265985763505)
11: (New_Zealand, 0.5014242983908542)
  0.797475 seconds (2.53 M allocations: 1.045 GiB, 9.35% gc time)
=#

# **90**
#=
julia> @time include("./88.jl")
1: (England, 0.0)
2: (Scotland, 0.3663732399867812)
3: (T20, 0.45535642583910907)
4: (Wales, 0.46629287761514726)
5: (Coope, 0.4867497657715091)
6: (Yorkshire, 0.5259334960472704)
7: (Australia, 0.5301869440377063)
8: (Leeds, 0.5345024659948912)
9: (Cardiff, 0.5370676668547449)
10: (match, 0.5375962084158777)
11: (Liverpool, 0.5377926492053353)
  0.593390 seconds (632.09 k allocations: 221.208 MiB, 20.46% gc time)
=#


text:
function get_eop(obj::EopParserLoader)
    return jcall(obj, "getEop", Collection, ())
end

function load_data(obj::EopParserLoader, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end

function still_accepts_data(obj::EopParserLoader)
    return jcall(obj, "stillAcceptsData", jboolean, ())
end


text:
export supported_platforms, platform_key, platform_dlext, valid_dl_path,
       arch, wordsize, triplet, Platform, Linux, MacOS, Windows

abstract type Platform end

struct Linux <: Platform
    arch::Symbol
    libc::Symbol

    function Linux(arch::Symbol, libc::Symbol=:glibc)
        if !in(arch, [:i686, :x86_64, :aarch64, :powerpc64le, :ppc64le, :armv7l])
            throw(ArgumentError("Unsupported architecture '$arch' for Linux"))
        end
        if libc !== :glibc && libc !== :musl
            throw(ArgumentError("Unsupported libc '$libc' for Linux"))
        end
        if arch === :ppc64le
            arch = :powerpc64le
        end
        new(arch, libc)
    end
end

struct MacOS <: Platform

    function MacOS(arch::Symbol)
        if arch !== :x86_64
            throw(ArgumentError("Unsupported architecture '$arch' for macOS"))
        end
        new()
    end
    MacOS() = new()
end

struct Windows <: Platform
    arch::Symbol

    function Windows(arch::Symbol)
        if arch !== :i686 && arch !== :x86_64
            throw(ArgumentError("Unsupported architecture '$arch' for Windows"))
        end
        new(arch)
    end
end

"""
    arch(platform)

Get the architecture for the given `Platform` object as a `Symbol`.

# Examples
```jldoctest
julia> arch(Linux(:aarch64))
:aarch64

julia> arch(MacOS())
:x86_64
```
"""
arch(p::Platform) = p.arch
arch(m::MacOS) = :x86_64

"""
    wordsize(platform)

Get the word size for the given `Platform` object.

# Examples
```jldoctest
julia> wordsize(Linux(:arm7vl))
32

julia> wordsize(MacOS())
64
```
"""
wordsize(l::Linux) = arch(l) === :i686 || arch(l) === :armv7l ? 32 : 64
wordsize(w::Windows) = arch(w) === :i686 ? 32 : 64
wordsize(m::MacOS) = 64

"""
    triplet(platform)

Get the target triplet for the given `Platform` object as a `String`.

# Examples
```jldoctest
julia> triplet(MacOS())
"x86_64-apple-darwin14"

julia> triplet(Windows(:i686))
"i686-w64-mingw32"

julia> triplet(Linux(:armv7l))
"arm-linux-gnueabihf"
```
"""
function triplet(l::Linux)
    c = l.libc === :glibc ? "gnu" : "musl" # Currently only glibc and musl are recognized
    if arch(l) === :armv7l
        string("arm-linux-", c, "eabihf")
    else
        string(arch(l), "-linux-", c)
    end
end
triplet(w::Windows) = string(arch(w), "-w64-mingw32")
triplet(m::MacOS) = "x86_64-apple-darwin14"

"""
    supported_platforms()

Return the list of supported platforms as an array of `Platform`s.
"""
function supported_platforms()
    return [
        Linux(:i686),
        Linux(:x86_64),
        Linux(:aarch64),
        Linux(:armv7l),
        Linux(:powerpc64le),
        MacOS(),
        Windows(:i686),
        Windows(:x86_64),
    ]
end

# Compat doesn't use the Base definitions for whatever terrible reason, so we'll overload
# both, ensuring the user gets our definitions regardless of whether they use Sys.is* or
# Compat.Sys.is*.
#if isdefined(Base.Sys, :isapple)
#    Base.Sys.isapple(p::Platform) = p isa MacOS
#    Base.Sys.islinux(p::Platform) = p isa Linux
#    Base.Sys.iswindows(p::Platform) = p isa Windows
# end
# Compat.Sys.isapple(p::Platform) = p isa MacOS
# Compat.Sys.islinux(p::Platform) = p isa Linux
# Compat.Sys.iswindows(p::Platform) = p isa Windows

"""
    platform_key(machine::AbstractString = Sys.MACHINE)

Returns the platform key for the current platform, or any other though the
the use of the `machine` parameter.
"""
function platform_key(machine::AbstractString = Sys.MACHINE)
    # First, off, if `machine` is literally one of the values of our mapping
    # above, just return the relevant key
    for key in supported_platforms()
        if machine == triplet(key)
            return key
        end
    end

    # Otherwise, try to parse the machine into one of our keys
    if startswith(machine, "x86_64-apple-darwin")
        return MacOS()
    end
    if ismatch(r"x86_64-(pc-)?(unknown-)?linux-gnu", machine)
        return Linux(:x86_64)
    end
    if ismatch(r"i\d86-(pc-)?(unknown-)?linux-gnu", machine)
        return Linux(:i686)
    end
    if ismatch(r"aarch64-(pc-)?(unknown-)?linux-gnu", machine)
        return Linux(:aarch64)
    end
    if ismatch(r"armv7l-(pc-)?(unknown-)?linux-gnueabihf", machine)
        return Linux(:armv7l)
    end
    if ismatch(r"powerpc64le-(pc-)?(unknown-)?linux-gnu", machine)
        return Linux(:powerpc64le)
    end

    throw(ArgumentError("Platform `$(machine)` is not an officially supported platform"))
end

"""
    platform_dlext(platform::Platform = platform_key())

Return the dynamic library extension for the given platform, defaulting to the
currently running platform.  E.g. returns "so" for a Linux-based platform,
"dll" for a Windows-based platform, etc...
"""
platform_dlext(l::Linux) = "so"
platform_dlext(m::MacOS) = "dylib"
platform_dlext(w::Windows) = "dll"
platform_dlext() = platform_dlext(platform_key())

"""
    valid_dl_path(path::AbstractString, platform::Platform)

Return `true` if the given `path` ends in a valid dynamic library filename.
E.g. returns `true` for a path like `"usr/lib/libfoo.so.3.5"`, but returns
`false` for a path like `"libbar.so.f.a"`.
"""
function valid_dl_path(path::AbstractString, platform::Platform)
    dlext_regexes = Dict(
        # On Linux, libraries look like `libnettle.so.6.3.0`
        "so" => r"^(.*).so(\.[\d]+){0,3}$",
        # On OSX, libraries look like `libnettle.6.3.dylib`
        "dylib" => r"^(.*).dylib$",
        # On Windows, libraries look like `libnettle-6.dylib`
        "dll" => r"^(.*).dll$"
    )

    # Given a platform, find the dlext regex that matches it
    dlregex = dlext_regexes[platform_dlext(platform)]

    # Return whether or not that regex matches the basename of the given path
    return ismatch(dlregex, basename(path))
end


text:
struct SparsePathway{I,J,S} <: Pathway
    pre::I      # presynaptic units
    syns::S     # synaptic state, of length sum(length.(targets))
    post::J     # postsynaptic units
    targets::Vector{Vector{Int}} # ragged arrays where targets[pre]
                                 # is a list of post indices.
end

"""
Creates a sparse pathway between `pre` and `post` with connection
probability `p`. The resulting `n` synapses are then instanciated
by calling `fsyn(n)`.
"""
function SparsePathway(pre::Group, post::Group, p::Real, fsyn)
    targets = [find(rand(length(post)) .< p) for _ in 1:length(pre)]
    syns = fsyn(sum(length.(targets)))
    SparsePathway(pre, syns, post, targets)
end

"""
Apply `f` on every connection in the pathway where `select` is true
for the presynaptic element.
"""
function Base.broadcast(f, path::SparsePathway, select=(_)->true)
    k = 0
    for i in eachindex(path.pre)
        if select(Elem(path.pre,i))
            @inbounds @simd for j in eachindex(path.targets[i])
                f(Elem(path.pre, i), Elem(path.syns, k+j),
                  Elem(path.post, path.targets[i][j]))
            end
        end
        k += length(path.targets[i])
    end
    nothing
end


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, NLsolve
	Random.seed!(0)
end;

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing5.08"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	a, b, c = 3, 5, 4
	dist = TriangularDist(a,b,c)
	n = 2000
	samples = rand(dist,n)
end

# ╔═╡ d91f4766-3a9b-11eb-2892-7b7925b084b2
begin
	m_k(k,data) = 1/n*sum(data.^k)
	mHats = [m_k(i,samples) for i in 1:3]
end

# ╔═╡ d91f7dd0-3a9b-11eb-313d-81f3bcee7e68
function equations(F, x)
    F[1] = 1/3*( x[1] + x[2] + x[3] ) - mHats[1]
    F[2] = 1/6*( x[1]^2 + x[2]^2 + x[3]^2 + x[1]*x[2] + x[1]*x[3] +
		 x[2]*x[3] ) - mHats[2]
    F[3] = 1/10*( x[1]^3 + x[2]^3 + x[3]^3 + x[1]^2*x[2] + x[1]^2*x[3] +
		 x[2]^2*x[1] + x[2]^2*x[3] + x[3]^2*x[1] + x[3]^2*x[2] +
		 x[1]*x[2]*x[3] ) - mHats[3]
end

# ╔═╡ d9200de0-3a9b-11eb-2204-db202244e154
begin
	nlOutput = nlsolve(equations, [ 0.1; 0.1; 0.1])
	sol = sort(nlOutput.zero)
	aHat, bHat, cHat = sol[1], sol[3], sol[2]
	Text("Found estimates for (a,b,c) = $((aHat, bHat, cHat))")
	nlOutput
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing5.08"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═d91f4766-3a9b-11eb-2892-7b7925b084b2
# ╠═d91f7dd0-3a9b-11eb-313d-81f3bcee7e68
# ╠═d9200de0-3a9b-11eb-2204-db202244e154
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
function DeInitInformation!(sys::System,comm::MPI.Comm,root::Int64)
    coords_all=SerializeArray(sys.atoms.coords_all)
    MPI.Gatherv!(SerializeArray(sys.atoms.coords),VBuffer(coords_all,sys.mpi.commsizes*3),root,comm)
    DeSerializeAllCoord!(sys,coords_all)

    vels_all=SerializeArray(sys.atoms.vels_all)
    MPI.Gatherv!(SerializeArray(sys.atoms.vels),VBuffer(vels_all,sys.mpi.commsizes*3),root,comm)
    DeSerializeAllVels!(sys,vels_all)

    bonds_all=SerializeArray(sys.atoms.bonds_all)
    MPI.Gatherv!(SerializeArray(sys.atoms.bonds),VBuffer(bonds_all,sys.mpi.bond_commsizes*3),root,comm)
    DeSerializeAllBonds!(sys,bonds_all)

    types_all=SerializeArray(sys.atoms.types_all)
    MPI.Gatherv!(SerializeArray(sys.atoms.types),VBuffer(types_all,sys.mpi.commsizes),root,comm)
    DeSerializeAllTypes!(sys,types_all)

    masses_all=SerializeArray(sys.atoms.masses_all)
    MPI.Gatherv!(SerializeArray(sys.atoms.masses),VBuffer(masses_all,sys.mpi.commsizes),root,comm)
    DeSerializeAllMasses!(sys,masses_all)
end

function RestartInformation(sys::System,comm::MPI.Comm,root::Int64)
    MPI.Barrier(comm)
    if sys.restart!=NoRestart()
        if sys.restart.freq==-1 && sys.current_step[1]==sys.steps+sys.start_step-1
            DeInitInformation!(sys,comm,root)
            if MPI.Comm_rank(comm)==root
                Restart(sys,sys.restart)
            end
        elseif sys.current_step[1]%sys.restart.freq==0 || sys.current_step[1]==sys.steps+sys.start_step-1
            DeInitInformation!(sys,comm,root)
            if MPI.Comm_rank(comm)==root
                Restart(sys,sys.restart)
            end
        end
    end
    MPI.Barrier(comm)
end

function Restart(sys::System,restart::JSONRestart)

    molecules=DeInitMolecules(sys)

    configs=Dict(
        "box"=>[sys.pbc.Lx,sys.pbc.Ly,sys.pbc.Lz],
        "molecules"=>molecules
        )
    #control_file=open("control.json","w")
    JSON.print(restart.io,configs,4)
end


text:
"""
Simulate various Dgps for dirBin0Rec0 various T N and filter
"""

#region import modules and init workers

using Pkg
Pkg.activate(".")
Pkg.instantiate()
using DrWatson
using JLD2
using Distributed
using SharedArrays
using ScoreDrivenERGM
using Logging
import ScoreDrivenERGM:StaticNets, DynNets
import ScoreDrivenERGM.DynNets:SdErgm,SdErgmDirBin0Rec0, simulate_and_estimate_parallel
using ScoreDrivenERGM.Utilities

begin
nWorkers = 10
addprocs(nWorkers - nprocs())
@sync @everywhere begin
    using Pkg
    Pkg.activate(".")
    Pkg.instantiate()
    using ScoreDrivenERGM
    import ScoreDrivenERGM:StaticNets, DynNets
    import ScoreDrivenERGM.DynNets:SdErgm,SdErgmDirBin0Rec0, simulate_and_estimate_parallel
    using ScoreDrivenERGM.Utilities

    model_mle = DynNets.SdErgmDirBin0Rec0_mle(scoreScalingType="FISH_D")

    indTvPar = trues(2)
    end
end
#endregion

dgpSetARlowlow, dgpSetARlow, dgpSetARmed, dgpSetARhigh, dgpSetSIN, dgpSetSDlow, dgpSetSD, dgpSetSDhigh = ScoreDrivenERGM.DynNets.list_example_dgp_settings(model_mle)

# define dictionary with all the different values for each setting of the simulation. The product of all settings will be executed using DrWatson.jl functionalities
c= Dict{String, Any}()
c["model"] =[DynNets.SdErgmDirBin0Rec0_mle(scoreScalingType="FISH_D"), DynNets.SdErgmDirBin0Rec0_pmle(scoreScalingType="FISH_D")]
c["model"][1].options["integrated"] = true
c["model"][2].options["integrated"] = true

c["T"] = [300, 3000]
c["N"] = [100]

dgpSetARlowlow.opt.B[1] = 1
c["dgpSettings"] = [dgpSetARlowlow]
c["nSample"] = 50

list = sort(dict_list(c), by=x->(string(x["model"])))

for d in list

    timeSim = @elapsed allObsT, allvEstSdResPar, allfVecT_filt, allParDgpT, allConvFlag, allftot_0, allfVecT_filt_SS =  ScoreDrivenERGM.DynNets.simulate_and_estimate_parallel(d["model"], d["dgpSettings"], d["T"], d["N"],  d["nSample"];)

    modelTag = string(d["model"])

    res1 =  (;modelTag, allObsT, allvEstSdResPar, allfVecT_filt, allParDgpT, allConvFlag, allftot_0) |>DrWatson.ntuple2dict |> DrWatson.tostringdict

    estDict = merge(res1, d)

    saveName = replace.( savename(d, "jld2";allowedtypes = (Real, String, Symbol, NamedTuple, Tuple, ScoreDrivenERGM.DynNets.SdErgm) ), r"[\"]" => "")

    timeSave = @elapsed save( datadir("sims", "dgp&Fil_est", saveName), estDict)

    Logging.@info("Time sim = $timeSim ,  time save = $timeSave ")

end

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
#####
##### Interface to Pascal Pons' Connect 4 Solver
##### https://github.com/PascalPons/connect4
#####

# Problem: no Connect4 module. We can change this:

module Solver

import ..Game, ..history, ..WHITE, ..NUM_CELLS
import AlphaZero: GI, GameInterface, Benchmark, AbstractPlayer, think

const DEFAULT_SOLVER_DIR = joinpath(@__DIR__, "solver", "connect4")

struct Player <: AbstractPlayer{Game}
  process :: Base.Process
  function Player(;
      solver_dir=DEFAULT_SOLVER_DIR,
      solver_name="c4solver",
      disable_stderr=true)
    cmd = Cmd(`./$solver_name`, dir=solver_dir)
    if disable_stderr
      cmd = pipeline(cmd, stderr=devnull)
    end
    p = open(cmd, "r+")
    return new(p)
  end
end

# Solver protocol
# - Standard input: one position per line
# - Standard output: space separated
#     position, score, number of explored node, computation time in μs

struct SolverOutput
  score :: Int
  num_explored_nodes :: Int64
  time :: Int64 # in μs
end

history_string(game) = reduce(*, map(string, history(game)))

function query_solver(p::Player, g)
  hstr = history_string(g)
  println(p.process, hstr)
  l = readline(p.process)
  args = map(split(l)[2:end]) do x
    parse(Int64, x)
  end
  return SolverOutput(args...)
end

function remaining_stones(game, player)
  @assert !isnothing(game.history)
  n = length(game.history)
  p = n ÷ 2
  (n % 2 == 1 && player == WHITE) && (p += 1)
  return NUM_CELLS ÷ 2 - p
end

function value(player, game)
  z = GI.white_reward(game)
  if isnothing(z)
    return query_solver(player, game).score
  elseif iszero(z)
    return 0
  else
    v = remaining_stones(game, game.winner) + 1
    z < 0 && (v = -v)
    GI.white_playing(game) || (v = -v)
    return v
  end
end

function qvalue(player, game, action)
  @assert isnothing(GI.white_reward(game))
  wp = GI.white_playing(game)
  game = copy(game)
  GI.play!(game, action)
  pswitch = wp != GI.white_playing(game)
  nextv = value(player, game)
  return pswitch ? - nextv : nextv
end

function think(p::Player, g, turn=nothing)
  as = GI.available_actions(g)
  qs = [qvalue(p, g, a) for a in as]
  maxq = maximum(qs)
  opt = findall(>=(maxq), qs)
  π = zeros(length(as))
  π[opt] .= 1 / length(opt)
  return as, π
end

Benchmark.PerfectPlayer(::Type{Game}) = Player

end


text:
function manialn(X::AbstractMatrix{T},Y::AbstractMatrix{T}) where T<:Real
    μ,dim=0.9,30
    n1,n2=size(X,1),size(Y,1)
    W₁,W₂=X.+1,Y.+1
    ℐ=Matrix(I,n1,n2)
    μ = μ*(sum(W₁)+sum(W₂))/(2*sum(ℐ))
    𝕎 = [W₁ μ*ℐ; μ*ℐ' W₂]
    L=diagm(vec(sum(abs.(𝕎),dims=1))).-𝕎
    # λ,V =KrylovKit.eigsolve(L,35,:SR,krylovdim=40)
    # V=hcat(V)
    λ,V = eigen(L)
    i=real(λ).>=1e-8
    V=real(V[:,i])
    dim=min(dim,size(V,2))
    V=V[:,1:dim]
    aln0=V[1:n1,:]
    aln1=V[n1+1:end,:]
    d = norm.((aln0.-aln1)[i,:] for i = 1:n1)
    # _,idx=findmax(dd)
    return d
end

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
abstract type AbstractPotential end

struct LennardJones <: AbstractPotential
    C12::Float64
    C6::Float64
end

function LennardJones(ϵ::Unitful.Energy, σ::Unitful.Length)
    s = auconvert(σ)
    e = auconvert(ϵ)
    A = 4*e*s^12
    B = 4*e*s^6
    return LennardJones(austrip(A),austrip(B))
end

function LennardJones(ϵ::Unitful.Temperature, σ)
    return LennardJones(ϵ*u"k", σ)
end

#(p::LennardJones)(r) = p.C12*r^-12 - p.C6*r^-6

#(p::LennardJones)(r::AbstractArray) = sum(p.C12.*r.^-12 .- p.C6.*r.^-6)

function (p::LennardJones)(r)
    @avx r6 = r.^-6
    return sum(muladd.(p.C12, r6, -p.C6) .* r6)
end


text:
module llh

using LinearAlgebra  #: cholesky, logdet, diag, inv, triu
using AxUtil
using NNlib: softmax
using ..Misc

export gmm_llh, gmm_llh_invLT

function log_gauss_llh(X, mu, sigma; bypass=false)
    if !bypass
        out = _log_gauss_llh(X, mu, sigma)
    else
        out = - ones(size(X, 2))*Inf
    end
    return out
end

function _log_gauss_llh(X, mu, sigma)
    d = size(X,1)
    invLT = Matrix(inv(cholesky(sigma).L))
    Z = invLT*(X .- mu)
    exponent = -0.5*sum(Z.^2, dims=1)  |> dropdim1
    lognormconst = -d*log(2*pi)/2 -0.5*logdet(sigma)  #.-0.5*(-2*sum(log.(diag(invLT))))
    return exponent .+ lognormconst
end

function _log_gauss_llh_invLT(X, mu, invLT)
    d = size(X,1)
    Z = invLT*(X .- mu)
    exponent = -0.5*sum(Z.^2, dims=1)  |> dropdim1
    lognormconst = -d*log(2*pi)/2 .-0.5*(-2*sum(log.(diag(invLT))))
    return exponent .+ lognormconst
end

function gmm_llh(X, pis, mus, sigmas; thrsh_comp=0.005)
    p, n = size(X)
    k = length(pis)
    inactive_ixs = pis[:] .< thrsh_comp

    P = zeros(k, n)
    for j = 1:k
        P[j,:] = log_gauss_llh(X, mus[j,:], sigmas[:,:,j],
            bypass=inactive_ixs[j]) .+ log(pis[j])
    end
    return AxUtil.Math.logsumexpcols(P)
end

function gmm_llh_invLT(X, pis, mus, invLTs::Array{T,1}; disp=false, thrsh_comp=0.005) where T <: AbstractMatrix
    p, n = size(X)
    k = length(pis)
    inactive_ixs = pis[:] .< thrsh_comp

    P = zeros(k, n)
    for j = 1:k
        if !inactive_ixs[j]
            P[j,:] = _log_gauss_llh_invLT(X, mus[j,:], invLTs[j]) .+ log(pis[j])
        else
            P[j,:] .= -Inf
        end
    end
    return AxUtil.Math.logsumexpcols(P)
end

function responsibilities(X, mus::AbstractArray{T,2}, sigmas::AbstractArray{T, 3}, pis::AbstractArray{T, 1}) where T <: AbstractFloat
    softmax(reduce(vcat, [log_gauss_llh(X, mus[j,:], sigmas[:,:,j]) .+ log(pis[j]) for j in 1:size(mus, 1)]'))
end

end


text:
using Revise
using gtpmip

sol = pure_symbolic_solve("test/pddl/domain.pddl", "test/pddl/hard_problem.pddl"; max_levels=1000)


text:
__precompile__()

module StaticArrays

import Base: @pure, @propagate_inbounds, getindex, setindex!, size, similar,
             length, convert, promote_op, map, map!, reduce, mapreduce,
             broadcast, broadcast!, conj, transpose, ctranspose, hcat, vcat,
             ones, zeros, eye, cross, vecdot, reshape, fill, fill!, det, inv,
             eig, trace, vecnorm, dot

export StaticArray, StaticVector, StaticMatrix
export SArray, SVector, SMatrix
export MArray, MVector, MMatrix
export FieldVector, MutableFieldVector

export @SVector, @SMatrix, @SArray
export @MVector, @MMatrix, @MArray

export similar_type

include("util.jl")

include("core.jl")
include("SVector.jl")
include("FieldVector.jl")
include("SMatrix.jl")
include("SArray.jl")
include("MVector.jl")
include("MMatrix.jl")
include("MArray.jl")

include("indexing.jl")
include("abstractarray.jl")
include("mapreduce.jl")
include("arraymath.jl")
include("linalg.jl")
include("matrix_multiply.jl")
include("deque.jl")
include("det.jl")
include("inv.jl")
include("eigen.jl")
include("cholesky.jl")

include("FixedSizeArrays.jl")

# TODO list
# ---------
#
# * tests
#
# * reshape() - accept Val? Currently uses `ReshapedArray`. Cool :)
#
# * permutedims() - accept Val? Or wait for `PermutedDimsArray` ?
#
# * Linear algebra - matrix functions (det, inv, eig, svd, qr, etc...)
#                    (currently, we use pointers to interact with LAPACK, etc)

end # module


text:
sprites = [
  [
    0b11110000,  # ****
    0b10010000,  # *  *
    0b10010000,  # *  *
    0b10010000,  # *  *
    0b11110000   # ****
  ]
  [
    0b00100000,  #  *
    0b01100000,  # **
    0b00100000,  #  *
    0b00100000,  #  *
    0b01110000   # ***
  ]
  [
    0b11110000,  # ****
    0b00010000,  #    *
    0b11110000,  # ****
    0b10000000,  # *
    0b11110000   # ****
  ]
  [
    0b11110000,  # ****
    0b00010000,  #    *
    0b11110000,  # ****
    0b00010000,  #    *
    0b11110000   # ****
  ]
  [
    0b10010000,  # *  *
    0b10010000,  # *  *
    0b11110000,  # ****
    0b00010000,  #    *
    0b00010000   #    *
  ]
  [
    0b11110000,  # ****
    0b10000000,  # *
    0b11110000,  # ****
    0b00010000,  #    *
    0b11110000   # ****
  ]
  [
    0b11110000,  # ****
    0b10000000,  # *
    0b11110000,  # ****
    0b10010000,  # *  *
    0b11110000   # ****
  ]
  [
    0b11110000,  # ****
    0b00010000,  #    *
    0b00100000,  #   *
    0b01000000,  #  *
    0b01000000   #  *
  ]
  [
    0b11110000,  # ****
    0b10010000,  # *  *
    0b11110000,  # ****
    0b10010000,  # *  *
    0b11110000   # ****
  ]
  [
    0b11110000,  # ****
    0b10010000,  # *  *
    0b11110000,  # ****
    0b00010000,  #    *
    0b11110000   # ****
  ]
  [
    0b11110000,  # ****
    0b10010000,  # *  *
    0b11110000,  # ****
    0b10010000,  # *  *
    0b10010000   # *  *
  ]
  [
    0b11100000,  # ***
    0b10010000,  # *  *
    0b11100000,  # ***
    0b10010000,  # *  *
    0b11100000   # ***
  ]
  [
    0b11110000,  # ****
    0b10000000,  # *
    0b10000000,  # *
    0b10000000,  # *
    0b11110000   # ****
  ]
  [
    0b11100000,  # ***
    0b10010000,  # *  *
    0b10010000,  # *  *
    0b10010000,  # *  *
    0b11100000   # ***
  ]
  [
    0b11110000,  # ****
    0b10000000,  # *
    0b11110000,  # ****
    0b10000000,  # *
    0b11110000   # ****
  ]
  [
    0b11110000,  # ****
    0b10000000,  # *
    0b11110000,  # ****
    0b10000000,  # *
    0b10000000   # *
  ]
]


text:
using KubeMeta, Test

@testset "differentiable controller" begin
    chan = Channel(32)
    cli = client()
    lw = TaskListWatcher(cli, chan)
    controller = DifferentiableController()
    informer = TaskInformer(lw, chan, controller)
    KubeMeta.run(informer)
end

max_stars_repo_name: pbarker/kubemeta

max_issues_repo_name: pbarker/kubemeta

max_forks_repo_name: pbarker/kubemeta


text:
const CITY_ALIASES = Dict(
    "LA" => "America/Los_Angeles",
    "NY" => "America/New_York",
    "Rio" => "America/Fortaleza",
    "Brussels" => "Europe/Brussels",
)

const DATE_TIME_PATTERNS = (
    dateformat"yyyy-mm-dd HH:MM",
    dateformat"yyyy-mm-dd HH:MMp", # AM/PM
)

function commander(c::Client, m::Message, ::Val{:time_zone})
    # @info "time_zone_commander called"
    # @info "Message content" m.content m.author.username m.author.discriminator
    startswith(m.content, COMMAND_PREFIX * "tz") || return nothing
    regex = Regex(COMMAND_PREFIX * raw"tz (help|convert|aliases) *(.*)$")
    matches = match(regex, m.content)
    if matches === nothing || matches.captures[1] == "help"
        help_commander(c, m, :time_zone)
        return nothing
    elseif matches.captures[1] == "aliases"
        help_time_zone_aliases(c, m)
    elseif matches.captures[1] == "convert"
        handle_time_zone_conversion(c, m, matches.captures[2])
    else
        @error "Bug: unreachable reached! :-)"
    end
end

function help_commander(c::Client, m::Message, ::Val{:time_zone})
    # @info "Sending help for message" m.id m.author
    return reply(
        c,
        m,
        """
        How to use the `tz` command:
        ```
        tz help
        tz convert <date> <time> <alias/timezone>
        tz aliases
        ```
        Some examples are listed below.
        You can use an alias or any official time zone names from https://en.wikipedia.org/wiki/List_of_tz_database_time_zones:
        ```
        2021-03-21 15:00 NY
        2021-03-21 01:00AM Rio
        2021-03-21 09:00PM America/Los_Angeles
        ```
        """,
    )
end

function help_time_zone_aliases(c::Client, m::Message)
    # @info "Sending aliases for message" m.id m.author
    aliases = join(sort(["$(rpad(k, 12)) => $v" for (k, v) in CITY_ALIASES]), "\n")
    return reply(
        c,
        m,
        """
        tz command's time zone aliases are:
        ```
        $aliases
        ```
        """,
    )
end

function handle_time_zone_conversion(c::Client, m::Message, dts::AbstractString)
    try
        zdt = parse_date_time_spec(dts)
        converted_times = convert_time(zdt)
        output = join([ct.display for ct in converted_times], "\n")
        reply(
            c,
            m,
            """
            Alright, the local time are listed as follows:
            ```
            $output
            ```
            """,
        )
    catch ex
        @show ex
        reply(c, m, "Sorry, I don't understand this format: `$dts`")
    end
end

# Examples:
# 2021-03-21 16:00 America/Los_Angeles
# 2021-03-21 9:00PM NY
function parse_date_time_spec(dts::AbstractString)
    tokens = split_date_time_tokens(dts)
    vanilla_dt = parse_date_time(tokens.datetime)
    timezone = parse_time_zone(tokens.city)
    return ZonedDateTime(vanilla_dt, timezone)
end

function split_date_time_tokens(dts::AbstractString)
    tokens = split(dts, " ")
    length(tokens) === 3 || error("Bad date time spec: $dts")
    return (datetime=join(tokens[1:2], " "), city=tokens[3])
end

function parse_date_time(dt::AbstractString)
    for pattern in DATE_TIME_PATTERNS
        local datetime
        try
            return DateTime(dt, pattern)
        catch ex
        end
    end
    return error("Unable to parse date/time: $dt")
end

function parse_time_zone(city::AbstractString)
    city_time_zone = get(CITY_ALIASES, city, city)
    return TimeZone(city_time_zone)
end

function convert_time(t::ZonedDateTime)
    display_format = "yyyy-mm-dd II:MMp"
    return [
        let
            datetime = astimezone(t, TimeZone(city))
            display = Dates.format(datetime, display_format) * " $alias"
            (; alias, datetime, display)
        end for (alias, city) in CITY_ALIASES
    ]
end


text:
using Documenter, GMT

makedocs(
    modules = [GMT],
    format = :html,
    sitename = "GMT",
    pages = [
        "Home" => "index.md",
        "Some examples" => "examples.md",
        "Draw rectangles examples" => "rectangles.md",
        "Draw frames examples" => "frames.md",
        "Manual" => [
            "usage.md",
            "monolitic.md",
            "modules.md",
        ],
        "The GMT types" => "types.md"
    ],
    html_prettyurls = true,
)

deploydocs(
    repo   = "github.com/GenericMappingTools/GMT.jl.git",
    target = "build",
    julia = "0.7",
    deps   = nothing,
    make   = nothing
)


text:
mutable struct MSetLPSolver <: MSetSolver
  solver
  m # Int
  n_arms # Int

  # Optimisation model (reused by solve_linear).
  model # ::Model
  x # ::Vector{Variable}

  function MSetLPSolver(solver)
    return new(solver, nothing, nothing, nothing, nothing)
  end
end

function build!(solver::MSetLPSolver, m::Int, n_arms::Int)
  solver.m = m
  solver.n_arms = n_arms
  solver.model = Model(solver.solver)
  solver.x = @variable(solver.model, [1:n_arms], binary=true)

  for i in 1:n_arms
    set_name(solver.x[i], "x_$i")
  end

  @constraint(solver.model, sum(solver.x) <= m)
end

has_lp_formulation(::MSetLPSolver) = true

function get_lp_formulation(solver::MSetLPSolver, rewards::Dict{T, Float64}) where T
  return solver.model,
    sum(rewards[i] * solver.x[i] for i in keys(rewards)),
    Dict{Int, JuMP.VariableRef}(i => solver.x[i] for i in keys(rewards))
end

function solve_linear(solver::MSetLPSolver, rewards::Dict{Int, Float64})
  m, obj, vars = get_lp_formulation(solver, rewards)
  @objective(m, Max, obj)

  # Maybe solve_budgeted_linear has already been called on this model.
  if :MSetLP in keys(m.ext)
    budget_constraint = m.ext[:MSetLP][:budget_constraint]
    set_normalized_rhs(budget_constraint, 0)
  end

  set_silent(m)
  optimize!(m)

  if termination_status(m) != MOI.OPTIMAL
    return Int[]
  end
  return Int[i for i in keys(rewards) if value(vars[i]) > 0.5]
end

function solve_budgeted_linear(solver::MSetLPSolver, rewards::Dict{Int, Float64}, weights::Dict{Int, Int}, budget::Int)
  m, obj, vars = get_lp_formulation(solver, rewards)
  @objective(m, Max, obj)

  # Add the budget constraint (or change the existing constraint).
  if :MSetLP in keys(m.ext)
    budget_constraint = m.ext[:MSetLP][:budget_constraint]
    set_normalized_rhs(budget_constraint, budget)
  else
    budget_constraint = @constraint(m, sum(weights[i] * solver.x[i] for i in keys(rewards)) >= budget)
    m.ext[:MSetLP] = Dict(:budget_constraint => budget_constraint)
  end

  set_silent(m)
  optimize!(m)

  if termination_status(m) != MOI.OPTIMAL
    return Int[]
  end
  return Int[i for i in keys(rewards) if value(vars[i]) > 0.5]
end


text:
"""
    csv2bin(;csvfile::AbstractString="", binfile::AbstractString="")
Convert a CSV file to Julia Binary file.
`csvfile` and `binfile` are specified such as Data.csv and Data.zst, respectively.
"""
function csv2bin(;csvfile::AbstractString="", binfile::AbstractString="")
    N = zeros(UInt32, 1)
    M = zeros(UInt32, 1)
    N[] = UInt32(nrow(csvfile=csvfile))
    M[] = UInt32(ncol(csvfile=csvfile))
    counter = 0
    open(binfile, "w") do file
        stream = ZstdCompressorStream(file)
        write(stream, N)
        write(stream, M)
        open(csvfile , "r") do f
            while !eof(f)
                counter += 1
                print("\r", counter)
                xx = readline(f)
                xx = split(xx, ",") # Imported as Character
                # Assume input data is Integer
                x = zeros(UInt32, length(xx))
                for i=1:length(x)
                    x[i] = floor(UInt32, parse(Float64, xx[i]))
                end
                write(stream, x)
            end
        end
        close(stream)
    end
    print("\n")
end

# no.of row of csv
function nrow(;csvfile::AbstractString="")
    nrow = 0
    open(csvfile, "r") do f
        while !eof(f)
            nrow += 1
            readline(f)
        end
    end
    nrow
end

# no. of columns of csv
function ncol(;csvfile::AbstractString="")
    counter = 0
    open(csvfile, "r") do f
        xx = readline(f)
        xx = split(xx, ",")
        length(xx)
    end
end


text:
# This file contains temporary methods which will likely be added to
# Arblib later on


text:
"""
    trace_game(p1::Player, p2::Player, s::State=State(); disp::Bool=true)

Using players `p1` and `p2`, simulate the game starting from state `s`. If
`disp` is true, the board will be printed at the start of each turn. Returns a
`Vector{State}` that records the history of the game
"""
function trace_game(p1::Player, p2::Player, s::State=State(); disp::Bool=true)
    players = (p1, p2)
    history = State[]
    while !(game_over(s)[1])
        push!(history, s)
        disp && show(s)
        actions_bits = legal_actions_bits(s)
        if actions_bits == 0
            s = do_pass(s)
            continue
        end
        action = select_action(players[s.player], s)
        s = do_action(s, action)
    end
    score(s), history
end

"""
     play_game(p1::Player, p2::Player, s::State=State(); disp::Bool=false)

Using players `p1` and `p2`, simulate the game starting from state `s`. If
`disp` is true, the board will be printed at the start of each turn. Returns
the final state
"""
function play_game(p1::Player, p2::Player, s::State=State(); disp::Bool=false)
    players = (p1, p2)
    while !(game_over(s)[1])
        actions_bits = legal_actions_bits(s)
        if actions_bits == 0
            s = do_pass(s)
            continue
        end
        action = select_action(players[s.player], s)
        s = do_action(s, action)
    end
    s
end


text:
using Test
using ClimateMachine
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Grids: polynomialorder
using ClimateMachine.BalanceLaws
using ClimateMachine.Ocean
using ClimateMachine.Ocean.SplitExplicit01
using ClimateMachine.Ocean.OceanProblems

using CLIMAParameters
using CLIMAParameters.Planet: grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

function config_simple_box(
    name,
    resolution,
    dimensions,
    boundary_conditions;
    dt_slow = 90.0 * 60.0,
    dt_fast = 240.0,
)

    problem = OceanGyre{FT}(
        dimensions...;
        τₒ = 0.1,
        λʳ = 10 // 86400,
        θᴱ = 10,
        BC = boundary_conditions,
    )

    add_fast_substeps = 2
    numImplSteps = 5
    numImplSteps > 0 ? ivdc_dt = dt_slow / FT(numImplSteps) : ivdc_dt = dt_slow
    model_3D = OceanModel{FT}(
        param_set,
        problem;
        cʰ = 1,
        κᶜ = FT(0.1),
        add_fast_substeps = add_fast_substeps,
        numImplSteps = numImplSteps,
        ivdc_dt = ivdc_dt,
    )

    N, Nˣ, Nʸ, Nᶻ = resolution
    resolution = (Nˣ, Nʸ, Nᶻ)

    config = ClimateMachine.OceanSplitExplicitConfiguration(
        name,
        N,
        resolution,
        param_set,
        model_3D;
        solver_type = SplitExplicitSolverType{FT}(dt_slow, dt_fast),
    )

    return config
end

function run_simple_box(driver_config, timespan; refDat = ())

    timestart, timeend = timespan
    solver_config = ClimateMachine.SolverConfiguration(
        timestart,
        timeend,
        driver_config,
        init_on_cpu = true,
        ode_dt = driver_config.solver_type.dt_slow,
    )

    ## Create a callback to report state statistics for main MPIStateArrays
    ## every ntFreq timesteps.
    nt_freq = 1 # floor(Int, 1 // 10 * solver_config.timeend / solver_config.dt)
    cb = ClimateMachine.StateCheck.sccreate(
        [
            (solver_config.Q, "oce Q_3D"),
            (solver_config.dg.state_auxiliary, "oce aux"),
            (solver_config.dg.modeldata.Q_2D, "baro Q_2D"),
            (solver_config.dg.modeldata.dg_2D.state_auxiliary, "baro aux"),
        ],
        nt_freq;
        prec = 12,
    )

    result = ClimateMachine.invoke!(solver_config; user_callbacks = [cb])

    ## Check results against reference if present
    ClimateMachine.StateCheck.scprintref(cb)
    if length(refDat) > 0
        @test ClimateMachine.StateCheck.scdocheck(cb, refDat)
    end
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using JuMP
using SCIP

function reopt_model(time_limit_seconds::Real)

    REopt = Model(with_optimizer(SCIP.Optimizer, display_verblevel=0, limits_time=time_limit_seconds))
	return REopt

end


text:
#=
Illustrates the usage of the BetaBinomial type

@author : Spencer Lyon <spencer.lyon@nyu.edu>
          Victoria Gregory <victoria.gregory@nyu.edu>

@date: 2014-08-04

=#

using QuantEcon
using Plots
pyplot()
using LaTeXStrings

n = 50
a_vals = [0.5, 1, 100]
b_vals = [0.5, 1, 100]

pdfs = []
labels = []
for (a, b) in zip(a_vals, b_vals)
    d = BetaBinomial(n, a, b)
    push!(pdfs, pdf(d))
    ab_label = LaTeXString("\$a=$a\$, \$b=$b\$")
    push!(labels, ab_label)
end

plot(0:n, pdfs, label=reshape(labels,1,length(labels)), markershape=:circle)


text:
import IJulia

# workaround #60:
if IJulia.Compat.Sys.isapple()
    ENV["PATH"] = JULIA_HOME*":"*ENV["PATH"]
end

IJulia.init(ARGS)

# import things that we want visible in IJulia but not in REPL's using IJulia
import IJulia: ans, In, Out, clear_history

pushdisplay(IJulia.InlineDisplay())

ccall(:jl_exit_on_sigint, Void, (Cint,), 0)

# the size of truncated output to show should not depend on the terminal
# where the kernel is launched, since the display is elsewhere
ENV["LINES"] = get(ENV, "LINES", 30)
ENV["COLUMNS"] = get(ENV, "COLUMNS", 80)

println(IJulia.orig_STDOUT[], "Starting kernel event loops.")
IJulia.watch_stdio()

# workaround JuliaLang/julia#4259
delete!(task_local_storage(),:SOURCE_PATH)

# workaround JuliaLang/julia#6765
eval(Base, :(is_interactive = true))

IJulia.waitloop()


text:
using CSV, DataFrames, Distributed, Dates, LinearAlgebra, Jevo, Distributions, DelimitedFiles, LambertW
using SharedArrays, TimerOutputs
# Get date to append to output file
date = Dates.format(Dates.today(), "yyyy_mm_dd")

# Get number of workers
if length(ARGS) == 1
    addprocs(parse(Int64, ARGS[1]))
elseif length(ARGS) > 1
    throw(ArgumentError("Only one command line argument (cores)."))
end

# Import packages needed for all workers
@everywhere  begin
    using Jevo
    using Distributions
    using DelimitedFiles
    using SharedArrays
    using LambertW
    using LinearAlgebra

end

# Parameters
@everywhere begin
    reps = 100
    steps = 2 * 10^7
    N = 100
    nu = 1/N
    ϵ = 2
    n = 4
    emat = ϵ * (ones(4, 4) - Matrix{Float64}(I, 4, 4))
    rho_array = [0, 0.1, 0.5, 1, 2, 5]
    l_array = [10]
    f0 = 200 / 2N
end

rescue = false

@everywhere l0_kappa(kappa, l, ϵ=2, n=4) = 1/2 * lambertw(2 * ϵ^2 * N * l * f0 * (n-1)/n^2 * exp(10)/(1+kappa))
#@everywhere fl(l_opt) = l0_kappa(0, 10)/l_opt^2 * n^2 / (n-1) * 1/ϵ
@everywhere fl(l_opt) = l0_kappa(0, 10)^2 / l_opt^2 * n^2 / (n-1) * 1/2

Gamma_results = SharedArray{Float64, 4}(length(rho_array), length(l_array), 1, reps)
RHO = deepcopy(Gamma_results)
L = deepcopy(Gamma_results)
F0 = deepcopy(Gamma_results)

# Function to run one simulation
@everywhere function run(N, f0, fl, rho, nu, l_0, emat, steps, rescue)
    # Initiate population
    pop = Jevo.mono_pop(N=N, l=l_0)
    Jevo.initiate!(pop, opt=true)
    f = Jevo.fermi_fitness(f0=f0, fl=fl(10)/2N)

    for i in 1:steps
        Jevo.bp_substitution!(pop, emat, f)
        if rand() < rho/N
            Jevo.driver_mutation!(pop)
        end
        if rand() < nu
            Jevo.l_substitution!(pop, emat, f)
        end
        # Recover lost sites
        if length(pop.seqs) < 7
            Jevo.initiate!(pop, opt=true)
        end
        if rescue && (Jevo.get_energy(pop, emat) > 3/4*l_0 - 5/2)
            Jevo.initiate!(pop, opt=true)
        end
    end
    return Jevo.get_energy(pop, emat), length(pop.seqs)
end

# Write Metadata
open(date*"_script1_sub_METADATA.txt", "a") do io
    write(io, "N=$N\n")
    write(io, "f0=$f0\n")
    write(io, "fl=$fl\n")
    write(io, "repetitions=$reps\n")
    write(io, "steps=$steps\n")
    write(io, "rho=$rho_array\n")
    write(io, "nu=$nu\n")
    write(io, "l=$l_array\n")
    write(io, "rescue=$rescue\n")
end

# Run simulation once for Julia
@sync @distributed for j in 1:reps
    for i in 1:length(rho_array)
        for l in 1:length(l_array)
            for r in 1:1
                f = Jevo.fermi_fitness(f0=f0, l=l_array[l])
                Gamma_results[i, l, r, j], L[i, l, r, j] = run(N, f0, fl, rho_array[i], nu, l_array[l], emat, steps, rescue)
                RHO[i, l, r, j] = rho_array[i]
                F0[i, l, r, j] = f0
            end
        end
    end
    println("Run $j done.")
end

df = DataFrame(gamma=[(Gamma_results...)...], l=[(L...)...], rho=[(RHO...)...], f0=[(F0...)...])
CSV.write(date * "_script1_sub.csv", df)


text:
module CitationGraphs

export CitationNode, CitationGraph, loadCitationGraph, saveCitationGraph

# ===========================================================================================
# struct CitationNode
# brief description: The node structure of a citation graph, the data of which are collected
#                    from crawling the website https://academic.microsoft.com , which is the
#                    public service website of MAG(Microsoft Academic Graph)
# fields:
#   id: The MAG(Microsoft Academic Graph) id of the paper at this node.
#       We can access the detail of the paper with the id by navigating the web link:
#       https://academic.microsoft.com/paper/$id
#   year: The year of publication of this paper.
#   title: The title of the paper.
#   labels: The labels from MAG(Microsoft Academic Graph).
#   refs: The references of this paper collected from MAG.
#         Please note that this field could be inaccurate: MAG suffers delay of info.
#         Many new papers with refs don't have refs listed in MAG.
#   cites: The citations to this paper from other papers (also collected from MAG).
#         Please note that this field is as inaccurate as refs due to the same reason.
struct CitationNode
    id::Int
    year::Int
    title::String
    labels::Vector{String}
    refs::Vector{Int}
    cites::Vector{Int}
end

# ===========================================================================================
# struct CitationGraph
# brief description: A data structure of citation graph, the data of which are collected from
#                    the website https://academic.microsoft.com , which is the public service
#                    website of MAG(Microsoft Academic Graph)
# fields:
#   nodes: The node of struct CitationNode stored in a dictionary with key = id
#   toBeAnalyzed: The list of nodes to be analyzed. In ijcai_clustering dataset, these nodes
#                 are those published in IJCAI.
struct CitationGraph
    nodes::Dict{Int, CitationNode}
    toBeAnalyzed::Vector{Int}
end

# ===========================================================================================
# function tidyTitle
# brief description: Make a title more tidy before using it. The procedure include the
#                    steps:
#                    (1) remove the spaces at the head and the tail of the title,
#                    (2) replace "<" with "<",
#                    (3) replace ">" with ">",
#                    (4) replace "&" with "&",
#                    (5) replace """ with "\"",
#                    (6) replace "'" with "'",
#                    (7) replace "&#[number];" with a corresponding unicode
# input:
#   title: The String (text) of a title.
# output:
#   A string of the tidied-up title.
function tidyTitle(title::Union{String,SubString{String}})::String
    # ---------------------------------------------------------------------------------------
    # step 1: remove the spaces at the head and the tail
    result = title
    while startswith(result, " ") || startswith(result, "\t")
        result = result[2:end]
    end
    while endswith(result, " ") || endswith(result, "\t")
        result = result[1:end-1]
    end

    # ---------------------------------------------------------------------------------------
    # step 2: replace "<" with "<"
    result = replace(result, "<"=>"<")

    # ---------------------------------------------------------------------------------------
    # step 3: replace ">" with ">",
    result = replace(result, ">"=>">")

    # ---------------------------------------------------------------------------------------
    # step 4: replace "&" with "&"
    result = replace(result, "&"=>"&")

    # ---------------------------------------------------------------------------------------
    # step 5: replace """ with "\""
    result = replace(result, """=>"\"")

    # ---------------------------------------------------------------------------------------
    # step 6: "'" with "'"
    result = replace(result, "'"=>"'")

    # ---------------------------------------------------------------------------------------
    # step 7: replace "&#[number];" with a corresponding unicode
    positions = findall(r"&#([XxA-Fa-f]|\d)+;",result)
    toReplace = Pair{String,String}[]
    for pos in positions
        s = result[pos[1]+2:pos[end]-1]
        if s[1] == 'X' || s[1] == 'x'
            push!(toReplace, "&#$s;"=>""*Char(parse(Int, "0$s")))
        else
            push!(toReplace, "&#$s;"=>""*Char(parse(Int, s)))
        end
    end
    for pattern in toReplace
        result = replace(result, pattern)
    end

    # ---------------------------------------------------------------------------------------
    # step 8: return the result
    result
end

# ===========================================================================================
# function loadCitationGraph
# brief description: Load data from three files (nodes, edges, labels) of a citation graph.
# input:
#   path: The name of the path where the three files are stored.
#   prefix: The prefix of the names of the three files. For example, the prefix of
#           ijcai-citation-graph-nodes.csv is ijcai.
# output:
#    The citation graph represented by the three files.
function loadCitationGraph(path::String, prefix::String)::CitationGraph
    # ---------------------------------------------------------------------------------------
    # step 1: prepare the data structure of the result
    result = CitationGraph(Dict{Int,CitationNode}(),Int[])

    # ---------------------------------------------------------------------------------------
    # step 2: assemble the file names for nodes, edges, and labels
    fileNameOfNodes = "$path/$prefix-citation-graph-nodes.csv"
    fileNameOfEdges = "$path/$prefix-citation-graph-edges.csv"
    fileNameOfLabels = "$path/$prefix-citation-graph-labels.csv"

    # ---------------------------------------------------------------------------------------
    # step 3: load the nodes
    # (3.1) open the file of Nodes
    fileOfNodes = open(fileNameOfNodes)

    # (3.2) examine the first line to check whether the file format is correct
    firstLine = readline(fileOfNodes)
    columnNames = replace.(split(firstLine, ","), " "=>"")
    @assert columnNames[1] == "#id" && columnNames[2] == "in-$prefix" &&
        columnNames[3] == "year" && columnNames[4] == "title"

    # (3.3) read the rest of lines
    for line in readlines(fileOfNodes)
        columns = split(line, ",")
        id = parse(Int, columns[1])
        toBeAnalyzed = parse(Bool, columns[2])
        year = parse(Int, columns[3])
        title = replace(tidyTitle(columns[4]), "[comma]"=>",")
        result.nodes[id] = CitationNode(id,year,title,String[],Int[],Int[])
        if toBeAnalyzed
            push!(result.toBeAnalyzed, id)
        end
    end

    # (3.4) close the file of Nodes
    close(fileOfNodes)

    # ---------------------------------------------------------------------------------------
    # step 4: load the edges
    # (4.1) open the file of edges
    fileOfEdges = open(fileNameOfEdges)

    # (4.2) examine the first line to check whether the file format is correct
    firstLine = readline(fileOfEdges)
    columnNames = replace.(split(firstLine, ","), " "=>"")
    @assert columnNames[1] == "#id" && columnNames[2] == "ref-id"

    # (4.3) read the rest of lines
    for line in readlines(fileOfEdges)
        columns = split(line, ",")
        id = parse(Int, columns[1])
        refID = parse(Int, columns[2])
        push!(result.nodes[id].refs, refID)
        push!(result.nodes[refID].cites, id)
    end

    # (4.4) close the file of edges
    close(fileOfEdges)

    # ---------------------------------------------------------------------------------------
    # step 5: load the labels
    # (5.1) open the file of labels
    fileOfLabels = open(fileNameOfLabels)

    # (5.2) examine the first line to check whether the file format is correct
    firstLine = readline(fileOfLabels)
    columnNames = replace.(split(firstLine, ","), " "=>"")
    @assert columnNames[1] == "#id" && columnNames[2] == "label"

    # (5.3) read the rest of lines
    for line in readlines(fileOfLabels)
        columns = split(line, ",")
        id = parse(Int, columns[1])
        label = columns[2]
        push!(result.nodes[id].labels, label)
    end

    # (5.4) close the file of labels
    close(fileOfLabels)

    # ---------------------------------------------------------------------------------------
    # step 6: return the result
    result
end

# ===========================================================================================
# function saveCitationGraph
# brief description: Save data to three files (nodes, edges, labels) of a citation graph.
# input:
#   path: The name of the path where the three files are stored.
#   prefix: The prefix of the names of the three files. For example, the prefix of
#           ijcai-citation-graph-nodes.csv is ijcai.
#   citationGraph: The citation graph represented by the three files.
# output:
#   nothing
function saveCitationGraph(path::String, prefix::String, citationGraph::CitationGraph)
    # ---------------------------------------------------------------------------------------
    # step 1: assemble the file names for nodes, edges, and labels
    fileNameOfNodes = "$path/$prefix-citation-graph-nodes.csv"
    fileNameOfEdges = "$path/$prefix-citation-graph-edges.csv"
    fileNameOfLabels = "$path/$prefix-citation-graph-labels.csv"

    # ---------------------------------------------------------------------------------------
    # step 2: save the nodes
    # (2.1) open the file of Nodes for writing
    fileOfNodes = open(fileNameOfNodes, "w")

    # (2.2) print the first line with the file format info
    println(fileOfNodes, "#id, in-$prefix, year, title")

    # (2.3) save data into the rest of lines
    idSetForAnalysis = Set(citationGraph.toBeAnalyzed)
    for (id, node) in citationGraph.nodes
        toBeAnalyzed = id ∈ idSetForAnalysis
        year = node.year
        title = tidyTitle(replace(node.title, ","=>"[comma]"))
        println(fileOfNodes, "$id, $toBeAnalyzed, $year, $title")
    end

    # (2.4) close the file of Nodes
    close(fileOfNodes)

    # ---------------------------------------------------------------------------------------
    # step 3: save the edges
    # (3.1) open the file of edges for writing
    fileOfEdges = open(fileNameOfEdges, "w")

    # (3.2) print the first line with the file format info
    println(fileOfEdges, "#id, ref-id")

    # (3.3) save data into the rest of lines
    edgeSet = Set{Tuple{Int,Int}}()
    for (id, node) in citationGraph.nodes
        for refID in node.refs
            push!(edgeSet, (id, refID))
        end
        for citeID in node.cites
            push!(edgeSet, (citeID, id))
        end
    end
    for (id, refID) in edgeSet
        println(fileOfEdges, "$id, $refID")
    end

    # (3.4) close the file of edges
    close(fileOfEdges)

    # ---------------------------------------------------------------------------------------
    # step 4: save the labels
    # (4.1) open the file of labels for writing
    fileOfLabels = open(fileNameOfLabels, "w")

    # (4.2) print the first line with the file format info
    println(fileOfLabels, "#id, label")

    # (4.3) save data to the rest of lines
    for id in idSetForAnalysis
        node = citationGraph.nodes[id]
        for label in node.labels
            println(fileOfLabels, "$id, $label")
        end
    end

    # (4.4) close the file of labels
    close(fileOfLabels)
end

end # module


text:
# -----------
# massflow

using Distributions
using NamedArrays

export massflow
export lost
export recovered
export entered
export recovery_ratio
export massflow_summary
export issource
export issink

const MassDict = Dict{Tech, NamedArray{Float64}}

issource{T <: AbstractTech}(t::T) = length(t.inputs) == 0
issink{T <: AbstractTech}(t::T) = length(t.outputs) == 0

"""
Compute massflows of a system. Optionally with Monte Carlo (MC) simulation.

Arguments:
- if `MC` is false, the transfer coefficients are the expected values of the Dirichlet distribution
- if `MC` is true, the transfer coefficients are sampled form Dirichlet distribution
- `scale_reliability` a factor to scale the `transC_reliability` of all Techs.
"""
function massflow{T <: Real}(sys::System, M_in::Dict{String, Dict{String, T}};
                             MC::Bool=false, scale_reliability::Real=1.0)

    # --- calculate flows
    flow_mats = Dict{String, AbstractArray}()
    M_out = MassDict()
    Ps = Dict{String, AbstractArray}()

    for substance in SUBSTANCE_NAMES

        # -- derive (random) adjacent matrix
        Pmean, rel_vect = get_adj_mat(sys, substance)

        if MC
            P = sample_P(Pmean, rel_vect)
        else
            P = Pmean
        end
        Ps[substance] = P

        # -- calulate flows
        technames = collect(keys(Pmean.dicts[1]))
        m_inputs = NamedArray(zeros(size(P,1)), (technames, ))

        for name in technames
            if haskey(M_in, name)
                m_inputs[name] = M_in[name][substance]
            end
        end

        flow_mats[substance] = calc_massflows(P, m_inputs)

    end

    # --- rearrange results

    for t in sys.techs

        outprod_names = ["$(p.name)" for p in t.outputs]
        if issink(t)
            cnames = ["recovered", "airloss", "soilloss", "waterloss"]
        else
            cnames = vcat(outprod_names, ["airloss", "soilloss", "waterloss"])
        end
        masses = NamedArray(zeros(NSUBSTANCE, length(cnames))-1,
                            (SUBSTANCE_NAMES, cnames),
                            ("substance", "product"))

        for substance in SUBSTANCE_NAMES

            flow = flow_mats[substance]
            P = Ps[substance]

            for prod in outprod_names
                ouput_con = collect(filter(c -> c[1]==Product(prod) && c[2] == t, sys.connections))[1]

                prod_tmp = Product(replace(prod, "transported", ""))

                if haskey(t.transC[substance], prod_tmp)
                    tc = t.transC[substance][prod_tmp]
                else
                    tc = t.transC[substance][Product("x")]
                end

                masses[substance, prod] = flow[t.name, ouput_con[3].name] *
                    tc / P[t.name, ouput_con[3].name]
            end

            for loss in ["airloss", "soilloss", "waterloss"]
                masses[substance, loss] =flow[t.name, "$(t.name)_$loss"]
            end

            if issink(t)
                masses[substance, "recovered"] = flow[t.name, "$(t.name)_recovered"]
            end

        end
        M_out[t] = masses
    end

    return M_out

end

function get_adj_mat(sys::System, substance::String)

    # -- build matrix
    technames = getfield.(collect(sys.techs), :name)
    lossnames = vcat([["$(n)_airloss", "$(n)_soilloss", "$(n)_waterloss"] for n in technames]...)
    sinkrecoverd = ["$(s.name)_recovered" for s in sys.techs if issink(s)]
    Pnames = vcat(technames, lossnames, sinkrecoverd)
    P = NamedArray(zeros(length(Pnames), length(Pnames)), (Pnames, Pnames), ("from", "to"))

    # fill in transfer coefficients
    for c in sys.connections
        prod, from_tech, to_tech = c
        # transC are the same for transported or non-transported products
        prod_tmp = Product(replace(String(prod.name), "transported", ""))

        if haskey(from_tech.transC[substance], prod_tmp)
            tc = from_tech.transC[substance][prod_tmp]
        else
            tc = from_tech.transC[substance][Product("x")]
        end
        P[from_tech.name, to_tech.name] += tc
    end

    # add all losses
    for t in sys.techs
        P[t.name, t.name * "_airloss"] = t.transC[substance][Product("airloss")]
        P[t.name, t.name * "_soilloss"] = t.transC[substance][Product("soilloss")]
        P[t.name, t.name * "_waterloss"] = t.transC[substance][Product("waterloss")]
    end

    # add 'recovered' for sinks
    for t in filter(issink, sys.techs)
        P[t.name, t.name * "_recovered"] = t.transC[substance][Product("recovered")]
    end

    # compile transC_reliability vector
    rel_vect = NamedArray(zeros(length(Pnames)), (Pnames, ))
    for t in sys.techs
        rel_vect[t.name] = t.transC_reliability[substance]
    end

    return P, rel_vect
end

# Sample a random adjacent matrix. Each row is Dirichlet distributed
function sample_P(P::AbstractArray, transC_reliability::AbstractArray)
    # m = zeros(P)
    P2 = 0.0*P
    for i in 1:size(P,1)
        i_zero = P[i,:] .<= 0.0
        alpha = P[i,.!i_zero] * transC_reliability[i]
        P2[i,.!i_zero] = rand( Dirichlet(alpha.array), 1 )
    end
    P2
end

function calc_massflows(P::AbstractArray, inp::AbstractVector)

    size(P,1) == size(P,2) || error("'P' must be a square matrix!")
    size(P,1) == length(inp) || error("'inp' must have $size(P,1) elements")

    names = P.dicts[1]
    P = P.array
    inp = inp.array

    inp = inp'
    m = inp*P * inv(eye(P) - P) # not an optimal implementation from a numerical point of view...

    flows = [(m+inp)[i]*P[i,j] for i=1:size(P,1), j=1:size(P,1)]

    return NamedArray(flows, (names, names), ("from", "to"))

end

# -----------
# summary functions

function lost(M_out::MassDict)
    sum(m[:,(end-2):end] for (t,m) in M_out)
end

function recovered(M_out::MassDict)
    sum(m[:,1:end-3] for (t,m) in M_out if issink(t))
end

function entered(M_in::Dict, sys::System)
    ent = NamedArray(zeros(NSUBSTANCE, 1), (SUBSTANCE_NAMES, ["entered"]), (:substance, :mass))
    in_dicts = [m for (t,m) in M_in if t in getfield.(collect(sys.techs), :name)]
    for substance in SUBSTANCE_NAMES
        ent[substance, "entered"] = sum(get.(in_dicts, substance, 0))
    end
    return ent
end

function recovery_ratio(M_out::MassDict, M_in::Dict, sys::System)
    mass_in = entered(M_in, sys)
    f = recovered(M_out) ./ mass_in
    f[mass_in .== 0.0] = 1.0  # define: 0/0 := 1.0
    setnames!(f, "ratio", 2, 1)
    return(f)
end

function functional_group_losses(M_out::MassDict, M_in::Dict, sys::System)

    functional_groups = [:U, :S, :C, :T, :D]
    losses = NamedArray(zeros(NSUBSTANCE, 1 + length(functional_groups)*3),
                        (SUBSTANCE_NAMES,
                        vcat("entered", [["$(f)_airloss", "$(f)_soilloss", "$(f)_waterloss"] for f in functional_groups]...)),
                        )

    losses[:,"entered"] = entered(M_in, sys)

    for f in functional_groups
        loss = zeros(4, 3)
        t_fg = filter(t -> t.functional_group == f, sys.techs)
        for t in t_fg
            loss += M_out[t][:,(end-2):end]
        end

        losses[:,["$(f)_airloss", "$(f)_soilloss", "$(f)_waterloss"]] = loss
    end
    losses
end

"""
Calculate summary statistics of a Monte Carlo massflow results
"""

function massflow_summary(sys::System, M_in::Dict; MC::Bool=true, n::Int=100,
                          scale_reliability::Real=1.0)

    # -- convert M_in

    sources = [t for t in sys.techs if issource(t)]
    for ts in sources
        haskey(M_in, ts.name) || error("Input masses are not defined for source '$(ts.name)'!")
    end

    summaries = Dict{String, NamedArray{Float64}}()

    #  -- compute masses
    ns = MC ? n : 1             # make only one run if MC == false
    m_outs = [massflow(sys, M_in, MC=MC, scale_reliability=scale_reliability) for i in 1:ns]

    ## quantiles to calculate
    qq = [0.1, 0.5, 0.9]

    # --  recovery ratio
    tmp = hcat((recovery_ratio(m, M_in, sys) for m in m_outs)...)
    rr = hcat(mean(tmp, 2),
              std(tmp, 2),
              NamedArray(mapslices(x -> quantile(x, qq), tmp.array, 2)))
    setnames!(rr, SUBSTANCE_NAMES, 1) #
    setnames!(rr, ["mean", "sd", ["q_$i" for i in qq]...], 2)
    rr.dimnames = (:substance, :stats)

    summaries["recovery_ratio"] = rr

    # --  recovered
    tmp = hcat((recovered(m) for m in m_outs)...)

    rm = hcat(mean(tmp, 2),
              std(tmp, 2),
              NamedArray(mapslices(x-> quantile(x, qq), tmp.array, 2)))
    setnames!(rm, SUBSTANCE_NAMES, 1)
    setnames!(rm, ["mean", "sd", ["q_$i" for i in qq]...], 2)
    rm.dimnames = (:substance, :stats)

    summaries["recovered"] = rm

    # --  losses
    tmp = cat(3, (lost(m) for m in m_outs)...)

    ll = NamedArray(cat(3,
                        mean(tmp.array, 3),
                        std(tmp.array, 3),
                        mapslices(x-> quantile(x, qq), tmp.array, 3)))

    setnames!(ll, SUBSTANCE_NAMES, 1)
    setnames!(ll, ["air loss", "soil loss", "water loss"],2)
    setnames!(ll, ["mean", "sd", ["q_$i" for i in qq]...], 3)
    ll.dimnames = (:substance, :losses, :stats)

    summaries["lost"] = ll

    # -- losses per functional group
    tmp = cat(3, (functional_group_losses(m, M_in, sys) for m in m_outs)...)

    ll = NamedArray(cat(3,
                        mean(tmp.array, 3),
                        std(tmp.array, 3),
                        mapslices(x-> quantile(x, qq), tmp.array, 3)))

    setnames!(ll, SUBSTANCE_NAMES, 1)
    functional_groups = [:U, :S, :C, :T, :D]
    nn = vcat("entered", [["$(f)_airloss", "$(f)_soilloss", "$(f)_waterloss"] for f in functional_groups]...)
    setnames!(ll, nn, 2)
    setnames!(ll, ["mean", "sd", ["q_$i" for i in qq]...], 3)
    ll.dimnames = (:substance, :losses, :stats)

    summaries["functional_group_losses"] = ll

    # --  entered
    summaries["entered"] = entered(M_in, sys)

    return(summaries)
end


text:
immutable MersenneFunctor{T} <: Func{1}
    mt::MersenneTwister
end
@inline call{T}(rf::MersenneFunctor{T}, i...) = rand(rf.mt, T)
immutable RandFunctor{T} <: Func{1}
    valuerange::T
end
@inline call(rf::RandFunctor, i...) = rand(rf.valuerange)
immutable ConstFunctor{T} <: Base.Func{1}
    args::T
end
call(f::ConstFunctor, i...) = f.args

immutable EyeFunc{T} <: Func{1} end
@inline call{T}(ef::Type{EyeFunc{T}}, i::Int, j::Int) = (i==j ? one(T) : zero(T))

immutable UnitFunctor <: Func{1}
    i::Int
    eltype::DataType
end
call(ef::UnitFunctor, i) = ef.i==i ? one(ef.eltype) : zero(ef.eltype)

immutable ConversionIndexFunctor{T, T2} <: Func{1}
    args1::T
    target::Type{T2}
end
call(f::ConversionIndexFunctor, i...) = f.target(f.args1[i...])

immutable IndexFunctor{T} <: Func{1}
    args1::T
end
call(f::IndexFunctor, i) = f.args1[i]

immutable IndexFunc{T} <: Base.Func{1}
    arg::T
    i::Int
end
Base.call{T}(a::IndexFunc{T}, j) = a.arg[a.i,j]

immutable CRowFunctor{M}
    mat::M
end
call(r::CRowFunctor, i::Int) = crow(r.mat, i)

immutable SetindexFunctor{T <: FixedArray, V, N} <: Func{1}
    target::T
    value::V
    index::NTuple{N, Int}
end

function call(sf::SetindexFunctor, i::Int...)
    sf.index == i && return eltype(sf.target)(sf.value)
    sf.target[i...]
end

immutable RowFunctor{M}
    mat::M
end
call(r::RowFunctor, i::Int) = row(r.mat, i)


text:
ifftc(x::Array{Complex{Float64},2})=fftshift(ifft(ifftshift(x)))*prod(size(x))
fftc(x::Array{Complex{Float64},2})=fftshift(fft(ifftshift(x)))


text:
# For visualization
import Graphs: attributes, vertex_index
# Export attributes, these are enumerates as properties for the variables and factors
# REF: http://www.graphviz.org/documentation/
function attributes(v::GraphsNode, g::T)::AttributeDict where T <:GenericIncidenceList
    AttributeDict(
        "label" => v.dfgNode.label,
        "color" => v.dfgNode isa DFGVariable ? "red" : "blue",
        "shape" => v.dfgNode isa DFGVariable ? "ellipse" : "box",
        "fillcolor" => v.dfgNode isa DFGVariable ? "red" : "blue"
        )
end

# This is insanely important - if we don't provide a valid index, the edges don't work correctly.
vertex_index(v::GraphsNode) = v.index

# Accessors
# getLabelDict(dfg::GraphsDFG) = dfg.labelDict
# getDescription(dfg::GraphsDFG) = dfg.description
# setDescription(dfg::GraphsDFG, description::String) = dfg.description = description
# getAddHistory(dfg::GraphsDFG) = dfg.addHistory
# getSolverParams(dfg::GraphsDFG) = dfg.solverParams
# function setSolverParams(dfg::GraphsDFG, solverParams::T) where T <: AbstractParams
#     dfg.solverParams = solverParams
# end

"""
    $(SIGNATURES)
Gets an empty and unique CloudGraphsDFG derived from an existing DFG.
"""
function _getDuplicatedEmptyDFG(dfg::GraphsDFG)::GraphsDFG
    newDfg = GraphsDFG{typeof(dfg.solverParams)}(;
        userId=dfg.userId, robotId=dfg.robotId, sessionId=dfg.sessionId,
        params=deepcopy(dfg.solverParams))
    newDfg.description ="(Copy of) $(dfg.description)"
    return newDfg
end

function exists(dfg::GraphsDFG, node::N) where N <: DFGNode
    return haskey(dfg.labelDict, node.label)
end
exists(dfg::GraphsDFG, nId::Symbol) = haskey(dfg.labelDict, nId)

function isVariable(dfg::GraphsDFG, sym::Symbol)
    return exists(dfg, sym) && dfg.g.vertices[dfg.labelDict[sym]].dfgNode isa DFGVariable
end

function isFactor(dfg::GraphsDFG, sym::Symbol)
    return exists(dfg, sym) && dfg.g.vertices[dfg.labelDict[sym]].dfgNode isa DFGFactor
end

function addVariable!(dfg::GraphsDFG, variable::DFGVariable)::DFGVariable
    if haskey(dfg.labelDict, variable.label)
        error("Variable '$(variable.label)' already exists in the factor graph")
    end
    dfg.nodeCounter += 1
    variable._dfgNodeParams._internalId = dfg.nodeCounter
    v = GraphsNode(dfg.nodeCounter, variable)
    Graphs.add_vertex!(dfg.g, v)
    push!(dfg.labelDict, variable.label=>variable._dfgNodeParams._internalId)
    # Track insertion
    push!(dfg.addHistory, variable.label)

    return variable
end

function addFactor!(dfg::GraphsDFG, factor::DFGFactor)::DFGFactor

    if haskey(dfg.labelDict, factor.label)
        error("Factor '$(factor.label)' already exists in the factor graph")
    end

    variableLabels = factor._variableOrderSymbols
    for v in variableLabels
        if !(v in keys(dfg.labelDict))
            error("Variable '$(v)' not found in graph when creating Factor '$(factor.label)'")
        end
    end

    dfg.nodeCounter += 1
    factor._dfgNodeParams._internalId = dfg.nodeCounter

    fNode = GraphsNode(dfg.nodeCounter, factor)
    f = Graphs.add_vertex!(dfg.g, fNode)
    # Add index
    push!(dfg.labelDict, factor.label=>factor._dfgNodeParams._internalId)
    # Add the edges...
    for varLbl in variableLabels
        variable = getVariable(dfg, varLbl)
        v = dfg.g.vertices[variable._dfgNodeParams._internalId]
        edge = Graphs.make_edge(dfg.g, v, f)
        Graphs.add_edge!(dfg.g, edge)
    end
    # Track insertion
    # push!(dfg.addHistory, factor.label)

    return factor
end

#moved to abstract
# function addFactor!(dfg::GraphsDFG, variableIds::Vector{Symbol}, factor::DFGFactor)::DFGFactor
#     variables = map(vId -> getVariable(dfg, vId), variableIds)
#     return addFactor!(dfg, variables, factor)
# end

# TODO: Confirm we can remove this.
# function getVariable(dfg::GraphsDFG, variableId::Int64)::DFGVariable
#     @warn "This may be slow, rather use by getVariable(dfg, label)"
#     #TODO: This may be slow (O(n)), can we make it better?
#     if !(variableId in values(dfg.labelDict))
#         error("Variable ID '$(variableId)' does not exist in the factor graph")
#     end
#     return dfg.g.vertices[variableId].dfgNode
# end

function getVariable(dfg::GraphsDFG, label::Union{Symbol, String})::DFGVariable
    if typeof(label) == String
        label = Symbol(label)
    end
    if !haskey(dfg.labelDict, label)
        error("Variable label '$(label)' does not exist in the factor graph")
    end
    node = dfg.g.vertices[dfg.labelDict[label]].dfgNode
    !isa(node, AbstractDFGVariable) && error("Node with label '$(label)' is not a variable")
    return node
end

function getFactor(dfg::GraphsDFG, factorId::Int64)::DFGFactor
    @warn "This may be slow, rather use by getFactor(dfg, label)"
    #TODO: This may be slow (O(n)), can we make it better?
    if !(factorId in values(dfg.labelDict))
        error("Factor ID '$(factorId)' does not exist in the factor graph")
    end
    return dfg.g.vertices[factorId].dfgNode
end

function getFactor(dfg::GraphsDFG, label::Union{Symbol, String})::DFGFactor
    if typeof(label) == String
        label = Symbol(label)
    end
    if !haskey(dfg.labelDict, label)
        error("Factor label '$(label)' does not exist in the factor graph")
    end
    node = dfg.g.vertices[dfg.labelDict[label]].dfgNode
    !isa(node, AbstractDFGFactor) && error("Node with label '$(label)' is not a factor")
    return node
end

function updateVariable!(dfg::GraphsDFG, variable::DFGVariable)::DFGVariable
    if !haskey(dfg.labelDict, variable.label)
        @warn "Variable label '$(variable.label)' does not exist in the factor graph, adding"
        return addVariable!(dfg, variable)
    end
    dfg.g.vertices[dfg.labelDict[variable.label]].dfgNode = variable
    return variable
end

function updateFactor!(dfg::GraphsDFG, factor::DFGFactor)::DFGFactor
    if !haskey(dfg.labelDict, factor.label)
        @warn "Factor label '$(factor.label)' does not exist in the factor graph, adding"
        return addFactor!(dfg, factor._variableOrderSymbols, factor)
    end
    dfg.g.vertices[dfg.labelDict[factor.label]].dfgNode = factor
    return factor
end

function deleteVariable!(dfg::GraphsDFG, label::Symbol)#::Tuple{AbstractDFGVariable, Vector{<:AbstractDFGFactor}}
    if !haskey(dfg.labelDict, label)
        error("Variable label '$(label)' does not exist in the factor graph")
    end

    deleteNeighbors = true # reserved, orphaned factors are not supported at this time
    if deleteNeighbors
        neigfacs = map(l->deleteFactor!(dfg, l), getNeighbors(dfg, label))
    end

    variable = dfg.g.vertices[dfg.labelDict[label]].dfgNode
    delete_vertex!(dfg.g.vertices[dfg.labelDict[label]], dfg.g)
    delete!(dfg.labelDict, label)
    return variable, neigfacs
end

function deleteFactor!(dfg::GraphsDFG, label::Symbol)::DFGFactor
    if !haskey(dfg.labelDict, label)
        error("Factor label '$(label)' does not exist in the factor graph")
    end
    factor = dfg.g.vertices[dfg.labelDict[label]].dfgNode
    delete_vertex!(dfg.g.vertices[dfg.labelDict[label]], dfg.g)
    delete!(dfg.labelDict, label)
    return factor
end

function getVariables(dfg::GraphsDFG,
                      regexFilter::Union{Nothing, Regex}=nothing;
                      tags::Vector{Symbol}=Symbol[],
                      solvable::Int=0)::Vector{DFGVariable}
    #
    variables = map(v -> v.dfgNode, filter(n -> (n.dfgNode isa DFGVariable) && (solvable != 0 ? solvable <= isSolvable(n.dfgNode) : true), Graphs.vertices(dfg.g)))
    # filter on solvable

    # filter on regex
    if regexFilter != nothing
        variables = filter(v -> occursin(regexFilter, String(v.label)), variables)
    end

    # filter on tags
    if length(tags) > 0
        mask = map(v -> length(intersect(v.tags, tags)) > 0, variables )
        return variables[mask]
    end
    return variables
end

function getFactors(dfg::GraphsDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0)::Vector{DFGFactor}
    factors = map(v -> v.dfgNode, filter(n -> (n.dfgNode isa DFGFactor) && (solvable != 0 ? solvable <= isSolvable(n.dfgNode) : true), Graphs.vertices(dfg.g)))

    if regexFilter != nothing
        factors = filter(f -> occursin(regexFilter, String(f.label)), factors)
    end

    if length(tags) > 0
        mask = map(v -> length(intersect(v.tags, tags)) > 0, factors )
        return factors[mask]
    end
    return factors
end

function isConnected(dfg::GraphsDFG)::Bool
    return length(Graphs.connected_components(dfg.g)) == 1
end

function getNeighbors(dfg::GraphsDFG, node::T; solvable::Int=0)::Vector{Symbol}  where T <: DFGNode
    if !haskey(dfg.labelDict, node.label)
        error("Variable/factor with label '$(node.label)' does not exist in the factor graph")
    end
    vert = dfg.g.vertices[dfg.labelDict[node.label]]
    neighbors = in_neighbors(vert, dfg.g) #Don't use out_neighbors! It enforces directiveness even if we don't want it
    # Additional filtering
    neighbors = solvable != 0 ? filter(v -> solvable <= isSolvable(v.dfgNode), neighbors) : neighbors
    # Variable sorting (order is important)
    if node isa DFGFactor
        order = intersect(node._variableOrderSymbols, map(v->v.dfgNode.label, neighbors))
        return order
    end

    return map(n -> n.dfgNode.label, neighbors)
end

function getNeighbors(dfg::GraphsDFG, label::Symbol; solvable::Int=0)::Vector{Symbol}  where T <: DFGNode
    if !haskey(dfg.labelDict, label)
        error("Variable/factor with label '$(label)' does not exist in the factor graph")
    end
    vert = dfg.g.vertices[dfg.labelDict[label]]
    neighbors = in_neighbors(vert, dfg.g) #Don't use out_neighbors! It enforces directiveness even if we don't want it
    # Additional filtering
    neighbors = solvable != 0 ? filter(v -> isSolvable(v.dfgNode) >= solvable, neighbors) : neighbors
    # Variable sorting when using a factor (function order is important)
    if vert.dfgNode isa DFGFactor
        vert.dfgNode._variableOrderSymbols
        order = intersect(vert.dfgNode._variableOrderSymbols, map(v->v.dfgNode.label, neighbors))
        return order
    end

    return map(n -> n.dfgNode.label, neighbors)
end

"""
    $(SIGNATURES)
Produces a dot-format of the graph for visualization.
"""
function toDot(dfg::GraphsDFG)::String
    return Graphs.to_dot(dfg.g)
end

function toDotFile(dfg::GraphsDFG, fileName::String="/tmp/dfg.dot")::Nothing
    open(fileName, "w") do fid
        write(fid, Graphs.to_dot(dfg.g))
    end
    return nothing
end


text:
using JuMP
using GLPK

financing = Model(with_optimizer(GLPK.Optimizer))

@variables(financing, begin
    0 <= u[1:5] <= 100
    0 <= v[1:3]
    0 <= w[1:5]
    m
end)

@objective(financing, Max, m) # Money at the end of June

@constraints(financing, begin
    u[1] + v[1] - w[1] == 150 # January
    u[2] + v[2] - w[2] - 1.01u[1] + 1.003w[1] == 100 # February
    u[3] + v[3] - w[3] - 1.01u[2] + 1.003w[2] == -200 # March
    u[4] - w[4] - 1.02v[1] - 1.01u[3] + 1.003w[3] == 200 # April
    u[5] - w[5] - 1.02v[2] - 1.01u[4] + 1.003w[4] == -50 # May
    -m - 1.02v[3] - 1.01u[5] + 1.003w[5] == -300 # June
end)

optimize!(financing)
@show objective_value(financing);

bid_values = [6 3 12 12 8 16]
bid_items = [[1], [2], [3 4], [1 3], [2 4], [1 3 4]]

auction = Model(with_optimizer(GLPK.Optimizer))
@variable(auction, y[1:6], Bin)
@objective(auction, Max, sum(y' .* bid_values))
for i in 1:6
    @constraint(auction, sum(y[j] for j in 1:6 if i in bid_items[j]) <= 1)
end

optimize!(auction)

@show objective_value(auction);
@show value.(y);

using Statistics # Useful for calculations
using Ipopt      # Ipopt since our objective is quadratic

stock_data = [
93.043 51.826 1.063;
84.585 52.823 0.938;
111.453 56.477 1.000;
99.525 49.805 0.938;
95.819 50.287 1.438;
114.708 51.521 1.700;
111.515 51.531 2.540;
113.211 48.664 2.390;
104.942 55.744 3.120;
99.827 47.916 2.980;
91.607 49.438 1.900;
107.937 51.336 1.750;
115.590 55.081 1.800;
]

# Calculating stock returns

stock_returns = Array{Float64}(undef, 12, 3)

for i in 1:12
    stock_returns[i, :] = (stock_data[i + 1, :] .- stock_data[i, :]) ./ stock_data[i, :]
end

# Calculating the expected value of monthly return

r = [Statistics.mean(stock_returns[:,1]) Statistics.mean(stock_returns[:,2]) Statistics.mean(stock_returns[:,3])]'

# Calculating the covariance matrix Q

Q = Statistics.cov(stock_returns)

# JuMP Model

portfolio = Model(with_optimizer(Ipopt.Optimizer, print_level=0))
@variable(portfolio, x[1:3] >= 0)
@objective(portfolio, Min, x' * Q * x)
@constraint(portfolio, sum(x) <= 1000)
@constraint(portfolio, sum(r .* x) >= 50)

optimize!(portfolio)

@show objective_value(portfolio);
@show value.(x);


text:
module test_aheui_basic

using Test
using Aheui

@test "" == 아희("")
@test "" == 아희("아희")
@test "" == 아희("밯망희")
@test "45" == 아희("발밞따망희")
@test "5" == 아희("반받다망희")
@test "Hello, world!\n" == 아희("""
밤밣따빠밣밟따뿌
빠맣파빨받밤뚜뭏
돋밬탕빠맣붏두붇
볻뫃박발뚷투뭏붖
뫃도뫃희멓뭏뭏붘
뫃봌토범더벌뿌뚜
뽑뽀멓멓더벓뻐뚠
뽀덩벐멓뻐덕더벅""")

# test
공간들 = Dict(:큐=>큐([], [], 순방향), :스택=>스택([], [], 순방향), :통로=>통로([], [], 순방향))
어딘가 = 가자(선택, 공간들, :스택)

가자(집어넣기, 어딘가, 2)
가자(집어넣기, 어딘가, 3)
@test [5] == 가자(덧셈, 어딘가, 0)

가자(집어넣기, 어딘가, 2)
가자(집어넣기, 어딘가, 3)
@test [5,6] == 가자(곱셈, 어딘가, 0)

가자(집어넣기, 어딘가, 2)
가자(집어넣기, 어딘가, 3)
@test [5,6,-1] == 가자(뺄셈, 어딘가, 0)

가자(집어넣기, 어딘가, 6)
가자(집어넣기, 어딘가, 3)
@test [5,6,-1,2.0] == 가자(나눗셈, 어딘가, 0)

가자(집어넣기, 어딘가, 5)
가자(집어넣기, 어딘가, 3)
@test [5,6,-1,2.0,2] == 가자(나머지, 어딘가, 0)

가자(집어넣기, 어딘가, 20)
@test 20 == 가자(뽑기, 어딘가, 0)

가자(집어넣기, 어딘가, 10)
@test [5,6,-1,2.0,2,10,10] == 가자(중복, 어딘가, 0)

가자(집어넣기, 어딘가, 5)
가자(집어넣기, 어딘가, 6)
@test (5,6) == 가자(바꿔치기, 어딘가, 0)

어딘가 = 가자(선택, 공간들, :큐)
가자(집어넣기, 어딘가, 7)

어딘가 = 가자(이동, 어딘가, 공간들, :스택)
어딘가 = 가자(선택, 공간들, :스택)

가자(집어넣기, 어딘가, 0)
어딘가 = 가자(조건, 어딘가)

end # module test_aheui_basic