text:
module TBHamiltonians

using StaticArrays, OffsetArrays, SparseArrays
using Distributed, LinearAlgebra, Arpack
using NearestNeighbors

export Orbital, Hopping, TBHamiltonian
export sort_hoppings!, dimension, hopping_type, norbs
export add_orbital!, add_hopping!

export BlochHamiltonian, BlochHamiltonian_noF

export hamiltoniank, hamiltoniank!

export KPath, KMesh, path, mpmesh
export EigenStates, KEigenStates, KEigenAmplitudes, KBandStates, KBandAmplitudes

export bands_full, bands_sparse

export latticebasis, reciprocalbasis, spanlattice, addlattice
export cdot, toSVector

include("path_mesh_structs.jl")
include("keigen_structs.jl")

include("path_mesh_utils.jl")
include("vector_utils.jl")
include("lattice_utils.jl")

include("tbham_structs.jl")
include("tbham_utils.jl")
include("tbham_builders.jl")
include("tbham_eval.jl")

include("blochham_structs.jl")
include("blochham_utils.jl")
include("blochham_eval.jl")

include("tbham_kspectral.jl")

include("presets.jl")

end #module


text:
using Test, MLToolkit.Datasets

@testset "Datasets" begin
    IGNORES = ["mnist", "cifar10"]

    for name in Datasets.DATASET_NAMES
        name in IGNORES && continue
        dataset = Dataset(name, 30_000)
    end

    @warn "Datasets in $IGNORES are not tested."
end

ext: jl


text:
using Flux
using Zygote
using Formatting
using IterTools
using LineSearches

using JuMP

import Base:ndims
# Captures the model state from model return and saves it.
# Returns just the output
function capture_state!(model_output, xt)
    # xt .= (model_output[2])
    copyto!(xt, model_output[2])
    return model_output[1]
end

function my_train!(train_data, val_data, model, opt; Epochs=200)

    ndims(xi::Zygote.Buffer) = ndims(copy(xi))  # ndims not defined for buffer?

    train_batches = size(nth(train_data, 1)[1][1], 2)
    x0 = init_state(model, train_batches)

    θ = Flux.trainable(model)
    ps = Flux.Params(θ)

    sample_loss(y1, y2) = norm(y1 - y2)^2
    loss(y1, y2) = mean(sample_loss.(y1, y2))

    train_loss_log = []

    for epoch in 1:Epochs
        # Reset the state
        x = x0
        new_x = Zygote.bufferfrom(x)

        batch_loss = []
        max_grads = []
        #  Loop through training data
        for (ut, yt) in train_data
            # ut, yt = nth(train_data, 1)
            # Recur state
            x = copy(new_x)
            new_x = Zygote.bufferfrom(x0)  # new buffer to store state

            #                                               fancy stuff to recur state
            train_loss, back = Zygote.pullback(() -> loss(yt, capture_state!(model(x, ut), new_x)), ps)

            append!(batch_loss, train_loss)
            # calculate gradients and update loss
            ∇J = back(one(train_loss))
            update!(opt, ps, ∇J)

            # print important training info...=
            grads = Iterators.Filter(!isnothing, ∇J[p] for p in ps)
            g_inf = maximum(norm.(grads))
            append!(max_grads, g_inf)
            printfmt("\r loss: {1:1.2E}", train_loss)
        end
        append!(train_loss_log, mean(batch_loss))
        printfmt("\nEpoch: {1:2d}\tTraining loss: {2:1.4f},\t|g|: {3:1.3E}\n", epoch, mean(batch_loss), maximum(max_grads))

    end
    return train_loss_log
end

function validate(model, val_data, stats; washout=100)
    sample_loss(y1, y2) = norm(y1 - y2)^2
    loss(y1, y2) = mean(sample_loss.(y1, y2))

    # Allows for case of multiple validation sets.
    res = map(enumerate(val_data)) do (id, (uv, yv))
        x0 = init_state(model, size(uv[1], 2))
        yest = model(x0, uv)[1]

        # MSE
        σy, μy = stats.σy, stats.μy
        mse = loss(σy .* yv[washout:end], σy .* yest[washout:end])

        # NRMSE
        mu = [mean(yv)]
        # nrmse = norm((yv - yest)[washout:end]) / norm(yv[washout:end] .- mu)
        square(x) = x.^2
        dy = (yv - yest)[washout:end]
        nmse = mean(square.(dy)) ./ mean(square.(yv[washout:end] .- mu))
        nrmse = sqrt.(nmse)

        res_i = Dict("u" => uv, "y" => yv, "yest" => yest, "mse" => mse,
                    "nrmse" => nrmse, "washout" => washout)

        return res_i
    end

    return res
end

mse(y1, y2) = mean(norm(y1[:, i] - y2[:, i]).^2 for i in 1:size(y1, 2))
# function train_observer!(model, xn, xt, input_data, opt; Epochs=200, regularizer=nothing, solve_tol=1E-5)
function train_observer!(model, data, opt; Epochs=200, regularizer=nothing, solve_tol=1E-5, min_lr=1E-7)
    θ = Flux.trainable(model)
    ps = Flux.Params(θ)

    mean_loss = [1E5]
    loss_std = []
    for epoch in 1:Epochs
        batch_loss = []
        for (xni, xi, ui) in data
            function calc_loss()
                xpred = model(xi, ui)[1]
                return mean(norm(xpred[:, i] - xni[:, i]).^2 for i in 1:size(xi, 2))
            end

            train_loss, back = Zygote.pullback(calc_loss, ps)

            # calculate gradients and update loss
            ∇J = back(one(train_loss))
            update!(opt, ps, ∇J)

            push!(batch_loss, train_loss)
            printfmt("Epoch: {1:2d}\tTraining loss: {2:1.4E} \t lr={3:1.1E}\n", epoch, train_loss, opt.eta)
        end

        # Print stats through epoch
        println("------------------------------------------------------------------------")
        printfmt("Epoch: {1:2d} \t mean loss: {2:1.4E}\t std: {3:1.4E}\n", epoch, mean(batch_loss), std(batch_loss))
        println("------------------------------------------------------------------------")
        push!(mean_loss, mean(batch_loss))
        push!(loss_std, std(batch_loss))

        # Check for decrease in loss.
        if mean_loss[end] >= mean_loss[end - 1]
            println("Reducing Learning rate")
            opt.eta *= 0.1
            if opt.eta <= min_lr  # terminate optim.
                return mean_loss, loss_std
            end
        end
    end
    return mean_loss, loss_std
end

function train_observer2!(model, f, g, nx, nu, opt; Epochs=200, regularizer=nothing, solve_tol=1E-5, batchsize=500)
    θ = Flux.trainable(model)
    ps = Flux.Params(θ)

    train_loss_log = []
    for epoch in 1:Epochs

        xt = 1.5 * rand(nx, batchsize) .- 0.2 |> device
        ut = 1.5 * rand(nu, batchsize) .- 0.2 |> device
        xn = f(xt, ut)
        yt = g(xt, ut)
        inputs = vcat(ut, yt)
        function calc_loss()
            xpred = model(xt, inputs)[1]
            return mean(norm(xpred[:, i] - xn[:, i]).^2 for i in 1:size(xt, 2))
        end

        train_loss, back = Zygote.pullback(calc_loss, ps)

        # calculate gradients and update loss
        ∇J = back(one(train_loss))
        update!(opt, ps, ∇J)

        append!(train_loss_log, train_loss)
        printfmt("Epoch: {1:2d}\tTraining loss: {2:1.4E}\n", epoch, train_loss)
        if train_loss < solve_tol
            break
        end
    end
end

function snlsdp!(train_data, val_data, model, opt; Epochs=500, β0=10000, patience=20, min_lr=1E-6)
    ndims(xi::Zygote.Buffer) = ndims(copy(xi))  # ndims not defined for buffer?

    batches = size(nth(train_data, 1)[1][1], 2)

    β = β0
    x0 = init_state(model, batches)
    θ = [Flux.trainable(model)..., x0]
    ps = Flux.Params(θ)

    best_model = deepcopy(model)

    square_el(M) = M.^2
    # Calc loss including barrier funcctions
    function calc_loss(utild, ytild)
        yest, x = collect(model(x0, utild))
        L = sum(sum(square_el.(yest .- ytild)) / length(ytild))
        LMIs = eval_LMIs(model)
        lmi_sizes = size.(LMIs, 1)  # normalize by lmi size

        ineqs = eval_inequalities(model)
        return L[1] - sum(log.(ineqs)) / length(ineqs) / β  - sum(logdet.(LMIs) ./ lmi_sizes) / β, x[end]
    end

    function check_constraints()
        # Check LMIs are all positive definite
        LMIs = eval_LMIs(model)
        for lmi in LMIs
            try # chol faster than eigs
                cholesky(lmi)
            catch e
                return false
            end
        end
        # Check inequalities
        ineqs = eval_inequalities(model)
        return all(ineqs .> 0)
    end

    best_loss = Inf

    no_decrease_counter = 0
    vloss_log = [Inf]
    tloss_log = [Inf]
    for epoch in 1:Epochs

        x = x0
        new_x = Zygote.bufferfrom(x)

        # Training loop
        tloss = []
        for (ut, yt) in train_data
            x = copy(new_x)
            new_x = Zygote.bufferfrom(x0)  # new buffer to store state

            train_loss, back = Zygote.pullback(() -> capture_state!(calc_loss(ut, yt), new_x), ps)
            old_θ = deepcopy(θ)

            ∇J = back(one(train_loss))
            update!(opt, ps, ∇J)  # step parameters

            grads = Iterators.Filter(!isnothing, ∇J[p] for p in ps)
            g_inf = maximum(norm.(grads))

            # backtracking line search
            Δθ = [p2 - p1 for (p1, p2) in zip(old_θ, ps)]
            α = 1
            ls_iters = 1
            valid = check_constraints()
            while ~valid
                α = α / 2
                for (pk, old_θk, Δθk) in zip(ps, old_θ, Δθ)
                    pk .= old_θk + α * Δθk
                end
                valid = check_constraints()

                ls_iters = ls_iters + 1
                if ls_iters > 200
                    break
                end
            end
            # printfmt("\rloss: {1:1.4E}\t ls_iters: {2:d} \t grad: {3:1.3E}", train_loss, ls_iters - 1, g_inf)
            printfmt("\rloss: {1:1.4E}\t ls_iters: {2:d} ", train_loss, ls_iters - 1)
            push!(tloss, train_loss)
        end
        push!(tloss_log, mean(tloss))

        # Check validation set performance
        vloss = []
        vx0 = init_state(model, size(nth(val_data, 1)[1][1], 2))
        for (uv, yv) in val_data
            yest = collect(model(vx0, uv))[1]
            L = sum(sum(square_el.(yest - yv)) / length(yv))
            push!(vloss, L[1])
        end

        if mean(vloss) < best_loss
            best_model = deepcopy(model)
            best_loss = mean(vloss)
        end

        push!(vloss_log, mean(vloss))
        printfmt("\nEpoch: {1:d}\t tloss: {2:1.4E}\t vloss: {3:1.4E}\t lr: {4:1.1E}\t β: {5:1.0E}\n", epoch, mean(tloss), mean(vloss), opt[2].eta, β)

        # Check for sufficient decrease - at least 0.1 percent improvement
        if vloss_log[end] <= minimum(vloss_log[1:end - 1]) - 0.001 * vloss_log[end]
            no_decrease_counter = 0
        else
            no_decrease_counter = no_decrease_counter + 1
        end

        if no_decrease_counter > patience
            no_decrease_counter = 0
            println("Decreasing learning rate and increase barrier parameter.")
            β = 100 * β
            opt[2].eta = opt[2].eta / 10

            if opt[2].eta < min_lr
                break
            end

        end
    end

    Flux.loadparams!(model, Flux.params(best_model))

    return tloss_log, vloss_log
end

max_stars_repo_name: imanchester/REN

max_issues_repo_name: imanchester/REN

max_forks_repo_name: imanchester/REN


text:
@testset "E321" begin
    msgs = lintstr("""
    function f(x)
        x + y
    end
    """)
    @test messageset(msgs) == Set([:E321])
    @test occursin("use of undeclared symbol", msgs[1].message)

    @test isempty(lintstr("""
    function f(x)
        @lintpragma("Ignore use of undeclared variable y")
        x + y
    end
    """))

    # from @pao
    @test_broken isempty(lintstr("""
    function f()
        @lintpragma("Ignore use of undeclared variable aone")
        addOne() = @withOneVar aone (aone + 1)
    end
    """))
end

s = """
function f(x)
    if x > 1
        local i = 1
        println(i)
        i = i + 1
    end
    i = 1
    i
end
"""
msgs = lintstr(s)
@test length(msgs) == 0

s = """
function f(x)
    try
        x > 1
        local i = 1
        println(i)
        i = i + 1
    end
    i = 1
    i
end
"""
msgs = lintstr(s)
@test_broken msgs[1].code == :I482
@test_broken occursin("used in a local scope", msgs[1].message)

s = """
function f()
    :(x + y)
end
"""
msgs = lintstr(s)
@test isempty(msgs)

if VERSION < v"0.5"
s = """
function f()
    open(readall, "tmp.txt")
end
"""
else
s = """
function f()
    read("tmp.txt", String)
end
"""
end
msgs = lintstr(s)
@test isempty(msgs)


text:
# String Cache
# ------------

using Printf: @printf

struct Record
    meta::UInt64
    data::String

    function Record(data::String)
        len = length(data) % UInt64
        meta = len
        if len > 0
            # mix the first and the last byte
            meta |= (codeunit(data, 1) % UInt64) << 56
            meta |= (codeunit(data, len) % UInt64) << 48
        end
        return new(meta, data)
    end
end

const EMPTY_RECORD = Record("")

mutable struct Stats
    hit::Int
    miss::Int
end

struct StringCache
    maxsize::Int
    records::Vector{Record}
    stats::Stats

    function StringCache(maxsize::Int)
        @assert maxsize > 0
        return new(maxsize, Record[], Stats(0, 0))
    end
end

function Base.show(io::IO, cache::StringCache)
    stats = cache.stats
    total = stats.hit + stats.miss
    @printf(
        io,
        "StringCache: maxsize = %d, #records = %d, cache hit = %.3f (%d/%d)",
        cache.maxsize,
        length(cache.records),
        stats.hit / total,
        stats.hit,
        total,
    )
end

