text:
"""
    DataFrameRow{<:AbstractDataFrame,<:AbstractIndex}

A view of one row of an `AbstractDataFrame`.

A `DataFrameRow` is returned by `getindex` or `view` functions when one row and a
selection of columns are requested, or when iterating the result
of the call to the [`eachrow`](@ref) function.

The `DataFrameRow` constructor can also be called directly:

```
DataFrameRow(parent::AbstractDataFrame, row::Integer, cols=:)
```

A `DataFrameRow` supports the iteration interface and can therefore be passed to
functions that expect a collection as an argument.

Indexing is one-dimensional like specifying a column of a `DataFrame`.
You can also access the data in a `DataFrameRow` using the `getproperty` and
`setproperty!` functions and convert it to a `NamedTuple` using the `copy` function.

It is possible to create a `DataFrameRow` with duplicate columns.
All such columns will have a reference to the same entry in the parent `DataFrame`.

If the selection of columns in a parent data frame is passed as `:` (a colon)
then `DataFrameRow` will always have all columns from the parent,
even if they are added or removed after its creation.

# Examples
```julia
df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
sdf1 = view(df, 2, :)
sdf2 = @view df[end, [:a]]
sdf3 = eachrow(df)[1]
sdf4 = DataFrameRow(df, 2, 1:2)
sdf5 = DataFrameRow(df, 1)
```
"""
struct DataFrameRow{D<:AbstractDataFrame,S<:AbstractIndex}
    df::D
    colindex::S
    row::Int

    @inline DataFrameRow(df::D, colindex::S, row::Union{Signed, Unsigned}) where
        {D<:AbstractDataFrame,S<:AbstractIndex} = new{D,S}(df, colindex, row)
end

Base.@propagate_inbounds function DataFrameRow(df::DataFrame, row::Integer, cols)
    @boundscheck if !checkindex(Bool, axes(df, 1), row)
        throw(BoundsError("attempt to access a data frame with $(nrow(df)) " *
                          "rows at index $row"))
    end
    DataFrameRow(df, SubIndex(index(df), cols), row)
end

Base.@propagate_inbounds DataFrameRow(df::DataFrame, row::Bool, cols) =
    throw(ArgumentError("invalid row index of type Bool"))

Base.@propagate_inbounds function DataFrameRow(sdf::SubDataFrame, row::Integer, cols)
    @boundscheck if !checkindex(Bool, axes(sdf, 1), row)
        throw(BoundsError("attempt to access a data frame with $(nrow(sdf)) " *
                          "rows at index $row"))
    end
    if index(sdf) isa Index # sdf was created using : as row selector
        colindex = SubIndex(index(sdf), cols)
    else
        colindex = SubIndex(index(parent(sdf)), parentcols(index(sdf), cols))
    end
    @inbounds DataFrameRow(parent(sdf), colindex, rows(sdf)[row])
end

Base.@propagate_inbounds DataFrameRow(df::SubDataFrame, row::Bool, cols) =
    throw(ArgumentError("invalid row index of type Bool"))

Base.@propagate_inbounds DataFrameRow(df::AbstractDataFrame, row::Integer) =
    DataFrameRow(df, row, :)

row(r::DataFrameRow) = getfield(r, :row)
Base.parent(r::DataFrameRow) = getfield(r, :df)
Base.parentindices(r::DataFrameRow) = (row(r), parentcols(index(r)))

Base.@propagate_inbounds Base.view(adf::AbstractDataFrame, rowind::Integer,
                                   colinds::Union{Colon, AbstractVector, Regex, Not, Between, All}) =
    DataFrameRow(adf, rowind, colinds)

Base.@propagate_inbounds Base.getindex(df::AbstractDataFrame, rowind::Integer,
                                       colinds::Union{AbstractVector, Regex, Not, Between, All}) =
    DataFrameRow(df, rowind, colinds)
Base.@propagate_inbounds Base.getindex(df::AbstractDataFrame, rowind::Integer, ::Colon) =
    DataFrameRow(df, rowind, :)
Base.@propagate_inbounds Base.getindex(r::DataFrameRow, idx::ColumnIndex) =
    parent(r)[row(r), parentcols(index(r), idx)]
Base.@propagate_inbounds Base.getindex(r::DataFrameRow, idxs::Union{AbstractVector, Regex, Not, Between, All}) =
    DataFrameRow(parent(r), row(r), parentcols(index(r), idxs))
Base.@propagate_inbounds Base.getindex(r::DataFrameRow, ::Colon) = r

for T in (:AbstractVector, :Regex, :Not, :Between, :All, :Colon)
    @eval function Base.setindex!(df::DataFrame,
                                  v::Union{DataFrameRow, NamedTuple, AbstractDict},
                                  row_ind::Integer,
                                  col_inds::$(T))
        idxs = index(df)[col_inds]
        if length(v) != length(idxs)
            throw(DimensionMismatch("$(length(idxs)) columns were selected but the assigned" *
                                    " collection contains $(length(v)) elements"))
        end

        if v isa AbstractDict
            for n in view(_names(df), idxs)
                if !haskey(v, n)
                    throw(ArgumentError("Column :$n not found in source dictionary"))
                end
            end
        elseif !all(((a, b),) -> a == b, zip(view(_names(df), idxs), keys(v)))
            mismatched = findall(view(_names(df), idxs) .!= collect(keys(v)))
            throw(ArgumentError("Selected column names do not match the names in assigned value in" *
                                " positions $(join(mismatched, ", ", " and "))"))
        end

        for (col, val) in pairs(v)
            df[row_ind, col] = val
        end
        return df
    end
end

Base.@propagate_inbounds Base.setindex!(r::DataFrameRow, value, idx) =
    setindex!(parent(r), value, row(r), parentcols(index(r), idx))

index(r::DataFrameRow) = getfield(r, :colindex)

Base.names(r::DataFrameRow) = _names(parent(r))[parentcols(index(r), :)]
_names(r::DataFrameRow) = view(_names(parent(r)), parentcols(index(r), :))

Base.haskey(r::DataFrameRow, key::Bool) =
    throw(ArgumentError("invalid key: $key of type Bool"))
Base.haskey(r::DataFrameRow, key::Integer) = 1 ≤ key ≤ size(r, 1)
function Base.haskey(r::DataFrameRow, key::Symbol)
    hasproperty(parent(r), key) || return false
    index(r) isa Index && return true
    # here index(r) is a SubIndex
    pos = index(parent(r))[key]
    remap = index(r).remap
    length(remap) == 0 && lazyremap!(index(r))
    checkbounds(Bool, remap, pos) || return false
    remap[pos] > 0
end

Base.getproperty(r::DataFrameRow, idx::Symbol) = getindex(r, idx)
Base.setproperty!(r::DataFrameRow, idx::Symbol, x::Any) = setindex!(r, x, idx)
# Private fields are never exposed since they can conflict with column names
Base.propertynames(r::DataFrameRow, private::Bool=false) = names(r)

Base.view(r::DataFrameRow, col::ColumnIndex) =
    view(parent(r)[!, parentcols(index(r), col)], row(r))
Base.view(r::DataFrameRow, cols::Union{AbstractVector, Regex, Not, Between, All}) =
    DataFrameRow(parent(r), row(r), parentcols(index(r), cols))
Base.view(r::DataFrameRow, ::Colon) = r

Base.size(r::DataFrameRow) = (length(index(r)),)
Base.size(r::DataFrameRow, i) = size(r)[i]
Base.length(r::DataFrameRow) = size(r, 1)
Base.ndims(r::DataFrameRow) = 1
Base.ndims(::Type{<:DataFrameRow}) = 1

Base.lastindex(r::DataFrameRow) = length(r)

Base.iterate(r::DataFrameRow) = iterate(r, 1)

function Base.iterate(r::DataFrameRow, st)
    st > length(r) && return nothing
    return (r[st], st + 1)
end

# Computing the element type requires going over all columns,
# so better let collect() do it only if necessary (widening)
Base.IteratorEltype(::Type{<:DataFrameRow}) = Base.EltypeUnknown()

function Base.convert(::Type{Vector}, dfr::DataFrameRow)
    df = parent(dfr)
    T = reduce(promote_type, (eltype(df[!, i]) for i in parentcols(index(dfr))))
    convert(Vector{T}, dfr)
end
Base.convert(::Type{Vector{T}}, dfr::DataFrameRow) where T =
    T[dfr[i] for i in 1:length(dfr)]
Base.Vector(dfr::DataFrameRow) = convert(Vector, dfr)
Base.Vector{T}(dfr::DataFrameRow) where T = convert(Vector{T}, dfr)

Base.keys(r::DataFrameRow) = Tuple(_names(r))
Base.values(r::DataFrameRow) =
    ntuple(col -> parent(r)[row(r), parentcols(index(r), col)], length(r))
Base.map(f, r::DataFrameRow, rs::DataFrameRow...) = map(f, copy(r), copy.(rs)...)
Base.get(dfr::DataFrameRow, key::ColumnIndex, default) =
    haskey(dfr, key) ? dfr[key] : default
Base.get(f::Base.Callable, dfr::DataFrameRow, key::ColumnIndex) =
    haskey(dfr, key) ? dfr[key] : f()
Base.broadcastable(::DataFrameRow) =
    throw(ArgumentError("broadcasting over `DataFrameRow`s is reserved"))

"""
    copy(dfr::DataFrameRow)

Convert a [`DataFrameRow`](@ref) to a `NamedTuple`.
"""
Base.copy(r::DataFrameRow) = NamedTuple{Tuple(keys(r))}(values(r))

# hash column element
Base.@propagate_inbounds hash_colel(v::AbstractArray, i, h::UInt = zero(UInt)) =
    hash(v[i], h)
Base.@propagate_inbounds function hash_colel(v::AbstractCategoricalArray, i,
                                             h::UInt = zero(UInt))
    ref = v.refs[i]
    if eltype(v) >: Missing && ref == 0
        hash(missing, h)
    else
        hash(CategoricalArrays.index(v.pool)[ref], h)
    end
end

# hash of DataFrame rows based on its values
# so that duplicate rows would have the same hash
# table columns are passed as a tuple of vectors to ensure type specialization
rowhash(cols::Tuple{AbstractVector}, r::Int, h::UInt = zero(UInt))::UInt =
    hash_colel(cols[1], r, h)
function rowhash(cols::Tuple{Vararg{AbstractVector}}, r::Int, h::UInt = zero(UInt))::UInt
    h = hash_colel(cols[1], r, h)
    rowhash(Base.tail(cols), r, h)
end

Base.hash(r::DataFrameRow, h::UInt = zero(UInt)) =
    rowhash(ntuple(col -> parent(r)[!, parentcols(index(r), col)], length(r)), row(r), h)

function Base.:(==)(r1::DataFrameRow, r2::DataFrameRow)
    if parent(r1) === parent(r2)
        parentcols(index(r1)) == parentcols(index(r2)) || return false
        row(r1) == row(r2) && return true
    else
        _names(r1) == _names(r2) || return false
    end
    all(((a, b),) -> a == b, zip(r1, r2))
end

function Base.isequal(r1::DataFrameRow, r2::DataFrameRow)
    if parent(r1) === parent(r2)
        parentcols(index(r1)) == parentcols(index(r2)) || return false
        row(r1) == row(r2) && return true
    else
        _names(r1) == _names(r2) || return false
    end
    all(((a, b),) -> isequal(a, b), zip(r1, r2))
end

# lexicographic ordering on DataFrame rows, missing > !missing
function Base.isless(r1::DataFrameRow, r2::DataFrameRow)
    length(r1) == length(r2) ||
        throw(ArgumentError("compared DataFrameRows must have the same number " *
                            "of columns (got $(length(r1)) and $(length(r2)))"))
    for (a,b) in zip(r1, r2)
        isequal(a, b) || return isless(a, b)
    end
    return false
end

function DataFrame(dfr::DataFrameRow)
    row, cols = parentindices(dfr)
    parent(dfr)[row:row, cols]
end

@noinline pushhelper!(x, r) = push!(x, x[r])

function Base.push!(df::DataFrame, dfr::DataFrameRow; columns::Symbol=:equal)
    if !(columns in (:equal, :intersect))
        throw(ArgumentError("`columns` keyword argument must be `:equal` or `:intersect`"))
    end
    nrows, ncols = size(df)
    targetrows = nrows + 1
    if ncols == 0
        for (n, v) in pairs(dfr)
            setproperty!(df, n, fill!(Tables.allocatecolumn(typeof(v), 1), v))
        end
        return df
    end

    try
        if parent(dfr) === df && index(dfr) isa Index
            # in this case we are sure that all we do is safe
            r = row(dfr)
            for col in _columns(df)
                # use a barrier function to improve performance
                pushhelper!(col, r)
            end
        else
            # DataFrameRow can contain duplicate columns and we disallow this
            # corner case when push!-ing
            # Only check for equal lengths, as an error will be thrown below if some names don't match
            if columns === :equal
                msg = "Number of columns of `row` does not match `DataFrame` column count."
                ncols == length(dfr) || throw(ArgumentError(msg))
            end
            for (col, nm) in zip(_columns(df), _names(df))
                push!(col, dfr[nm])
            end
        end
        for col in _columns(df)
            @assert length(col) == targetrows
        end
    catch err
        for col in _columns(df)
            resize!(col, nrows)
        end
        rethrow(err)
    end
    df
