text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script to compile data for some descriptive statistics
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics, StatsBase, StatsPlots, Dates

# other scripts
dir_home = "x:/VIVES/1-Personal/Florian/git/2022_01_comext/src/"
include(dir_home * "functions.jl")

# location of data input/output (io)
dir_io = "C:/Users/u0148308/data/comext/"
dir_dropbox = "C:/Users/u0148308/Dropbox/BREXIT/"

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Raw data
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

df_VLAIO = CSV.read(dir_dropbox * "rawdata/" * "df_VLAIO" * ".csv", DataFrame)
transform!(df_VLAIO, ["PRODUCT_NC", "PRODUCT_NC_LAB1", "PRODUCT_NC_LAB2", "DECLARANT_ISO", "TRADE_TYPE", "PARTNER_ISO", "FLOW"] .=> ByRow(string), renamecols=false)
transform!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> convert(Union{Missing, Float64}, x)), renamecols=false)

# remove missings
df_VLAIO.PRODUCT_NC = lpad.(string.(df_VLAIO.PRODUCT_NC), 8, '0')
subset!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x))) # lose 1-16027346/21667553 ~26% observations

# MAD adjustment
outlier_cutoff = 3
transform!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)
gdf = groupby(df_VLAIO, "PRODUCT_NC")
df_VLAIO = transform(gdf, :UNIT_PRICE => MAD_method => :MAD)
subset!(df_VLAIO, :MAD => ByRow(x -> x < outlier_cutoff)) # lose further 1-13110297/16027346 ~18% observations (~40% of total)
df_VLAIO = df_VLAIO[:, Not([:UNIT_PRICE, :MAD])]

# add EU/WORLD again
EU27 = ["Roemenië", "Griekenland", "Oostenrijk", "Polen", "Duitsland", "Spanje", "Hongarije", "Slovakije", "Italië", "Nederland",
       "Frankrijk", "Letland", "Kroatië", "Cyprus", "Malta", "Litouwen", "Slovenië", "Estland", "Portugal", "Finland", "Tsjechië",
       "Luxemburg", "Zweden", "Denemarken", "Bulgarije", "Ierland", "België"]

# slightly modified functions
function append_WORLD(df::DataFrame)

    # sum over PARTNER
    cols_grouping = ["DECLARANT_ISO", "PRODUCT_NC", "PRODUCT_NC_LAB1", "PRODUCT_NC_LAB2", "TRADE_TYPE", "FLOW", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

    df_WORLD.TRADE_TYPE .= "total"
    df_WORLD.PARTNER_ISO .= "WORLD"

    df = vcat(df, df_WORLD)

    return df
end

function append_EU(df::DataFrame, EU::Vector{String})

    # subset and sum over EU ctrys
    df_EU = subset(df, :PARTNER_ISO => ByRow(x -> x in EU))
    cols_grouping = ["DECLARANT_ISO", "PRODUCT_NC", "PRODUCT_NC_LAB1", "PRODUCT_NC_LAB2", "TRADE_TYPE", "FLOW", "PERIOD"]
    gdf = groupby(df_EU, cols_grouping)
    df_EU = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

    df_EU.TRADE_TYPE .= "intra"
    df_EU.PARTNER_ISO .= "EU"

    df = vcat(df, df_EU)

    return df
end

df_VLAIO = append_WORLD(df_VLAIO)
df_VLAIO = append_EU(df_VLAIO, EU27)

# double check if no zeros/missing introduced by WORLD/EU
transform!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> ifelse(iszero(x) | ismissing(x), missing, x)), renamecols=false)
subset!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x))) # lose 1-15228494/15228494 ~0% observations

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

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD"]
gdf = groupby(df_VLAIO[:, [cols_grouping; "PRODUCT_NC"]], cols_grouping)

years = string.(2014:2021)
months = lpad.(1:12, 2, '0')
period = repeat(years, inner=12) .* repeat(months, outer=length(years))
period = parse.(Int64, period) # convert to Int64

intersection_exports = unique(df_VLAIO.PRODUCT_NC) # simply start with all products
intersection_imports = unique(df_VLAIO.PRODUCT_NC) # simply start with all products

for (i, time) in enumerate(period[1:end-2])

    # exports
    prod_GB_t1 = gdf[(DECLARANT_ISO = "Vlaanderen", PARTNER_ISO = "Verenigd Koninkrijk", FLOW = "exports", PERIOD = period[i])].PRODUCT_NC
    prod_GB_t2 = gdf[(DECLARANT_ISO = "Vlaanderen", PARTNER_ISO = "Verenigd Koninkrijk", FLOW = "exports", PERIOD = period[i+1])].PRODUCT_NC

    prod_NL_t1 = gdf[(DECLARANT_ISO = "Vlaanderen", PARTNER_ISO = "Nederland", FLOW = "exports", PERIOD = period[i])].PRODUCT_NC
    prod_NL_t2 = gdf[(DECLARANT_ISO = "Vlaanderen", PARTNER_ISO = "Nederland", FLOW = "exports", PERIOD = period[i+1])].PRODUCT_NC

    intersection_exports = intersect(intersection_exports, prod_GB_t1, prod_GB_t2, prod_NL_t1, prod_NL_t2)

    # imports
    prod_GB_t1 = gdf[(DECLARANT_ISO = "Vlaanderen", PARTNER_ISO = "Verenigd Koninkrijk", FLOW = "imports", PERIOD = period[i])].PRODUCT_NC
    prod_GB_t2 = gdf[(DECLARANT_ISO = "Vlaanderen", PARTNER_ISO = "Verenigd Koninkrijk", FLOW = "imports", PERIOD = period[i+1])].PRODUCT_NC

    prod_NL_t1 = gdf[(DECLARANT_ISO = "Vlaanderen", PARTNER_ISO = "Nederland", FLOW = "imports", PERIOD = period[i])].PRODUCT_NC
    prod_NL_t2 = gdf[(DECLARANT_ISO = "Vlaanderen", PARTNER_ISO = "Nederland", FLOW = "imports", PERIOD = period[i+1])].PRODUCT_NC

    intersection_imports = intersect(intersection_exports, prod_GB_t1, prod_GB_t2, prod_NL_t1, prod_NL_t2)

end

ext: jl


text:
include("4.5.1.1 (a+b sec)^n.jl")
include("4.5.1.2 (d sec)^n (a+b sec)^m.jl")
include("4.5.1.3 (d sin)^n (a+b sec)^m.jl")
include("4.5.1.4 (d tan)^n (a+b sec)^m.jl")
include("4.5.10 (c+d x)^m (a+b sec)^n.jl")
include("4.5.11 (e x)^m (a+b sec(c+d x^n))^p.jl")
include("4.5.2.1 (a+b sec)^m (c+d sec)^n.jl")
include("4.5.2.2 (g sec)^p (a+b sec)^m (c+d sec)^n.jl")
include("4.5.3.1 (a+b sec)^m (d sec)^n (A+B sec).jl")
include("4.5.4.1 (a+b sec)^m (A+B sec+C sec^2).jl")
include("4.5.4.2 (a+b sec)^m (d sec)^n (A+B sec+C sec^2).jl")
include("4.5.7 (d trig)^m (a+b (c sec)^n)^p.jl")
include("4.5.9 trig^m (a+b sec^n+c sec^(2 n))^p.jl")

max_stars_repo_path: src/rules/4 Trig functions/4.5 Secant/.jl


text:
@testset "Prints" begin
    @testset "routetable" begin
        io = IOBuffer()
        app = App()
        add_post!(app, "test") do req
            body = req.body
            return body
        end

        WebAPI.routetable(io, app)
        iostr = String(take!(io))

        str = """
          Method        Route
          ==============================================
          POST          test
        """

        @test iostr == str
    end

    @testset "App show" begin
        io = IOBuffer()
        app = App()
        add_post!(app, "test") do req
            body = req.body
            return body
        end

        show(io, app)
        iostr = String(take!(io))

        str = """
        App with reqparser = JSONParser and routes:\n
          Method        Route
          ==============================================
          POST          test
        """

        @test iostr == str
    end

    @testset "apprunning" begin
        # localhost
        io = IOBuffer()
        app = App()
        add_post!(app, "test") do req
            body = req.body
            return body
        end

        WebAPI.apprunning(io, app, Sockets.localhost, 80)
        iostr = String(take!(io))

        str = """

          App running at:
          http://localhost:80/\n
          Method        Route
          ==============================================
          POST          test
        """

        @test iostr == str

        # other ips
        io = IOBuffer()
        app = App()
        add_post!(app, "test") do req
            body = req.body
            return body
        end

        WebAPI.apprunning(io, app, ip"192.168.0.1", 80)
        iostr = String(take!(io))

        str = """

          App running at:
          http://192.168.0.1:80/\n
          Method        Route
          ==============================================
          POST          test
        """

        @test iostr == str
    end
end


text:
# Network inference algorithms and the InferredNetwork type. The algorithms MI, CLR, PUC and
# PIDC are explained in http://biorxiv.org/content/early/2017/04/26/082099 - along with terms
# such as specific information, proportional unique contribution, context, etc.

# Network inference algorithms
abstract type AbstractNetworkInference end
mutable struct CorrelationNetworkInference <: AbstractNetworkInference
    type::String
    signed::Bool
    values::Union{Array, Nothing}
end
struct MINetworkInference <: AbstractNetworkInference end
struct CLRNetworkInference <: AbstractNetworkInference end
struct PUCNetworkInference <: AbstractNetworkInference end
struct PIDCNetworkInference <: AbstractNetworkInference end

# correlation trait
get_correlation(::CorrelationNetworkInference) = true
get_correlation(::MINetworkInference) = false
get_correlation(::CLRNetworkInference) = false
get_correlation(::PUCNetworkInference) = false
get_correlation(::PIDCNetworkInference) = false

# Context trait
apply_context(::CorrelationNetworkInference) = false
apply_context(::MINetworkInference) = false
apply_context(::CLRNetworkInference) = true
apply_context(::PUCNetworkInference) = false
apply_context(::PIDCNetworkInference) = true

# PUC trait
get_puc(::CorrelationNetworkInference) = false
get_puc(::MINetworkInference) = false
get_puc(::CLRNetworkInference) = false
get_puc(::PUCNetworkInference) = true
get_puc(::PIDCNetworkInference) = true

# For sorting the edges
get_weight(edge::Edge) = edge.weight

# Gets the joint probability distribution for two Nodes.
function get_joint_probabilities(node1, node2, estimator)

    frequencies = get_frequencies_from_bin_ids(
        node1.binned_values,
        node2.binned_values,
        node1.number_of_bins,
        node2.number_of_bins
    )

    probabilities = get_probabilities(estimator, frequencies)
    # probabilities is already property of a node, but doing this gets correct array shapes.
    # Also, for MI and CLR, it means that we don't assume that the marginal probabilities for
    # a node are always the same, no matter what the second node is, meaning that we can use
    # estimators other than maximum likelihood. (We still can't do this for PUC and PIDC,
    # because we do make that assumption for 3-node joint distributions, in get_puc.)
    probabilities1 = sum(probabilities, dims = 2)
    probabilities2 = sum(probabilities, dims = 1)

    return (probabilities, probabilities1, probabilities2)

end

# get correlation between all pairs of Nodes
function get_correlation_scores(inference::AbstractNetworkInference, nodes, number_of_nodes)

    function get_correlation(i, j, correlation_scores)

        if inference.type == "Pearson"
            corr = cor(inference.values[i,:], inference.values[j,:])
        elseif inference.type == "Spearman"
            corr = corspearman(inference.values[i,:], inference.values[j,:])
        end

        if !inference.signed
            corr = abs(corr)
        end

        correlation_scores[i, j] = corr
        correlation_scores[j, i] = corr

    end

    correlation_scores = SharedArray{Float64}(number_of_nodes, number_of_nodes)

    @sync @distributed for i in 1:number_of_nodes
        for j in i+1:number_of_nodes
            get_correlation(i, j, correlation_scores)
        end
    end

    correlation_scores

end

# Gets the mutual information between all pairs of Nodes.
function get_mi_scores(nodes, number_of_nodes, estimator, base)

    function get_mi(node1, node2, i, j, base, mi_scores)
        probabilities, probabilities1, probabilities2 = get_joint_probabilities(node1, node2, estimator)
        mi = apply_mutual_information_formula(probabilities, probabilities1, probabilities2, base)
        mi_scores[i, j] = mi
        mi_scores[j, i] = mi
    end

    mi_scores = SharedArray{Float64}(number_of_nodes, number_of_nodes)

    @sync @distributed for i in 1 : number_of_nodes
        for j in i+1 : number_of_nodes
            get_mi(nodes[i], nodes[j], i, j, base, mi_scores)
        end
    end

    return mi_scores

end

# Gets the proportional unique contribution between all pairs of Nodes.
function get_puc_scores(nodes, number_of_nodes, estimator, base)

    function get_mi_and_si(node1, node2, base) # Mutual information and specific information
        probabilities, probabilities1, probabilities2 = get_joint_probabilities(node1, node2, estimator)
        mi = apply_mutual_information_formula(probabilities, probabilities1, probabilities2, base)
        si1 = apply_specific_information_formula(probabilities, probabilities1, probabilities2, 1, base)
        si2 = apply_specific_information_formula(probabilities, probabilities2, probabilities1, 2, base)
        return (mi, si1, si2)
    end

    function get_node_pairs(node1, node2, i, j, base)
        mi, si1, si2 = get_mi_and_si(node1, node2, base)
        node_pairs[i, j] = NodePair(mi, si1)
        node_pairs[j, i] = NodePair(mi, si2)
    end

    function increment_puc_scores(x, z, mi, redundancy, puc_scores)
        puc_score = (mi - redundancy) / mi
        puc_score = isfinite(puc_score) ? puc_score : zero(puc_score)
        puc_scores[x, z] += puc_score
        puc_scores[z, x] += puc_score
    end

    function get_puc(target, source1_target, source2_target, x, y, z, puc_scores)
        redundancy = apply_redundancy_formula(
            target.probabilities,
            source1_target.si,
            source2_target.si,
            base
        )
        increment_puc_scores(x, z, source1_target.mi, redundancy, puc_scores)
        increment_puc_scores(y, z, source2_target.mi, redundancy, puc_scores)
    end

    node_pairs = Array{NodePair}(undef, number_of_nodes, number_of_nodes)
    puc_scores = SharedArray{Float64}(number_of_nodes, number_of_nodes)

    for i in 1 : number_of_nodes
        for j in i+1 : number_of_nodes
            get_node_pairs(nodes[i], nodes[j], i, j, base)
        end
    end

    @sync @distributed for i in 1 : number_of_nodes
        for j in i+1 : number_of_nodes
            for k in j+1 : number_of_nodes
                get_puc(nodes[k], node_pairs[i, k], node_pairs[j, k], i, j, k, puc_scores)
                get_puc(nodes[j], node_pairs[i, j], node_pairs[k, j], i, k, j, puc_scores)
                get_puc(nodes[i], node_pairs[j, i], node_pairs[k, i], j, k, i, puc_scores)
            end
        end
    end

    return puc_scores

