text:
@testset "clean facet attributes" begin
    # consistent `x` and `y` labels
    df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c"], 100), j=rand(["d", "e", "f"], 100))
    plt = data(df) * mapping(:x, :y, col=:i, row=:j)
    aes = compute_axes_grid(Figure(), plt)
    @test consistent_xlabels(aes)
    @test consistent_ylabels(aes)
    @test colwise_consistent_xlabels(aes)
    @test rowwise_consistent_ylabels(aes)
    @test clean_facet_attributes(aes, hidexdecorations=automatic) ==
        (linkxaxes=:all, linkyaxes=:all, hidexdecorations=true, hideydecorations=true)
    @test clean_facet_attributes(aes, hidexdecorations=automatic, linkxaxes=false) ==
        (linkxaxes=:none, linkyaxes=:all, hidexdecorations=false, hideydecorations=true)
    @test clean_facet_attributes(aes, hideydecorations=false, linkyaxes=:minimal) ==
        (linkxaxes=:all, linkyaxes=:rowwise, hidexdecorations=true, hideydecorations=false)

    @test @test_logs(
        (:warn, "Replaced invalid keyword linkxaxes = :rowwise by automatic. " *
            "Valid values are :all, :colwise, :minimal, :none, true, false, or automatic."),
        clean_facet_attributes(aes, hidexdecorations=automatic, linkxaxes=:rowwise)
    ) == (linkxaxes=:all, linkyaxes=:all, hidexdecorations=true, hideydecorations=true)
    @test @test_logs(
        (:warn, "Replaced invalid keyword hidexdecorations = nothing by automatic. " *
            "Valid values are true, false, or automatic."),
        clean_facet_attributes(aes, hidexdecorations=nothing, linkxaxes=:colwise)
    ) == (linkxaxes = :colwise, linkyaxes = :all, hidexdecorations = true, hideydecorations = true)

    # consistent `x` and directionally consistent `y` labels
    df = (x=rand(100), y1=rand(100), y2=rand(100), i=rand(["a", "b", "c"], 100))
    plt = data(df) * mapping(:x, [:y1, :y2], col=:i, row=dims(1))
    aes = compute_axes_grid(Figure(), plt)
    @test consistent_xlabels(aes)
    @test !consistent_ylabels(aes)
    @test colwise_consistent_xlabels(aes)
    @test rowwise_consistent_ylabels(aes)
    @test clean_facet_attributes(aes, hidexdecorations=automatic) ==
        (linkxaxes=:all, linkyaxes=:rowwise, hidexdecorations=true, hideydecorations=true)
    @test clean_facet_attributes(aes, hidexdecorations=automatic, linkxaxes=false) ==
        (linkxaxes=:none, linkyaxes=:rowwise, hidexdecorations=false, hideydecorations=true)
    @test clean_facet_attributes(aes, hideydecorations=false, linkxaxes=:minimal) ==
        (linkxaxes=:colwise, linkyaxes=:rowwise, hidexdecorations=true, hideydecorations=false)

    # facet wrap
    df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c", "d"], 100))
    plt = data(df) * mapping(:x, :y, layout=:i)
    aes = compute_axes_grid(Figure(), plt)
    @test consistent_xlabels(aes)
    @test consistent_ylabels(aes)
    @test colwise_consistent_xlabels(aes)
    @test rowwise_consistent_ylabels(aes)
    @test clean_facet_attributes(aes, hidexdecorations=automatic) ==
        (linkxaxes=:all, linkyaxes=:all, hidexdecorations=true, hideydecorations=true)
    @test clean_facet_attributes(aes, hidexdecorations=automatic, linkxaxes=false) ==
        (linkxaxes=:none, linkyaxes=:all, hidexdecorations=false, hideydecorations=true)
    @test clean_facet_attributes(aes, hideydecorations=false, linkyaxes=:minimal) ==
        (linkxaxes=:all, linkyaxes=:rowwise, hidexdecorations=true, hideydecorations=false)
end

@testset "facet labels" begin
    df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c"], 100), j=rand(["d", "e", "f"], 100))
    plt = data(df) * mapping(:x, :y, col=:i, row=:j)
    fig = Figure()
    aes = compute_axes_grid(fig, plt)
    scales = first(aes).categoricalscales
    ax = first(aes).axis

    labels = col_labels!(fig, aes, scales[:col])
    @test labels[1].text[] == "a"
    @test labels[2].text[] == "b"
    @test labels[3].text[] == "c"
    for (i, label) in enumerate(labels)
        @test label in contents(fig[1, i, Top()])
        @test label.rotation[] == 0
        @test label.padding[] == (0, 0, ax.titlegap[], 0)
        @test label.color == ax.titlecolor
        @test label.font == ax.titlefont
        @test label.textsize == ax.titlesize
        @test label.visible == ax.titlevisible
    end

    labels = row_labels!(fig, aes, scales[:row])
    @test labels[1].text[] == "d"
    @test labels[2].text[] == "e"
    @test labels[3].text[] == "f"
    for (i, label) in enumerate(labels)
        @test label in contents(fig[i, 3, Right()])
        @test label.rotation[] == -π/2
        @test label.padding[] == (ax.titlegap[], 0, 0, 0)
        @test label.color == ax.titlecolor
        @test label.font == ax.titlefont
        @test label.textsize == ax.titlesize
        @test label.visible == ax.titlevisible
    end

    df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c", "d"], 100))
    plt = data(df) * mapping(:x, :y, layout=:i)
    fig = Figure()
    aes = compute_axes_grid(fig, plt)
    scales = first(aes).categoricalscales
    ax = first(aes).axis

    labels = panel_labels!(fig, aes, scales[:layout])
    @test labels[1].text[] == "a"
    @test labels[2].text[] == "b"
    @test labels[3].text[] == "c"
    @test labels[4].text[] == "d"
    for (I, label) in enumerate(labels)
        i, j = fldmod1(I, 2)
        @test label in contents(fig[i, j, Top()])
        @test label.rotation[] == 0
        @test label.padding[] == (0, 0, ax.titlegap[], 0)
        @test label.color == ax.titlecolor
        @test label.font == ax.titlefont
        @test label.textsize == ax.titlesize
        @test label.visible == ax.titlevisible
    end
end

@testset "spanned labels" begin
    df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c"], 100), j=rand(["d", "e", "f"], 100))
    plt = data(df) * mapping(:x => "xlabel", :y => "ylabel", col=:i, row=:j)
    fig = Figure()
    aes = compute_axes_grid(fig, plt)
    ax = first(aes).axis

    label = span_xlabel!(fig, aes)
    @test label.rotation[] == 0.0
    @test label.color == ax.xlabelcolor
    @test label.font == ax.xlabelfont
    @test label.textsize == ax.xlabelsize
    @test label.text[] == "xlabel"
    @test label in contents(fig[3, :, Bottom()])

    label = span_ylabel!(fig, aes)
    @test label.rotation[] == π/2
    @test label.color == ax.ylabelcolor
    @test label.font == ax.ylabelfont
    @test label.textsize == ax.ylabelsize
    @test label.text[] == "ylabel"
    @test label in contents(fig[:, 1, Left()])
end

@testset "linking" begin
    df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c"], 100), j=rand(["d", "e", "f"], 100))
    plt = data(df) * mapping(:x => "xlabel", :y => "ylabel", col=:i, row=:j)
    fg = draw(plt; facet=(linkxaxes=:all, linkyaxes=:rowwise))
    aes = fg.grid
    axs = [ae.axis for ae in aes]
    for c in CartesianIndices(axs)
        i, j = Tuple(c)
        ax = axs[i, j]
        @test ax.xaxislinks == setdiff(axs, [ax])
        @test ax.yaxislinks == setdiff(axs[i, :], [ax])
    end

    fg = draw(plt; facet=(linkxaxes=:colwise, linkyaxes=:all))
    aes = fg.grid
    axs = [ae.axis for ae in aes]
    for c in CartesianIndices(axs)
        i, j = Tuple(c)
        ax = axs[i, j]
        @test ax.xaxislinks == setdiff(axs[:, j], [ax])
        @test ax.yaxislinks == setdiff(axs, [ax])
    end

    fg = draw(plt; facet=(linkxaxes=false, linkyaxes=false))
    aes = fg.grid
    axs = [ae.axis for ae in aes]
    for c in CartesianIndices(axs)
        i, j = Tuple(c)
        ax = axs[i, j]
        @test ax.xaxislinks == Axis[]
        @test ax.yaxislinks == Axis[]
    end

    df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c", "d"], 100))
    plt = data(df) * mapping(:x, :y, layout=:i)
    fg = draw(plt; facet=(linkxaxes=:all, linkyaxes=:rowwise))
    aes = fg.grid
    axs = [ae.axis for ae in aes]
    for c in CartesianIndices(axs)
        i, j = Tuple(c)
        ax = axs[i, j]
        @test ax.xaxislinks == setdiff(axs, [ax])
        @test ax.yaxislinks == setdiff(axs[i, :], [ax])
    end

    df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c", "d", "e"], 100))
    plt = data(df) * mapping(:x, :y, layout=:i)
    fg = draw(plt; facet=(linkxaxes=:all, linkyaxes=:rowwise))
    aes = fg.grid
    axs = [ae.axis for ae in aes]
    for c in CartesianIndices(axs)
        i, j = Tuple(c)
        ax = axs[i, j]
        @test ax.xaxislinks == setdiff(axs, [ax])
        @test ax.yaxislinks == setdiff(axs[i, :], [ax])
    end
    @test isempty(contents(fg.figure[2, 3]))
end

@testset "hidedecorations" begin
    x = rand(500)
    y = rand(500)
    i = rand(["a", "b", "c"], 500)
    j = rand(["d", "e", "f"], 500)
    idxs = @. !((i == "b") & (j == "e")) & !((i == "c") & (j == "f"))
    df = (x=x[idxs], y=y[idxs], i=i[idxs], j=j[idxs])
    plt = data(df) * mapping(:x => "xlabel", :y => "ylabel", col=:i, row=:j)
    axis = (
        xticksvisible=true,
        xminorticksvisible=true,
        xgridvisible=true,
        xminorgridvisible=true,
        yticksvisible=true,
        yminorticksvisible=true,
        ygridvisible=true,
        yminorgridvisible=true,
    )
    fg = draw(plt; axis)
    for i in 1:3, j in 1:3
        ax = content(fg.figure[i, j])
        hidex = i != 3
        hidey = j != 1
        @test ax.xticklabelsvisible[] == !hidex
        @test ax.xticksvisible[] == !hidex
        @test ax.xminorticksvisible[] == !hidex
        @test ax.xgridvisible[] == true
        @test ax.xminorgridvisible[] == true

        @test ax.yticklabelsvisible[] == !hidey
        @test ax.yticksvisible[] == !hidey
        @test ax.yminorticksvisible[] == !hidey
        @test ax.ygridvisible[] == true
        @test ax.yminorgridvisible[] == true
    end

    df = (x=rand(100), y=rand(100), i=rand(["a", "b", "c", "d", "e"], 100))
    plt = data(df) * mapping(:x, :y, layout=:i)
    fg = draw(plt; facet=(hideydecorations=false,), axis)
    aes = fg.grid
    axs = [ae.axis for ae in aes]
    for c in CartesianIndices(axs)
        ax = axs[c]
        i, j = Tuple(c)
        aboveempty = i == 1 && j == 3 # subplot above an empty subplot
        hidex = (i != 2) && !aboveempty
        hidey = false
        @test ax.xticklabelsvisible[] == !hidex
        @test ax.xticksvisible[] == !hidex
        @test ax.xminorticksvisible[] == !hidex
        @test ax.xgridvisible[] == true
        @test ax.xminorgridvisible[] == true

        @test ax.yticklabelsvisible[] == !hidey
        @test ax.yticksvisible[] == !hidey
        @test ax.yminorticksvisible[] == !hidey
        @test ax.ygridvisible[] == true
        @test ax.yminorgridvisible[] == true

        alignmode = aboveempty ? Mixed(bottom=Protrusion(0)) : Inside()
        @test ax.alignmode[] == alignmode
    end
    @test isempty(contents(fg.figure[2, 3]))
end


text:
#using Pkg
#Pkg.API.develop(Pkg.PackageSpec(name="Sherlock", path="/home/caseykneale/.julia/dev/Sherlock"))

using Sherlock
sherlock_UI()

using MonteCarloMeasurements
#Pkg.add("Simple")
#using SimpleWeightedGraphs

d = Detective( Sherlock )
typetype_edges( d )
functiontype_edges( d )

functions( d )
types( d )
abstracttypes( d )
undefined( d )

#sherlockplot(d)
#Plots.default(size = (1000,1000))
inquire( d, :Detective )

colors = [ k => (v, inquire( d, v ))  for (k,v) in d.tag ]
d.tag
d.graph

#inquire( d, Symbol("Thunk") )
#magnify( d, :Detective )
#magnify( d, :NodeManager )
#using Plots, GraphRecipes, Blink, Interact
sherlock_UI()

using Plots
l = @layout [a ; b ]

p1 = scatter( [ 0,1,2,3,4 ],[ 0,0,0,0,0 ], markershape = :rect,
            xlim = [ -0.45, 4.45 ], ylim = [ -0.01, 0.01 ],
            legend = false, axis = nothing, ticks = nothing,
            border = :none, markersize = 30,
            size = (640,96),
            markercolor = [ :lightblue, :lightgreen, :tomato,
                            :orchid1, :grey]
            )
