text:
mutable struct RTree{ND,BD} <: AbstractRTree
    nodes::Vector{RNode{ND}}
    branches::Vector{Branch{BD}}

    root::RNode{ND}

    RTree(
        nodes::Vector{RNode{ND}},
        branches::Vector{Branch{BD}}
        ) where {ND,BD} =
            new{ND,BD}(nodes, branches)
end

getroot(tree::RTree) = tree.root

"""
    getnodes(tree)

Return an iterable collection with the nodes in a `tree`
"""
getnodes(tree::RTree) = tree.nodes

getnode(tree::RTree, id) = tree.nodes[id]

nnode(tree::RTree) = length(tree.nodes)

ninner(tree::RTree) = mapreduce(isinner, +, getnodes(tree))

nouter(tree::RTree) = mapreduce(isouter, +, getnodes(tree))

Base.show(io::IO, tree::T) where {T<:AbstractTree} =
    print(io, T, ": ", nnode(tree), " nodes")

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

"""
    createtree([ND, BD,]n)

Create a tree with `n` nodes of data type `ND` and branches of data type `BD`.

`ND` and `BD` are `Dict` by default.
"""
function createtree(::Type{ND}, ::Type{BD}, n) where {ND, BD}
    nodes = [RNode{ND}(i) for i in 1:n]
    branches = [Branch{BD}(i) for i in 1:n]
    for c in nodes, br in branches
        c.branch = br
    end

    return RTree(nodes, branches)
end
createtree(n::Int) = createtree(Dict, Dict, n)

"""
    addnode!(tree, n=1)

Add `n` nodes to a `tree`. Return the first new node.
"""
function addnode!(tree::RTree{ND,BD}, n=1) where {ND, BD}
    i = nnode(tree)
    push!(tree.nodes, (RNode{ND}(j) for j in (i+1):(i+n))...)
    push!(tree.branches, (Branch{BD}(j) for j in (i+1):(i+n))...)
    for c in tree.nodes[i+1:i+n], br in tree.branches[i+1:i+n]
        c.branch = br
    end

    return getnode(tree, i+1)
end

max_stars_repo_name: eascarrunz/RootedTrees

max_issues_repo_name: eascarrunz/RootedTrees

max_forks_repo_name: eascarrunz/RootedTrees


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "libx11"
version = v"1.6.9"

