text:
module ReindexTests

using Test
using FillArrays
using Gridap.Arrays
using Gridap.TensorValues

a = [1,2,3]
b = lazy_map(Reindex(a),[3,2,1])
for i=1:3
    b[i]=i
end
@test b == [1,2,3]
@test a == [3,2,1]

a = [[1,2,4,5],[2,4,6,7],[4,3,5,1],[2,3]]
b = [3,1,2]

c = lazy_map(Reindex(a),b)

d = a[b]
test_array(c,d)

a = Fill(30.0,10)
b = [3,1,2]

c = lazy_map(Reindex(a),b)
@test isa(c,Fill)

d = a[b]
test_array(c,d)

a = Fill(1.0,5,5)
b = [13 23; 15 25]

c = lazy_map(Reindex(a),b)
@test isa(c,Fill)

d = a[b]
test_array(c,d)

a = CompressedArray([30,40,10,20,30],[1,2,3,5,3,1,4,2])
b = [3,1,2]

c = lazy_map(Reindex(a),b)
@test isa(c,CompressedArray)

d = a[b]
test_array(c,d)

x = [1,2,3,5,3,1,4,2]
a = lazy_map(-,x)
b = [3,1,2]

c = lazy_map(Reindex(a),b)

@test isa(c,LazyArray)

d = [a[bi] for bi in b]
test_array(c,d)

d = a[b]
test_array(c,d)

# x = [[1,2,4,5],[2,4,6,7],[4,3,5,1],[2,3]]
# a = Table(x)
# b = [3,1,2]
# c = lazy_map(Reindex(a),b)
# d = a[b]
# test_array(c,d)

p = VectorValue(2.0,1.0)
data = [2,3,1,3,4,4,3,2,5,4,3,4]
ptrs = [1,4,4,7,13]
gid_to_val = [p, 2*p, 3*p, -p, p]
lid_to_gid = Table(data,ptrs)
ca = lazy_map(Broadcasting(Reindex(gid_to_val)),lid_to_gid)
a = [ gid_to_val[gids] for gids in lid_to_gid ]
test_array( ca, a )

ids = [3,1,2]
r = lazy_map(Reindex(ca),ids)
@test isa(r,LazyArray)
test_array(r,ca[ids])

i_to_v=[[1, 2], [5, 6], [1, 5], [2, 6], [2, 3], [6, 7], [3, 7], [3, 4], [7, 8], [4, 8], [9, 10], [5, 9], [6, 10], [10, 11], [7, 11], [11, 12], [8, 12], [13, 14], [9, 13], [10, 14], [14, 15], [11, 15], [15, 16], [12, 16]]
j_to_i=Int[]
lid_to_gid=lazy_map(Reindex(i_to_v),j_to_i)
gid_to_val=VectorValue{2,Float64}[(0.0, 0.25), (0.25, 0.25), (0.5, 0.25), (0.75, 0.25), (0.0, 0.5), (0.25, 0.5), (0.5, 0.5), (0.75, 0.5), (0.0, 0.75), (0.25, 0.75), (0.5, 0.75), (0.75, 0.75), (0.0, 1.0), (0.25, 1.0), (0.5, 1.0), (0.75, 1.0)]
f2 = lazy_map(Broadcasting(Reindex(gid_to_val)),lid_to_gid)

test_array(f2,Vector{VectorValue{2,Float64}}[])

values = Float64[]
indices = Int[]
c = lazy_map(Reindex(values),indices)
test_array(c,Float64[])

end # module


text:
using jInv.Mesh
using jInv.LinearSolvers
using Test
using KrylovMethods
using Multigrid
using LinearAlgebra
using SparseArrays

println("===  Example 2D DivSigGrad ====");

domain = [0.0, 1.0, 0.0, 1.0];
n      = [50,50];
Mr     = getRegularMesh(domain,n)
G      = getNodalGradientMatrix(Mr);
m      = sparse(Diagonal(exp.(randn(size(G,1)))));
Ar     = G'*m*G;
Ar     = Ar + 1e-1*norm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));
N      = size(Ar,2);
B      = Ar*rand(N,4);

levels      = 5;
numCores 	= 8;
maxIter     = 15;
relativeTol = 1e-2;
relaxType   = "SPAI";
relaxParam  = 1.0;
relaxPre 	= 2;
relaxPost   = 2;
cycleType   ='V';
coarseSolveType = "Julia";

MG = getMGparam(Float64,Int64,levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType);

sSAPCG   = getSA_AMGsolver(MG, "PCG",sym=1,out=1);
X,  = solveLinearSystem(Ar,B,sSAPCG);
@test norm(Ar*X-B)/norm(B) < sSAPCG.tol

clear!(sSAPCG);
sSABiCG   = copySolver(sSAPCG);
sSABiCG.Krylov = "BiCGSTAB"
X,  = solveLinearSystem(Ar,B,sSABiCG);
@test norm(Ar*X-B)/norm(B) < sSABiCG.tol

ext: jl


text:
using Distributed, StatsBase, Distributions, PyPlot
using TransitionalMCMC

lb  = -15
ub  = 15

fT(x) = logpdf(Uniform(lb,ub), x)
sample_fT(Nsamples) = rand(Uniform(lb,ub), Nsamples, 1)

function log_fD_T(x)
    return log(pdf(Normal(0,1), x[1]) + pdf(Normal(5,0.2), x[1]))
end

Nsamples = 2000
samps, acc = tmcmc(log_fD_T, fT, sample_fT, Nsamples)

plt.hist(samps,50)


text:
# ---
# title: Crop
# cover: crop.gif
# ---

# Subset image using `Crop` and `CropNative`

using Augmentor
using ImageShow, ImageCore
using OffsetArrays

img_in = testpattern(RGB, ratio=0.5)
img_out = augment(img_in, Crop(20:75,25:120))

mosaicview(img_in, img_out; fillvalue=colorant"white", nrow=1)

# If the input image is plain arrays without offset indices, then `Crop` and `CropNative` is equivalent.

augment(img_in, Crop(20:75,25:120)) == augment(img_in, CropNative(20:75,25:120))

# Whether you should use `Crop` or `CropNative` depends on if you want to take the index offset
# of the input image into consideration.

imgo_in = OffsetArray(img_in, -50, -50)
imgo_out = augment(imgo_in, Crop(20:75,25:120))
imgo_out_native = augment(imgo_in, CropNative(20:75,25:120))

(
    imgo_in[(first.(axes(imgo_in)) .+ (20, 25))...] == imgo_out[1, 1],
    imgo_in[20, 25] == imgo_out_native[1, 1]
)

# A typical scenario that you may want to use `CropNative` is when you have affine operations, e.g.,
# `Rotate` and `ShearX`.

mosaicview(
    augment(img_in, Rotate(30) |> Crop(axes(img_in))),
    augment(img_in, Rotate(30) |> CropNative(axes(img_in))),

    augment(img_in, ShearX(10) |> Crop(axes(img_in))),
    augment(img_in, ShearX(10) |> CropNative(axes(img_in)));

    fillvalue=colorant"white", rowmajor=true, nrow=2, npad=10
)

# ## Reference

#md # ```@docs
#md # Crop
#md # CropNative
#md # ```

## save covers #src
using ImageMagick #src
include(joinpath("..", "assets", "utilities.jl")) #src
cover = make_gif(img_in, Crop(20:75,25:120), 2) #src
ImageMagick.save("crop.gif", cover; fps=1) #src


text:
using SubspacePreservingSparsification
using SubspacePreservingSparsification: near_zero_row_col
using MAT
using DelimitedFiles
using LinearAlgebra

sq = "square10"

function default_ssa_compute_bin_no_null(M)
    M = sparsify(M, 0.4, 2, 50)
    return M
end

function default_ssa_compute_no_bin_no_null(M)
    M = sparsify(M, 0.4, 2, 0)
    return M
end

function my_default_ssa_compute_bin_null(M)
    M = sparsify(M, 0.4, 2, 50, true)
    return M
end

function my_default_compute_sparse_pattern(M)
    _, rnull, lnull = pinv_qr(M)

    # sparsity pattern
    num_near_zero_rows, num_near_zero_cols = near_zero_row_col(M)
    min_per_row = max(0, min(size(rnull)[2] - num_near_zero_cols, size(M)[2]))
    min_per_col = max(0, min(size(lnull)[2] - num_near_zero_rows, size(M)[1]))
    M = sparsity_pattern(M, 0.4, 2, min_per_row, min_per_col)
    return M
end

function my_default_compute_bin_pattern(M)
    _, rnull, lnull = pinv_qr(M)

    # sparsity pattern
    num_near_zero_rows, num_near_zero_cols = near_zero_row_col(M)
    min_per_row = max(0, min(size(rnull)[2] - num_near_zero_cols, size(M)[2]))
    min_per_col = max(0, min(size(lnull)[2] - num_near_zero_rows, size(M)[1]))
    M_pat = sparsity_pattern(M, 0.4, 2, min_per_row, min_per_col)
    M = bin_sparse_matrix!(M, M_pat, 50)
    return M
end

function init_script()
    ROOT_DIR = realpath(dirname(@__FILE__) * "/..")
    cd(ROOT_DIR)
end

init_script()

trans_names = readdlm(joinpath("test", "test_data", sq, "trans_names.txt"))

transformations = Dict(
    "pinv_rrqr" => pinv_qr,
    "default_ssa_compute_bin_no_null" => default_ssa_compute_bin_no_null,
    "default_ssa_compute_no_bin_no_null" => default_ssa_compute_no_bin_no_null,
    "default_ssa_compute_bin_null" => my_default_ssa_compute_bin_null,
    "sparse_pattern" => my_default_compute_sparse_pattern,
    "bin_pattern" => my_default_compute_bin_pattern,
)

# load untransformed matrix
file = matopen(joinpath("test", "test_data", sq, "untransformed_matrices", "with_nullspace.mat"))
A = read(file, "A")
close(file)

for j_trans_name in trans_names
    trans = transformations[j_trans_name]
    trans_A = @inferred trans(A)
    if typeof(trans_A) <: Tuple{Any, Any, Any}
        trans_A = trans_A[1]
    end
    if typeof(trans_A) != Matrix{Float64}
        trans_A = Matrix{Float64}(trans_A)
    end
    local file = matopen(joinpath("test", "test_data", sq,
    "transformed_matrices", "matlab", "with_nullspace", j_trans_name*".mat"))
    trans_A_mat = read(file, "trans_A")  # TODO: no sparse matrix
    close(file)
    @test trans_A ≈ trans_A_mat
end


text:
 open("src/RandomTestMatrix.jl") do file
    s = false

    while !eof(file)
        l = readline(file)
        if length(l) >= 3 && l[1:3] == "\"\"\""
           s = !s
           l = l[3:end]
        end
        if s
            print(readline(file));
        end
    end
end


text:
@testset "Encoding video across all supported colortypes" begin
    for el in [UInt8, RGB{N0f8}]
        @testset "Encoding $el imagestack" begin
            n = 100
            imgstack = map(x -> rand(el, 100, 100), 1:n)
            encoder_options = (color_range = 2, crf = 0, preset = "medium")
            VideoIO.save(tempvidpath, imgstack, framerate = 30, encoder_options = encoder_options)
            @test stat(tempvidpath).size > 100
            @test VideoIO.openvideo(VideoIO.counttotalframes, tempvidpath) == n
        end
    end
end

@testset "Simultaneous encoding and muxing" begin
    n = 100
    encoder_options = (color_range = 2,)
    container_private_options = (movflags = "+write_colr",)
    for el in [Gray{N0f8}, Gray{N6f10}, RGB{N0f8}, RGB{N6f10}]
        codec_name = el <: RGB ? "libx264rgb" : "libx264" # the former is necessary for lossless RGB
        for scanline_arg in [true, false]
            for sz in [100, 128] # 100 tests where julia<>ffmpeg imgbuf size doesn't match, 128 when it does
                @testset "Encoding $el imagestack, scanline_major = $scanline_arg, size = $sz" begin
                    img_stack = map(x -> rand(el, sz, sz), 1:n)
                    encoder_private_options = (crf = 0, preset = "medium")
                    VideoIO.save(
                        tempvidpath,
                        img_stack;
                        codec_name = codec_name,
                        encoder_private_options = encoder_private_options,
                        encoder_options = encoder_options,
                        container_private_options = container_private_options,
                        scanline_major = scanline_arg,
                    )
                    @test stat(tempvidpath).size > 100
                    f = VideoIO.openvideo(tempvidpath, target_format = VideoIO.get_transfer_pix_fmt(el))
                    try
                        notempty = !eof(f)
                        @test notempty
                        if notempty
                            img = read(f)
                            test_img = scanline_arg ? parent(img) : img
                            i = 1
                            if el in [Gray{N0f8}, RGB{N0f8}]
                                @test test_img == img_stack[i]
                            else
                                @test_broken test_img == img_stack[i]
                            end
                            while !eof(f) && i < n
                                read!(f, img)
                                i += 1
                                if el in [Gray{N0f8}, RGB{N0f8}]
                                    @test test_img == img_stack[i]
                                else
                                    @test_broken test_img == img_stack[i]
                                end
                            end
                            @test i == n
                        end
                    finally
                        close(f)
                    end
                end
            end
        end
    end
end