annotate!(p1, [0.0], [0.00], "Function" )
annotate!(p1, [1.0], [0.00], "Type"     )
annotate!(p1, [2.0], [0.00], "Abstract" )
annotate!(p1, [3.0], [0.00], "Untyped"  )
annotate!(p1, [4.0], [0.00], "???"      )

p2 = scatter(randn(100,2))
plot(p1, p2, layout = l)

using GraphRecipes, Plots, Interact
type_in = getfield( Sherlock, :Detective )
plot( type_in, method=:tree, fontsize=10, nodeshape=:rect,
    nodesize = 0.07, nodecolor = :lightgray, color = :black )

# cool = Observable{Any}([])
# cool[] = plot( type_in, method=:tree, fontsize=10, nodeshape=:rect)
#

ext: jl


text:
"""
    cfg_net(;[agent_count=100, growth_rate=4, new_follows=4])

Define configuration parameters for the network.

# Example
```julia-repl
julia>using ABM4OSN

julia>cfg_net()
(agent_count = 100, growth_rate = 4, new_follows = 4)
```

# Arguments
- `agent_count`: How many agents are in the simulation
- `growth_rate`: How many agents to add each tick
- `new_follows`: How many new agent an agent follows each tick

See also: [Config](@ref), [cfg_sim](@ref), [cfg_ot](@ref), [cfg_ag](@ref), [cfg_feed](@ref), [cfg_mech](@ref)
"""
function cfg_net(
    ;
    agent_count::Int64=100,
    growth_rate::Int64=4,
    new_follows::Int64=4
)
    return (
        agent_count=agent_count,
        growth_rate=growth_rate,
        new_follows=new_follows
    )
end

"""
    cfg_sim(;[n_iter=100, max_inactive_ticks=2])

Define configuration parameters for the simulation run.

# Example
```julia-repl
julia>using ABM4OSN

julia>cfg_sim()
(n_iter = 100, max_inactive_ticks = 2, logging = false)
```

# Arguments
- `n_iter`: Number of iterations
- `max_inactive_ticks`: After how many ticks an agent's status is set to inactive
- `Logging`: Determines if extended logging on simulation run is performed

See also: [Config](@ref), [cfg_net](@ref), [cfg_ot](@ref), [cfg_ag](@ref), [cfg_feed](@ref), [cfg_mech](@ref)
"""
function cfg_sim(
    ;
    n_iter::Int64=100,
    max_inactive_ticks::Int64=2,
    logging::Bool=false,
    repcount::Int64=1
)
    return (
        n_iter=n_iter,
        max_inactive_ticks=max_inactive_ticks,
        logging=logging,
        repcount=repcount
    )
end

"""
    cfg_ot(;[like=0.2, share=0.3, backfire=0.4, check_unease=0.3, follow=0.2, unfollow=0.4])

Define thresholds for opinion differences.

# Example
```julia-repl
julia using ABM4OSN

julia>cfg_ot()
(like = 0.2, share = 0.3, backfire = 0.5, check_unease = 0.4, follow = 0.3, unfollow = 0.5)
```

# Arguments
- `like`: Opinion difference threshold for likes
- `share`:  Opinion difference threshold for shares
- `backfire`: Opinion difference for the backfire effect
- `check_unease`: Opinion difference threshold for lowering check regularity
- `follow`: Opinion difference threshold for follows
- `unfollow`: Opinion difference threshold for unfollows

See also: [Config](@ref), [cfg_net](@ref), [cfg_sim](@ref), [cfg_ag](@ref), [cfg_feed](@ref), [cfg_mech](@ref)
"""
function cfg_ot(
    ;
    like::Float64=0.2,
    share::Float64=0.3,
    backfire::Float64=0.5,
    check_unease::Float64=0.4,
    follow::Float64=0.3,
    unfollow::Float64=0.5
)
    return (
        like=like,
        share=share,
        backfire=backfire,
        check_unease=check_unease,
        follow=follow,
        unfollow=unfollow
    )
end

"""
    cfg_ag(;[own_opinion_weight=0.95, check_decrease=0.9, inclin_interact_lambda=log(25), unfollow_rate=0.05, min_input_count=0, mean_desired_input_count=100])

Define agent parameters.

# Example
```julia-repl
julia>using ABM4OSN

julia>cfg_ag()
(own_opinion_weight = 0.99, check_decrease = 0.9, inclin_interact_lambda = 3.2188758248682006, unfollow_rate = 0.2, mean_desired_input_count = 100)
```

# Arguments
- `own_opinion_weight`: Weight of own opinion in relation to other opinions
- `check_decrease`: Decrease factor for check regularity
- `inclin_interact_lambda`: Lambda for exponential distribution for inclination to interact
- `unfollow_rate`: Fraction of agents to unfollow each tick
- `min_input_count`: Minimum input count for an agent to stay in the network
- `mean_desired_input_count`: Mean of desired input count over all agents

See also: [Config](@ref), [cfg_net](@ref), [cfg_sim](@ref), [cfg_ot](@ref), [cfg_feed](@ref), [cfg_mech](@ref)
"""
function cfg_ag(
    ;
    own_opinion_weight::Float64=0.95,
    check_decrease::Float64=0.9,
    inclin_interact_lambda::Float64=log(25),
    unfollow_rate::Float64=0.05,
    min_input_count::Int64=0,
    mean_desired_input_count::Int64=100
)
    return (
        own_opinion_weight=own_opinion_weight,
        check_decrease=check_decrease,
        inclin_interact_lambda=inclin_interact_lambda,
        unfollow_rate=unfollow_rate,
        min_input_count=min_input_count,
        mean_desired_input_count=mean_desired_input_count
    )
end

"""
    cfg_feed(;[feed_size=15, post_decay=0.5])

Define feed parameters.

# Example
```julia-repl
julia>using ABM4OSN

julia>cfg_feed()
(feed_size = 15, post_decay = 0.5)
```

# Arguments
- `feed_size`: length of post feed
- `post_decay`: decay factor for posts in each tick

See also: [Config](@ref), [cfg_net](@ref), [cfg_sim](@ref), [cfg_ot](@ref), [cfg_ag](@ref), [cfg_mech](@ref)
"""
function cfg_feed(
    ;
    feed_size::Int64=15,
    post_decay::Float64=0.5
)
    return (
        feed_size=feed_size,
        post_decay=post_decay
    )
end

"""
    cfg_mech(;[dynamic_net = false, like=true, dislike=true, share=true])

Define mechanics of the social network.

# Example
```julia-repl
julia>using ABM4OSN

julia>cfg_mech()
(dynamic_net = false, like = true, dislike = false, share = true)
```

# Arguments
- `dynamic_net`: Whether or not to include functionality 'dynamic net'
- `like`: Whether or not to include functionality 'like'
- `dislike`: Whether or not to include functionality 'dislike'
- `share`: Whether or not to include functionality 'share'

See also: [Config](@ref), [cfg_net](@ref), [cfg_sim](@ref), [cfg_ot](@ref), [cfg_ag](@ref), [cfg_feed](@ref)
"""
function cfg_mech(
    ;
    dynamic_net::Bool=false,
    like::Bool=true,
    dislike::Bool=false,
    share::Bool=true
)
    return (
        dynamic_net=dynamic_net,
        like=like,
        dislike=dislike,
        share=share
    )
end

"""
    Config(;[network=cfg_net(), simulation=cfg_sim(), opinion_threshs=cfg_ot(), agent_props=cfg_ag(), feed_props=cfg_feed(), mechanics=cfg_mech()])

Provide configuration parameters.

# Example
```julia-repl
julia>using ABM4OSN

julia> Config()
Config((agent_count = 100, growth_rate = 4, new_follows = 4), (n_iter = 100, max_inactive_ticks = 2), (like = 0.2, share = 0.3, backfire = 0.5, check_unease = 0.4, follow = 0.3, unfollow = 0.5), (own_opinion_weight = 0.95, check_decrease = 0.9, inclin_interact_lambda = 3.2188758248682006, unfollow_rate = 0.05, min_input_count = 0, mean_desired_input_count = 100), (feed_size = 15, post_decay = 0.5), (dynamic_net = false, like = true, dislike = false, share = true))
```

# Arguments
- `network`: Tuple of network parameters as created by `cfg_net`
- `simulation`: Tuple of simulation parameters as created by `cfg_sim`
- `opinion_threshs`: Tuple of opinion difference thresholds as created by `cfg_ot`
- `agent_props`: Tuple of agent parameters as created by `cfg_ag`
- `feed_props`: Tuple of feed parameters as created by `cfg_feed`
- `mechanics`: Tuple of mechanics configurations as created by `cfg_mech`

See also: [cfg_net](@ref), [cfg_sim](@ref), [cfg_ot](@ref), [cfg_ag](@ref), [cfg_feed](@ref), [cfg_mechanics](@ref)
"""
struct Config

    network::NamedTuple{
        (
            :agent_count,
            :growth_rate,
            :new_follows
        ),
        NTuple{3, Int64}
    }
    simulation::NamedTuple{
        (
            :n_iter,
            :max_inactive_ticks,
            :logging,
            :repcount
        ),
        <:Tuple{Int64, Int64, Bool, Int64}
    }
    opinion_threshs::NamedTuple{
        (
            :like,
            :share,
            :backfire,
            :check_unease,
            :follow,
            :unfollow
        ),
        NTuple{6, Float64}
    }
    agent_props::NamedTuple{
        (
            :own_opinion_weight,
            :check_decrease,
            :inclin_interact_lambda,
            :unfollow_rate,
            :min_input_count,
            :mean_desired_input_count
        ),
        <:Tuple{Float64, Float64, Float64, Float64, Int64, Int64}
    }
    feed_props::NamedTuple{
        (
            :feed_size,
            :post_decay
        ),
        <:Tuple{Int64, Float64}
    }
    mechanics::NamedTuple{
        (
            :dynamic_net,
            :like,
            :dislike,
            :share
        ),
        NTuple{4, Bool}
    }

    function Config(
        ;
        network=cfg_net(),
        simulation=cfg_sim(),
        opinion_threshs=cfg_ot(),
        agent_props=cfg_ag(),
        feed_props=cfg_feed(),
        mechanics=cfg_mech()
    )
        new(
            network,
            simulation,
            opinion_threshs,
            agent_props,
            feed_props,
            mechanics
        )
    end

end

# suppress output of include()
;


text:
function rankedge{V}(graph::AbstractGraph{V})
    weights = zeros(num_edges(graph))
    for e in edges(graph)
        e_idx = edge_index(e, graph)
        u = source(e, graph)
        v = target(e, graph)
        u_neighbors = Set{V}(out_neighbors(u, graph))
        v_neighbors = Set{V}(out_neighbors(v, graph))
        common_neighbors = intersect(u_neighbors, v_neighbors)
        union_neighbors = union(u_neighbors, v_neighbors)
        union_neighbors = setdiff(union_neighbors, Set{V}({u,v}))
        no_common_neighbors = setdiff(union_neighbors, common_neighbors)
        weights[e_idx] = (length(common_neighbors) + 1) / (length(no_common_neighbors) + 1)
    end
    weights
end

function rankedge1{V}(graph::AbstractGraph{V})
    weights = zeros(num_edges(graph))
    for e in edges(graph)
        e_idx = edge_index(e, graph)
        u = source(e, graph)
        v = target(e, graph)
        u_neighbors = Set{V}(out_neighbors(u, graph))
        v_neighbors = Set{V}(out_neighbors(v, graph))
        common_neighbors = intersect(u_neighbors, v_neighbors)
        union_neighbors = union(u_neighbors, v_neighbors)
        union_neighbors = setdiff(union_neighbors, Set{V}({u,v}))
        no_common_neighbors = setdiff(union_neighbors, common_neighbors)
        weights[e_idx] = (length(common_neighbors) + 1) / (length(no_common_neighbors) + length(common_neighbors) + 1)
    end
    weights
end

function rankedgenocommon{V}(graph::AbstractGraph{V})
    weights = zeros(num_edges(graph))
    for e in edges(graph)
        e_idx = edge_index(e, graph)
        u = source(e, graph)
        v = target(e, graph)
        u_neighbors = Set{V}(out_neighbors(u, graph))
        v_neighbors = Set{V}(out_neighbors(v, graph))
        common_neighbors = intersect(u_neighbors, v_neighbors)
        union_neighbors = union(u_neighbors, v_neighbors)
        union_neighbors = setdiff(union_neighbors, Set{V}({u,v}))
        no_common_neighbors = setdiff(union_neighbors, common_neighbors)
        weights[e_idx] = length(no_common_neighbors)
    end
    weights
end

function rankedgebynei{V}(graph::AbstractGraph{V})
    weights = zeros(num_edges(graph))
    for e in edges(graph)
        e_idx = edge_index(e, graph)
        u = source(e, graph)
        v = target(e, graph)
        u_neighbors = Set{V}(out_neighbors(u, graph))
        v_neighbors = Set{V}(out_neighbors(v, graph))
        common_neighbors = intersect(u_neighbors, v_neighbors)
        weights[e_idx] = length(common_neighbors)
    end
    weights
end


