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

using Markdown
using InteractiveUtils

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

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

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

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
trueVar, trueStd = 1/12, sqrt(1/12)

# ╔═╡ 81b26bd0-3a9a-11eb-0a09-db9eed1b0609
function estVar(n)
    sample = rand(n)
    sum((sample .- 0.5).^2)/n
end

# ╔═╡ 81b2aba4-3a9a-11eb-1bd6-a938171cf060
begin
	N = 10^7
	df = DataFrame()
	for n in 5:5:30
		biasVar = mean([estVar(n) for _ in 1:N]) - trueVar
		biasStd = mean([sqrt(estVar(n)) for _ in 1:N]) - trueStd
		append!(df, DataFrame(:n => n, :var_bias => round(biasVar, digits=5),
				:std_bias => round(biasStd, digits=5)))
	end
	df
end

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

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═81b26bd0-3a9a-11eb-0a09-db9eed1b0609
# ╠═81b2aba4-3a9a-11eb-1bd6-a938171cf060
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
"""
Module Internet

Provides functions for navigating in the internet by voice.

# Functions

###### Searching
- [`search`](@ref)

To see a description of a function type `?<functionname>`.

See also: [`Email`](@ref)
"""
module Internet

using PyCall
import DefaultApplication
import ..JustSayIt: @voiceargs, Keyboard, Key

## CONSTANTS

const SEARCH_ENGINE = "https://google.com"
const SEARCH_DOC = "start search engine, enter \"type text\" mode to obtain search words from speech and then trigger the search when the keyword \"search\" is spoken."

## FUNCTIONS

@doc """
    internet `search`

Navigate the internet, performing one of the following actions:
- `search`: $SEARCH_DOC
"""
internet
@enum Action search
@voiceargs action=>(valid_input_auto=true) function internet(action::Action)
    if (action == search) search_internet()
    else @info "unknown action" #NOTE: this should never happen.
    end
end

@doc SEARCH_DOC
function search_internet()
    DefaultApplication.open(SEARCH_ENGINE)
    Keyboard.type(Keyboard.text; end_keyword="search")
    Keyboard.press_keys(Key.enter)
end

end # module Internet

ext: jl


text:
module LLRBVisualize

using LLRBTrees
import Base: replace, display

export buildstring, draw

include("GraphViz.jl")

function buildstring{K,V}(tree::LLRBTree{K,V})

    treestring = """
    digraph tree {

        }
        """
    linenow = 16

    if(isdefined(tree, :root))
        root=tree.root
        if(!isleaf(root))

            if(!isleftleaf(root))
                treestring, linenow = buildstringR(treestring, root.left, root, linenow)
            else
                treestring, linenow = leafstring(treestring, root, linenow)
            end

            if(!isrightleaf(root))
                treestring, linenow = buildstringR(treestring, root.right, root, linenow)
            else
                treestring, linenow = leafstring(treestring, root, linenow)
            end

        end

    end

    return treestring
end

function leafstring{K,V}(treestring::ASCIIString, parent::TreeNode{K,V}, linenow::Int)

    linkstring = """null$(linenow) [shape=point];
    $(parent.value.value) -> null$(linenow);
    """
    treestring = replace(treestring, linenow, linenow+1, linkstring)
    linenow += length(linkstring)

    return treestring, linenow
end

function buildstringR{K,V}(treestring::ASCIIString, node::TreeNode{K,V}, parent::TreeNode{K,V}, linenow::Int)

    nodestring = string(node.value.value)*";\n"
    linkstring = string(parent.value.value)*" -> "*nodestring
    if(node.isRed)
        len=length(linkstring)
        linkstring = replace(linkstring, len-2, len-1, " [color=red]")
    end
    treestring = replace(treestring, linenow, linenow+1, linkstring)
    linenow += length(linkstring)

    if(!isleftleaf(node))
        treestring, linenow = buildstringR(treestring, node.left, node, linenow)
    else
        treestring, linenow = leafstring(treestring, node, linenow)
    end

    if(!isrightleaf(node))
        treestring, linenow = buildstringR(treestring, node.right, node, linenow)
    else
        treestring, linenow = leafstring(treestring, node, linenow)
    end

    return treestring, linenow
end

function display{K,V}(tree::LLRBTree{K,V})
    GraphViz.Graph(buildstring(tree))
end

#LLRBVisualize end
end


text:
# Unit testing for MxArray

using MATLAB, Compat, Base.Test
using MATLAB: nfields

m = 5
n = 6

# test basic types in 1D & 2D

macro mx_test_basic_types(ty, testfun)
    quote
        a = mxarray($(ty), n)
        @test elsize(a) == sizeof($(ty))
        @test eltype(a) === $(ty)
        @test nrows(a) == n
        @test ncols(a) == 1
        @test nelems(a) == n
        @test ndims(a) == 2
        @test size(a) == (n, 1)
        @test size(a, 1) == n
        @test size(a, 2) == 1
        @test size(a, 3) == 1
        @test !is_complex(a)
        @test $(testfun)(a)
        delete(a)

        b = mxarray($(ty), m, n)
        @test elsize(b) == sizeof($(ty))
        @test eltype(b) === $(ty)
        @test nrows(b) == m
        @test ncols(b) == n
        @test nelems(b) == m * n
        @test ndims(b) == 2
        @test size(b) == (m, n)
        @test size(b, 1) == m
        @test size(b, 2) == n
        @test size(b, 3) == 1
        @test !is_complex(b)
        @test $(testfun)(b)
        delete(b)
    end
end

# empty array

a = mxempty()
@test nrows(a) == 0
@test ncols(a) == 0
@test nelems(a) == 0
@test ndims(a) == 2
@test eltype(a) == Float64
@test is_empty(a)

# basic arrays

@mx_test_basic_types Float64 is_double
@mx_test_basic_types Float32 is_single
@mx_test_basic_types Int64   is_int64
@mx_test_basic_types Uint64  is_uint64
@mx_test_basic_types Int32   is_int32
@mx_test_basic_types Uint32  is_uint32
@mx_test_basic_types Int16   is_int16
@mx_test_basic_types Uint16  is_uint16
@mx_test_basic_types Int8    is_int8
@mx_test_basic_types Uint8   is_uint8
@mx_test_basic_types Bool    is_logical

# complex arrays

macro mx_test_complex_type(ty, testfun)
    quote
        b = mxarray(Complex{$(ty)}, m, n)
        @test elsize(b) == sizeof($(ty))
        @test eltype(b) === $(ty)
        @test nrows(b) == m
        @test ncols(b) == n
        @test nelems(b) == m * n
        @test ndims(b) == 2
        @test size(b) == (m, n)
        @test size(b, 1) == m
        @test size(b, 2) == n
        @test size(b, 3) == 1
        @test is_complex(b)
        @test $(testfun)(b)
        delete(b)
    end
end
@mx_test_complex_type Float64 is_double
@mx_test_complex_type Float32 is_single

# test creating multi-dimensional arrays

a = mxarray(Float64, (6, 5, 4))
@test elsize(a) == sizeof(Float64)
@test eltype(a) === Float64
@test size(a) == (6, 5, 4)
@test size(a, 1) == 6
@test size(a, 2) == 5
@test size(a, 3) == 4
@test size(a, 4) == 1
@test nelems(a) == 6 * 5 * 4
@test is_numeric(a)
@test !is_sparse(a)

a = mxarray(Bool, (6, 5, 4))
@test elsize(a) == 1
@test eltype(a) === Bool
@test size(a) == (6, 5, 4)
@test size(a, 1) == 6
@test size(a, 2) == 5
@test size(a, 3) == 4
@test size(a, 4) == 1
@test nelems(a) == 6 * 5 * 4
@test is_logical(a)
@test !is_sparse(a)

# scalars

a_mx = mxarray(3.25)
@test eltype(a_mx) == Float64
@test size(a_mx) == (1, 1)
@test jscalar(a_mx) == 3.25
delete(a_mx)

a_mx = mxarray(@compat Int32(12))
@test eltype(a_mx) == Int32
@test size(a_mx) == (1, 1)
@test jscalar(a_mx) == @compat Int32(12)
delete(a_mx)

a_mx = mxarray(true)
@test eltype(a_mx) == Bool
@test size(a_mx) == (1, 1)
@test jscalar(a_mx)
delete(a_mx)

a_mx = mxarray(false)
@test eltype(a_mx) == Bool
@test size(a_mx) == (1, 1)
@test !jscalar(a_mx)
delete(a_mx)

a_mx = mxarray(3.25 + 4im)
@test eltype(a_mx) == Float64
@test size(a_mx) == (1, 1)
@test jscalar(a_mx) == 3.25 + 4im
delete(a_mx)

# conversion between Julia and MATLAB numeric arrays

a = rand(5, 6)
a_mx = mxarray(a)
a2 = jarray(a_mx)
@test isequal(a, a2)
delete(a_mx)

a = rand(5)
a_mx = mxarray(a)
a2 = jvector(a_mx)
@test isequal(a, a2)
delete(a_mx)

a_t = reshape(a, 1, 5)
a_mx = mxarray(a_t)
a2 = jvector(a_mx)
@test isequal(a, a2)
delete(a_mx)

a = 1:5
a_mx = mxarray(a)
a2 = jvector(a_mx)
@test isequal([1:5;], a2)
delete(a_mx)

a = rand(5, 6) + rand(5, 6)*im
a_mx = mxarray(a)
a2 = jarray(a_mx)
@test isequal(a, a2)
delete(a_mx)

# sparse matrices

a = sprand(8, 9, 0.2)
a_mx = mxarray(a)
@test is_double(a_mx)
@test is_sparse(a_mx)
@test nrows(a_mx) == 8
@test ncols(a_mx) == 9

a2 = jsparse(a_mx)
@test size(a2) == (8, 9)
@test countnz(a2) == countnz(a)
@test isequal(a2, a)
delete(a_mx)

a = sparse(convert(Array{Bool}, rand(8, 9) .< 0.3))
a_mx = mxarray(a)
@test is_logical(a_mx)
@test is_sparse(a_mx)
@test nrows(a_mx) == 8
@test ncols(a_mx) == 9

a2 = jsparse(a_mx)
@test size(a2) == (8, 9)
@test countnz(a2) == countnz(a)
@test isequal(a2, a)
delete(a_mx)

# strings

s = "MATLAB.jl"
s_mx = mxarray(s)
@test classid(s_mx) == MATLAB.mxCHAR_CLASS
@test nrows(s_mx) == 1
@test ncols(s_mx) == length(s)
@test nelems(s_mx) == length(s)
@test ndims(s_mx) == 2
@test is_char(s_mx)

s2 = jstring(s_mx)
@test s == s2
delete(s_mx)

s = ""
s_mx = mxarray(s)
@test classid(s_mx) == MATLAB.mxCHAR_CLASS
@test is_char(s_mx)
@test is_empty(s_mx)

s2 = jstring(s_mx)
@test s == s2
delete(s_mx)

# cell arrays

a = mxcellarray(10)
@test nrows(a) == 10
@test ncols(a) == 1
@test nelems(a) == 10
@test classid(a) == MATLAB.mxCELL_CLASS
@test is_cell(a)
delete(a)

a = mxcellarray(4, 5)
@test nrows(a) == 4
@test ncols(a) == 5
@test nelems(a) == 20
@test classid(a) == MATLAB.mxCELL_CLASS
@test is_cell(a)
delete(a)

a = mxcellarray((3, 4, 5))
@test size(a) == (3, 4, 5)
@test nelems(a) == 60
@test classid(a) == MATLAB.mxCELL_CLASS
@test is_cell(a)
delete(a)

s = ["abc", "efg"]
s_mx = mxcellarray(s)
@test jstring(get_cell(s_mx, 1)) == "abc"
@test jstring(get_cell(s_mx, 2)) == "efg"
delete(s_mx)

# struct

a = mxstruct("abc", "efg", "xyz")
@test is_struct(a)
@test nfields(a) == 3
@test nrows(a) == 1
@test ncols(a) == 1
@test nelems(a) == 1
@test ndims(a) == 2

@test get_fieldname(a, 1) == "abc"
@test get_fieldname(a, 2) == "efg"
@test get_fieldname(a, 3) == "xyz"
delete(a)

s = @compat Dict{Any,Any}("name"=>"MATLAB", "version"=>12.0, "data"=>[1,2,3])
a = mxstruct(s)
@test is_struct(a)
@test nfields(a) == 3
@test jstring(get_field(a, "name")) == "MATLAB"
@test jscalar(get_field(a, "version")) == 12.0
@test isequal(jvector(get_field(a, "data")), [1,2,3])
delete(a)

type TestType
    name::String
    version::Float64
    data::Vector{Int}
end
t = TestType("MATLAB", 12.0, [1,2,3])
a = mxstruct(t)
@test is_struct(a)
@test nfields(a) == 3
@test jstring(get_field(a, "name")) == "MATLAB"
@test jscalar(get_field(a, "version")) == 12.0
@test isequal(jvector(get_field(a, "data")), [1,2,3])
delete(a)

a = mxstructarray([TestType("MATLAB", 12.0, [1,2,3]),
    TestType("Julia", 0.2, [4,5,6])])