end

# Applies context to the raw edge weights.
function get_weights(inference, scores, number_of_nodes, nodes)

    # In their respective original implementations, CLR and PIDC applied network context in slightly
    # different ways. Those differences are respected here; in informal tests, they have not been
    # found to make much of a difference.
    function get_weight(::PIDCNetworkInference, i, j, scores, weights, nodes)
        score = scores[i, j]
        scores_i = vcat(scores[1:i-1, i], scores[i+1:end, i])
        scores_j = vcat(scores[1:j-1, j], scores[j+1:end, j])
        try
            weights[i, j] = cdf(fit(Gamma, scores_i), score) + cdf(fit(Gamma, scores_j), score)
        catch
            println(string("Gamma distribution failed for ", nodes[i].label, " and ", nodes[j].label, "; used normal instead."))
            apply_clr_context(i, j, score, scores_i, scores_j, weights)
        end
    end

    function get_weight(::CLRNetworkInference, i, j, scores, weights, nodes)
        score = scores[i, j]
        scores_i = vcat(scores[1:i-1, i], scores[i+1:end, i])
        scores_j = vcat(scores[1:j-1, j], scores[j+1:end, j])
        apply_clr_context(i, j, score, scores_i, scores_j, weights)
    end

    function apply_clr_context(i, j, score, scores_i, scores_j, weights)
        difference_i = score - mean(scores_i)
        difference_j = score - mean(scores_j)
        weights[i, j] = sqrt(
            (var(scores_i) == 0 || difference_i < 0 ? 0 : difference_i^2 / var(scores_i)) +
            (var(scores_j) == 0 || difference_j < 0 ? 0 : difference_j^2 / var(scores_j))
        )
    end

    weights = SharedArray{Float64}(number_of_nodes, number_of_nodes)

    @sync @distributed for i in 1 : number_of_nodes
        for j in i+1 : number_of_nodes
            get_weight(inference, i, j, scores, weights, nodes)
        end
    end

    return weights

end

"""
InferredNetwork type. Represents a weighted, fully connected network, where an
edges's weight indicates the relative confidence of that edge existing in the true
network.

Fields:
* `nodes`: array of all the nodes, in an arbitrary order
* `edges`: array of all the edges, in descending order of weight
"""
struct InferredNetwork
    nodes::Array{Node}
    edges::Array{Edge}
end

# Constructs an InferredNetwork given a network inference algorithm and an array of
# Nodes.
#
# Keyword arguments:
# - estimator: algorithm for estimating the probability distribution
# (The "maximum_likelihood" estimator is recommended for PUC and PIDC, because speedups
# are made here, based on the assumption that the marginal probability distribution for
# a node, from the joint distribution with any two other nodes is always the same. If
# the joint distributions are estimated using other estimators, this assumption is
# violated for PUC and PIDC in get_puc and get_joint_probabilities.)
# - base: base for the information measures
function InferredNetwork(inference::AbstractNetworkInference, nodes::Array{Node}; estimator = "maximum_likelihood", base = 2)

    # Constants and containers
    number_of_nodes = length(nodes)
    edges = Array{Edge}(undef, binomial(number_of_nodes, 2))

    # Get the raw scores
    if get_correlation(inference)
        scores = get_correlation_scores(inference, nodes, number_of_nodes)
    elseif get_puc(inference)
        scores = get_puc_scores(nodes, number_of_nodes, estimator, base)
    else
        scores = get_mi_scores(nodes, number_of_nodes, estimator, base)
    end

    # Apply context if necessary
    if apply_context(inference)
        weights = get_weights(inference, scores, number_of_nodes, nodes)
    else
        weights = scores
    end

    # Get edges from scores
    index = 0
    for i in 1 : number_of_nodes
        node1 = nodes[i]
        for j in i+1 : number_of_nodes
            index += 1
            node2 = nodes[j]
            edges[index] = Edge(
                [node1, node2],
                weights[i, j]
            )
        end
    end
    sort!(edges; by = get_weight, rev = true)

    return InferredNetwork(nodes, edges)

end


text:
"""
A Robust Online Sequential Extreme Learning Machine.
Lecture Notes in Computer Science.

Minh-Tuelm T Hoelmg, Hieu T Huynh, Nguyen H Vo, elmd Yonggwelm Won.

Advelmces in Neural Networks – ISNN 2007, Chapter 126, 1077-1086. Berlin,
Heidelberg.

http://link.springer.com/10.1007/978-3-540-72383-7_126

"""
mutable struct ROSELM{TN} <: AbstractSLFN where TN <: Linear
    p::Int  # Number of training points
    q::Int  # Dimensionality of function domain
    s::Int  # number of neurons
    c::Float64
    maxit::Int
    μx::Vector{Float64}
    σx::Vector{Float64}
    Wt::Matrix{Float64}  # trelmspose of W matrix
    d::Vector{Float64}
    v::Vector{Float64}
    neuron_type::TN

    # internal state
    M::Matrix{Float64}
    p_tot::Int

    function ROSELM(p::Int, q::Int, s::Int, c::Float64, maxit::Int,
                    neuron_type::TN, μx, σx) where TN
        Wt = Array{Float64}(undef, q, s)
        d = Array{Float64}(undef, s)
        v = Array{Float64}(undef, s)
        M = zeros(Float64, 0, 0)
        new{TN}(p, q, s, c, maxit, μx, σx, Wt, d, v, neuron_type, M, 0)
    end
end

function ROSELM(x::AbstractArray, y::AbstractVector;
                neuron_type::TN=Linear(Tanh()),
                s::Int=size(x, 1), maxit::Int=1000,
                c::Float64=2.5) where TN <: Linear
    p = size(x, 1)
    q = size(x, 2)

    @assert size(y, 1) == p "x and y must have same number of observations"

    s = min(s, p)
    xn, μx, σx = standardize(x[:, :])
    out = ROSELM(p, q, s, c, maxit, neuron_type, μx, σx)
    fit!(out, xn, y)
    out
end

## API methods
function fit!(elm::ROSELM, x::AbstractArray, y::AbstractVector)
    if elm.p_tot == 0
        local S
        for i in 1:elm.maxit
            # initialization phase -- try to find invertible S matrix
            randn!(elm.Wt)
            elm.Wt ./= elm.c
            elm.d = -diag(x * elm.Wt)
            S = hidden_out(elm, x)

            if !(rank(S) < elm.s)
                break
            end
        end

        # elm.v = S \ y
        elm.M = pinv(S'S)
        elm.v = elm.M*S'*y
    else
        # post-init phase
        S = hidden_out(elm, x)
        elm.M -= elm.M * S' * inv(I + S*elm.M*S') * S * elm.M
        elm.v += elm.M * S' * (y - S*elm.v)
    end

    elm.p_tot += size(x, 1)
    elm

end

function Base.show(io::IO, elm::ROSELM{TA}) where TA
    s =
    """
    ROSELM with
      - $(TA) Activation function
      - $(elm.q) input dimension(s)
      - $(elm.s) neuron(s)
      - $(elm.p) training point(s)
    """
    print(io, s)
end

function test_me(::Type{ROSELM})
    x = range(0, stop=1, length=20) + 0.1*randn(20)
    y = sin.(6*x[:])
    roselm = ROSELM(x, y)

    xt = range(0, stop=1, length=300)
    yt = sin.(6*xt)

    @show maximum(abs.(roselm(xt) - yt))

    for i in 1:50
        x2 = rand(15) + 0.1 * randn(15)
        y2 = sin.(6*x2)
        fit!(roselm, x2, y2)
        @printf "%i\t%2.4e\n" i maximum(abs.(roselm(xt) - yt))
    end
end


text:
# Bukdu v0.4.4
using Bukdu # ApplicationController Conn render HTML routes get Router
using Documenter.Utilities.DOM: @tags

struct Calculator <: ApplicationController
    conn::Conn
end

function index(c::Calculator)
    @tags h3 div a

    x = get(c.params, :x, 0)
    y = get(c.params, :y, 0)

    rand_x = rand(1:10)
    rand_y = rand(1:10)

    render(HTML, div(
        h3("Calculator"),
        div(
            string(x, +, y, '=', x+y),
        ),
        div(
            a[:href => "/?x=$rand_x&y=$rand_y"]("/?x=$rand_x&y=$rand_y")
        )
    ))
end

if PROGRAM_FILE == basename(@__FILE__)

routes() do
    get("/", Calculator, index)
end

Bukdu.start(8080)

Router.call(get, "/") #
# CLI.routes()

Base.JLOptions().isinteractive==0 && wait()

# Bukdu.stop()

end # if


text:
require("test.jl")

using Scallop

my_tests = ["test/backtest.jl",
            "test/book.jl",
            "test/chart.jl",
            "test/metrics.jl",
            "test/underlying.jl",
            "test/walkforward.jl"]

print_with_color(:cyan, "Running tests: ")
println("")

for my_test in my_tests
    print_with_color(:magenta, "**   ")
    print_with_color(:blue, "$my_test")
    println("")
    include(my_test)
end


text:
# Simple rewrite of Conv without biases, usefull for convolutional layer followed by BatchNorm
using Flux: @functor, glorot_uniform, expand

struct ConvNB{N,M,F,A}
  σ::F
  weight::A
  stride::NTuple{N,Int}
  pad::NTuple{M,Int}
  dilation::NTuple{N,Int}
end

function ConvNB(w::AbstractArray{T,N}, σ = identity; stride = 1, pad = 0,
                    dilation = 1) where {T,N}
  stride = expand(Val(N-2), stride)
  pad = expand(Val(2*(N-2)), pad)
  dilation = expand(Val(N-2), dilation)
  return ConvNB(σ, w, stride, pad, dilation)
end

function ConvNB(k::NTuple{N,Integer}, ch::Pair{<:Integer,<:Integer}, σ = identity;
              init = glorot_uniform,  stride = 1, pad = 0, dilation = 1) where N
  ConvNB(init(k..., ch...), σ, stride = stride, pad = pad, dilation = dilation)
end

@functor ConvNB

function (c::ConvNB)(x::AbstractArray)
  # TODO: breaks gpu broadcast :(
  # ndims(x) == ndims(c.weight)-1 && return squeezebatch(c(reshape(x, size(x)..., 1)))
  σ = c.σ
  cdims = DenseConvDims(x, c.weight; stride=c.stride, padding=c.pad, dilation=c.dilation)
  σ.(conv(x, c.weight, cdims))
end

function Base.show(io::IO, l::ConvNB)
  print(io, "ConvNB(", size(l.weight)[1:ndims(l.weight)-2])
  print(io, ", ", size(l.weight, ndims(l.weight)-1), "=>", size(l.weight, ndims(l.weight)))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end