text:
"""
    Simulation(config, init_state, final_state, agent_log, post_log, graph_list)

Provide data structure for a simulation.

# Examples
```julia-repl
julia>using ABM4OSN

julia>Simulation()
Simulation{Config, Any, Any, DataFrame, Any, Array{AbstractGraph}}
```

# Arguments
- `config`: Config object as provided by Config
- `init_state`: Initial graph and agent_list of Simulation
- `final_state`: Final graph and agent_list of Simulation
- `agent_log`: DataFrame that logs the agent states at each simulation tick
- `post_log`: During simulation holds an array of all posts, for output converted to DataFrame
- `graph_list`: Array of graphs which are generated at 10%-intermediate steps of simulation

See also: [Post](@ref), [generate_opinion](@ref), [generate_inlinc_interact](@ref), [generate_check_regularity](@ref), [generate_desired_input_count](@ref)
"""
mutable struct Simulation

	name::String
	batch_name::String
	repnr::Int64
	config::Config
	rng::MersenneTwister
    init_state::Any
    final_state::Any
    agent_log::DataFrame
    post_log::Any
    graph_list::Array{AbstractGraph}

    function Simulation(config=Config(), batch_name="")
        new(
			"",
			batch_name,
			0,
            config,
			MersenneTwister(),
            (nothing, nothing),
            (nothing, nothing),
            DataFrame(),
            DataFrame(),
            Array{AbstractGraph, 1}(undef, 0)
        )
    end

end

# Base.show(io::IO, s::Simulation) = print(
#     io, "Simulation{Config, Any, Any, DataFrame, Any, Array{AbstractGraph}}"
# )

"""
    tick!(state, post_list, tick_nr, config)

Runs a single tick of the simulation and returns the updated state and simulation logs.

# Arguments
- `state`: a tuple of the current graph and agent_list
- `post_list`: List of all published posts in network
- `tick_nr`: Number of current simulation tick
- `config`: Config object as provided by Config()

See also: [log_network](@ref), [simulate!](@ref), [Config](@ref)
"""
function tick!(
    state::Tuple{AbstractGraph, AbstractArray},
    post_list::AbstractArray,
    tick_nr::Int64,
	config::Config
)

    agent_list = state[2]

    for agent_idx in shuffle(1:length(agent_list))

        this_agent = agent_list[agent_idx]

        if this_agent.active && (rand() < this_agent.check_regularity)

            this_agent.inactive_ticks = 0
            update_feed!(state, agent_idx, config)
            update_perceiv_publ_opinion!(state, agent_idx)
            update_opinion!(state, agent_idx, config)

            if config.mechanics.like
                like!(state, agent_idx, config)
            end

            if config.mechanics.dislike
                dislike!(state, agent_idx, config)
            end

            if config.mechanics.share
                share!(state, agent_idx, config)
            end

            if config.mechanics.dynamic_net
                drop_input!(state, agent_idx, config)
				if (
					indegree(state[1], agent_idx)
					< this_agent.desired_input_count
				)
					add_input!(state, agent_idx, post_list, config)
				end
				update_check_regularity!(state, agent_idx, config)
			else
				update_input!(state, agent_idx, config)
			end

            inclin_interact = deepcopy(this_agent.inclin_interact)
            while inclin_interact > 0
                if rand() < inclin_interact
                    publish_post!(state, post_list, agent_idx, tick_nr)
                end
                inclin_interact -= 1.0
            end

        elseif this_agent.active
            this_agent.inactive_ticks += 1
            if (
                this_agent.inactive_ticks > config.simulation.max_inactive_ticks
                && config.mechanics.dynamic_net
            )
                set_inactive!(state, agent_idx, post_list)
            end
        end

    end

    if config.mechanics.dynamic_net
        update_network!(state, config)
    end

    return log_network(state, tick_nr)

end

"""
    run!(simulation=Simulation(); [name="_"])

Creates the initial state, performs and logs simulation ticks and returns the collected data

# Arguments
- `simulation`: Simulation object that provides data structure for simulation results and a config
- `name`: Name of current simulation

See also: [rerun_single!](@ref), [log_network](@ref), [tick!](@ref), [Config](@ref)
"""
function run!(
    simulation::Simulation=Simulation()
    ;
    name::String="_"
)

	if name * ".jld2" in readdir("results")
		raw = load(joinpath("results", name * ".jld2"))
		rep = raw[first(keys(raw))]
	else
		rep = Simulation[]
	end

	for current_rep in 1:simulation.config.simulation.repcount

		simulation.name = name
		simulation.repnr = length(rep) + 1
		simulation.rng = Random.seed!(sum(codeunits(name)) + simulation.repnr)

		agent_list = create_agents(simulation.config)
	    simulation.init_state = (
	        create_network(agent_list, simulation.config), agent_list
	    )
	    state = deepcopy(simulation.init_state)
	    post_log = Array{Post, 1}(undef, 0)

		if simulation.config.simulation.logging
		    simulation.graph_list = Array{AbstractGraph, 1}([simulation.init_state[1]])
		    agent_log = DataFrame(
		        TickNr = Int64[],
		        AgentID = Int64[],
		        Opinion = Float64[],
		        PerceivPublOpinion = Float64[],
		        CheckRegularity = Float64[],
		        InclinInteract = Float64[],
		        DesiredInputCount = Int64[],
		        InactiveTicks = Int64[],
		        Indegree = Int64[],
		        Outdegree = Int64[],
		        ActiveState = Bool[]
		    )
		end

	    for i in 1:simulation.config.simulation.n_iter

			if simulation.config.simulation.logging
				append!(
					agent_log,
					tick!(state, post_log, i, simulation.config))
			else
				tick!(state, post_log, i, simulation.config)
			end

	        if i % ceil(simulation.config.simulation.n_iter / 10) == 0

	            print(".")

				if (
					simulation.config.simulation.logging
					&& simulation.config.mechanics.dynamic_net
				)
					current_network = deepcopy(state[1])
		            rem_vertices!(
		                current_network,
		                [agent.id for agent in state[2] if !agent.active]
		            )
		            push!(simulation.graph_list, current_network)
				end
	        end
	    end

	    simulation.final_state = state
		if simulation.config.simulation.logging
			simulation.agent_log = agent_log
		    simulation.post_log = DataFrame(
		        Opinion = [p.opinion for p in post_log],
		        Weight = [p.weight for p in post_log],
		        Source_Agent = [p.source_agent for p in post_log],
		        Published_At = [p.published_at for p in post_log],
		        Seen = [p.seen_by for p in post_log],
		        Likes = [p.like_count for p in post_log],
		        Dislikes = [p.dislike_count for p in post_log],
		        Reposts = [p.share_count for p in post_log]
		    )
		end

		push!(rep, deepcopy(simulation))

		if !in("results", readdir())
	        mkdir("results")
	    end
    	save(joinpath("results", name * ".jld2"), name, rep)
	end

    print(
        "\n---\nFinished simulation run with the following specifications:\n
        $(simulation.config)\n---\n"
    )

    return rep

end

"""
    run_batch(config_list; batch_desc)

Creates the initial state, performs and logs simulation ticks and returns the collected data

# Arguments
- `config_list`: List of `Config` objects as provided by `Config()`
- `batch_name`: Name of current simulation batch

See also: [run!](@ref), [rerun_single!](@ref), [tick!](@ref), [Config](@ref)
"""
function run_batch(
    config_list::Array{Config, 1}
    ;
    resume_at::Int64=1,
    stop_at::Int64=length(config_list),
    batch_name::String = ""
)
    for i in resume_at:stop_at
        run_nr = lpad(string(i),length(string(length(config_list))),"0")
        run!(
            Simulation(config_list[i], batch_name=batch_name),
            name = (batch_name * "_run$run_nr")
        )
    end
end

"""
    rerun_single!(path)

Repeats a single simulation run with detailed simulation tracking

# Arguments
- `simulation_imported`: Simulation object of the run that shall be repeated
- `logging`: Bool for activating the detailed logging

See also: [run_batch](@ref), [run!](@ref), [tick!](@ref), [Config](@ref)
"""
function rerun_single!(
    simulation_imported::Simulation;
    logging::Bool = true
    )
    simulation = deepcopy(simulation_imported)

    simulation.config = Config(
        network = simulation.config.network,
        simulation = cfg_sim(
			n_iter = simulation.config.simulation.n_iter,
			max_inactive_ticks = simulation.config.simulation.max_inactive_ticks,
            repcount = simulation.config.simulation.repcount,
            logging = logging
        ),
        opinion_threshs = simulation.config.opinion_threshs,
        agent_props = simulation.config.agent_props,
		feed_props = simulation.config.feed_props,
		mechanics = simulation.config.mechanics
    )

	Random.seed!(simulation.rng.seed)
	agent_list = create_agents(simulation.config)
	simulation.init_state = (
		create_network(agent_list, simulation.config), agent_list
	)
	state = deepcopy(simulation.init_state)
	post_log = Array{Post, 1}(undef, 0)

	if simulation.config.simulation.logging
		simulation.graph_list = Array{AbstractGraph, 1}([simulation.init_state[1]])
		agent_log = DataFrame(
			TickNr = Int64[],
			AgentID = Int64[],
			Opinion = Float64[],
			PerceivPublOpinion = Float64[],
			CheckRegularity = Float64[],
			InclinInteract = Float64[],
			DesiredInputCount = Int64[],
			InactiveTicks = Int64[],
			Indegree = Int64[],
			Outdegree = Int64[],
			ActiveState = Bool[]
		)
	end

	for i in 1:simulation.config.simulation.n_iter

		if simulation.config.simulation.logging
			append!(
				agent_log,
				tick!(state, post_log, i, simulation.config))
		else
			tick!(state, post_log, i, simulation.config)
		end

		if i % ceil(simulation.config.simulation.n_iter / 10) == 0

			print(".")

			if (
				simulation.config.simulation.logging
				&& simulation.config.mechanics.dynamic_net
			)
				current_network = deepcopy(state[1])
				rem_vertices!(
					current_network,
					[agent.id for agent in state[2] if !agent.active]
				)
				push!(simulation.graph_list, current_network)
			end
		end
	end

	simulation.final_state = state
	if simulation.config.simulation.logging
		simulation.agent_log = agent_log
		simulation.post_log = DataFrame(
			Opinion = [p.opinion for p in post_log],
			Weight = [p.weight for p in post_log],
			Source_Agent = [p.source_agent for p in post_log],
			Published_At = [p.published_at for p in post_log],
			Seen = [p.seen_by for p in post_log],
			Likes = [p.like_count for p in post_log],
			Dislikes = [p.dislike_count for p in post_log],
			Reposts = [p.share_count for p in post_log]
		)
	end

	if !in("results", readdir())
        mkdir("results")
    end
    save(joinpath("results", simulation.name * "_singlerun.jld2"), simulation.name, simulation)

	return simulation
end

# suppress output of include()
;


text:
"""
Calibrating a component with Optim.jl
"""

using Agtor
using Optim, CSV, DataFrames, Dates
using Statistics

gs_start = (5, 1)

function run_model(zone)
    time_sequence = zone.climate.time_steps

    # Example annual water allocations
    allocs = (surface_water=150.0, groundwater=40.0)

    @inbounds for (idx, dt_i) in enumerate(time_sequence)
        run_timestep!(zone.manager, zone, idx, dt_i)

        # Resetting allocations for example run
        if monthday(dt_i) == gs_start
            update_available_water!(zone, allocs)
        end
    end

    zone_results, field_results = collect_results(zone)

    return zone_results, field_results
end

"""
Generate an objective function
"""
function create_obj_func(tgt_spec, pars)

    # Loading climate data
    climate_data = "$(data_dir)climate/calib_climate_data.csv"
    use_threads = Threads.nthreads() > 1
    climate_seq = DataFrame!(CSV.File(climate_data, threaded=use_threads, dateformat="dd-mm-YYYY"))

    # Create example zone
    climate = Climate(climate_seq)
    farmer = BaseManager("Calibration")
    zone = create(tgt_spec, climate)

    # Load observed wheat production
    hist_wheat_prod = "$(data_dir)crops/Vic_NC_wheat_production.csv"
    obs_wheat = DataFrame!(CSV.File(hist_wheat_prod, threaded=use_threads, dateformat="YYYY", comment="#", types=Dict(1=>Dates.Date)))
    # obs_wheat[!, :Year] = obs_wheat[!, :Year] .- Dates.Year(1)

    obs_wheat = obs_wheat[obs_wheat[!, :Year] .<= Dates.Date("2015-01-01"), :]
    obs_wheat[!, :Average] = (obs_wheat[!, Symbol(" Wheat produced (t)")] ./ obs_wheat[!, Symbol(" Wheat area sown (ha)")])

    function run_func(params=nothing)
        if !isnothing(params)
            tmp = (; zip(map(Symbol, pars[!, :name]), params)...)
            tmp_zone = update_model!(zone, tmp)
        else
            tmp_zone = zone
        end

        zone_results, field_results = run_model(farmer, tmp_zone)

        f1 = field_results["field1"]

        mod_years = map(x -> Year(x), f1.Date)
        avail_years = (mod_years .>= Dates.Year(1990)) .& (mod_years .<= Dates.Year(2015))

        subset = f1[avail_years, :]
        avg_yield = (subset[!, :irrigated_yield] + subset[!, :dryland_yield]) ./ (subset[!, :irrigated_area] + subset[!, :dryland_area])

        rmse = sqrt(mean((avg_yield .- obs_wheat.Average).^2))

        return rmse
    end

    return run_func