@test is_struct(a)
@test nfields(a) == 3
@test jstring(get_field(a, 1, "name")) == "MATLAB"
@test jscalar(get_field(a, 1, "version")) == 12.0
@test isequal(jvector(get_field(a, 1, "data")), [1,2,3])
@test jstring(get_field(a, 2, "name")) == "Julia"
@test jscalar(get_field(a, 2, "version")) == 0.2
@test isequal(jvector(get_field(a, 2, "data")), [4,5,6])
delete(a)

# bi-directional conversions

x = mxarray(12.0)
y = jvariable(x)
delete(x)
@test isa(y, Float64)
@test y == 12.0

a = rand(5)
x = mxarray(a)
y = jvariable(x)
delete(x)
@test isa(y, Vector{Float64})
@test isequal(y, a)

a = rand(3, 4)
x = mxarray(a)
y = jvariable(x)
delete(x)
@test isa(y, Matrix{Float64})
@test isequal(y, a)

a = rand(3, 4, 5)
x = mxarray(a)
y = jvariable(x)
delete(x)
@test isa(y, Array{Float64, 3})
@test isequal(y, a)

a = "MATLAB"
x = mxarray(a)
y = jvariable(x)
delete(x)
@test isa(y, ASCIIString)
@test y == a

a = ["abc", 3, "efg"]
x = mxarray(a)
y = jvariable(x)
delete(x)
@test isa(y, Vector{Any})
@test length(y) == 3
@test y[1] == a[1]
@test y[2] == a[2]
@test y[3] == a[3]

a = @compat Dict{Any,Any}("abc"=>10.0, "efg"=>[1, 2, 3], "xyz"=>"MATLAB")
x = mxarray(a)
y = jvariable(x)
delete(x)
@test isa(y, Dict{String, Any})

@test y["abc"] == 10.0
@test isequal(y["efg"], [1, 2, 3])
@test y["xyz"] == "MATLAB"

gc()


text:
struct Markdown <: Format end

function parse(io::IO, ::Markdown; parser = default_md_parser())
    ast = parser(io)
    return convert(XTree, ast)
end

function default_md_parser()
    # adapted from https://github.com/MichaelHatherly/Publish.jl/blob/master/src/utilities.jl
    cm = CommonMark
    parser = cm.enable!(cm.Parser(), [
        ## CommonMark-provided.
        cm.AdmonitionRule(),
        cm.AttributeRule(),
        cm.AutoIdentifierRule(),
        cm.CitationRule(),
        cm.DollarMathRule(),
        cm.FootnoteRule(),
        cm.FrontMatterRule(toml=TOML.parse),
        cm.MathRule(),
        cm.RawContentRule(),
        cm.TableRule(),
        cm.TypographyRule(),
    ])
    return parser
end

extensionformat(::Val{:md}) = Markdown()
formatextension(::Markdown) = "md"

## Parsing helpers

function mdchildren(node::CommonMark.Node)
    if !isdefined(node.first_child, :t)
        return CommonMark.Node[]
    end

    child = node.first_child
    childs = [child]

    while child != node.last_child
        child = child.nxt
        push!(childs, child)
    end
    return childs
end

function mdchildrenattrs(node::CommonMark.Node)
    allcs = mdchildren(node)
    cs = CommonMark.Node[]
    attrs = Dict{Symbol, String}[]
    as = Dict{Symbol, String}()

    for (i, c) in enumerate(allcs)
        if c.t isa Attributes
            as = Dict{Symbol, Any}((Symbol(k), v) for (k, v) in c.t.dict)
            if haskey(as, :class)
                as[:class] = join(as[:class], ';')
            end
        else
            push!(cs, c)
            push!(attrs, as)
            as = Dict{Symbol, String}()
        end
    end
    return cs, attrs
end

function childrenxtrees(node::Node)
    cs, attrs = mdchildrenattrs(node)
    return XTree[convert(XTree, c, as) for (c, as) in zip(cs, attrs)]
end

Base.convert(::Type{XTree}, node::Node, attrs::Dict = Dict{Symbol, Any}()) = convert(XTree, node, node.t, attrs)

const BLOCK_TO_TAG = Dict(
    Document => :body,
    Item => :li,
    List => :ul,
    Paragraph => :p,
    Text => :span,
    Emph => :em,
    SoftBreak => :br,
    ThematicBreak => :hr,
    BlockQuote => :blockquote,
    Admonition => :admonition,
    Citation => :citation,
    CommonMark.Strong => :strong,
    CommonMark.Table => :table,
    CommonMark.TableHeader => :span,
    CommonMark.TableRow => :tr,
    CommonMark.TableCell => :td,
    CommonMark.TableBody => :div,
    CommonMark.FrontMatter => :fm,
)

function Base.convert(::Type{XTree}, node::Node, c::AbstractContainer, attrs = Dict{Symbol, String}())
    tag = BLOCK_TO_TAG[typeof(c)]
    # TODO: respect `Attributes`
    return XNode(tag, attrs, childrenxtrees(node))
end

# For some `AbstractContainer`s, the behavior is customized.

Base.convert(::Type{XTree}, node::Node, ::Text, attrs) = XLeaf(node.literal)
Base.convert(::Type{XTree}, node::Node, ::Code, attrs) = XNode(:code, [XLeaf(node.literal)])

function Base.convert(::Type{XTree}, node::Node, i::Image, attrs) XNode(:code, [XLeaf(node.literal)])
    return XNode(
        :img,
        Dict(:src => i.destination, :alt => i.title),
        childrenxtrees(node)
    )
end

function Base.convert(::Type{XTree}, node::Node, c::CodeBlock, attrs)
    return XNode(
        :pre,
        merge(attrs, Dict(:lang => c.info)),
        [XNode(:code, [XLeaf(node.literal)])],)
end

function Base.convert(::Type{XTree}, node::Node, l::Link, attrs)
    return XNode(
        :a,
        Dict(:href => l.destination, :title => l.title),
        childrenxtrees(node))
end

function Base.convert(::Type{XTree}, node::Node, c::Heading, attrs)
    tag = Symbol("h$(c.level)")
    return XNode(tag, attrs, childrenxtrees(node))
end

# TODO: add conversion for Table Nodes


text:
import StateSpaceReconstruction:
	AbstractEmbedding,
	LinearlyInvariantEmbedding

@recipe function f(E::AbstractEmbedding;
        vars = [1, 2, 3],
        mc = :black, ms = 1, lc = :black, lw = 1)
    if size(E.points, 1) < 3
        error("Cannot plot reconstruction with < 3 variables.")
    end
    legend --> false
    xlabel --> "x1"
    ylabel --> "x2"
    zlabel --> "x3"
    markercolor --> mc
    markersize --> ms
    ([E.points[i, :] for i in vars]...,)
end

@recipe function f(E::LinearlyInvariantEmbedding;
        vars = [1, 2, 3],
        mc = :black, ms = 1, lc = :black, lw = 1)
    if size(E.points, 1) < 3
        error("Cannot plot reconstruction with < 3 variables.")
    end
    legend --> false
    xlabel --> "x1"
    ylabel --> "x2"
    zlabel --> "x3"
    markercolor --> mc
    markersize --> ms
    ([E.points[i, :] for i in vars]...,)
end


text:
module TestCore

using Mutabilities
using Test

@testset "readonly" begin
    @testset for x in Any[[1, 2], [1 2; 3 4]]
        @test parent(readonly(x)) === x  # (no copy)
        @test parent(readonly(move!(x))) === x  # (no copy)
        @test readonly(x) == x
        @test_throws Exception readonly(x)[1] = 1
        @test_throws Exception push!(readonly(x), 1)

        z = readonly(x)
        x[1] = 111
        @test z[1] == 111

        if ndims(x) == 1
            push!(x, 333)
            @test z[end] == 333
        end
    end
end

@testset "freeze" begin
    @testset for x in Any[[1, 2], [1 2; 3 4]]
        @test parent(freeze(x)) !== x  # (copy)
        @test parent(freeze(move!(x))) === x  # (no copy)
        @test freeze(x) == x
        @test_throws Exception freeze(x)[1] = 1
        @test_throws Exception push!(freeze(x), 1)

        z = freeze(x)
        x[1] = 111
        @test z[1] != 111

        if ndims(x) == 1
            push!(x, 333)
            @test z[end] != 333
        end
    end
end

@testset "freezevalue" begin
    @testset for x in Any[[1, 2], [1 2; 3 4]]
        @test parent(freezevalue(x)) !== x  # (copy)
        @test parent(freezevalue(move!(x))) === x  # (no copy)
        @test freezevalue(x) == x
        @test_throws Exception freezevalue(x)[1] = 1
        if ndims(x) == 1
            @test push!(freezevalue(x), 333) == vcat(x, [333])
        end

        z = freezevalue(x)
        x[1] = 111
        @test z[1] != 111

        if ndims(x) == 1
            push!(x, 333)
            @test z[end] != 333
        end
    end
end

@testset "freezeindex" begin
    @testset for x in Any[[1, 2], [1 2; 3 4]]
        if x isa Matrix
            @test freezeindex(x) === x
        else
            @test parent(freezeindex(x)) !== x  # (copy)
        end
        @test parent(freezeindex(move!(x))) === x  # (no copy)
        @test freezeindex(x) == x
        @test_throws Exception push!(freezeindex(x), 333)

        x0 = copy(x)
        z = freezeindex(x)
        z[1] = 111
        @test z[1] == 111
        @test z != x0
    end
end

@testset "melt" begin
    @testset for freezer in [readonly, freeze, freezeindex, freezevalue]
        @testset for x in Any[[1, 2], [1 2; 3 4]]
            z = freezer(x)

            @test melt(z) !== x
            @test melt(move!(z)) === parent(z)
            if freezer === readonly || (freezer == freezeindex && x isa Matrix)
                @test melt(move!(z)) === x
            else
                @test melt(move!(z)) !== x
            end

            a = melt(z)
            a[1] = 123
            @test z[1] != 123

            b = melt(move!(z))
            b[1] = 123
            @test z[1] == 123
        end
    end
end

@testset "immutable input" begin
    @testset for x in Any[1:2, (1, 2), (a = (1, 2), c = 3)]
        @test readonly(x) === x
        @test freeze(x) === x
        @test freezevalue(x) === x
        @test freezeindex(x) === x
        @test readonly(move!(x)) === x
        @test freeze(move!(x)) === x
        @test freezevalue(move!(x)) === x
        @test freezeindex(move!(x)) === x
    end
end

@testset "showarg" begin
    @testset for freezer in [readonly, freeze, freezeindex, freezevalue]
        z = freezer(Symbol[])
        msg0 = sprint(Base.showarg, z, false)
        msg1 = sprint(Base.showarg, z, true)
        @test startswith(msg0, "$freezer(")
        @test startswith(msg1, "$freezer(")
        @test startswith(msg1, msg0)
        @test endswith(msg1, " with eltype Symbol")
    end
end

end  # module


text:
@testset "Lines2D" begin

  expected = Gray{N0f8}.([0 0 1 0 0
                          0 1 1 1 0
                          0 0 1 0 0
                          0 0 0 0 0
                          0 0 0 0 0])
  img = fill(Gray{N0f8}(0), 5, 5)
  draw!(img, Cross(Point(3,2), 1), Gray{N0f8}(1))
  @test img == expected

end


text:
@testset "indent" begin
    @test indent("", 4) == ""
    @test indent("\n", 4) == "\n"
    @test indent("  \n  ", 2) == "  \n  "

    @test indent("a", 4) == "    a"
    @test indent("a\n\nb", 4) == "    a\n\n    b"
end


text:
function get_inits(::MIMIXGaussian, params::Dict{Symbol, Any}, data::Dict{Symbol, Any})
    θ_init = iclr(proportionalize(data[:Y]))
    ω_init = params[:ω] == "ones" ? ones(data[:L], data[:p]) : zeros(data[:L], data[:p])
    π_init = params[:π] == "even" ? [0.5 for _ in 1:data[:p]] : params[:π]
    b_init = params[:b] == "zeros" ? zeros(data[:L], data[:p]) : params[:b]
    g_init = params[:g] == "zeros" ? zeros(data[:L], data[:q]) : params[:g]
    b_var_init = params[:b_var] == "ones" ? ones(data[:p]) : params[:b_var]
    g_var_init = params[:g_var] == "ones" ? ones(data[:num_blocking_factors]) : params[:g_var]

    inits = Dict{Symbol, Any}(
        :Y => data[:Y],
        :θ => θ_init,
        :θ_var => vec(var(θ_init, dims=1)),
        :μ => vec(mean(θ_init, dims=1)),
        :μ_var => params[:μ_var],
        :Λ => Matrix(1.0I, data[:L], data[:K]),
        :λ_var => ones(data[:L]),
        :F => zeros(data[:N], data[:L]),
        :b_full => b_init,
        :g => g_init,
        :ω => ω_init,
        :π => π_init,
        :b_var => b_var_init,
        :g_var => g_var_init
    )

    return inits
end

