text:
"""
This program is to discover and plot the communities of a network by leiden algorithm
NOTE: the leiden algorithm is implemented by the python package leidenalg, so
        before doing community discovery, Conda, PyCall have to be installed as
        follows:

            import Pkg
            Pkg.add("Conda")
            Pkg.add("PyCall")
            Pkg.build("PyCall")

            using Conda

            Conda.pip_interop(true)
            #Conda.pip("install", "scipy")
            #Conda.pip("install", "numpy")
            Conda.pip("install", "leidenalg")

Contributors: Xiaoshan Nian
Date: August, 2020
Email: cen@njust.edu.cn
Github: https://github.com/yottoo/JuliaCommunity
"""

module JuliaCommunity

using LightGraphs, SimpleWeightedGraphs, GraphPlot
using PyCall
using DataFrames, CSV
using Statistics, StatsBase, Random
using Parameters, ProgressMeter
using Gadfly, Cairo, Compose

export JuliaCommunityInstance,
        discover_communities,
        optimise_resolution,
        save,
        plot_community,
        compute_cluster_coef

const leiden = pyimport("leidenalg")
const ig = pyimport("igraph")

@with_kw mutable struct PartitionMethod
    louvain::String = "louvain"
    CPM::String = "CPM"
    modularity::String = "modularity"
end

@with_kw mutable struct JuliaCommunityInstance
    task_series::String = ""

    methods::PartitionMethod = PartitionMethod()
    method::String = "CPM"    # CPM, Modularity (for leiden algorithm) and Louvain
    """
    NETWORK: the network data.
        from: id of the node from
        to: id of the node to
        weight: edge weight (if edge_weighted is false, this column could be ignored)
    NODES: the data for nodes, with ID required
        id: node id
    NODE_IMPORTANCES: the importances of vertices
        id: node id
        importance
    """
    network::DataFrame = DataFrame()
    graph = nothing
    igraph = nothing
    nodes::DataFrame = DataFrame()
    max_node_id::Int = 1
    node_importance_field::String = "importance"
    node_label_field::String = "id"

    edge_weighted::Bool = true
    node_weighted::Bool = false
    is_directed::Bool = true
    γ::Float16 = 0.001

    """
    COMMUNITIES: communities discovered.
        c: community id (start from 1)
        size: community size
        cluster_cof: clustering coefficent of the community
    MEMBERSHIPS: memberships indicating the nodes belonging to communities
        id: node id
        c: community id
    """
    communities::DataFrame = DataFrame()
    membership_vector::Array{Int} = []
    memberships::DataFrame = DataFrame()
    n_community::Int = 0

    modularity::Float64 = 0
    quality::Float64 = 0
end

function JuliaCommunityInstance(network::DataFrame;
        nodes::DataFrame=DataFrame(),
        node_label_field::String="id",
        node_importance_field="importance",
        edge_weighted::Bool=true,
        node_weighted::Bool=false,
        is_directed::Bool=true,
        method::String="CPM",
        to_summarise_graph::Bool=true,
        task_series::String="",)
    filter!(row -> row.from > 0 && row.to > 0, network)
    jc = JuliaCommunityInstance()
    jc.task_series = replace("_" * replace(task_series, "-" => "_"), "__" => "_")
    set_method(jc, method=method)
    print("\nLoading the network and nodes data...")
    jc.network = network
    check_network(jc)
    jc.nodes = nodes
    check_nodes(jc)
    jc.max_node_id = maximum(jc.nodes.id)
    jc.node_importance_field = node_importance_field
    jc.node_label_field = node_label_field
    jc.edge_weighted = edge_weighted
    jc.node_weighted = node_weighted
    jc.is_directed = is_directed

    if edge_weighted
        filter!(row -> row.weight > 0, jc.network)
    #= ============================================
    else if !is_directed
        jc.network = vcat(jc.network, jc.network)
        unique(jc.network)
    ============================================ =#
    end

    print("\n\tBuilding the graph from the network...")
    if jc.is_directed && jc.edge_weighted
        jc.graph = SimpleWeightedDiGraph(jc.network.from, jc.network.to, jc.network.weight)
        jc.igraph = ig.Graph(zip(jc.network.from .- 1, jc.network.to .- 1), directed=true, edge_attrs=Dict("weight" => jc.network.weight))
    elseif jc.is_directed
        jc.graph = SimpleDiGraph(jc.max_node_id)
        for i in 1:nrow(jc.network) add_edge!(jc.graph, jc.network.from[i], jc.network.to[i]) end
        jc.igraph = ig.Graph(zip(jc.network.from .- 1, jc.network.to .- 1), directed=true)
    elseif jc.edge_weighted
        jc.graph = SimpleWeightedGraph(jc.network.from, jc.network.to, jc.network.weight)
        jc.igraph = ig.Graph(zip(jc.network.from .- 1, jc.network.to .- 1), directed=false, edge_attrs=Dict("weight" => jc.network.weight))
    else
        jc.graph = SimpleGraph(jc.max_node_id)
        for i in 1:nrow(jc.network) add_edge!(jc.graph, jc.network.from[i], jc.network.to[i]) end
        jc.igraph = ig.Graph(zip(jc.network.from .- 1, jc.network.to .- 1), directed=false)
    end

    if to_summarise_graph summarise_graph(jc) end

    jc
end

function check_network(jc::JuliaCommunityInstance)
    columns = names(jc.network)
    if isnothing(findfirst(name -> name == "from", columns)) throw(error("The network data must have a 'from' column.")) end
    if isnothing(findfirst(name -> name == "to", columns)) throw(error("The network data must have a 'to' column.")) end
    if jc.edge_weighted && isnothing(findfirst(name -> name == "weight", columns)) throw(error("The network data must have a 'weight' column.")) end
end

function check_nodes(jc::JuliaCommunityInstance)
    columns = names(jc.nodes)
    if isnothing(findfirst(name -> name == "id", columns)) throw(error("The nodes data must have a 'id' column.")) end
end

function set_method(jc::JuliaCommunityInstance; method::String="CPM")
    jc.method = method
    if (jc.method != jc.methods.louvain && jc.method != jc.methods.CPM && jc.method != jc.methods.modularity)
        throw(error("The partition method has to be louvain, CMP or modularity."))
    end
end

function summarise_graph(jc::JuliaCommunityInstance)
    print("\n\tSummary of the graph built on the network: ")
    print("\n\t\tNodes: $(nv(jc.graph))[$(nrow(jc.nodes))]")
    print("\n\t\tEdges: $(ne(jc.graph))")
    print("\n\t\tDensity: $(density(jc.graph))")
    print("\n\t\tDiameter: $(jc.igraph.diameter())")
    print("\n\t\tAverage path length: $(jc.igraph.average_path_length(directed=jc.is_directed))")
    print("\n\t\tClustering coefficient: $(global_clustering_coefficient(jc.graph))")
    print("\n\t\tModularity: $(modularity(jc.graph, fill(1, nv(jc.graph)), γ=1.0))\n")
end

function plot_network(jc::JuliaCommunityInstance; fig_path::String="fig", mute::Bool=false, node_size_smoother::Float64=0.5, edge_width_smoother::Float64=0.5, line_type::String="straight", arrow_length_frac::Float64=0.015, arrow_angle::Float64=π / 8)
    if !mute print("\n\nPlotting the network graph......\n\t")    end

    g = jc.graph

    if nv(g) > 5000 throw(error("It makes no sense to plot a huge network with more than 5000 vertices.")) end

    labels = jc.nodes[:, jc.node_label_field]
    nodesize = jc.node_weighted ? jc.nodes[:, jc.node_importance_field].^node_size_smoother : fill(1, nv(g))
    # labels[findall(size -> size < ceil(ne(jc.graph) / 100), nodesize)] .= ""
    # nodesize = log.(nodesize .+ 1)
    # edge_weights = jc.network[:weight].^0.72
    edge_weight = jc.edge_weighted ? jc.network.weight.^edge_width_smoother : fill(1, ne(g))
    # g = Graph(adjacency_matrix(g))
    node_colorss = ["orange", "purple", "turquoise", "green", "red", "blue", "violet", "olive", "tan", "magenta", "cyan", "pink", "gold"]
    edge_colors = ["navajowhite3", "coral2", "orange3", "yellow3", "yellowgreen", "turquoise", "lightskyblue", "mediumpurple1", "hotpink2", "tan3", "grey64"]
    shuffle!(node_colorss)
    # node_colors = colors[1 .+ Int.((length(colors) - 1) .* ceil.((nodesize .-  minimum(nodesize)) ./ (maximum(nodesize) - minimum(nodesize))))]
    # edge_colors = colors[1 .+ Int.((length(colors) - 1) .* ceil.((weights .-  minimum(weights)) ./ (maximum(weights) - minimum(weights))))]

    #= ========================================================================
    Further partition the target community and render each sub partition with a same color
    ======================================================================== =#
    igraph = jc.igraph

    categories = leiden.find_partition(igraph, leiden.ModularityVertexPartition)
    # categories = leiden.find_partition(igraph, leiden.CPMVertexPartition, resolution_parameter= 1 / nv(g))
    node_colors = node_colorss[(categories.membership .+ 1) .% length(node_colorss) .+ 1]
    # layout = (args...) -> spring_layout(args...; C = 12)   # where C influences the desired distance between nodes.

    run_label = "$(jc.method)-$(jc.γ)$(jc.task_series)"

    plot = nothing

    if jc.is_directed
        plot = gplot(g, nodesize=nodesize, nodelabel=labels,
                    nodelabeldist=0.2, nodelabelangleoffset=π / 4,
                    nodelabelsize=nodesize, edgelinewidth=edge_weight,
                    nodefillc=node_colors,
                    nodelabelc=node_colors,
                    linetype=line_type,
                    edgestrokec=edge_colors[rand(1:length(edge_colors), 1)[1]], arrowlengthfrac=arrow_length_frac, arrowangleoffset=arrow_angle
                    );
                    # or linetype = "curve" or "straight"
                    # edgestrokec = edge_colors,
                    # edgestrokec=colors[rand(1:length(colors), 1)[1]]
    else
        plot = gplot(g, nodesize=nodesize, nodelabel=labels,
                    nodelabeldist=0.2, nodelabelangleoffset=π / 4,
                    nodelabelsize=nodesize, edgelinewidth=edge_weight,
                    nodefillc=node_colors,
                    nodelabelc=node_colors,
                    linetype=line_type,
                    edgestrokec=edge_colors[rand(1:length(edge_colors), 1)[1]]);
    end
    if !ispath(fig_path) mkpath(fig_path) end
    filename = "$fig_path/network-graph-$run_label.svg"
    draw(SVG(filename, 20cm, 16cm), plot);
    open_file(filename)
end

"""
caculate_centralities:  Caculate the centralities of the network.
    Note: for weighted networks, eigenvector centrality measure may not lead to convergency.
"""
function caculate_centralities(jc::JuliaCommunityInstance; to_save::Bool=true)
    print("\n\tCaculating the centralities of the target network......\n")

    centralities = DataFrame(id=collect(vertices(jc.graph)))
    insertcols!(centralities, "btw_cent"  => betweenness_centrality(jc.graph),
        "close_cent"  => closeness_centrality(jc.graph),
        "in_cent"  => indegree_centrality(jc.graph),
        "out_cent"  => outdegree_centrality(jc.graph),
        "eigen_cent"  => eigenvector_centrality(jc.graph),
        "katz_cent"  => katz_centrality(jc.graph, 0.3),
        "stress_cent"  => stress_centrality(jc.graph),
        "radia_cent"  => radiality_centrality(jc.graph))

    insertcols!(centralities, "btw_cent_rank" => competerank(centralities[:, "btw_cent"], rev=true),
        "close_cent_rank" => competerank(centralities[:, "close_cent"], rev=true),
        "in_cent_rank" => competerank(centralities[:, "in_cent"], rev=true),
        "out_cent_rank" => competerank(centralities[:, "out_cent"], rev=true),
        "eigen_cent_rank" => competerank(centralities[:, "eigen_cent"], rev=true),
        "katz_cent"  * "_rank" => competerank(centralities[:, "katz_cent"], rev=true),
        "stress_cent_rank" => competerank(centralities[:, "stress_cent"], rev=true),
        "radia_cent"  * "_rank" => competerank(centralities[:, "radia_cent"], rev=true))

    centralities = leftjoin(jc.nodes, centralities, on=:id)

    if to_save
        save_path = "data/centralities$(jc.task_series).csv"
        centralities_merged =  ispath(save_path) ? hcat(CSV.File(save_path) |> DataFrame, centralities) : centralities
        CSV.write(save_path, centralities_merged)
    end

    centralities
end

"""
_louvain:    the louvain algorithm implemented by leiden package
"""
function _louvain(g)
    optimiser = leiden.Optimiser()
    partitions = leiden.ModularityVertexPartition(g)
    partitions_agg = partitions.aggregate_partition()
    while optimiser.move_nodes(partitions) > 0
        partitions.from_coarse_partition(partitions_agg)
        partitions_agg = partitions_agg.aggregate_partition()
    end
    partitions
end

"""
discover_communities:    discover the communities by leiden algorithm (both CMP and modularity method)
    and louvain algorithm.
"""
function discover_communities(jc::JuliaCommunityInstance; mute::Bool=false)
    if !mute print("\nDiscovering the communities for the built network......") end

    partitions = nothing
    if jc.method == jc.methods.CPM
        partitions = leiden.find_partition(jc.igraph, leiden.CPMVertexPartition, resolution_parameter=jc.γ)
    elseif jc.method == jc.methods.modularity
        partitions = leiden.find_partition(jc.igraph, leiden.ModularityVertexPartition)
    elseif jc.method == jc.methods.louvain
        partitions = _louvain(jc.igraph)
    end

    if isnothing(partitions) return end

    jc.n_community = length(partitions)
    if !mute println("\t\t$(jc.n_community) communities have been discovered.") end

    # partitions = leiden.find_partition(g, leiden.ModularityVertexPartition, resolution_parameter=0.2)
    # The following code could not be done to Leiden community members
    # for i in 1:length(partitions) partitions[i] .+= 1 end
    # partitions.membership .+= 1
    jc.membership_vector = partitions.membership .+ 1
    # println(partitions.membership)

    jc.modularity = modularity(jc.graph, jc.membership_vector, γ=1.0)
    jc.quality =  jc.method == jc.methods.CPM ? partitions.quality() : jc.modularity

    jc.memberships = DataFrame(id=1:maximum(vcat(jc.network.from, jc.network.to)), c=jc.membership_vector)
    # jc.memberships = DataFrame(id=1:nv(g), c=jc.membership_vector)
    jc.communities = DataFrame(c=1:jc.n_community, size=length.(partitions))
end

