text:
#==============================================================================

          Solving the HJB with a state constraint using
                the explicit method

               Translated from matlab code on Ben Moll's website:
                              http://www.princeton.edu/~moll/HACTproject.htm

==============================================================================#

using LinearAlgebra, SparseArrays, Plots

γ = 2 # parameter from CRRA utility
r = 0.03 # the interest rate
ρ = 0.05 # the discount rate

z1=.1
z2 =.2
z= [z1 z2]

λ1 = 0.02
λ2 = 0.03

λ = [λ1 λ2]

H = 500 # number of points in the grid space
amin=-0.02
amax=2

a = LinRange(amin,amax,H)
da = (amax-amin)/(H-1)

aa = [a a]
zz = ones(H,1)*z

maxit = 20000 # the maximum number of iterations we allow the finite differencing algorithm
crit = 10^(-6) # our critical value

dVf=zeros(H,2)
dVb = zeros(H,2)
c= zeros(H,2)

# Inital Guess
v0=zeros(H,2)
v0[:,1] = (z[1] .+ r.*a).^(1-γ)/(1-γ)/ρ
v0[:,2] = (z[2] .+ r.*a).^(1-γ)/(1-γ)/ρ
global v=v0

dist=[]

# The finite differeing loop

for n in 1:maxit
    global V=v
    # Forward differencing
    dVf[1:H-1,:]=(V[2:H,:]-V[1:H-1,:])/da
    dVf[H,:] = zeros(1,2) # will never be used
    # Backward differencing
    dVb[2:H,:]=(V[2:H,:]-V[1:H-1,:])/da
    dVb[1,:]= (z .+r.*amin).^(-γ)

    # Find consumption and savings with the forward difference
    cf = dVf.^(-1/γ)
    sf = zz +r.*aa -cf

    # Find consumption and savings with the backward difference
    cb= dVb.^(-1/γ)
    sb = zz +r.*aa -cb

    # Find consumption and savings at steady state
    c0 = zz + r.*aa
    dV0 = c0.^(-γ)

    # Now implement the upwind scheme in order to select the best differencing method
    If = sf .> 0
    Ib = sb .< 0
    I0 = (1 .- If - Ib)

    # Make sure the backward difference is using at the grid max
    Ib[H,:] = ones(1,2)
    If[H,:] = zeros(1,2)

    # State constraint at amin is automatically implemented

    global dV_Upwind = dVf.*If + dVb.*Ib + dV0.*I0

    global c = dV_Upwind.^(-1/γ)
    global V_switch = [V[:,2] V[:,1]]

    V_change = c.^(1-γ)/(1-γ) + dV_Upwind.*(zz + r.*aa -c) + ones(H,1)*λ.*(V_switch-V)-ρ.*V

    Δ = .9*da/(findmax(z2 .+r.*a)[1])
    global v = v  + Δ*V_change

    push!(dist, findmax(abs.(V_change))[1])

    if dist[n] < crit
        println("Value Function converged, Iteration=")
        println(n)
        break
    end
end

# Graphs
V_err = c.^(1-γ)/(1-γ) + dV_Upwind.*(zz +r.*aa -c) + ones(H,1)*λ.*(V_switch-V)-ρ.*V
adot = zz + r.*aa -c

plot(dist[:], xlabel="Iteration", ylabel="\$\\lvert V^{n+1}-V^{n}\\rvert\$")

plot(a,V_err, xlabel="k",ylabel="Error in HJB Equation", xlims=(amin,amax), legend=false)

plot(a,V, xlabel="a", ylabel="\$V_{i}(a)\$", xlims=(amin,amax),legend=false)

plot(a,c,xlabel="a", ylabel="\$c_{i}(a)\$", xlims=(amin,amax), legend=false)

plot(a, adot, xlabel="a", ylabel="\$s_{i}(a)\$", xlims=(amin,amax))

# Approximation at the borrowing constraint
u1=(z1+r*amin)^(-γ)
u2 = c[1,2]^(-γ)
u11 = -γ*(z1+r*amin)^(-γ-1)
ν = sqrt(-2*((ρ-r)*u1 + λ1*(u1-u2))/u11)
s_approx = -ν*(a.-amin).^(1/2)

plot(a, adot[:,1], xlabel="a", ylabel="\$s_{i}(a)\$", xlims=(amin,amax), label="\$s_{1}(a)\$")
plot!(a, adot[:,2], label="\$s_{2}(a)\$")
plot!(a, zeros(H,1), line=:dash, label="", color=:black)
plot!(a, s_approx, label="Approximatation")


text:
# # Strings

## load packages
using PlotPlants
#------------------------------------------------------------------------------

# ## Latex strings
# `PlotPlants` uses `PyPlot` to plot figures and render math symbols. To make
#     it more convenient to generate symbols in Latex format,
#     [`latex_symbol`](@ref) and [`latex_unit`](@ref) are provided. Also, to
#     accomendate different text renders (normal render and Latex render),
#     different options are used. For example, when normal render is used,
#     `mathrm` option is used to make strings non-italic; when Latex render is
#     used, an extra `text` option (besides `mathrm`) can be used. The
#     advantage of `text` option is that the text can be bold and italic at the
#     same time.

# ### Latex symbol
for str in [latex_symbol("G"),
            latex_symbol("LAI"),
            latex_symbol("G"; sub="sub", sup="sup"),
            latex_symbol("G"; presub="sub", presup="sup"),
            latex_symbol("G"; option="text"),
            latex_symbol("LAI"; option="text"),
            latex_symbol("G"; sub="sub", sup="sup", option="text"),
            latex_symbol("G"; presub="sub", presup="sup", option="text")]
    @show str;
end
#------------------------------------------------------------------------------

# ### pre-defined symbols (unicode and Latex)
for str in [parse_symbol("alpha"),
            parse_symbol("Omega"),
            parse_symbol("alpha", true),
            parse_symbol("Omega", true),
            parse_symbol("degc"),
            parse_symbol("degc", true)]
    @show str;
end
#------------------------------------------------------------------------------

# ### pre-defined units (unicode and Latex)
for str in [latex_unit("A"),
            latex_unit("undefined"),
            latex_unit("A"; paren=false),
            latex_unit("A"; squared=true),
            latex_unit("A"; space=true),
            latex_unit("A"; unicode=true)]
    @show str;
end
#------------------------------------------------------------------------------

# ## Text
# ### Text font
# Currently, `PlotPlants` supports sans (san-serif) and serif fonts.
use_sans();
use_serif();
#------------------------------------------------------------------------------

# ### Text render
# Currently, `PlotPlants` supports Latex and normal text renders.
use_sans_tex();
use_sans_text();
use_serif_tex();
use_serif_text();
#------------------------------------------------------------------------------

ext: jl


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã ``2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Computes the min depths by rule and action for a grammar
"""
module MinDepths

export MinDepthByRule, MinDepthByAction
export min_depth_rule, min_depth_actions

using ExprSearch
using GrammaticalEvolution

typealias MinDepthByRule Dict{Symbol,Int64}
typealias MinDepthByAction Dict{Symbol,Vector{Int64}}

"""
Compute minimum depth for each rule
"""
function min_depth_rule(grammar::Grammar)
    d = MinDepthByRule()
    changed = Dict{Symbol,Bool}()
    for (k,v) in grammar.rules
        d[k] = typemax(Int64)/2
        changed[k] = true
    end
    while any(values(changed))
        for (k,rule) in grammar.rules
           d_k = min_depth_rule(d, rule)
           changed[k] = d[k] != d_k
           d[k] = d_k
       end
    end
    d
end

min_depth_rule(d::MinDepthByRule, rule::ReferencedRule) = d[rule.symbol]
min_depth_rule(d::MinDepthByRule, rule::Union{RangeRule,Symbol}) = 0 #terminals
min_depth_rule(d::MinDepthByRule, x::Any) = 0 #terminals such as constants?
function min_depth_rule(d::MinDepthByRule, rule::OrRule)
    1 + minimum(map(r->min_depth_rule(d,r), rule.values))
end
function min_depth_rule(d::MinDepthByRule, rule::ExprRule)
    a = filter(r->isa(r,ReferencedRule), rule.args)
    1 + maximum(map(r->min_depth_rule(d,r), a))
end
function min_depth_rule(d::MinDepthByRule, rule::AndRule)
    a = filter(r->isa(r,ReferencedRule), rule.values)
    1 + maximum(map(r->min_depth_rule(d,r), a))
end

"""
Compute minimum depth per action of decision rule
"""
function min_depth_actions(grammar::Grammar)
    d = min_depth_rule(grammar)
    da = min_depth_actions(d, grammar)
    da
end
function min_depth_actions(d::MinDepthByRule, grammar::Grammar)
    da = MinDepthByAction()
    for (k,rule) in grammar.rules
        da[k] = min_depth_actions(d, rule)
    end
    da
end
min_depth_actions(d::MinDepthByRule, rule::ReferencedRule) = Int64[d[rule.symbol]]
min_depth_actions(d::MinDepthByRule, rule::Union{Symbol,Terminal}) = zeros(Int64, 1)
min_depth_actions(d::MinDepthByRule, rule::RangeRule) = zeros(Int64, length(rule.range))
function min_depth_actions(d::MinDepthByRule, rule::OrRule)
    1 + Int64[min_depth_rule(d, v) for v in rule.values]
end
function min_depth_actions(d::MinDepthByRule, rule::Union{AndRule,ExprRule})
    Int64[min_depth_rule(d, rule)]
end

end #module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
import Observables: observe

is_independant(x) = true
is_widget(x) = false
function value_range end
function update_value!(x, value) end

function extract_widgets(dom_root)
    result = []
    walk_dom(dom_root) do x
        if is_widget(x)
            push!(result, x)
        end
    end
    return result
end

struct Dependant <: WidgetsBase.AbstractWidget{Any}
    value
end
# Implement interface for dependant
jsrender(session::Session, x::Dependant) = jsrender(session, x.value)
is_independant(x::Dependant) = false
is_widget(::Dependant) = true
value_range(x::Dependant) = value_range(x.value)
update_value!(x::Dependant, value) = update_value!(x.value, value)

# Implement interface for slider!
is_widget(::Slider) = true
value_range(slider::Slider) = slider.range[]
update_value!(slider::Slider, value) = (slider.value[] = value)

function record_values(f, session, widget)
    empty!(session.message_queue)
    try
        f()
        messages = filter(session.message_queue) do msg
            # filter out the event that triggers updating the obs
            # we actually update on js already
            !(msg[:msg_type] == UpdateObservable &&
                msg[:id] == observe(widget).id)
        end
        return Dict(:msg_type => FusedMessage, :payload => messages)
    catch e
        Base.showerror(stderr, e)
    end
end

function record_state_map(session::Session, handler::Function)
    return record_state_map(session, handler(session, nothing))
end

function while_disconnected(f, session::Session)
    connection = session.connection[]
    session.connection[] = nothing
    f()
    session.connection[] = connection
end

function record_states(session::Session, dom::Hyperscript.Node)
    widgets = extract_widgets(dom)
    rendered = jsrender(session, dom)
    # We'll mess with the message_queue to record the statemap
    # So we copy the current message queue and disconnect the session!

    msgs = copy(session.message_queue)
    empty!(session.message_queue)
    independent = filter(is_independant, widgets)
    independent_states = Dict{String, Any}()
    while_disconnected(session) do
        for widget in independent
            state = Dict{Any, Dict{Symbol,Any}}()
            for value in value_range(widget)
                state[value] = record_values(session, widget) do
                    update_value!(widget, value)
                end
            end
            independent_states[observe(widget).id] = state
        end
    end

    append!(session.message_queue, msgs)
    observable_ids = Observables.obsid.(observe.(independent))

    statemap_script = js"""
        const statemap = $(independent_states)
        const observables = $(observable_ids)
        observables.forEach(id => {
            console.log(id)
            JSServe.on_update(id, function (val) {
                // messages to send for this state of that observable
                const messages = statemap[id][val]
                // not all states trigger events
                // so some states won't have any messages recorded
                if (messages){
                    JSServe.process_message(messages)
                }
            })
        })
    """
    evaljs(session, statemap_script)
    return rendered
end

"""
export_standalone(dom_handler, folder::String;
        absolute_urls=false, content_delivery_url="file://" * folder * "/",
    )

