text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

@testset "test_coordinates" begin
     img = Raytracing.HDRimage(7, 4)

     @test Raytracing.valid_coordinates(img, 0, 0)
     @test Raytracing.valid_coordinates(img, 3, 2)
     @test Raytracing.valid_coordinates(img, 6, 3)

     @test !Raytracing.valid_coordinates(img, 6, 4)
     @test !Raytracing.valid_coordinates(img, 7, 3)
     @test !Raytracing.valid_coordinates(img, -1, 0)
     @test !Raytracing.valid_coordinates(img, 0, -1)

end

@testset "test_pixel_offset" begin
     rgb_matrix= [ RGB( 3i/255, (3i+1)/255, (3i+2)/255 ) for i in 0:5]
     img = Raytracing.HDRimage(3, 2, rgb_matrix)

     @test_throws AssertionError Raytracing.pixel_offset(img, 3, 0)
     @test_throws AssertionError Raytracing.pixel_offset(img, 2, 2)

     @test img.rgb_m[Raytracing.pixel_offset(img, 0, 0)] ≈ RGB(0/255, 1/255, 2/255)
     @test img.rgb_m[Raytracing.pixel_offset(img, 2, 0)] ≈ RGB(6/255, 7/255, 8/255)
     @test img.rgb_m[Raytracing.pixel_offset(img, 2, 1)] ≈ RGB(15/255, 16/255, 17/255)

end

@testset "test_get_pixel" begin
     rgb_matrix= [ RGB( 3i/255, (3i+1)/255, (3i+2)/255 ) for i in 0:5]
     img = Raytracing.HDRimage(3, 2, rgb_matrix)

     @test_throws AssertionError Raytracing.get_pixel(img, 3, 0)
     @test_throws AssertionError Raytracing.get_pixel(img, 2, 2)

     @test Raytracing.get_pixel(img, 0, 0)  ≈ RGB(0/255, 1/255, 2/255)
     @test Raytracing.get_pixel(img, 2, 0) ≈ RGB(6/255, 7/255, 8/255)
     @test Raytracing.get_pixel(img, 2, 1) ≈ RGB(15/255, 16/255, 17/255)

end

@testset "test_set_pixel" begin
     rgb_matrix= [ RGB( 3i/255, (3i+1)/255, (3i+2)/255 ) for i in 0:5]
     img = Raytracing.HDRimage(3, 2, rgb_matrix)
     C = RGB(0.,0.,0.)

     @test_throws AssertionError Raytracing.set_pixel(img, 3, 0, C)
     @test_throws AssertionError Raytracing.set_pixel(img, 2, 2, C)

     Raytracing.set_pixel(img, 2, 0, C)
     @test Raytracing.get_pixel(img, 2, 0) ≈ RGB(0., 0., 0.)
     @test Raytracing.get_pixel(img, 0, 0)  ≈ RGB(0/255, 1/255, 2/255)
     @test Raytracing.get_pixel(img, 2, 1) ≈ RGB(15/255, 16/255, 17/255)

end

@testset "test_write_pfm" begin
     # Inizializzo e definisco una HDRimage
     img = Raytracing.HDRimage(3, 2)
     Raytracing.set_pixel(img, 0, 0, RGB(1.0e1, 2.0e1, 3.0e1)) # Each component is
     Raytracing.set_pixel(img, 1, 0, RGB(4.0e1, 5.0e1, 6.0e1)) # different from any
     Raytracing.set_pixel(img, 2, 0, RGB(7.0e1, 8.0e1, 9.0e1)) # other: important in
     Raytracing.set_pixel(img, 0, 1, RGB(1.0e2, 2.0e2, 3.0e2)) # tests!
     Raytracing.set_pixel(img, 1, 1, RGB(4.0e2, 5.0e2, 6.0e2))
     Raytracing.set_pixel(img, 2, 1, RGB(7.0e2, 8.0e2, 9.0e2))

     # Leggo file di rifeirmento
     inpf = open("reference_le.pfm", "r") do file	# read( ) legge già di base i bytes grezzi,
               read(file)							# opzioni da poter decidere sono solo "r" e "w" - specificando UInt8 legge solo il primo carattere
          end

     # Variabile stream
     io = IOBuffer(UInt8[], read=true, write=true)	# read, write = true opzionali, lo sono già di default
     Raytracing.write(io, img)

     # Stesso array del file, controllo sulla corretta lettura
     reference_bytes = [
          0x50, 0x46, 0x0a, 0x33, 0x20, 0x32, 0x0a, 0x2d, 0x31, 0x2e, 0x30, 0x0a,
          0x00, 0x00, 0xc8, 0x42, 0x00, 0x00, 0x48, 0x43, 0x00, 0x00, 0x96, 0x43,
          0x00, 0x00, 0xc8, 0x43, 0x00, 0x00, 0xfa, 0x43, 0x00, 0x00, 0x16, 0x44,
          0x00, 0x00, 0x2f, 0x44, 0x00, 0x00, 0x48, 0x44, 0x00, 0x00, 0x61, 0x44,
          0x00, 0x00, 0x20, 0x41, 0x00, 0x00, 0xa0, 0x41, 0x00, 0x00, 0xf0, 0x41,
          0x00, 0x00, 0x20, 0x42, 0x00, 0x00, 0x48, 0x42, 0x00, 0x00, 0x70, 0x42,
          0x00, 0x00, 0x8c, 0x42, 0x00, 0x00, 0xa0, 0x42, 0x00, 0x00, 0xb4, 0x42
     ]

     # Si salva il contenuto di io: se c'è bisogno di accedere al contenuto più di una volta
     # la seconda lettura fallisce poiché operazione stream (il puntatore non torna al principio)
     value = take!(io)		# sbagliato usare read(io): è già inizializzato come IOBuffer()

     # Si controlla che il contenuto di io, quello del file letto e dell'array inizializzato
     # corrispondano tra loro
     @test value == inpf
     @test value == reference_bytes
     @test reference_bytes == inpf

end

@testset "test_read_line" begin
     line = IOBuffer(b"hello\nworld   again \n ...")
     @test Raytracing.read_line(line) == "hello"
     @test Raytracing.read_line(line) == "world   again "
     @test Raytracing.read_line(line) == " ..."
end

@testset "test_parse_img_size" begin
     @test Raytracing.parse_img_size("3 2") == (3, 2)
     @test Raytracing.parse_img_size("1. 2.") == (1,2)
     @test_throws Raytracing.InvalidPfmFileFormat var = Raytracing.parse_img_size("3")
     @test_throws Raytracing.InvalidPfmFileFormat var = Raytracing.parse_img_size("3.14 4")
     @test_throws Raytracing.InvalidPfmFileFormat var = Raytracing.parse_img_size("2 -1")
end

@testset "test_read_float" begin
     # creo matrice come test precedenti, ma con 1 byte in meno per testare errore a fine lettura
     reference_bytes2 = IOBuffer([
          0x50, 0x46, 0x0a, 0x33, 0x20, 0x32, 0x0a, 0x2d, 0x31, 0x2e, 0x30, 0x0a,
          0x00, 0x00, 0xc8, 0x42, 0x00, 0x00, 0x48, 0x43, 0x00, 0x00, 0x96, 0x43,
          0x00, 0x00, 0xc8, 0x43, 0x00, 0x00, 0xfa, 0x43, 0x00, 0x00, 0x16, 0x44,
          0x00, 0x00, 0x2f, 0x44, 0x00, 0x00, 0x48, 0x44, 0x00, 0x00, 0x61, 0x44,
          0x00, 0x00, 0x20, 0x41, 0x00, 0x00, 0xa0, 0x41, 0x00, 0x00, 0xf0, 0x41,
          0x00, 0x00, 0x20, 0x42, 0x00, 0x00, 0x48, 0x42, 0x00, 0x00, 0x70, 0x42,
          0x00, 0x00, 0x8c, 0x42, 0x00, 0x00, 0xa0, 0x42, 0x00, 0x00, 0xb4#, 0x42
     ])

     # leggo l'HEAD
     for i in 0:2
          Raytracing.read_line(reference_bytes2)
     end

     # testo la lettura corretta di tutti i colori (ordine basso → alto nella scrittura)
     @test Raytracing.read_float(reference_bytes2, -1.0) == 1.0e2
     @test Raytracing.read_float(reference_bytes2, -1.0) == 2.0e2
     @test Raytracing.read_float(reference_bytes2, -1.0) == 3.0e2
     @test Raytracing.read_float(reference_bytes2, -1.0) == 4.0e2
     @test Raytracing.read_float(reference_bytes2, -1.0) == 5.0e2
     @test Raytracing.read_float(reference_bytes2, -1.0) == 6.0e2
     @test Raytracing.read_float(reference_bytes2, -1.0) == 7.0e2
     @test Raytracing.read_float(reference_bytes2, -1.0) == 8.0e2
     @test Raytracing.read_float(reference_bytes2, -1.0) == 9.0e2
     @test Raytracing.read_float(reference_bytes2, -1.0) == 1.0e1
     @test Raytracing.read_float(reference_bytes2, -1.0) == 2.0e1
     @test Raytracing.read_float(reference_bytes2, -1.0) == 3.0e1
     @test Raytracing.read_float(reference_bytes2, -1.0) == 4.0e1
     @test Raytracing.read_float(reference_bytes2, -1.0) == 5.0e1
     @test Raytracing.read_float(reference_bytes2, -1.0) == 6.0e1
     @test Raytracing.read_float(reference_bytes2, -1.0) == 7.0e1
     @test Raytracing.read_float(reference_bytes2, -1.0) == 8.0e1
     # errore nella lettura dell'ultimo byte: ne manca 1 per fare un Float32
     @test_throws Raytracing.InvalidPfmFileFormat var = Raytracing.read_float(reference_bytes2, -1.0)

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

     # creo matrice come test precedenti, ma con 1 byte in più per testare errore a fine lettura
     reference_bytes2 = IOBuffer([
          0x50, 0x46, 0x0a, 0x33, 0x20, 0x32, 0x0a, 0x2d, 0x31, 0x2e, 0x30, 0x0a,
          0x00, 0x00, 0xc8, 0x42, 0x00, 0x00, 0x48, 0x43, 0x00, 0x00, 0x96, 0x43,
          0x00, 0x00, 0xc8, 0x43, 0x00, 0x00, 0xfa, 0x43, 0x00, 0x00, 0x16, 0x44,
          0x00, 0x00, 0x2f, 0x44, 0x00, 0x00, 0x48, 0x44, 0x00, 0x00, 0x61, 0x44,
          0x00, 0x00, 0x20, 0x41, 0x00, 0x00, 0xa0, 0x41, 0x00, 0x00, 0xf0, 0x41,
          0x00, 0x00, 0x20, 0x42, 0x00, 0x00, 0x48, 0x42, 0x00, 0x00, 0x70, 0x42,
          0x00, 0x00, 0x8c, 0x42, 0x00, 0x00, 0xa0, 0x42, 0x00, 0x00, 0xb4, 0x42, 0x42
     ])

     # leggo l'HEAD
     for i in 0:2
          Raytracing.read_line(reference_bytes2)
     end
     for i in 1:18
          Raytracing.read_float(reference_bytes2, -1.0)
     end

     # errore nella lettura dell'ultimo byte: ne mancano 3 per fare un Float32
     @test_throws Raytracing.InvalidPfmFileFormat var = Raytracing.read_float(reference_bytes2, -1.0)
end

@testset "test_parse_endianness" begin
     @test Raytracing.parse_endianness("1.0") == 1.0
     @test Raytracing.parse_endianness("+1.0") == 1.0
     @test Raytracing.parse_endianness("-1.0") == -1.0
     @test_throws Raytracing.InvalidPfmFileFormat var = Raytracing.parse_endianness("1.5")
     @test_throws Raytracing.InvalidPfmFileFormat var = Raytracing.parse_endianness("2")
     @test_throws Raytracing.InvalidPfmFileFormat var = Raytracing.parse_endianness("abc")
end

@testset "test_read_pfm" begin
     img_le = open("reference_le.pfm", "r") do file
               Raytracing.read(file, Raytracing.HDRimage)
               end
     img_be = open("reference_be.pfm", "r") do file
               Raytracing.read(file, Raytracing.HDRimage)
               end

     @test img_le.width == 3
     @test img_le.height == 2
     @test Raytracing.get_pixel(img_le, 0, 0) == RGB(1.0e1, 2.0e1, 3.0e1)
     @test Raytracing.get_pixel(img_le, 1, 0) == RGB(4.0e1, 5.0e1, 6.0e1)
     @test Raytracing.get_pixel(img_le, 2, 0) == RGB(7.0e1, 8.0e1, 9.0e1)
     @test Raytracing.get_pixel(img_le, 0, 1) == RGB(1.0e2, 2.0e2, 3.0e2)
     @test Raytracing.get_pixel(img_le, 1, 1) == RGB(4.0e2, 5.0e2, 6.0e2)
     @test Raytracing.get_pixel(img_le, 2, 1) == RGB(7.0e2, 8.0e2, 9.0e2)

     @test img_be.width == 3
     @test img_be.height == 2
     @test Raytracing.get_pixel(img_be, 0, 0) == RGB(1.0e1, 2.0e1, 3.0e1)
     @test Raytracing.get_pixel(img_be, 1, 0) == RGB(4.0e1, 5.0e1, 6.0e1)
     @test Raytracing.get_pixel(img_be, 2, 0) == RGB(7.0e1, 8.0e1, 9.0e1)
     @test Raytracing.get_pixel(img_be, 0, 1) == RGB(1.0e2, 2.0e2, 3.0e2)
     @test Raytracing.get_pixel(img_be, 1, 1) == RGB(4.0e2, 5.0e2, 6.0e2)
     @test Raytracing.get_pixel(img_be, 2, 1) == RGB(7.0e2, 8.0e2, 9.0e2)