function optimise_resolution(jc::JuliaCommunityInstance; γ_from::Float64=0.0001, γ_end::Float64=0.01, γ_step::Float64=0.0001)
    println("\n")
    qualities = DataFrame(resolution=[], n_community=[], modularity=[], quality=[])
    jc_copy = deepcopy(jc)
    progress = Progress(length(γ_from:γ_step:γ_end), desc="Finding the best resolution γ for the Leiden-based community discovery algorithm: ")
    for γ  in γ_from:γ_step:γ_end
        jc_copy.γ = γ
        discover_communities(jc_copy, mute = true)
        push!(qualities, (γ, jc_copy.n_community, jc_copy.modularity, jc_copy.quality))
        next!(progress)
        #print("\n\t\tResolution: $γ: $(jc_copy.n_community) commxunities discovered; Modularity: $(jc_copy.modularity); CPM Quality: $(jc_copy.quality).")
    end
    CSV.write("data/community_discover_optimisation-$(jc.method)$(jc.task_series).csv", qualities)

    p_modularity = plot(
        layer(qualities, x=:resolution, y=:modularity, Geom.line, Geom.point, Theme(default_color="blue")),
        Guide.xticks(ticks=γ_from:γ_step * 2:γ_end),
        Guide.xlabel("resolution γ"),
        Guide.ylabel("modularity"),
        Theme(major_label_font_size=10pt),
        Scale.y_continuous(format=:plain)
    )
    fig_modularity = "fig/community_discover_optimisation-modularities-$(jc.method)$(jc.task_series).svg"
    draw(SVG(fig_modularity, 24cm, 16cm), p_modularity);
    open_file(fig_modularity)

    p_quality = plot(
        layer(qualities, x=:resolution, y=:quality, Geom.line, Geom.point, Theme(default_color="blue")),
        Guide.xticks(ticks=γ_from:γ_step * 2:γ_end),
        Guide.xlabel("resolution γ"),
        Guide.ylabel("quality"),
        Theme(major_label_font_size=10pt),
        Scale.y_continuous(format=:plain)
    )
    fig_quality = "fig/community_discover_optimisation-qualities-$(jc.method)$(jc.task_series).svg"
    draw(SVG(fig_quality, 24cm, 16cm), p_quality);
    open_file(fig_quality)

    p_n_communities = plot(
        layer(qualities, x=:resolution, y=:n_community, Geom.line, Geom.point, Theme(default_color="blue")),
        Guide.xticks(ticks=γ_from:γ_step * 2:γ_end),
        Guide.xlabel("resolution γ"),
        Guide.ylabel("number of communities"),
        Theme(major_label_font_size=10pt),
        Scale.y_continuous(format=:plain)
    )
    fig_n_communities = "fig/community_discover_optimisation-n-communities-$(jc.method)$(jc.task_series).svg"
    draw(SVG(fig_n_communities, 24cm, 16cm), p_n_communities);
    open_file(fig_n_communities)
end

function save(jc::JuliaCommunityInstance; file::String="data/communities")
    print("\n\nSaving the communities discovery outcomes......\n")
    run_label = "$(jc.method)-$(jc.γ)$(jc.task_series)"
    CSV.write("$file-$run_label.csv", jc.communities)
    CSV.write("$file-memberships-$run_label.csv", jc.memberships)
end

function build_community_graph(jc::JuliaCommunityInstance, c::Int; to_save_data::Bool=true, mute::Bool=false)
    if !mute print("\tBuilding the graph for community $c......\n\t") end
    community = select(filter(row -> row.c == c, jc.memberships), Not(:c))
    if nrow(community) < 2 return (graph = nothing, edges = 1) end
    sort!(community, :id)
    community_size = nrow(community)

    nodes = DataFrame(_label=community.id, new_id=1:nrow(community))
    if jc.node_weighted
        nodes = select(innerjoin(nodes, select(jc.nodes, "id", jc.node_importance_field), on=:_label => :id), "_label", "new_id", jc.node_importance_field => "importance")
    end

    if jc.edge_weighted
        sub_network = rename(select(innerjoin(nodes, jc.network, on=:_label => :from), [:new_id, :to, :weight]), :new_id => "from")
        sub_network = rename(select(innerjoin(nodes, sub_network, on=:_label => :to), [:new_id, :from, :weight]), :new_id => "to")
    else
        sub_network = rename(select(innerjoin(nodes, jc.network, on=:_label => :from), [:new_id, :to]), :new_id => "from")
        sub_network = rename(select(innerjoin(nodes, sub_network, on=:_label => :to), [:new_id, :from]), :new_id => "to")
    end

    # There may be communities in which all nodes are not connected to others
    if nrow(sub_network) <= 0 return (graph = nothing, edges = 0) end

    if jc.node_label_field != "id"
        nodes = innerjoin(select(jc.nodes, :id, jc.node_label_field), nodes, on=:id => :_label)
        rename!(nodes, jc.node_label_field => "label")
    else
        rename!(nodes, :_label => "label")
    end

    graph = nothing
    if jc.is_directed && jc.edge_weighted
        graph = SimpleWeightedDiGraph(sub_network.from, sub_network.to, sub_network.weight)
    elseif jc.is_directed
        graph = SimpleDiGraph(community_size)
        for i in 1:nrow(sub_network) add_edge!(graph, sub_network.from[i], sub_network.to[i]) end
    elseif jc.edge_weighted
        graph = SimpleWeightedGraph(sub_network.from, sub_network.to, sub_network.weight)
    else
        graph = SimpleGraph(community_size)
        for i in 1:nrow(sub_network) add_edge!(graph, sub_network.from[i], sub_network.to[i]) end
    end

    if to_save_data
        CSV.write("data/graph_network_$(jc.method)-$(jc.γ)$(jc.task_series)-$c.csv", sub_network)
        if jc.node_weighted CSV.write("data/graph_nodes_$(jc.method)-$(jc.γ)$(jc.task_series)-$c.csv", nodes) end
    end

    if jc.node_weighted
        (graph = graph, node_labels = nodes.label, network = sub_network, community_size = community_size, node_weights = nodes.importance)
    else
        (graph = graph, node_labels = nodes.label, network = sub_network, community_size = community_size)
    end
end

function plot_community(jc::JuliaCommunityInstance, c::Int; fig_path::String="fig", to_save_data::Bool=true, mute::Bool=false, node_size_smoother::Float64=0.5, edge_width_smoother::Float64=0.5, line_type="straight", arrow_length_frac::Float64=0.012, arrow_angle::Float64=π / 18)
    if !mute print("\n\nPlotting the community $c......\n\t")    end

    community_graph = build_community_graph(jc, c, to_save_data=to_save_data, mute=mute)
    if isnothing(community_graph) return end
    g = community_graph.graph
    if isnothing(g) return  end
    labels = community_graph.node_labels
    nodesize = jc.node_weighted ? community_graph.node_weights.^node_size_smoother : fill(1, nv(g))
    # labels[findall(size -> size < ceil(community.edges / 100), nodesize)] .= ""
    # nodesize = log.(nodesize .+ 1)
    # edge_weights = log.(community.network[:co_f] .+ 1)
    # edge_weights = community_graph.network[:weight].^0.72
    edge_weight = jc.edge_weighted ? community_graph.network.weight.^edge_width_smoother : fill(1, ne(g))
    # g = Graph(adjacency_matrix(g))
    node_colorss = ["orange", "purple", "turquoise", "green", "red", "blue", "violet", "olive", "tan", "magenta", "cyan", "pink", "gold"]
    edge_colors = ["navajowhite3", "coral2", "orange3", "yellow3", "yellowgreen", "turquoise", "lightskyblue", "mediumpurple1", "hotpink2", "tan3", "grey64"]
    shuffle!(node_colorss)
    # node_colors = colors[1 .+ Int.((length(colors) - 1) .* ceil.((nodesize .-  minimum(nodesize)) ./ (maximum(nodesize) - minimum(nodesize))))]
    # edge_colors = colors[1 .+ Int.((length(colors) - 1) .* ceil.((weights .-  minimum(weights)) ./ (maximum(weights) - minimum(weights))))]

    #= ========================================================================
    Further partition the target community and render each sub partition with a same color
    ======================================================================== =#
    network = community_graph.network
    igraph = nothing
    if jc.is_directed && jc.edge_weighted
        igraph = ig.Graph(zip(network.from .- 1, network.to .- 1), directed=true, edge_attrs=Dict("weight" => network.weight))
    elseif jc.is_directed
        igraph = ig.Graph(zip(network.from .- 1, network.to .- 1), directed=true)
    elseif jc.edge_weighted
        igraph = ig.Graph(zip(network.from .- 1, network.to .- 1), directed=false, edge_attrs=Dict("weight" => network.weight))
    else
        igraph = ig.Graph(zip(network.from .- 1, network.to .- 1), directed=false)
    end

    categories = leiden.find_partition(igraph, leiden.ModularityVertexPartition)
    # categories = leiden.find_partition(igraph, leiden.CPMVertexPartition, resolution_parameter= 1 / nv(g))
    node_colors = node_colorss[(categories.membership .+ 1) .% length(node_colorss) .+ 1]
    # layout = (args...) -> spring_layout(args...; C = 12)   # where C influences the desired distance between nodes.

    run_label = "$(jc.method)-$(jc.γ)$(jc.task_series)"

    plot = nothing

    if jc.is_directed
        plot = gplot(g, nodesize=nodesize, nodelabel=labels,
                    nodelabeldist=0.2, nodelabelangleoffset=π / 4,
                    nodelabelsize=nodesize, edgelinewidth=edge_weight,
                    nodefillc=node_colors,
                    nodelabelc=node_colors,
                    linetype=line_type,
                    edgestrokec=edge_colors[rand(1:length(edge_colors))], arrowlengthfrac=arrow_length_frac, arrowangleoffset=arrow_angle
                    );
                    # or linetype = "curve" or "straight"
                    # edgestrokec = edge_colors,
                    # edgestrokec=colors[rand(1:length(colors), 1)[1]]
    else
        plot = gplot(g, nodesize=nodesize, nodelabel=labels,
                    nodelabeldist=0.2, nodelabelangleoffset=π / 4,
                    nodelabelsize=nodesize, edgelinewidth=edge_weight,
                    nodefillc=node_colors,
                    nodelabelc=node_colors,
                    linetype=line_type,
                    edgestrokec=edge_colors[rand(1:length(edge_colors))]);
    end
    save_path = "$fig_path/community-$run_label"
    if !ispath(save_path) mkpath(save_path) end
    draw(SVG("$save_path/$c.svg", 20cm, 16cm), plot);
end

function compute_cluster_coef(jc::JuliaCommunityInstance)
    print("\nComputing the clustering coefficients for the partitioned communities......\n\t")

    # If multiple threads are used to speed the caculation, then we have to bind the thread id with the corresponding index of the array.
    cluster_cofs = zeros(jc.n_community)

    Threads.@threads for c in 1:jc.n_community
        community = build_community_graph(jc, c, to_save_data=false, mute=true)
        g = community.graph
        if isnothing(g)
            cluster_cofs[c] = community.edges
        else
            cluster_cofs[c] = global_clustering_coefficient(g)
        end
    end

    insertcols!(jc.communities, :cluster_cof => cluster_cofs)
    avg_cluster_cof = mean(cluster_cofs)

    print("\tThe total clustering coefficient of the network  all the communities is $(global_clustering_coefficient(jc.graph)); The average clustering coefficient for all the communities is $avg_cluster_cof.")
end

function open_file(filename)
    if Sys.isapple()
        run(`open $(filename)`)
    elseif Sys.islinux() || Sys.isbsd()
        run(`xdg-open $(filename)`)
    elseif Sys.iswindows()
        run(`$(ENV["COMSPEC"]) /c start $(filename)`)
    else
        @warn "Showing plots is not supported on OS $(string(Sys.KERNEL))"
    end
end

end

max_stars_repo_name: yottoo/JuliaCommunity

max_issues_repo_name: yottoo/JuliaCommunity

max_forks_repo_name: yottoo/JuliaCommunity


text:
#=
type and methods for handling signals generating by consuming data exhaust from indicators
=#

struct Signal
    switch::Expr
    fun::Function
    function Signal(switch::Expr)
        @assert typeof(eval(switch.args[1])) <: Function
        f::Function = eval(switch.args[1])
        a::Symbol = switch.args[2]
        b::Symbol = switch.args[3]
        #pair = eval(switch.args[1])::Function => Tuple{Symbol,Symbol}(switch.args[2]::Symbol, switch.args[3]::Symbol)
        #fun(x::TS)::BitVector = pair.first(x[pair.second[1]].values[:]::Vector, x[pair.second[2]].values[:]::Vector)
        function fun(x::TS)::BitVector
            #fld1::Symbol = switch.args[2]
            #fld2::Symbol = switch.args[3]
            vec1::Vector = x[a].values[:]
            vec2::Vector = x[b].values[:]
            comp::Function = eval(switch.args[1])
            out::BitVector = comp(vec1, vec2)
            return out
        end
        return new(switch, fun)
    end
end

function prep_signal(signal::Signal, indicator_data::TS)::Expr
    local switch = copy(signal.switch)
    for i in 2:length(switch.args)
        switch.args[i] = indicator_data[switch.args[i]]
    end
    return switch
end

macro signal(logic::Expr)
    return Signal(logic)
end

↑(x, y) = Indicators.crossover(x, y)
↓(x, y) = Indicators.crossunder(x, y)

ext: jl


text:
"""
    update_bgc_params(tmp::Ditc)
Update parameter values based on a `Dict` provided by user

Keyword Arguments
=================
- `tmp` is a `Dict` containing the parameters needed to be upadated
"""
function update_bgc_params(tmp::Dict)
    parameters = bgc_params_default()
    tmp_keys = collect(keys(tmp))
    pkeys = collect(keys(parameters))
    for key in tmp_keys
        if length(findall(x->x==key, pkeys))==0
            throw(ArgumentError("PARAM: parameter not found $key"))
        else
            parameters[key] = tmp[key]
        end
    end
    return parameters
end

"""
    update_phyt_params(tmp::Dict, mode::AbstractMode)
Update parameter values based on a `Dict` provided by user

Keyword Arguments
=================
- `tmp` is a `Dict` containing the parameters needed to be upadated
- `mode` is the mode of phytoplankton physiology resolved in the model
"""
function update_phyt_params(tmp::Dict, mode::AbstractMode)
    parameters = phyt_params_default(1,mode)
    tmp_keys = collect(keys(tmp))
    pkeys = collect(keys(parameters))
    for key in tmp_keys
        if length(findall(x->x==key, pkeys))==0
            throw(ArgumentError("PARAM: parameter not found $key"))
        else
            parameters[key] = tmp[key]
        end
    end
    return parameters