function get_model(::MIMIXGaussian, monitor::Dict{Symbol, Any}, hyper::Dict{Symbol, Any})

    model = Model(
        # High-dimensional counts
        Y = Stochastic(2,
            (m, ϕ, N) -> MultivariateDistribution[
                Multinomial(m[i], ϕ[i, :]) for i in 1:N
            ],
            monitor[:Y]
        ),

        ϕ = Logical(2,
            (θ) -> clr(θ),
            monitor[:ϕ]
        ),

        # Sample-specific random effects

        θ = Stochastic(2,
            (θ_mean, θ_var, N) -> MultivariateDistribution[
                MvNormal(θ_mean[i, :], sqrt.(θ_var)) for i in 1:N
            ],
            monitor[:θ]
        ),

        θ_mean = Logical(2,
            (μ, ΛF) -> ΛF .+ transpose(μ),
            monitor[:θ_mean]
        ),

        θ_var = Stochastic(1,
            () -> InverseGamma(1, 1),
            monitor[:θ_var]
        ),

        # Population mean

        μ = Stochastic(1,
            (K, μ_var) -> MvNormal(K, sqrt(μ_var)),
            monitor[:μ]
        ),

        μ_var = Stochastic(
            () -> InverseGamma(1, 1),
            monitor[:μ_var]
        ),

        # Full-dimensional fixed effect estimates (monitor these for OTU-level inference!)

        β = Logical(2,
            (Λ, b) -> transpose(Λ) * b,
            monitor[:β]
        ),

        # Factor analysis w/ loadings matrix Λ (L by K) and scores F (N by L)

        ΛF = Logical(2,
            (Λ, F) -> F * Λ,
            false
        ),

        Λ = Stochastic(2,
            (λ_var, L, K) -> MultivariateDistribution[
                MvNormal(zeros(K), sqrt(λ_var[l])) for l in 1:L
            ],
            monitor[:Λ]
        ),

        λ_var = Stochastic(1,
            () -> InverseGamma(1, 1),
            monitor[:λ_var]
        ),

        F = Stochastic(2,
            (F_mean, N) -> MultivariateDistribution[
                MvNormal(F_mean[i, :], 1.0) for i in 1:N
            ],
            false
        ),

        F_mean = Logical(2,
            (Xb, sumZg) -> Xb + sumZg,
            false
        ),

        # Blocking factor random effects

        sumZg = Logical(2,
            (Zg) -> transpose(squeezesum(Zg, 3)),
            false
        ),

        Zg = Logical(3,
            (Z, g) -> g[:, Z],
            false
        ),

        g = Stochastic(2,
            (q, L, g_var, blocking_factor) -> UnivariateDistribution[
                Normal(0.0, sqrt(g_var[blocking_factor[r]])) for l in 1:L, r in 1:q
            ],
            monitor[:g]
        ),

        g_var = Stochastic(1,
            () -> InverseGamma(1, 1),
            monitor[:g_var]
        ),

        # Fixed treatment effects w/ spike-and-slab variable selection

        Xb = Logical(2,
            (X, b) -> X * transpose(b),
            false
        ),

        b = Logical(2,
            (b_full, ω) -> ω .* b_full,
            monitor[:b]
        ),

        b_full = Stochastic(2,
            (p, L, b_var) -> MultivariateDistribution[
                MvNormal(zeros(p), sqrt.(b_var)) for l in 1:L
            ],
            monitor[:b_full]
        ),

        b_var = Stochastic(1,
            () -> InverseGamma(1, 1),
            monitor[:b_var]
        ),

        ω = Stochastic(2,
            (p, L, π) -> UnivariateDistribution[
                Bernoulli(π[j]) for l in 1:L, j in 1:p
            ],
            monitor[:ω]
        ),

        π = Stochastic(1,
            (L) -> Beta(1, L),
            monitor[:π]
        )
    )

    samplers = [
        Sampler(:θ, (Y, m, θ, θ_mean, θ_var, N, K) ->
            begin
                acc = 0
                for i in 1:N
                    eps = rand(Exponential(hyper[:hmc_epsilon]))
                    yi = Y[i, :]
                    mi = m[i]
                    u1 = θ[i, :]
                    θi_mean = θ_mean[i, :]
                    logf0, grad0 = logf1, grad1 = logfgrad(u1, θi_mean, θ_var, yi, mi)
                    v0 = v1 = randn(length(u1))
                    v1 += 0.5eps * grad0
                    for step in 1:hyper[:hmc_steps]
                        u1 += eps * v1
                        logf1, grad1 = logfgrad(u1, θi_mean, θ_var, yi, mi)
                        v1 += eps * grad1
                    end
                    v1 -= 0.5eps * grad1
                    v1 *= -1.0
                    Kv0 = 0.5sum(abs2, v0)
                    Kv1 = 0.5sum(abs2, v1)
                    if rand() < exp((logf1 - Kv1) - (logf0 - Kv0))
                        acc += 1
                        θ[i, :] = u1
                    end
                end
                if hyper[:hmc_verbose]
                    println("HMC accept rate: $(round(100acc/N, 1))%, ")
                end
                θ
            end
        ),

        Sampler(:μ, (θ, θ_var, μ, μ_var, ΛF, N, K) ->
            begin
                Sigma = 1 ./ (N ./ θ_var .+ 1 / μ_var)
                mu = Sigma .* vec(sum(θ - ΛF, dims=1)) ./ θ_var
                for k in 1:K
                    μ[k] = rand(Normal(mu[k], sqrt(Sigma[k])))
                end
                μ
            end
        ),

        Sampler(:Λ, (θ, θ_var, μ, Λ, λ_var, F, K, L) ->
            begin
                FtF = transpose(F) * F
                B = (transpose(F) * (θ .- transpose(μ))) ./ transpose(θ_var)
                Sigma = zeros(L, L)
                for k in 1:K
                    Sigma[:, :] = inv(cholesky(Hermitian(Diagonal(1 ./ λ_var) + FtF ./ θ_var[k])))
                    Λ[:, k] = rand(MvNormal(Sigma * B[:, k], Sigma))
                end
                Λ
            end
        ),

        Sampler(:F, (θ, θ_var, μ, Λ, F, F_mean, b_var, L, N) ->
            begin
                A = transpose(Λ) ./ θ_var
                Sigma = inv(cholesky(Hermitian(Λ * A + Matrix(1.0I, L, L))))
                mu = Sigma * transpose((θ .- transpose(μ)) * A + F_mean)
                for i in 1:N
                    F[i, :] = rand(MvNormal(mu[:, i], Sigma))
                end
                F
            end
        ),

        Sampler(:g, (F, Xb, g, g_var, Z, sumZg, q, L) ->
            begin
                A = F - Xb - sumZg
                for j in 1:size(Z, 2)
                    z = Z[:, j]
                    A += transpose(g[:, z])
                    for r in unique(z)
                        in_block = z .== r
                        B = vec(sum(A[in_block, :], dims=1))
                        Sigma = 1 / (sum(in_block) + 1 / g_var[j])
                        mu = Sigma .* B
                        for l in 1:L
                            g[l, r] = rand(Normal(mu[l], sqrt(Sigma)))
                        end
                    end
                    A -= transpose(g[:, z])
                end
                g
            end
        ),

        Sampler(:b_full, (F, sumZg, b_full, b_var, X, ω, L, p) ->
            begin
                A = F - sumZg
                for l in 1:L
                    Xtilde = X * Diagonal(ω[l, :])
                    Sigma = inv(cholesky(Hermitian(transpose(Xtilde) * Xtilde + Diagonal(1 ./ b_var))))
                    mu = Sigma * (transpose(Xtilde) * A[:, l])
                    b_full[l, :] = rand(MvNormal(mu, Sigma))
                end
                b_full
            end
        ),

        Sampler(:ω, (ω, F, sumZg, X, Xb, b_full, b, π, p, L) ->
            begin
                D = F - sumZg - Xb
                for l in 1:L
                    d = D[:, l]
                    for j in 1:p
                        x = X[:, j]
                        d += x .* b[l, j]
                        lprob_in = log(π[j]) - 0.5sum(abs2, d - x .* b_full[l, j])
                        lprob_out = log(1 - π[j]) - 0.5sum(abs2, d)
                        diff = min(lprob_in - lprob_out, 10)
                        prob_in = exp(diff) / (1 + exp(diff))
                        ω[l, j] = rand(Bernoulli(prob_in))
                        d -= x .* (ω[l, j] * b_full[l, j])
                    end
                end
                ω
            end
        ),

        Sampler(:π, (π, ω, L, p) ->
            begin
                S = vec(sum(ω, dims=1))
                c0, d0 = params(π.distr)
                c = c0 .+ S
                d = d0 + L .- S
                for j in 1:p
                    π[j] = rand(Beta(c[j], d[j]))
                end
                π
            end
        ),

        Sampler(:μ_var, (μ, μ_var, K) ->
            begin
                u = 0.5K + shape(μ_var.distr)
                v = 0.5sum(abs2, μ) + scale(μ_var.distr)
                rand(InverseGamma(u, v))
            end
        ),

        Sampler(:g_var, (g, g_var, L, q, blocking_factor) ->
            begin
                num_blocking_factor = length(g_var)
                u = fill(shape(g_var.distr), num_blocking_factor)
                v = fill(scale(g_var.distr), num_blocking_factor)
                G = sum(abs2, g, dims=1)
                for r in 1:q
                    u[blocking_factor[r]] += 0.5L
                    v[blocking_factor[r]] += 0.5G[r]
                end
                for i in 1:num_blocking_factor
                    g_var[i] = rand(InverseGamma(u[i], v[i]))
                end
                g_var
            end
        ),

        Sampler(:b_var, (b, b_var, ω, p) ->
            begin
                u = 0.5vec(sum(ω, dims=1)) .+ shape(b_var.distr)
                v = 0.5vec(sum(abs2, b, dims=1)) .+ scale(b_var.distr)
                for j in 1:p
                    b_var[j] = rand(InverseGamma(u[j], v[j]))
                end
                b_var
            end
        ),

        Sampler(:θ_var, (θ, θ_mean, θ_var, N, K) ->
            begin
                u = 0.5N + shape(θ_var.distr)
                v = 0.5vec(sum(abs2, θ - θ_mean, dims=1)) .+ scale(θ_var.distr)
                for k in 1:K
                    θ_var[k] = rand(InverseGamma(u, v[k]))
                end
                θ_var
            end
        ),

        Sampler(:λ_var, (λ_var, Λ, K, L) ->
        begin
            u = 0.5K + shape(λ_var.distr)
            v = 0.5vec(sum(abs2, Λ, dims=2)) .+ scale(λ_var.distr)
            [rand(InverseGamma(u, v[l])) for l in 1:L]
        end
    )
    ]

    setsamplers!(model, samplers)
    return model
end

max_stars_repo_name: reich-group/MIMIX

max_issues_repo_name: reich-group/MIMIX

max_forks_repo_name: reich-group/MIMIX


text:
abstract Simplexer

immutable AffineSimplexer <: Simplexer
    a::Float64
    b::Float64
end
AffineSimplexer(;a = 0.025, b = 0.5) = AffineSimplexer(a, b)

function simplexer{T, N}(S::AffineSimplexer, initial_x::Array{T, N})
    n = length(initial_x)
    initial_simplex = Array{T, N}[copy(initial_x) for i = 1:n+1]
    for j = 1:n
        initial_simplex[j+1][j] = (1+S.b) * initial_simplex[j+1][j] + S.a
    end
    initial_simplex
end

abstract NMParameters

immutable AdaptiveParameters <: NMParameters
    α::Float64
    β::Float64
    γ::Float64
    δ::Float64
end

AdaptiveParameters(;  α = 1.0, β = 1.0, γ = 0.75 , δ = 1.0) = AdaptiveParameters(α, β, γ, δ)
parameters(P::AdaptiveParameters, n::Integer) = (P.α, P.β + 2/n, P.γ - 1/2n, P.δ - 1/n)

immutable FixedParameters <: NMParameters
    α::Float64
    β::Float64
    γ::Float64
    δ::Float64
end

FixedParameters(; α = 1.0, β = 2.0, γ = 0.5, δ = 0.5) = FixedParameters(α, β, γ, δ)
parameters(P::FixedParameters, n::Integer) = (P.α, P.β, P.γ, P.δ)

immutable NelderMead{Ts <: Simplexer, Tp <: NMParameters} <: Optimizer
    initial_simplex::Ts
    parameters::Tp
end

function NelderMead(; kwargs...)
    KW = Dict(kwargs)
    if haskey(KW, :a) || haskey(KW, :g) || haskey(KW, :b)
        a, g, b = 1.0, 2.0, 0.5
        haskey(KW, :a) && (a = KW[:a])
        haskey(KW, :g) && (g = KW[:g])
        haskey(KW, :b) && (b = KW[:b])
        return NelderMead(a, g, b)
    elseif haskey(KW, :initial_simplex) || haskey(KW, :parameters)
        initial_simplex, parameters = AffineSimplexer(), AdaptiveParameters()
        haskey(KW, :initial_simplex) && (initial_simplex = KW[:initial_simplex])
        haskey(KW, :parameters) && (parameters = KW[:parameters])
        return NelderMead(initial_simplex, parameters)
    else
        return NelderMead(AffineSimplexer(), AdaptiveParameters())
    end
end

# centroid except h-th vertex
function centroid!{T}(c::Array{T}, simplex, h=0)
    n = length(c)
    fill!(c, zero(T))
    @inbounds for i in 1:n+1
        if i != h
            xi = simplex[i]
            @simd for j in 1:n
                c[j] += xi[j]
            end
        end
    end
    scale!(c, 1/n)