Exports the app defined by `dom_handler` with all its assets to `folder`.
Will write the main html out into `folder/index.html`.
Overwrites all existing files!
If this gets served behind a proxy, set `absolute_urls=true` and
set `content_delivery_url` to your proxy url.
If `clear_folder=true` all files in `folder` will get deleted before exporting again!
"""
function export_standalone(dom_handler, folder::String;
        clear_folder=false, write_index_html=true,
        absolute_urls=false, content_delivery_url="file://" * folder * "/",
    )
    if clear_folder
        for file in readdir(folder)
            rm(joinpath(folder, file), force=true, recursive=true)
        end
    end
    serializer = UrlSerializer(false, folder, absolute_urls, content_delivery_url)
    # set id to "", since we dont needed, and like this we get nicer file names
    session = Session(url_serializer=serializer)
    html_dom = Base.invokelatest(dom_handler, session, (target="/",))
    html_str = dom2html(session, html_dom)
    if write_index_html
        open(joinpath(folder, "index.html"), "w") do io
            println(io, html_str)
        end
        return html_str, session
    else
        return html_str, session
    end
end


text:
using PkgBenchmark
using HMMBase

# path = dirname(dirname(pathof(HMMBase)))
path = pwd()
println("Benchmarking $path...")

timestamp = Int64(time() * 1e9)
results = benchmarkpkg(path, resultfile = "benchmark_$timestamp.json")
export_markdown("benchmark_$timestamp.md", results)


text:
# SucVector
# ---------
#
# Almost the same as CompactBitVector, but the layout of data is different.
# Four chunks of 64 bits are stored in a block and blocks are stored in a vector.
# As a result, all data are interleaved and thus the CPU cache may be used
# efficiently.  The idea is taken from:
# https://github.com/herumi/cybozulib/blob/master/include/cybozu/sucvector.hpp
# Note that in v0.3 blocks will not be packed in a vector.

immutable Block
    # large block
    large::UInt32
    # small blocks
    #   the first small block is used for 8-bit extension of the large block
    #   hence, 40 (= 32 + 8) bits are available in total
    smalls::NTuple{4,UInt8}
    # bit chunks (64bits × 4 = 256bits)
    chunks::NTuple{4,UInt64}
end

const bits_per_chunk =  64
const bits_per_block = 256

function Block(chunks::NTuple{4,UInt64}, offset::Int)
    a =     convert(UInt8, count_ones(chunks[1]))
    b = a + convert(UInt8, count_ones(chunks[2]))
    c = b + convert(UInt8, count_ones(chunks[3]))
    Block(offset & ~UInt32(0), (offset >>> 32, a, b, c), chunks)
end

@inline function block_id(i)
    j = Int(i - 1)
    (j >>> 8) + 1, (j & 0b11111111) + 1
end

@inline function chunk_id(i)
    j = Int(i - 1)
    (j >>> 6) + 1, (j & 0b00111111) + 1
end

"""
Static indexable bit vector.

`SucVector` requires 5/4 bits per bit in order to keep bits and the number of bits
(i.e. a `SucVector` object becomes ~25% bigger when converted from a `BitVector` object).
Bits are split into blocks: four chunks of 64 bits are stored in a block and blocks are stored in a vector.
Blocks are contiguous in momery so that memory cache works efficiently.

The number of bits is limited up to 2^40 bits. This is 128GiB so would be enough in most cases.

Let `n` be the length of a `SucVector`, the asymptotic query times are

* getindex: `O(1)`
* rank: `O(1)`
* select: `O(log n)`
"""
type SucVector <: AbstractIndexableBitVector
    blocks::Vector{Block}
    len::Int
end

SucVector() = SucVector(Block[], 0)

Base.copy(bv::SucVector) = SucVector(copy(bv.blocks), bv.len)

# data size of payload in bytes
sizeof(bv::SucVector) = sizeof(bv.blocks)

function convert(::Type{SucVector}, vec::AbstractVector{Bool})
    len = length(vec)
    @assert len ≤ 2^40
    n_blocks = cld(len, bits_per_block)
    blocks = Vector{Block}(n_blocks)
    offset = 0
    for i in 1:n_blocks
        chunks = read_4chunks(vec, (i - 1) * bits_per_block + 1)
        blocks[i] = Block(chunks, offset)
        for j in 1:4
            offset += count_ones(chunks[j])
        end
    end
    return SucVector(blocks, len)
end

function read_chunk(src, from::Int)
    @assert bits_per_chunk == sizeof(UInt64) * 8
    chunk = UInt64(0)
    for k in 0:bits_per_chunk-1
        i = from + k
        chunk >>= 1
        if i ≤ endof(src) && src[i]
            chunk |= UInt64(1) << 63
        end
    end
    return chunk
end

function read_4chunks(src::AbstractVector{Bool}, from::Int)
    a = read_chunk(src, from)
    b = read_chunk(src, from + bits_per_chunk * 1)
    c = read_chunk(src, from + bits_per_chunk * 2)
    d = read_chunk(src, from + bits_per_chunk * 3)
    return a, b, c, d
end

function read_4chunks(src::BitVector, from::Int)
    @assert rem(from, bits_per_block) == 1
    # NOTE: this method depends on the internal data layout of BitVector
    # (but much faster).
    if length(src) >= from + bits_per_block - 1
        i = div(from - 1, bits_per_chunk) + 1
        a = src.chunks[i]
        b = src.chunks[i+1]
        c = src.chunks[i+2]
        d = src.chunks[i+3]
    else
        a = read_chunk(src, from)
        b = read_chunk(src, from + bits_per_chunk * 1)
        c = read_chunk(src, from + bits_per_chunk * 2)
        d = read_chunk(src, from + bits_per_chunk * 3)
    end
    return a, b, c, d
end

length(v::SucVector) = v.len

@inline function getindex(v::SucVector, i::Integer)
    checkbounds(v, i)
    return unsafe_getindex(v, i)
end

@inline function unsafe_getindex(v::SucVector, i::Integer)
    q, r = block_id(i)
    @inbounds block = v.blocks[q]
    q, r = chunk_id(r)
    @inbounds chunk = block.chunks[q]
    return (chunk >> (r - 1)) & 1 == 1
end

@inline function rank1(bv::SucVector, i::Int)
    i = clamp(i, 0, bv.len)
    if i == 0
        return 0
    end
    q, r = block_id(i)
    @inbounds begin
        block = bv.blocks[q]
        # large block
        ret = Int(block.large) + Int(block.smalls[1]) << 32
        # small block
        q, r = chunk_id(r)
        ret += ifelse(q == 1, 0x00, block.smalls[q])
        # remaining bits
        chunk = block.chunks[q]
    end
    ret += count_ones(chunk & rmask(UInt64, r))
    return ret
end

@inline rank1(bv::SucVector, i::Integer) = rank1(bv, Int(i))

# run getindex and rank1 in one shot
function accrank1(bv::SucVector, i)
    q, r = block_id(i)
    @inbounds begin
        block = bv.blocks[q]
        # large block
        rnk1 = Int(block.large) + Int(block.smalls[1]) << 32
        q, r = chunk_id(r)
        # small block
        rnk1 += ifelse(q == 1, 0x00, block.smalls[q])
        chunk = block.chunks[q]
        rnk1 += count_ones(chunk & rmask(UInt64, r))
    end
    return (chunk >> (r - 1)) & 1 == 1, rnk1
end


text:
"""
$(TYPEDEF)