end


text:
using ArgParse

function parse_cmdline()
    settings = ArgParseSettings()
    @add_arg_table settings begin
        "--gc"
            arg_type = Bool
            default = true
            help = "collect garbage between iterations"
        "--ci"
            arg_type = FloatingPoint
            default = 0.95
            help = "bootstrap confidence interval"
        "--list"
            action = :store_true
            help = "print only a list of benchmark names"
        "--report"
            arg_type = String
            default = nothing
            help = "file to write report"
        "--quiet"
            action = :store_true
            help = "print less output"
        "--resamples"
            arg_type = Integer
            default = 100 * 1000
            help = "number of botstrap resamples to perform"
        "--results"
            arg_type = String
            default = nothing
            help = "file to write raw results to"
        "--samples"
            arg_type = Integer
            default = 100
            help = "number of samples to collect"
        "--summary"
            arg_type = String
            default = nothing
            help = "produce a summary csv file of all results"
        "--no-measurements"
            action = :store_true
            help = "don't do any measurements"
        "--junit"
            arg_type = String
            default = nothing
            help = "produce a junit report file of all results"
        end
    parse_args(settings)
end

function main()
    parsed_args = parse_cmdline()
    cfg = config(parsed_args)
end


text:
sz = (4,5)
#filename = "/tmp/test-9.nc"
# The mode "c" stands for creating a new file (clobber)

for format in [:netcdf4, :netcdf4_classic, :netcdf3_classic, :netcdf3_64bit_offset]
    filenamefmt = tempname()
    NCDatasets.NCDataset(filenamefmt,"c"; format = format) do ds

        # define the dimension "lon" and "lat"
        NCDatasets.defDim(ds,"lon",sz[1])
        NCDatasets.defDim(ds,"lat",sz[2])

        # variables
        vfmt = NCDatasets.defVar(ds,"var-Float32",Float32,("lon","lat"))

        # write array
        vfmt[:,:] = fill(Float32(123),size(vfmt))

        # check content
        @test all(vfmt[:,:][:] .== 123)
    end
    rm(filenamefmt)
end

@test_throws NCDatasets.NetCDFError NCDatasets.NCDataset(
    tempname(),"c";
    format = :netcdf3000_perfect_format)


text:
@testset "Arithmetics" begin
    f0,f1,f2,f3,f4 = 0.125,-0.5,1.75,2.5,0.0

    # basic arithmetic operations
    @testset for T in (Posit8, Posit16, Posit32, Posit24_2)
        @test T(f0 + f1 * f2 - f3) == T(f0) + T(f1) * T(f2) - T(f3)
        @test T(f0 / f1) == T(f0) / T(f1)
        @test T(sqrt(f0) * f1 / f3) == sqrt(T(f0)) * T(f1) / T(f3)
        @test T(f1) ^ 2 == T(f1 ^ 2)
        @test T(f3) ^ 5 == T(f3 ^ 5)
    end

    # triggering complex infinity and conversion to inf
    @test Inf == Float64(Posit24_2(1.0)/Posit24_2(0.0))
end


text:
using ZigZagBoomerang
using LinearAlgebra
using Random
using SparseArrays
using Distributions
using Statistics
using ZigZagBoomerang: idot, idot_moving!
const ZZB = ZigZagBoomerang
using ReverseDiff #ReverseDiffSparse?
using Test
using FileIO
using GLM

println("Sparse logistic regression")

# Design matrix
Random.seed!(2)
sparsity(A, d = 3) = round(nnz(A)/length(A), digits=d)
include("./../sparsedesign.jl")
# create mock design matrix with 2 categorical explanatory variables
# and their interaction effects and 2 continuous explanatory variable
A = sparse_design((20,20), 2, 4*20)
n, p = size(A)
@show n, p
Γ = A'*A
@show sparsity(A), sparsity(Γ)

const m = 1 # m Bernoulli samples for each set of covariates
println("Av. number of regressors per column ", mean(sum(A .!= 0, dims=1)), ", row ", mean(sum(A .!= 0, dims=2)))