end

# run_func(farmer, wheat_zone, obs_wheat)

data_dir = "test/data/"
calib_zone = load_spec("$(data_dir)zones/CalibZone.yml")
tgt_spec = calib_zone[:CalibZone]
agparams = collect_agparams!(tgt_spec, [])

# Find relevant crop parameters for French-Schultz calibration
crop_params = agparams[occursin.("_uncalibrated_wheat", agparams[!, "name"]) .&
                       (occursin.("wue_coef", agparams[!, "name"]) .|
                        occursin.("et_coef", agparams[!, "name"]) .|
                        occursin.("ssm_coef", agparams[!, "name"]) .|
                        occursin.("rainfall_threshold", agparams[!, "name"])), :]

lower = crop_params.min_val
upper = crop_params.max_val
initial_x = crop_params.default

obj_func = create_obj_func(tgt_spec, crop_params)

res = optimize(obj_func, lower, upper, initial_x, SAMIN(), Optim.Options(iterations=10^4))
# 8.770169e-01 (0.877) w/ SAMIN(), no convergence - took 3055 seconds (~51 mins)

@info res
summary(res)

@info Optim.minimizer(res)
@info (; zip(map(Symbol, crop_params[!, :name]), Optim.minimizer(res))...)

# ┌ Info:  * Status: success
# │
# │  * Candidate solution
# │     Final objective value:     8.896884e-01
# │
# │  * Found with
# │     Algorithm:     Fminbox with L-BFGS
# │
# │  * Convergence measures
# │     |x - x'|               = 0.00e+00 ≤ 0.0e+00
# │     |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00
# │     |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
# │     |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
# │     |g(x)|                 = 2.09e-04 ≰ 1.0e-08
# │
# │  * Work counters
# │     Seconds run:   8549  (vs limit Inf)
# │     Iterations:    7
# │     f(x) calls:    1957
# └     ∇f(x) calls:   1957
# [ Info: [379.7167523959479, 100.00000000000001, 7.755910986915787]
# [ Info: (Zone-CalibZone└──CropField-field1└──Crop-uncalibrated_wheat~rainfall_threshold = 379.7167523959479, Zone-CalibZone└──CropField-field1└──Crop-uncalibrated_wheat~et_coef = 100.00000000000001, Zone-CalibZone└──CropField-field1└──Crop-uncalibrated_wheat~wue_coef = 7.755910986915787)

# Fminbox with L-BFGS
# ┌ Info:  * Status: success
# │
# │  * Candidate solution
# │     Final objective value:     8.896766e-01
# │
# │  * Found with
# │     Algorithm:     Fminbox with L-BFGS
# │
# │  * Convergence measures
# │     |x - x'|               = 0.00e+00 ≤ 0.0e+00
# │     |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00
# │     |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
# │     |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
# │     |g(x)|                 = 2.09e-04 ≰ 1.0e-08
# │
# │  * Work counters
# │     Seconds run:   6642  (vs limit Inf)
# │     Iterations:    7
# │     f(x) calls:    2507
# └     ∇f(x) calls:   2507
# [ Info: [380.4056132641329, 100.00000000000001, 7.703705039894828, 0.31532465148791083]
# [ Info: (Zone-CalibZone└──CropField-field1└──Crop-uncalibrated_wheat~rainfall_threshold = 380.4056132641329, Zone-CalibZone└──CropField-field1└──Crop-uncalibrated_wheat~et_coef = 100.00000000000001, Zone-CalibZone└──CropField-field1└──Crop-uncalibrated_wheat~wue_coef = 7.703705039894828, Zone-CalibZone└──CropField-field1└──Crop-uncalibrated_wheat~ssm_coef = 0.31532465148791083)

# ┌ Info:  * Status: success
# │
# │  * Candidate solution
# │     Final objective value:     6.609913e-01
# │
# │  * Found with
# │     Algorithm:     Fminbox with L-BFGS
# │
# │  * Convergence measures
# │     |x - x'|               = 0.00e+00 ≤ 0.0e+00
# │     |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00
# │     |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
# │     |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
# │     |g(x)|                 = 1.07e-01 ≰ 1.0e-08
# │
# │  * Work counters
# │     Seconds run:   353  (vs limit Inf)
# │     Iterations:    6
# │     f(x) calls:    2667
# └     ∇f(x) calls:   2667
# [ Info: [371.4033677239416, 100.00000000000004, 8.97837123677021, 0.39999999999999986]
# [ Info: (Zone-CalibZone└──CropField-field1└──Crop-uncalibrated_wheat~rainfall_threshold = 371.4033677239416,
#    Zone-CalibZone└──CropField-field1└──Crop-uncalibrated_wheat~et_coef = 100.00000000000004,
#    Zone-CalibZone└──CropField-field1└──Crop-uncalibrated_wheat~wue_coef = 8.97837123677021,
#    Zone-CalibZone└──CropField-field1└──Crop-uncalibrated_wheat~ssm_coef = 0.39999999999999986)

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
#=

    This is a port of the R2 model SurveyResultUnbias2.cs

    With the datafiles
    - survey_result_unbias2_population.csv
    - survey_result_unbias2_gender.csv
    - survey_result_unbias2_person_gender.csv
    - survey_result_unbias2_answer.csv

    Note: there is only 5 observations in the answer file, but reading the files makes it look
          more general...

    Result from the R2 model:
    ```
    <0>[0] Mean: 0.590981              bias[0]
    <0>[0] Variance: 0.0427666

    <0>[1] Mean: 0.472561              bias[1]
    <0>[1] Variance: 0.0494528

    <1>[0] Mean: 591.075               votes[0]
    <1>[0] Variance: 43109

    <1>[1] Mean: 944.156               votes[1]
    <1>[1] Variance: 197880
    ```

    This model:

Summary Statistics
  parameters        mean        std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol     Float64    Float64    Float64   Float64     Float64   Float64       Float64

     bias[1]      0.5994     0.2008     0.0020    0.0025   4269.9355    0.9999      385.8259
     bias[2]      0.5022     0.2250     0.0023    0.0035   3757.7959    1.0002      339.5496
    votes[1]    599.4513   201.1412     2.0114    2.5556   4243.2179    0.9999      383.4118
    votes[2]   1004.9705   450.3128     4.5031    6.9752   3788.1458    1.0002      342.2920

Quantiles
  parameters       2.5%      25.0%       50.0%       75.0%       97.5%
      Symbol    Float64    Float64     Float64     Float64     Float64

     bias[1]     0.1927     0.4562      0.6150      0.7582      0.9293
     bias[2]     0.1006     0.3207      0.5015      0.6832      0.9061
    votes[1]   191.9750   456.0000    614.0000    754.2500    931.0000
    votes[2]   201.0000   643.0000   1002.0000   1369.0000   1812.0000

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function survey_result_unbias2(population, person_gender, answer)
    pop_len = length(population)
    answer_len = length(answer)

    bias ~ filldist(Beta(1,1),pop_len)
    votes = tzeros(pop_len)
    for i in 1:pop_len
        votes[i] ~ Binomial(population[i],bias[i])
    end

    ansBias =  Vector{Real}(undef, answer_len) # tzeros(person_gender_len)
    for i in 1:answer_len
        ansBias[i] = bias[person_gender[i]+1]
        answer[i] ~ Bernoulli(ansBias[i])
    end

end

# 1000 2000
population = parse.(Int64,split(readline("survey_result_unbias2_population.csv"),","))
# 0 1
# gender = parse.(Int64,split(readline("survey_result_unbias2_gender.csv"),","))
# 0 1 0 1 0
person_gender = parse.(Int64,split(readline("survey_result_unbias2_person_gender.csv"),","))
# 1 0 1 1 0
answer = parse.(Int64,split(readline("survey_result_unbias2_answer.csv"),","))

model = survey_result_unbias2(population, person_gender, answer)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 1_000)
# chns = sample(model, IS(), 1_000)

# chns = sample(model, NUTS(), 1_000)
# chns = sample(model, HMC(0.1,6), 1_000)

display(chns)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# COV_EXCL_START
module hadoop
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  module common
    const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
    include("HAServiceProtocol_pb.jl")
    include("ProtobufRpcEngine_pb.jl")
    include("RpcHeader_pb.jl")
    include("Security_pb.jl")
    include("IpcConnectionContext_pb.jl")
  end
  module hdfs
    const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
    include("hdfs_pb.jl")
    include("datatransfer_pb.jl")
    include("acl_pb.jl")
    include("xattr_pb.jl")
    include("encryption_pb.jl")
    include("inotify_pb.jl")
    include("ClientNamenodeProtocol_pb.jl")
  end
  module yarn
    const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
    include("yarn_protos_pb.jl")
    include("yarn_service_protos_pb.jl")
    include("applicationmaster_protocol_pb.jl")
    include("applicationclient_protocol_pb.jl")
    include("containermanagement_protocol_pb.jl")
    include("application_history_client_pb.jl")
  end
end
# COV_EXCL_STOP


text:
module ImageHistogramTest

export
    calc_histogram_norm_factor,
    imhistogramGray,
    imhistogramRGB,
    imhistogramRGB3d_16bit,
    imhistogramRGB3d_new2,
    normalize_histogram,
    plot_imhi,
    plot_imhi_GrayRGB,
    plot_imhi_3D

using Images, Colors
using Plots #  necessary to be able include specials function for ploting histograms in 2D & 3D
# using Gnuplot; # this is currently better for 3D plots, but one can load one package only.
# ToDO: check if one can use 2 plot packages if they live in separate submodules. Mmaybe move into submodule 'imhi_plot3D' and 2D into submodule 'imhi_plot2D'
# Worst case:  make individual modules, one for 2D, one for 3D plus a file for common code which will be included from the module files.

#####
# TWA 2018-03-22 : CustomUnitRanges is not usable with julia 0.6.2 for me ; ==> use OffsetArrays
# using CustomUnitRanges: filename_for_zerorange
# include(filename_for_zerorange)

# ==> use OffsetArrays
#  BUT
#     package Plots cannot handle OffsetArrays
#     thus the need for a simple function to convert to common julia array
# flame on
#     and all this just because julia has 1 offset arrays - grrrr
# flame off
#

#=    ##### begin{OffsetArray off}
using OffsetArrays
t_oa = OffsetArray{Int32}(0:255); # defines an array with index 0 .. 255; values are not initialized
summary(t_oa); # gives some info about t_oa
t_oa = OffsetArray(Int32, 0:255);
summary(t_oa); # gives some info about t_oa
t_oa = OffsetArray(Int16, 0:255);
summary(t_oa); # gives some info about t_oa
t_oa = OffsetArray(UInt16, 0:255);
summary(t_oa); # gives some info about t_oa
t_oa = OffsetArray(zeros(Int16, 256), 0:255)
summary(t_oa); # gives some info about t_oa

#####
## From OffsetArrays package Readme:
#   If your package makes use of OffsetArrays, you can also add the following internal convenience definitions:
#   These versions should work for all types.

_size(A::AbstractArray) = map(length, indices(A))
_size(A) = size(A)

_length(A::AbstractArray) = length(linearindices(A))
_length(A) = length(A)

#####
# convert OffsetArray to Julia array
# flame on
#  julia has one major design mistake - array indices start at 1 instead of zero
# flame off
#

function conv_Offset_to_julia_Array(vma::OffsetArray)
    vma[0:end];
end
cOtjA(x) = conv_Offset_to_julia_Array(x); # short name for above
=#     ##### end{OffsetArray off}

#####
# convert sRGB (non-linear) to linearRGB (also used for sRGB to XYZ)
function linearRGB_from_sRGB(v)
    # v <= 0.04045 ? v/12.92 : ((v + 0.055) / 1.055)^2.4
    v <= 0.04045 ? v/12.92 : ((v + 0.055) / 1.055)^2.4
end

#####
# convert linearRGB to sRGB (non-linear) (also used for XYZ to sRGB)
function sRGB_from_linearRGB(v)
    # v <= 0.0031308 ? 12.92v : 1.055v^(1/2.4) - 0.055
    v <= 0.0031308 ? 12.92v : 1.055v^(inv(2.4)) - 0.055
end

#####
# histogram of 2D-array of image data
# only for internal use within module
#=
function _imhistogram(array2D::T) where T<:AbstractArray
    # I don't like fortran order - 1-offset arrays - driving on the left side
    ihist = ZeroRange(256);
    for j=1:size(array2D,2)
        for i=1:size(array2D,1)
            channel = Int(floor(array2D[i,j]));
            ihist[channel] += 1;
        end # i
    end # j

    return ihist;
end # function histogram of gray image
=#
#=
function _imhistogram(array2D::T) where T<:AbstractArray
    # I don't like fortran order - 1-offset arrays - driving on the left side
    ihist = ZeroRange(256);
    ihist[1] = 1; #  dies with 'ERROR: LoadError: indexing not defined for ImageHistogram.ZeroRange{Int64}'
    #=
    for p in array2D
        channel = Int32(floor(p));
        ihist[channel] += 1;
    end # j
    =#
    return ihist;