end

centroid(simplex, h) = centroid!(similar(simplex[1]), simplex, h)

nmobjective(y::Vector, m::Integer, n::Integer) = sqrt(var(y) * (m / n))

macro nmtrace()
    quote
        if tracing
            dt = Dict()
            if o.extended_trace
                dt["centroid"] = x_centroid
                dt["step_type"] = step_type
            end
            update!(tr,
                    iteration,
                    f_lowest,
                    f_x,
                    dt,
                    o.store_trace,
                    o.show_trace,
                    o.show_every,
                    o.callback)
        end
    end
end

function print_header(mo::NelderMead, options::OptimizationOptions)
    if options.show_trace
        @printf "Iter     Function value    √(Σ(yᵢ-ȳ)²)/n \n"
        @printf "------   --------------    --------------\n"
    end
end

function Base.show(io::IO, t::OptimizationTrace{NelderMead})
    @printf io "Iter     Function value    √(Σ(yᵢ-ȳ)²)/n \n"
    @printf io "------   --------------    --------------\n"
    for state in t.states
        show(io, state)
    end
    return
end

function Base.show(io::IO, t::OptimizationState{NelderMead})
    @printf io "%6d   %14e    %14e\n" t.iteration t.value t.g_norm
    if !isempty(t.metadata)
        for (key, value) in t.metadata
            @printf io " * %s: %s\n" key value
        end
    end
    return
end

function optimize{T}(f::Function,
                     initial_x::Vector{T},
                     mo::NelderMead,
                     o::OptimizationOptions)

    # Print header if show_trace is set
    print_header(mo, o)

    # Set up a simplex of points
    m = length(initial_x)
    if m == 1
        error("Use optimize(f, scalar, scalar) for 1D problems")
    end
    n = m + 1
    simplex = simplexer(mo.initial_simplex, initial_x)
    f_simplex = zeros(T, n)
    @inbounds for i in 1:length(simplex)
        f_simplex[i] = f(simplex[i])
    end

    # Get the indeces that correspond to the ordering of the f values
    # at the vertices. i_order[1] is the index in the simplex of the vertex
    # with the lowest function value, and i_order[end] is the index in the
    # simplex of the vertex with the highest function value
    i_order = sortperm(f_simplex)
    x_centroid = centroid(simplex,  i_order[n])

    # Count function calls
    f_calls = n

    # Setup parameters
    α, β, γ, δ = parameters(mo.parameters, m)
    # Count iterations
    iteration = 0

    step_type = "initial"

    # Maintain a trace
    f_x_previous, f_x = NaN, nmobjective(f_simplex, m, n)
    f_lowest = f_simplex[i_order[1]]
    tr = OptimizationTrace{typeof(mo)}()
    tracing = o.show_trace || o.store_trace || o.extended_trace || o.callback != nothing
    @nmtrace

    # Cache x_centroid, y_bar, x_reflect, x_expand, x_lowest, x_highest
    x_reflect = Array(T, m)
    x_expand = Array(T, m)
    x_cache = Array(T, m)

    x_lowest = Array(T, m)
    x_second_highest = Array(T, m)
    x_highest = Array(T, m)

    # Iterate until convergence or exhaustion
    x_converged = false
    f_converged = false
    g_converged = false

    while !g_converged && !f_converged && iteration < o.iterations
        # Augment the iteration counter
        shrink = false
        iteration += 1

        centroid!(x_centroid, simplex,  i_order[n])
        copy!(x_lowest, simplex[i_order[1]])
        copy!(x_second_highest, simplex[i_order[m]])
        copy!(x_highest, simplex[ i_order[n]])

        f_lowest = f_simplex[i_order[1]]
        f_second_highest = f_simplex[i_order[m]]
        f_highest = f_simplex[ i_order[n]]
        # Compute a reflection
        @inbounds for j in 1:m
            x_reflect[j] = x_centroid[j] + α * (x_centroid[j]-x_highest[j])
        end

        f_reflect = f(x_reflect)
        f_calls += 1
        if f_reflect < f_lowest
            # Compute an expansion
            @inbounds for j in 1:m
                x_cache[j] = x_centroid[j] + β *(x_reflect[j] - x_centroid[j])
            end
            f_expand = f(x_cache)
            f_calls += 1

            if f_expand < f_reflect
                copy!(simplex[ i_order[n]], x_cache)
                @inbounds f_simplex[ i_order[n]] = f_expand
                step_type = "expansion"
            else
                copy!(simplex[ i_order[n]], x_reflect)
                @inbounds f_simplex[ i_order[n]] = f_reflect
                step_type = "reflection"
            end
            # shift all order indeces, and wrap the last one around to the first
            i_highest = i_order[n]
            @inbounds for i = n:-1:2
                i_order[i] = i_order[i-1]
            end
            i_order[1] = i_highest
        elseif f_reflect < f_second_highest
            copy!(simplex[ i_order[n]], x_reflect)
            @inbounds f_simplex[ i_order[n]] = f_reflect
            step_type = "reflection"
            sortperm!(i_order, f_simplex)
        else
            if f_reflect < f_highest
                # Outside contraction
                @simd for j in 1:m
                    @inbounds x_cache[j] = x_centroid[j] + γ * (x_reflect[j]-x_centroid[j])
                end
                f_outside_contraction = f(x_cache)
                if f_outside_contraction < f_reflect
                    copy!(simplex[ i_order[n]], x_cache)
                    @inbounds f_simplex[ i_order[n]] = f_outside_contraction
                    step_type = "outside contraction"
                    sortperm!(i_order, f_simplex)

                else
                    shrink = true
                end
            else # f_reflect > f_highest
                # Inside constraction
                @simd for j in 1:m
                    @inbounds x_cache[j] = x_centroid[j] - γ *(x_reflect[j] - x_centroid[j])
                end
                f_inside_contraction = f(x_cache)
                if f_inside_contraction < f_highest
                    copy!(simplex[ i_order[n]], x_cache)
                    @inbounds f_simplex[ i_order[n]] = f_inside_contraction
                    step_type = "inside contraction"
                    sortperm!(i_order, f_simplex)
                else
                    shrink = true
                end
            end
        end

        if shrink
            for i = 2:n
                ord = i_order[i]
                copy!(simplex[ord], x_lowest + δ*(simplex[ord]-x_lowest))
                f_simplex[ord] = f(simplex[ord])
            end
            step_type = "shrink"
            sortperm!(i_order, f_simplex)
        end

        f_x_previous, f_x = f_x, nmobjective(f_simplex, m, n)

        @nmtrace
        if f_x <= o.g_tol
            g_converged = true
        end
    end

    sortperm!(i_order, f_simplex)
    x_centroid_min = centroid(simplex,  i_order[n])
    f_centroid_min = f(x_centroid)
    f_calls += 1
    f_min, i_f_min = findmin(f_simplex)
    x_min = simplex[i_f_min]
    if f_centroid_min < f_min
        x_min = x_centroid_min
        f_min = f_centroid_min
    end

    return MultivariateOptimizationResults("Nelder-Mead",
                                           initial_x,
                                           x_min,
                                           Float64(f_min),
                                           iteration,
                                           iteration == o.iterations,
                                           x_converged,
                                           NaN,
                                           f_converged,
                                           NaN,
                                           g_converged,
                                           o.g_tol,
                                           tr,
                                           f_calls,
                                           0)
end


text:
#
# unidimensional.jl --
#
# Unidimensional interpolation (the result may however be multi-dimensional).
#
#------------------------------------------------------------------------------
#
# This file is part of the LinearInterpolators package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2016-2021, Éric Thiébaut.
#

# All code is in a module to "hide" private methods.
module UnidimensionalInterpolators

using InterpolationKernels

using LazyAlgebra
using LazyAlgebra.Foundations
import LazyAlgebra: apply, apply!, vcreate

using ..LinearInterpolators
using ..LinearInterpolators: limits, getcoefs
import ..LinearInterpolators.Meta

#------------------------------------------------------------------------------
# Out-of place versions (coordinates cannot be a function).

# FIXME: This is type-piracy because the only non-standard type, `Kernel`, comes
#        from another package.

function apply(ker::Kernel,
               x::AbstractArray,
               src::AbstractVector)
    return apply(Direct, ker, x, src)
end

function apply(::Type{Direct}, ker::Kernel,
               x::AbstractArray,
               src::AbstractVector)
    # The element type of the result is the type T of the product of the
    # interpolation weights by the elements of the source array.
    T = promote_type(eltype(ker), eltype(src))
    return apply!(Array{T}(undef, size(x)), Direct, ker, x, src)
end

function apply(::Type{Adjoint}, ker::Kernel,
               x::AbstractArray,
               src::AbstractArray, len::Integer)
    # The element type of the result is the type T of the product of the
    # interpolation weights by the elements of the source array.
    T = promote_type(eltype(ker), eltype(src))
    return apply!(Array{T}(undef, len), Adjoint, ker, x, src)
end

#------------------------------------------------------------------------------
# In-place wrappers.

"""
    apply!(dst, ker, x, src)
Interpolate `src` using kernel `ker` and samples `x`,
writing result to `dst`.
"""
function apply!(dst::AbstractArray{<:Any,N},
                ker::Kernel,
                x::Union{Function,AbstractArray{<:Any,N}},
                src::AbstractVector) where {N}
    apply!(1, Direct, ker, x, src, 0, dst)
end

function apply!(dst::AbstractArray{<:Any,N},
                ::Type{Direct},
                ker::Kernel,
                x::Union{Function,AbstractArray{<:Any,N}},
                src::AbstractVector) where {N}
    apply!(1, Direct, ker, x, src, 0, dst)
end

"""
    apply!(dst, Adjoint, ker, x, src)
Adjoint of interpolation using kernel `ker` and samples `x` from `src`,
writing result to `dst`.
"""
function apply!(dst::AbstractVector,
                ::Type{Adjoint},
                ker::Kernel,
                x::Union{Function,AbstractArray{<:Any,N}},
                src::AbstractArray{<:Any,N}) where {N}
    apply!(1, Adjoint, ker, x, src, 0, dst)
end

#------------------------------------------------------------------------------
# In-place direct operation.

function __generate_interp(n::Integer, ker::Symbol, lim::Meta.Arg,
                           pos::Meta.Arg, arr::Symbol)
    J = [Symbol(:j_,s) for s in 1:n]
    W = [Symbol(:w_,s) for s in 1:n]
    code = Meta.generate_getcoefs(J, W, ker, lim, pos)
    expr = Meta.generate_interp_expr(arr, J, W)
    return (code, expr)
end

@generated function apply!(α::Number,
                           ::Type{Direct},
                           ker::Kernel{<:Any,S},
                           x::AbstractArray{<:Any,N},
                           src::AbstractVector,
                           β::Number,
                           dst::AbstractArray{<:Any,N}) where {S,N}
    code, expr = __generate_interp(S, :ker, :lim, :pos, :src)
    quote
        @assert size(dst) == size(x)
        lim = limits(ker, length(src))
        if α == 0
            vscale!(dst, β)
        elseif β == 0
            if α == 1
                @inbounds for i in eachindex(dst, x)
                    pos = x[i]
                    $code
                    dst[i] = $expr
                end
            else
                alpha = promote_multiplier(α, eltype(ker), eltype(src))
                @inbounds for i in eachindex(dst, x)
                    pos = x[i]
                    $code
                    dst[i] = $expr*alpha
                end

            end
        else
            alpha = promote_multiplier(α, eltype(ker), eltype(src))
            beta  = promote_multiplier(β, eltype(dst))
            @inbounds for i in eachindex(dst, x)
                pos = x[i]
                $code
                dst[i] = $expr*alpha + beta*dst[i]
            end
        end
        return dst
    end
end

@generated function apply!(α::Number,
                           ::Type{Direct},
                           ker::Kernel{<:Any,S},
                           f::Function,
                           src::AbstractVector,
                           β::Number,
                           dst::AbstractArray{<:Any,N}) where {S,N}
    code, expr = __generate_interp(S, :ker, :lim, :pos, :src)
    quote
        lim = limits(ker, length(src))
        if α == 0
            vscale!(dst, β)
        elseif β == 0
            if α == 1
                @inbounds for i in eachindex(dst)
                    pos = f(i)
                    $code
                    dst[i] = $expr
                end
            else
                alpha = promote_multiplier(α, eltype(ker), eltype(src))
                @inbounds for i in eachindex(dst)
                    pos = f(i)
                    $code
                    dst[i] = $expr*alpha
                end
            end
        else
            alpha = promote_multiplier(α, eltype(ker), eltype(src))
            beta  = promote_multiplier(β, eltype(dst))
            @inbounds for i in eachindex(dst)
                pos = f(i)
                $code
                dst[i] = $expr*alpha + beta*dst[i]
            end
        end
        return dst
    end
end

#------------------------------------------------------------------------------
# In-place adjoint operation.

function __generate_interp_adj(n::Integer, ker::Symbol, lim::Meta.Arg,
                               pos::Meta.Arg, dst::Symbol, val::Symbol)
    J = [Symbol(:j_,s) for s in 1:n]
    W = [Symbol(:w_,s) for s in 1:n]
    return (Meta.generate_getcoefs(J, W, ker, lim, pos),
            [:($dst[$(J[i])] += $(W[i])*$val) for i in 1:n]...)