# Collection of sources required to build libx11
sources = [
    "https://xorg.freedesktop.org/releases/individual/util/util-macros-1.19.2.tar.bz2" =>
    "d7e43376ad220411499a79735020f9d145fdc159284867e99467e0d771f3e712",

    "https://xorg.freedesktop.org/archive/individual/proto/xorgproto-2019.2.tar.bz2" =>
    "46ecd0156c561d41e8aa87ce79340910cdf38373b759e737fcbba5df508e7b8e",

    "https://xcb.freedesktop.org/dist/xcb-proto-1.13.tar.bz2" =>
    "7b98721e669be80284e9bbfeab02d2d0d54cd11172b72271e47a2fe875e2bde1",

    "https://xorg.freedesktop.org/releases/individual/lib/libXdmcp-1.1.3.tar.bz2" =>
    "20523b44aaa513e17c009e873ad7bbc301507a3224c232610ce2e099011c6529",

    "https://xorg.freedesktop.org/releases/individual/lib/libXau-1.0.9.tar.bz2" =>
    "ccf8cbf0dbf676faa2ea0a6d64bcc3b6746064722b606c8c52917ed00dcb73ec",

    "https://xcb.freedesktop.org/dist/libxcb-1.13.1.tar.bz2" =>
    "a89fb7af7a11f43d2ce84a844a4b38df688c092bf4b67683aef179cdf2a647c4",

    "https://xcb.freedesktop.org/dist/xcb-util-0.4.0.tar.bz2" =>
    "46e49469cb3b594af1d33176cd7565def2be3fa8be4371d62271fabb5eae50e9",

    "https://xcb.freedesktop.org/dist/xcb-util-wm-0.4.1.tar.bz2" =>
    "28bf8179640eaa89276d2b0f1ce4285103d136be6c98262b6151aaee1d3c2a3f",

    "https://xcb.freedesktop.org/dist/xcb-util-image-0.4.0.tar.bz2" =>
    "2db96a37d78831d643538dd1b595d7d712e04bdccf8896a5e18ce0f398ea2ffc",

    "https://xcb.freedesktop.org/dist/xcb-util-keysyms-0.4.0.tar.bz2" =>
    "0ef8490ff1dede52b7de533158547f8b454b241aa3e4dcca369507f66f216dd9",

    "https://xcb.freedesktop.org/dist/xcb-util-renderutil-0.3.9.tar.bz2" =>
    "c6e97e48fb1286d6394dddb1c1732f00227c70bd1bedb7d1acabefdd340bea5b",

    "https://xorg.freedesktop.org/releases/individual/lib/libX11-1.6.9.tar.bz2" =>
    "9cc7e8d000d6193fa5af580d50d689380b8287052270f5bb26a5fb6b58b2bed1",

    "https://xorg.freedesktop.org/releases/individual/lib/libXext-1.3.4.tar.bz2" =>
    "59ad6fcce98deaecc14d39a672cf218ca37aba617c9a0f691cac3bcd28edf82b",

    "https://xorg.freedesktop.org//releases/individual/lib/xtrans-1.4.0.tar.bz2" =>
    "377c4491593c417946efcd2c7600d1e62639f7a8bbca391887e2c4679807d773",

    "https://xkbcommon.org/download/libxkbcommon-0.9.1.tar.xz" =>
    "d4c6aabf0a5c1fc616f8a6a65c8a818c03773b9a87da9fbc434da5acd1199be0"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd util-macros-1.19.2/
./configure --prefix=$prefix --host=$target
make install

cd ../xorgproto-2019.2/
./configure --prefix=$prefix --host=$target
make install
export PKG_CONFIG_PATH=$PKG_CONFIG_PATH:$prefix/share/pkgconfig

cd ../xcb-proto-1.13/
./configure --prefix=$prefix --host=$target
make install
cd ../libXdmcp-1.1.3/
./configure --prefix=$prefix --host=$target --disable-static
make
make install

cd ../libXau-1.0.9/
./configure --prefix=$prefix --host=$target
make
make install

cd ../libxcb-1.13.1/
curl https://git.archlinux.org/svntogit/packages.git/plain/trunk/libxcb-1.1-no-pthread-stubs.patch?h=packages/libxcb --output libxcb-1.1-no-pthread-stubs.patch
patch -Np1 -i libxcb-1.1-no-pthread-stubs.patch
autoreconf -I$prefix/share/aclocal -vfi
./configure --prefix=$prefix --host=$target --disable-static
make
make install

cd ../xcb-util-0.4.0
./configure --prefix=$prefix --host=$target --disable-static
make
make install

cd ../xcb-util-wm-0.4.1
./configure --prefix=$prefix --host=$target --disable-static
make
make install

cd ../xcb-util-image-0.4.0
./configure --prefix=$prefix --host=$target --disable-static
make
make install

cd ../xcb-util-keysyms-0.4.0
./configure --prefix=$prefix --host=$target --disable-static
make
make install

cd ../xcb-util-renderutil-0.3.9
./configure --prefix=$prefix --host=$target --disable-static
make
make install

cd ../xtrans-1.4.0/
ls
./configure --prefix=$prefix --host=$target
make
make install

cd ../libX11-1.6.9/
./configure --prefix=$prefix --host=$target --disable-static
make
make install

cd ../libXext-1.3.4/
./configure --prefix=$prefix --host=$target --disable-static
make
make install

cd ../libxkbcommon-0.9.1
apk add meson
meson setup build -Dprefix=/workspace/destdir -Denable-wayland=false -Denable-docs=false
ninja -C build
ninja -C build install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libX11", :libx11),
    LibraryProduct(prefix, "libxcb-glx", Symbol("libxcb-glx")),
    LibraryProduct(prefix, "libxcb-damage", Symbol("libxcb-damage")),
    LibraryProduct(prefix, "libxcb-dri2", Symbol("libxcb-dri2")),
    LibraryProduct(prefix, "libxcb-xkb", Symbol("libxcb-xkb")),
    LibraryProduct(prefix, "libxcb", :libxcb),
    LibraryProduct(prefix, "libxcb-dpms", Symbol("libxcb-dpms")),
    LibraryProduct(prefix, "libXau", :libxau),
    LibraryProduct(prefix, "libxcb-sync", Symbol("libxcb-sync")),
    LibraryProduct(prefix, "libxcb-record", Symbol("libxcb-record")),
    LibraryProduct(prefix, "libxcb-present", Symbol("libxcb-present")),
    LibraryProduct(prefix, "libxcb-shape", Symbol("libxcb-shape")),
    LibraryProduct(prefix, "libxcb-res", Symbol("libxcb-res")),
    LibraryProduct(prefix, "libxcb-xf86dri", Symbol("libxcb-xf86dri")),
    LibraryProduct(prefix, "libxcb-xtest", Symbol("libxcb-xtest")),
    LibraryProduct(prefix, "libxcb-dri3", Symbol("libxcb-dri3")),
    LibraryProduct(prefix, "libxcb-shm", Symbol("libxcb-shm")),
    LibraryProduct(prefix, "libxcb-screensaver", Symbol("libxcb-screensaver")),
    LibraryProduct(prefix, "libxcb-xfixes", Symbol("libxcb-xfixes")),
    LibraryProduct(prefix, "libxcb-composite", Symbol("libxcb-composite")),
    LibraryProduct(prefix, "libxcb-xinput", Symbol("libxcb-xinput")),
    LibraryProduct(prefix, "libxcb-icccm", Symbol("libxcb-icccm")),
    LibraryProduct(prefix, "libXdmcp", :libxdcmp),
    LibraryProduct(prefix, "libxcb-xvmc", Symbol("libxcb-xvmc")),
    LibraryProduct(prefix, "libXext", :libxext),
    LibraryProduct(prefix, "libxcb-xinerama", Symbol("libxcb-xinerama")),
    LibraryProduct(prefix, "libxcb-randr", Symbol("libxcb-randr")),
    LibraryProduct(prefix, "libxcb-xv", Symbol("libxcb-xv")),
    LibraryProduct(prefix, "libxcb-render", Symbol("libxcb-render")),
    LibraryProduct(prefix, "libX11-xcb", Symbol("libx11-xcb")),
    LibraryProduct(prefix, "libxkbcommon", Symbol("libxkbcommon")),
    LibraryProduct(prefix, "libxkbcommon-x11", Symbol("libxkbcommon-x11"))
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

ext: jl


text:
using MinimaxAI
using Test

# playBones()

playBones()

MinimaxAI.playConnect4()

# @testset "MinimaxAI.jl" begin
#     @test my_f(2,2) == 10
# end

# @testset "MinimaxAI.jl" begin
#     @test ... == ...
# end


text:
function gamma_posterior(chain)
    gamma_rank = denserank(chain.gamma)
    gamma_unique = sort(unique(chain.gamma))
    gamma_card = length(unique(gamma_unique))
    gamma_freq = counts(gamma_rank, gamma_card)
    pγ = gamma_freq ./ length(chain.gamma)
    N = length(gamma_unique)
    G = length(gamma_unique[1])
    df = DataFrame(
        [zeros(Bool, N) for k in 2:G],
        [Symbol("group$i") for i in 2:G]
    )
    for i in 1:N, j in 2:G
        df[i, j - 1] = gamma_unique[i][j]
    end
    df[!, :prob] = pγ
    return df
end


text:
cd(@__DIR__); include("setups/simpletree.jl")
gr(dpi = 200)

## build SGWT frame
SGWT = sgwt_frame(Matrix(W); nf = 6)
SGWT = reshape(SGWT, (N, :))

important_idx = sortperm((SGWT' * f).^2; rev = true)
plot(layout = Plots.grid(3, 7), size = (1400, 1500))
for i = 1:21
    ind = important_idx[i]
    j = Int(floor(ind / N))
    x = ind - j * N
    w = SGWT[:, ind]
    w ./= norm(w, 2)
    scatter!(X[:, 1], X[:, 2], marker_z = w, ms = 5, c = :viridis, subplot = i, mswidth = 0)
    plot!(frame = :none, cbar = false, grid = false, legend = false, subplot = i)
    plot!(title = latexstring("\\psi^{\\mathrm{SGWT}}_{", j, ",", x, "}"), titlefont = 18, subplot = i)
end
plt = current()
savefig(plt, "../figs/simpletree_f_sgwt_top21.png")


text:
function temperature()
    temp = 0
    num = 0
    side = Ost
    flag = false
    while !flag
        while !isborder(r, side)
            if ismarker(r)
                temp += temperature(r)
                num += 1
            end
            move!(r, side)
        end
        side = inverse_side(side)
        flag = isborder(r, Nord)
        if !flag
            move!(r, Nord)
        end
    end
    if num != 0
        println(temp / num)
    else
        println(0)
    end
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
# Automatically generated using Clang.jl

# const CSFML_NETWORK_API = CSFML_API_IMPORT

@cenum sfFtpTransferMode::UInt32 begin
    sfFtpBinary = 0
    sfFtpAscii = 1
    sfFtpEbcdic = 2
end

@cenum sfFtpStatus::UInt32 begin
    sfFtpRestartMarkerReply = 110
    sfFtpServiceReadySoon = 120
    sfFtpDataConnectionAlreadyOpened = 125
    sfFtpOpeningDataConnection = 150
    sfFtpOk = 200
    sfFtpPointlessCommand = 202
    sfFtpSystemStatus = 211
    sfFtpDirectoryStatus = 212
    sfFtpFileStatus = 213
    sfFtpHelpMessage = 214
    sfFtpSystemType = 215
    sfFtpServiceReady = 220
    sfFtpClosingConnection = 221
    sfFtpDataConnectionOpened = 225
    sfFtpClosingDataConnection = 226
    sfFtpEnteringPassiveMode = 227
    sfFtpLoggedIn = 230
    sfFtpFileActionOk = 250
    sfFtpDirectoryOk = 257
    sfFtpNeedPassword = 331
    sfFtpNeedAccountToLogIn = 332
    sfFtpNeedInformation = 350
    sfFtpServiceUnavailable = 421
    sfFtpDataConnectionUnavailable = 425
    sfFtpTransferAborted = 426
    sfFtpFileActionAborted = 450
    sfFtpLocalError = 451
    sfFtpInsufficientStorageSpace = 452
    sfFtpCommandUnknown = 500
    sfFtpParametersUnknown = 501
    sfFtpCommandNotImplemented = 502
    sfFtpBadCommandSequence = 503
    sfFtpParameterNotImplemented = 504
    sfFtpNotLoggedIn = 530
    sfFtpNeedAccountToStore = 532
    sfFtpFileUnavailable = 550
    sfFtpPageTypeUnknown = 551
    sfFtpNotEnoughMemory = 552
    sfFtpFilenameNotAllowed = 553
    sfFtpInvalidResponse = 1000
    sfFtpConnectionFailed = 1001
    sfFtpConnectionClosed = 1002
    sfFtpInvalidFile = 1003
end

@cenum sfHttpMethod::UInt32 begin
    sfHttpGet = 0
    sfHttpPost = 1
    sfHttpHead = 2
    sfHttpPut = 3
    sfHttpDelete = 4
end

@cenum sfHttpStatus::UInt32 begin
    sfHttpOk = 200
    sfHttpCreated = 201
    sfHttpAccepted = 202
    sfHttpNoContent = 204
    sfHttpResetContent = 205
    sfHttpPartialContent = 206
    sfHttpMultipleChoices = 300
    sfHttpMovedPermanently = 301
    sfHttpMovedTemporarily = 302
    sfHttpNotModified = 304
    sfHttpBadRequest = 400
    sfHttpUnauthorized = 401
    sfHttpForbidden = 403
    sfHttpNotFound = 404
    sfHttpRangeNotSatisfiable = 407
    sfHttpInternalServerError = 500
    sfHttpNotImplemented = 501
    sfHttpBadGateway = 502
    sfHttpServiceNotAvailable = 503
    sfHttpGatewayTimeout = 504
    sfHttpVersionNotSupported = 505
    sfHttpInvalidResponse = 1000
    sfHttpConnectionFailed = 1001
end

struct sfIpAddress
    address::NTuple{16, UInt8}
end

@cenum sfSocketStatus::UInt32 begin
    sfSocketDone = 0
    sfSocketNotReady = 1
    sfSocketPartial = 2
    sfSocketDisconnected = 3
    sfSocketError = 4
end

const sfFtpDirectoryResponse = Cvoid
const sfFtpListingResponse = Cvoid
const sfFtpResponse = Cvoid
const sfFtp = Cvoid
const sfHttpRequest = Cvoid
const sfHttpResponse = Cvoid
const sfHttp = Cvoid
const sfPacket = Cvoid
const sfSocketSelector = Cvoid
const sfTcpListener = Cvoid
const sfTcpSocket = Cvoid
const sfUdpSocket = Cvoid


text:
using JLD
using ProgressMeter

export window_bed_file, streaming_cor, cov2cor!, streaming_edgeimpact

# creates a BitArray with a 1 for every bin overlapped by the bed file regions
function window_bed_file(stream, contigs; binSize=1000)
    numBins = ceil(Int64, sum(contigs.sizes) / binSize)
    chrOffsets = Dict{ASCIIString,Int64}()
    for i in 1:contigs.count
        chrOffsets[contigs.names[i]] = contigs.offsets[i]
    end

    # mark all bins that are touched with 1
    binValues = falses(numBins)
    for line in eachline(stream)
        parts = split(line, '\t')
        if haskey(chrOffsets, parts[1])
            startPos = ceil(Int64, (chrOffsets[parts[1]]+parse(Int64, parts[2]))/binSize)
            endPos = ceil(Int64, (chrOffsets[parts[1]]+parse(Int64, parts[3]))/binSize)
            for i in startPos:endPos
                binValues[i] = true
            end
        end
    end

    binValues
end

"Designed to compute the correlation matrix from a ChromNet data bundle."
function streaming_cor(dataBundles; chunkSize=10000, quiet=false)

    # load data arrays
    P = 0
    datas = Any[]
    N = read(dataBundles[1], "N")
    for bundle in dataBundles
        push!(datas, readmmap(bundle["data"]))
        P += read(bundle, "P")
        @assert N == read(bundle, "N") "Data bundles have differing numbers of samples!"
    end

    # get the means
    means = zeros(P)
    lastPos = 0
    for bundle in dataBundles
        bundleP = read(bundle, "P")
        means[lastPos+1:lastPos+bundleP] = read(bundle, "means")
        lastPos = lastPos+bundleP
    end

    # loop through all the chunks
    XtX = zeros(Float64, P, P)
    chunk = Array(Float64, P, chunkSize)
    numChunks = round(Int64, ceil(N/chunkSize))
    p = Progress(N, 0.1, "Computing correlation: ", 34, :green, STDERR)
    for i in 1:chunkSize:N
        if i+chunkSize-1 > N # account for the last unevenly sized chunk
            chunk = Array(Float64, P, N-i+1)
        end
        streaming_cor_chunk!(XtX, datas, (i, min(i+chunkSize-1,N)), chunk, means)
        if !quiet update!(p, min(i+chunkSize-1,N)) end
    end

    cov2cor!(XtX)
    XtX + XtX' - eye(P)
end
function streaming_cor_chunk!(XtX, datas, range, chunk, means)
    w = size(chunk)[2]

    lastPos = 0
    for data in datas
        bundleP = size(data)[1]
        chunk[lastPos+1:lastPos+bundleP,:] = data[:,range[1]:range[1]+w-1]
        lastPos = lastPos+bundleP
    end

    for j in 1:size(chunk)[2]
        @simd for i in 1:size(chunk)[1]
            @inbounds chunk[i,j] -= means[i]
        end
    end
    BLAS.syrk!('U', 'N', 1.0, chunk, 1.0, XtX)
end

function streaming_edgeimpact(dataBundles, C, rows, cols; chunkSize=10000, quiet=false)
    IC = inv(C)

    # load data arrays
    P = 0
    datas = Any[]
    N = read(dataBundles[1], "N")
    for bundle in dataBundles
        push!(datas, readmmap(bundle["data"]))
        P += read(bundle, "P")
        @assert N == read(bundle, "N") "Data bundles have differing numbers of samples!"
    end

    # get the means and norms
    means = zeros(P)
    norms = zeros(P)
    lastPos = 0
    for bundle in dataBundles
        bundleP = read(bundle, "P")
        means[lastPos+1:lastPos+bundleP] = read(bundle, "means")
        norms[lastPos+1:lastPos+bundleP] = read(bundle, "norms")
        lastPos = lastPos+bundleP
    end

    # loop through all the chunks
    out = Float64[]
    vtmp = zeros(Float64, P)
    ctmp = zeros(Float64, P)
    outChunk = zeros(Float64, chunkSize)
    chunk = Array(Float64, P, chunkSize)
    chunkTmp = Array(Float64, P, chunkSize)
    p = Progress(N, 0.1, "Computing sample relevance: ", 34, :green, STDERR)
    for i in 1:chunkSize:N
        if i+chunkSize-1 > N # account for the last unevenly sized chunk
            outChunk = Array(Float64, N-i+1)
            chunk = Array(Float64, P, N-i+1)
            chunkTmp = Array(Float64, P, N-i+1)
        end
        streaming_edgeimpact_chunk!(outChunk, datas, (i, min(i+chunkSize-1,N)), IC, rows, cols, chunk, chunkTmp, means, norms, vtmp, ctmp)
        append!(out, outChunk)
        if !quiet update!(p, min(i+chunkSize-1,N)) end
    end
    out
end
function streaming_edgeimpact_chunk!(out, datas, range, IC, rows, cols, chunk, chunkTmp, means, norms, vtmp, ctmp)
    P,w = size(chunk)
    @assert size(IC) == (P,P)
    @assert length(means) == P
    @assert length(norms) == P
    @assert length(vtmp) == P
    @assert length(ctmp) == P
    @assert length(out) == w

    # load the data chunk and normalize it
    lastPos = 0
    for data in datas
        bundleP = size(data)[1]
        chunk[lastPos+1:lastPos+bundleP,:] = data[:,range[1]:range[1]+w-1] # converting to a float array is also important to get BLAS speed
        lastPos = lastPos+bundleP
    end
    for j in 1:size(chunk)[2]
        @simd for i in 1:P
            @inbounds chunk[i,j] -= means[i]
            @inbounds chunk[i,j] /= norms[i]
        end
    end

    # precompute the matrix-vector multiplication for all samples
    BLAS.symm!('L', 'U', 1.0, IC, chunk, 0.0, chunkTmp)

    # compute the sum of the changes in the given edges for each sample
    for ind in 1:length(out)
        copy!(ctmp, 1, chunk, P*(ind-1)+1, P)
        copy!(vtmp, 1, chunkTmp, P*(ind-1)+1, P)
        s = 1 - dot(ctmp,vtmp)
        totalChange = 0.0
        for j in cols
            @inbounds dj = sqrt(1 - ctmp[j].^2)
            for i in rows
                @inbounds di = sqrt(1 - ctmp[i].^2)
                @inbounds totalChange += di*dj*(IC[i,j] + vtmp[i]*vtmp[j]/s) - IC[i,j]
            end
        end
        @inbounds out[ind] = totalChange
    end
    out
end

# "Designed to compute the correlation matrix from a ChromNet data bundle."
# function streaming_cor(dataFile, memData=nothing; chunkSize=10000, quiet=false)
#
#     # open data file
#     f = jldopen(dataFile)
#     bundleP = read(f, "P")
#     N = read(f, "N")
#     means = read(f, "means")
#     data = readmmap(f["data"])
#     means = [read(f, "means"); mean(memData, 2)]
#     close(f)
#
#     if memData == nothing
#         memData = zeros(0,N)
#     else
#         @assert size(memData)[2] == N "Passed in-memory data size does not match bundle ($(size(memData)[2]) != $N)"
#     end
#     P = size(memData)[1] + bundleP
#
#     # loop through all the chunks
#     XtX = zeros(Float64, P, P)
#     chunk = Array(Float64, P, chunkSize)
#     numChunks = round(Int64, ceil(N/chunkSize))
#     p = Progress(N, 0.1, "Computing correlation: ", 34, :green, STDERR)
#     for i in 1:chunkSize:N
#         if i+chunkSize-1 > N # account for the last unevenly sized chunk
#             chunk = Array(Float64, P, N-i+1)
#         end
#         streaming_cor_chunk!(XtX, data, memData, (i, min(i+chunkSize-1,N)), chunk, means)
#         if !quiet update!(p, min(i+chunkSize-1,N)) end
#     end
#
#     cov2cor!(XtX)
#     XtX + XtX' - eye(P)
# end
# function streaming_cor_chunk!(XtX, data, memData, range, chunk, means)
#     memP = size(memData)[1]
#     bundleP = size(chunk)[1] - memP
#     w = size(chunk)[2]
#
#     chunk[1:bundleP,:] = data[:,range[1]:range[1]+w-1]
#     chunk[bundleP+1:end,:] = memData[:,range[1]:range[1]+w-1]
#     #copy!(chunk, 1, data, P*(range[1]-1)+1, P*w) # converting to a float array is also important to get BLAS speed
#     for j in 1:size(chunk)[2]
#         @simd for i in 1:size(chunk)[1]
#             @inbounds chunk[i,j] -= means[i]
#         end
#     end
#     BLAS.syrk!('U', 'N', 1.0, chunk, 1.0, XtX)
# end

function cov2cor!(M)
    for i in 1:size(M)[1]
        val = sqrt(M[i,i])
        if val > 0.0
            M[i,:] ./= val
            M[:,i] ./= val
        end
    end
    M
end


text:
export isunivariate

###############################################################################
# other stuff, trivia and non-trivia
###############################################################################

function Nemo.PolynomialRing(R::Nemo.Ring, n::Int, s::String="x";
  cached::Bool = false, ordering::Symbol = :lex)
  return Nemo.PolynomialRing(R, ["$s$i" for i=1:n], cached = cached,
                         ordering = ordering)
end

#TODO: only makes sense if f is univ (uses only one var)
function (Rx::FmpzPolyRing)(f::fmpq_mpoly)
  fp = Rx()
  R = base_ring(Rx)
  d = denominator(f)
  @assert d == 1
  for t = terms(f)
    e = total_degree(t)
    @assert length(t) == 1
    c = coeff(t, 1)
    setcoeff!(fp, e, numerator(c*d))
  end
  return fp
end

function (Rx::FmpqPolyRing)(f::fmpq_mpoly)
  fp = Rx()
  R = base_ring(Rx)
  for t = terms(f)
    e = total_degree(t)
    @assert length(t) == 1
    c = coeff(t, 1)
    setcoeff!(fp, e, c)
  end
  return fp
end

function (Rx::GFPPolyRing)(f::fmpq_mpoly)
  fp = Rx()
  R = base_ring(Rx)
  d = denominator(f)
  for t = terms(f)
    e = total_degree(t)
    @assert length(t) == 1
    c = coeff(t, 1)
    setcoeff!(fp, e, R(numerator(c*d)))
  end
  return fp * inv(R(d))
end

function Hecke.lead(f::AbstractAlgebra.MPolyElem)
  iszero(f) && error("zero poly")
  return coeff(f, 1)
end

function Hecke.coefficients(f::AbstractAlgebra.MPolyElem)
  return Hecke.PolyCoeffs(f)
end

function Base.iterate(PC::Hecke.PolyCoeffs{<:AbstractAlgebra.MPolyElem}, st::Int = 0)
  st += 1
  if st > length(PC.f)
    return nothing
  else
    return coeff(PC.f, st), st
  end
end

Base.length(M::Hecke.PolyCoeffs{<:AbstractAlgebra.MPolyElem}) = length(M.f)

function mul!(res::fmpq_mpoly, a::fmpq_mpoly, c::fmpz)
  ccall((:fmpq_mpoly_scalar_mul_fmpz, libflint), Nothing,
    (Ref{fmpq_mpoly}, Ref{fmpq_mpoly}, Ref{fmpz}, Ref{FmpqMPolyRing}), res, a, c, parent(a))
  return nothing
end

@doc Markdown.doc"""
    isunivariate(f::Generic.MPoly{T}) where T <: NumFieldElem -> Bool, PolyElem{T}

Tests if $f$ involves only one variable. If so, return a corresponding univariate polynomial.
"""
function isunivariate(f::Generic.MPoly{T}) where T
  kx, x = PolynomialRing(base_ring(f), "x", cached = false)
  if ngens(parent(f)) == 1
    f1 = kx()
    for i = 1:f.length
      setcoeff!(f1, Int(f.exps[1, i]), f.coeffs[i])
    end
    return true, f1
  end
  if f.length == 0
    @assert iszero(f)
    return true, kx(0)
  end
  n = ngens(parent(f))
  i = 1
  while i <= n && iszero(f.exps[i, :])
    i += 1
  end
  j = n
  while j >= 1 && iszero(f.exps[j, :])
    j -= 1
  end
  if i != j
    return false, x
  end
  f1 = kx()
  for j = 1:f.length
    setcoeff!(f1, Int(f.exps[i, j]), f.coeffs[j])
  end
  return true, f1
end


text:
@testset "Exceptions" begin
    # Test exception messages
    tz = VariableTimeZone(
        "Imaginary/Zone",
        [Transition(DateTime(1800,1,1), FixedTimeZone("IST",0,0))],
        DateTime(1980,1,1),
    )

    @test sprint(showerror, AmbiguousTimeError(DateTime(2015,1,1), tz)) ==
        "AmbiguousTimeError: Local DateTime 2015-01-01T00:00:00 is ambiguous within Imaginary/Zone"
    @test sprint(showerror, NonExistentTimeError(DateTime(2015,1,1), tz)) ==
        "NonExistentTimeError: Local DateTime 2015-01-01T00:00:00 does not exist within Imaginary/Zone"
    @test sprint(showerror, UnhandledTimeError(tz)) ==
        "UnhandledTimeError: TimeZone Imaginary/Zone does not handle dates on or after 1980-01-01T00:00:00 UTC"
end


text:
using Literate
using Plots
using Graphs, MetaGraphs
using GraphRecipes
using NetworkLayout
using Logging
using JSON3
using PowerSystems
using PowerSystemCaseBuilder
using PowerSimulations
using InteractiveUtils
using Base64
using REPL
using Markdown
using JSONTables
import Literate

export D3PowerSystemTypes
export D3NetworkLayout

disable_logging() = configure_logging(console_level = Logging.AboveMaxLevel, filename = nothing)

function lx_textinput(p)
    rpath = Xranklin.unixify(strip(p[1]))
    lc = Xranklin.cur_lc()
    Xranklin.attach(lc, rpath)
    Xranklin.html(open(f->read(f, String), rpath), lc)
end

#################################################################

struct D3LoadProfile
    names
    data
end

function render(loadprofile::D3LoadProfile)
    t = read(joinpath(@__DIR__, "./templates/d3loadprofile.html"), String)
    r = replace(t, "{{data}}" => JSON3.write(Dict([n => JSON3.read(arraytable(d)) for (n, d) in zip(loadprofile.names, loadprofile.data)]...)))
    return r
end

Base.showable(::MIME"juliavscode/html", ::D3LoadProfile) = true
Base.show(io::IO, ::MIME"juliavscode/html", t::D3LoadProfile) = print(io, render(t))

#################################################################

struct D3PowerSystemTypes
    datatype::DataType
end

function get_field_doc(t, x)
    try
        REPL.fielddoc(t, x)
    catch
        try
            REPL.fielddoc(Base.Docs.Binding(PowerSystems, Symbol(t)), x)
        catch
            Markdown.MD()
        end
    end
end
function children!(data, t)
    # @show t
    data["name"] = string(t)

    if !isabstracttype(t)
        data["fieldnames"] = [x for x in string.(fieldnames(t))]
        data["fieldtypes"] = [x for x in string.(fieldtypes(t))]
        data["fielddocs"] = [get_field_doc(t, x) for x in fieldnames(t)]
        data["fielddocs"] = map(data["fielddocs"]) do d
            if startswith(string(d), "`$(string(t))`")
                Markdown.MD()
            else
                Markdown.MD(d)
            end
        end
        data["fielddocs"] = [strip(Markdown.string(x)) for x in data["fielddocs"]]
        data["value"] = length(fieldnames(t))
    else
        data["children"] = []
        for subtype in subtypes(t)
            c = Dict()
            push!(data["children"], c)
            children!(c, subtype)
        end
    end
end

function render(types::D3PowerSystemTypes, id = "graph")
    t = read(joinpath(@__DIR__, "./templates/d3powersystemtypes.html"), String)
    data = Dict()
    children!(data, types.datatype)
    r = replace(t, "{{data}}" => JSON3.write(data), "{{id}}" => id)
    return r
end

function hfun_d3types_component()
    return render(D3PowerSystemTypes(PowerSystems.Component), "component")
end

function hfun_d3types_timeseries()
    return render(D3PowerSystemTypes(PowerSystems.TimeSeriesData), "timeseries")
end

Base.showable(::MIME"juliavscode/html", ::D3PowerSystemTypes) = true
Base.show(io::IO, ::MIME"juliavscode/html", t::D3PowerSystemTypes) = print(io, render(t))

#################################################################

function create_plot(g)
    n = first(size(g))
    graphplot(
        g,
        curves = false,
        names = rpad.(lpad.(string.(1:n), 2), 2),
        nodesize = 0.2,
        nodeshape = :circle,
        method = :sfdp,
        # layout_kw=Dict(:tol=>0.000001,  :iterations=>10000)
    )
end

function create_graph(sys::PowerSystems.System, res::Union{Nothing,ProblemResults} = nothing)
    buses = collect(get_components(Bus, sys))
    n_bus = length(buses)
    gens = collect(get_components(Generator, sys))
    n_gen = length(gens)
    loads = collect(get_components(PowerLoad, sys))
    n_load = length(loads)
    graph = MetaGraph(SimpleGraph(n_bus + n_gen + n_load))
    from_to_buses = map(get_components(Branch, sys)) do branch
        (
            only(findall(x -> x.name === get_arc(branch).from.name, buses)),
            only(findall(x -> x.name === get_arc(branch).to.name, buses))
        )
    end
    for (f, t) in from_to_buses
        add_edge!(graph, f, t)
    end
    for (i, bus) in enumerate(buses)
        set_prop!(graph, i, :group, get_bustype(bus).value)
        set_prop!(graph, i, :name, get_name(bus))
    end
    for (i, gen) in enumerate(gens)
        gen_node = i + n_bus
        add_edge!(graph, gen_node, only(findall(x -> x.name === get_bus(gen).name, buses)))
        set_prop!(graph, gen_node, :group, 0)
        set_prop!(graph, gen_node, :name, get_name(gen))
    end
    for (i, load) in enumerate(loads)
        load_node = i + n_bus + n_gen
        add_edge!(graph, load_node, only(findall(x -> x.name === get_bus(load).name, buses)))
        set_prop!(graph, load_node, :group, 6)
        set_prop!(graph, load_node, :name, get_name(load))
        set_prop!(graph, load_node, :timeseries, arraytable(get_time_series_array(Deterministic, load, "max_active_power")))
    end
    graph
end
struct D3NetworkLayout
    graph
    marked
end

D3NetworkLayout(graph) = D3NetworkLayout(graph, [])

function render(l::D3NetworkLayout)
    t = read(joinpath(@__DIR__, "./templates/d3networklayout.html"), String)
    data = Dict("nodes" => [], "links" => [])
    graph = l.graph
    for node in vertices(graph)
        d = Dict(
            "id" => node,
            "name" => get_prop(graph, node, :name),
            "group" => get_prop(graph, node, :group),
        )
        if d["group"] == 6
            d["timeseries"] = JSON3.read(get_prop(graph, node, :timeseries))
        end
        push!(data["nodes"], d)
    end
    for (i, edge) in enumerate(edges(graph))
        push!(data["links"], Dict("id" => i, "source" => src(edge), "target" => dst(edge)))
    end
    r = replace(t, "{{time_series_data}}" => JSON3.write([]), "{{graph}}" => JSON3.write(data), "{{marked}}" => JSON3.write(l.marked))
    return r
end

Base.showable(::MIME"juliavscode/html", ::D3NetworkLayout) = true
Base.show(io::IO, ::MIME"juliavscode/html", l::D3NetworkLayout) = print(io, render(l))

function hfun_d3networklayout()
    system = build_system(PSITestSystems, "modified_RTS_GMLC_DA_sys");
    graph = create_graph(system)
    return render(D3NetworkLayout(graph))
end

max_stars_repo_name: NREL-SIIP/SIIP-Tutorial

max_issues_repo_name: NREL-SIIP/SIIP-Tutorial

max_forks_repo_name: NREL-SIIP/SIIP-Tutorial


text:
abstract NUTSState{F<:VariateForm} <: HMCSamplerState{F}

type MuvNUTSState <: NUTSState{Multivariate}
  pstateplus::ParameterState{Continuous, Multivariate}
  pstateminus::ParameterState{Continuous, Multivariate}
  pstateprime::ParameterState{Continuous, Multivariate}
  pstatedprime::ParameterState{Continuous, Multivariate}
  tune::MCTunerState
  ratio::Real
  a::Real
  momentum::RealVector
  momentumplus::RealVector
  momentumminus::RealVector
  momentumprime::RealVector
  oldhamiltonian::Real
  newhamiltonian::Real
  u::Real
  v::Integer
  j::Integer
  n::Integer
  nprime::Integer
  ndprime::Integer
  s::Bool
  sprime::Bool
  sdprime::Bool
  update::Bool
  count::Integer

  function MuvNUTSState(
    pstateplus::ParameterState{Continuous, Multivariate},
    pstateminus::ParameterState{Continuous, Multivariate},
    pstateprime::ParameterState{Continuous, Multivariate},
    pstatedprime::ParameterState{Continuous, Multivariate},
    tune::MCTunerState,
    ratio::Real,
    a::Real,
    momentum::RealVector,
    momentumplus::RealVector,
    momentumminus::RealVector,
    momentumprime::RealVector,
    oldhamiltonian::Real,
    newhamiltonian::Real,
    u::Real,
    v::Integer,
    j::Integer,
    n::Integer,
    nprime::Integer,
    ndprime::Integer,
    s::Bool,
    sprime::Bool,
    sdprime::Bool,
    update::Bool,
    count::Integer
  )
    if !isnan(ratio)
      @assert ratio > 0 "Acceptance ratio should be positive"
    end
    if !isnan(a)
      @assert 0 <= a <= 1 "Acceptance probability should be in [0, 1]"
    end
    @assert count >= 0 "Number of iterations (count) should be non-negative"
    new(
      pstateplus,
      pstateminus,
      pstateprime,
      pstatedprime,
      tune,
      ratio,
      a,
      momentum,
      momentumplus,
      momentumminus,
      momentumprime,
      oldhamiltonian,
      newhamiltonian,
      u,
      v,
      j,
      n,
      nprime,
      ndprime,
      s,
      sprime,
      sdprime,
      update,
      count
    )
  end
end

MuvNUTSState(pstate::ParameterState{Continuous, Multivariate}, tune::MCTunerState=BasicMCTune()) =
  MuvNUTSState(
    pstate,
    pstate,
    pstate,
    pstate,
    tune,
    NaN,
    NaN,
    Array(eltype(pstate), pstate.size),
    Array(eltype(pstate), pstate.size),
    Array(eltype(pstate), pstate.size),
    Array(eltype(pstate), pstate.size),
    NaN,
    NaN,
    NaN,
    0,
    0,
    0,
    0,
    0,
    true,
    true,
    true,
    true,
    0
  )

immutable NUTS <: HMCSampler
  leapstep::Real
  maxδ::Integer
  maxndoublings::Integer

  function NUTS(leapstep::Real, maxδ::Integer, maxndoublings::Integer)
    @assert leapstep > 0 "Leapfrog step is not positive"
    @assert maxδ > 0 "maxδ is not positive"
    @assert maxndoublings > 0 "Maximum number of doublings is not positive"
    new(leapstep, maxδ, maxndoublings)
  end
end

NUTS(leapstep::Real=0.1; maxδ::Integer=1000, maxndoublings::Integer=5) = NUTS(leapstep, maxδ, maxndoublings)

function initialize!{F<:VariateForm}(
  pstate::ParameterState{Continuous, F},
  parameter::Parameter{Continuous, F},
  sampler::NUTS,
  outopts::Dict
)
  parameter.uptogradlogtarget!(pstate)
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"
  @assert all(isfinite(pstate.gradlogtarget)) "Gradient of log-target not finite: initial values out of support"

  if !isempty(outopts[:diagnostics])
    pstate.diagnostickeys = copy(outopts[:diagnostics])
    pstate.diagnosticvalues = Array(Any, length(pstate.diagnostickeys))
  end
end

sampler_state(
  parameter::Parameter{Continuous, Multivariate},
  sampler::NUTS,
  tuner::MCTuner,
  pstate::ParameterState{Continuous, Multivariate},
  vstate::VariableStateVector
) =
  MuvNUTSState(generate_empty(pstate), tuner_state(parameter, sampler, tuner))

uturn(positionplus::RealVector, positionminus::RealVector, momentumplus::RealVector, momentumminus::RealVector) =
  dot(positionplus-positionminus, momentumplus) < 0. || dot(positionplus-positionminus, momentumminus) < 0.

function build_tree!(
  sstate::NUTSState{Multivariate},
  pstate::ParameterState{Continuous, Multivariate},
  momentum::RealVector,
  u::Real,
  v::Real,
  j::Integer,
  step::Real,
  logtarget!::Function,
  gradlogtarget!::Function,
  sampler::NUTS
)
  if j == 0 # Base case: take one leapfrog step in the direction v
    local hamiltonianprime::Real

    leapfrog!(sstate.pstateprime, sstate.momentumprime, pstate, momentum, v*sstate.tune.step, gradlogtarget!)
    logtarget!(sstate.pstateprime)
    hamiltonianprime = hamiltonian(sstate.pstateprime.logtarget, sstate.momentumprime)

    sstate.nprime = Int(u <= hamiltonianprime)
    sstate.sprime = u < sampler.maxδ+hamiltonianprime

    return sstate.pstateprime,
      sstate.momentumprime,
      sstate.pstateprime,
      sstate.momentumprime,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime
  else # Recursion: implicitly build the left and right subtrees
    sstate.pstateminus,
    sstate.momentumminus,
    sstate.pstateplus,
    sstate.momentumplus,
    sstate.pstateprime,
    sstate.nprime,
    sstate.sprime =
      build_tree!(sstate, pstate, momentum, u, v, j-1, step, logtarget!, gradlogtarget!, sampler)

    if sstate.sprime
      if v == -1
        sstate.pstateminus, sstate.momentumminus, _, _, sstate.pstatedprime, sstate.ndprime, sstate.sdprime =
          build_tree!(
            sstate,
            sstate.pstateminus,
            sstate.momentumminus,
            u,
            v,
            j-1,
            step,
            logtarget!,
            gradlogtarget!,
            sampler
          )
      else
        _, _, sstate.pstateplus, sstate.momentumplus, sstate.pstatedprime, sstate.ndprime, sstate.sdprime =
          build_tree!(
            sstate,
            sstate.pstateplus,
            sstate.momentumplus,
            u,
            v,
            j-1,
            step,
            logtarget!,
            gradlogtarget!,
            sampler
          )
      end

      if rand() <= sstate.ndprime/(sstate.ndprime+sstate.nprime)
        sstate.pstateprime.value = copy(sstate.pstatedprime.value)
        sstate.pstateprime.gradlogtarget = copy(sstate.pstatedprime.gradlogtarget)
      end

      sstate.nprime += sstate.ndprime

      sstate.sprime =
        sstate.sdprime &&
        !uturn(sstate.pstateplus.value, sstate.pstateminus.value, sstate.momentumplus, sstate.momentumminus)
    end

    return sstate.pstateminus,
      sstate.momentumminus,
      sstate.pstateplus,
      sstate.momentumplus,
      sstate.pstateprime,
      sstate.nprime,
      sstate.sprime
  end
end


text:
function isobj(filepath)
    base, ext = splitext(filepath)
    ext == ".obj" ? true : false
end

"""
    loadobj(shapepath::String; scale=1, message=true) -> nodes, faces
"""
function loadobj(shapepath::String; scale=1, static=true, message=true)

    nodes = Vector{Float64}[]
    faces = Vector{Int64}[]

    open(shapepath, "r") do f
        for line in eachline(f)
            line == "" && continue  # Skip blank lines
            data = split(line)

            if data[1] == "v"
                node = parse.(Float64, data[2:4])  # (x, y, z) coordinates [km]
                push!(nodes, node)
            elseif data[1] == "f"
                face = parse.(Int64, data[2:4])  # indices of three vertices
                push!(faces, face)
            end
        end
    end

    nodes *= scale  # if scale is 1000, converted [km] to [m]

    if static == true
        nodes = [SVector{3,Float64}(node) for node in nodes]
        faces = [SVector{3,Int64}(face) for face in faces]
    end

    if message == true
        println("+-----------------------------+")
        println("|        Load OBJ file        |")
        println("+-----------------------------+")
        println(" Nodes: ", length(nodes))
        println(" Faces: ", length(faces))
    end

    return nodes, faces
end

"""
    saveobj(path::String, nodes, faces)
"""
function saveobj(path::String, nodes, faces)
    open(path, "w") do f
        println(f, "# Generated by obj.jl written by M. Kanamaru")
        println(f, "# Vertices: ", length(nodes))
        println(f, "# Faces:    ", length(faces))
        # println(f)

        for node in nodes
            println(f, "v  ", node[1], "  ", node[2], "  ", node[3])
        end

        println(f)

        for face in faces
            println(f, "f  ", face[1], "  ", face[2], "  ", face[3])
        end
    end
end

"""
Make a perturbed shape model
"""
function perturbnodes!(nodes, noise_weight)

    for node in nodes
        noise = rand()*2 - 1  # random number in [-1, 1]
        noise *= noise_weight

        dr = normalize(node) * noise / 1000  # [km]
        node .+= dr
    end
    nodes
end

perturbnodes(nodes, noise_weight) = perturbnodes!(deepcopy(nodes), noise_weight)

# ****************************************************************

getFaceCenters(nodes, faces) = [getcenter(nodes[face]) for face in faces]
getFaceNormals(nodes, faces) = [getnormal(nodes[face]) for face in faces]
getFaceAreas(nodes, faces) = [getarea(nodes[face]) for face in faces]

max_stars_repo_name: MasanoriKanamaru/Astroshaper

max_issues_repo_name: MasanoriKanamaru/Astroshaper

max_forks_repo_name: MasanoriKanamaru/Astroshaper


text:
"""
This function reads the RTS-GMLC timeseries data and parses it into the format
compatible with AgentSimulation
"""

function read_rts(data_dir::String,
                  test_system_dir::String,
                  base_dir::String,
                  annual_growth::AxisArrays.AxisArray{Float64, 2},
                  start_year::Int64,
                  n_rep_days::Int64)

    test_system_load = DataFrames.DataFrame(CSV.File(joinpath(test_system_dir, "RTS_Data", "timeseries_data_files", "Load", "DAY_AHEAD_regional_Load.csv")))

    base_year = test_system_load[1, "Year"]

    test_sys_num_hours = DataFrames.nrow(test_system_load)
    if test_sys_num_hours >= 8760
        test_sys_hour_weight = ones(test_sys_num_hours)
    else
        test_sys_hour_weight = ones(test_sys_num_hours) * 8760 / test_sys_num_hours
    end

    zone_numbers = names(test_system_load)[5:end]
    zones = ["zone_$(i)" for i in zone_numbers]

    test_system_reserve_up = test_system_load[:, 1:4]
    test_system_reserve_down = test_system_load[:, 1:4]

    for zone in zone_numbers
        reserve_up_df = DataFrames.DataFrame(CSV.File(joinpath(test_system_dir, "RTS_Data", "timeseries_data_files", "Reserves", "DAY_AHEAD_regional_Reg_Up_R$(zone).csv")))

        test_system_reserve_up[:, zone] = reserve_up_df[:, "Reg_Up_R$(zone)"]

        reserve_down_df = DataFrames.DataFrame(CSV.File(joinpath(test_system_dir, "RTS_Data", "timeseries_data_files", "Reserves", "DAY_AHEAD_regional_Reg_Down_R$(zone).csv")))
        test_system_reserve_down[:, zone] = zeros(DataFrames.nrow(test_system_reserve_down))
        for d in 1:Int(DataFrames.nrow(test_system_reserve_down)/24)
            for h in 1:24
                test_system_reserve_down[(d - 1) * 24 + h, zone] = reserve_down_df[d, Symbol(h)]
            end
        end
    end

    scaled_test_sys_load = deepcopy(test_system_load)
    scaled_test_sys_reserve_up = deepcopy(test_system_reserve_up)
    scaled_test_sys_reserve_down = deepcopy(test_system_reserve_down)

    scaled_test_sys_load[:, "Year"] = fill(start_year, DataFrames.nrow(scaled_test_sys_load))
    remove_leap_day!(scaled_test_sys_load, start_year)

    scaled_test_sys_reserve_up[:, "Year"] = fill(start_year, DataFrames.nrow(scaled_test_sys_reserve_up))
    remove_leap_day!(scaled_test_sys_reserve_up, start_year)

    scaled_test_sys_reserve_down[:, "Year"] = fill(start_year, DataFrames.nrow(scaled_test_sys_reserve_down))
    remove_leap_day!(scaled_test_sys_reserve_down, start_year)

    for y in 1:(start_year - base_year)
        for zone in zone_numbers
            scaled_test_sys_load[:, "$(zone)"] =  scaled_test_sys_load[:, "$(zone)"] * (1 + annual_growth[y, Symbol("load_zone_$(zone)")])
            scaled_test_sys_reserve_up[:, "$(zone)"] =  scaled_test_sys_reserve_up[:, "$(zone)"] * (1 + annual_growth[y, Symbol("load_zone_$(zone)")])
            scaled_test_sys_reserve_down[:, "$(zone)"] =  scaled_test_sys_reserve_down[:, "$(zone)"] * (1 + annual_growth[y, Symbol("load_zone_$(zone)")])
        end
    end

    net_load_df = scaled_test_sys_load[:, 1:4]

    for zone in zone_numbers
        net_load_df[:, "load_zone_$(zone)"] = scaled_test_sys_load[:, zone]
    end

    existing_generator_data = DataFrames.DataFrame(CSV.File(joinpath(test_system_dir, "RTS_Data", "SourceData", "gen.csv")))

    wind_timeseries_file = joinpath(test_system_dir, "RTS_Data", "timeseries_data_files", "WIND", "DAY_AHEAD_wind.csv")
    if isfile(wind_timeseries_file)
        wind_timeseries_data = DataFrames.DataFrame(CSV.File(wind_timeseries_file))
        remove_leap_day!(wind_timeseries_data, start_year)
    end

    pv_timeseries_file = joinpath(test_system_dir, "RTS_Data", "timeseries_data_files", "PV", "DAY_AHEAD_pv.csv")
    if isfile(pv_timeseries_file)
        pv_timeseries_data = DataFrames.DataFrame(CSV.File(pv_timeseries_file))
        remove_leap_day!(pv_timeseries_data, start_year)
    end

    rtpv_timeseries_file = joinpath(test_system_dir, "RTS_Data", "timeseries_data_files", "RTPV", "DAY_AHEAD_rtpv.csv")
    if isfile(rtpv_timeseries_file)
        rtpv_timeseries_data = DataFrames.DataFrame(CSV.File(rtpv_timeseries_file))
        remove_leap_day!(rtpv_timeseries_data, start_year)
    end

    hydro_timeseries_file = joinpath(test_system_dir, "RTS_Data", "timeseries_data_files", "Hydro", "DAY_AHEAD_hydro.csv")
    if isfile(hydro_timeseries_file)
        hydro_timeseries_data = DataFrames.DataFrame(CSV.File(hydro_timeseries_file))
        remove_leap_day!(hydro_timeseries_data, start_year)
    end

    gen_availability_df = scaled_test_sys_load[:, 1:4]

    for i in 1:DataFrames.nrow(existing_generator_data)
        if existing_generator_data[i, "Unit Type"] == "WIND"
            gen_availability_df[:, existing_generator_data[i, "GEN UID"]] = wind_timeseries_data[:, existing_generator_data[i, "GEN UID"]] / existing_generator_data[i, "PMax MW"]
            net_load_df[:, existing_generator_data[i, "GEN UID"]] = wind_timeseries_data[:, existing_generator_data[i, "GEN UID"]]
        elseif existing_generator_data[i, "Unit Type"] == "PV"
            gen_availability_df[:, existing_generator_data[i, "GEN UID"]] = pv_timeseries_data[:, existing_generator_data[i, "GEN UID"]] / existing_generator_data[i, "PMax MW"]
            net_load_df[:, existing_generator_data[i, "GEN UID"]] = pv_timeseries_data[:, existing_generator_data[i, "GEN UID"]]
        elseif existing_generator_data[i, "Unit Type"] == "RTPV"
            gen_availability_df[:, existing_generator_data[i, "GEN UID"]] = rtpv_timeseries_data[:, existing_generator_data[i, "GEN UID"]] / existing_generator_data[i, "PMax MW"]
        elseif existing_generator_data[i, "Unit Type"] == "HYDRO"
            gen_availability_df[:, existing_generator_data[i, "GEN UID"]] = hydro_timeseries_data[:, existing_generator_data[i, "GEN UID"]] / existing_generator_data[i, "PMax MW"]
        else
            gen_availability_df[:, existing_generator_data[i, "GEN UID"]] = ones(DataFrames.nrow(gen_availability_df))
        end
    end

    write_data(joinpath(data_dir, "timeseries_data_files", "Availability"), "DAY_AHEAD_availability.csv", gen_availability_df)
    write_data(joinpath(data_dir, "timeseries_data_files", "Net Load Data"), "load_n_vg_data.csv", net_load_df)

    representative_days = find_representative_days(data_dir, base_dir, n_rep_days)
    rep_load_data = filter(row -> in(Dates.Date(row[:Year], row[:Month], row[:Day]), keys(representative_days)), scaled_test_sys_load)

    rep_hour_weight = zeros(DataFrames.nrow(rep_load_data))

    for i in 1:length(rep_hour_weight)
        rep_hour_weight[i] = representative_days[Dates.Date(rep_load_data[i, :Year], rep_load_data[i, :Month], rep_load_data[i, :Day])]
    end

    write_data(joinpath(data_dir, "timeseries_data_files", "Load"), "load_0.csv", scaled_test_sys_load)
    write_data(joinpath(data_dir, "timeseries_data_files", "Load"), "rep_load_0.csv", rep_load_data)

    system_peak_load = maximum(sum(scaled_test_sys_load[:, zone] for zone in zone_numbers))

    rep_reserve_up_data = filter(row -> in(Dates.Date(row[:Year], row[:Month], row[:Day]), keys(representative_days)), scaled_test_sys_reserve_up)
    write_data(joinpath(data_dir, "timeseries_data_files", "Reserves"), "reserve_up_0.csv", scaled_test_sys_reserve_up)
    write_data(joinpath(data_dir, "timeseries_data_files", "Reserves"), "rep_reserve_up_0.csv", rep_reserve_up_data)

    rep_reserve_down_data = filter(row -> in(Dates.Date(row[:Year], row[:Month], row[:Day]), keys(representative_days)), scaled_test_sys_reserve_down)
    write_data(joinpath(data_dir, "timeseries_data_files", "Reserves"), "reserve_down_0.csv", scaled_test_sys_reserve_down)
    write_data(joinpath(data_dir, "timeseries_data_files", "Reserves"), "rep_reserve_down_0.csv", rep_reserve_down_data)

    reserve_params_df = read_data(joinpath(test_system_dir, "RTS_Data", "SourceData", "reserves.csv"))

    reserve_up_eligible_categories = String[]
    reserve_down_eligible_categories = String[]

    for zone in zone_numbers
        push!(reserve_up_eligible_categories, reserve_params_df[findfirst(x-> x == "Reg_Up_R$(zone)", reserve_params_df[:, "Reserve Product"]), "Eligible Device SubCategories"])
        push!(reserve_down_eligible_categories, reserve_params_df[findfirst(x-> x == "Reg_Down_R$(zone)", reserve_params_df[:, "Reserve Product"]), "Eligible Device SubCategories"])
    end

    reserve_up_mkt_df = read_data(joinpath(data_dir, "markets_data", "reserve_up_mkt_param.csv"))
    reserve_down_mkt_df = read_data(joinpath(data_dir, "markets_data", "reserve_down_mkt_param.csv"))

    reserve_up_mkt_df[!, "eligible categories"] = reserve_up_eligible_categories
    reserve_down_mkt_df[!, "eligible categories"] = reserve_down_eligible_categories

    write_data(joinpath(data_dir, "markets_data"), "reserve_up_mkt_param.csv", reserve_up_mkt_df)
    write_data(joinpath(data_dir, "markets_data"), "reserve_down_mkt_param.csv", reserve_down_mkt_df)

    # Create zonal lines

    zonal_lines = ZonalLine[]

    branches = read_data(joinpath(test_system_dir, "RTS_Data", "SourceData", "branch.csv"))
    dc_branches = read_data(joinpath(test_system_dir, "RTS_Data", "SourceData", "dc_branch.csv"))

    for b in 1:DataFrames.nrow(branches)
        from_bus = "$(branches[b, "From Bus"])"
        from_zone = "zone_$(first(from_bus, 1))"

        to_bus = "$(branches[b, "To Bus"])"
        to_zone = "zone_$(first(to_bus, 1))"

        similar_line = filter(l -> (in(from_zone, [get_from_zone(l), get_to_zone(l)]) && in(to_zone, [get_from_zone(l), get_to_zone(l)])), zonal_lines)

        if length(similar_line) < 1
            push!(zonal_lines, ZonalLine(branches[b, "UID"], from_zone, to_zone, branches[b, "Cont Rating"]))
        else
            set_active_power_limit!(similar_line[1], get_active_power_limit(similar_line[1]) + branches[b, "Cont Rating"])
        end
    end

    for b in 1:DataFrames.nrow(dc_branches)
        from_bus = "$(dc_branches[b, "From Bus"])"
        from_zone = "zone_$(first(from_bus, 1))"

        to_bus = "$(dc_branches[b, "To Bus"])"
        to_zone = "zone_$(first(to_bus, 1))"

        similar_line = filter(l -> (in(from_zone, [get_from_zone(l), get_to_zone(l)]) && in(to_zone, [get_from_zone(l), get_to_zone(l)])), zonal_lines)

        if length(similar_line) < 1
            push!(zonal_lines, ZonalLine(dc_branches[b, "UID"], from_zone, to_zone, dc_branches[b, "MW Load"]))
        else
            set_active_power_limit!(similar_line[1], get_active_power_limit(similar_line[1]) + dc_branches[b, "MW Load"])
        end
    end

    return zones, representative_days, rep_hour_weight, system_peak_load, test_sys_hour_weight, zonal_lines
end


text:
using Test, LogSynth

@test length(SkipListDistribution(5)) == 0

let s = SkipListDistribution(5)
    push!(s, 1.0)
    push!(s, 2.0)
    push!(s, 1.0)
    @test length(s) == 3
    @test LogSynth.total(s) == 4.0
    @test LogSynth.getindexvector(s, 2)[end] == 2
    @test LogSynth.getindexvector(s, 0.24999)[end] == 1
    @test LogSynth.getindexvector(s, 0.25)[end] == 2
    @test LogSynth.getindexvector(s, 0.5)[end] == 2
    @test LogSynth.getindexvector(s, 0.75)[end] == 3

    @test LogSynth.inc!(s, 0.5, 1.0) == 2
    @test length(s) == 3
    @test LogSynth.total(s) == 5.0
    @test LogSynth.getindexvector(s, 0.19999)[end] == 1
    @test LogSynth.getindexvector(s, 0.2)[end] == 2
    @test LogSynth.getindexvector(s, 0.7999)[end] == 2
    @test LogSynth.getindexvector(s, 0.8)[end] == 3
end

let s = SkipListDistribution(5)
    for i in 1:100_000
        push!(s, 1.0)
    end

    # verify that the lengths of the different levels are about right
    # 5σ gives us about 1.5 x 10^-6 chance of random failure
    p = 1.0
    for i in s.height:-1:1
        σ = sqrt(p * (1-p) * length(s))
        @test length(s.child[i]) ≈ length(s) * p atol = 5σ
        p *= s.skipProbability
    end
end


text:
# #-------------------------------------------------------------------------------------------------------
# #-------------------------------------------------------------------------------------------------------
# # This file contains functions that are used for various climate projection and SC-CH4 calculations.
# #-------------------------------------------------------------------------------------------------------
# #-------------------------------------------------------------------------------------------------------

#######################################################################################################################
# CREATE FOLDER STRUCTURE TO STORE RESULTS.
#######################################################################################################################
# Description: This function creates a directory structure to store all of the model output in the "Results" folder.
#              The main folder will have the user assigned "results_name".
#
# Function Arguments:
#
#       results_name = A vector of global carbon tax values to be optimized.
#----------------------------------------------------------------------------------------------------------------------

function build_result_folders(results_name::String)

    # Set SNEASY+CH4 and IAM model names.
    climate_model_names = ["s_fair", "s_fund", "s_hector", "s_magicc"]
    iam_model_names     = ["dice", "fund"]

    results_root_path = joinpath(@__DIR__, "..", "..", "results")

    # Loop over the various IAM and climate model combinations.
    for climate_model in climate_model_names

        # Create folders to store calibrated posterior parameter samples for baseline and wider prior sensitivity runs.
        mkpath(joinpath(results_root_path, results_name, "calibrated_parameters", climate_model))
        mkpath(joinpath(results_root_path, results_name, "calibrated_parameters", "wider_priors", climate_model))

        # Create folders to store climate projections for the different scenarios.
        mkpath(joinpath(results_root_path, results_name, "climate_projections", "baseline_run", climate_model))
        mkpath(joinpath(results_root_path, results_name, "climate_projections", "outdated_ch4_forcing", climate_model))
        mkpath(joinpath(results_root_path, results_name, "climate_projections", "remove_correlations", climate_model))
        mkpath(joinpath(results_root_path, results_name, "climate_projections", "us_climate_sensitivity", climate_model))
        mkpath(joinpath(results_root_path, results_name, "climate_projections", "rcp26", climate_model))
        mkpath(joinpath(results_root_path, results_name, "climate_projections", "wider_priors", climate_model))

        # Create folders to store SC-CH4 estimates for DICE and FUND for baseline and alternative scenario runs.
        for iam_model in iam_model_names
            mkpath(joinpath(results_root_path, results_name, "scch4_estimates", "baseline_run", iam_model, climate_model))
            mkpath(joinpath(results_root_path, results_name, "scch4_estimates", "outdated_ch4_forcing", iam_model, climate_model))
            mkpath(joinpath(results_root_path, results_name, "scch4_estimates", "remove_correlations", iam_model, climate_model))
            mkpath(joinpath(results_root_path, results_name, "scch4_estimates", "us_climate_sensitivity", iam_model, climate_model))
            mkpath(joinpath(results_root_path, results_name, "scch4_estimates", "rcp26", iam_model, climate_model))
            mkpath(joinpath(results_root_path, results_name, "scch4_estimates", "wider_priors", iam_model, climate_model))
        end
    end

    # Create folder to store Bayesian model averaging (BMA) weights within the calibration folder.
    mkpath(joinpath(results_root_path, results_name, "calibrated_parameters", "bma_weights"))

    # Create folder to save equity-weighted SC-CH4 estimates for FUND + S-MAGICC.
    mkpath(joinpath(results_root_path, results_name, "scch4_estimates", "equity_weighting", "fund", "s_magicc"))

    # Create folder to store .jpg and .pdf figures.
    mkpath(joinpath(results_root_path, results_name, "figures", "jpg_figures"))
    mkpath(joinpath(results_root_path, results_name, "figures", "pdf_figures"))
end

#######################################################################################################################
# CREATE A FUNCTION TO UPDATE PARAMETERS FOR SPECIFIC METHANE CYCLE MODEL.
#######################################################################################################################
# Description: This function will update the uncertain parameters specific to the methane cycle being used.
#
# Function Arguments:
#
#       climate_model = A symbol identifying the specific version of SNEASY+CH4 (options are :sneasy_fair, :sneasy_fund,
#                       :sneasy_hector, and :sneasy_magicc).
#----------------------------------------------------------------------------------------------------------------------

function create_update_ch4_function(climate_model::Symbol)

    # Create a function to update the uncertain CH₄ parameters, given the selected CH₄ model.
    update_ch4_params! =

        if climate_model == :sneasy_fair

            function(mi, params::Array{Float64,1})
                CH₄_0         = params[12]
                N₂O_0         = params[13]
                rf_scale_CH₄  = params[17]
                τ_troposphere = params[22]
                CH₄_natural   = params[23]

                update_param!(mi, :natural_emiss_CH₄, ones(length(dim_keys(mi.md, :time))) .* CH₄_natural)
                update_param!(mi, :CH₄_0,             CH₄_0)
                update_param!(mi, :τ_CH₄,             τ_troposphere)
                update_param!(mi, :N₂O_0,             N₂O_0)
                update_param!(mi, :scale_CH₄,         rf_scale_CH₄)
                return
            end

        elseif climate_model == :sneasy_fund

            function(mi, params::Array{Float64,1})
                CH₄_0         = params[12]
                N₂O_0         = params[13]
                rf_scale_CH₄  = params[17]
                τ_troposphere = params[22]

                update_param!(mi, :ch4pre,    CH₄_0)
                update_param!(mi, :acch4_0,   CH₄_0)
                update_param!(mi, :lifech4,   τ_troposphere)
                update_param!(mi, :CH₄_0,     CH₄_0)
                update_param!(mi, :N₂O_0,     N₂O_0)
                update_param!(mi, :scale_CH₄, rf_scale_CH₄)
                return
            end

        elseif climate_model == :sneasy_hector

            function(mi, params::Array{Float64,1})
                CH₄_0          = params[12]
                N₂O_0          = params[13]
                rf_scale_CH₄   = params[17]
                τ_troposphere  = params[22]
                CH₄_natural    = params[23]
                τ_soil         = params[24]
                τ_stratosphere = params[25]

                update_param!(mi, :TOH0,      τ_troposphere)
                update_param!(mi, :M0,        CH₄_0)
                update_param!(mi, :CH4N,      CH₄_natural)
                update_param!(mi, :Tsoil,     τ_soil)
                update_param!(mi, :Tstrat,    τ_stratosphere)
                update_param!(mi, :CH₄_0,     CH₄_0)
                update_param!(mi, :N₂O_0,     N₂O_0)
                update_param!(mi, :scale_CH₄, rf_scale_CH₄)
                return
            end

        elseif climate_model == :sneasy_magicc

            function(mi, params::Array{Float64,1})
                CH₄_0          = params[12]
                N₂O_0          = params[13]
                rf_scale_CH₄   = params[17]
                τ_troposphere  = params[22]
                CH₄_natural    = params[23]
                τ_soil         = params[24]
                τ_stratosphere = params[25]

                update_param!(mi, :CH₄_0,       CH₄_0)
                update_param!(mi, :CH4_natural, CH₄_natural)
                update_param!(mi, :TAUSOIL,     τ_soil)
                update_param!(mi, :TAUSTRAT,    τ_stratosphere)
                update_param!(mi, :TAUINIT,     τ_troposphere)
                update_param!(mi, :N₂O_0,       N₂O_0)
                update_param!(mi, :scale_CH₄,   rf_scale_CH₄)
                return
            end

        else
            error("Incorrect climate model selected. Options include :sneasy_fair, :sneasy_fund, :sneasy_hector, or :sneasy_magicc.")
        end

    # Return the paramter updated function for the specific version of SNEASY+CH4.
    return update_ch4_params!
end

#######################################################################################################################
# CREATE A FUNCTION TO ACCESS PROJECTED METHANE CONCENTRATIONS FOR SPECIFIC METHANE CYCLE MODEL.
#######################################################################################################################
# Description: This function will access the CH₄ concentration projections given the particular version of
#              SNEASY+CH4 (component and variable names differ across model versions).
#
# Function Arguments:
#
#       climate_model = A symbol identifying the specific version of SNEASY+CH4 (options are :sneasy_fair, :sneasy_fund,
#                       :sneasy_hector, and :sneasy_magicc).
#----------------------------------------------------------------------------------------------------------------------

function create_get_ch4_results_function(climate_model::Symbol)

    # Create a function to get CH₄ concentration projections, given the selected CH₄ model.
    get_ch4_results! =

        if climate_model == :sneasy_fair || climate_model == :sneasy_magicc

            function(m)
                return m[:ch4_cycle, :CH₄]
            end

        elseif climate_model == :sneasy_fund

            function(m)
                return m[:climatech4cycle, :acch4]
            end

        elseif climate_model == :sneasy_hector

            function(m)
                return m[:ch4_cycle, :CH4]
            end
        end

    return get_ch4_results!
end

#######################################################################################################################
# CALCULATE DISCOUNT FACTORS
#######################################################################################################################
# Description: This function calculates discount factors given user-specified discounting settings.
#
# Function Arguments:
#
#       df               = A pre-allocated array to store discount factors.
#       pc_consumption   = Per capita consumption values.
#       ρ                = Pure rate of time preference.
#       η                = Elasticity of marginal utility of consumption.
#       years            = Vector of years that the model is being run for.
#       pulse_year       = The year that the CH₄ emission pulse is added.
#       pulse_year_index = The timestep index corresponding to the CH₄ emission pulse year.
#----------------------------------------------------------------------------------------------------------------------

function calculate_discount_factors!(df::Array{Float64,1}, pc_consumption::Array{Float64,1}, ρ::Float64, η::Float64, years::Array{Int,1}, pulse_year::Int, pulse_year_index::Int)
    for t=pulse_year_index:length(years)
        df[t] = (pc_consumption[pulse_year_index] / pc_consumption[t])^η * 1.0 / (1.0 + ρ)^(years[t]-pulse_year)
    end
end

#######################################################################################################################
# CALCULATE CO₂ RADIATIVE FORCING SCALAR
#######################################################################################################################
# Description: This function calculates a scaling coefficient so CO₂ radiative forcing is consistent with user supplied
#              value for the radiative forcing due to doubling CO₂ (following Smith et al. (2018), Geosci. Model Dev.).
#
# Function Arguments:
#
#       F2x   = Radiative forcing due to doubling carbon dioxide (Wm⁻²).
#       CO₂_0 = Initial (pre-industrial) carbon dioxide concentration (ppm).
#       N₂O_0 = Initial (pre-industrial) nitrous oxide concentration (ppb).
#----------------------------------------------------------------------------------------------------------------------

function co2_rf_scale(F2x::Float64, CO₂_0::Float64, N₂O_0::Float64)

    # Calcualte forcing from doubling of CO₂.
    co2_doubling = (-2.4e-7 * CO₂_0^2 + 7.2e-4 * CO₂_0 - 2.1e-4 * N₂O_0 + 5.36) * log(2)

    # Calculate scaling factor, given user supplied F2x parameter.
    CO₂_scale = F2x / co2_doubling

    return CO₂_scale
end

#######################################################################################################################
# SIMULATE STATIONARY AR(1) PROCESS WITH TIME VARYING OBSERVATION ERRORS.
#######################################################################################################################
# Description: This function simulates a stationary AR(1) process (given time-varying observation errors supplied with
#              each calibration data set) to superimpose noise onto the climate model projections.
#
# Function Arguments:
#
#       n = Number of time periods (years) the model is being run for.
#       σ = Calibrated standard deviation.
#       ρ = Calibrated autocorrelation coefficient.
#       ϵ = Time-varying observation errors.
#----------------------------------------------------------------------------------------------------------------------

function simulate_ar1_noise(n::Int, σ::Float64, ρ::Float64, ϵ::Array{Float64,1})

    # Define AR(1) stationary process variance.
    σ_process = σ^2/(1-ρ^2)

    # Initialize AR(1) covariance matrix (just for convenience).
    H = abs.(collect(1:n)' .- collect(1:n))

    # Calculate residual covariance matrix (sum of AR(1) process variance and observation error variances).
    # Note: This follows Supplementary Information Equation (10) in Ruckert et al. (2017).
    cov_matrix = σ_process * ρ .^ H + Diagonal(ϵ.^2)

    # Return a mean-zero AR(1) noise sample accounting for time-varying observation error.
    return rand(MvNormal(cov_matrix))
end

#######################################################################################################################
# SIMULATE STATIONARY CAR(1) PROCESS WITH TIME VARYING OBSERVATION ERRORS.
#######################################################################################################################
# Description: This function simulates a stationary CAR(1) process (given time-varying observation errors supplied with
#              each calibration data set) to superimpose noise onto the climate model projections.
#
# Function Arguments:
#
#       n              = Number of time periods (years) the model is being run for.
#       α₀             = Calibrated term describing correlation memory of CAR(1) process.
#       σ²_white_noise = Calibrated continuous white noise process variance term.
#       ϵ              = Time-varying observation errors.
#----------------------------------------------------------------------------------------------------------------------

function simulate_car1_noise(n, α₀, σ²_white_noise, ϵ)

    # Indices for full time horizon.
    indices = collect(1:n)

    # Initialize covariance matrix for irregularly spaced data with relationships decaying exponentially.
    H = exp.(-α₀ .* abs.(indices' .- indices))

    # Define the variance of x(t), a continous stochastic time-series.
    σ² = σ²_white_noise / (2*α₀)

    # Calculate residual covariance matrix (sum of CAR(1) process variance and observation error variances).
    cov_matrix = σ² .* H + Diagonal(ϵ.^2)

    # Return a mean-zero CAR(1) noise sample accounting for time-varying observation error.
    return rand(MvNormal(cov_matrix))
end

#######################################################################################################################
# REPLICATE TIME-VARYING OBSERVATION ERRORS FOR PERIODS WITHOUT COVERAGE.
#######################################################################################################################
# Description: This function creates a time-series of observation errors for the entire model time horizon. For years
#              without observation error estimates, the error remains constant at the average of the ten nearest error
#              values in time.
#
# Function Arguments:
#
#       start_year = The first year to run the climate model.
#       end_year   = The final year to run the climate model.
#       error_data = A vector of time-varying observation errors supplied with each calibration data set.
#----------------------------------------------------------------------------------------------------------------------

function replicate_errors(start_year::Int, end_year::Int, error_data)

    # Initialize model years and measurement error vector.
    model_years = collect(start_year:end_year)
    errors = zeros(length(model_years))

    # Find indices for periods that have observation errors.
    err_indices = findall(x-> !ismissing(x), error_data)

    # Replicate errors for all periods prior to start of observations based on average of 10 nearest errors in time.
    errors[1:(err_indices[1]-1)] .= mean(error_data[err_indices[1:10]])

    # Add errors for periods with observations.
    errors[err_indices[1:end]] = error_data[err_indices[1:end]]

    # Replicate errors for all periods after observation data based on average of 10 nearest errors in time.
    errors[(err_indices[end]+1):end] .= mean(error_data[err_indices[(end-9):end]])

    return errors
end

#######################################################################################################################
# LINEARLY INTERPOLATE DICE RESULTS TO ANNUAL VALUES
#######################################################################################################################
# Description: This function uses linear interpolation to create an annual time series from DICE results (which have
#              five year timesteps).
#
# Function Arguments:
#
#       data    = The DICE results to be interpolated
#       spacing = Length of model time steps (5 for DICE).
#----------------------------------------------------------------------------------------------------------------------

function dice_interpolate(data, spacing)

    # Create an interpolation object for the data (assume first and last points are end points, e.g. no interpolation beyond support).
    interp_linear = interpolate(data, BSpline(Linear()), OnGrid())

    # Create points to interpolate for (based on spacing term).
    interp_points = collect(1:(1/spacing):length(data))

    # Carry out interpolation.
    return interp_linear[interp_points]
end

#######################################################################################################################
# CREATE CLIMATE PROJECTION CREDIBLE INTERVALS
#######################################################################################################################
# Description: This function calculates the upper and lower credible interval ranges for the climate projections carried
#              out with each posterior parameter sample. It does so for two different % levels.
#
# Function Arguments:
#
#       years          = The years the model projection is carried out for.
#       model_result   = An array of model projections (each row is a new projection, each column is a new year).
#       conf_1_percent = First percentage value to calculate credible interval over (i.e. 0.95 corresponds to 95% credible interval).
#       conf_2_percent = Second percentage value to calculate credible interval over.
#----------------------------------------------------------------------------------------------------------------------

function get_confidence_interval(years, model_result, conf_1_percent, conf_2_percent)

    # Set intervals for quantile function and calculate total number of years in results.
    α1      = 1-conf_1_percent
    α2      = 1-conf_2_percent
    n_years = length(years)

    # Initialize dataframe of results with a column of years, mean results, and missing values.
    ci_results = DataFrame(Year=years, Mean=vec(mean(model_result, dims=1)), Lower1=zeros(Union{Missing,Float64}, n_years), Upper1=zeros(Union{Missing,Float64}, n_years), Lower2=zeros(Union{Missing,Float64}, n_years), Upper2=zeros(Union{Missing,Float64}, n_years))

    # Calculate credible intervals for each year (CI for years with 'missing' values also set to 'missing').
    for i in 1:n_years
        if all(x-> x !== missing, model_result[:,i])
            ci_results.Lower1[i] = quantile(model_result[:,i], α1/2)
            ci_results.Upper1[i] = quantile(model_result[:,i], 1-α1/2)
            ci_results.Lower2[i] = quantile(model_result[:,i], α2/2)
            ci_results.Upper2[i] = quantile(model_result[:,i], 1-α2/2)
        else
            ci_results.Lower1[i] = missing
            ci_results.Upper1[i] = missing
            ci_results.Lower2[i] = missing
            ci_results.Upper2[i] = missing
        end
    end

    # Rename columns to have names specific to user-provided credible interval percentages.
    rename!(ci_results, :Lower1 => Symbol(join(["LowerConf" conf_1_percent], '_')))
    rename!(ci_results, :Upper1 => Symbol(join(["UpperConf" conf_1_percent], '_')))
    rename!(ci_results, :Lower2 => Symbol(join(["LowerConf" conf_2_percent], '_')))
    rename!(ci_results, :Upper2 => Symbol(join(["UpperConf" conf_2_percent], '_')))

    return ci_results
end


text:
using runtests.jl

function runtest()
  # do something
end


text:
module HSA

include("LibHSARuntime.jl")

# Forward prefixed names
hsa_names = map(string, names(LibHSARuntime))
for name in filter(n->!(getproperty(LibHSARuntime, Symbol(n)) isa Function), hsa_names)
    newname = if endswith(name, "_s")
        continue
    elseif endswith(name, "_t")
        Symbol(replace(titlecase(replace(name[begin:end-2], "hsa_"=>"")), "_"=>"",
                                                                          "Amd"=>"AMD",
                                                                          "Isa"=>"ISA"))
    else
        Symbol(replace(name, "HSA_"=>""))
    end
    @eval const $newname = LibHSARuntime.$(Symbol(name))
end
for func in filter(n->startswith(n, "hsa_") && getproperty(LibHSARuntime, Symbol(n)) isa Function, hsa_names)
    newfunc = Symbol(replace(func, "hsa_"=>""))
    @eval const $newfunc = LibHSARuntime.$(Symbol(func))
end

end


text:
module PrecompileFuncs
    using FluxDash
    using CSV, DataFrames
    using FluxDash
    using FluxDash.Flux
    using Random
    function precompile_nn_train()
        df = CSV.read("FluxDash/datasets/real_estate.csv", DataFrame)

        n_row = nrow(df)
        df = df[shuffle(1:n_row)[:], :] #shuffle data

        labels = names(df)
        in_labels = labels[end]
        out_labels = labels[end]

        X_train, X_test, y_train, y_test = FluxDash.format_nn_data(
            df, in_labels, out_labels; training_percent=75
        )

        chn = Base.Channel{Tuple{Float64, Int, Chain}}(Inf)

        hidden_outs = [2]
        hidden_activations =["relu"]
        FluxDash.train_nn(
            X_train, y_train, in_labels, out_labels, hidden_outs,
            chn, hidden_activations; training_percent = 75,
            epoc = 4, opt = "adam", eta = 0.1
        )
        return "success"
    end
end

max_stars_repo_name: DEICIDIA/fluxdash

max_issues_repo_name: DEICIDIA/fluxdash

max_forks_repo_name: DEICIDIA/fluxdash


text:
@testset "arithmetic operations" begin
    x = RevisedAffineForm(1, [1, 2, 3], 1)
    y = RevisedAffineForm(2, [3, 4, 5], 0.1)

    @test +x == x
    @test -x == RevisedAffineForm(-1, [-1, -2, -3], 1)

    @test x + y == RevisedAffineForm(3, [4, 6, 8], 1.1)
    @test x - y == RevisedAffineForm(-1, [-2, -2, -2], 1.1)

    @test x + 1 == RevisedAffineForm(2, [1, 2, 3], 1)
    @test x - 1 == RevisedAffineForm(0, [1, 2, 3], 1)
    @test 1 + x == RevisedAffineForm(2, [1, 2, 3], 1)
    @test 1 - x == RevisedAffineForm(0, [-1, -2, -3], 1)
    @test -2 * x == RevisedAffineForm(-2, [-2, -4, -6], 2)
    @test x * 2 == RevisedAffineForm(2, [2, 4, 6], 2)
    @test x / 2 == RevisedAffineForm(0.5, [0.5, 1, 1.5], 0.5)
end

@testset "correcty rounded sum and dot product" begin
    a = [Interval(0.1), Interval(0.2), Interval(0.3)]
    b = [Interval(0.4), Interval(0.6), Interval(0.7)]
    am = mid.(a)
    bm = mid.(b)

    si = sum(a)
    di = dot(a, b)
    @test RAA._sum_up(am) == sup(si)
    @test RAA._sum_down(am) == inf(si)
    @test RAA._dot_up(am, bm) == sup(di)
    @test RAA._dot_down(am, bm) == inf(di)

    @test RAA._sum_up(abs, -am) == sup(si)
    @test RAA._sum_down(abs, -am) == inf(si)
    @test RAA._dot_up(abs, -am, -bm) == sup(di)
    @test RAA._dot_down(abs, -am, -bm) == inf(di)
end


text:
### ndft functions ###

function ndft(plan::NFFTPlan{D}, f::AbstractArray{T,D}) where {D,T}
    plan.N == size(f) || throw(DimensionMismatch("Data is not consistent with NFFTPlan"))

    g = zeros(T, plan.M)

    for l=1:prod(plan.N)
        idx = CartesianIndices(plan.N)[l]

        for k=1:plan.M
            arg = zero(T)
            for d=1:D
                arg += plan.x[d,k] * ( idx[d] - 1 - plan.N[d] / 2 )
            end
            g[k] += f[l] * cis(-2*pi*arg)
        end
    end

    return g
end

function ndft(x, f::AbstractArray{T,D}, rest...; kwargs...) where {D,T}
    p = NFFTPlan(x, size(f), rest...; kwargs...)
    return ndft(p, f)
end

function ndft_adjoint(plan::NFFTPlan{D}, fHat::AbstractArray{T,1}) where {D,T}
    plan.M == length(fHat) || throw(DimensionMismatch("Data is not consistent with NFFTPlan"))

    g = zeros(T, plan.N)

    for l=1:prod(plan.N)
        idx = CartesianIndices(plan.N)[l]

        for k=1:plan.M
            arg = zero(T)
            for d=1:D
                arg += plan.x[d,k] * ( idx[d] - 1 - plan.N[d] / 2 )
            end
            g[l] += fHat[k] * cis(2*pi*arg)
        end
    end

    return g
end

function ndft_adjoint(x, N, fHat::AbstractVector{T}, rest...; kwargs...) where T
    p = NFFTPlan(x, N, rest...; kwargs...)
    return ndft_adjoint(p, fHat)
end