end # function histogram of gray image
=#
function _imhistogram(array2D::T, ncolors) where T<:AbstractArray
    # I don't like fortran order - 1-offset arrays - driving on the left side
    ncolorsteps = ncolors - 1;
    ihist = zeros(Int32, ncolors); # zeros(UInt16, ncolors) or zeros(Int32, ncolors);
    # ihist = zeros(OffsetArray(UInt16, 0:ncolorsteps));
    for p in array2D
        # channel = Int32(floor(p * ncolorsteps));
        #channel = floor(UInt16, p * ncolorsteps); # using julia-0.6.2 : OK with RGB ; FAIL with Gray
        channel = UInt16(floor(p * ncolorsteps)); # using julia-0.6.2 : OK with Gray and RGB

        ihist[channel+1] += 1; # Bugfix : use 'channel+1' as channel can be 0 if no color. max(channel) is <= 255 for 'N0f8'
        #ihist[channel] += 1; # using OffsetArray the channel of black is a valid index
    end # for p

    return ihist;
end # function histogram of gray image

#####
# Capacity of a FixedPointNumber
#  Helper function, used in several functions; internal; not for export
#  Colors of images are encoded as FixedPointNumber of type N0f8, N6f10, N4f12, ...
#  The part after 'f' gives the number of bits 'nb' used per color.
#  The number of shades per coler are given by 2^nb.
function _capacity_of_FPNtype(array_elem)
    # ToDo: currently assume 8-bit per pixel per color
    # typname = string("t", zero(Gray(imarray[1]))); # string("t", zero(Gray.(imarray)));
    typname = string("t", zero(Gray(array_elem)));
    if searchindex(typname, "N0f8") > 0
        ncolorshades = 2^8;
        typref = N0f8;
    elseif searchindex(typname, "N6f10") > 0
        ncolorshades = 2^10;
        typref = N6f10;
    elseif searchindex(typname, "N4f12") > 0
        ncolorshades = 2^12;
        typref = N4f12;
    elseif searchindex(typname, "N2f14") > 0
        ncolorshades = 2^14;
        typref = N2f14;
    elseif searchindex(typname, "N0f16") > 0
        ncolorshades = 2^16;
        typref = N0f16;
    else
        println("WARNING: unknown FixedPointNumber type used");
        ncolorshades = 369; # artifical value to make the result worse
        typref = N0f16;
    end # if searchindex(typname, "N0f8") > 0

    return ncolorshades, typref;
end # function _capacity_of_FPNtype(array_elem)

#####
# histogram of gray image
# ToDo: simply convert Gray Image to 2d-array and call '_histogram'
#    Gray.(img_col256)+0
#
#geht net: function imhistogram(imarray::AbstractArray{C}) where C<:AbstractGray
function imhistogramGray(imarray::AbstractArray)
    (ncolors, typref) = _capacity_of_FPNtype(Gray(imarray[1]));

    histo_gray = _imhistogram(Gray.(imarray), ncolors);

    return histo_gray;
end # function histogram of gray image

#####
# histogram of rgb image
# returns 3 histograms
#geht net: function imhistogram(imarray::AbstractArray{C}) where C<:AbstractRGB
function imhistogramRGB(imarray::AbstractArray)
    (ncolors, typref) = _capacity_of_FPNtype(red(imarray[1]));

    histo_red   = _imhistogram(red.(imarray), ncolors);
    histo_green = _imhistogram(green.(imarray), ncolors);
    histo_blue  = _imhistogram(blue.(imarray), ncolors);

    return histo_red, histo_green, histo_blue;
end # function histogram of RGB image

function imhistogramRGB3d_16bit(imarray::AbstractArray)
    # previous name of this func:  imhistogramRGB3d_old
    # use this function to support images with more than 8bits per color channel til package 'Colors' supports RGB48.
    # But the function will always support 8bits.  At least for speed tests, ...

    (ncolors, typref) = _capacity_of_FPNtype(red(imarray[1]));

    # 1st version
#    red_vector   = reshape(red.(imarray), :) .* myfactor; # reshape(red.(imarray), :, 1) * myfactor;
#    green_vector = reshape(green.(imarray), :) .* myfactor; # reshape(green.(imarray), :, 1) * myfactor;
#    blue_vector  = reshape(blue.(imarray), :) .* myfactor; # reshape(blue.(imarray), :, 1) * myfactor;
#    col_vector   = reshape(RGB24.(imarray), :); # reshape(RGB24.(imarray), :, 1);

    # make as 1 dim array
    red_vector   = reshape(red.(imarray), :); # reshape(red.(imarray), :, 1) * myfactor;
    green_vector = reshape(green.(imarray), :); # reshape(green.(imarray), :, 1) * myfactor;
    blue_vector  = reshape(blue.(imarray), :); # reshape(blue.(imarray), :, 1) * myfactor;

    # convert to ...
    if ncolors == 2^8 # searchindex(typname, "N0f8") > 0
        # ... 8-bit int
        red_ivector   = reinterpret.(N0f8.(red_vector));
        green_ivector = reinterpret.(N0f8.(green_vector));
        blue_ivector  = reinterpret.(N0f8.(blue_vector));

        col_ivector = (red_ivector .% UInt32) .<< 16 .+ (green_ivector .% UInt32) .<< 8 .+ (blue_ivector .% UInt32);

        # now do unique sort to reduce the number of color values
        # and prepare the data for the color cube
        col_cube_iv   = unique(sort!(col_ivector));
        red_cube_iv   = (col_cube_iv .& 0x00ff0000) .>> 16;
        green_cube_iv = (col_cube_iv .& 0x0000ff00) .>> 8;
        blue_cube_iv  = (col_cube_iv .& 0x000000ff);
    else # all colors here use 16bit for storage, but they may use not the full range.
        # ... 16-bit int ; currently raw images have 10, 12 or 14 bits per color, some image tools work with 16 bit
        # Depending on the used camera, one can have 10bit, 12bit or 14bit pictures.

        if ncolors == 2^10 # N6f10
            red_ivector   = reinterpret.(N6f10.(red_vector));
            green_ivector = reinterpret.(N6f10.(green_vector));
            blue_ivector  = reinterpret.(N6f10.(blue_vector));
        elseif ncolors == 2^12 # N4f12
            red_ivector   = reinterpret.(N4f12.(red_vector));
            green_ivector = reinterpret.(N4f12.(green_vector));
            blue_ivector  = reinterpret.(N4f12.(blue_vector));
        elseif ncolors == 2^14 # N2f14
            red_ivector   = reinterpret.(N2f14.(red_vector));
            green_ivector = reinterpret.(N2f14.(green_vector));
            blue_ivector  = reinterpret.(N2f14.(blue_vector));
        elseif ncolors == 2^16 # N0f16
            red_ivector   = reinterpret.(N0f16.(red_vector));
            green_ivector = reinterpret.(N0f16.(green_vector));
            blue_ivector  = reinterpret.(N0f16.(blue_vector));
        end # if ncolors == 2^10 # N6f10

        col_ivector = (red_ivector .% UInt64) .<< 32 .+ (green_ivector .% UInt64) .<< 16 .+ (blue_ivector .% UInt64);

        # now do unique sort to reduce the number of color values
        # and prepare the data for the color cube
        col_cube_iv   = unique(sort!(col_ivector));
        red_cube_iv   = (col_cube_iv .& 0x0000ffff00000000) .>> 32;
        green_cube_iv = (col_cube_iv .& 0x00000000ffff0000) .>> 16;
        blue_cube_iv  = (col_cube_iv .& 0x000000000000ffff);
    end # if ncolors == 2^8 # searchindex(typname, "N0f8") > 0

    # now this can be plotted as 3D using gnuplot from julia
    # with the command
    #   @gp(splot=true, redv[:], greenv[:], bluev[:], colv[:], "with points pt 13 ps 0.7 lc rgb variable")
    # since Gnuplot.jl v0.2.0
    #   @gsp(redv[7:10:end], greenv[7:10:end], bluev[7:10:end], gen_pcv(colv[7:10:end]), "with points pt 13 ps 0.7 lc rgb variable", xrange=(0,255), yrange=(0,255), zrange=(0,255), xlabel="red", ylabel="green", zlabel="blue", "set border -1", "set tics in mirror", "set grid", "set zticks out mirror", "set grid ztics", "set xyplane at 0.0")

    # to plot this with Plots
    # keep red_cube_iv, green_cube_iv, blue_cube_iv as is
    # convert col_cube_iv with
    #    colv_new=(RGB24.(N0f8.(UInt8.(redv)/255), N0f8.(UInt8.(greenv)/255), N0f8.(UInt8.(bluev)/255)))
    # then this can be plotted with
    #    scatter3d(redv[140000:1:end], greenv[140000:1:end], bluev[140000:1:end], color=colv_new[140000:1:end], ms=1, msw=0, marker=1, xlims=(0,255), ylims=(0,255), bg=:blue)

#    return red_ivector, green_ivector, blue_ivector, col_ivector;
    return red_cube_iv, green_cube_iv, blue_cube_iv, col_cube_iv;
end # function imhistogramRGB3d_16bit of RGB image

function imhistogramRGB3d_new2(imarray::AbstractArray)
    # previous name of this func:  imhistogramRGB3d_new2
    (ncolors, typref) = _capacity_of_FPNtype(red(imarray[1]));

    if (ncolors > 256)
        println("Sorry, this function supports 8-bits per color channel only !!!");
        println("Please call / use function 'imhistogramRGB3d_16bit(imarray::AbstractArray)' !!!");

        # return red_cube_iv, green_cube_iv, blue_cube_iv, col_cube_iv;
        return [0], [0], [0], [0];
    end # if (ncolors > 256)

    # 1st version
#    red_vector   = reshape(red.(imarray), :) .* myfactor; # reshape(red.(imarray), :, 1) * myfactor;
#    green_vector = reshape(green.(imarray), :) .* myfactor; # reshape(green.(imarray), :, 1) * myfactor;
#    blue_vector  = reshape(blue.(imarray), :) .* myfactor; # reshape(blue.(imarray), :, 1) * myfactor;
#    col_vector   = reshape(RGB24.(imarray), :); # reshape(RGB24.(imarray), :, 1);

    # make as 1 dim array
    #red_vector   = reshape(red.(imarray), :); # reshape(red.(imarray), :, 1) * myfactor;
    #green_vector = reshape(green.(imarray), :); # reshape(green.(imarray), :, 1) * myfactor;
    #blue_vector  = reshape(blue.(imarray), :); # reshape(blue.(imarray), :, 1) * myfactor;

    ### the next uses too much mem
    #col_vector   = reshape(hex.(RGB24.(imarray)), :); # println(col_vector);
    #col_cube_us   = unique(sort!(col_vector));
    #col_cube_iv   = reinterpret.(RGB24, parse.(UInt32, string.("0x", col_cube_us)));

    ### another try => Bingo!
    # calling sort() with the special isless() related to RGB24 reduced used mem and time
    # and then use sort!() for another improvement
    col_vector  = reshape(RGB24.(imarray), :); # println(col_vector);
    col_cube_iv = unique(sort!(col_vector, lt = (x, y)->isless(x.color, y.color)));

    red_cube_iv   = red.(col_cube_iv);
    green_cube_iv = green.(col_cube_iv);
    blue_cube_iv  = blue.(col_cube_iv);

    # now this can be plotted as 3D using gnuplot from julia
    # with the command
    #   @gp(splot=true, redv[:], greenv[:], bluev[:], colv[:], "with points pt 13 ps 0.7 lc rgb variable")
    #
    #   much better
    # @gp(splot=true, redv[:], greenv[:], bluev[:], gen_pcv(colv[:]), "with points pt 13 ps 0.7 lc rgb variable", xlabel="red", ylabel="green", zlabel="blue", "set border -1", "set tics in mirror", "set grid", "set zticks out mirror", "set grid ztics", "set xyplane at 0.0")
    #
    # @gp(splot=true, redv[1:10:end], greenv[1:10:end], bluev[1:10:end], gen_pcv(colv[1:10:end]), "with points pt 13 ps 0.7 lc rgb variable", xrange=(0,255), yrange=(0,255), zrange=(0,255), xlabel="red", ylabel="green", zlabel="blue", "set border -1", "set tics in mirror", "set grid", "set zticks out mirror", "set grid ztics", "set xyplane at 0.0")
    #
    # befor plotting with gnuplot, some crafting is necessary
    # 1) convert colv using
    #   gen_pcv(cv24_a)=(pcv24=zeros(length(cv24_a));for i = 1:endof(cv24_a); pcv24[i]=cv24_a[i].color; end;return pcv24)
    # 2) convert red, green, blue parts
    #   redv=redv*255.0;
    #   greenv=greenv*255.0;
    #   bluev=bluev*255.0;
    #
    # damn - still no 3d-plot with @gp.  somewhere there was a significant change:
    #     !!!! the syntax is not accepted any more !!!!
    #     !!!! splot wants a function instead of data !!!!
    #
    # !!!!! Intermediate Solution no longer necessary; new Gnuplot.jl v-0.2.0 !!!!!
    # manually checkout commit  1665b78a54a0b67ce6b61c2b8ebfe0f409a47ae4  from Gnuplot.jl
    #
    # manually checkout commit  56b64fcef797cc337eb2d589edf5e95a9abd37f5  from Gnuplot.jl
    #   works also
    #
    # starting with commit 6a4e95f377e3da76026731a99f2c6a918cd87fe8 from Gnuplot.jl 'splot' stops working
    #
    # !!!!! !!!!!
    #

    # 1st own 3 plot
    # @gp("set hidden3d", "set grid", "splot 'data3d.txt' using 1:2:3 with points pt 13 ps 0.7 lc rgb variable")
    # @gp("set grid; set xlabel 'red'; set ylabel 'green'; set zlabel 'blue'", "splot 'data3d.txt' using 1:2:3 with points pt 13 ps 0.7 lc rgb variable")
    # still not possible to splot data. Worked last year !!!

    # since Gnuplot.jl v0.2.0
    #  @gsp(redv[7:10:end], greenv[7:10:end], bluev[7:10:end], gen_pcv(colv[7:10:end]), "with points pt 13 ps 0.7 lc rgb variable", xrange=(0,255), yrange=(0,255), zrange=(0,255), xlabel="red", ylabel="green", zlabel="blue", "set border -1", "set tics in mirror", "set grid", "set zticks out mirror", "set grid ztics", "set xyplane at 0.0")