end

@generated function apply!(α::Number,
                           ::Type{Adjoint},
                           ker::Kernel{<:Any,S},
                           x::AbstractArray{<:Any,N},
                           src::AbstractArray{<:Any,N},
                           β::Number,
                           dst::AbstractVector) where {S,N}
    code = __generate_interp_adj(S, :ker, :lim, :pos, :dst, :val)
    quote
        @assert size(src) == size(x)
        if β != 0
            vscale!(dst, β)
        else
            vzero!(dst)
        end
        lim = limits(ker, length(dst))
        if α == 1
            @inbounds for i in eachindex(src, x)
                pos = x[i]
                val = src[i]
                $(code...)
            end
        elseif α != 0
            alpha = promote_multiplier(α, eltype(ker), eltype(src))
            @inbounds for i in eachindex(src, x)
                pos = x[i]
                val = alpha*src[i]
                $(code...)
            end
        end
        return dst
    end
end

@generated function apply!(α::Number,
                           ::Type{Adjoint},
                           ker::Kernel{<:Any,S},
                           f::Function,
                           src::AbstractArray{<:Any,N},
                           β::Number,
                           dst::AbstractVector) where {S,N}
    code = __generate_interp_adj(S, :ker, :lim, :pos, :dst, :val)
    quote
        if β != 0
            vscale!(dst, β)
        else
            vzero!(dst)
        end
        lim = limits(ker, length(dst))
        if α == 1
            @inbounds for i in eachindex(src)
                pos = f(i)
                val = src[i]
                $(code...)
            end
        elseif α != 0
            alpha = promote_multiplier(α, eltype(ker), eltype(src))
            @inbounds for i in eachindex(src)
                pos = f(i)
                val = alpha*src[i]
                $(code...)
            end
        end
        return dst
    end
end

end # module


text:
module GibbsSeaWater

using Compat

deps = joinpath(dirname(dirname(pathof(GibbsSeaWater))), "deps", "deps.jl")
isfile(deps) ? include(deps) : error("GibbsSeaWater is not properly installed")

include("gen_gswteos_h.jl")
include("gen_gswteos10.jl")
include("wrappers.jl")

"""
wrap a gsw function with `nin` input arguments and `nout` output arguments defined
as references. All types are `Cdouble`.

Expanded functions look like these:
https://github.com/TEOS-10/GibbsSeaWater.jl/blob/9375647d97fe7cd54e11135a18e0c5447ddbed8c/src/GibbsSeaWater.jl#L17
"""
macro wr(name, nin, nout)
    return Expr(:function,
                Expr(:call,
                     Symbol(:gsw_, name),
                     (Expr(:(::), Symbol(:i, i), :Cdouble) for i = 1:nin)...),
                Expr(:block,
                     ( Expr(:(=), Symbol(:p, i), :(Ref{Cdouble}())) for i = 1:nout)...,
                     Expr(:call,
                          Symbol(:gsw_, name),
                          (Symbol(:i, i) for i = 1:nin)..., (Symbol(:p, i) for i = 1:nout)...),
                     Expr(:return,
                          Expr(:tuple,
                               (Expr(:ref, Symbol(:p, i)) for i = 1:nout)...
                               )
                          )
                     )
                )
end

#   function name                                nin nout
@wr specvol_alpha_beta                           3   3
@wr specvol_first_derivatives                    3   3
@wr specvol_second_derivatives                   3   5
@wr specvol_first_derivatives_wrt_enthalpy       3   2
@wr specvol_second_derivatives_wrt_enthalpy      3   3
@wr rho_alpha_beta                               3   3
@wr rho_first_derivatives                        3   3
@wr rho_second_derivatives                       3   5
@wr rho_first_derivatives_wrt_enthalpy           3   2
@wr rho_second_derivatives_wrt_enthalpy          3   3

export gsw_add_barrier
export gsw_add_mean
export gsw_adiabatic_lapse_rate_from_ct
export gsw_adiabatic_lapse_rate_ice
export gsw_alpha
export gsw_alpha_on_beta
export gsw_alpha_wrt_t_exact
export gsw_alpha_wrt_t_ice
export gsw_beta_const_t_exact
export gsw_beta
export gsw_cabbeling
export gsw_c_from_sp
export gsw_chem_potential_water_ice
export gsw_chem_potential_water_t_exact
export gsw_cp_ice
export gsw_cp_t_exact
export gsw_ct_first_derivatives
export gsw_ct_first_derivatives_wrt_t_exact
export gsw_ct_freezing
export gsw_ct_freezing_first_derivatives
export gsw_ct_freezing_first_derivatives_poly
export gsw_ct_freezing_poly
export gsw_ct_from_enthalpy
export gsw_ct_from_enthalpy_exact
export gsw_ct_from_entropy
export gsw_ct_from_pt
export gsw_ct_from_rho
export gsw_ct_from_t
export gsw_ct_maxdensity
export gsw_ct_second_derivatives
export gsw_deltasa_atlas
export gsw_deltasa_from_sp
export gsw_dilution_coefficient_t_exact
export gsw_dynamic_enthalpy
export gsw_enthalpy_ct_exact
export gsw_enthalpy_diff
export gsw_enthalpy
export gsw_enthalpy_first_derivatives_ct_exact
export gsw_enthalpy_first_derivatives
export gsw_enthalpy_ice
export gsw_enthalpy_second_derivatives_ct_exact
export gsw_enthalpy_second_derivatives
export gsw_enthalpy_sso_0
export gsw_enthalpy_t_exact
export gsw_entropy_first_derivatives
export gsw_entropy_from_ct
export gsw_entropy_from_pt
export gsw_entropy_from_t
export gsw_entropy_ice
export gsw_entropy_part
export gsw_entropy_part_zerop
export gsw_entropy_second_derivatives
export gsw_fdelta
export gsw_frazil_properties
export gsw_frazil_properties_potential
export gsw_frazil_properties_potential_poly
export gsw_frazil_ratios_adiabatic
export gsw_frazil_ratios_adiabatic_poly
export gsw_geo_strf_dyn_height
export gsw_geo_strf_dyn_height_1
export gsw_geo_strf_dyn_height_pc
export gsw_gibbs_ice
export gsw_gibbs_ice_part_t
export gsw_gibbs_ice_pt0
export gsw_gibbs_ice_pt0_pt0
export gsw_gibbs
export gsw_gibbs_pt0_pt0
export gsw_grav
export gsw_helmholtz_energy_ice
export gsw_hill_ratio_at_sp2
export gsw_ice_fraction_to_freeze_seawater
export gsw_internal_energy
export gsw_internal_energy_ice
export gsw_ipv_vs_fnsquared_ratio
export gsw_kappa_const_t_ice
export gsw_kappa
export gsw_kappa_ice
export gsw_kappa_t_exact
export gsw_latentheat_evap_ct
export gsw_latentheat_evap_t
export gsw_latentheat_melting
export gsw_linear_interp_sa_ct
export gsw_melting_ice_equilibrium_sa_ct_ratio
export gsw_melting_ice_equilibrium_sa_ct_ratio_poly
export gsw_melting_ice_into_seawater
export gsw_melting_ice_sa_ct_ratio
export gsw_melting_ice_sa_ct_ratio_poly
export gsw_melting_seaice_equilibrium_sa_ct_ratio
export gsw_melting_seaice_equilibrium_sa_ct_ratio_poly
export gsw_melting_seaice_into_seawater
export gsw_melting_seaice_sa_ct_ratio
export gsw_melting_seaice_sa_ct_ratio_poly
export gsw_nsquared
export gsw_o2sol
export gsw_o2sol_sp_pt
export gsw_pot_enthalpy_from_pt_ice
export gsw_pot_enthalpy_from_pt_ice_poly
export gsw_pot_enthalpy_ice_freezing
export gsw_pot_enthalpy_ice_freezing_first_derivatives
export gsw_pot_enthalpy_ice_freezing_first_derivatives_poly
export gsw_pot_enthalpy_ice_freezing_poly
export gsw_pot_rho_t_exact
export gsw_pressure_coefficient_ice
export gsw_pressure_freezing_ct
export gsw_pt0_cold_ice_poly
export gsw_pt0_from_t
export gsw_pt0_from_t_ice
export gsw_pt_first_derivatives
export gsw_pt_from_ct
export gsw_pt_from_entropy
export gsw_pt_from_pot_enthalpy_ice
export gsw_pt_from_pot_enthalpy_ice_poly_dh
export gsw_pt_from_pot_enthalpy_ice_poly
export gsw_pt_from_t
export gsw_pt_from_t_ice
export gsw_pt_second_derivatives
export gsw_rho_alpha_beta
export gsw_rho
export gsw_rho_first_derivatives
export gsw_rho_first_derivatives_wrt_enthalpy
export gsw_rho_ice
export gsw_rho_second_derivatives
export gsw_rho_second_derivatives_wrt_enthalpy
export gsw_rho_t_exact
export gsw_rr68_interp_sa_ct
export gsw_saar
export gsw_sa_freezing_estimate
export gsw_sa_freezing_from_ct
export gsw_sa_freezing_from_ct_poly
export gsw_sa_freezing_from_t
export gsw_sa_freezing_from_t_poly
export gsw_sa_from_rho
export gsw_sa_from_sp_baltic
export gsw_sa_from_sp
export gsw_sa_from_sstar
export gsw_sa_p_inrange
export gsw_seaice_fraction_to_freeze_seawater
export gsw_sigma0
export gsw_sigma1
export gsw_sigma2
export gsw_sigma3
export gsw_sigma4
export gsw_sound_speed
export gsw_sound_speed_ice
export gsw_sound_speed_t_exact
export gsw_specvol_alpha_beta
export gsw_specvol_anom_standard
export gsw_specvol
export gsw_specvol_first_derivatives
export gsw_specvol_first_derivatives_wrt_enthalpy
export gsw_specvol_ice
export gsw_specvol_second_derivatives
export gsw_specvol_second_derivatives_wrt_enthalpy
export gsw_specvol_sso_0
export gsw_specvol_t_exact
export gsw_sp_from_c
export gsw_sp_from_sa_baltic
export gsw_sp_from_sa
export gsw_sp_from_sk
export gsw_sp_from_sr
export gsw_sp_from_sstar
export gsw_spiciness0
export gsw_spiciness1
export gsw_spiciness2
export gsw_sr_from_sp
export gsw_sstar_from_sa
export gsw_sstar_from_sp
export gsw_t_deriv_chem_potential_water_t_exact
export gsw_t_freezing
export gsw_t_freezing_first_derivatives_poly
export gsw_t_freezing_first_derivatives
export gsw_t_freezing_poly
export gsw_t_from_ct
export gsw_t_from_pt0_ice
export gsw_thermobaric
export gsw_turner_rsubrho
export gsw_util_indx
export gsw_util_interp1q_int
export gsw_util_linear_interp
export gsw_util_sort_real
export gsw_util_xinterp1
export gsw_util_pchip_interp
export gsw_z_from_p
export gsw_p_from_z

end # module


text:
@testset "174.dungeon-game.jl" begin
    @test calculate_minimum_hp([[-2,-3,3],[-5,-10,1],[10,30,-5]]) == 7
    @test calculate_minimum_hp([[9, 0, 0, -6, 5], [6, 2, 6, 8, -10], [-1, -5, 0, 6, 1], [-1, 8, -6, -5, -10]]) == 1
end


text:
################################################################################
#
#          NfOrd/FracIdl.jl : Elements of orders of number fields
#
# This file is part of hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

export parent, order, basis_mat, basis_mat_inv, basis, norm,
       ring_of_multipliers, ==

################################################################################
#
#  Parent
#
################################################################################

parent(a::NfOrdFracIdl) = a.parent

################################################################################
#
#  Order
#
################################################################################

order(a::NfOrdFracIdlSet) = a.order

order(a::NfOrdFracIdl) = order(parent(a))

################################################################################
#
#  (Inverse) basis matrix
#
################################################################################

doc"""
***
    basis_mat(I::NfOrdFracIdl) -> FakeFmpqMat

> Returns the basis matrix of $I$ with respect to the basis of the order.
"""
function basis_mat(a::NfOrdFracIdl)
  return deepcopy(a.basis_mat)
end

doc"""
***
    basis_mat_inv(I::NfOrdFracIdl) -> FakeFmpqMat

> Returns the inverse of the basis matrix of $I$.
"""
function basis_mat_inv(a::NfOrdFracIdl)
  if isdefined(a, :basis_mat_inv)
    return deepcopy(a.basis_mat_inv)
  else
    a.basis_mat_inv = inv(basis_mat(a))
    return deepcopy(a.basis_mat_inv)
  end
end

################################################################################
#
#  Basis
#
################################################################################

doc"""
***
    basis(I::NfOrdFracIdl) -> Array{nf_elem, 1}

> Returns the $\mathbf Z$-basis of $I$.
"""
function basis(a::NfOrdFracIdl)
  B = basis_mat(a)
  d = degree(order(a))
  O = order(a)
  K = nf(O)
  Oba = O.basis_nf
  res = Array(nf_elem, d)

  for i in 1:d
    z = K()
    for j in 1:d
      z = z + B.num[i, j]*Oba[j]
    end
    z = divexact(z, B.den)
    res[i] = z
  end

  return res
