text:
using ImageCore, Colors, FixedPointNumbers, OffsetArrays, Test

if VERSION >= v"1.2.0-DEV.229"
    sumsz(img) = Base.dims2string(size(img)) * ' '
else
    sumsz(img) = ""
end

@testset "show" begin
    thismodule = string(@__MODULE__)
    if thismodule != "Main"
        prefixF = "FixedPointNumbers."
        prefixC = "ColorTypes."
    else
        prefixF = prefixC = ""
    end
    rgb32 = rand(RGB{Float32}, 3, 5)
    v = view(rgb32, 2:3, :)
    @test summary(v) == "2×5 view(::Array{RGB{Float32},2}, 2:3, :) with eltype $(prefixC)RGB{Float32}"
    a = channelview(rgb32)
    @test summary(a) == "3×3×5 reinterpret(Float32, ::Array{RGB{Float32},3})"
    num64 = rand(3,5)
    b = colorview(RGB, num64)
    @test summary(b) == "5-element reshape(reinterpret(RGB{Float64}, ::Array{Float64,2}), 5) with eltype $(prefixC)RGB{Float64}"
    rgb8 = rand(RGB{N0f8}, 3, 5)
    c = rawview(channelview(rgb8))
    @test summary(c) == "3×3×5 rawview(reinterpret(N0f8, ::Array{RGB{N0f8},3})) with eltype UInt8"
    @test summary(rgb8) == "3×5 Array{RGB{N0f8},2} with eltype $(prefixC)RGB{$(prefixF)Normed{UInt8,8}}"
    rand8 = rand(UInt8, 3, 5)
    d = normedview(permuteddimsview(rand8, (2,1)))
    @test summary(d) == "5×3 normedview(N0f8, PermutedDimsArray(::Array{UInt8,2}, (2, 1))) with eltype $(prefixF)Normed{UInt8,8}"
    e = permuteddimsview(normedview(rand8), (2,1))
    @test summary(e) == "5×3 PermutedDimsArray(reinterpret(N0f8, ::Array{UInt8,2}), (2, 1)) with eltype $(prefixF)Normed{UInt8,8}"
    f = permuteddimsview(normedview(N0f16, rand(UInt16, 3, 5)), (2,1))
    @test summary(f) == "5×3 PermutedDimsArray(reinterpret(N0f16, ::Array{UInt16,2}), (2, 1)) with eltype $(prefixF)Normed{UInt16,16}"
    g = channelview(rgb8)
    @test summary(g) == "3×3×5 reinterpret(N0f8, ::Array{RGB{N0f8},3})"
    h = OffsetArray(rgb8, -1:1, -2:2)
    @test summary(h) == "$(sumsz(h))OffsetArray(::Array{RGB{N0f8},2}, -1:1, -2:2) with eltype $(prefixC)RGB{$(prefixF)Normed{UInt8,8}} with indices -1:1×-2:2"
    i = channelview(h)
    @test summary(i) == "$(sumsz(i))reinterpret(N0f8, OffsetArray(::Array{RGB{N0f8},3}, 1:1, -1:1, -2:2)) with indices 1:3×-1:1×-2:2"
    c = channelview(rand(RGB{N0f8}, 2))
    o = OffsetArray(c, -1:1, 0:1)
    @test summary(o) == "$(sumsz(o))OffsetArray(reinterpret(N0f8, ::Array{RGB{N0f8},2}), -1:1, 0:1) with eltype $(prefixF)Normed{UInt8,8} with indices -1:1×0:1"
    # Issue #45
    a = collect(tuple())
    @test summary(a) == "0-element Array{Union{},1}"
    b = view(a, :)
    @test summary(b) == "0-element view(::Array{Union{},1}, :) with eltype Union{}"
end

nothing


text:
"""
`scannerlist(lst...)`

Define scanner list. Scanner list is used by the `SD1` command (see [`SD1cmd`](@ref) and [`SD1`](@ref)) to define which scanners are available to the user.

The arguments of the function are tuples composed of the following parameters:

 * `scn` Integer or Integer range specifying the scanners used
 * `npp` Number of pressure ports in the scanner
 * `lrn` Logical Range number. I don't really understand this. Just use `1`

# Examples

```julia-repl
julia> scannerlist((1:2, 64, 1), (3:4, 32, 1), (5,16,2))
5-element Vector{Tuple{Int64, Int64, Int64}}:
 (1, 64, 1)
 (2, 64, 1)
 (3, 32, 1)
 (4, 32, 1)
 (5, 16, 2)
```

"""
function scannerlist(lst...; npp=64, lrn=1)

    nscanners = 0
    scnlst = Tuple{Int,Int,Int}[]

    for s in lst
        if isa(s, Tuple)
            scn = s[1]
            if length(s) >= 3
                npp1 = Int(s[2])
                lrn1 = Int(s[3])
            elseif length(s) == 2
                npp1 = Int(s[2])
                lrn1 = lrn
            elseif length(s) == 1
                npp1 = npp
                lrn1 = lrn
            end
        else
            scn = s
            npp1 = npp
            lrn1 = lrn
        end
        nscanners += length(scn)
        if nscanners > 8
            throw(BoundsError(nscanners, "Maximum of 8 scanners is possible"))
        end
        for scanner in scn
            if scanner < 1 || scanner > 8
                throw(BoundsError(scanner, "Only scanners 1-8 are possible"))
            end
            if npp1 ∉ (16,32,64)
                throw(BoundsError(npp, "Only ESP with 16, 32 or 64 possible"))
            end
            push!(scnlst, (scanner, npp1, lrn1))
        end
    end

    return scnlst
end

"""
`availablechans(scnlst)`
`availablechans(scnlst...)`

Return the number of available channels in scanner list.

"""
availablechans(scnlst) = sum(s[2] for s in scnlst)
availablechans(scnlst...) = sum(s[2] for s in scnlst...)

"""
`SD1cmd(scnlst; crs="111")`

Implements the `SD1` command which configures the scanners connected to the DTC Initium.
This function actually creates a `String` containing the command to be sent to the DTC Inicit. Usually, to actually send the command, this function is called by the [`SD1`](@ref) function.

The `scnlst` parameter should be created with the [`scannerlist`](@ref)  function.

# Examples

```julia-repl
julia> scnlst = scannerlist((1:8,64,1))  # Use eight scanners with 64 pressure ports
8-element Vector{Tuple{Int64, Int64, Int64}}:
 (1, 64, 1)
 (2, 64, 1)
 (3, 64, 1)
 (4, 64, 1)
 (5, 64, 1)
 (6, 64, 1)
 (7, 64, 1)
 (8, 64, 1)

julia> SD1cmd(scnlst)
"SD1 111 (1,64,1) (2,64,1) (3,64,1) (4,64,1) (5,64,1) (6,64,1) (7,64,1) (8,64,1);"

```
"""
function SD1cmd(scnlst; crs="111")

    cmd = "SD1 $crs"

    for (scn,npp,lrn) in scnlst
        cmd *= " ($scn,$npp,$lrn)"
    end
    cmd *= ";"
    return cmd
end

"""
`setdaqparams(;stbl=1, nfr=64, nms=1, msd=100, trm=0, scm=1, ocf=2)`

Creates a dictionary with data acquisition setup parameters to be used by the [`SD2`](@ref) command that configures data acquisition. This function creates arguments that can be used the the `SD2cmd` function.

Parameters
 * `stbl` Chooses a DA Setup table. Should be an integer 1-5
 * `nfr` Integer that specifies the number of samples that should be averaged to create a measurement set.
 * `nms` Default number of pressure measurement sets that should be acquired. 0 if continuous data acquisition is desired.
 * `msd`  Delay between each pressure set. If 0, the scanners will acquire pressure as fast as possible
 * `trm` Trigger mode.
     - 0 no trigger
     - 1 Initial measurement is triggered
     - 2 Each pressure measurement is triggered
 * `scm` Defines scam modes. Just use 1. See the manual for further information
 * `ocf` Output format for the measurement. Just use 2 in general
     - 1 Raw measurements
     - 2 Receive pressure measurements in EU units
     - 3 EU press-sets and ez-set and temp-sets

# Example

```julia-repl
julia> setdaqparams()
Dict{Symbol, Int64} with 9 entries:
  :frd   => 0
  :scm   => 1
  :ocf   => 2
  :nfrez => 64
  :nms   => 1
  :trm   => 0
  :nfr   => 64
  :msd   => 100
  :stbl  => 1

```

"""
function setdaqparams(;stbl=1, nfr=64, nms=1, msd=100, trm=0, scm=1, ocf=2)
    if !(1 ≤ stbl ≤ 5)
        throw(BoundsError(stbl, "Scan table limited to 1-5!"))
    end

    if length(nfr) == 2
        nfrez = nfr[2]
        nfr = nfr[1]
    else
        nfrez = nfr
    end

    if !(1 ≤ nfr ≤ 127)
        throw(BoundsError(nfr, "nfr should range from 1-127"))
    end

    if !(1 ≤ nfrez ≤ 127)
        throw(BoundsError(nfr, "nfrez should range from 1-127"))
    end

    if !(0 ≤ nms ≤ 65000)
        throw(BoundsError(nfr, "nms should range from 0-65000"))
    end

    if !(0 ≤ msd ≤ 600_000)
        throw(BoundsError(nfr, "msd should range from 0-600_000"))
    end

    if !(0 ≤ trm ≤ 2)
        throw(BoundsError(nfr, "trm should be 0, 1 or 2"))
    end
    if !(0 ≤ scm ≤ 1)
        throw(BoundsError(nfr, "scm should be 0 or 1"))
    end
    if !(1 ≤ ocf ≤ 3)
        throw(BoundsError(ocf, "trm should be 1, 2 or 3"))
    end

    frd = 0  # Unused

    return Dict(:stbl=>stbl, :nfr=>nfr, :nfrez=>nfrez, :frd=>frd, :nms=>nms,
                :msd=>msd, :trm=>trm, :scm=>scm, :ocf=>ocf)

end

"""
`SD2cmd(p; crs="111")`

This function will actually create the command `String` that should be sent to the DTC Initium. The parameter `p` is usually created by function [`setdaqparams`](@ref). The user will usually use function [`SD2`](@ref) that actually communicates with the DTC Initium.

# Example

```julia-repl
julia> p = setdaqparams()
Dict{Symbol, Int64} with 9 entries:
  :frd   => 0
  :scm   => 1
  :ocf   => 2
  :nfrez => 64
  :nms   => 1
  :trm   => 0
  :nfr   => 64
  :msd   => 100
  :stbl  => 1

julia> SD2cmd(p)
"SD2 111 1 (64-64 0) (1 100) (0 1) 2;"

```
"""
function SD2cmd(p; crs="111")

    #cmd = "SD2 $crs $(p[:stbl]) ($(p[:nfr])-$(p[:nfrez]) $(p[:frd])) ($(p[:nms]) $(p[:msd])) ($(p[:trm]) $(p[:scm])) $(p[:ocf]);"
    cmd = "SD2 $crs $(p[:stbl]) ($(p[:nfr]) $(p[:frd])) ($(p[:nms]) $(p[:msd])) ($(p[:trm]) $(p[:scm])) $(p[:ocf]);"

end

"""
`SD3cmd(stbl, plst; crs="111")`

Creates the String that can be sent to the DTC Initium with information specifying
which pressure ports should be acquired.

The syntax of the port list `plst` is described in the user manual but there are auxiliary functions that help with this task. See [`PortRange`](@ref), [`strport`](@ref), [`portlist`](@ref), [`strportlist`](@ref) for further information.

Usually, this function is called by function [`SD3`](@ref) and [`addinputs`](@ref).
"""
SD3cmd(stbl, plst; crs="111") = "SD3 $crs $stbl,$plst;"

"""
`SD5cmd(stbl, actx; crs="111")`

Implements the `SD5` command. This command does a lot of things but for now, only the control option is implemented (`stbl`=-1).

The default mode of operation uses `actx`=1. In this case, on each measurement, both the pressure and the temperature are measured simultaneously. If high speed is necessary, the temperature can be measured only once (at the beginning) and this is done setting `actx=0`. Positive values greater than 1 of `actx` specify the intermittent rate at which
temp-sets should be acquired.

Again, this function only creates the String of the command. This function is usually called by [`SD5`](@ref).

"""
function SD5cmd(stbl, actx; crs="111")

    if stbl != -1
        throw(DomainError(stbl, "Control form of SD5 implemented only: stbl = -1!"))
    end
    if actx < 0
        throw(DomainError(actx, "Non-negative integer is acceptable only!"))
    end

    return "SD5 $crs -1 $actx;"
end

"""
`PC4cmd(unx, fct=0; lrn=1)`

Create command to change pressure units. There are predefined units for `unx` in 1-12. If `unx` is 0 or 13, a unit conversion factor from psi is employed.

This function is usually called by [`PC4`](@ref).
"""
function PC4cmd(unx, fct=0; lrn=1)

    if unx==0 || unx==13
        if fct ≤ 0
            throw(DomainError(unx, "For specifying unit conversion factors, unx should be either 0 or 13! and fct should be positive"))
        end
    elseif 1 ≤ unx ≤ 12
        fct = 0
    else
        throw(DomainError(unx, "unx should be 0-13!"))
    end
    return "PC4 $lrn $unx $fct;"
end

"""
`CV1cmd(valpos, puldur)`

Create String command to set calibration valve position for all scanners.
"""
function CV1cmd(valpos, puldur)

    if valpos != 0 || valpos != 1
        throw(DomainError(valpos, "Valve position should be either 0 (RUNPOS) or 1 (CALPOS)"))
    end

    if !(0 ≤ puldur ≤ 199)
        throw(DomainError(puldur, "Pulse duration should be 0-199"))
    end

    return "CV1 $valpos, $puldur;"
end

"""
`CP1cmd(puldur)`

Create String command to set pneumatic pressure calibration valve mode.
"""
function CP1cmd(puldur)
    if puldur < 0 || puldur > 30
        throw(DomainError(puldur, "PULSE duration should be between 0 and 30"))
    end
    return "CP1 $puldur;"