end


text:
# A semi-separable matrix is defined by
#
# S = triu(uvᵀ, bu+1) + tril(pqᵀ, -bl-1),
#
# where bu, bl >= 0. u and v are n × ru matrices and p and q are n × rl matrices.
# See [1] for the definition.

# [1] Chandrasekaran and Gu, Fast and stable algorithms for banded plus
# Semiseparable systems of linear equations, SIMAX, 25 (2003), pp. 373-384.

## Constructors

struct SemiSeparableMatrix{T} <: AbstractMatrix{T}
    L::LowRankMatrix{T}
    U::LowRankMatrix{T}
    bl::Int
    bu::Int

    function SemiSeparableMatrix(L::LowRankMatrix{T}, U::LowRankMatrix{T}, bl, bu) where T
        Lm, Ln = size(L)
        Um, Un = size(U)
        @assert Um == Un == Lm == Ln && Un >= bu+1 && Lm >= bl+1
        new{T}(L, U, bl, bu)
    end
end

function Base.convert(::Type{Matrix}, S::SemiSeparableMatrix)
    return triu(Matrix(S.U), S.bu+1) + tril(Matrix(S.L), -S.bl-1)
end

size(S::SemiSeparableMatrix) = size(S.L)

function getindex(S::SemiSeparableMatrix{T}, k::Int, j::Int)  where T
    k-j ≥ S.bl && return S.L[k,j]
    j-k ≥ S.bu && return S.U[k,j]
    return zero(T)
end

ext: jl


text:
"""
    import_data(filename:: String ; transf = "DEFAULT")

extracts the desired data from a file *.EMTPFT.
# Arguments
- `transf = "DEFAULT"` : is the label for the transformation if Only 1 configuration was computed.
"""
function import_data(filename:: String ; transf = "DEFAULT")
    df            = DataFrame(readdlm(filename,' ')) # this line will have to be modified with DataFrames 0.11
    dfPabs,dfPrad = extract_data( df , transf)
    return dfPabs,dfPrad
end

"""
    import_data(filename:: Vector{String} ; transf = "DEFAULT")

puts everything in a single dataframe if given a list of files.
# Arguments
- `transf = "DEFAULT"` : is the label for the transformation if Only 1 configuration was computed.
"""
function import_data(filename:: Vector{String} ; transf = "DEFAULT")
    dfPabsv = DataFrame[]
    dfPradv = DataFrame[]
    for i = 1:length(filename)
        dfPabs_temp, dfPrad_temp = import_data(filename[i] ; transf=transf)
        dfPabsv = push!(dfPabsv,dfPabs_temp)
        dfPradv = push!(dfPradv,dfPrad_temp)
    end
    dfPabs = extract_data(dfPabsv , transf)
    dfPrad = extract_data(dfPradv , transf)
    return dfPabs,dfPrad
end

function extract_data(df :: DataFrame , transf)
    df1 = df[df[:x1] .== transf, :]
    nrows = size(df1, 1)
    ncols = size(df1, 2)

    Pabs  = zeros(Float64,Int(nrows/4),5)
    Prad  = zeros(Float64,Int(nrows/4),5)

    colnames    = [:x1 , :x2, :x3, :x4, :x5]
    newcolnames = [:Freq , :c11, :c12, :c21 ,:c22]

    count = 1
    for i = 1:4:nrows
        Pabs[count,1] = df1[i,2]
        Prad[count,1] = df1[i,2]
        for j=0:3
            Pabs[count,j+2]   = df1[i+j,4]
            Prad[count,j+2]   = df1[i+j,5]
        end
        count +=1
    end
    dfPabs = DataFrame(Pabs)
    dfPrad = DataFrame(Prad)

    rename!(dfPabs,colnames,newcolnames)
    rename!(dfPrad,colnames,newcolnames)

    return dfPabs, dfPrad
end

function extract_data(df :: Vector{DataFrame} , transf)
    df_temp = df[1]
    for i =2:length(df)
        append!(df_temp,df[i])
    end
    sort!(df_temp, cols = :Freq)
    return df_temp
end

function convert_to_array(dfPabs,dfPrad,columnname)
    dfPabs[:Freq] = w0.*dfPabs[:Freq]
    dfPrad[:Freq] = w0.*dfPrad[:Freq]
    wv   = dfPabs[:Freq]
    Prad = dfPrad[:,columnname]
    Pabs = dfPabs[:,columnname]
    return wv, Pabs, Prad
end

function write_data(flname,dict)
   for k in keys(dict)
        propname = k
        writetable(flname*"_"*propname*".dat", dict[k])
    end
end

"""
    simulation_data(filetype :: SIFlux,...)

Return and save (optional) the spectral flux and transfer function.

# Arguments
- `trans = "DEFAULT"` : the number of the transformation, trans = "DEFAULT" if no transformation is given
- `filename ::  Union{String,Vector{String}}` : can be single filename or an array of filenames for parallel computations
- `columnname` : column name. Example : [:c12]. Only use one name at the time.
"""
function simulation_data(filetype :: SIFlux, filename :: Union{String,Vector{String}}, columnname ,T1,T2, trans = "DEFAULT"; savefile = (false," "))
    dfPabs,dfPrad  = import_data(filename ; transf = transf)
    wv, Pabs, Prad = convert_to_array(dfPabs,dfPrad,columnname)

    qtrans = zeros(Float64,length(wv))
    qtrans = transfer_w.(T1,T2,wv,Prad)
    Qtrans = [wv qtrans]

    if savefile[1] == true
        dfQtrans = DataFrame(Qtrans)
        dc = Dict("Prad" => dfPrad ,"Pabs" => dfPabs , "Qtrans" => dfQtrans)
        write_data(savefile[2], dc)
    end
    return wv, Pabs, Prad, qtrans
end

"""
    simulation_data(filetype :: TotalFlux{:vsT},...)

Return and save (optional) total heat_transfer as a function temparature of one of the bodies and transfer function.
# Arguments
- `Tmin` : minimum value of temperature of body 2
- `Tmax` : maximum value of temperature of body 2
- `T1 = 0.0` : temperature of body 1. Set by default at 0.0K
"""
function simulation_data(filetype :: TotalFlux{:vsT}, filename:: Union{String,Vector{String}}, columnname,Tmin,Tmax, trans= "DEFAULT";T1=0.0, savefile = (false," "))
    dfPabs,dfPrad  = import_data(filename ; transf = transf)
    wv, Pabs, Prad = convert_to_array(dfPabs,dfPrad,columnname)

    Tempv = collect(linspace(Tmin,Tmax,100))
    q_tot = zeros(Float64,100)
    tt(T) = total_transfer(T1,T,wv,Prad)
    q_tot = tt.(Tempv)
    Qtrans = [Tempv q_tot]

    if savefile[1] == true
        dfQtrans = DataFrame(Qtrans)
        dc = Dict("Prad" => dfPrad ,"Pabs" => dfPabs , "Qtrans_vs_T_" => dfQtrans)
        write_data(savefile[2], dc)
    end
    return wv, Pabs, Prad, Tempv, Qtrans
end

"""
    simulation_data(filetype :: TotalFlux{:vsd},...)

Return and save (optional) total heat_transfer as a function separation distance and transfer function for each separation distance.
# Arguments
- `T1` : Temperature of body 1
- `T2` : Temperature of body 2
- `trans` : vector containing the transformation label
"""
function simulation_data(filetype :: TotalFlux{:vsd}, filename :: Union{String,Vector{String}}, columnname ,T1,T2, trans; savefile = (false," "))
    q_tot   = zeros(Float64,length(collect(trans)))

    for i in trans
        dfPabs,dfPrad  = import_data(filename ; transf = Float64(trans[i]))
        wv, Pabs, Prad = convert_to_array(dfPabs,dfPrad,columnname)
        q_tot[i]       = total_transfer(T1,T2,wv,Prad)
        if savefile[1] == true
            dc = Dict("Prad_trans="*"_$i" => dfPrad ,"Pabs_trans="*"_$i" => dfPabs)
            write_data(savefile[2], dc)
        end
    end

    Qtrans = [trans q_tot]
    if savefile[1] == true
        dfQtrans = DataFrame(Qtrans)
        dc = Dict("Qtrans_vs_dist_" => dfQtrans)
        write_data(savefile[2], dc)
    end
    return trans, q_tot
end

"""
    simulation_data(filetype :: TotalFlux{:vsth},...)

Return and save (optional) total heat_transfer as a function cylinder thickness.
# Arguments
- `T1` : temperature of body 1
- `T2` : temperature of body 2
"""
function simulation_data(filetype :: TotalFlux{:vsth}, filename:: Vector{String}, columnname,T1,T2, trans= "DEFAULT"; savefile = (false," "))

    lengthlist = length(filename)
    q_tot      = zeros(Float64,lengthlist)
    thv        = zeros(Float64,lengthlist)

    for i in 1:lengthlist
        dfPabs,dfPrad  = import_data(filename[i] ; transf = trans)
        wv, Pabs, Prad = convert_to_array(dfPabs,dfPrad,columnname)
        q_tot[i]       = total_transfer(T1,T2,wv,Prad)
        indexb         = searchindex(filename[i],"=")
        indexe         = searchindex(filename[i],"microns")
        thv[i]         = parse(Float64,filename[i][indexb+1:indexe-1])

        if savefile[1] == true
            dc = Dict("Prad_th="*"_$i" => dfPrad ,"Pabs_th="*"_$i" => dfPabs)
            write_data(savefile[2], dc)
        end

    end

    Qtrans = [thv q_tot]

    if savefile[1] == true
        dfQtrans = DataFrame(Qtrans)
        dc = Dict("Qtrans_vs_th_" => dfQtrans)
        write_data(savefile[2], dc)
    end
    return thv, q_tot
end


text:
using HDF5
using DocOpt
using HPAT