end

################################################################################
#
#  Norm
#
################################################################################

doc"""
***
    norm(I::NfOrdFracIdl) -> fmpq

> Returns the norm of $I$
"""
function norm(a::NfOrdFracIdl)
  if isdefined(a, :norm)
    return deepcopy(a.norm)
  else
    a.norm = det(basis_mat(a))
    return deepcopy(a.norm)
  end
end

################################################################################
#
#  Equality
#
################################################################################

doc"""
***
    ==(x::NfOrdFracIdl, y::NfOrdFracIdl) -> Bool

> Returns whether $x$ and $y$ are equal.
"""
==(x::NfOrdFracIdl, y::NfOrdFracIdl) = basis_mat(x) == basis_mat(y)


text:
#Presumably I will add more ways to initialize cells...

####################################################
# Functions to set initial cell locations
####################################################

#Creates a cluster of cells in the middle of the grid
function GrowCells(gr::SimpleGraph, M::ModelParameters)

    #initialize matrix of cell indicies (σ)
    cellMembership = zeros(Int, M.graphDimension...)

    #Find the center node
    centerIdx = CartesianIndex(M.graphDimension.÷2)
    nodeIdx = LinearIndices(M.graphDimension)
    centerNode = nodeIdx[centerIdx]

    #Determine how far nodes are from the center
    nodeDis = gdistances(gr, centerNode)

    #How many nodes need to be initialized?
    totalNodes = M.cellCounts ⋅ M.cellVolumes #dot product to sum and multiply

    #Get a sorted permutation of the distance
    sortedDis = sortperm(nodeDis)

    #Assign 1:totalNodes to be filled with cells and the rest medium
    networkIdx = sortedDis[1:totalNodes]

    #Partition the identified nodes by the number of cells needed
    if sum(M.cellCounts) == 1 #There is only one cell (no need to partition)
        cellMembership[networkIdx] .= 1
    else
        cellMembership[networkIdx] = Metis.partition(gr[networkIdx], sum(M.cellCounts))
    end

    return cellMembership
end


text:
# StepSimulator
# maintained by @zsunberg

mutable struct StepSimulator <: Simulator
    rng::AbstractRNG
    max_steps::Union{Nothing,Any}
    spec
end
function StepSimulator(spec; rng=Random.GLOBAL_RNG, max_steps=nothing)
    return StepSimulator(rng, max_steps, spec)
end

function simulate(sim::StepSimulator, mdp::MDP{S}, policy::Policy, init_state::S=initialstate(mdp, sim.rng)) where {S}
    symtuple = convert_spec(sim.spec, MDP)
    max_steps = something(sim.max_steps, typemax(Int64))
    return MDPSimIterator(symtuple, mdp, policy, sim.rng, init_state, max_steps)
end

function simulate(sim::StepSimulator, pomdp::POMDP, policy::Policy, bu::Updater=updater(policy))
    dist = initialstate_distribution(pomdp)
    return simulate(sim, pomdp, policy, bu, dist)
end

function simulate(sim::StepSimulator, pomdp::POMDP, policy::Policy, bu::Updater, dist::Any, is=initialstate(pomdp, sim.rng))
    initial_belief = initialize_belief(bu, dist)
    symtuple = convert_spec(sim.spec, POMDP)
    max_steps = something(sim.max_steps, typemax(Int64))
    return POMDPSimIterator(symtuple, pomdp, policy, bu, sim.rng, initial_belief, is, max_steps)
end

struct MDPSimIterator{SPEC, M<:MDP, P<:Policy, RNG<:AbstractRNG, S}
    mdp::M
    policy::P
    rng::RNG
    init_state::S
    max_steps::Int
end

function MDPSimIterator(spec::Union{Tuple, Symbol}, mdp::MDP, policy::Policy, rng::AbstractRNG, init_state, max_steps::Int)
    return MDPSimIterator{spec, typeof(mdp), typeof(policy), typeof(rng), typeof(init_state)}(mdp, policy, rng, init_state, max_steps)
end

Base.IteratorSize(::Type{<:MDPSimIterator}) = Base.SizeUnknown()

function Base.iterate(it::MDPSimIterator, is::Tuple{Int, S}=(1, it.init_state)) where S
    if isterminal(it.mdp, is[2]) || is[1] > it.max_steps
        return nothing
    end
    t = is[1]
    s = is[2]
    a, ai = action_info(it.policy, s)
    on = outputnames(DDNStructure(it.mdp))
    out = gen(DDNOut(on), it.mdp, s, a, it.rng)
    nt = merge(namedtuple(on, out), (t=t, s=s, a=a, action_info=ai))
    return (out_tuple(it, nt), (t+1, nt.sp))
end

struct POMDPSimIterator{SPEC, M<:POMDP, P<:Policy, U<:Updater, RNG<:AbstractRNG, B, S}
    pomdp::M
    policy::P
    updater::U
    rng::RNG
    init_belief::B
    init_state::S
    max_steps::Int
end
function POMDPSimIterator(spec::Union{Tuple,Symbol}, pomdp::POMDP, policy::Policy, up::Updater, rng::AbstractRNG, init_belief, init_state, max_steps::Int)
    return POMDPSimIterator{spec,
                            typeof(pomdp),
                            typeof(policy),
                            typeof(up),
                            typeof(rng),
                            typeof(init_belief),
                            typeof(init_state)}(pomdp,
                                                policy,
                                                up,
                                                rng,
                                                init_belief,
                                                init_state,
                                                max_steps)
end

Base.IteratorSize(::Type{<:POMDPSimIterator}) = Base.SizeUnknown()

function Base.iterate(it::POMDPSimIterator, is::Tuple{Int,S,B} = (1, it.init_state, it.init_belief)) where {S,B}
    if isterminal(it.pomdp, is[2]) || is[1] > it.max_steps
        return nothing
    end
    t = is[1]
    s = is[2]
    b = is[3]
    a, ai = action_info(it.policy, b)
    on = outputnames(DDNStructure(it.pomdp))
    out = gen(DDNOut(on), it.pomdp, s, a, it.rng)
    outnt = namedtuple(on, out)
    bp, ui = update_info(it.updater, b, a, outnt.o)
    nt = merge(outnt, (t=t, b=b, s=s, a=a, action_info=ai, bp=bp, update_info=ui))
    return (out_tuple(it, nt), (t+1, nt.sp, nt.bp))
end

@generated function out_tuple(it::Union{MDPSimIterator{spec}, POMDPSimIterator{spec}}, all::NamedTuple) where spec
    # the only reason this is generated is to check for :ai and :ui - can get rid of in v0.4
    newspec = Meta.quot(fixdeps(spec))
    quote
        if isa($newspec, Tuple)
            return NamedTupleTools.select(all, $newspec)
        else
            @assert isa(spec, Symbol) "Invalid specification: $spec is not a Symbol or Tuple."
            return all[spec]
        end
    end
end

# XXX can get rid of in v0.4
function fixdeps(tpl::Tuple)
    fixed = []
    for s in tpl
        if s == :ai && !(:action_info in tpl)
            @warn("Automatically switching :ai to :action_info. To disable this switch (e.g. if you have an :ai node in your DDN), also include :action_info in your output spec.")
            push!(fixed, :action_info)
        elseif s == :ui && !(:update_info in tpl)
            @warn("Automatically switching :ui to :update_info. To disable this switch (e.g. if you have an :ui node in your DDN), also include :update_info in your output spec.")
            push!(fixed, :update_info)
        else
            push!(fixed, s)
        end
    end
    return tuple(fixed...)
end
# if there's just one symbol, don't worry about checking for :ai and :ui
fixdeps(s::Symbol) = s

convert_spec(spec, T::Type{M}) where {M<:POMDP} = convert_spec(spec, union(Set(nodenames(DDNStructure(T))), Set(tuple(:bp, :b, :action_info, :update_info, :t))))
convert_spec(spec, T::Type{M}) where {M<:MDP} = convert_spec(spec, union(Set(nodenames(DDNStructure(T))), Set(tuple(:action_info, :t))))

function convert_spec(spec, recognized::Set{Symbol})
    conv = convert_spec(spec)
    convtpl = isa(conv, Tuple) ? conv : tuple(conv)
    for s in convtpl
        if s == :ai && !(:action_info in convtpl)
            @warn("Using :ai to access the action info in a history is deprecated. Use :action_info instead.") # XXX get rid of in v0.4 or greater
        elseif s == :ui && !(:update_info in convtpl)
            @warn("Using :ui to access the update info in a history is deprecated. Use :update_info instead.") # XXX get rid of in v0.4 or greater
        elseif !(s in recognized)
            @warn("uncrecognized symbol $s in step iteration specification $spec.")
        end
    end
    return conv
end

function convert_spec(spec::String)
    syms = spec |> x->strip(x,['(',')']) |> x->split(x,',') |> x->strip.(x) |> x->Symbol.(x)
    if length(syms) == 1
        return Symbol(first(syms))
    else
        return tuple(syms...)
    end
end

function convert_spec(spec::Tuple)
    for s in spec
        @assert isa(s, Symbol)
    end
    return spec
end

convert_spec(spec::Symbol) = spec

default_spec(m::MDP) = tuple(nodenames(DDNStructure(m))..., :t, :action_info)
default_spec(m::POMDP) = tuple(nodenames(DDNStructure(m))..., :t, :action_info, :b, :bp, :update_info)

"""
    stepthrough(problem, policy, [spec])
    stepthrough(problem, policy, [spec], [rng=rng], [max_steps=max_steps])
    stepthrough(mdp::MDP, policy::Policy, [init_state], [spec]; [kwargs...])
    stepthrough(pomdp::POMDP, policy::Policy, [up::Updater, [initial_belief, [initial_state]]], [spec]; [kwargs...])

Create a simulation iterator. This is intended to be used with for loop syntax to output the results of each step *as the simulation is being run*.

Example:

    pomdp = BabyPOMDP()
    policy = RandomPolicy(pomdp)

    for (s, a, o, r) in stepthrough(pomdp, policy, "s,a,o,r", max_steps=10)
        println("in state \$s")
        println("took action \$o")
        println("received observation \$o and reward \$r")
    end

The optional `spec` argument can be a string, tuple of symbols, or single symbol and follows the same pattern as [`eachstep`](@ref) called on a `SimHistory` object.

Under the hood, this function creates a `StepSimulator` with `spec` and returns a `[PO]MDPSimIterator` by calling simulate with all of the arguments except `spec`. All keyword arguments are passed to the `StepSimulator` constructor.
"""
function stepthrough end # for documentation

function stepthrough(mdp::MDP,
                     policy::Policy,
                     spec::Union{String, Tuple, Symbol}=default_spec(mdp);
                     kwargs...)
    sim = StepSimulator(spec; kwargs...)
    return simulate(sim, mdp, policy)
end

function stepthrough(mdp::MDP{S},
                     policy::Policy,
                     init_state::S,
                     spec::Union{String, Tuple, Symbol}=default_spec(mdp);
                     kwargs...) where {S}
    sim = StepSimulator(spec; kwargs...)
    return simulate(sim, mdp, policy, init_state)
end

function stepthrough(pomdp::POMDP, policy::Policy, args...; kwargs...)
    spec_included=false
    if !isempty(args) && isa(last(args), Union{String, Tuple, Symbol})
        spec = last(args)
        spec_included = true
        if spec isa statetype(pomdp) && length(args) == 3
            error("Ambiguity between `initial_state` and `spec` arguments in stepthrough. Please explicitly specify the initial state and spec.")
        end
    else
        spec = default_spec(pomdp)
    end
    sim = StepSimulator(spec; kwargs...)
    return simulate(sim, pomdp, policy, args[1:end-spec_included]...)
end


text:
@with_kw struct TradeRegulations{T <: AbstractFloat}
    lot::T = 0.1                  # Trade lot in MWh
    dayaheadfee::T = 0.04         # Fee for day ahead trades (EUR/MWh)
    intradayfee::T = 0.11         # Fee for intraday trades (EUR/MWh)
    blocklimit::T = 500.0         # Maximum volume (MWh) in block orders
    blockminlength::Int = 3       # Minimum amount of consecutive hours in block orders
    nblockorders::Int = 50        # Maximum number of block orders in portfolio
    imbalancelimit::T = 200.0     # Maximum volume (MWh) that can be settled in the intraday market
    lowerorderlimit::T = -500.0   # Lower technical order price limit (EUR)
    upperorderlimit::T = 3000.0   # Upper technical order price limit (EUR)
end
NordPoolRegulations() = TradeRegulations{Float64}()

struct DayAheadData{T <: AbstractFloat} <: AbstractModelData
    hydrodata::HydroPlantCollection{T,2}
    water_value::PolyhedralWaterValue{T}
    regulations::TradeRegulations{T}
    intraday_trading::Bool
    penalty_percentage::Float64
    use_blockbids::Bool
    bidlevels::Vector{Vector{T}}

    function DayAheadData(plantdata::HydroPlantCollection{T,2},
                          water_value::PolyhedralWaterValue{T},
                          regulations::TradeRegulations{T},
                          bidlevels::Vector{Vector{T}};
                          intraday_trading::Bool = true,
                          penalty_percentage::Float64 = 0.0,
                          use_blockbids::Bool = true) where T <: AbstractFloat
        length(bidlevels) == 24 || error("Supply exactly 24 bidlevel sets")
        all(length.(bidlevels) .== length(bidlevels[1])) || error("All bidlevel sets must be of the same length.")
        return new{T}(plantdata, water_value, regulations, intraday_trading, penalty_percentage, use_blockbids, bidlevels)
    end