Constrain only the yaw angle of a Pose2, generally used for Gyrocompass, Magnetometer, Dual-GNSS heading type measurements, or any other similar construct.
"""
mutable struct PartialPriorYawPose2{T <: IIF.SamplableBelief} <: IIF.AbstractPrior
    Z::T
    partial::Tuple{Int}
    PartialPriorYawPose2{T}() where T = new{T}()
    PartialPriorYawPose2{T}(x::T) where {T <: IIF.SamplableBelief}  = new{T}(x, (3,))
end
PartialPriorYawPose2(x::T) where {T <: IIF.SamplableBelief} = PartialPriorYawPose2{T}(x)

function getSample(cf::CalcFactor{<:PartialPriorYawPose2})

  Z = cf.factor.Z
  M = getManifold(cf.factor)
  p = identity_element(M)

  Xc = [0,0,rand(Z)]

  X = hat(M, p, Xc)
  points = exp(M, p, X)
  return points
end

getManifold(::PartialPriorYawPose2) = SpecialEuclidean(2)

## Serialization support

"""
$(TYPEDEF)
"""
mutable struct PackedPartialPriorYawPose2 <: IIF.PackedInferenceType
    Z::String
    PackedPartialPriorYawPose2() = new()
    PackedPartialPriorYawPose2(x::T) where {T <: AbstractString}  = new(x)
end

function convert(::Type{PackedPartialPriorYawPose2}, d::PartialPriorYawPose2)
  PackedPartialPriorYawPose2(string(d.Z))
end
function convert(::Type{PartialPriorYawPose2}, d::PackedPartialPriorYawPose2)
  PartialPriorYawPose2(extractdistribution(d.Z))
end


text:
# Download and install RADMC-3D
binary_file = "master.zip"

# Change to the src directory
cd("src")

# Download and extract the RADMC-3D executable
download("https://github.com/dullemond/radmc3d-0.41/archive/$binary_file", binary_file)

# Remove the directory, if a stale version exists from a previous build
if ispath("radmc-3d")
    rm("radmc-3d", recursive=true)
end
# Unzipping the file
run(`unzip $binary_file`)

# make the binary
cd("radmc3d-0.41-master/src")
run(`make`)


text:
const use_return  = haskey(ENV, "USE_RETURN" ) ? parse(Bool, ENV["USE_RETURN"] ) : false
const USE_GPU     = haskey(ENV, "USE_GPU"    ) ? parse(Bool, ENV["USE_GPU"]    ) : false
const do_viz      = haskey(ENV, "DO_VIZ"     ) ? parse(Bool, ENV["DO_VIZ"]     ) : false
const do_save     = haskey(ENV, "DO_SAVE"    ) ? parse(Bool, ENV["DO_SAVE"]    ) : false
const do_save_viz = haskey(ENV, "DO_SAVE_VIZ") ? parse(Bool, ENV["DO_SAVE_VIZ"]) : false
const nx          = haskey(ENV, "NX"         ) ? parse(Int , ENV["NX"]         ) : 64
const ny          = haskey(ENV, "NY"         ) ? parse(Int , ENV["NY"]         ) : 64
const nz          = haskey(ENV, "NZ"         ) ? parse(Int , ENV["NZ"]         ) : 64
###
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 3)
else
    @init_parallel_stencil(Threads, Float64, 3)
end
using ImplicitGlobalGrid, Plots, Printf, LinearAlgebra, MAT
import MPI

norm_g(A) = (sum2_l = sum(A.^2); sqrt(MPI.Allreduce(sum2_l, MPI.SUM, MPI.COMM_WORLD)))

@views inn(A) = A[2:end-1,2:end-1,2:end-1]

macro innH3(ix,iy,iz)       esc(:( H[$ix+1,$iy+1,$iz+1]*H[$ix+1,$iy+1,$iz+1]*H[$ix+1,$iy+1,$iz+1] )) end
macro av_xi_H3(ix,iy,iz)    esc(:( 0.5*(H[$ix,$iy+1,$iz+1]+H[$ix+1,$iy+1,$iz+1]) * 0.5*(H[$ix,$iy+1,$iz+1]+H[$ix+1,$iy+1,$iz+1]) * 0.5*(H[$ix,$iy+1,$iz+1]+H[$ix+1,$iy+1,$iz+1]) )) end
macro av_yi_H3(ix,iy,iz)    esc(:( 0.5*(H[$ix+1,$iy,$iz+1]+H[$ix+1,$iy+1,$iz+1]) * 0.5*(H[$ix+1,$iy,$iz+1]+H[$ix+1,$iy+1,$iz+1]) * 0.5*(H[$ix+1,$iy,$iz+1]+H[$ix+1,$iy+1,$iz+1]) )) end
macro av_zi_H3(ix,iy,iz)    esc(:( 0.5*(H[$ix+1,$iy+1,$iz]+H[$ix+1,$iy+1,$iz+1]) * 0.5*(H[$ix+1,$iy+1,$iz]+H[$ix+1,$iy+1,$iz+1]) * 0.5*(H[$ix+1,$iy+1,$iz]+H[$ix+1,$iy+1,$iz+1]) )) end
macro av_xi_Re(ix,iy,iz)    esc(:( π + sqrt(π*π + max_lxyz2 / @av_xi_H3($ix,$iy,$iz) * _dt) )) end
macro av_yi_Re(ix,iy,iz)    esc(:( π + sqrt(π*π + max_lxyz2 / @av_yi_H3($ix,$iy,$iz) * _dt) )) end
macro av_zi_Re(ix,iy,iz)    esc(:( π + sqrt(π*π + max_lxyz2 / @av_zi_H3($ix,$iy,$iz) * _dt) )) end
macro Re(ix,iy,iz)          esc(:( π + sqrt(π*π + max_lxyz2 / @innH3($ix,$iy,$iz)    * _dt) )) end
macro av_xi_θr_dτ(ix,iy,iz) esc(:( max_lxyz / Vpdτ / @av_xi_Re($ix,$iy,$iz) * Resc )) end
macro av_yi_θr_dτ(ix,iy,iz) esc(:( max_lxyz / Vpdτ / @av_yi_Re($ix,$iy,$iz) * Resc )) end
macro av_zi_θr_dτ(ix,iy,iz) esc(:( max_lxyz / Vpdτ / @av_zi_Re($ix,$iy,$iz) * Resc )) end
macro dτ_ρ(ix,iy,iz)        esc(:( Vpdτ * max_lxyz / @innH3($ix,$iy,$iz) / @Re($ix,$iy,$iz) * Resc )) end

@parallel_indices (ix,iy,iz) function compute_flux!(qHx, qHy, qHz, H, Vpdτ, Resc, _dt, max_lxyz, max_lxyz2, _dx, _dy, _dz)
    if (ix<=size(qHx,1) && iy<=size(qHx,2) && iz<=size(qHx,3)) qHx[ix,iy,iz]  = (qHx[ix,iy,iz] * @av_xi_θr_dτ(ix,iy,iz) - @av_xi_H3(ix,iy,iz) * _dx * (H[ix+1,iy+1,iz+1] - H[ix,iy+1,iz+1]) ) / (1.0 + @av_xi_θr_dτ(ix,iy,iz))  end
    if (ix<=size(qHy,1) && iy<=size(qHy,2) && iz<=size(qHy,3)) qHy[ix,iy,iz]  = (qHy[ix,iy,iz] * @av_yi_θr_dτ(ix,iy,iz) - @av_yi_H3(ix,iy,iz) * _dy * (H[ix+1,iy+1,iz+1] - H[ix+1,iy,iz+1]) ) / (1.0 + @av_yi_θr_dτ(ix,iy,iz))  end
    if (ix<=size(qHz,1) && iy<=size(qHz,2) && iz<=size(qHz,3)) qHz[ix,iy,iz]  = (qHz[ix,iy,iz] * @av_zi_θr_dτ(ix,iy,iz) - @av_zi_H3(ix,iy,iz) * _dz * (H[ix+1,iy+1,iz+1] - H[ix+1,iy+1,iz]) ) / (1.0 + @av_zi_θr_dτ(ix,iy,iz))  end
    return
end

@parallel_indices (ix,iy,iz) function compute_update!(H, Hold, qHx, qHy, qHz, Vpdτ, Resc, _dt, max_lxyz, max_lxyz2, _dx, _dy, _dz, size_innH_1, size_innH_2, size_innH_3)
    if (ix<=size_innH_1 && iy<=size_innH_2 && iz<=size_innH_3)  H[ix+1,iy+1,iz+1] = (H[ix+1,iy+1,iz+1] + @dτ_ρ(ix,iy,iz) * (_dt * Hold[ix+1,iy+1,iz+1] - (_dx * (qHx[ix+1,iy,iz] - qHx[ix,iy,iz]) + _dy * (qHy[ix,iy+1,iz] - qHy[ix,iy,iz])  + _dz * (qHz[ix,iy,iz+1] - qHz[ix,iy,iz])) )) / (1.0 + _dt * @dτ_ρ(ix,iy,iz))  end
    return
end

@parallel_indices (ix,iy,iz) function compute_flux_res!(qHx2, qHy2, qHz2, H, _dx, _dy, _dz)
    if (ix<=size(qHx2,1) && iy<=size(qHx2,2) && iz<=size(qHx2,3))  qHx2[ix,iy,iz] = -@av_xi_H3(ix,iy,iz) * _dx * (H[ix+1,iy+1,iz+1] - H[ix,iy+1,iz+1])  end
    if (ix<=size(qHy2,1) && iy<=size(qHy2,2) && iz<=size(qHy2,3))  qHy2[ix,iy,iz] = -@av_yi_H3(ix,iy,iz) * _dy * (H[ix+1,iy+1,iz+1] - H[ix+1,iy,iz+1])  end
    if (ix<=size(qHz2,1) && iy<=size(qHz2,2) && iz<=size(qHz2,3))  qHz2[ix,iy,iz] = -@av_zi_H3(ix,iy,iz) * _dz * (H[ix+1,iy+1,iz+1] - H[ix+1,iy+1,iz])  end
    return
end

@parallel_indices (ix,iy,iz) function check_res!(ResH, H, Hold, qHx2, qHy2, qHz2, _dt, _dx, _dy, _dz)
    if (ix<=size(ResH,1) && iy<=size(ResH,2) && iz<=size(ResH,3))  ResH[ix,iy,iz] = -_dt * (H[ix+1,iy+1,iz+1] - Hold[ix+1,iy+1,iz+1]) - (_dx * (qHx[ix+1,iy,iz] - qHx[ix,iy,iz]) + _dy * (qHy[ix,iy+1,iz] - qHy[ix,iy,iz])  + _dz * (qHz[ix,iy,iz+1] - qHz[ix,iy,iz]))  end
    return
end

@parallel_indices (ix,iy,iz) function assign!(Hold, H)
    if (ix<=size(H,1) && iy<=size(H,2) && iz<=size(H,3))  Hold[ix,iy,iz] = H[ix,iy,iz]  end
    return
end

@views function diffusion_3D_()
    # Physics
    lx, ly, lz = 10.0, 10.0, 10.0 # domain size
    ttot       = 0.4              # total simulation time
    dt         = 0.2              # physical time step
    # Numerics
    tol        = 1e-8             # tolerance
    itMax      = 1e3              # max number of iterations
    nout       = 2000             # tol check
    CFL        = 1 / sqrt(3)      # CFL number
    Resc       = 1 / 1.2          # iteration parameter scaling
    me, dims, nprocs = init_global_grid(nx, ny, nz) # MPI initialisation
    @static if USE_GPU select_device() end    # select one GPU per MPI local rank (if >1 GPU per node)
    b_width    = (16, 4, 4)       # boundary width for comm/comp overlap
    # Derived numerics
    dx, dy, dz = lx/nx_g(), ly/ny_g(), lz/nz_g() # cell sizes
    Vpdτ       = CFL * min(dx, dy, dz)
    max_lxyz   = max(lx, ly, lz)
    max_lxyz2  = max_lxyz^2
    xc, yc, zc = LinRange(dx/2, lx-dx/2, nx), LinRange(dy/2, ly-dy/2, ny), LinRange(dz/2, lz-dz/2, nz)
    _dx, _dy, _dz, _dt = 1.0/dx, 1.0/dy, 1.0/dz, 1.0/dt
    # Array allocation
    qHx        = @zeros(nx-1,ny-2,nz-2)
    qHy        = @zeros(nx-2,ny-1,nz-2)
    qHz        = @zeros(nx-2,ny-2,nz-1)
    qHx2       = @zeros(nx-1,ny-2,nz-2)
    qHy2       = @zeros(nx-2,ny-1,nz-2)
    qHz2       = @zeros(nx-2,ny-2,nz-1)
    ResH       = @zeros(nx-2,ny-2,nz-2)
    # Initial condition
    H0         = zeros(nx,ny,nz)
    H0         = Data.Array([exp(-(x_g(ix,dx,H0)-0.5*lx+dx/2)*(x_g(ix,dx,H0)-0.5*lx+dx/2) - (y_g(iy,dy,H0)-0.5*ly+dy/2)*(y_g(iy,dy,H0)-0.5*ly+dy/2) - (z_g(iz,dz,H0)-0.5*lz+dz/2)*(z_g(iz,dz,H0)-0.5*lz+dz/2)) for ix=1:size(H0,1), iy=1:size(H0,2), iz=1:size(H0,3)])
    Hold       = @ones(nx,ny,nz) .* H0
    H          = @ones(nx,ny,nz) .* H0
    size_innH_1, size_innH_2, size_innH_3 = size(H,1)-2, size(H,2)-2, size(H,3)-2
    len_ResH_g = ((nx-2-2)*dims[1]+2)*((ny-2-2)*dims[2]+2)*((nz-2-2)*dims[3]+2)
    if do_viz || do_save_viz
        if (me==0) ENV["GKSwstype"]="nul"; if do_viz !ispath("../../figures") && mkdir("../../figures") end; end
        nx_v, ny_v, nz_v = (nx-2)*dims[1], (ny-2)*dims[2], (nz-2)*dims[3]
        if (nx_v*ny_v*nz_v*sizeof(Data.Number) > 0.8*Sys.free_memory()) error("Not enough memory for visualization.") end
        H_v    = zeros(nx_v, ny_v, nz_v) # global array for visu
        H_inn  = zeros(nx-2, ny-2, nz-2) # no halo local array for visu
        z_sl   = Int(ceil(nz_g()/2))     # Central z-slice
        Xi_g, Yi_g = dx+dx/2:dx:lx-dx-dx/2, dy+dy/2:dy:ly-dy-dy/2 # inner points only
    end
    t = 0.0; it = 0; ittot = 0; nt = Int(ceil(ttot/dt)); niter = 0
    # Physical time loop
    while it<nt
        iter = 0; err = 2*tol
        # Pseudo-transient iteration
        while err>tol && iter<itMax
            if (it==1 && iter==0) tic(); niter = 0 end
            @parallel compute_flux!(qHx, qHy, qHz, H, Vpdτ, Resc, _dt, max_lxyz, max_lxyz2, _dx, _dy, _dz)
            @hide_communication b_width begin # communication/computation overlap
                @parallel compute_update!(H, Hold, qHx, qHy, qHz, Vpdτ, Resc, _dt, max_lxyz, max_lxyz2, _dx, _dy, _dz, size_innH_1, size_innH_2, size_innH_3)
                update_halo!(H)
            end
            iter += 1;  niter += 1
            if iter % nout == 0
                @parallel compute_flux_res!(qHx2, qHy2, qHz2, H, _dx, _dy, _dz)
                @parallel check_res!(ResH, H, Hold, qHx2, qHy2, qHz2, _dt, _dx, _dy, _dz)
                err = norm_g(ResH) / sqrt(len_ResH_g)
            end
        end
        ittot += iter; it += 1; t += dt
        @parallel assign!(Hold, H)
        if isnan(err) error("NaN") end
    end
    t_toc = toc()
    A_eff = (2*4+2)/1e9*nx*ny*nz*sizeof(Data.Number) # Effective main memory access per iteration [GB]
    t_it  = t_toc/niter                              # Execution time per iteration [s]
    T_eff = A_eff/t_it                               # Effective memory throughput [GB/s]
    if (me==0) @printf("PERF: Time = %1.3f sec, T_eff = %1.2f GB/s (niter = %d)\n", t_toc, round(T_eff, sigdigits=3), niter) end
    if (me==0) @printf("Total time = %1.2f, time steps = %d, nx = %d, iterations tot = %d \n", round(ttot, sigdigits=2), it, nx_g(), ittot) end
    # Visualise
    if do_viz || do_save_viz
        H_inn .= inn(H); gather!(H_inn, H_v)
        if me==0 && do_viz
            heatmap(Xi_g, Yi_g, H_v[:,:,z_sl]', dpi=150, aspect_ratio=1, framestyle=:box, xlims=(Xi_g[1],Xi_g[end]), ylims=(Yi_g[1],Yi_g[end]), xlabel="lx", ylabel="ly", c=:viridis, clims=(0,1), title="nonlinear diffusion (nt=$it, iters=$ittot)")
            savefig("../../figures/diff_3D_nonlin_multixpu_perf_$(nx_g()).png")
        end
    end
    if me==0 && do_save
        !ispath("../../output") && mkdir("../../output")
        open("../../output/out_diff_3D_nonlin_multixpu_perf.txt","a") do io
            println(io, "$(nprocs) $(nx_g()) $(ny_g()) $(nz_g()) $(ittot) $(t_toc) $(A_eff) $(t_it) $(T_eff)")
        end
    end
    if me==0 && do_save_viz
        !ispath("../../out_visu") && mkdir("../../out_visu")
        matwrite("../../out_visu/diff_3D_nonlin_multixpu_perf.mat", Dict("H_3D"=> Array(H_v), "xc_3D"=> Array(xc), "yc_3D"=> Array(yc), "zc_3D"=> Array(zc)); compress = true)
    end
    finalize_global_grid()
    return xc, yc, zc, H
end

if use_return
    xc, yc, zc, H = diffusion_3D_();
else
    diffusion_3D = begin diffusion_3D_(); return; end
end


text:
@testset "test_nothing_match" begin
    @test @match 1 begin
        nothing => false
        _ => true
    end
end


text:
# Fit a mixture model using EM. Input argument L is the n x m
# conditional likelihood matrix, where n is the number of samples and
# m is the number of mixture components; optional input argument w is
# the initial estimate of the mixture weights.
function mixEM(L; w = ones(size(L,2))/size(L,2), maxiter = 10000,
               tol = 1e-6, eps = 1e-8)

  # Get the number of rows (n) and columns (k) of the conditional
  # likelihood matrix.
  n = size(L,1);
  k = size(L,2);

  # Initialize storage for outputs obj, maxd and timing.
  obj    = zeros(maxiter);
  maxd   = zeros(maxiter);
  timing = zeros(maxiter);

  # Compute the objective function value at the initial iterate.
  obj[1] = -sum(log.(L * w .+ eps));

  # Repeat until convergence criterion is met, or until the maximum
  # number of iterations is reached.
  numiter = 0
  for i = 2:maxiter
    numiter = i;

    # Save the current estimate of the mixture weights.
    out, timing[i] = @timed begin w0 = w;

      # E STEP
      # ------
      # Compute the posterior probabilities.
      P = L * sparse(Diagonal(w));
      P = P ./ repeat(sum(P,dims = 2) .+ eps,1,k);

      # M STEP
      # ------
      # Update the mixture weights.
      w = mean(P,dims = 1)'[:];

      # COMPUTE OBJECTIVE
      # -----------------
      obj[i] = -sum(log.(L * w .+ eps));
    end

    # CHECK CONVERGENCE
    # -----------------
    # Convergence is reached when the maximum difference between the
    # mixture weights at two successive iterations is less than the
    # specified tolerance, or when objective increases.
    maxd[i] = maximum(abs.(w - w0));
    if maxd[i] < tol
      break
    end
  end

  # Return the mixture weights and other optimization info.
  return w, obj[1:numiter], maxd[1:numiter], timing[1:numiter]
end

max_stars_repo_name: stephenslab/mixsqp-paper

max_issues_repo_name: stephenslab/mixsqp-paper

max_forks_repo_name: stephenslab/mixsqp-paper


text:
using Colors, FixedPointNumbers, Base.Test

# test utility function weighted_color_mean
parametric2 = [GrayA,AGray32,AGray]
parametric3 = ColorTypes.parametric3
parametric4A = setdiff(subtypes(ColorAlpha),[GrayA])
parametricA4 = setdiff(subtypes(AlphaColor),[AGray32,AGray])
colortypes = vcat(parametric3,parametric4A,parametricA4)
colorElementTypes = [Float16,Float32,Float64,BigFloat,N0f8,N6f10,N4f12,N2f14,N0f16]

for T in colorElementTypes
    c1 = Gray{T}(1)
    c2 = Gray{T}(0)
    @test weighted_color_mean(0.5,c1,c2) == Gray{T}(0.5)
end

for C in parametric2
    for T in colorElementTypes
        C == AGray32 && (T=N0f8)
	c1 = C(T(0),T(1))
        c2 = C(T(1),T(0))
	@test weighted_color_mean(0.5,c1,c2) == C(T(1)-T(0.5),T(0.5))
    end
end

for C in colortypes
    for T in (Float16,Float32,Float64,BigFloat)
	if issubtype(C,Color)
            c1 = C(T(1),T(1),T(0))
            c2 = C(T(0),T(1),T(1))
	    @test weighted_color_mean(0.5,c1,c2) == C(T(0.5),T(0.5)+T(1)-T(0.5),T(1)-T(0.5))
	else
            C == ARGB32 && (T=N0f8)
            c1 = C(T(1),T(1),T(0),T(1))
            c2 = C(T(0),T(1),T(1),T(0))
	    @test weighted_color_mean(0.5,c1,c2) == C(T(0.5),T(0.5)+T(1)-T(0.5),T(1)-T(0.5),T(0.5))
	end
    end
end

# test utility function linspace
# linspace uses weighted_color_mean which is extensively tested.
# Therefore it suffices to test the function using gray colors.
for T in colorElementTypes
    c1 = Gray(T(1))
    c2 = Gray(T(0))
    linc1c2 = linspace(c1,c2,43)
    @test length(linc1c2) == 43
    @test linc1c2[1] == c1
    @test linc1c2[end] == c2
    @test linc1c2[22] == Gray(T(0.5))
    @test typeof(linc1c2) == Array{Gray{T},1}
end


text:
@testset "molgraph" begin

    @testset "graphmol" begin
        mol = graphmol(SDFileAtom, SDFileBond)
        @test atomcount(mol) == 0
        @test bondcount(mol) == 0
        edges = [(1, 2)]
        atoms = [SDFileAtom(),SDFileAtom(),SDFileAtom()]
        bonds = [SDFileBond()]
        mol = graphmol(edges, atoms, bonds)
        newbond = addbond!(mol, 2, 3, SDFileBond())
        @test newbond == 2
        setatom!(mol, 1, SDFileAtom(:O))
        setbond!(mol, 1, SDFileBond(2))
        @test getatom(mol, 1).symbol == :O
        @test getbond(mol, 1).order == 2
        mol.cache[:test] = [1, 1, 1]
        cp = graphmol(mol)
        @test !haskey(cp.cache, :test)
        cl = clone(mol)
        @test length(cl.cache[:test]) == 3
        @test length(mol.cache[:test]) == 3
    end

    @testset "querymol" begin
        q = querymol(SmartsAtom, SmartsBond)
        addatom!(q, SmartsAtom(:Any => true))
        addatom!(q, SmartsAtom(:Any => true))
        addbond!(q, 1, 2, SmartsBond(:Any => true))
        @test atomcount(q) == 2
        @test bondcount(q) == 1
        @test hasbond(q, 1, 2)
    end

end


text:
function histogram(data::Coord;bins::Int=10,norm::Real=1.0,args...)
    # validation
    bins < 1 && throw(DomainError(bins, "at least one bin is required"))
    norm < 0 && throw(DomainError(norm, "norm must be a positive number."))

    x, y = hist(data,bins)
    y = norm*y/(step(x)*sum(y))  # make area under histogram equal to norm

    bar(x, y; boxwidth="0.9 relative", fillstyle="solid 0.5",args...)
end


text:
using DSGE, ModelConstructors, Dates, OrderedCollections, Test, CSV, DataFrames, Random, FileIO

regenerate_output = false

# Test automatic enforcement of ZLB as a temporary ZLB with no regime-switching
custom_settings = Dict{Symbol, Setting}(
    :data_vintage             => Setting(:data_vintage, "160812"),
    :cond_vintage             => Setting(:cond_vintage, "160812"),
    :cond_id                  => Setting(:cond_id, 0),
    :use_population_forecast  => Setting(:use_population_forecast, true),
    :date_presample_start     => Setting(:date_presample_start, Date(1959, 9, 30)),
    :date_forecast_start      => Setting(:date_forecast_start, DSGE.quartertodate("2016-Q3")),
    :date_conditional_end     => Setting(:date_conditional_end, DSGE.quartertodate("2016-Q3")),
    :n_mon_anticipated_shocks => Setting(:n_mon_anticipated_shocks, 6))

dfs = Dict()
path = dirname(@__FILE__)
dfs[:full] = load("$path/../reference/regime_switch_data.jld2", "full")
dfs[:none] = load("$path/../reference/regime_switch_data.jld2", "none")
dfs[:full][end, :obs_nominalrate] = .5
dfs[:none][end, :obs_nominalrate] = 0.5
dfs[:full][end, :obs_gdp] = -.5

m = Model1002("ss10", custom_settings = custom_settings)
m <= Setting(:rate_expectations_source, :ois)
m.settings[:regime_switching] = Setting(:regime_switching, false)
m.settings[:forecast_ndraws] = Setting(:forecast_ndraws, 3)
m.settings[:forecast_block_size] = Setting(:forecast_block_size, 2)
m.settings[:forecast_jstep] = Setting(:forecast_jstep, 1)
output_vars = [:forecastobs, :bddforecastobs]

para = repeat(Matrix(map(x -> x.value, m.parameters)'), 3, 1)
out = Dict()
Random.seed!(1793 * 1000)
for cond_type in [:full, :none]
    out[cond_type] = Dict()
    forecast_one(m, :mode, cond_type, output_vars, verbose = :none, params = para[1, :], df = dfs[cond_type],
                 zlb_method = :temporary_altpolicy)
    forecast_one(m, :full, cond_type, output_vars, verbose = :none, params = para, df = dfs[cond_type],
                 zlb_method = :temporary_altpolicy)
    for input_type in [:mode, :full]
        out[cond_type][input_type] = Dict()
        outfn = get_forecast_output_files(m, input_type, cond_type, output_vars)
        out[cond_type][input_type][:bddforecastobs] = load(outfn[:bddforecastobs], "arr")
        out[cond_type][input_type][:forecastobs]    = load(outfn[:forecastobs], "arr")
    end
end

@testset "Automatic enforcement of ZLB as temporary alternative policy (no regime-switching)" begin
    for cond_type in [:none, :full]
        for input_type in [:mode, :full]
            if input_type == :mode
                @test all(out[cond_type][input_type][:bddforecastobs][m.observables[:obs_nominalrate], :] .> -1e-14)
            else
                @test all(out[cond_type][input_type][:bddforecastobs][:, m.observables[:obs_nominalrate], :] .> -1e-14)
            end
        end
    end
end

## Regime switching and full-distriution
# Initialize model objects
Random.seed!(1793)
m = Model1002("ss60")
if (VERSION >= v"1.3")
    df_full = DataFrame!(CSV.File(joinpath(dirname(@__FILE__), "../reference/uncertain_altpolicy_data.csv")))
else
    df_full = DataFrame(CSV.read(joinpath(dirname(@__FILE__), "../reference/uncertain_altpolicy_data.csv")))
end
m <= Setting(:forecast_horizons, 30)
m <= Setting(:cond_full_names, [:obs_gdp, :obs_corepce, :obs_spread, # Have to add anticipated rates to conditional data
                                :obs_nominalrate, :obs_longrate,
                                :obs_nominalrate1, :obs_nominalrate2, :obs_nominalrate3,
                                :obs_nominalrate4, :obs_nominalrate5, :obs_nominalrate6])
m <= Setting(:cond_semi_names, [:obs_spread,
                                :obs_nominalrate, :obs_longrate,
                                :obs_nominalrate1, :obs_nominalrate2, :obs_nominalrate3,
                                :obs_nominalrate4, :obs_nominalrate5, :obs_nominalrate6])
m <= Setting(:date_forecast_start, Date(2020, 6, 30))
m <= Setting(:date_conditional_end, Date(2020, 6, 30))
m <= Setting(:forecast_ndraws, 10)
m <= Setting(:forecast_block_size, get_setting(m, :forecast_ndraws))
m <= Setting(:forecast_jstep, 1)
m <= Setting(:use_parallel_workers, false)
setup_regime_switching_inds!(m; cond_type = :full)

# For parameters, use modal ones, but also draw from prior for the COVID-19 shocks
mparas = repeat(ModelConstructors.get_values(m.parameters)', get_setting(m, :forecast_ndraws), 1)
j = length(m.parameters) # To help index for regime-switching parameters
biidc_keys = Dict()
ziid_keys = Dict()
φ_keys = Dict()
for p in m.parameters # Find the keys of the COVID-19 shocks
    if haskey(p.regimes, :value)
        for k in keys(p.regimes[:value])
            if k > 1
                global j += 1
                if p.key == :σ_φ
                    φ_keys[k] = j
                elseif p.key == :σ_biidc
                    biidc_keys[k] = j
                elseif p.key == :σ_ziid
                    ziid_keys[k] = j
                end
            end
        end
    end
end

for i in 1:size(mparas, 1)
    mparas[i, φ_keys[2]] = rand(get(m.parameters[m.keys[:σ_φ]].prior)) # draw from prior for second-regime value
    mparas[i, biidc_keys[2]] = rand(get(m.parameters[m.keys[:σ_biidc]].prior)) # draw from prior for second-regime value
    mparas[i, ziid_keys[2]] = rand(get(m.parameters[m.keys[:σ_ziid]].prior)) # draw from prior for second-regime value
    mparas[i, φ_keys[3]] = rand(get(m.parameters[m.keys[:σ_φ]].prior)) # draw from prior for third-regime value
    mparas[i, biidc_keys[3]] = rand(get(m.parameters[m.keys[:σ_biidc]].prior)) # draw from prior for third-regime value
    mparas[i, ziid_keys[3]] = rand(get(m.parameters[m.keys[:σ_ziid]].prior)) # draw from prior for third-regime value
end

# Function for adding new regimes to parameters
function covid_set_regime_vals(m::AbstractDSGEModel, n::Int)
    if n > 4
        for p in m.parameters
            if haskey(p.regimes, :value)
                for i in 5:n
                    ModelConstructors.set_regime_val!(p, i, ModelConstructors.regime_val(p, 4))
                end
            end
        end
    end
end

# Run forecast!
output_vars = [:forecastobs, :bddforecastobs]
forecast_one(m, :full, :full, output_vars, verbose = :none, params = mparas, df = df_full,
             zlb_method = :temporary_altpolicy, set_regime_vals_altpolicy = covid_set_regime_vals)

# Either test against saved output or re-generate output
output_files = get_forecast_output_files(m, :full, :full, output_vars)
if regenerate_output
    using JLD2, FileIO
    if (VERSION >= v"1.5")
        JLD2.jldopen(joinpath(dirname(@__FILE__), "../reference/automatic_tempalt_zlb_fulldist_v1p5.jld2"), true, true, true, IOStream) do file
            for (k, v) in output_files
                write(file, string(k), load(v, "arr"))
            end
        end
    else
        JLD2.jldopen(joinpath(dirname(@__FILE__), "../reference/automatic_tempalt_zlb_fulldist.jld2"), true, true, true, IOStream) do file
            for (k, v) in output_files
                write(file, string(k), load(v, "arr"))
            end
        end
    end
else
    refdata = (VERSION >= v"1.5") ?
        load(joinpath(dirname(@__FILE__), "../reference/automatic_tempalt_zlb_fulldist_v1p5.jld2")) :
        load(joinpath(dirname(@__FILE__), "../reference/automatic_tempalt_zlb_fulldist.jld2"))
    @testset "Automatic enforcement of ZLB as a temporary alternative policy during full-distribution forecast" begin
        for (k, v) in output_files
            @test @test_matrix_approx_eq refdata[string(k)] load(v, "arr")
            if k == :bddforecastobs
                @test all(refdata[string(k)][:, m.observables[:obs_nominalrate], :] .> -1e-14)
            else
                @test !all(refdata[string(k)][:, m.observables[:obs_nominalrate], :] .> -1e-14)
            end
        end
    end
    for v in values(output_files)
        rm(v)
    end
end


text:
using StatsBase
using StatsBase: PValue
using Test, Random

v1 = [1.45666, -23.14, 1.56734e-13]
v2 = ["Good", "Great", "Bad"]
v3 = [1, 56, 2]
v4 = [-12.56, 0.1326, 2.68e-16]
v5 = [0.12, 0.3467, 1.345e-16]
@test sprint(show, CoefTable(Any[v1, v2, v3, v4, v5],
    ["Estimate", "Comments", "df", "t", "p"],
    ["x1", "x2", "x3"], 5, 4)) == """