function generate_files(path,N)
    D = 10
    A = rand(Float32,D+1,N)
    hdf5_file = path*"logistic_regression.hdf5"
    csv_file = path*"logistic_regression.csv"
    if isfile(hdf5_file)
        rm(hdf5_file)
    end
    if isfile(csv_file)
        rm(csv_file)
    end
    h5write(hdf5_file,"/points", A[2:end,:])
    h5write(hdf5_file,"/responses", reshape(A[1,:],N))
    writecsv(csv_file,A')
end

function main()
    doc = """generate_logistic_regression.jl

generate random input for logistic regression example.

Usage:
  generate_logistic_regression.jl -h | --help
  generate_logistic_regression.jl [--instances=<instances>] [--path=<path>]

Options:
  -h --help                  Show this screen.
  --instances=<instances>    Specify number of instances; defaults to 2000000.
  --path=<path>              Specify output path for generated files; defaults to HPAT's default data path.
"""

    arguments = docopt(doc)

    if (arguments["--instances"] != nothing)
        instances = parse(Int, arguments["--instances"])
    else
        instances = 2*10^6
    end

    if (arguments["--path"] != nothing)
        path = arguments["--path"]
    else
        path = HPAT.getDefaultDataPath()
    end

    generate_files(path, instances)
end

main()


text:
function get_transform(obj::MODProvider, arg0::AbsoluteDate)
    return jcall(obj, "getTransform", Transform, (AbsoluteDate,), arg0)
end

function get_transform(obj::MODProvider, arg0::FieldAbsoluteDate)
    return jcall(obj, "getTransform", FieldTransform, (FieldAbsoluteDate,), arg0)
end


text:
module Nodes

import Base: show, length

export Node

# NOTE: we define an abstract type here to allow for mutually recursive types
#       it is expected that T = block for all purposes within the code
#       this is solely to deal with a limitation of julia
# NOTE: we make this mutable so that independent calls to the constructor
#       with equivalent parameters ensures the return of a brand new object
"""
	mutable struct Node{T}
		block  :: T
		strand :: Bool
	end

Node represents a portion of a sequence path that passes through a single `block`.
`strand` stores whether we pass along the forward strand of `block` (if true) or reverse (if false).
"""
mutable struct Node{T}
    block  :: T
    strand :: Bool
end

# --------------------------------
# constructors

"""
	Node{T}(b::T; strand=true)

Create a Node that passed through block `b`.
Default to forward strand orientation.
"""
Node(b::T; strand=true) where T = Node{T}(b,strand)

# --------------------------------
# operators

# Base.show(io::IO, n::Node) = Base.show(io, UInt64(pointer_from_objref(n)))
Base.show(io::IO, n::Node) = Base.show(io,
    (block = n.block, strand = n.strand)
)

"""
	length(n::Node) = length(n.block, n)

Return the length of sequence stored within node `n`
"""
length(n::Node) = length(n.block, n)

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
# A SortedMultiDict is a wrapper around balancedTree.
## Unlike SortedDict, a key in SortedMultiDict can
## refer to multiple data entries.

mutable struct SortedMultiDict{K, D, Ord <: Ordering}
    bt::BalancedTree23{K,D,Ord}

    ## Base constructors

    """
        SortedMultiDict{K,V,Ord}(o)

    Construct an empty sorted multidict in which type parameters are
    explicitly listed; ordering object is explicitly specified. (See
    below for discussion of ordering.) An empty SortedMultiDict may also
    be constructed via `SortedMultiDict(K[], V[], o)` where the `o`
    argument is optional.
    """
    SortedMultiDict{K,D,Ord}(o::Ord) where {K,D,Ord} = new{K,D,Ord}(BalancedTree23{K,D,Ord}(o))
    function SortedMultiDict{K,D,Ord}(o::Ord, kv) where {K,D,Ord}
        smd = new{K,D,Ord}(BalancedTree23{K,D,Ord}(o))

        if eltype(kv) <: Pair
            # It's (possibly?) more efficient to access the first and second
            # elements of Pairs directly, rather than destructure
            for p in kv
                insert!(smd, p.first, p.second)
            end
        else
            for (k,v) in kv
                insert!(smd, k, v)
            end
        end
        return smd

    end
end

"""
    SortedMultiDict()

Construct an empty `SortedMultiDict` with key type `Any` and value type
`Any`. Ordering defaults to `Forward` ordering.

**Note that a key type of `Any` or any other abstract type will lead
to slow performance.**
"""
SortedMultiDict() = SortedMultiDict{Any,Any,ForwardOrdering}(Forward)

"""
    SortedMultiDict(o)

Construct an empty `SortedMultiDict` with key type `Any` and value type
`Any`, ordered using `o`.

**Note that a key type of `Any` or any other abstract type will lead
to slow performance.**
"""
SortedMultiDict(o::O) where {O<:Ordering} = SortedMultiDict{Any,Any,O}(o)

# Construction from Pairs
"""
    SortedMultiDict(k1=>v1, k2=>v2, ...)

Arguments are key-value pairs for insertion into the multidict. The
keys must be of the same type as one another; the values must also
be of one type.
"""
SortedMultiDict(ps::Pair...) = SortedMultiDict(Forward, ps)

"""
    SortedMultiDict(o, k1=>v1, k2=>v2, ...)

The first argument `o` is an ordering object. The remaining
arguments are key-value pairs for insertion into the multidict. The
keys must be of the same type as one another; the values must also
be of one type.
"""
SortedMultiDict(o::Ordering, ps::Pair...) = SortedMultiDict(o, ps)
SortedMultiDict{K,D}(ps::Pair...) where {K,D} = SortedMultiDict{K,D,ForwardOrdering}(Forward, ps)
SortedMultiDict{K,D}(o::Ord, ps::Pair...) where {K,D,Ord<:Ordering} = SortedMultiDict{K,D,Ord}(o, ps)

# Construction from AbstractDicts
SortedMultiDict(o::Ord, d::AbstractDict{K,D}) where {K,D,Ord<:Ordering} = SortedMultiDict{K,D,Ord}(o, d)

## Construction from iteratables of Pairs/Tuples

# Construction specifying Key/Value types
# e.g., SortedMultiDict{Int,Float64}([1=>1, 2=>2.0])
"""
    SortedMultiDict{K,D}(iter)

Takes an arbitrary iterable object of key=>value pairs with
key type `K` and value type `D`. The default Forward ordering is used.
"""
SortedMultiDict{K,D}(kv) where {K,D} = SortedMultiDict{K,D}(Forward, kv)

"""
    SortedMultiDict{K,D}(o, iter)

Takes an arbitrary iterable object of key=>value pairs with
key type `K` and value type `D`. The ordering object `o` is explicitly given.
"""
function SortedMultiDict{K,D}(o::Ord, kv) where {K,D,Ord<:Ordering}
    try
        SortedMultiDict{K,D,Ord}(o, kv)
    catch e
        if not_iterator_of_pairs(kv)
            throw(ArgumentError("SortedMultiDict(kv): kv needs to be an iterator of tuples or pairs"))
        else
            rethrow(e)
        end
    end
end

# Construction inferring Key/Value types from input
# e.g. SortedMultiDict{}

SortedMultiDict(o1::Ordering, o2::Ordering) = throw(ArgumentError("SortedMultiDict with two parameters must be called with an Ordering and an interable of pairs"))
SortedMultiDict(kv, o::Ordering=Forward) = SortedMultiDict(o, kv)
function SortedMultiDict(o::Ordering, kv)
    try
        _sorted_multidict_with_eltype(o, kv, eltype(kv))
    catch e
        if not_iterator_of_pairs(kv)
            throw(ArgumentError("SortedMultiDict(kv): kv needs to be an iterator of tuples or pairs"))
        else
            rethrow(e)
        end
    end
end

_sorted_multidict_with_eltype(o::Ord, ps, ::Type{Pair{K,D}}) where {K,D,Ord} = SortedMultiDict{  K,  D,Ord}(o, ps)
_sorted_multidict_with_eltype(o::Ord, kv, ::Type{Tuple{K,D}}) where {K,D,Ord} = SortedMultiDict{  K,  D,Ord}(o, kv)
_sorted_multidict_with_eltype(o::Ord, ps, ::Type{Pair{K}}  ) where {K,  Ord} = SortedMultiDict{  K,Any,Ord}(o, ps)
_sorted_multidict_with_eltype(o::Ord, kv, ::Type            ) where {    Ord} = SortedMultiDict{Any,Any,Ord}(o, kv)

## TODO: It seems impossible (or at least very challenging) to create the eltype below.
##       If deemed possible, please create a test and uncomment this definition.
# _sorted_multi_dict_with_eltype{  D,Ord}(o::Ord, ps, ::Type{Pair{K,D} where K}) = SortedMultiDict{Any,  D,Ord}(o, ps)

const SMDSemiToken = IntSemiToken

const SMDToken = Tuple{SortedMultiDict, IntSemiToken}

## This function inserts an item into the tree.
## It returns a token that
## points to the newly inserted item.

"""
    insert!(sc, k)

Argument `sc` is a SortedDict or SortedMultiDict, `k` is a key and
`v` is the corresponding value. This inserts the `(k,v)` pair into
the container. If the key is already present in a SortedDict, this
overwrites the old value. In the case of SortedMultiDict, no
overwriting takes place (since SortedMultiDict allows the same key
to associate with multiple values). In the case of SortedDict, the
return value is a pair whose first entry is boolean and indicates
whether the insertion was new (i.e., the key was not previously
present) and the second entry is the semitoken of the new entry. In
the case of SortedMultiDict, a semitoken is returned (but no
boolean). Time: O(*c* log *n*)
"""
@inline function insert!(m::SortedMultiDict{K,D,Ord}, k_, d_) where {K, D, Ord <: Ordering}
    b, i = insert!(m.bt, convert(K,k_), convert(D,d_), true)
    IntSemiToken(i)
end

## push! is an alternative to insert!; it returns the container.
"""
    push!(sc, k=>v)

Argument `sc` is a SortedDict or SortedMultiDict and `k=>v` is a
key-value pair. This inserts the key-value pair into the container.
If the key is already present, this overwrites the old value. The
return value is `sc`. Time: O(*c* log *n*)
"""
@inline function push!(m::SortedMultiDict{K,D}, pr::Pair) where {K,D}
    insert!(m.bt, convert(K,pr[1]), convert(D,pr[2]), true)
    m
end

## First and last return the first and last (key,data) pairs
## in the SortedMultiDict.  It is an error to invoke them on an
## empty SortedMultiDict.

"""
    first(sc)

Argument `sc` is a SortedDict, SortedMultiDict or SortedSet. This
function returns the first item (a `k=>v` pair for SortedDict and
SortedMultiDict or a key for SortedSet) according to the sorted
order in the container. Thus, `first(sc)` is equivalent to
`deref((sc,startof(sc)))`. It is an error to call this function on
an empty container. Time: O(log *n*)
"""
@inline function first(m::SortedMultiDict)
    i = beginloc(m.bt)
    i == 2 && throw(BoundsError())
    return Pair(m.bt.data[i].k, m.bt.data[i].d)
end

"""
    last(sc)

Argument `sc` is a SortedDict, SortedMultiDict or SortedSet. This
function returns the last item (a `k=>v` pair for SortedDict and
SortedMultiDict or a key for SortedSet) according to the sorted
order in the container. Thus, `last(sc)` is equivalent to
`deref((sc,lastindex(sc)))`. It is an error to call this function on an
empty container. Time: O(log *n*)
"""
@inline function last(m::SortedMultiDict)
    i = endloc(m.bt)
    i == 1 && throw(BoundsError())
    return Pair(m.bt.data[i].k, m.bt.data[i].d)
end

function searchequalrange(m::SortedMultiDict, k_)
    k = convert(keytype(m),k_)
    i1 = findkeyless(m.bt, k)
    i2, exactfound = findkey(m.bt, k)
    if exactfound
        i1a = nextloc0(m.bt, i1)
        return IntSemiToken(i1a), IntSemiToken(i2)
    else
        return IntSemiToken(2), IntSemiToken(1)
    end
end

## '(k,d) in m' checks whether a key-data pair is in
## a sorted multidict.  This requires a loop over
## all data items whose key is equal to k.

function in_(k_, d_, m::SortedMultiDict)
    k = convert(keytype(m), k_)
    d = convert(valtype(m), d_)
    i1 = findkeyless(m.bt, k)
    i2,exactfound = findkey(m.bt,k)
    !exactfound && return false
    ord = m.bt.ord
    while true
        i1 = nextloc0(m.bt, i1)
        @assert(eq(ord, m.bt.data[i1].k, k))
        m.bt.data[i1].d == d && return true
        i1 == i2 && return false
    end
end

"""
    eltype(sc)

Returns the (key,value) type (a 2-entry pair, i.e., `Pair{K,V}`) for
SortedDict and SortedMultiDict. Returns the key type for SortedSet.
This function may also be applied to the type itself. Time: O(1)
"""
@inline eltype(m::SortedMultiDict{K,D,Ord}) where {K,D,Ord <: Ordering} =  Pair{K,D}
@inline eltype(::Type{SortedMultiDict{K,D,Ord}}) where {K,D,Ord <: Ordering} =  Pair{K,D}

"""
    in(p, sc)

Returns true if `p` is in `sc`. In the case that `sc` is a
SortedDict or SortedMultiDict, `p` is a key=>value pair. In the
case that `sc` is a SortedSet, `p` should be a key. Time: O(*c* log
*n* + *d*) for SortedDict and SortedSet, where *d* stands for the
time to compare two values. In the case of SortedMultiDict, the time
is O(*c* log *n* + *dl*), and *l* stands for the number of entries
that have the key of the given pair. (So therefore this call is
inefficient if the same key addresses a large number of values, and
an alternative should be considered.)
"""
@inline in(pr::Pair, m::SortedMultiDict) =
    in_(pr[1], pr[2], m)
@inline in(::Tuple{Any,Any}, ::SortedMultiDict) =
    throw(ArgumentError("'(k,v) in sortedmultidict' not supported in Julia 0.4 or 0.5.  See documentation"))

"""
    keytype(sc)

Returns the key type for SortedDict, SortedMultiDict and SortedSet.
This function may also be applied to the type itself. Time: O(1)
"""
@inline keytype(m::SortedMultiDict{K,D,Ord}) where {K,D,Ord <: Ordering} = K
@inline keytype(::Type{SortedMultiDict{K,D,Ord}}) where {K,D,Ord <: Ordering} = K

"""
    valtype(sc)

Returns the value type for SortedDict and SortedMultiDict. This
function may also be applied to the type itself. Time: O(1)
"""
@inline valtype(m::SortedMultiDict{K,D,Ord}) where {K,D,Ord <: Ordering} = D
@inline valtype(::Type{SortedMultiDict{K,D,Ord}}) where {K,D,Ord <: Ordering} = D

"""
    ordtype(sc)

Returns the order type for SortedDict, SortedMultiDict and
SortedSet. This function may also be applied to the type itself.
Time: O(1)
"""
@inline ordtype(m::SortedMultiDict{K,D,Ord}) where {K,D,Ord <: Ordering} = Ord
@inline ordtype(::Type{SortedMultiDict{K,D,Ord}}) where {K,D,Ord <: Ordering} = Ord

"""
    orderobject(sc)

Returns the order object used to construct the container. Time: O(1)
"""
@inline orderobject(m::SortedMultiDict) = m.bt.ord

"""
    haskey(sc,k)

Returns true if key `k` is present for SortedDict, SortedMultiDict
or SortedSet `sc`. For SortedSet, `haskey(sc,k)` is a synonym for
`in(k,sc)`. For SortedDict and SortedMultiDict, `haskey(sc,k)` is
equivalent to `in(k,keys(sc))`. Time: O(*c* log *n*)
"""
@inline function haskey(m::SortedMultiDict, k_)
    i, exactfound = findkey(m.bt,convert(keytype(m),k_))
    exactfound
end

## Check if two SortedMultiDicts are equal in the sense of containing
## the same (K,D) pairs in the same order.  This sense of equality does not mean
## that semitokens valid for one are also valid for the other.
"""
    isequal(sc1,sc2)

Checks if two containers are equal in the sense that they contain
the same items; the keys are compared using the `eq` method, while
the values are compared with the `isequal` function. In the case of
SortedMultiDict, equality requires that the values associated with a
particular key have same order (that is, the same insertion order).
Note that `isequal` in this sense does not imply any correspondence
between semitokens for items in `sc1` with those for `sc2`. If the
equality-testing method associated with the keys and values implies
hash-equivalence in the case of SortedDict, then `isequal` of the
entire containers implies hash-equivalence of the containers. Time:
O(*cn* + *n* log *n*)
"""
function isequal(m1::SortedMultiDict, m2::SortedMultiDict)
    ord = orderobject(m1)
    if !isequal(ord, orderobject(m2)) || !isequal(eltype(m1), eltype(m2))
        throw(ArgumentError("Cannot use isequal for two SortedMultiDicts unless their element types and ordering objects are equal"))
    end
    p1 = startof(m1)
    p2 = startof(m2)
    while true
        if p1 == pastendsemitoken(m1)
            return p2 == pastendsemitoken(m2)
        end
        if p2 == pastendsemitoken(m2)
            return false
        end
        k1,d1 = deref((m1,p1))
        k2,d2 = deref((m2,p2))
        if !eq(ord,k1,k2) || !isequal(d1,d2)
            return false
        end
        p1 = advance((m1,p1))
        p2 = advance((m2,p2))
    end
end

const SDorAbstractDict = Union{AbstractDict,SortedMultiDict}

function mergetwo!(m::SortedMultiDict{K,D,Ord},
                   m2::SDorAbstractDict) where {K,D,Ord <: Ordering}
    for (k,v) in m2
        insert!(m.bt, convert(K,k), convert(D,v), true)
    end
end

"""
    packcopy(sc)

This returns a copy of `sc` in which the data is packed. When
deletions take place, the previously allocated memory is not
returned. This function can be used to reclaim memory after many
deletions. Time: O(*cn* log *n*)
"""
function packcopy(m::SortedMultiDict{K,D,Ord}) where {K,D,Ord <: Ordering}
    w = SortedMultiDict{K,D}(orderobject(m))
    mergetwo!(w,m)
    w
end

"""
    packdeepcopy(sc)

This returns a packed copy of `sc` in which the keys and values are
deep-copied. This function can be used to reclaim memory after many
deletions. Time: O(*cn* log *n*)
"""
function packdeepcopy(m::SortedMultiDict{K,D,Ord}) where {K,D,Ord <: Ordering}
    w = SortedMultiDict{K,D}(orderobject(m))
    for (k,v) in m
        insert!(w.bt, deepcopy(k), deepcopy(v), true)
    end
    w
end

"""
    merge!(sc, sc1...)

This updates `sc` by merging SortedDicts or SortedMultiDicts `sc1`,
etc. into `sc`. These must all must have the same key-value types.
In the case of keys duplicated among the arguments, the rightmost
argument that owns the key gets its value stored for SortedDict. In
the case of SortedMultiDict all the key-value pairs are stored, and
for overlapping keys the ordering is left-to-right. This function is
not available for SortedSet, but the `union!` function (see below)
provides equivalent functionality. Time: O(*cN* log *N*), where *N*
is the total size of all the arguments.
"""
function merge!(m::SortedMultiDict{K,D,Ord},
                others::SDorAbstractDict...) where {K,D,Ord <: Ordering}
    for o in others
        mergetwo!(m,o)
    end
end

"""
    merge(sc1, sc2...)

This returns a SortedDict or SortedMultiDict that results from
merging SortedDicts or SortedMultiDicts `sc1`, `sc2`, etc., which
all must have the same key-value-ordering types. In the case of keys
duplicated among the arguments, the rightmost argument that owns the
key gets its value stored for SortedDict. In the case of
SortedMultiDict all the key-value pairs are stored, and for keys
shared between `sc1` and `sc2` the ordering is left-to-right. This
function is not available for SortedSet, but the `union` function
(see below) provides equivalent functionality. Time: O(*cN* log
*N*), where *N* is the total size of all the arguments.
"""
function merge(m::SortedMultiDict{K,D,Ord},
               others::SDorAbstractDict...) where {K,D,Ord <: Ordering}
    result = packcopy(m)
    merge!(result, others...)
    result
end

function Base.show(io::IO, m::SortedMultiDict{K,D,Ord}) where {K,D,Ord <: Ordering}
    print(io, "SortedMultiDict(")
    print(io, orderobject(m), ",")
    l = length(m)
    for (count,(k,v)) in enumerate(m)
        print(io, k, " => ", v)
        if count < l
            print(io, ", ")
        end
    end
    print(io, ")")
end

"""
    empty(sc)

Returns a new `SortedDict`, `SortedMultiDict`, or `SortedSet` of the same
type and with the same ordering as `sc` but with no entries (i.e.,
empty). Time: O(1)
"""
empty(m::SortedMultiDict{K,D,Ord}) where {K,D,Ord<:Ordering} =
   SortedMultiDict{K,D}(orderobject(m))
@deprecate similar(m::SortedMultiDict) empty(m)

isordered(::Type{T}) where {T<:SortedMultiDict} = true


text:
@testset "integration tests" begin
    julia_repo = "bcbi-test/julia"

    global_kwargs = Dict(
        :auth => BumpStdlibs.get_input_from_environment(:auth, "BUMPSTDLIBS_TOKEN_FOR_TESTS"),
        :close_old_pull_requests => true,
        :close_old_pull_requests_older_than => Dates.Minute(5),
        :pr_branch_suffix => Random.randstring(4),
        :push_if_no_changes => true,
        :stdlibs_to_include => "Pkg",
    )

    @testset "all stdlibs" begin
        @testset "push_if_no_changes = true" begin
            kwargs = deepcopy(global_kwargs)
            kwargs[:push_if_no_changes] = true
            kwargs[:stdlibs_to_include] = "all"
            result = bump_stdlibs(julia_repo; kwargs...)
            @test result isa Nothing
        end
        @testset "push_if_no_changes = false" begin
            kwargs = deepcopy(global_kwargs)
            kwargs[:push_if_no_changes] = false
            kwargs[:stdlibs_to_include] = "all"
            result = bump_stdlibs(julia_repo; kwargs...)
            @test result isa Nothing
        end
    end

    @testset "only specified stdlibs" begin
        kwargs = deepcopy(global_kwargs)
        result = bump_stdlibs(julia_repo; kwargs...)
        @test result isa Nothing
    end

    @testset "override the default branch" begin
        kwargs = deepcopy(global_kwargs)
        kwargs[:julia_repo_default_branch] = "master"
        result = bump_stdlibs(julia_repo; kwargs...)
        @test result isa Nothing
    end
end


text:
function derivative(A::LinearInterpolation{<:AbstractVector}, t::Number)
  idx = searchsortedfirst(A.t, t)
  if A.t[idx] >= t
    idx -= 1
  end
  idx == 0 ? idx += 1 : nothing
  θ = 1 / (A.t[idx+1] - A.t[idx])
  (A.u[idx+1] - A.u[idx]) / (A.t[idx+1] - A.t[idx])
end

function derivative(A::LinearInterpolation{<:AbstractMatrix}, t::Number)
  idx = searchsortedfirst(A.t, t)
  if A.t[idx] >= t
    idx -= 1
  end
  idx == 0 ? idx += 1 : nothing
  θ = 1 / (A.t[idx+1] - A.t[idx])
  @views @. (A.u[:, idx+1] - A.u[:, idx]) / (A.t[idx+1] - A.t[idx])
end

function derivative(A::QuadraticInterpolation{<:AbstractVector}, t::Number)
  idx = searchsortedfirst(A.t, t)
  if A.t[idx] >= t
    idx -= 1
  end
  idx == 0 ? idx += 1 : nothing
  if idx == length(A.t) - 1
    i₀ = idx - 1; i₁ = idx; i₂ = i₁ + 1;
  else
    i₀ = idx; i₁ = i₀ + 1; i₂ = i₁ + 1;
  end
  dl₀ = (2t - A.t[i₁] - A.t[i₂]) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂]))
  dl₁ = (2t - A.t[i₀] - A.t[i₂]) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂]))
  dl₂ = (2t - A.t[i₀] - A.t[i₁]) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁]))
  A.u[i₀] * dl₀ + A.u[i₁] * dl₁ + A.u[i₂] * dl₂