end


text:
Base.size(d :: BoundsDeviation) = Base.size(lbounds(d));

get_data_values(d :: BoundsDeviation) = nothing;
lbounds(d :: BoundsDeviation) = d.lbV;
ubounds(d :: BoundsDeviation) = d.ubV;
inside_bounds(d :: BoundsDeviation) = all(d.ubV .>= d.modelV .>= d.lbV);
set_model_values(d :: BoundsDeviation, modelV) = (d.modelV = deepcopy(modelV));

"""
    $(SIGNATURES)

Scalar deviation from one Deviation object.

Optionally includes `scalarWt` factor.
"""
function scalar_dev(d :: BoundsDeviation{F1}; inclScalarWt :: Bool = true) where F1
    if inside_bounds(d)
        scalarDev = zero(F1);
    else
        modelV = get_model_values(d);
        devV = d.wtV .* max.(zero(F1), modelV .- ubounds(d)) .+
            d.wtV .* max.(zero(F1), lbounds(d) .- modelV);
        scalarDev = sum(devV);
        if inclScalarWt
            scalarDev *= d.scalarWt;
        end
    end
    scalarStr = sprintf1(d.fmtStr, scalarDev);

    return scalarDev, scalarStr
end

"""
	$(SIGNATURES)

Show a bounds deviation. As a table.
"""
function bounds_show_fct(d :: BoundsDeviation{F1}; showModel :: Bool = true, fPath :: String = "") where F1
    io = open_show_path(d, fPath = fPath);

    # Dimensions of data matrix
    lbV = lbounds(d);
    ubV = ubounds(d);
    modelV = get_model_values(d);
    nd = ndims(lbV);
    if nd > 2
        @warn "Showing deviation not implemented for ndims = $nd"
    elseif nd == 2
        (nr, nc) = size(lbV);
    elseif nd == 1
        nr = length(lbV);
        nc = 1;
    else
        error("Empty deviation")
    end

    println(io, "$(d.name):  Lb / Model / Ub")
    for ir = 1 : nr
        print(io, "\t $ir: ");
        for ic = 1 : nc
            if showModel
                mStr = sprintf1(d.fmtStr, modelV[ir, ic]);
            else
                mStr = "  --  ";
            end
            lbStr = sprintf1(d.fmtStr, lbV[ir, ic]);
            ubStr = sprintf1(d.fmtStr, ubV[ir, ic]);
            print(io, "\t $lbStr / $mStr / $ubStr");
        end
        print(io, "\n");
    end

    close_show_path(d, io);
    return nothing
end

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

ext: jl


text:
"Your optional docstring here"
function distance(a, b)
    if length(a) != length(b)
        throw((a, b), ArgumentError("Unequal strand lengths!"))
    end
    hamming_distance = 0

    for i in 1:length(a)
        if a[i] != b[i]
            hamming_distance += 1
        end
    end

    return hamming_distance
end

max_stars_repo_name: jfcliang/jl-exercism-solutions

max_issues_repo_name: jfcliang/jl-exercism-solutions

max_forks_repo_name: jfcliang/jl-exercism-solutions


text:
# This file was generated, do not modify it. # hide
best_ensemble = fitted_params(tuned_ensemble).best_model
@show best_ensemble.atom.K
@show best_ensemble.bagging_fraction


text:
# Tests for pca.jl

test_ens_prefix = testfile("AKE_ensemble", "AKE.pdb")

@testset "PCA" begin
    # Values tested against Bio3D values in R
    ensemble = readensemble(test_ens_prefix, 10) # These structures are already aligned
    pca = PCA(ensemble)
    @test isapprox(pca.evals[1] / pca.evals[2], 5.26960, rtol=0.001)
    @test isapprox(abs(pca.evecs[1, 1]), 0.01500, rtol=0.001)
    @test isapprox(abs(pca.evecs[642, 1]), 0.00201, rtol=0.001)
    @test isapprox(abs(pca.evecs[1, 2]), 0.00381, rtol=0.001)

    ensemble = readensemble(test_ens_prefix, 10) # These structures are already aligned
    pca = PCA(ensemble)
    pcs = projectensemble(ensemble, pca)
    @test isapprox(abs(pcs[1, 1]), 20.32696, rtol=0.01)
    @test isapprox(abs(pcs[9, 1]), 3.73002, rtol=0.01)
    @test isapprox(abs(pcs[1, 10]), 41.06011, rtol=0.01)

    ensemble = readensemble(test_ens_prefix, 10) # These structures are already aligned
    pca = PCA(ensemble)
    pcs = projectstructure(ensemble.atoms, pca)
    @test isapprox(abs(pcs[1]), 20.32696, rtol=0.01)
    @test isapprox(abs(pcs[9]), 3.73002, rtol=0.01)

    atom_one = Atom("CA", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("CA", "ALA", 'A', 21, [0.0, 0.0, 0.0], "C")
    atom_thr = Atom("CA", "ALA", 'A', 22, [0.0, 0.0, 0.0], "C")
    atoms = [atom_one, atom_two, atom_thr]
    coords_one = [
        0.0 -2.0 0.0;
        0.0 0.0 0.0;
        0.0 0.0 0.0;
    ]
    coords_two = [
        1.0 0.0 3.0;
        3.0 0.0 sqrt(3);
        0.0 0.0 0.0;
    ]
    coords_thr = [
        2.0 2.0 3.0;
        0.0 0.0 -sqrt(3);
        0.0 0.0 0.0;
    ]
    struc_one = ModelledStructure(0.0, coords_one)
    struc_two = ModelledStructure(0.0, coords_two)
    struc_thr = ModelledStructure(0.0, coords_thr)
    strucs = [struc_one, struc_two, struc_thr]
    ensemble = ModelledEnsemble(atoms, strucs)
    flucs_test = fluctuations(ensemble)
    flucs_real = [8/3, 8/3, 4.0]
    @test isapprox(flucs_test, flucs_real)
end


text:
module mqpit1
using Elfel
using Elfel.RefShapes: npts, param_coords, weights
using Elfel.FElements: FE, FEData, refshape
using Elfel.FElements: bfun, bfungradpar, FEH1_T3
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns
using MeshSteward: Mesh, load, baseincrel
using Elfel.QPIterators: QPIterator, bfun, bfungradpar, weight
using Test
# using BenchmarkTools
function test()
    fe = FEH1_T3()
    mesh = load(Mesh(), "qmesh.mesh")
    fesp = FESpace(Float64, mesh, fe)

    qpit = QPIterator(fesp, (kind = :default,))
    @test length(qpit) == 1
    # @test isapprox(param_coords(quadrule(idom)), [0.3333333333333333 0.3333333333333333])
    # @test isapprox(weights(quadrule(idom)), [0.5])

    refNs = [[0.3333333333333334, 0.3333333333333333, 0.3333333333333333]]
    k = 1
    for qp in qpit
        Ns = bfun(qp)
        @test isapprox(Ns, refNs[k])
        k = k + 1
    end

    qpit = QPIterator(fesp, (kind = :default, npts = 3))
    @test length(qpit) == 3
    # @test isapprox(param_coords(quadrule(idom)), [0.3333333333333333 0.3333333333333333])
    # @test isapprox(weights(quadrule(idom)), [0.5])

    refNs = [[0.1666666666666667, 0.6666666666666666, 0.16666666666666666],
    [0.16666666666666674, 0.16666666666666666, 0.6666666666666666],
    [0.6666666666666667, 0.16666666666666666, 0.16666666666666666]]
    k = 1
    for qp in qpit
        Ns = bfun(qp)
        @test isapprox(Ns, refNs[k])
        k = k + 1
    end
    k = 1
    for qp in qpit
        Ns = bfun(qp)
        @test isapprox(Ns, refNs[k])
        k = k + 1
    end

    refgNs = [[-1.0 -1.0], [1.0 0.0], [0.0 1.0]]
    for qp in qpit
        gradNs = bfungradpar(qp)
        for i in 1:length(gradNs)
            @test isapprox(gradNs[i], refgNs[i])
        end
    end

    for qp in qpit
        w = weight(qp)
        @test w == 0.16666666666666666
    end
    # Ns, gradNparams = bfundata(idom)
    # @test length(Ns[1]) == 3

    # geom = geomattr(femesh)
    # conn = connectivity(femesh)

    # it = FEIterator(fesp)
    # @show c = iterate(it)
    # Ns, gradNparams = bfundata(idom)
    # J = jac(geom, c, gradNparams[1])
    # # @btime J = jac($geom, retrieve($ir, 1), $gradNparams[1])
    # @test isapprox(J, [0.5 0.5; 0.0 0.3333333333333333])
end
end
using .mqpit1
mqpit1.test()


text:
# This file is part of Kpax3. License is MIT.

function test_mcmc_merge_init()
  # merge cluster 2 and cluster 3
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, maxclust=2, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  R = [13; 13; 13; 42; 42; 76]
  k = length(unique(R)) - 1

  ij = [4; 6]
  S = 1
  u = 5

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  Kpax3.initsupportmerge!(ij, k, data, priorC, support)

  wi = zeros(Float64, 4, m)
  for col in 1:m
    wi[1, col] = priorC.logγ[1] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[1, col], priorC.B[1, col])
    wi[2, col] = priorC.logγ[2] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[2, col], priorC.B[2, col])
    wi[3, col] = priorC.logγ[3] + priorC.logω[k][1] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[3, col], priorC.B[3, col])
    wi[4, col] = priorC.logγ[3] + priorC.logω[k][2] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[4, col], priorC.B[4, col])
  end

  ci = Float64[log(sum(exp.(support.wi.w[:, b]))) for b in 1:m]

  wj = zeros(Float64, 4, m)
  for col in 1:m
    wj[1, col] = priorC.logγ[1] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[1, col], priorC.B[1, col])
    wj[2, col] = priorC.logγ[2] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[2, col], priorC.B[2, col])
    wj[3, col] = priorC.logγ[3] + priorC.logω[k][1] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[3, col], priorC.B[3, col])
    wj[4, col] = priorC.logγ[3] + priorC.logω[k][2] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[4, col], priorC.B[4, col])
  end

  cj = Float64[log(sum(exp.(support.wj.w[:, b]))) for b in 1:m]

  @test support.vi == 1
  @test support.ni == float(data[:, ij[1]])
  @test support.ui == [ij[1]; 0; 0; 0; 0; 0]

  @test maximum(abs.(support.wi.w - wi)) <= ε
  @test maximum(abs.(support.wi.c - ci)) <= ε
  @test support.wi.z == zeros(Float64, 4, m)

  @test support.vj == 1
  @test support.nj == float(data[:, ij[2]])
  @test support.uj == [ij[2]; 0; 0; 0; 0; 0]

  @test maximum(abs.(support.wj.w - wj)) <= ε
  @test maximum(abs.(support.wj.c - cj)) <= ε
  @test support.wj.z == zeros(Float64, 4, m)

  nothing
end

test_mcmc_merge_init()