@testset "Monochrome rescaling" begin
    nw = nh = 100
    s = VideoIO.GrayTransform()
    s.srcframe.color_range = VideoIO.AVCOL_RANGE_JPEG
    s.dstframe.color_range = VideoIO.AVCOL_RANGE_MPEG
    s.srcframe.format = VideoIO.AV_PIX_FMT_GRAY8
    s.dstframe.format = VideoIO.AV_PIX_FMT_GRAY8
    s.src_depth = s.dst_depth = 8
    f, src_t, dst_t = VideoIO.make_scale_function(s)
    @test f(0x00) == 16
    @test f(0xff) == 235
    s.srcframe.format = s.dstframe.format = VideoIO.AV_PIX_FMT_GRAY10LE
    s.src_depth = s.dst_depth = 10
    f, src_t, dst_t = VideoIO.make_scale_function(s)
    @test f(0x0000) == 64
    @test f(UInt16(1023)) == 940

    # Test that range conversion is working properly
    img_full_range = reinterpret(UInt16, test_tone(N6f10, nw, nh))
    writer = VideoIO.open_video_out(
        tempvidpath,
        img_full_range;
        target_pix_fmt = VideoIO.AV_PIX_FMT_GRAY10LE,
        scanline_major = true,
    )
    @test VideoIO.get_codec_name(writer) != "None"
    try
        VideoIO.write(writer, img_full_range)
        bwidth = nw * 2
        buff = Vector{UInt8}(undef, bwidth * nh)
        # Input frame should be full range
        copy_imgbuf_to_buf!(
            buff,
            bwidth,
            nh,
            writer.frame_graph.srcframe.data[1],
            writer.frame_graph.srcframe.linesize[1],
        )
        raw_vals = reinterpret(UInt16, buff)
        @test extrema(raw_vals) == (0x0000, 0x03ff)
        # Output frame should be limited range
        copy_imgbuf_to_buf!(
            buff,
            bwidth,
            nh,
            writer.frame_graph.dstframe.data[1],
            writer.frame_graph.dstframe.linesize[1],
        )

        @test extrema(raw_vals) == (0x0040, 0x03ac)
    finally
        VideoIO.close_video_out!(writer)
    end
    @test_throws ErrorException VideoIO.write(writer, img_full_range)
end

@testset "Encoding monochrome videos" begin
    encoder_private_options = (crf = 0, preset = "fast")
    nw = nh = 100
    nf = 5
    for elt in (N0f8, N6f10)
        if elt == N0f8
            limited_min = 16
            limited_max = 235
            full_min = 0
            full_max = 255
            target_fmt = VideoIO.AV_PIX_FMT_GRAY8
        else
            limited_min = 64
            limited_max = 940
            full_min = 0
            full_max = 1023
            target_fmt = VideoIO.AV_PIX_FMT_GRAY10LE
        end
        img_stack_full_range = make_test_tones(elt, nw, nh, nf)
        # Test that full-range input is automatically converted to limited range
        VideoIO.save(
            tempvidpath,
            img_stack_full_range,
            target_pix_fmt = target_fmt,
            encoder_private_options = encoder_private_options,
        )

        minp, maxp = get_raw_luma_extrema(elt, tempvidpath, nw, nh)
        @test minp > full_min
        @test maxp < full_max

        # Test that this conversion is NOT done if output video is full range
        VideoIO.save(
            tempvidpath,
            img_stack_full_range,
            target_pix_fmt = target_fmt,
            encoder_private_options = encoder_private_options,
            encoder_options = (color_range = 2,),
        )
        minp, maxp = get_raw_luma_extrema(elt, tempvidpath, nw, nh)
        @test minp == full_min
        @test maxp == full_max

        # Test that you can override this automatic conversion when writing videos
        img_stack_limited_range = make_test_tones(elt, nw, nh, nf, limited_min, limited_max)
        VideoIO.save(
            tempvidpath,
            img_stack_limited_range,
            target_pix_fmt = target_fmt,
            encoder_private_options = encoder_private_options,
            input_colorspace_details = VideoIO.VioColorspaceDetails(),
        )
        minp, maxp = get_raw_luma_extrema(elt, tempvidpath, nw, nh)
        @test minp > full_min # Actual N6f10 values are messed up during encoding
        @test maxp < full_max # Actual N6f10 values are messed up during encoding
    end
end

@testset "Encoding video with rational frame rates" begin
    n = 100
    fr = 59 // 2 # 29.5
    target_dur = 3.39
    @testset "Encoding with frame rate $(float(fr))" begin
        imgstack = map(x -> rand(UInt8, 100, 100), 1:n)
        encoder_options = (color_range = 2, crf = 0, preset = "medium")
        VideoIO.save(tempvidpath, imgstack, framerate = fr, encoder_options = encoder_options)
        @test stat(tempvidpath).size > 100
        measured_dur_str = VideoIO.FFMPEG.exe(
            `-v error -show_entries format=duration -of default=noprint_wrappers=1:nokey=1 $(tempvidpath)`,
            command = VideoIO.FFMPEG.ffprobe,
            collect = true,
        )
        @test parse(Float64, measured_dur_str[1]) == target_dur
    end
end

@testset "Encoding video with float frame rates" begin
    n = 100
    fr = 29.5 # 59 // 2
    target_dur = 3.39
    @testset "Encoding with frame rate $(float(fr))" begin
        imgstack = map(x -> rand(UInt8, 100, 100), 1:n)
        encoder_options = (color_range = 2, crf = 0, preset = "medium")
        VideoIO.save(tempvidpath, imgstack, framerate = fr, encoder_options = encoder_options)
        @test stat(tempvidpath).size > 100
        measured_dur_str = VideoIO.FFMPEG.exe(
            `-v error -show_entries format=duration -of default=noprint_wrappers=1:nokey=1 $(tempvidpath)`,
            command = VideoIO.FFMPEG.ffprobe,
            collect = true,
        )
        @test parse(Float64, measured_dur_str[1]) == target_dur
    end
end


text:
# The rules below are inspired by those found here:
# https://github.com/ebnn/draftcheck/blob/master/draftcheck/rules.py
rule("centering",r"\\begin{center}", "Use \\centering instead of \\begin{center}")
rule("footnote_spaces",r"\s+\\footnote{", "Do not precede footnotes with spaces.")
rule("citations_before_period",r"\.\\cite{", "Place citations before periods with a non-breaking space."),
rule("url", r"(?<!\\url{)(\bhttps?://)[^\s.]+\.[-A-Za-z0-9+&@#/%?=~_|!:,.;]+", "Wrap URLs with the \\url command.")
rule("left_and_right_quotes",r"\s\'.+?\'[\s\.,]", "Use left and right quotation marks ` and ' rather than '")
rule("citation_nouns",r"\b(:?in|as|on|by)[ ~]\\cite{", "Avoid using citations as nouns.")
rule("non-breaking_space_citation", r"[^~]\\ref{", "Place a single, non-breaking space '~' before references.")
rule("escape_percent",r"\d+%", "Escape percentages with backslash.")
rule("punctuation", r"\s[,;.!?]", "Do not precede punctuation characters with spaces.")
rule("mutiply",r"\d+\s?x\d+", "When multiplying, do not use 'x' and instead use \$\\times\$.")
rule("duplicate",r"\b([a-z]+)\s+\1\b(?![^{]*})", "Remove duplicated word.")
rule("ellipses",r"\.\.\.", "Typeset ellipses by \\dots, not ...")
rule("non-breaking_space_periods",r"\.\\cite{", "Place citations before periods with a non-breaking space.")
rule("double_quotes",r"\"", "Use left and right quotation marks `` and '' rather than \".")
rule("begin_equation",r"^\$\$", "Use \\[ or \\begin{equation} instead of \$\$")
rule("endash", r"\d\s?-\s?\d", "Use endash '--' for numeric ranges instead of hyphens.")
rule("footnotes_after_punc", r"\\footnote{.+?}[,;.?]", "Place footnotes after punctuation marks.")
rule("angle_brac",r"<[^\s](.+?)[^\s]>", "Use \\langle and \\rangle instead of '<' and '>' for angle brackets.")
rule("use_cite", r"\\cite{.+?}\s?\\cite{", "Use \\cite{...,...} for multiple citations")
rule("space_num_unit",r"\d(m|A|kg|s|K|mol|cd)\b", "Place a non-breaking space between a number and its unit.")
rule("backslash_space",r"\b(e\.g\.|i\.e\.)\s+", "Place a '\\ ' (backslash space) after the period of an abbreviation.")
rule("def_command",r"\\def\\[a-z]+{", "Do not use the \\def command. Use \\newcommand instead.")
rule("sloppy",r"\\sloppy", "Avoid the \\sloppy command.")
rule("ordinal_numbers",r"1st|2nd|3rd", "Spell out ordinal numbers (1st, 2nd, etc.) in words.")
rule("single_digit_spell",r"[a-z]+ \d [a-z]+", "Spell out single digit numbers in words.")
rule("single_space_period",r"\.  [A-Z]", "Prefer single space over double space after a period.")


text:
using ElectronicStructure

## using PyCall will trigger loading pyscf-specific code

using PyCall

using ElectronicStructure: one_electron_integrals,
   MolecularData, PySCF

##  Define geometries: specification of atoms and their positions
##  for various molecules.

geoms = (
    Geometry(Atom(:H, (0., 0., 0.)), Atom(:H, (0., 0., 0.7414))),

    Geometry(Atom(:Li, (0., 0., 0.)), Atom(:H, (0., 0., 1.4))),

    Geometry(Atom(:O, (0., 0., 0.)), Atom(:H, (0.757, 0.586, 0.)),
             Atom(:H, (-0.757, 0.586, 0.)))
)

## Choose one of the geometries
geom = geoms[1]

## Construct specification of electronic structure problem
mol_spec = MolecularSpec(geometry=geom)

## Do calculations and populate MolecularData with results
mol_pyscf = MolecularData(PySCF, mol_spec)

## Create interaction operator from one and two body integrals and constant.
## This does just a bit of manipulation of mol_data; converting space orbitals
## into space-and-spin orbitals.
## This is the same as the operator by the same name in OpenFermion.
iop = InteractionOperator(mol_pyscf; spin_order=:block, index_order=:default)

nothing


text:
using Statistics
using DataFrames
using CSV
using Plots; pyplot(fmt="png")
PyPlot.matplotlib.rc("mathtext", fontset="cm");
PyPlot.matplotlib.rc("text", usetex=true);
PyPlot.matplotlib.rc("font", family="serif", size=12);
using LaTeXStrings
using LsqFit

path = "/home/user/thesis/data/scaling/";
graph_type = "Network";
evolution_method = "stochastic_edge_acceptance";
data_file = string(path, graph_type, "_", evolution_method, ".csv");

savepath = "/home/user/thesis/latex/images"
save_bool = true
dpi = 144
if save_bool
	dpi = 300
end

# %% loading and computing df

df_raw = CSV.read(data_file) |> DataFrame;

Int.(df_raw[:, :N] .* df_raw[:, :r_1])
maximum(df_raw[:, :r_1])
cols = [
	:N,
	:r_0,
	:r_1,
	:m_0,
	:m_1,
	:Δr,
	:Δm,
];

df = aggregate(df_raw[:, cols], :N, [mean, std, median]);

# %% plotting Δr

x_min  = 15;
Δ_x    = 1;
x_max  = 24;
border = 0.25;
y_min  = 0.68;
Δ_y    = 0.01;
y_max  = 0.74;

plot_ = plot(legend=:right, dpi=dpi);
plot!(log2.(df[:, :N]), df[:, :r_0_mean], yerror=df[:, :r_0_std],
	label=L"t_0/N",
	line=(1, :black, 0.8),
	marker=(5, :dodgerblue, :circle, stroke(:dodgerblue)),
	xaxis=(latexstring("\\log_2 (N)"), (x_min-border, x_max+border), x_min:Δ_x:x_max),
	yaxis=(latexstring("t/N"), (y_min, y_max), y_min:Δ_y:y_max)
);
plot!(log2.(df[:, :N]), df[:, :r_1_mean], yerror=df[:, :r_1_std],
	label=L"t_1/N",
	line=(1, :black, 0.8),
	marker=(5, :firebrick1, :circle, stroke(:firebrick1)),
);

if save_bool
	savefig(plot_, joinpath(savepath, "r_scaling.png"))
end

df[:, :r_0_mean]
err = round((df[end, :r_1_mean] .- df[end, :r_0_mean]) / 2, digits=3)
round(df[end, :r_0_mean] + err, digits=3)

plot_

# %% plotting and fitting Δr

@. function f_linear(x, α)
	return α[1] * x + α[2]
end

α     = [-0.1, -1.5];
α_exp = [1., -1e6, 0.];
x     = log2.(df[:, :N]);
y     = log2.(df[:, :Δr_mean]);
σ     = df[:, :Δr_std];
fit   = curve_fit(f_linear, x, y, α);
y_fit = f_linear(x, fit.param);
α     = fit.param
σ     = stderror(fit)

plot_ = plot(dpi=dpi);
plot!(x, y_fit,
	label=latexstring("\\mathrm{Fit} \\sim $(round(α[1], digits=3)) ($(Int(round(σ[1] * 1e3, digits=0)))) \\cdot \\log_2(N)"),
	line=(1, :black, 0.8)
);
scatter!(x, y,
	label="Ensemble Average",
	markerstrokecolor=:black,
	marker=(6, :dodgerblue, :circle, stroke(:dodgerblue)),
	xaxis=(L"\log_2 (N)", 15:1:24),
	yaxis=(L"\log_2 (\Delta r)", (-8, -4.5))
);

if save_bool
	savefig(plot_, joinpath(savepath, "delta_r_scaling.png"))
end

plot_

# %%

@. function f_linear(x, α)
	return α[1] * x + α[2]
end

α     = [-0.1, -1.5];
α_exp = [1., -1e6, 0.];
x     = log2.(df[:, :N]);
y     = log2.(df[:, :Δm_mean]);
σ     = df[:, :Δm_std];
fit   = curve_fit(f_linear, x, y, α);
y_fit = f_linear(x, fit.param);
α     = fit.param
σ     = stderror(fit)

plot_ = plot(dpi=dpi);
plot!(x, y_fit,
	label=latexstring("\\mathrm{Fit} \\sim $(round(α[1], digits=3)) ($(Int(round(σ[1] * 1e3, digits=0)))) \\cdot \\log_2(N)"),
	line=(1, :black, 0.8)
);