(a::ConvNB{<:Any,<:Any,W})(x::AbstractArray{T}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
  invoke(a, Tuple{AbstractArray}, x)

(a::ConvNB{<:Any,<:Any,W})(x::AbstractArray{<:Real}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
  a(T.(x))


text:
#####
##### Exceptions
#####

struct SlackError <: Exception
    error::String
end

#####
##### Attachments
#####

# Upload dispatch strategy:
# We need a local file to upload.
#
# entrypoint is `local_file(item)`.
# Pairs are destructured to to the two argument `local_file(name, object)` methods.
# Non-pairs are assumed to be file paths and are dispached to `local_file(path)`.
local_file(item::Pair; kw...) = local_file(first(item), last(item); kw...)
local_file(file::AbstractString; kw...) = file

# If the `path` is not a `AbstractString`, we assume it isn't a local path
# (maybe it's an S3Path, etc.). So we conservatively write a new local file to upload.
# Thus, we only generically require `basename` and `read` to be supported for path types.
local_file(file; kw...) = local_file(basename(file), read(file); kw...)

# Two-argument `local_file`. Second argument is assumed to be an object to write, not a filepath
# (since that is only allowed in the 1-arg case). Bytes and strings are written to a file and everything else uses `FileIO.save`.

function local_file(name, object; dir=mktempdir())
    local_path = joinpath(dir, name)
    if object isa Union{Vector{UInt8},<:AbstractString}
        write(local_path, object)
    else
        save(local_path, object)
    end
    return local_path
end

function upload_file(local_path::AbstractString; extra_args=String[])
    api = "https://slack.com/api/files.upload"

    token = get(ENV, "SLACK_TOKEN", nothing)
    if token === nothing
        @warn "No Slack token provided; file not sent." api local_path
        return nothing
    else
        @debug "Uploading slack file" api local_path
    end

    auth = "Authorization: Bearer $(token)"

    response = @maybecatch begin
        JSON3.read(@mock readchomp(`curl -s -F file=@$(local_path) $(extra_args) -H $auth $api`))
    end "Error when attempting to upload file to Slack"

    response === nothing && return nothing

    if haskey(response, :ok) === true && response[:ok] != true
        return @maybecatch begin
            err = haskey(response, :error) ? string(response.error) :
                  "No error field returned"
            throw(SlackError(err))
        end "Error reported by Slack API"
    end

    @debug "Slack responded" response
    return response
end

#####
##### Messages
#####

function send_message(thread::SlackThread, text::AbstractString; options...)
    data = Dict{String,Any}(string(k) => v for (k, v) in pairs(options))
    data["channel"] = thread.channel
    data["text"] = text
    if thread.ts !== nothing
        data["thread_ts"] = thread.ts
    end
    data_str = JSON3.write(data)
    api = "https://slack.com/api/chat.postMessage"

    token = get(ENV, "SLACK_TOKEN", nothing)
    if token === nothing
        @warn "No Slack token provided; message not sent." data api
        return nothing
    elseif thread.channel === nothing
        @warn "No Slack channel configured; message not sent." data api
        return nothing
    else
        @debug "Sending slack message" data api
    end
    auth = "Authorization: Bearer $(token)"

    response = @maybecatch begin
        JSON3.read(@mock readchomp(`curl -s -X POST -H $auth -H 'Content-type: application/json; charset=utf-8' --data $(data_str) $api`))
    end "Error when attempting to send message to Slack thread"

    response === nothing && return nothing
    @debug "Slack responded" response

    if haskey(response, :ok) === true && response[:ok] === false
        return @maybecatch begin
            err = haskey(response, :error) ? string(response.error) :
                  "No error field returned"
            throw(SlackError(err))
        end "Error reported by Slack API"
    end

    if thread.ts === nothing && hasproperty(response, :ts) === true
        thread.ts = response.ts
    end
    return response
end


text:
export quick_sort!

function quick_sort!(arr::Vector, h::Int64, t::Int64)
    len = t-h+1
    if len < 2
        return arr
    end

    n, i = h, h
    pivot = arr[t]
    for n = h:t
        if arr[n] <= pivot
            arr[i], arr[n] = arr[n], arr[i]
            i += 1
        end
    end

    if i > t
        i = t
    end

    quick_sort!(arr, h, i-1)
    quick_sort!(arr, i, t)

    return arr
end


text:
_err_str(err) = sprint(showerror, err, catch_backtrace())

"""
    print the err text
"""
_printerr(io::IO, err) = print(io, _err_str(err))
_printerr(err) = _printerr(stdout, err)

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
"""

```
closest(x,y)
```

Computes the minimum distance between two sets of atoms and returns the indexes of the atoms
and their distance. Both vector of atoms or vectors of coordinates can be used as input.

### Examples

```julia-repl
julia> model = wget("1BSX");

julia> protein = select(model,"protein");

julia> ligand = select(model,"resname T3");

julia> closest(ligand,protein)
(43, 3684, 2.7775834820937417)

julia> ligand[43]
    4037   O1      T3     B        2      512  -22.568   81.625    3.159 36.59  1.00     1       -      4041

julia> closest(ligand[43],protein)
(1, 3684, 2.7775834820937417)

julia> x = coor(protein)
3994-element Vector{SVector{3, Float64}}:
 [52.884, 24.022, 35.587]
 [52.916, 24.598, 36.993]
 ⋮
 [-46.887, 86.925, 13.235]
 [-47.164, 83.593, 15.25]

julia> closest(ligand,x)
(43, 3684, 2.7775834820937417)

```

"""
function closest(x::AbstractVector{T1}, y::AbstractVector{T2}) where {T1,T2 <: Union{SVector, Atom}}
    imin = -1
    jmin = -1
    dmin = +Inf
    for (i,xatom) in pairs(x)
        for (j,yatom) in pairs(y)
            d = distance(xatom, yatom)
            if d < dmin
                imin = i
                jmin = j
                dmin = d
            end
        end
    end
    return imin, jmin, dmin
end
closest(x::Atom, y::AbstractVector{<:Union{SVector, Atom}}) = closest(SVector(x),y)
closest(x::AbstractVector{<:Union{SVector, Atom}}, y::Atom) = closest(x,SVector(y))
closest(x::AbstractVector{<:Union{SVector, Atom}}, y::SVector{3,<:Real}) = closest(x,SVector{1,typeof(y)}(y))
closest(x::SVector{3,<:Real},y::AbstractVector{<:Union{SVector, Atom}}) = closest(SVector{1,typeof(x)}(x),y)
closest(x::Atom, y::Atom) = closest(SVector(x),SVector(y))
closest(x::SVector{3,<:Real}, y::SVector{3,<:Real}) = closest(SVector{1,typeof(x)}(x),SVector{1,typeof(y)}(y))
closest(x::Atom, y::SVector{3,<:Real}) = closest(SVector(x),SVector{1,typeof(y)}(y))
closest(x::SVector{3,<:Real}, y::Atom) = closest(SVector{1,typeof(x)}(x),SVector(y))
closest(x::Residue,y::Residue) = closest(x.atoms[x.range],y.atoms[y.range])

"""

```
distance(x,y)
```

Computes the minimum distance between two sets of atoms, between an atom and a set of atoms, or simply
the distance between two atoms. The input may be a vector of `Atom`s, or the
coordinates that are output of the `coor` function.

### Examples

```julia-repl
julia> model = wget("1BSX");

julia> protein = select(model,"protein");

julia> ligand = select(model,"resname T3");

julia> distance(protein,ligand)
2.7775834820937417

julia> distance(protein[1],ligand[3])
36.453551075306784

julia> distance(coor(ligand),protein)
2.7775834820937417

```

"""
distance(x::SVector, y::SVector) = norm(x - y)
distance(x::Atom, y::Atom) = norm(coor(x) - coor(y))
distance(x::Atom, y::SVector) = norm(coor(x) - y)
distance(x::SVector, y::Atom) = norm(x - coor(y))
distance(x,y) = closest(x,y)[3]

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
using FirstOrderLogic
using Test

@testset "parser & types " begin
  # test parser on simple first order formulae
  @test isa(fol"f(R)", PredicateTerm)
  @test isa(fol"~f(R)", NegationTerm)
  @test isa(fol"a(X) & b(Y)", AndTerm)
  @test isa(fol"a(X)|b(Y)", OrTerm)
  @test isa(fol"![X,Y] : a(f(X),Y)", AQuantifierTerm)
  @test isa(fol"?[X,Y,Z] : (a(f(X)) & b(c(Y,Z)))", EQuantifierTerm)
  @test isa(fol"d(Z) & (a(f(X)) & b(c(Y,Z)))", AndTerm)

  # test parser on annotated first order formulae
  @test isa(tptp"fof(union,axiom,
                 ! [X,A,B] :
                 ( member(X,union(A,B))
                 <=> ( member(X,A)
                 | member(X,B) ) ))."[2], AQuantifierTerm)
  @test length(parseTPTP("fofTest.p")) == 7

  # test parser on FOOL formula (with types and Boolean variables reifying formulas)
  @test isa(fol"! [A: animal] : ? [H: human] : H = owner_of(A) ", AQuantifierTerm)
  @test isa(fol" ! [X: $o]: (X | ~X)", AQuantifierTerm)
  @test length(parseTPTP("tffTest.p")) == 4

  # validate the environment created by tffTypeTest.p
  tmp = parseTPTP("tffTypeTest.p")
  en = FirstOrderLogic.e
  @test haskey(en,"A") && haskey(en,"C") && haskey(en,"D") && haskey(en,"H") &&
    haskey(en,"dog") && haskey(en,"cat") && haskey(en,"human") && haskey(en,"animal") &&
    haskey(en,"garfield") && haskey(en,"odie") && haskey(en,"jon") &&
    haskey(en,"chased") && haskey(en,"hates") && haskey(en,"owner_of") &&
    haskey(en,"cat_to_animal") && haskey(en,"dog_to_animal")

  # test the ingested SententialTerms for their type
  wellTyped = Vector{Bool}(undef, 0)
  for s in filter(x->isa(x[2],FirstOrderLogic.SententialTerm), tmp)
    push!(wellTyped, isWellTyped(s[2],FirstOrderLogic.e))
  end
  @test all(wellTyped)
end # parser tests

@testset "transformations" begin
  # conversion of FOOL to FOL

  fool = tptp"tff(imply_defn,axiom, ! [X: $o,Y: $o]: (imply(X,Y) <=> (~X | Y)) )."[end]
  fol = simplify(FOOL2FOL!(fool, FirstOrderLogic.e))
  @test isa(fol, AQuantifierTerm{AndTerm{OrTerm{SententialTerm}}})

  fool = parseTPTP("tffTest.p")[end][end]
  fol = simplify(FOOL2FOL!(fool[end][end], FirstOrderLogic.e))
  @test isa(fol, AndTerm{AQuantifierTerm})

  #skolemization
  f = fol"?[Z]:(?[X]: q(X, Z) | ?[X]: p(X)) => ~(?[X] : p(X) & ![X] : ?[Z]: q(Z, X))"

end # transformation tests


text:
module PortNumbers

import Base: string, convert

export PortClass, EstablishedPort, RootPort, UserPort, AnyPort, PortNumber
export portmin, portmax
export PortNumber
export string, convert
export make_port_number

const ESTABLISHED_MIN_DEFAULT = 1
const ESTABLISHED_MAX_DEFAULT = 256
const ROOT_MIN_DEFAULT = 1
const ROOT_MAX_DEFAULT = 1023
const USER_MIN_DEFAULT = 1024
const USER_MAX_DEFAULT = 65535
const ANY_MIN_DEFAULT = 1
const ANY_MAX_DEFAULT = 65535

const ESTABLISHED_MIN_DEFAULT = 1
const ESTABLISHED_MAX_DEFAULT = 256
const ROOT_MIN_DEFAULT = 1
const ROOT_MAX_DEFAULT = 1023
const USER_MIN_DEFAULT = 1024
const USER_MAX_DEFAULT = 65535
const ANY_MIM_DEFAULT = 1
const ANY_MAX_DEFAULT = 65535

est_min = ESTABLISHED_MIN_DEFAULT
est_max = ESTABLISHED_MAX_DEFAULT

root_min = ROOT_MIN_DEFAULT
root_max = ROOT_MAX_DEFAULT

user_min = USER_MIN_DEFAULT
user_max = USER_MAX_DEFAULT

any_min = ANY_MIN_DEFAULT
any_max = ANY_MAX_DEFAULT

"""
### PortClass

An abstract type used as the super type
 for the various port number classes.
"""
abstract PortClass

abstract AnyPort <: PortClass
minport( ::Type{AnyPort} ) = any_min
maxport( ::Type{AnyPort} ) = any_max

abstract EstablishedPort <: PortClass
minport( ::Type{EstablishedPort} ) = est_min
maxport( ::Type{EstablishedPort} ) = est_max

abstract RootPort <: PortClass
minport( ::Type{RootPort} ) = root_min
maxport( ::Type{RootPort} ) = root_max

abstract UserPort <: PortClass
minport( ::Type{UserPort} ) = user_min
maxport( ::Type{UserPort} ) = user_max

immutable PortNumber{T<:PortClass}
    value::Integer
    function PortNumber( x::Integer )
        @assert x >= minport( T )
        @assert x <= maxport( T )
        new( x )
    end
end

convert{T<:Integer}( ::Type{T}, x::PortNumber ) = T( x.value )
string( x::PortNumber ) = Int( x ) |> string
convert{T <: AbstractString }( ::Type{T}, x::PortNumber ) = string( x ) |> T

"""
### `make_port_number( x::Integer )`

Return a port number object.

#### Example

`*julia> make_port_number( 22 )`

  --> `PortNumber{EstablishedPort}( 22 )`

"""
function make_port_number( x::Integer )
    @assert x >= minport( AnyPort )
    @assert x <= maxport( AnyPort )
    if x <= maxport( EstablishedPort )
        return PortNumber{EstablishedPort}( x )
    elseif x <= maxport( RootPort )
        return PortNumber{RootPort}( x )
    else
        return PortNumber{UserPort}( x )
    end
end

end # module

max_stars_repo_name: sabjohnso/PortNumbers

max_issues_repo_name: sabjohnso/PortNumbers

max_forks_repo_name: sabjohnso/PortNumbers


text:
using SparseArrays
using NPZ

"""
    tosparse(obj) -> SparseMatrixCSC

Returns the sparse matrix corresponding to the PyObject obj.
obj can be a scipy.sparse matrix or a qutip Qobj
"""
function tosparse(obj::PyObject)
    # This is just a wrapper to the ugly py"..." syntax
    (I, J, V, m, n) = pystuff.sparse_to_ijv(obj)
    return sparse(Int64.(I), Int64.(J), V, Int64(m), Int64(n))
end

"""
"""
function sparse_fromfile(filename)
    data = npzread(filename)
    return sparse(Int64.(data["arr_0"]), Int64.(data["arr_1"]), data["arr_2"], Int64(data["arr_3"]), Int64(data["arr_4"]))
end

"""
    liouvillian(N) -> SparseMatrixCSC

Return a dephasing Liouvillian for N spins in the Dicke basis
"""
function liouvillian(N::T) where T<:Integer
    sys = piqs.Dicke(N)
    sys.dephasing = 4.

    return tosparse(sys.liouvillian())
end

"""
    css(N) -> SparseMatrixCSC

Return a coherent spin state in the Dicke basis for N spins
"""
function css(N::T) where T<:Integer
    return tosparse(piqs.css(N))
end

"""
    jspin(N) -> Tuple{SparseMatrixCSC}

Return the tuple (Jx, Jy, Jz) of the spin operators in the
Dicke basis. The returned operators are sparse arrays.
"""
function jspin(N::T) where T<:Integer
    return tosparse.(piqs.jspin(N))
end

"""
    j_min(N)

Return the minimum value of j for N spins
"""
function j_min(N::T) where T<:Integer
    if N % 2 == 0
        return 0
    else
        return 0.5
    end
end

"""
    j_vals(N)

Returns a list of values of j in decreasing order
"""
function j_vals(N::T) where T<:Integer
    j = (N/2):-1:j_min(N)
    return j
end

"""
    block_sizes(N)

Return a list of block sizes for the density matrix of a N-spin
system in the Dicke basis.
"""
function block_sizes(N::T) where T<:Integer
   return Int.(2*j_vals(N) .+ 1)
end


text:
abstract ImplicitStack
immutable ImplicitIndexStack{S} <: ImplicitStack
    stack::Vector{S}
end
Base.copy(s::ImplicitIndexStack) = typeof(s)(copy(s.stack))
immutable ImplicitNodeStack{T, S} <: ImplicitStack
    node_stack::Vector{T}
    idx_stack::ImplicitIndexStack{S}
end
Base.copy(s::ImplicitNodeStack) = typeof(s)(copy(s.node_stack), copy(s.idx_stack))
Base.isempty(s::ImplicitNodeStack) = isempty(s.node_stack)
Base.isempty(s::ImplicitIndexStack) = isempty(s.stack)
immutable ImplicitChildStates{T, S}
    tree::T
    stack::S
end
Base.iteratorsize{T<:ImplicitChildStates}(::Type{T}) = Base.SizeUnknown()
children(states::ImplicitChildStates) = children(states.tree, states.stack)
children(stack::ImplicitNodeStack) = children(stack.node_stack[end])

parentstate(tree, state::ImplicitNodeStack) =
    ImplicitNodeStack(state.node_stack[1:end-1], parentstate(tree, state.idx_stack))
parentstate(tree, state::ImplicitIndexStack) = ImplicitIndexStack(state.stack[1:end-1])
childstates(tree, state) = childstates(tree, state, treekind(tree))
childstates(tree, state::ImplicitStack) = ImplicitChildStates(tree, state)
children(tree, state::ImplicitNodeStack) = isempty(state) ?
    children(tree) : children(tree, state.node_stack[end])
children(tree, state::ImplicitIndexStack) = isempty(state) ?
    children(tree) : children(tree, tree[state.stack[end]])

childstates(s::ImplicitChildStates) = isempty(s.stack) ?
  childstates(s.tree, s.tree) : childstates(s.tree,
    isa(s.stack, ImplicitNodeStack) ? s.stack.node_stack[end] : s.stack.stack[end])
nextind(s::ImplicitChildStates, ind) = nextind(childstates(s), ind)
start(s::ImplicitChildStates) = start(childstates(s))
function next(s::ImplicitChildStates, ind)
    ni = next(childstates(s), ind)[2]
    (update_state!(copy(s.stack), childstates(s), ind, treekind(s.tree)), ni)
end
done(s::ImplicitChildStates, ind) = done(childstates(s), ind)

update_state!(ns, cs, idx, _) = update_state!(ns, cs, idx)
update_state!(ns::ImplicitIndexStack, cs, idx) = (push!(ns.stack, idx); ns)
update_state!(ns::ImplicitIndexStack, cs, idx, ::IndexedTree) = (push!(ns.stack, next(cs, idx)[1]); ns)
function update_state!(ns::ImplicitNodeStack, cs, idx)
    push!(ns.node_stack, next(cs, idx)[1])
    update_state!(ns.idx_stack, cs, idx)
    ns
end

joinstate(state::ImplicitNodeStack, new_state::ImplicitNodeStack) =
    ImplicitNodeStack([state.node_stack; new_state.node_stack],
        joinstate(state.idx_stack, new_state.idx_stack))
joinstate(state::ImplicitIndexStack, new_state::ImplicitIndexStack) =
    ImplicitIndexStack([state.stack; new_state.stack])

isroot(tree, state::ImplicitStack) = isempty(state)


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

# ToDo: Add literature references to AdaptiveMetropolisTuning docstring.

"""
    AdaptiveMetropolisTuning(...)

Adaptive MCMC tuning strategy for Metropolis-Hastings samplers.

Adapts the proposal function based on the acceptance ratio and covariance
of the previous samples.

Fields:

* `λ`: Controls the weight given to new covariance information in adapting
  the proposal distribution. Defaults to `0.5`.

* `α`: Metropolis-Hastings acceptance ratio target, tuning will try to
  adapt the proposal distribution to bring the acceptance ratio inside this
  interval. Defaults to `IntervalSets.ClosedInterval(0.15, 0.35)`

* `β`: Controls how much the spread of the proposal distribution is
  widened/narrowed depending on the current MH acceptance ratio.

* `c`: Interval for allowed scale/spread of the proposal distribution.
  Defaults to `ClosedInterval(1e-4, 1e2)`.

Constructors:

```julia
AdaptiveMetropolisTuning(
    λ::Real,
    α::IntervalSets.ClosedInterval{<:Real},
    β::Real,
    c::IntervalSets.ClosedInterval{<:Real}
)
```
"""
@with_kw struct AdaptiveMetropolisTuning <: AbstractMCMCTuningStrategy
    λ::Float64 = 0.5
    α::IntervalSets.ClosedInterval{Float64} = ClosedInterval(0.15, 0.35)
    β::Float64 = 1.5
    c::IntervalSets.ClosedInterval{Float64} = ClosedInterval(1e-4, 1e2)
end

export AdaptiveMetropolisTuning

# Deprecate:
AbstractMCMCTuningStrategy(algorithm::MetropolisHastings) = AdaptiveMetropolisTuning()

(config::AdaptiveMetropolisTuning)(chain::MHIterator) = ProposalCovTuner(config, chain)

mutable struct ProposalCovTuner{
    S<:MCMCBasicStats
} <: AbstractMCMCTuner
    config::AdaptiveMetropolisTuning
    stats::S
    iteration::Int
    scale::Float64
end

function ProposalCovTuner(
    config::AdaptiveMetropolisTuning,
    chain::MHIterator
)
    m = nparams(chain)
    scale = 2.38^2 / m
    ProposalCovTuner(config, MCMCBasicStats(chain), 1, scale)
end

isviable(tuner::ProposalCovTuner, chain::MHIterator) = nsamples(chain) >= 2

function tuning_init!(tuner::ProposalCovTuner, chain::MHIterator)
    Σ_unscaled = cov(getprior(getposterior(chain)))
    Σ = Σ_unscaled * tuner.scale

    next_cycle!(chain)
    chain.proposaldist = set_cov(chain.proposaldist, Σ)

    nothing
end

function tuning_update!(tuner::ProposalCovTuner, chain::MHIterator)
    config = tuner.config

    α_min = minimum(config.α)
    α_max = maximum(config.α)

    c_min = minimum(config.c)
    c_max = maximum(config.c)

    β = config.β

    t = tuner.iteration
    λ = config.λ
    c = tuner.scale
    Σ_old = Matrix(get_cov(chain.proposaldist))

    S = convert(Array, tuner.stats.param_stats.cov)
    a_t = 1 / t^λ
    new_Σ_unscal = (1 - a_t) * (Σ_old/c) + a_t * S

    α = eff_acceptance_ratio(chain)

    max_log_posterior = tuner.stats.logtf_stats.maximum

    if α_min <= α <= α_max
        chain.info = MCMCIteratorInfo(chain.info, tuned = true)
        @debug "MCMC chain $(chain.info.id) tuned, acceptance ratio = $(Float32(α)), proposal scale = $(Float32(c)), max. log posterior = $(Float32(max_log_posterior))"
    else
        chain.info = MCMCIteratorInfo(chain.info, tuned = false)
        @debug "MCMC chain $(chain.info.id) *not* tuned, acceptance ratio = $(Float32(α)), proposal scale = $(Float32(c)), max. log posterior = $(Float32(max_log_posterior))"

        if α > α_max && c < c_max
            tuner.scale = c * β
        elseif α < α_min && c > c_min
            tuner.scale = c / β
        end
    end

    Σ_new_raw = new_Σ_unscal * tuner.scale
    Σ_new = PDMat(cholesky(Positive, Hermitian(Σ_new_raw)))

    next_cycle!(chain)
    chain.proposaldist = set_cov(chain.proposaldist, Σ_new)
    tuner.iteration += 1

    nothing
end

function run_tuning_cycle!(
    callbacks,
    tuners::AbstractVector{<:ProposalCovTuner},
    chains::AbstractVector{<:MHIterator};
    kwargs...
)
    run_tuning_iterations!(callbacks, tuners, chains; kwargs...)
    tuning_update!.(tuners, chains)
    nothing
end

function run_tuning_iterations!(
    callbacks,
    tuners::AbstractVector{<:ProposalCovTuner},
    chains::AbstractVector{<:MHIterator};
    max_nsamples::Int64 = Int64(1000),
    max_nsteps::Int64 = Int64(10000),
    max_time::Float64 = Inf
)
    user_callbacks = mcmc_callback_vector(callbacks, eachindex(chains))

    combined_callbacks = broadcast(tuners, user_callbacks) do tuner, user_callback
        (level, chain) -> begin
            if level == 1
                get_samples!(tuner.stats, chain, true)
            end
            user_callback(level, chain)
        end
    end

    mcmc_iterate!(combined_callbacks, chains, max_nsamples = max_nsamples, max_nsteps = max_nsteps, max_time = max_time)
    nothing
end


text:
# extra_file.jl

my_f(x, y) = 2x + 3y

derivative(x, y) = ForwardDiff.derivative(x -> my_f(x, y), x)


text:
module Drosophila

using GZip
using LinearAlgebra, NMF
using Statistics, StatsBase
using FileIO, JLD2

include("../src/util.jl")
include("../src/infer.jl")
include("../src/normalize.jl")
include("../src/pointcloud.jl")

# ------------------------------------------------------------------------
# globals

const DATA = "proc/drosophila"

# ------------------------------------------------------------------------
# data i/o

module DataIO

include("../src/scrna.jl")

function counts(io::IO)
	genes = String[]
	lines = [
		let
			entry = split(line)
			push!(genes, entry[1])
			parse.(Int, entry[2:end])
		end for line in eachline(io)
	]
	count = Matrix(reduce(hcat,lines)')
	return count, genes
end

function counts(dir::String, subdir::String)
	count = reduce(∪,
		scRNA.load("$dir/$d") for d ∈ readdir(dir) if occursin(subdir,d)
	)

    # remove any non melongaster cells (experiment specific)
	markers  = (
		dvir = scRNA.searchloci(count, "Dvir_"),
	)

	count = scRNA.filtercell(count) do cell, _
		sum(cell[markers.dvir]) < .10*sum(cell)
	end

    # remove any low count genes or non melongaster specific genes
	count = scRNA.filtergene(count) do count, gene
		!occursin("Dvir_", gene) && sum(count) > 20
	end

    # remove yolk and pole cells
	markers  = (
		yolk = scRNA.locus(count, "sisA", "CG8129", "Corp", "CG8195", "CNT1", "ZnT77C"),
		pole = scRNA.locus(count, "pgc"),
	)

	count = scRNA.filtercell(count) do cell, _
		(sum(cell[markers.yolk]) < 10
	  && sum(cell[markers.pole]) < 3
	  && sum(cell) > 1e4)
	end

    # remove any remaining lowly expressed genes
	count = scRNA.filtergene(count) do count, _
		sum(count) > 20
	end

	return count.data, count.gene
end

end

# ------------------------------------------------------------------------
# figures

module Figures

using Printf
using LaTeXStrings
using ColorSchemes
using Makie, CairoMakie
using Statistics, StatsBase
using SpecialFunctions: erf

using ..Normalize, ..PointCloud

const FIG = "figs/drosophila"

function mkpath()
    isdir(FIG) || Base.mkpath(FIG)
end

function heteroskedastic(raw)
    fig = Figure(font="Latin Modern Math", fontsize=26)

	ax1 = Axis(
		fig[1,1],
		xscale=Makie.log10,
		xlabel="cell sequencing depth",
        # xticks=[10^4.0, 10^4.5, 10^5.0], [L"10^4", L"10^{4.5}", L"10^5"]),
		ylabel="cumulative",
        yticks=([0, 0.25, 0.5, 0.75, 1.0], [L"0.0", L"0.25", L"0.50", L"0.75", L"1.00"])
	)

	ax2 = Axis(
		fig[1,2],
		xlabel="gene sequencing depth",
		xscale=Makie.log10,
        xticks=([1e1, 1e2, 1e3, 1e4, 1e5, 1e6], [L"10^1", L"10^2", L"10^3", L"10^4", L"10^5", L"10^6"]),
		ylabel="cumulative",
        yticks=([0, 0.25, 0.5, 0.75, 1.0], [L"0.0", L"0.25", L"0.50", L"0.75", L"1.00"])
	)

	depth = vec(sum(raw,dims=1))
	level = vec(sum(raw,dims=2))

	lines!(ax1, sort(depth), range(0,1,size(raw,2)),
		linewidth=2,
		color=:black
	)
    xlims!(ax1, (1e4, 3e5))

	lines!(ax2, sort(level), range(0,1,size(raw,1)),
		linewidth=2,
		color=:black
	)

    save("$FIG/heteroskedastic.png", fig, px_per_unit=2)
end

function overdispersed(raw)
	depth = sum(raw,dims=1)
    χ = Normalize.logmean(depth|>vec) .* (raw ./ depth)

	# figure 1
    fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		xscale=Makie.log10,
		yscale=Makie.log10,
		xlabel="mean / gene",
        xticks=([1e-2, 1e0, 1e2, 1e4], [L"10^{-2}", L"10^0", L"10^2", L"10^4"]),
		ylabel="variance / gene",
        yticks=([1e-2, 1e0, 1e2, 1e4], [L"10^{-2}", L"10^0", L"10^2", L"10^4"]),
    )

	scatter!(ax1,
		vec(mean(χ,dims=2)), vec(var(χ,dims=2)),
		color=(:black,0.03),
	)
	lines!(ax1,
		1e-2:1e4, 1e-2:1e4,
		linewidth=2,
		color=:red,
		linestyle=:dashdot,
		label="poisson",
	)
	ylims!(ax1, 1e-2, 1e4)
	axislegend(ax1, position = :rb)
    save("$FIG/overdispersed_mean_vs_variance.png", fig, px_per_unit=2)

	# figure 2
    fig = Figure(font="Latin Modern Math", fontsize=26)
	ax2 = Axis(fig[1,1],
		xscale=Makie.pseudolog10,
		xlabel="expected number of zeros (poisson)",
        xticklabelrotation=π/5,
		yscale=Makie.pseudolog10,
		ylabel="measured number of zeros",
    )

	g = vec(mean(χ,dims=2))
	z = vec(sum(χ.==0, dims=2))

	scatter!(ax2,
		size(χ,2)*exp.(-g),
		z,
		color=(:black,0.05),
	)
	lines!(ax2,
		1:1500, 1:1500,
		linewidth=2,
		color=:red,
		linestyle=:dashdot,
		label="poisson",
	)
	axislegend(ax2, position = :rb)

    save("$FIG/overdispersed_zeros.png", fig, px_per_unit=2)
end

function nb_uncertainty(raw, model)
    χ = [ gene |> vec |> Normalize.logmean for gene in eachrow(raw) ]

    # figure 1
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		xscale=Makie.log10,
		yscale=Makie.log10,
		xlabel=L"\chi",
        xticks=([1e-2, 1e0, 1e2], [L"10^{-2}", L"10^0", L"10^2"]),
		ylabel=L"Tr\left[\delta\Theta^2\right]",
        yticks=([1e0, 1e5, 1e10], [L"10^0", L"10^5", L"10^{10}"])
	)

    scatter!(ax1, χ, sqrt.(model.δΘ₁.^2 .+ model.δΘ₂.^2 .+ model.δΘ₃.^2),
		color=model.likelihood,
		colormap=:inferno
	)
	hlines!(ax1, [1e1], color=:turquoise, linewidth=4, linestyle=:dashdot)
    save("$FIG/nb_total_uncertainty_vs_expression.png", fig, px_per_unit=2)

    # figure 2
    fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		xscale=Makie.pseudolog10,
		yscale=Makie.log10,
		xlabel=L"\Theta_1",
        xticks=([-5, 0, +5], [L"-5", L"0", L"5"]),
		ylabel=L"\delta\Theta_L1",
        yticks=([1e-4, 1e-3, 1e-2, 1e-1, 1e0], [L"10^{-4}", L"10^{-3}", L"10^{-2}", L"10^{-1}", L"10^{0}"]),
	)
	scatter!(ax1, model.Θ₁, model.δΘ₁, color=log.(χ), colormap=:inferno)
    save("$FIG/nb_1_uncertainty_vs_expression.png", fig, px_per_unit=2)

    # figure 3
    fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		xscale=Makie.pseudolog10,
		yscale=Makie.log10,
		xlabel=L"\Theta_2",
        xticks=([-1, 0, 1, 2, 3], [L"-1", L"0", L"1", L"2", L"3"]),
		ylabel=L"\delta\Theta_2",
        yticks=([1e-3, 1e-2, 1e-1, 1e0], [L"10^{-3}", L"10^{-2}", L"10^{-1}", L"10^{0}"]),
	)
	scatter!(ax1, model.Θ₂, model.δΘ₂, color=log.(χ), colormap=:inferno)
    save("$FIG/nb_2_uncertainty_vs_expression.png", fig, px_per_unit=2)

    # figure 4
    fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		xscale=Makie.pseudolog10,
		yscale=Makie.log10,
		xlabel=L"\Theta_3",
        xticks=([1e-2, 1e0, 1e2, 1e4], [L"10^{-2}", L"10^{0}", L"10^{2}", L"10^{4}"]),
		ylabel=L"\delta\Theta_3",
        yticks=([1e-5, 1e0, 1e5, 1e10], [L"10^{-5}", L"10^{0}", L"10^{5}", L"10^{10}"]),
	)
	scatter!(ax1, model.Θ₃, model.δΘ₃, color=log.(χ), colormap=:inferno)
    save("$FIG/nb_3_uncertainty_vs_expression.png", fig, px_per_unit=2)