At = SparseMatrixCSC(A')
γ0 = 0.01
# Data from the model
wc = 0.5
xtrue = 2.0.*(rand(p) .> wc)  # either 2 or 0
sigmoid(x) = inv(one(x) + exp(-x))
lsigmoid(x) = -log(one(x) + exp(-x))
y_ = Matrix(rand(n, m) .< sigmoid.(A*xtrue))
y = sum(y_, dims=2)[:]
ny = m .- y

### prior and potential
# prior precision
# dot(x,x) might be more efficient if we know that many of those elements are 0.
# A*x might be more efficient if we know that many of those elements are 0.
ϕ(x, A, y::Vector) = γ0*dot(x,x)/2 - sum(y .* lsigmoid.(A*x) + (m .- y) .* lsigmoid.(-A*x))

# Sparse gradient
# helper functions for sparse gradient
function fdot(A::SparseMatrixCSC, At::SparseMatrixCSC, f, j, x, y)
   rows = rowvals(A)
   vals = nonzeros(A)
   s = zero(eltype(A))
   @inbounds for i in nzrange(A, j)
       s += vals[i]*y[rows[i]]*f(idot(At, rows[i], x))
   end
   s
end

sigmoidn(x) = sigmoid(-x)
nsigmoid(x) = -sigmoid(x)
# helper function for sparse gradient estimate through subsampling
function fdotr(A::SparseMatrixCSC, At::SparseMatrixCSC, j, x, μ, y, ny, k)
   rows = rowvals(A)
   vals = nonzeros(A)
   s = zero(eltype(A))
   r = nzrange(A, j)
   sampler = Random.SamplerRangeNDL(r)
   l = length(r)
   @inbounds for _ in 1:k
       i = rand(sampler)
       u = idot(At, rows[i], x)
       s += l/k*vals[i]*y[rows[i]]*sigmoidn(u)
       s += l/k*vals[i]*ny[rows[i]]*nsigmoid(u)
       u0 = idot(At, rows[i], μ)
       s -= l/k*vals[i]*y[rows[i]]*sigmoidn(u0)
       s -= l/k*vals[i]*ny[rows[i]]*nsigmoid(u0)   end
   s
end

# the same, but advancing only the coordinates necessary with `idot_moving!`
function fdot_moving(A::SparseMatrixCSC, At::SparseMatrixCSC, j, t, x, θ, t′, F, μ, y, ny, k)
   rows = rowvals(A)
   vals = nonzeros(A)
   s = zero(eltype(A))
   r = nzrange(A, j)
   sampler = Random.SamplerRangeNDL(r)
   l = length(r)
   @inbounds for _ in 1:k
       i = rand(sampler)
       u = idot_moving!(At, rows[i], t, x, θ, t′, F)
       s += l/k*vals[i]*y[rows[i]]*sigmoidn(u)
       s += l/k*vals[i]*ny[rows[i]]*nsigmoid(u)
       u0 = idot(At, rows[i], μ)
       s -= l/k*vals[i]*y[rows[i]]*sigmoidn(u0)
       s -= l/k*vals[i]*ny[rows[i]]*nsigmoid(u0)
   end
   s
end

# Gradient (found using http://www.matrixcalculus.org/ S. Laue, M. Mitterreiter, and J. Giesen. A Simple and Efficient Tensor Calculus, AAAI 2020.)
∇ϕ(x, A, At, y) = γ0*x - A'*(y .* sigmoidn.(A*x)) - A'*((m .- y).*nsigmoid.(A*x))

# Element i of the gradient exploiting sparsity
∇ϕ(x, i, A, At, y, ny = m .- y) = γ0*x[i] - fdot(A, At, sigmoidn, i, x, y) - fdot(A, At, nsigmoid, i, x, ny)

#res = glm(A, y/m, Binomial(m), LogitLink())
#x0 = coef(res)
#norm(xtrue - x0)

# Some Newton steps towards the mode as starting point
println("Newton steps")
x0 = 0.1rand(p)
@time for i in 1:30
    global x0
    H = hcat((sparse(ReverseDiff.gradient(x -> ∇ϕ(x, i, A, At, y, ny), x0)) for i in 1:p)...)
    x0 = x0 - (Symmetric(0.00I + H))\[∇ϕ(x0, i, A, At, y, ny) for i in 1:p]
end
norm(xtrue - x0)
fullgradx0 = ∇ϕ(x0, A, At, y)
# Element i of the gradient exploiting sparsity and random subsampling
∇ϕr(x, i, A, At, μ, y, ny = m .- y, k = 5) = γ0*x[i] - fdotr(A, At, i, x, μ, y, ny, k) - fullgradx0[i]

# The same, but exploiting that only dependencies in subsamples need to be evaluated
∇ϕmoving(t, x, θ, i, t′, F, A, At, μ, y, ny = m .- y, k = 5) = γ0*x[i] - fdot_moving(A, At, i, t, x, θ, t′, F, μ, y, ny, k) - fullgradx0[i]
# Tests, to be sure
#@test norm(ReverseDiff.gradient(x->ϕ(x, A, y), xtrue) - ∇ϕ(xtrue, A, At, y)) < 1e-7
#@test norm(ReverseDiff.gradient(x->ϕ(x, A, y), xtrue) - [∇ϕ(xtrue, i, A, At, y) for i in 1:p]) < 1e-7

# Note that the A'A has the same sparsity structure as the Hessian of ϕ
# We can use it to tune the doubly local ZigZig with spdmp sparsifying the
# precision matrix estimate
Γ0 = A'*A
nonzeros(Γ0) .= 1

Γ = hcat((sparse(ReverseDiff.gradient(x -> ∇ϕ(x, i, A, At, y), x0)) for i in 1:p)...)
μ = copy(x0)

Γdrop = droptol!(copy(Γ), 1e-2)
@assert nnz(diag(Γdrop)) == p
println("Precision Γ is ", 100*nnz(Γ)/length(Γ), "% filled")
println("Dropped precision is ", 100*nnz(Γdrop)/length(Γdrop), "% filled")

#defining tailored bounds
struct MyBoundLog
    c::Float64
end

c = [MyBoundLog(0.0) for i in 1:p]
rows = rowvals(A)
vals = nonzeros(A)
for i in 1:p
    C_i = 0.0
    Ni = length(nzrange(A, i))
    for j in nzrange(A, i)
        row = rows[j]
        C_i = max(C_i, Ni*0.25*abs(vals[j])*norm(A[row,:]))
    end
    c[i] = MyBoundLog(C_i)
end
# a and b
function ZZB.ab(G, i, x, θ, c::Vector{MyBoundLog}, F::ZigZag, args...)
    # NB: if ξ and θ are sparse vectors, this implemention is not efficient
    max(0, θ[i]*(fullgradx0[i] + γ0*x[i])) + abs(θ[i])*c[i].c*norm(x - μ), abs(θ[i])*c[i].c*norm(θ)
end

# Random initial values
t0 = 0.0
# Define ZigZag
σ = sqrt.(diag(inv(Matrix(Γ))))
#σ = (Vector(diag(Γ))).^(-0.5) # cheaper
μ = copy(x0)
Z = ZigZag(Γ, μ, σ, ρ=0.5, λref=0.00)
Zdiag = ZigZag(sparse(Diagonal(Γ)), μ, σ, ρ=0.5, λref=0.00)
Zdrop = ZigZag(Γdrop, μ, σ, ρ=0.5, λref=0.00)
θ0 = rand([-1.0,1.0], p) .* σ

# Run sparse ZigZag for T time units and collect trajectory
println("Distance starting point")
@show norm(x0 - xtrue)

println("Run spdmp")
#traj, u, (acc,num), c = @time spdmp(∇ϕr, t0, x0, θ0, T, c, Z, A, At, μ, y, m .- y, 12; adapt=true, factor=5)
# Probability of not being 0 with spike and slab prior
w = 1- wc
#kappa given the prior (here Gaussian) and w
κ = (γ0/sqrt(2π))/(1/w -1)

t0, x0, T = 0.0, randn(p), 1000.0
su = true #strong upperbounds
adapt = false

traj, (t, x, θ), (acc, num), c = @time ZZB.sspdmp(∇ϕr, t0, x0, θ0, T, c, Z, κ,
                                                A, At, μ, y, m .- y, 12;
                                                strong_upperbounds = su ,
                                                adapt = adapt)

error("")
#traj, u, (acc,num), c = @time spdmp(∇ϕ, t0, x0, θ0, T, c, Z, A, At, y, m .- y; adapt=true)
@show acc/num
# @show extrema(c ./ c0)

dt = T/4000
x̂ = mean(x for (t,x) in discretize(traj, dt))
X = Float64[]
ts = Float64[]
for (t,x) in discretize(traj, dt)
    append!(ts, t)
    append!(X, x)
end
X = reshape(X, p, length(X)÷p)

[xtrue x̂ abs.(xtrue - x̂)]

@show norm(xtrue - μ)
@show norm(xtrue - x̂)
error("")
println("Plot")
using CairoMakie, AbstractPlotting
using Colors
using GoldenSequences

ps = [4, 22, 50, p-3, p-2, p-1]
p0 = length(ps)
cs = map(x->RGB(x...), (Iterators.take(GoldenSequence(3), p0)))
pis = []
for i in 1:p0
    p_i = lines(ts, X[ps[i], :], color=cs[i])
    # lines!(p_i, [0, T], [xtrue[ps[i]], xtrue[ps[i]]], linewidth=2.0)
    # ylabel!(p_i, "var$(ps[i])")
    # xlabel!(p_i, "t")
    push!(pis, p_i)
end
p1 = title(hbox(pis...), "Sparse design logistic regression n=$(m*n), p=$p", textsize=20)
pis[6]
save(joinpath("figures","logistic$(typeof(Z).name).png"), p1)

p2 = title(vbox(hbox([text("$p", show_axis=false, textsize=10) for p in ps]...),
    [hbox([i==j ? lines(ts, X[i,:]) : lines(X[i, :], X[j, :]) for i in ps]...) for j in ps]...), "$ps")


text:
@testset "Coupling constraints" begin
    cp = convert(CoreModelCoupled, test_LP())
    @test size(cp.lm.S) == (4, 3)
    @test size(stoichiometry(convert(CoreModel, cp))) == (4, 3)
    new_cp = add_coupling_constraints(cp, stoichiometry(cp)[end, :], -1.0, 1.0)
    @test n_coupling_constraints(cp) + 1 == n_coupling_constraints(new_cp)

    new_cp =
        add_coupling_constraints(cp, stoichiometry(cp)[1:2, :], [-1.0; -1.0], [1.0; 1.0])
    @test n_coupling_constraints(cp) + 2 == n_coupling_constraints(new_cp)

    n_c = n_coupling_constraints(cp)
    add_coupling_constraints!(cp, stoichiometry(cp)[end, :], -1.0, 1.0)
    @test n_c + 1 == n_coupling_constraints(cp)
    add_coupling_constraints!(cp, stoichiometry(cp)[1:2, :], [-1.0; -1.0], [1.0; 1.0])
    @test n_c + 3 == n_coupling_constraints(cp)

    n_c = n_coupling_constraints(cp)
    remove_coupling_constraints!(cp, 1)
    @test n_c - 1 == n_coupling_constraints(cp)
    remove_coupling_constraints!(cp, [1, 2])
    @test n_c - 3 == n_coupling_constraints(cp)
    @test n_coupling_constraints(cp) == 0

    cp = test_coupledLP()
    n_c = n_coupling_constraints(cp)
    new_cp = remove_coupling_constraints(cp, 1)
    @test size(coupling(cp)) == (n_c, n_reactions(cp))
    @test n_c - 1 == n_coupling_constraints(new_cp)
    @test n_coupling_constraints(cp) == n_c
    new_cp = remove_coupling_constraints(cp, [1, 2])
    @test n_c - 2 == n_coupling_constraints(new_cp)
    new_cp = remove_coupling_constraints(cp, Vector(1:n_coupling_constraints(cp)))
    @test n_coupling_constraints(new_cp) == 0
    @test n_coupling_constraints(cp) == n_c

    cp = test_coupledLP()
    change_coupling_bounds!(cp, [3, 1], cl = [-10.0, -20], cu = [10.0, 20])
    cl, cu = coupling_bounds(cp)
    @test cl[[1, 3]] == [-20, -10]
    @test cu[[1, 3]] == [20, 10]
    change_coupling_bounds!(cp, [1000, 1001], cl = [-50.0, -60.0])
    cl, cu = coupling_bounds(cp)
    @test cl[[1000, 1001]] == [-50.0, -60.0]
end

@testset "Add reactions" begin
    cp = convert(CoreModelCoupled, test_LP())
    cp = add_coupling_constraints(cp, stoichiometry(cp)[end, :], -1.0, 1.0)

    new_cp = add_reactions(cp, 2.0 * ones(4), 3 .* ones(4), 2.0, -1.0, 1.0)
    @test new_cp isa CoreModelCoupled
    @test cp.C == new_cp.C[:, 1:end-1]
    @test cp.cl == new_cp.cl
    @test cp.cu == new_cp.cu

    new_cp = add_reactions(
        cp,
        2.0 * ones(4),
        3 .* ones(4),
        2.0,
        -1.0,
        1.0,
        "r4",
        ["m$i" for i = 1:4],
    )
    @test cp.C == new_cp.C[:, 1:end-1]
    @test cp.cl == new_cp.cl
    @test cp.cu == new_cp.cu

    new_cp = add_reactions(
        cp,
        2.0 * ones(4, 10),
        3 .* ones(4),
        2 .* ones(10),
        -ones(10),
        ones(10),
    )
    @test cp.C == new_cp.C[:, 1:end-10]
    @test cp.cl == new_cp.cl
    @test cp.cu == new_cp.cu

    new_cp = add_reactions(
        cp,
        2.0 * ones(4, 10),
        3 .* ones(4),
        2 .* ones(10),
        -ones(10),
        ones(10),
        ["r$i" for i = 1:10],
        ["m$i" for i = 1:4],
    )
    @test cp.C == new_cp.C[:, 1:end-7] # 3 reactions were already present
    @test cp.cl == new_cp.cl
    @test cp.cu == new_cp.cu

    new_cp =
        add_reactions(cp, 2.0 * sprand(4000, 0.5), 3 .* sprand(4000, 0.5), 2.0, -1.0, 1.0)
    @test cp.C == new_cp.C[:, 1:end-1]
    @test cp.cl == new_cp.cl
    @test cp.cu == new_cp.cu

    cm = CoreModel(
        2.0 * ones(4, 10),
        3 .* ones(4),
        2 .* ones(10),
        -ones(10),
        ones(10),
        ["r$i" for i = 1:10],
        ["m$i" for i = 1:4],
    )
    new_cp = add_reactions(cp, cm)
    @test cp.C == new_cp.C[:, 1:end-7] # 3 reactions were already present
    @test cp.cl == new_cp.cl
    @test cp.cu == new_cp.cu
end

@testset "Remove reactions" begin
    cp = convert(CoreModelCoupled, test_LP())
    cp = add_coupling_constraints(cp, 1.0 .* collect(1:n_reactions(cp)), -1.0, 1.0)

    new_cp = remove_reactions(cp, [3, 2])
    @test new_cp isa CoreModelCoupled
    @test new_cp.C[:] == cp.C[:, 1]
    @test new_cp.cl == cp.cl
    @test new_cp.cu == cp.cu

    new_cp = remove_reaction(cp, 2)
    @test new_cp.C == cp.C[:, [1, 3]]
    @test new_cp.cl == cp.cl
    @test new_cp.cu == cp.cu

    new_cp = remove_reaction(cp, "r1")
    @test new_cp.C == cp.C[:, 2:3]
    @test new_cp.cl == cp.cl
    @test new_cp.cu == cp.cu

    new_cp = remove_reactions(cp, ["r1", "r3"])
    @test new_cp.C[:] == cp.C[:, 2]
    @test new_cp.cl == cp.cl
    @test new_cp.cu == cp.cu

    new_cp = remove_reactions(cp, [1, 4])
    @test new_cp.C == cp.C[:, 2:3]

    new_cp = remove_reactions(cp, [1, 1, 2])
    @test new_cp.C[:] == cp.C[:, 3]
end

@testset "Change bounds" begin
    cp = convert(CoreModelCoupled, test_LP())
    @test cp isa CoreModelCoupled

    change_bound!(cp, 1, lower = -10, upper = 10)
    @test cp.lm.xl[1] == -10
    @test cp.lm.xu[1] == 10
    change_bounds!(cp, [1, 2]; lower = [-11, -12.2], upper = [11, 23.0])
    @test cp.lm.xl[2] == -12.2
    @test cp.lm.xu[1] == 11

    change_bound!(cp, "r1", lower = -101, upper = 101)
    @test cp.lm.xl[1] == -101
    @test cp.lm.xu[1] == 101
    change_bounds!(cp, ["r1", "r2"]; lower = [-113, -12.23], upper = [113, 233.0])
    @test cp.lm.xl[2] == -12.23
    @test cp.lm.xu[1] == 113

    new_model = change_bound(cp, 1, lower = -10, upper = 10)
    @test new_model.lm.xl[1] == -10
    @test new_model.lm.xu[1] == 10
    new_model = change_bounds(cp, [1, 2]; lower = [-11, -12.2], upper = [11, 23.0])
    @test new_model.lm.xl[2] == -12.2
    @test new_model.lm.xu[1] == 11

    new_model = change_bound(cp, "r1", lower = -101, upper = 101)
    @test new_model.lm.xl[1] == -101
    @test new_model.lm.xu[1] == 101
    new_model =
        change_bounds(cp, ["r1", "r2"]; lower = [-113, -12.23], upper = [113, 233.0])
    @test new_model.lm.xl[2] == -12.23
    @test new_model.lm.xu[1] == 113

end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

function Signature(ptr::Ptr{SignatureStruct})
    sig   = unsafe_load(ptr)::SignatureStruct
    name  = unsafe_string(sig.name)
    email = unsafe_string(sig.email)
    time   = sig.when.time
    offset = sig.when.offset
    return Signature(name, email, time, offset)
end
Signature(sig::GitSignature) = Signature(sig.ptr)

function Signature(name::AbstractString, email::AbstractString)
    sig_ptr_ptr = Ref{Ptr{SignatureStruct}}(C_NULL)
    @check ccall((:git_signature_now, :libgit2), Cint,
                 (Ptr{Ptr{SignatureStruct}}, Cstring, Cstring), sig_ptr_ptr, name, email)
    sig = GitSignature(sig_ptr_ptr[])
    s = Signature(sig.ptr)
    finalize(sig)
    return s
end

function Signature(repo::GitRepo)
    sig = default_signature(repo)
    s = Signature(sig.ptr)
    finalize(sig)
    return s
end

function Base.convert(::Type{GitSignature}, sig::Signature)
    sig_ptr_ptr = Ref{Ptr{SignatureStruct}}(C_NULL)
    @check ccall((:git_signature_new, :libgit2), Cint,
                 (Ptr{Ptr{SignatureStruct}}, Cstring, Cstring, Int64, Cint),
                 sig_ptr_ptr, sig.name, sig.email, sig.time, sig.time_offset)
    return GitSignature(sig_ptr_ptr[])
end

"""Return signature object. Free it after use."""
function default_signature(repo::GitRepo)
    sig_ptr_ptr = Ref{Ptr{SignatureStruct}}(C_NULL)
    @check ccall((:git_signature_default, :libgit2), Cint,
                 (Ptr{Ptr{SignatureStruct}}, Ptr{Void}), sig_ptr_ptr, repo.ptr)
    return GitSignature(sig_ptr_ptr[])
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/iup/include/iupim.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static is_windows()? (const libiupim_ = "iupim") : (const libiupim_ = "libiupim")  # Name of shared lib.

typealias Ihandle Void

function IupLoadImage(file_name::String)
	ccall((:IupLoadImage, libiupim_), Ptr{Ihandle}, (Ptr{UInt8},), file_name)
end
function IupSaveImage(ih::Ptr{Ihandle}, file_name::String, format::String)
	ccall((:IupSaveImage, libiupim_), Cint, (Ptr{Ihandle}, Ptr{UInt8}, Ptr{UInt8}), ih, file_name, format)
end

function IupGetNativeHandleImage(handle::Ptr{Void})
	ccall((:IupGetNativeHandleImage, libiupim_), Ptr{imImage}, (Ptr{Void},), handle)
end
function IupGetImageNativeHandle(image::Ptr{imImage})
	ccall((:IupGetImageNativeHandle, libiupim_), Ptr{Void}, (Ptr{imImage},), image)
end
function IupImageFromImImage(image::Ptr{imImage})
	ccall((:IupImageFromImImage, libiupim_), Ptr{Ihandle}, (Ptr{imImage},), image)
end


text:
const DB_ROOT = joinpath(artifact"refractiveindex.info", "refractiveindex.info-database-2021-07-18", "database")

struct MaterialEntry
    name::String
    path::String
    book_category::String
    page_category::String
end

function Base.show(io::IO, ::MIME"text/plain", g::MaterialEntry)
    println(io, typeof(g))
    println(io, "  name: ", g.name)
    println(io, "  book_category: ", g.book_category)
    println(io, "  page_category: ", g.page_category)
    print(io,   "  path: ", g.path)
end

const DB = Dict{Tuple{String,String,String},MaterialEntry}()

function __init__()

    lib = YAML.load_file(joinpath(DB_ROOT, "library.yml"), dicttype = Dict{String,Any})

    last_book_divider, last_page_divider = "", ""
    for shelf in lib
        shelfname = shelf["SHELF"]
        for book in shelf["content"]
            if haskey(book, "DIVIDER")
                last_book_divider = book["DIVIDER"]
                continue
            end
            bookname = book["BOOK"]
            for page in book["content"]
                if haskey(page, "DIVIDER")
                    last_page_divider = page["DIVIDER"]
                    continue
                end
                pagename = string(page["PAGE"])
                DB[(shelfname, bookname, pagename)] = MaterialEntry(
                    page["name"],
                    page["data"],
                    last_book_divider,
                    last_page_divider
                )
            end
        end
    end
end

const FORMULAS = Dict{String,Symbol}(
    "formula 1" => :Sellmeier,
    "formula 2" => :Sellmeier2,
    "formula 3" => :Polynomial,
    "formula 4" => :RIInfo,
    "formula 5" => :Cauchy,
    "formula 6" => :Gases,
    "formula 7" => :Herzberger,
    "formula 8" => :Retro,
    "formula 9" => :Exotic,
)

function str2tuple(str)
    arr = Base.parse.(Float64, split(str))
    ntuple(i -> arr[i], length(arr))
end

function parse(dict)

    reference = get(dict, "REFERENCES", "")
    comment = get(dict, "COMMENTS", "")
    specs = get(dict, "SPECS", Dict{Any,Any}())

    meta = Metadata(reference, comment, specs)

    data = dict["DATA"]
    N = length(data)

    if N == 1

        if data[1]["type"] in keys(FORMULAS)

            coeffs = str2tuple(data[1]["coefficients"])
            λrange = str2tuple(data[1]["wavelength_range"])
            n = eval(FORMULAS[data[1]["type"]])(λrange, coeffs)

            return FormulaN(meta, n)

        elseif data[1]["type"] == "tabulated n"

            raw = readdlm(IOBuffer(data[1]["data"]), ' ', Float64)
            λ = convert(Vector{Float64}, raw[:,1])
            n = convert(Vector{Float64}, raw[:,2])

            return TabulatedN(meta, λ, n)

        elseif data[1]["type"] == "tabulated nk"

            raw = readdlm(IOBuffer(data[1]["data"]), ' ', Float64)
            λ = convert(Vector{Float64}, raw[:,1])
            n = convert(Vector{Float64}, raw[:,2])
            k = convert(Vector{Float64}, raw[:,3])

            return TabulatedNK(meta, λ, n, k)

        end

    elseif N == 2

        if data[1]["type"] in keys(FORMULAS) && data[2]["type"] == "tabulated k"

            coeffs = str2tuple(data[1]["coefficients"])
            λrange = str2tuple(data[1]["wavelength_range"])
            n = eval(FORMULAS[data[1]["type"]])(λrange, coeffs)

            raw = readdlm(IOBuffer(data[2]["data"]), ' ', Float64)
            λ = convert(Vector{Float64}, raw[:,1])
            k = convert(Vector{Float64}, raw[:,2])

            return FormulaNTabulatedK(meta, n, λ, k)

        elseif data[1]["type"] == "tabulated n" && data[2]["type"] == "tabulated k"

            raw = readdlm(IOBuffer(data[1]["data"]), ' ', Float64)
            λ_n = convert(Vector{Float64}, raw[:,1])
            n = convert(Vector{Float64}, raw[:,2])

            raw = readdlm(IOBuffer(data[2]["data"]), ' ', Float64)
            λ_k = convert(Vector{Float64}, raw[:,1])
            k = convert(Vector{Float64}, raw[:,2])

            return TabulatedNTabulatedK(meta, λ_n, n, λ_k, k)
        end
    end

    error("unsupported data type")

end

"""
    get_material(shelf, book, page)

Load the refractive index data for the material corresponding to the specified
shelf, book, and page within the [refractiveindex.info](https://refractiveindex.info/) database. The data
can be queried by calling the returned object at a given wavelength.

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

julia> MgLiTaO3 = get_material("other", "Mg-LiTaO3", "Moutzouris-o")
"Mg-LiTaO3 (Moutzouris et al. 2011: n(o) 0.450-1.551 µm; 8 mol.% Mg)"

julia> MgLiTaO3(0.45) # default unit is microns
2.2373000025056826

julia> using Unitful

julia> MgLiTaO3(450u"nm") # auto-conversion from generic Unitful.jl length units
2.2373000025056826

julia> MgLiTaO3(450e-9, "m") # strings can be used to specify units (parsing is cached)
2.2373000025056826
```
"""
function get_material(shelf, book, page)
    item = DB[(shelf, book, page)]
    return load_file(joinpath(DB_ROOT, "data", item.path))
end

"""
    load_file(path)

Load the refractive index data for the material corresponding to a YAML
file within the [refractiveindex.info](https://refractiveindex.info/) database. The data
can be queried by calling the returned object at a given wavelength.

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

julia> path = joinpath(RefractiveIndexDatabase.DB_ROOT, "data", "main", "Ar", "Peck-15C.yml");

julia> Ar = load_file(path)
RefractiveIndexDatabase.FormulaN

julia> Ar(532, "nm")
1.0002679711455778
```
"""
function load_file(path)

    dict = YAML.load_file(path)
    return parse(dict)

end

"""
    load_url(path)

Load the refractive index data for the material corresponding to the specified
url within the [refractiveindex.info](https://refractiveindex.info/) database. The data
can be queried by calling the returned object at a given wavelength.

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

julia> Ar = load_url("https://refractiveindex.info/database/data/main/Ar/Peck-15C.yml")
RefractiveIndexDatabase.FormulaN

julia> Ar(532, "nm")
1.0002679711455778
```
"""
function load_url(path)

    r = request("GET", path)
    dict = YAML.load(String(r.body))
    return parse(dict)

end

"""
    search(shelf="", book="", page="")

Returns with a set of keys for all materials in the database that match the given query.

```julia-repl
julia> RefractiveIndexDatabase.search()
Set{Tuple{String, String, String}} with 2915 elements:
  ("glass", "HIKARI-F", "F1")
  ("main", "Au", "Ciesielski")
  ("glass", "SUMITA-PG", "K-PG375")
  ("glass", "SUMITA-VC", "K-VC179")
  ("glass", "SUMITA-SFLD", "K-SFLD11")
  ("main", "SrF2", "Li")
  ("glass", "SCHOTT-SK", "N-SK11")
  ("main", "BaB2O4", "Zhang-o")
  ("glass", "OHARA-BAH", "BAH54")
  ("glass", "HIKARI-LaSF", "E-LASFH9")
  ("main", "ZnTe", "Marple")
  ("main", "Si", "Green-2008")
  ("glass", "HOYA-TaF", "M-TAF31")
  ("main", "C", "Djurisic-e")
  ("organic", "propanol", "Sani")
  ("glass", "BARBERINI-sun-protection", "D6420")
  ("3d", "crystals", "germanium")
  ("glass", "SCHOTT-SK", "P-SK60")
  ("glass", "SUMITA-GIR", "K-GIR79")
  ("glass", "SCHOTT-F", "F2HT")
  ("glass", "OHARA-LAM", "LAM59")
  ("main", "He", "Borzsonyi")
  ⋮
```
```julia-repl
julia> RefractiveIndexDatabase.search(book="Ag", shelf="main")
Set{Tuple{String, String, String}} with 31 elements:
  ("main", "AgGaSe2", "Harasaki-o")
  ("main", "Ag", "Johnson")
  ("main", "Ag", "Rakic-BB")
  ("main", "AgCl", "Tilton")
  ("main", "Ag", "Jiang")
  ("main", "Ag", "Choi")
  ("main", "AgGaS2", "Boyd-o")
  ("main", "Ag3AsS3", "Hulme-e")
  ("main", "AgGaS2", "Kato-e")
  ("main", "AgGaS2", "Kato-o")
  ("main", "Ag", "Ciesielski-Ge")
  ("main", "Ag", "McPeak")
  ("main", "AgGaSe2", "Boyd-e")
  ("main", "Ag", "Ciesielski-Ni")
  ("main", "Ag", "Werner")
  ("main", "AgGaS2", "Boyd-e")
  ("main", "AgBr", "Schröter")
  ("main", "AgGaS2", "Takaoka-e")
  ("main", "Ag", "Werner-DFT")
  ("main", "Ag", "Yang")
  ("main", "Ag", "Stahrenberg")
  ("main", "Ag", "Windt")
  ⋮
```
"""
function search(; shelf::String = "", book::String = "", page::String = "")

    if isempty(shelf) && isempty(book) && isempty(page)
        return Set(keys(DB))
    end

    return Set(filter((key) -> isin(shelf, key[1]) && isin(book, key[2]) && isin(page, key[3]), keys(DB)))
end

isin(query::String, key::String) = isempty(query) || occursin(query, key)

"""
    info(key)

Return with the metadata for the material corresponding to the given key.

```julia-repl
julia> info("main", "Ag", "Johnson")
RefractiveIndexDatabase.MaterialEntry
  name: Johnson and Christy 1972: n,k 0.188-1.94 µm
  book_category: Ag - Silver
  page_category: Experimental data: bulk, thick film
  path: main/Ag/Johnson.yml
```
"""
function info(key...)
    return DB[key]
end


text:
# # A sphere in a homogenous field

# For a three axial ellipsoid in a constant field, we have an analytic solution expressed in terms of demagnetization coefficients $n_i$:
# ```math
# H_i = \frac{H_{0i}}{(1 + (\mu-1)n_i)}
# ```
# whre $\vec H$ is the interior field, and $\vec H_{0}$ is the constant field. That makes ellipsoid a perfect subject to test the numerics.

using QuadGK

function ellipsoid_demagnetization_coefficients(a,b,c)

    UP_BOUND = 1000

    Ru2(u) = (u+a^2)*(u+b^2)*(u+c^2)

    nx = 1/2 * a*b*c * quadgk(s -> 1/(s+a^2)/sqrt(Ru2(s)), 0, UP_BOUND)[1]
    ny = 1/2 * a*b*c * quadgk(s -> 1/(s+b^2)/sqrt(Ru2(s)), 0, UP_BOUND)[1]
    nz = 1/2 * a*b*c * quadgk(s -> 1/(s+c^2)/sqrt(Ru2(s)), 0, UP_BOUND)[1]

    return [nx, ny, nz]
end

function EllipsoidField(a,b,c,mu,H0)

    H0x, H0y, H0z = H0
    nx, ny, nz = ellipsoid_demagnetization_coefficients(a,b,c)

    Hix = H0x/(1 + (mu-1)*nx)
    Hiy = H0y/(1 + (mu-1)*ny)
    Hiz = H0z/(1 + (mu-1)*nz)

    return [Hix,Hiy,Hiz]
end;

# We define that we have a field of strength $H_0 = e_z$. A corresponding potential for that is $\psi_0 = z$, which we need to pass to the solver.

using LaplaceBIE
using LinearAlgebra

H0 = [0,0,1]

ψ0(x) = dot(x,H0)
∇ψ0(x) = H0;

# For simplicity, we consider a sphere with relative permittivity $\epsilon=15$. The mesh for the sphere we generate from subdivisions of icosahedron given in sphere.jl. We also need normal vectors which we get from the vertex positions with high accuracy

ϵ = 15

include("sphere.jl")
msh = unitsphere(3)
vertices, faces = msh.vertices, msh.faces
n = vertices;

# At the moment users need to take the derivative of the potential themselves which can be done in multiple ways. It is also possible that a user might wish to calculate the field due to the surface current for which `normalderivatives` method could be useful.

# Now we can proceed with calculation. To calculate surface potential everywhere on the surface one executes a method
ψ = surfacepotential(vertices,n,faces,ϵ,ψ0);
# which solves a regularized boundary integral equation with BLAS. That then, for example, can be used to calculate the energy of the field.

# Usually, however, one wants to know the field on the surface for force calculations due to $(M \cdot \nabla)H$. The library offers to do a finite differentiation on the calculated potential for calculating tangential field components. That can be easily achieved with method `tangentderivatives`:

P∇ψ = tangentderivatives(vertices,n,faces,ψ);

# And lastly to calculate the normal derivatives, we can use Biot-Savarat law treating the tangential field as raising from a surface current:

n∇ψ = normalderivatives(vertices,n,faces,P∇ψ,ϵ,∇ψ0);

# We can compare numerics with analytics easily by the use of azimuthal symmetry in the z-direction. The normal field is equal to potential since vertex positions coincide with normals which is thus not shown twice. And so the comparison:

Hin = EllipsoidField(1,1,1,ϵ,H0)

x = [v[3] for v in vertices]
sp = sortperm(x)

for xkey in sp[2:43:end]
    psit = dot(Hin,vertices[xkey])
    println("$(vertices[xkey][3]) \t $psit \t $(ψ[xkey]) \t $(n∇ψ[xkey])")
end


text:
using Polynomials
@testset "Sastre" begin
    α = big(1.2)

    ## Test approximation of exponential and number of multiplications
    for (k, dk, Ak, m) in [
        (0, 1, 0.0001, :ps_degopt),
        (1, 2, 0.001, :ps_degopt),
        (2, 4, 0.1, :ps_degopt),
        (3, 8, 1.1, :y1s),
        (4, 12, 2.9, :y1s),
        (6, 16, 5.1, :h2m),
        (8, 30, 16.5, :z1ps),
    ] # Where do I find these limits?
        for i = 1:2
            @testset "k = $k, and i=$i" begin
                if i == 1
                    (graph, cref) = graph_sastre_exp(k, m)
                else
                    (graph, cref) = graph_sastre_exp(k, :auto)
                end
                err1 = eval_graph(graph, big(Ak)) - exp(big(Ak))
                err2 = eval_graph(graph, big(Ak) / α) - exp(big(Ak) / α)
                @test -log(abs(err2 / err1)) / log(α) > dk + 1
                @test sum(values(graph.operations) .== :mult) == k
            end
        end
    end

    ## Test the degree-9 solver graph_sastre_poly(b) by comparison to tabulated
    ## values.
    # 1: Table 4 against equations (16)-(32) for exponential
    (graph1, cref1) = graph_sastre_exp(3, :y1s)
    (graph2, cref2) = graph_sastre_poly(1 ./ factorial.(0:8))
    @test sum(values(graph2.operations) .== :mult) == 3
    @test all(cref1 .== cref2)
    for c in cref1
        @test graph1.coeffs[c[1]][c[2]] ≈ graph2.coeffs[c[1]][c[2]]
    end

    # 2: Table 11 against Table 10 and equations (16)-(32)
    (graph_exp, _) = graph_sastre_exp(6, :h2m)
    A = randn(100, 100) / 20 * 5.3
    # Naive implementation from coefficients in Table 10
    b8 = 2.186201576339059e-7
    bp8 = 2.186201576339059e-7
    b7 = 9.839057366529322e-7
    bp7 = 2.514016785489562e-6
    b6 = 1.058964584814256e-5
    bp6 = 3.056479369585950e-5
    b5 = 1.554700173279057e-4
    bp5 = 3.197607034851565e-4
    b4 = 2.256892506343887e-3
    bp4 = 2.585006547542889e-3
    b3 = 2.358987357109499e-2
    bp3 = 1.619043970183846e-2
    b2 = 1.673139636901279e-1
    bp2 = 8.092036376147299e-2
    b1 = 7.723603212944010e-1
    bp1 = 3.229486011362677e-1
    b0 = 3.096467971936040
    bp0 = 0.0
    β0 = 1
    (graph, _) = graph_sastre_poly([b0, b1, b2, b3, b4, b5, b6, b7, b8])
    (graphp, _) =
        graph_sastre_poly([bp0, bp1, bp2, bp3, bp4, bp5, bp6, bp7, bp8])
    v1 = eval_graph(graph, A)
    vp1 = eval_graph(graphp, A)
    @test v1 * vp1 + β0 * I ≈ eval_graph(graph_exp, A)

    v1 = eval_graph(graph, big(5.3))
    vp1 = eval_graph(graphp, big(5.3))
    err1 = (v1 * vp1 + β0 * I) - exp(big(5.3))

    v1 = eval_graph(graph, big(5.3) / α)
    vp1 = eval_graph(graphp, big(5.3) / α)
    err2 = (v1 * vp1 + β0 * I) - exp(big(5.3) / α)
    @test -log(abs(err2 / err1)) / log(α) > 16 + 1

    ## Test conversion of (62)-(65) to Degopt
    # 1: By using numbers from (57)-(59)
    c16 = 4.018761610201036e-4
    c8  = 2.116367017255747e0
    c15 = 2.945531440279683e-3
    c7  = -5.792361707073261e0
    c14 = 8.712167566050691e-2
    c6  = -1.491449188999246e-1
    c13 = 4.017568440673568e-1
    c5  = 1.040801735231354e1
    c12 = -6.352311335612147e-2
    c4  = -6.331712455883370e1
    c11 = 2.684264296504340e-1
    c3  = 3.484665863364574e-1
    c10 = 1.857143141426026e1
    c2  = -1.224230230553340e-1
    c9  = 2.381070373870987e-1
    c1  = 1

    c = [
        [[c15, c16], [0.0, 0, 0]], # (57)
        [[1.0], [0, c13, c14],
         [1.0], [c11, 0, c12],
         [c10], [0.0, 0, 0]], # (58)
        [[0.0, 1], [0, c8, c9],
         [c7, 1], [0, c6, 0],
         [c4, c5], [c1, c2, c3]], # (59)
    ]

    (graph, cref) = graph_sastre_yks_degopt(2, 2, c)
    err1 = eval_graph(graph, big(3.1)) - exp(big(3.1))
    err2 = eval_graph(graph, big(3.1) / α) - exp(big(3.1) / α)
    @test -log(abs(err2 / err1)) / log(α) > 15 + 1

    x = Polynomial("x")
    CC = coeffs(eval_graph(Compgraph(Any, graph), x))
    t = (CC[17] - 1 / factorial(big(16))) * factorial(big(16)) - (-0.454) # (60)
    @test (abs(t) < 1e-3) && (abs(t) > 1e-4) # (60)
    @test all(abs.((CC.*factorial.(big.(0:16)))[1:16] .- 1) .<= 2e-15)

    # 2: By constructing an onw polynomial of the form y_33
    c1  = 1.1e-2
    c2  = 1.2e-2
    c3  = 1.3e-2
    c4  = 1.4e-2
    c5  = 1.5e-2
    c6  = 1.6e-2
    c7  = 1.7e-2
    c8  = 2.05e-2
    c9  = 2.1e-2
    c10 = 2.15e-2
    c11 = 2.2e-2
    c12 = 2.25e-2
    c13 = 2.3e-2
    c14 = 2.35e-2
    c15 = 2.4e-2
    c16 = 2.45e-2
    c17 = 2.5e-2
    c18 = 2.55e-2
    c19 = 2.6e-2
    c20 = 2.65e-2
    c21 = 2.7e-2
    c22 = 2.75e-2
    c23 = 3.05e-2
    c24 = 3.1e-2
    c25 = 3.15e-2
    c26 = 3.2e-2
    c27 = 3.25e-2
    c28 = 3.3e-2
    c29 = 3.35e-2
    c30 = 3.4e-2
    c31 = 3.45e-2
    c32 = 3.5e-2
    c33 = 3.55e-2
    c34 = 3.6e-2
    c35 = 3.65e-2
    c36 = 3.7e-2
    c37 = 3.75e-2
    c38 = 3.8e-2
    c39 = 3.85e-2
    c40 = 3.9e-2
    c41 = 4.05e-2
    c42 = 4.1e-2
    c43 = 4.15e-2
    c44 = 4.2e-2
    c45 = 4.25e-2
    c46 = 4.3e-2
    c47 = 4.35e-2
    c48 = 4.4e-2
    c49 = 4.45e-2
    c50 = 4.5e-2
    c51 = 4.55e-2
    c52 = 4.6e-2
    c53 = 4.65e-2
    c54 = 4.7e-2
    c55 = 4.8e-2
    c56 = 4.85e-2
    c57 = 4.9e-2
    c58 = 4.95e-2
    c59 = 5e-2
    c60 = 5.05e-2
    c61 = 5.10e-2

    c = [
        [[c1, c2, c3], [c4, c5, c6, c7]],
        [
            [c8],
            [c9, c10, c11, c12],
            [c13],
            [c14, c15, c16, c17],
            [c18],
            [c19, c20, c21, c22],
        ],
        [
            [c23, c24],
            [c25, c26, c27, c28],
            [c29, c30],
            [c31, c32, c33, c34],
            [c35, c36],
            [c37, c38, c39, c40],
        ],
        [
            [c41, c42, c43],
            [c44, c45, c46, c47],
            [c48, c49, c50],
            [c51, c52, c53, c54],
            [c55, c56, c57],
            [c58, c59, c60, c61],
        ],
    ]

    II = Matrix(1.0I, 100, 100)
    A2 = A * A
    A3 = A2 * A
    y03 =
        A3 * (c1 * A + c2 * A2 + c3 * A3) + c4 * II + c5 * A + c6 * A2 + c7 * A3
    (graph, cref) = graph_sastre_yks_degopt(0, 3, c[1:1])
    @test eval_graph(graph, A) ≈ y03

    y13 =
        (c8 * y03 + c9 * II + c10 * A + c11 * A2 + c12 * A3) *
        (c13 * y03 + c14 * II + c15 * A + c16 * A2 + c17 * A3) +
        c18 * y03 +
        c19 * II +
        c20 * A +
        c21 * A2 +
        c22 * A3
    (graph, cref) = graph_sastre_yks_degopt(1, 3, c[1:2])
    @test eval_graph(graph, A) ≈ y13

    y23 =
        (c23 * y03 + c24 * y13 + c25 * II + c26 * A + c27 * A2 + c28 * A3) *
        (c29 * y03 + c30 * y13 + c31 * II + c32 * A + c33 * A2 + c34 * A3) +
        c35 * y03 +
        c36 * y13 +
        c37 * II +
        c38 * A +
        c39 * A2 +
        c40 * A3
    (graph, cref) = graph_sastre_yks_degopt(2, 3, c[1:3])
    @test eval_graph(graph, A) ≈ y23

    y33 =
        (
            c41 * y03 +
            c42 * y13 +
            c43 * y23 +
            c44 * II +
            c45 * A +
            c46 * A2 +
            c47 * A3
        ) * (
            c48 * y03 +
            c49 * y13 +
            c50 * y23 +
            c51 * II +
            c52 * A +
            c53 * A2 +
            c54 * A3
        ) +
        c55 * y03 +
        c56 * y13 +
        c57 * y23 +
        c58 * II +
        c59 * A +
        c60 * A2 +
        c61 * A3
    (graph, cref) = graph_sastre_yks_degopt(3, 3, c)
    @test eval_graph(graph, A) ≈ y33
end


text:
# Import package and any files here ↴
include("Muffins.jl")

# Make sure test.txt is in the same folder as test.jl
lines = open("test.txt") do file
    readlines(file)
end

# Array of failed cases: [(muffins, students), method, generated answer, actual answer]
incorrect = []

# Array of cases with no solutions over all methods: [(muffins, students), methods]
nosol = []

# Array of skipped cases: [(muffins, students), method]
skipped = []

for case in lines
    println(case)
    (m, s, methods, alpha) = tuple(split(case, " ")...)

    m = parse(Int64, m)
    s = parse(Int64, s)

    alphaF = split(alpha, "/")
    alpha = parse(Int64, alphaF[1]) // parse(Int64, alphaF[2])

    for method in split(methods, ",")
        res = 0

        if method == "FC"
            res = Muffins.fc(m, s, output = 0)
        end

        if method == "HALF"
            res = Muffins.half(m, s, output = 0)
        end

        if method == "INT"
            res = Muffins.int(m, s, output = 0)
        end

        if method == "MID"
            res = Muffins.mid(m, s, output = 0)
        end

        if startswith(method, "EBM")
            res = Muffins.ebm(m, s, output = 0)
        end

        if startswith(method, "HBM")
            res = Muffins.hbm(m, s, output = 0)
        end

        if method == "GAP"
            res = Muffins.gap(m, s, output = 0)
        end

        if res != 0
            if res != alpha
                append!(incorrect, [[(m, s), method, res, alpha]])
            end
        else
            append!(skipped, [[(m, s), method]])
        end
    end

    if Muffins.muffins(m, s, output = 0) != alpha
        append!(nosol, [[(m, s), methods]])
    end

end

# Customize output here ↴


text:
using ColorTypes

const Palette = Vector{RGB24}

rgb24(x::UInt32) = RGB24(((x >> 16) & 0xFF) / 255.0, ((x >> 8) & 0xFF) / 255.0, (x & 0xFF) / 255.0)

const DefaultPalette = vcat(
  rgb24(0x666666), rgb24(0x002A88), rgb24(0x1412A7), rgb24(0x3B00A4),
  rgb24(0x5C007E), rgb24(0x6E0040), rgb24(0x6C0600), rgb24(0x561D00),
  rgb24(0x333500), rgb24(0x0B4800), rgb24(0x005200), rgb24(0x004F08),
  rgb24(0x00404D), rgb24(0x000000), rgb24(0x000000), rgb24(0x000000),
  rgb24(0xADADAD), rgb24(0x155FD9), rgb24(0x4240FF), rgb24(0x7527FE),
  rgb24(0xA01ACC), rgb24(0xB71E7B), rgb24(0xB53120), rgb24(0x994E00),
  rgb24(0x6B6D00), rgb24(0x388700), rgb24(0x0C9300), rgb24(0x008F32),
  rgb24(0x007C8D), rgb24(0x000000), rgb24(0x000000), rgb24(0x000000),
  rgb24(0xFFFEFF), rgb24(0x64B0FF), rgb24(0x9290FF), rgb24(0xC676FF),
  rgb24(0xF36AFF), rgb24(0xFE6ECC), rgb24(0xFE8170), rgb24(0xEA9E22),
  rgb24(0xBCBE00), rgb24(0x88D800), rgb24(0x5CE430), rgb24(0x45E082),
  rgb24(0x48CDDE), rgb24(0x4F4F4F), rgb24(0x000000), rgb24(0x000000),
  rgb24(0xFFFEFF), rgb24(0xC0DFFF), rgb24(0xD3D2FF), rgb24(0xE8C8FF),
  rgb24(0xFBC2FF), rgb24(0xFEC4EA), rgb24(0xFECCC5), rgb24(0xF7D8A5),
  rgb24(0xE4E594), rgb24(0xCFEF96), rgb24(0xBDF4AB), rgb24(0xB3F3CC),
  rgb24(0xB5EBF2), rgb24(0xB8B8B8), rgb24(0x000000), rgb24(0x000000)
)

const width = 256
const height = 240
const numcolors = 64

# Frame containing pixel IDs.
const Frame = Matrix{UInt8}
# Rendered frame containing pixel RGB values.
const Screen = Matrix{RGB24}

Frame() = zeros(UInt8, height, width)
Screen() = fill(RGB24(0x00), (height, width))

# Render the frame as a screen with RGB pixel values using the given palette.
function render!(f::Frame, s::Screen, p::Palette)
  @assert size(f) == (height, width)
  @assert size(s) == (height, width)
  @assert size(p) == (numcolors,)
  for i = 1:width, j = 1:height
    index = f[j, i] + 1
    s[j, i] = p[index]
  end
end

render!(f::Frame, s::Screen) = render!(f, s, DefaultPalette)

function Screen(f::Frame, p::Palette)
  s = Screen()
  render!(f, s, p)
  s
end

Screen(f::Frame) = Screen(f, DefaultPalette)


text:
module FeatureSelectors

@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), "```julia" => "```jldoctest")
end FeatureSelectors