───────────────────────────────────────────────
         Estimate  Comments  df       t       p
───────────────────────────────────────────────
x1    1.45666         Good    1  -12.56  0.1200
x2  -23.14            Great  56    0.13  0.3467
x3    1.56734e-13     Bad     2    0.00  <1e-15
───────────────────────────────────────────────"""

Random.seed!(10)
m = rand(3,4)
@test sprint(show, CoefTable(m, ["Estimate", "Stderror", "df", "p"], [], 4)) == """
──────────────────────────────────────────
     Estimate   Stderror        df       p
──────────────────────────────────────────
[1]  0.112582  0.0566454  0.381813  0.8198
[2]  0.368314  0.120781   0.815104  0.6699
[3]  0.344454  0.179574   0.242208  0.4531
──────────────────────────────────────────"""

@test sprint(show, PValue(1.0)) == "1.0000"
@test sprint(show, PValue(1e-1)) == "0.1000"
if VERSION > v"1.6.0-DEV"
    @test sprint(show, PValue(1e-5)) == "<1e-04"
else
    @test sprint(show, PValue(1e-5)) == "<1e-4"
end
@test sprint(show, PValue(NaN)) == "NaN"
@test_throws ErrorException PValue(-0.1)
@test_throws ErrorException PValue(1.1)
@test PValue(PValue(0.05)) === PValue(0.05)

@test sprint(showerror, ConvergenceException(10)) == "failure to converge after 10 iterations."

@test sprint(showerror, ConvergenceException(10, 0.2, 0.1)) ==
    "failure to converge after 10 iterations. Last change (0.2) was greater than tolerance (0.1)."

@test sprint(showerror, ConvergenceException(10, 0.2, 0.1, "Try changing maxIter.")) ==
    "failure to converge after 10 iterations. Last change (0.2) was greater than tolerance (0.1). Try changing maxIter."

err = @test_throws ArgumentError ConvergenceException(10,.1,.2)
@test err.value.msg == "Change must be greater than tol."


text:
# https://github.com/nalimilan/FreqTables.jl #
##############################################

import Pkg
Pkg.instantiate() # downloads all dependencies for the current project
Pkg.add("FreqTables")

using DataFrames, CSV, FreqTables, DelimitedFiles

df = CSV.read("../externals/Core.Math.Data/data/Pejcic_318.csv", copycols = true)

println(df.MJESTO)

af = freqtable(df.MJESTO)
println(af)

rf = prop(af)
println(rf)

af_2 = freqtable(df.SPOL, df.RAZRED)
println(af_2)


text:
using Mustache
using Test

@testset " inverted " begin

	## Falsey sections should have their contents rendered.
tpl = """\"{{^boolean}}This should be rendered.{{/boolean}}\""""

	@test Mustache.render(tpl, Dict{Any,Any}("boolean"=>false)) == """\"This should be rendered.\""""

	## Truthy sections should have their contents omitted.