# LRU caching
function allocate!(cache::StringCache, p::Ptr{UInt8}, length::UInt64)
    meta = length | ((unsafe_load(p) % UInt64) << 56) | ((unsafe_load(p, length) % UInt64) << 48)
    records = cache.records
    n = Base.length(records)
    stats = cache.stats
    @inbounds for i in 1:n
        r = records[i]
        if r.meta == meta && memcmp(p, pointer(r.data), length) == 0
            stats.hit += 1
            for j in i:-1:2
                records[j] = records[j-1]
            end
            records[1] = r
            return r.data
        end
    end
    # not found
    stats.miss += 1
    string = unsafe_string(p, length)
    newrecord = Record(string)
    if n < cache.maxsize
        push!(records, EMPTY_RECORD)
        n += 1
    end
    @inbounds for j in n:-1:2
        records[j] = records[j-1]
    end
    records[1] = newrecord
    return string
end

@inline function memcmp(p1::Ptr, p2::Ptr, length::UInt64)
    return ccall(:memcmp, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), p1, p2, length)
end


text:
using FStrings
using Test

@testset "FStrings.jl" begin
    # Standard use cases
    x = 2
    y = 3
    @test f"{y}, {x}, {x+y}," == "3, 2, 5,"
    @test f"x={x:.2f}" == "x=2.00"
    @test f"π ≈ {π:.2f}" == "π ≈ 3.14"
    @test f"π ≈ {π:.4f}" == "π ≈ 3.1416"
    @test f"10π ≈ 0x{31:04x}" == "10π ≈ 0x001f"

    # Test trivial case
    @test f"" == ""

    # Test special characters
    # f_str macro receives escaped strings, which need to be unescaped:
    @test f"\n" == "\n"

    # `%` need to be escaped to `%%` before passed to `Printf.@sprintf`:
    @test f"%" == "%"
    @test f"""{"%"}""" == "%"

    # Potential extension to binary format
    # @test f"{123:b}" == "1111011"
    # @test f"{123:5b}" == "11011"
    # @test f"{123:8b}" == "1111011"
    # @test f"{123:08b}" == "01111011"

    # Use case questionable:
    @test f"{[x^2 for x in 1:3]}" == "[1, 4, 9]"
    z = 11:13
    @test f"{z[2:end]}" == "12:13"

    # Assert that proper errors are thrown, if format specifier is malformed.
    @test_throws ArgumentError f"{0x2:02b}"  # Unfortunately bitfields are not supported
end


text:
function update_Z_v2!(s::State, c::Constants, d::Data, tuners::Tuners,
                      sb_ibp::Bool; use_repulsive::Bool=false)
  if 0.1 > rand()
    # update Z marginalizing over lam and gam
    update_Z_marg_lamgam!(s, c, d, sb_ibp, use_repulsive=use_repulsive)
  else
    update_Z!(s, c, d, sb_ibp, use_repulsive=use_repulsive)
  end
end

function update_Z_marg_lamgam!(j::Int, k::Int,
                               A::Vector{Vector{Float64}},
                               B0::Vector{Matrix{Float64}},
                               B1::Vector{Matrix{Float64}},
                               s::State, c::Constants, d::Data, sb_ibp::Bool;
                               use_repulsive::Bool=false)
  v = sb_ibp ? cumprod(s.v) : s.v
  Z0 = deepcopy(s.Z)
  Z0[j, k] = false
  lp0 = MCMC.log1m(v[k]) + log_dmix_nolamgam(Z0, A, B0, B1, s, c, d)

  Z1 = deepcopy(s.Z)
  Z1[j, k] = true
  lp1 = log(v[k]) + log_dmix_nolamgam(Z1, A, B0, B1, s, c, d)

  if use_repulsive
    # NOTE: log_penalty_repFAM = log repulsive term
    lp0 += log_penalty_repFAM(k, Z0, c.log_repulsive_fn)
    lp1 += log_penalty_repFAM(k, Z1, c.log_repulsive_fn)
  end

  if lp0 == lp1  # especially if lp0 == lp1 == (Inf or -Inf)
    lp0 = 0.0
    lp1 = 0.0
  end
  p1_post = 1 / (1 + exp(lp0 - lp1))
  new_Zjk_is_one = p1_post > rand()
  s.Z[j, k] = new_Zjk_is_one
end

function update_Z_marg_lamgam!(s::State, c::Constants, d::Data, sb_ibp::Bool;
                               use_repulsive::Bool=false)
  # Precompute mus0, mus1, sig
  mus0 = mus(false, s, c, d)
  mus1 = mus(true, s, c, d)
  sig = sqrt.(s.sig2 * c.temper)

  # Precompute A, B0, B1

  # size(A[i]) = N[i]
  A = [if s.eps[i] > 0
         [logdnoisy(i, n, s, c, d) for n in 1:d.N[i]]
       else
         zeros(d.N[i])
       end for i in 1:d.I]

  # size(B0[i]) = (N[i], J)
  B0 = [[logdmixture(s.y_imputed[i][n, j], mus0, sig[i], s.eta[false][i, j, :])
         for n in 1:d.N[i], j in 1:d.J] for i in 1:d.I]

  # size(B1[i]) = (N[i], J)
  B1 = [[logdmixture(s.y_imputed[i][n, j], mus1, sig[i], s.eta[true][i, j, :])
         for n in 1:d.N[i], j in 1:d.J]
        for i in 1:d.I]

  # NOTE: The updates of Z_{j,k} cannot be made asynchronous naively.
  for j in 1:d.J
    for k in 1:c.K
      update_Z_marg_lamgam!(j, k, A, B0, B1, s, c, d, sb_ibp,
                            use_repulsive=use_repulsive)
    end
  end
end


text:
# LinkStream tests
ls = LinkStream("test", Intervals([(0.0, 10.0)]), Set(), Dict())
add_node!(ls,"a")
@test ls.V == Set(["a"])
@test ls.E == Dict()
add_node!(ls,"b")
@test ls.V == Set(["a","b"])
@test ls.E == Dict()

l = Link("a to b", Intervals([(1.5,3.0)]), "a", "b", 1.0)
add_link!(ls,l)
@test ls.V == Set(["a","b"])
@test ls.E == Dict("a"=>Dict("b"=>l))
@test 2.0 ∈ times(ls,"a","b")
@test 2.0 ∈ times(ls,"b","a")
@test 1.4999 ∉ times(ls,"a","b")
@test 1.4999 ∉ times(ls,"b","a")
@test duration(ls) == 10.0
@test node_duration(ls) == 10.0
@test link_duration(ls) == 1.5
@test nodes(ls,1.6)==ls.V
@test nodes(ls,1.0)==ls.V
@test links(ls,1.6)==[l]
@test links(ls,1.0)==[]

l2 = Link("b to c", Intervals([(2.5,5.0)]), "b", "c", 1.0)
add_link!(ls,l2)
@test ls.V == Set(["a","b","c"])
@test ls.E == Dict("a"=>Dict("b"=>l),"b"=>Dict("c"=>l2))
@test 2.55 ∈ times(ls,"a","b")
@test 2.55 ∈ times(ls,"b","a")
@test 2.55 ∈ times(ls,"b","c")
@test 2.55 ∈ times(ls,"c","b")
@test 1.4999 ∉ times(ls,"b","c")
@test duration(ls) == 10.0
@test node_duration(ls) == 10.0
@test link_duration(ls) == 4.0/3.0
@test links(ls,1.499) == []
@test Set(links(ls,2.66)) == Set([l,l2])
@test links(ls,2.499) == [l]
@test links(ls,4.0) == [l2]
@test links(ls,5.1) == []

record!(ls,2.0,2.6,"a","c")
@test ls.V == Set(["a","b","c"])
@test 2.1 ∈ times(ls,"a","c")
@test 2.55 ∈ times(ls,"b","c")
@test 1.4999 ∉ times(ls,"a","c")
@test duration(ls) == 10.0
@test node_duration(ls) == 10.0
@test link_duration(ls) == 4.6/3.0

ls = LinkStream("test-figure-1",Intervals([(0.0,10.0)]))
load!(ls,"./test_data/link_stream_1.txt","abuv")
@test ls.T == Intervals([(0.0,10.0)])
@test ls.V == Set(["a","b","c","d"])
@test 1.55 ∈ times(ls,"a","b")
@test node_duration(ls) == 10.0
@test link_duration(ls) == 23.0/6.0
@test isapprox(density(ls),0.38;atol=0.01)
@test isapprox(density(ls,1.0),2.3333,atol=0.0001)
@test Set(nodes(ls,8.0)) == Set(["a","b","c","d"])
@test length((links(ls,8.0)))==4.0
@test links(ls,"c","d")==links(ls,"d","c")
@test links(ls,"c","d")==Link("c to d", Intervals(Tuple{Float64,Float64}[(6.0, 9.0)]), "c", "d", 1.0)
@test length(links(ls,"d"))==2.0
@test length(links(ls,"a"))==2.0
@test length(times(ls,"b","d"))==3.0
@test length(times(ls,"d","b"))==3.0
@test length(times(ls,"b","c"))==7.0
@test length(times(ls,"c","b"))==7.0
@test coverage(ls)==1.0
@test uniformity(ls)==1.0
@test node_clustering(ls,"a")==1.0
@test node_clustering(ls,"b")==0.375
@test node_clustering(ls,"c")==0.6
@test node_clustering(ls,"d")==0.5
@test node_clustering(ls)==0.61875
N_c=neighborhood(ls,"c")
@test keys(N_c)==Set(["a","b","d"])
@test N_c["a"].presence==Intervals([(2.0,5.0)])
@test N_c["b"].presence==Intervals([(1.0,8.0)])
@test N_c["d"].presence==Intervals([(6.0,9.0)])
@test degree(ls,"c")==1.3
@test times(ls)==[0.0,1.0,2.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0]


text:
precompile(Tuple{typeof(OhMyREPL.Prompt.insert_keybindings), Any})
precompile(Tuple{typeof(OhMyREPL.__init__)})


text:
# This file was generated, do not modify it. # hide
mach.model.min_samples_split  = 10
fit!(mach, rows=train_rows) # re-train with new hyper-parameter

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# converts exhaustive optical data to optical gradients
using LinearAlgebra
using GaussianDistributions
using SARA
using SARA: TemperatureProfile, stripe_to_global, get_relevant_T, iSARA
using Base.Threads: @threads
using SARA
using JSON
using HDF5

# path = "SARA/ScienceAdvances2021/inner/"
path = "../inner/"
include(path * "inner_load_data.jl")
using SARA: estimate_noise

path = "SARA/ScienceAdvances2021/inner/"
file = "Bi2O3_19F44_01_inner_loop_data.json"
f = JSON.parsefile(path * "data/" * file)
position, optical, rescaling_parameters, T_max, log10_τ = load_data(f)

# omat = [optical[k][j][i] for i in eachindex(optical[1][1]), j in eachindex(optical[1]), k in eachindex(optical)]
# println(maximum(omat))
########### estimating noise variance of normalized coefficient
l = .2 # from inner loop benchmarks
i = 3
x = position[i]
y = optical[i][1]
θ = rescaling_parameters[i]
k = iSARA(l, θ)
σ = estimate_noise(x, y, k) # make dependent on iSARA kernel
println(σ)

# using Plots
# plotly()
# C = Gaussian(k) | (x, y, σ^2)
# plot(x, mean(C, x), ribbon = 2std(C).(x))
# scatter!(x, y)

############## convert to global data
# dT = 20 # uncertainty in peak temperature in C
# dx = 10/1000 # uncertainty in position in mm
dT = 25 # uncertainty in peak temperature in C
dx = 50/1000 # uncertainty in position in mm
# dx = 100/1000 # uncertainty in position in mm
P = TemperatureProfile(dT, dx)
nout = 32 # we can still change this after the creation of the data
# constant_offset = Val(false)
# T_proportions = (.75, .99) # as a proportion of T_max, what data to convert to gradients
constant_offset_bool = true
constant_offset = Val(constant_offset_bool)
T_offset = [200, 10] # number of degrees from T_max we are generating data for
relevant_T = get_relevant_T(constant_offset, T_offset..., nout)

input_noise_bool = false
input_noise = Val(input_noise_bool)
nstripes = length(optical)
gradients = zeros(nout, nstripes)
var_gradients = zeros(nout, nstripes)
@threads for i in 1:nstripes
    println(i)
    kernel = iSARA(l, rescaling_parameters[i])
    gradients[:, i], var_gradients[:, i] = stripe_to_global(position[i], optical[i],
                                            σ, kernel,
                                            P, T_max[i], log10_τ[i], relevant_T, input_noise)
end

temperatures = zeros(nout, nstripes)
dwelltimes = zeros(nout, nstripes)
for i in 1:nstripes
    dwelltimes[:, i] .= log10_τ[i]
    temperatures[:, i] = relevant_T(T_max[i], log10_τ[i])
end

# savefile = "Bi2O3_19F44_01_outer_loop_data.h5"
# path = "../outer/data/"
# path = "SARA/ScienceAdvances2021/outer/data/"
path = ""
# how to read file name: l_2 means kernel lengthscale was .2,
# input_noise means input noise propagation was switched on
savefile = "Bi2O3_19F44_01_outer_optical_gradients_l_2_input_noise_0.h5"
f = h5open(path * savefile, "w")
f["dwelltimes"] = dwelltimes
f["temperatures"] = temperatures
f["gradients"] = gradients
f["var(gradients)"] = var_gradients
f["input noise"] = input_noise_bool
f["nstripes"] = nstripes
f["constant offset"] = constant_offset_bool
f["T_offset"] = T_offset
f["nout"] = nout
f["l"] = l
f["optical noise variance"] = σ^2
f["temperature uncertainty parameters"] = [dT, dx]
close(f)

# creates gradient map
# k = oSARA()
# G = Gaussian(k)
# temperatures = vec(temperatures)
# dwelltimes = vec(dwelltimes)
# x = tuple.(temperatures, dwelltimes)
# y = vec(target)
# σ² = vec(uncertainty) .+ .1
# C = G | (x, y, Diagonal(σ²))

# using Plots
# pyplot()
# T = 200:20.:1400
# τ = 2:.2:4
# surface(T, τ, (x,y)->mean(C)((x,y)))
# gui()


text:
# $Id: util.jl,v 1.1 2018/04/15 00:59:35 manabu Exp manabu $

function dereference(t::Any)
    t <: Cxx.CppRef ? Cxx.CppValue{Cxx.CxxQualType{(t.parameters[1]),Cxx.NullCVR},N} where N : t