end

function nb_badfits(raw, model)
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(
		fig[1,1],
		xlabel="gene expression",
		ylabel="cumulative",
        yticks=([0, 0.25, 0.5, 0.75, 1.0], [L"0.0", L"0.25", L"0.50", L"0.75", L"1.00"])
	)

    δΘ  = sqrt.(model.δΘ₁.^2 .+ model.δΘ₂.^2 .+ model.δΘ₃.^2)
    bad = findall(δΘ .> 1e1)

    for g in sample(bad, min(100, length(bad)); replace=false)
		lines!(ax1,
			sort(vec(raw[g,:])), range(0,1,size(raw,2)),
			color=(:black,.05),
		)
	end
    save("$FIG/nb_badfits.png", fig, px_per_unit=2)
end

function nb_params(raw, model)
    δΘ = sqrt.(model.δΘ₁.^2 .+ model.δΘ₂.^2 .+ model.δΘ₃.^2)
	ι  = findall(δΘ .< 1e1)

    # figure 1
    fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(
		fig[1,1],
		xlabel=L"\Theta_2",
		ylabel="cumulative",
        yticks=([0, 0.25, 0.5, 0.75, 1.0], [L"0.0", L"0.25", L"0.50", L"0.75", L"1.00"])
	)
	lines!(ax1,
		sort(model.Θ₂[ι]), range(0,1,length(ι)),
		color=:black,
		linewidth=2,
	)
    save("$FIG/nb_param2.png", fig, px_per_unit=2)

    # figure 2
    fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(
		fig[1,1],
		xlabel=L"\Theta_3",
        xscale=Makie.log10,
        xticks=([1e-2, 1e-1, 1e0, 1e1], [L"10^{-2}", L"10^{-1}", L"10^0", L"10^1"]),
		ylabel="cumulative",
        yticks=([0, 0.25, 0.5, 0.75, 1.0], [L"0.0", L"0.25", L"0.50", L"0.75", L"1.00"])
	)
	lines!(ax1,
		sort(model.Θ₃[ι]), range(0,1,length(ι)),
		color=:black,
		linewidth=2,
	)
    save("$FIG/nb_param3.png", fig, px_per_unit=2)