end

function derivative(A::QuadraticInterpolation{<:AbstractMatrix}, t::Number)
  idx = searchsortedfirst(A.t, t)
  if A.t[idx] >= t
      idx -= 1
  end
  idx == 0 ? idx += 1 : nothing
  if idx == length(A.t) - 1
    i₀ = idx - 1; i₁ = idx; i₂ = i₁ + 1;
  else
    i₀ = idx; i₁ = i₀ + 1; i₂ = i₁ + 1;
  end
  dl₀ = (2t - A.t[i₁] - A.t[i₂]) / ((A.t[i₀] - A.t[i₁]) * (A.t[i₀] - A.t[i₂]))
  dl₁ = (2t - A.t[i₀] - A.t[i₂]) / ((A.t[i₁] - A.t[i₀]) * (A.t[i₁] - A.t[i₂]))
  dl₂ = (2t - A.t[i₀] - A.t[i₁]) / ((A.t[i₂] - A.t[i₀]) * (A.t[i₂] - A.t[i₁]))
  @views @. A.u[:, i₀] * dl₀ + A.u[:, i₁] * dl₁ + A.u[:, i₂] * dl₂
end

function derivative(A::LagrangeInterpolation{<:AbstractVector}, t::Number)
  idxs = findRequiredIdxs(A, t)
  if A.t[idxs[1]] == t
    return zero(A.u[idxs[1]])
  end
  G = zero(A.u[1]); F = zero(A.t[1])
  DG = zero(A.u[1]); DF = zero(A.t[1])
  tmp = G
  for i = 1:length(idxs)
    if isnan(A.bcache[idxs[i]])
      mult = one(A.t[1])
      for j = 1:(i - 1)
        mult *= (A.t[idxs[i]] - A.t[idxs[j]])
      end
      for j = (i+1):length(idxs)
        mult *= (A.t[idxs[i]] - A.t[idxs[j]])
      end
      A.bcache[idxs[i]] = mult
    else
      mult = A.bcache[idxs[i]]
    end
    wi = inv(mult)
    tti = t - A.t[idxs[i]]
    tmp = wi / (t - A.t[idxs[i]])
    g = tmp * A.u[idxs[i]]
    G += g
    DG -= g / (t - A.t[idxs[i]])
    F += tmp
    DF -= tmp / (t - A.t[idxs[i]])
  end
  (DG * F - G * DF) / (F ^ 2)
end

function derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number)
  idxs = findRequiredIdxs(A, t)
  if A.t[idxs[1]] == t
    return zero(A.u[:, idxs[1]])
  end
  G = zero(A.u[:, 1]); F = zero(A.t[1])
  DG = zero(A.u[:, 1]); DF = zero(A.t[1])
  tmp = G
  for i = 1:length(idxs)
    if isnan(A.bcache[idxs[i]])
      mult = one(A.t[1])
      for j = 1:(i - 1)
        mult *= (A.t[idxs[i]] - A.t[idxs[j]])
      end
      for j = (i+1):length(idxs)
        mult *= (A.t[idxs[i]] - A.t[idxs[j]])
      end
      A.bcache[idxs[i]] = mult
    else
      mult = A.bcache[idxs[i]]
    end
    wi = inv(mult)
    tti = t - A.t[idxs[i]]
    tmp = wi / (t - A.t[idxs[i]])
    g = tmp * A.u[:, idxs[i]]
    @. G += g
    @. DG -= g / (t - A.t[idxs[i]])
    F += tmp
    DF -= tmp / (t - A.t[idxs[i]])
  end
  @. (DG * F - G * DF) / (F ^ 2)
end

function derivative(A::AkimaInterpolation{<:AbstractVector}, t::Number)
  i = searchsortedlast(A.t, t)
  i == 0 && return zero(A.u[1])
  i == length(A.t) && return zero(A.u[end])
  wj = t - A.t[i]
  @evalpoly wj A.b[i] 2A.c[i] 3A.d[i]
end

function derivative(A::ConstantInterpolation{<:AbstractVector}, t::Number)
  return isempty(searchsorted(A.t, t)) ? zero(A.u[1]) : eltype(A.u)(NaN)
end

function derivative(A::ConstantInterpolation{<:AbstractMatrix}, t::Number)
  return isempty(searchsorted(A.t, t)) ? zero(A.u[:, 1]) : eltype(A.u)(NaN) .* A.u[:, 1]