end

"""
`CP2cmd(stbtim)`

Create String command to set the delay time the Initium will wait for pressure to stabilize before reading any data.
"""
function CP2cmd(stbtim)
    if stbtim < 1 || stbtim > 199
        throw(DomainError(stbtim, "Calibration stabilization time should be 1-199 seconds!"))
    end
    return "CP2 $stbtim;"
end

"""
`CA2cmd(lrn=1)`

Creates command String to perform zero pressure calibration.
"""
CA2cmd(lrn=1) = "CA2 $lrn;"

"Create OP2 - output DA setup table's coefficients - command"
OP2cmd(stbl, ports; crs="111") = "OP2 $crs $(-stbl) $ports;"

"Create OP3 - output DA setup table's coefficients (all) - command"
OP3cmd(stbl, ports; crs="111") = "OP3 $crs $stbl $ports;"

"Create OP5 - output scan list of setup table - command"
OP5cmd(stbl; crs="111") = "OP5 $crs $stbl;"

"""
`AD0cmd()`

Creates String command that stops high-speed data acquisition. The command [`AD0`](@ref) is usually used since it actually communicates with the DTC Initium.
"""
AD0cmd() = "AD0;"

"""
`AD2cmd(stbl, nms)`

Creates command to start data acquisition run.

"""
function AD2cmd(stbl, nms)

    if !(1 ≤ stbl ≤ 5)
        throw(DomainError(stbl, "stbl should be 1-5"))
    end

    if !(0 ≤ nms ≤ 65_000)
        throw(DomainError(nms, "nms should be 0-65000"))
    end

    return "AD2 $stbl $nms;"
end

function AD2cmd(stbl)
    if !(1 ≤ stbl ≤ 5)
        throw(DomainError(stbl, "stbl should be 1-5"))
    end
    return "AD2 $stbl;"

end

"Create LA1 - Look at scanner status - command"
LA1cmd(port; crs="111") = "LA1 $crs $port;"

"Create LA4 - Look at Firmware version - command"
LA4cmd(;crs="111") = "LA4 $crs;"

"""
`sendcommand!(io, cmd, buf, nbytes)`

Sends a command to socket.
"""
function sendcommand!(io, cmd, buf, nbytes)
    println(io, cmd)
    readbytes!(io, buf, nbytes)
end

ext: jl


text:
using LightXML
using JuMP, GLPK
using ElasticArrays
using DataStructures
using OffsetArrays

global counter = 0
global counter1 = 0
global counter2 = 0
const EPS = 1e-9
const MAX_INT = 9223372036854775807
const DIGITS = 9

include("data_structure.jl")
include("utils.jl")
include("solomon.jl")
include("ESPPRC.jl")
include("pulse.jl")
include("labeling.jl")
include("relaxed_master_problem.jl")
include("branch_and_bound.jl")


text:
struct TextXMLElement <: MyXMLElement
    el::XMLElement
end

function name(el::TextXMLElement)
    return name(el.el)
end

function make_xml(::TextXMLElement)
    # do nothing
end

function get_loggables(el::TextXMLElement)
    return [el.el]
end

function get_priors(el::TextXMLElement)
    return [el]
end

function add_ref_el(el::XMLElement, txt::TextXMLElement)
    add_child(el, txt.el)
end

function add_ref_els(el::XMLElement, txts::Vector{TextXMLElement})
    for txt in txts
        add_ref_el(el, txt)
    end
end

function get_id(el::TextXMLElement)
    id = attribute(el.el, bn.ID)
    if isnothing(id)
        id = attribute(el.el, bn.IDREF)
    end
    if isnothing(id)
        error("no id or idref")
    end

    return id
end

################################################################################
## Build from existing xml
################################################################################

function BEASTXMLElement(xml::String)
    xdoc = parse_file(xml)
    xroot = root(xdoc)
    if name(xroot) != bn.BEAST
        error("Not a BEAST xml.")
    end

    bx = BEASTXMLElement()
    taxa_parsed = false
    for element in child_elements(xroot)
        el = parse_element(element, taxa_parsed = taxa_parsed)
        if typeof(el) <: EmptyDataXMLElement
            taxa_parsed = true
        end
        if !isnothing(el)
            add_child(bx, el)
        end
    end
    return bx
end

# const SPECIAL_PARSERS = Dict(bn.TAXA => parse_taxa)

const DO_NOT_PARSE = ["coalescentSimulator", bn.TREE_MODEL, "rescaledTree", bn.NEWICK, "newick"]

function parse_element(el::XMLElement; taxa_parsed = false)
    nm = name(el)
    if nm == bn.TAXA && !taxa_parsed
        return parse_taxa(el)
    elseif nm in DO_NOT_PARSE
        return nothing
    elseif nm == bn.OPERATORS
        return parse_operators(el)
    elseif nm == bn.MCMC
        return parse_mcmc(el)
    end
    return TextXMLElement(el)
end

struct EmptyDataXMLElement <: MyXMLElement
    taxa::Vector{String}
end

function parse_taxa(el::XMLElement)
    taxa_els = el[bn.TAXON]
    taxa = [get_id(x) for x in taxa_els]
    return EmptyDataXMLElement(taxa)
end

function parse_operators(el::XMLElement)
    ops = TextXMLElement[]
    for el in child_elements(el)
        push!(ops, TextXMLElement(el))
    end
    return OperatorsXMLElement(ops)
end

struct ParsedMCMCXMLElement <: MyXMLElement
    priors::Vector{TextXMLElement}
    likelihoods::Vector{TextXMLElement}
    file_log::TextXMLElement
    tree_log::TextXMLElement
    file_loggables::Vector{TextXMLElement}
end

function parse_mcmc(el::XMLElement)
    j_el = find_element(el, "joint")
    prior_el = find_element(j_el, bn.PRIOR)
    priors = [TextXMLElement(cel) for cel in child_elements(prior_el)]

    like_el = find_element(j_el, bn.LIKELIHOOD)
    likes = [TextXMLElement(cel) for cel in collect(child_elements(like_el))]

    file_log = nothing
    log_els = el[bn.LOG]
    file_loggables = nothing
    for lg in log_els
        id = get_id(lg)
        if id == "fileLog"
            set_attribute(lg, bn.ID, "fileLogSequence")
            file_log = TextXMLElement(lg)
            c_els = child_elements(lg)
            file_loggables = TextXMLElement[]
            ind = 1
            for el in c_els
                if ind > 3 #first thre are alwasy posterior, likelihood, prior
                    push!(file_loggables, TextXMLElement(el))
                end
                ind += 1
            end
        end
    end

    tree_log = TextXMLElement(find_element(el, "logTree"))
    j_el = find_element(tree_log.el, "joint")
    unlink(j_el)
    free(j_el)
    post_el = new_child(tree_log.el, bn.POSTERIOR)
    set_attribute(post_el, bn.IDREF, bn.POSTERIOR)

    return ParsedMCMCXMLElement(priors, likes, file_log, tree_log, file_loggables)
end

################################################################################
## Merge xml
################################################################################

# This will only work with some very specific circumnstances
#   1. Both have some el <: AbstractDataXMLElement

function findfirst_name(bx::BEASTXMLElement, nm::String)
    return findfirst(x -> name(x) == nm, bx.components)
end

function findall_name(bx::BEASTXMLElement, nm::String)
    return findall(x -> name(x) == nm, bx.components)
end

function merge_xml!(traitxml::BEASTXMLElement, seqxml::BEASTXMLElement; separate_logs::Bool = true)

    # merge taxa
    tx_trait = find_element(traitxml, DataXMLElement)
    # @show tx_trait
    tx_seq = find_element(seqxml, EmptyDataXMLElement)
    taxa_trait = tx_trait.taxa
    taxa_seq = tx_seq.taxa
    if Set(taxa_trait) != Set(taxa_seq)
        only_trait = setdiff(taxa_trait, taxa_seq)
        only_seq = setdiff(taxa_seq, taxa_trait)
        if length(only_seq) > 0
            error("The following taxa are present in the sequence xml but not the data:\n\t" *
                  join(only_seq, ' '))
        else
            error("The following taxa are present in the data but not the sequence xml:\n\t" *
                  join(only_trait, ' '))
        end
    end

    seq_start = findfirst(x -> typeof(x) <: EmptyDataXMLElement, seqxml.components) + 1

    seq_ops_ind = findfirst(x -> typeof(x) <: OperatorsXMLElement, seqxml.components)
    seq_stop = seq_ops_ind - 1

    trait_ind = findfirst(x -> typeof(x) <: TreeModelXMLElement, traitxml.components) + 1

    for i = seq_start:seq_stop
        add_child(traitxml, seqxml.components[i], trait_ind)
        trait_ind += 1
    end

    trait_ops_ind = findfirst(x -> typeof(x) <: OperatorsXMLElement, traitxml.components)
    merge_operators!(traitxml.components[trait_ops_ind],
                     seqxml.components[seq_ops_ind])

    seqmc = find_element(seqxml, ParsedMCMCXMLElement)
    traitmc = find_element(traitxml, MCMCXMLElement)
    merge_mcmc!(traitmc, seqmc, separate_logs = separate_logs)

    # display(traitxml)
    # add_loggables(traitxml, LoggablesXMLElement(seqmc.file_loggables))
    # display(traitxml)
    # error()

    newick_el = find_element(traitxml, NewickXMLElement)
    newick_el.fix_tree = false
end


text:
# Only run coverage from linux release build on travis.
get(ENV, "TRAVIS_OS_NAME", "")       == "linux" || exit()
get(ENV, "TRAVIS_JULIA_VERSION", "") == "0.6"   || exit()

using Coverage

Codecov.submit(Codecov.process_folder())


text:
struct System
    criteria        :: Vector{String}
    levels          :: Vector{Vector{String}}
    points          :: Dict{Pair{String,String},Int}
end
System(criteria, levels) = System(criteria, levels, Dict())
System() = System([], [])

criteria(s::System) = s.criteria
levels(s::System) = s.levels

shape(s::System) = length.(levels(s))

function addcriterion(s::System, c)
    push!(criteria(s), c)
    push!(levels(s), [])
end

addlevel(s::System, l) = push!(last(levels(s)), l)

setpts(s::System, key, pts) = (s.points[key] = pts)
getpts(s::System, key) = get(s.points, key, nothing)

_key(s, i, j) = criteria(s)[i] => levels(s)[i][j]

setpts(s::System, i, j, pts) = setpts(s, _key(s, i, j), pts)
getpts(s::System, i, j) = getpts(s, _key(s, i, j))


text:
using Documenter
using RegisterCore

makedocs(
    sitename = "RegisterCore",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    modules = [RegisterCore],
    pages = ["index.md", "api.md"]
)

deploydocs(
    repo = "github.com/HolyLab/RegisterCore.jl.git"
)


text:
using Statistics
using BlackBoxOptim
using NLopt
using Distributions
using StatsBase

## Objective functions year

function objfun(g...; m=mchain_op)
    # @unpack nsim, nfast = Par
    gamma = g
    N, nfast = size(m)

    logunif=log.(rand(nfast))
    geta=m[:,1]
    hM=zeros(N)
    h=zeros(N,nfast)
    # ar=0
    for r=1:nfast
        val=0
        gtry=m[:,r]
        for j=1:N
            val=+gtry[j].*gamma[j]-geta[j].*gamma[j]
        end
        geta=logunif[r] < val[1] ? gtry : geta
        # ar=+ logunif[r] < val[1] ? 1 : 0
        hM=+geta
        h[:,r]=geta
    end
    # acc_r=round((ar/nfast)*100, digits=2)
    # println("Acceptance rate is : ",acc_r)
    # hM=hM./nfast
    # var=(h.-hM)*(h.-hM)'
    #
    # lambda, QM = eigen(var)
    # sel = lambda.>0
    # An = QM[:,sel]
    # dvs = An'*(h.-hM)
    # vs = An'*var*An
    # omega = inv(vs)
    # obj=dvs'*omega*dvs
    #
    # # obj=hM'*inv(var)*hM
    #
    # return obj[1]
    hM=hM./nfast
    # var=(h.-hM)*(h.-hM)'
    bn=Integer(round(sqrt(N)))
    bt=Integer(N-bn+1)
    d=zeros(N,bt)
    for b=1:bt
        ind=b:bn+b-1
        hMb=sum(h[:,ind], dims=2)./bn
        d[:,b] = (hMb.-hM)
    end
    dv= d*d'
    vhat=(bn/(bt-1)).*dv
    lambda, QM = eigen(vhat)
    sel = lambda.>0
    An = QM[:,sel]
    dvs = An'*(hM)
    vs = An'*vhat*An
    omega = inv(vs)
    t=dvs'*omega*dvs
    t[1]
end

function objfun_bbo(gamma0; m=mchain_op)
    # @unpack nfast = Par
    gamma = gamma0
    N, nfast = size(m)
    # nfast = size(m,1)

    logunif=log.(rand(nfast))
    geta=m[:,1]
    hM=zeros(N)
    h=zeros(N,nfast)
    # ar=0
    for r=1:nfast
        val=0
        gtry=m[:,r]
        for j=1:N
            val=+gtry[j].*gamma[j]-geta[j].*gamma[j]
        end
        geta=logunif[r] < val[1] ? gtry : geta
        # ar=+ logunif[r] < val[1] ? 1 : 0
        h[:,r]=geta
        hM=+geta
    end
    # acc_r=round((ar/nfast)*100, digits=2)
    # println("Acceptance rate is : ",acc_r)
    # hM=hM./nfast
    # var=(h.-hM)*(h.-hM)'
    #
    # lambda, QM = eigen(var)
    # sel = lambda.>0
    # An = QM[:,sel]
    # dvs = An'*(h.-hM)
    # vs = An'*var*An
    # omega = inv(vs)
    # obj=dvs'*omega*dvs
    #
    # # obj=hM'*inv(var)*hM
    #
    # return obj[1]

    hM=hM./nfast
    # var=(h.-hM)*(h.-hM)'
    bn=Integer(round(sqrt(N)))
    bt=Integer(N-bn+1)
    d=zeros(N,bt)
    for b=1:bt
        ind=b:bn+b-1
        hMb=sum(h[:,ind], dims=2)./bn
        d[:,b] = (hMb.-hM)
    end
    dv= d*d'
    vhat=(bn/(bt-1)).*dv
    lambda, QM = eigen(vhat)
    sel = lambda.>0
    An = QM[:,sel]
    dvs = An'*(hM)
    vs = An'*vhat*An
    omega = inv(vs)
    t=dvs'*omega*dvs
    t[1]