end

function bootstrap(count, model, fits, index)
	χ = [ gene |> vec |> Normalize.logmean for gene in eachrow(count) ]

	function kernel(i)
		μ = [ median(b[i]) for b in fits ]
		δ = [ var(b[i]) for b in fits ]
		return μ, δ
	end

	Θ₁, δΘ₁ = kernel(1)
	Θ₂, δΘ₂ = kernel(2)
	Θ₃, δΘ₃ = kernel(3)

    # figure 1
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(
		fig[1,1],
		xlabel="MLE Θ₁",
		ylabel="bootstrap Θ₁",
	)
	scatter!(ax1, model.Θ₁[index], Θ₁, color=log.(χ[index]), colormap=:inferno)
	ax2 = Axis(
		fig[1,2],
		xlabel="MLE δΘ₁",
		xscale=Makie.log10,
        xticks=([1e-3, 1e-2, 1e-1, 1e0], [L"10^{-3}", L"10^{-2}", L"10^{-1}", L"10^{0}"]),
		ylabel="bootstrap δΘ₁",
		yscale=Makie.log10,
        yticks=([1e-3, 1e-2, 1e-1, 1e0], [L"10^{-3}", L"10^{-2}", L"10^{-1}", L"10^{0}"]),
	)
	scatter!(ax2, model.δΘ₁[index], δΘ₁, color=log.(χ[index]), colormap=:inferno)
    save("$FIG/bootstrap_1.png", fig, px_per_unit=2)

	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(
		fig[1,1],
		xlabel="MLE Θ₂",
		ylabel="bootstrap Θ₂",
	)
	scatter!(ax1, model.Θ₂[index], Θ₂, color=log.(χ[index]), colormap=:inferno)

	ax2 = Axis(
		fig[1,2],
		xlabel="MLE δΘ₂",
        xticks=([1e-3, 1e-2, 1e-1, 1e0], [L"10^{-3}", L"10^{-2}", L"10^{-1}", L"10^{0}"]),
		xscale=Makie.log10,
		ylabel="bootstrap δΘ₂",
		yscale=Makie.log10,
        yticks=([1e-3, 1e-2, 1e-1, 1e0], [L"10^{-3}", L"10^{-2}", L"10^{-1}", L"10^{0}"]),
	)
	scatter!(ax2, model.δΘ₂[index], δΘ₂, color=log.(χ[index]), colormap=:inferno)
    save("$FIG/bootstrap_2.png", fig, px_per_unit=2)

	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(
		fig[1,1],
		xlabel="MLE Θ₃",
		ylabel="bootstrap Θ₃",
		xscale=Makie.log10,
        xticks=([1e-1, 1e0, 1e1, 1e2], [L"10^{-1}", L"10^{0}", L"10^{1}", L"10^2"]),
		yscale=Makie.log10,
        yticks=([1e-2, 1e-1, 1e0, 1e1, 1e2], [L"10^{-2}", L"10^{-1}", L"10^{0}", L"10^{1}", L"10^2"]),
	)
	scatter!(ax1, model.Θ₃[index], Θ₃, color=log.(χ[index]), colormap=:inferno)
	ylims!(ax1, (1e-2, 1e2))

	ax2 = Axis(
		fig[1,2],
		xlabel="MLE δΘ₃",
        xticks=([1e-4, 1e-2, 1e0], [L"10^{-4}", L"10^{-2}", L"10^{0}"]),
		xscale=Makie.log10,
		ylabel="bootstrap δΘ₃",
		yscale=Makie.log10,
        yticks=([1e-4, 1e-2, 1e0, 1e2], [L"10^{-4}", L"10^{-2}", L"10^{0}", L"10^{2}"]),
	)
	scatter!(ax2, model.δΘ₃[index], δΘ₃, color=log.(χ[index]), colormap=:inferno)
	ylims!(ax2, (1e-5, 1e2))
    save("$FIG/bootstrap_3.png", fig, px_per_unit=2)
end

function estimated_rank(count, Λ)
    fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(
		fig[1,1],
		xlabel="singular value",
		xscale=Makie.log10,
        xticks=([1e2, 10^(2.5), 1e3, 10^(3.5)], [L"10^{2.0}", L"10^{2.5}", L"10^{3.0}", L"10^{3.5}"]),
		ylabel="rank",
		yscale=Makie.log10,
        yticks=([1e0, 1e1, 1e2, 1e3], [L"10^{0}", L"10^{1}", L"10^{2}", L"10^{3}"]),
	)
	λ = sum(sqrt.(size(count)))
	k = sum(Λ.S .>= λ)

	lines!(ax1, Λ.S, 1:length(Λ.S),
		color=:black,
		linewidth=2,
		label="empirical spectrum",
	)
	vlines!(ax1, [λ],
		color=:red,
		linestyle=:dashdot,
		linewidth=2,
		label="MP max value",
	)
	axislegend("Rank≈$(k)", position = :rt)
    save("$FIG/rank_estimate.png", fig, px_per_unit=2)

    fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(
		fig[1,1],
		xlabel="principal component",
		xscale=Makie.log10,
        xticks=([1e0, 1e1, 1e2, 1e3], [L"10^0", L"10^1", L"10^2", L"10^3"]),
		ylabel="participation ratio",
		yscale=Makie.log10,
        yticks=([10^2.6, 10^2.8, 10^3.0, 10^3.2, 10^3.4], [L"10^{2.6}", L"10^{2.8}", L"10^{3.0}", L"10^{3.2}", L"10^{3.4}"]),
	)

    IPR = ( sum(Λ.U.^4, dims=1) ./ sum(Λ.U.^2,dims=1) ) |> vec

    lines!(ax1, 1 ./ IPR,
           color=:black,
           linewidth=2
    )
    save("$FIG/participation_ratio.png", fig, px_per_unit=2)
end

function gamma_fits(model, count)
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		xlabel="theoretical quantile",
        xticks=([0, 0.25, 0.5, 0.75, 1.0], [L"0.0", L"0.25", L"0.50", L"0.75", L"1.00"]),
		ylabel="empirical quantile",
        yticks=([0, 0.25, 0.5, 0.75, 1.0], [L"0.0", L"0.25", L"0.50", L"0.75", L"1.00"]),
	)

    Z = .5*(1 .+ erf.( (model.residual) ./ sqrt(2) ))
    χ = [ Normalize.logmean(gene |> vec) for gene in eachrow(count) ]

	for i in 1:10:size(model.residual,1)
        χ[i] ≥ 1 || continue

        x = sort(Z[i,:] |> vec)
		lines!(ax1, sort(x), range(0,1,length(x)), color=(:black,.05))
	end

    save("$FIG/gamma_qq.png", fig, px_per_unit=2)
end

function scaling(D, N; δ=2)
    ϕ, Rs = PointCloud.scaling(D, N)

	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		xlabel="ball size",
		xscale=Makie.log10,
        xticks=([10^1, 10^1.5, 10^2.0, 10^2.5], [L"10^{1.0}", L"10^{1.5}", L"10^{2.0}", L"10^{2.5}"]),
		ylabel="points enclosed",
		yscale=Makie.log10,
        yticks=([10^1, 10^1.5, 10^2.0, 10^2.5, 10^3.0], [L"10^{1.0}", L"10^{1.5}", L"10^{2.0}", L"10^{2.5}", L"10^{3.0}"]),
	)
	for i in 1:δ:size(ϕ,1)
	    lines!(ax1, Rs, vec(ϕ[i,:]),
	           color=(:black,.01),
	           label="",
	    )
	end
    lines!(ax1, Rs, 6*size(ϕ,1)*(Rs ./ maximum(Rs)).^3,
           color=:red,
           linestyle=:dashdot,
           linewidth=2,
           label="3d scaling",
    )
	ylims!(ax1, (10, size(D,1)))

    save("$FIG/pointcloud_scaling.png", fig, px_per_unit=2)
end

function locality(Dg, Dx)
    x = PointCloud.upper_tri(Dx)
	g = PointCloud.upper_tri(Dg)

	bins = collect(range(50,350,length=10))

	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
		xlabel="expression distance (AU)",
		ylabel="cumulative",
		yticks=([0, 0.25, 0.5, 0.75, 1.0], [L"0.0", L"0.25", L"0.50", L"0.75", L"1.00"])
	)

	color = cgrad(:matter, length(bins)-1, categorical=true)
	for i in 1:length(bins)-1
		lb, ub = bins[i:i+1]
		x0 = (lb + ub) / 2

		data = g[lb .≤ x .≤ ub]
		lines!(ax1, sort(data), range(0,1,length(data)),
            color=color[i],
            linewidth=2,
			label=@sprintf "%.2f μm" x0)
	end
	axislegend("est. physical distance", position = :rb)
    save("$FIG/pointcloud_locality.png", fig, px_per_unit=2)
end

function inverse(mapping, database, residual; len=100)
    cmean   = zeros(len)
    cstd    = zeros(len)
    entropy = Array{Array{Float64,1},1}(undef,len)

    βs = range(1,2000,len)
    for (i,β) in enumerate(βs)
        ψ = minimum(size(residual))*mapping.invert(β)
        ψ = ψ ./ sum(ψ,dims=1)

        c = [cor( (ψ*(residual[g,:] |> vec)) , database.expression.real[:,i] |> vec) for (i,g) in enumerate(mapping.index) if g !== nothing]

        cmean[i]   = mean(c)
        cstd[i]    = std(c)
        entropy[i] = -sum(ψ .* log2.(ψ), dims=1) |> vec
    end

    # figure 1
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
        xlabel="β",
        ylabel="residual",
    )
    lines!(ax1, βs, cmean, color=:black, linewidth=2)
    errorbars!(ax1, βs, cmean, cstd./4, color=:red, linewidth=2)
    save("$FIG/bdtnp_fit.png", fig, px_per_unit=2)

    # figure 2
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
        xlabel="entropy (bits)",
        ylabel="cumulative",
        yticks=([0, 0.25, 0.5, 0.75, 1.0], [L"0.0", L"0.25", L"0.50", L"0.75", L"1.00"])
    )

    δ = 10
    colors = cgrad(:matter,length(entropy[1:δ:end]), categorical=true)
    for (i,S) in enumerate(entropy[1:δ:end])
        lines!(ax1, sort(S)./1.25, range(0,1,length(S));
               linewidth=2,
               color=colors[i],
               label=@sprintf "%.2f" βs[δ*(i-1)+1])
    end
    axislegend("β", position = :rb)
    save("$FIG/bdtnp_entropy.png", fig, px_per_unit=2)