end

function reference(t::Any)
    if typeof(t) == UnionAll
        t0 = t.body
    else
        t0 = t
    end
    t0 <: Cxx.CppValue ? Cxx.CppRef{((t0.parameters[1]).parameters[1]), Cxx.NulCVR} : t
end

function pointer(t::Any)
    t <: Cxx.CppValue ? Cxx.CppPtr{t.body.parameters[1],Cxx.NullCVR} : t
end

max_stars_repo_name: moyhig/CxxROS

max_issues_repo_name: moyhig/CxxROS

max_forks_repo_name: moyhig/CxxROS


text:
mutable struct GlobalCache
  time::Float64
  cache::Vector{DocObj}
end
const CACHE = GlobalCache(0.0, DocObj[])

CACHETIMEOUT = 30 # s

MAX_RETURN_SIZE = 20 # how many results to return at most

function searchdocs(needle::AbstractString; loaded::Bool = true, mod::Module = Main,
                    maxreturns::Int = MAX_RETURN_SIZE, exportedonly::Bool = false,
                    name_only::Bool = false)
  loaded ? dynamicsearch(needle, mod, exportedonly, maxreturns, name_only) :
           dynamicsearch(needle, mod, exportedonly, maxreturns, name_only, loaddocsdb())
end

# TODO:
# We may want something like `CodeTools.getmodule` here, so that we can accept `mod` as `String`:
# - then we can correctly score bindings even in unloaded packages
# - it would make `isdefined` checks below more robust -- currently it won't work when e.g.
#   we try to find `Atom.JunoDebugger.isdebugging` given `mod == Atom`
function dynamicsearch(needle::AbstractString, mod::Module = Main,
                       exportedonly::Bool = false, maxreturns::Int = MAX_RETURN_SIZE,
                       name_only::Bool = false, docs::Vector{DocObj} = alldocs(mod))
  isempty(docs) && return DocObj[]
  scores = zeros(size(docs))
  modstr = string(mod)
  Threads.@threads for i in eachindex(docs)
    scores[i] = score(needle, docs[i], modstr, name_only)
  end
  perm = sortperm(scores, rev=true)
  out = [(scores[p], docs[p]) for p in perm]

  f = if exportedonly
    if mod == Main
      x -> x[2].exported
    else
      let mod = mod, modstr = modstr
        x -> begin
          # filters out unexported bindings
          x[2].exported &&
          # filters bindings that can be reached from `mod`
          (
            isdefined(mod, Symbol(x[2].mod)) ||
            modstr == x[2].mod # needed since submodules are not defined in themselves
          )
        end
      end
    end
  else
    if mod == Main
      x -> true
    else
      let mod = mod, modstr = modstr
        x -> begin
          # filters bindings that can be reached from `mod`
          isdefined(mod, Symbol(x[2].mod)) ||
          modstr == x[2].mod # needed since submodules are not defined in themselves
        end
      end
    end
  end
  filter!(f, out)

  return out[1:min(length(out), maxreturns)]
end

function modulebindings(mod, exportedonly = false, binds = Dict{Module, Set{Symbol}}(), seenmods = Set{Module}())
  # This does fairly stupid things, but whatever. Works for now.
  for mod in Base.loaded_modules_array()
    mod in seenmods && continue
    push!(seenmods, mod)
    modulebindings(mod, exportedonly, binds, seenmods)
  end

  for name in names(mod, all=!exportedonly, imported=!exportedonly)
    startswith(string(name), '#') && continue
    if isdefined(mod, name) && !Base.isdeprecated(mod, name)
      obj = getfield(mod, name)
      !haskey(binds, mod) && (binds[mod] = Set{Symbol}())
      push!(binds[mod], name)
      if (obj isa Module) && !(obj in seenmods)
        push!(seenmods, obj)
        modulebindings(obj, exportedonly, binds, seenmods)
      end
    end
  end
  return binds
end

"""
    alldocs(topmod = Main) -> Vector{DocObj}

Find all docstrings in all currently loaded Modules.
"""
function alldocs(topmod = Main)::Vector{DocObj}
  time() - CACHE.time < CACHETIMEOUT && return CACHE.cache

  results = DocObj[]
  # all bindings
  modbinds = modulebindings(topmod, false)
  # exported bindings only
  exported = modulebindings(topmod, true)

  # loop over all loaded modules
  for mod in keys(modbinds)
    parentmod = parentmodule(mod)
    meta = Docs.meta(mod)

    # loop over all names handled by the docsystem
    for b in keys(meta)
      # kick everything out that is handled by the docsystem
      haskey(modbinds, mod) && delete!(modbinds[mod], b.var)
      haskey(exported, mod) && delete!(exported[mod], b.var)

      expb = (haskey(exported, mod) && (b.var in exported[mod])) ||
             (haskey(exported, parentmod) && (b.var in exported[parentmod]))

      multidoc = meta[b]
      for sig in multidoc.order
        d = multidoc.docs[sig]
        md = Markdown.parse(join(d.text, ' '))
        text = stripmd(md)
        path = d.data[:path] == nothing ? "<unknown>" : d.data[:path]
        dobj = DocObj(string(b.var), string(b.mod), string(determinetype(b.mod, b.var)),
                      # sig,
                      text, md, path, d.data[:linenumber], expb)
        push!(results, dobj)
      end
    end

    # resolve everything that is not caught by the docsystem
    for name in modbinds[mod]
      b = Docs.Binding(mod, name)

      # figure out how to do this properly...
      expb = (haskey(exported, mod) && (name in exported[mod])) ||
             (haskey(exported, parentmod) && (name in exported[parentmod]))

      if isdefined(mod, name) && !Base.isdeprecated(mod, name) && name != :Vararg
        # HACK: For now we don't need this -> free 50% speedup.
        # bind = getfield(mod, name)
        # meths = methods(bind)
        # if !isempty(meths)
        #   for m in meths
        #     dobj = DocObj(string(name), string(mod), string(determinetype(mod, name)),
        #                   "", Hiccup.div(), m.file, m.line, expb)
        #     push!(results, dobj)
        #   end
        # else
        #   dobj = DocObj(string(name), string(mod), string(determinetype(mod, name)),
        #                 "", Markdown.parse(""), "<unknown>", 0, expb)
        #   push!(results, dobj)
        # end
        dobj = DocObj(string(name), string(mod), string(determinetype(mod, name)),
                      "", Markdown.parse(""), "<unknown>", 0, expb)
        push!(results, dobj)
      end
    end
  end
  append!(results, keywords())
  results = unique(results)

  # update cache
  CACHE.time = time()
  CACHE.cache = results

  return results
end

function keywords()
  out = DocObj[]
  for k in keys(Docs.keywords)
    d = Docs.keywords[k]
    md = Markdown.parse(join(d.text, ' '))
    text = stripmd(md)
    dobj = DocObj(string(k), "Base", "Keyword", text, md, "", 0, true)
    push!(out, dobj)
  end
  return out
end

function determinetype(mod, var)
  (isdefined(mod, var) && !Base.isdeprecated(mod, var)) || return ""

  b = getfield(mod, var)

  b isa Function && return "Function"
  b isa UnionAll && return "DataType"

  string(typeof(b))
end