end
## Bootstrap testing

function my_dist(m, g; P::Par=Param)
    @unpack nboost = P
    # nboost = 20000
    N, nfast = size(m)
    gamma=g
    logunif=log.(rand(nfast))
    geta=m[:,1]
    hM=zeros(N)
    h=zeros(N,nfast)
    d=zeros(nboost)
    ar=0.0
    for r=1:nfast
        val=0.0
        gtry=m[:,r]
        for j=1:N
            val=+gtry[j].*gamma[j]-geta[j].*gamma[j]
        end
        geta=logunif[r] < val[1] ? gtry : geta

        logunif[r] < val[1] && (ar=+1.0)

        h[:,r]=geta
        hM=+geta
    end
    acc_r=100.0*(ar/nfast)
    # println("Acceptance rate is : ", acc_r)
    @show (acc_r)

    hM=hM./nfast
    var=(h.-hM)*(h.-hM)'

    for b=1:nboost
        ind=rand(1:nfast,nfast)
        hMb=sum(h[:,ind], dims=2)./nfast
        varb=(h[:,ind].-hMb)*(h[:,ind].-hMb)'
        lambda, QM = eigen(varb)
        sel = lambda.>0
        An = QM[:,sel]
        dvs = An'*(hMb.-hM)
        vs = An'*varb*An
        omega = inv(vs)
        t=dvs'*omega*dvs
        d[b]=nfast.*t[1]
    end
    # d=[nfast*(h[:,i]-hM)'*inv(var)*(h[:,i]-hM) for i=1:nfast]
    f=ecdf(d)
    f, d
end

function sim_err(m, g; P::Par=Param)
    @unpack nboost = P
    # nboost = 20000
    N, nfast = size(m)
    gamma=g
    logunif=log.(rand(nfast))
    geta=m[:,1]
    hM=zeros(N)
    h=zeros(N,nfast)
    ar=0.0
    for r=1:nfast
        val=0.0
        gtry=m[:,r]
        for j=1:N
            val=+gtry[j].*gamma[j]-geta[j].*gamma[j]
        end
        geta=logunif[r] < val[1] ? gtry : geta

        # logunif[r] < val[1] && (ar=+1.0)

        h[:,r]=geta
        hM=+geta
    end
    # acc_r=100.0*(ar/nfast)
    # # println("Acceptance rate is : ", acc_r)
    # @show (acc_r)

    hM=hM./nfast
    # var=(h.-hM)*(h.-hM)'
    bn=Integer(round(sqrt(N)))
    bt=Integer(N-bn+1)
    d=zeros(N,bt)
    for b=1:bt
        ind=b:bn+b-1
        hMb=sum(h[:,ind], dims=2)./bn
        d[:,b] = (hMb.-hM)
    end
    dv= d*d'
    vhat=(bn/(bt-1)).*dv
    lambda, QM = eigen(vhat)
    sel = lambda.>0
    An = QM[:,sel]
    dvs = An'*(hM)
    vs = An'*vhat*An
    omega = inv(vs)
    t=dvs'*omega*dvs
    t[1]
end

function acc_rate(m, g; P::Par=Param, f=0)
    # @unpack nboost = P
    # nboost = 20000
    d=ndims(m)
    if d==3
        T, N, nfast = size(m)
        geta=m[:,:,1]
    else
        T=1
        N, nfast = size(m)
        geta=m[:,1]
    end

    gamma=g
    unif=rand(nfast)

    # hM=zeros(N)
    # h=zeros(N,nfast)
    # d=zeros(nboost)
    ar=0.0
    for t = 1:T, r=1:nfast
        val=0.0
        d==3 ? gtry=m[t,:,r] : gtry=m[:,r]
        for j=1:N
            val=+gtry[j].*gamma[j]-geta[j].*gamma[j]
        end
        # geta=logunif[r] < val[1] ? gtry : geta
        if f==0
            # log.(unif[r]) < val[1] && (ar=+1.0)
            ars=log(unif[r]) < val[1] ? 1.0 : 0.0
        else
            ars=unif[r] < 1.4^(val[1]) ? 1.0 : 0.0
        end
        ar=+ars
        # ar=+logunif[r] < val[1] ? 1.0 : 0.0
        # h[:,r]=geta
        # hM=+geta
    end
    acc_r=100.0*(ar/nfast)
    # println("Acceptance rate is : ", acc_r)
    @show (acc_r)
    acc_r
end

## Testing function

function RP_Cournot_test(m ; Par::Par=Param)
    @unpack nfast = Par
    N, nsim = size(m)
    # Generate random index
    indfast=rand(1:nsim, nfast)
    indfast[1]=1 # Keeping the first element fixed.
    mchain_op=m[:,indfast]

    # First optimization BlackBox
    res = bboptimize(x->objfun_bbo(x, m=mchain_op); SearchRange = (-10e300,10e300), NumDimensions = N, TraceMode=:silent, MaxTime=400.0)
    ## Getting results
    # minbb = best_fitness(res)
    γ0 = best_candidate(res)

    # mycdf, dist = my_dist(mchain_op,γ0)
    # TS = minbb*nfast
    # p_val = mycdf(TS)
    println("First optimization from ", myid()," done.")
    # Second optimization NLopt
    fop=nothing
    fop=Model(NLopt.Optimizer)
    set_optimizer_attribute(fop, "algorithm", :LN_BOBYQA)
    set_optimizer_attribute(fop, "xtol_rel", 1e-6)
    # set_optimizer_attribute(fop, "maxtime", 300)
    @variable(fop, gamma0[1:N])

    register(fop, :objfun, N, (x...)->objfun(x..., m=mchain_op); autodiff=true )
    @NLobjective(fop, Min, objfun(gamma0...))

    for i=1:N
        set_start_value(gamma0[i],γ0[i])
    end

    ## Solving in NLopt
    JuMP.optimize!(fop)

    ## Getting Results
    # termination_status(fop)
    # primal_status(fop)

    γ1=value.(gamma0)
    # obj_v=objective_value(fop)
    # TS=obj_v*nfast
    # mycdf, dist = my_dist(mchain_op,γ1)
    # p_val = mycdf(TS)
    println("Second optimization from ", myid()," done.")
    ## Second optimization for refinement
    for i=1:N
        set_start_value(gamma0[i],γ1[i])
    end
    ## Solving in NLopt
    JuMP.optimize!(fop)

    ## Getting Results
    γ2=value.(gamma0)
    obj_v=objective_value(fop)
    TS=obj_v
    # mycdf, dist = my_dist(mchain_op,γ2)
    # p_val = 1-mycdf(TS)
    p_val = ccdf(Chisq(N),TS)
    println("Third optimization from ", myid()," done.")
    # fop=nothing
    acc_r = MyModules.acc_rate(mchain_op,γ2)
    return TS, p_val , acc_r
end

## Objective functions semester

# function objfun_s(gamma0...; mchain=mchain_op)
#     # @unpack nsim, nfast = Par
#     gamma = Vector(gamma0[1])
#     T,N,nfast = size(mchain)
#
#     logunif=log.(rand(T,nfast))
#     geta=mchain[:,:,1]
#     hM=zeros(T,N)
#     for t=1:T, r=1:nfast
#         val=0
#         gtry=mchain[t,:,r]
#         for j=1:N
#             val=+gtry[j].*gamma[j]-geta[t,j].*gamma[j]
#         end
#         geta[t,:]=logunif[t,r] < val[1] ? gtry : geta[t,:]
#         hM[t,:]=+geta[t,:]./nfast
#     end
#     mhM=sum(hM, dims=1)./T
#     omega=(hM.-mhM)'*(hM.-mhM)
#
#     obj=mhM*inv(omega)*mhM'
#
#     return obj[1]
# end
#
# function objfun_s_bbo(gamma0; mchain=mchain_op)
#     # @unpack nfast = Par
#     gamma = gamma0
#     T,N,nfast = size(mchain)
#
#     logunif=log.(rand(T,nfast))
#     geta=mchain[:,:,1]
#     hM=zeros(T,N)
#     for t=1:T, r=1:nfast
#         val=0
#         gtry=mchain[t,:,r]
#         for j=1:N
#             val=+gtry[j].*gamma[j]-geta[t,j].*gamma[j]
#         end
#         geta[t,:]=logunif[t,r] < val[1] ? gtry : geta[t,:]
#         hM[t,:]=+geta[t,:]./nfast
#     end
#     mhM=sum(hM, dims=1)./T
#     omega=(hM.-mhM)'*(hM.-mhM)
#
#     obj=mhM*inv(omega)*mhM'
#
#     return obj[1]
# end

# Testing function
#
# function RP_Cournot_test_sem(mchain ; Par::Par=Param)
#     @unpack nfast = Par
#     T, N, nsim = size(mchain)
#     # Generate random index
#     indfast=rand(1:nsim, nfast)
#     indfast[1]=1 # Keeping the first element fixed.
#     mchain_op=mchain[:,:,indfast]
#
#     # First optimization BlackBox
#     res = bboptimize(x->objfun_s_bbo(x, mchain=mchain_op); SearRange = (-10e300,10e300), NumDimensions = N, TraceMode=:silent)
#     ## Getting results
#     minbb = best_fitness(res)
#     γ0 = best_candidate(res)
#     # TS_bb = minbb*N
#     # p_val_bb = ccdf(Chisq(N),TS_bb)
#     println("First optimization from ", myid()," done.")
#     ## Second optimization NLopt
#     fop=nothing
#     fop=Model(NLopt.Optimizer)
#     set_optimizer_attribute(fop, "algorithm", :LN_BOBYQA)
#     set_optimizer_attribute(fop, "xtol_rel", 1e-8)
#     @variable(fop, gamma0[1:N])
#
#     register(fop, :objfun_s, N, (x...)->objfun_s(x..., mchain=mchain_op); autodiff=true )
#     @NLobjective(fop, Min, objfun_s(gamma0...))
#
#     for i=1:N
#         set_start_value(gamma0[i],γ0[i])
#     end
#
#     ## Solving in NLopt
#     JuMP.optimize!(fop)
#
#     ## Getting Results
#     # termination_status(fop)
#     # primal_status(fop)
#
#     γ1=value.(gamma0)
#     obj_v=objective_value(fop)
#     TS=T*obj_v
#     p_val = ccdf(Chisq(N),TS)
#     println("Second optimization from ", myid()," done.")
#     ## Second optimization for refinement
#     for i=1:N
#         set_start_value(gamma0[i],γ1[i])
#     end
#     ## Solving in NLopt
#     JuMP.optimize!(fop)
#
#     ## Getting Results
#     obj_v=objective_value(fop)
#     TS=T*obj_v
#     p_val = ccdf(Chisq(N),TS)
#     println("Third optimization from ", myid()," done.")
#     fop=nothing
#
#     return TS, p_val
# end
##

function RP_Cournot_test_par(mchain ; Par::Par=Param)
    @unpack nfast = Par
    T, N, nsim = size(mchain)
    # Generate random index
    indfast=rand(1:nsim, nfast)
    indfast[1]=1 # Keeping the first element fixed.
    mchain_op=mchain[:,:,indfast]

    # First optimization BlackBox
    res = bboptimize(x->objfun_bbo_par(x, mchain=mchain_op); SearchRange = (-10e300,10e300), NumDimensions = N, TraceMode=:silent, MaxTime=400.0)
    ## Getting results
    minbb = best_fitness(res)
    γ0 = best_candidate(res)
    # TS_bb = minbb*N
    # p_val_bb = ccdf(Chisq(N),TS_bb)
    println("First optimization from ", myid()," done.")
    ## Second optimization NLopt
    fop=nothing
    fop=Model(NLopt.Optimizer)
    set_optimizer_attribute(fop, "algorithm", :LN_BOBYQA)
    set_optimizer_attribute(fop, "xtol_rel", 1e-6)
    set_optimizer_attribute(fop, "maxtime", 300)
    @variable(fop, gamma0[1:N])

    register(fop, :objfun_par, N, (x...)->objfun_par(x..., mchain=mchain_op); autodiff=true )
    @NLobjective(fop, Min, objfun_par(gamma0...))

    for i=1:N
        set_start_value(gamma0[i],γ0[i])
    end

    ## Solving in NLopt
    JuMP.optimize!(fop)

    ## Getting Results
    # termination_status(fop)
    # primal_status(fop)

    γ1=value.(gamma0)
    obj_v=objective_value(fop)
    TS=T*obj_v
    p_val = ccdf(Chisq(N),TS)
    println("Second optimization from ", myid()," done.")
    ## Second optimization for refinement
    for i=1:N
        set_start_value(gamma0[i],γ1[i])
    end
    ## Solving in NLopt
    JuMP.optimize!(fop)

    ## Getting Results
    γ2=value.(gamma0)
    obj_v=objective_value(fop)
    TS=T*obj_v
    p_val = ccdf(Chisq(N),TS)
    println("Third optimization from ", myid()," done.")
    fop=nothing
    acc_r = MyModules.acc_rate(mchain_op, γ2)
    return TS, p_val, acc_r
end