δ = 0.302
τ = 2.213
exp_upper_bound = -δ * (τ - 2) / (τ - 1);
α_upper_bound = [exp_upper_bound, -exp_upper_bound * x[1] + y[1]]
plot!(x, f_linear(x, α_upper_bound),
	label=latexstring("\\mathrm{Upper Bound} \\sim $(round(exp_upper_bound, digits=3)) \\cdot \\log_2(N)"),
	line=(1, :black, :dash)
);

scatter!(x, y,
	label="Ensemble Average",
	markerstrokecolor=:black,
	marker=(6, :dodgerblue, :circle, stroke(:dodgerblue)),
	xaxis=(L"\log_2 (N)", 15:1:24),
	yaxis=(L"\log_2 (\Delta m)", (-3.8, -2), -3.8:0.2:-2)
);

if save_bool
	savefig(plot_, joinpath(savepath, "delta_m_scaling.png"))
end

plot_

max_stars_repo_name: cameronperot/explosive-percolation

max_issues_repo_name: cameronperot/explosive-percolation

max_forks_repo_name: cameronperot/explosive-percolation


text:
#using Plots
#
#function plotinvariants(H::HybridSystem, args...; kws...)
#    plot(first(H.invariants), args...; kws...)
#    for inv in Iterators.drop(H.invariants, 1)
#        plot!(inv, args...; kws...)
#    end
#end
#function plotinvariants!(H::HybridSystem, args...; kws...)
#    for inv in H.invariants
#        plot!(inv, args...; kws...)
#    end
#end

using FillArrays

ConeLyap(p::P, Q::Matrix{S}, b::Vector{S}, c, H) where {T, P<:AbstractPolynomial{T}, S} = ConeLyap{T, P, S}(p, Q, b, c, H)
JuMP.getvalue(p::ConeLyap) = ConeLyap(getvalue(p.p), getvalue(p.Q), getvalue(p.b), p.c, p.H)
ellipsoid(p::ConeLyap{T, P, JuMP.Variable}) where {T, P<:AbstractPolynomial{T}} = ellipsoid(getvalue(p))

function Ellipsoid(ell::LiftedEllipsoid)
    # P is
    # λ * [c'Qc-1  -c'Q
    #         -Qc   Q]
    # Let P be [β b'; b B]
    # We have
    # β = λ c'Qc - λ
    # b = -λ Qc <=> Q^{-1/2}b = -λ Q^{1/2} c
    # hence
    # λ c'Qc = β + λ
    # λ^2 c'Qc = b'Q^{-1}b = λ b'B^{-1}b <=> λ c'Qc = b'B^{-1}b
    # Hence λ = b'B^{-1}b - β
    n = LinAlg.checksquare(ell.P)-1
    ix = 1+(1:n)
    β = ell.P[1, 1]
    b = ell.P[1, ix]
    B = ell.P[ix, ix]
    λ = dot(b, B \ b) - β
    c = -(B \ b)
    Q = B / λ
    Ellipsoid(Q, c)
end

@recipe function f(ell::LiftedEllipsoid)
    Ellipsoid(ell)
end

function ellipsoid(p::ConeLyap)
    n = size(p.Q, 1)
    P = [-1. p.b'
         p.b  p.Q]
    HPH = p.H * P * p.H
    LiftedEllipsoid(inv(HPH))
end

function algebraiclift(s::DiscreteLinearControlSystem{T, MT, FullSpace}) where {T, MT}
    n = statedim(s)
    z = find(i -> iszero(sum(abs.(s.B[i,:]))), 1:n)
    # TODO ty - 1//2y^3 + 3//1xy + 2//1yhe affine space may not be parallel to classical axis
    DiscreteLinearAlgebraicSystem(s.A[z, :], (eye(n))[z, :])
end
algebraiclift(s::DiscreteIdentitySystem) = s
algebraiclift(S::AbstractVector) = algebraiclift.(S)
algebraiclift(S::Fill) = Fill(algebraiclift(first(S)), length(S))
function algebraiclift(h::HybridSystem)
    HybridSystem(h.automaton, algebraiclift(h.modes), h.invariants, h.guards, algebraiclift(h.resetmaps), h.switchings)
end

function r(A::Matrix{T}, c::Vector{T}=zeros(T, size(A, 1))) where T
    [one(T) zeros(T, 1, size(A, 2))
     c A]
end

function householder(x)
    y = copy(x)
    t = LinAlg.reflector!(y)
    v = [1; y[2:end]]
    eye(length(x)) - t * v * v'
end