end

# QuadraticSpline Interpolation
function derivative(A::QuadraticSpline{<:AbstractVector{<:Number}}, t::Number)
  i = searchsortedfirst(A.t, t)
  i == 1 ? i += 1 : nothing
  σ = 1//2 * (A.z[i] - A.z[i - 1]) / (A.t[i] - A.t[i - 1])
  A.z[i-1] + 2σ * (t - A.t[i-1])
end

# CubicSpline Interpolation
function derivative(A::CubicSpline{<:AbstractVector{<:Number}}, t::Number)
  i = searchsortedfirst(A.t, t)
  isnothing(i) ? i = length(A.t) - 1 : i -= 1
  i == 0 ? i += 1 : nothing
  dI = -3A.z[i] * (A.t[i + 1] - t)^2 / (6A.h[i + 1]) + 3A.z[i + 1] * (t - A.t[i])^2 / (6A.h[i + 1])
  dC = A.u[i + 1] / A.h[i + 1] - A.z[i + 1] * A.h[i + 1] / 6
  dD = -(A.u[i] / A.h[i + 1] - A.z[i] * A.h[i + 1] / 6)
  dI + dC + dD
end

function derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number)
  # change t into param [0 1]
  idx = searchsortedlast(A.t,t)
  idx == length(A.t) ? idx -= 1 : nothing
  n = length(A.t)
  scale = (A.p[idx+1] - A.p[idx]) / (A.t[idx+1] - A.t[idx])
  t_ = A.p[idx] + (t - A.t[idx]) * scale
  N = DataInterpolations.spline_coefficients(n, A.d-1, A.k, t_)
  ducum = zero(eltype(A.u))
  for i = 1:(n - 1)
    ducum += N[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
  end
  ducum * A.d * scale
end

# BSpline Curve Approx
function derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number)
  # change t into param [0 1]
  idx = searchsortedlast(A.t,t)
  idx == 0 ? idx += 1 : nothing
  scale = (A.p[idx+1] - A.p[idx]) / (A.t[idx+1] - A.t[idx])
  t_ = A.p[idx] + (t - A.t[idx]) * scale
  N = spline_coefficients(A.h, A.d-1, A.k, t_)
  ducum = zero(eltype(A.u))
  for i = 1:(A.h - 1)
    ducum += N[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
  end
  ducum * A.d * scale
end


text:
using Test
using LinearAlgebra
using QCMaterial

import PyCall: pyimport
import Random

@testset "ucc.uccgsd" begin
    Random.seed!(1)
    scipy_opt = pyimport("scipy.optimize")

    nsite = 2
    n_qubit = 2*nsite
    U = 1.0
    t = -0.01
    μ = 0

    ham = FermionOperator()
    for i in 1:nsite
        up = up_index(i)
        down = down_index(i)
        ham += FermionOperator("$(up)^ $(down)^ $(up) $(down)", -U)
    end

    for i in 1:nsite-1
        ham += FermionOperator("$(up_index(i+1))^ $(up_index(i))", t)
        ham += FermionOperator("$(up_index(i))^ $(up_index(i+1))", t)
        ham += FermionOperator("$(down_index(i+1))^ $(down_index(i))", t)
        ham += FermionOperator("$(down_index(i))^ $(down_index(i+1))", t)
    end

    for i in 1:nsite
        up = up_index(i)
        down = down_index(i)
        ham += FermionOperator("$(up)^  $(up) ", -μ)
        ham += FermionOperator("$(down)^ $(down)", -μ)
    end

    n_electron = 2
    @assert mod(n_electron, 2) == 0
    sparse_mat = get_number_preserving_sparse_operator(ham, n_qubit, n_electron);

    enes_ed = eigvals(sparse_mat.toarray());

    ham_jw = jordan_wigner(ham)

    circuit = uccgsd(n_qubit, orbital_rot=true)

    function cost(theta_list)
        state = create_hf_state(n_qubit, n_electron)
        update_circuit_param!(circuit, theta_list)
        update_quantum_state!(circuit, state)
        get_expectation_value(ham_jw, state)
    end

    theta_init = rand(num_theta(circuit))
    cost_history = Float64[]
    init_theta_list = theta_init
    push!(cost_history, cost(init_theta_list))

    method = "BFGS"
    options = Dict("disp" => true, "maxiter" => 200, "gtol" => 1e-5)
    callback(x) = push!(cost_history, cost(x))
    opt = scipy_opt.minimize(cost, init_theta_list, method=method, callback=callback)

    EigVal_min = minimum(enes_ed)
    @test abs(EigVal_min-cost_history[end]) < 1e-6
end


text:
"""
    get_node(graph, label)

Returns the node corresponding to `label`.
"""
get_node(graph::DispatchGraph, node::T) where T<:DispatchNode = node

get_node(graph::DispatchGraph, label::T) where T<:AbstractString = begin
    found = Set{DispatchNode}()
    for node in nodes(graph)
        if has_label(node)
            get_label(node) == label && push!(found, node)
        end
    end
    length(found) > 1 && throw(ErrorException("Labels in dispatch graph are not unique."))
    length(found) < 1 && throw(ErrorException("No nodes with label $label found."))
    return pop!(found)
end

get_node(graph::DispatchGraph, node::T) where T = begin
    throw(ArgumentError("A node identifier can be either a " *
                        "::DispatchNode or ::AbstractString."))
end

"""
    load_hashchain(cachedir [; compression=DEFAULT_COMPRESSION])

Loads the hashchain file found in the directory `cachedir`. Before
loading, the `compression` value is checked against the one stored
in the hashchain file (both have to match). If the file does not exist,
it is created.
"""
function load_hashchain(cachedir::String=DEFAULT_CACHE_DIR;
                        compression::String=DEFAULT_COMPRESSION)
    cachedir = abspath(expanduser(cachedir))
    file = joinpath(cachedir, DEFAULT_HASHCHAIN_FILENAME)
    cachedir_outputs = joinpath(cachedir, DEFAULT_HASHCACHE_DIR)
    if !ispath(cachedir_outputs)
        @debug "Creating the cache directory..."
        mkpath(cachedir_outputs)
    end

    local hashchain
    if !isfile(file)
        @debug "Creating a new hashchain file $file..."
        hashchain = Dict{String, Any}()
        store_hashchain(hashchain, cachedir, compression=compression)
    else
        local data
        open(file, "r") do fid  # read the whole JSON hashchain file
            data = JSON.parse(read(fid, String))
        end
        if compression != data["compression"]
            throw(ErrorException("Compression mismatch: $compression vs. "*
                                 "$(data["compression"])"))
        end
        hashchain = data["hashchain"]
        # Clean up hashchain based on what exists already on disk
        # i.e. remove keys not found on disk
        on_disk_hashes = map(filename->split(filename, ".")[1],
                             filter!(!isfile, readdir(cachedir_outputs)))
        keys_to_delete = setdiff(keys(hashchain), on_disk_hashes)
        for key in keys_to_delete
            delete!(hashchain, key)
        end
        store_hashchain(hashchain, cachedir, compression=compression)
    end
    return hashchain
end

"""
    store_hashchain(hashchain, cachedir=DEFAULT_CACHE_DIR [; compression=DEFAULT_COMPRESSION, version=1])

Stores the `hashchain` object in a file named `DEFAULT_HASHCHAIN_FILENAME`,
in the directory `cachedir`. The values of `compression` and `version` are
stored as well in the file.
"""
function store_hashchain(hashchain::Dict{String, Any},
                         cachedir::String=DEFAULT_CACHE_DIR;
                         compression::String=DEFAULT_COMPRESSION,
                         version::Int=1)
    cachedir = abspath(expanduser(cachedir))
    if !ispath(cachedir)
        @debug "Creating the cache directory..."
        mkpath(cachedir)
    end
    file = joinpath(cachedir, DEFAULT_HASHCHAIN_FILENAME)
    hashchain = Dict("version" => version,
                     "compression" => compression,
                     "hashchain" => hashchain)
    open(file, "w+") do fid
        write(fid, JSON.json(hashchain, 4))
    end
end

"""
    get_compressor(compression, action)

Return a `TranscodingStreams` compatible compressor or decompressor
based on the values of `compression` and `action`.
"""
function get_compressor(compression::AbstractString, action::AbstractString)
    # Checks
    if !(compression in ["bz2", "bzip2", "gz", "gzip", "none"])
        throw(ErrorException("Unknown compression option,"*
                             " aborting."))
    end
    if !(action in ["compress", "decompress"])
        throw(ErrorException("The action can only be \"compress\" or \"decompress\"."))
    end
    # Get compressor/decompressor
    if compression == "bz2" || compression == "bzip2"
        compressor = ifelse(action == "compress",
                            Bzip2CompressorStream,
                            Bzip2DecompressorStream)
    elseif compression == "gz" || compression == "gzip"
        compressor = ifelse(action == "compress",
                            GzipCompressorStream,
                            GzipDecompressorStream)
    elseif compression == "none"
        compressor = NoopStream  # no compression
    end
    return compressor
end

"""
    root_nodes(graph::DispatchGraph) ->

Return an iterable of all nodes in the graph with no input edges.
"""
function root_nodes(graph::DispatchGraph)
    imap(n->graph.nodes[n], filter(1:nv(graph.graph)) do node_index
        indegree(graph.graph, node_index) == 0
    end)
end


text:
using FredData
using Test

function with_key_env(f::Function, key::AbstractString)
    withenv(FredData.KEY_ENV_NAME => key) do
        f()
    end
end

function with_key_file(f::Function, key::AbstractString)
    mktempdir() do tmpdir
        open(joinpath(tmpdir, FredData.KEY_FILE_NAME), "w") do f
            write(f, key)
        end
        withenv(FredData.KEY_ENV_NAME => nothing,
                "USERPROFILE" => tmpdir,
                "HOME" => tmpdir, ) do
            f()
        end
    end
end

function with_key_none(f::Function)
    mktempdir() do tmpdir
        withenv(FredData.KEY_ENV_NAME => nothing,
                "USERPROFILE" => tmpdir,
                "HOME" => tmpdir, ) do
            f()
        end
    end
end

@testset "Client creation with key" begin
    fake_key = repeat("0", FredData.API_KEY_LENGTH)
    fake_key1 = repeat("1", FredData.API_KEY_LENGTH)

    # pass key directly
    f1 = Fred(fake_key)
    @test f1.key == fake_key

    # detect from ENV
    with_key_env(fake_key) do
        f2 = Fred()
        @test f2.key == fake_key
    end

    # detect from ~/.freddatarc
    # from libuv::uv_os_homedir, we find we need to set USERPROFILE for windows
    # and HOME for *nix
    with_key_file(fake_key1) do
        f3 = Fred()
        @test f3.key == fake_key1
    end
end

@testset "Client creation fails with invalid/missing key" begin
    bad_key_short = repeat("c",
                           convert(Int, round(FredData.API_KEY_LENGTH/2)))
    @test_throws Exception Fred(bad_key_short)

    # no key anywhere
    with_key_none() do
        @test_throws Exception Fred()
    end
end


text:
using FountainCodes
using Random, StatsBase, Distributions
using LinearAlgebra, SparseArrays

function r10_example(nsrc=10, nenc=20)
    code = R10(nsrc)
    println("Example for $code w. $nsrc source symbols and $nenc coded symbols")
    src = rand(GF256, nsrc)
    println("source symbols:")
    println(src)
    println()

    # compute the intermediate symbols
    A = constraint_matrix(code, 0:nsrc-1)
    pre = zeros(GF256, size(A, 2))
    pre[end-nsrc+1:end] .= src
    int = decode(A, pre)
    println("intermediate symbols:")
    println(int)
    println()

    # randomly selected nenc unique encodes symbol identifiers, corresponding to the nenc encoded
    # symbols received by the decoder (the upper limit is chosen arbitrarily)
    Xs = sample(0:1000, nenc, replace=false)

    # compute the encoded symbols
    G = generator_matrix(code, Xs)
    enc = G'*int
    println("encoded symbols:")
    println(enc)
    println()

    # decode the intermediate symbols from the encoded symbols
    A = constraint_matrix(code, Xs)
    pre = zeros(GF256, size(A, 2))
    pre[end-nenc+1:end] .= enc
    dec_int = decode(A, pre)
    println("decoded intermediate symbols:")
    println(dec_int)
    println()

    # generate the source symbols from the decoded intermediate symbols
    G = generator_matrix(code, 0:nsrc-1)
    dec_src = G'*dec_int
    println("decoded source symbols:")
    println(dec_src)
end

function rq_example(nsrc=10, nenc=20)
    code = RQ(nsrc)
    nsrc = code.Kp # w. zero-padding up to the next supported number of source symbols
    println("Example for $code w. $nsrc source symbols and $nenc coded symbols")
    src = rand(GF256, nsrc)
    println("source symbols:")
    println(src)
    println()

    # compute the intermediate symbols
    A = constraint_matrix(code, 0:nsrc-1)
    pre = zeros(GF256, size(A, 2))
    pre[end-nsrc+1:end] .= src
    int = decode(A, pre)
    println("intermediate symbols:")
    println(int)
    println()

    # randomly selected nenc unique encodes symbol identifiers, corresponding to the nenc encoded
    # symbols received by the decoder (the upper limit is chosen arbitrarily)
    Xs = sample(0:1000, nenc, replace=false)

    # compute the encoded symbols
    G = generator_matrix(code, Xs)
    enc = G'*int
    println("encoded symbols:")
    println(enc)
    println()

    # decode the intermediate symbols from the encoded symbols
    A = constraint_matrix(code, Xs)
    pre = zeros(GF256, size(A, 2))
    pre[end-nenc+1:end] .= enc
    dec_int = decode(A, pre)
    println("decoded intermediate symbols:")
    println(dec_int)
    println()

    # generate the source symbols from the decoded intermediate symbols
    G = generator_matrix(code, 0:nsrc-1)
    dec_src = G'*dec_int
    println("decoded source symbols:")
    println(dec_src)
end

function lt_example(nsrc=10, nenc=20)
    dd = Soliton(nsrc, nsrc, 1e-6)
    code = LT{GF256}(nsrc, dd) # w. generator matrix coefficients of type GF256
    println("Example for $code w. $nsrc source symbols and $nenc coded symbols")
    src = rand(GF256, nsrc)
    println("source symbols:")
    println(src)
    println()

    # randomly selected nenc unique encodes symbol identifiers, corresponding to the nenc encoded
    # symbols received by the decoder (the upper limit is chosen arbitrarily)
    Xs = sample(0:1000, nenc, replace=false)

    # compute the encoded symbols
    G = generator_matrix(code, Xs)
    enc = G'*src
    println("encoded symbols:")
    println(enc)
    println()

    # decode the source symbols from the encoded symbols
    dec = decode(G, enc)
    println("decoded intermediate symbols:")
    println(dec)
end

using DelimitedFiles

function ldpc_example(nerasures=20)
    # read the parity check matrix
    # (transposed, since the decoder expects columns to correspond to constraints and rows to symbols)
    # despite it being a binary matrix, we're using symbols of type GF256 since there are no
    # optimizations specific to the binary field
    H = SparseMatrixCSC{GF256,Int}(readdlm("./examples/H_612_1224.txt")')
    # H = SparseMatrixCSC{GF256,Int}(readdlm("./examples/H_2400_4800.txt")')
    nenc, nconstraints = size(H)

    # since we do not consider encoding, we use the all-zeros vector, which is a valid codeword in
    # any linear code
    enc = zeros(GF256, nenc)

    # randomly select nerasures symbols to erase
    Xs = collect(1:nenc)
    shuffle!(Xs)
    Xs_rec = sort!(Xs[1:nenc-nerasures])
    Xs_era = sort!(Xs[end-nerasures+1:end])
    H_rec = H[Xs_rec, :] # constraints corresponding to received symbols
    H_era = H[Xs_era, :] # constraints corresponding to erased symbols
    enc_rec = enc[Xs_rec]
    enc_era = enc[Xs_era]
    println("erased encoded symbols:")
    println(enc_era)
    println()

    # note that H'*enc = H_rec'*enc_rec + H_era'*enc_era = 0, i.e., we can solve for enc_era from
    # the system of equations H_era'*enc_era = - H_rec'*enc_rec = H_rec'*enc_rec,
    # where the second equality is due to the fact that -GF256(x) = GF256(x)
    # finally, assuming that the code is systematic, the source symbols can be recovered from the
    # decoded enc vector by selecting the subset of its indices corresponding to the systematic entries
    rhs = H_rec'*enc_rec
    enc_era_dec = decode(H_era, rhs)
    println("decoded erased symbols:")
    println(enc_era_dec)
end

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
# handles google authentication

# temp loads while testing
using JSON3
using HTTP
using JSON
using Base64
import Base: show
import MbedTLS
using Dates
using Sockets
using Lazy

const api_scope = "https://www.googleapis.com/auth/analytics.readonly"
const api_endpoint = "https://analyticsreporting.googleapis.com/v4/reports:batchGet"

struct Servicekey
    type::String
    project_id::String
    private_key_id::String
    private_key::String
    client_email::String
    client_id::String
    auth_uri::String
    token_uri::String
    auth_provider_x509_cert_url::String
    client_x509_cert_url::String
end

StructTypes.StructType(::Type{Servicekey}) = StructTypes.Struct()

mutable struct Credentials
    access_token::String
    refresh_token::String
end

# functions prototypes
function load_keyfile(file_location::AbstractString) end
function Base.show(io::IO, key::Servicekey) end
function refresh_access_token(creds::Credentials) end

# service key functions
function load_keyfile(file_location::AbstractString)
    open(file_location, "r") do f
        JSON3.read(f, Servicekey)
    end
end

# extend base show function for Servicekey objects
function Base.show(io::IO, key::Servicekey)
    print(io, "Service key data:\n")
    print(io, "type:\t", key.type, "\n")
    print(io, "project_id:\t", key.project_id, "\n")
    print(io, "private_key_id:\t", key.private_key_id, "\n")
    print(io, "private_key:\t", key.private_key, "\n")
    print(io, "client_email:\t", key.client_email, "\n")
    print(io, "client_id:\t", key.client_id, "\n")
    print(io, "auth_uri:\t", key.auth_uri, "\n")
    print(io, "token_uri:\t", key.token_uri, "\n")
    print(io, "auth_provider_x509_cert_url:\t", key.auth_provider_x509_cert_url, "\n")
    print(io, "client_x509_cert_url:\t", key.client_x509_cert_url, "\n")
end

# ------------------ Shamelessly lifted from GitHub.jl (thanks for saving my sanity!)
const ENTROPY = Ref{MbedTLS.Entropy}()
const RNG     = Ref{MbedTLS.CtrDrbg}()

ENTROPY[] = MbedTLS.Entropy()
RNG[]     = MbedTLS.CtrDrbg()
MbedTLS.seed!(RNG[], ENTROPY[])

#######################
# Authorization Types #
#######################

abstract type Authorization end

# TODO: SecureString on 0.7
struct OAuth2 <: Authorization
    token::String
end

struct UsernamePassAuth <: Authorization
    username::String
    password::String
end

struct AnonymousAuth <: Authorization end

struct JWTAuth <: Authorization
    JWT::String
end

####################
# JWT Construction #
####################

function base64_to_base64url(string)
    replace(replace(replace(string, "=" => ""), '+' => '-'), '/' => '_')
end

function JWTAuth(app_id::Int, key::MbedTLS.PKContext; servicekey::Servicekey, iat = now(Dates.UTC), exp_mins = 1)
    algo = base64_to_base64url(base64encode(JSON.json(Dict(
        "alg" => "RS256",
        "typ" => "JWT"
    ))))
    data = base64_to_base64url(base64encode(JSON.json(Dict(
        "iss"      => servicekey.client_email,
        "scope"    => "https://www.googleapis.com/auth/analytics.readonly",
        "aud"      => servicekey.token_uri,
        "exp" => trunc(Int64, Dates.datetime2unix(iat+Dates.Minute(exp_mins))),
        "iat" => trunc(Int64, Dates.datetime2unix(iat))

    ))))
    signature = base64_to_base64url(base64encode(MbedTLS.sign(key, MbedTLS.MD_SHA256,
        MbedTLS.digest(MbedTLS.MD_SHA256, string(algo,'.',data)), RNG[])))
    JWTAuth(string(algo,'.',data,'.',signature))
end

function JWTAuth(app_id::Int, privkey::String; kwargs...)
    JWTAuth(app_id, MbedTLS.parse_keyfile(privkey); kwargs...)
end

#########################
# Header Authentication #
#########################

authenticate_headers!(headers, auth::AnonymousAuth) = headers

function authenticate_headers!(headers, auth::OAuth2)
    headers["Authorization"] = "token $(auth.token)"
    return headers
end

function authenticate_headers!(headers, auth::JWTAuth)
    headers["Authorization"] = "Bearer $(auth.JWT)"
    return headers
end

function authenticate_headers!(headers, auth::UsernamePassAuth)
    headers["Authorization"] = "Basic $(base64encode(string(auth.username, ':', auth.password)))"
    return headers
end

###################
# Pretty Printing #
###################

function Base.show(io::IO, a::OAuth2)
    token_str = a.token[1:6] * repeat("*", length(a.token) - 6)
    print(io, "GitHub.OAuth2($token_str)")
end

token_endpoint = "https://oauth2.googleapis.com/token"
grant = HTTP.escapeuri("urn:ietf:params:oauth:grant-type:jwt-bearer")
auth = JWTAuth(0, "examples/test.pem"; servicekey = servicekey).JWT
# payload = "grant_type=urn%3Aietf%3Aparams%3Aoauth%3Agrant-type%3Ajwt-bearer&assertion=$auth"
payload = Dict(
    "grant_type" => "urn:ietf:params:oauth:grant-type:jwt-bearer",
    "assertion" => auth
)
raw_response = HTTP.request("POST", token_endpoint, ["Content-Type" => "application/json"], JSON.json(payload))
response = String(raw_response.body)
access_token_blob = JSON.parse(response)

# test request
access_token = access_token_blob["access_token"]
token_type = access_token_blob["token_type"]
payload = Dict(
        "reportRequests" => [Dict(
            "viewId"        => "43047246",
            "dateRanges"    => [Dict(
                                    "startDate" => "7daysAgo",
                                    "endDate" => "today"),],
          "metrics"         => [Dict("expression" => "ga:sessions"),],
          "dimensions" => [Dict("name" => "ga:country"),]
        ),]
      )

batchGet_endpoint = "https://analyticsreporting.googleapis.com/v4/reports:batchGet"
raw_response = HTTP.request("POST",
    batchGet_endpoint,
    ["Authorization" => "$token_type $access_token"],
    JSON.json(payload)
)
response = String(raw_response.body)
JSON.parse(response)


text:
"""
Distributionally robust extension
"""

abstract type AmbiguitySet end
abstract type WassersteinSet <: AmbiguitySet end

mutable struct DroData{WassersteinSet}
    lp_norm::Float64
    ϵ::Float64
end

"""
    set

Set the parameters for Wasserstein ambiguity set

# Arguments
- `dsp`: DSPProblem struct
- `p`: norm of the Wasserstein distance
- `ϵ`: the size of the Wasserstein ball
"""
function set(::Type{WassersteinSet}, p, ϵ)
    dspenv.dro = DroData{WassersteinSet}(p, ϵ)
end

function loadDroData(dsp::DSPProblem, dro::Nothing)
end

function loadDroData(dsp::DSPProblem, dro::DroData{WassersteinSet})
    @dsp_ccall("setWassersteinAmbiguitySet", Cvoid,
        (Ptr{Cvoid}, Cdouble, Cdouble),
        dsp.p, dro.lp_norm, dro.ϵ)
end


text:
using SafeTestsets

using Test, Documenter, ThreeBodyProblem
# doctest(ThreeBodyProblem)

@safetestset "dynamics.jl tests" begin
    include("dynamics_tests.jl")
end

@safetestset "frames.jl tests" begin
    include("frames_tests.jl")
end

@safetestset "orbitalelements.jl tests" begin
    include("orbitalelements_tests.jl")
end

@safetestset "parameters.jl tests" begin
    include("parameters_tests.jl")
end

@safetestset "plot.jl tests" begin
    include("plot_tests.jl")
end

@safetestset "special.jl tests" begin
    include("special_tests.jl")
end

@safetestset "util.jl tests" begin
    include("util_tests.jl")
end


text:
import Nemo.sub!, Base.gcd
export induce_rational_reconstruction, induce_crt, root, roots,
       ismonic, radical_extension, isradical_extension,
       iskummer_extension, cyclotomic_field, wildanger_field,
       compositum

add_verbose_scope(:PolyFactor)
add_assert_scope(:PolyFactor)
add_verbose_scope(:CompactPresentation)
add_assert_scope(:CompactPresentation)

if Int==Int32
  global const p_start = next_prime(2^30)
else
  global const p_start = next_prime(2^60)
end

include("NfAbs/NfAbs.jl")
include("NfAbs/Conjugates.jl")
include("NfAbs/CompactRepresentation.jl")
include("NfAbs/Elem.jl")
include("NfAbs/NonSimple.jl")
include("NfAbs/Poly.jl")
include("NfAbs/Simplify.jl")
include("NfAbs/DiscriminantBounds.jl")
include("NfAbs/Hilbert.jl")
include("NfAbs/NormRelation.jl")
include("NfAbs/PolyFact.jl")
include("NfAbs/MPolyGcd.jl")
include("NfAbs/ConjugatesNS.jl")


text:
include("utils.jl")

build_so("003")

x = "abcdef"

y = ccall(
    (:foo, "binary/lib003.so"),
    Void,
    (Ptr{Uint8}, Cint),
    x, length(x)
)

@assert x == "bcdefg"

@assert y === nothing


text:
# #-------------------------------------------------------------------------------------------------------
# #-------------------------------------------------------------------------------------------------------
# # This file contains functions that are used for the model calibrations.
# #-------------------------------------------------------------------------------------------------------
# #-------------------------------------------------------------------------------------------------------

#######################################################################################################################
# LOAD AND CLEAN UP DATA USED FOR MODEL CALIBRATION.
#######################################################################################################################
# Description: This function loads, cleans up, and merges all of the calibration data into a single dataframe.
#
# Function Arguments:
#
#       model_end = The final year to include in the calibration data set (defaults to start in 1765).
#----------------------------------------------------------------------------------------------------------------------

function load_calibration_data(model_end::Int)

    # Set model start year to 1765 by default.
    model_start = 1765

    # Create column of calibration years and calculate indicies for calibration time period (relative to 1765-2018).
    df = DataFrame(year = collect(1765:2018))
    model_calibration_indices = findall((in)(collect(model_start:model_end)), collect(1765:2018))

    #-------------------------------------------------------------------
    # HadCRUT4 temperature data (anomalies relative to 1861-1880 mean).
    #-------------------------------------------------------------------

    # Load raw temperature data.
    raw_temp_data = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "calibration_data", "global_temp_hadcrut4.csv"), skiplines_begin=24))

    # Find indices to normalize temperature data to 1861-1880 mean.
    hadcrut_norm_indices = findall((in)(collect(1861:1880)), collect(1850:2017))

    # Normalize temperature data to 1861-1880 mean.
    norm_temp_data  = DataFrame(year=raw_temp_data[!,:year], hadcrut_temperature_obs = raw_temp_data[!,:median] .- mean(raw_temp_data[hadcrut_norm_indices, :median]))

    # Join data on year.
    df = join(df, norm_temp_data, on=:year, kind=:outer)

    # Read in HadCRUT4 1σ errors and rename column.
    raw_temp_errors  = DataFrame(load(joinpath(@__DIR__,  "..", "..", "data", "calibration_data", "global_temp_hadcrut4_1sigma_uncertainty.csv"), skiplines_begin=21))
    rename!(raw_temp_errors, :one_sigma_all => :hadcrut_temperature_sigma)

    # Join data on year.
    df = join(df, raw_temp_errors[!, [:year, :hadcrut_temperature_sigma]], on=:year, kind=:outer)

    #---------------------------------------------------------------------------------
    # Annual Global Ocean Heat Content (0-3000 m).
    #---------------------------------------------------------------------------------

    # Load ocean heat content (0-3000m) observations and errors.
    ocean_heat_raw = DataFrame(load(joinpath(@__DIR__,  "..", "..", "data", "calibration_data", "ocean_heat_gouretski_3000m.csv"), colnames=["year", "ocean_heat_obs", "ocean_heat_sigma"], skiplines_begin=3))

    # Join data on year.
    df = join(df, ocean_heat_raw, on=:year, kind=:outer)

    #--------------------------------------------------------
    # Mauna Loa Instrumental Atmospheric CO₂ Concentrations.
    #--------------------------------------------------------

    # Load Mauna Loa CO₂ observations and errors, and rename columns.
    raw_mauna_loa_co2_data  = DataFrame(load(joinpath(@__DIR__,  "..", "..", "data", "calibration_data", "co2_mauna_loa.csv"), skiplines_begin=59))
    rename!(raw_mauna_loa_co2_data, :mean => :maunaloa_co2_obs, :unc => :maunaloa_co2_sigma)

    # Join data on year.
    df = join(df, raw_mauna_loa_co2_data, on=:year, kind=:outer)

    #-----------------------------------------------------
    # Law Dome Ice Core Atmospheric CO₂ Concentrations.
    #-----------------------------------------------------

    # Load Law Dome CO₂ observations and errors, and rename columns.
    raw_law_dome_co2_data = DataFrame(load(joinpath(@__DIR__,  "..", "..", "data", "calibration_data", "law_dome_co2.csv"), skiplines_begin=4))
    rename!(raw_law_dome_co2_data, :co2_ice => :lawdome_co2_obs, :one_sigma_error => :lawdome_co2_sigma)

    # Join data on year.
    df = join(df, raw_law_dome_co2_data, on=:year, kind=:outer)

    #---------------------------------------------------------------------------------
    # Decadal Ocean Carbon Fluxes
    #---------------------------------------------------------------------------------

    # Observations and errors from McNeil et al. (2003).
    ocean_co2_flux_data = DataFrame(year=[1985, 1995], oceanco2_flux_obs=[1.6, 2.0], oceanco2_flux_sigma=[0.4, 0.4])

    # Join data on year.
    df = join(df, ocean_co2_flux_data, on=:year, kind=:outer)

    #---------------------------------------------------------------------------------
    # NOAA Instrumental Atmospheric CH₄ Concentrations.
    #---------------------------------------------------------------------------------

    # Load NOAA CH₄ observations and errors, and rename columns.
    raw_noaa_ch4_data = DataFrame(load(joinpath(@__DIR__,  "..", "..", "data", "calibration_data", "ch4_noaa.csv"), skiplines_begin=58))
    rename!(raw_noaa_ch4_data, :mean => :noaa_ch4_obs, :unc => :noaa_ch4_sigma)

    # Join data on year.
    df = join(df, raw_noaa_ch4_data, on=:year, kind=:outer)

    #---------------------------------------------------------------------------------
    # Law Dome Ice Core Atmospheric CH₄ Concentrations.
    #---------------------------------------------------------------------------------

    # Load Law Dome CH₄ observations and dates.
    raw_law_dome_ch4_data = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "calibration_data", "ch4_ice_core_etheridge_1998.csv"), skiplines_begin=14))

    # Get all listed years and number of unique years (multiple observations may occur for a single year).
    years = trunc.(Int, raw_law_dome_ch4_data.air_age)
    unique_years = unique(raw_law_dome_ch4_data.air_age)

    # Create array to hold annual CH₄ values.
    n_years = length(unique_years)
    law_dome_ch4_annual_avg = DataFrame(year=zeros(Int64, n_years), ch4=zeros(n_years))

    # For years with multiple observations, calculate mean CH₄ value.
    for t in 1:n_years
        #Find indices where there are multiple observations across a single year.
        index = findall(unique_years[t] .== years)
        #Assign year and mean value for array with one value per year.
        law_dome_ch4_annual_avg[t,:year] = unique_years[t]
        law_dome_ch4_annual_avg[t,:ch4]  = mean(raw_law_dome_ch4_data[index, :CH4])
    end

    # Load CH₄ ice firn observations and dates (also Etheridge 1998).
    raw_ch4_firn_data = DataFrame(load(joinpath(@__DIR__,  "..", "..", "data", "calibration_data", "ch4_firn_etheridge_1998.csv"), skiplines_begin=14))

    # Get years (averaging data together for multiple observations in the same year).
    years = trunc.(Int, raw_ch4_firn_data[:,:air_age])
    unique_years = unique(years)

    #Create array to hold annual CH₄ ice firn values.
    n_years = length(unique_years)
    firn_ch4_annual_avg = DataFrame(year=zeros(Int64, n_years), ch4=zeros(n_years))

    #Calculate mean value for each unique year.
    for t in 1:n_years
        #Find indices where there are multiple observations across a single year.
        index = findall(unique_years[t] .== years)
        #Assign year and mean value for array with one value per year.
        firn_ch4_annual_avg[t,:year] = unique_years[t]
        firn_ch4_annual_avg[t,:ch4] = mean(raw_ch4_firn_data[index, :CH4])
    end

    # Merge data for ice core (1852-1977) and firn (1978-1981) that occur during calibration period. After that, calibration uses NOAA flask data.
    core_indices = findall((in)(collect(1852:1977)), law_dome_ch4_annual_avg[:,:year])
    firn_indices = findall((in)(collect(1978:1981)), firn_ch4_annual_avg[:,:year])
    ch4_ice = vcat(law_dome_ch4_annual_avg[core_indices,:], firn_ch4_annual_avg[firn_indices,:])

    # Do a linear interpolation between interpolar (N-S) differences as described in Etheridge et al. (1998).
    # Uses reported Greenland-Antarctica difference of ≈ 41.9 ppb for 1842 and 143 ppb for 1980s (assumed centered on 1985).
    start_year = 1842
    end_year   = 1985
    n_years    = length(start_year:end_year)

    # Initialize array to hold merged, re-scaled global CH₄ data.
    interp_vals = zeros(n_years, 2)
    interp_vals[:,1] = collect(start_year:end_year)

    # Set initial value.
    interp_vals[1,2] = 41.9

    #Calculate annual amount to add for each year in interpolation.
    adder = (143-41.9) / (length(start_year:end_year)-1)

    # Calculate interpolar values to be added to ice core data over time if we had annually consecutive observations.
    for t in 2:n_years
        interp_vals[t,2] = interp_vals[t-1,2] + adder
    end

    #Find indices of annual interpolated interpolar differences that correspond to ice core CH₄ observations.
    interp_index = findall((in)(ch4_ice[!,:year]), interp_vals[:,1])

    # Approximate global CH4 concentration as Antartica ice core values plus 37% of interpolar difference following Etheridge et al. (1998).
    ch4_ice[!,:global_ch4] = ch4_ice[!,:ch4] .+ interp_vals[interp_index, 2] .* 0.37

    # Set observation error for CH₄ ice core data as 15 ppb.
    # Etheridge notes 5ppb (1σ) measurement error and global calculation (using interhemispheric mixing assumptions) could introduce an error of ≈ 10ppb.
    ch4ice_error = ones(size(ch4_ice)[1]) .* 15.0
    final_law_dome_ch4_data = DataFrame(year=ch4_ice[!,:year],  lawdome_ch4_obs = ch4_ice[!,:global_ch4], lawdome_ch4_sigma = ch4ice_error)

    # Join re-scaled global CH₄ observations with rest of calibration data.
    df = join(df, final_law_dome_ch4_data, on=:year, kind=:outer)

    #---------------------------------------------------------------------------------
    # Finalize Joint Calibration Data Set.
    #---------------------------------------------------------------------------------

    # Sort all calibration data by year.
    sort!(df, :year)

    # Crop data to appropriate calibration years and return.
    return df[model_calibration_indices, :]