function RP_Cournot_test_par_alt(mchain ; Par::Par=Param)
    @unpack nfast = Par
    T, N, nsim = size(mchain)
    # Generate random index
    indfast=rand(1:nsim, nfast)
    indfast[1]=1 # Keeping the first element fixed.
    mchain_op=mchain[:,:,indfast]

    # First optimization BlackBox
    res = bboptimize(x->objfun_bbo_par_alt(x, mchain=mchain_op); SearchRange = (-10e300,10e300), NumDimensions = N, TraceMode=:silent)
    ## Getting results
    minbb = best_fitness(res)
    γ0 = best_candidate(res)
    # TS_bb = minbb*N
    # p_val_bb = ccdf(Chisq(N),TS_bb)
    println("First optimization from ", myid()," done.")
    ## Second optimization NLopt
    fop=nothing
    fop=Model(NLopt.Optimizer)
    set_optimizer_attribute(fop, "algorithm", :LN_BOBYQA)
    set_optimizer_attribute(fop, "xtol_rel", 1e-6)
    @variable(fop, gamma0[1:N])

    register(fop, :objfun_par_alt, N, (x...)->objfun_par_alt(x..., mchain=mchain_op); autodiff=true )
    @NLobjective(fop, Min, objfun_par_alt(gamma0...))

    for i=1:N
        set_start_value(gamma0[i],γ0[i])
    end

    ## Solving in NLopt
    JuMP.optimize!(fop)

    ## Getting Results
    # termination_status(fop)
    # primal_status(fop)

    γ1=value.(gamma0)
    obj_v=objective_value(fop)
    TS=T*obj_v
    p_val = ccdf(Chisq(N),TS)
    println("Second optimization from ", myid()," done.")
    ## Second optimization for refinement
    for i=1:N
        set_start_value(gamma0[i],γ1[i])
    end
    ## Solving in NLopt
    JuMP.optimize!(fop)

    ## Getting Results
    γ2=value.(gamma0)
    obj_v=objective_value(fop)
    TS=T*obj_v
    p_val = ccdf(Chisq(N),TS)
    println("Third optimization from ", myid()," done.")
    fop=nothing
    acc_r = MyModules.acc_rate(mchain_op, γ2,f=1)
    return TS, p_val, acc_r
end

max_stars_repo_name: hans-mtz/SCR

max_issues_repo_name: hans-mtz/SCR

max_forks_repo_name: hans-mtz/SCR


text:
@info "Generating examples."

using CairoMakie
CairoMakie.activate!(type = "png")

scatter(rand(100), rand(100))

#import Literate
#import BAT
#
#bat_path = dirname(dirname(pathof(BAT)))
#
#tutorial_src = joinpath(bat_path, "docs", "src", "tutorial_lit.jl")
#
#Literate.notebook(tutorial_src, ".", execute = false, name = "bat_tutorial", documenter = true, credit = true)
#Literate.script(tutorial_src, ".", keep_comments = false, name = "bat_tutorial", documenter = true, credit = false)

@info "Examples created."

max_stars_repo_name: greimel/julia-binder

max_issues_repo_name: greimel/julia-binder

max_forks_repo_name: greimel/julia-binder


text:
function part1(input)
    p1, p2 = parse.(Int, last.(split.(readlines(input), ":")))
    s1, s2 = 0, 0
    die = 1
    while true
        p1 = mod1(p1 + 3die + 3, 10)
        die += 3
        s1 += p1
        s1 >= 1000 && return s2 * (die - 1)
        p2 = mod1(p2 + 3die + 3, 10)
        die += 3
        s2 += p2
        s2 >= 1000 && return s1 * (die - 1)
    end
end

function part2(input)
    p1, p2 = parse.(Int, last.(split.(readlines(input), ":")))
    x1 = Dict((p1, 0) => 1)
    x2 = Dict((p2, 0) => 1)
    sum1 = 1
    sum2 = 2
    y1 = empty(x1)
    y2 = empty(x2)
    w1, w2 = 0, 0
    dice = ((3, 1), (4, 3), (5, 6), (6, 7), (7, 6), (8, 3), (9, 1))
    while !isempty(x1)
        empty!(y1)
        sum1 = 0
        for ((p1, s1), n) in x1
            for (d, m) in dice
                pp1 = mod1(p1 + d, 10)
                ss1 = s1 + pp1
                if ss1 >= 21
                    w1 += m * n * sum2
                else
                    y1[(pp1, ss1)] = get(y1, (pp1, ss1), 0) + m * n
                    sum1 += m * n
                end
            end
        end
        x1, y1 = y1, x1
        empty!(y2)
        sum2 = 0
        for ((p2, s2), n) in x2
            for (d, m) in dice
                pp2 = mod1(p2 + d, 10)
                ss2 = s2 + pp2
                if ss2 >= 21
                    w2 += m * n * sum1
                else
                    y2[(pp2, ss2)] = get(y2, (pp2, ss2), 0) + m * n
                    sum2 += m * n
                end
            end
        end
        x2, y2 = y2, x2
    end
    return max(w1, w2)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libavtp_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libavtp")
JLLWrappers.@generate_main_file("libavtp", UUID("c33e1e51-3cb7-5762-8964-d0cfbd951039"))
end  # module libavtp_jll


text:
"""
    eom_mediated_link(u, p, t) -> SVector{9}

Equations of motion for a three-subsystem dynamical system where `X` and `Y` are
driven by `Z`. At the default value of the coupling constant `c = 0`, the
responses `X` and `Y` are already synchronized to the driver `Z`.
The frequencies `ωx`, `ωy` and ``ωz` have near unit ratios by default.

# References
Krakovská, Anna, et al. "Comparison of six methods for the detection of causality in a bivariate time series." Physical Review E 97.4 (2018): 042207
"""
function eom_mediated_link(u, p, t)
    ωx, ωy, ωz, k, l, m, c = (p...,)
    x₁, x₂, x₃, y₁, y₂, y₃, z₁, z₂, z₃ = (u...,)

    dx₁ = -ωx*x₂ - x₃ + c*(z₁ - x₁)
	dx₂ = ωx*x₁ + k*x₂
	dx₃ = l + x₃*(x₁ - m)

	dy₁ = -ωy*y₂ - y₃ + c*(z₁ - y₁)
	dy₂ = ωy*y₁ + k*y₂
	dy₃ = l + y₃*(y₁ - m)

	dz₁ = -ωz*z₂ - z₃
	dz₂ = ωz*z₁ + k*z₂
	dz₃ = l + z₃*(z₁ - m)

    SVector{9}(dx₁, dx₂, dx₃, dy₁, dy₂, dy₃, dz₁, dz₂, dz₃)
end

function mediated_link(u₀, ωx, ωy, ωz, k, l, m, c)
    p = [ωx, ωy, ωz, k, l, m, c]
    ContinuousDynamicalSystem(eom_mediated_link, u₀, p)
end

"""
    mediated_link(;u₀ = rand(9), ωx = 1, ωy = 1.015, ωz = 0.985,
        k = 0.15, l = 0.2, m = 10.0, c = 0.06) -> ContinuousDynamicalSystem

Equations of motion for a three-subsystem dynamical system where `X` and `Y` are
driven by `Z`. At the default value of the coupling constant `c = 0`, the
responses `X` and `Y` are already synchronized to the driver `Z`.

The dynamics is generated by the vector field:

```math
\\begin{aligned}
dx_1 &= -\\omega_x x_2 - x_3 + c*(z_1 - x_1) \\
dx_2 &= \\omega_x x_1 + k*x_2 \\
dx_3 &= l + x_3(x_1 - m) \\
dy_1 &= -\\omega_y y_2 - y_3 + c*(z_1 - y_1) \\
dy_2 &= \\omega_y y_1 + k*y_2 \\
dy_3 &= l + y_3(y_1 - m) \\
dz_1 &= -\\omega_z z_2 - z_3 \\
dz_2 &= \\omega_z z_1 + k*z_2 \\
dz_3 &= l + z_3(z_1 - m) \\
\end{aligned}
```

# References
Krakovská, Anna, et al. "Comparison of six methods for the detection of causality in a bivariate time series." Physical Review E 97.4 (2018): 042207
"""
mediated_link(;u₀ = rand(9), ωx = 1, ωy = 1.015, ωz = 0.985,
            k = 0.15, l = 0.2, m = 10.0, c = 0.06) =
    mediated_link(u₀, ωx, ωy, ωz, k, l, m, c)


text:
function putmarkers1!(r::Robot,side::HorizonSide,num_steps::Int)
    if isborder(r,side)==false
        move!(r,side)
    elseif isborder(r,side)==true
        k=0
        while isborder(r,side)==true
            move!(r,Nord)
            k+=1
        end
        move!(r,side)
        if isborder(r,Sud)==true
            move!(r,side)
        end
        for _ in 1:k
            move!(r,Sud)
        end
    end
    go_and_obhod(r,side,num_steps)
    putmarker!(r)
    if isborder(r,Nord)==false
        move!(r,Nord)
    end
end


text:
"""
            Kernels, Kernelvalues and -derivatives for sph mapping.

    Author: Ludwig Böss
    Contact: lboess@usm.lmu.de
    Created: 2018-12-12

"""

#using Distributed

abstract type SPHKernel end

"""
            Cubic
"""
# struct
struct Cubic <: SPHKernel
    n_neighbours::Int64
    norm_2D::Float64
    norm_3D::Float64
    function Cubic(n_neighbours::Int64=64)
        new(n_neighbours, 8.0/π, 8.0/π)
    end
end

# kernel values
@inline function kernel_value_2D(kernel::Cubic, u::Float64, h_inv::Float64)

    n = kernel.norm_3D * h_inv^2

    if u < 0.5
        return ( 1.0 - 6.0 * (1.0 - u) * u^2) * n
    elseif u < 1.0
        return ( 2.0 * (1.0 - u) * (1.0 - u) * (1.0 - u)) * n
    else
        return 0.
    end

end

@inline function kernel_value_3D(kernel::Cubic, u::Float64, h_inv::Float64)

    n = kernel.norm_3D * h_inv^3

    if u < 0.5
        return ( 1.0 + 6.0 * (u - 1.0) * u^2) * n
    elseif u < 1.0
        return ( 2.0 * (1.0 - u) * (1.0 - u) * (1.0 - u)) * n
    else
        return 0.
    end

end

"""
            Quintic
"""
# struct
struct Quintic <: SPHKernel
    n_neighbours::Int64
    norm_2D::Float64
    norm_3D::Float64
    function Quintic(n_neighbours::Int64=200)
        new(n_neighbours, 15309.0/(478.0*π), 2187.0/(40.0*pi))
    end
end

# kernel values
@inline function kernel_value_2D(kernel::Quintic, u::Float64, h_inv::Float64)

    n = kernel.norm_2D * h_inv^2

    if u < 1.0/3.0
        u_m1  = ( 1.0 - u )
        u_m23 = ( 2.0/3.0 - u )
        u_m13 = ( 1.0/3.0 - u )
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 - 6.0 * u_m23*u_m23*u_m23*u_m23*u_m23 + 15.0 * u_m13*u_m13*u_m13*u_m13*u_m13 ) * n
    elseif u < 2.0/3.0
        u_m1  = ( 1.0 - u )
        u_m23 = ( 2.0/3.0 - u )
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 - 6.0 * u_m23*u_m23*u_m23*u_m23*u_m23 ) * n
    elseif u < 1.0
        u_m1  = ( 1.0 - u )
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 ) * n
    else
        return 0.
    end

end

@inline function kernel_value_3D(kernel::Quintic, u::Float64, h_inv::Float64)

    n = kernel.norm_3D * h_inv^3

    if u < 1.0/3.0
        u_m1  = ( 1.0 - u )
        u_m23 = ( 2.0/3.0 - u )
        u_m13 = ( 1.0/3.0 - u )
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 - 6.0 * u_m23*u_m23*u_m23*u_m23*u_m23 + 15.0 * u_m13*u_m13*u_m13*u_m13*u_m13 ) * n
    elseif u < 2.0/3.0
        u_m1  = ( 1.0 - u )
        u_m23 = ( 2.0/3.0 - u )
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 - 6.0 * u_m23*u_m23*u_m23*u_m23*u_m23 ) * n
    elseif u < 1.0
        u_m1  = ( 1.0 - u )
        return ( u_m1*u_m1*u_m1*u_m1*u_m1 ) * n
    else
        return 0.
    end

end

"""
            Wendland C4
"""
# struct
struct WendlandC4 <: SPHKernel
    n_neighbours::Int64
    norm_2D::Float64
    norm_3D::Float64
    function WendlandC4(n_neighbours::Int64=200)
        new(n_neighbours,9.0/π, 495.0/(32.0 * π))
    end
end

# kernel values
@inline function kernel_value_2D(kernel::WendlandC4, u::Float64, h_inv::Float64)

    if u < 1.0
        n = kernel.norm_2D * h_inv^2
        u_m1 = 1.0 - u
        u_m1_2 = u_m1 * u_m1  # (1.0 - u)^2
        u_m1_4 = u_m1 * u_m1  # (1.0 - u)^4
        return ( u_m1_2*u_m1_4 * ( 1.0 + 6u + 35.0/3.0 * u^2 ) ) * n
    else
        return 0.
    end

end

@inline function kernel_value_3D(kernel::WendlandC4, u::Float64, h_inv::Float64)

    if u < 1.0
        n = kernel.norm_3D * h_inv^3
        u_m1 = 1.0 - u
        u_m1_2 = u_m1 * u_m1  # (1.0 - u)^2
        u_m1_4 = u_m1 * u_m1  # (1.0 - u)^4
        return ( u_m1_2*u_m1_4 * ( 1.0 + 6u + 35.0/3.0 * u^2 ) ) * n
    else
        return 0.
    end

end

"""
            Wendland C6
"""
# struct
struct WendlandC6 <: SPHKernel
    n_neighbours::Int64
    norm_2D::Float64
    norm_3D::Float64
    function WendlandC6(n_neighbours::Int64=295)
        new(n_neighbours, 78.0/(7.0*π), 1365.0/(64.0*π))
    end
end