end

function DayAheadData(plantfilename::String, watervalue_filename, bidlevelsets::AbstractMatrix; intraday_trading = true, penalty_percentage = 0.0, use_blockbids = true)
    size(bidlevelsets, 2) == 24 || error(" ")
    plantdata = HydroPlantCollection(plantfilename)
    water_value = PolyhedralWaterValue(watervalue_filename)
    regulations = NordPoolRegulations()
    bidlevels = [bidlevelsets[:,i] for i = 1:size(bidlevelsets, 2)]
    for i in eachindex(bidlevels)
        prepend!(bidlevels[i], regulations.lowerorderlimit)
        push!(bidlevels[i], regulations.upperorderlimit)
    end
    DayAheadData(plantdata, water_value, regulations, bidlevels, intraday_trading = intraday_trading, penalty_percentage = penalty_percentage, use_blockbids = use_blockbids)
end

function DayAheadData(plantfilename::String, watervalue_filename, bidlevels::AbstractVector; intraday_trading = true, use_blockbids = true)
    plantdata = HydroPlantCollection(plantfilename)
    water_value = PolyhedralWaterValue(watervalue_filename)
    regulations = NordPoolRegulations()
    prepend!(bidlevels, regulations.lowerorderlimit)
    push!(bidlevels, regulations.upperorderlimit)
    DayAheadData(plantdata, water_value, regulations, fill(bidlevels, 24), intraday_trading = intraday_trading, use_blockbids = use_blockbids)
end


text:
@doc raw"""
    SoftMaxLikelihood(num_class::Int) -> MultiClassLikelihood

## Arguments
- `num_class::Int` : Total number of classes

    SoftMaxLikelihood(labels::AbstractVector) -> MultiClassLikelihood

## Arguments
- `labels::AbstractVector` : List of classes labels

Multiclass likelihood with [Softmax transformation](https://en.wikipedia.org/wiki/Softmax_function):

```math
p(y=i|\{f_k\}_{k=1}^K) = \frac{\exp(f_i)}{\sum_{k=1}^K\exp(f_k)}
```

There is no possible augmentation for this likelihood
"""
SoftMaxLikelihood(x) = MultiClassLikelihood(SoftMaxLink(), x)

implemented(::MultiClassLikelihood{<:SoftMaxLink}, ::MCIntegrationVI) = true

Base.show(io::IO, ::SoftMaxLink) = print(io, "SoftMax Link")

function grad_samples!(
    model::AbstractGPModel{T,<:MultiClassLikelihood{<:SoftMaxLink}},
    samples::AbstractMatrix{T},
    opt_state,
    y,
    index,
) where {T}
    grad_μ = zeros(T, n_latent(model))
    grad_Σ = zeros(T, n_latent(model))
    num_sample = size(samples, 1)
    samples .= mapslices(StatsFuns.softmax, samples; dims=2)
    @inbounds for i in 1:num_sample
        s = samples[i, y][1]
        @views g_μ = grad_softmax(samples[i, :], y) / s
        grad_μ += g_μ
        @views h = diaghessian_softmax(samples[i, :], y) / s
        grad_Σ += h - abs2.(g_μ)
    end
    for k in 1:n_latent(model)
        opt_state[k].ν[index] = -grad_μ[k] / num_sample
        opt_state[k].λ[index] = grad_Σ[k] / num_sample
    end
end

function log_like_samples(
    ::AbstractGPModel{T,<:MultiClassLikelihood{<:SoftMaxLink}},
    samples::AbstractMatrix,
    y::BitVector,
) where {T}
    num_sample = size(samples, 1)
    return mapslices(logsumexp, samples; dims=2) / num_sample
end

function grad_softmax(s::AbstractVector{<:Real}, y)
    return (y - s) * s[y][1]
end

function diaghessian_softmax(s::AbstractVector{<:Real}, y)
    return s[y][1] * (abs2.(y - s) - s .* (1 .- s))
end

function hessian_softmax(s::AbstractVector{T}, y) where {T}
    m = length(s)
    i = findfirst(y)
    hessian = zeros(T, m, m)
    for j in 1:m
        for k in 1:m
            hessian[j, k] =
                s[i] * ((δ(i, k) - s[k]) * (δ(i, j) - s[j]) - s[j] * (δ(j, k) - s[k]))
        end
    end
    return hessian
end


text:
"""
  IsingJ

  ver 0.2.2 Released Sep 2020
  This code is distributed under the constitution
  of MIT LICENSE.

  module Initial

  2019/04/24  Released by PearCandy                 """

module Params
    mutable struct Variables
        site_size::Int64
        exch_j::Float64
        lower_temp::Float64
        upper_temp::Float64
        ext_field::Float64
        step_num::Int64
        sweep_num::Int64
        warm_up::Int64
    end
end

 """ Arguments
     N_site    : number of site
     j0_val    : temperature
     low_t     : minimum temperature
     high_t    : maximum temperature
     ext       : external field
     inum_t    : number of temperature steps
     sweep_num : sweep number
     warm_up   : number of "warming-up" iteration  """

site_size=64
exch_j=1.0
lower_temp=0.1
upper_temp=5.0
ext_field=0.0
step_num=30
sweep_num=2000
warm_up=2000

#using .Params
#params=Params.Variables(site_size,exch_j,lower_temp,upper_temp,
#                        ext_field,step_num,sweep_num,warm_up)


text:
"""
    src

A general module name. In this demo, it contains:
2. DiscreteMarkovChain, a type to define a discrete Markov chain
3. stationary_distribution, a function to solve stationary distribution of a `DiscreteMarkovChain`
4. approx_ar1, approximates an AR(1) process with a discrete Markov Chain
5. golden_section, Golden Section search
6. golden_section_improved, improved 3-point Golden Section search
7. DiscreteMarkovChain, discrete state Markov Chain
"""
module src
    import Distributions
    # ----------------------
    export approx_ar1
    export DiscreteMarkovChain, stationary_distribution
    export DiscreteTimeAR1
    export golden_section, golden_section_improved

# ------------------ TYPE: DISCRETE MARKOV CHAIN
"""
    DiscreteMarkovChain( P::Matrix{T} where T <: Real ; w::Vector = Array(1:size(T,1)) )

immutable type for discrete markov chains,
including:
1. K::Int, the number of states
2. w::Vector, state values, can be anything
2. P::Matrix{T} where T <: Real, transition matrix, square matrix, P[i,j], prob from state i to j
"""
mutable struct DiscreteMarkovChain <: Any
    K::Int  # number of states
    w::Vector # state values, can be anything
    P::Matrix{T} where T <: Real  # transition matrix
    function DiscreteMarkovChain( P::Matrix{T} where T <: Real ; w::Vector = Array(1:size(P,1)) )
        (size(P)[1] == size(P)[2])  ?  nothing  :  throw(ErrorException("non-square transition matrix received"))
        all(isapprox.( sum(P,dims = 2), 1.0))  ?  nothing  : throw(ErrorException("row-sum == 1 transition matrix received!"))
        # adjust P, making row-sum EXACT 1
        local newP = copy(P)
        for x in 1:size(P,1)
            newP[x,:] ./= sum(newP[x,:])
        end
        return new( size(P,1), w, newP )
    end
end
# -----
"""
    stationary_distribution( mc::DiscreteMarkovChain )

get the stationary distribution of a given discrete Markov Chain.
returns a `Vector{Float64}`, elements are in the order of `mc.w`.
"""
function stationary_distribution( mc::DiscreteMarkovChain )
    local tmpP = copy(mc.P)
    while true
        (sum(abs.(tmpP .- tmpP^2)) < 1E-8) ? break : nothing
        tmpP ^= 2
    end # while
    return tmpP[1,:]::Vector{Float64}
end # stationary_distribution

abstract type AbstractARProcess <: Any end
# ------------------------------- AR(1) process
"""
    DiscreteTimeAR1( ρ::Real, σ::Real ; Z̄::Real = 0.0, Z0::Union{Real,Distributions.UnivariateDistribution} = 0.0 )

An AR(1) process defined as: `` Z_{t+1} = (1-\\rho)\\cdot \\bar{Z} + \\rho Z_t + \\epsilon_t ``.
where `ρ` is the auto-reg coefficient, `Z̄` is the mean value, `σ` is the std of the error term ``\\epsilon_t``
which follows i.i.d. `N(0,σ)`, and `Z0` is the initial value/distribution of this process.
It can be a constant, or an instance of one of the subtypes of `Distributions.UnivariateDistribution`, e.g. `Normal()`.
"""
struct DiscreteTimeAR1 <: AbstractARProcess
    ρ::Real
    Z̄::Real
    Z0::Union{Real,Distributions.UnivariateDistribution}
    σ::Real
    # ---------------
    function DiscreteTimeAR1( ρ::Real, σ::Real ; Z̄::Real = 0.0, Z0::Union{Real,Distributions.UnivariateDistribution} = 0.0 )
        @assert(isfinite(ρ), "rho must be finite"); @assert(isfinite(Z̄), "the mean value must be finite")
        @assert(0<σ<Inf, "sigma must be non-negative and finite")
        isa(Z0, Real) ?  @assert(isfinite(Z0), "Z0 must be finite") : nothing
        # -----
        return new(ρ,Z̄,Z0,σ)
    end # AR1()
end # AR1

# ------------------ APPROXIMATE AR(1) WITH MARKOV CHAIN
"""
    approx_ar1(ρ::Real, σ::Real ; m::Int = 9, λ::Int = 3, Z̄::Real = 0.0 )

approximates a zero-mean AR(1) process with a discrete finite-state Markov Chain
whose states are equally divided.
The zero-mean AR(1) process is defined as:
`` Z_{t+1} = (1-\\rho)\\cdot \\bar{Z} + \\rho Z_t + \\epsilon_t ``,
where `ρ` is the auto-regression coefficient, `σ` is the standard error of ``\\epsilon_t``,
`m` is the number of the new Markov chain's states, `Z̄` is the mean of `Zt`.
`λ>0` is the absolute maximum deviation of the Markov chain's states,
it is defined as `max(z) = λσ`, where `z` is the state value of the new Markov Chain.

returns an instance of `DiscreteMarkovChain`.

This algorithm is from `(Tauchen, 1986)` and developed by `(Heer & Maussner, 2015)`.
In most cases, `m=9` and `λ=3` are enough to approximate an AR(1).
"""
function approx_ar1(ρ::Real, σ::Real ; m::Int = 9, λ::Real = 3, Z̄::Real = 0.0 )
    # assertions
    abs(ρ) < 1  ?  nothing : @warn("the AR(1) process is not stationary becuase abs(ρ) >=1")
    @assert(0 < σ < Inf, "the standard error must be greater than 0 and finite")
    @assert(1 < m, "there must be at least two states to make a Markov chain")
    @assert(0 < λ < Inf, "the absolute deviation must be greater than 0 and finite")
    @assert(isfinite(Z̄), "the mean value of Zt must be finite")
    # 1. states
    local w::Vector = zeros(m)
    w[1] = -λ * σ / √(1-ρ^2)  # minimum state value (lower bound, finite)
    local steplen::Real = -2 * w[1] / (m-1)  # step length between two states
    for x in 2:m; w[x] = w[x-1] + steplen; end;  # fill
    # 2. transition matrix
    local P::Matrix = zeros(m,m)
    local tmpNorm = Distributions.Normal() # std norm distrib
    for x in 1:m
        P[x,1] = Distributions.cdf(tmpNorm, (w[1]-ρ*w[x])/σ + 0.5 * steplen / σ )
        if m>2
            for y in 2:(m-1)
                P[x,y] = Distributions.cdf(tmpNorm, (w[y]-ρ*w[x])/σ + 0.5 * steplen / σ ) -
                    Distributions.cdf(tmpNorm, (w[y]-ρ*w[x])/σ - 0.5 * steplen / σ )
            end # for
        end # if
        P[x,m] = 1 - sum(P[x,1:(m-1)])
    end # x
    return DiscreteMarkovChain(P,w = w .+ Z̄)::DiscreteMarkovChain # add mean value to state
end # approx_ar1
# -------------------
approx_ar1(ar1::DiscreteTimeAR1 ; m::Int = 9, λ::Real = 3) = begin
    return approx_ar1(ar1.ρ, ar1.σ, m = m, λ = λ, Z̄ = ar1.Z̄)::DiscreteMarkovChain
end # approx_ar1