function test_mcmc_merge_updatei()
  # move the first unit (u) to cluster 2 (test inverse split operator)
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, maxclust=2, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  R = [13; 13; 13; 42; 42; 76]
  k = length(unique(R)) - 1

  ij = [4; 6]
  S = 1
  u = 5

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  Kpax3.initsupportmerge!(ij, k, data, priorC, support)

  lcp = zeros(Float64, 2)
  for b in 1:m
    lcp[1] += Kpax3.computeclusteriseqprobs!(data[b, u], b, priorC, support)
    lcp[2] += Kpax3.computeclusterjseqprobs!(data[b, u], b, priorC, support)
  end
  Kpax3.updateclusteri!(u, data, support)

  wi = zeros(Float64, 4, m)
  for col in 1:m
    wi[1, col] = priorC.logγ[1] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[1, col], priorC.B[1, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[1, col], priorC.B[1, col])
    wi[2, col] = priorC.logγ[2] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[2, col], priorC.B[2, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[2, col], priorC.B[2, col])
    wi[3, col] = priorC.logγ[3] + priorC.logω[k][1] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[3, col], priorC.B[3, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[3, col], priorC.B[3, col])
    wi[4, col] = priorC.logγ[3] + priorC.logω[k][2] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[4, col], priorC.B[4, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[4, col], priorC.B[4, col])
  end

  zi = copy(wi)

  ci = Float64[log(sum(exp.(support.wi.w[:, b]))) for b in 1:m]

  wj = zeros(Float64, 4, m)
  zj = zeros(Float64, 4, m)
  for col in 1:m
    wj[1, col] = priorC.logγ[1] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[1, col], priorC.B[1, col])
    wj[2, col] = priorC.logγ[2] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[2, col], priorC.B[2, col])
    wj[3, col] = priorC.logγ[3] + priorC.logω[k][1] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[3, col], priorC.B[3, col])
    wj[4, col] = priorC.logγ[3] + priorC.logω[k][2] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[4, col], priorC.B[4, col])

    zj[1, col] = wj[1, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[1, col], priorC.B[1, col])
    zj[2, col] = wj[2, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[2, col], priorC.B[2, col])
    zj[3, col] = wj[3, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[3, col], priorC.B[3, col])
    zj[4, col] = wj[4, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[4, col], priorC.B[4, col])

  end

  cj = Float64[log(sum(exp.(support.wj.w[:, b]))) for b in 1:m]

  @test support.vi == 2
  @test support.ni == float(data[:, ij[1]]) + float(data[:, u])
  @test support.ui == [ij[1]; u; 0; 0; 0; 0]

  @test maximum(abs.(support.wi.w - wi)) <= ε
  @test maximum(abs.(support.wi.c - ci)) <= ε
  @test maximum(abs.(support.wi.z - zi)) <= ε

  @test support.vj == 1
  @test support.nj == float(data[:, ij[2]])
  @test support.uj == [ij[2]; 0; 0; 0; 0; 0]

  @test maximum(abs.(support.wj.w - wj)) <= ε
  @test maximum(abs.(support.wj.c - cj)) <= ε
  @test maximum(abs.(support.wj.z - zj)) <= ε

  nothing
end

test_mcmc_merge_updatei()

function test_mcmc_merge_updatej()
  # move the first unit (u) to cluster 2 (test inverse split operator)
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, maxclust=2, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  R = [13; 13; 13; 42; 42; 76]
  k = length(unique(R)) - 1

  ij = [4; 6]
  S = 1
  u = 5

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  Kpax3.initsupportmerge!(ij, k, data, priorC, support)

  lcp = zeros(Float64, 2)
  for b in 1:m
    lcp[1] += Kpax3.computeclusteriseqprobs!(data[b, u], b, priorC, support)
    lcp[2] += Kpax3.computeclusterjseqprobs!(data[b, u], b, priorC, support)
  end
  Kpax3.updateclusterj!(u, data, support)

  wi = zeros(Float64, 4, m)
  zi = zeros(Float64, 4, m)
  for col in 1:m
    wi[1, col] = priorC.logγ[1] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[1, col], priorC.B[1, col])
    wi[2, col] = priorC.logγ[2] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[2, col], priorC.B[2, col])
    wi[3, col] = priorC.logγ[3] + priorC.logω[k][1] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[3, col], priorC.B[3, col])
    wi[4, col] = priorC.logγ[3] + priorC.logω[k][2] + Kpax3.logmarglik(data[col, ij[1]], 1, priorC.A[4, col], priorC.B[4, col])

    zi[1, col] = wi[1, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[1, col], priorC.B[1, col])
    zi[2, col] = wi[2, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[2, col], priorC.B[2, col])
    zi[3, col] = wi[3, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[3, col], priorC.B[3, col])
    zi[4, col] = wi[4, col] + Kpax3.logcondmarglik(data[col, u], data[col, ij[1]], 1, priorC.A[4, col], priorC.B[4, col])
  end

  ci = Float64[log(sum(exp.(support.wi.w[:, b]))) for b in 1:m]

  wj = zeros(Float64, 4, m)
  for col in 1:m
    wj[1, col] = priorC.logγ[1] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[1, col], priorC.B[1, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[1, col], priorC.B[1, col])
    wj[2, col] = priorC.logγ[2] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[2, col], priorC.B[2, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[2, col], priorC.B[2, col])
    wj[3, col] = priorC.logγ[3] + priorC.logω[k][1] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[3, col], priorC.B[3, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[3, col], priorC.B[3, col])
    wj[4, col] = priorC.logγ[3] + priorC.logω[k][2] + Kpax3.logmarglik(data[col, ij[2]], 1, priorC.A[4, col], priorC.B[4, col]) + Kpax3.logcondmarglik(data[col, u], data[col, ij[2]], 1, priorC.A[4, col], priorC.B[4, col])
  end
  zj = copy(wj)

  cj = Float64[log(sum(exp.(support.wj.w[:, b]))) for b in 1:m]

  @test support.vi == 1
  @test support.ni == float(data[:, ij[1]])
  @test support.ui == [ij[1]; 0; 0; 0; 0; 0]

  @test maximum(abs.(support.wi.w - wi)) <= ε
  @test maximum(abs.(support.wi.c - ci)) <= ε
  @test maximum(abs.(support.wi.z - zi)) <= ε

  @test support.vj == 2
  @test support.nj == float(data[:, ij[2]]) + float(data[:, u])
  @test support.uj == [ij[2]; u; 0; 0; 0; 0]

  @test maximum(abs.(support.wj.w - wj)) <= ε
  @test maximum(abs.(support.wj.c - cj)) <= ε
  @test maximum(abs.(support.wj.z - zj)) <= ε

  nothing
end

test_mcmc_merge_updatej()


text:
# Common functions compatible for use with physical vectors

function dot_product(A::AbstractArray{T1,N},B::AbstractArray{T2,N}) where {T1<:Number,T2<:Number,N}

    size(A) === size(B) || throw(ArgumentError("Vectors must have the same shape"))
    u = zeros(promote_type(T1,T2),size(A)[1:end-1]...)
    dot_product!(u,A,B)
    return u
end

function dot_product!(u::AbstractArray{T1,N}, A::AbstractArray{T2,N2},B::AbstractArray{T3,N2}) where {T1<:Number,T2<:Number,T3<:Number,N,N2}

    (size(A) === size(B)) || throw(ArgumentError("Vectors must have the same shape"))
    (N === N2-1) || throw(ArgumentError("Output should be scalar matrix, one dimension less than input"))

    T = promote_type(T1,T2,T3)
    @turbo for I in CartesianIndices(u)
        x_temp = zero(T)
        for i in 1:N
            x_temp += A[I,i]*B[I,i]
        end
        u[I] = x_temp
    end
    return u
end

function cross_product(A::AbstractArray{T1,4},B::AbstractArray{T2,4}) where {T1<:Number,T2<:Number}

    size(A) === size(B) || throw(ArgumentError("Vectors must have the same shape"))
    s = size(A)
    u = zeros(promote_type(T1,T2),s...)
    cross_product!(u,A,B)
    return u
end

function cross_product!(u::AbstractArray{T1,4},A::AbstractArray{T2,4},B::AbstractArray{T3,4}) where {T1<:Number,T2<:Number,T3<:Number}

    (size(A) === size(B) && size(A) === size(u))|| throw(ArgumentError("Vectors must have the same shape"))
    s = size(u)

    @turbo for r in 1:s[3], q in 1:s[2], p in 1:s[1]
        u[p,q,r,1] = A[p,q,r,2]*B[p,q,r,3] - A[p,q,r,3]*B[p,q,r,2]
        u[p,q,r,2] = A[p,q,r,3]*B[p,q,r,1] - A[p,q,r,1]*B[p,q,r,3]
        u[p,q,r,3] = A[p,q,r,1]*B[p,q,r,2] - A[p,q,r,2]*B[p,q,r,1]
    end
end

function square_norm(A::AbstractArray{T,N}) where {T<:Number,N}
    u = zeros(T,size(A)[1:end-1]...)
    square_norm!(u,A)
    return u
end

function square_norm!(u::AbstractArray{T1,N1},A::AbstractArray{T2,N2}) where {T1<:Number,T2<:Number,N1,N2}
    dot_product!(u,A,A)
    @. u = u^0.5
end


text:
function update!(configuration::Configuration)
    for key in keys(configuration.parameters)
        configuration.value[key] = configuration[key].value
    end
    configuration
end

function update!{T <: Parameter}(configuration::Configuration, parameters::Array{T})
    for param in parameters
        configuration[param.name] = deepcopy(param)
    end
    update!(configuration)
    configuration
end

function update!{T <: Parameter}(configuration::Configuration, parameters::Dict{String, T})
    key_set = collect(keys(parameters))
    for key in key_set
        configuration[key] = deepcopy(parameters[key])
    end
    update!(configuration)
    configuration
end


text:
using XYZ
using StaticArrays
using BenchmarkTools, Compat

a = collect(1:16)
a = reshape(a, (2,2,4))
XYZ.offsets(a)

bbox = XYZ.NBoundingBox(xmin=542900.0, xmax=544100.0, ymin=9679900.0, ymax=9681077.6, zmin=1.06, zmax=329.09)
dims = SVector(5, 5, 5)
coords = MVector(1.0,2.0,3.0)
coords = MVector(1,-1,2)
@btime invalid_coords(dims, coords)
@btime sub2ind(dims, coords[1], coords[2], coords[3])


text:
using Documenter, MagneticResonanceSignals

makedocs(;
    modules=[MagneticResonanceSignals],
    format=Documenter.HTML(),
    pages=[
        "Overview" => "index.md",
        "How-to Guide" => "guide.md",
        "Reference" => "reference.md",
    ],
    repo="https://github.com/TRIImaging/MagneticResonanceSignals.jl/blob/{commit}{path}#L{line}",
    sitename="MagneticResonanceSignals.jl",
    authors="Chris Foster <chris42f@gmail.com>"
)

deploydocs(;
    repo="github.com/TRIImaging/MagneticResonanceSignals.jl",
    push_preview=true
)

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
using JsonGrinder
using Test
using Documenter

function buf_printtree(data; kwargs...)
    buf = IOBuffer()
	printtree(buf, data; kwargs...)
	String(take!(buf))
end

const ≃ = isequal

@testset "JsonGrinder.jl" begin

	# this must be first, because other tests include more imports which break string asserts on types
	@testset "Doctests" begin
	    # there are breaking changes in printing of things between julia 1.5 and 1.6, so we run tests only on 1.6
	    @static if VERSION >= v"1.6.0-"
        	DocMeta.setdocmeta!(JsonGrinder, :DocTestSetup, :(using JsonGrinder); recursive=true)
        	doctest(JsonGrinder)
    	end
	end

	@testset "Extractors" begin
		include("extractors.jl")
		include("multirepresentation.jl")
	end
	@testset "Pipeline" begin
		include("pipelinetest.jl")
	end
	@testset "Schema" begin
		include("schema.jl")
	end
	@testset "HierarchicalUtils" begin
		include("hierarchical_utils_extractors.jl")
		include("hierarchical_utils_schema.jl")
	end
	include("show_html.jl")
	include("util.jl")
	include("json3.jl")

end


text:
# Julia benchmarks from R's data.table
# https://github.com/Rdatatable/data.table/wiki/Benchmarks-:-Grouping

using Random
using DataFrames, DataFramesMeta
using CategoricalArrays

N=10_0000; K=100
Random.seed!(1)

# Array version

DA = DataFrame(
  id1 = rand([Symbol("id", i) for i=1:K], N),          # large groups (char)
  id2 = rand([Symbol("id", i) for i=1:K], N),          # large groups (char)
  id3 = rand([Symbol("id", i) for i=1:N÷K], N),        # small groups (char)
  id4 = rand(1:K, N),                          # large groups (int)
  id5 = rand(1:K, N),                          # large groups (int)
  id6 = rand(1:N÷K, N),                        # small groups (int)
  v1 =  rand(1:5, N),                          # int in range [1,5]
  v2 =  rand(1:5, N),                          # int in range [1,5]
  v3 =  rand(N)                                # numeric e.g. 23.5749
);

# CategoricalArray version

DCA = DataFrame(
  id1 = CategoricalArray(rand([Symbol("id", i) for i=1:K], N)),     # large groups (char)
  id2 = CategoricalArray(rand([Symbol("id", i) for i=1:K], N)),     # large groups (char)
  id3 = CategoricalArray(rand([Symbol("id", i) for i=1:N÷K], N)),   # small groups (char)
  id4 = CategoricalArray(rand(1:K, N)),                          # large groups (int)
  id5 = CategoricalArray(rand(1:K, N)),                          # large groups (int)
  id6 = CategoricalArray(rand(1:N÷K, N)),                        # small groups (int)
  v1 =  rand(1:5, N),                          # int in range [1,5]
  v2 =  rand(1:5, N),                          # int in range [1,5]
  v3 =  rand(N)                                # numeric e.g. 23.5749
);

# Array{Union{T, Missing}} version