end

#######################################################################################################################
# BRIDGE SAMPLING FOR MODEL MARGINAL LIKELIHOODS.
#######################################################################################################################
# Description: This function carries out a single bridge sampling iteration for calculating the marginal likelihoods
#              of the four versions of SNEASY+CH4. It closely follows the code from "Neglecting model structural
#              uncertainty underestimates upper tails of flood hazard (Wong et al., 2018).
#
# Function Arguments:
#
#       norm_constant               = Log normalizing constant estimate (passed in as the bridge sampling estimate from previous iteration).
#       mcmc_posterior              = Vector of log-posterior estimates evaluated for each calibrated posterior parameter sample (from MCMC calibration).
#       mcmc_importance_density     = Density of posterior parameter samples (from MCMC calibration) evaluated on the importance distribution.
#       importance_sample_posterior = Vector of log-posterior estimates evaluated for each sample from the importance distribution.
#       importance_sample_density   = Density of importance distribution samples evaluated on the importance distribution.
#----------------------------------------------------------------------------------------------------------------------

function bridge_sampling_iteration(norm_constant, mcmc_posterior, mcmc_importance_density, importance_sample_posterior, importance_sample_density)

    # Filter out non-finite values and normalize finite values based on normalizing constant (i.e. previous bridge sampling estimate).
    mcmc_posterior_norm = mcmc_posterior[isfinite.(mcmc_posterior)] .- norm_constant
    importance_sample_posterior_norm = importance_sample_posterior[isfinite.(importance_sample_posterior)] .- norm_constant

    mcmc_importance_density = mcmc_importance_density[isfinite.(mcmc_posterior)]
    importance_sample_density = importance_sample_density[isfinite.(importance_sample_posterior)]

    # Get number of samples with finite values.
    n_mcmc = length(mcmc_posterior_norm)
    n_importance = length(importance_sample_posterior_norm)

    # Compute updated estimates for numerator and denominator means.
    # Note: A small subset of cases produce non-finite values that lead to an NaN result. First calculate temporary variables with these potential NaNs and then filter them out.
    importance_temp = exp.(importance_sample_posterior_norm) ./ (n_importance .* exp.(importance_sample_density) .+ n_mcmc .* exp.(importance_sample_posterior_norm))
    mcmc_temp       = exp.(mcmc_importance_density) ./ (n_importance .* exp.(mcmc_importance_density) .+ n_mcmc .* exp.(mcmc_posterior_norm))

    # Calcualte updated mean values for non-NaN results.
    importance_mean =  mean(importance_temp[findall(x -> !isnan(x), importance_temp)])
    mcmc_mean       =  mean(mcmc_temp[findall(x -> !isnan(x), mcmc_temp)])

    # Return updated bridge sampling estimate.
    return (norm_constant + log(importance_mean) - log(mcmc_mean))