# kernel values
@inline function kernel_value_2D(kernel::WendlandC6, u::Float64, h_inv::Float64)

    if u < 1.0
        n = kernel.norm_2D * h_inv^2
        u_m1 = 1.0 - u
        u_m1 = u_m1 * u_m1  # (1.0 - u)^2
        u_m1 = u_m1 * u_m1  # (1.0 - u)^4
        u_m1 = u_m1 * u_m1  # (1.0 - u)^8
        u2 = u*u
        return ( u_m1 * ( 1.0 + 8u + 25u2 + 32u2*u )) * n
    else
       return 0.0
   end

end

@inline function kernel_value_3D(kernel::WendlandC6, u::Float64, h_inv::Float64)

    if u < 1.0
        n = kernel.norm_3D * h_inv^3
        u_m1 = 1.0 - u
        u_m1 = u_m1 * u_m1  # (1.0 - u)^2
        u_m1 = u_m1 * u_m1  # (1.0 - u)^4
        u_m1 = u_m1 * u_m1  # (1.0 - u)^8
        u2 = u*u
        return ( u_m1 * ( 1.0 + 8u + 25u2 + 32u2*u )) * n
    else
       return 0.0
   end

end

# @inline function kernel_deriv(kernel::Cubic, u::Float64, h::Float64)
#
#     norm = 8.0/π
#     n = norm/h^4
#
#     if u < 0.5
#         return ( u * (18.0 * u - 12.0 )) * n
#     elseif u < 1.0
#         return ( -6.0 * (1.0 - u) * (1.0 - u) ) * n
#     else
#         return 0.
#     end
#
# end

# @inline function kernel_deriv(kernel::Quintic, u::Float64, h::Float64)
#
#     norm = ( 2187.0 / ( 40. * π))
#     n = norm/h^4
#
#     if u < 1.0/3.0
#         return ( -5.0 * ( 1.0 - u )^4 + 30.0 * ( 2.0/3.0 - u )^4  - 75.0 * ( 1.0/3.0 - u )^4 ) * n
#     elseif u < 2.0/3.0
#         return ( -5.0 * ( 1.0 - u )^4 + 30.0 * ( 2.0/3.0 - u )^4 - 75.0 ) * n
#     elseif u < 1.0
#         return ( -5.0 * ( 1.0 - u )^4 ) * n
#     else
#         return 0.
#     end
#
# end

# @inline function kernel_deriv(kernel::WendlandC4, u, h)
#
#     norm = 495.0/(32. * π)
#     n = norm/h^4
#
#     if u < 1.0
#         return ( -288.0/3.0 * ( 1. - u )^5 * u^2 - 56.0/3.0 * u * ( 1. - u)^5 ) * n
#     else
#         return 0.
#     end
#
# end

# @inline function kernel_deriv(kernel::WendlandC6, u::Float64, h::Float64)
#
#     norm = 1365.0/(64.0*π)
#     n = norm/h^4
#
#     if u < 1.0
#         return ( -22. * (1.0 - u)^7 * u * ( 16. * u^2 + 7. * u + 1. )) * n
#     else
#         return 0.
#     end
#
# end

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
const outputname = joinpath(@__DIR__, "currency-data.jl")

# First, check if currency-data.jl already exists
isfile(outputname) && exit()

using JSON3

const src = "https://pkgstore.datahub.io/core/country-codes/country-codes_json/data/471a2e653140ecdd7243cdcacfd66608/country-codes_json.json"

const inputname = joinpath(@__DIR__, "country-codes.json")

const currency_list = Dict{String,Tuple{Int,Int,String}}()

const SymCurr = Symbol("ISO4217-currency_alphabetic_code")
const SymUnit = Symbol("ISO4217-currency_minor_unit")
const SymName = Symbol("ISO4217-currency_name")
const SymCode = Symbol("ISO4217-currency_numeric_code")

function genfile(io)
    for country in country_list
        (abbrlist = country[SymCurr]) === nothing && continue
        (unitlist = country[SymUnit]) === nothing && continue
        (namelist = country[SymName]) === nothing && continue
        (codelist = country[SymCode]) === nothing && continue
        currencies = split(abbrlist, ',')
        units = split(string(unitlist), ',')
        names = split(namelist, ',')
        codes = split(string(codelist), ',')

        for (curr, unit, code, name) in zip(currencies, units, codes, names)
            length(curr) != 3 && continue
            haskey(currency_list, curr) && continue
            currency_list[curr] = (parse(Int, unit), parse(Int, code), string(name))
        end
    end
    println(io, "const _currency_data = Dict(")
    for (curr, val) in currency_list
        println(io, "    :$curr => (Currency{:$curr}, $(val[1]), $(lpad(val[2], 4)), \"$(val[3])\"),")
    end
    println(io, ")\n")
end

# Only download the file from datahub.io if not already present
if !isfile(inputname)
    println("Downloading currency data: ", src)
    download(src, inputname)
end

const country_list = open(io -> JSON3.read(io), inputname)

open(io -> genfile(io), outputname, "w")

max_stars_repo_name: TotalVerb/Currencies

max_issues_repo_name: TotalVerb/Currencies

max_forks_repo_name: TotalVerb/Currencies


text:
# constants.jl
#  - Physical constants for Orbital 3 Density Functional Theory

# Constants:
const hbar = 1.054E-34
const h =    6.62606957E-34
const kb =   8.6173324E-5 # in units of eV
const ε_0 =  8.854E-12 #Units: C2N−1m−2, permittivity of free space
const q = eV = 1.60217657E-19
const c = 3E8 # Speed of light in this universe
const k_e = 8.98755E+9 # Coulomb's constant
const me = 9.10938291E-31 # Mass of electron (kg)
const rbohr = 5.2917721067E-11 # Bohr radius (m)

const Ha=27.2111eV
const Ry=Ha/2

const Å=1E-10 # Angstrom


text:
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

type CreateEncryptionZoneRequestProto
    src::AbstractString
    keyName::AbstractString
    CreateEncryptionZoneRequestProto() = (o=new(); fillunset(o); o)
end #type CreateEncryptionZoneRequestProto
const __req_CreateEncryptionZoneRequestProto = Symbol[:src]
meta(t::Type{CreateEncryptionZoneRequestProto}) = meta(t, __req_CreateEncryptionZoneRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::CreateEncryptionZoneRequestProto) = ProtoBuf.protohash(v)
isequal(v1::CreateEncryptionZoneRequestProto, v2::CreateEncryptionZoneRequestProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::CreateEncryptionZoneRequestProto, v2::CreateEncryptionZoneRequestProto) = ProtoBuf.protoeq(v1, v2)

type CreateEncryptionZoneResponseProto
    CreateEncryptionZoneResponseProto() = (o=new(); fillunset(o); o)
end #type CreateEncryptionZoneResponseProto
hash(v::CreateEncryptionZoneResponseProto) = ProtoBuf.protohash(v)
isequal(v1::CreateEncryptionZoneResponseProto, v2::CreateEncryptionZoneResponseProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::CreateEncryptionZoneResponseProto, v2::CreateEncryptionZoneResponseProto) = ProtoBuf.protoeq(v1, v2)

type ListEncryptionZonesRequestProto
    id::Int64
    ListEncryptionZonesRequestProto() = (o=new(); fillunset(o); o)
end #type ListEncryptionZonesRequestProto
const __req_ListEncryptionZonesRequestProto = Symbol[:id]
meta(t::Type{ListEncryptionZonesRequestProto}) = meta(t, __req_ListEncryptionZonesRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::ListEncryptionZonesRequestProto) = ProtoBuf.protohash(v)
isequal(v1::ListEncryptionZonesRequestProto, v2::ListEncryptionZonesRequestProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::ListEncryptionZonesRequestProto, v2::ListEncryptionZonesRequestProto) = ProtoBuf.protoeq(v1, v2)

type EncryptionZoneProto
    id::Int64
    path::AbstractString
    suite::Int32
    cryptoProtocolVersion::Int32
    keyName::AbstractString
    EncryptionZoneProto() = (o=new(); fillunset(o); o)
end #type EncryptionZoneProto
const __req_EncryptionZoneProto = Symbol[:id,:path,:suite,:cryptoProtocolVersion,:keyName]
meta(t::Type{EncryptionZoneProto}) = meta(t, __req_EncryptionZoneProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::EncryptionZoneProto) = ProtoBuf.protohash(v)
isequal(v1::EncryptionZoneProto, v2::EncryptionZoneProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::EncryptionZoneProto, v2::EncryptionZoneProto) = ProtoBuf.protoeq(v1, v2)

type ListEncryptionZonesResponseProto
    zones::Array{EncryptionZoneProto,1}
    hasMore::Bool
    ListEncryptionZonesResponseProto() = (o=new(); fillunset(o); o)
end #type ListEncryptionZonesResponseProto
const __req_ListEncryptionZonesResponseProto = Symbol[:hasMore]
meta(t::Type{ListEncryptionZonesResponseProto}) = meta(t, __req_ListEncryptionZonesResponseProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::ListEncryptionZonesResponseProto) = ProtoBuf.protohash(v)
isequal(v1::ListEncryptionZonesResponseProto, v2::ListEncryptionZonesResponseProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::ListEncryptionZonesResponseProto, v2::ListEncryptionZonesResponseProto) = ProtoBuf.protoeq(v1, v2)

type GetEZForPathRequestProto
    src::AbstractString
    GetEZForPathRequestProto() = (o=new(); fillunset(o); o)
end #type GetEZForPathRequestProto
const __req_GetEZForPathRequestProto = Symbol[:src]
meta(t::Type{GetEZForPathRequestProto}) = meta(t, __req_GetEZForPathRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::GetEZForPathRequestProto) = ProtoBuf.protohash(v)
isequal(v1::GetEZForPathRequestProto, v2::GetEZForPathRequestProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::GetEZForPathRequestProto, v2::GetEZForPathRequestProto) = ProtoBuf.protoeq(v1, v2)

type GetEZForPathResponseProto
    zone::EncryptionZoneProto
    GetEZForPathResponseProto() = (o=new(); fillunset(o); o)
end #type GetEZForPathResponseProto
hash(v::GetEZForPathResponseProto) = ProtoBuf.protohash(v)
isequal(v1::GetEZForPathResponseProto, v2::GetEZForPathResponseProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::GetEZForPathResponseProto, v2::GetEZForPathResponseProto) = ProtoBuf.protoeq(v1, v2)

export CreateEncryptionZoneRequestProto, CreateEncryptionZoneResponseProto, ListEncryptionZonesRequestProto, EncryptionZoneProto, ListEncryptionZonesResponseProto, GetEZForPathRequestProto, GetEZForPathResponseProto


text:
using ArgParse
using DataFrames
using Statistics

import Colors

function parse_segfree_commandline(args::Union{Nothing, Array{String, 1}}=nothing)
    s = ArgParseSettings(prog="baysor preview")
    @add_arg_table! s begin
        "--k", "-k"
            help = "Number of neighbors for segmentation-free pseudo-cells. It's suggested to set it to the expected minimal number of molecules per cell."
        "--config", "-c"
            help = "TOML file with config"
        "--x-column", "-x"
            help = "Name of x column. Overrides the config value."
        "--y-column", "-y"
            help = "Name of gene column. Overrides the config value."
        "--z-column", "-z"
            help = "Name of gene column. Overrides the config value."
        "--gene-column", "-g"
            help = "Name of gene column. Overrides the config value."
        "--exclude-genes"
            help = "Comma-separated list of genes to ignore during segmentation"
        "--output", "-o"
            help = "Name of the output file or path to the output directory"
            default = "preview.html"

        "coordinates"
            help = "CSV file with coordinates of transcripts and gene type"
            required = true
    end

    return (args === nothing) ? parse_args(s) : parse_args(args, s)
end

function parse_segfree_configs(args::Union{Nothing, Array{String, 1}}=nothing)
    r = parse_segfree_commandline(args)
    if r === nothing
        return nothing
    end

    if r["config"] !== nothing
        extend_params_with_config!(r, parse_toml_config(r["config"]))
    else
        extend_params_with_config!(r, get_default_config())
    end

    for k in ["gene-column", "x-column", "y-column", "z-column"]
        r[k] = Symbol(r[k])
    end

    return r
end

function run_cli_segfree(args::Union{Nothing, Array{String, 1}}=nothing)
    args = parse_segfree_configs(args)

    # Set up logger

    log_file = open(append_suffix(args["output"], "segfree_log.log"), "w")
    Base.CoreLogging.global_logger(DoubleLogger(log_file, stdout; force_flush=true))

    # Run preview

    @info get_baysor_run_str()
    @info "Loading data..."
    args["min-molecules-per-gene"] = 0
    df_spatial, gene_names = load_df(args)

    @info "Loaded $(size(df_spatial, 1)) transcripts"
    if size(df_spatial, 1) != size(unique(df_spatial), 1)
        @warn "$(size(df_spatial, 1) - size(unique(df_spatial), 1)) records are duplicates. You may need to filter them beforehand."
    end

    @error "NOT IMPLEMENTED" # TODO

    close(log_file)

    return 0
end

max_stars_repo_name: repushko/Baysor

max_issues_repo_name: repushko/Baysor

max_forks_repo_name: repushko/Baysor


text:
#osaka.jl

"""
    Implementation of Osaka's extended treatment of Feynman's variational technique applied to the Polaron model; generalising it from the case at O^{∘}K to the case at finite temperature.

    See Osaka 1959:
    Progress of Theoretical Physics, Vol. 22, No.3, September 1959

    And Hellwarth 1999:
    https://doi.org/10.1103/PhysRevB.60.299

"""

using Optim
using QuadGK

# Equation 62b in Hellwarth. Equation 20 in Osaka.
A(v, w, β) = 3 / β * (log(v / w) - 1 / 2 * log(2 * π * β) - log(sinh(v * β / 2) / sinh(w * β / 2)))

# Equation 62d in Hellwarth.
Y(x, v, β) = 1 / (1 - exp(-v * β)) * (1 + exp(-v * β) - exp(-v * x) - exp(v * (x - β)))