#    return red_ivector, green_ivector, blue_ivector, col_ivector;
    return red_cube_iv, green_cube_iv, blue_cube_iv, col_cube_iv;
end # function imhistogramRGB3d_new2 of RGB image

#####
# Helper used in function plot_imhi_3D().  Needed to pass the RGB24 colors to Gnuplot.
# It extracts the RGB24 color value and stores it in the output array.
# This is understood as color to use for the related marker in plots
gen_pcv(cv24_a)=(pcv24=zeros(length(cv24_a));for i = 1:endof(cv24_a); pcv24[i]=cv24_a[i].color; end;return pcv24)

#####
# How to do the 3D plot.
# But encapsulation into one function does not work.
# Thus, it is not possible to use 2 different plotting tools.
# Try using sub-modules and try 'using' and 'import':  one for 2D another one for 3D
# Else make an appropriate example.
#
function plot_imhi_3D(imarray::AbstractArray; bg::Int = 0, range_step::Int = 50)
#    using Gnuplot; # ToDO: maybe move into submodule 'imhi_plot3D' and 2D into submodule 'imhi_plot2D'

    # default to dark theme; for 3D it looks funny, but still better to read
    # Hmm PlotTheme broken since 2018-04-21 ~22:00
#    if (bg == 0)
#        theme(:dark); # good for electronic media
#        GrayShade=:lightgray;
#    else
#        theme(:default); # good for paper
#        GrayShade=:gray71;
#    end

    redv, greenv, bluev, colv = ImageHistogramTest.imhistogramRGB3d_new2(imarray);

    # convert FixedPontNumbers to plain numbers. They are used as coordinates in 3D.
    # ToDo: extend to respect more than 8 bits per color
    max_color_value = 2^8 - 1.0; # need a float as result
    redv *= max_color_value; greenv *= max_color_value; bluev *= max_color_value;

    # the Gnuplot tools want simple RGB24 values encoded as 32bit ints.
    # colv = gen_pcv(colv);

    # do the plot with Plots
    # scatter3d(redv[1:50:end], greenv[1:50:end], bluev[1:50:end], color=colv[1:50:end], markersize=3, marker=:cross)
    # scatter3d(redv[1:range_step:end], greenv[1:range_step:end], bluev[1:range_step:end], color=colv[1:range_step:end], title="Color Cube with 3D Histogram, draft quality", grid = :all, markersize=3, marker=:cross)

    # With the hints of 'mkborregaard' (JuliaPlots member); thanks
    scatter3d(redv[1:range_step:end], greenv[1:range_step:end], bluev[1:range_step:end], color=colv[1:range_step:end], title="Color Cube with 3D RGB-Histogram, draft quality", xlims=(0,255), ylims=(0,255), zlims=(0,255), grid = :all, markerstrokewidth = 0, markersize = 0.1, bg = :black, camera = (60, 60), marker=:circle)  # disable border of marker with 'markerstrokewidth = 0'

    # do the plot with Gnuplot manually in case the arrays are larger than 3000 - 4000 elements each.
    # @gp(splot=true, redv[1:10:end], greenv[1:10:end], bluev[1:10:end], gen_pcv(colv[1:10:end]), "with points pt 13 ps 0.7 lc rgb variable", xrange=(0,255), yrange=(0,255), zrange=(0,255), xlabel="red", ylabel="green", zlabel="blue", "set border -1", "set tics in mirror", "set grid", "set zticks out mirror", "set grid ztics", "set xyplane at 0.0")
    # since Gnuplot.jl v0.2.0
    #  @gsp(redv[7:10:end], greenv[7:10:end], bluev[7:10:end], gen_pcv(colv[7:10:end]), "with points pt 13 ps 0.7 lc rgb variable", xrange=(0,255), yrange=(0,255), zrange=(0,255), xlabel="red", ylabel="green", zlabel="blue", "set border -1", "set tics in mirror", "set grid", "set zticks out mirror", "set grid ztics", "set xyplane at 0.0")

end # function plot_imhi_3D(imarray::AbstractArray; how::Int = 2, bg::Int = 0)

function calc_histogram_norm_factor(ImHi1::AbstractArray, ImHi2::AbstractArray = [0], ImHi3::AbstractArray = [0], ImHi4::AbstractArray = [0])
    # calculat the norm-factor for the histograms.
    # expect up to 4 histograms : Gray, Red, Green, Blue
    # the order does not matter, thus name it 1 2 3 4
    # at minimum, one histogram - the 1st -  with data must be given via 'ImHi1'.

    if length(ImHi1) != 0
        norm_factor_overall = max(maximum(ImHi1), maximum(ImHi2), maximum(ImHi3), maximum(ImHi4));
    else
        norm_factor_overall = 1;
    end

    return norm_factor_overall;
end # function calc_histogram_norm_factor(...)

#=
# No.  In most cases input and output type differ
function normalize_histogram!(imhisto::AbstractArray, ih_nf = 0.0)
    # normalize the histogram
    #   ih_nf := 0 => max of all channels is set to 1
    #   ih_nf > 0 => use as scale factor

    # modification of arguments does not work this way
    # TWA ToDo : Bug or Feature ???
    # norm_factor = maximum(imhisto);
    # imhisto /= norm_factor;
    #
    # only this way
    if (ih_nf != 0)
        norm_factor = ih_nf;
    else
        norm_factor = typeof(ih_nf)(maximum(imhisto));
    end

    for i = 1:endof(imhisto) # ToDo: change julia-0.6.2 notation to julia-0.7.0 notation if the latter is public.
        imhisto[i] /= norm_factor;
    end

    return imhisto;
end # function normalize_histogram!(imhisto_raw::AbstractArray, ih_nf = 0.0)
=#

function normalize_histogram(imhisto_raw::AbstractArray, ih_nf::Real = 0.0)
    # normalize the histogram
    #   ih_nf := 0 => max of all channels is set to 1
    #   ih_nf > 0 => use as scale factor
    #
    # the output is always an array of reals. Thus in most cases different from input.
    # => no normalize_histogram!()-function

    if (ih_nf != zero(ih_nf))
        norm_factor = ih_nf;
    else
        norm_factor = typeof(ih_nf)(maximum(imhisto_raw));
    end

    imhisto_normed = imhisto_raw / norm_factor;

    return imhisto_normed;
end # function normalize_histogram(imhisto_raw::AbstractArray, ih_nf::Real = 0.0)

function combine_color_slots(imhisto_in::AbstractArray, nslots::Int = 0)
    # imhisto_in :== input of a histogram. It is a 1D array
    # nslots :== how much color slots shall be combined
    #      0 : auto (do nothing for histograms of 8 bit colors having a length of 256)
    #          if larger than 1024 scale down.
    #     >0 : combine. Best if a power of 2.

    if (mod(nslots, 2) != 0)
        println("Number of color slots to combine is best if a pwer of 2");
        println("You feed me with something else, be sure you know what you are doing");
    end

    imhisto_out = zeros(imhisto_in);
    len_ih_in = length(imhisto_in);

    if (nslots == 0)
        if (len_ih_in > 1024) # 2^10
            nslots = div(len_ih_in, 1024);
        else
            nslots = 1;
        end
    else
        nslots = abs(nslots);
        if (nslots > len_ih_in)
            println("The given 'nslots' is larger than the length of the histogram => using 1 instead");
        end
    end

    for i = 1:len_ih_in
        imhisto_out[div(i, nslots)] += imhisto_in[i];
    end
end # function combine_color_slots(imhisto_in::AbstractArray, nslots::Int = 0)

# ToDo : => no, array tooooo large
# my_rgb24int = Int32.(floor.(red.(colv)*255)) .<< 16 .+ Int32.(floor.(green.(colv)*255)) .<< 8 .+ Int32.(floor.(blue.(colv)*255))

function plot_imhi(;ihGray_cooked::AbstractArray = [0], ihR_cooked::AbstractArray = [0],
                   ihG_cooked::AbstractArray = [0], ihB_cooked::AbstractArray = [0],
                   how::Int = 4, bg::Int = 0)
    # default to dark theme
    if (bg == 0)
        theme(:dark); # good for electronic media
        GrayShade=:lightgray;
    else
        theme(:default); # good for paper
        GrayShade=:gray71;
    end

    if (how == 1)
        plot(ihGray_cooked, color=GrayShade, w=3, line=:sticks);
        plot!(ihR_cooked, color=:red, w=3);
        plot!(ihG_cooked, color=:green, w=3);
        plot!(ihB_cooked, color=:blue, w=3);
    elseif (how == 2)
        pl_Gray = plot(ihGray_cooked, color=GrayShade, w=3, line=:sticks);
        pl_RGB = plot(ihR_cooked, color=:red, w=3);
        pl_RGB = plot!(ihG_cooked, color=:green, w=3);
        pl_RGB = plot!(ihB_cooked, color=:blue, w=3);

        plot(pl_Gray, pl_RGB, layout=(2,1), legend=false);
    elseif (how == 3)
        nf = calc_histogram_norm_factor(ihR_cooked, ihG_cooked, ihB_cooked); # ToDo : make it optional

        pl_R = plot(ihR_cooked, ylims=(0,nf), color=:red, w=3);
        pl_G = plot(ihG_cooked, ylims=(0,nf), color=:green, w=3);
        pl_B = plot(ihB_cooked, ylims=(0,nf), color=:blue, w=3);

        plot(pl_R, pl_G, pl_B, layout=(3,1), legend=false);
    elseif (how == 4)
        nf = calc_histogram_norm_factor(ihGray_cooked, ihR_cooked, ihG_cooked, ihB_cooked); # ToDo : make it optional

        pl_Gray = plot(ihGray_cooked, ylims=(0,nf), color=GrayShade, w=3, line=:sticks);
        pl_R = plot(ihR_cooked, ylims=(0,nf), color=:red, w=3);
        pl_G = plot(ihG_cooked, ylims=(0,nf), color=:green, w=3);
        pl_B = plot(ihB_cooked, ylims=(0,nf), color=:blue, w=3);

        plot(pl_Gray, pl_R, pl_G, pl_B, layout=(2,2), legend=false);
    end # if (how == 1)
end # function plot_imhi(...)

function plot_imhi_GrayRGB(imarray::AbstractArray; how::Int = 2, bg::Int = 0)
    # bg := 0 => dark theme as background <=> see PlotThemes package for details
    # bg := 1 => default theme as background
    #
    # how := 1 => plot 4 graphs in 1 diagram
    #     := 2 => 2 subplots : Gray + RGB
    #     := 3 => 3 subplots : Red, Green, Blue
    #     := 4 => 4 subplots : Gray, Red, Green, Blue
    #
    # ToDo:  change how & background to keyword args.
    #

    ihR, ihG, ihB = ImageHistogramTest.imhistogramRGB(imarray);
    ihGray        = ImageHistogramTest.imhistogramGray(imarray);

    plot_imhi(ihGray_cooked=ihGray, ihR_cooked=ihR, ihG_cooked=ihG, ihB_cooked=ihB, how=how, bg=bg)

    #=
    # das folgende braucht noch 'using Plots'
    # Plot examples:
    # or using as subplot
    # plot_red = plot(ihr, line=:red, w=2);
    # plot(plot_red, plot_green, plot_blue, plot_gray, layout=(2,2), legend=false)

    # WATCH OUT !!  WATCH OUT !!
    # after switching to use OffsetArrays
    # plot needs another syntax
    #     plot(ihR[0:end], line=:red)
    #

    ploting into a 3d RGB-cube:
    using Images, TestImages ; import ImageHistogram ; img_col256 = testimage("lena_color_256");
    redv, greenv, bluev, colv = ImageHistogram.imhistogramRGB3d(img_col256);
    redv1=redv[1:32:65536]; greenv1=greenv[1:32:65536]; bluev1=bluev[1:32:65536]; colv1=colv[1:32:65536];

    plot using Plots:
    scatter3d(redv1, greenv1, bluev1, xlabel="red", ylabel="green", zlabel="blue"; color=colv1)
    scatter3d(redv[1:50:end], greenv[1:50:end], bluev[1:50:end], color=colv[1:50:end], markersize=3, marker=:cross)

    plot using GR:
    setmarkersize=1; setmarkertype(GR.MARKERTYPE_DOT)
    scatter3(redv1, greenv1, bluev1, xlabel="red", ylabel="green", zlabel="blue"; color=colv1)

    the previous is not as good as using gnuplot
    some reasons:
    Plots and GR do not support x-, y-, z-lables as expected
    take too long
    ? grid ?
    ? user defined view angle ?

    =#