using Statistics
using StatsBase
using DataFrames
using HypothesisTests
using RDatasets

include("utils.jl")
include("UnivariateFeatureSelector.jl")

#! format: off
export
    UnivariateFeatureSelector

export
    select_features,
    # Measurements
    pearson_correlation,
    f_test,
    chisq_test,
    # Utils
    one_hot_encode,
    calculate_feature_importance

#! format: on
end # module


text:
include("Jcampdx.jl")

export BrukerFile, BrukerFileMeas, BrukerFileCalib, BrukerFileFast, latin1toutf8,
       sfPath, rawDataLengthConsistent

function latin1toutf8(str::AbstractString)
  buff = Char[]
  for c in Vector{UInt8}(str)
    push!(buff,c)
  end
  string(buff...)
end

function latin1toutf8(str::Nothing)
  error("Can not convert noting to UTF8!")
end

abstract type BrukerFile <: MPIFile end

mutable struct BrukerFileMeas <: BrukerFile
  path::String
  params::JcampdxFile
  paramsProc::JcampdxFile
  methodRead::Bool
  acqpRead::Bool
  visupars_globalRead::Bool
  recoRead::Bool
  methrecoRead::Bool
  visuparsRead::Bool
  mpiParRead::Bool
  maxEntriesAcqp::Int
  keylistAcqp::Vector{String}
  keylistMethod::Vector{String}
  pretendToBeSinglePatch::Bool