tpl = """\"{{^boolean}}This should not be rendered.{{/boolean}}\""""

	@test Mustache.render(tpl, Dict{Any,Any}("boolean"=>true)) == """\"\""""

	## Objects and hashes should behave like truthy values.
tpl = """\"{{^context}}Hi {{name}}.{{/context}}\""""

	@test Mustache.render(tpl, Dict{Any,Any}("context"=>Dict{Any,Any}("name"=>"Joe"))) == """\"\""""

	## Lists should behave like truthy values.
tpl = """\"{{^list}}{{n}}{{/list}}\""""

	@test Mustache.render(tpl, Dict{Any,Any}("list"=>Dict{Any,Any}[Dict("n"=>1), Dict("n"=>2), Dict("n"=>3)])) == """\"\""""

	## Empty lists should behave like falsey values.
tpl = """\"{{^list}}Yay lists!{{/list}}\""""

	@test Mustache.render(tpl, Dict{Any,Any}("list"=>Any[])) == """\"Yay lists!\""""

	## Multiple inverted sections per template should be permitted.
tpl = """{{^bool}}
* first
{{/bool}}
* {{two}}
{{^bool}}
* third
{{/bool}}
"""

	@test Mustache.render(tpl, Dict{Any,Any}("two"=>"second","bool"=>false)) == """* first
* second
* third
"""

	## Nested falsey sections should have their contents rendered.
tpl = """| A {{^bool}}B {{^bool}}C{{/bool}} D{{/bool}} E |"""

	@test Mustache.render(tpl, Dict{Any,Any}("bool"=>false)) == """| A B C D E |"""

	## Nested truthy sections should be omitted.
tpl = """| A {{^bool}}B {{^bool}}C{{/bool}} D{{/bool}} E |"""

	@test Mustache.render(tpl, Dict{Any,Any}("bool"=>true)) == """| A  E |"""

	## Failed context lookups should be considered falsey.
tpl = """[{{^missing}}Cannot find key 'missing'!{{/missing}}]"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """[Cannot find key 'missing'!]"""

	## Dotted names should be valid for Inverted Section tags.
tpl = """\"{{^a.b.c}}Not Here{{/a.b.c}}\" == \"\""""

	@test Mustache.render(tpl, Dict{Any,Any}("a"=>Dict{Any,Any}("b"=>Dict{Any,Any}("c"=>true)))) == """\"\" == \"\""""

	## Dotted names should be valid for Inverted Section tags.
tpl = """\"{{^a.b.c}}Not Here{{/a.b.c}}\" == \"Not Here\""""

	@test Mustache.render(tpl, Dict{Any,Any}("a"=>Dict{Any,Any}("b"=>Dict{Any,Any}("c"=>false)))) == """\"Not Here\" == \"Not Here\""""

	## Dotted names that cannot be resolved should be considered falsey.
tpl = """\"{{^a.b.c}}Not Here{{/a.b.c}}\" == \"Not Here\""""

	@test Mustache.render(tpl, Dict{Any,Any}("a"=>Dict{Any,Any}())) == """\"Not Here\" == \"Not Here\""""

	## Inverted sections should not alter surrounding whitespace.
tpl = """ | {{^boolean}}	|	{{/boolean}} |
"""

	@test Mustache.render(tpl, Dict{Any,Any}("boolean"=>false)) == """ | 	|	 |
"""

	## Inverted should not alter internal whitespace.
tpl = """ | {{^boolean}} {{! Important Whitespace }}
 {{/boolean}} |
"""

        @test Mustache.render(tpl, Dict{Any,Any}("boolean"=>false)) == " |  \n  |\n"
#""" |\\s
#  |
#"""

	## Single-line sections should not alter surrounding whitespace.
tpl = """ {{^boolean}}NO{{/boolean}}
 {{^boolean}}WAY{{/boolean}}
"""

	@test Mustache.render(tpl, Dict{Any,Any}("boolean"=>false)) == """ NO
 WAY
"""

	## Standalone lines should be removed from the template.
tpl = """| This Is
{{^boolean}}
|
{{/boolean}}
| A Line
"""

	@test Mustache.render(tpl, Dict{Any,Any}("boolean"=>false)) == """| This Is
|
| A Line
"""

	## Standalone indented lines should be removed from the template.
tpl = """| This Is
  {{^boolean}}
|
  {{/boolean}}
| A Line
"""

	@test Mustache.render(tpl, Dict{Any,Any}("boolean"=>false)) == """| This Is
|
| A Line
"""

	## "\r\n" should be considered a newline for standalone tags.
tpl = """|
{{^boolean}}
{{/boolean}}
|"""

	@test Mustache.render(tpl, Dict{Any,Any}("boolean"=>false)) == """|
|"""

	## Standalone tags should not require a newline to precede them.
tpl = """  {{^boolean}}
^{{/boolean}}
/"""

	@test Mustache.render(tpl, Dict{Any,Any}("boolean"=>false)) == """^
/"""

	## Standalone tags should not require a newline to follow them.
tpl = """^{{^boolean}}
/
  {{/boolean}}"""

	@test Mustache.render(tpl, Dict{Any,Any}("boolean"=>false)) == """^
/
"""

	## Superfluous in-tag whitespace should be ignored.
tpl = """|{{^ boolean }}={{/ boolean }}|"""

	@test Mustache.render(tpl, Dict{Any,Any}("boolean"=>false)) == """|=|"""
end


text:
struct Shape
    points::Vector{Point}
    simplified_points::Vector{Point}
    centroid::Point
    centered_points::Vector{Point}
    subpaths::Vector{Vector{Point}}
    num_acute_angles::Int
    num_obtuse_angles::Int
    num_right_angles::Int
end

function EmptyShape()
    return Shape(Point[], Point[], O, Point[], Vector{Vector{Point}}(), 0, 0, 0)
end

function Shape(num_points::Int, points_of_subpaths::Vector{Int})
    points = Vector{Point}(undef, num_points)
    subpaths = Vector{Vector{Point}}()
    for points_in_subpath in points_of_subpaths
        push!(subpaths, Vector{Point}(undef, points_in_subpath))
    end
    return Shape(points, Point[], O, Point[], subpaths, 0, 0, 0)
end

function Shape(points, subpaths; compute_info = true)
    simplified_points = simplify(points)
    if compute_info
        points, aa, oa, ra = get_angles(simplified_points)
    else
        aa = 0
        oa = 0
        ra = 0
    end
    centroid = polycentroid(points)
    shape =
        Shape(points, simplified_points, centroid, points .- centroid, subpaths, aa, oa, ra)
    return shape
end

function Base.isempty(shape::Shape)
    return length(shape.points) == 0
end

function get_angles(p)
    # TODO: let's assume it's a closed path
    num_acute_angles = 0
    num_obtuse_angles = 0
    num_right_angles = 0
    simplified = Point[]

    for i in 1:length(p)
        pA = p[i]
        if i < length(p) - 1
            u = p[i] - p[i + 1]
            v = p[i + 2] - p[i + 1]
            pB = p[i + 1]
            pC = p[i + 2]
        elseif i == length(p) - 1
            u = p[i] - p[i + 1]
            v = p[1] - p[i + 1]
            pB = p[i + 1]
            pC = p[1]
        else
            u = p[i] - p[1]
            v = p[2] - p[1]
            pB = p[1]
            pC = p[2]
        end
        ang = rad2deg(acos(clamp(
            dotproduct(u, v) / (sqrt(u.x^2 + u.y^2) * sqrt(v.x^2 + v.y^2)),
            -1,
            1,
        )))
        if ang < 179.8
            push!(simplified, pB)
        end

        if ang < 160
            if ang > 100 # obtuse
                num_obtuse_angles += 1
            elseif ang < 80 # acute
                num_acute_angles += 1
            else # right
                num_right_angles += 1
            end
        end
    end
    return simplified, num_acute_angles, num_obtuse_angles, num_right_angles
end

function get_similarity(shapeA::Shape, shapeB::Shape)
    if isempty(shapeA) || isempty(shapeB)
        return 0.0
    end
    score_points = 10.0
    score_holes = 500.0
    score_acute_angles = 100.0
    score_right_angles = 100.0
    score_obtuse_angles = 100.0
    score_point_diff = 100.0
    score_centered_point_diff = 100.0

    score = 0.0

    # number of points
    nA = length(shapeA.simplified_points)
    nB = length(shapeB.simplified_points)
    perc = nA < nB ? nA / nB : nB / nA
    perc = isnan(perc) ? 1.0 : perc
    # println("Point score: ", perc*score_points)
    score += perc * score_points

    # number of holes
    nA = length(shapeA.subpaths)
    nB = length(shapeB.subpaths)
    perc = nA < nB ? nA / nB : nB / nA
    perc = isnan(perc) ? 1.0 : perc
    # println("Hole score: ", perc*score_holes)
    score += perc * score_holes

    # number of angles
    num_angles = shapeA.num_acute_angles
    num_angles += shapeB.num_acute_angles
    num_angles += shapeA.num_obtuse_angles
    num_angles += shapeB.num_obtuse_angles
    num_angles += shapeA.num_right_angles
    num_angles += shapeB.num_right_angles

    # TODO: think about how to compare this better
    # Probably dividing by the whole number of angles is dumb
    # Or one might need a score for the difference in number of angles
    # acute
    nA = shapeA.num_acute_angles
    nB = shapeB.num_acute_angles
    perc = 1 - abs(nA - nB) / num_angles
    perc = isnan(perc) ? 1.0 : perc
    score += perc * score_acute_angles
    # println("Acute score: ", perc*score_acute_angles)

    # obtuse
    nA = shapeA.num_obtuse_angles
    nB = shapeB.num_obtuse_angles
    perc = 1 - abs(nA - nB) / num_angles
    perc = isnan(perc) ? 1.0 : perc
    score += perc * score_obtuse_angles
    # println("Obtuse score: ", perc*score_obtuse_angles)

    # right
    nA = shapeA.num_right_angles
    nB = shapeB.num_right_angles
    perc = 1 - abs(nA - nB) / num_angles
    perc = isnan(perc) ? 1.0 : perc
    score += perc * score_right_angles
    # println("right score: ", perc*score_right_angles)

    # difference in centered_points
    pointsA, pointsB = match_num_points(shapeA.centered_points, shapeB.centered_points)
    smallest_i, smallest_distance = compute_shortest_morphing_dist(pointsA, pointsB)
    # TODO: maybe there is a more reasonable denominator
    smallest_distance /= length(pointsA)
    wmin, wmax = extrema(p -> p.x, pointsA)
    hmin, hmax = extrema(p -> p.x, pointsA)
    w = wmax - wmin
    h = hmax - hmin
    perc = clamp(1 - smallest_distance / sqrt(w^2 + h^2), 0, 1)
    score += perc * score_centered_point_diff

    # difference in movement
    pointsA, pointsB = match_num_points(shapeA.points, shapeB.points)
    smallest_i, smallest_distance = compute_shortest_morphing_dist(pointsA, pointsB)
    video = CURRENT_VIDEO[1]
    smallest_distance /= length(pointsA)
    perc = clamp(1 - smallest_distance / sqrt(video.width^2 + video.height^2), 0, 1)
    score += perc * score_centered_point_diff

    return score