DMA = DataFrame(
  id1 = Array{Union{Symbol, Missing}}(rand([Symbol("id", i) for i=1:K], N)),          # large groups (char)
  id2 = Array{Union{Symbol, Missing}}(rand([Symbol("id", i) for i=1:K], N)),          # large groups (char)
  id3 = Array{Union{Symbol, Missing}}(rand([Symbol("id", i) for i=1:N÷K], N)),        # small groups (char)
  id4 = Array{Union{Int, Missing}}(rand(1:K, N)),                          # large groups (int)
  id5 = Array{Union{Int, Missing}}(rand(1:K, N)),                          # large groups (int)
  id6 = Array{Union{Int, Missing}}(rand(1:N÷K, N)),                        # small groups (int)
  v1 =  Array{Union{Int, Missing}}(rand(1:5, N)),                          # int in range [1,5]
  v2 =  Array{Union{Int, Missing}}(rand(1:5, N)),                          # int in range [1,5]
  v3 =  Array{Union{Float64, Missing}}(rand(N))                            # numeric e.g. 23.5749
);

function dt_timings(D)
    @time @by(D, :id1, :sv =sum(:v1));
    @time @by(D, :id1, :sv =sum(:v1));
    @time @by(D, [:id1, :id2], :sv =sum(:v1));
    @time @by(D, [:id1, :id2], :sv =sum(:v1));
    @time @by(D, :id3, :sv = sum(:v1), :mv3 = mean(:v3));
    @time @by(D, :id3, :sv = sum(:v1), :mv3 = mean(:v3));
    @time aggregate(D[!, [4;7:9]], :id4, mean);
    @time aggregate(D[!, [4;7:9]], :id4, mean);
    @time aggregate(D[!, [6;7:9]], :id6, sum);
    @time aggregate(D[!, [6;7:9]], :id6, sum);
    return
end

dt_timings(DA)
dt_timings(DCA)
dt_timings(DMA)

@profile @by(DA, :id1, :sv =sum(:v1));


text:
using Downloads

bathymetry_path = joinpath(@__DIR__, "bathymetry-360x150-latitude-75.0.jld2")
boundary_conditions_path = joinpath(@__DIR__, "boundary_conditions-1degree.jld2")
initial_condition_path = joinpath(@__DIR__, "initial_conditions-1degree.jld2")

# TODO: convert to DataDeps

download_bathymetry(path=bathymetry_path) =
    Downloads.download("https://www.dropbox.com/s/axyzt88g0nr9dbc/bathymetry-360x150-latitude-75.0.jld2", path)

download_boundary_conditions(path=boundary_conditions_path) =
    Downloads.download("https://www.dropbox.com/s/7sbrq5fcnuhtvqp/boundary_conditions-1degree.jld2", path)

download_initial_condition(path=initial_condition_path) =
    Downloads.download("https://www.dropbox.com/s/hbbzeg301hwm2xd/initial_conditions-1degree.jld2", path)


text:
using Documenter, MCMCBenchmarks

function add_image(folder,file)
    mkpath("docs/build/benchmarks/images/"*folder)
    cp("Examples/"*folder*"/results/"*file,
       "docs/build/benchmarks/images/"*folder*"/"*file)
end

makedocs(
    modules = [MCMCBenchmarks],
    checkdocs = :exports,
    format = :html,
    clean = true,
    authors = "Christopher R. Fisher, Rob J Goedman",
    sitename = "StatisticalRethinkingJulia/MCMCBenchmarks.jl",
    pages = Any[
        "Home"=>"index.md",
        "Purpose"=>"purpose.md",
        "Design"=>"design.md",
        "Functions"=>"functions.md",
        "Example"=>"example.md",
        "Benchmark Results"=>"benchmarks.md"
        ]
)

folders = ["Gaussian","Gaussian","Gaussian","SDT","SDT","SDT",
    "Linear_Regression","Linear_Regression","Linear_Regression",
    "LBA","LBA","LBA","Hierarchical_Poisson","Hierarchical_Poisson",
    "Hierarchical_Poisson","Autodiff","Autodiff","Autodiff"]

files = ["summary_time.png","summary_allocations.png","density_mu_ess.png",
    "summary_time.png","summary_allocations.png","density_d_ess.png",
    "summary_time.png","summary_allocations.png","summary_B0_ess.png",
    "summary_time.png","summary_allocations.png","density_A_ess.png",
    "summary_time.png","summary_allocations.png","density_a0_ess.png",
    "summary_time.png","summary_allocations.png","density_a0_ess.png"]

[add_image(folder,file) for (folder,file) in zip(folders,files)]

deploydocs(
    repo = "github.com/StatisticalRethinkingJulia/MCMCBenchmarks.jl.git",
)


text:
function extract_frames(video, frame_folder)
    path = joinpath(frame_folder, "frames%04d.png")
    FFMPEG.ffmpeg_exe(`-loglevel quiet -i $video -y $path`)
end

function compare_media(a, b; sigma=[1,1], eps=0.02)
    file, ext = splitext(a)
    if ext in (".png", ".jpg", ".jpeg", ".JPEG", ".JPG")
        imga = load(a)
        imgb = load(b)
        if size(imga) != size(imgb)
            @warn "images don't have the same size, difference will be Inf"
            return Inf
        end
        return Images.test_approx_eq_sigma_eps(imga, imgb, sigma, Inf)
    elseif ext in (".mp4", ".gif")
        mktempdir() do folder
            afolder = joinpath(folder, "a")
            bfolder = joinpath(folder, "b")
            mkpath(afolder); mkpath(bfolder)
            extract_frames(a, afolder)
            extract_frames(b, bfolder)
            aframes = joinpath.(afolder, readdir(afolder))
            bframes = joinpath.(bfolder, readdir(bfolder))
            if length(aframes) > 10
                # we don't want to compare too many frames since it's time costly
                # so we just compare 10 random frames if more than 10
                samples = rand(1:length(aframes), 10)
                aframes = aframes[samples]
                bframes = bframes[samples]
            end
            # test by maximum diff
            return mean(compare_media.(aframes, bframes; sigma=sigma, eps=eps))
        end
    else
        error("Unknown media extension: $ext")
    end
end

function compare(test_files::Vector{String}, reference_dir::String; missing_refimages=String[], scores=Dict{String,Float64}())
    for test_path in test_files
        ref_path = joinpath(reference_dir, basename(test_path))
        if isdir(test_path)
            if !isdir(ref_path)
                push!(missing_refimages, test_path)
            else
                compare(joinpath.(test_path, readdir(test_path)), ref_path; missing_refimages=missing_refimages, scores=scores)
            end
        elseif isfile(test_path)
            if !isfile(ref_path)
                push!(missing_refimages, test_path)
            else
                diff = compare_media(test_path, ref_path)
                name = replace(ref_path, reference_dir => "")[2:end]
                @info(@sprintf "%1.4f == %s\n" diff name)
                scores[name] = diff
            end
        end
    end
    return missing_refimages, scores
end

# missing_imgs, scores = compare(recording_dir, joinpath(@__DIR__, "refimages"))

function reference_tests(recorded; ref_images = ReferenceTests.download_refimages(), difference=0.03)
    @testset "Reference Image Tests" begin
        missing_files, scores = ReferenceTests.compare(joinpath.(recorded, readdir(recorded)), ref_images)
        @testset "$name" for (name, score) in scores
            @test score < difference
        end
    end
end

function record_tests(db=load_database(); recording_dir=joinpath(@__DIR__, "..", "recorded"))
    recorded_files = String[]
    @testset "Record tests" begin
        rm(recording_dir, recursive=true, force=true)
        mkdir(recording_dir)
        AbstractPlotting.inline!(true)
        no_backend = AbstractPlotting.current_backend[] === missing
        for (source_location, entry) in db
            AbstractPlotting.set_theme!(resolution=(500, 500))
            # we currently can't record anything without a backend!
            if no_backend && ((:Record in entry.used_functions) || (:Stepper in entry.used_functions))
                continue
            end
            RNG.seed_rng!()
            result = Base.invokelatest(entry.func)
            # only save if we have a backend for saving
            uname = unique_name(entry)
            if !no_backend
                save_result(joinpath(recording_dir, uname), result)
            end
            push!(recorded_files, uname)
            @info("Tested: $(nice_title(entry))")
            @test true
        end
    end
    return recorded_files, recording_dir
end

function run_tests(db=load_database(); ref_images = ReferenceTests.download_refimages(),
                    recording_dir=joinpath(@__DIR__, "..", "recorded"), difference=0.03)
    files, dir = record_tests(db, recording_dir=recording_dir)
    reference_tests(recording_dir, ref_images=ref_images, difference=difference)
end


text:
@testset "H-W test" begin
    # <https://en.wikipedia.org/wiki/Hardy–Weinberg_principle#Example_.7F.27.22.60UNIQ--postMath-0000001E-QINU.60.22.27.7F_test_for_deviation>
    #@code_warntype hwe(1469, 138, 5)
    @test VCFTools.hwe(1469, 138, 5) ≈ 0.406002899
end

@testset "openvcf" begin
    # download and extract test file if not exist
    isfile("test.08Jun17.d8b.vcf.gz") || download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
        joinpath(dirname(pathof(VCFTools)), "..", "test/test.08Jun17.d8b.vcf.gz"))
    @test_throws ArgumentError openvcf("test.08Jun17.d8b.vcf.gz", "r+")
    @test_throws ArgumentError openvcf("test.08Jun17.d8b.vcf.gz", "w+")
    @test_throws ArgumentError openvcf("test.08Jun17.d8b.vcf.gz", "c")
    # VCF file name has to end with vcf or vcf.gz
    @test_throws ArgumentError openvcf("test.08Jun17.d8b")
end

@testset "nrecords(vcf)" begin
    # download and extract test file if not exist
    isfile("test.08Jun17.d8b.vcf.gz") ||   download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
        joinpath(dirname(pathof(VCFTools)), "..", "test/test.08Jun17.d8b.vcf.gz"))
    @test nrecords("test.08Jun17.d8b.vcf.gz") == 1356
end

@testset "nsamples(vcf)" begin
    # download test file if not exist
    isfile("test.08Jun17.d8b.vcf.gz") ||  download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
        joinpath(dirname(pathof(VCFTools)), "..", "test/test.08Jun17.d8b.vcf.gz"))
    @test nsamples("test.08Jun17.d8b.vcf.gz") == 191
end

@testset "gtstats(record)" begin
    record = VCF.Record("20\t14370\trs6054257\tG\tA\t29\tPASS\tNS=3;DP=14;AF=0.5;DB;H2\tGT:GQ:DP:HQ\t0|0:48:1:51,51\t1|0:48:8:51,51")
    #@code_warntype gtstats(record)
    @test gtstats(record)[1:end-1] == (1, 1, 0, 3, 1, 0.75, 0.25, 0, true, 0.25)
    @test gtstats(record)[end] ≈ 0.82414089574
    # with missing data
    record = VCF.Record("20\t14370\trs6054257\tG\tA\t29\tPASS\tNS=3;DP=14;AF=0.5;DB;H2\tGT:GQ:DP:HQ\t.\t1|0:48:8:51,51")
    @test gtstats(record)[1:end-1] == (0, 1, 0, 1, 1, 0.5, 0.5, 1, false, 0.5)
    @test gtstats(record)[end] ≈ 0.31731050786291415
end

@testset "gtstats(vcf)" begin
    #@code_warntype gtstats("test.08Jun17.d8b.vcf.gz", STDOUT)
    #@time records, samples, lines = gtstats("test.08Jun17.d8b.vcf.gz")
    #@time gtstats("test.08Jun17.d8b.vcf")
    # output tuple: (records, samples, lines, missings_by_sample, missings_by_record, maf_by_record, minorallele_by_record)

    # download and extract test file and unzip
    isfile("test.08Jun17.d8b.vcf.gz") || download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
        joinpath(dirname(pathof(VCFTools)), "..", "test/test.08Jun17.d8b.vcf.gz"))
    write("test.08Jun17.d8b.vcf", read(GzipDecompressorStream(open("test.08Jun17.d8b.vcf.gz", "r"))))

    @testset "input: text file, output: text file" begin
        @time out = gtstats("test.08Jun17.d8b.vcf", "gtstats.out.txt")
        @test out[1:3] == (1356, 191, 1356)
        outtxt = readdlm("gtstats.out.txt")
        @test all(outtxt[:, 9]  .== out[5]) # missings_by_record
        @test all(outtxt[:, 14] .== out[6]) # maf_by_record
    end

    @testset "input: text file, output: gz file" begin
        @time out = gtstats("test.08Jun17.d8b.vcf", "gtstats.out.gz")
        @test out[1:3] == (1356, 191, 1356)
        outtxt = readdlm(GzipDecompressorStream(open("gtstats.out.gz", "r")))
        @test all(outtxt[:, 9]  .== out[5]) # missings_by_record
        @test all(outtxt[:, 14] .== out[6]) # maf_by_record
    end

    @testset "input: gz file, output: text file" begin
        @time out = gtstats("test.08Jun17.d8b.vcf.gz", "gtstats.out.txt")
        @test out[1:3] == (1356, 191, 1356)
        outtxt = readdlm("gtstats.out.txt")
        @test all(outtxt[:, 9]  .== out[5]) # missings_by_record
        @test all(outtxt[:, 14] .== out[6]) # maf_by_record
    end

    @testset "input: gz file, output: gz file" begin
        @time out = gtstats("test.08Jun17.d8b.vcf.gz", "gtstats.out.gz")
        @test out[1:3] == (1356, 191, 1356)
        outtxt = readdlm(GzipDecompressorStream(open("gtstats.out.gz", "r")))
        @test all(outtxt[:, 9]  .== out[5]) # missings_by_record
        @test all(outtxt[:, 14] .== out[6]) # maf_by_record
    end

    #@time records, samples, lines = gtstats("test.08Jun17.d8b.vcf.gz", "gtstats.out.gz")
    #@time records, samples, lines = gtstats("chr22.1kg.ref.phase1_release_v3.20101123.vcf.gz", "gtstats.out.gz") # about 180 seconds
    #@show maf_by_record