# ----------------------- Linear interpolation & expolation
"""
    linear_function_interpolation( x0::Real, xsamples::Vector{T} where T <: Real, fsamples::Vector{T} where T <: Real )

linear interpolation on one-parameter function `f(x)` which is defined by a vector of sample `x` points
and a vector of sample `f(x)` corresponding to sample `x`;
`xsamples` should be increasingly sorted.
returns the function value evaluated at the given point `x0`.
if `x0<findmin(x)` or `x0>findmax(x)`, this function uses the linear expolation
of minimum/maximum sample `x` and the second minimum/maximum sample `x`.
"""
function linear_function_interpolation( x0::Real, xsamples::Vector{T} where T <: Real, fsamples::Vector{T} where T <: Real )
    issorted(xsamples) ? nothing : throw(ErrorException("xsamples should be pre-sorted"))
    local N = length(fsamples); local xmax = xsamples[end]; local xmin = xsamples[1]
    local x0_location = (x0 - xmin) / (xmax - xmin) * (N - 1) + 1  # locating x0
    local nearest_xsample_loc = floor(Int,x0_location)  # the nearest sample x point (on the left side) to interpolate x0
    local dist_from_nearest_xsample = x0_location - nearest_xsample_loc # the fractional distance between x0 and the nearest x point
    # NOTE: relationship: x[near] < x0 < x[near+1]

    # case: if touching the upper bound
    if (x0 > xmax) | (nearest_xsample_loc == N)  # to prevent BoundsError at index N
        res = fsamples[end]
        return res::Float64
    # case: if touching the lower bound
    elseif (x0 <= xmin) | (nearest_xsample_loc == 1)
        res = fsamples[1] - (1 - x0_location) * ( fsamples[2] - fsamples[1] )
        return res::Float64
    # case: normal case
    else
        res = (1 - dist_from_nearest_xsample) * fsamples[nearest_xsample_loc] +
            dist_from_nearest_xsample * fsamples[nearest_xsample_loc + 1]
        return res::Float64
    end # if

end # linear_function_interpolation

# ------------------------------ conventional Golden Section for 1-parameter function
"""
    golden_section(f::Function, lb::Real, ub::Real ; atol::Real = 1E-8, maxiter::Int = 5000)

conventional Golden Section search for MINIMIZATION problem.
`f` must only receive one `Real` type location parameter;
`lb` is the finite lower bound of searching; `ub` is the finite upper bound of searching;
`atol` is the absolute error to converge,
and `maxiter` sets the maximum loops to prevent dead loop.
returns a `Tuple` (in order) consisting of optimal x, optimal `f` value, a flag indicating if converged,
and in which round the algorithm ends.

testing code:
```
g(x::Real) = (x-1)^2
BenchmarkTools.@benchmark tmp = src.golden_section(g,-1.0,3.0,atol=1E-12)
```

performance tips:
```
BenchmarkTools.Trial:
  memory estimate:  96 bytes
  allocs estimate:  1
  --------------
  minimum time:     144.933 ns (0.00% GC)
  median time:      150.576 ns (0.00% GC)
  mean time:        181.232 ns (7.65% GC)
  maximum time:     71.892 μs (99.47% GC)
  --------------
  samples:          10000
  evals/sample:     852
```
"""
function golden_section(f::Function, lb::Real, ub::Real ; atol::Real = 1E-8, maxiter::Int = 5000)
    @assert(-Inf < lb < ub < Inf, "lb must be less than ub, both must be finite")
    local goldnum = 0.618033988749895  # the golden number
    local bounds = [lb,ub]  # changable bounds to update
    local left_x = bounds[1] + (1-goldnum) * (bounds[2] - bounds[1])  # left trial
    local right_x = bounds[1] + goldnum * (bounds[2]-bounds[1])  # right trial
    # initial evaluation (NOTE: using type assertion to test if `f` only returns one `Real` value)
    local left_fval::Real = f(left_x)
    local right_fval::Real = f(right_x)
    # test if both trials are defined & finite (because NaN & Inf cannot be operated rationally)
    if isnan(left_fval) | isinf(left_fval) | isnan(left_fval) | isinf(left_fval)
        throw(DomainError("NaN or Inf function value(s) found in the initial evaluation of Golden section search"))
    end # if
    # begin search
    for j in 1:maxiter
        # check convergency
        if abs(bounds[2]-bounds[1]) < atol
            local finalx::Real = (bounds[1] + bounds[2]) / 2
            return ( finalx, f(finalx), true, j )::Tuple
        end # if
        # if not converge, go on to update
        if left_fval > right_fval
            bounds[1] = left_x; bounds[2] = bounds[2]; left_x = right_x
            right_x = bounds[1] + goldnum * (bounds[2] - bounds[1])
            left_fval = f(left_x)
            right_fval = f(right_x)
        else # i.e. if left_fval <= right_fvals
            bounds[1] = bounds[1]; bounds[2] = right_x; right_x = left_x
            left_x = bounds[1] + (1-goldnum) * (bounds[2] - bounds[1])
            left_fval = f(left_x)
            right_fval = f(right_x)
        end # if
    end # for j
    # if loop normally ends, it means that the algorithm did not converge
    local finalx::Real = (left_x+right_x)/2
    return ( finalx, f(finalx), false, maxiter )::Tuple
end # golden_section

# ---------------------------- improved Golden Section search
"""
    golden_section_improved(f::Function, LLB::Float64, LB::Float64, RB::Float64 ;  TOL::Float64=1E-08, ITER=300)

improved Golden Section search, using a third point x_mid between (x_low,x_high) to improve performance.
designed for one-dim function MAXIMIZATION problem

Par:
    1. f [annonymous func]: target to search (a maximization problem), it MUST HAVE ONLY ONE PARAMETER f(x)
    1. LLB [num]: the (initial) very left bound of searching (initial)
    1. LB [num]: the (initial) left bound of searching
    1. RB [num]: the (initial) right bound of searching
    1. TOL [num]: tolerance, 1E-4 or 1E-5 is enough
    1. ITER [int]: maximum iteration times in searching

Ret:
    1. Xmin [num]: the solution

Depend:
    1. func
"""
function golden_section_improved(f::Function, LLB::Real, LB::Real, RB::Real ;  TOL::Float64=1E-08)
    # -------- DATA PROCESS ----------
    # NOTE: in Julia, the golden numebr is an integrated const (before v0.6), but we redefine a golden number for case in v1.0 and later
    GoldenNumber = 1.618033988749895
    r1 = GoldenNumber - 1; r2 = 1 - r1
    x0 = LLB; x3 = RB # initia-lize bounds
    if abs(RB-LB)<=abs(LB-LLB)
        x1 = LB; x2 = LB+r2*(RB-LB)
    else
        x2 = LB; x1 = LB-r2*(LB-LLB)
    end
    # initialization of function value (and turns a maximization to a minimization)
    f1 = - f(x1); f2 = - f(x2)
    # Searching
    while true
        if f2<f1
            x0=x1 # Update the very lower bound
            x1=x2; x2=r1*x1+r2*x3 # new left golden position
            f1=f2; f2= - f(x2)
        else
            x3=x2; x2=x1; x1=r1*x2+r2*x0; f2=f1; f1= - f(x1)
        end # if
        if abs(x3-x0) > TOL * (abs(x1)+abs(x2))
            break
        end # if
    end # while
    # Post-convergence
    Xmin = f1<=f2 ? x1 : x2

    return Xmin::Float64
end

end # src
#

max_stars_repo_name: Clpr/EconToolBox

max_issues_repo_name: Clpr/EconToolBox

max_forks_repo_name: Clpr/EconToolBox


text:
function trydiffΔcradius(Δc_partition_radius_candidates::Vector{T},
    molecule_names, base_path_JLD, Δcs_max_mixture, hz2ppmfunc, ppm2hzfunc,
    fs, SW, λ0, ν_0ppm, early_exit_part_size, Δcs_max, tol_coherence,
    α_relative_threshold,
    dummy_SSFID::SST)::Tuple{Vector{NMRSpectraSimulator.FIDModelType{T,SST}}, T} where {T <: Real, SST}

    @assert early_exit_part_size > 0
    @assert all(Δc_partition_radius_candidates .> zero(T))

    Δcs_max_mixture = collect( Δcs_max for i = 1:length(molecule_names))

    for Δc_partition_radius in Δc_partition_radius_candidates[1:end-1]

        mixture_params = NMRSpectraSimulator.setupmixtureproxies(molecule_names,
            base_path_JLD, Δcs_max_mixture, hz2ppmfunc, ppm2hzfunc, fs, SW,
            λ0, ν_0ppm, dummy_SSFID;
            tol_coherence = tol_coherence,
            α_relative_threshold = α_relative_threshold,
            Δc_partition_radius = Δc_partition_radius)
        As = mixture_params

        if all( all(displaypartitionsizes(As[n]) .<= early_exit_part_size) for n = 1:length(As) )
            return mixture_params, Δc_partition_radius
        end
    end

    mixture_params = NMRSpectraSimulator.setupmixtureproxies(molecule_names,
    base_path_JLD, Δcs_max_mixture, hz2ppmfunc, ppm2hzfunc, fs, SW, λ0,
    ν_0ppm, dummy_SSFID;
    tol_coherence = tol_coherence,
    α_relative_threshold = α_relative_threshold,
    Δc_partition_radius = Δc_partition_radius_candidates[end])

    return mixture_params, Δc_partition_radius_candidates[end]
end


text:
include("plotting.jl")
include("algConnectivity.jl")
out_fldr = "projectResults/empirical_results_sg/"
fnames = ["initial","final"]

for fname in fnames
    f = open(string(joinpath(out_fldr,fname),".dot"),"r")

    g = SimpleGraphs.SimpleGraph()

    for i in enumerate(eachline(f))
        if i[1] == 1 || eof(f)
            continue
        else
            inew = [parse(Int,x) for x in split(i[2],"--")]
            if length(inew) == 1
                add!(g,inew[1])
            else
                add!(g,inew[1],inew[2])
            end
        end
    end
    close(f)
    clf()
    ac = AlgConnectivity(g,false)
    plot_deg_hist(g,out_fldr,string("d_hist_",fname),ac)
end


text:
using ToggledLH
using Test

foo(x) = 2x;

@testset "ToggledLH" begin
    toggle(false);
    x = 3;
    @toggled (global x = 4)
    @test x == 3
    @toggled (1 == 1.0) (global x = 5)
    @test x == 3

    z = 0;
    @toggled z = foo(3)
    @test z == 0

    @tog_assert (1 == 2.0) "Nothing happens here"

    toggle(true)
    x = 3;
    @toggled (1 == 1.0) (global x = 5)
    @test x == 5
    @toggled (global x = 4)
    @test x == 4
    @toggled (1 == 2.0) (global x = 7)
    @test x == 4

    z = 0;
    @toggled z = foo(3)
    @test z == foo(3)

    z = 0;
    @toggled (1 == 1.0) z = foo(3)
    @test z == foo(3)
    z = 0;
    @toggled (2 == 1.0) z = foo(3)
    @test z == 0

    @toggled (1 == 1.0) begin
        global x = 11
    end
    @test x == 11
    @test_throws AssertionError (@tog_assert 1 == 2.0)
end

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


text:
using DataFrames
using PooledArrays
using BenchmarkTools
using Random

Random.seed!(1)

grouping_benchmarks = BenchmarkGroup()

for n in (1000, 10_000_000),
    k in (10, 10_000),
    x in (rand(1:k, n),
          rand([1:k; missing], n),
          rand(string.(1:k), n),
          rand([string.(1:k); missing], n),
          PooledArray(rand(1:k, n)),
          PooledArray(rand([missing; 1:k], n)))
        df = DataFrame(x=x)

        grouping_benchmarks[n, k, typeof(x), "skipmissing=false"] =
            @benchmarkable groupby($df, :x)

        # Skipping missing values
        grouping_benchmarks[n, k, typeof(x), "skipmissing=true"] =
            @benchmarkable groupby($df, :x, skipmissing=true)

        for l in (10, 10_000),
            y in (rand(1:l, n),
                  rand([1:l; missing], n),
                  rand(string.(1:l), n),
                  rand([string.(1:l); missing], n),
                  PooledArray(rand(1:l, n)),
                  PooledArray(rand([missing; 1:l], n)))
            df.y = y

            grouping_benchmarks[n, (k, l), (typeof(x), typeof(y)),
                                "skipmissing=false"] =
                @benchmarkable groupby($df, [:x, :y])

            # Skipping missing values
            grouping_benchmarks[n, (k, l), (typeof(x), typeof(y)),
                                "skipmissing=true"] =
                @benchmarkable groupby($df, [:x, :y], skipmissing=true)
        end

        if df.x isa PooledArray
            # Empty group which requires adjusting group indices
            replace!(df.x, levels(df.x)[1] => levels(df.x)[2])
            grouping_benchmarks[n, k, typeof(x), "empty group"] =
                @benchmarkable groupby($df, :x)
        end
end

# If a cache of tuned parameters already exists, use it, otherwise, tune and cache
# the benchmark parameters. Reusing cached parameters is faster and more reliable
# than re-tuning `suite` every time the file is included.
paramspath = joinpath(dirname(@__FILE__), "params.json")

if isfile(paramspath)
    loadparams!(grouping_benchmarks, BenchmarkTools.load(paramspath)[1], :evals);
else
    tune!(grouping_benchmarks)
    BenchmarkTools.save(paramspath, params(grouping_benchmarks));
end

grouping_results = run(grouping_benchmarks, verbose=true)
# using Serialization
# serialize("grouping_results.jls", grouping_results)
# leaves(judge(median(grouping_results1), median(grouping_results2)))
# leaves(regressions(judge(median(grouping_results1), median(grouping_results2))))
# leaves(improvements(judge(median(grouping_results1), median(grouping_results2))))