end

mutable struct BrukerFileCalib <: BrukerFile
  path::String
  params::JcampdxFile
  paramsProc::JcampdxFile
  methodRead::Bool
  acqpRead::Bool
  visupars_globalRead::Bool
  recoRead::Bool
  methrecoRead::Bool
  visuparsRead::Bool
  mpiParRead::Bool
  maxEntriesAcqp::Int
  keylistAcqp::Vector{String}
  keylistMethod::Vector{String}
  pretendToBeSinglePatch::Bool
end

function _iscalib(path::AbstractString)
    calib = false
    acqpPath = joinpath(path,"acqp")
    # if there is a file ismeasurement in path treat data in path as measurement
    if isfile(joinpath(path,"ismeasurement"))
        return calib
    end
    # else use the information provided within the acqp file
    if isfile(acqpPath)
        open(acqpPath, "r") do io
            for line in eachline(io)
                if !isnothing(findfirst("MPICalibration",line))
                    calib = true
                    break
                end
            end
        end
    end
    return calib
end

function BrukerFile(path::String; isCalib=_iscalib(path), fastMode=false, pretendToBeSinglePatch=false)
  if fastMode
    maxEntriesAcqp = 400
	  keylistAcqp = ["ACQ_scan_name", "ACQ_size", "ACQ_jobs","ACQ_MPI_drive_field_strength",
				 "ACQ_MPI_selection_field_gradient","NA", "ACQ_operator",
				 "ACQ_time", "ACQ_ReceiverSelect"]
	  keylistMethod = ["PVM_MPI_Bandwidth", "MPI_RepetitionsPerStep",
				 "PVM_MPI_NrCalibrationScans","MPI_NSteps",
				 "PVM_MPI_NrBackgroundMeasurementCalibrationAdditionalScans",
				 "PVM_MPI_ChannelSelect", "PVM_MPI_Tracer", "PVM_MPI_DriveFieldCycle", "PVM_Matrix", "PVM_Fov"]
  else
    maxEntriesAcqp = 2000
	  keylistAcqp = String[]
	  keylistMethod = String[]
	end
  params = JcampdxFile()
  paramsProc = JcampdxFile()

  if isCalib
    return BrukerFileCalib(path, params, paramsProc, false, false, false,
               false, false, false, false, maxEntriesAcqp, keylistAcqp,
               keylistMethod, pretendToBeSinglePatch)
  else
    return BrukerFileMeas(path, params, paramsProc, false, false, false,
               false, false, false, false, maxEntriesAcqp, keylistAcqp,
               keylistMethod, pretendToBeSinglePatch)
  end
end

function BrukerFile()
  params = JcampdxFile()
  paramsProc = JcampdxFile()
  return BrukerFileMeas("", params, paramsProc, false, false, false,
             false, false, false, false, 1, String[], String[], false)
end

BrukerFileFast(path) = BrukerFile(path, fastMode=true)

function getindex(b::BrukerFile, parameter)#::String
  if !b.acqpRead && ( parameter=="NA" || parameter[1:3] == "ACQ" )
    acqppath = joinpath(b.path, "acqp")
    read(b.params, acqppath, b.keylistAcqp, maxEntries=b.maxEntriesAcqp)
    b.acqpRead = true
  elseif !b.methodRead && length(parameter) >= 3 &&
         (parameter[1:3] == "PVM" || parameter[1:3] == "MPI")
    methodpath = joinpath(b.path, "method")
    read(b.params, methodpath, b.keylistMethod)
    b.methodRead = true
  elseif !b.visupars_globalRead && length(parameter) >= 4 &&
         parameter[1:4] == "Visu"
    visupath = joinpath(b.path, "visu_pars")
    if isfile(visupath)
      keylist = ["VisuStudyId","VisuStudyNumber","VisuExperimentNumber",
		 "VisuSubjectId","VisuSubjectName","VisuStudyDate", "VisuUid","VisuStudyUid"]
      read(b.params, visupath,keylist)
      b.visupars_globalRead = true
    end
  elseif !b.mpiParRead && length(parameter) >= 6 &&
         parameter[1:6] == "CONFIG"
    mpiParPath = joinpath(b.path, "mpi.par")
    if isfile(mpiParPath)
      read(b.params, mpiParPath)
      b.mpiParRead = true
    end
  end

  if haskey(b.params, parameter)
    return b.params[parameter]
  else
    return ""
  end