end


text:
"""
    colorbar(parent, plotobject, kwargs...)

Add a colorbar based on the attribute chose to color the plot. plotobject must be a plot of
an MTG colored by an attribute. Use Makie.Colorbar for any other use case instead.
"""
function colorbar(parent, plotobject; kwargs...)
    color = plotobject.attributes.color[]
    mtg = plotobject.converted[1][]

    if !(typeof(mtg) <: MultiScaleTreeGraph.Node)
        error("This is not a plot of an MTG. Use Makie.Colorbar instead.")
    end

    if !(color in get_attributes(mtg))
        error(
            "The plot must be colored by an MTG attribute for making a colorbar.",
            "Use Makie.Colorbar instead."
        )
    end

    range_val = attribute_range(plotobject.converted[1][], plotobject.attributes.color[])

    Makie.Colorbar(
        parent,
        label = string(plotobject.attributes.color[]),
        colormap = plotobject.attributes.colormap,
        limits = Float64.(range_val);
        kwargs...
    )
end


text:
# ---
# title: 1504. Count Submatrices With All Ones
# id: problem1504
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/count-submatrices-with-all-ones/description/>
# hidden: true
# ---
#
# Given a `rows * columns` matrix `mat` of ones and zeros, return how many
# **submatrices** have all ones.
#
#
#
# **Example 1:**
#
#
#
#     Input: mat = [[1,0,1],
#                   [1,1,0],
#                   [1,1,0]]
#     Output: 13
#     Explanation: There are **6** rectangles of side 1x1.
#     There are **2** rectangles of side 1x2.
#     There are **3** rectangles of side 2x1.
#     There is **1** rectangle of side 2x2.
#     There is **1** rectangle of side 3x1.
#     Total number of rectangles = 6 + 2 + 3 + 1 + 1 = **13.**
#
#
# **Example 2:**
#
#
#
#     Input: mat = [[0,1,1,0],
#                   [0,1,1,1],
#                   [1,1,1,0]]
#     Output: 24
#     Explanation:
#     There are **8** rectangles of side 1x1.
#     There are **5** rectangles of side 1x2.
#     There are **2** rectangles of side 1x3.
#     There are **4** rectangles of side 2x1.
#     There are **2** rectangles of side 2x2.
#     There are **2** rectangles of side 3x1.
#     There is **1** rectangle of side 3x2.
#     Total number of rectangles = 8 + 5 + 2 + 4 + 2 + 2 + 1 = 24 **.**
#
#
# **Example 3:**
#
#
#
#     Input: mat = [[1,1,1,1,1,1]]
#     Output: 21
#
#
# **Example 4:**
#
#
#
#     Input: mat = [[1,0,1],[0,1,0],[1,0,1]]
#     Output: 5
#
#
#
#
# **Constraints:**
#
#   * `1 <= rows <= 150`
#   * `1 <= columns <= 150`
#   * `0 <= mat[i][j] <= 1`
#
#
## @lc code=start
using LeetCode

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


text:
include("applications.jl")

using Fire

@main function sample_cluster()
    m = model(Val(:cluster); nbit=3, B=10)
    @show gensample(m, X)
end

"""
    julia j1j2.jl train [--symmetry <su2|u1|general>] [--depth <Int>]

Train a 4 x 4 frustrated Heisenberg model with J2 = 0.5.
Available ansatz symmetries includes `general`, `u1` and `su2`.
"""
@main function train(;symmetry::Symbol=:su2, depth::Int=5)
    USE_CUDA || @warn "You are not using GPU (35 x speed up), this training may take life long. Turn on the switch in file `applications.jl` if you have a GPU that supports CUDA!"
    model = simple_model_j1j2(4, 4)
    ansatz = simple_ansatz(16, symmetry, depth; load_params=false)

    run_train(ansatz, model; SAVE_ID=Symbol(symmetry,:_d,depth), niter=500, start_point=0)
end

"""
    julia j1j2.jl measure <energy|fidelity|szsz> [--symmetry <su2|u1|general>] [--depth <Int>]

Load a pre-trained ansatz for 4 x 4 frustrated Heisenberg model with J2 = 0.5, tasks includes
* energy, sample the energy.
* szsz, calculate the <sz(i)*sz(j)> correlation matrix.

Also, we can obtain some exact quantities in simulation for analysis
* fidelity, calculated the exact fidelity
* energy_exact, calculated the exact energy

Pre-trained ansatz includes
* --symmetry su2, --depth <1-5>
* --symmetry u1, --depth 5
* --symmetry general, --depth 5
"""
@main function measure(task::String; symmetry::Symbol=:su2, depth::Int=5)
    model = simple_model_j1j2(4, 4)
    ansatz = simple_ansatz(16, symmetry, depth; load_params=true)
    nbit = nbit_simulated(ansatz)

    if task == "energy"
        eng_sample = energy(ansatz, model)/nbit
        println("Sampled value of energy/site = $eng_sample")
    elseif task == "energy_exact"
        eng = energy_exact(ansatz, model)/nbit
        println("Exact value of energy/site = $eng")
    elseif task == "fidelity"
        gs = ground_state(model)
        fid = fidelity_exact(ansatz, gs)[]
        println("Exact value of fidelity = $fid")
    elseif task == "szsz"
        correlation_matrix(ansatz; SAVE_ID=symmetry)
    else
        throw(ArgumentError("$task is not defined!"))
    end
end

max_stars_repo_name: GiggleLiu/QuantumMPS

max_issues_repo_name: GiggleLiu/QuantumMPS

max_forks_repo_name: GiggleLiu/QuantumMPS


text:
"""
Redirects stdout to a file.

# Example

```julia
import Util

Util.redirect_stdout_to_file("log.txt") do
  for i in 1:20
    sleep(.1)
    print("#")
    flush(stdout)
  end
  println()
end
```
"""
function redirect_stdout_to_file(f::Function, path::String)
  open(path, "w") do io
    redirect_stdout(io) do
      f()
    end
  end
end


text:
module Utils
include("testtasks.jl")
include("mockservices.jl")
include("testservices.jl")
end # module Utils

max_forks_repo_name: seung-lab/SimpleTasks


text:
type FormatHints
    width         :: Int  # column width, not the format width
    scale         :: Real
    precision     :: Int
    commas        :: Bool
    stripzeros    :: Bool
    parens        :: Bool
    rednegative   :: Bool # print in red when negative?
    hidezero      :: Bool
    alternative   :: Bool
    mixedfraction :: Bool
    suffix        :: UTF8String
    autoscale     :: Symbol
    conversion    :: ASCIIString
end

function FormatHints{T<:Integer}( ::Type{T} )
    FormatHints( 8, 1, 0, true, false, false, true, true, false, false, "", :none, "d" )
end
function FormatHints{T<:Unsigned}( ::Type{T} )
    FormatHints( 8, 1, 0, true, false, false, true, true, false, false, "", :none, "x" )
end
function FormatHints{T<:AbstractFloat}( ::Type{T} )
    FormatHints( 10, 1.0, 2, true, false, false, true, true, false, false, "", :none, "f" )
end
function FormatHints{T<:Rational}( ::Type{T} )
    FormatHints( 12, 1, 0, false, false, false, true, true, false, true, "", :none, "s" )
end
function FormatHints( ::Type{Date} )
    FormatHints( 10, 1, 0, false, false, false, false, false, false, false,"", :none,
       "yyyy-mm-dd" )
end
function FormatHints( ::Type{DateTime} )
    FormatHints( 20, 1, 0, false, false, false, false, false, false, false,"", :none,
       "yyyy-mm-dd HH:MM:SS" )
end
function FormatHints( ::Type{} )
    FormatHints( 14, 1, 0, false, false, false, true, true, false, false, "", :none, "s" )
end

function applyformat{T<:Number}( v::T, fmt::FormatHints )
    if fmt.hidezero && v == 0
        ""
    else
        format( v * fmt.scale,
            precision     = fmt.precision,
            commas        = fmt.commas,
            stripzeros    = fmt.stripzeros,
            parens        = fmt.parens,
            alternative   = fmt.alternative,
            mixedfraction = fmt.mixedfraction,
            suffix        = fmt.suffix,
            autoscale     = fmt.autoscale,
            conversion    = fmt.conversion
            )
    end
end

function applyformat( v::Union{Date,DateTime}, fmt::FormatHints )
    Dates.format( v, fmt.conversion )
end

function applyformat{T<:AbstractString}( v::T, fmt::FormatHints )
    @lintpragma( "Ignore unused fmt" )
    return v
end

function applyformat( v::AbstractArray, fmt::FormatHints )
    strs = UTF8String[]
    for s in v
        push!( strs, applyformat( s, fmt ) )
    end
    join( strs, "," )
end

function applyformat( v, fmt::FormatHints )
    @lintpragma( "Ignore unused fmt" )
    return string( v )
end


text:
using AbstractGPs
using Distributions
using FiniteDifferences
using LinearAlgebra
using OILMMs
using Random
using Stheno
using Test
using Zygote

using AbstractGPs: AbstractGP, FiniteGP
using OILMMs: denoised_marginals, rand_latent
using Stheno: GaussianProcessProbabilisticProgramme

# Helper functionality, doesn't actually run any tests.
include("test_util.jl")

@testset "OILMMs.jl" begin
    include("util.jl")
    include("oilmm.jl")
    include("missing_data.jl")
end


text:
using UnFortran
using Base.Test

filename = joinpath(mktempdir(), "test")
v = (1,2,3,1pi)
m,n = 5,4
a = rand(m,n)

fort_open(filename, "w") do file
    write(file, v...)
    write(file, m, n)
    write(file, a)
end

fort_open(filename) do file
    @test (read(file, Int64, Int64, Int64, Float64)...) == v
    mi,ni = read(file, Int64, Int64)
    @test (mi,ni) == (m,n)
    ai = Array(Float64, m, n)
    read!(file, ai)
    @test ai == a
end

# Test reading record of wrong size
fort_open(filename) do file
    @test_throws ErrorException read(file, Int64)
end


text:
"""
    save_input_mat(case)

Read the original case.mat file, replace relevant parameters from Case struct,
save a new input.mat file with parameters as they're passed to solver.
"""
function save_input_mat(
    case::Case, storage::Storage, inputfolder::String, outputfolder::String
)
    # MATPOWER column indices
    gen_PMAX = 9
    gen_PMIN = 10
    gen_RAMP_30 = 19
    gencost_MODEL = 1
    gencost_NCOST = 4
    gencost_COST = 5

    # Read original
    case_mat_file = MAT.matopen(joinpath(inputfolder, "case.mat"))
    mpc = read(case_mat_file, "mpc")
    mdi = Dict("mpc" => mpc)

    # Save modifications to gen
    mpc["gen"][:, gen_PMIN] = case.gen_pmin
    mpc["gen"][:, gen_RAMP_30] = case.gen_ramp30
    # Save modifications to gencost table
    mpc["gencost"] = case.gencost
    mpc["gencost_orig"] = case.gencost_orig

    # Save storage details in mpc.[gencost, genfuel, gencost] and in 'Storage'
    if size(storage.gen, 1) > 0
        num_storage = size(storage.gen, 1)
        # Add fuel types for storage 'generators'. ESS = energy storage system
        mpc["genfuel"] = [mpc["genfuel"]; repeat(["ess"], num_storage)]
        # Save storage modifications to gen table
        mpc["gen"] = [mpc["gen"]; storage.gen]
        # Save storage modifications to gencost
        #@show case.gencost[:, gencost_MODEL]
        segment_indices = findall(case.gencost[:, gencost_MODEL] .== 1)
        segment_orders = case.gencost[segment_indices, gencost_NCOST]
        num_segments = convert(Int, maximum(segment_orders)) - 1
        storage_gencost = zeros(num_storage, (6 + 2 * num_segments))
        # Storage is specified by two points, PMIN and PMAX, both with cost 0
        storage_gencost[:, gencost_MODEL] .= 1
        storage_gencost[:, gencost_NCOST] .= 2
        storage_gencost[:, gencost_COST] = storage.gen[:, gen_PMIN]
        storage_gencost[:, gencost_COST + 2] = storage.gen[:, gen_PMAX]
        mpc["gencost"] = [mpc["gencost"]; storage_gencost]
        # Save addition of 'iess' field (index, energy storage systems) to mpc
        num_gen = length(case.genid)
        # Mimic the array that MATPOWER/MOST would create
        iess = collect((num_gen + 1):(num_gen + num_storage))
        mpc["iess"] = iess
        # Build new struct for data in 'Storage'
        input_mat_storage = Dict(
            String(s) => storage.sd_table[!, s] for s in DataFrames.names(storage.sd_table)
        )
        mdi["Storage"] = input_mat_storage
    end

    output_path = joinpath(outputfolder, "input.mat")
    MAT.matwrite(output_path, Dict("mdi" => mdi); compress=true)
    return nothing