end # function plot_imhi_GrayRGB

end # module ImageHistoryTest


text:
using MomentClosure
using MomentClosure: define_μ, define_M
#using SymbolicUtils: expand
using Symbolics: expand
using Test

M = define_M(2, 4)
κ = MomentClosure.cumulants_to_central_moments(2,4)
@test isequal(κ[1,2], M[1,2])
expr = M[3,1] - 3*M[1,1]*M[2,0]
@test isequal(κ[3,1], expr)

μ = define_μ(2, 3)
κ = MomentClosure.cumulants_to_raw_moments(2,3)
@test isequal(κ[1,0], μ[1,0])
expr = μ[2,1] + 2*μ[0,1]*μ[1,0]^2 - μ[0,1]*μ[2,0] - 2*μ[1,0]*μ[1,1]
# if this one fails most likely simplification functions have broke
@test isequal(expand(κ[2,1]),expr)

raw_to_central = MomentClosure.raw_to_central_moments(2,3)
expr = 2*μ[1,0]^3 + μ[3,0] - 3*μ[1,0]*μ[2,0]
@test isequal(raw_to_central[3,0], expr)
expr = μ[2,1] + 2*μ[0,1]*μ[1,0]^2 - 2*μ[1,0]*μ[1,1] - μ[0,1]*μ[2,0]
@test isequal(raw_to_central[2,1], expr)

central_to_raw = MomentClosure.central_to_raw_moments(2,3)
expr = M[2,0] + μ[1,0]^2
@test isequal(central_to_raw[2,0], expr)
expr = 2*M[1,1]*μ[0,1] + M[0,2]*μ[1,0] + M[1,2] + μ[1,0]*μ[0,1]^2
@test isequal(expand(central_to_raw[1,2]), expr)


text:
function truncate(y, negtrunc, postrunc)

# Truncate input if its too big or too small
    if (y < negtrunc) # If y is less than the value negtrunc
       truncy = negtrunc;
    elseif (y > postrunc)
       truncy= postrunc; # If y is greater than the value postrunc
    else
       truncy = y;
    end

    return truncy

end

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
# using RCall
# @rimport rcommon

include("imports.jl")
using Distributed
rmprocs(filter(w -> w > 1, workers()))
addprocs(16)

@everywhere include("imports.jl")
@everywhere Random.seed!(0)

@everywhere eye(T::Type, n::Int) = Matrix{T}(LinearAlgebra.I, n, n)

@everywhere function gen_distribution(P)
  S = rand(InverseWishart(P, eye(Float64, P)))
  m = randn(P) * 10
  (m, S)
end

# True distribution has this covariance
# @everywhere P = 5  # NOTE: works great
@everywhere P = 20 # NOTE: works
# @everywhere P = 50 # NOTE: dimension too high?
@everywhere m, S = gen_distribution(P)
@everywhere mvn = MvNormal(m, S)

@everywhere function update(init::T, tau, propCov, iters) where T
  samps = T[]
  append!(samps, [copy(init)])

  logprob(z) = logpdf(mvn, z) / tau

  x = copy(init)
  for _ in 1:iters
    x = ParTemp.metropolisMV(x, logprob, propCov)[1]
    append!(samps, [copy(x)])
  end

  return samps
end

taus = 1.1 .^ (0:15)
# taus = .99 .^ (15:-1:0)
iters = 100
epochs = 100
buffsize = 1000

@time samps = let
  out = nothing
  state = [rand(P) for _ in taus]
  propCov = [1e-6 * eye(Float64, P) for _ in taus]
  states = [ParTemp.Buffer(Vector{Float64}, buffsize) for _ in taus]

  for i in 1:epochs
    print("\rIter $i / $epochs")
    out = pmap(update,
               state,
               taus,
               propCov,
               [iters for _ in taus])

    for t in 1:length(taus)
      ParTemp.append!(states[t], out[t])
    end

    propCov = [cov(states[t].x) + 1e-6 * eye(Float64, P)
               for t in 1:length(taus)]
    state = map(o -> o[end], out)

    # Swap chains
    for t in length(taus):-1:2
      j = t
      i = j - 1
      logprob(k) = logpdf(mvn, state[k])
      logprobswap = (1/taus[j] - 1/taus[i]) * (logprob(i) - logprob(j))
      if logprobswap > log(rand())
        println("Swapping ($j, $i)")
        tmp = copy(state[j])
        state[j] = copy(state[i])
        state[i] = tmp
      end
    end
  end
  println()

  states[1].x
end

begin
  tail = Int(length(samps) / 2)
  K = 5
  mean(samps[(end-tail):end])
  println("Est. Cov:")
  show(stdout, "text/plain", cov(samps[(end-tail):end])[1:K, 1:K])
  # show(stdout, "text/plain", cor(hcat(samps[(end-tail):end]...)')[1:K, 1:K])

  println()
  println()

  println("True Cov:")
  show(stdout, "text/plain", cov(mvn)[1:K, 1:K])
  # show(stdout, "text/plain", cor(mvn)[1:K, 1:K])
  # rcommon.plotPosts(Matrix(hcat(samps[(end-tail):end]...)')[:, 1:K]);

  println("mean(mvn), mean(samps)")
  show(stdout, "text/plain", [mean(mvn) mean(samps)])

  nothing
end;


text:
#!/usr/bin/env julia

using Luxor

using Test, Random

function drawbbox(apoly)
    gsave()
    setline(0.3)
    setdash("dotted")
    box(BoundingBox(apoly), :stroke)
    grestore()
end

function test1(p, x, y)
    prettypoly(p, :fill,
        vertexlabels = (n, l) -> (text(string(n, " of ", l))),
        # all these commands are executed for each vertex of the polygon
        () ->
        begin
            scale(0.5, 0.5)
            randomhue()
            opacity = Luxor.get_current_alpha()
            setopacity(rand(5:10)/10)
            ngon(0, 0, 15, rand(3:12), 0, :fill)
            setopacity(opacity)
            sethue("black")
        end,
    )
end

function test2(p)
    prettypoly(p,
        :fill,
        () ->
        begin
          ## all these commands are executed for each vertex of the polygon
            randomhue()
            scale(0.5, 0.5)
            rotate(pi/2)
            prettypoly(p, :fill,
                () ->
                begin
                # and all these commands are executed for each vertex of that polygon
                randomhue()
                scale(0.25, 0.25)
                rotate(pi/2)
                prettypoly(p, :fill)
                end
            )
        end
    )
end

function get_png_files()
    folder = dirname(dirname(pathof(Luxor))) * "/test/"
    imagelist = filter(f -> !startswith(f, ".") && endswith(f, "png"), readdir(folder))
    imagepathlist = map(fn -> folder * "/" * fn, imagelist)
    return imagepathlist
end

function test3(p)
    imagelist = get_png_files()
    shuffle!(imagelist)
    img = readpng(imagelist[1])
    w = img.width
    h = img.height
    prettypoly(p,
        :fill,
        () ->
        begin
            ## all these commands are executed for each vertex of the polygon
            placeimage(readpng(imagelist[rand(1:end)]), -w/2, -h/2)
        end
        )
end

function draw_lots_of_polys(pagewidth, pageheight)
    setopacity(0.75)
    pagetiles = Tiler(pagewidth, pageheight, 9, 9, margin=20)
    for (pos, n) in pagetiles
        if rand(Bool)
            p = randompointarray(rand(-140:-10), rand(-140:-10), rand(10:140), rand(10:140), rand(5:12))
        else
            p = ngon(0, 0, rand(50:100), rand(3:12), vertices=true)
        end
        gsave()
        translate(pos)
        setline(1)
        randomhue()
        if n % 3  == 0
            test1(p, pos.x, pos.y)
        elseif n % 3 == 1
            test2(p)
        elseif n % 3 == 2
            test3(p)
        end
        drawbbox(p)
        grestore()
    end
end

function test_pretty_poly(fname)
    width, height = 3000, 3000
    Drawing(width, height, fname)
    origin()
    background("ivory")
    # test empty prettypoly
    prettypoly(Point[], :stroke)
    draw_lots_of_polys(width, height)
    @test finish() == true
end

fname = "pretty-poly-test.pdf"
test_pretty_poly(fname)
println("...finished test: output in $(fname)")


text:
function freqresp(sys::LTISystem, w::Real)
    # Create imaginary freq vector s
    if iscontinuous(sys)
        s = im*w
    else
        s = cis(w*sys.Ts)
    end
    evalfr(sys, s)
end

"""
    sys_fr = freqresp(sys, w)

Evaluate the frequency response of a linear system

`w -> C*((iw*im -A)^-1)*B + D`

of system `sys` over the frequency vector `w`.
"""
@autovec () function freqresp(sys::LTISystem, w_vec::AbstractVector{<:Real})
    te = sys.timeevol
    ny,nu = noutputs(sys), ninputs(sys)
    [evalfr(sys[i,j], _freq(w, te))[] for w in w_vec, i in 1:ny, j in 1:nu]
end

_freq(w, ::Continuous) = complex(0, w)
_freq(w, te::Discrete) = cis(w*te.Ts)

@autovec () function freqresp(sys::AbstractStateSpace, w_vec::AbstractVector{W}) where W <: Real
    ny, nu = size(sys)
    T = promote_type(Complex{real(eltype(sys.A))}, Complex{W})
    if sys.nx == 0 # Only D-matrix
        return PermutedDimsArray(repeat(T.(sys.D), 1, 1, length(w_vec)), (3,1,2))
    end
    local F
    try
        F = hessenberg(sys.A)
    catch e
        # For matrix types that do not have a hessenberg implementation, we call the standard version of freqresp.
        e isa MethodError && return freqresp_nohess(sys, w_vec)
        rethrow()
    end
    Q = Matrix(F.Q)
    A = F.H
    C = sys.C*Q
    B = Q\sys.B
    D = sys.D

    te = sys.timeevol
    R = Array{T, 3}(undef, ny, nu, length(w_vec))
    Bc = similar(B, T) # for storage
    for i in eachindex(w_vec)
        Ri = @views R[:,:,i]
        copyto!(Ri,D) # start with the D-matrix
        isinf(w_vec[i]) && continue
        copyto!(Bc,B) # initialize storage to B
        w = -_freq(w_vec[i], te)
        ldiv!(A, Bc, shift = w) # B += (A - w*I)\B # solve (A-wI)X = B, storing result in B
        mul!(Ri, C, Bc, -1, 1) # use of 5-arg mul to subtract from D already in Ri. - rather than + since (A - w*I) instead of (w*I - A)
    end
    PermutedDimsArray(R, (3,1,2)) # PermutedDimsArray doesn't allocate to perform the permutation
end

"""
    freqresp_nohess(sys::AbstractStateSpace, w_vec::AbstractVector{<:Real})

Compute the frequency response of `sys` without forming a Hessenberg factorization.
This function is called automatically if the Hessenberg factorization fails.
"""
freqresp_nohess
@autovec () function freqresp_nohess(sys::AbstractStateSpace, w_vec::AbstractVector{W}) where W <: Real
    ny, nu = size(sys)
    nx = sys.nx
    T = promote_type(Complex{real(eltype(sys.A))}, Complex{W})
    if nx == 0 # Only D-matrix
        return PermutedDimsArray(repeat(T.(sys.D), 1, 1, length(w_vec)), (3,1,2))
    end
    A,B,C,D = ssdata(sys)
    te = sys.timeevol
    R = Array{T, 3}(undef, ny, nu, length(w_vec))
    Ac = (A+one(T)*I) # for storage
    Adiag = diagind(A)
    for i in eachindex(w_vec)
        Ri = @views R[:,:,i]
        copyto!(Ri,D) # start with the D-matrix
        isinf(w_vec[i]) && continue
        w = _freq(w_vec[i], te)
        @views copyto!(Ac[Adiag],A[Adiag]) # reset storage to A
        @views Ac[Adiag] .-= w # Ac = A - w*I
        Bc = Ac \ B # Bc = (A - w*I)\B # avoid inplace to handle sparse matrices etc.
        mul!(Ri, C, Bc, -1, 1) # use of 5-arg mul to subtract from D already in Ri. - rather than + since (A - w*I) instead of (w*I - A)
    end
    PermutedDimsArray(R, (3,1,2)) # PermutedDimsArray doesn't allocate to perform the permutation
end

function _evalfr_return_type(sys::AbstractStateSpace, s::Number)
    T0 = numeric_type(sys)
    temp_product = one(T0)*one(typeof(s))
    typeof(temp_product/temp_product)
end

"""
`evalfr(sys, x)` Evaluate the transfer function of the LTI system sys
at the complex number s=x (continuous-time) or z=x (discrete-time).

For many values of `x`, use `freqresp` instead.
"""
function evalfr(sys::AbstractStateSpace, s::Number)
    T = _evalfr_return_type(sys, s)
    try
        R = s*I - sys.A
        sys.D + sys.C*((R\sys.B))
    catch e
        @warn "Got exception $e, returning Inf" maxlog=1
        fill(convert(T, Inf), size(sys))
    end