end

function getindex(b::BrukerFile, parameter, procno::Int64)#::String
  if !b.recoRead && lowercase( parameter[1:4] ) == "reco"
    recopath = joinpath(b.path, "pdata", string(procno), "reco")
    read(b.paramsProc, recopath, maxEntries=13)
    b.recoRead = true
  elseif !b.methrecoRead && parameter[1:3] == "PVM"
    methrecopath = joinpath(b.path, "pdata", string(procno), "methreco")
    read(b.paramsProc, methrecopath)
    b.methrecoRead = true
  elseif !b.visuparsRead && parameter[1:4] == "Visu"
    visuparspath = joinpath(b.path, "pdata", string(procno), "visu_pars")
    if isfile(visuparspath)
      read(b.paramsProc, visuparspath)
      b.visuparsRead = true
    end
  end

  return b.paramsProc[parameter]
end

# Helper
activeChannels(b::BrukerFile) = [parse(Int64,s) for s=b["PVM_MPI_ActiveChannels"]]
selectedChannels(b::BrukerFile) = b["PVM_MPI_ChannelSelect"] .== "Yes"
selectedReceivers(b::BrukerFile) = b["ACQ_ReceiverSelect"] .== "Yes"

# general parameters
version(b::BrukerFile) = nothing
uuid(b::BrukerFile) = nothing
time(b::BrukerFile) = nothing

# study parameters
studyName(b::BrukerFile) = latin1toutf8(b["VisuStudyId"])
# old study name
studyNameOld(b::BrukerFile) = string(latin1toutf8(b["VisuSubjectId"])*latin1toutf8(b["VisuSubjectName"]),"_",
                                  latin1toutf8(b["VisuStudyId"]),"_",
                                  b["VisuStudyNumber"])
studyNumber(b::BrukerFile) = parse(Int64,b["VisuStudyNumber"])
function studyUuid(b::BrukerFile)
  rng = MersenneTwister(hash(b["VisuStudyUid"])) # use VisuStudyUid as seed to generate uuid4
  return uuid4(rng)
end
studyDescription(b::BrukerFile) = "n.a."
function studyTime(b::BrukerFile)
  m = match(r"<(.+)\+",b["VisuStudyDate"])
  timeString = replace(m.captures[1],"," => ".")
  return DateTime( timeString )
end

# study parameters
experimentName(b::BrukerFile) = latin1toutf8(b["ACQ_scan_name"])
experimentNumber(b::BrukerFile) = parse(Int64,b["VisuExperimentNumber"])
function experimentUuid(b::BrukerFile)
  rng = MersenneTwister(hash(b["VisuUid"])) # use VisuUid as seed to generate uuid4
  return uuid4(rng)
end
experimentDescription(b::BrukerFile) = latin1toutf8(b["ACQ_scan_name"])
function experimentSubject(b::BrukerFile)
  id = latin1toutf8(b["VisuSubjectId"])
  name = latin1toutf8(b["VisuSubjectName"])
  if id == name
    return id
  else
    return id*name
  end
end
experimentIsSimulation(b::BrukerFile) = false
experimentIsCalibration(b::BrukerFile) = _iscalib(b.path)
experimentHasProcessing(b::BrukerFile) = experimentIsCalibration(b)
experimentHasReconstruction(b::BrukerFile) = false # fixme later
experimentHasMeasurement(b::BrukerFile) = true

# tracer parameters
tracerName(b::BrukerFile) = [b["PVM_MPI_Tracer"]]
tracerBatch(b::BrukerFile) = [b["PVM_MPI_TracerBatch"]]
tracerVolume(b::BrukerFile) = [parse(Float64,b["PVM_MPI_TracerVolume"])*1e-6]
tracerConcentration(b::BrukerFile) = [parse(Float64,b["PVM_MPI_TracerConcentration"])]
tracerSolute(b::BrukerFile) = ["Fe"]
function tracerInjectionTime(b::BrukerFile)
  initialFrames = b["MPI_InitialFrames"]
  if initialFrames == ""
    return [acqStartTime(b)]
  else
    return [acqStartTime(b) + Dates.Millisecond(
       round(Int64,parse(Int64, initialFrames)*dfCycle(b)*1000 ) )]
  end
end
tracerVendor(b::BrukerFile) = ["n.a."]

# scanner parameters
scannerFacility(b::BrukerFile) = latin1toutf8(b["ACQ_institution"])
scannerOperator(b::BrukerFile) = latin1toutf8(b["ACQ_operator"])
scannerManufacturer(b::BrukerFile) = "Bruker/Philips"
scannerName(b::BrukerFile) = b["ACQ_station"]
scannerTopology(b::BrukerFile) = "FFP"

# acquisition parameters
function acqStartTime(b::BrukerFile)
  if b["ACQ_time"]==""
    m = match(r"<(.+)\+","<0000-01-01T00:00:00,000+000>")
    timeString = replace(m.captures[1],"," => ".")
    return DateTime( timeString )
  else
    m = match(r"<(.+)\+",b["ACQ_time"])
    timeString = replace(m.captures[1],"," => ".")
    return DateTime( timeString )
  end
end
function acqNumFrames(b::BrukerFileMeas)
  M = Int64(b["ACQ_jobs"][1][8])
  return div(M,acqNumPeriodsPerFrame(b))
end

function acqNumFrames(b::BrukerFileCalib)
  M = parse(Int64,b["PVM_MPI_NrCalibrationScans"])
  A_ = b["PVM_MPI_NrBackgroundMeasurementCalibrationAdditionalScans"]
  A = (A_ == "") ? 0 : parse(Int64, A_)
  return div(M-A,acqNumPeriodsPerFrame(b))
end

function acqNumPatches(b::BrukerFile)
  if b.pretendToBeSinglePatch
    return 1
  else
    M = b["MPI_NSteps"]
    return (M == "") ? 1 : parse(Int64,M)
  end
end
function acqNumPeriodsPerFrame(b::BrukerFile)
  if b.pretendToBeSinglePatch
    return 1
  else
    M = b["MPI_RepetitionsPerStep"]
    N = acqNumPatches(b)
    return (M == "") ? N : N*parse(Int64,M)
  end
end

acqNumAverages(b::BrukerFileMeas) = parse(Int,b["NA"])
acqNumAverages(b::BrukerFileCalib) = parse(Int,b["NA"])*numSubPeriods(b)

function acqNumBGFrames(b::BrukerFile)
  n = b["PVM_MPI_NrBackgroundMeasurementCalibrationAllScans"]
  a = b["PVM_MPI_NrBackgroundMeasurementCalibrationAdditionalScans"]
  if n == ""
    n = "0"
  end
  if a == ""
    a = "0"
  end

  return parse(Int64,n)-parse(Int64,a)
end

function acqGradient(b::BrukerFile)
  G1::Float64 = parse(Float64,b["ACQ_MPI_selection_field_gradient"])
  G2 = Matrix(Diagonal([-0.5;-0.5;1.0])) .* G1

  G = zeros(3,3,1,acqNumPeriodsPerFrame(b))
  G[:,:,1,:] .= G2
  return G
end