end

"""Given a Results object and a filename, save a matfile with results data."""
function save_results(results::Results, filename::String; demand_scaling::Number=1.0)
    mdo_save = Dict(
        "results" => Dict("f" => results.f),
        "demand_scaling" => demand_scaling,
        "flow" => Dict(
            "mpc" => Dict(
                "bus" => Dict("LAM_P" => results.lmp),
                "gen" => Dict("PG" => results.pg),
                "branch" => Dict(
                    "PF" => results.pf,
                    "MU_SF" => results.congu,
                    "MU_ST" => results.congl,
                ),
            ),
        ),
    )
    # For DC lines, storage power/energy, load_shed, and flexible demand,
    # save only if nonempty
    if size(results.pf_dcline) != (0, 0)
        mdo_save["flow"]["mpc"]["dcline"] = Dict("PF_dcline" => results.pf_dcline)
    end
    if size(results.storage_pg) != (0, 0)
        mdo_save["flow"]["mpc"]["storage"] = Dict(
            "PG" => results.storage_pg, "Energy" => results.storage_e
        )
    end
    if size(results.load_shed) != (0, 0)
        mdo_save["flow"]["mpc"]["load_shed"] = Dict("load_shed" => results.load_shed)
    end
    if size(results.load_shift_dn) != (0, 0)
        mdo_save["flow"]["mpc"]["flexible_demand"] = Dict(
            "load_shift_up" => results.load_shift_up,
            "load_shift_dn" => results.load_shift_dn,
        )
    end
    return MAT.matwrite(filename, Dict("mdo_save" => mdo_save); compress=true)
end

"""
    redirect_stdout_stderr("stdout.log", "stderr.err") do
        run_scenario(; kwargs...)
    end

While executing a function, redirect stdout and stderr to files.
"""
function redirect_stdout_stderr(dofunc, stdoutfile, stderrfile)
    open(stdoutfile, "a") do out
        open(stderrfile, "a") do err
            redirect_stdout(out) do
                redirect_stderr(err) do
                    dofunc()
                end
            end
        end
    end
end


text:
module AbstractType

export analyticalSolution

export GoverningDiffEqn, Analysis, Approach
export Dimensions, Material
export BoundaryCondition, Source
export ResultSpec

export Err_InputCombo

abstract type GoverningDiffEqn  end
abstract type Analysis          end
abstract type Approach          end
abstract type Dimensions        end
abstract type Material          end
abstract type BoundaryCondition end
abstract type Source            end
abstract type ResultSpec        end

const Err_InputCombo = "Analytical solution not yet implemented for requested input combination"

function analyticalSolution(g::GoverningDiffEqn, an::Analysis, ap::Approach,
                            d::Dimensions, m::Material,
                            b::BoundaryCondition, s::Source,
                            r::ResultSpec)
  error(Err_InputCombo)
end

end


text:
push!(LOAD_PATH,"../src")

using Test
using TypesParticles
using TypesBasis
using BuildBasis
using BuildOverlap
using InputParser

"""
		run_buildoverlap_test()

This is a function to test the user types in BuildOverlap.jl. In asserts overlap integral calculationfor two hydrogen atoms is the same as the output from Prof. James John's matlab code.

""" function run_buildoverlap_test()
    @testset "Overlap" begin
        system,basisset = getcalcsetup("../examples/HatreeFock.in")
        basisfunc,numelec = buildbasisfunc(system,basisset);
        Scalc = buildelecoverlap(numelec,basisfunc,basisset)
        Stest=  [1.00000e0 0.65829e0 0.77339e0 0.60892e0;
                 0.65829e0 1.00000e0 0.60892e0 0.95331e0;
                 0.77339e0 0.60892e0 1.00000e0 0.65829e0;
                 0.60892e0 0.95331e0 0.65829e0 1.00000e0]
        #Loop over and check each entry
        tolerance = 1.0e-5;
        for i=1:4
            for j=1:4
                @test begin
                    sign(Scalc[i,j]) == sign(Stest[i,j]);
                    diff = Scalc[i,j]-Stest[i,j];
                    diff <= tolerance;
                end
	    end
        end

    end #testset
end #run_overlap_test()


text:
function find_ksym(addr::UInt64)
    start, stop = UInt64(0), addr
    rgx = r"([0-9abcdef]*) ([a-zA-Z]) ([0-9a-zA-Z_\-]*)"
    last_addr = start
    last_kind = "?"
    last_name = ""
    for line in readlines(open("/proc/kallsyms", "r"))
        m = match(rgx, line)
        @assert m !== nothing
        start_addr, kind, name = m.captures
        start_addr = parse(UInt64, "0x"*start_addr)
        if start_addr > stop
            return last_addr, last_kind, last_name
        elseif start_addr == stop
            return addr, kind, name
        end
        last_addr = addr
        last_kind = kind
        last_name = name
    end
end
function stack_to_string(nt::NTuple{N,UInt64}) where N
    iob = IOBuffer()
    for i in 1:N
        addr = nt[i]
        if addr == UInt64(0)
            break
        end
        println(iob, "$(find_ksym(addr)[3])")
    end
    String(take!(iob))
end


text:
function list_tarball(
    tar::IO;
    strict::Bool = true,
    buf::Vector{UInt8} = Vector{UInt8}(undef, 512),
)
    headers = Header[]
    while !eof(tar)
        hdr = read_header(tar, buf=buf)
        hdr === nothing && break
        strict && check_header(hdr)
        push!(headers, hdr)
        skip(tar, 512 * ((hdr.size + 511) ÷ 512))
    end
    return headers
end

function extract_tarball(
    tar::IO,
    root::String;
    buf::Vector{UInt8} = Vector{UInt8}(undef, 512),
)
    links = Set{String}()
    while !eof(tar)
        hdr = read_header(tar, buf=buf)
        hdr === nothing && break
        check_header(hdr)
        # normalize path and check for symlink attacks
        path = ""
        parts = String[]
        for part in split(hdr.path, '/')
            (isempty(part) || part == ".") && continue
            path in links && error("""
            refusing to extract path with symlink prefix, possible attack
             * symlink prefix: $(repr(path))
             * path: $(repr(hdr.path))
            """)
            path = isempty(path) ? part : "$path/$part"
            push!(parts, part)
        end
        if hdr.type == :symlink
            push!(links, path)
        end
        sys_path = joinpath(root, parts...)
        # create the path
        if hdr.type == :directory
            mkpath(sys_path)
        else
            if ispath(sys_path)
                # delete and replace path
                rm(sys_path, force=true, recursive=true)
            else
                dir = dirname(sys_path)
                # ensure `dir` is a directory
                st = stat(dir)
                if !isdir(st)
                    ispath(st) && rm(dir, force=true, recursive=true)
                    mkpath(dir)
                end
            end
            hdr.type == :file && read_data(tar, sys_path, size=hdr.size)
            hdr.type == :symlink && symlink(hdr.link, sys_path)
        end
        if hdr.type != :symlink
            chmod(sys_path, hdr.mode)
        elseif Sys.isbsd()
            # BSD system support symlink permissions, so try setting them...
            ret = ccall(:lchmod, Cint, (Cstring, Base.Cmode_t), sys_path, hdr.mode)
            systemerror(:lchmod, ret != 0)
        end
    end
end

function read_header(io::IO; buf::Vector{UInt8} = Vector{UInt8}(undef, 512))
    hdr = read_standard_header(io, buf=buf)
    hdr === nothing && return nothing
    hdr.type in (:x, :g) || return hdr
    size = path = link = nothing
    metadata = read_extended_metadata(io, hdr.size, buf=buf)
    for (key, value) in metadata
        if key == "size"
            size = tryparse(UInt64, value)
            size === nothing &&
                error("invalid extended header size value: $(repr(value))")
        elseif key == "path"
            path = value
        elseif key == "linkpath"
            link = value
        end
    end
    hdr = read_standard_header(io, buf=buf)
    hdr === nothing && error("premature end of tar file")
    return Header(
        something(path, hdr.path),
        hdr.type, hdr.mode,
        something(size, hdr.size),
        something(link, hdr.link),
    )
end

using Base.Checked: mul_with_overflow, add_with_overflow

function read_extended_metadata(
    io::IO,
    size::Integer;
    buf::Vector{UInt8} = Vector{UInt8}(undef, 512),
)
    n = readbytes!(io, buf, size)
    n < size && "premature end of tar file"
    skip(io, mod(512 - n, 512)) # advance to end of block
    malformed() = error("malformed extended header metadata: $(repr(String(buf)))")
    metadata = Pair{String,String}[]
    i = 0
    while i < n
        j, m = i, 0
        while j ≤ n
            byte = buf[j += 1]
            byte == UInt8(' ') && break
            UInt8('0') ≤ byte ≤ UInt8('9') || malformed()
            m, fm = mul_with_overflow(m, 10)
            m, fa = add_with_overflow(m, Int(byte - UInt8('0')))
            fm | fa && error("extended header record size too large: $(repr(buf))")
        end
        k, l = j, i + m
        while k ≤ l
            byte = buf[k += 1]
            byte == UInt8('=') && break
        end
        # buf[i+1:j-1] is the length in decimal
        # buf[j+1:k-1] is the key string
        # buf[k+1:l-1] is the value string
        # buf[i] is end of previous
        # buf[j] is ` ` (space)
        # buf[k] is `=` (equals)
        # buf[l] is `\n` (newline)
        i+1 < j < k < l || malformed()
        @assert buf[j] == UInt8(' ')
        @assert buf[k] == UInt8('=')
        buf[l] == UInt('\n') || malformed()
        @views push!(metadata, String(buf[j+1:k-1]) => String(buf[k+1:l-1]))
        i = l
    end
    return metadata
end

function read_standard_header(io::IO; buf::Vector{UInt8} = Vector{UInt8}(undef, 512))
    resize!(buf, 512)
    read!(io, buf)
    all(iszero, buf) && return nothing
    n = length(buf)
    n == 0 && error("premature end of tar file")
    n < 512 && error("incomplete trailing block with length $n < 512")
    @assert n == 512
    name    = read_header_str(buf, 0, 100)
    mode    = read_header_int(buf, 100, 8)
    size    = buf[124+1] & 0x80 == 0 ?
              read_header_int(buf, 124, 12) :
              read_header_bin(buf, 124, 12)
    chksum  = read_header_int(buf, 148, 8)
    type    = read_header_chr(buf, 156)
    link    = read_header_str(buf, 157, 100)
    magic   = read_header_str(buf, 257, 6)
    version = read_header_str(buf, 263, 2)
    prefix  = read_header_str(buf, 345, 155)
    # check various fields
    buf[index_range(148, 8)] .= ' ' # fill checksum field with spaces
    buf_sum = sum(buf)
    chksum == buf_sum ||
        error("incorrect header checksum = $chksum; should be $buf_sum\n$(repr(String(buf)))")
    occursin(r"^ustar\s*$", magic) ||
        error("unknown magic string for tar file: $(repr(magic))")
    occursin(r"^0* *$", version) ||
        error("unkonwn version string for tar file: $(repr(version))")
    path = isempty(prefix) ? name : "$prefix/$name"
    return Header(path, to_symbolic_type(type), mode, size, link)
end

index_range(offset::Int, length::Int) = offset .+ (1:length)

read_header_chr(buf::Vector{UInt8}, offset::Int) = Char(buf[offset+1])

function read_header_str(buf::Vector{UInt8}, offset::Int, length::Int)
    r = index_range(offset, length)
    for i in r
        byte = buf[i]
        byte == 0 && return String(buf[first(r):i-1])
    end
    return String(buf[r])
end

function read_header_int(buf::Vector{UInt8}, offset::Int, length::Int)
    n = UInt64(0)
    for i in index_range(offset, length)
        byte = buf[i]
        byte in (0x00, UInt8(' ')) && break
        UInt8('0') <= byte <= UInt8('7') ||
            error("invalid octal digit: $(repr(Char(byte)))")
        n <<= 3
        n |= byte - 0x30
    end
    return n
end

function read_header_bin(buf::Vector{UInt8}, offset::Int, length::Int)
    n = UInt64(0)
    for i in index_range(offset, length)
        n <<= 8
        n |= buf[i]
    end
    return n
end

function read_data(
    tar::IO,
    file::IO;
    size::Integer,
    buf::Vector{UInt8} = Vector{UInt8}(undef, 512),
)
    resize!(buf, 512)
    while size > 0
        r = readbytes!(tar, buf)
        r < 512 && eof(io) && error("premature end of tar file")
        size < 512 && resize!(buf, size)
        size -= write(file, buf)
    end
    resize!(buf, 512)
    @assert size == 0
    return