end

function create_shapes(polys)
    shapes = Vector{Shape}()

    is_last_subpath = false
    current_points = Point[]
    current_subpaths = Vector{Vector{Point}}()
    for poly in polys
        if ispolyclockwise(poly) && !is_last_subpath
            empty!(current_subpaths)
            if !isempty(current_points)
                shape = Shape(current_points, Vector{Vector{Point}}())
                push!(shapes, shape)
            end
            current_points = poly
            is_last_subpath = false
        elseif ispolyclockwise(poly)
            shape = Shape(current_points, copy(current_subpaths))
            push!(shapes, shape)
            is_last_subpath = false
            current_points = poly
            empty!(current_subpaths)
        else # is a hole
            push!(current_subpaths, poly)
            is_last_subpath = true
        end
    end
    shape = Shape(current_points, copy(current_subpaths))
    push!(shapes, shape)
    return shapes
end

function print_basic(s::Shape)
    println("Shape: #Points: $(length(s.points))")
    println("Angles: #Acute: $(s.num_acute_angles) #Obtuse: $(s.num_obtuse_angles) #Right: $(s.num_right_angles)")
    println("#Holes: $(length(s.subpaths))")
end

function prepare_to_interpolate(from_shape, to_shape)
    # match number of points for outer polygon
    from_outer, to_outer = match_num_points(from_shape.points, to_shape.points)

    # for holes
    from_holes = Vector{Vector{Point}}()
    to_holes = Vector{Vector{Point}}()

    for (from_hole, to_hole) in zip(from_shape.subpaths, to_shape.subpaths)
        from_hole_matched, to_hole_matched = match_num_points(from_hole, to_hole)
        push!(from_holes, from_hole_matched)
        push!(to_holes, to_hole_matched)
    end

    # rotate outer polygon
    rotate_i, _ = compute_shortest_morphing_dist(from_outer, to_outer)
    new_from_outer = circshift(from_outer, -rotate_i + 1)

    new_from_holes = Vector{Vector{Point}}()

    # rotate holes
    for (from_hole, to_hole) in zip(from_holes, to_holes)
        rotate_i, _ = compute_shortest_morphing_dist(from_hole, to_hole)
        push!(new_from_holes, circshift(from_hole, -rotate_i + 1))
    end

    from = Shape(new_from_outer, new_from_holes; compute_info = false)
    to = Shape(to_outer, to_holes; compute_info = false)

    return from, to
end

function interpolate_shape!(inter_shape, from, to, t)
    # outer
    for (i, p1, p2) in zip(1:length(from.points), from.points, to.points)
        new_point = p1 + t * (p2 - p1)
        inter_shape.points[i] = new_point
    end

    # for holes
    for hi in 1:length(from.subpaths)
        for (i, p1, p2) in
            zip(1:length(from.subpaths[hi]), from.subpaths[hi], to.subpaths[hi])
            new_point = p1 + t * (p2 - p1)
            inter_shape.subpaths[hi][i] = new_point
        end
    end
end

function draw_shape(shape, draw_action)
    newpath()
    poly(shape.points, :path; close = true)
    for hi in 1:length(shape.subpaths)
        newsubpath()
        poly(shape.subpaths[hi], :path; close = true)
    end
    do_action(draw_action)
end


text:
bstr = [path * "/SampleFiles/Bottle/B024*",
        path * "/SampleFiles/Bottle/B2021901700*",
        path * "/SampleFiles/Bottle/B203*",
        path * "/SampleFiles/Bottle/B20319115Rainfallmm"]

printstyled("  Bottle (UNAVCO strain data)\n", color=:light_green)
for i = 1:4
  S = verified_read_data("bottle", bstr[i], nx_new=14400, nx_add=14400, vl=true)
  fill_pbo!(S)
  if i == 4
    @test S.name[1] == "quarry203bwa2007"
    @test isapprox(S.fs[1], 1.0/(30*60))
    @test isapprox(S.loc[1].el, 814.4)
    @test S.units[1] == "mm"
  end
end

printstyled("    source logging\n", color=:light_green)
redirect_stdout(out) do
  S = SeisIO.read_bottle(bstr[4], 14400, 14400, true, true, 0)
  show_src(S, 1)
  show_src(S[1])
  show_src(S)
end


text:
nums = Int32[]
for line in eachline()
    push!(nums, parse(Int32, line))
end

for num1 in nums
    for num2 in nums
        for num3 in nums
            if num1 + num2 + num3 == 2020
                println(num1 * num2 * num3)
                return
            end
        end
    end
end


text:
####################################################################
# Basic Data Structures
####################################################################

#------PDAGs------#

#PDAG is parametric type because A could be a BitArray or a SubArray
"""
Define a Partially Directed Acyclic Graph (PDAG).

PDAGs are defined by the number of vertices and a Boolean adjacency matrix
"""
mutable struct PDAG{S<:AbstractMatrix{Bool}}
    nv::Int #number of vertices
    A::S #adjacency matrix to represent graph
end

"""
    PDAG(x::Integer)
Create an empty graph with `x` vertices.
"""
PDAG(nv::Integer) = PDAG(nv, falses(nv,nv))

"""
    PDAG(numNodes, numEdges; seed=-1)
Generate a random PDAG with `numNodes` vertices and `numEdges` edges.
"""
function PDAG(numNodes::Integer, numEdges::Integer; seed::Int=-1)

    #Calculate the maximum number of edges
    maxEdges = numNodes*(numNodes - 1) ÷ 2
    @assert(numEdges <= maxEdges, "Maximum number of edges for this graph is $maxEdges")

    #If no random seed is given, use the default RNG
    rng = seed==-1 ? Xoshiro() : Xoshiro(seed);

    #Generate a empty graph
    g = PDAG(numNodes)

    #Continually add edges until the correct number is reached
    edgesAdded = 0
    while edgesAdded < numEdges
        src = rand(rng, 1:numNodes)
        dest = rand(rng, 1:numNodes)
        orientation = rand(rng, Bool)

        #Do not allow self loops and do not over-write edges
        if src != dest && !isadjacent(g, src, dest)
            orientation ? addedge!(g, src, dest, directed=false) : addedge!(g, src, dest)
            edgesAdded += 1
        end
    end

    #Return the random PDAG
    return g
end

#------Edges------#

#A data type to hold edge Information
# parent → child
"""
    Edge(parent, child, directed::Bool)
    Edge(Tuple(parent,child), directed::Bool)
    Edge(CartesianIndex(parent,child), directed::Bool)
A data type to hold edge information.

It is also possible to convert a `Tuple` or `CartesianIndex` into an Edge.
"""
struct Edge
    parent::Int
    child::Int
    directed::Bool
end

#converting edges from other types
Edge(idx::Tuple{Int, Int}, directed) = Edge(idx...,directed)
Edge(idx::CartesianIndex{2}, directed) = Edge(Tuple(idx),directed)

"A data type used to iterate through the edges of the graph"
struct IterEdge{G<:PDAG}
    g::G
end

####################################################################
# Overload Base functions
####################################################################

#Used to subset the graph
#nodes contains a list of nodes you want to keep in the graph
#no need to make a copy of the graph using views!
view(g::PDAG, nodes) = PDAG(length(nodes), view(g.A, nodes, nodes))

#An element of an IterEdge is an Edge (gives functions like collect() information about the element type)
eltype(::Type{IterEdge{PDAG{T}}}) where T = Edge

#knowing the number of edges will allow us to pre-allocate a correctly sized collection of edges
#ne() is a function defined below to get the number of edges
length(iterEdge::IterEdge) = ne(iterEdge.g)

#overload to print PDAG (requires nv() and ne() which are functions defined below)
function show(io::IO, g::PDAG)
    print(io, "{$(nv(g)), $(ne(g))} PDAG")
end

#overload to print Edge type
function show(io::IO, edge::Edge)
    if edge.directed
        print(io,"$(edge.parent) → $(edge.child)")
    else
        print(io,"$(edge.parent) - $(edge.child)")
    end
end

#check if two graphs are equal by compares adjacency matrices
==(g₁::PDAG, g₂::PDAG) = g₁.A == g₂.A

####################################################################
# Number of Vertices and Edges
####################################################################

"Return the number of vertices"
nv(g::PDAG) = g.nv

"Return the number of edges"
function ne(g::PDAG)
    A = g.A
    totalEdges = 0
    #basically count each edge twice because some edges appear only on upper/lower diagonal but undirected edges appear on both
    for i in eachindex(A)
        if A[i] && A'[i]  #undirected
            totalEdges += 1
        elseif A[i] #directed
            totalEdges += 2
        end
    end

    return totalEdges ÷ 2 #integer division
end

####################################################################
# is* functions
####################################################################

"""
    isneighbor(g::PDAG, x, y)
    isneighbor(g::PDAG, edge::Edge)
Test if `x` and `y` are connected by a undirected edge in the graph `g`.

Can also perform the same test given an `edge`.
"""
isneighbor(g::PDAG, x, y) = @inbounds g.A[x,y] && g.A[y,x]
isneighbor(g::PDAG, edge::Edge) = isneighbor(g, edge.parent, edge.child)

"""
    isadjacent(g::PDAG, x, y)
    isadjacent(g::PDAG, edge::Edge)
Test if `x` and `y` are connected by a any edge in the graph `g`.

Can also perform the same test given an `edge`.
"""
isadjacent(g::PDAG, x, y) = @inbounds g.A[x,y] || g.A[y,x]
isadjacent(g::PDAG, edge::Edge) = isadjacent(g, edge.parent, edge.child)

"""
    isparent(g::PDAG, x, y)
    isparent(g::PDAG, edge::Edge)
Test if `x` is a parent of `y` in the graph `g`.

Can also perform the same test given an `edge`.
"""
isparent(g::PDAG, x, y) = @inbounds g.A[x,y] && !g.A[y,x]
isparent(g::PDAG, edge::Edge) = isparent(g, edge.parent, edge.child)

"""
    ischild(g::PDAG, x, y)
    ischild(g::PDAG, edge::Edge)
Test if `x` is a child of `y` in the graph `g`.

Can also perform the same test given an `edge`.
"""
ischild(g::PDAG, x, y) = @inbounds !g.A[x,y] && g.A[y,x]
ischild(g::PDAG, edge::Edge) = ischild(g, edge.parent, edge.child)

"""
    isdescendent(g::PDAG, x, y)
Return `true` if `x`←`y` OR `x`-`y` in the graph `g`.
"""
isdescendent(g::PDAG, x, y) = @inbounds g.A[y,x]

"""
    isclique(g::PDAG, nodes)
Return `true` if all vertices in `nodes` are connected to each other in the graph `g`.
"""
function isclique(g::PDAG, nodes)

    for (v₁, v₂) in allpairs(nodes)
        if !isadjacent(g, v₁, v₂)
            return false
        end
    end

    return true
end

#Do the nodesRemoved block all semi-directed paths between src and dest?
"""
    isblocked(g::PDAG, src, dest, nodesRemoved)
Return `true` if there is no semi-directed path between `src` and `dest` in the graph `g`.

A set of vertices (`nodesRemoved`) can be removed from the graph before searching for a semi-directed path.

A semi-directed path between `src` and `dest` is a list of edges in `g` where every edge is either undirected or points toward `dest`.

    src → x₁ - x₂ → dest ✓
    src → x₁ ← x₂ - dest ✖
"""
function isblocked(g::PDAG, src, dest, nodesRemoved)

    #Keep track of all the nodes visited
    visited = zeros(Bool, nv(g))

    # mark excluded vertices as visited
    for vᵢ in nodesRemoved
        visited[vᵢ] = true
    end

    #If src or dest were in nodesRemoved, the path is blocked
    (visited[src] || visited[dest]) && return true

    #if the src and dest are the same, path is itself
    src == dest && return false

    #check if scr or dest have no neighbors
    if countNeighbors(g, src)==0 || countNeighbors(g, dest)==0
        return true
    end

    queue = [src]
    visited[src] = true
    while !isempty(queue)
        currentNode = popfirst!(queue) # get new element from queue
        for vᵢ in descendent(g,currentNode)
            vᵢ == dest && return false
            if !visited[vᵢ]
                push!(queue, vᵢ) # push onto queue
                visited[vᵢ] = true
            end
        end
    end
    return true
end

"Return true if `edge` is directed"
isdirected(edge::Edge) = edge.directed

####################################################################
# Modify Edges
####################################################################

"""
    addedge!(g::PDAG, x, y; directed=true)
    addedge!(g::PDAG, edge::Edge)
Add the edge `x`→`y` to the graph `g`.

Can also perform the same addition given an `edge`.

By default, this function will direct the new edge from `x` to `y` unless the `directed` keyword is changed to `false`.
"""
function addedge!(g::PDAG, x, y; directed=true)
    g.A[x,y] = true
    g.A[y,x] = !directed
    return nothing