# Osaka's version: Within Equation 16.
# Y(x, v, β) = 1 - exp(-v * x) + (1 - coth(v * β / 2)) * (cosh(v * x) - 1)

# Integrand of Equation 62c in Hellwarth.
B_integrand(x, v, w, β) =  (exp(β - x) + exp(x)) / sqrt(w^2 * x * (1 - x / β) + Y(x, v, β) * (v^2 - w^2) / v)

# Osaka's version: Within Equation 16.
# B_integrand(x, v, w, β) =  (exp(-x)) / sqrt(abs(w^2 * x * (1 - x / β) + Y(x, v, β) * (v^2 - w^2) / v))

# Equation 62c. NB: Extra exp(β) correction term from Osaka to correct Hellwarth typo.
B(v, w, β, α) = α * v / (sqrt(π) * (exp(β) - 1)) * QuadGK.quadgk(x -> B_integrand(x, v, w, β), 0, β / 2)[1]

# Osaka's version: Equation 16.
# B(v, w, β, α) = α * v * exp(β) / (sqrt(π) * (exp(β) - 1)) * QuadGK.quadgk(x -> B_integrand(x, v, w, β), 0, β)[1]

# Equation 62e in Hellwarth. Equation 17 in Osaka.
C(v, w, β) = 3 / 4 * (v^2 - w^2) / v * (coth(v * β / 2) - 2 / (v * β))

# Equation 62a in Hellwarth. In paragraph below Equation 22 in Osaka; has extra 1/β due to different definition of A, B & C.
osaka_free_energy(v, w, β, α) = -(A(v, w, β) + B(v, w, β, α) + C(v, w, β))


text:
tol = 1e-10

@testset "Grid estimator" begin
    points_2D = rand(2, 200)
    points_3D = rand(3, 400)
    E_2D = invariantize(StateSpaceReconstruction.cembed(points_2D))
    E_3D = invariantize(StateSpaceReconstruction.cembed(points_3D))
    ϵ = 3
    bins_visited_by_orbit_2D = assign_bin_labels(E_2D, ϵ)
    bins_visited_by_orbit_3D = assign_bin_labels(E_3D, ϵ)

    bininfo_2D = get_binvisits(bins_visited_by_orbit_2D)
    bininfo_3D = get_binvisits(bins_visited_by_orbit_3D)

    TO_2D = estimate_transferoperator_from_binvisits(bininfo_2D)
    TO_3D = estimate_transferoperator_from_binvisits(bininfo_3D)

    invm_2D = invariantmeasure(TO_2D)
    invm_3D = invariantmeasure(TO_3D)

    @test all(invm_2D.dist .>= -tol)
    @test sum(invm_2D.dist) <= 1 + tol || sum(invm_2D.dist) ≈ 1
end


text:
using CmdStan

ProjDir = @__DIR__
cd(ProjDir)

bernoulli_model = "
  functions{
    #include model_specific_funcs.stan
    #include shared_funcs.stan // a comment
    //#include shared_funcs.stan // a comment
  }
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    model_specific_function();
    theta ~ beta(my_function(),1);
    y ~ bernoulli(theta);
  }
";

observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

stanmodel = Stanmodel(Sample(save_warmup=false, num_warmup=1000,
  num_samples=2000, thin=1), name="bernoulli", model=bernoulli_model,
  printsummary=false, tmpdir=mktempdir());

rc, a3d, cnames = stan(stanmodel, observeddata, ProjDir);

if rc == 0
  sdf  = read_summary(stanmodel)
  @test sdf[sdf.parameters .== :theta, :mean][1] ≈ 0.34 rtol=0.1
end


text:
using Rotations
using BenchmarkTools
import Base.Iterators: product
import Random

const T = Float64

const suite = BenchmarkGroup()
Random.seed!(1)

noneuler = suite["Non-Euler conversions"] = BenchmarkGroup()
rotationtypes = [RotMatrix3{T}, Quat{T}, UnitQuaternion{T}, SPQuat{T}, MRP{T}, AngleAxis{T}, RodriguesVec{T}]
for (from, to) in product(rotationtypes, rotationtypes)
    if from != to
        name = "$(string(from)) -> $(string(to))"
        # use eval here because of https://github.com/JuliaCI/BenchmarkTools.jl/issues/50#issuecomment-318673288
        noneuler[name] = eval(:(@benchmarkable $to(rot) setup = rot = rand($from)))
    end
end

euler = suite["Euler conversions"] = BenchmarkGroup()
eulertypes = [
    RotX{T}, RotY{T}, RotZ{T},
    RotXY{T}, RotYX{T}, RotZX{T}, RotXZ{T}, RotYZ{T}, RotZY{T},
    RotXYX{T}, RotYXY{T}, RotZXZ{T}, RotXZX{T}, RotYZY{T}, RotZYZ{T},
    RotXYZ{T}, RotYXZ{T}, RotZXY{T}, RotXZY{T}, RotYZX{T}, RotZYX{T}]
for from in eulertypes
    to = RotMatrix3{T}
    name = "$(string(from)) -> $(string(to))"
    # use eval here because of https://github.com/JuliaCI/BenchmarkTools.jl/issues/50#issuecomment-318673288
    euler[name] = eval(:(@benchmarkable $to(rot) setup = rot = rand($from)))
end

composition = suite["Composition"] = BenchmarkGroup()
# use eval here because of https://github.com/JuliaCI/BenchmarkTools.jl/issues/50#issuecomment-318673288
composition["RotMatrix{3} * RotMatrix{3}"] = eval(:(@benchmarkable r1 * r2 setup = begin
    r1 = rand(RotMatrix3{T})
    r2 = rand(RotMatrix3{T})
end))

paramspath = joinpath(@__DIR__, "benchmarkparams.json")
if isfile(paramspath)
    loadparams!(suite, BenchmarkTools.load(paramspath)[1], :evals, :samples);
else
    tune!(suite, verbose = true)
    BenchmarkTools.save(paramspath, BenchmarkTools.params(suite))
end

results = run(suite, verbose=true)
println()
for (groupname, groupresults) in results
    println("Group: $groupname")
    max_name_length = maximum(length, keys(groupresults))
    for (name, trial) in groupresults
        if minimum(trial).allocs != 0
            @warn("$name  allocates!")
        end
        println(rpad(name, max_name_length), " ", BenchmarkTools.prettytime(minimum(trial).time))
    end
    println()
end


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

using Markdown
using InteractiveUtils

# ╔═╡ dc9d1ef7-23e2-4bc3-ab47-e9ca80af5811
begin
	using PrettyTables
	using PlutoUI
end

# ╔═╡ bd0378cb-ca70-4803-b5fc-7bd08ad1b9f3
md"""
#### Practice Problem: Multiple Chemical Reactions in a Single CSTR

Component $A$ is converted to $B$ (desired product) and $C$ (undesired product) in single well-mixed continuous stirred tank reactor (CSTR) with a single input and output and the following reaction chemistry:

$$\begin{eqnarray}
\hat{r}_{1}: A &\rightleftharpoons& B\\
\hat{r}_{2}: A &\rightleftharpoons& C\\
\end{eqnarray}$$

The substrate $A$ (component 1) is supplied to the reactor at 10 mol/hr, where the input stream has a volumetric flow rate of $\dot{F}_{1}$ = 2 L/hr. The fractional conversion of $A$ was measured and found to be equal to 70%. Lastly, the extent of the desired reaction was estimated to $\dot{\epsilon}_{1}$ = 1.5 mol/hr.

__Assumptions__:
* The CSTR is well-mixed, operating at steady-state, and has a constant volume V = 14L
* The CSTR has a single input stream and a single output stream
* Neither of the chemical reactions occur in the streams
* There is no $B$ or $C$ in the input stream

__Compute__:
* The chemical reaction rate(s) $\hat{r}_{1}$ and $\hat{r}_{2}$
* The mol flow rates of all components in the CSTR exit stream
* If $\hat{r}_{1}$ and $\hat{r}_{2}$ are irreversible, and first-order in $A$, compute the forward rate constants $k_{1}$ and $k_{2}$ for reaction 1 and 2, respectively.

"""

# ╔═╡ 6b07d7ae-8b5d-4c41-8184-9329cdaa94a6
begin

	# Setup problem parameters -
	F_dot = 2.0 	# units: L/hr
	f1 = 0.70 		# units: dimensionless
	ϵ₁_dot = 1.5 	# units: mol/hr
	n11_dot = 10.0 	# units: mol/hr
	V = 14.0 		# units: L

	# setup stoichiometric matrix for "easy" computation -
	S = [

		# ϵ₁_dot ϵ₂_dot
		-1.0 -1.0 	; # 1 A
		1.0 0.0 	; # 2 B
		0.0 1.0 	; # 3 C
	];

	# setup input flow rate -
	n_dot_in = [
		n11_dot ; 	# 1 A
		0.0 	; 	# 2 B
		0.0 	; 	# 3 C
	];

	# show -
	nothing
end

# ╔═╡ 7e8aaae7-555f-4b70-bae7-4adc1d1c9923
md"""
#### Compute: the chemical reaction rate(s) $\hat{r}_{1}$ and $\hat{r}_{2}$

To compute the reaction rates $\hat{r}_{1}$ and $\hat{r}_{2}$, we first must estimate the open extent for the second reaction $\dot{\epsilon}_{2}$
by relating the $\dot{\epsilon}_{2}$ to the fractional conversion of A. To start, we know that for component $i$ the species mole balance is given by:

$$\dot{n}_{i1} - \dot{n}_{i2} + \sum_{j=1}^{\mathcal{R}}\sigma_{ij}\dot{\epsilon}_{j} = 0$$

where $\epsilon_{j} = \hat{r}_{j}V$. Rearranging the species mole balance for $\dot{n}_{i2}$:

$$\dot{n}_{i2} = \dot{n}_{i1} + \sum_{j=1}^{\mathcal{R}}\sigma_{ij}\dot{\epsilon}_{j}$$

However, the fractional conversion is given by:

$$\dot{n}_{i2} = \dot{n}_{i1} - f_{i}\dot{n}_{i1}$$

Thus, we can eliminate $\dot{n}_{i2}$ from the species mol balance ans solve for the fractional conversion $f_{i}$:

$$f_{i} = -\frac{1}{\dot{n_{i,1}}}\left(\sum_{j=1}^{\mathcal{R}}\sigma_{ij}\dot{\epsilon}_{j}\right)\qquad{i=1}$$

For the currenr problem, we know the fractional concersion of A (i=1), and we have two reactions where $\dot{\epsilon}_{2}$ is known, thus
we can solve for $\dot{\epsilon}_{2}$:

$$\dot{\epsilon}_{2} = \dot{n}_{11}f_{1} - \dot{\epsilon}_{1}$$

Once we have $\dot{\epsilon}_{1}$ and $\dot{\epsilon}_{2}$, we can compute the rates: $\hat{r}_{j}=\dot{\epsilon}_{j}V^{-1}$.

"""

# ╔═╡ fddd2b82-8eba-4517-bdfd-0953bd10c259
begin

	# compute -
	ϵ₂_dot = n11_dot*f1  - ϵ₁_dot;

	# display -
	with_terminal() do
		# compute ϵ₂_dot:
		println("ϵ₂_dot = $(ϵ₂_dot) mol/hr")
	end
end

# ╔═╡ 86964693-9968-4524-a2bc-8c3a8f202cc1
begin

	# compute the rates -
	r1_hat = (1/V)*ϵ₁_dot
	r2_hat = (1/V)*ϵ₂_dot

	# display -
	with_terminal() do
		println("r1_hat = $(r1_hat) mol/L-hr and r2_hat = $(r2_hat) mol/L-hr")
	end
end

# ╔═╡ 34ddef24-462a-4a1d-8a61-9c760f278883
md"""

#### Compute: the mole flow rates of all components in the CSTR exit stream

We know have both $\dot{\epsilon}_{1}$ and $\dot{\epsilon}_{2}$, so we can use the steady-state species mole balance for component $i=1,2,3$

$$\dot{n}_{i1} - \dot{n}_{i2} + \sum_{j=1}^{\mathcal{R}}\sigma_{ij}\dot{\epsilon}_{j} = 0\qquad{i=1,2,3}$$

to compute the composition of the exit stream.

"""

# ╔═╡ fff53850-a627-4349-bef2-c5aaf8828c3d
begin

	# setup extent vector -
	ϵ_vector = [ϵ₁_dot ; ϵ₂_dot];

	# compute n_dot coming out of reactor -
	n_dot_out = n_dot_in + S*ϵ_vector

	# setup data array =
	state_table_array = Array{Any,2}(undef, 3,3)

	# A -
	state_table_array[1,1] = "A"
	state_table_array[1,2] = n_dot_in[1]
	state_table_array[1,3] = n_dot_out[1]

	# B -
	state_table_array[2,1] = "B"
	state_table_array[2,2] = n_dot_in[2]
	state_table_array[2,3] = n_dot_out[2]

	# C -
	state_table_array[3,1] = "C"
	state_table_array[3,2] = n_dot_in[3]
	state_table_array[3,3] = n_dot_out[3]

	# setup a pretty table =
	# header row -
	state_table_header_row = (["Species","n_dot_in","n_dot_out"],["","mol/hr","mol/hr"]);

	with_terminal() do
		# write the table -
		pretty_table(state_table_array; header=state_table_header_row)
	end
end

# ╔═╡ 8b5a451f-5d31-49f5-91fb-3e24a00b87ac
md"""
#### Compute the forward rate constants $k^{+}_{1}$ and $k^{+}_{2}$

To compute the rate constants, we can solve the rate expression $\hat{r}_{j} = k^{+}_{j}\left[A\right]$ once we have the concentration of A.
To compute concentration, we know that in stream $s$

$$\left[C_{i,s}\right] = \dot{n}_{i,s}\left(\dot{F}_{s}\right)^{-1}$$

where $\left[C_{i,s}\right]$ denotes the concentration of component $i$ in stream $s$, and $\dot{F}_{s}$ denotes the volumetric flow rate of stream $s$. Because the reactor is well-mixed, and is constant volume, the concentration of $A$ (denoted in index notation as $C_{12}$) is given by:

$$C_{12} = \dot{n}_{12}\left(\dot{F}\right)^{-1}$$

This, the rate constant(s) are given by:

$$k_{j}^{+} = \left(\frac{\dot{F}}{\dot{n}_{12}}\right)\hat{r}_{j}\qquad{j=1,2}$$

"""