using DynamicPolynomials
using MultivariatePolynomials
using PolyJuMP
using SumOfSquares
using JuMP
using LightGraphs
function ATrp(p, x, A)
    B = r(A)'
    y = x[1:size(B, 2)]
    p(x => r(A)' * y)
end
function lhs(p, x, Re::Fill)
    # If it is not constant, I would be comparing dummy variables of different meaning
    ATrp(p, x, first(Re).A)
end

function getp(m::Model, c, x, z::AbstractVariable)
    y = [z; x]
    n = length(x)
    β = 1.#@variable m lowerbound=0.
    b = @variable m [1:n]
    #@constraint m b .== 0
    Q = @variable m [1:n, 1:n] Symmetric
    @constraint m x' * Q * x in DSOSCone()
    H = householder([1; c]) # We add 1, for z
    P = [-β b'
         b  Q]
    HPH = H * P * H
    p = y' * HPH * y
    ConeLyap(p, Q, b, c, H)
    #@constraint m sum(Q) == 1 # dehomogenize
    #@variable m L[1:n, 1:n]
    #@variable m λinv[1:(n-1)] >= 0
    #@SDconstraint m [Q  L
    #                 L' diagm([λinv; -1])] ⪰ 0
    #ConeLyap(x' * Q * x, Q, L, λinv)
end
function getis(s::HybridSystem, solver, c)
    g = s.automaton.G
    n = nv(g)
    m = SOSModel(solver=solver)
    @polyvar x[1:2] z
    y = [z; x]
    l = [getp(m, c[u], x, z) for u in 1:n]
    #X = monomials(x, 2)
    #@variable m p[1:n] Poly(X)
    #@variable m vol
    #@objective m Max vol

    @objective m Max sum(p -> trace(p.Q), l)

    λouts = Vector{Vector{JuMP.Variable}}(n)
    #λouts = Vector{Vector{Float64}}(n)
    for u in 1:n
        # Constraint 1
        N = length(out_neighbors(g, u))
        λout = @variable m [1:N] lowerbound=0
        λouts[u] = λout
        #@constraint m sum(λin) == 1
        Σ = symbol.(s.automaton, Edge.(u, out_neighbors(g, u)))
        expr = lhs(l[u].p, y, s.resetmaps[Σ])
        for (j, v) in enumerate(out_neighbors(g, u))
            E = s.resetmaps[Σ[j]].E
            #expr -= λout[j] * p[v](x => r(E)' * x)
            newp = ATrp(l[v].p, y, E)
            expr -= λout[j] * newp
        end
        @constraint m expr in DSOSCone()
        # Constraint 2
        #@SDconstraint m differentiate(p[u], x, 2) >= 0
        # Constraint 3
        for hs in ineqs(s.invariants[u])
            @constraint m l[u].p(y => [-hs.β; hs.a]) <= 0
        end
    end

    status = solve(m)

    @show getobjectivevalue(m)

    for u in 1:n
        @show getvalue.(λouts[u])
    end

    @show status
    @assert status == :Optimal
    ellipsoid.(l)
#    if uc
#        [Ellipsoid(inv(Q[u].value), c[u] + 2*reshape(C[u].value, d)) for u in vertices(g)]
#    else
#        [Ellipsoid(inv(Q[u].value), c[u]) for u in vertices(g)]
#    end
end


text:
using Molly
using DelimitedFiles

data_dir = normpath(dirname(pathof(Molly)), "..", "data")
ff_dir = joinpath(data_dir, "force_fields")
openmm_dir = joinpath(data_dir, "openmm_6mrr")

ff = OpenMMForceField(joinpath.(ff_dir, ["ff99SBildn.xml", "tip3p_standard.xml", "his.xml"])...)

atoms, atoms_data, specific_inter_lists, general_inters, neighbor_finder, coords, box_size = setupsystem(
    joinpath(data_dir, "6mrr_equil.pdb"), ff)

n_steps = 500
timestep = 0.0005u"ps"
temp = 300.0u"K"
velocities = SVector{3}.(eachrow(readdlm(joinpath(openmm_dir, "velocities_300K.txt"))))u"nm * ps^-1"

s = Simulation(
    simulator=VelocityVerlet(),
    atoms=atoms,
    atoms_data=atoms_data,
    specific_inter_lists=specific_inter_lists,
    general_inters=general_inters,
    coords=coords,
    velocities=velocities,
    temperature=temp,
    box_size=box_size,
    neighbor_finder=neighbor_finder,
    thermostat=AndersenThermostat(1.0u"ps"),
    timestep=timestep,
    n_steps=n_steps,
)

parallel = false
simulate!(s, 5; parallel=parallel)
@time simulate!(s, n_steps; parallel=parallel)

parallel = true
simulate!(s, 5; parallel=parallel)
@time simulate!(s, n_steps; parallel=parallel)


text:
using DelimitedFiles

# Some file utils
function getFileExtension(filename)
    return filename[findlast(isequal('.'),filename):end]
end
function getFileHeader(filename)
    return filename[findlast(isequal('\\'),filename)+1:findlast(isequal('.'),filename)-1]
end

# Loads dataset points from a .txt file from the given path.
# Returns an array of points representing the dataset
function loadFromDelimTxt(filepath::String, k::Int64, idHeader::String)
    id = 0
    arrayPoints = Vector{Point}()
    io = open(filepath, "r")
    result = readdlm(filepath, '\t', Float64, '\n')

    for i in 1:size(result,1)
        push!(arrayPoints,Point(string(idHeader,id), result[i,:], k))
        id+=1
    end
    close(io)
    return arrayPoints
end

# Generate a .txt file containing a dataset starting from a matrix representing points. Used when input dataset is not provided
function generatePoints(filepath, matr)
    io = open(filepath, "w")

    for row in eachrow(matr)
        write(io, string(row[1], "\t", row[2]), '\n')
    end
    close(io)
end


text:
function backtrack_constrained(ϕ, α::Real, αmax::Real, αImax::Real,
                               Lcoefsα::Tuple{<:Real,<:Real,<:Real}, c1::Real = 0.5,
                               ρ::Real=oftype(α, 0.5),
                               αminfrac::Real = sqrt(eps(one(α)));
                               show_linesearch::Bool=false)
    # TODO: Specify that all elements should be of the same type T <: Real?
    # TODO: What does αI do??
    α, αI = min(α, 0.999*αmax), min(α, 0.999*αImax)
    αmin = αminfrac * α
    L0, L1, L2 = Lcoefsα
    if show_linesearch
        println("L0 = $L0, L1 = $L1, L2 = $L2")
    end
    while α >= αmin
        val = ϕ((α, αI))
        δ = evalgrad(L1, α, αI)
        if show_linesearch
            println("α = $α, αI = $αI, value: ($L0, $val, $(L0+δ))")
        end
        if isfinite(val) && val - (L0 + δ) <= c1*abs(val-L0)
            return α, αI
        end
        α *= ρ
        αI *= ρ
    end
    ϕ((zero(α), zero(αI)))  # to ensure that state gets set appropriately
    return zero(α), zero(αI)
end

function backtrack_constrained_grad(ϕ, α::Real, αmax::Real, Lcoefsα::Tuple{<:Real,<:Real,<:Real},
                                    c1::Real = 0.9, c2::Real = 0.9, ρ::Real=oftype(α, 0.5),
                                    αminfrac::Real = sqrt(eps(one(α))); show_linesearch::Bool=false)
    # TODO: Specify that all elements should be of the same type T <: Real?
    # TODO: Should c1 be 0.9 or 0.5 default?
    # TODO: Should ρ be 0.9 or 0.5 default?
    α = min(α, 0.999*αmax)
    αmin = αminfrac * α
    L0, L1, L2 = Lcoefsα
    if show_linesearch
        println("L0 = $L0, L1 = $L1, L2 = $L2")
    end
    while α >= αmin
        val, slopeα = ϕ(α)
        δval = L1*α
        δslope = L2*α
        if show_linesearch
            println("α = $α, value: ($L0, $val, $(L0+δval)), slope: ($L1, $slopeα, $(L1+δslope))")
        end
        if isfinite(val) && val - (L0 + δval) <= c1*abs(val-L0) &&
            (slopeα < c2*abs(L1) ||
             slopeα - (L1 + δslope) .<= c2*abs.(slopeα-L1))
            return α
        end
        α *= ρ
    end
    ϕ(zero(α))  # to ensure that state gets set appropriately
    return zero(α)
end

# Evaluate for a step parametrized as [α, α, αI, α]
function evalgrad(slopeα, α, αI)
    α*(slopeα[1] + slopeα[2] + slopeα[4]) + αI*slopeα[3]
end

# TODO: Never used anywhere? Intended for a linesearch that depends on ϕ''?
function mulhess(Hα, α, αI)
    αv = [α, α, αI, α]
    Hα*αv
end

# TODO: Never used anywhere? Intended for a linesearch that depends on ϕ''?
function evalhess(Hα, α, αI)
    αv = [α, α, αI, α]
    dot(αv, Hα*αv)
end


text:
#md # # Advanced Tutorial
#md # Notebook [download](advanced_tutorial.ipynb) [nbviewer](@__NBVIEWER_ROOT_URL__/advanced_tutorial.ipynb) [source](advanced_tutorial.jl)

# ## Introduction

# In this tutorial, we will fit the [coal mining disaster dataset](coal_mining_data.tsv) with a Gaussian process modulated
# Poisson process.

# ## Prepare the environment

# We start by importing:
# * MGVI that we will use for posterior fit
# * Distributions and FFTW to define the statistical model
# * Optim to pass Optim.Options to MGVI and to find Maximum A-Posteriori fit that we will use for comparison
# * StatsBase for histogram preparation from the data and also for error bands visualization
# * Plots.jl for visualization

using MGVI

using Distributions
using DelimitedFiles
using Random
using Optim
using StatsBase

using Plots

using FFTW

import ForwardDiff
#-
Random.seed!(84612);

# ## Load data

# Dataset is attached to the repository and contains intervals in days between
# disasters happend at British coal mines between March 1851 and March 1962.
# We split the entire time range into intervals of 365 days, then the number of events
# in each interval costitutes the measurement we are going to model.

function read_coal_mining_data(filepath, binsize)
    init_year = empty
    data = empty
    open(filepath) do io
        raw = readline(io)
        while ! occursin("init_date", raw)
            raw = readline(io)
        end

        init_year = parse(Float64, split(split(strip(raw[2:end]), "\t")[2], "-")[1])
        data = readdlm(io, '\t', Int, '\n', comments=true)[:]
    end
    dates_fract_years = init_year .+ cumsum(data)/365
    left_edge = dates_fract_years[1]
    num_bins = ((dates_fract_years[end] - left_edge) ÷ binsize)
    right_edge = left_edge + binsize*num_bins
    fit(Histogram, dates_fract_years, left_edge:binsize:right_edge).weights
end

coal_mine_disaster_data = read_coal_mining_data(joinpath(@__DIR__, "coal_mining_data.tsv"), 1);

# ## Global parameters and the grid

# Now we define several model properties:
# * `DATA_DIM` is just the size of the dataset
# * `DATA_XLIM` specifies the time range of the data
# * `GP_GRAIN_FACTOR` determines the numbers of finer bins into which the data bin is split.
#    This is useful when there are several datasets defined on different grids.
# * `GP_PADDING` adds empty paddings to the dataset. We use Fourier transform to sample from the Gaussian process
#    with a finite correlation length. `GP_PADDING` helps us to ensure that periodic boundary conditions
#    imposed by Fourier transform won't affect the data region.

DATA_DIM = size(coal_mine_disaster_data, 1);

data = coal_mine_disaster_data;

DATA_XLIM = [1851., 1962.];

GP_GRAIN_FACTOR = 3;
GP_PADDING = 80;
#-
function produce_bins()
    data_binsize = (DATA_XLIM[2] - DATA_XLIM[1])/DATA_DIM
    gp_binsize = data_binsize/GP_GRAIN_FACTOR
    gp_dim = Integer(((DATA_XLIM[2] - DATA_XLIM[1]) + 2*GP_PADDING) ÷ gp_binsize)
    gp_left_bin_offset = gp_right_bin_offset = (gp_dim - DATA_DIM) ÷ 2
    if (2*gp_left_bin_offset + DATA_DIM*GP_GRAIN_FACTOR) % 2 == 1
        gp_left_bin_offset += 1
    end
    gp_left_xlim = DATA_XLIM[1] - gp_left_bin_offset*gp_binsize
    gp_right_xlim = DATA_XLIM[2] + gp_right_bin_offset*gp_binsize
    gp_left_xs = collect(gp_left_xlim + gp_binsize/2:gp_binsize:DATA_XLIM[1])
    gp_right_xs = collect(DATA_XLIM[2] + gp_binsize/2:gp_binsize:gp_right_xlim)
    gp_data_xs = collect(DATA_XLIM[1] + gp_binsize/2:gp_binsize:DATA_XLIM[2])
    gp_xs = [gp_left_xs; gp_data_xs; gp_right_xs]
    data_idxs = collect(gp_left_bin_offset+1:GP_GRAIN_FACTOR:gp_left_bin_offset+DATA_DIM*GP_GRAIN_FACTOR)
    gp_xs, gp_binsize, data_idxs
end;

# Based on the defined model properties, we generate the grid. GP grid is the fine-grained grid
# with offsets added to the data range.
# * `_GP_XS` represent bin centers of such a fine-grained grid
# * `_GP_BINSIZE` is the width of the bin (that is 1/`GP_GRAIN_FACTOR` of data bin size)
# * `_DATA_IDXS` - integer indices of the left edges of the data bins

_GP_XS, _GP_BINSIZE, _DATA_IDXS = produce_bins();
_GP_DIM = length(_GP_XS);

# ## Model parameters

# The Gaussian process in this tutorial is modeled in the Fourier space with zero mean
# and two hyperparameters defining properties of its kernel. To sample from this
# Gaussian process, we also need a parameter per bin that will represent the particular
# realization of the GP in the bin.

function assemble_paridx(;kwargs...)
    pos = 0
    res = []
    for (k, v) in kwargs
        new_start, new_stop = v.start+pos, v.stop+pos
        push!(res, (k, (v.start+pos):(v.stop+pos)))
        pos = new_stop
    end
    (;res...)
end;

# MGVI is an iterative procedure, so we will need to introduce an initial guess for the state of the model.
# We create one vector of the size of the count of all parameters `starting_point` and a NamedTuple
# `PARDIX` that assigns names to the sub-regions in the vector of parameters. In the correct case:
# * `gp_hyper` two hyperparameters of the Gaussian process stored in the first two cells of the parameter vector
# * `gp_latent` `_GP_DIM` parameters used to define the particular realization of the gaussian process,
#    stored at indices between `3` to `2 + _GP_DIM`.
#
# Function `assemble_paridx` is responsible for constructing such a NamedTuple from the parameter specification.

PARIDX = assemble_paridx(gp_hyper=1:2, gp_latent=1:_GP_DIM);

starting_point = randn(last(PARIDX).stop);

# ## Model implementation

k = collect(0:(_GP_DIM)÷2 -1);

# Gaussian process covariance in the Fourier space is represented with a diagonal matrix. Values
# on the diagonal follow squared exponential function with parameters depending on priors.
# A kernel that is diagonal and mirrored around the center represents periodic and translation invariant function
# in the coordinate space. This property restricts covariance to have a finite correlation length in the coordinate
# space.
#
# MGVI assumes that all priors are distributed as standard normals `N(0, 1)`; thus,
# to modify the shapes of the priors, we explicitly rescale them at the model implementation phase.
#
# We also exponentiate each prior before using it for tuning squared exponential shape. In this way,
# we impose only positive values of the hyperparameters of the kernel.
#
# Actually, for the sake of numeric stability we model already square root of the covariance.
# This can be traced by missing `sqrt` in the next level, where we sample from the Gaussian process.

function kernel_model(p)
    kernel_A_c, kernel_l_c = p[PARIDX.gp_hyper]
    kernel_A = 60*exp(kernel_A_c*0.9)*GP_GRAIN_FACTOR
    kernel_l = 0.025*exp(kernel_l_c/15)/(GP_GRAIN_FACTOR^0.3)
    positive_modes = kernel_A .* sqrt(2 * π * kernel_l) .* exp.( -π^2 .* k.^2 .* kernel_l^2)
    negative_modes = positive_modes[end:-1:1]
    [positive_modes; negative_modes]
end;

# As a Fourier transform we use Discrete Hartley Transform that ensures that Fourier
# coefficients of the real valued function are real valued.

ht = FFTW.plan_r2r(zeros(_GP_DIM), FFTW.DHT);

# Before we proceed, let's have a brief look at the kernel's shape. Below
# we plot the kernel in the coordinate space `K(r) = K(x2 - x1)` as a function of time in years
# between two points. The further we go along `x`-axis the larger is the time interval and the
# smaller is the covariance between two points in time.

function plot_kernel_model(p, width; plot_args=(;))
    xs = collect(1:Int(floor(width/_GP_BINSIZE)))
    plot!(xs .* _GP_BINSIZE, (ht * kernel_model(p))[xs] ./ _GP_DIM, label=nothing, linewidth=2.5; plot_args...)
end

plot()
plot_kernel_model(starting_point, 20)

# To make it even more visual, we also plot the structure of the covariance matrix as a heatmap.
# We see that the finite correlation length shows up as a band around the diagonal. We also
# see small artifacts in the antidiagonal corners. They come from the assumption that the
# kernel is periodic.

function plot_kernel_matrix(p)
    xkernel = ht * kernel_model(p) ./ _GP_DIM
    res = reduce(hcat, [circshift(xkernel, i) for i in _GP_DIM-1:-1:0])'
    heatmap!(res)
end

plot()
plot_kernel_matrix(starting_point)

# After we defined the square root of the kernel function (`kernel_model`),
# we just follow the regular procedure of sampling from the normal distribution.
# Since the covariance matrix in the Fourier space is diagonal, Gaussian variables
# in each bin are independent of each other. Thus, sampling ends up rescaling
# the `gp_latent` part of the prior vector responsible for the Gaussian process state.
#
# After we produced a sample of Gaussian random values following the kernel model,
# we apply Fourier transform to return back to the coordinate space.

function gp_sample(p)
    flat_gp = kernel_model(p) .* p[PARIDX.gp_latent]
    (ht * flat_gp) ./ _GP_DIM
end;

# Together with the implementation of `gp_sample` we also need
# it's version for the `Dual`s. It is a little patch that makes
# application of the Hartley transform differentiatiable.

function gp_sample(dp::Vector{ForwardDiff.Dual{T, V, N}}) where {T,V,N}
    flat_gp_duals = kernel_model(dp) .* dp[PARIDX.gp_latent]
    val_res = ht*ForwardDiff.value.(flat_gp_duals) ./ _GP_DIM
    psize = size(ForwardDiff.partials(flat_gp_duals[1]), 1)
    ps = x -> ForwardDiff.partials.(flat_gp_duals, x)
    val_ps = map((x -> ht*ps(x) ./ _GP_DIM), 1:psize)
    ForwardDiff.Dual{T}.(val_res, val_ps...)
end;

# Gaussian process realization is meant to serve as a Poisson rate of the Poisson
# process. Since the Gaussian process is not restricted to positive values, we
# exponentiate its values to forcefully make the function positive.

function poisson_gp_link(fs)
    exp.(fs)
end;

# Now when we have a function representing the Poisson rate density,
# we have to integrate it over each data bin to define the Poisson rate in these bins.
# Function `agg_lambdas` does precisely that. When `GP_GRAIN_FACTOR = 1`, this function
# just multiplies the value of the Gaussian process in the bin by the `_GP_BINSIZE`.
# When we have more GP bins per data bin (`GP_GRAIN_FACTOR > 1`), then we apply
# rectangular quadrature to integrate over the bin.

function _forward_agg(data, idxs, steps_forward)
    [sum(data[i:i+steps_forward-1]) for i in idxs]
end;

function agg_lambdas(lambdas)
    gps = _forward_agg(lambdas, _DATA_IDXS, GP_GRAIN_FACTOR) .* _GP_BINSIZE
    xs = _GP_XS[_DATA_IDXS .+ (GP_GRAIN_FACTOR ÷ 2)]
    xs, gps
end;

# Finally, we arrive to the model definition assembled from the building blocks we defined above:
# * `gp_sample` sample from the Gaussian process with defined `kernel_model` covariance
# * `poisson_gp_link` ensures Gaussian process is positive
# * `agg_lambdas` integrates Gaussian process over each data bin to turn it into Poisson rate for each bin
# * `model` maps parameters into the product of the Poisson distributions counting events in each bin.

function model(params)
    fs = gp_sample(params)
    fine_lambdas = poisson_gp_link(fs)
    _, lambdas = agg_lambdas(fine_lambdas)
    Product(Poisson.(lambdas))
end;

# ## Visualization utilities

function agg_full_lambdas(lambdas)
    left_idxs = 1:GP_GRAIN_FACTOR:(_DATA_IDXS[1]-GP_GRAIN_FACTOR)
    left_gp = _forward_agg(lambdas, left_idxs, GP_GRAIN_FACTOR) .* _GP_BINSIZE
    left_xs = _GP_XS[left_idxs .+ (GP_GRAIN_FACTOR ÷ 2)]
    right_idxs = (_DATA_IDXS[end]+1):GP_GRAIN_FACTOR:(size(lambdas, 1) - GP_GRAIN_FACTOR)
    right_gp = _forward_agg(lambdas, right_idxs, GP_GRAIN_FACTOR) .* _GP_BINSIZE
    right_xs = _GP_XS[right_idxs .+ (GP_GRAIN_FACTOR ÷ 2)]
    middle_xs, middle_gp = agg_lambdas(lambdas)
    full_xs = [left_xs; middle_xs; right_xs]
    full_gp = [left_gp; middle_gp; right_gp]
    full_xs, full_gp
end;

function _mean(p; full=false)
    agg_func = if (!full) agg_lambdas else agg_full_lambdas end
    xs, gps = agg_func(poisson_gp_link(gp_sample(p)))
    xs, gps
end;

function plot_mean(p, label="mean"; plot_args=(;), full=false)
    plot!(_mean(p; full=full)..., label=label, linewidth=2; plot_args...)
end;

function plot_prior_samples(num_samples)
    for _ in 1:num_samples
        p = randn(last(PARIDX).stop)
        plot_mean(p, nothing)
    end
end;

function plot_kernel_prior_samples(num_samples, width)
    for _ in 1:num_samples
        p = randn(last(PARIDX).stop)
        plot_kernel_model(p, width)
    end
    plot!()
end;

function plot_data(; scatter_args=(;), smooth_args=(;))
    scatter!(_GP_XS[_DATA_IDXS .+ (GP_GRAIN_FACTOR ÷ 2)], data, la=0, markersize=2., markerstrokewidth=0, label="data"; scatter_args...)
    smooth_step = 4
    smooth_xs = _GP_XS[_DATA_IDXS .+ (GP_GRAIN_FACTOR ÷ 2)][(smooth_step+1):(end-smooth_step)]
    smooth_data = [sum(data[i-smooth_step:i+smooth_step])/(2*smooth_step+1) for i in (smooth_step+1):(size(data, 1)-smooth_step)]
    plot!(smooth_xs, smooth_data, linewidth=2, linealpha=1, ls=:dash, label="smooth data"; smooth_args...)
end;

function plot_mgvi_samples(params)
    for sample in eachcol(params.samples)
        if any(isnan.(sample))
            print("nan found in samples", "\n")
            continue
        end
        plot!(_mean(Vector(sample))..., linealpha=0.5, linewidth=1, label=nothing)
    end
    plot!()
end;

function plot_kernel_mgvi_samples(params, width)
    for sample in eachcol(params.samples)
        if any(isnan.(sample))
            print("nan found in samples", "\n")
            continue
        end
        plot_kernel_model(sample, width; plot_args=(linealpha=0.5, linewidth=1, label=nothing))
    end
    plot!()
end;

function produce_posterior_samples(p, num_residuals)
    batch_size = 10

    if num_residuals <= 2*batch_size
        batch_size = num_residuals ÷ 2
    end

    est_res_sampler = MGVI._create_residual_sampler(model, p;
                                                    residual_sampler=ImplicitResidualSampler,
                                                    jacobian_func=FwdRevADJacobianFunc,
                                                    residual_sampler_options=(;cg_params=(;abstol=1E-2)))
    batches = []
    for _ in 1:(num_residuals ÷ batch_size ÷ 2)
        batch_residual_samples = MGVI.rand(Random.GLOBAL_RNG, est_res_sampler, batch_size)
        push!(batches, p .+ batch_residual_samples)
        push!(batches, p .- batch_residual_samples)
    end
    reduce(hcat, batches)
end

function _extract_quantile(sorted_gp_realizations, p)
    map(s -> quantile(s, p; sorted=true), eachrow(sorted_gp_realizations))
end;

function plot_posterior_bands(p, num_samples; full=false)
    bands = [(0.997, :red), (0.955, :goldenrod1), (0.683, :green)]
    samples = produce_posterior_samples(p, num_samples)
    xs, first_gp = _mean(samples[1:end, 1]; full=full)
    gp_realizations = reduce(hcat, [_mean(Vector(sample); full=full)[2] for sample in eachcol(samples[1:end, 2:end])]; init=first_gp)
    for (i, one_x_sample) in enumerate(eachrow(gp_realizations))
        gp_realizations[i, 1:end] .= sort(Vector(one_x_sample))
    end
    for (band, color) in bands
        quant_l = _extract_quantile(gp_realizations, (1-band)/2)
        quant_u = _extract_quantile(gp_realizations, (1+band)/2)
        plot!(xs, quant_l; fillrange=quant_u, fillcolor=color, linealpha=0, label=band)
    end
    sample_median = _extract_quantile(gp_realizations, 0.5)
    plot!(xs, sample_median; linewidth=2, linecolor=:grey25, label="median")
end;

# ## Visualization and fitting

# We start by plotting the Gaussian process's dynamic range by
# sampling many possible realizations of it unconditionally on the data.
# We expect the set of lines to populate the regions where there are data.

plot()
plot_prior_samples(200)
plot_data()
plot!(ylim=[0, 8])

# We also plot prior samples for the kernel in the coordinate space. The plot below
# shows that the kernel is flexible in the amplitude while the correlation length
# is quite strongly predefined:

plot()
plot_kernel_prior_samples(200, 20)

# Now we see that the Gaussian process potentially is able to fit the data,
# we plot the initial guess (`starting_point`) to see where we start from.
# This plot shows:
# * data points
# * smoothed data with a moving average of 9 years
# * Poisson rate for each bin

plot()
plot_mean(starting_point, "starting_point")
plot_data()

# We also want to introduce the `full` plot that shows not only the data region
# but includes the region with the padding we added with `GP_PADDING`. We will use
# this plot to make sure that periodic boundary conditions don't interfere with
# the data.

plot()
plot_mean(starting_point, "full gp"; full=true)
plot_mean(starting_point, "starting_point")
plot_data()

# Let's make a first iteration of the MGVI. We limited Optim option to 1 iteration on purpose
# to let MGVI coverge slowly, so that we'll see nice convergence curve.

first_iteration = mgvi_kl_optimize_step(Random.GLOBAL_RNG,
                                        model, data,
                                        starting_point;
                                        num_residuals=3,
                                        jacobian_func=FwdRevADJacobianFunc,
                                        residual_sampler=ImplicitResidualSampler,
                                        optim_options=Optim.Options(iterations=1, show_trace=false),
                                        residual_sampler_options=(;cg_params=(;abstol=1E-2,verbose=false)));

# We again plot data and the Poisson rate. Then we again show the Gaussian process with padding.
# After one iteration the Poisson rate doesn't seem to get much closer to the data.

plot()
plot_mean(first_iteration.result, "first_iteration")
plot_data()
#-
plot()
plot_data()
plot_mean(first_iteration.result, "full gp"; full=true)
plot_mean(first_iteration.result, "first_iteration")

# We also would like to have a look at the kernel. Below we plot it together
# with the MGVI samples that represent the possible variation of the kernel
# shape around the mean:

plot()
plot_kernel_model(first_iteration.result, 20; plot_args=(;label="kernel model"))
plot_kernel_mgvi_samples(first_iteration, 20)

# In order to visualize convergence we prepare few functions to compute
# average posterior likelihood, store it and plot it.

function compute_avg_likelihood(model, samples, data)
    tot = 0
    for sample in eachcol(samples)
        tot += -MGVI.posterior_loglike(model, sample, data)
    end
    tot/size(samples, 2)
end;

function show_avg_likelihood(series)
    scatter!(1:size(series, 1), series, label="-loglike")
end;

# Now we do 30 more iterations of the MGVI and store the average likelihood after each step.
# We feed the fitted result of the previous step as an input to the next iteration.

next_iteration = first_iteration;
avg_likelihood_series = [];
push!(avg_likelihood_series, compute_avg_likelihood(model, next_iteration.samples, data));
for i in 1:30
    tmp_iteration = mgvi_kl_optimize_step(Random.GLOBAL_RNG,
                                          model, data,
                                          next_iteration.result;
                                          num_residuals=8,
                                          jacobian_func=FwdRevADJacobianFunc,
                                          residual_sampler=ImplicitResidualSampler,
                                          optim_options=Optim.Options(iterations=1, show_trace=false),
                                          residual_sampler_options=(;cg_params=(;abstol=1E-2,verbose=false)))
    global next_iteration = tmp_iteration
    push!(avg_likelihood_series, compute_avg_likelihood(model, next_iteration.samples, data))
end;

# Firstly, let's have a look at the convergence plots. We see that MGVI converged after 10 iterations
# while being limited to very poor Optim performance.

plot(yscale=:log)
show_avg_likelihood(avg_likelihood_series)

# Below we plot the result of the fit. Together with data and Poisson rate we also plot
# MGVI residuals. They are samples from the Gaussian posterior, sampled respecting the posterior
# covariance structure. Thus MGVI residual samples are deviations from the MGVI fit that represent
# how confident we are about the prediction.

plot(ylim=[0,8])
plot_mgvi_samples(next_iteration)
plot_mean(next_iteration.result, "many_iterations", plot_args=(color=:deepskyblue2, linewidth=3.5))
plot_data(scatter_args=(;color=:blue2, marker_size=3.5), smooth_args=(;color=:deeppink3, linewidth=3))

# To present credibility intervals we also plot credibility bands. We sample 400 residual samples
# from MGVI and then plot quantiles for each data bin. This should give us a feeling of how
# the MGVI fit is compatible with the data.

plot(ylim=[0,8])
plot_posterior_bands(next_iteration.result, 400)
plot_mean(next_iteration.result, "many_iterations", plot_args=(color=:deepskyblue2, linewidth=3.5))
plot_data(scatter_args=(;color=:blue2, marker_size=3.5), smooth_args=(;color=:deeppink3, linewidth=3))

# We also make sure boundary conditions do not interfere with the data. Here is the Gaussian process
# plot with the paddings included:

plot()
plot_data()
plot_mean(next_iteration.result; full=true)
plot_mean(next_iteration.result, "many_iterations")

# Let's have a look at the kernel again. We expect the variation of samples
# to become narrower:

plot()
plot_kernel_model(next_iteration.result, 20; plot_args=(;label="kernel model"))
plot_kernel_mgvi_samples(next_iteration, 20)

# ## Maximum A-Posteriori estimation

# We build a MAP as a cross check of MGVI results. We just optimize posterior likelihood with Optim
# without any particular tuning:

max_posterior = Optim.optimize(x -> -MGVI.posterior_loglike(model, x, data), starting_point, LBFGS(), Optim.Options(show_trace=false, g_tol=1E-10, iterations=300));

# We observe that the bump in the middle (around 1910) is caught by MAP while it is less pronounced in the MGVI fit.
# MAP also has finer structure around 1875 and 1835.

plot()
plot_mean(Optim.minimizer(max_posterior), "map")
plot_mean(next_iteration.result, "mgvi mean")
plot_data()

# We also can see the difference at the left edge of the data region. While MGVI smoothed the data,
# MAP predicted a consequent peak:

plot()
plot_data()
plot_mean(Optim.minimizer(max_posterior), "full gp"; full=true)
plot_mean(next_iteration.result, "mgvi full gp"; full=true)


text:
function rectangle(w, h)
    x = -w/2
    y = -h/2
    return Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
end

prior_bounds = [800 800; 500 1000; 500 1000; 800 800;
                1500 1500; 1500 1500; 1000 1500; 1000 1500;
                500 1200; 1000 1000; 1000 1000; 500 1200
]

function add_prior_region!(prior_bounds;label=true)

    rect = prior_bounds*1E-3
    if label==true
        plot!(rectangle(rect...), opacity=.5, c=:gray, linecolor=:white, label="Prior", legend=:outertop)
    else
        plot!(rectangle(rect...), opacity=.5, c=:gray, linecolor=:white, label=false, legend=:outertop)
    end

end

function add_trueloc!(true_loc, label, char)

    dis_lab = get_dis_lab(char, phi, chi)
    true_loc_im = dis_to_image(true_loc, dis_lab)
    num_dislocs = Int(length(true_loc)/2)
    for i = 1:num_dislocs
        scatter!([true_loc_im[2*i-1]*1E-3],[true_loc_im[2*i]*1E-3], c=:cyan, markersize = 8, markershape = :star7, legend=:outertop, legendfontsize=18, label=label)
    end

end

function dfxm_image_plot_base(im, char, true_loc, save_str; width=450, true_loc_label = false, tickfontsize=16, prior=false, prior_label=false)

    heatmap(y_im_vec, z_im_vec, im, clim=(0,maximum(im)), aspect_ratio=:equal, bbox=:tight)

    if prior
        add_prior_region!(prior_bounds[char,:], label=prior_label)
    end

    add_trueloc!(true_loc,  true_loc_label, char)

    xlims!(y_im_vec[1],y_im_vec[end])
    ylims!(z_im_vec[1],z_im_vec[end])

    xticks!([-1,0,1])
    yticks!([-4,-2,0,2,4])

    xlabel!(L"y_{\ell} \ (\mu \mathrm{m})", xtickfontsize=tickfontsize, xguidefontsize=18, legendfontsize=18, margin=5mm)
    ylabel!(L"x_{\ell} \ (\mu \mathrm{m})", ytickfontsize=tickfontsize, yguidefontsize=18)

    plot!(size=(width,width))
    Plots.savefig(save_str*".pdf")

end

max_stars_repo_name: michael-c-brennan/DFXMTools

max_issues_repo_name: michael-c-brennan/DFXMTools

max_forks_repo_name: michael-c-brennan/DFXMTools


text:
# Version that takes `op("S=1/2", "X", Index(2, "i"))`.
# Shortcuts searching for the correct tag, so is a bit faster.
ITensors.op(on::AbstractString,
            st::AbstractString,
            s::Index...; kwargs...) = op(OpName(on), SiteType(st), s...; kwargs...)


text:
using Pkg

Pkg.add(PackageSpec(url="https://github.com/danphenderson/TRS.jl.git", rev="8f5af784ee98d72df050fb03deadd88cd41d1209"))


text:
@testset "mean_functions" begin
    rng = MersenneTwister(123456)
    N, D = 5, 3
    x1 = randn(rng, N)
    xD_colvecs = ColVecs(randn(rng, D, N))
    xD_rowvecs = RowVecs(randn(rng, N, D))

    zero_mean_testcase = (; mean_function=ZeroMean(), calc_expected=_ -> zeros(N))

    c = randn(rng)
    const_mean_testcase = (; mean_function=ConstMean(c), calc_expected=_ -> fill(c, N))

    foo_mean = x -> sum(abs2, x)
    custom_mean_testcase = (;
        mean_function=CustomMean(foo_mean), calc_expected=x -> map(foo_mean, x)
    )

    @testset "$(typeof(testcase.mean_function))" for testcase in [
        zero_mean_testcase, const_mean_testcase, custom_mean_testcase
    ]
        for x in [x1, xD_colvecs, xD_rowvecs]
            m = testcase.mean_function
            @test AbstractGPs._map_meanfunction(m, x) == testcase.calc_expected(x)
            differentiable_mean_function_tests(rng, m, x)
        end
    end
end

max_stars_repo_name: JuliaGaussianProcesses/AbstractGP

max_issues_repo_name: JuliaGaussianProcesses/AbstractGP

max_forks_repo_name: JuliaGaussianProcesses/AbstractGP


text:
using CSV, SparseArrays, DataFrames
using Recommenders: rows2sparse, reindex_id_column, make_u2i_dataset
using Test

@testset "CSV to sparse matrix" begin
    csv_str = """
    col1,col2,col3
    1,1,0.1
    1,2,-1
    2,2,5
    """
    csv = CSV.File(IOBuffer(csv_str))
    mat_sparse = rows2sparse(csv, col_user = :col1, col_item = :col2, col_rating = :col3)
    expected = sparse([1, 1, 2], [1, 2, 2], [0.1, -1, 5])
    @test mat_sparse == expected
end

@testset "DataFrame to sparse matrix" begin
    df = DataFrame(col1 = [1, 1, 2], col2 = [1, 2, 2], col3 = [0.1, -1, 5])
    mat_sparse = rows2sparse(df, col_user = :col1, col_item = :col2, col_rating = :col3)
    expected = sparse([1, 1, 2], [1, 2, 2], [0.1, -1, 5])
    @test mat_sparse == expected
end

@testset "Re-index id column" begin
    df = DataFrame(id = [1, 4, 2, 1])
    table, id2index = reindex_id_column(df, :id)

    expected_table = DataFrame(id = [1, 2, 3, 1])
    expected_mapper = Dict(1 => 1, 4 => 2, 2 => 3)

    @test table == expected_table
    @test id2index == expected_mapper
end

@testset "Make u2i dataset" begin
    df = DataFrame(userid = [1, 1, 2, 3], itemid = [3, 4, 4, 1])
    xs, ys = make_u2i_dataset(df)
    sorted_idx = sortperm(xs)
    xs = xs[sorted_idx]
    ys = ys[sorted_idx]

    @test xs == [1, 2, 3]
    @test ys == [[3, 4], [4], [1]]
end


text:
export save, load

# import JLD2
using DelimitedFiles

"""
    save(Π, "file.csv") = writecsv("file.csv", Π)
    save(Π, "file.json")
    save(Π, "file.jld") = JLD.save(file, "Π", Π)
Save `Π::Weighted` to disk, in one of several formats:

* CSV: weights are rightmost column, as in `[Π.array' Π.weights]`.
  Options aren't saved, and `load()` will try to guess them from the numbers.
  Mathematica: `arr = Import["file.csv"][[;; , 1 ;; -2]]; wei = Import["file.csv"][[;; , -1]]`

* JSON: saves a dictionary, but `dict["array"]` is a nested Vector{Vector{Any}} for now,
  which `load()` converts to Float64. Includes `dict["opt"] == String(Π.opt)`.
  Mathematica: `arr = Import["file.json"][[1, 2]]; wei = Import["file.json"][[3, 2]]`
  ... and then maybe `Style @@@ Thread[{arr, PointSize /@ (0.3*Sqrt[wei]), Opacity[0.5]}] // ListPlot` ?

* JLD: built-in HDF5 format binary saving, fast & neat but perhaps fragile.
  Removed for now!

    Π = load("file.csv") = readcsv("file.csv", Π)
    Π = load("file.json")
    Π = load("file.jld") = JLD.load(file, "Π")
Reverse the above.
"""
function save(x::Weighted, file::String; verbose=true)
    if endswith(file, ".csv")
        writecsv(file, x)
    elseif endswith(file, ".json")
        writejson(file, x)
    # elseif endswith(file, ".jld")
    #     JLD.save(file, "Π", x)
    else
        error("load doesn't understand file extension of $file")
    end
end

save(file::String, x::Weighted; kw...) = save(x, file; kw...)

@doc @doc(save)
function load(file::String)
    if endswith(file, ".csv")
        readcsv(file, Weighted)
    elseif endswith(file, ".json")
        readjson(file)
    # elseif endswith(file, ".jld")
    #     JLD.load(file, "Π")
    else
        error("save doesn't understand file extension of $file")
    end
end

# writecsv(io, x::Weighted; kw...) = writedlm(io, [x.array' x.weights], ','; kw...)
writecsv(io, x::Weighted; kw...) = writedlm(io, [x.array' |> copy x.weights], ','; kw...) ## copy for Flux bug with Adjoint

function readcsv(io, T::Type{Weighted}; kw...)

    mat = readdlm(io, ','; kw...)
    array = mat[:, 1:end-1]' |> copy
    weights = mat[:,end]

    norm = sum(weights)≈1
    if minimum(array)>=0 && maximum(array)<=1
        clamp,lo,hi = true,0,1
    else
        clamp,lo,hi = false,-Inf,Inf
    end
    if endswith(string(io), ".csv")
        aname = string(io)[1:end-4] |> Symbol
    else
        aname = string(io) |> Symbol
    end

    Weighted(array, weights, WeightOpt(norm=norm, clamp=clamp,lo=lo,hi=hi, aname=aname))
end

# Base.writecsv(io, x::Weighted; kw...) = writecsv(io, (x.array, x.weights, [x.opt.norm, x.opt.clamp, x.opt.lo, x.opt.lo, x.opt.names]); kw...)

import JSON

function writejson(io, x::Weighted)
    d = Dict("array" => x.array, "weights" => x.weights, "opt" => string(x.opt))
    write(io, JSON.json(d))
end

function readjson(io)
    d = JSON.parse(read(io, String))

    array = reduce(hcat, d["array"]) |> Matrix{Float64}
    weights = d["weights"] |> Vector{Float64}
    opt = eval(Meta.parse(d["opt"]))

    Weighted(array, weights, opt)
end


text:
module Day15

using DataStructures

import ..data_dir # from parent module
input = read(joinpath(data_dir, "day15"), String)

export part1, part2

function string_block_to_matrix(string_block::AbstractString)
    as_arrays = split.(split(string_block, "\n"), "")
    as_int_arrays = [parse.(Int64, array) for array in as_arrays]
    # Transposing changes the type!
    Matrix(transpose(hcat(as_int_arrays...)))
end

function out_of_bounds(heatmap::Matrix{Int64}, position::Tuple{Int64,Int64})
    row, column = position
    heatmap_height, heatmap_width = size(heatmap)
    row < 1 || row > heatmap_height || column < 1 || column > heatmap_width
end

function neighbours(heatmap::Matrix{Int64}, position::Tuple{Int64,Int64})
    row, column = position
    candidates = [
        (row - 1, column), (row + 1, column),
        (row, column - 1), (row, column + 1)
    ]
    filter(candidate -> !out_of_bounds(heatmap, candidate), candidates)
end

function dijkstra(heatmap::Matrix{Int64})
    source = (1, 1) # top-left corner
    target = size(heatmap) # bottom-right corner

    risk_from_source = BinaryMinHeap{Tuple{Int64,Int64,Int64}}()
    push!(risk_from_source, (0, source...))

    counted = Vector{Tuple{Int64,Int64}}()

    while true
        risk, row, column = pop!(risk_from_source)
        vertex = (row, column)
        if vertex in counted
            continue
        end
        push!(counted, vertex)

        if vertex == target
            return risk
        end

        uncounted_neighbours = filter(v -> v ∉ counted, neighbours(heatmap, vertex))

        for neighbour in uncounted_neighbours
            neighbour_risk = risk + heatmap[neighbour...]
            push!(risk_from_source, (neighbour_risk, neighbour...))
        end
    end
end

function part1(input = input)
    heatmap = string_block_to_matrix(input)
    dijkstra(heatmap)
end

function increment(heatmap::Matrix{Int64}, times::Int64 = 1)
    for i = 1:times
        heatmap = @. max(mod(heatmap + 1, 10), 1)
    end
    heatmap
end

function expand_heatmap(heatmap::Matrix{Int64}, times::Int64)
    horizontal_expansion_vector = Vector{Matrix{Int64}}()
    for i = 1:times
        push!(horizontal_expansion_vector, increment(heatmap, i - 1))
    end
    horizontal_expansion = hcat(horizontal_expansion_vector...)

    vertical_expansion_vector = Vector{Matrix{Int64}}()
    for i = 1:times
        push!(vertical_expansion_vector, increment(horizontal_expansion, i - 1))
    end
    vertical_expansion = vcat(vertical_expansion_vector...)

    vertical_expansion
end

function part2(input = input)
    heatmap = string_block_to_matrix(input)
    expanded_heatmap = expand_heatmap(heatmap, 5)
    dijkstra(expanded_heatmap)
end

end #module


text:
module MultimodeCavity

include("potential_box.jl")

module quantum
include("system_quantum.jl")
include("timeevolution_quantum.jl")
include("steadystate_quantum.jl")
end

module classical
include("system_classical.jl")
include("timeevolution_classical.jl")
include("steadystate_classical.jl")
end

end # module


text:
module AbstractInstances

_singletonsymbol(T::UnionAll) = _singletonsymbol(T.body)
_singletonsymbol(T::DataType) = Symbol("Singleton", string(T))

_new_singleton(T::UnionAll) = _new_singleton(T.body)
function _new_singleton(T::DataType)
    S = _singletonsymbol(T)
    return quote
        struct $S <: $T end
    end
end

_concretesymbol(T::UnionAll) = _concretesymbol(T.body)
_concretesymbol(T::DataType) = Symbol("Concrete", string(T))

_new_concrete(T::UnionAll) = _new_concrete(T.body)
function _new_concrete(T::DataType)
    S = _concretesymbol(T)
    return quote
        mutable struct $S <: $T end
    end
end

"""
    singletontype(T::Type)

Return a concrete immutable type with no fields which is a subtype of
`T`.  The immutability of the returned type ensures that instances of
the type are singletons.
"""
function singletontype(T::Type)
    @assert isabstracttype(T)
    S = _singletonsymbol(T)
    if !(@eval @isdefined $S)
        @eval $(_new_singleton(T))
    end
    return @eval $S
end

"""
    singleton(T::Type)

Create an instance of type, `singletontype(T)`.
"""
singleton(T::Type) = @eval singletontype($T)()

"""
    concretetype(T::Type)

Return a concrete mutable type with no fields which is a subtype of
`T`.  The mutability of the returned type ensures that instances of
the type are unique.
"""
function concretetype(T::Type)
    @assert isabstracttype(T)
    S = _concretesymbol(T)
    if !(@eval @isdefined $S)
        @eval $(_new_concrete(T))
    end
    return @eval $S
end

"""
    oftype(T::Type)

Create a unique instance of type, `concretetype(T)`.
"""
oftype(T::Type) = @eval concretetype($T)()

end # module


text:
branches(::Type{AbstractBinomialTree}) = 2

mutable struct BinomialTreeCommon
  x0::Float64
  driftPerStep::Float64
  dt::Float64
  columns::Int
end

get_x0(t::AbstractBinomialTree) = t.common.x0
get_drift_per_step(t::AbstractBinomialTree) = t.common.driftPerStep
get_dt(t::AbstractBinomialTree) = t.common.dt

function BinomialTreeCommon(process::StochasticProcess1D, endTime::Float64, steps::Int)
  x0 = get_x0(process)
  dt = endTime/steps
  driftPerStep = drift(process, 0.0, x0) * dt
  columns = steps + 1
  BinomialTreeCommon(x0, driftPerStep, dt, columns)
end

struct Joshi4 <: BinomialTreeType end
struct LeisenReimer <: BinomialTreeType end
struct Tian <: BinomialTreeType end

mutable struct BinomialTree{T <: BinomialTreeType} <: AbstractBinomialTree
  up::Float64
  down::Float64
  pu::Float64
  pd::Float64
  common::BinomialTreeCommon
  treeType::T
end

const Joshi4BinomialTree = BinomialTree{Joshi4}
const LeisenReimerBinomialTree = BinomialTree{LeisenReimer}
const TianBinomialTree = BinomialTree{Tian}

function Tian(process::StochasticProcess1D, endTime::Float64, steps::Int, ::Float64)
  common = BinomialTreeCommon(process, endTime, steps)
  q = exp(variance(process, 0.0, common.x0, common.dt))
  r = exp(common.driftPerStep) * sqrt(q)

  up = 0.5 * r * q * (q + 1.0 + sqrt(q * q + 2.0 * q - 3.0))
  down = 0.5 * r * q * (q + 1.0 - sqrt(q * q + 2.0 * q - 3.0))

  pu = (r - down) / (up - down)
  pd = 1.0 - pu

  pu <= 1.0 || error("negative probability")
  pu >= 0.0 || error("negative probability")

  return BinomialTree{Tian}(up, down, pu, pd, common, Tian())
end

function LeisenReimer(process::StochasticProcess1D, endTime::Float64, steps::Int, strike::Float64)
  common = BinomialTreeCommon(process, endTime, steps)
  strike > 0.0 || error("strike must be positive")
  oddSteps = isodd(steps) ? steps : steps + 1
  variance_ = variance(process, 0.0, common.x0, endTime)
  ermqdt = exp(common.driftPerStep + 0.5 * variance_ / oddSteps)
  d2 = (log(common.x0 / strike) + common.driftPerStep * oddSteps) / sqrt(variance_)

  pu = peizer_pratt_method_2_inversion(d2, oddSteps)
  pd = 1.0 - pu
  pdash = peizer_pratt_method_2_inversion(d2 + sqrt(variance_), oddSteps)

  up = ermqdt * pdash / pu
  down = (ermqdt - pu * up) / (1.0 - pu)

  return BinomialTree{LeisenReimer}(up, down, pu, pd, common, LeisenReimer())
end

function compute_up_prob(k::Float64, dj::Float64)
  alpha = dj / sqrt(8.0)
  # alpha2 = alpha * alpha
  alpha3 = alpha ^ 3
  alpha5 = alpha ^ 5
  alpha7 = alpha ^ 7
  beta = -0.375 * alpha - alpha3
  gamma = (5.0 / 6.0) * alpha5 + (13.0/12.0) * alpha3 + (25.0 / 128.0) * alpha
  delta = -0.1025 * alpha - 0.9285 * alpha3 - 1.43 * alpha5 - 0.5 * alpha7
  p = 0.5
  rootk = sqrt(k)
  p += alpha / rootk
  p += beta / (k * rootk)
  p += gamma / (k * k * rootk)
  # delete next line to get results for j-three tree
  p += delta / (k * k * k * rootk)
  return p
end

function Joshi4(process::StochasticProcess1D, endTime::Float64, steps::Int, strike::Float64)
  common = BinomialTreeCommon(process, endTime, steps)
  strike > 0.0 || error("strike must be positive")
  oddSteps = isodd(steps) ? steps : steps + 1
  variance_ = variance(process, 0.0, common.x0, endTime)
  ermqdt = exp(common.driftPerStep + 0.5 * variance_ / oddSteps)
  d2 = (log(common.x0 / strike) + common.driftPerStep * oddSteps) / sqrt(variance_)

  pu = compute_up_prob((oddSteps - 1.0) / 2.0, d2)
  pd = 1.0 - pu
  pdash = compute_up_prob((oddSteps - 1.0) / 2.0, d2 + sqrt(variance_))

  up = ermqdt * pdash / pu
  down = (ermqdt - pu * up) / (1.0 - pu)

  return BinomialTree{Joshi4}(up, down, pu, pd, common, Joshi4())
end

struct JarrowRudd <: EqualProbabilitiesBinomialTreeType end
struct AdditiveEQP <: EqualProbabilitiesBinomialTreeType end

mutable struct EqualProbabilitiesBinomialTree{T <: EqualProbabilitiesBinomialTreeType} <: AbstractBinomialTree
  up::Float64
  common::BinomialTreeCommon
  treeType::T
end

const JarrowRuddBinomialTree = EqualProbabilitiesBinomialTree{JarrowRudd}
const AdditiveEQPBinomialTree = EqualProbabilitiesBinomialTree{AdditiveEQP}

# overloaded to build this alias
function JarrowRudd(process::StochasticProcess1D, endTime::Float64, steps::Int, ::Float64)
  common = BinomialTreeCommon(process, endTime, steps)
  up = std_deviation(process, 0.0, common.x0, common.dt)

  return EqualProbabilitiesBinomialTree{JarrowRudd}(up, common, JarrowRudd())
end

function AdditiveEQP(process::StochasticProcess1D, endTime::Float64, steps::Int, ::Float64)
  common = BinomialTreeCommon(process, endTime, steps)
  up = -0.5 * common.driftPerStep + 0.5 * sqrt(4.0 * variance(process, 0.0, common.x0, common.dt) - 3.0 * common.driftPerStep * common.driftPerStep)

  return EqualProbabilitiesBinomialTree{AdditiveEQP}(up, common, AdditiveEQP())
end

struct CoxRossRubinstein <: EqualJumpsBinomialTreeType end
struct Trigeorgis <: EqualJumpsBinomialTreeType end

mutable struct EqualJumpsBinomialTree{T <: EqualJumpsBinomialTreeType} <: AbstractBinomialTree
  dx::Float64
  pu::Float64
  pd::Float64
  common::BinomialTreeCommon
  treeType::T
end

const CoxRossRubinsteinBinomialTree = EqualJumpsBinomialTree{CoxRossRubinstein}
const TrigeorgisBinomialTree = EqualJumpsBinomialTree{Trigeorgis}

function CoxRossRubinstein(process::StochasticProcess1D, endTime::Float64, steps::Int, ::Float64)
  common = BinomialTreeCommon(process, endTime, steps)
  dx = std_deviation(process, 0.0, common.x0, common.dt)
  pu = 0.5 + 0.5 * common.driftPerStep / dx
  pd = 1.0 - pu

  pu <= 1.0 || error("negative probability")
  pu >= 0.0 || error("negative probability")

  return EqualJumpsBinomialTree{CoxRossRubinstein}(dx, pu, pd, common, CoxRossRubinstein())
end

function Trigeorgis(process::StochasticProcess1D, endTime::Float64, steps::Int, ::Float64)
  common = BinomialTreeCommon(process, endTime, steps)
  dx = sqrt(variance(process, 0.0, common.x0, common.dt) + common.driftPerStep * common.driftPerStep)
  pu = 0.5 + 0.5 * common.driftPerStep / dx
  pd = 1.0 - pu

  pu <= 1.0 || error("negative probability")
  pu >= 0.0 || error("negative probability")

  return EqualJumpsBinomialTree{Trigeorgis}(dx, pu, pd, common, Trigeorgis())
end

get_size(tree::AbstractBinomialTree, idx::Int) = idx

descendant(tree::AbstractBinomialTree, ::Int, idx::Int, branch::Int) = (idx - 1) + (branch - 1) + 1

function get_underlying(tree::EqualProbabilitiesBinomialTree, i::Int, idx::Int)
  j = 2 * (idx - 1) - (i - 1)
  # exploiting the forward value tree centering
  return get_x0(tree) * exp((i-1) * get_drift_per_step(tree) + j * tree.up)
end
function get_underlying(tree::EqualJumpsBinomialTree, i::Int, idx::Int)
  j = 2 * (idx - 1) - (i - 1)
  # exploting equal jump and the x0 tree centering
  return get_x0(tree) * exp(j * tree.dx)
end
get_underlying(tree::BinomialTree, i::Int, idx::Int) = get_x0(tree) * ^(tree.down, (i - 1) - (idx - 1)) * ^(tree.up, (idx - 1))

probability(tree::EqualProbabilitiesBinomialTree, ::Int, ::Int, ::Int) = 0.5
probability(tree::Union{EqualJumpsBinomialTree, BinomialTree}, ::Int, ::Int, branch::Int) = branch == 2 ? tree.pu : tree.pd


text:
#==============================================================================#
# Translate.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/translate-2017-07-01.normal.json
#==============================================================================#

__precompile__()

module Translate

using AWSCore

"""
    using AWSSDK.Translate.translate_text
    translate_text([::AWSConfig], arguments::Dict)
    translate_text([::AWSConfig]; Text=, SourceLanguageCode=, TargetLanguageCode=)

    using AWSCore.Services.translate
    translate([::AWSConfig], "TranslateText", arguments::Dict)
    translate([::AWSConfig], "TranslateText", Text=, SourceLanguageCode=, TargetLanguageCode=)

# TranslateText Operation

Translates input text from the source language to the target language. You can translate between English (en) and one of the following languages, or between one of the following languages and English.

*   Arabic (ar)

*   Chinese (Simplified) (zh)

*   French (fr)

*   German (de)

*   Portuguese (pt)

*   Spanish (es)

To have Amazon Translate determine the source language of your text, you can specify `auto` in the `SourceLanguageCode` field. If you specify `auto`, Amazon Translate will call Amazon Comprehend to determine the source language.

# Arguments

## `Text = ::String` -- *Required*
The text to translate.

## `SourceLanguageCode = ::String` -- *Required*
One of the supported language codes for the source text. If the `TargetLanguageCode` is not "en", the `SourceLanguageCode` must be "en".

To have Amazon Translate determine the source language of your text, you can specify `auto` in the `SourceLanguageCode` field. If you specify `auto`, Amazon Translate will call Amazon Comprehend to determine the source language.

## `TargetLanguageCode = ::String` -- *Required*
One of the supported language codes for the target text. If the `SourceLanguageCode` is not "en", the `TargetLanguageCode` must be "en".

# Returns

`TranslateTextResponse`

# Exceptions

`InvalidRequestException`, `TextSizeLimitExceededException`, `TooManyRequestsException`, `UnsupportedLanguagePairException`, `DetectedLanguageLowConfidenceException`, `InternalServerException` or `ServiceUnavailableException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/translate-2017-07-01/TranslateText)
"""
@inline translate_text(aws::AWSConfig=default_aws_config(); args...) = translate_text(aws, args)

@inline translate_text(aws::AWSConfig, args) = AWSCore.Services.translate(aws, "TranslateText", args)

@inline translate_text(args) = translate_text(default_aws_config(), args)

end # module Translate

#==============================================================================#
# End of file
#==============================================================================#


text:
module Millboard

include("types.jl")
include("table_modes.jl")

export table
include("table.jl")

end # module


text:
export pairwise_marginal, pairwise_MI

using LinearAlgebra: diagind, diag
using CUDA

xlogx(x) =
    iszero(x) ? zero(x) : x * log(x)

xlogy(x, y) =
    iszero(x) && !isnan(y) ? zero(x) : x * log(y)

#############################
# Mutual Information on Binary Data
#############################

function pairwise_marginal(data::Union{Matrix{Bool}, CuMatrix{Bool}, BitMatrix};
        weights::Union{Vector, CuVector, Nothing}=nothing,
        pseudocount=0.0,
        Float=Float64)
    N = isnothing(weights) ? size(data, 1) : sum(weights)
    D = size(data, 2)
    base = N + pseudocount
    if data isa CuMatrix{Bool} && !isnothing(weights)
        weights = CuVector(weights)
    end

    not_data = Float.(.!data)
    data = Float.(data)
    pxy = similar(data, D, D, 4)

    if isnothing(weights)
        pxy[:, :, 4] = transpose(data) * data
        pxy[:, :, 3] = transpose(data) * not_data
        pxy[:, :, 2] = transpose(not_data) * data
    else
        pxy[:, :, 4] = transpose(data) * (data .* weights)
        pxy[:, :, 3] = transpose(data) * (not_data .* weights)
        pxy[:, :, 2] = transpose(not_data) * (data .* weights)
    end
    pxy[:, :, 1] = N .- pxy[:, :, 4] .- pxy[:, :, 3] .- pxy[:, :, 2]

    joint_count = pseudocount / 4
    pxy .+= joint_count

    # diagonal is the marginal of p(x)
    diag = diagind(pxy[:, :, 1])
    @view(pxy[:,:, 1])[diag] .+= joint_count
    @view(pxy[:,:, 2])[diag] .= zero(Float)
    @view(pxy[:,:, 3])[diag] .= zero(Float)
    @view(pxy[:,:, 4])[diag] .+= joint_count

    if data isa CuMatrix
        CUDA.unsafe_free!(data)
        CUDA.unsafe_free!(not_data)
    end
    pxy ./= base
end

function pairwise_MI(data::Union{Matrix{Bool}, CuMatrix{Bool}, BitMatrix};
        weights=nothing,
        pseudocount=0.0,
        Float=Float64)
    pxy = pairwise_marginal(data; weights, pseudocount, Float)
    p0 = diag(pxy[:, :, 1])
    p1 = diag(pxy[:, :, 4])
    D = size(data, 2)
    pxpy = Float.(similar(data, D, D, 4))
    pxpy[:,:,1] = p0 * p0'
    pxpy[:,:,2] = p0 * p1'
    pxpy[:,:,3] = p1 * p0'
    pxpy[:,:,4] = p1 * p1'
    pxy_log_pxy = @. xlogx(pxy)
    pxy_log_pxpy = @. xlogy(pxy, pxpy)
    MI = dropdims(sum(pxy_log_pxy - pxy_log_pxpy, dims=3), dims=3)
    if data isa CuMatrix{Bool}
        CUDA.unsafe_free!(pxy)
        CUDA.unsafe_free!(pxpy)
        CUDA.unsafe_free!(pxy_log_pxy)
        CUDA.unsafe_free!(pxy_log_pxpy)
    end
    MI
end

#############################
# Mutual Information on Categorical Data
#############################

function pairwise_marginal(data::Matrix;
                           weights::Union{Vector, Nothing} = nothing,
                           num_cats = maximum(data) - minimum(data) + 1,
                           pseudocount = zero(Float32),
                           Float = Float32)

    @assert minimum(data) >= 0 "Categorical data labels are assumed to be indexed starting 0"

    num_samples = size(data, 1)
    num_vars = size(data, 2)

    if weights === nothing
        weights = ones(Float32, num_samples)
    else
        pseudocount = pseudocount * sum(weights) / num_samples
    end

    num_vars = size(data, 2)
    pair_margs = Array{Float}(undef, num_vars, num_vars, num_cats, num_cats)
    Z = sum(weights) + pseudocount
    single_smooth = Float(pseudocount / num_cats)
    pair_smooth = Float(single_smooth / num_cats)

    for i = 1:num_vars
        # note: multithreading needs to be on inner loop for thread-safe copying across diagonal
        Threads.@threads for j = 1:num_vars
            if i<=j
                if i!=j
                    @inbounds pair_margs[i,j,:,:] .= pair_smooth
                else
                    @inbounds pair_margs[i,j,:,:] .= zero(eltype(pair_margs))
                    for l = 1:num_cats
                        @inbounds pair_margs[i,j,l,l] = single_smooth
                    end
                end
                @simd for k = 1:size(data,1) # @avx here gives incorrect results
                    @inbounds pair_margs[i,j,data[k,i]+1,data[k,j]+1] += weights[k]
                end
                @inbounds pair_margs[i,j,:,:] ./= Z
            else
                @inbounds pair_margs[i,j,:,:] .= transpose(@view pair_margs[j,i,:,:])
            end
            nothing
        end
    end

    return pair_margs
end

function balance_threads(num_items, num_examples, config; mine, maxe)
    block_threads = config.threads
    # make sure the number of example threads is a multiple of 32
    num_item_batches = cld(num_items, maxe)
    num_blocks = cld(num_item_batches * num_examples, block_threads)
    if num_blocks < config.blocks
        max_num_item_batch = cld(num_items, mine)
        max_num_blocks = cld(max_num_item_batch * num_examples, block_threads)
        num_blocks = min(config.blocks, max_num_blocks)
        num_item_batches = (num_blocks * block_threads) ÷ num_examples
    end
    item_work = cld(num_items, num_item_batches)
    @assert item_work*block_threads*num_blocks >= num_examples*num_items
    block_threads, num_blocks, num_examples, item_work
end

function pairwise_marginal_kernel(pair_margs, data, weights, num_ex_threads::Int32, vpair_work::Int32)
    threadid = ((blockIdx().x - one(Int32)) * blockDim().x) + threadIdx().x

    vpair_batch, ex_id = fldmod1(threadid, num_ex_threads)

    num_vars = size(pair_margs, 1)
    vpair_start = one(Int32) + (vpair_batch - one(Int32)) * vpair_work
    vpair_end = min(vpair_start + vpair_work - one(Int32), num_vars ^ Int32(2))

    @inbounds if ex_id <= size(data, 1)
        for vpair_id = vpair_start : vpair_end
            v1, v2 = fldmod1(vpair_id, num_vars)
            if v1 <= v2
                d1 = data[ex_id, v1]
                d2 = data[ex_id, v2]
                CUDA.@atomic pair_margs[v1,v2,d1+1,d2+1] += weights[ex_id]
            end
        end
    end
    nothing
end

function pairwise_marginal(data::CuMatrix;
                           weights::Union{CuVector, Vector, Nothing} = nothing,
                           num_cats = maximum(data) - minimum(data) + 1,
                           pseudocount = zero(Float32),
                           Float = Float32)

    @assert eltype(data) != Bool
    @assert minimum(data) >= 0 "Categorical data labels are assumed to be indexed starting 0"

    num_examples = size(data, 1)
    num_vars = size(data, 2)
    num_var_pairs = num_vars * num_vars

    if weights === nothing
        weights = ones(Float, num_examples)
    else
        pseudocount = pseudocount * sum(weights) / num_examples
    end

    num_vars = size(data,2)
    pair_margs = zeros(Float, num_vars, num_vars, num_cats, num_cats)
    Z = sum(weights) + pseudocount
    single_smooth = Float(pseudocount / num_cats)
    pair_smooth = Float(single_smooth / num_cats)

    # init pair_margs
    for i = 1 : num_vars
        for j = 1 : num_vars
            if i <= j
                if i != j
                    @inbounds pair_margs[i,j,:,:] .= pair_smooth
                else
                    @inbounds pair_margs[i,j,:,:] .= zero(eltype(pair_margs))
                    for l = 1 : num_cats
                        @inbounds pair_margs[i,j,l,l] = single_smooth
                    end
                end
            end
            nothing
        end
    end

    pair_margs = cu(pair_margs)
    weights = cu(weights)

    dummy_args = (pair_margs, data, weights, Int32(1), Int32(1))
    kernel = @cuda name="pairwise_marginal" launch=false pairwise_marginal_kernel(dummy_args...)
    config = launch_configuration(kernel.fun)

    threads, blocks, num_example_threads, vpair_work =
        balance_threads(num_var_pairs, num_examples, config; mine=2, maxe=32)

    args = (pair_margs, data, weights, Int32(num_example_threads), Int32(vpair_work))
    kernel(args...; threads, blocks)

    pair_margs ./= Z
    pair_margs = Array(pair_margs)

    for i = 1 : num_vars
        for j = 1 : num_vars
            if i > j
                @inbounds pair_margs[i,j,:,:] .= transpose(@view pair_margs[j,i,:,:])
            end
        end
    end

    pair_margs
end

function pairwise_MI(data::Union{Matrix,CuMatrix};
                     weights::Union{Vector, Nothing} = nothing,
                     num_cats = maximum(data) - minimum(data) + 1,
                     pseudocount = zero(Float32),
                     Float = Float32)
    num_samples = size(data, 1)
    num_vars = size(data, 2)

    if weights === nothing
        weights = ones(Int32, num_samples)
    else
        pseudocount = pseudocount * sum(weights) / num_samples
    end

    # `joint_cont[i, j, k, w]' is the total weight of samples whose i- and j-th variable are k and w, respectively
    joint_cont = pairwise_marginal(data; weights, num_cats, pseudocount)

    # `marginal_cont[i, j]' is the total weight of sample whose i-th variable is j
    marginal_cont = zeros(Float, num_vars, num_cats)
    for i = 1:num_vars
        for j = 1:num_cats
            @inbounds marginal_cont[i,j] = joint_cont[i,i,j,j]
        end
    end

    # Compute mutual information
    MI = zeros(Float, num_vars, num_vars)
    for var1_idx = 1 : num_vars
        for var2_idx = var1_idx : num_vars
            @inbounds MI[var1_idx, var2_idx] = sum(
                xlogx.(joint_cont[var1_idx, var2_idx, :, :]) .-
                xlogy.(joint_cont[var1_idx, var2_idx, :, :], marginal_cont[var1_idx, :] .* marginal_cont[var2_idx, :]'))
        end
    end

    for var1_idx = 2 : num_vars
        for var2_idx = 1 : var1_idx - 1
            @inbounds MI[var1_idx, var2_idx] = MI[var2_idx, var1_idx]
        end
    end

    MI
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
(((x, y),) -> x + y).([(1, 2), (3, 4)])

# %%
f((x, y)) = x + y
f.([(1, 2), (3, 4)])

# %%
(t -> +(t...)).([(1, 2), (3, 4)])

# %%

max_stars_repo_path: 0019/Anonymous function of a tuple.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0019/Anonymous function of a tuple.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0019/Anonymous function of a tuple.jl

max_forks_repo_name: genkuroki/public


text:
#-- PRIVATE METHODS BELOW ------------------------------------------------------------------------------- #
function _build_sbml_species_section(problem_dictionary::Dict{String,Any})::Array{String,1}

    # initialize -
    buffer = Array{String,1}()

    # setup list of species -
    # start tag -
    +(buffer, "<listOfSpecies>\n")

    # grab the list of species -
    list_of_species_objects = problem_dictionary["list_of_species"]
    for (index,species) in enumerate(list_of_species_objects)

        # get species information -
        species_symbol = species.species_symbol

        # tmp line -
        tmp_line = "\t<species id=\"$(species_symbol)\" name=\"$(species_symbol)\" compartment=\"default\" boundaryCondition=\"false\"/>\n"

        # push -
        +(buffer, tmp_line)
    end

    # end tag -
    +(buffer, "</listOfSpecies>\n")

    # return -
    return buffer
end

function _build_sbml_reaction_section(problem_dictionary::Dict{String,Any})::Array{String,1}

    # inner functions -
    function _process_reaction_string(reaction_phrase::String)::Array{VLMetabolicSpeciesReference,1}

        # Split around + -
        species_ref_array = Array{VLMetabolicSpeciesReference,1}()
        fragment_array = split(reaction_phrase,"+")
        for fragment in fragment_array

            if (contains(fragment,"*"))
                local_fragment_array = split(fragment,"*");
                coefficient = local_fragment_array[1]
                species_symbol = local_fragment_array[2];

                # build -
                species_ref_object = VLMetabolicSpeciesReference()
                species_ref_object.species_symbol = species_symbol
                species_ref_object.stoichiometry = coefficient

                # cache -
                push!(species_ref_array, species_ref_object)
            else

                # Build -
                species_symbol = fragment;
                coefficient = "1.0"

                # build -
                species_ref_object = VLMetabolicSpeciesReference()
                species_ref_object.species_symbol = species_symbol
                species_ref_object.stoichiometry = coefficient

                # cache -
                push!(species_ref_array, species_ref_object)
            end
        end

        # return -
        return species_ref_array
    end

    # initialize -
    buffer = Array{String,1}()

    # setup list of reactions -
    # start tag -
    +(buffer, "<listOfReactions>\n")

    # grab the array of reactions -
    reaction_order_array = problem_dictionary["metabolic_reaction_order_array"]
    reaction_dictionary = problem_dictionary["metabolic_reaction_dictionary"]
    for (index, reaction_key) in enumerate(reaction_order_array)

        # check do we have this key?
        if (haskey(reaction_dictionary,reaction_key) == true)

            # grab this reaction object -
            reaction_object = reaction_dictionary[reaction_key]

            # get stuff from the reaction object -
            reaction_name = reaction_object.reaction_name
            reversible = reaction_object.reversible
            left_phrase = reaction_object.left_phrase
            right_phrase = reaction_object.right_phrase

            # tmp line -
            reaction_open_line = "\t<reaction id=\"$(reaction_name)\" name=\"$(reaction_name)\" reversible=\"$(reversible)\">\n"

            # process the left phrase -
            left_species_ref_list = _process_reaction_string(left_phrase)
            list_of_reactants_buffer = Array{String,1}()
            +(list_of_reactants_buffer, "\t\t<listOfReactants>\n")
            for (index,specief_ref) in enumerate(left_species_ref_list)

                species_symbol = specief_ref.species_symbol
                stoichiometry = specief_ref.stoichiometry

                # build tmp line -
                tmp_line = "\t\t\t<speciesReference species=\"$(species_symbol)\" stoichiometry=\"$(stoichiometry)\"/>\n"

                # add -
                +(list_of_reactants_buffer,tmp_line)
            end
            +(list_of_reactants_buffer, "\t\t</listOfReactants>\n")

            # process the right phrase -
            right_species_ref_list = _process_reaction_string(right_phrase)
            list_of_products_buffer = Array{String,1}()
            +(list_of_products_buffer, "\t\t<listOfProducts>\n")
            for (index,specief_ref) in enumerate(right_species_ref_list)

                species_symbol = specief_ref.species_symbol
                stoichiometry = specief_ref.stoichiometry

                # build tmp line -
                tmp_line = "\t\t\t<speciesReference species=\"$(species_symbol)\" stoichiometry=\"$(stoichiometry)\"/>\n"

                # add -
                +(list_of_products_buffer,tmp_line)
            end
            +(list_of_products_buffer, "\t\t</listOfProducts>\n")

            # close -
            reaction_close_line = "\t</reaction>\n"

            # grab -
            +(buffer, reaction_open_line)
            +(buffer, list_of_reactants_buffer)
            +(buffer, list_of_products_buffer)
            +(buffer, reaction_close_line)
        else
            error("$(reaction_key) not found in the problem dictionary")
        end
  end

    # end tag -
    +(buffer,"</listOfReactions>\n")

    # return -
    return buffer
end

function _build_sbml_header_section(problem_dictionary::Dict{String,Any})::Array{String,1}

    # initialize -
    buffer = Array{String,1}()

    # add header lines -
    +(buffer,"<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
    +(buffer,"<sbml xmlns=\"http://www.sbml.org/sbml/level2\" level=\"2\" version=\"1\">\n")
    +(buffer,"<model id=\"default_model\" name=\"default_model\">\n")

    # return -
    return buffer
end

function _build_sbml_footer_section(problem_dictionary::Dict{String,Any})::Array{String,1}

    # initialize -
    buffer = Array{String,1}()

    # add footer lines -
    +(buffer, "</model>\n")
    +(buffer, "</sbml>\n")

    # return -
    return buffer
end

function _build_smbl_compartments_section(problem_dictionary::Dict{String,Any})::Array{String,1}

    # initialize -
    buffer = Array{String,1}()

    # <listOfCompartments>
	#     <compartment id="Extra_organism"/>
    # </listOfCompartments>

    # add footer lines -
    +(buffer, "<listOfCompartments>\n")
    +(buffer, "\t<compartment id=\"default\"/>\n")
    +(buffer, "</listOfCompartments>\n")

    # return -
    return buffer
end
# ------------------------------------------------------------------------------------------------------ #

#-- PUBLIC METHODS BELOW ------------------------------------------------------------------------------- #
function build_sbml_model_program_component(problem_dictionary::Dict{String,Any}, filename::String)::VLProgramComponent

    # initalize -
    sbml_buffer = Array{String,1}()
    program_component = VLProgramComponent()

    # build the header buffer -
    header_buffer = _build_sbml_header_section(problem_dictionary)

    # build the list of compartments buffer -
    compartments_buffer = _build_smbl_compartments_section(problem_dictionary)

    # build the species buffer -
    species_buffer = _build_sbml_species_section(problem_dictionary)

    # build the reaction buffer -
    reaction_buffer = _build_sbml_reaction_section(problem_dictionary)

    # build the footer buffer -
    footer_buffer = _build_sbml_footer_section(problem_dictionary)

    # combine the different sections together -
    +(sbml_buffer, header_buffer)

    # compartments -
    +(sbml_buffer, compartments_buffer)

    # add the species -
    +(sbml_buffer,species_buffer)

    # add the reactions -
    +(sbml_buffer, reaction_buffer)

    # add the footer -
    +(sbml_buffer, footer_buffer)

    # collapse -
    flat_buffer = ""
    [flat_buffer *= line for line in sbml_buffer]

    # add data to program_component -
    program_component.filename = filename;
    program_component.buffer = flat_buffer
    program_component.type = :buffer

    # return -
    return program_component
end
# ------------------------------------------------------------------------------------------------------ #

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
# Discretisations based on pointwise samples

"""
    abstract type Sampling{N} <: Discretisation{N}

An interface that only requires D!, `sample(d)` and `interpolant` to be defined.
"""
abstract type Sampling{N} <: Discretisation{N} end

sample(f, d::Sampling) = f.(daxes(d)...)

argand(f, d::Sampling) = f.(argand(d))
# TODO figure out why Complex doesn't work
# argand(d::Discretisation) = sample(Complex, d)
argand(d::Sampling) = sample((x, y) -> x + 1im * y, d)

"""
    interpolate(d, u, r)

Return the value of u interpolated at r
"""
interpolate(d::Sampling, u, r) = dot(interpolant(d, r), u)

"""
    interpolant(d, r)

Return the dual vector that interpolates a discretised function at r.
"""
function interpolant(d::Sampling, r) end

dif!(y, d::Sampling{2}, a, u; axis) = D!(y, d, a, u, 1, axis)
dif2!(y, d::Sampling{2}, a, u; axis) = D!(y, d, a, u, 2, axis)

function daxes(d::Sampling{2})
    x = d.h / 2 * (1-d.n:2:d.n-1)
    x = reshape(x, d.n, 1)
    (x, transpose(x))
end

function primitive_operators(d::Sampling{2})
    # TODO confirm that using scratch storage this way works
    x, y = daxes(d)

    function Δ!(w, a::Number, u)
        for axis = 1:2
            dif2!(w, d, a, u; axis)
        end
    end

    # φ is normalised as a number, so the density is φ/h^N
    U!(w, a, u, r) = (@. w += a * r / d.h^2 * u; w)

    function J!(w, a::Number, u)
        dif!(w, d, 1im * a * y, u, axis = 1)
        dif!(w, d, -1im * a * x, u, axis = 2)
    end

    Δ!, U!, J!
end