end

#If an edge type if given, call the previous function
addedge!(g::PDAG, edge::Edge) = addedge!(g::PDAG, edge.parent, edge.child, directed=edge.directed)

"""
    remedge!(g::PDAG, x, y)
    remedge!(g::PDAG, edge::Edge)
Remove any edge between `x` and `y` in the graph `g`.

Can also perform the same deletion given an `edge`.

No warning is given if `edge` is not in the graph.
"""
function remedge!(g::PDAG, x, y)
    g.A[x,y] = false
    g.A[y,x] = false
    return nothing
end

#remove edge (Edge type)
remedge!(g::PDAG, edge::Edge) = remedge!(g, edge.parent, edge.child)

"""
    orientedge!(g::PDAG, x, y)

Update the edge `x`-`y` to `x`→`y` in the graph `g`.
"""
function orientedge!(g::PDAG, x, y)
    g.A[y,x] = false
    return nothing
end

####################################################################
# Neighborhood functions
####################################################################

#Should these be non-allocating iterators?
neighborsGeneral(isFunction::Function, g::PDAG, x) = Iterators.filter(vᵢ -> isFunction(g,vᵢ,x), vertices(g))

#Loop through all the vertices and perform some neighborhood test on the given vertex
#Output all vertices that pass the given test
# function neighborsGeneral(isFunction::Function, g::PDAG, x)
#     neighborList = Int64[]
#     for vᵢ in vertices(g)
#         #check depends on desired list of neighbors
#         #e.g. isFunction == isParent when we're searching for vᵢ→x
#         if isFunction(g,vᵢ,x)
#             push!(neighborList,vᵢ)
#         end
#     end

#     return neighborList
# end

#All vertices connected to x
neighbors(g::PDAG, x) = neighborsGeneral(isadjacent, g, x)

#All vertices pointing to x
neighbors_in(g::PDAG, x) = neighborsGeneral(isparent, g, x)

#All vertices pointing away from x
neighbors_out(g::PDAG, x) = neighborsGeneral(ischild, g, x)

#All vertices connected to x by an undirected edge
neighbors_undirect(g::PDAG, x) = neighborsGeneral(isneighbor, g, x)

#these are just aliases for the functions above
parents(g::PDAG, x) = neighbors_in(g, x)
children(g::PDAG, x) = neighbors_out(g, x)

descendent(g::PDAG, x) = neighborsGeneral(isdescendent, g, x)

####################################################################
# Counting Functions
####################################################################

#Very similar to the neighborsGeneral function except we only keep a count of the valid nodes
function countGeneral(isFunction::Function, g::PDAG,x)
    counter=0
    for vᵢ in vertices(g)
        if isFunction(g,vᵢ,x)
            counter += 1
        end
    end
    return counter
end

#All vertices connected to x
countNeighbors(g::PDAG, x) = countGeneral(isadjacent, g, x)

#All vertices pointing to x
countNeighbors_in(g::PDAG, x) = countGeneral(isparent, g, x)

#All vertices pointing away from x
countNeighbors_out(g::PDAG, x) = countGeneral(ischild, g, x)

#All vertices connected to x by an undirected edge
countNeighbors_undirect(g::PDAG, x) = countGeneral(isneighbor, g, x)

#these are just aliases for the functions above
countParents(g::PDAG, x) = countNeighbors_in(g, x)
countChildren(g::PDAG, x) = countNeighbors_out(g, x)

####################################################################
# all* functions
####################################################################