# ╔═╡ 86b16793-3647-4bd6-870e-3ca2be3035ae
begin

	# concetration -
	A_concentration = (n_dot_out[1]/F_dot)

	# compute -
	k1f = r1_hat*(1/A_concentration)
	k2f = r2_hat*(1/A_concentration)

	# display -
	with_terminal() do
		println("k1f = $(k1f) 1/hr and k2f = $(k2f) 1/hr")
	end
end

# ╔═╡ 4bab103c-3d4b-11ec-1ceb-9566ca26ce6b
html"""
<style>
main {
    max-width: 1200px;
    width: 85%;
    margin: auto;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"

[compat]
PlutoUI = "~0.7.19"
PrettyTables = "~1.2.3"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "0bc60e3006ad95b4bb7497698dd7c6d649b9bc06"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "e071adf21e165ea0d904b595544a8e514c8bb42c"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.19"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "d940010be611ee9d67064fe559edbb305f8cc0eb"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.2.3"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─bd0378cb-ca70-4803-b5fc-7bd08ad1b9f3
# ╠═6b07d7ae-8b5d-4c41-8184-9329cdaa94a6
# ╟─7e8aaae7-555f-4b70-bae7-4adc1d1c9923
# ╠═fddd2b82-8eba-4517-bdfd-0953bd10c259
# ╠═86964693-9968-4524-a2bc-8c3a8f202cc1
# ╟─34ddef24-462a-4a1d-8a61-9c760f278883
# ╟─fff53850-a627-4349-bef2-c5aaf8828c3d
# ╟─8b5a451f-5d31-49f5-91fb-3e24a00b87ac
# ╠═86b16793-3647-4bd6-870e-3ca2be3035ae
# ╠═dc9d1ef7-23e2-4bc3-ab47-e9ca80af5811
# ╟─4bab103c-3d4b-11ec-1ceb-9566ca26ce6b
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
include("cont_2d_worlds.jl")

module FourRoomsContParams

const BASE_WALLS = [0 0 0 0 0 1 0 0 0 0 0;
                    0 0 0 0 0 1 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 1 0 0 0 0 0;
                    0 0 0 0 0 1 0 0 0 0 0;
                    1 0 1 1 1 1 0 0 0 0 0;
                    0 0 0 0 0 1 1 1 0 1 1;
                    0 0 0 0 0 1 0 0 0 0 0;
                    0 0 0 0 0 1 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 1 0 0 0 0 0;]

const GOAL_LOCS =  [0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;
                    0 0 0 0 0 0 0 0 0 0 0;]

const ROOM_TOP_LEFT = 1
const ROOM_TOP_RIGHT = 2
const ROOM_BOTTOM_LEFT = 3
const ROOM_BOTTOM_RIGHT = 4

const REW_FUNCS = Dict{Int, Function}()

end

FourRooms(max_action_noise=0.1,
          drift_noise=0.001;
          normalized=false) =
              ContGridWorld(FourRoomsContParams.BASE_WALLS[end:-1:1, :],
                            FourRoomsContParams.GOAL_LOCS[end:-1:1, :],
                            FourRoomsContParams.REW_FUNCS,
                            max_action_noise,
                            drift_noise,
                            normalized)

module TMazeContParams

import ..BasicRewFunc

const BASE_WALLS = [1 1 1 1 1 1 1 1 1 1 1;
                    1 0 0 1 1 1 1 1 0 0 1;
                    1 0 0 0 0 0 0 0 0 0 1;
                    1 0 0 0 0 0 0 0 0 0 1;
                    1 0 0 0 0 0 0 0 0 0 1;
                    1 0 0 1 0 0 0 1 0 0 1;
                    1 1 1 1 0 0 0 1 1 1 1;
                    1 1 1 1 0 0 0 1 1 1 1;
                    1 1 1 1 0 0 0 1 1 1 1;
                    1 1 1 1 0 0 0 1 1 1 1;
                    1 1 1 1 0 0 0 1 1 1 1;]

const GOAL_LOCS = [-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1;
                   -1  1  1 -1 -1 -1 -1 -1  3  3 -1;
                   -1  0  0  0  0  0  0  0  0  0 -1;
                   -1  0  0  0  0  0  0  0  0  0 -1;
                   -1  0  0  0  0  0  0  0  0  0 -1;
                   -1  2  2 -1  0  0  0 -1  4  4 -1;
                   -1 -1 -1 -1  0  0  0 -1 -1 -1 -1;
                   -1 -1 -1 -1  0  0  0 -1 -1 -1 -1;
                   -1 -1 -1 -1  0  0  0 -1 -1 -1 -1;
                   -1 -1 -1 -1  0  0  0 -1 -1 -1 -1;
                   -1 -1 -1 -1  0  0  0 -1 -1 -1 -1;]

const REW_FUNCS = Dict([i=>BasicRewFunc(i) for i in 1:4]...)

end

TMaze(max_action_noise=0.1, drift_noise=0.001; normalized=false) =
    ContGridWorld(TMazeContParams.BASE_WALLS[end:-1:1, :],
                  TMazeContParams.GOAL_LOCS[end:-1:1, :],
                  TMazeContParams.REW_FUNCS,
                  max_action_noise,
                  drift_noise,
                  normalized)

module OpenWorldContParams

import ..BasicRewFunc
import ..Distributions

const REW_FUNCS = Dict([i=>BasicRewFunc(i) for i in 1:4]...)

function center_start_func(env, rng)
    sze = size(env.walls)
    x_rng = Distributions.Uniform((round(sze[1]/2)) - 0.5, round(sze[1]/2) + 0.5)
    y_rng = Distributions.Uniform((round(sze[2]/2)) - 0.5, round(sze[2]/2) + 0.5)

    x = rand(rng, x_rng)
    y = rand(rng, y_rng)
    [x, y]
end

end

OpenWorld(width, height, max_action_noise=0.1, drift_noise=0.001; per_step_rew=0.0, cumulant_schedule=nothing, normalized=true, start_type=:none) = begin
    base_walls = zeros(Int, width, height)
    goal_locs = zeros(Int, width, height)
    goal_locs[1, 1] = 1
    goal_locs[1, end] = 3
    goal_locs[end, 1] = 2
    goal_locs[end, end] = 4

    start = if start_type == :uniform
        nothing
    elseif start_type == :center
        OpenWorldContParams.center_start_func
    else
        throw("Only starting dist availble is center (center area) or uniform (uniform random)")
    end

    ContGridWorld(base_walls,
                  goal_locs[end:-1:1, :],
                  cumulant_schedule,
                  start,
                  max_action_noise,
                  drift_noise,
                  normalized,
                  per_step_rew)
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
## Data structure
"""
    struct EDDataDescriptor

    EDDataDescriptor(model[, diagtype][, pathdata]; suffix, prefix)

The important bits of information needed to specify for performing exact diagonalization on spin models:
 - `model`: See [`EDModels`](@ref)
 - `diagtype`: See [`DiagonalizationType`](@ref)
 - `pathdata`: `SaveLocation` for the results
For `load`ing data only the first field is required!