end

function read_data(
    tar::IO,
    file::String;
    size::Integer,
    buf::Vector{UInt8} = Vector{UInt8}(undef, 512),
)
    open(file, write=true) do file′
        read_data(tar, file′, size=size, buf=buf)
    end
end


text:
immutable CoordinateSystem
    center::XYZ
    rotation_matrix::Mat3x3
end

function CoordinateSystem(pc::XYZ,pz::XYZ,pxz::XYZ)
    vecz = normalize(pz - pc)
    vecy = normalize(cross(vecz,pxz-pc))
    vecx = normalize(cross(vecy,vecz))
    CoordinateSystem(pc,Mat3x3([Array(vecx) Array(vecy) Array(vecz)]))
end

rotate_from_global(csys::CoordinateSystem,p::XYZ) = csys.rotation_matrix' * p
rotate_to_global(csys::CoordinateSystem,p::XYZ) = csys.rotation_matrix * p

rotate_from_global(csys::CoordinateSystem,p::Mat3x3) = csys.rotation_matrix' * p * csys.rotation_matrix
rotate_to_global(csys::CoordinateSystem,p::Mat3x3) = csys.rotation_matrix * p * csys.rotation_matrix'

from_global(csys::CoordinateSystem,p::XYZ) = rotate_from_global(csys,p-csys.center)
to_global(csys::CoordinateSystem,p::XYZ) = rotate_to_global(csys,p) + csys.center


text:
push!(LOAD_PATH, "@stdlib")
import Pkg; Pkg.add("Conda");
using Conda
using ADCME

@info "Install Boost"
CONDA = get_conda()
run(`$CONDA install boost`)

@info "Install AMGCL"
UNZIP = joinpath(ADCME.BINDIR, "unzip")
if !isdir("$(@__DIR__)/amgcl")
    download("https://github.com/ddemidov/amgcl/archive/master.zip", "$(@__DIR__)/amgcl.zip")
    run(`$UNZIP -o $(@__DIR__)/amgcl.zip -d $(@__DIR__)`)
    mv("$(@__DIR__)/amgcl-master","$(@__DIR__)/amgcl", force=true)
    rm("$(@__DIR__)/amgcl.zip")
end

@info "Build Custom Operators"
change_directory("CustomOps/build")
require_file("build.ninja") do
    ADCME.cmake()
end
ADCME.make()


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: chime_sdk_meetings
using AWS.Compat
using AWS.UUIDs