text:
# Julia wrapper for header: /usr/include/scip/cons_cardinality.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrCardinality(scip)
    ccall((:SCIPincludeConshdlrCardinality, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsCardinality(scip, cons, name, nvars, vars, cardval, indvars, weights, initial, separate, enforce, check, propagate, _local, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsCardinality, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nvars, vars, cardval, indvars, weights, initial, separate, enforce, check, propagate, _local, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicCardinality(scip, cons, name, nvars, vars, cardval, indvars, weights)
    ccall((:SCIPcreateConsBasicCardinality, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}), scip, cons, name, nvars, vars, cardval, indvars, weights)
end

function SCIPchgCardvalCardinality(scip, cons, cardval)
    ccall((:SCIPchgCardvalCardinality, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Cint), scip, cons, cardval)
end

function SCIPaddVarCardinality(scip, cons, var, indvar, weight)
    ccall((:SCIPaddVarCardinality, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var, indvar, weight)
end

function SCIPappendVarCardinality(scip, cons, var, indvar)
    ccall((:SCIPappendVarCardinality, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Ptr{SCIP_VAR}), scip, cons, var, indvar)
end

function SCIPgetNVarsCardinality(scip, cons)
    ccall((:SCIPgetNVarsCardinality, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetVarsCardinality(scip, cons)
    ccall((:SCIPgetVarsCardinality, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetCardvalCardinality(scip, cons)
    ccall((:SCIPgetCardvalCardinality, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetWeightsCardinality(scip, cons)
    ccall((:SCIPgetWeightsCardinality, libscip), Ptr{Cdouble}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end


text:
"""
Given a file path or paths to abundance tables (eg humann2 or metaphlan2),
create abundance table. Table is presumed to have samples in columns and
features in rows. First column is taken as feature IDs.
"""
function import_abundance_table(file::AbstractString; delim::Char='\t')
    @info "Importing abundance table" file
    df = CSV.File(file, delim=delim) |> DataFrame
    rename!(df, names(df)[1] => :col1)
    return df
end

function import_abundance_tables(files::Array{<:AbstractString, 1}; delim::Char='\t')
    @info "Importing abundance tables"
    fulltable = DataFrame(col1=String[])
    for t in files
        df = import_abundance_table(t, delim=delim)
        fulltable = join(fulltable, df, on=:col1, kind=:outer)
    end

    # replace all missing values (from joins) with 0.
    fulltable = map(c -> eltype(c) <: Union{<:Number, Missing} ? collect(Missings.replace(c, 0)) : c, eachcol(fulltable))
    return fulltable
end

function clean_abundance_tables(files::Array{String, 1};
                        delim::Char='\t',
                        col1::Symbol=:taxon,
                        suffix::String="_taxonomic_profile")
    @info "Cleaning"
    t = import_abundance_tables(files, delim=delim)

    rename!(t, :col1 => col1)
    rename!(n-> Symbol(replace(String(n), suffix => "")), t)

    return t
end

function rm_strat!(df::DataFrame; col::Union{Int, Symbol}=1)
    filter!(x->!occursin(r"\|", x[1]), df)
end

function permanova(dm::Array{<:Real,2}, metadata::AbstractVector, nperm::Int=999; filter=fill(true, length(metadata)))
    size(dm,1) != size(dm,2) && throw(ArgumentError("dm must be symetrical distance matrix"))
    size(dm,2) != length(metadata) && throw(ArgumentError("Metadata does not match the size of distance matrix"))

    r_dm = dm[filter, filter]
    metadata = metadata[filter]
    @rput r_dm
    @rput metadata

    R"""
    library(vegan)

    p <- adonis(r_dm ~ metadata,
            method = "bray", permutations = $nperm)
    """

    @rget p

    return p[:aov_tab]
end

#==============
PanPhlAn Utils
==============#

function panphlan_calcs(df::DataFrame)
    abun = abundancetable(df)
    dm = pairwise(Jaccard(), abun, dims=2)
    rowdm = pairwise(Jaccard(), abun, dims=1)
    col_clust = hclust(dm, :single)
    row_clust = hclust(rowdm, :single)
    optimalorder!(col_clust, dm)
    optimalorder!(row_clust, rowdm)

    mds = fit(MDS, dm, distances=true)

    return abun, dm, col_clust, row_clust, mds
end


text:
printstyled("  seismogram differentiation/integration\n", color=:light_green)

convert_file = path*"/SampleFiles/SEIS/2019_M7.1_Ridgecrest.seis"
S = randSeisData(3, s=1.0)
S.units = ["m/s2", "m", "m/s"]
U = deepcopy(S)
S1 = deepcopy(S)

printstyled("  conversion to m/s\n", color=:light_green)
redirect_stdout(out) do
  convert_seis!(S, v=3)
  @test S.units[1] == S.units[2] == S.units[3] == "m/s"
  @test S.x[3] == U.x[3]
end

printstyled("  conversion to m/s2\n", color=:light_green)
redirect_stdout(out) do
  convert_seis!(S, v=2, units_out="m/s2")
  @test S.units[1] == S.units[2] == S.units[3] == "m/s2"
  @test isapprox(S.x[1], U.x[1])
end

printstyled("  conversion to m\n", color=:light_green)
redirect_stdout(out) do
  convert_seis!(S1, v=2, units_out="m")
  @test S1.units[1] == S1.units[2] == S1.units[3] == "m"
  @test isapprox(S1.x[2], U.x[2])
end

printstyled("  conversion from m to m/s2 and back to m\n", color=:light_green)
redirect_stdout(out) do
  for j = 1:100
    S = randSeisData(3, s=1.0)
    S.units = ["m", "m", "m"]

    for i = 1:3
      Nx = min(length(S.x[i]), 20000)
      S.t[i] = [1 S.t[i][1,2]; Nx 0]
      S.x[i] = rand(Float64, Nx)
    end
    U = deepcopy(S)

    convert_seis!(S, units_out="m/s2")
    @test S.units[1] == S.units[2] == S.units[3] == "m/s2"
    convert_seis!(S, units_out="m")
    @test S.units[1] == S.units[2] == S.units[3] == "m"

    for i = 1:3
      @test isapprox(S.x[i], U.x[i])
    end
  end
end

printstyled("    at Float32 precision\n", color=:light_green)
redirect_stdout(out) do
  S = rseis(convert_file)[1]
  U = deepcopy(S)
  convert_seis!(S, units_out="m/s2", v=2)
  convert_seis!(S)
  for i = 1:16
    @test isapprox(S.x[i], U.x[i])
  end
end

printstyled("    at Float64 precision\n", color=:light_green)
redirect_stdout(out) do
  S = rseis(convert_file)[1]
  for i = 1:S.n
    S.x[i] = Float64.(S.x[i])
  end
  detrend!(S)

  U = deepcopy(S)
  convert_seis!(S, units_out="m/s2")
  for i = 17:19
     @test S.x[i] == U.x[i]
  end
  convert_seis!(S)
  for i = 1:16
    @test isapprox(S.x[i], U.x[i])
    @test isapprox(S.t[i], U.t[i])
  end

  convert_seis!(S, units_out="m")
  convert_seis!(S, units_out="m/s")
  convert_seis!(S, units_out="m")
  convert_seis!(S, units_out="m/s2")
  T = convert_seis(S, units_out="m")

  # Test on a SeisChannel
  C = deepcopy(U[1])
  convert_seis!(C, units_out="m", v=1)
  convert_seis!(C, units_out="m/s", v=1)
  convert_seis!(C, units_out="m/s2", v=1)
  convert_seis!(C, units_out="m", v=1)
  D = convert_seis(C, units_out="m/s2", v=1)
end


text:
using Mortgage

include("mortgageArguments.jl")

function main(x...)
	if !isempty(x)
		args = split(x[1])
	else
		if isempty(ARGS)
			println("Please enter command line arguments for mortgage.jl, or hit \"enter\" for defaults:")
			args = split(readline(STDIN))
		else
			args = ARGS
		end
	end
	p = parse_commandline(args)
	if p["output"] == nothing
		printsummary(mortgagecalc(p["principal"], p["rate"], p["amortization"], p["frequency"], p["compounding"], p["startdate"], p["term"], p["payment"])...)
	else
		printtable(mortgagecalc(p["principal"], p["rate"], p["amortization"], p["frequency"], p["compounding"], p["startdate"], p["term"], p["payment"])...)
	end
end

main()


text:
struct System{N}
    matrix_entries::SparseMatrixCSC{GraphBasedSystems.Entry, Int64}
    vector_entries::Vector{Entry}
    diagonal_inverses::Vector{Entry}
    acyclic_children::Vector{Vector{Int64}} # Contains direct children that are not part of a cycle
    cyclic_children::Vector{Vector{Int64}}  # Contains direct and indirect children that are part of a cycle
    parents::Vector{Vector{Int64}}          # Contains direct and cycle-opening parents
    dfs_list::SVector{N,Int64}
    graph::SimpleGraph{Int64}
    dfs_graph::SimpleDiGraph{Int64}

    function System{T}(A, dims; force_static = false) where T
        N = length(dims)
        static = force_static || all(dims.<=10)

        full_graph = Graph(A)

        matrix_entries = spzeros(Entry,N,N)

        for (i,dimi) in enumerate(dims)
            for (j,dimj) in enumerate(dims)
                if i == j
                    matrix_entries[i,j] = Entry{T}(dimi, dimj, static = static)
                elseif j ∈ all_neighbors(full_graph,i)
                    matrix_entries[i,j] = Entry{T}(dimi, dimj, static = static)
                    matrix_entries[j,i] = Entry{T}(dimj, dimi, static = static)
                end
            end
        end

        vector_entries = [Entry{T}(dim, static = static) for dim in dims]
        diagonal_inverses = [Entry{T}(dim, dim, static = static) for dim in dims]

        graphs, roots = split_adjacency(A)
        dfs_list = Int64[]
        acyclic_children = [Int64[] for i=1:N]
        cycles = [Vector{Int64}[] for i=1:N]
        parents = [Int64[] for i=1:N]
        edgelist = LightGraphs.SimpleEdge{Int64}[]

        for (i,graph) in enumerate(graphs)
            root = roots[i]
            dfs_graph = dfs_tree(graph, root)
            sub_dfs_list, cycle_closures = dfs(graph, root)
            append!(dfs_list, sub_dfs_list)

            cycle_dfs_graph = copy(dfs_graph)
            cycle_dfs_graph_reverse = copy(dfs_graph)
            for cycle_closure in cycle_closures
                add_edge!(cycle_dfs_graph, cycle_closure...)
                add_edge!(cycle_dfs_graph_reverse, reverse(cycle_closure)...)
            end
            append!(edgelist, collect(edges(cycle_dfs_graph_reverse)))

            for v in sub_dfs_list
                acyclic_children[v] = neighbors(dfs_graph, v)
                parents[v] = neighbors(reverse(dfs_graph), v)
            end

            for cyclic_members in simplecycles(cycle_dfs_graph)
                v, cyclic_children = cycle_parent_children(cyclic_members, parents)
                cyclic_children = reverse(cyclic_children)
                lump_cycles!(cycles[v], cyclic_children)

                acyclic_children[v] = setdiff(acyclic_children[v], cyclic_children)
                for c in cyclic_children
                    matrix_entries[v,c] = Entry{T}(dims[v], dims[c], static = static)
                    matrix_entries[c,v] = Entry{T}(dims[c], dims[v], static = static)

                    v ∉ parents[c] && push!(parents[c],v)
                end
            end
        end

        full_dfs_graph = SimpleDiGraph(edgelist)
        cyclic_children = [unique(vcat(cycles[i]...)) for i=1:N]

        new{N}(matrix_entries, vector_entries, diagonal_inverses, acyclic_children, cyclic_children, parents, dfs_list, full_graph, full_dfs_graph)
    end

    System(A, dims; force_static = false) = System{Float64}(A, dims; force_static = force_static)
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, system::System{N}) where {N}
    println(io, "System with "*string(N)*" nodes.")
    SparseArrays._show_with_braille_patterns(io, system.matrix_entries)
end

@inline children(system, v) = outneighbors(system.dfs_graph, v)
@inline connections(system, v) = neighbors(system.graph, v)
@inline parents(system, v) = inneighbors(system.dfs_graph, v)

# There probably exists a smarter way of getting the dense matrix from the spares one
function full_matrix(system::System{N}) where N
    dims = [length(system.vector_entries[i].value) for i=1:N]

    range = [1:dims[1]]

    for (i,dim) in enumerate(collect(Iterators.rest(dims, 2)))
        push!(range,sum(dims[1:i])+1:sum(dims[1:i])+dim)
    end
    A = zeros(sum(dims),sum(dims))

    for (i,row) in enumerate(system.matrix_entries.rowval)
        col = findfirst(x->i<x,system.matrix_entries.colptr)-1
            A[range[row],range[col]] = system.matrix_entries[row,col].value
    end
    return A
end

full_vector(system) = vcat(getfield.(system.vector_entries,:value)...)


text:
#!/usr/bin/julia

# Trizen
# Date: 12 February 2017
# https://github.com/trizen

# Implementation of the sigma function (sum of the divisors).

# See also:
#    https://en.wikipedia.org/wiki/Divisor_function

using Primes

function divisor_sum(n, x=1)
    sigma = 1
    for (p,e) in factor(n)
        s = 1
        for j in 1:e
            s += p^(j*x)
        end
        sigma *= s
    end
    sigma
end

println(divisor_sum(1234))        # 1854
println(divisor_sum(99000))       # 365040
println(divisor_sum(99000, 2))    # 15359172920

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
using ParamLevelSet
using jInv.Mesh
using SparseArrays
using LinearAlgebra

# SpMatBuilder = ParamLevelSet.SpMatBuilder

# function getAlpha(k::Int64,m::Vector{Float64},numParamOfRBF::Int64 = 5)
# 	return m[numParamOfRBF*(k-1) + 1];
# end

# function getBeta(k::Int64,m::Vector{Float64})
# 	return m[5*(k-1) + 2];
# end

# function getX(k::Int64,m::Vector{Float64},numParamOfRBF::Int64 = 5)
# 	return m[numParamOfRBF*(k-1).+((numParamOfRBF-2):numParamOfRBF)];
# end

# function getXt(k::Int64,m::Vector{Float64},numParamOfRBF::Int64 = 5)
# 	return (m[numParamOfRBF*(k-1)+numParamOfRBF-2],m[numParamOfRBF*(k-1)+numParamOfRBF-1],m[numParamOfRBF*(k-1)+numParamOfRBF]);
# end

# function setX(k::Int64,m::Vector{Float64},x::Array)
# 	m[5*(k-1)+(3:5)] = x;
# end

# function computeFunc!(m::Vector{Float64},k::Int64,Xc,h::Vector{Float64},n::Vector{Int64},
#     u::Vector{Float64},xmin::Vector{Float64})
# (x1,x2,x3) = getXt(k,m);
# i = round.(Int64,(getX(k,m) - xmin)./h .+ 0.5);
# # ## here we take a reasonable box around the bf.
# # ## The constant depends on the type of RBF that we take.
# beta = getBeta(k,m);
# betaSq = beta^2;
# thres = 0.81/betaSq;
# boxL = ceil.(Int64,0.9./(h*abs(beta)));
# imax = min.(i + boxL,n);
# imin = max.(i - boxL,[1;1;1]);
# alpha = getAlpha(k,m);

# @inbounds for l = imin[3]:imax[3]
# @inbounds iishift3 = (l-1)*n[1]*n[2];
# @inbounds for j = imin[2]:imax[2]
#    @inbounds iishift2 = iishift3 + (j-1)*n[1];
#    @inbounds for q = imin[1]:imax[1]
#        ii = iishift2 + q;
#        @inbounds y = Xc[ii,1]-x1; y*=y;
#        nX = y;
#        @inbounds y = Xc[ii,2]-x2; y*=y;
#        nX+=y;
#        @inbounds y = Xc[ii,3]-x3; y*=y;
#        nX+=y;
#        if (nX <= thres) # 0.77~0.875^2
#            nX*=betaSq;
#            argii = radiust(nX);
#            argii = psi1(argii);
#            argii *= alpha;
#            @inbounds u[ii] += argii;
#        end
#    end
# end
# end
# end

# function updateJacobian!(k::Int64,m::Vector{Float64},dsu::Vector{Float64},h::Vector{Float64},n::Vector{Int64},Jbuilder::SpMatBuilder,u::Vector{Float64},iifunc::Function,Xc,xmin::Vector{Float64})
# (x1,x2,x3) = getXt(k,m);
# alpha = getAlpha(k,m);
# beta = getBeta(k,m);
# betaSQ = beta*beta;
# invBeta = (1.0/beta);
# thres = 0.81/betaSQ;

# alphaBetaSq = alpha*betaSQ;
# offset = convert(Int64,(k-1)*5 + 1);
# md = 1e-3*maximum(abs.(dsu));
# boxL = ceil.(Int64,0.9./(h*abs(beta)));
# i = round.(Int64,([x1;x2;x3] - xmin)./h .+ 0.5);
# imax = min.(i + boxL,n);
# imin = max.(i - boxL,[1;1;1]);

# @inbounds for l = imin[3]:imax[3]
# @inbounds iishift3 = (l-1)*n[1]*n[2];
# @inbounds for j = imin[2]:imax[2]
#    @inbounds iishift2 = iishift3 + (j-1)*n[1];
#    @inbounds for q = imin[1]:imax[1]
#        ii = iishift2 + q;
#        temp = dsu[ii];
#        if temp >= md
#            @inbounds y1 = x1 - Xc[ii,1]; nX =y1*y1;
#            @inbounds y2 = x2 - Xc[ii,2]; nX+=y2*y2;
#            @inbounds y3 = x3 - Xc[ii,3]; nX+=y3*y3;
#            if (nX <= thres) # 0.77~0.875^2
#                radii = radiust(nX*betaSQ);
#                psi,dpsi = dpsi1_t(radii);
#                psi*=temp;
#                temp*= alphaBetaSq;
#                temp/= radii;
#                temp*=dpsi;
#                nX *= temp;
#                nX *= invBeta;
#                y1*=temp; y2*=temp; y3*=temp;
#                ii = iifunc(ii);
#                setNext!(Jbuilder,ii,offset,psi);
#                setNext!(Jbuilder,ii,offset+1,nX);
#                setNext!(Jbuilder,ii,offset+2,y1);
#                setNext!(Jbuilder,ii,offset+3,y2);
#                setNext!(Jbuilder,ii,offset+4,y3);
#            end
#        end
#    end
# end
# end
# end

# function MeshFreeParamLevelSetModelFunc(Mesh::RegularMesh, m::Vector; computeJacobian = 1, bf::Int64 = 1,
#     sigma::Function = (m,n)->(n[:] .= 1.0),
#     Xc = convert(Array{Float64,2}, getCellCenteredGrid(Mesh)),u::Vector = zeros(size(Xc,1)),dsu::Vector = zeros(size(Xc,1)),
#     Jbuilder::SpMatBuilder{Int64,Float64} = getSpMatBuilder(Int64,Float64,size(Xc,1),length(m),computeJacobian*10*size(Xc,1)),
#     iifunc::Function=identity, numParamOfRBF::Int64 = 5,tree = BallTree(Matrix(Xc')))

#     if length(u) != size(Xc,1)
#         error("preallocated u is of wrong size");
#     end

#     u[:] .= 0.0;
#     reset!(Jbuilder);
#     nRBFs = div(length(m),numParamOfRBF);
#     ## in run 1 we calculate u. in run 2 we calculate J after we know the derivative of the heaviside.
#     beta_arr = Array{Float64}(undef,nRBFs);

#    for k=1:nRBFs
#        MeshFreeComputeFunc!(m,k,u,Xc,tree);
#    end

#    sigma(u,dsu);
#     if computeJacobian == 1
#        for k=1:nRBFs
#            MeshFreeUpdateJacobian!(k,m,dsu,Jbuilder,u,iifunc,Xc,tree);
#        end
#     end
#     return u,Jbuilder;
# end

# function eval_f_grad()
#     sigma::Function = (m,n)->(n[:] .= 1.0),
#     Xc = convert(Array{Float64,2}, getCellCenteredGrid(Mesh)),u::Vector = zeros(size(Xc,1)),dsu::Vector = zeros(size(Xc,1)),
#     # Jbuilder::SpMatBuilder{Int64,Float64} = getSpMatBuilder(Int64,Float64,size(Xc,1),length(m),computeJacobian*10*size(Xc,1)),
#     # iifunc::Function=identity, numParamOfRBF::Int64 = 5,tree = BallTree(Matrix(Xc')))

#     if length(u) != size(Xc,1)
#         error("preallocated u is of wrong size");
#     end

#     u[:] .= 0.0;
#     reset!(Jbuilder);
#     nRBFs = div(length(m),numParamOfRBF);
#     ## in run 1 we calculate u. in run 2 we calculate J after we know the derivative of the heaviside.
#     beta_arr = Array{Float64}(undef,nRBFs);

#    for k=1:nRBFs
#        MeshFreeComputeFunc!(m,k,u,Xc,tree);
#    end

#    sigma(u,dsu);
#     if computeJacobian == 1
#        for k=1:nRBFs
#            MeshFreeUpdateJacobian!(k,m,dsu,Jbuilder,u,iifunc,Xc,tree);
#        end
#     end
#     return u,Jbuilder;
# end

function init_pals(vol_size, nbasis=4)
	half = vol_size ./ 2
    # Mesh = getRegularMesh([-half[1];half[1];-half[2];half[2];-half[3];half[3]],[vol_size[1],vol_size[2],vol_size[3]]);
    Mesh = getRegularMesh([1;vol_size[1];1;vol_size[2];1;vol_size[3]],[vol_size[1],vol_size[2],vol_size[3]]);

	alpha = [1.5;2.5;-2.0;-1.0];
	beta = [2.5;2.0;-1.5;2.5]*0.001
	Xs = [0.5 0.5 0.5; 2.0 2.0 2.0; 1.2 2.3 1.5; 2.2 1.5 2.0] .* vol_size[1] .* 0.3
	m = wrapTheta(alpha,beta,Xs);

	return Mesh, m
end

function recon2d_slices_pals!(u::Array{T, 3}, A, p, niter::Int=1, nbasis::Int=4) where {T <: AbstractFloat}

    global Jt, J0, mesh, m, Hessian, grad, r

    mesh, m = init_pals(size(u), nbasis)
    sigmaH = getDefaultHeaviside();
    phat = similar(p)
    r = similar(p)
    At = A'
    stepsize = 0.01
    nslice = size(p,2)

    HW = size(u,1)*size(u,2)
    szpslice = size(p,1)*size(p,3)
    J = spzeros(prod(size(p)), nbasis*5) # A * J0

    for i=1:niter
        # u0, JBuilder = eval_f_grad(mesh, m; computeJacobian=1, sigma=sigmaH)
        u0, JBuilder = MeshFreeParamLevelSetModelFunc(mesh, m; computeJacobian = 1)

        J0 = getSparseMatrix(JBuilder)

        # construct J
        for ss=1:nslice
            J[szpslice*(ss-1)+1:szpslice*ss, :] .= A * J0[HW*(ss-1)+1:HW*ss, :]

            u_slice = view(u, :, :, ss)
            Av = A*vec(u_slice)

            # phat_slice = reshape(Av, size(p,1), size(p,3))
            p_vec_ = reshape(Av, size(p,1), size(p,3))
            @views r[:,ss,:] = p_vec_ - p[:,ss,:]
        end

        Jt = J'
        Hessian = Jt*J
        Hinv = pinv(Matrix(Hessian))
        grad = Jt * vec(r)
        m .-= stepsize * Hinv * grad

        fres = sum(r .^ 2) / length(r)
        @show fres, m
        # println("residual: $fres, $m")
    end

    # u0, JBuilder = eval_f_grad(mesh,m;computeJacobian=1, sigma=sigmaH)
    u0, JBuilder = MeshFreeParamLevelSetModelFunc(mesh, m; computeJacobian = 1)
    u .= reshape(u0, size(u))
    return u, m
end

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
# Basic arithmetic
Base.:*(a::Number, x::SparseArray) =
    mul!(similar(x, Base.promote_eltypeof(a, x)), a, x)
Base.:*(x::SparseArray, a::Number) =
    mul!(similar(x, Base.promote_eltypeof(a, x)), x, a)
Base.:\(a::Number, x::SparseArray) =
    mul!(similar(x, Base.promote_eltypeof(a, x)), inv(a), x)
Base.:/(x::SparseArray, a::Number) =
    mul!(similar(x, Base.promote_eltypeof(a, x)), x, inv(a))
Base.:+(x::SparseArray, y::SparseArray) =
    (T = Base.promote_eltypeof(x, y); axpy!(+one(T), y, copy!(similar(x, T), x)))
Base.:-(x::SparseArray, y::SparseArray) =
    (T = Base.promote_eltypeof(x, y); axpy!(-one(T), y, copy!(similar(x, T), x)))

Base.:-(x::SparseArray) = LinearAlgebra.lmul!(-one(eltype(x)), copy(x))

Base.zero(x::SparseArray) = similar(x)
Base.iszero(x::SparseArray) = nonzero_length(x) == 0

function Base.one(x::SparseArray{<:Any,2})
    m, n = size(x)
    m == n ||
        throw(DimensionMismatch("multiplicative identity defined only for square matrices"))

    u = similar(x)
    @inbounds for i = 1:m
        u[i,i] = one(eltype(x))
    end
    return u
end

# comparison
function Base.:(==)(x::SparseArray, y::SparseArray)
    keys = collect(nonzero_keys(x))
    intersect!(keys, nonzero_keys(y))
    if !(length(keys) == length(nonzero_keys(x)) == length(nonzero_keys(y)))
        return false
    end
    for I in keys
        x[I] == y[I] || return false
    end
    return true
end

# Vector space functions
#------------------------
function Base.conj!(x::SparseArray)
    conj!(x.data.vals)
    return x
end
function LinearAlgebra.mul!(dst::SparseArray, a::Number, src::SparseArray)
    _zero!(dst)
    for (k, v) in nonzero_pairs(src)
        dst[k] = a*v
    end
    return dst
end
function LinearAlgebra.mul!(dst::SparseArray, src::SparseArray, a::Number)
    _zero!(dst)
    for (k, v) in nonzero_pairs(src)
        dst[k] = v*a
    end
    return dst
end
function LinearAlgebra.lmul!(a::Number, x::SparseArray)
    lmul!(a, x.data.vals)
    # typical occupation in a dict is about 30% from experimental testing
    # the benefits of scaling all values (e.g. SIMD) largely outweight the extra work
    return x
end
function LinearAlgebra.rmul!(x::SparseArray, a::Number)
    rmul!(x.data.vals, a)
    return x
end
function LinearAlgebra.ldiv!(a::Number, x::SparseArray)
    ldiv!(a, x.data.vals)
    return x
end
function LinearAlgebra.rdiv!(x::SparseArray, a::Number)
    rdiv!(x.data.vals, a)
    return x
end
function LinearAlgebra.axpby!(α, x::SparseArray, β, y::SparseArray)
    β == one(β) || (iszero(β) ? _zero!(y) : LinearAlgebra.lmul!(β, y))
    for (k, v) in nonzero_pairs(x)
        increaseindex!(y, α*v, k)
    end
    return y
end
function LinearAlgebra.axpy!(α, x::SparseArray, y::SparseArray)
    for (k, v) in nonzero_pairs(x)
        increaseindex!(y, α*v, k)
    end
    return y
end

function LinearAlgebra.norm(x::SparseArray, p::Real = 2)
    norm(nonzero_values(x), p)
end

function LinearAlgebra.dot(x::SparseArray, y::SparseArray)
    size(x) == size(y) || throw(DimensionMismatch("dot arguments have different size"))
    s = dot(zero(eltype(x)), zero(eltype(y)))
    if nonzero_length(x) >= nonzero_length(y)
        @inbounds for I in nonzero_keys(x)
            s += dot(x[I], y[I])
        end
    else
        @inbounds for I in nonzero_keys(y)
            s += dot(x[I], y[I])
        end
    end
    return s
end

# permutedims
Base.permutedims!(dst::SparseArray, src::SparseArray, p) =
    add!(one(eltype(dst)), src, :N, zero(eltype(dst)), dst, tuple(p...))

# matrix functions
const SV{T} = SparseArray{T,1}
const SM{T} = SparseArray{T, 2}
const ASM{T} = Union{SparseArray{T, 2},
                    Transpose{T, <:SparseArray{T,2}},
                    Adjoint{T, <:SparseArray{T,2}}}

LinearAlgebra.mul!(C::SM, A::ASM, B::ASM) = mul!(C, A, B, one(eltype(C)), zero(eltype(C)))
function LinearAlgebra.mul!(C::SM, A::ASM, B::ASM, α::Number, β::Number)
    CA = A isa Adjoint ? :C : :N
    CB = B isa Adjoint ? :C : :N
    oindA = A isa Union{Adjoint,Transpose} ? (2,) : (1,)
    cindA = A isa Union{Adjoint,Transpose} ? (1,) : (2,)
    oindB = B isa Union{Adjoint,Transpose} ? (1,) : (2,)
    cindB = B isa Union{Adjoint,Transpose} ? (2,) : (1,)

    AA = A isa Union{Adjoint,Transpose} ? parent(A) : A
    BB = B isa Union{Adjoint,Transpose} ? parent(B) : B

    contract!(α, AA, CA, BB, CB, β, C, oindA, cindA, oindB, cindB, (1, 2))
end

function LinearAlgebra.adjoint!(C::SM, A::SM)
    add!(one(eltype(C)), A, :C, zero(eltype(C)), C, (2, 1))
    return C
end
function LinearAlgebra.transpose!(C::SM, A::SM)
    add!(one(eltype(C)), A, :N, zero(eltype(C)), C, (2, 1))
    return C
end


text:
struct aType
  a::Integer
end

function (::Type{aType})(arg1::Int, arg2::Int)
  println("thang")
end

outputty = aType(1, 3)

println(outputty.a)

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
using ..ChemistryFeaturization.AbstractType: AbstractCodec
import ..ChemistryFeaturization.encode
import ..ChemistryFeaturization.decode
using Flux: onecold

export encode, decode, output_shape

"""
    OneHotOneCold(encode_f, decode_f, nbins, logspaced)

AbstractCodec type which uses a dummy variable (as defined in statistical literature), i.e., which employs
one-hot encoding and a one-cold decoding scheme.
"""
struct OneHotOneCold <: AbstractCodec
    categorical::Bool
    bins::Vector
end

output_shape(ohoc::OneHotOneCold) =
    ohoc.categorical ? length(ohoc.bins) : length(ohoc.bins) - 1

"A flexible version of Flux.onehot that can handle both categorical and continuous-valued encoding."
function encode(ohoc::OneHotOneCold, val)
    local bin_index, onehot_vec
    if ohoc.categorical
        onehot_vec = [0.0 for i = 1:length(ohoc.bins)]
        bin_index = findfirst(isequal(val), ohoc.bins)
    else
        @assert eltype(ohoc.bins) <: Number "Your bins aren't numbers...are you sure you didn't mean for this feature to be categorical?"
        @assert ohoc.bins[1] <= val <= ohoc.bins[end] "The value $val is outside the range of bins $ohoc.bins"
        onehot_vec = [0.0 for i = 1:(length(ohoc.bins)-1)]
        bin_index = searchsorted(ohoc.bins, val).stop
        if bin_index == length(ohoc.bins) # got the max value
            bin_index = bin_index - 1
        elseif isapprox(val, ohoc.bins[1]) # sometimes will get 0 if this doesn't get checked
            bin_index = 1
        end
    end
    onehot_vec[bin_index] = 1.0
    return onehot_vec
end

encode(ohoc::OneHotOneCold, vals::Vector) = encode.(Ref(ohoc), vals)

function decode(ohoc::OneHotOneCold, encoded)
    @assert size(encoded)[1] == output_shape(ohoc)
    local decoded
    if ohoc.categorical # return value
        decoded = onecold(encoded, ohoc.bins)
    else # return bounds (TODO: should this be a tuple or a vector..? I like tuple for distinguishing from encoded vectors, but it doesn't play so nice with broadcasting...)
        decoded = (onecold(encoded, ohoc.bins[1:end-1]), onecold(encoded, ohoc.bins[2:end]))
    end
    return decoded
end


text:
# all functions must receive client and message arguments!
using Discord

const DATA = (
    userpoints => Dict{"String", "Int"}(),
    usermoney => Dict{"String", "Int"}(),
    currentusersworking => String[]
) # not too certain how to store and load data. Databases?

function replypong(c::Client, m::Message)
    reply(c, m, "pong")
end

function replyping(c::Client, m::Message)
    reply(c, m, "ping")
end

function replytopleaderboard(c::Client, m::Message)
    # Unfold data and create Embed to show rich content
    # https://posts.xh4h.com/Discord.jl/dev/types/#Discord.Embed
    # Some tools
    # https://leovoel.github.io/embed-visualizer/
    # https://cog-creators.github.io/discord-embed-sandbox/
    reply(c, m, DATA.userpoints)
end

function workcontroller(c::Client, m::Message)
    # check if user has been working
    # if so, stop work
    # else, start work
end

max_stars_repo_name: dzonimn/HoJDiscordBot

max_issues_repo_name: dzonimn/HoJDiscordBot

max_forks_repo_name: dzonimn/HoJDiscordBot


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

using Markdown
using InteractiveUtils

# ╔═╡ 33616748-f796-11ea-1c1d-2dae00de8c6e
using Pkg, DrWatson

# ╔═╡ 33619e36-f796-11ea-0b36-bda68b9a5e2f
begin
	@quickactivate "StatisticalRethinkingStan"
	using Random, DynamicHMC, Distributions
	using TransformVariables, LogDensityProblems
	using StatisticalRethinking
end

# ╔═╡ 97cf3b66-f795-11ea-18b2-d9e9689a0b20
md"# Intro-dhmc-1d.jl"

# ╔═╡ 336211e8-f796-11ea-3d91-fbe93bfeafb3
md"### snippet 9.3"

# ╔═╡ 33705476-f796-11ea-0972-edddcf5fc6de
md"##### Construct the logdensity problem."

# ╔═╡ 3370fc94-f796-11ea-02d3-1f9c5c37c18f
struct clip_9_3_model{TY <: AbstractVector, TX <: AbstractVector}
    "Observations."
    y::TY
    "Covariate"
    x::TX
end

# ╔═╡ 337c7ccc-f796-11ea-372d-03b0e499af83
md"##### Make the type callable with the parameters *as a single argument*."

# ╔═╡ 337d0f3e-f796-11ea-1afe-978a61c7cf51
function (problem:: clip_9_3_model)(θ)
    @unpack y, x, = problem    # extract the data
    @unpack muy, mux = θ     # works on the named tuple too
    ll = 0.0
    ll += loglikelihood(Normal(mux, 1), x)
    ll += loglikelihood(Normal(muy, 1), y)
    ll += logpdf(Normal(0, 1), mux)
    ll += logpdf(Normal(0, 1), muy)
    ll
end

# ╔═╡ 3387db1c-f796-11ea-3c9b-912fb0729b97
md"##### Instantiate the model with data and inits."

# ╔═╡ 33888ca4-f796-11ea-3137-cf3d79ddeb9b
Random.seed!(1234591);

# ╔═╡ c40363d8-f7c2-11ea-33ab-bdf7b5f51496
d = rand(Normal(0,1), 10)

# ╔═╡ 3394ced0-f796-11ea-026a-c5e2f92cd86c
begin
	N = 100
	x = rand(Normal(0, 1), N)
	y = rand(Normal(0, 1), N)

	p = clip_9_3_model(y, x)
	θ = (muy = 0.0, mux=0.0)
	p(θ)
end;

# ╔═╡ 339bdb26-f796-11ea-263a-c52b4b4acce2
md"##### Write a function to return properly dimensioned transformation."

# ╔═╡ 339d8368-f796-11ea-28d5-cf374f141754
problem_transformation(p::clip_9_3_model) =
    as((muy = asℝ, mux = asℝ));

# ╔═╡ 33aa58b8-f796-11ea-0f79-173450a7276c
md"##### Wrap the problem with a transformation, then use Flux for the gradient."

# ╔═╡ 33abb8a2-f796-11ea-3bc9-319850bbeb9b
begin
	P = TransformedLogDensity(problem_transformation(p), p)
	∇P = ADgradient(:ForwardDiff, P);
end

# ╔═╡ 33b2a428-f796-11ea-01e4-93e8d05c759d
md"##### Generate and show fig 9.3."

# ╔═╡ 33bb5ece-f796-11ea-1ed2-1f9c93c1c9ac
function generate_n_samples(model, grad;
  epsilon = 0.03, # Step size
  L = 11, # No of leapfrog steps
  pr = 0.3, # Plot range
  n_samples = 4, # No of samples
  q = [-0.1, 0.2]) # Initial position

  fig = plot(ylab="muy", xlab="mux", xlim=(-pr, pr), ylim=(-pr, pr), leg=false)

  for i in 1:n_samples
    q, ptraj, qtraj, accept, dH = StatisticalRethinking.HMC(model, grad, 0.03, 11, q)
    if n_samples < 10
      for j in 1:L
        K0 = sum(ptraj[j, :] .^ 2) / 2
        plot!(fig, [qtraj[j, 1], qtraj[j+1, 1]], [qtraj[j, 2], qtraj[j+1, 2]], leg=false)
        scatter!(fig, [(qtraj[j, 1], qtraj[j, 2])])
      end
    end
  end
  fig
end

# ╔═╡ 33c2b99e-f796-11ea-0fde-89c09665dda0
fig = generate_n_samples(p, ∇P; pr=0.5)

# ╔═╡ 33cbd4ac-f796-11ea-2dbe-510c35a57460
md"## End of clip-09-03d.jl"

# ╔═╡ Cell order:
# ╟─97cf3b66-f795-11ea-18b2-d9e9689a0b20
# ╠═33616748-f796-11ea-1c1d-2dae00de8c6e
# ╠═33619e36-f796-11ea-0b36-bda68b9a5e2f
# ╟─336211e8-f796-11ea-3d91-fbe93bfeafb3
# ╟─33705476-f796-11ea-0972-edddcf5fc6de
# ╠═3370fc94-f796-11ea-02d3-1f9c5c37c18f
# ╟─337c7ccc-f796-11ea-372d-03b0e499af83
# ╠═337d0f3e-f796-11ea-1afe-978a61c7cf51
# ╟─3387db1c-f796-11ea-3c9b-912fb0729b97
# ╠═33888ca4-f796-11ea-3137-cf3d79ddeb9b
# ╠═c40363d8-f7c2-11ea-33ab-bdf7b5f51496
# ╠═3394ced0-f796-11ea-026a-c5e2f92cd86c
# ╟─339bdb26-f796-11ea-263a-c52b4b4acce2
# ╠═339d8368-f796-11ea-28d5-cf374f141754
# ╟─33aa58b8-f796-11ea-0f79-173450a7276c
# ╠═33abb8a2-f796-11ea-3bc9-319850bbeb9b
# ╟─33b2a428-f796-11ea-01e4-93e8d05c759d
# ╠═33bb5ece-f796-11ea-1ed2-1f9c93c1c9ac
# ╠═33c2b99e-f796-11ea-0fde-89c09665dda0
# ╟─33cbd4ac-f796-11ea-2dbe-510c35a57460


text:
module MessageRequest

export body_is_a_stream, body_was_streamed, setuseragent!, resource

import ..Layer, ..request
using ..IOExtras
using URIs
using ..Messages
import ..Messages: bodylength
import ..Headers
import ..Form, ..content_type

"""
"request-target" per https://tools.ietf.org/html/rfc7230#section-5.3
"""
resource(uri::URI) = string( isempty(uri.path)     ? "/" :     uri.path,
                            !isempty(uri.query)    ? "?" : "", uri.query,
                            !isempty(uri.fragment) ? "#" : "", uri.fragment)

"""
    request(MessageLayer, method, ::URI, headers, body) -> HTTP.Response

Construct a [`Request`](@ref) object and set mandatory headers.
"""
struct MessageLayer{Next <: Layer} <: Layer{Next} end
export MessageLayer

function request(::Type{MessageLayer{Next}},
                 method::String, url::URI, headers::Headers, body;
                 http_version=v"1.1",
                 target=resource(url),
                 parent=nothing, iofunction=nothing, kw...) where Next

    defaultheader!(headers, "Host" => url.host)
    if USER_AGENT[] !== nothing
        defaultheader!(headers, "User-Agent" => USER_AGENT[])
    end

    if !hasheader(headers, "Content-Length") &&
       !hasheader(headers, "Transfer-Encoding") &&
       !hasheader(headers, "Upgrade")
        l = bodylength(body)
        if l != unknown_length
            setheader(headers, "Content-Length" => string(l))
        elseif method == "GET" && iofunction isa Function
            setheader(headers, "Content-Length" => "0")
        end
    end
    if !hasheader(headers, "Content-Type") && body isa Form && method == "POST"
        # "Content-Type" => "multipart/form-data; boundary=..."
        setheader(headers, content_type(body))
    end

    req = Request(method, target, headers, bodybytes(body);
                  parent=parent, version=http_version)

    return request(Next, url, req, body; iofunction=iofunction, kw...)
end

const USER_AGENT = Ref{Union{String, Nothing}}("HTTP.jl/$VERSION")

"""
    setuseragent!(x::Union{String, Nothing})

Set the default User-Agent string to be used in each HTTP request.
Can be manually overridden by passing an explicit `User-Agent` header.
Setting `nothing` will prevent the default `User-Agent` header from being passed.
"""
function setuseragent!(x::Union{String, Nothing})
    USER_AGENT[] = x
    return
end

bodylength(body) = unknown_length
bodylength(body::AbstractVector{UInt8}) = length(body)
bodylength(body::AbstractString) = sizeof(body)
bodylength(body::Form) = length(body)
bodylength(body::Vector{T}) where T <: AbstractString = sum(sizeof, body)
bodylength(body::Vector{T}) where T <: AbstractArray{UInt8,1} = sum(length, body)
bodylength(body::IOBuffer) = bytesavailable(body)
bodylength(body::Vector{IOBuffer}) = sum(bytesavailable, body)

const body_is_a_stream = UInt8[]
const body_was_streamed = bytes("[Message Body was streamed]")
bodybytes(body) = body_is_a_stream
bodybytes(body::Vector{UInt8}) = body
bodybytes(body::IOBuffer) = read(body)
bodybytes(body::AbstractVector{UInt8}) = Vector{UInt8}(body)
bodybytes(body::AbstractString) = bytes(body)
bodybytes(body::Vector) = length(body) == 1 ? bodybytes(body[1]) :
                                              body_is_a_stream

end # module MessageRequest


text:
cd(raw"C:\Users\francis.smart.ctr\RandomnessAndTalent.jl")
Pkg.activate(".")

module RandomnessAndTalent

using DataFrames, Gadfly

xgrid = 200
ygrid = 200

ni = 1000
ix = rand(1:xgrid, ni)
iy = rand(1:ygrid, ni)

talentave = .6
talentsd  = .1

italent = max.(0, min.(1, talentave .+ talentsd.*rand(Normal(), ni)))

maxt = round(maximum(italent), digits=2)
mint = round(minimum(italent), digits=2)

print("Maximum Talent:$maxt \nMinimum Talent:$mint")

ngood = 250

goodx = rand(1:xgrid, ngood)
goody = rand(1:ygrid, ngood)

nbad  = 250
badx = rand(1:xgrid, ngood)
bady = rand(1:ygrid, ngood)

timeperiods    = 80
initialcapital = 100.0

people = DataFrame(x = ix,    y = iy,    type = "individual")
good   = DataFrame(x = goodx, y = goody, type = "good")
bad    = DataFrame(x = badx,  y = bady,  type = "bad")

gridset = vcat(people, good, bad)

set_default_plot_size(20cm, 16cm)  # hide

plot(gridset, x = "x", y = "y", color = "type", Geom.point, Theme(point_size = 2pt),
     Scale.x_continuous(minvalue=0, maxvalue=xgrid),
     Scale.y_continuous(minvalue=0, maxvalue=ygrid))

function randomwalk!(dataframe, rstep=1, xborder = xgrid, yborder = ygrid)
  roundboundary(xmax, x) = round.(max.(0, min.(x, xmax)))

  dataframe.x =  roundboundary(xborder, dataframe.x + rstep.*rand(Normal(), size(dataframe, 1)))
  dataframe.y =  roundboundary(yborder, dataframe.y + rstep.*rand(Normal(), size(dataframe, 1)))

end

function nextphase(dataframe, good = dataframe(x=-1, y=-1), bad = dataframe(x=-1, y=-1), verbose = false)
  t = maximum(dataframe.t)
  recentdata = dataframe[dataframe.t .== t, :]
  recentdata.t .+= 1
  recentdata.good .= 0
  recentdata.realised .= 0
  recentdata.bad .= 0

  #Loop through each individual and check if a good event happened
  for i in 1:size(recentdata, 1)
    # Check if good event could happen
    if sum((good.x .== recentdata.x[i]) .& (good.y .== recentdata.y[i])) > 0
        verbose && println("$i possibly get fortuitous event")
        recentdata.good[i] = 1
        if rand(1,1)[1] < recentdata.talent[i]
            recentdata.capital[i] *= 2
            recentdata.realised[i] = 1
            verbose && println("$i get fortuitous event")

        end
    end

    # Check if bad even happens
    if sum((bad.x .== recentdata.x[i]) .& (bad.y .== recentdata.y[i])) > 0
        verbose && println("$i is unlucky")
        recentdata.capital[i] *= .5
        recentdata.bad[i] = 1
    end
  end

  dataframe = vcat(dataframe, recentdata)
end # module

peopleset = DataFrame(i = 1:ni, t = 0, capital = initialcapital, talent = italent,
  x = ix, y = iy, good=0, realised=0, bad=0)

for i in 1:timeperiods;
    global peopleset
    peopleset = nextphase(peopleset, good, bad);
    randomwalk!(good)
    randomwalk!(bad)
end

retirement = peopleset[peopleset.t .== timeperiods,:]

top = retirement[retirement.capital .> 1000, :].i

plot(peopleset[[i in top for i in peopleset.i],:], x = "t", y = "capital", color = "i",
  Geom.point, Geom.line, Scale.y_log)

plot(peopleset[[i in top for i in peopleset.i],:], x = "t", y = "capital", color = "i",
  Geom.point, Geom.line)

peopleset[[i in top for i in peopleset.i] .& (peopleset.t .== 0), :talent]

plot(peopleset[peopleset.t .== 80,:], x = "talent", y = "capital", Geom.point, Scale.y_log)

eventcounts = by(peopleset, :i,
  (:good, :bad, :realised) => x -> (good = sum(x[1]), bad = sum(x[2]), realised = sum(x[3])))

peopleset[peopleset.good .== 1, :]

mean(eventcounts.good)
mean(eventcounts.bad)
mean(eventcounts.realised)

end #module


text:
module JuliaInterface

# include and implement the interface of this module
include("interface/IJuliaInterface.jl")

function getResults(aStruct::OneStruct)
    Addition.doAdd(aStruct.a, aStruct.b)
end

end

max_stars_repo_name: sonolac/JuliaInterface

max_issues_repo_name: sonolac/JuliaInterface

max_forks_repo_name: sonolac/JuliaInterface


text:
# ---
# title: 295. Find Median from Data Stream
# id: problem295
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Heap, Design
# link: <https://leetcode.com/problems/find-median-from-data-stream/description/>
# hidden: true
# ---
#
# Median is the middle value in an ordered integer list. If the size of the list
# is even, there is no middle value. So the median is the mean of the two middle
# value.
#
# For example,
#
# `[2,3,4]`, the median is `3`
#
# `[2,3]`, the median is `(2 + 3) / 2 = 2.5`
#
# Design a data structure that supports the following two operations:
#
#   * void addNum(int num) - Add a integer number from the data stream to the data structure.
#   * double findMedian() - Return the median of all elements so far.
#
#
#
# **Example:**
#
#
#
#     addNum(1)
#     addNum(2)
#     findMedian() -> 1.5
#     addNum(3)
#     findMedian() -> 2
#
#
#
#
# **Follow up:**
#
#   1. If all integer numbers from the stream are between 0 and 100, how would you optimize it?
#   2. If 99% of all integer numbers from the stream are between 0 and 100, how would you optimize it?
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Implement code blocks 3.20 through 3.26 in Statistical Rethinking, 2nd edition.

using Random
using Distributions: Binomial, pdf, rand
using StatsBase: countmap, proportionmap, sample, Weights
using Plots: histogram

Random.seed!(17)

# Code block 3.20
distro = Binomial(2, 0.7)
probas = pdf(distro)

# Code block 3.21
rand(distro, 1)

# Code block 3.22
rand(distro, 10)

# Code block 3.23
draw = rand(distro, 100_000)
countmap(draw)
proportionmap(draw)

# Code block 3.24
distro = Binomial(9, 0.7)
draw = rand(distro, 100_000)
display(histogram(draw, xlabel="Dummy water count", ylabel="Frequency"))

# Code block 3.25
distro = Binomial(9, 0.6)
draw = rand(distro, 10_000)

# Code block 3.26 (also block 3.27 in the practice questions)

# This one is a little confusing - it's not clear what `samples` refers to in the book's
# code snippet. I *think* these are samples from the posterior distribution of the
# parameter P(theta | data), and we use these to then get the posterior predictive
# distribution P(y | data).
# See snippets 3.11 and 3.27 as examples.
grid_size = 1_000
num_samples = 10_000

param_grid = range(0, 1, length=grid_size)
prior = ones((grid_size))
likelihood = [pdf(Binomial(9, p), 6) for p in param_grid]
posterior = (likelihood .* prior)
posterior = posterior ./ sum(posterior)

param_samples = sample(
    param_grid, Weights(posterior), num_samples; replace=true, ordered=false
)

posterior_draw = [rand(Binomial(9, p)) for p in param_samples]  # this seems inefficient

# This figure is the bottom plot in Figure 3.6.
display(
    histogram(
        posterior_draw, xlabel="Water obs count", ylabel="Posterior Predictve Frequency"
    )
)


text:
using PooksoftOptionsKit
using PyPlot
using Dates

# setup the call test -
stock_price = 54.17
implied_volatility = 0.67
DTE = (31.0/365.0)
number_of_levels = 4
risk_free_rate = 0.15 # this is in percent
dividend_rate = 0.0 # this is in percent

# setup options calculation -
amd_option_parameters = PSOptionKitPricingParameters(implied_volatility, DTE, number_of_levels, risk_free_rate, dividend_rate)

# setup an amd call option -
assetSet = Set{PSAbstractAsset}()
callOptionContract = PSCallOptionContract("AMD", Date(2020,6,25), 60.0, 0.0, 1; sense=:buy, contractMultiplier=1.0)
push!(assetSet, callOptionContract)

# build the pricing tree -
tree = build_binary_price_tree(assetSet, amd_option_parameters, stock_price)

# compute -
result = option_contract_price(tree, amd_option_parameters)


text:
function Measure_Fluxes!(highway, t, sections, T, flux_local0, flux_local1, flux_local2)

    for (i, section) in enumerate(sections)
      for j = section+4:-1:section
				if highway[j].speed > 0
            add_Fluxes!(flux_local0, highway[j], i, t, section, T, Int8(0))
            add_Fluxes!(flux_local1, highway[j], i, t, section, T, Int8(1))
            add_Fluxes!(flux_local2, highway[j], i, t, section, T)
					end
      end
    end
end

function Measure_Speeds!(highway, t, sections, T, speed_local0, speed_local1, speed_local2)

    for (i, section) in enumerate(sections)
      for j = section+4:-1:section
				if highway[j].speed > 0
            add_Speeds!(speed_local0[i, t], highway[j], section, Int8(0))
            add_Speeds!(speed_local1[i, t], highway[j], section, Int8(1))
            add_Speeds!(speed_local2[i, t], highway[j], section)
					end
      end
    end
end

function add_Fluxes!(flux, v, i, t, x_section, T, tipo::Int8 = Int8(2))
  if v.position-v.speed < x_section
		if tipo == Int8(0)
			if v.tipo > Int8(0)
				flux[i, t] += 1/T
			end
		else
			if v.tipo == tipo
				flux[i, t] += 1/T
			end
		end
  end
end

function add_Speeds!(array, v, x_section, tipo::Int8 = Int8(2))

	if v.position-v.speed < x_section
		## discriminative count
    if v.tipo == tipo #tipo != Int8(0) &&
			push!(array, v.speed)
    end

		## tipo == 0 to count vehicles no matter the kind
    if tipo == Int8(0) && v.tipo > Int8(0)
			push!(array, v.speed)
    end
	end
end

function Measure_Densities!(mean_density, D_mean_density, mean_speed, D_mean_speed, mean_flux, D_mean_flux)

  idx_s = find(mean_speed .== 0)
  idx_f = find(mean_flux .== 0)

  # Normal case
  mean_density = mean_flux./mean_speed
  D_mean_density = (D_mean_flux.*mean_speed - mean_flux.*D_mean_speed)./(mean_speed.^2)

  # Vehicles present but no flux case
  mean_density[idx_s] = 1
  D_mean_density[idx_s] = 0

  # No vehicles case
  mean_density[intersect(idx_s, idx_f)] = 0
  D_mean_density[intersect(idx_s, idx_f)] = 0
end

function Measure!(Lanes::Int64, S::Int64, N::Int64, A::Int64, mean_speed, mean_flux,
									D_mean_speed, D_mean_flux, speeds, fluxes; sense::Int64 = 1)

	mean_speed_n = zeros(S, Lanes, A)
	mean_flux_n = zeros(S, Lanes, A)

	# ensamble averaging
	for t = 1:A, k = 1:Lanes, s = 1:S
    K = (sense == 1 ? k : k+1)
		idx_n = find([f for f in fluxes[s, K, t, :]])
    mean_speed_n[s, k, t] = (length(idx_n) != 0 ? mean(speeds[s, k, t, idx_n]) : 0)
    mean_flux_n[s, k, t] = (length(idx_n) != 0 ? mean(fluxes[s, K, t, idx_n]) : 0)
	end

	# temporal averaging
	for k = 1:Lanes, s = 1:S
		idx_t = find([f for f in mean_flux_n[s, k, :]])

		mean_flux[s, k] = (length(idx_t) != 0 ? mean(mean_flux_n[s, k, idx_t]) : 0)
		mean_speed[s, k] = (length(idx_t) != 0 ? mean(mean_speed_n[s, k, idx_t]) : 0)

		D_mean_flux[s, k] = (length(idx_t) > 1 ? std(mean_flux_n[s, k, idx_t]) : 0)
		D_mean_speed[s, k] = (length(idx_t) > 1 ? std(mean_speed_n[s, k, idx_t]) : 0)
	end
	K = idx_n = idx_t = mean_speed_n = mean_flux_n = 0
end

function Measure_s!(Lanes::Int64, S::Int64, N::Int64, A::Int64, mean_speed, mean_flux,
									D_mean_speed, D_mean_flux, speeds, fluxes; sense::Int64 = 1)

	mean_speed_n = zeros(S, Lanes, A)
	mean_flux_n = zeros(S, Lanes, A)

	# ensamble averaging
	for t = 1:A, k = 1:Lanes, s = 1:S
		idx_n = find(view(fluxes, s, k, t, :))
		if sense == 1
			mean_speed_n[s, k, t] = mean(speeds[s, k, t, idx_n])
			mean_flux_n[s, k, t] = mean(fluxes[s, k, t, idx_n])
		else
			mean_speed_n[s, k, t] = mean(speeds[s, k+1, t, idx_n])
			mean_flux_n[s, k, t] = mean(fluxes[s, k+1, t, idx_n])
		end
	end

	# temporal averaging
	for k = 1:Lanes, s = 1:S
		idx_t = find(view(mean_flux_n, s, k, :))
    idx_s = find(view(mean_speed_n, s, k, :))

		mean_flux[s, k] = mean(mean_flux_n[s, k, idx_t])
		mean_speed[s, k] = 1./mean(1./mean_speed_n[s, k, idx_s])

		D_mean_flux[s, k] = std(mean_flux_n[s, k, idx_t])
		D_mean_speed[s, k] = 1./std(1./mean_speed_n[s, k, idx_s])
	end
	idx_n = idx_t = mean_speed_n = mean_flux_n = 0
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using UnitCommitment

@testset "convert" begin
    @testset "EGRET solution" begin
        solution = UnitCommitment.read_egret_solution("fixtures/egret_output.json.gz")
        for attr in ["Is on", "Production (MW)", "Production cost (\$)"]
            @test attr in keys(solution)
            @test "115_STEAM_1" in keys(solution[attr])
            @test length(solution[attr]["115_STEAM_1"]) == 48
        end
        @test solution["Production cost (\$)"]["315_CT_6"][15:20] == [0., 0., 884.44, 1470.71, 1470.71, 884.44]
        @test solution["Startup cost (\$)"]["315_CT_6"][15:20] == [0., 0., 5665.23, 0., 0., 0.]
        @test length(keys(solution["Is on"])) == 154
    end
end


text:
struct GTO
    ζ::Float64
    L::Vec3{Int}
    norm::Float64
end

function GTO(
    ζ::Float64,
    L::AbstractVector{Int},
)
    norm = normalizaton(ζ, Vec3(L))
    return GTO(ζ, L, norm)
end

struct CGTO{N}
    d::SVector{N, Float64}
    gtos::SVector{N, GTO}
    norm::Float64
end

function CGTO(
    d::AbstractVector{Float64},
    ζ::AbstractVector{Float64},
    L::AbstractVector{Int},
)
    N = length(d)
    gtos = Vector{GTO}(undef, N)
    for i = 1:N
        gtos[i] = GTO(ζ[i], L)
    end
    norm = normalizaton(d, gtos)
    return CGTO{N}(d, gtos, norm)
end

Base.iterate(cgto::CGTO, state=1) = state > length(cgto) ? nothing : (cgto[state], state+1)
Base.eltype(::CGTO) = GTO
Base.length(cgto::CGTO) = length(cgto.d)
Base.size(cgto::CGTO) = length(cgto.d)

Base.getindex(cgto::CGTO, i::Int) = cgto.gtos[i]
Base.firstindex(cgto::CGTO) = 1
Base.lastindex(cgto::CGTO) = length(cgto)

struct Intor <: AbstractIntor
    nuclei::Nuclei
    cgtos::Vector{CGTO}
    idx::Vector{Int}
end

function Intor(nuclei::Nuclei, basis::Basis)
    cgtos = CGTO[]
    idx = Int[]
    for i = 1:length(nuclei)
        number = nuclei.numbers[i]
        shells_idx = findfirst(==(number), basis.numbers)
        shells = basis.shells[shells_idx]
        for (j, shell) in enumerate(shells)
            exponents = shell.exponents
            angular_momentum = shell.angular_momentum
            for coefficients in shell.coefficients
                for L in ANGULAR_MOMENTUM[angular_momentum + 1]
                    push!(cgtos, CGTO(coefficients, exponents, collect(L)))
                    push!(idx, i)
                end
            end
        end
    end
    return Intor(nuclei, cgtos, idx)
end

Intor(mole::Mole) = Intor(mole.nuclei, mole.basis)

function (intor::Intor)(intor_name::String)
    func = getfield(@__MODULE__, Symbol(intor_name))
    if startswith(intor_name, "int1e")
        return getints2c(func, intor.nuclei, intor.cgtos, intor.idx)
    elseif startswith(intor_name, "int2e")
        return getints4c(func, intor.nuclei, intor.cgtos, intor.idx)
    end
end

function (intor::Intor)(intor_name::String, args...)
    func = getfield(@__MODULE__, Symbol(intor_name))
    return getints2c(func, intor.nuclei, intor.cgtos, intor.idx, args...)
end

function normalizaton(
    a::Float64,
    L::Int,
    )
    return sqrt(sqrt(2*a/π) * (8*a)^L * factorial(L) / factorial(2*L))
end

function normalizaton(
    a::Float64,
    L::Vec3{Int},
    )
    return normalizaton(a, L[1]) * normalizaton(a, L[2]) * normalizaton(a, L[3])
end

normalizaton(gto::GTO) = normalizaton(gto.ζ, gto.L)

function normalizaton(
    d::AbstractVector{Float64},
    gtos::AbstractVector{GTO},
    )
    N = length(d)
    ret = 0.0
    for i = 1:N, j = 1:N
        ret += d[i] * d[j] * int1e_ovlp(gtos[i], gtos[j], Vec3(zeros(3)))
    end
    return 1 / sqrt(ret)
end

function int1e_ovlp(
    A::GTO,
    B::GTO,
    RAB::Vec3{Float64},
    )
    return A.norm * B.norm * int1e_ovlp(A.ζ, A.L, B.ζ, B.L, RAB)
end

int1e_ovlp(A::CGTO, RA::Vec3{Float64}, B::CGTO, RB::Vec3{Float64}) = contract(int1e_ovlp, A, RA, B, RB)

function int1e_kin(
    A::GTO,
    B::GTO,
    RAB::Vec3{Float64},
    )
    return A.norm * B.norm * int1e_kin(A.ζ, A.L, B.ζ, B.L, RAB)
end

int1e_kin(A::CGTO, RA::Vec3{Float64}, B::CGTO, RB::Vec3{Float64}) = contract(int1e_kin, A, RA, B, RB)

function int1e_rinv(
    A::GTO,
    RA::Vec3{Float64},
    B::GTO,
    RB::Vec3{Float64},
    RC::Vec3{Float64},
    )
    RAB = RA - RB
    RP = gaussian_product_center(A.ζ, RA, B.ζ, RB)
    RPC = RP - RC
    return A.norm * B.norm * int1e_rinv(A.ζ, A.L, B.ζ, B.L, RAB, RPC)
end

int1e_rinv(A::CGTO, RA::Vec3{Float64}, B::CGTO, RB::Vec3{Float64}, RC::Vec3{Float64}) = contract(int1e_rinv, A, RA, B, RB, RC)

int1e_nuc(A::CGTO, RA::Vec3{Float64}, B::CGTO, RB::Vec3{Float64}, Q::Float64, RC::Vec3{Float64}) = -Q * contract(int1e_rinv, A, RA, B, RB, RC)

function int2e(
    A::GTO,
    RA::Vec3{Float64},
    B::GTO,
    RB::Vec3{Float64},
    C::GTO,
    RC::Vec3{Float64},
    D::GTO,
    RD::Vec3{Float64},
    )
    RAB = RA - RB
    RCD = RC - RD
    RP = gaussian_product_center(A.ζ, RA, B.ζ, RB)
    RQ = gaussian_product_center(C.ζ, RC, D.ζ, RD)
    RPQ = RP - RQ
    return A.norm * B.norm * C.norm * D.norm * int2e(A.ζ, A.L, B.ζ, B.L, C.ζ, C.L, D.ζ, D.L, RAB, RCD, RPQ)
end

function int2e(
    A::CGTO,
    RA::Vec3{Float64},
    B::CGTO,
    RB::Vec3{Float64},
    C::CGTO,
    RC::Vec3{Float64},
    D::CGTO,
    RD::Vec3{Float64},
    )
    contract(int2e, A, RA, B, RB, C, RC, D, RD)
end

function contract(f::Function, A::CGTO, RA::Vec3{Float64}, B::CGTO, RB::Vec3{Float64})
    RAB = RA - RB
    ret = 0.0
    for i = 1:length(A), j = 1:length(B)
        ret += A.d[i] * B.d[j] * f(A[i], B[j], RAB)
    end
    return A.norm * B.norm * ret
end

function contract(f::Function, A::CGTO, RA::Vec3{Float64}, B::CGTO, RB::Vec3{Float64}, RC::Vec3{Float64})
    ret = 0.0
    for i = 1:length(A), j = 1:length(B)
        ret += A.d[i] * B.d[j] * f(A[i], RA, B[j], RB, RC)
    end
    return A.norm * B.norm * ret
end

function contract(f::Function, A::CGTO, RA::Vec3{Float64}, B::CGTO, RB::Vec3{Float64}, C::CGTO, RC::Vec3{Float64}, D::CGTO, RD::Vec3{Float64})
    ret = 0.0
    for i = 1:length(A), j = 1:length(B), k = 1:length(C), l = 1:length(D)
        ret += A.d[i] * B.d[j] * C.d[k] * D.d[l] *
            f(A[i], RA, B[j], RB, C[k], RC, D[l], RD)
    end
    return A.norm * B.norm * C.norm * D.norm * ret
end

function getints2c(
    func::Function,
    nuclei::Nuclei,
    cgtos::Vector{CGTO},
    idx::Vector{Int},
    )
    n = length(cgtos)
    M = zeros(n, n)
    for i=1:n, j=i:n
        A = cgtos[i]
        RA = nuclei.positions[idx[i]]
        B = cgtos[j]
        RB = nuclei.positions[idx[j]]
        M[i, j] = func(A, RA, B, RB)
    end
    symmetrize!(M)
    return M
end

function getints2c(
    func::typeof(int1e_nuc),
    nuclei::Nuclei,
    cgtos::Vector{CGTO},
    idx::Vector{Int},
    )
    n = length(cgtos)
    M = zeros(n, n)
    for i=1:n, j=i:n
        A = cgtos[i]
        RA = nuclei.positions[idx[i]]
        B = cgtos[j]
        RB = nuclei.positions[idx[j]]
        for (Q, RC) in zip(nuclei.charges, nuclei.positions)
            M[i,j] += func(A, RA, B, RB, Q, RC)
        end
    end
    symmetrize!(M)
    return M
end

function getints2c(
    func::typeof(int1e_nuc),
    nuclei::Nuclei,
    cgtos::Vector{CGTO},
    idx::Vector{Int},
    pg::AbstractParticleGroup,
    )
    n = length(cgtos)
    M = zeros(n, n)
    for i=1:n, j=i:n
        A = cgtos[i]
        RA = nuclei.positions[idx[i]]
        B = cgtos[j]
        RB = nuclei.positions[idx[j]]
        for (Q, RC) in zip(pg.charges, pg.positions)
            M[i,j] += func(A, RA, B, RB, Q, RC)
        end
    end
    symmetrize!(M)
    return M
end

function getints2c(
    func::typeof(int1e_rinv),
    nuclei::Nuclei,
    cgtos::Vector{CGTO},
    idx::Vector{Int},
    point::Vec3{Float64},
    )
    n = length(cgtos)
    M = zeros(n, n)
    for i=1:n, j=i:n
        A = cgtos[i]
        RA = nuclei.positions[idx[i]]
        B = cgtos[j]
        RB = nuclei.positions[idx[j]]
        M[i,j] += func(A, RA, B, RB, point)
    end
    symmetrize!(M)
    return M
end

function getints4c(
    func::Function,
    nuclei::Nuclei,
    cgtos::Vector{CGTO},
    idx::Vector{Int},
    )
    n = length(cgtos)
    T = zeros(n, n, n, n)
    for i = 1:n, j = i:n
        for k = i:n, l = k:n
            A = cgtos[i]
            RA = nuclei.positions[idx[i]]
            B = cgtos[j]
            RB = nuclei.positions[idx[j]]
            C = cgtos[k]
            RC = nuclei.positions[idx[k]]
            D = cgtos[l]
            RD = nuclei.positions[idx[l]]
            T[i,j,k,l] = func(A, RA, B, RB, C, RC, D, RD)
        end
    end
    symmetrize!(T)
    return T
end


text:
to_length(n::Integer) = n
to_length(t::Tuple{Vararg{Integer}}) = prod(t)
to_length(t::Tuple{Vararg{AbstractUnitRange}}) = prod(map(length, t))

default_initializer(::Type{T}, sz) where {T} = similar(arrayof(T), to_length(sz))

"""
    collect_columns(itr)

Collect an iterable as a `Columns` object if it iterates `Tuples` or `NamedTuples`, as a normal
`Array` otherwise.

# Examples

    s = [(1,2), (3,4)]
    collect_columns(s)

    s2 = Iterators.filter(isodd, 1:8)
    collect_columns(s2)
"""
collect_columns(v) = vec(collect_structarray(v, initializer = default_initializer))

_append!!(v, itr) = append!!(v, itr)
_append!!(v::StructArray{NamedTuple{(),Tuple{}}}, itr) = collect_columns(itr)

collect_columns_flattened(itr) = foldl(_append!!, itr, init = Columns())


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libX11"], :libX11),
    LibraryProduct(prefix, ["libX11-xcb"], :libX11_xcb),
    LibraryProduct(prefix, ["libXau"], :libXau),
    LibraryProduct(prefix, ["libxcb-composite"], :libxcb_composite),
    LibraryProduct(prefix, ["libxcb-damage"], :libxcb_damage),
    LibraryProduct(prefix, ["libxcb-dpms"], :libxcb_dpms),
    LibraryProduct(prefix, ["libxcb-dri2"], :libxcb_dri2),
    LibraryProduct(prefix, ["libxcb-dri3"], :libxcb_dri3),
    LibraryProduct(prefix, ["libxcb-glx"], :libxcb_glx),
    LibraryProduct(prefix, ["libxcb-present"], :libxcb_present),
    LibraryProduct(prefix, ["libxcb-randr"], :libxcb_randr),
    LibraryProduct(prefix, ["libxcb-record"], :libxcb_record),
    LibraryProduct(prefix, ["libxcb-render"], :libxcb_render),
    LibraryProduct(prefix, ["libxcb-res"], :libxcb_res),
    LibraryProduct(prefix, ["libxcb-screensaver"], :libxcb_screensaver),
    LibraryProduct(prefix, ["libxcb-shape"], :libxcb_shape),
    LibraryProduct(prefix, ["libxcb-shm"], :libxcb_shm),
    LibraryProduct(prefix, ["libxcb"], :libxcb),
    LibraryProduct(prefix, ["libxcb-sync"], :libxcb_sync),
    LibraryProduct(prefix, ["libxcb-xf86dri"], :libxcb_xf86dri),
    LibraryProduct(prefix, ["libxcb-xfixes"], :libxcb_xfixes),
    LibraryProduct(prefix, ["libxcb-xinerama"], :libxcb_xinerama),
    LibraryProduct(prefix, ["libxcb-xinput"], :libxcb_xinput),
    LibraryProduct(prefix, ["libxcb-xkb"], :libxcb_xkb),
    LibraryProduct(prefix, ["libxcb-xtest"], :libxcb_xtest),
    LibraryProduct(prefix, ["libxcb-xvmc"], :libxcb_xvmc),
    LibraryProduct(prefix, ["libxcb-xv"], :libxcb_xv),
    LibraryProduct(prefix, ["libXcomposite"], :libXcomposite),
    LibraryProduct(prefix, ["libXcursor"], :libXcursor),
    LibraryProduct(prefix, ["libXdamage"], :libXdamage),
    LibraryProduct(prefix, ["libXext"], :libXext),
    LibraryProduct(prefix, ["libXfixes"], :libXfixes),
    LibraryProduct(prefix, ["libXinerama"], :libXinerama),
    LibraryProduct(prefix, ["libXi"], :libXi),
    LibraryProduct(prefix, ["libXrandr"], :libXrandr),
    LibraryProduct(prefix, ["libXrender"], :libXrender),
    LibraryProduct(prefix, ["libxshmfence"], :libxshmfence),
    LibraryProduct(prefix, ["libXtst"], :libXtst),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaBinaryWrappers/X11_jll.jl/releases/download/X11-v1.6.8+4"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc) => ("$bin_prefix/X11.v1.6.8.aarch64-linux-gnu.tar.gz", "3e0d549012b0f547e2be29434a38a9fc539bf3907563c35bb1b5b671f180eb2f"),
    Linux(:aarch64, libc=:musl) => ("$bin_prefix/X11.v1.6.8.aarch64-linux-musl.tar.gz", "5e23e284f230561e9511b39a787389f1572481adefeeff67aaf8b0324b20473a"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf) => ("$bin_prefix/X11.v1.6.8.arm-linux-gnueabihf.tar.gz", "005f53458fedc6a62539fc3e8251be1db6ea9959ff664dde67a6e65894b78ef5"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf) => ("$bin_prefix/X11.v1.6.8.arm-linux-musleabihf.tar.gz", "8b23a814b230e2daf224fe50e9a04b36d0cdc330a0b46d93f6aae35d5f3bc319"),
    Linux(:i686, libc=:glibc) => ("$bin_prefix/X11.v1.6.8.i686-linux-gnu.tar.gz", "9e453780811b27923660c43aa96c7f46a7c946ac34b626fca4f0de920857b7ec"),
    Linux(:i686, libc=:musl) => ("$bin_prefix/X11.v1.6.8.i686-linux-musl.tar.gz", "70ac7542eed7d4bad979eaa84ecdb50ed735a435b3cf024b677255efc421f51e"),
    Linux(:powerpc64le, libc=:glibc) => ("$bin_prefix/X11.v1.6.8.powerpc64le-linux-gnu.tar.gz", "b3974379065f2dd3c9ef16ffc88ae57c7b6a2f9367ec2d055f527d92ceaac5ec"),
    Linux(:x86_64, libc=:glibc) => ("$bin_prefix/X11.v1.6.8.x86_64-linux-gnu.tar.gz", "31ccdbc12db080ea40d6c04b4197a26135dba7af63181b36ed1590e21f8b3ce1"),
    Linux(:x86_64, libc=:musl) => ("$bin_prefix/X11.v1.6.8.x86_64-linux-musl.tar.gz", "f8f576c34046f84b10df0c4116ba602ac43c55409386d6c34a17c97923d4a788"),
    FreeBSD(:x86_64) => ("$bin_prefix/X11.v1.6.8.x86_64-unknown-freebsd11.1.tar.gz", "164ad28c9573225f59dd987eaf4580c8b2855a560cbe1ba443907ed08fa2f63c"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
using Test
using Expronicon

@testset "one line expression" begin
    @test sprint_expr(:(name::type)) == "name::type"
    @test sprint_expr(:("abc $name")) == "\"abc \$name\""
    @test sprint_expr(:(compare_expr(lhs, rhs::Int))) == "compare_expr(lhs, rhs::Int)"
    @test sprint_expr(:(compare_expr(lhs, rhs::Int; a = "c"))) == "compare_expr(lhs, rhs::Int; a = \"c\")"
    @test sprint_expr(:(return a, b, c)) == "return a, b, c"
    @test sprint_expr(:(a = b)) == "a = b"
    @test sprint_expr(Expr(:kw, :a, :b)) == "a = b"
    @test sprint_expr(:(!x)) == "!x"
    @test sprint_expr(:(x + 1)) == "x + 1"
    @test sprint_expr(:(x * 1)) == "x * 1"
    str = sprint_expr(:(f(x) = x))
    @test occursin("f(x) = x", str)
    str = sprint_expr(:(x->2x))
    @test occursin("x -> 2 * x", str)
    @test sprint_expr(:(Type{T <: Real})) == "Type{T <: Real}"
end

print_expr(:(let x, y
    x + 1
    y + 1
end))

print_expr(:(if x > 0
    x + 1
end))

print_expr(quote
x < 0
end)

print_expr(:(if x > 0
    x + 1
elseif x > 1
    x + 2
elseif x < 0
    x + 3
else
    x + 4
end))

print_expr(:(function foo(x, y::T; z::Int = 1) where {N, T <: Real}
    x + 1
end))

ex = @expr struct Foo <: Super
    x::Int
    Foo(x::Int) = new(x)
    Foo(x::Int) = new(x)
end

print_expr(ex)

ex = @expr mutable struct Goo <: Super
    x::Int
    Foo(x::Int) = new(x)
    Foo(x::Int) = new(x)
end

print_expr(ex)

ex = :(function compare_expr(lhs, rhs)
@switch (lhs, rhs) begin
    @case (::Symbol, ::Symbol)
        lhs === rhs
    @case (Expr(:curly, name, lhs_vars...), Expr(:curly, &name, rhs_vars...))
        all(map(compare_vars, lhs_vars, rhs_vars))
    @case (Expr(:where, lbody, lparams...), Expr(:where, rbody, rparams...))
        compare_expr(lbody, rbody) &&
            all(map(compare_vars, lparams, rparams))
    @case (Expr(head, largs...), Expr(&head, rargs...))
            isempty(largs) && isempty(rargs) ||
        (length(largs) == length(rargs) && all(map(compare_expr, largs, rargs)))
    # ignore LineNumberNode
    @case (::LineNumberNode, ::LineNumberNode)
        true
    @case _
        lhs == rhs
end
end)

print_expr(ex)

ex = :(try
    1 + 1
catch e
    rethrow(ex)
end)

print_expr(ex)

ex = :(try
    1 + 1
finally
    rethrow(ex)
end)

print_expr(ex)

ex = :(try
    1 + 1
catch e
    rethrow(ex)
finally
    1 + 2
end)

print_expr(ex)

def = @expr JLFunction function foo(x, y)
    1 + 1
end

print_expr(def)

def = @expr JLKwStruct struct Moo
    x::Int = 1
end

print_expr(def)

ex = @expr for i in 1:10, j in 1:10
    M[i, j] += 1
end

print_expr(ex)

ex = @expr function foo(i, j)
    for i in 1:10
        M[i, j] += 1
    end
end
print_expr(ex)

ex = quote
    """
        foo(i, j)

    test function. test function.
    test function. test function.
    """
    function foo(i, j)
        for i in 1:10
            M[i, j] += 1
        end
    end
end

print_expr(ex)

ex = quote
    """
        foo(i, j)

    test function. test function.
    test function. test function.
    """
    function foo(i, j)
        function goo(x)
            return x
        end
    end
end

print_expr(ex)

function broutine2x2_m_kernel_expr(idx_1, idx_2)
    return quote
        ST1 = U11 * st[b, $idx_1] + U12 * st[b, $idx_2]
        ST2 = U21 * st[b, $idx_1] + U22 * st[b, $idx_2]

        st[b, $idx_1] = ST1
        st[b, $idx_2] = ST2
    end
end

ex = quote
    Base.Cartesian.@nexprs 16 k->begin
    $(broutine2x2_m_kernel_expr(:(i+k), :(i+step_1+k)))
    end
end

print_expr(ex)

ex = quote
    try
        ex = include_string(mod, "quote $code end", path)
        mod.eval(mod.eval(ex))
        return nothing
    finally
        if prev === nothing
            delete!(tls, :SOURCE_PATH)
        else
            tls[:SOURCE_PATH] = prev
        end
    end
end

print_expr(ex)