Note: This datastructure forwards properties of `model`.
"""
struct EDDataDescriptor <: SimLib.AbstractDataDescriptor
    model::EDModels.Model
    diagtype::Maybe{DiagonalizationType}
    pathdata::SaveLocation
    # function EDDataDescriptor(model, diagtype, pathdata::SaveLocation)
    #     geometry ∈ SimLib.GEOMETRIES || error("Unknown geometry: $geom")
    #     scale_fields ∈ [:none, :ensemble, :shot]
    #     new(geometry, dimension, system_size, α, shots, ismissing(ρs) ? missing : unique!(sort(collect(ρs))), ismissing(fields) ? missing : unique!(sort(collect(fields))), scale_fields, basis, diagtype, pathdata::SaveLocation)
    # end
    function EDDataDescriptor(model, diagtype=missing, pathdata=SaveLocation(); suffix=pathdata.suffix, prefix=pathdata.prefix)
        new(model, diagtype, SaveLocation(;suffix, prefix))
    end
end

# constructor for easier loading
EDDataDescriptor(model, pathdata::SaveLocation; kwargs...) = EDDataDescriptor(model, missing, pathdata; kwargs...)
EDDataDescriptor(eddd::EDDataDescriptor) = eddd

# forward properties to model
function Base.getproperty(edd::EDDataDescriptor, s::Symbol)
    if hasfield(typeof(edd), s)
        return getfield(edd, s)
    else
        return getproperty(getfield(edd, :model), s)
    end
end

function Base.:(==)(d1::EDDataDescriptor, d2::EDDataDescriptor)
    all(getfield(d1, f) == getfield(d2, f) for f in [:model, :diagtype])
end


text:
function test_sparse_vanilla_same(m)
    vanilla_solver = ValueIterationSolver(max_iterations=500, belres=1e-8)
    sparse_solver = SparseValueIterationSolver(max_iterations=500, belres=1e-8)
    vanilla_policy = solve(vanilla_solver, m)
    sparse_policy = solve(sparse_solver, m)
    @test vanilla_policy.policy == sparse_policy.policy
    @test isapprox(vanilla_policy.util, sparse_policy.util; atol=1e-3)
    for s in states(m)
        si = stateindex(m, s)
        for a in actions(m, s)
            ai = actionindex(m, a)
            @test isapprox(vanilla_policy.qmat[si, ai], sparse_policy.qmat[si, ai]; atol=1e-3)
        end
    end
end

m1 = SimpleGridWorld(size=(10,10))
test_sparse_vanilla_same(m1)

m2 = SpecialGridWorld()
test_sparse_vanilla_same(m2)


text:
using SumProductNetworks

# Load an automatic differentiation (AD) package, this case we use ForwardDiff.
# Note that using a reverse mode AD is more efficient.
using ForwardDiff

# Loading some utility functions.
using Random, LinearAlgebra

# Generate some training data.
N = 1000
D = 4
x = rand(1:2, N, D)

# Helper function to construct an SPN structure.
function buildSPN(T::Type{<:Real}, x; K = 2)
    N, D = size(x)
    root = FiniteSumNode{T}(;D = D);

    for k in 1:K
        add!(root, FiniteProductNode(D = D), map(T, log(1.0/K)))
        for d in 1:D
            add!(root[k], FiniteSumNode{T}(D = D))
            for l in 1:2
                add!(root[k][d], IndicatorNode(l, d), map(T, log(0.5)))
            end
        end
    end
    updatescope!(root)
    return SumProductNetwork(root)
end

# Number of children under a sum node.
K = 2

# Construct the SPN.
spn = buildSPN(Real, x, K = K)

# Collect all sum nodes.
snodes = filter(n -> isa(n, SumNode), values(spn))

# Helper function used by ForwardDiff.
function f(θ)

    N, D = size(x)
    c = 1
    for i in 1:length(snodes)
        K = length(snodes[i])

        # Extraction of parameters.
        ϕ = θ[c:(c+K-1)]

        # Copy parameters to SPN node.
        snodes[i].logweights[:] = log.(projectToPositiveSimplex!(ϕ))

        c += length(snodes[i])
    end

    # Return llh.
    return -mean(logpdf(spn.root, x))
end

# Create an initial guess for all parameters and project them to the positive simplex.
q = mapreduce(n -> projectToPositiveSimplex!(rand(length(n))), vcat, snodes)

@info "Number of parameters to optimize: $(length(q))"

# Configure the AD.
chunk = ForwardDiff.Chunk(10)
config = ForwardDiff.GradientConfig(f, q, chunk)

# Construct AD function.
g(θ) = ForwardDiff.gradient!(similar(θ), f, θ, config)

# Number of iterations.
τ = 100

# log likelihood values.
ℓ = zeros(τ)

@info "Starting optimization for $(τ) iterations..."
@info "Initial log likelihood: $(-f(q))"

for i in 1:τ

    # Save LLH.
    ℓ[i] = -f(q)
    # Gradient step.
    δq = g(q)
    q .-= 0.1 * δq
end

@info "Log likelihood of optimized model: $(last(ℓ))"


text:
module TestMeasures

# using Revise
using Test
using MLJBase
import Distributions
using CategoricalArrays
import Random.seed!
seed!(1234)

@testset "built-in regressor measures" begin

    y    = [1, 2, 3, 4]
    yhat = [4, 3, 2, 1]
    w = [1, 2, 4, 3]
    @test isapprox(mav(yhat, y), 2)
    @test isapprox(mav(yhat, y, w), 9/5)
    @test isapprox(rms(yhat, y), sqrt(5))
    @test isapprox(rms(yhat, y, w), sqrt(21/5))
    @test isapprox(mean(l1(yhat, y)), 2)
    @test isapprox(mean(l1(yhat, y, w)), 9/5)
    @test isapprox(mean(l2(yhat, y)), 5)
    @test isapprox(mean(l2(yhat, y, w)), 21/5)

    yhat = y .+ 1
    @test isapprox(rmsl(yhat, y),
                   sqrt((log(1/2)^2 + log(2/3)^2 + log(3/4)^2 + log(4/5)^2)/4))
    @test isapprox(rmslp1(yhat, y),
                   sqrt((log(2/3)^2 + log(3/4)^2 + log(4/5)^2 + log(5/6)^2)/4))
    @test isapprox(rmsp(yhat, y), sqrt((1 + 1/4 + 1/9 + 1/16)/4))

end

@testset "built-in classifier measures" begin

    y    = categorical(collect("asdfasdfaaassdd"))
    yhat = categorical(collect("asdfaadfaasssdf"))
    w = 1:15
    @test misclassification_rate(yhat, y) ≈ 0.2
    @test misclassification_rate(yhat, y, w) ≈ 4/15
    y = categorical(collect("abb"))
    L = [y[1], y[2]]
    d1 = UnivariateFinite(L, [0.1, 0.9])
    d2 = UnivariateFinite(L, [0.4, 0.6])
    d3 = UnivariateFinite(L, [0.2, 0.8])
    yhat = [d1, d2, d3]
    @test mean(cross_entropy(yhat, y)) ≈ -(log(0.1) + log(0.6) + log(0.8))/3
    scores = BrierScore()(yhat, y)
    @test scores ≈ [-1.62, -0.32, -0.08]
    wscores = BrierScore()(yhat, y, [1, 2, 7])
    @test wscores ≈ scores .* [0.3, 0.6, 2.1]

end

@testset "MLJBase.value" begin

    yhat = rand(5)
    X = (weight=rand(5), x1 = rand(5))
    y = rand(5)
    w = rand(5)

    @test MLJBase.value(mav, yhat, nothing, y, nothing) ≈ mav(yhat, y)
    @test MLJBase.value(mav, yhat, nothing, y, w) ≈ mav(yhat, y, w)

    spooky(yhat, y) = abs.(yhat - y) |> mean
    @test MLJBase.value(spooky, yhat, nothing, y, nothing) ≈ mav(yhat, y)

    cool(yhat, y, w) = abs.(yhat - y) .* w ./ (sum(w)/length(y)) |> mean
    MLJBase.supports_weights(::Type{typeof(cool)}) = true
    @test MLJBase.value(cool, yhat, nothing, y, w) ≈ mav(yhat, y, w)

    funky(yhat, X, y) = X.weight .* abs.(yhat - y) ./ (sum(X.weight)/length(y)) |> mean
    MLJBase.is_feature_dependent(::Type{typeof(funky)}) = true
    @test MLJBase.value(funky, yhat, X, y, nothing) ≈ mav(yhat, y, X.weight)

    weird(yhat, X, y, w) = w .* X.weight .* abs.(yhat - y) ./ sum(w .* X.weight) |> sum
    MLJBase.is_feature_dependent(::Type{typeof(weird)}) = true
    MLJBase.supports_weights(::Type{typeof(weird)}) = true
    @test MLJBase.value(weird, yhat, X, y, w) ≈ mav(yhat, y, X.weight .* w)

end

# for when ROC is added as working dependency:
# y = ["n", "p", "n", "p", "n", "p"]
# yhat = [0.1, 0.2, 0.3, 0.6, 0.7, 0.8]
# @test auc("p")(yhat, y) ≈ 2/3

end
true


text:
using DiffEqPhysics
using StochasticDiffEq

T = 120.0 # °K
T0 = 90
kb = 8.3144598e-3 # kJ/(K*mol)
ϵ = T * kb
σ = 0.34 # nm
ρ = 1374/1.6747# Da/nm^3
m = 39.95# Da
N = 216
L = (m*N/ρ)^(1/3)#10.229σ
R = 0.5*L
v_dev = sqrt(kb * T / m)
bodies = generate_bodies_in_cell_nodes(N, m, v_dev, L)

τ = 0.5e-3 # ps or 1e-12 s
t1 = 0.0
t2 = 10000τ

parameters = LennardJonesParameters(ϵ, σ, R)
lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => parameters));
#thermostat = AndersenThermostat(T0, 1000)
#thermostat = BerendsenThermostat(T0, 20τ)
thermostat = LangevinThermostat(T0, 10.0)
pbc = CubicPeriodicBoundaryConditions(L)
simulation = NBodySimulation(lj_system, (t1, t2), pbc, thermostat, kb);
#result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)
result = @time run_simulation_sde(simulation, EM(),  dt=τ)

time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
Nactual = length(bodies)
timesteps = round(length(result.solution.t))

using Plots
t = t1:τ:result.solution.t[end-1]
temper = temperature.(result, t)
pl=plot(t, temper, ylim=[0,500], xlabel="t, ps", ylabel = "T, °K", label="Temperature, °K", linewidth=2)
plot!(pl, t, T0*ones(length(t)), label = "$T0 °K", linewidth=2)
plot!(pl, title="Langevin thermostat at gamma=$(thermostat.γ)  1/ps")

using JLD
save("d:/$(typeof(thermostat)) thermostat for water $T0 _$time_now.jld", "t",t, "temper", temper, "T0", T0, "thermostat", thermostat)


text:
#__precompile__()

module HDR

using Images, ImageView, Plots, ImageFiltering, FileIO

# Currently N0f16 and others are NOT supported
HDRSequence = Vector{Tuple{Matrix{RGB{N0f8}}, Float64}}

include("tonemapper.jl")
include("synthesizer.jl")
include("reader.jl")

include("test_synthesizer.jl")

export read_hdrgen_file
export image_synthesis
export hdr_sequence_downsample
export pseudo_color_visualizer

end # module HDR


text:
# ---------  GradRule simple
# Graduate after given `t` and with at least `hGradV`.

Base.@kwdef mutable struct GradRuleSwitchesSimple <: GradRuleSwitches
    "Students can graduate at the END of this year or thereafter"
    tFirstGrad :: TimeInt = 4
    "Number of courses required for graduation (in DATA courses)"
    minNcForGrad :: ncInt = 40
    "Number of 2 year colleges (cannot graduate from those)"
    n2Year :: CollInt = 1
    "Number of colleges"
    nColleges :: CollInt = 4
    # Min and max grad prob when graduation is feasible
    gradProbMin :: Double = 0.05
    gradProbMax :: Double = 0.95
    # Grad prob depends on levels or gains
    useHcLevel :: Bool = true
    byCollege :: Bool = false
end

mutable struct GradRuleSetSimple <: GradRuleSet
    @grad_rule_set_common
    # Min h to graduate for each college
    hGradV :: IncreasingVector{Double}
    switches :: GradRuleSwitchesSimple
end

## For one college
struct GradRuleSimple <: GradRule
    # Can graduate after at least this many years (at the END of this year)
    tFirstGrad :: TimeInt
    # Min h to graduate
    hGrad :: Double
    # Min no of courses attempted (cumulative) (at the END of the year)
    nMin :: ncInt
    # Min and max grad prob when graduation is feasible
    gradProbMin :: Double
    gradProbMax :: Double
    # Grad prob depends on levels or gains
    useHcLevel :: Bool
end

## ----------- Simple GradRule

# There is no hMin property here. Ignore.
function set_hmin!(switches :: GradRuleSwitchesSimple, hMin :: Double) end

ModelParams.has_pvector(gs :: GradRuleSetSimple) = false;

# function make_test_gradrule(gs :: GradRuleSwitchesSimple)
#     return GradRuleSimple(CollInt(4), 2.0, ncInt(20),
#         grad_prob_min(gs), grad_prob_max(gs), use_hc_level(gs))
# end

"""
	$(SIGNATURES)

For all colleges.
Assumes college 1 is a 2 year college. Others are 4 year colleges.
"""
function make_grad_rule_set(objId :: ObjectId,
    grSwitches :: GradRuleSwitchesSimple)

    nc = grSwitches.nColleges;
    hGradV = init_hgrad(objId, grSwitches);

    pvec = ParamVector(objId,  []);
    return GradRuleSetSimple(objId, pvec, hGradV, grSwitches)
end

"""
	$(SIGNATURES)

Make grad rule for one college
"""
function make_grad_rule(gs :: GradRuleSetSimple, iCollege :: Integer)
    return GradRuleSimple(t_first_grad(gs, iCollege),
        h_grad(gs, iCollege),  min_courses_to_grad(gs, iCollege),
        grad_prob_min(gs), grad_prob_max(gs),
        use_hc_level(gs))
end

## Graduation probability at end of t
# One state
function grad_prob(g :: GradRuleSimple, t :: Integer, h :: F1, n :: I1,
    h0 :: F1) where {I1 <: Integer, F1 <: AbstractFloat}

    if can_graduate(g, t)
        if (h >= g.hGrad) && (n >= g.nMin)
            gProb = grad_prob_max(g);
        else
            gProb = grad_prob_min(g);
        end
    else
        gProb = zero(F1);
    end
    return gProb
end

# --------------


text:
module G1
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    1.8558518988717e-6*l1.^5 + 9.70087471350808e-6*l1.^4.*l2 - 1.89926015002423e-5*l1.^4.*l3 - 5.70920347683978e-8*l1.^4 + 8.93480174177124e-6*l1.^3.*l2.^2 - 5.67541019579733e-5*l1.^3.*l2.*l3 - 4.60335359066924e-8*l1.^3.*l2 + 6.64276883744224e-5*l1.^3.*l3.^2 + 2.74798041182713e-7*l1.^3.*l3 - 2.44419795743875e-10*l1.^3 - 1.56450323689193e-5*l1.^2.*l2.^3 + 1.99678170677335e-5*l1.^2.*l2.^2.*l3 + 3.09604520022297e-7*l1.^2.*l2.^2 + 6.54477101135341e-5*l1.^2.*l2.*l3.^2 - 4.80212518198379e-7*l1.^2.*l2.*l3 - 4.07104029469703e-10*l1.^2.*l2 - 8.84038175553207e-5*l1.^2.*l3.^3 - 1.7313330904159e-7*l1.^2.*l3.^2 + 1.14159743254739e-9*l1.^2.*l3 - 3.20459750258103e-13*l1.^2 - 1.62995586742055e-5*l1.*l2.^4 + 9.63865835476286e-5*l1.*l2.^3.*l3 - 5.76159669483844e-8*l1.*l2.^3 - 0.000164232243744917*l1.*l2.^2.*l3.^2 - 4.47135751195212e-7*l1.*l2.^2.*l3 + 7.88991063020458e-10*l1.*l2.^2 + 6.5456473828544e-5*l1.*l2.*l3.^3 + 9.27226967263222e-7*l1.*l2.*l3.^2 - 7.60394380347721e-10*l1.*l2.*l3 - 7.1886514515434e-13*l1.*l2 + 2.80178737068122e-5*l1.*l3.^4 - 1.9291801143875e-7*l1.*l3.^3 - 7.64342748517079e-10*l1.*l3.^2 + 1.36072328692881e-12*l1.*l3 + 1.77973142550235e-5*l1 - 3.71608982966903e-6*l2.^5 + 3.4860347811527e-5*l2.^4.*l3 - 8.99561440808747e-8*l2.^4 - 0.000117823809431305*l2.^3.*l3.^2 + 4.16948693822733e-7*l2.^3.*l3 + 3.53672039292904e-11*l2.^3 + 0.000172372563165961*l2.^2.*l3.^3 - 4.0073008496488e-7*l2.^2.*l3.^2 - 8.95304236140046e-10*l2.^2.*l3 + 4.99582038836009e-13*l2.^2 - 0.000102352765801864*l2.*l3.^4 - 4.26318851821059e-8*l2.*l3.^3 + 1.27404682191573e-9*l2.*l3.^2 - 2.77898404065488e-13*l2.*l3 - 2.25988265609486e-5*l2 + 1.47914348941021e-5*l3.^5 + 5.88810194343027e-8*l3.^4 - 1.68437205927552e-10*l3.^3 - 5.43101594146752e-13*l3.^2 + 4.6819529309252e-6*l3 + 1.83176148000071e-8
end
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.8.2
#   kernelspec:
#     display_name: Julia 1.6.1
#     language: julia
#     name: julia-1.6
# ---

using Pkg
Pkg.add(["Combinatorics", "Plots"])

using Combinatorics
using LinearAlgebra
using Plots

# +
function polynomial_features(X, degree)

    n_samples, n_features = size(X)

    combs = [with_replacement_combinations(1:n_features, i) for i in 0:degree+1]

    combinations = []
    for sublist in combs
        for j in sublist
             push!(combinations, j)
        end
    end
    n_output_features = length(combinations)
    X_new = zeros(n_samples, n_output_features)

    for (i, index_combs) in enumerate(combinations)
        X_new[:, i] = prod(X[:, index_combs], dims=2)
    end

    return X_new

end
# -

# # Regression polynomial non locale

# +
f(x, y) = - sin(x) .+ 0.5 * cos(y)
n, p = 1000, 2
degree = 3
X = randn(n, p)
Y_true = f.(X[:,1],X[:,2])
Y = Y_true .+ 0.2 .* randn(n)
X_new = polynomial_features(X, degree)
beta = X_new \ Y

X_grid = hcat(vcat([0.1i for i in -20:20, j in -20:20]...),
              vcat([0.1j for i in -20:20, j in -20:20]...))
Y_true = f.(X_grid[:,1], X_grid[:,2])
X_grid_new = polynomial_features( X_grid, degree)
Y_grid = X_grid_new * beta
plot(X_grid[:,1], X_grid[:,2], Y_grid, st=:surface, alpha=0.5 )
scatter!(X[:,1], X[:,2], Y, ms = 1)

# +

scatter( X_grid[:,1], X_grid[:,2], Y_true)
scatter!(X_grid[:,1], X_grid[:,2], Y_grid)

# +
abstract type AbstractKernel end

struct Epanechnikov <: AbstractKernel end

"""
    epanechnikov(t)

Compute the Epanechnikov kernel.
"""
function (::Epanechnikov)(t)

    kernel = zero(t)
    idx = findall(abs.(t) .<= 1)
    kernel[idx] = 0.75 .* (1 .- t[idx].^2)
    return kernel

end

# -

# # Regression locale

# +
kernel = Epanechnikov()
width = 1

Y_pred = Float64[]
for x in eachrow(X_grid)
    dist = sqrt.(sum((X .- x').^2, dims=2)) ./ width
    w = vec(kernel(dist))
    idx = findall(abs.(w) .> 1e-10)
    beta = pinv(X_new[idx,:]' * Diagonal(w[idx]) * X_new[idx,:]) * (X_new[idx,:]' * Diagonal(w[idx]) * Y[idx])
    P = polynomial_features(x', degree)

    push!(Y_pred, P * beta...)
end

scatter( X_grid[:,1], X_grid[:,2], Y_pred)
# -