"""
    batch_create_attendee(attendees, meeting_id)
    batch_create_attendee(attendees, meeting_id, params::Dict{String,<:Any})

Creates up to 100 attendees for an active Amazon Chime SDK meeting. For more information
about the Amazon Chime SDK, see Using the Amazon Chime SDK in the Amazon Chime Developer
Guide.

# Arguments
- `attendees`: The attendee information, including attendees' IDs and join tokens.
- `meeting_id`: The Amazon Chime SDK ID of the meeting to which you're adding attendees.

"""
function batch_create_attendee(
    Attendees, MeetingId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return chime_sdk_meetings(
        "POST",
        "/meetings/$(MeetingId)/attendees?operation=batch-create",
        Dict{String,Any}("Attendees" => Attendees);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function batch_create_attendee(
    Attendees,
    MeetingId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return chime_sdk_meetings(
        "POST",
        "/meetings/$(MeetingId)/attendees?operation=batch-create",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("Attendees" => Attendees), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_attendee(external_user_id, meeting_id)
    create_attendee(external_user_id, meeting_id, params::Dict{String,<:Any})

 Creates a new attendee for an active Amazon Chime SDK meeting. For more information about
the Amazon Chime SDK, see Using the Amazon Chime SDK in the Amazon Chime Developer Guide.

# Arguments
- `external_user_id`: The Amazon Chime SDK external user ID. An idempotency token. Links
  the attendee to an identity managed by a builder application.
- `meeting_id`: The unique ID of the meeting.

"""
function create_attendee(
    ExternalUserId, MeetingId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return chime_sdk_meetings(
        "POST",
        "/meetings/$(MeetingId)/attendees",
        Dict{String,Any}("ExternalUserId" => ExternalUserId);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_attendee(
    ExternalUserId,
    MeetingId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return chime_sdk_meetings(
        "POST",
        "/meetings/$(MeetingId)/attendees",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ExternalUserId" => ExternalUserId), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_meeting(client_request_token, external_meeting_id, media_region)
    create_meeting(client_request_token, external_meeting_id, media_region, params::Dict{String,<:Any})

Creates a new Amazon Chime SDK meeting in the specified media Region with no initial
attendees. For more information about specifying media Regions, see Amazon Chime SDK Media
Regions in the Amazon Chime Developer Guide. For more information about the Amazon Chime
SDK, see Using the Amazon Chime SDK in the Amazon Chime Developer Guide.

# Arguments
- `client_request_token`: The unique identifier for the client request. Use a different
  token for different meetings.
- `external_meeting_id`: The external meeting ID.
- `media_region`: The Region in which to create the meeting.  Available values: af-south-1
  , ap-northeast-1 , ap-northeast-2 , ap-south-1 , ap-southeast-1 , ap-southeast-2 ,
  ca-central-1 , eu-central-1 , eu-north-1 , eu-south-1 , eu-west-1 , eu-west-2 , eu-west-3 ,
  sa-east-1 , us-east-1 , us-east-2 , us-west-1 , us-west-2 .

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MeetingFeatures"`: Lists the audio and video features enabled for a meeting, such as
  echo reduction.
- `"MeetingHostId"`: Reserved.
- `"NotificationsConfiguration"`: The configuration for resource targets to receive
  notifications when meeting and attendee events occur.
"""
function create_meeting(
    ClientRequestToken,
    ExternalMeetingId,
    MediaRegion;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return chime_sdk_meetings(
        "POST",
        "/meetings",
        Dict{String,Any}(
            "ClientRequestToken" => ClientRequestToken,
            "ExternalMeetingId" => ExternalMeetingId,
            "MediaRegion" => MediaRegion,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_meeting(
    ClientRequestToken,
    ExternalMeetingId,
    MediaRegion,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return chime_sdk_meetings(
        "POST",
        "/meetings",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ClientRequestToken" => ClientRequestToken,
                    "ExternalMeetingId" => ExternalMeetingId,
                    "MediaRegion" => MediaRegion,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_meeting_with_attendees(attendees, client_request_token, external_meeting_id, media_region)
    create_meeting_with_attendees(attendees, client_request_token, external_meeting_id, media_region, params::Dict{String,<:Any})

 Creates a new Amazon Chime SDK meeting in the specified media Region, with attendees. For
more information about specifying media Regions, see Amazon Chime SDK Media Regions in the
Amazon Chime Developer Guide. For more information about the Amazon Chime SDK, see Using
the Amazon Chime SDK in the Amazon Chime Developer Guide.

# Arguments
- `attendees`: The attendee information, including attendees' IDs and join tokens.
- `client_request_token`: The unique identifier for the client request. Use a different
  token for different meetings.
- `external_meeting_id`: The external meeting ID.
- `media_region`: The Region in which to create the meeting.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MeetingFeatures"`: Lists the audio and video features enabled for a meeting, such as
  echo reduction.
- `"MeetingHostId"`: Reserved.
- `"NotificationsConfiguration"`: The configuration for resource targets to receive
  notifications when meeting and attendee events occur.
"""
function create_meeting_with_attendees(
    Attendees,
    ClientRequestToken,
    ExternalMeetingId,
    MediaRegion;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return chime_sdk_meetings(
        "POST",
        "/meetings?operation=create-attendees",
        Dict{String,Any}(
            "Attendees" => Attendees,
            "ClientRequestToken" => ClientRequestToken,
            "ExternalMeetingId" => ExternalMeetingId,
            "MediaRegion" => MediaRegion,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_meeting_with_attendees(
    Attendees,
    ClientRequestToken,
    ExternalMeetingId,
    MediaRegion,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return chime_sdk_meetings(
        "POST",
        "/meetings?operation=create-attendees",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "Attendees" => Attendees,
                    "ClientRequestToken" => ClientRequestToken,
                    "ExternalMeetingId" => ExternalMeetingId,
                    "MediaRegion" => MediaRegion,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_attendee(attendee_id, meeting_id)
    delete_attendee(attendee_id, meeting_id, params::Dict{String,<:Any})

Deletes an attendee from the specified Amazon Chime SDK meeting and deletes their
JoinToken. Attendees are automatically deleted when a Amazon Chime SDK meeting is deleted.
For more information about the Amazon Chime SDK, see Using the Amazon Chime SDK in the
Amazon Chime Developer Guide.

# Arguments
- `attendee_id`: The Amazon Chime SDK attendee ID.
- `meeting_id`: The Amazon Chime SDK meeting ID.

"""
function delete_attendee(
    AttendeeId, MeetingId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return chime_sdk_meetings(
        "DELETE",
        "/meetings/$(MeetingId)/attendees/$(AttendeeId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_attendee(
    AttendeeId,
    MeetingId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return chime_sdk_meetings(
        "DELETE",
        "/meetings/$(MeetingId)/attendees/$(AttendeeId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_meeting(meeting_id)
    delete_meeting(meeting_id, params::Dict{String,<:Any})

Deletes the specified Amazon Chime SDK meeting. The operation deletes all attendees,
disconnects all clients, and prevents new clients from joining the meeting. For more
information about the Amazon Chime SDK, see Using the Amazon Chime SDK in the Amazon Chime
Developer Guide.

# Arguments
- `meeting_id`: The Amazon Chime SDK meeting ID.

"""
function delete_meeting(MeetingId; aws_config::AbstractAWSConfig=global_aws_config())
    return chime_sdk_meetings(
        "DELETE",
        "/meetings/$(MeetingId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_meeting(
    MeetingId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return chime_sdk_meetings(
        "DELETE",
        "/meetings/$(MeetingId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_attendee(attendee_id, meeting_id)
    get_attendee(attendee_id, meeting_id, params::Dict{String,<:Any})

 Gets the Amazon Chime SDK attendee details for a specified meeting ID and attendee ID. For
more information about the Amazon Chime SDK, see Using the Amazon Chime SDK in the Amazon
Chime Developer Guide.

# Arguments
- `attendee_id`: The Amazon Chime SDK attendee ID.
- `meeting_id`: The Amazon Chime SDK meeting ID.

"""
function get_attendee(
    AttendeeId, MeetingId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return chime_sdk_meetings(
        "GET",
        "/meetings/$(MeetingId)/attendees/$(AttendeeId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_attendee(
    AttendeeId,
    MeetingId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return chime_sdk_meetings(
        "GET",
        "/meetings/$(MeetingId)/attendees/$(AttendeeId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_meeting(meeting_id)
    get_meeting(meeting_id, params::Dict{String,<:Any})

Gets the Amazon Chime SDK meeting details for the specified meeting ID. For more
information about the Amazon Chime SDK, see Using the Amazon Chime SDK in the Amazon Chime
Developer Guide.

# Arguments
- `meeting_id`: The Amazon Chime SDK meeting ID.

"""
function get_meeting(MeetingId; aws_config::AbstractAWSConfig=global_aws_config())
    return chime_sdk_meetings(
        "GET",
        "/meetings/$(MeetingId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_meeting(
    MeetingId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return chime_sdk_meetings(
        "GET",
        "/meetings/$(MeetingId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_attendees(meeting_id)
    list_attendees(meeting_id, params::Dict{String,<:Any})

 Lists the attendees for the specified Amazon Chime SDK meeting. For more information about
the Amazon Chime SDK, see Using the Amazon Chime SDK in the Amazon Chime Developer Guide.

# Arguments
- `meeting_id`: The Amazon Chime SDK meeting ID.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"max-results"`: The maximum number of results to return in a single call.
- `"next-token"`: The token to use to retrieve the next page of results.
"""
function list_attendees(MeetingId; aws_config::AbstractAWSConfig=global_aws_config())
    return chime_sdk_meetings(
        "GET",
        "/meetings/$(MeetingId)/attendees";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_attendees(
    MeetingId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return chime_sdk_meetings(
        "GET",
        "/meetings/$(MeetingId)/attendees",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    start_meeting_transcription(meeting_id, transcription_configuration)
    start_meeting_transcription(meeting_id, transcription_configuration, params::Dict{String,<:Any})

Starts transcription for the specified meetingId.

# Arguments
- `meeting_id`: The unique ID of the meeting being transcribed.
- `transcription_configuration`: The configuration for the current transcription operation.
  Must contain EngineTranscribeSettings or EngineTranscribeMedicalSettings.

"""
function start_meeting_transcription(
    MeetingId, TranscriptionConfiguration; aws_config::AbstractAWSConfig=global_aws_config()
)
    return chime_sdk_meetings(
        "POST",
        "/meetings/$(MeetingId)/transcription?operation=start",
        Dict{String,Any}("TranscriptionConfiguration" => TranscriptionConfiguration);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function start_meeting_transcription(
    MeetingId,
    TranscriptionConfiguration,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return chime_sdk_meetings(
        "POST",
        "/meetings/$(MeetingId)/transcription?operation=start",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "TranscriptionConfiguration" => TranscriptionConfiguration
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    stop_meeting_transcription(meeting_id)
    stop_meeting_transcription(meeting_id, params::Dict{String,<:Any})

Stops transcription for the specified meetingId.

# Arguments
- `meeting_id`: The unique ID of the meeting for which you stop transcription.

"""
function stop_meeting_transcription(
    MeetingId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return chime_sdk_meetings(
        "POST",
        "/meetings/$(MeetingId)/transcription?operation=stop";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function stop_meeting_transcription(
    MeetingId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return chime_sdk_meetings(
        "POST",
        "/meetings/$(MeetingId)/transcription?operation=stop",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
@bsr_exe zgenconf
@bsr_exe zgenterm

@bsr_exe bsr_breit
@bsr_mpi_exe bsr_breit_mpi
@bsr_exe ci_bnk
@bsr_exe bsr_mat
@bsr_mpi_exe bsr_mat_mpi
@bsr_exe bsr_hd

@bsr_exe mult
@bsr_exe bsr_dmat
@bsr_exe bsr_phot
@bsr_mpi_exe bsr_phot_mpi

function dipole(cfg1, cfg2, t1, t2, trans="E1", args...; kwargs...)
    mult(cfg1, cfg2, trans)
    bsr_dmat(cfg1, cfg2, t1, t2, args...; kwargs...)
end

function photo(photo_inp::Matrix, klsp::Integer)
    open("bsr_phot.inp", "w") do file
        write(file, "-------\n")
        for i = 1:size(photo_inp,1)
            if photo_inp[i,1] != ""
                write(file,
                      string(photo_inp[i,1]), " = ",
                      string(photo_inp[i,2]), "\n")
            end
        end
        write(file, "-------\n")
    end
    debug("bsr_phot.inp:\n$(photo_inp)")

    bsr_phot(klsp = klsp)
end

export dipole, photo


text:
using Revise

using ImageViewer, MLDatasets, Plots

x = FashionMNIST.traintensor(1);
plot(image(x); axis = nothing, border = :none)

plot(
    plot(image(x; flip = true); title = "flip = true"),
    plot(image(x; flip = false); title = "flip = false");
    axis = nothing,
    border = :none,
)

x1 = FashionMNIST.traintensor(1);
x2 = CIFAR10.traintensor(2);
plot(
    plot(image(x1)),
    plot(image(x2));
    axis = nothing,
    border = :none
)

x = FashionMNIST.traintensor(1:10);
plot(plot.(image(x, [1,2]))...; axis = nothing, border = :none)


text:
facts("is_square_stochastic with sparse transition") do
    @fact is_square_stochastic(speye(10)) => true
    @fact is_square_stochastic([ speye(10) for _ = 1:3 ]) => true
    @fact is_square_stochastic([ 2*speye(10) for _ = 1:3 ]) => false
end

facts("square stochastic with Int transition") do
    @fact is_square_stochastic(eye(Int, 10)) => true
    @fact is_square_stochastic(2*eye(Int, 10)) => false
end

facts("ismdp with dense transition and vector reward") do
    P = cat(3, eye(10), eye(10))
    @fact ismdp(P, rand(10)) => true
    @fact ismdp(P, rand(11)) => false
end

facts("ismdp with dense transition and array reward") do
    P = cat(3, eye(10), eye(10))
    @fact ismdp(P, rand(10, 10, 2)) => true
    @fact ismdp(P, rand(10, 11, 2)) => false
    @fact ismdp(P, rand(9, 10, 2)) => false
    @fact ismdp(P, rand(10, 10, 3)) => false
end

facts("ismdp with sparse transition and reward") do
    P = [ speye(10) for _ = 1:3 ]
    R = sprand(10, 3, 1/3)
    @fact ismdp(P, R) => true
    P = SparseMatrixCSC[ speye(10), speye(11), speye(10) ]
    R = sprand(10, 3, 1/3)
    @fact ismdp(P, R) => false
    P = [ speye(10) for _ = 1:2 ]
    R = sprand(10, 3, 1/3)
    @fact ismdp(P, R) => false
end

facts("ismdp with sparse transition and dense reward") do
    P = [ speye(10) for _ = 1:3 ]
    @fact ismdp(P, rand(10)) => true
    @fact ismdp(P, rand(11)) => false
    @fact ismdp(P, rand(10, 3)) => true
    @fact ismdp(P, rand(11, 3)) => false
    @fact ismdp(P, rand(11, 4)) => false
    @fact ismdp(P, rand(10, 10, 3)) => true
    @fact ismdp(P, rand(10, 11, 3)) => false
    @fact ismdp(P, rand(11, 10, 3)) => false
    @fact ismdp(P, rand(10, 10, 2)) => false
end

facts("ismdp with dense transition and sparse reward") do
    P = cat(3, eye(10), eye(10), eye(10))
    @fact ismdp(P, sprand(10, 3, 1/3)) => true
    @fact ismdp(P, sprand(11, 3, 1/3)) => false
    @fact ismdp(P, sprand(10, 4, 1/3)) => false
end

facts("ismdp with transition probability and reward types") do
    P = TransitionProbability(cat(3, eye(10), eye(10)))
    R = Reward(rand(10, 2))
    @fact ismdp(P, R) => true
end

facts("ismdp with MDP type") do
    P, R = MDPs.Examples.small()
    @fact ismdp(MDP(P, R)) => true
end


text:
####################
## vine testcases ##
####################

@doc doc"""
Create `Vine` testcases.

1. D-Vine, 6 dimensional

2. 6 dimensional

3. 6 dimensional

4. 7 dimensional

5. 14 dimensional
"""->
function testvine(number::Int)
    if number == 1
        ## dvine
        kk = [0 2 2 2 2 2;
              1 0 3 3 3 3;
              2 2 0 4 4 4;
              3 3 3 0 5 5;
              4 4 4 4 0 6;
              5 5 5 5 5 0]
        vn = Copulas.Vine(kk)
    elseif number == 2
        kk = [0 3 3 3 3 3;
              3 0 3 5 5 5;
              1 2 0 2 2 2;
              5 5 5 0 5 6;
              2 2 2 4 0 4;
              4 4 4 4 4 0]
        vn = Copulas.Vine(kk)
    elseif number == 3
        kk = [0 2 5 4 2 5;
              1 0 3 1 5 3;
              5 2 0 5 2 6;
              1 1 1 0 1 1;
              2 2 2 2 0 2;
              3 3 3 3 3 0]
        vn = Copulas.Vine(kk)
    elseif number == 4
        p1 = Copulas.CTreePaths(1, [2 3 4], [2 3 6 7], [2 3 5])
        p2 = Copulas.CTreePaths(2, [1], [3 6 7], [3 4], [5])
        p3 = Copulas.CTreePaths(3, [2 1], [2 5], [4], [6 7])
        p4 = Copulas.CTreePaths(4, [3 2 1 5 6 7])
        p5 = Copulas.CTreePaths(5, [2 3 1 4], [2 3 1 6 7])
        p6 = Copulas.CTreePaths(6, [3 2 1 5 4], [7])
        p7 = Copulas.CTreePaths(7, [6 3 2 1 5 4])
        vn = Copulas.Vine([p1, p2, p3, p4, p5, p6, p7])
    elseif number == 5
        ## factors 5 and 6 with 2 lags (5 with 1 and 3)
        ## 7, 8, 9 and 10 linked to factor 5
        ## 11, 12, 13 and 14 linked to factor 6
        ## no simulation with 1, 2, 3 and 4 given
        trs = [Copulas.CTreePaths(1, [3,5,6,4,2,14,13,12,11],
                          [3,5,6,7,8,9,10]);
               Copulas.CTreePaths(2, [4,6,5,3,1,7,8,9,10],
                          [4,6,5,14,13,12,11]);
               Copulas.CTreePaths(3, [1], [5,6,4,2,14,13,12,11], [5,6,7],
                          [5,6,8], [5,6,9], [5,6,10]);
               Copulas.CTreePaths(4, [2], [6,5,3,1,7,8,9,10], [6,5,11],
                          [6,5,12], [6,5,13], [6,5,14]);
               Copulas.CTreePaths(5, [3,1], [6,4,2], [6,11], [6,12],
                          [6,13], [6,14], [7], [8], [9], [10]);
               Copulas.CTreePaths(6, [4,2], [5,3,1], [5,7], [5,8], [5,9],
                          [5,10],
                          [11], [12], [13], [14]);
               Copulas.CTreePaths(7, [5,6,3,1,4,2,14,13,12,11],
                          [5,6,3,8,9,10]);
               Copulas.CTreePaths(8, [5,6,3,7,1,4,2,14,13,12,11],
                          [5,6,3,9,10]);
               Copulas.CTreePaths(9, [5,6,3,8,7,1,4,2,14,13,12,11],
                          [5,6,3,10]);
               Copulas.CTreePaths(10, [5,6,3,9,8,7,1,4,2,14,13,12,11]);
               Copulas.CTreePaths(11, [6,5,4,12,13,14,2,3,1,7,8,9,10]);
               Copulas.CTreePaths(12, [6,5,4,11],
                          [6,5,4,13,14,2,3,1,7,8,9,10]);
               Copulas.CTreePaths(13, [6,5,4,12,11],
                          [6,5,4,14,2,3,1,7,8,9,10]);
               Copulas.CTreePaths(14, [6,5,4,2,3,1,7,8,9,10],
                          [6,5,4,13,12,11])]

        vn = Copulas.Vine(trs)
    elseif number == 6
        ## factors 5 and 6 with 2 lags (5 with 1 and 3)
        ## 7, 8, 9 and 10 linked to factor 5
        ## 11, 12, 13 and 14 linked to factor 6
        ## no simulation with 1, 2, 3 and 4 given
        ## faster DVine structure than testvine number 5

        trs = [Copulas.CTreePaths(1, [3,5,6,4,2,14,13,12,11],
                          [3,5,6,7,8,9,10]);
               Copulas.CTreePaths(2, [4,6,5,3,1,7,8,9,10],
                          [4,6,5,14,13,12,11]);
               Copulas.CTreePaths(3, [1], [5,6,4,2,14,13,12,11],
                          [5,6,7,8,9,10]);
               Copulas.CTreePaths(4, [2], [6,5,3,1,7,8,9,10],
                          [6,5,14,13,12,11]);
               Copulas.CTreePaths(5, [3,1], [6,4,2], [6,11], [6,12], [6,13],
                          [6,14],
                          [7], [8], [9], [10]);
               Copulas.CTreePaths(6, [4,2], [5,3,1], [5,7], [5,8], [5,9],
                          [5,10],
                          [11], [12], [13], [14]);
               Copulas.CTreePaths(7, [5,6,3,1,4,2,14,13,12,11], [5,6,8,9,10]);
               Copulas.CTreePaths(8, [5,6,7,3,1,4,2,14,13,12,11], [5,6,9,10]);
               Copulas.CTreePaths(9, [5,6,8,7,3,1,4,2,14,13,12,11], [5,6,10]);
               Copulas.CTreePaths(10, [5,6,9,8,7,3,1,4,2,14,13,12,11]);
               Copulas.CTreePaths(11, [6,5,12,13,14,4,2,3,1,7,8,9,10]);
               Copulas.CTreePaths(12, [6,5,11], [6,5,13,14,4,2,3,1,7,8,9,10]);
               Copulas.CTreePaths(13, [6,5,12,11], [6,5,14,4,2,3,1,7,8,9,10]);
               Copulas.CTreePaths(14, [6,5,4,2,3,1,7,8,9,10], [6,5,13,12,11])]

        vn = Copulas.Vine(trs)
    end
    return vn
end

function testpcc(number::Int)
    if number == 1
        vn = Copulas.setFirstLayer([1 2; 2 3; 3 4])
        vnFinal = Copulas.vineFrom2(vn)
        cops = [Copulas.GaussianPC_Cpp([0.5]),
                Copulas.GaussianPC_Cpp([0.2]),
                Copulas.GaussianPC_Cpp([-0.6]),
                Copulas.ClaytonPC_Cpp([3.5]),
                Copulas.ClaytonPC_Cpp([2.5]),
                Copulas.tPC_Cpp([0.5, 4.2])
                ]
        pcc = Copulas.PCC(cops, vnFinal)

    elseif number == 2
        vn = Copulas.setFirstLayer([1 2; 2 3; 3 4])
        vnFinal = Copulas.vineFrom2(vn)
        pcc = Copulas.PCC(vn)
        pcc[1, 2] = Copulas.GaussianPC_Cpp([0.5])
        pcc[3, 2] = Copulas.GaussianPC_Cpp([0.2])
        pcc[3, 4] = Copulas.ClaytonPC_Cpp([3.2])
    else
        error("No pair copula test case defined for number $number.")
    end
    return pcc
end