end

#######################################################################################################################
# CALCULATE BAYESIAN MODEL AVERAGING (BMA) WEIGHTS USING BRIDGE SAMPLING.
#######################################################################################################################
# Description: This function uses bridge sampling to calculate the marginal likelihoods of the four versions of
#              SNEASY+CH4. It uses a multivariate normal approximation to the joint posterior distribtion as the
#              importance distribution and closely follows the code from "Neglecting model structural uncertainty
#              underestimates upper tails of flood hazard (Wong et al., 2018).
#
# Function Arguments:
#
#       parameters_fairch4      = Calibrated parameters for SNEASY-FAIR (each row is a new sample from joint posterior distribution, each column is a different parameter).
#       parameters_fundch4      = Calibrated parameters for SNEASY-FUND (each row is a new sample from joint posterior distribution, each column is a different parameter).
#       parameters_hectorch4    = Calibrated parameters for SNEASY-Hector (each row is a new sample from joint posterior distribution, each column is a different parameter).
#       parameters_magiccch4    = Calibrated parameters for SNEASY-MAGICC (each row is a new sample from joint posterior distribution, each column is a different parameter).
#       log_posterior_fairch4   = Function to calculate log-posterior for SNEASY-FAIR.
#       log_posterior_fundch4   = Function to calculate log-posterior for SNEASY-FUND.
#       log_posterior_hectorch4 = Function to calculate log-posterior for SNEASY-Hector.
#       log_posterior_magiccch4 = Function to calculate log-posterior for SNEASY-MAGICC.
#----------------------------------------------------------------------------------------------------------------------