end

function isomap(inverse, database, count, Z, D)
    # figure 1
	fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis(fig[1,1],
        xlabel="latent dimension",
        ylabel="distance correlation",
    )
    ξ = PointCloud.isomap(count,10;k=8);
    c = [
		let
			Dest = PointCloud.Distances.euclidean(ξ[:,1:i]')
			cor(PointCloud.upper_tri(Dest), PointCloud.upper_tri(D))
		end	for i in 1:size(ξ,2)
	]
    lines!(ax1, c, color=:black, linewidth=2)
    save("$FIG/isomap_dimension.png", fig, px_per_unit=2)

    ξ = PointCloud.isomap(Z,10;k=10);
    ψ = inverse.invert(150)
    X = ψ'*database.position

    # figure 2
    fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis3(fig[1,1],
		elevation=π/10,
		azimuth=5π/4,
        title="AP"
	)
	scatter!(ax1,
		ξ[:,1],ξ[:,2],ξ[:,3],
		color=X[:,1],
        colormap=:inferno,
        markersize=25
	)
	hidedecorations!(ax1)
    save("$FIG/isomap_AP.png", fig, px_per_unit=2)

    # figure 3
    fig = Figure(font="Latin Modern Math", fontsize=26)
	ax1 = Axis3(fig[1,1],
		elevation=π/10,
		azimuth=5π/4,
        title="DV"
	)
	scatter!(ax1,
		ξ[:,1],ξ[:,2],ξ[:,3],
        color=atan.(X[:,2],X[:,3]),
        colormap=:inferno,
        markersize=25
	)
	hidedecorations!(ax1)
    save("$FIG/isomap_DV.png", fig, px_per_unit=2)
end

const genelist = String["eve","odd","twi","ftz","zen"]

function gene_patterns(inverse, database, gene, count; genes=genelist)
    X = database.position
    ψ = inverse.invert(150)
    ψ = ψ ./ sum(ψ,dims=1)

    G = count*ψ'

    isdir("$FIG/gene") || Base.mkpath("$FIG/gene")

    for g in genes
        i = findfirst(gene .== g)

        fig = Figure(font="Latin Modern Math", fontsize=26)
        ax1 = Axis3(fig[1,1],
            elevation=0,
            azimuth=π/2,
        )
        scatter!(ax1,
            -X[:,1],X[:,2],X[:,3],
            color=G[i,:] |> vec,
            colormap=:inferno,
            markersize=15,
        )
        hidedecorations!(ax1)

        save("$FIG/gene/$g.png", fig, px_per_unit=2)
    end
end

end

# ------------------------------------------------------------------------
# main point of entry

function finalize(count, gene, model; cutoff=10)
    uncertainty = sqrt.(model.δΘ₁.^2 .+ model.δΘ₂.^2 .+ model.δΘ₃.^2)
	keep = findall(
		(uncertainty .<= cutoff) .&
		(vec(sum(count .== 0, dims=2)) .<= .99*size(count,2))
	)

	return count[keep,:], gene[keep], (
		likelihood = model.likelihood[keep],
		residual   = model.residual[keep, :],

		Θ₁ = model.Θ₁[keep],
		Θ₂ = model.Θ₂[keep],
		Θ₃ = model.Θ₃[keep],

		δΘ₁ = model.δΘ₁[keep],
		δΘ₂ = model.δΘ₂[keep],
		δΘ₃ = model.δΘ₃[keep],
	)
end

function bootstrap(count)
	χ = [ cell |> vec |> Normalize.logmean for cell in eachcol(count) ]
	ι = 1:10:size(count,1)

	data = Array{Tuple}(undef,length(ι))
	Threads.@threads for (n,i) in collect(enumerate(ι))
		data[n] = Normalize.bootstrap(vec(count[i,:]), χ)
	end

	return data, ι
end

function rescale(count, model)
	σ² = count.*(count.+model.Θ₃) ./ (1 .+ model.Θ₃)
	u², v², _ = Utility.sinkhorn(σ²)

	μ  = Diagonal(.√u²)*count*Diagonal(.√v²)
	return μ, .√u², .√v²
end

function main(root::String; subdir=nothing)
    isdir(DATA) || mkpath(DATA)

    # raw empirics
    count, gene = if subdir === nothing
        if isfile("$DATA/count_dvex.jld2")
            load("$DATA/count_dvex.jld2", "count", "gene")
        else
            c, g = open(DataIO.counts, root)
            save("$DATA/count_dvex.jld2", Dict("count"=>c, "gene"=>g))
            c, g
        end
    else
        if isfile("$DATA/count_raw.jld2")
            load("$DATA/count_raw.jld2", "count", "gene")
        else
            c, g = DataIO.counts(root, subdir)
            save("$DATA/count_raw.jld2", Dict("count"=>c, "gene"=>g))
            c, g
        end
    end

    Figures.mkpath()
    Figures.heteroskedastic(count)
    Figures.overdispersed(count)

    # negative binomial estimation
    model = if isfile("$DATA/nbmodel.jld2")
        load("$DATA/nbmodel.jld2", "model")
    else
        m = Normalize.glm(count)
        save("$DATA/nbmodel.jld2", Dict("model"=>m))
        m
    end

    Figures.nb_uncertainty(count, model)
    Figures.nb_badfits(count, model)
    Figures.nb_params(count, model)

    count, gene, model = finalize(count, gene, model)

    # bootstrap
    fits, index = if isfile("$DATA/bootstrap.jld2")
        load("$DATA/bootstrap.jld2", "fits", "index")
    else
        f, i = bootstrap(count)
        save("$DATA/bootstrap.jld2", Dict("fits"=>f, "index"=>i))
        f, i
    end
    Figures.bootstrap(count, model, fits, index)

    # rank of matrix
    norm, u, v, Λ = if isfile("$DATA/normalize.jld2")
        load("$DATA/normalized.jld2", "count", "u", "v", "Λ")
    else
        norm, u, v = rescale(count, model)
        Λ = svd(norm)
        save("$DATA/normalized.jld2", Dict("count"=>norm, "u"=>u, "v"=>v, "Λ"=>Λ))
        norm, u, v, Λ
    end

    Figures.estimated_rank(norm, Λ)

    # linear dimensional reduction
    μ, gene = if isfile("$DATA/nnmf.jld2")
        load("$DATA/nnmf.jld2", "matrix", "gene")
    else
        r = nnmf(norm,40;alg=:multmse,maxiter=200,init=:nndsvda)
        m = r.W*r.H
        save("$DATA/nnmf.jld2", Dict("matrix"=>m, "gene"=>gene))
        m, gene
    end

    # estimate differential expression
    model = if isfile("$DATA/gammamodel.jld2")
        load("$DATA/gammamodel.jld2", "model")
    else
        m = Normalize.glm(μ.+1e-8; stochastic=Normalize.gamma, ϵ=1e-8);
        save("$DATA/gammamodel.jld2", Dict("model"=>m))
        m
    end
    Figures.gamma_fits(model, count)

    # spatial inference
    Z = model.residual
    D = PointCloud.geodesics(μ,3)

    inverse, database = if isfile("$DATA/inverse.jld2")
        load("$DATA/inverse.jld2", "inverse", "database")
    else
        ν, ω = if isfile("$DATA/bdtnp_params.jld2")
            # load("$DATA/bdtnp_params.jld2", "ν", "ω")
            nothing, nothing
        else
            nothing, nothing
        end
        db = Inference.virtualembryo(directory="/home/nolln/mnt/data/drosophila/dvex");
        I  = Inference.inversion(norm./sum(norm,dims=1), gene; refdb=db, ν=ν, ω=ω);
        save("$DATA/inverse.jld2", Dict("inverse"=>I, "transform"=>I.invert(150), "database"=>db))
        I, db
    end

    Figures.scaling(D, 100)
    Figures.inverse(inverse, database, norm./sum(norm,dims=1))
    Figures.gene_patterns(inverse, database, gene, norm./sum(norm,dims=1))
    Figures.isomap(inverse, database, μ, Z, D)

    ξ = PointCloud.isomap(Z,10;k=10)
    if !isfile("$DATA/isomap.jld2")
        save("$DATA/isomap.jld2", Dict(
            "latent"=>ξ,
            "distance"=>D,
        ))
    end

    ψ = inverse.invert(150)
    ψ = ψ ./ sum(ψ,dims=1)
    D = ψ'*PointCloud.geodesics(database.position',6)*ψ
    Figures.locality(PointCloud.Distances.euclidean(ξ[:,1:5]'), D)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main("/home/nolln/mnt/data/drosophila/raw"; subdir="rep")
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
module Cmdl

"""
Cmdl.dict( ARGS )
"""
function dict( args::Array{String,1}; extra::Bool=false, debug::Bool=false )
 rv=Dict()
 lastkey=""
 let args = ( !extra && (x = findfirst( x->x=="--", args))>0 ) ? args[1:x-1] : args
  st = start(args)
  while !done(args,st)
   (arg,st)=next(args,st)
   debug && info(arg)
   if (m=match( r"^-(?<k>\w)(?<v>.+)", arg)) != nothing ||
      (m=match( r"^--(?<k>[\w\-]+)\=(?<v>.+)", arg)) != nothing
        push!( get!(rv, m[:k], []), m[:v] )
        lastkey=m[:k]
   elseif (m=match( r"^-(?<k>\w)$", arg)) != nothing ||
          (m=match( r"^--(?<k>[\w\-]+)$", arg)) != nothing
              get!(rv, m[:k],[])
              lastkey=m[:k]
   elseif !isempty( (vv=split( arg, r"\s+");) )
        append!( get!(rv, lastkey, []), vv )
   end
  end
 end
 rv
end

"""
tu = Cmdl.dict( ARGS ) |>
 opt( "file", mustbe=isfile, msg="--file must be a real file", musthave=true ) |>
 opt( "o", to=Int, mustbe=o->0<o<3) |>
 opt( "fi", musthave=1) |>
 opt( "asd",to=Int, musthave=2) |>
 opt( "a", to=Int) |>
 opt( "new", to=Int, mustbe=x->x>5) |>
 unexpected( warning=true)

Boolean options:

julia -e'1' -q -i --  --opt1 --opt2=true --opt3=false

julia> using Cmdl

julia> opts = Cmdl.dict( ARGS ) |>
        opt("opt1",to=Bool ) |>
        opt("opt2",to=Bool) |>
        opt("opt3",to=Bool) |>
        unexpected(warning=true)
Dict{Any,Any} with 3 entries:
  "opt3" => false
  "opt1" => true
  "opt2" => true
"""
function opt( n::AbstractString; to::Type=AbstractString, mustbe::Function=x->true, msg::AbstractString="$n !",
    musthave::Union{Bool,Int}=false )
 function a( tu::Tuple )
  cmddict::Dict = tu[1]
  found::Dict = (length(tu)>1) ? tu[2] : Dict()
  nvals = Base.get( cmddict, n, [] )
  if to==Bool
    v = false
    if haskey( cmddict, n)
      if length(nvals)>1 error("Option \"$n\" may have only one value, which is one of 'true','false' or nothing")
      elseif length(nvals)==1
        v = nvals[1]=="true"?true :
            nvals[1]=="false"?false:
            error("Option \"$n\" may have only 'true'|'false'")
      else # has key, but values array is isempty
        v = true
      end
    end
    get!( found, n, v)
  else
    vv=[]
    if typeof(musthave)<:Bool
      musthave && isempty(nvals) && error("Command line option \"$n\" must be defined: $msg")
    else
      musthave!=length(nvals) && error("Command line option \"$n\" must have $musthave value(s)")
    end
    for s in nvals
     v = (typeof(s)<:to) ? s : parse(to, s)
     !mustbe(v) && error( "$msg : $n $v")
     push!(vv,v)
    end
    append!( get!( found, n, []), vv)
  end
  return (cmddict, found)
 end

 a( cmddict::Dict) = a( (cmddict,Dict()) )
end
export opt

unexpected(;warning=false) = x->unexpected(x,warning=warning)

function unexpected(tu; warning=false)
 cmddict, needdict = tu
 for opt in keys( cmddict )
  !haskey( needdict, opt) && "unexpected option \"$opt\" in command line args"|> e->warning?warn(e):error(e)
 end
 needdict
end
export unexpected

#=
function args( wants...; extra=true )
 argd = dict(ARGS, extra=extra)
 rv = Dict()
 for w in wants
  if typeof(w)<:AbstractString
   get( argd, w, "")
  end
 end # не доделано
end
=#

#=
immutable WantOpt
 wantname::AbstractString
 wanttypestr::AbstractString
 wanttype::Type
 converting::Function
 re::Regex

 function WantOpt(n,s,t,f)
    if s=="b"
    r = Regex("-{1,2}$n")
    end
    if s=="i"
    r = Regex("-{1,2}$n.*=(\\d+)")
    end
    if s=="f"
    r = Regex("-{1,2}$n.*=(\\d+\\.?\\d?)")
    end
    if s=="s"
    r = Regex("-{1,2}$n.*=(.+)")
    end

    new(n,s,t,f,r)
 end
end

immutable GetOpt
 w::WantOpt
 rv
 from::AbstractString
end

function wantopt(w::AbstractString)
    (n,s,t) = name_str_type(w)
    f = Base.get(str_convert,s,nothing)
    WantOpt(n::AbstractString,s::AbstractString,t::Type,f::Function)
end

function wantopt{S<:AbstractString}(w::Pair{S,Function})
    (n,s,t) = name_str_type(w[1])
    WantOpt(n,s,t,w[2])
end

function name_str_type(w::AbstractString)
    m = match(r"(?P<wantname>.+?)\=(?P<wanttype>\w)", w)
    if m!=nothing
        t = Base.get(str_type, m[:wanttype], AbstractString)::Type
        return (m[:wantname], m[:wanttype], t)
    end

    m = match(r"(?P<wantname>.+)", w)
    if m!=nothing
        t = Base.get(str_type, "b", AbstractString)::Type
        return (m[:wantname],"b", t)
    end

    ("", "", AbstractString)
end

str_type = Dict( "i"=>Int, "f"=>Float64, "b"=>Bool )
str_convert = Dict(
    "i"=>s->parse(Int,s),
    "f"=>s->parse(Float64,s),
    "b"=>s->!isempty(s),
    "s"=>s->s
)

"""
    If you want to see full return structure, use find()
    It work same as get, except, it return full structure
     for wanted key(s).
"""
function find(wants...)
 wantopts = []
 for w in wants
    push!(wantopts, wantopt(w))
 end
 rv = findopt(wantopts)
end

findopt(ww::Vector{WantOpt}, given_args::Array) = map(w->findopt(w,ARGS) , wantopts)

function findopt(w::WantOpt, given_args::Array)
    rv = []
    for arg in given_args
    m = match(w.re, arg)
    if m!=nothing
        if isempty(m.captures)
            push!(rv, GetOpt(w,true,arg))
        else
            push!(rv, GetOpt(w, m.captures[1]|>w.converting, arg))
        end
    end
    end
    rv
end

"""
    If you wait multiple values for the same command line key,
     f.e. your script runned as:

    myscript -key1=123 -key1=234 -key2=345

    then:

    myarray = get(\"key1=i\")::Array

    (will empty array if not found key)
    or

    (arr1, arr2) = get(\"key1=i\",\"key2=i\")
"""
get() = nothing
function get(args...)
 rv = map( find(args...) ) do found_arr
      map(found_arr) do found
        found.rv
      end
 end
 l = length(rv)
 l==0 ? nothing :
 l==1 ? rv[1] :
 rv
end

"""
    # If your script runned as: myscript --mykey1=kuku -mykey2=lala
    # then you can get (string or nothing):

    val1 = first(\"mykey1=s\")::Union{AbstractString,Void}

    # two values at a time:

    (val1, val2) = first(\"mykey1=s\",\"mykey2=s\")

    first(),get(),find() arguments - wanted key descriptions
    in format \"keyname\" or \"keyname=T\",
    where T may be one of:
        s - AbstractString
        i - Int
        f - Float64
        b - Bool (default, if \"=T\" not present in format)

    Each wanted commandline argument will converted into required type.

    first() will raise error if got more then one value for ther one key.
    For multiple key values use get().

"""
function first(args...)
 rv = []
 for arg in args
    rv_arr = find(arg)[1]
    l = length(rv_arr)
    topush = l == 0 ? nothing :
         l == 1 ? Base.first(rv_arr).rv :
             begin
                error("""
Excpected single command line argument via first(\"$arg\").
Found $l matched arguments: $(join(map(a->a.from,rv_arr),",")).
        """)
        Base.first(rv_arr).rv
        end
    push!(rv, topush)
 end

 l = length(rv)
 l==0 ? nothing :
 l==1 ? rv[1] :
 rv
end
=#

end # module

max_stars_repo_name: closescreen/Cmdl

max_issues_repo_name: closescreen/Cmdl

max_forks_repo_name: closescreen/Cmdl


text:
function accum(m)
	tmp = 1
	for i in find(m)
		tmp *= abs(m[i])
	end
	tmp
end

function get_size(a, b, c)
	P = 1
	P *= accum(a)
	P *= accum(b)
	P *= accum(c)
	P = size(a)[1] * size(a)[2] + log10(abs(P))
	P = round(Int, P)
	println("size: ", P)
	P
end

function set_prec(n, L)
	println("prec: ", 2 * log2(n) + 2 * L)
end

function get_maxiter(b, c, n, R2, eps)
	p = sqrt(sum(x for x in b))
	m = sqrt(sum(x^2 for x in c))
	k = 2*(n+1)^2 * (abs(log(eps)) + log(1/p*m*R2))
	k = round(Int, k)
	k
end

a = Float64[1 0; 0 2; 3 2]
b = Float64[180, 150, 300]
c = Float64[-3, -5]

L = get_size(a, b, c)

n = size(a)[2]

set_prec(n, L)

x = zeros(2, 1)
R2 = Int128(n ^ 2) * (2 ^ (2 * L))

H = zeros(n, n)

for i in 1:n
	for j in 1:n
		if i == j
			H[i, j] = R2
		end
	end
end

eps = 1e-3
maxiter = get_maxiter(b, c, n, R2, eps)

function check(x, a, b)
	max = -Inf
	idx = 0
	for i in 1:size(a)[1]
		lhs = sum(a[i,:] .* x)
		rhs = b[i]
		#println(lhs, " >= ", rhs)
		if lhs >= rhs
			diff = lhs - rhs
			if diff > max
				max = diff
				idx = i
			end
		end
	end
	idx
end

for iter in 1:maxiter
	#println("\niter: ", iter)

	g = c

	idx = check(x, a, b)
	if idx != 0
		g = a[idx,:]
	end

	#println("g: ", g)

	Hg = H * g
	#println("Hg: ", Hg)

	Hgg = sum(x for x in (Hg .* g))
	#println("Hgg: ", Hgg)

	if Hgg < 0
		print("iters: ", iter)
		break
	end

	dx = - 1 / ((n + 1) * sqrt(Hgg)) * Hg

	x = x + dx
	H = n^2/(n^2-1) * (H - 2/(n+1) * (Hg * (g' * H)) / Hgg)

	#println("dx: ", dx)
	#println("x: ", x)
	#println("H: ", H)

	if iter == 2
		#break
	end
end

println("/", maxiter)

println("x: ", x)
println("f(x): ", c[1]*x[1] + c[2]*x[2])

max_stars_repo_name: vyacheslav-bezborodov/libem

max_issues_repo_name: vyacheslav-bezborodov/libem

max_forks_repo_name: vyacheslav-bezborodov/libem


text:
include("make_gui_html.jl")
include("make_gui_block.jl")
include("make_gui_widget.jl")
include("make_gui_bootloader.jl")

make_gui_list(cur_widgets::Widget...) = [cur_widgets...]

include("_manipulate.jl")


text:
abstract type AbstractConstraintInfo end

get_component_name(d::AbstractConstraintInfo) = d.component_name

abstract type AbstractRangeConstraintInfo <: AbstractConstraintInfo end
abstract type AbstractRampConstraintInfo <: AbstractConstraintInfo end
abstract type AbstractStartConstraintInfo <: AbstractConstraintInfo end

""" Data Container to construct range constraints"""
struct DeviceRangeConstraintInfo <: AbstractRangeConstraintInfo
    component_name::String
    limits::MinMax
    additional_terms_ub::Vector{Symbol}
    additional_terms_lb::Vector{Symbol}
end

function DeviceRangeConstraintInfo(name::String, limits::MinMax)
    return DeviceRangeConstraintInfo(name, limits, Vector{Symbol}(), Vector{Symbol}())
end

function DeviceRangeConstraintInfo(name::String)
    return DeviceRangeConstraintInfo(
        name,
        (min = -Inf, max = Inf),
        Vector{Symbol}(),
        Vector{Symbol}(),
    )
end

get_component_name(d::DeviceRangeConstraintInfo) = d.component_name

struct DeviceTimeSeriesConstraintInfo
    bus_number::Int
    multiplier::Float64
    timeseries::Vector{Float64}
    range::DeviceRangeConstraintInfo
    function DeviceTimeSeriesConstraintInfo(
        bus_number,
        multiplier,
        timeseries,
        range_constraint_info,
    )
        return new(bus_number, multiplier, timeseries, range_constraint_info)
    end
end

get_component_name(d::DeviceTimeSeriesConstraintInfo) = get_component_name(d.range)
get_limits(d::DeviceTimeSeriesConstraintInfo) = d.range.limits
get_timeseries(d::DeviceTimeSeriesConstraintInfo) = d.timeseries

function DeviceTimeSeriesConstraintInfo(
    device::PSY.Device,
    multiplier_function::Function,
    ts_vector::Vector{Float64},
    get_constraint_values::Union{Function, Nothing} = nothing,
)
    name = PSY.get_name(device)
    bus_number = PSY.get_number(PSY.get_bus(device))
    multiplier = multiplier_function(device)
    if isnothing(get_constraint_values)
        range_constraint_info = DeviceRangeConstraintInfo(name)
    else
        range_constraint_info =
            DeviceRangeConstraintInfo(name, get_constraint_values(device))
    end
    return DeviceTimeSeriesConstraintInfo(
        bus_number,
        multiplier,
        ts_vector,
        range_constraint_info,
    )
end

struct DeviceMultiStartRangeConstraintsInfo <: AbstractRangeConstraintInfo
    component_name::String
    limits::PSY.Min_Max
    lag_ramp_limits::NamedTuple{(:startup, :shutdown), Tuple{Float64, Float64}}
    additional_terms_ub::Vector{Symbol}# [:R1, :R2]
    additional_terms_lb::Vector{Symbol}
end

function DeviceMultiStartRangeConstraintsInfo(
    component_name::String,
    limits::PSY.Min_Max,
    lag_ramp_limits::NamedTuple{(:startup, :shutdown), Tuple{Float64, Float64}},
)
    return DeviceMultiStartRangeConstraintsInfo(
        component_name,
        limits,
        lag_ramp_limits,
        Vector{Symbol}(),
        Vector{Symbol}(),
    )
end

function DeviceMultiStartRangeConstraintsInfo(name::String)
    return DeviceMultiStartRangeConstraintsInfo(
        component_name,
        (min = -Inf, max = Inf),
        (startup = Inf, shutdown = Inf),
        Vector{Symbol}(),
        Vector{Symbol}(),
    )
end

struct DeviceRampConstraintInfo <: AbstractRampConstraintInfo
    component_name::String
    limits::PSI.MinMax
    ic_power::InitialCondition
    ramp_limits::PSI.UpDown
    additional_terms_ub::Vector{Symbol}
    additional_terms_lb::Vector{Symbol}
end

function DeviceRampConstraintInfo(
    component_name::String,
    limits::PSY.Min_Max,
    ic_power::InitialCondition,
    ramp_limits::PSI.UpDown,
)
    return DeviceRampConstraintInfo(
        component_name,
        limits,
        ic_power,
        ramp_limits,
        Vector{Symbol}(),
        Vector{Symbol}(),
    )
end

function DeviceRampConstraintInfo(name::String, ic_power::InitialCondition)
    return DeviceRampConstraintInfo(
        component_name,
        (min = -Inf, max = Inf),
        ic_power,
        (up = Inf, down = Inf),
    )
end

get_limits(d::DeviceRampConstraintInfo) = d.limits
get_ic_power(d::DeviceRampConstraintInfo) = d.ic_power
get_ramp_limits(d::DeviceRampConstraintInfo) = d.ramp_limits
get_additional_terms_ub(d::DeviceRampConstraintInfo) = d.additional_terms_ub
get_additional_terms_lb(d::DeviceRampConstraintInfo) = d.additional_terms_lb

struct DeviceStartUpConstraintInfo <: AbstractStartConstraintInfo
    component_name::String
    time_limits::StartUpStages
    startup_types::Int
end

struct DeviceStartTypesConstraintInfo <: AbstractStartConstraintInfo
    component_name::String
    startup_types::Int
end


text:
#############################################################
#Basic types and constructors

mutable struct AE # basic autoencoder type (MLP)
	nl::Int64
	x::Array{Array{Float64, 1}, 1}
	ax::Array{Array{Float64, 1}, 1}
	e::Array{Array{Float64, 1}, 1}
	w::Array{Array{Float64, 2}, 1}
	b::Array{Array{Float64, 1}, 1}
end

function AE(ns::Array{Int64, 1})
	nl = length(ns)
	AE(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			[randn(ns[i+1], ns[i])/(10*sqrt(ns[i])) for i in 1:nl - 1],
			[zeros(ns[i])/10 for i in 2:nl])
end

#Autoencoder learnt by sparse coding (and with tied weights)
#Minimizes lasso problem (see, e.g. Mairal and Bach 2009)
mutable struct AE_sparse
	nl::Int64
	x::Array{Array{Float64, 1}, 1}
	ax::Array{Array{Float64, 1}, 1}
	e::Array{Array{Float64, 1}, 1}
	w::Array{Array{Float64, 2}, 1}
	b::Array{Array{Float64, 1}, 1}
	A::Array{Array{Float64, 2}, 1}
	B::Array{Array{Float64, 2}, 1}
end

#init for sparse coding AE (no biases)
function AE_sparse(ns::Array{Int64, 1})
	nl = length(ns)
	AE_sparse(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			[randn(ns[i+1], ns[i])/(10*sqrt(ns[i])) for i in 1:nl - 1],
			[5 .* ones(ns[i]) for i in 2:nl],
			[zeros(ns[i], ns[i]) for i in 2:nl],#[rand(ns[i], ns[i])/5 for i in 2:nl],
			[zeros(ns[i], ns[i+1]) for i in 1:nl - 1])
end

mutable struct AE_sparse_patchy
	n_AE_sparse_patches::Int64
	n_out_per_population::Int64
	patchsize::Int64
	stride::Int64
	x::Array{Array{Float64, 1}, 1}
	AE_sparse_patches::Array{AE_sparse, 1}
end

function AE_sparse_patchy(n_in::Int64, n_out_per_population::Int64, n_of_populations::Int64;
					patchsize = 10, stride = 3, im_size = 28, colorimage = true)
	#n_AE_sparse_patches = Int(((im_size - patchsize)/stride + 1)^2)
	if colorimage
		ns = [3 * Int(patchsize^2), n_out_per_population]
	else
		ns = [Int(patchsize^2), n_out_per_population]
	end
	AE_sparse_patchy(n_of_populations,
		n_out_per_population,
		patchsize,
		stride,
		[zeros(n_in),zeros(n_of_populations * n_out_per_population)],
		[AE_sparse(ns) for i in 1:n_of_populations])
end

#same type as AE, different initialization
function Classifier(ns::Array{Int64, 1})
	nl = length(ns)
	AE(nl - 1,
		    [zeros(ns[i]) for i in 1:nl],
			[zeros(ns[i]) for i in 2:nl],
			[zeros(ns[i]) for i in 2:nl],
			[randn(ns[i+1], ns[i])/(10*sqrt(ns[i])) for i in 1:nl - 1],
			[rand(ns[i])/10 for i in 2:nl])
end

mutable struct SAE_sparse #stacked autoencoder
	nae::Int64
	aes::Array{AE_sparse, 1}
	reps::Array{Array{Array{Float64, 2}, 1}, 1} #AE -> train/test -> 2d array
end
function SAE_sparse(nAE::Int64,layers::Array{Array{Int64, 1}, 1},n_samples::Array{Int64, 1})
	SAE_sparse(nAE,
	[AE_sparse(layers[i]) for i in 1:nAE],
	[[zeros(layers[i][end],n_samples[j])  for j in [1,2]] for i in 1:nAE])
end

mutable struct SAE_sparse_patchy # stacked auto encoder
	nae::Int64
	aes::Array{AE_sparse_patchy, 1}
	reps::Array{Array{Array{Float64, 2}, 1}, 1} #AE -> train/test -> 2d array
end
function SAE_sparse_patchy(n_AE_sparse_patches::Int64,
		n_in::Int64, n_out_per_population::Int64;
		patchsize = 10, stride = 3, im_size = 28,
		n_samples = [60000,10000])
	SAE_sparse_patchy(1,
	[AE_sparse_patchy(n_in, n_out_per_population, n_AE_sparse_patches; patchsize = patchsize, stride = stride, im_size = im_size) for i in 1:nAE],
	[[zeros(n_out_per_population * n_AE_sparse_patches,n_samples[j]) for j in [1,2]]])
end

mutable struct SAEc_sparse #stacked sparse autoencoder with back-end classifier
	sae::SAE_sparse
	classifier::Array{AE, 1}
end

#############################################################
#Configuration Type

mutable struct configClass
	task::String
	hidden_sizes::Array{Int64, 1}
	n_inits::Int64
	iterations::Int64
	learningrates::Array{Float64, 1}
	weight_decay::Float64
	nonlinearity::Array{Function, 1}
	nonlinearity_diff::Array{Function, 1}
	nonlin_diff::Bool
end
mutable struct configAE
	task::String
	hidden_sizes::Array{Int64, 1}
  	n_inits::Int64
	iterations::Int64
	learningrates::Array{Float64, 1}
	weight_decay::Float64
	noise_type::String
	noise::Float64
	nonlinearity::Array{Function, 1}
	nonlinearity_diff::Array{Function, 1}
	nonlin_diff::Bool
end
mutable struct configAE_sparse
	task::String
	hidden_sizes::Array{Int64, 1}
  	n_inits::Int64
	iterations::Int64
	learningrates::Array{Float64, 1}
	lambda::Array{Float64, 1}
	weight_decay::Float64
	noise_type::String
	noise::Float64
	nonlinearity::Array{Function, 1}
	nonlinearity_diff::Array{Function, 1}
	nonlin_diff::Bool
end
mutable struct configSAE_sparse
	task::String
	n_samples::Array{Int64, 1}
	nAE::Int64
	hidden_sizes::Array{Array{Int64, 1}, 1}
  	n_inits::Int64
	iterations::Array{Int64, 1}
	learningrates::Array{Array{Float64, 1}, 1}
	lambdas::Array{Array{Float64, 1}, 1}
	weight_decay::Float64
	noise_type::String
	noise::Array{Float64, 1}
	nonlinearity::Array{Array{Function, 1}, 1}
	nonlinearity_diff::Array{Array{Function, 1}, 1}
	nonlin_diff::Bool
end
mutable struct configSAEc_sparse
	configSAE_sparse::configSAE_sparse
	configClass::configClass
end
function init_nonlin(nonlinarray, nonlindiffarray)
	nl = length(nonlinarray)
	nonlinearity = Array{Function, 1}()
	nonlinearity_diff = Array{Function, 1}()
	for i in 1:nl
		push!(nonlinearity,nonlinarray[i])
		push!(nonlinearity_diff,nonlindiffarray[i])
	end
	return nonlinearity, nonlinearity_diff
end


text:
# two main type systems - static and dynamic

#Integers
a = 1;
println(typeof(a));

println(Sys.WORD_SIZE);
println(Int);
println(UInt);

# Types of Integers

println(typemin(Int8));
println(typemax(Int8));
println(typemin(UInt8));
println(typemax(UInt8));

println(typemin(Int16));
println(typemax(Int16));
println(typemin(UInt16));
println(typemax(UInt16));

println(typemin(Int32));
println(typemax(Int32));
println(typemin(UInt32));
println(typemax(UInt32));

println(typemin(Int64));
println(typemax(Int64));
println(typemin(UInt64));
println(typemax(UInt64));

println(typemin(Int128));
println(typemax(Int128));
println(typemin(UInt128));
println(typemax(UInt128));

println(typemin(Bool));
println(false == 0);

println(typemax(Bool));

# Number Bases

b = 4; #decimal
b = 0b0100; # binary
b = 0o010; # octal
b = 0x1f; # hexadecimal

# overflow handled by wrap around

a = typemax(Int64);
a = a + 1;
println(a);

# BigInt -> larger than Int128

b = BigInt(typemax(Int128) + 1);
println(b);

# floating point numbers (IEEE 754)
c = 1.2
println(typeof(c));

println(typemin(Float16));
println(typemax(Float16));

Float16; # half precision 16-bit float
Float32; # single precision 32-bit float
Float64; # double precision 64-bit float

# Simple floats

f = 1.0;
f = 1.;
f = 0.1;
f = .1;

f = 22.5f-3; # single precision
println(typeof(f));

f = 22.5e-3; # double precision
println(typeof(f));

f = 0x1fp-1; # p - base2 exponent for 64-bit hex floats

# Special floats
Inf;
-Inf;
NaN;

println(Inf / Inf);

# BigFloat for big computations

a = parse(BigFloat, "0.00000000000000000000000000000000000000001");
println(a);

# static type a variable using ::
local a::Int8;

# This throws TypeError
# local a::Int8 = 128;

# isa -> check if an object is a particular type
a = 2;
println(isa(a, Int));

a = 2.0;
println(isa(a, Int));

# supertype -> see type's supertype or parent type

println(supertype(Float64));
println(supertype(AbstractFloat));
println(supertype(Real));
println(supertype(Number));
println(supertype(Any));

# Type Conversion

x = 12;
println(typeof(x));
println(typeof(convert(UInt8, x)));

# Type Promotion

println(promote(x, 2.3));

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
module JuliaSet

export GenerateJuliaSet

# package code goes here

# Populates a matrix A with the values to construct a Julia set image
function GenerateJuliaSet(R::Function, x::Array{Float64,1}, y::Array{Float64,1}, n_iter::Int64 = 200, escape_tol::Real = 2)

	nx = length(x)
	ny = length(y)
	A = zeros(nx, ny)
	for i in 1:nx
		for j in 1:ny
			z = x[i] + im * y[j]
			for k in 1:n_iter
				z = R(z)
				if (abs(z) > escape_tol)
					A[i,j] = k
				end
			end
			if A[i,j] == 0
				A[i,j] = escape_tol + 1 # Should be n_iter+1, but escape_tol gives prettier pictures
			end
		end
	end

	return A
end

end # module

max_stars_repo_name: restrin/JuliaSet

max_issues_repo_name: restrin/JuliaSet

max_forks_repo_name: restrin/JuliaSet


text:
#================
Preprocessing
     Module
================#
@doc """
      |====== Lathe.preprocess =====\n
      |____________/ Generalized Processing ___________\n
      |_____preprocess.TrainTestSplit(array)\n
      |_____preprocess.SortSplit(array)\n
      |_____preprocess.UniformSplit(array)\n
      |____________/ Feature Scaling ___________\n
      |_____preprocess.Rescalar(array)\n
      |_____preprocess.ArbitraryRescale(array)\n
      |_____preprocess.MeanNormalization(array)\n
      |_____preprocess.StandardScalar(array)\n
      |____________/ Categorical Encoding ___________\n
      |_____preprocess.OneHotEncode(array)\n

       """ ->
module preprocess
using Random
using DataFrames
using Lathe
#===============
Generalized
    Data
        Processing
===============#
# Train-Test-Split-----

function _dfTrainTestSplit(df,at = 0.75)
    sample = randsubseq(1:size(df,1), at)
    trainingset = df[sample, :]
    notsample = [i for i in 1:size(df,1) if isempty(searchsorted(sample, i))]
    testset = df[notsample, :]
    return(trainingset,testset)
end
function _ArraySplit(data, at = 0.7)
    n = length(data)
    idx = Random.shuffle(1:n)
    train_idx = view(idx, 1:floor(Int, at*n))
    test_idx = view(idx, (floor(Int, at*n)+1):n)
    data[train_idx,:], data[test_idx,:]
    return(test_idx,train_idx)
end
@doc """
      TrainTestSplit takes either a DataFrame or an Array and splits it according to the at parameter.\n
      --------------------\n
      [data] <- Iterable dictionary, dataframe, or Array.\n
      a <- Percentage value used to determine a point to split the data.\n
      -------------------\n
       """
TrainTestSplit(data::Array, at::Float64) = _ArraySplit(data,at)
TrainTestSplit(data::DataFrame, at::Float64) = _dfTrainTestSplit(data,at)
# Sort-Split -------------
@doc """
      SortSplit sorts the data from least to greatest, and then splits it,
      ideal for quartile calculations.\n
      --------------------\n
      array = [5,10,15]\n
      top25, lower75 = Lathe.preprocess.SortSplit(array,at = 0.75,rev = false)\n
      -------------------\n
      PARAMETERS:\n
      at:: Percentage value used to determine a point to split the data.\n
      rev:: Reverse, false by default, determines whether to sort least to
      greatest, or greatest to least.\n
       """
function SortSplit(data, at = 0.25, rev=false)
  n = length(data)
  sort!(data, rev=rev)  # Sort in-place
  train_idx = view(data, 1:floor(Int, at*n))
  test_idx = view(data, (floor(Int, at*n)+1):n)
  return(test_idx,train_idx)
end
# Unshuffled Split ----
@doc """
      Uniform Split does the exact same thing as ArraySplit(), but observations
      are returned split, but unsorted and unshuffled.\n
      --------------------\n
      array = [5,10,15]\n
      test, train = Lathe.preprocess.UniformSplit(array,at = 0.75)\n
      -------------------\n
      PARAMETERS:\n
      at:: Percentage value used to determine a point to split the data.
       """
function UniformSplit(data, at = 0.7)
    n = length(data)
    idx = data
    train_idx = view(idx, 1:floor(Int, at*n))
    test_idx = view(idx, (floor(Int, at*n)+1):n)
    data[train_idx,:], data[test_idx,:]
    return(test_idx,train_idx)
end
#=======
Numerical
    Scaling
=======#
# ---- Rescalar (Standard Deviation) ---
@doc """
      Rescalar scales a feature based on the minimum and maximum of the array.\n
      --------------------\n
      array = [5,10,15]\n
      scaled_feature = Lathe.preprocess.Rescalar(array)\n
      --------------------\n
      ==Functions==\n
      predict(xt) <- Returns a prediction from the model based on the xtrain value passed (xt)
       """ ->
function Rescalar(array)
    min = minimum(array)
    max = maximum(array)
    predict(array) = [i = (i-min) / (max - min) for i in array]
    (var) -> (predict)
end
# ---- Arbitrary Rescalar ----
@doc """
      Arbitrary Rescaling scales a feature based on the minimum and maximum
       of the array.\n
      --------------------\n
      array = [5,10,15]\n
      scaled_feature = Lathe.preprocess.Rescalar(array)\n
      --------------------\n
      ==Functions==\n
      predict(xt) <- Returns a prediction from the model based on the xtrain value passed (xt)
       """ ->
function ArbitraryRescale(array)
    a = minimum(array)
    b = maximum(array)
    predict(array) = [x = a + ((i-a*i)*(b-a)) / (b-a) for x in array]
    (var) -> (predict)
end
# ---- Mean Normalization ----
@doc """
      Mean Normalization normalizes the data based on the mean.\n
      --------------------\n
      array = [5,10,15]\n
      scaled_feature = Lathe.preprocess.MeanNormalization(array)\n
      --------------------\n
      ==Functions==\n
      predict(xt) <- Returns a prediction from the model based on the xtrain value passed (xt)
       """ ->
function MeanNormalization(array)
    avg = Lathe.stats.mean(array)
    a = minimum(array)
    b = maximum(array)
    predict(array) = [i = (i-avg) / (b-a) for i in array]
    (var) -> (predict)
end
# ---- Quartile Normalization ----
function QuartileNormalization(array)
    q1 = firstquar(array)
    q2 = thirdquar(array)

end
# ---- Z Normalization ----
@doc """
      Standard Scalar z-score normalizes a feature.\n
      --------------------\n
      array = [5,10,15]\n
      scaled_feature = Lathe.preprocess.StandardScalar(array)\n
      --------------------\n
      ==Functions==\n
      predict(xt) <- Returns a prediction from the model based on the xtrain value passed (xt)
       """ ->
function StandardScalar(array)
    q = Lathe.stats.std(array)
    avg = Lathe.stats.mean(array)
    predict(array) = [i = (i-avg) / q for i in array]
    (var) -> (predict)
end
# ---- Unit L-Scale normalize ----
@doc """
      FUNCTION NOT YET WRITTEN\n
      Unit L Scaling uses eigen values to normalize the data.\n
      --------------------\n
      array = [5,10,15]\n
      scaled_feature = Lathe.preprocess.UnitLScale(array)\n
       """ ->
function UnitLScale(array)

end
#==========
Categorical
    Encoding
==========#
# <---- One Hot Encoder ---->

@doc """
      One hot encoder replaces a single feature with sub arrays containing
      boolean values (1 or 0) for each individual category.\n
      --------------------\n
      df = DataFrame(:A => ['w','b','w'], :B => [5, 10, 15])\n
      scaled_feature = Lathe.preprocess.OneHotEncode(df,:A)\n
       """
function OneHotEncode(df,symb)
    copy = df
    for c in unique(copy[!,symb])
    copy[!,Symbol(c)] = copy[!,symb] .== c
    end
    predict(copy) = copy
end
# <---- Invert Encoder ---->
#==
@doc """
      FUNCTION NOT YET WRITTEN\n
      Invert Encoder (Not written.)\n
      --------------------\n
      array = [5,10,15]\n
      scaled_feature = Lathe.preprocess.OneHotEncode(array)\n
       """ ->
       ==#
function InvertEncode(array)

end
#-----------------------------
end