#Get every undirected edge in the entire graph g
allundirected(g::PDAG) = Edge.(findall(tril(g.A .& g.A')),false)

#Get every directed edge in the entire graph g
alldirected(g::PDAG) = [Edge(i,true) for i in CartesianIndices(g.A) if g.A[i] && !g.A'[i]]

#combine and sort all directed and undirected edges
alledges(g::PDAG) = sort([allundirected(g); alldirected(g)], by=x->(x.parent, x.child))

#Generate an iterator equivalent to combinations(x,2) but with base packages
#The combinations are not in order, but that doesn't matter in this case
#benefit: packages have more methods defined for the base iterators
allpairs(v) = Iterators.filter(i -> isless(i...), Iterators.product(v,v))

####################################################################
# Iterators
####################################################################

#Given an initial (i,j) index, output the next non-zero index in adjacency matrix
function iterate(iterEdge::IterEdge, state=(1,1))
    g = iterEdge.g
    n = nv(g)
    (row,col) = state

    #loop starting at the current state (row,col)
    for i = row:n
        start = row==i ? (col+1) : (i+1)
        for j = start:n

            #set the next potential state
            state = (i,j)

            #Check for an edge
            if g.A[i,j] && g.A'[i,j] #undirected
                return (Edge(state,false), state)
            elseif g.A[i,j] #forward edge
                return (Edge(state,true), state)
            elseif g.A'[i,j] #backward edge
                return (Edge(reverse(state),true), state)
            end
        end
    end

    return nothing
end

"""
    edges(g::PDAG)

Return an iterator to generate all edges within the graph `g`. Use
"""
edges(g::PDAG) = IterEdge(g)

"""
    vertices(g::PDAG)

Return an iterator to generate all vertices within the graph `g`.
"""
vertices(g::PDAG) = Base.OneTo(nv(g))

####################################################################
# Misc
####################################################################

# for x-y, get undirected neighbors of y and any neighbor of x
#calcNAyx(g::PDAG, y::Integer, x::Integer) = setdiff(neighbors_undirect(g,y) ∩ neighbors(g,x), x)
function calcNAyx(g::PDAG, y, x)
    neighborList = Int64[]

    #loop through all vertices except for x
    for vᵢ in Iterators.filter(!isequal(x),vertices(g))
        #look for neighbors of y and adjacencies to x
        if isneighbor(g,vᵢ,y) && isadjacent(g,vᵢ,x)
            push!(neighborList,vᵢ)
        end
    end

    return neighborList
end
calcNAyx(g::PDAG, edge::Edge) = calcNAyx(g, edge.child, edge.parent)

#for x-y, undirected neighbors of y not connected to x
#calcT(g::PDAG, y::Integer, x::Integer) = setdiff(neighbors_undirect(g,y), neighbors(g,x), x)
function calcT(g::PDAG, y, x)
    neighborList = Int64[]

    #loop through all vertices except for x
    for vᵢ in Iterators.filter(!isequal(x),vertices(g))
        #look for neighbors of y and adjacencies to x
        if isneighbor(g,vᵢ,y) && !isadjacent(g,vᵢ,x)
            push!(neighborList,vᵢ)
        end
    end

    return neighborList
end
calcT(g::PDAG, edge::Edge) = calcT(g, edge.child, edge.parent)

#Single pass algoirthm to calculate the average degree
function degreeAverage(g::PDAG)

    ave = 0.0

    for (i,vᵢ) in enumerate(vertices(g))
        ave += (countNeighbors(g, vᵢ) - ave) / i
    end

    return ave
end

####################################################################
# Save/load graphs
####################################################################

"""
    saveGraph(fileName, g::PDAG)

Save the output of fges to a text file.
"""
function saveGraph(fileName,g::PDAG)
    open(fileName,"w") do io
        for edge in edges(g)
            println(io,edge)
        end
    end
end

"""
    loadGraph(fileName)

Load a saved output of fges from a text file.
"""
function loadGraph(fileName)
    #find the largest number in the file
    numFeatures = maximum([maximum(parse.(Int, split(line," ")[[1,3]])) for line in eachline(fileName)])

    g = PDAG(numFeatures)

    for line in eachline(fileName)
        sublines = split(line," ")

        v₁ = parse(Int, sublines[1])
        v₂ = parse(Int, sublines[3])
        edgeDirected = sublines[2] == "→" ? true : false

        addedge!(g, v₁, v₂, directed=edgeDirected)
    end

    return g
end

####################################################################
# Convert graph to table
####################################################################

#Well that was easy, thanks Tables.jl
edgetable(g::PDAG) = DataFrame(edges(g))


text:
# functions to set parameters and to creat inital values

function q_set_para_initial()
    # set parameters

    # model type and shock
    type = 0 # 0 - rep agent ;; do not change
    action_type = 5 # type of action function - 0 constant /  1 linear / 2 polynomial.... see setup_functions.jl
    market_type = 0 # not relevant for representative agent model
    shock_type = 2 # 0 nothing , 1 postitive productivity, 2 negative productivity
    t_shock = 20 #0 default; otherwise period at which shock happens -- if bigger than periods nothing happens -- set shock duration below
    shock_duration = 5  # 1 to something , number of periods shock is lasting -- do not set to zero
    r_seed = 123

    # simulation parameters
    n_agents = 5000  # number of agents -- default 5000
    periods = 50  # number of periods in each generation --- default 50 for shock, 200 without shock
    generations = 100 # number of gnerations default 1000
    reset_capital = 1   # if =1 capital reset for each new generation to populations average
    sim_set = Dict("n_agents" => n_agents, "periods" => periods, "generations" => generations, "r_seed" => r_seed,
    "type"=> type, "action_type" => action_type, "market_type" => market_type , "t_shock" => t_shock, "shock_duration" => shock_duration, "shock_type" => shock_type )

    # choice and learning parameters
    scale = 10   # inside decsion function - maps (0;1) parameters to (-scale;+scale) used inside sigmoid function, default 10
    mutation_prob = .01  # probablity that new learning parmeter is mutated, default .01
    decide_set = Dict("scale" => scale, "mutation_prob" => mutation_prob )

    # economy parameters
    prod_level = 1   # productivity level, default 1
    prod_trans_type = 1 # type of productivity shock: 0 discrete ; 1 normal  -- for standard growth model: 1
    product_trans_p = [.9, 0.1]  # low to high, high to high --- or for transition type 1  roh and sigma --  roh persitence -- sigma variance of new draw, default .9 and .1
    alpha = .25  # parameter for for production function
    beta = 1   # discount rate utility
    delta = 0.5      # depreciation capital  0 full depreciation 1 nothing depreciates
    price_para = [0, .1, 0.1]  # not relevant for representative agent model
    print  = 0 # not relevant for representative agent model
    econ_set = Dict("prod_level"=> prod_level, "prod_trans_type"=> prod_trans_type, "productivity_transition"=> product_trans_p, "alpha" => alpha, "beta" => beta, "delta" => delta, "reset_capital" => reset_capital, "price_para" => price_para, "print_money" => print )

    return sim_set, decide_set, econ_set
end

function q_set_accounting(sim_set)
    # set up array to keep track of simulation results
    type = sim_set["type"]
    action_type = sim_set["type"]
    market_type =sim_set["type"]
    # accounting
    # economy features recordecd every period
    # output, capital, consumption, save_rate
    mean_output = zeros(0)
    mean_consume = zeros(0)
    mean_capital = zeros(0)
    mean_s_goods = zeros(0)
    mean_s_cash = zeros(0)
    mean_price = zeros(0)
    mean_cash = zeros(0)

    account_econ = [mean_output, mean_consume, mean_capital, mean_s_goods, mean_s_cash, mean_price, mean_cash]

    # parameters for decsions recorded every generation
    account_para = zeros(0,n_decide_para)

    return account_econ, account_para
end

function q_set_initial_values(type, n_agents)
    # set initial values
    productivity = ones(n_agents)
    capital = ones(n_agents)
    cash = ones(n_agents) * 20
    price = ones(n_agents)
    s_cash = ones(n_agents) * .5   # save rate cash
    s_goods = ones(n_agents) * .5
    if type == 1
        state = [capital, productivity, cash, price, s_cash, s_goods]
    elseif type == 2
        state = [productivity, cash, price, s_cash]
    else
        state = [capital, productivity]
    end
    decide_para = q_draw_para(n_agents, n_decide_para)
    initial_set = Dict( "state" => state, "decide_para" => decide_para)
return initial_set
end


text:
# GO

TESTING = false

blank = 1
black = 2
white = 3
stone = [".", "⚈", "⚆"]

size = 9

state = ones(Int8, size*size)
move = black
passes = 0
last_captured = []

function msg(s)
    if !TESTING
        println(s)
    end
end

function set_game(blacks, whites, moveof)
    new_game()
    global move = moveof
    for m in blacks
        col = Int(m[1]) - 96
        row = size - parse(Int8, m[2])
        state[row*size+col] = black
    end
    for m in whites
        col = Int(m[1]) - 96
        row = size - parse(Int8, m[2])
        state[row*size+col] = white
    end
end

function set_move(moveof)
    global move = moveof
end

function new_game()
    global state = ones(Int8, size*size)
    global move = black
    global passes = 0
end

function display()
    println()
    println("  a  b  c  d  e  f  g  h  i")
    for j = 0:size-1
        print(size-j, " ")
        for i = 1:size
            print(stone[state[j*size+i]])

            if i != size
                print("  ")
            end
        end
        print(" ", size-j)
        println()
        print("  ")
        if j != size-1
            for i = 1:size
                print("   ")
            end
            println()
        end
    end
    println("a  b  c  d  e  f  g  h  i")
end

function is_in_bounds(row, col)
    return col>0 && col<=size && row>=0 && row<size
end

function is_surrounded(row, col, moveof, vg = Nothing)
    if vg === Nothing vg = zeros(size*size) end
    oppo = moveof==black ? white : black

    idx = row*size+col
    vg[idx] = 1 # this cell has been visited

    # quick check: this cell is surrounded
    if (col+1>size || state[row*size+col+1]==oppo) &&
        (col-1<=0 || state[row*size+col-1]==oppo) &&
        (row+1>=size || state[(row+1)*size+col]==oppo) &&
        (row-1<0 || state[(row-1)*size+col]==oppo)
        return true
    end

    # recursively check adjacent cells
    adj = [[-1,0], [1,0], [0, -1], [0, 1]]
    for a in adj
        acol = col+a[1]
        arow = row+a[2]
        if !is_in_bounds(arow, acol) continue end

        # dont check already visited ones
        if vg[arow*size+acol] == 1 continue end

        # found a liberty?
        if state[arow*size+acol] == blank return false end

        # opposition stone?
        if state[arow*size+acol] == oppo continue end

        # same colored stone, so recurse
        if !is_surrounded(arow, acol, moveof, vg) return false end
    end

    return true
end

function is_capturing(row, col)
    oppo = move==black ? white : black
    vg = zeros(size*size)
    vg[row*size+col] = 1 # do not check this cell
    # recursively check adjacent cells
    adj = [[-1,0], [1,0], [0, -1], [0, 1]]
    for a in adj
        acol = col+a[1]
        arow = row+a[2]
        if !is_in_bounds(arow, acol) continue end
        if is_surrounded(arow, acol, oppo, vg) return true end
    end
    return false
end

function is_legal(row, col)
    if !is_in_bounds(row, col) return false end
    if state[row*size+col] != blank return false end
    if is_surrounded(row, col, move)
        if is_capturing(row, col)
            if length(last_captured) == 1 && last_captured[1] == row*size+col
                msg("KO not allowed!")
                return false
            else return true end
        else
            msg("Move is illigal (captured). Also not a capturing move!")
            return false
        end
    end

    return true
end

function stone_on(m)
    col = Int(m[1]) - 96
    row = size - parse(Int8, m[2])
    return state[row*size+col]
end

function update_board()
    empty!(last_captured)

    # mark stones to remove
    for row = 0:size-1
        for col = 1:size
            if state[row*size+col]==move && is_surrounded(row, col, move)
                push!(last_captured, row*size+col)
            end
        end
    end

    # remove stones
    for i in last_captured
        state[i] = blank
    end
end

function process_input(m)
    col = Int(m[1]) - 96
    row = size - parse(Int8, m[2])
    if is_legal(row, col)
        state[row*size+col] = move

        global move = move==black ? white : black
        global passes = 0

        update_board()

        return true
    end
    return false
end

function search_pocket(p, row, col, vg)
    idx = row*size+col
    if state[idx]==blank
        push!(p, idx)
    else
        return state[idx]
    end

    vg[idx] = 1
    m = state[idx]
    adj = [[-1,0], [1,0], [0, -1], [0, 1]]
    for a in adj
        acol = col+a[1]
        arow = row+a[2]
        if !is_in_bounds(arow, acol) continue end
        if vg[arow*size+acol] == 1 continue end
        ms = search_pocket(p, arow, acol, vg)
        if ms!=m && ms!=blank && m!=blank
            return 4
        end
        if ms!=blank
            m = ms
        end
    end
    return m
end

function calculate_score()
    score = [0, 0, 0, 0] # blank/black/white scores
    vg = zeros(size*size)
    for row = 0:size-1
        for col = 1:size
            idx = row*size+col
            if vg[idx]==1 continue end
            if state[idx] == blank
                pocket = []
                m = search_pocket(pocket, row, col, vg)
                score[m] += length(pocket)
            else
                score[state[idx]] += 1
            end
        end
    end
    return score
end

function user_input()
    score = calculate_score()
    println(stone[black], " :", score[black], ". ", stone[white], " :", score[white])
    print("place ", stone[move], " >> ")
    user = lowercase(strip(readline()))
    if user == "end"
        return false
    elseif user == "pass"
        global passes += 1
        if passes==2
            return false
        end
        global move = move==black ? white : black
    elseif length(user) == 2
        if occursin(user[1], "abcdefghi") && occursin(user[2], "123456789")
            process_input(user)
        elseif occursin(user[2], "abcdefghi") && occursin(user[1], "123456789")
            process_input(user[2]*user[1])
        end
    end

    return true
end

if !TESTING
    println("GO")
    new_game()
    display()
    while user_input()
        display()
    end
    score = calculate_score()
    if score[black]>score[white]
        println("Blacks wins!")
    elseif score[white]>score[black]
        println("White win!")
    else
        println("Draw!")
    end
end

max_stars_repo_name: aole/Games-in-Julia

max_issues_repo_name: aole/Games-in-Julia

max_forks_repo_name: aole/Games-in-Julia


text:
using HorizonSideRobots
include("BackPath.jl")
include("Structs.jl")

function paint_perimeter(robot::Robot)
    back_path = BackPath(robot)
    r = PutmarkersRobot(robot)
    sides = [Nord, Ost, Sud, West]
    for side in sides
        movements!(r, side)
    end
    back!(robot, back_path)
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

@doc doc"""
    abstract type AbstractPosteriorDensity <: AbstractDensity end

Abstract super-type for posterior probability densities.
"""
abstract type AbstractPosteriorDensity <: AbstractDensity end
export AbstractPosteriorDensity

@doc doc"""
    getlikelihood(posterior::AbstractPosteriorDensity)::AbstractDensity

*BAT-internal, not part of stable public API.*

The likelihood density of `posterior`. The likelihood may or may not be
normalized.
"""
function getlikelihood end

@doc doc"""
    getprior(posterior::AbstractPosteriorDensity)::AbstractDensity

*BAT-internal, not part of stable public API.*

The prior density of `posterior`. The prior may or may not be normalized.
"""
function getprior end

function eval_logval(
    density::AbstractPosteriorDensity,
    v::Any,
    T::Type{<:Real} = density_logval_type(v);
    use_bounds::Bool = true,
    strict::Bool = false
)
    v_shaped = reshape_variate(varshape(density), v)
    if use_bounds && !variate_is_inbounds(density, v_shaped, strict)
        return log_zero_density(T)
    end

    prior_logval = eval_logval(
        getprior(density), v_shaped,
        use_bounds = false, strict = false
    )

    # Don't evaluate likelihood if prior probability is zero. Prevents
    # failures when algorithms try to explore parameter space outside of
    # definition of likelihood (as long as prior is chosen correctly).
    if !is_log_zero(prior_logval, T)
        likelihood_logval = eval_logval(
            getlikelihood(density), v_shaped,
            use_bounds = false, strict = false
        )
        convert(T, likelihood_logval + prior_logval)
    else
        log_zero_density(T)
    end
end

function eval_logval_unchecked(density::AbstractPosteriorDensity, v::Any)
    eval_logval(
        density, v,
        use_bounds = false, strict = false
    )
end

function var_bounds(density::AbstractPosteriorDensity)
    li_bounds = var_bounds(getlikelihood(density))
    pr_bounds = var_bounds(getprior(density))
    if ismissing(li_bounds)
        pr_bounds
    else
        li_bounds ∩ pr_bounds
    end
end

@doc doc"""
    PosteriorDensity{
        Li<:AbstractDensity,
        Pr<:DistLikeDensity,
        ...
    } <: AbstractPosteriorDensity

A representation of a PosteriorDensity, based a likelihood and prior, a
representation of the posterior density is cached internally. The densities
be accessed via

```julia
getlikelihood(posterior::PosteriorDensity)::Li
getprior(posterior::PosteriorDensity)::Pr
```

Constructors:

```julia
PosteriorDensity(likelihood::AbstractDensity, prior::AbstractDensity)
PosteriorDensity(likelihood::Any, prior::Any)
PosteriorDensity(log_likelihood::Function, prior::Any)
```
"""
struct PosteriorDensity{
    L<:AbstractDensity,
    P<:AbstractDensity,
    B<:AbstractVarBounds,
    S<:AbstractValueShape
} <: AbstractPosteriorDensity
    likelihood::L
    prior::P
    parbounds::B
    parshapes::S
end

export PosteriorDensity

function PosteriorDensity(likelihood::Any, prior::Any)
    li = convert(AbstractDensity, likelihood)
    pr = convert(AbstractDensity, prior)

    parbounds = _posterior_parbounds(
        var_bounds(li),
        var_bounds(pr)
    )

    parshapes = _posterior_parshapes(
        varshape(li),
        varshape(pr)
    )

    PosteriorDensity(li, pr, parbounds, parshapes)
end

getlikelihood(posterior::PosteriorDensity) = posterior.likelihood

getprior(posterior::PosteriorDensity) = posterior.prior

var_bounds(posterior::PosteriorDensity) = posterior.parbounds

ValueShapes.varshape(posterior::PosteriorDensity) = posterior.parshapes

function _posterior_parshapes(li_ps::AbstractValueShape, pr_ps::AbstractValueShape)
    if li_ps == pr_ps
        li_ps
    else
        throw(ArgumentError("Variable shapes of likelihood and prior are incompatible"))
    end
end

function _posterior_parshapes(li_ps::AbstractValueShape, pr_ps::ArrayShape{T,1}) where T
    n = totalndof(li_ps)
    if n == totalndof(pr_ps)
        @assert size(pr_ps) == (n,)
        li_ps
    else
        throw(ArgumentError("Likelihood and prior have different number of free parameters"))
    end
end

_posterior_parshapes(li_ps::ArrayShape{T,1}, pr_ps::AbstractValueShape) where T =
    _posterior_parshapes(pr_ps, li_ps)

function _posterior_parshapes(li_ps::ArrayShape{T,1}, pr_ps::ArrayShape{U,1}) where {T,U}
    if size(li_ps) == size(pr_ps)
        li_ps
    else
        throw(ArgumentError("Likelihood and prior have different number of free parameters"))
    end
end

_posterior_parshapes(li_ps::Missing, pr_ps::AbstractValueShape) = pr_ps

_posterior_parshapes(li_ps::Union{AbstractValueShape,Missing}, pr_ps::Missing) =
    throw(ArgumentError("Parameter shapes of prior must not be missing"))

_posterior_parbounds(li_bounds::AbstractVarBounds, pr_bounds::AbstractVarBounds) =
     li_bounds ∩ pr_bounds

_posterior_parbounds(li_bounds::Missing, pr_bounds::AbstractVarBounds) = pr_bounds


text:
export SourceTime

"""
    SourceTime()

QuakeML-compliant seismic source-time parameterization.

Field     | Type    | Meaning         | SeisIO conventions/behavior
--------: |:------- |:--------------  | :----------
 desc     | String  | description     |
 dur      | Float64 | duration        |
 rise     | Float64 | rise time       |
 decay    | Float64 | decay time      |

"""
mutable struct SourceTime
  desc  ::String    # description
  dur   ::Real      # duration
  rise  ::Real      # rise time
  decay ::Real      # decay time

  function SourceTime(
                      desc  ::String,    # description
                      dur   ::Real,      # duration
                      rise  ::Real,      # rise time
                      decay ::Real,      # decay time
                    )
    return new(desc, dur, rise, decay)
  end
end

SourceTime(;
          desc    ::String            = "",
          dur     ::Real              = zero(Float64),
          rise    ::Real              = zero(Float64),
          decay   ::Real              = zero(Float64),
          ) = SourceTime(desc, dur, rise, decay)

isempty(ST::SourceTime) = min(getfield(ST, :dur) == zero(Float64),
                              getfield(ST, :rise) == zero(Float64),
                              getfield(ST, :decay) == zero(Float64),
                              isempty(getfield(ST, :desc)))

function hash(ST::SourceTime)
  h = hash(zero(UInt64))
  for f in fieldnames(SourceTime)
    h = hash(getfield(ST, f), h)
  end
  return h
end

function isequal(S::SourceTime, U::SourceTime)
  q::Bool = isequal(getfield(S, :desc), getfield(U, :desc))
  for f in (:dur, :rise, :decay)
    q = min(q, getfield(S,f) == getfield(U,f))
  end
  return q
end
==(S::SourceTime, U::SourceTime) = isequal(S, U)

sizeof(ST::SourceTime) = 56 + sizeof(getfield(ST, :desc))

function write(io::IO, ST::SourceTime)
  desc = getfield(ST, :desc)
  write(io, Int64(sizeof(desc)))
  write(io, desc)
  for f in (:dur, :rise, :decay)
    write(io, getfield(ST, f))
  end
  return nothing
end

function read(io::IO, ::Type{SourceTime})
  ST = SourceTime()
  L = fastread(io, Int64)
  setfield!(ST, :desc, String(fastread(io, L)))
  for f in (:dur, :rise, :decay)
    setfield!(ST, f, fastread(io, Float64))
  end
  return ST
end

function show(io::IO, ST::SourceTime)
  if get(io, :compact, false) == false
    println(io, "SourceTime with fields:")
    for f in (:desc, :dur, :rise, :decay)
      fn = lpad(String(f), 5, " ")
      println(io, fn, ": ", getfield(ST, f))
    end
  else
    c = :compact => true
    st_str = string("dur ", repr(getfield(ST, :dur), context=c),
                    ", rise ", repr(getfield(ST, :rise), context=c),
                    ", decay ", repr(getfield(ST, :decay), context=c))

    print(io, str_trunc(st_str, max(80, displaysize(io)[2]-2) - show_os))
  end
  return nothing
end