function calculate_bma_weights(parameters_fairch4::Array{Float64,2}, parameters_fundch4::Array{Float64,2}, parameters_hectorch4::Array{Float64,2}, parameters_magiccch4::Array{Float64,2}, log_posterior_fairch4, log_posterior_fundch4, log_posterior_hectorch4, log_posterior_magiccch4)

    # Calculate number of posterior parameter samples.
    n_samples = size(parameters_fairch4, 1)

    # Pre-allocate arrays to store log-posterior values resulting from sampled posterior and importance pdf model parameters.
    posterior_mcmc_fair   = zeros(n_samples)
    posterior_mcmc_fund   = zeros(n_samples)
    posterior_mcmc_hector = zeros(n_samples)
    posterior_mcmc_magicc = zeros(n_samples)

    posterior_importance_fair   = zeros(n_samples)
    posterior_importance_fund   = zeros(n_samples)
    posterior_importance_hector = zeros(n_samples)
    posterior_importance_magicc = zeros(n_samples)

    # Calculate log-posterior values resulting from calibrated posterior model parameters.
    for i = 1:n_samples
        posterior_mcmc_fair[i]   = log_posterior_fairch4(parameters_fairch4[i,:])
        posterior_mcmc_fund[i]   = log_posterior_fundch4(parameters_fundch4[i,:])
        posterior_mcmc_hector[i] = log_posterior_hectorch4(parameters_hectorch4[i,:])
        posterior_mcmc_magicc[i] = log_posterior_magiccch4(parameters_magiccch4[i,:])
    end

    # Calculate posterior parameter mean and covariance matrix, then fit to a multivariate normal distribution (importance pdf).
    posterior_mean_fair = vec(mean(parameters_fairch4, dims=1))
    posterior_cov_fair  = cov(parameters_fairch4)
    fair_importance_pdf = MvNormal(posterior_mean_fair, posterior_cov_fair)

    posterior_mean_fund = vec(mean(parameters_fundch4, dims=1))
    posterior_cov_fund  = cov(parameters_fundch4)
    fund_importance_pdf = MvNormal(posterior_mean_fund, posterior_cov_fund)

    posterior_mean_hector = vec(mean(parameters_hectorch4, dims=1))
    posterior_cov_hector  = cov(parameters_hectorch4)
    hector_importance_pdf = MvNormal(posterior_mean_hector, posterior_cov_hector)

    posterior_mean_magicc = vec(mean(parameters_magiccch4, dims=1))
    posterior_cov_magicc  = cov(parameters_magiccch4)
    magicc_importance_pdf = MvNormal(posterior_mean_magicc, posterior_cov_magicc)

    # Calculate importance densities using calibrated posterior parameters.
    importance_mcmc_density_fair   = logpdf(fair_importance_pdf,   transpose(parameters_fairch4))
    importance_mcmc_density_fund   = logpdf(fund_importance_pdf,   transpose(parameters_fundch4))
    importance_mcmc_density_hector = logpdf(hector_importance_pdf, transpose(parameters_hectorch4))
    importance_mcmc_density_magicc = logpdf(magicc_importance_pdf, transpose(parameters_magiccch4))

    # Sample parameters from fitted importance pdf.
    importance_sample_fair   = rand(fair_importance_pdf,   n_samples)
    importance_sample_fund   = rand(fund_importance_pdf,   n_samples)
    importance_sample_hector = rand(hector_importance_pdf, n_samples)
    importance_sample_magicc = rand(magicc_importance_pdf, n_samples)

    # Calculate importance densities using samples from importance pdf.
    importance_sample_density_fair   = logpdf(fair_importance_pdf,   importance_sample_fair)
    importance_sample_density_fund   = logpdf(fund_importance_pdf,   importance_sample_fund)
    importance_sample_density_hector = logpdf(hector_importance_pdf, importance_sample_hector)
    importance_sample_density_magicc = logpdf(magicc_importance_pdf, importance_sample_magicc)

    # Calculate log-posterior values of samples from importance pdf.
    # Note: A small subset of samples from importance pdf produce non-physical model results, so also set them as -Inf to be filtered out.

    # SNEASY-FAIR
    for i = 1:length(posterior_importance_fair)
        try
            posterior_importance_fair[i] = log_posterior_fairch4(importance_sample_fair[:,i])
        catch
            posterior_importance_fair[i] = -Inf
        end
    end

    # SNEASY-FUND
    for i = 1:length(posterior_importance_fund)
        try
            posterior_importance_fund[i] = log_posterior_fundch4(importance_sample_fund[:,i])
        catch
            posterior_importance_fund[i] = -Inf
        end
    end

    # SNEASY-Hector
    for i = 1:length(posterior_importance_hector)
        try
            posterior_importance_hector[i] = log_posterior_hectorch4(importance_sample_hector[:,i])
        catch
            posterior_importance_hector[i] = -Inf
        end
    end

    # SNEASY-MAGICC
    for i = 1:length(posterior_importance_magicc)
        try
            posterior_importance_magicc[i] = log_posterior_magiccch4(importance_sample_magicc[:,i])
        catch
            posterior_importance_magicc[i] = -Inf
        end
    end

    # Set convergence tolerance and maximum number of bridge-sampling iterations for all models.
    max_iterations = 10_000
    tolerance = 1e-10

    # Set initial conditions to begin iterating over for all models (following Wong et al., 2018, set initial value by "averaging the ratios on a log scale.").
    x_fair   = [0.0, -mean(importance_mcmc_density_fair   - posterior_mcmc_fair)]
    x_fund   = [0.0, -mean(importance_mcmc_density_fund   - posterior_mcmc_fund)]
    x_hector = [0.0, -mean(importance_mcmc_density_hector - posterior_mcmc_hector)]
    x_magicc = [0.0, -mean(importance_mcmc_density_magicc - posterior_mcmc_magicc)]

    #-------------------------------------------------------------
    # Use bridge sampling to estimate model marginal likelihoods.
    #-------------------------------------------------------------

    # Iterate over SNEASY-FAIR values until tolerance or maximum iterations reached.
    let iterations = 0
        while abs(x_fair[2] - x_fair[1]) >= tolerance && iterations <= max_iterations
            # Set x[1] to previous iteration's estimated value.
            x_fair[1] = x_fair[2]
            # Calculate new x[2] value for current iteration.
            x_fair[2] = bridge_sampling_iteration(x_fair[1], posterior_mcmc_fair, importance_mcmc_density_fair, posterior_importance_fair, importance_sample_density_fair)
            # Add one to iteration.
            iterations += 1
        end
    end

    # Iterate over SNEASY-FUND values until tolerance or maximum iterations reached.
    let iterations = 0
        while abs(x_fund[2] - x_fund[1]) >= tolerance && iterations <= max_iterations
            # Set x[1] to previous iteration's estimated value.
            x_fund[1] = x_fund[2]
            # Calculate new x[2] value for current iteration.
            x_fund[2] = bridge_sampling_iteration(x_fund[1], posterior_mcmc_fund, importance_mcmc_density_fund, posterior_importance_fund, importance_sample_density_fund)
            # Add one to iteration.
            iterations += 1
        end
    end

    # Iterate over SNEASY-Hector values until tolerance or maximum iterations reached.
    let iterations = 0
        while abs(x_hector[2] - x_hector[1]) >= tolerance && iterations <= max_iterations
            # Set x[1] to previous iteration's estimated value.
            x_hector[1] = x_hector[2]
            # Calculate new x[2] value for current iteration.
            x_hector[2] = bridge_sampling_iteration(x_hector[1], posterior_mcmc_hector, importance_mcmc_density_hector, posterior_importance_hector, importance_sample_density_hector)
            # Add one to iteration.
            iterations += 1
        end
    end

    # Iterate over SNEASY-MAGICC values until tolerance or maximum iterations reached.
    let iterations = 0
        while abs(x_magicc[2] - x_magicc[1]) >= tolerance && iterations <= max_iterations
            # Set x[1] to previous iteration's estimated value.
            x_magicc[1] = x_magicc[2]
            # Calculate new x[2] value for current iteration.
            x_magicc[2] = bridge_sampling_iteration(x_magicc[1], posterior_mcmc_magicc, importance_mcmc_density_magicc, posterior_importance_magicc, importance_sample_density_magicc)
            # Add one to iteration.
            iterations += 1
        end
    end

    # Collect marginal likelihoods of each model estimated from bridge sampling approach.
    model_marginal_likelihood = [x_fair[2], x_fund[2], x_hector[2], x_magicc[2]]

    # Normalize model likelihoods relative to largest value.
    normalized_model_marginal_likelihood = model_marginal_likelihood .- maximum(model_marginal_likelihood)

    # Use model marginal likelihoods from bridge sampling to calculate BMA weights.
    bma_weights = exp.(normalized_model_marginal_likelihood) ./ sum(exp.(normalized_model_marginal_likelihood))

    return bma_weights
end

#######################################################################################################################
# CALCULATE AR(1) LOG-LIKELIHOOD.
########################################################################################################################
# Description: This function calculates the AR(1) log-likelihood in terms of the data-model residuls when accounting for
#              time-varying observation errors. It follows "The Effects of Time-Varying Observation Errors on Semi-Empirical
#              Sea-Level Projections" (Ruckert et al., 2017) DOI 10.1007/s10584-016-1858-z.
#
# Function Arguments:
#
#       residuals = A vector of data-model residuals.
#       σ         = AR(1) innovation standard deviation.
#       ρ         = AR(1) autocorrelation term.
#       ϵ         = A vector of time-varying observation error estimates (from calibration data sets).
#----------------------------------------------------------------------------------------------------------------------

function hetero_logl_ar1(residuals::Array{Float64,1}, σ::Float64, ρ::Float64, ϵ::Array{Union{Float64, Missings.Missing},1})

    # Calculate length of residuals.
    n=length(residuals)

    # 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 the log-likelihood.
    return logpdf(MvNormal(cov_matrix), residuals)
end

#######################################################################################################################
# CALCULATE CAR(1) LOG-LIKELIHOOD.
########################################################################################################################
# Description: This function calculates the continuous time autoregressive, or CAR(1), log-likelihood for irregularly
#              spaced data in terms of the data-model residuls when accounting for time-varying observation errors. It
#              builds off of "The Effects of Time-Varying Observation Errors on Semi-Empirical Sea-Level Projections"
#              (Ruckert et al., 2017) DOI 10.1007/s10584-016-1858-z and "The Analysis of Irregularly Observed Stochastic
#              Astronomical Time-Series – I. Basics of Linear Stochastic Differential Equations" (Koen, 2005)
#              doi.org/10.1111/j.1365-2966.2005.09213.x
#
# Function Arguments:
#
#       residuals      = A vector of data-model residuals.
#       indices        = Index positions of observations relative to model time horizon (i.e. the first model time period = 1, the second = 2, etc.).
#       σ²_white_noise = Variance of the continuous white noise process.
#       α₀             = Parameter describing correlation memory of CAR(1) process.
#       ϵ              = A vector of time-varying observation error estimates (from calibration data sets).
#----------------------------------------------------------------------------------------------------------------------

function hetero_logl_car1(residuals::Array{Float64,1}, indices::Array{Int,1}, σ²_white_noise::Float64, α₀::Float64, ϵ::Array{Union{Float64, Missings.Missing},1})

    # Calculate length of residuals.
    n=length(residuals)

    # 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 the log-likelihood.
    return logpdf(MvNormal(cov_matrix), residuals)
end