function acqOffsetField(b::BrukerFile)
  if b["MPI_FocusFieldX"] != ""
    off = repeat(1e-3 * cat([-parse(Float64,a) for a in b["MPI_FocusFieldX"]],
                 [-parse(Float64,a) for a in b["MPI_FocusFieldY"]],
                 [parse(Float64,a) for a in b["MPI_FocusFieldZ"]],dims=2)',inner=(1,acqNumPeriodsPerPatch(b)))
  elseif b["CONFIG_MPI_FF_calibration"] != ""
    voltage = [parse(Float64,s) for s in b["ACQ_MPI_frame_list"]]
    voltage = reshape(voltage,4,:)
    voltage = repeat(voltage,inner=(1,acqNumPeriodsPerPatch(b)))
    calibFac = 1.0 / 100000 ./ parse.(Float64,b["CONFIG_MPI_FF_calibration"])
    off = [-1, -1, 1].*Float64[voltage[d,j]*calibFac[d-1] for d=2:4, j=1:acqNumPeriodsPerFrame(b)]
  else # legacy
    voltage = [parse(Float64,s) for s in b["ACQ_MPI_frame_list"]]
    voltage = reshape(voltage,4,:)
    voltage = repeat(voltage,inner=(1,acqNumPeriodsPerPatch(b)))
    calibFac = [2.5 / 49.45, 0.5 * (-2.5)*0.008/22.73, 0.5*2.5*0.008/22.73, -1.5*0.0094/13.2963]
    off = [-1, -1, 1].*Float64[voltage[d,j]*calibFac[d] for d=2:4, j=1:acqNumPeriodsPerFrame(b)]
  end
  return reshape(off, 3, 1, :)
end

# drive-field parameters
dfNumChannels(b::BrukerFile) = sum( selectedReceivers(b)[1:3] .== true )
   #sum( dfStrength(b)[1,:,1] .> 0) #TODO Not sure about this
function dfStrength(b::BrukerFile)
  str::Vector{String} = b["ACQ_MPI_drive_field_strength"]
  df = parse.(Float64,str) * 1e-3
  dfr = zeros(1,length(df),acqNumPeriodsPerFrame(b))
  dfr[1,:,:] .= df
  return dfr
end
dfPhase(b::BrukerFile) = dfStrength(b) .*0 .+  1.5707963267948966 # Bruker specific!
dfBaseFrequency(b::BrukerFile) = 2.5e6
dfCustomWaveform(b::BrukerFile) = nothing
dfDivider(b::BrukerFile) = reshape([102; 96; 99],:,1)
dfWaveform(b::BrukerFile) = "sine"
dfCycle(b::BrukerFile) = parse(Float64,b["PVM_MPI_DriveFieldCycle"]) / 1000 / numSubPeriods(b)

# receiver parameters
rxNumChannels(b::BrukerFile) = sum( selectedReceivers(b)[1:3] .== true )
rxBandwidth(b::BrukerFile) = parse(Float64,b["PVM_MPI_Bandwidth"])*1e6
rxNumSamplingPoints(b::BrukerFile) = div(parse(Int64,b["ACQ_size"][1]),numSubPeriods(b))
function rxTransferFunction(b::BrukerFile)
  filename = b["ACQ_comment"]
  if isfile(filename)
    return sampleTF(TransferFunction(filename), b)
  else
    return nothing
  end
end
function rxHasTransferFunction(b::BrukerFile)
  filename = b["ACQ_comment"]
  return isfile(filename)
end
rxInductionFactor(b::BrukerFile) = nothing
rxUnit(b::BrukerFile) = "a.u."
rxDataConversionFactor(b::BrukerFileMeas) =
                 repeat([1.0/acqNumAverages(b), 0.0], outer=(1,rxNumChannels(b)))
rxDataConversionFactor(b::BrukerFileCalib) =
                 repeat([1.0, 0.0], outer=(1,rxNumChannels(b)))

function rawDataLengthConsistent(b::BrukerFile)
  dataFilename = joinpath(b.path,"rawdata.job0")
  dType = acqNumAverages(b) == 1 ? Int16 : Int32

  # We derive numFrames from ACQ_jobs, since calibration files
  # are a bit longer than our acqNumFrames function reports
  # Bruker is padding the file for later processing
  numFrames = Int64(b["ACQ_jobs"][1][8])

  N = rxNumSamplingPoints(b)*numSubPeriods(b)*rxNumChannels(b)*
      acqNumPeriodsPerFrame(b)*numFrames*sizeof(dType)

  M = filesize(dataFilename)
  if N != M
    @warn "raw data length inconsistent N != M" N M
  end
  return N == M
end

function measData(b::BrukerFileMeas, frames=1:acqNumFrames(b), periods=1:acqNumPeriodsPerFrame(b),
                  receivers=1:rxNumChannels(b))

  dataFilename = joinpath(b.path,"rawdata.job0")
  dType = acqNumAverages(b) == 1 ? Int16 : Int32

  s = open(dataFilename)

  if numSubPeriods(b) == 1
    raw = Mmap.mmap(s, Array{dType,4},
             (rxNumSamplingPoints(b),rxNumChannels(b),acqNumPeriodsPerFrame(b),acqNumFrames(b)))
  else
    raw = Mmap.mmap(s, Array{dType,5},
             (rxNumSamplingPoints(b),numSubPeriods(b),rxNumChannels(b),acqNumPeriodsPerFrame(b),acqNumFrames(b)))
    raw = dropdims(sum(raw,dims=2),dims=2)
  end
  data = raw[:,receivers,periods,frames]
  close(s)

  return reshape(data, rxNumSamplingPoints(b), length(receivers),length(periods),length(frames))
end

function measData(b::BrukerFileCalib, frames=1:acqNumFrames(b), periods=1:acqNumPeriodsPerFrame(b),
                  receivers=1:rxNumChannels(b))

  sfFilename = joinpath(b.path,"pdata", "1", "systemMatrix")
  nFreq = div(rxNumSamplingPoints(b)*numSubPeriods(b),2)+1

  s = open(sfFilename)
  data = Mmap.mmap(s, Array{ComplexF64,4}, (prod(calibSize(b)),nFreq,rxNumChannels(b),1))
  #S = data[:,:,:,:]
  S = map(ComplexF32, data)
  close(s)
  rmul!(S,1.0/acqNumAverages(b))

  bgFilename = joinpath(b.path,"pdata", "1", "background")

  s = open(bgFilename)
  data = Mmap.mmap(s, Array{ComplexF64,4}, (acqNumBGFrames(b),nFreq,rxNumChannels(b),1))
  #bgdata = data[:,:,:,:]
  bgdata = map(ComplexF32, data)
  close(s)
  rmul!(bgdata,1.0/acqNumAverages(b))
  S_ = cat(S,bgdata,dims=1)
  if numSubPeriods(b) == 1
    return S_
  else
    return S_[:,1:numSubPeriods(b):end,:,:]
  end
end

function measDataTDPeriods(b::BrukerFile, periods=1:acqNumPeriods(b),
                  receivers=1:rxNumChannels(b))

  dataFilename = joinpath(b.path,"rawdata.job0")
  dType = acqNumAverages(b) == 1 ? Int16 : Int32

  s = open(dataFilename)
  if numSubPeriods(b) == 1
    raw = Mmap.mmap(s, Array{dType,3},
    (rxNumSamplingPoints(b),rxNumChannels(b),acqNumPeriods(b)))
  else
    raw = Mmap.mmap(s, Array{dType,4},
    (rxNumSamplingPoints(b),numSubPeriods(b),rxNumChannels(b),acqNumPeriods(b)))
    raw = dropdims(sum(raw,2),dims=2)
  end
  data = raw[:,receivers,periods]
  close(s)

  return reshape(data, rxNumSamplingPoints(b), length(receivers),length(periods))
end

systemMatrixWithBG(b::BrukerFileCalib) = measData(b)

# This is a special variant used for matrix compression
function systemMatrixWithBG(b::BrukerFileCalib, freq)
  sfFilename = joinpath(b.path,"pdata", "1", "systemMatrix")
  nFreq = div(rxNumSamplingPoints(b)*numSubPeriods(b),2)+1

  s = open(sfFilename)
  data = Mmap.mmap(s, Array{ComplexF64,4}, (prod(calibSize(b)),nFreq,rxNumChannels(b),1))
  #S = data[:,:,:,:]
  S = map(ComplexF32, data[:,freq,:,:])
  close(s)
  rmul!(S,1.0/acqNumAverages(b))

  bgFilename = joinpath(b.path,"pdata", "1", "background")

  s = open(bgFilename)
  data = Mmap.mmap(s, Array{ComplexF64,4}, (acqNumBGFrames(b),nFreq,rxNumChannels(b),1))
  #bgdata = data[:,:,:,:]
  bgdata = map(ComplexF32, data[:,freq,:,:])
  close(s)
  rmul!(bgdata,1.0/acqNumAverages(b))
  S_ = cat(S,bgdata,dims=1)
  if numSubPeriods(b) == 1
    return S_
  else
    return S_[:,1:numSubPeriods(b):end,:,:]
  end
end

function systemMatrix(b::BrukerFileCalib, rows, bgCorrection=true)

  localSFFilename = bgCorrection ? "systemMatrixBG" : "systemMatrix"
  sfFilename = joinpath(b.path,"pdata", "1", localSFFilename)
  nFreq = div(rxNumSamplingPoints(b)*numSubPeriods(b),2)+1

  if numSubPeriods(b) > 1
    rows_ = collect(rows)
    NFreq = rxNumFrequencies(b)
    NRx = rxNumChannels(b)
    stepsize = numSubPeriods(b)
    for k=1:length(rows)
      freq = mod1(rows[k],NFreq)
      rec = div(rows[k]-1,NFreq)
      rows_[k] = (freq-1)*stepsize+1 + rec*nFreq
    end
  else
    rows_ = rows
  end

  s = open(sfFilename)
  data = Mmap.mmap(s, Array{ComplexF64,2}, (prod(calibSize(b)),nFreq*rxNumChannels(b)))
  S = data[:,rows_]
  close(s)
  rmul!(S, 1.0/acqNumAverages(b))
  return S
end

measIsCalibProcessed(b::BrukerFile) = isfile(joinpath(b.path,"pdata", "1", "systemMatrix"))

measIsFourierTransformed(b::BrukerFileMeas) = false
measIsFourierTransformed(b::BrukerFileCalib) = measIsCalibProcessed(b)
measIsTFCorrected(b::BrukerFile) = false
measIsBGCorrected(b::BrukerFileMeas) = false
# We have it, but by default we pretend that it is not applied
measIsBGCorrected(b::BrukerFileCalib) = false

measIsFastFrameAxis(b::BrukerFileMeas) = false
measIsFastFrameAxis(b::BrukerFileCalib) = true

measIsFramePermutation(b::BrukerFileMeas) = false
measIsFramePermutation(b::BrukerFileCalib) = measIsCalibProcessed(b)

function measIsBGFrame(b::BrukerFileMeas)
  if !experimentIsCalibration(b)
    # If the file is not a calibration file we cannot say if any particular scans
    # were BG scans
    return zeros(Bool, acqNumFrames(b))
  else
    # In case of a calibration file we know the particular indices corresponding
    # to BG measurements
    isBG = zeros(Bool, acqNumFrames(b))
    increment = parse(Int,b["PVM_MPI_BackgroundMeasurementCalibrationIncrement"])+1
    isBG[1:increment:end] .= true

    return isBG
  end
end

# If the file is considered to be a calibration file, we will load
# the measurement in a processed form. In that case the BG measurements
# will be put at the end of the frame dimension.
function measIsBGFrame(b::BrukerFileCalib)
  if measIsCalibProcessed(b)
    return cat(zeros(Bool,acqNumFGFrames(b)),ones(Bool,acqNumBGFrames(b)),dims=1)
  else
    isBG = zeros(Bool, acqNumFrames(b))
    increment = parse(Int,b["PVM_MPI_BackgroundMeasurementCalibrationIncrement"])+1
    isBG[1:increment:end] .= true

    return isBG
  end
end

# measurements are not permuted
measFramePermutation(b::BrukerFileMeas) = nothing
# calibration scans are permuted
function measFramePermutation(b::BrukerFileCalib)
  if measIsCalibProcessed(b)
    # The following is a trick to obtain the permutation applied to the measurements
    # in a calibration measurement.
    bMeas = BrukerFile(b.path, isCalib=false)

    return fullFramePermutation(bMeas)
  else
    return nothing
  end
end

fullFramePermutation(f::BrukerFile) = fullFramePermutation(f, true)

measIsSpectralLeakageCorrected(b::BrukerFile) = get(b.params, "ACQ_MPI_spectral_cleaningl", "No") != "No"
measIsFrequencySelection(b::BrukerFile) = false
measIsSparsityTransformed(b::BrukerFile) = false

# calibrations
function calibSNR(b::BrukerFile)
  snrFilename = joinpath(b.path,"pdata", "1", "snr")

  if !isfile(snrFilename)
    return nothing
  end

  nFreq = div(rxNumSamplingPoints(b)*numSubPeriods(b),2)+1
  s = open(snrFilename)
  data = Mmap.mmap(s, Array{Float64,3}, (nFreq,rxNumChannels(b),1))
  snr = data[:,:,:]
  close(s)

  if numSubPeriods(b) == 1
    return snr
  else
    return snr[1:numSubPeriods(b):end,:,:]
  end
end
calibFov(b::BrukerFile) = [parse(Float64,s) for s = b["PVM_Fov"] ] * 1e-3
calibFovCenter(b::BrukerFile) =
          [parse(Float64,s) for s = b["PVM_MPI_FovCenter"] ] * 1e-3
calibSize(b::BrukerFile) = [parse(Int64,s) for s in b["PVM_Matrix"]]
calibOrder(b::BrukerFile) = "xyz"
calibPositions(b::BrukerFile) = nothing
calibOffsetField(b::BrukerFile) = nothing
calibDeltaSampleSize(b::BrukerFile) = [0.0, 0.0, 0.0] #TODO
calibMethod(b::BrukerFile) = "robot"
calibIsMeanderingGrid(b::BrukerFile) = true

# additional functions that should be implemented by an MPIFile
filepath(b::BrukerFile) = b.path

# special additional methods
function sfPath(b::BrukerFile)
  tmp = b["PVM_MPI_FilenameSystemMatrix",1]
  m = match(r"^(.+)\/pdata",tmp)
  return string(m.captures[1])
end

### The following is for field measurements from Alex Webers method
numCurrentSettings(b::BrukerFile) = parse(Int64,b["MPI_NrCurrentSettings"])
function currentSetting(b::BrukerFile)
  c = Float64[]
  for s in b["MPI_CurrentSetting"]
    append!(c,s)
  end
  return reshape(c,4,div(length(c),4))
end
ballRadius(b::BrukerFile) = parse(Float64,b["MPI_BallRadius"])
numLatitude(b::BrukerFile) = parse(Int64,b["MPI_NrLatitude"])
numMeridian(b::BrukerFile) = parse(Int64,b["MPI_NrMeridian"])

"""
This function is used as a workaround that no 1D and 2D sequences
were functional. Therefore we acquired 3D data and set the corresponding
drive-field amplitude to zero. Stepsize is then the number of period that
fit into the full 3D sequence.
"""
function numSubPeriods(f::BrukerFile)
  df = dfStrength(f)
  active_divider = copy(dfDivider(f))
  selected_channels = selectedChannels(f)

  for d=1:3
    active_divider[d] = (df[d] >= 0.0000001 && selected_channels[d]) ?
                         active_divider[d] : 1
  end
  floor(Int,(lcm(dfDivider(f)[selected_channels]) / lcm(active_divider)))
end

##### Reco

function recoData(f::BrukerFile)
  recoFilename = joinpath(f.path,"pdata", "1", "2dseq")
  N = recoSize(f)

  #if f["RECO_wordtype",1] != "_16BIT_SGN_INT"
  #  @error "Not yet implemented!"
  #end

  I = open(recoFilename,"r") do fd
    read!(fd,Array{Int16,3}(undef,1,prod(N),1))
  end
  return map(Float32,I)
end

recoResolution(f::BrukerFile) = push!(parse.(Float64,f["PVM_SpatResol"])./1000,
                                parse(Float64,f["ACQ_slice_thick"])./1000)

recoFov(f::BrukerFile) = recoResolution(f).*recoSize(f)

recoFovCenter(f::BrukerFile) = zeros(3)
recoSize(f::BrukerFile) = push!(parse.(Int,f["RECO_size",1]),
                                parse(Int,f["RecoObjectsPerRepetition",1]))
#recoOrder(f::BrukerFile) = f["/reconstruction/order"]
#recoPositions(f::BrukerFile) = f["/reconstruction/positions"]

###############################
# delta sample functions
###############################

export deltaSampleConcentration, deltaSampleVolume

function deltaSampleConcentration(b::BrukerFile)
 tmp = b["PVM_MPI_TracerConcentration"]
 if tmp != nothing
   return parse(Float64, tmp)
 else
   return 1.0
 end
end

deltaSampleConcentration(b::Array{T,1}) where {T<:BrukerFile} =
    map(deltaSampleConcentration, b)

function deltaSampleVolume(b::BrukerFile)
 V = parse(Float64, b["PVM_MPI_TracerVolume"] )*1e-6 # mu l
 return V
end


text:
import Pkg
Pkg.add("Debugger")
Pkg.add("Parameters")
Pkg.add("ResumableFunctions")

max_stars_repo_name: timm/misc

max_issues_repo_name: timm/misc

max_forks_repo_name: timm/misc


text:
# minhash.jl
# ==========
#
# Distance measures for MinHash sketches
#
# see DOI: 10.1186/s13059-016-0997-x
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/GeneticVariation.jl/blob/master/LICENSE.md

struct Jaccard end

"""
[MASH distances](http://doi.org/10.1186/s13059-016-0997-x), based on MinHash
sketches of genome sequences provide rapid genome-scale sequence comparisons
when sequence distance (not specific mutations) are all that's required.

A MinHash sketch is made by taking the `s` smallest hash values for kmers of
length `k` for a given sequence. The genome distance for two genomes is then
essentially the [Jaccard index](https://en.wikipedia.org/wiki/Jaccard_index)
of the minhashes, with some additional modification to account for the size of
the kmers used.
"""
struct MASH end

@inline function distance(::Type{Jaccard}, sketch1::MinHashSketch, sketch2::MinHashSketch)
    sketch1.kmersize == sketch2.kmersize || error("sketches must have same kmer length")
    length(sketch1) == length(sketch2) || error("sketches must be the same size")

    matches = 0
    sketchlen = length(sketch1)
    i = 1
    j = 1

    while i <= sketchlen && j <= sketchlen
        if sketch1.sketch[i] == sketch2.sketch[j]
            matches += 1
            i += 1
            j += 1
        elseif sketch1.sketch[i] < sketch2.sketch[j]
            while i <= sketchlen && sketch1.sketch[i] < sketch2.sketch[j]
                i += 1
            end
        elseif sketch2.sketch[j] < sketch1.sketch[i]
            while j <= sketchlen && sketch2.sketch[j] < sketch1.sketch[i]
                j += 1
            end
        end
    end

    if matches == sketchlen
        return 1.0
    else
        return matches / (2 * sketchlen - matches)
    end
end

@inline function distance(::Type{MASH}, sketch1::MinHashSketch, sketch2::MinHashSketch)
    j = distance(Jaccard, sketch1, sketch2)
    k = sketch1.kmersize
    return -1/k * log(2j / (1+j))
end

@inline mash(sketch1::MinHashSketch, sketch2::MinHashSketch) = distance(MASH, sketch1, sketch2)
@inline jaccard(sketch1::MinHashSketch, sketch2::MinHashSketch) = distance(Jaccard, sketch1, sketch2)


text:
#
# PostgreSQL access.
#

include("pgsql/libpq.jl")
include("pgsql/driver.jl")
include("pgsql/introspect.jl")
include("pgsql/attach.jl")

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
@testset "@fastkw" begin

end


text:
# # Instructions for Assignment {{{ASSIGNMENT}}}
#
# ## Introduction
#
# Some description of the Assignment

# ## Question 1
#
# Instructions

question1 = ""

# ## Question 2
#
# Instructions

question2 = ""

# ## Question 3
#
# Instructions

question3 = ""

# ## Question 4
#
# Instructions

question4 = ""

# ## Question 5
#
# Instructions

question5 = ""


text:
using LinearAlgebra

println("k_local(Elastic Modulus,thickness,height,length,number of nodes,row one, row two, column one, column two)")

function k_local(E,t,h,l,n,r_1,r_2,c_1,c_2)
    k_global = zeros(n,n)
    A = t*h
    k = ((E*A)/l)
    k_global[r_1,c_1] = k
    k_global[r_1,c_2] = -k
    k_global[r_2,c_1] = -k
    k_global[r_2,c_2] = k
    return k_global
end


text:
# This file creates the Manpower module, holding all the types, definitions, and
#   methods for the manpower planning project.

module ManpowerPlanning
    using SimJulia
    using ResumableFunctions
    using Plots
    using LightGraphs
    using MetaGraphs
    using GraphIO
    using EzXML
    using GraphPlot
    using Gadfly
    using Distributions
    using Polynomials
    using FileIO
    using SQLite
    using DataFrames
    using StatsBase
    using XLSX
    # using ExcelWrapper

    types = [
#        "historyEntry",
#        "history",
        "attrition",
        "personnelAttribute",
        "state",
        "compoundState",
        "condition",
        "transition",
#        "personnel",
#        "personnelDatabase",
#        "prerequisite",
#        "prerequisiteGroup",
        "recruitment",
        "retirement",
        "simulationReport",
        "manpowerSimulation"
    ]

    rootPath = Base.source_path()
    rootPath = rootPath isa Void ? "" : dirname( rootPath )
    funcPath = joinpath( rootPath, "Functions" )
    typePath = joinpath( rootPath, "Types" )

    # The types
    map( mpType -> include( joinpath( typePath, mpType * ".jl" ) ), types )

    # Some union type aliases
    include( joinpath( typePath, "typeAliases.jl" ) )

    # The functions
    map( mpType -> include( joinpath( funcPath, mpType * ".jl" ) ), types )
end  # module ManpowerPlanning