end

function evalfr(G::TransferFunction{<:TimeEvolution,<:SisoTf}, s::Number)
    map(m -> evalfr(m,s), G.matrix)
end

"""
`F(s)`, `F(omega, true)`, `F(z, false)`

Notation for frequency response evaluation.
- F(s) evaluates the continuous-time transfer function F at s.
- F(omega,true) evaluates the discrete-time transfer function F at exp(im*Ts*omega)
- F(z,false) evaluates the discrete-time transfer function F at z
"""
function (sys::TransferFunction)(s)
    evalfr(sys,s)
end

function (sys::TransferFunction)(z_or_omega::Number, map_to_unit_circle::Bool)
    isdiscrete(sys) || throw(ArgumentError("It only makes no sense to call this function with discrete systems"))
    if map_to_unit_circle
        isreal(z_or_omega) ? evalfr(sys,exp(im*z_or_omega.*sys.Ts)) : error("To map to the unit circle, omega should be real")
    else
        evalfr(sys,z_or_omega)
    end
end

function (sys::TransferFunction)(z_or_omegas::AbstractVector, map_to_unit_circle::Bool)
    isdiscrete(sys) || throw(ArgumentError("It only makes no sense to call this function with discrete systems"))
    vals = sys.(z_or_omegas, map_to_unit_circle)# evalfr.(sys,exp.(evalpoints))
    # Reshape from vector of evalfr matrizes, to (in,out,freq) Array
    nu,ny = size(vals[1])
    [v[i,j]  for v in vals, i in 1:nu, j in 1:ny]
end

"""`mag, phase, w = bode(sys[, w])`

Compute the magnitude and phase parts of the frequency response of system `sys`
at frequencies `w`

`mag` and `phase` has size `(length(w), ny, nu)`"""
@autovec (1, 2) function bode(sys::LTISystem, w::AbstractVector)
    resp = freqresp(sys, w)
    return abs.(resp), rad2deg.(unwrap!(angle.(resp),1)), w
end
@autovec (1, 2) bode(sys::LTISystem) = bode(sys, _default_freq_vector(sys, Val{:bode}()))

"""`re, im, w = nyquist(sys[, w])`

Compute the real and imaginary parts of the frequency response of system `sys`
at frequencies `w`

`re` and `im` has size `(length(w), ny, nu)`"""
@autovec (1, 2) function nyquist(sys::LTISystem, w::AbstractVector)
    resp = freqresp(sys, w)
    return real(resp), imag(resp), w
end
@autovec (1, 2) nyquist(sys::LTISystem) = nyquist(sys, _default_freq_vector(sys, Val{:nyquist}()))

"""`sv, w = sigma(sys[, w])`

Compute the singular values `sv` of the frequency response of system `sys` at
frequencies `w`

`sv` has size `(length(w), max(ny, nu))`"""
@autovec (1) function sigma(sys::LTISystem, w::AbstractVector)
    resp = freqresp(sys, w)
    sv = dropdims(mapslices(svdvals, resp, dims=(2,3)),dims=3)
    return sv, w
end
@autovec (1) sigma(sys::LTISystem) = sigma(sys, _default_freq_vector(sys, Val{:sigma}()))

function _default_freq_vector(systems::Vector{<:LTISystem}, plot)
    min_pt_per_dec = 60
    min_pt_total = 200
    bounds = map(sys -> _bounds_and_features(sys, plot)[1], systems)
    w1 = minimum(minimum.(bounds))
    w2 = maximum(maximum.(bounds))

    nw = round(Int, max(min_pt_total, min_pt_per_dec*(w2 - w1)))
    return exp10.(range(w1, stop=w2, length=nw))
end
_default_freq_vector(sys::LTISystem, plot) = _default_freq_vector(
        [sys], plot)

function _bounds_and_features(sys::LTISystem, plot::Val)
    # Get zeros and poles for each channel
    if !isa(plot, Val{:sigma})
        zs, ps = zpkdata(sys)
        # Compose vector of all zs, ps, positive conjugates only.
        zpType = promote_type(eltype(eltype(zs)), eltype(eltype(ps)))
        zp = vcat(zpType[], zs..., ps...) # Emty vector to avoid type unstable vcat()
        zp = zp[imag(zp) .>= 0.0]
    else
         # For sigma plots, use the MIMO poles and zeros
         zp = [tzeros(sys); poles(sys)]
    end
    # Get the frequencies of the features, ignoring low frequency dynamics
    fzp = log10.(abs.(zp))
    fzp = fzp[fzp .> -4]
    fzp = sort!(fzp)
    # Determine the bounds on the frequency vector
    if !isempty(fzp)
        w1 = floor(fzp[1] - 1.2)
        w2 = ceil(fzp[end] + 1.2)
        # Expand the range for nyquist plots
        if plot isa Val{:nyquist}
            w1 -= 0.0
            w2 += 1.0
        end
    else
        w1 = 0.0
        w2 = 2.0
    end
    if isdiscrete(sys) # Do not draw above Nyquist freq for disc. systems
        w2 = min(w2, log10(π/sys.Ts))
    end
    return [w1, w2], zp
end


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

JLLWrappers.@generate_main_file_header("msolve")
JLLWrappers.@generate_main_file("msolve", UUID("6d01cc9a-e8f6-580e-8c54-544227e08205"))
end  # module msolve_jll


text:
#Version 2
function build_nabla2_x(LF)
  D2jl = LF.LFx.D2jl
  Nx = LF.Nx
  Ny = LF.Ny
  Nz = LF.Nz
  Npoints = Nx*Ny*Nz

  # initialize with
  rowval = zeros(Int64,Nx*Npoints)
  colptr = zeros(Int64,Npoints+1)
  colptr[1] = 1
  nzval = zeros(Float64,Nx*Npoints)

  idx = 1
  for colLoc = 1:Nx
    colGbl_start = (colLoc-1)*(Ny*Nz)+1
    colGbl_stop = colLoc*Ny*Nz
    #println("")
    #println("colGbl_start = ", colGbl_start)
    #println("colGbl_stop  = ", colGbl_stop)
    for colGbl = colGbl_start:colGbl_stop
      for rowLoc = 1:Nx
        rowGbl = (rowLoc-1)*Ny*Nz + colGbl - (colLoc-1)*Ny*Nz
        rowval[idx] = rowGbl
        nzval[idx] = D2jl[rowLoc,colLoc]
        idx = idx + 1
        #println(colLoc, " ", rowGbl)
        #append!( rowval, [rowGbl] )
        #append!( nzval, [D2jl[rowLoc,colLoc]] )
      end
      colptr[colGbl+1] = colptr[colGbl] + Nx
    end
  end # colLoc
  #@printf("Finished looping\n")

  #println(size(rowval[2:end]))
  #println(size(nzval[2:end]))
  #println(size(colptr))
  #println(colptr)
  return SparseMatrixCSC( Npoints, Npoints, colptr, rowval, nzval )
end

# Version 1
function build_nabla2_x(LF)
  D2jl = LF.LFx.D2jl
  Nx = LF.Nx
  Ny = LF.Ny
  Nz = LF.Nz
  Npoints = Nx*Ny*Nz

  # initialize with
  rowval = [0]  # start with dummy value 0
  colptr = zeros(Int64,Npoints+1)
  colptr[1] = 1
  nzval = [0.0]  # start with dummy value 0.0

  colLoc = 1
  for colGbl = 1:Ny*Nz
    for rowLoc = 1:Nx
      rowGbl = (rowLoc-1)*Ny*Nz + colGbl - 0*Ny*Nz
      #println(colLoc, " ", rowGbl)
      append!( rowval, [rowGbl] )
      append!( nzval, [D2jl[rowLoc,colLoc]] )
    end
    colptr[colGbl+1] = colptr[colGbl] + Nx
  end

  #println("")
  colLoc = 2
  colGblStart = (colLoc-1)*(Ny*Nz+1)
  colGblStop = colLoc*Ny*Nz
  for colGbl =Ny*Nz+1:2*Ny*Nz
    for rowLoc = 1:Nx
      rowGbl = (rowLoc-1)*Ny*Nz + colGbl - 1*Ny*Nz
      #println(colLoc, " ", rowGbl)
      append!( rowval, [rowGbl] )
      append!( nzval, [D2jl[rowLoc,colLoc]] )
    end
    colptr[colGbl+1] = colptr[colGbl] + Nx
  end

  #println("")
  colLoc = 3
  for colGbl = 2*Ny*Nz+1:3*Ny*Nz
    for rowLoc = 1:Nx
      rowGbl = (rowLoc-1)*Ny*Nz + colGbl - 2*Ny*Nz
      #println(colLoc, " ", rowGbl)
      append!( rowval, [rowGbl] )
      append!( nzval, [D2jl[rowLoc,colLoc]] )
    end
    colptr[colGbl+1] = colptr[colGbl] + Nx
  end

  println(size(rowval[2:end]))
  println(size(nzval[2:end]))
  println(size(colptr))
  println(colptr)
  return SparseMatrixCSC( Npoints, Npoints, colptr, rowval[2:end], nzval[2:end] )
end

function build_nabla2_x(LF)
  D2jl = LF.LFx.D2jl
  Nx = LF.Nx
  Ny = LF.Ny
  Nz = LF.Nz
  Npoints = Nx*Ny*Nz

  # initialize with
  rowval = [0]  # start with dummy value 0
  colptr = zeros(Int,Npoints)
  colptr[1] = 1
  nzval = [0.0]  # start with dummy value 0.0

  for jx=1:Nx
    j = (jx-1)*Ny*Nz + 1
    nnzc = 0
    for ix=1:Nx
      i = (ix-1)*Ny*Nz + 1
      for dd = 1:Ny*Nz
        nnzc = nnzc + 1
        append!( rowval, [i+dd-1])
        append!( nzval, [D2jl[ix,jx]] )
      end
      append!( colptr, [colptr[end] + nnzc])
    end
  end

  @printf("Converting to sparse matrix ..\n");
  println(size(colptr))
  println(rowval[2:end])
  #return sparse( rows[2:end], cols[2:end], vals[2:end] )
  #return SparseMatrixCSC( Npoints, Npoints, colptr, rowval[2:end], nzval[2:end] )
  exit()
end


text:
"""
Add a bias/intercept term to a vector on the fly without creating or copying data:

- `BiasVector(rand(10))` is roughly equivalent to `vcat(rand(10), 1.0)`
"""
immutable BiasVector{A<:AVecF} <: AVecF
  vec::A
end

Base.length(v::BiasVector) = length(v.vec) + 1
Base.size(v::BiasVector) = (length(v),)
Base.getindex(v::BiasVector, i::Integer) = i > length(v.vec) ? 1.0 : getindex(v.vec, i)
Base.setindex!(v::BiasVector, x::Real, i::Integer) = setindex!(v.vec, x, i)

"""
Adda bias/intercept term to a matrix on the fly without creating or copying data:

- `BiasMatrix(rand(10,5))` is roughly equivalent to `hcat(rand(10,5), ones(10))`
"""
immutable BiasMatrix{A<:AMatF} <: AMatF
  mat::A
end

Base.length(m::BiasMatrix) = length(m.mat) + nrows(m.mat)
Base.size(m::BiasMatrix) = (nrows(m.mat), ncols(m.mat)+1)
Base.getindex(m::BiasMatrix, i::Integer, j::Integer) = j > ncols(m.mat) ? 1.0 : getindex(m.mat, i, j)
Base.setindex!(m::BiasMatrix, x::Real, i::Integer, j::Integer) = setindex!(m.mat, x, i, j)

#-----------------------------------------------------------------------------------#
"""
Add second-order interaction terms on the fly without creating or copying data:

- `TwoWayInteractionVector(rand(p))` "adds" the `binomial(p, 2)` interaction terms
"""
immutable TwoWayInteractionVector{A<:AVecF} <: AVecF
    vec::A
end
Base.length(v::TwoWayInteractionVector) = length(v.vec) + binomial(length(v.vec), 2)
Base.size(v::TwoWayInteractionVector) = (length(v),)
function Base.getindex(v::TwoWayInteractionVector, i::Integer)
    first = 0
    second = i
    n = length(v.vec)
    if i <= n
        return v.vec[i]
    else
        while second > n
            second -= n
            first += 1
            n -= 1
        end
        return v.vec[first] * v.vec[first + second]
    end
end

"""
Add second-order interaction terms on the fly without creating or copying data:

- `TwoWayInteractionMatrix(rand(n, p))` "adds" the `binomial(p, 2)` interaction terms
to each row
"""
immutable TwoWayInteractionMatrix{A<:AMatF} <: AMatF
    mat::A
end
function Base.length(m::TwoWayInteractionMatrix)
    n, p = size(m.mat)
    n * (p + binomial(p, 2))
end
function Base.size(m::TwoWayInteractionMatrix)
    n, p = size(m. mat)
    (n, p + binomial(p, 2))
end
function Base.getindex(m::TwoWayInteractionMatrix, i::Integer, j::Integer)
    first = 0
    second = j
    n, p = size(m.mat)
    if j <= p
        return m.mat[i, j]
    else
        while second > p
            second -= p
            first += 1
            p -= 1
        end
        return m.mat[i, first] * m.mat[i, first + second]
    end
end