text:
"""
    L0_reg(x::BEDFile, y, k, kernfile::String)

If supplied a `BEDFile` `x` and an OpenCL kernel file `kernfile` as an String, then `L0_reg` will attempt to accelerate the calculation of the dense gradient `x' * (y - x*b)` with a GPU device.
This variant introduces a host of extra arguments for the GPU encapsulated in an optional `PlinkGPUVariables` argument `v`.
The optional argument `v` facilitates the calculation of a regularization path by `iht_path`.
"""
function L0_reg{T <: Float, V <: DenseVector}(
    x        :: BEDFile{T},
    y        :: V,
    k        :: Int,
    kernfile :: String;
    pids     :: Vector{Int}       = procs(x),
    w     :: IHTVariables{T,V} = IHTVariables(x, y, k),
    mask_n   :: Vector{Int}       = ones(Int, size(y)),
    v        :: PLINK.PlinkGPUVariables{T} = PLINK.PlinkGPUVariables(w.df, x, y, kernfile, mask_n),
    tol      :: T    = convert(T, 1e-4),
    max_iter :: Int  = 100,
    max_step :: Int  = 50,
    quiet    :: Bool = true
)

    # start timer
    tic()

    # first handle errors
    k        >= 0      || throw(ArgumentError("Value of k must be nonnegative!\n"))
    max_iter >= 0      || throw(ArgumentError("Value of max_iter must be nonnegative!\n"))
    max_step >= 0      || throw(ArgumentError("Value of max_step must be nonnegative!\n"))
    tol      >  eps(T) || throw(ArgumentError("Value of global tol must exceed machine precision!\n"))
    n = length(y)
    sum((mask_n .== 1) $ (mask_n .== 0)) == n || throw(ArgumentError("Argument mask_n can only contain 1s and 0s"))

    # initialize return values
    mm_iter   = 0       # number of iterations of L0_reg
    mm_time   = zero(T) # compute time *within* L0_reg
    next_loss = zero(T) # loss function value

    # initialize floats
    current_loss = convert(T, Inf)    # tracks previous objective function value
    the_norm     = zero(T)            # norm(b - b0)
    scaled_norm  = zero(T)            # the_norm / (norm(b0) + 1)
    mu           = zero(T)            # Landweber step size, 0 < tau < 2/rho_max^2

    # initialize integers
    i       = 0         # used for iterations in loops
    mu_step = 0         # counts number of backtracking steps for mu

    # initialize booleans
    converged = false   # scaled_norm < tol?

    # update Xb, r, and gradient
    if sum(w.idx) == 0
        fill!(w.xb, zero(T))
        copy!(sdata(w.r), sdata(y))
        mask!(w.r, mask_n, 0, zero(T))
    else
        A_mul_B!(w.xb, x, w.b, w.idx, k, mask_n)
        difference!(w.r, y, w.xb)
        mask!(w.r, mask_n, 0, zero(T))
    end

    # calculate the gradient using the GPU
    At_mul_B!(w.df, x, w.r, mask_n, v)

    # update loss
    next_loss = convert(T, Inf)

    # formatted output to monitor algorithm progress
    !quiet && print_header()

    # main loop
    for mm_iter = 1:max_iter

        # notify and break if maximum iterations are reached.
        if mm_iter >= max_iter

            # alert about hitting maximum iterations
            !quiet && print_maxiter(max_iter, loss)

            # send elements below tol to zero
            threshold!(b, tol)

            # stop timer
            mm_time = toq()

            # these are output variables for function
            # wrap them into a Dict and return
            return IHTResults(mm_time, next_loss, mm_iter, copy(w.b))

            return output
        end

        # save values from previous iterate
        copy!(sdata(w.b0), sdata(w.b))   # b0 = b
        copy!(sdata(w.xb0), sdata(w.xb)) # Xb0 = Xb
        loss = next_loss

        # now perform IHT step
        (mu, mu_step) = iht!(w, x, y, k, nstep=max_step, iter=mm_iter, pids=pids)

        # update residuals
        difference!(w.r, y, w.xb)
        mask!(w.r, mask_n, 0, zero(T))

        # use updated residuals to recompute the gradient on the GPU
        At_mul_B!(w.df, x, w.r, mask_n, v)

        # update objective
        next_loss = sumabs2(sdata(w.r)) / 2

        # guard against numerical instabilities
        # ensure that objective is finite
        # if not, throw error
        check_finiteness(next_loss)

        # track convergence
        the_norm    = chebyshev(w.b, w.b0)
        scaled_norm = (the_norm / ( norm(w.b0,Inf) + 1)) :: T
        converged   = scaled_norm < tol

        # output algorithm progress
        quiet || @printf("%d\t%d\t%3.7f\t%3.7f\t%3.7f\n", mm_iter, mu_step, mu, the_norm, next_loss)

        # check for convergence
        # if converged and in feasible set, then algorithm converged before maximum iteration
        # perform final computations and output return variables
        if converged

            # send elements below tol to zero
            threshold!(w.b, tol)

            # stop time
            mm_time = toq()

            # announce convergence
            !quiet && print_convergence(mm_iter, next_loss, mm_time)

            # these are output variables for function
            return IHTResults(mm_time, next_loss, mm_iter, copy(w.b))
        end

        # algorithm is unconverged at this point, so check descent property
        # if objective increases, then abort
        if next_loss > current_loss + tol
            !quiet && print_descent_error(mm_iter, loss, next_loss)
            throw(ErrorException("Descent failure!"))
        end
    end # end main loop
end # end function

"""
    iht_path(x::BEDFile, y, k, kernfile::String)

If supplied a `BEDFile` `x` and an OpenCL kernel file `kernfile` as an String, then `iht_path` will attempt to accelerate the calculation of the dense gradient `x' * (y - x*b)` in `L0_reg` with a GPU device.
"""
function iht_path{T <: Float, V <: DenseVector}(
    x        :: BEDFile{T},
    y        :: DenseVector{T},
    path     :: DenseVector{Int},
    kernfile :: String;
    pids     :: DenseVector{Int}  = procs(x),
    w        :: IHTVariables{T,V} = IHTVariables(x, y, 1),
    mask_n   :: Vector{Int}       = ones(Int,length(y)),
    v        :: PLINK.PlinkGPUVariables{T} = PLINK.PlinkGPUVariables(w.df, x, y, kernfile, mask_n),
    tol      :: T    = convert(T, 1e-4),
    max_iter :: Int  = 100,
    max_step :: Int  = 50,
    quiet    :: Bool = true
)
    # size of problem?
    n,p = size(x)

    # how many models will we compute?
    num_models = length(path)

    # also preallocate matrix to store betas
    betas = spzeros(T,p,num_models) # a matrix to store calculated models

    # compute the path
    @inbounds for i = 1:num_models

        # model size?
        q = path[i]

        # monitor progress
        quiet || print_with_color(:blue, "Computing model size $q.\n\n")

        # these arrays change in size from iteration to iteration
        # we must allocate them for every new model size
        update_variables!(w, x, q)

        # store projection of beta onto largest k nonzeroes in magnitude
        project_k!(w.b, q)

        # now compute current model
        output = L0_reg(x, y, q, kernfile, w=w, v=v, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet, pids=pids, mask_n=mask_n)

        # ensure that we correctly index the nonzeroes in b
        update_indices!(w.idx, output.beta)
        fill!(w.idx0, false)

        # put model into sparse matrix of betas
        betas[:,i] = sparsevec(output.beta)
    end

    return betas
end

"""
    one_fold(x::BEDFile, y, path, kernfile::String, folds, fold)

If supplied a `BEDFile` `x` and an OpenCL kernel file `kernfile` as an String, then `one_fold` will attempt to accelerate the calculation of the dense gradient `x' * (y - x*b)` in `L0_reg` with a GPU device.
"""
function one_fold{T <: Float}(
    x        :: BEDFile{T},
    y        :: DenseVector{T},
    path     :: DenseVector{Int},
    kernfile :: String,
    folds    :: DenseVector{Int},
    fold     :: Int;
    pids     :: Vector{Int} = procs(x),
    tol      :: T    = convert(T, 1e-4),
    max_iter :: Int  = 100,
    max_step :: Int  = 50,
    header   :: Bool = false,
    quiet    :: Bool = true
)
    # dimensions of problem
    n,p = size(x)

    # get list of available GPU devices
    # var device gets pointer to device indexed by variable devidx
    #device = cl.devices(:gpu)[devidx]

    # make vector of indices for folds
    test_idx = folds .== fold

    # train_idx is the vector that indexes the TRAINING set
    train_idx = !test_idx

    # how many indices are in test set?
    test_size = sum(test_idx)

    # GPU code requires Int variant of training indices, so do explicit conversion
    train_idx = convert(Vector{Int}, train_idx) :: Vector{Int}
    test_idx  = convert(Vector{Int}, test_idx)  :: Vector{Int}

    # compute the regularization path on the training set
    betas = iht_path(x, y, path, kernfile, max_iter=max_iter, quiet=quiet, max_step=max_step, mask_n=train_idx, pids=pids, tol=tol)

    # tidy up
    #gc()

    # preallocate vector for output
    myerrors = zeros(T, length(path))

    # allocate an index vector for b
    indices = falses(p)

    # allocate temporary arrays for the test set
    xb = SharedArray(T, n, pids=pids) :: SharedVector{T}
    b  = SharedArray(T, p, pids=pids) :: SharedVector{T}
    r  = SharedArray(T, n, pids=pids) :: SharedVector{T}

    # compute the mean out-of-sample error for the TEST set
    # do this for every computed model in regularization path
    for i = 1:size(betas,2)

        # pull ith model in dense vector format
        b2 = full(vec(betas[:,i]))

        # copy it into SharedArray b
        copy!(b,b2)

        # indices stores Boolean indexes of nonzeroes in b
        update_indices!(indices, b)

        # compute estimated response Xb with $(path[i]) nonzeroes
        #A_mul_B!(xb, x, b, indices, path[i], test_idx, pids=pids)
        A_mul_B!(xb, x, b, indices, path[i], test_idx)

        # compute residuals
        difference!(r, y, xb)

        # mask data from training set
        # training set consists of data NOT in fold:
        # r[folds .!= fold] = zero(Float64)
        mask!(r, test_idx, 0, zero(T))

        # compute out-of-sample error as squared residual averaged over size of test set
        myerrors[i] = sumabs2(r) / test_size / 2
    end

    return myerrors :: Vector{T}
end

# subroutine to calculate the approximate memory load of one fold on the GPU
function onefold_device_memload(x::BEDFile, wg_size::Int, y_chunks::Int; prec64::Bool = true)

    # floating point bytes multiplier depends on precision
    # prec64 = true -> use double precision (8 bytes per float)
    # prec64 = false -> use single precision (4 bytes per float)
    bytemult = ifelse(prec64, 8, 4)

    # get dimensions of problem
    n = x.n
    p = size(x,2)

    # start counting bytes
    # genotype array has bitstype Int8, one byte per component
    nx = length(x.x)

    # residual
    nr = bytemult * n

    # gradient, means, precisions
    ndf = 3 * bytemult * p

    # reduction vector
    nrv = bytemult * p * y_chunks

    # bitmask, which has bitstype Int32 (4 bytes per component)
    nbm = 4 * n

    # local memory footprint, just in case
    loc = bytemult * wg_size

    # total memory in Mb, rounded up to nearest int
    totmem = nx + ndf + nrv + nbm + loc

    return ceil(Int, totmem / 1024^2)
end

# subroutine to compute the number of folds that will fit on the GPU at one time
function compute_max_gpu_load(x::BEDFile, wg_size::Int, device::cl.Device; prec64::Bool = true)

    # number of chunks in residual
    y_chunks = div(x.n, wg_size) + (x.n % wg_size != 0 ? 1 : 0)

    # total available memory on current device
    gpu_memtot = ceil(Int, device[:global_mem_size] / 1024^2)

    # memory load of one CV fold on current device
    onefold_mem = onefold_device_memload(x,wg_size,y_chunks, prec64=prec64)

    # how many folds could we fit on the current GPU?
    max_folds = div(gpu_memtot, onefold_mem)

    return max_folds
end

"""
    pfold(xfile, xtfile, x2file,yfile, meanfile, precfile,path,kernfile,folds,q [, pids=procs(), devindices=ones(Int,q])

This function is the parallel execution kernel in `cv_iht()`. It is not meant to be called outside of `cv_iht()`.
It will distribute `q` crossvalidation folds across the processes supplied by the optional argument `pids` and call `one_fold()` for each fold.
Each fold will use the GPU device indexed by its corresponding component of the optional argument `devindices` to compute a regularization path given by `path`.
`pfold()` collects the vectors of MSEs returned by calling `one_fold()` for each process, reduces them, and returns their average across all folds.
"""
function pfold(
    T          :: Type,
    xfile      :: String,
    xtfile     :: String,
    x2file     :: String,
    yfile      :: String,
    meanfile   :: String,
    precfile   :: String,
    path       :: DenseVector{Int},
    kernfile   :: String,
    folds      :: DenseVector{Int},
    q          :: Int;
    devindices :: Vector{Int} = ones(Int,q),
    pids       :: Vector{Int} = procs(),
    max_iter   :: Int  = 100,
    max_step   :: Int  = 50,
    quiet      :: Bool = true,
    header     :: Bool = false
)

    # ensure correct type
    T <: Float || throw(ArgumentError("Argument T must be either Float32 or Float64"))

    # how many CPU processes can pfold use?
    np = length(pids)

    # report on CPU processes
    quiet || println("pfold: np = ", np)
    quiet || println("pids = ", pids)

    # set up function to share state (indices of folds)
    i = 1
    nextidx() = (idx=i; i+=1; idx)

    # preallocate cell array for results
    results = SharedArray(T, (length(path),q), pids=pids) :: SharedMatrix{T}

    # master process will distribute tasks to workers
    # master synchronizes results at end before returning
    @sync begin

        # loop over all workers
        for worker in pids

            # exclude process that launched pfold, unless only one process is available
            if worker != myid() || np == 1

                # asynchronously distribute tasks
                @async begin
                    while true

                        # grab next fold
                        current_fold = nextidx()

                        # if current fold exceeds total number of folds then exit loop
                        current_fold > q && break

                        # grab index of GPU device
                        devidx = devindices[current_fold]

                        # report distribution of fold to worker and device
                        quiet || print_with_color(:blue, "Computing fold $current_fold on worker $worker and device $devidx.\n\n")

                        # launch job on worker
                        # worker loads data from file paths and then computes the errors in one fold
                        r = remotecall_fetch(worker) do
                                pids = [worker]
                                x = BEDFile(T, xfile, xtfile, x2file, meanfile, precfile, pids=pids, header=header) :: BEDFile{T}
                                y = SharedArray(abspath(yfile), T, (x.geno.n,), pids=pids) :: SharedVector{T}

                                one_fold(x, y, path, kernfile, folds, current_fold, max_iter=max_iter, max_step=max_step, quiet=quiet, pids=pids)
                        end # end remotecall_fetch()
                        setindex!(results, r, :, current_fold)
                    end # end while
                end # end @async
            end # end if
        end # end for
    end # end @sync

    # return reduction (row-wise sum) over results
    return (vec(sum(results, 2) ./ q)) :: Vector{T}
end

# default type for pfold is Float64
pfold(xfile::String, xtfile::String, x2file::String, yfile::String, meanfile::String, precfile::String, path::DenseVector{Int}, kernfile::String, folds::DenseVector{Int}, q::Int; devindices::Vector{Int}=ones(Int,q), pids::Vector{Int}=procs(), max_iter::Int=100, max_step::Int =50, quiet::Bool=true, header::Bool=false) = pfold(Float64, xfile, xtfile, x2file, yfile, meanfile, precfile, path, kernfile, folds, q, devindices=devindices, pids=pids, max_iter=max_iter, max_step=max_step, quiet=quiet, header=header)

"""
    cv_iht(xfile, xtfile, x2file, yfile, meanfile, precfile, kernfile)
This variant of `cv_iht()` uses a GPU to perform `q`-fold crossvalidation with a `BEDFile` object loaded by `xfile`, `xtfile`, and `x2file`,
with column means stored in `meanfile` and column precisions stored in `precfile`.
The continuous response is stored in the binary file `yfile`.
The calculations employ GPU acceleration by calling OpenCL kernels from `kernfile` with workgroup size `wg_size`.
"""
function cv_iht(
    T        :: Type,
    xfile    :: String,
    xtfile   :: String,
    x2file   :: String,
    yfile    :: String,
    meanfile :: String,
    precfile :: String,
    kernfile :: String;
    q        :: Int = cv_get_num_folds(3,5),
    path     :: DenseVector{Int} = begin
           # find p from the corresponding BIM file, then make path
            bimfile = xfile[1:(endof(xfile)-3)] * "bim"
            p       = countlines(bimfile)
            collect(1:min(20,p))
            end,
    folds    :: DenseVector{Int} = begin
           # find n from the corresponding FAM file, then make folds
            famfile = xfile[1:(endof(xfile)-3)] * "fam"
            n       = countlines(famfile)
            cv_get_folds(n, q)
            end,
    pids     :: Vector{Int} = procs(),
    tol      :: Float = convert(T, 1e-4),
    max_iter :: Int   = 100,
    max_step :: Int   = 50,
    wg_size  :: Int   = 512,
    quiet    :: Bool  = true,
    header   :: Bool  = false
)
    # enforce type
    T <: Float || throw(ArgumentError("Argument T must be either Float32 or Float64"))

    # how many elements are in the path?
    num_models = length(path)

    # how many GPU devices are available to us?
    devs = cl.devices(:gpu)
    ndev = length(devs)

    # how many folds can we fit on a GPU at once?
    # count one less per GPU device, just in case
#   max_folds = zeros(Int, ndev)
#   @inbounds for i = 1:ndev
#       max_folds[i] = max(compute_max_gpu_load(x, wg_size, devs[i], prec64 = true) - 1, 0)
#   end

    # how many rounds of folds do we need to schedule?
#   fold_rounds = zeros(Int, ndev)
#   @inbounds for i = 1:ndev
#       fold_rounds[i] = div(q, max_folds[i]) + (q% max_folds[i] != 0 ? 1 : 0)
#   end

    # assign index of a GPU device for each fold
    # default is first GPU device (devidx = 1)
    devindices = ones(Int, q)
#   @inbounds for i = 1:q
#       devindices[i] += i % ndev
#   end

    # want to compute a path for each fold
    # the folds are computed asynchronously
    # only use the worker processes
    mses = pfold(T, xfile, xtfile, x2file, yfile, meanfile, precfile, path, kernfile, folds, q, max_iter=max_iter, max_step=max_step, quiet=quiet, devindices=devindices, pids=pids, header=header)

    # what is the best model size?
    k = path[indmin(mses)] :: Int

    # print results
    !quiet && print_cv_results(mses, path, k)

    # load data on *all* processes
    x = BEDFile(T, xfile, xtfile, x2file, meanfile, precfile, header=header, pids=pids) :: BEDFile{T}
    y = SharedArray(abspath(yfile), T, (x.geno.n,), pids=pids) :: SharedVector{T}
    w = IHTVariables(x, y, k)
    m = ones(Int, size(y))
    v = PLINK.PlinkGPUVariables(w.df, x, y, kernfile, m)

    # first use L0_reg to extract model
    output = L0_reg(x, y, k, kernfile, w=w, mask_n=m, v=v, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet, pids=pids)

    # which components of beta are nonzero?
    inferred_model = output.beta .!= zero(T)
    bidx = find(inferred_model)

    # allocate the submatrix of x corresponding to the inferred model
    x_inferred = zeros(T, x.geno.n, sum(inferred_model))
    decompress_genotypes!(x_inferred, x, inferred_model)

    # refit best model
    b, bidx = refit_iht(x_inferred, y, k, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet)

    bids = prednames(x)[bidx]
    return IHTCrossvalidationResults{T}(mses, sdata(path), b, bidx, k, bids)
end

# default type for cv_iht is Float64
cv_iht(xfile::String, x2file::String, yfile::String, meanfile::String, precfile::String, kernfile::String; q::Int = cv_get_num_folds(3,5), path::DenseVector{Int} = begin bimfile=xfile[1:(endof(xfile)-3)] * "bim"; p=countlines(bimfile); collect(1:min(20,p)) end, folds::DenseVector{Int} = begin famfile=xfile[1:(endof(xfile)-3)] * "fam"; n=countlines(famfile); cv_get_folds(n, q) end, pids::Vector{Int}=procs(), tol::Float64=1e-4, max_iter::Int=100, max_step::Int=50, quiet::Bool=true, header::Bool=false) = cv_iht(Float64, xfile, x2file, yfile, meanfile, precfile, kernfile, path=path, folds=folds, q=q, pids=pids, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet, header=header)

"""
    pfold(xfile, x2file, yfile, path, kernfile, folds, q [, pids=procs(), devindices=ones(Int,q])

`pfold` can also be called without the transposed genotype files or the mean/precision files.
In this case, it will attempt to precompute them before calling `one_fold`.
"""
function pfold(
    T          :: Type,
    xfile      :: String,
    x2file     :: String,
    yfile      :: String,
    path       :: DenseVector{Int},
    kernfile   :: String,
    folds      :: DenseVector{Int},
    q          :: Int;
    devindices :: Vector{Int} = ones(Int,q),
    pids       :: Vector{Int} = procs(),
    max_iter   :: Int  = 100,
    max_step   :: Int  = 50,
    quiet      :: Bool = true,
    header     :: Bool = false
)

    # ensure correct type
    T <: Float || throw(ArgumentError("Argument T must be either Float32 or Float64"))

    # how many CPU processes can pfold use?
    np = length(pids)

    # report on CPU processes
    quiet || println("pfold: np = ", np)
    quiet || println("pids = ", pids)

    # set up function to share state (indices of folds)
    i = 1
    nextidx() = (idx=i; i+=1; idx)

    # preallocate cell array for results
    results = SharedArray(T, (length(path),q), pids=pids) :: SharedMatrix{T}

    # master process will distribute tasks to workers
    # master synchronizes results at end before returning
    @sync begin

        # loop over all workers
        for worker in pids

            # exclude process that launched pfold, unless only one process is available
            if worker != myid() || np == 1

                # asynchronously distribute tasks
                @async begin
                    while true

                        # grab next fold
                        current_fold = nextidx()

                        # if current fold exceeds total number of folds then exit loop
                        current_fold > q && break

                        # grab index of GPU device
                        devidx = devindices[current_fold]

                        # report distribution of fold to worker and device
                        quiet || print_with_color(:blue, "Computing fold $current_fold on worker $worker and device $devidx.\n\n")

                        # launch job on worker
                        # worker loads data from file paths and then computes the errors in one fold
                        r = remotecall_fetch(worker) do
                                pids = [worker]
                                x = BEDFile(T, xfile, x2file, pids=pids, header=header) :: BEDFile{T}
                                y = SharedArray(abspath(yfile), T, (x.geno.n,), pids=pids) :: SharedVector{T}
                                one_fold(x, y, path, kernfile, folds, current_fold, max_iter=max_iter, max_step=max_step, quiet=quiet, pids=pids)
                        end # end remotecall_fetch()
                        setindex!(results, r, :, current_fold)
                    end # end while
                end # end @async
            end # end if
        end # end for
    end # end @sync

    # return reduction (row-wise sum) over results
    return (vec(sum(results, 2) ./ q)) :: Vector{T}
end

# default type for pfold is Float64
pfold(xfile::String, x2file::String, yfile::String, path::DenseVector{Int}, kernfile::String, folds::DenseVector{Int}, q::Int; devindices::Vector{Int}=ones(Int,q), pids::Vector{Int}=procs(), max_iter::Int=100, max_step::Int =50, quiet::Bool=true, header::Bool=false) = pfold(Float64, xfile, x2file, yfile, path, kernfile, folds, q, devindices=devindices, pids=pids, max_iter=max_iter, max_step=max_step, quiet=quiet, header=header)

"""
    cv_iht(xfile,x2file,yfile,path,kernfile,folds,q [, pids=procs()])

If `cv_iht` is called without filepaths to the transposed genotype data or the means and precisions, then it will attempt to precompute them.
"""
function cv_iht(
    T        :: Type,
    xfile    :: String,
    x2file   :: String,
    yfile    :: String,
    kernfile :: String;
    q        :: Int = cv_get_num_folds(3,5),
    path     :: DenseVector{Int} = begin
           # find p from the corresponding BIM file, then make path
            bimfile = xfile[1:(endof(xfile)-3)] * "bim"
            p       = countlines(bimfile)
            collect(1:min(20,p))
            end,
    folds    :: DenseVector{Int} = begin
           # find n from the corresponding FAM file, then make folds
            famfile = xfile[1:(endof(xfile)-3)] * "fam"
            n       = countlines(famfile)
            cv_get_folds(n, q)
            end,
    pids     :: Vector{Int} = procs(),
    tol      :: Float = convert(T, 1e-4),
    max_iter :: Int   = 100,
    max_step :: Int   = 50,
    wg_size  :: Int   = 512,
    quiet    :: Bool  = true,
    header   :: Bool  = false
)
    # enforce type
    T <: Float || throw(ArgumentError("Argument T must be either Float32 or Float64"))

    # how many elements are in the path?
    num_models = length(path)

    # how many GPU devices are available to us?
    devs = cl.devices(:gpu)
    ndev = length(devs)

    # how many folds can we fit on a GPU at once?
    # count one less per GPU device, just in case
#   max_folds = zeros(Int, ndev)
#   @inbounds for i = 1:ndev
#       max_folds[i] = max(compute_max_gpu_load(x, wg_size, devs[i], prec64 = true) - 1, 0)
#   end

    # how many rounds of folds do we need to schedule?
#   fold_rounds = zeros(Int, ndev)
#   @inbounds for i = 1:ndev
#       fold_rounds[i] = div(q, max_folds[i]) + (q% max_folds[i] != 0 ? 1 : 0)
#   end

    # assign index of a GPU device for each fold
    # default is first GPU device (devidx = 1)
    devindices = ones(Int, q)
#   @inbounds for i = 1:q
#       devindices[i] += i % ndev
#   end

    # want to compute a path for each fold
    # the folds are computed asynchronously
    # only use the worker processes
    mses = pfold(T, xfile, x2file, yfile, path, kernfile, folds, q, max_iter=max_iter, max_step=max_step, quiet=quiet, devindices=devindices, pids=pids, header=header)

    # what is the best model size?
    k = path[indmin(mses)] :: Int

    # print results
    !quiet && print_cv_results(mses, path, k)

    # recompute ideal model
    # first load data on *all* processes
    x = BEDFile(T, xfile, xtfile, x2file, header=header, pids=pids) :: BEDFile{T}
    y = SharedArray(abspath(yfile), T, (x.geno.n,), pids=pids) :: SharedVector{T}
    w = IHTVariables(x, y, k)
    m = ones(Int, size(y))
    v = PLINK.PlinkGPUVariables(w.df, x, y, kernfile, m)

    # first use L0_reg to extract model
    output = L0_reg(x, y, k, kernfile, w=w, mask_n=m, v=v, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet, pids=pids)

    # which components of beta are nonzero?
    inferred_model = output.beta .!= zero(T)
    bidx = find(inferred_model)

    # allocate the submatrix of x corresponding to the inferred model
    x_inferred = zeros(T, x.geno.n, sum(inferred_model))
    decompress_genotypes!(x_inferred, x, inferred_model)

    # refit the best model
    b, bidx = refit_iht(x_inferred, y, k, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet)

    bids = prednames(x)[bidx]
    return IHTCrossvalidationResults{T}(mses, sdata(path), b, bidx, k, bids)
end

# default type for cv_iht is Float64
cv_iht(xfile::String, x2file::String, yfile::String, kernfile::String; q::Int = cv_get_num_folds(3,5), path::DenseVector{Int} = begin bimfile=xfile[1:(endof(xfile)-3)] * "bim"; p=countlines(bimfile); collect(1:min(20,p)) end, folds::DenseVector{Int} = begin famfile=xfile[1:(endof(xfile)-3)] * "fam"; n=countlines(famfile); cv_get_folds(n, q) end, pids::Vector{Int}=procs(), tol::Float64=1e-4, max_iter::Int=100, max_step::Int=50, quiet::Bool=true, header::Bool=false) = cv_iht(Float64, xfile, x2file, yfile, kernfile, path=path, folds=folds, q=q, pids=pids, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet, header=header)


text:
module BibParser

export parse_file

include("utils.jl")

# BibTeX module
include("bibtex.jl")
import .BibTeX

# CSL-JSON module
include("csl.jl")
import .CSL

"""
    parse_file(path::String; parser::Symbol = :BibTeX)
Parse a bibliography file. Default to BibTeX format. No other options available yet (CSL-JSON coming soon).
"""
function parse_file(path::String; parser::Symbol = :BibTeX)
    if parser == :BibTeX
        BibTeX.parse_file(path)
    end
end

"""
    parse_entry(entry::String; parser::Symbol = :BibTeX)
Parse a string entry. Default to BibTeX format. No other options available yet (CSL-JSON coming soon).
"""
function parse_entry(entry::String; parser::Symbol = :BibTeX)
    if parser == :BibTeX
        BibTeX.parse_bibtex(entry)
    end
end

end # module

ext: jl


text:
function pooled_pava_isotonic_regression!(y::Vector{Float64}, weights::Vector{Float64})

    n = length(y)
    if n <= 1
        return y
    end
    if n != length(weights)
        throw(DimensionMismatch("Lengths of values and weights mismatch"))
    end

    j = 1
    S = Dict(0 => 0, 1 => 1)
    ydash = Dict(1 => y[1])
    wdash = Dict(1 => weights[1])
    @inbounds begin
        for i in 2 : n
            j = j + 1
            ydash[j] = y[i]
            wdash[j] = weights[i]
            while j > 1 && ydash[j] < ydash[j-1]
                ydash[j-1] = (wdash[j] * ydash[j] + wdash[j-1] * ydash[j-1]) /
                   (wdash[j] + wdash[j-1])
                wdash[j-1] = wdash[j] + wdash[j-1]
                j = j-1
            end
            S[j] = i
        end
        for k in 1 : j
            for l in S[k-1] + 1 : S[k]
                y[l] = ydash[k]
            end
        end
    end
    return y
end

pooled_pava_isotonic_regression!(y::Vector{Float64}) = pooled_pava_isotonic_regression!(y, ones(size(y, 1)))

# non-mutating versions
pooled_pava_isotonic_regression(y::Vector{Float64}, weights::Vector{Float64}) = pooled_pava_isotonic_regression!(copy(y), weights)
pooled_pava_isotonic_regression(y::Vector{Float64}) = pooled_pava_isotonic_regression(y, ones(size(y,1)))


text:
#!/usr/bin/env julia

using Test, SafeTestsets

@time begin
@time @safetestset "Project.toml append values to package.json" begin include("project_toml_append_values_test.jl") end
@time @safetestset "Webpack integration to existing project" begin include("project_integration_test.jl") end
end


text:
module train
using DataFrames
using ArgParse

##function export list
#export train

function import_data()
    s = ArgParseSettings()
    @add_arg_table s begin
	"--source", "-s"
	help = "the relative path to the training set. defaults to ../data/train.csv"
	arg_type = String
	default = "../data/train.csv"
    end

    passengers = readtable(s.source)
    (num, len) = size(passengers)
    println("%d passengers", num)
end

end

max_stars_repo_name: kfuruya/titanic

max_issues_repo_name: kfuruya/titanic

max_forks_repo_name: kfuruya/titanic


text:
<!DOCTYPE html>

    <title>{{ title ?? 'Angel' }}</title>

    <link href="https://fonts.googleapis.com/css?family=Lato:100" rel="stylesheet" type="text/css">
    <link rel="stylesheet" href="/css/site.css">
    <link rel="icon" href="/images/favicon.png">

        <block name="content"></block>

max_stars_repo_name: misterdev/angel

max_issues_repo_name: misterdev/angel

max_forks_repo_name: misterdev/angel


text:
"""
    AbstractSpecialPolynomial{T,X}

An abstract type to distinguish the different polynomial types in this package.

The concrete types specify different bases for the space of polynomials of degree `n` or less.

This package includes:

* several classic orthogonal polynomials.
* Newton and Lagrange interpolating polynomials
* Bernstein polynomials

As many of the methods for the base `Polynomials` class are directly coded if possible, but quite a few
depend on conversion to the base `Polynomial` type (which uses the standard polynomial basis).

"""
abstract type AbstractSpecialPolynomial{T,X} <: Polynomials.AbstractPolynomial{T,X} end

# polynomial like Vector{T} with variable
Base.eltype(::Type{<:AbstractSpecialPolynomial{T}}) where {T} = T
Base.eltype(::Type{<:AbstractSpecialPolynomial}) = Float64

# to strip off type parameters. See Polynomials.constructorof
# We want ⟒(P{α,T,N}) = P{α} where {T, N}x
@generated function constructorof(::Type{T}) where {T}
    getfield(parentmodule(T), nameof(T))
end
⟒(P::Type{<:AbstractSpecialPolynomial}) = constructorof(P)
⟒(::Type{<:Polynomial}) = Polynomial

## Display

function unicode_subscript(io, j)
    a = ("⁻", "", "", "₀", "₁", "₂", "₃", "₄", "₅", "₆", "₇", "₈", "₉")
    for i in string(j)
        print(io, a[Int(i) - 44])
    end
end

function Polynomials.showterm(
    io::IO,
    ::Type{P},
    pj::T,
    var,
    j,
    first::Bool,
    mimetype,
) where {N,T,P<:AbstractSpecialPolynomial}
    iszero(pj) && return false
    !first && print(io, " ")

    if Polynomials.hasneg(T) && Polynomials.isneg(pj)
        print(io, "- ")
        pj = -pj
    else
        !first && print(io, "+ ")
    end

    Polynomials.printcoefficient(io, pj, 1, mimetype) # j != 0; want to have ()
    print(io, "⋅")
    print(io, basis_symbol(P))

    unicode_subscript(io, j)
    print(io, "($(var))")
    return true
end

# Conversion is *possible* when multiplication  in `P` is defined
# but otherwise,  it must be handled specially
# function Base.convert(::Type{P}, q::Q) where {P <: AbstractSpecialPolynomial, Q <: AbstractSpecialPolynomial}
#     convert(P, convert(Polynomial, q))
# end

# return the domain as a tuple, not an interval object
function Base.extrema(::Type{P}) where {P<:AbstractSpecialPolynomial}
    dom = domain(P)
    first(dom), last(dom)
end
Base.extrema(p::P) where {P<:AbstractSpecialPolynomial} = extrema(P)

## Polynomial operations

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

function Polynomials.derivative(p::P, order::Integer=1) where {P<:AbstractSpecialPolynomial}
    T = eltype(one(P))
    q = convert(Polynomial{T}, p)
    convert(⟒(P), derivative(q, order))
end

function Polynomials.integrate(p::P) where {P<:AbstractSpecialPolynomial}
    q = convert(Polynomial, p)
    integrate(q)
end

function Base.divrem(num::P, den::P) where {P<:AbstractSpecialPolynomial}
    p1 = convert(Polynomial, num)
    p2 = convert(Polynomial, den)
    q, r = divrem(p1, p2)
    convert.(P, (q, r))
end

function Polynomials.companion(p::P) where {P<:AbstractSpecialPolynomial}
    companion(convert(Polynomial, p))
end

function Polynomials.vander(
    ::Type{P},
    x::AbstractVector{T},
    n::Integer,
) where {P<:AbstractSpecialPolynomial,T}
    N = length(x) - 1
    R = typeof(one(T) / one(T))
    V = zeros(R, N + 1, n + 1)

    for j in 0:n
        bj = Polynomials.basis(P, j)
        for i in 0:N
            V[i + 1, j + 1] = bj(x[i + 1])
        end
    end

    V
end


text:
module DummyPolyModule

using LinearAlgebra
using MathOptInterface
const MOI = MathOptInterface
using JuMP
using PolyJuMP
using MultivariatePolynomials
import MultivariateBases
const MB = MultivariateBases
using SemialgebraicSets

struct NonNeg{BT <: MB.AbstractPolynomialBasis,
              DT <: SemialgebraicSets.AbstractSemialgebraicSet,
              MT <: MultivariatePolynomials.AbstractMonomial,
              MVT <: AbstractVector{MT}} <: MOI.AbstractVectorSet
    basis::Type{BT}
    domain::DT
    monomials::MVT
    kwargs
end
MOI.dimension(set::NonNeg) = length(set.monomials)
Base.copy(set::NonNeg) = set

struct DummyNonNeg <: PolyJuMP.PolynomialSet end

JuMP.reshape_set(::NonNeg, ::PolyJuMP.PolynomialShape) = DummyNonNeg()

struct DummyNonNegBridge{T,F} <: MOIB.Constraint.AbstractBridge end
function PolyJuMP.bridges(
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{<:NonNeg},
)
    return [DummyNonNegBridge]
end
function MOI.Bridges.added_constrained_variable_types(::Type{<:DummyNonNegBridge})
    return Tuple{Type}[]
end
function MOI.Bridges.added_constraint_types(::Type{DummyNonNegBridge{T,F}}) where {T, F}
    return Tuple{Type,Type}[(F, MOI.EqualTo{T})]
end
function MOI.Bridges.Constraint.concrete_bridge_type(
    ::Type{<:DummyNonNegBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    ::Type{<:NonNeg}) where {T}
    # This tests that `T` matches the coefficient type of `F`
    # as otherwise it would error.
    G = MOI.Utilities.promote_operation(-, T, F, MOI.VectorOfVariables)
    H = MOI.Utilities.scalar_type(G)
    return DummyNonNegBridge{T,H}
end

function JuMP.moi_set(cone::DummyNonNeg,
                      monos::AbstractVector{<:AbstractMonomial};
                      domain::AbstractSemialgebraicSet=FullSpace(),
                      basis=MB.MonomialBasis, kwargs...)
    return NonNeg(basis, domain, monos, kwargs)
end

function JuMP.build_constraint(_error::Function, p::AbstractPolynomialLike,
                               s::DummyNonNeg; kwargs...)
    coefs = PolyJuMP.non_constant_coefficients(p)
    monos = monomials(p)
    set = JuMP.moi_set(s, monos; kwargs...)
    return PolyJuMP.bridgeable(
        JuMP.VectorConstraint(coefs, set, PolyJuMP.PolynomialShape(monos)),
        JuMP.moi_function_type(typeof(coefs)),
        typeof(set),
    )
end

struct MatrixPolynomialShape{MT <: AbstractMonomial,
                             MVT <: AbstractVector{MT}} <: JuMP.AbstractShape
    side_dimension::Int
    monomials::Matrix{MVT}
end

function JuMP.reshape_vector(x::Vector,
                             shape::MatrixPolynomialShape{MT}) where {MT}
    n = shape.side_dimension
    p = Matrix{polynomialtype(MT, eltype(x))}(undef, n, n)
    k = 0
    for j in 1:n
        for i in 1:n
            m = length(shape.monomials[i, j])
            p[i, j] = polynomial(x[k .+ (1:m)], shape.monomials[i, j])
            k += m
        end
    end
    @assert length(x) == k
    return p
end

struct PosDefMatrix{BT <: MB.AbstractPolynomialBasis,
                    DT <: SemialgebraicSets.AbstractSemialgebraicSet,
                    MT <: MultivariatePolynomials.AbstractMonomial,
                    MVT <: AbstractVector{MT}} <: MOI.AbstractVectorSet
    basis::Type{BT}
    domain::DT
    monomials::Matrix{MVT}
    kwargs
end
MOI.dimension(set::PosDefMatrix) = sum(length, set.monomials)
Base.copy(set::PosDefMatrix) = set

struct DummyPosDefMatrix <: PolyJuMP.PolynomialSet end

function JuMP.reshape_set(::PosDefMatrix, ::MatrixPolynomialShape)
    return DummyPosDefMatrix()
end
function JuMP.moi_set(::DummyPosDefMatrix,
                      monos::Matrix{<:AbstractVector{<:AbstractMonomial}};
                      domain::AbstractSemialgebraicSet=FullSpace(),
                      basis=MB.MonomialBasis, kwargs...)
    return PosDefMatrix(basis, domain, monos, kwargs)
end

function JuMP.build_constraint(_error::Function,
                               p::Matrix{<:AbstractPolynomialLike},
                               s::DummyPosDefMatrix; kwargs...)
    n = LinearAlgebra.checksquare(p)
    # TODO we should use `non_constant_coefficients_type` once it exists
    coefs = coefficienttype(p[1, 1])[]
    for j in 1:n
        for i in 1:n
            append!(coefs, coefficients(p[i, j]))
        end
    end
    monos = monomials.(p)
    set = JuMP.moi_set(s, monos; kwargs...)
    return JuMP.VectorConstraint(coefs, set, MatrixPolynomialShape(n, monos))
end

# TODO the function in JuMP should not require the eltype to be
#      `AbstractJuMPScalar` so that we don't have to define this
# These methods are just copy-paste from JuMP/src/print.jl
function JuMP.function_string(::MIME"text/plain",
                              A::AbstractMatrix{<:AbstractPolynomialLike})
    str = sprint(show, MIME"text/plain"(), A)
    lines = split(str, '\n')
    # We drop the first line with the signature "m×n Array{...}:"
    lines = lines[2:end]
    # We replace the first space by an opening `[`
    lines[1] = '[' * lines[1][2:end]
    for i in 1:length(lines)
        lines[i] = lines[i] * (i == length(lines) ? ']' : ';')
    end
    return join(lines, '\n')
end
function JuMP.function_string(print_mode::MIME"text/latex",
                              A::AbstractMatrix{<:AbstractPolynomialLike})
    str = sprint(show, MIME"text/plain"(), A)
    str = "\\begin{bmatrix}\n"
    for i in 1:size(A, 1)
        line = ""
        for j in 1:size(A, 2)
            if j != 1
                line *= " & "
            end
            if A isa Symmetric && i > j
                line *= "\\cdot"
            else
                line *= JuMP.function_string(print_mode, A[i, j])
            end
        end
        str *= line * "\\\\\n"
    end
    return str * "\\end{bmatrix}"
end

function setdefaults!(data::PolyJuMP.Data)
    PolyJuMP.setdefault!(data, PolyJuMP.NonNegPoly, DummyNonNeg)
    PolyJuMP.setdefault!(data, PolyJuMP.PosDefPolyMatrix, DummyPosDefMatrix)
end

end


text:
# ------------------ CMAX agent ------------------------
struct MountainCarCMAXAgent
    mountaincar::MountainCar
    cmax_planner::MountainCarCMAXPlanner
end

function run(agent::MountainCarCMAXAgent; max_steps = 1e5)
    state = init(agent.mountaincar; cont = true)
    num_steps = 0
    while !checkGoal(agent.mountaincar, state) && num_steps < max_steps
        num_steps += 1
        best_action, info = act(agent.cmax_planner, cont_state_to_disc(agent.mountaincar, state))
        new_state, cost = step(agent.mountaincar, state, best_action, true_params)
        # Check if sim prediction matches true next state
        sim_state, _ = step(
            agent.cmax_planner.planner.mountaincar,
            state,
            best_action,
            agent.cmax_planner.planner.params,
        )
        if sim_state != new_state
            # println("Discrepancy detected")
            addDiscrepancy!(
                agent.cmax_planner,
                cont_state_to_disc(agent.mountaincar, state),
                best_action
            )
            # replan
            _, info = act(agent.cmax_planner, cont_state_to_disc(agent.mountaincar, state))
        end
        updateResiduals!(agent.cmax_planner, info)
        state = new_state
    end
    if num_steps < max_steps
        println("Reached goal in ", num_steps, " steps")
    else
        println("Did not reach goal")
    end
    num_steps
end


text: ROOT_FILEPATH = "/home/member/gail-driver/"


text:
function SensitivityMultiShot{T1<:Complex,T2<:Int,T3<:AbstractFloat}(IN::Array{T1,1},adj::Bool;R::SparseMatrixCSC{T2,T2}=sparse(ones(T2,1,1)),H::Base.SparseArrays.UMFPACK.UmfpackLU{T1,T2}=lufact(sparse(ones(T1,1,1))),U::Array{T1,2}=im*ones(T1,1,1),A::SparseMatrixCSC{T1,T2}=sparse(ones(T1,1,1)),w::T3=1.)

# This function is the forward and adjoint sensitivity/jacobian operator
# for multiple shots.
#
# INPUTS:     IN        - Data residuals for adjoint; Model perturbation for forward
#             adj       - true for adjoint; false for forward
#             R         - Restriction operator
#             H         - LU factorization of the Helmholtz operator
#             U         - Predicted forward-propagated wavefield for each shot (stored as a matrix with dimensions nz*nx,ns)
#             A         - Diagonal matrix containing the complex-valued attenuation coefficients
#             w         - Angular frequency
#
# OUTPUTS:    OUT       - Model perturbation for adjoint; Data residuals for forward

    ng = size(R,1)
    ns = size(U,2)
    n = size(A,2)

    if adj==false
        OUT = zeros(T1,ng*ns)
        for ishot = 1:ns
            OUT[(ishot-1)*ng+1:ishot*ng] = Sensitivity1Shot(IN,adj,R=R,H=H,U=U[:,ishot],A=A,w=w)
        end
    end

    if adj==true
        OUT = zeros(T1,n)
        for ishot = 1:ns
            OUT += Sensitivity1Shot(IN[(ishot-1)*ng+1:ishot*ng],adj,R=R,H=H,U=U[:,ishot],A=A,w=w)
        end
    end

    return OUT

end

function Sensitivity1Shot{T1<:Complex,T2<:Int,T3<:AbstractFloat}(IN::Array{T1,1},adj::Bool;R::SparseMatrixCSC{T2,T2}=sparse(ones(T2,1,1)),H::Base.SparseArrays.UMFPACK.UmfpackLU{T1,T2}=lufact(sparse(ones(T1,1,1))),U::Array{T1,1}=im*ones(T1,1),A::SparseMatrixCSC{T1,T2}=sparse(ones(T1,1,1)),w::T3=1.)

# This function is the forward and adjoint sensitivity/jacobian operator
# for a single shot.
#
# INPUTS:     IN        - Data residual for adjoint; Model perturbation for forward
#             adj       - true for adjoint; false for forward
#             R         - Restriction operator
#             H         - LU factorization of the Helmholtz operator
#             U         - Predicted forward-propagated wavefield
#             A         - Diagonal matrix containing the complex-valued attenuation coefficients
#             w         - Angular frequency
#
# OUTPUTS:    OUT       - Model perturbation for adjoint; Data residual for forward

    if adj==false
        OUT = -w^2*R*(H\(spdiagm(U)*A*IN))
    end

    if adj==true
        OUT = -w^2*A'*spdiagm(U)'*(H'\(R'*IN))
    end

    return OUT

end


text:
# Hello World
# variable substitution with the $ operator
word = "world"
output = "Hello $word !"
println(output)

# tab autocomplete of LateX characters
ℵ = "world"
output = "Hello $ℵ, again!"
println(output)

# this extends to function evaluations as well
α = 2
println("The sine of $α is $(sin(α))")

# sidenote:
# so if printing to a file you can directly direference the evalutioan
# without passing it to an array which can then be printed.

typeof(α)
typeof(sin)
methods(sin)

# built in type promotion
check = int8(10)
typemax(Int8)
typeof(check)
check = check *1232
typeof(check)

# defining numerical functions
# looks and feels very much like real math

f(x) = 2*x + x^2 + 5
f(3.4)

g(x, y) = 34 + 2*x^2 + y^3
g(2, 3)

# passing function results to arrays
arr = Float64[]
for x in 1.3:0.01:3.9 push!(arr, f(x)) end
println(arr)
for x in 1.3:0.01:9.9 push!(arr, f(x)) end
println(arr)

arr2 = [f(i) for i in [1:2:10, 15]]
println(arr2)

using ASCIIPlots
lineplot([1.3:0.01:3.9], [1.3:0.01:3.9])
lineplot([1:2:10], [f(i) for i in [1:2:10]])

# working with imaginary numbers
a = 23
b = (1 + 34im)
c = a+b

typeof(a)
typeof(b)
typeof(c)

# now c can be dereferenced into real and imaginary parts
d = c.re
e = c.im

typeof(d)
typeof(e)

d = c*conj(c)
typeof(d)

max_stars_repo_name: thehackerwithin/swinburne

max_issues_repo_name: thehackerwithin/swinburne

max_forks_repo_name: thehackerwithin/swinburne


text:
@testset "Formulas" begin

    using StreamModels: sort_terms!, parse!

    @testset "Formula parsing" begin

        @testset "Associative property" begin
            @test parse!(:(a+(b+c))) == parse!(:(a+b+c))
            @test parse!(:((a+b)+c)) == parse!(:(a+b+c))
            @test parse!(:(a&(b&c))) == parse!(:(a&b&c))
            @test parse!(:((a&b)&c)) == parse!(:(a&b&c))
        end

        @testset "Distributive property" begin
            @test parse!(:(a & (b+c))) == parse!(:(a&b + a&c))
            @test parse!(:((a+b) & c)) == parse!(:(a&c + b&c))
            @test parse!(:((a+b) & (c+d))) == parse!(:(a&c + a&d + b&c + b&d))
            @test parse!(:(a & (b+c) & d)) == parse!(:(a&b&d + a&c&d))
        end

        @testset "Expand * to main effects + interactions" begin
            @test parse!(:(a*b)) == parse!(:(a+b+a&b))
            @test sort_terms!(parse!(:(a*b*c))) == parse!(:(a+b+c+a&b+a&c+b&c+a&b&c))
            @test parse!(:(a + b*c)) == parse!(:(a + b + c + b&c))
            @test parse!(:(a*b + c)) == parse!(:(a + b + a&b + c))
        end

    end

    # @testset "@formula" begin

    #     @test (@formula(y ~ x)) == Formula(:y, :x)
    #     @test (@formula(y ~ 1)) == Formula(:y, 1)
    #     @test (@formula(y ~ 1 + x)) == Formula(:y, :(1+x))

    #     @test (@formula( ~ x)) == Formula(nothing, :x)
    #     @test (@formula( ~ 1)) == Formula(nothing, 1)
    #     @test (@formula( ~ 1 + x)) == Formula(nothing, :(1+x))

    # end

end


text:
# This file should probably be moved to the general directory and used by
# classical MC as well

################################################################################
### Utility updates
################################################################################

"""
    AbstractUpdate

A update should be a struct inhereting from AbstractUpdate or one of its
abstract children `AbstractLocalUpdate`, `AbstractGlobalUpdate` or
`AbstractParallelUpdate`.

```
struct MyGlobalUpdate <: MonteCarlo.AbstractGlobalUpdate
    ...
end
```

It should implement the methods

```
function MonteCarlo.update(u::MyGlobalUpdate, mc, model)
    mc.temp_conf = ...
    return global_update(mc, model, mc.temp_conf)
end
MonteCarlo.name(::MyGlobalUpdate) = "MyGlobalUpdate"
```

The latter is used for printing statistics of the scheduler. The former defines
what the implemented update does. Usually this means creating a new
configuration and handing it to `global_update` to do a standard Metropolis
update. Note that you can and should use `mc.temp_conf` as temporary storage.

Behind the scenes, the scheduler will wrap `MyGlobalUpdate` in
`AcceptanceStatistics` which collects the number of requested and accepted
updates. It is expected that you return `0` if the update is denied or `1` if it
is accepted (as does the `global_update` returned above).
"""
abstract type AbstractUpdate end
abstract type AbstractLocalUpdate <: AbstractUpdate end

init!(mc, update::AbstractUpdate) = nothing
is_full_sweep(update::AbstractUpdate) = true

"""
    NoUpdate([mc, model])

An update that does nothing. Used internally to keep the adaptive scheduler
running if all (other) adaptive updates are discarded.
"""
struct NoUpdate <: AbstractUpdate end
NoUpdate(mc, model) = NoUpdate()
function update(u::NoUpdate, args...)
    # we count this as "denied" global update
    return 0
end
name(::NoUpdate) = "NoUpdate"
is_full_sweep(update::NoUpdate) = false

"""
    Adaptive()

A placeholder for adaptive updates in the `AdaptiveScheduler`.
"""
struct Adaptive <: AbstractUpdate end
name(::Adaptive) = "Adaptive"

# Should this inherit from AbstractUpdate?
# This is required for AdaptiveScheduler, used by both
mutable struct AcceptanceStatistics{Update <: AbstractUpdate}
    accepted::Float64
    total::Int
    update::Update
end
AcceptanceStatistics(update) = AcceptanceStatistics(0.0, 0, update)
AcceptanceStatistics(wrapped::AcceptanceStatistics) = wrapped
AcceptanceStatistics(proxy::Adaptive) = proxy
name(w::AcceptanceStatistics) = name(w.update)
function update(w::AcceptanceStatistics, mc, m)
    accepted = update(w.update, mc, m)
    w.total += 1
    w.accepted += accepted
    return accepted
end
function Base.show(io::IO, u::AcceptanceStatistics)
    @printf(io,
        "%s with %i/%i = %0.1f%c accepted",
        name(u), u.accepted, u.total, u.accepted/max(1, u.accepted), '%'
    )
end
init!(mc, update::AcceptanceStatistics) = init!(mc, update.update)
is_full_sweep(update::AcceptanceStatistics) = is_full_sweep(update.update)

################################################################################
### Generic
################################################################################

updates(s::AbstractUpdateScheduler) = s.sequence

function init_scheduler!(mc, scheduler::AbstractUpdateScheduler)
    for update in updates(scheduler)
        init!(mc, update)
    end
    nothing
end

################################################################################
### SimpleScheduler
################################################################################

"""
    SimpleScheduler(updates...)

Schedules updates in the order specified through `updates...`. Note that local
updates (`LocalSweep([N=1])`) always needs to be part of those updates, as they
define what a sweep is.

Example:
```
scheduler = SimpleScheduler(
    LocalSweep(10),
    GlobalFlip(),
    GlobalShuffle(),
    LocalSweep(10),
    ReplicaExchange(3)
)
```
This defines a scheduler that performs 10 local sweeps, then attempts a global
flip of the configuration, then a global shuffle, another 10 local sweeps and
a replica exchange with worker 3. After the replica exchange the sequence
repeats from the start.
"""
mutable struct SimpleScheduler{ST} <: AbstractUpdateScheduler
    sequence::ST
    idx::Int64

    # For loading
    function SimpleScheduler(sequence::ST, idx::Int) where {ST}
        new{ST}(sequence, idx)
    end

    function SimpleScheduler(updates...)
        # Without local sweeps we don't have a sweep increment
        updates = map(AcceptanceStatistics, Tuple(vcat(updates...)))
        if !any(u -> u isa AcceptanceStatistics{<: AbstractLocalUpdate}, updates)
            error("The scheduler requires local updates, but none were passed.")
        end
        any(isa.(updates, Adaptive)) && error("SimpleScheduler cannot process Adaptive. Use AdaptiveScheduler.")
        obj = new{typeof(updates)}()
        obj.sequence = updates
        obj.idx = 0
        obj
    end
end

function update(s::SimpleScheduler, mc::DQMC, model)
    while true
        s.idx = mod1(s.idx + 1, length(s.sequence))
        update(s.sequence[s.idx], mc, model)
        is_full_sweep(s.sequence[s.idx]) && break
    end
    mc.last_sweep += 1
    return nothing
end

function show_statistics(io::IO, s::SimpleScheduler, prefix="")
    println(io, prefix, "Update statistics (since start):")

    cum_accepted  = mapreduce(u -> u isa AcceptanceStatistics ? u.accepted : 0, +, s.sequence)
    cum_total  = mapreduce(u -> u isa AcceptanceStatistics ? u.total : 0, +, s.sequence)

    accumulated = flatten_sequence_statistics(s.sequence)
    show_accumulated_sequence(io, accumulated, prefix * "\t", cum_total)

    p = @sprintf("%2.1f", 100cum_accepted/max(1, cum_total))
    s = rpad("Total", 20) * " " * lpad(p, 5) * "% accepted" *
        "   (" * string(round(Int, cum_accepted)) * " / $cum_total)"

    println(io, prefix, '\t', '-' ^ length(s))
    println(io, "$prefix\t", s)

    nothing
end

function Base.show(io::IO, s::SimpleScheduler)
    sequence = combine_sequence_to_string(s.sequence)
    print(io, "SimpleScheduler(): $sequence -> (repeat)")
end

function save_scheduler(file::JLDFile, s::SimpleScheduler, entryname::String="/Scheduler")
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/tag", "SimpleScheduler")
    write(file, entryname * "/sequence", s.sequence)
    write(file, entryname * "/idx", s.idx)
    nothing
end
function _load(data, ::Val{:SimpleScheduler})
    SimpleScheduler(data["sequence"], data["idx"])
end

to_tag(::Type{<: SimpleScheduler}) = Val(:SimpleScheduler)

################################################################################
### AdaptiveScheduler
################################################################################

"""
    AdaptiveScheduler(sequence::Tuple, pool::Tuple; kwargs...])

Creates an adaptive scheduler which prioritizes updates that succeed more
frequently.

`sequence` defines the static sequence of global updates, with
`Adaptive()` as a proxy for adaptive updates. When polling the scheduler
`Adaptive()` will be replaced by one of the updates in `pool` based on their
success rate.

Example:
```
scheduler = AdaptiveScheduler(
    (Adaptive(), LocalSweep(10), ReplicaExchange(3)),
    (GlobalFlip(), GlobalShuffle())
)
```
This setup will perform one of `GlobalFlip()` and `GlobalShuffle()` followed by
10 local sweeps, a `ReplicaExchange(3)` update, repeating after that. The
replica exchange will always be a `ReplicaExchange()` while the ealier global
update will be picked based on their relative sampling rates which are sluggishly
derived from acceptance rates.

If those sampling rates fall below a certain threshhold they may be set to 0.
If all updates drop to 0 sampling rate, `NoUpdate()` is performed instead.
(`NoUpdate()` is always added and has a static probability of `1e-10`.)

### Keyword Arguments:
- `minimum_sampling_rate = 0.01`: This defines the threshold under which the
sampling rate is set to 0.
- `grace_period = 99`: This sets a minimum number of times an update needs to
be called before its sampling rate is adjusted.
- `adaptive_rate = 9.0`: Controls how fast the sampling rate is adjusted to the
acceptance rate. More means slower. This follows the formula
`(adaptive_rate * sampling_rate + accepted/total) / (adaptive_rate + 1)`
"""
mutable struct AdaptiveScheduler{PT, ST} <: AbstractUpdateScheduler
    minimum_sampling_rate::Float64
    grace_period::Int
    adaptive_rate::Float64

    adaptive_pool::PT
    sampling_rates::Vector{Float64}
    sequence::ST
    idx::Int

    # for loading
    function AdaptiveScheduler(msr, gp, ar, ap::PT, sr, s::ST, i) where {PT, ST}
        new{PT, ST}(msr, gp, ar, ap, sr, s, i)
    end

    function AdaptiveScheduler(sequence, pool;
            minimum_sampling_rate = 0.01, grace_period = 99, adaptive_rate = 9.0
        )
        # Without local updates we have no sweep increment
        sequence = map(AcceptanceStatistics, vcat(sequence...))
        if !any(u -> u isa AcceptanceStatistics{<: AbstractLocalUpdate}, sequence)
            error("The scheduler requires local updates, but none were passed.")
        end

        # sampling rate x after i steps with constant acceptance rate p:
        # x_i = (c/c+1)^i x_0 + a * \sum_{i} c^{i-1} / (c+1)^i
        # where c is the adaptive rate and x_0 is the initial sampling rate.
        # Minimum number of samples to discard:
        i_min = ceil(Int, log(
            adaptive_rate / (adaptive_rate+1),
            minimum_sampling_rate / 0.5
        ))
        @debug("Minimum number of samples for discard: $grace_period + $i_min")

        if !any(x -> x isa NoUpdate, pool)
            pool = tuple(vcat(pool...)..., NoUpdate())
        end

        # Wrap in AcceptanceStatistics
        adaptive_pool = map(AcceptanceStatistics, pool)

        adaptive_pool = Tuple(adaptive_pool)
        sequence = Tuple(sequence)

        obj = new{typeof(adaptive_pool), typeof(sequence)}()
        obj.minimum_sampling_rate = minimum_sampling_rate
        obj.grace_period = grace_period
        obj.adaptive_rate = adaptive_rate
        obj.adaptive_pool = adaptive_pool
        obj.sampling_rates = [u isa AcceptanceStatistics{NoUpdate} ? 1e-10 : 0.5 for u in adaptive_pool]
        obj.sequence = sequence
        obj.idx = 0

        obj
    end
end

function update(s::AdaptiveScheduler, mc::DQMC, model)
    while true
        s.idx = mod1(s.idx + 1, length(s.sequence))

        if s.sequence[s.idx] === Adaptive()
            # Find appropriate (i.e. with probability matching sampling rates)
            # update from adaptive pool
            total_weight = sum(s.sampling_rates)
            target = rand() * total_weight
            current_weight = 0.0
            idx = 1
            while idx < length(s.sampling_rates)
                current_weight += s.sampling_rates[idx]
                if target <= current_weight
                    break
                else
                    idx += 1
                end
            end

            # Apply the update and adjust sampling rate
            updater = s.adaptive_pool[idx]
            update(updater, mc, model)
            if !(updater isa AcceptanceStatistics{NoUpdate}) && updater.total > s.grace_period
                s.sampling_rates[idx] = (
                    s.adaptive_rate * s.sampling_rates[idx] +
                    updater.accepted / updater.total
                ) / (s.adaptive_rate + 1)

                # Hit miniomum threshold - this can no longer be accepted.
                if s.sampling_rates[idx] < s.minimum_sampling_rate
                    s.sampling_rates[idx] = 0.0
                end
            end

        else
            # Some sort of non-adaptive update, just perform it.
            updater = s.sequence[s.idx]
            update(updater, mc, model)
        end

        is_full_sweep(updater) && break
    end

    mc.last_sweep += 1

    return nothing
end

function show_statistics(io::IO, s::AdaptiveScheduler, prefix="")
    println(io, prefix, "Update statistics (since start):")

    cum_accepted  = mapreduce(u -> u isa AcceptanceStatistics ? u.accepted : 0, +, s.sequence)
    cum_accepted += mapreduce(u -> u isa AcceptanceStatistics ? u.accepted : 0, +, s.adaptive_pool)
    cum_total  = mapreduce(u -> u isa AcceptanceStatistics ? u.total : 0, +, s.sequence)
    cum_total += mapreduce(u -> u isa AcceptanceStatistics ? u.total : 0, +, s.adaptive_pool)

    accumulated = flatten_sequence_statistics(s.adaptive_pool)
    flatten_sequence_statistics(s.sequence, accumulated)
    show_accumulated_sequence(io, accumulated, prefix * "\t", cum_total)

    p = @sprintf("%2.1f", 100cum_accepted/max(1, cum_total))
    s = rpad("Total", 20) * " " * lpad(p, 5) * "% accepted" *
        "   (" * string(round(Int, cum_accepted)) * " / $cum_total)"

    println(io, prefix, '\t', '-' ^ length(s))
    println(io, "$prefix\t", s)

    nothing
end

function Base.show(io::IO, s::AdaptiveScheduler)
    sequence = combine_sequence_to_string(s.sequence)
    total_weight = sum(s.sampling_rates)
    pool = mapreduce((a, b) -> "$a, $b", s.sampling_rates, s.adaptive_pool) do sr, u
        @sprintf("%0.0f%s %s", 100sr / total_weight, "%", name(u))
    end
    println(io, "AdaptiveScheduler():")
    println(io, "\t$sequence -> (repeat)")
    print(io, "\twith Adaptive() = ($pool)")
end

function save_scheduler(file::JLDFile, s::AdaptiveScheduler, entryname::String="/Scheduler")
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/tag", "AdaptiveScheduler")
    write(file, entryname * "/sequence", s.sequence)
    write(file, entryname * "/sampling_rates", s.sampling_rates)
    write(file, entryname * "/pool", s.adaptive_pool)
    write(file, entryname * "/minimum_sampling_rate", s.minimum_sampling_rate)
    write(file, entryname * "/grace_period", s.grace_period)
    write(file, entryname * "/adaptive_rate", s.adaptive_rate)
    write(file, entryname * "/idx", s.idx)
    nothing
end
function _load(data, ::Val{:AdaptiveScheduler})
    s = AdaptiveScheduler(data["sequence"], data["pool"])
    s.sampling_rates = data["sampling_rates"]
    s.minimum_sampling_rate = data["minimum_sampling_rate"]
    s.grace_period = data["grace_period"]
    s.adaptive_rate = data["adaptive_rate"]
    s.idx = data["idx"]
    s
end

to_tag(::Type{<: AdaptiveScheduler}) = Val(:AdaptiveScheduler)

updates(s::AdaptiveScheduler) = (s.sequence..., s.adaptive_pool...)

################################################################################
### Printing utillities
################################################################################

function flatten_sequence_statistics(sequence, accumulated::Dict = Dict{String, Tuple{Float64, Int}}())
    for update in sequence
        if update isa Adaptive || !(update isa AcceptanceStatistics)
            continue
        end

        key = name(update)
        acc = update.accepted
        total = update.total

        if haskey(accumulated, key)
            accumulated[key] = accumulated[key] .+ (acc, total)
        else
            accumulated[key] = (acc, total)
        end
    end

    accumulated
end

function show_accumulated_sequence(io::IO, accumulated, prefix = "", max_total=0)
    N = length(string(max_total))
    sorted = sort!(collect(accumulated), by = x -> x[2][1] / max(1, x[2][2]))
    for (key, (acc, total)) in sorted
        total == 0 && continue
        p = @sprintf("%2.1f", 100acc / max(1, total))
        println(
            io, prefix, rpad(key, 20), " ", lpad(p, 5), "% accepted",  "   (",
            lpad(string(round(Int, acc)), N), " / ", lpad(string(total), N), ")"
        )
    end

    nothing
end

function show_sequence(io::IO, sequence, prefix = "", max_total=0)
    N = length(string(max_total))
    for update in sequence
        if update isa AcceptanceStatistics{NoUpdate} ||
            update isa Adaptive || !(update isa AcceptanceStatistics)
            continue
        end

        p = @sprintf("%2.1f", 100update.accepted/max(1, update.total))
        println(
            io, prefix, rpad(name(update), 20), " ", lpad(p, 5), "% accepted", "(",
            lpad(string(round(Int, update.accepted)), N), " / ",
            lpad(string(update.total), N), ")"
        )
    end

    nothing
end

function combine_sequence_to_string(sequence)
    last_update = first(sequence)
    counter = 1
    output = ""
    for update in sequence[2:end]
        if update == last_update
            counter += 1
        else
            _name = name(last_update)
            if counter == 1
                output *= " -> $_name"
            else
                output *= " -> $(_name)($counter)"
            end
            counter = 1
            last_update = update
        end
    end

    _name = name(last(sequence))
    if counter == 1
        output *= " -> $_name"
    else
        output *= " -> $(_name)($counter)"
    end

    output
end

# Literally just for tests
for T in (SimpleScheduler, AdaptiveScheduler, AcceptanceStatistics)
    @eval begin
        function Base.:(==)(a::$T, b::$T)
            for field in fieldnames($T)
                getfield(a, field) == getfield(b, field) || return false
            end
            return true
        end
    end
end


text:
# Invertible CNN layer from Dinh et al. (2017)/Kingma and Dhariwal (2018)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test, Random

# Random seed
Random.seed!(11)

###################################################################################################
# Test invertibility

# Input
nx = 24
ny = 24
k = 4
n_in = 2
n_hidden = 4
batchsize = 1

# Input images
Xa = randn(Float32, nx, ny, Int(k/2), batchsize)
Xb = randn(Float32, nx, ny, Int(k/2), batchsize)
Xa0 = randn(Float32, nx, ny, Int(k/2), batchsize)
Xb0 = randn(Float32, nx, ny, Int(k/2), batchsize)
dXa = Xa - Xa0
dXb = Xb - Xb0

# 1x1 convolution and residual blocks
RB = ResidualBlock(nx, ny, n_in, n_hidden, batchsize; fan=true)
L = CouplingLayerBasic(RB; logdet=true)

###################################################################################################
# Invertibility tests

# Layer w/ logdet
Ya, Yb, logdet = L.forward(Xa, Xb)
Xa_, Xb_ = L.inverse(Ya, Yb)
@test isapprox(norm(Xa - Xa_)/norm(Xa), 0f0; atol=1e-2)
@test isapprox(norm(Xb - Xb_)/norm(Xb), 0f0; atol=1e-2)

Ya, Yb, logdet = L.forward(Xa, Xb)
Xa_, Xb_ = L.backward(Ya.*0f0, Yb.*0f0, Ya, Yb)[3:4]
@test isapprox(norm(Xa - Xa_)/norm(Xa), 0f0; atol=1e-2)
@test isapprox(norm(Xb - Xb_)/norm(Xb), 0f0; atol=1e-2)

Ya, Yb = L.inverse(Xa, Xb)
Xa_, Xb_, logdet = L.forward(Ya, Yb)
@test isapprox(norm(Xa - Xa_)/norm(Xa), 0f0; atol=1e-2)
@test isapprox(norm(Xb - Xb_)/norm(Xb), 0f0; atol=1e-2)

# Layer w/o logdet
L.logdet = false
Ya, Yb = L.forward(Xa, Xb)
Xa_, Xb_ = L.inverse(Ya, Yb)
@test isapprox(norm(Xa - Xa_)/norm(Xa), 0f0; atol=1e-2)
@test isapprox(norm(Xb - Xb_)/norm(Xb), 0f0; atol=1e-2)

Ya, Yb = L.forward(Xa, Xb)
Xa_, Xb_ = L.backward(Ya.*0f0, Yb.*0f0, Ya, Yb)[3:4]
@test isapprox(norm(Xa - Xa_)/norm(Xa), 0f0; atol=1e-2)
@test isapprox(norm(Xb - Xb_)/norm(Xb), 0f0; atol=1e-2)

Ya, Yb = L.inverse(Xa, Xb)
Xa_, Xb_ = L.forward(Ya, Yb)
@test isapprox(norm(Xa - Xa_)/norm(Xa), 0f0; atol=1e-2)
@test isapprox(norm(Xb - Xb_)/norm(Xb), 0f0; atol=1e-2)

# Reverse layer w/ logdet
L.logdet = true
L_rev = reverse(L)
Ya, Yb, logdet = L_rev.forward(Xa, Xb)
Xa_, Xb_ = L_rev.inverse(Ya, Yb)
@test isapprox(norm(Xa - Xa_)/norm(Xa), 0f0; atol=1e-2)
@test isapprox(norm(Xb - Xb_)/norm(Xb), 0f0; atol=1e-2)

Ya, Yb, logdet = L_rev.forward(Xa, Xb)
Xa_, Xb_ = L_rev.backward(Ya.*0f0, Yb.*0f0, Ya, Yb)[3:4]
@test isapprox(norm(Xa - Xa_)/norm(Xa), 0f0; atol=1e-2)
@test isapprox(norm(Xb - Xb_)/norm(Xb), 0f0; atol=1e-2)

Ya, Yb = L_rev.inverse(Xa, Xb)
Xa_, Xb_, logdet = L_rev.forward(Ya, Yb)
@test isapprox(norm(Xa - Xa_)/norm(Xa), 0f0; atol=1e-2)
@test isapprox(norm(Xb - Xb_)/norm(Xb), 0f0; atol=1e-2)

###################################################################################################
# Gradient tests

# Loss Function
function loss(L, Xa, Xb, Ya, Yb)
    Ya_, Yb_, logdet = L.forward(Xa, Xb)
    f = mse(tensor_cat(Ya_, Yb_), tensor_cat(Ya, Yb)) - logdet
    ΔY = ∇mse(tensor_cat(Ya_, Yb_), tensor_cat(Ya, Yb))
    ΔYa, ΔYb = tensor_split(ΔY)
    ΔXa, ΔXb = L.backward(ΔYa, ΔYb, Ya_, Yb_)[1:2]

    # Pass back gradients w.r.t. input X and from the residual block and 1x1 conv. layer
    return f, ΔXa, ΔXb, L.RB.W1.grad, L.RB.W2.grad, L.RB.W3.grad
end

# Invertible layers
RB0 = ResidualBlock(nx, ny, n_in, n_hidden, batchsize; fan=true)
L = CouplingLayerBasic(RB; logdet=true)
L01 = CouplingLayerBasic(RB; logdet=true)
L02 = CouplingLayerBasic(RB0; logdet=true)

# Gradient test w.r.t. input X0
Ya, Yb = L.forward(Xa, Xb)[1:2]
f0, ΔXa, ΔXb  = loss(L, Xa0, Xb0, Ya, Yb)[1:3]
h = 0.1f0
maxiter = 6
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test coupling layer\n")
for j=1:maxiter
    f = loss(L, Xa0 + h*dXa, Xb0 + h*dXb, Ya, Yb)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dXa, ΔXa) - h*dot(dXb, ΔXb))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test w.r.t. weights of residual block
Ya, Yb = L.forward(Xa, Xb)[1:2]
Lini = deepcopy(L02)
dW1 = L.RB.W1.data - L02.RB.W1.data
dW2 = L.RB.W2.data - L02.RB.W2.data
dW3 = L.RB.W3.data - L02.RB.W3.data

f0, ΔXa, ΔXb, ΔW1, ΔW2, ΔW3 = loss(L02, Xa, Xb, Ya, Yb)
h = 0.1f0
maxiter = 4
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

print("\nGradient test coupling layer\n")
for j=1:maxiter
    L02.RB.W1.data = Lini.RB.W1.data + h*dW1
    L02.RB.W2.data = Lini.RB.W2.data + h*dW2
    L02.RB.W3.data = Lini.RB.W3.data + h*dW3
    f = loss(L02, Xa, Xb, Ya, Yb)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(dW1, ΔW1) - h*dot(dW2, ΔW2) - h*dot(dW3, ΔW3))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)

###################################################################################################
# Gradient tests

# Loss Function
function loss(L, Xa, Xb, Ya, Yb)
    Ya_, Yb_, logdet = L.forward(Xa, Xb)
    f = mse(tensor_cat(Ya_, Yb_), tensor_cat(Ya, Yb)) - logdet
    ΔY = ∇mse(tensor_cat(Ya_, Yb_), tensor_cat(Ya, Yb))
    ΔYa, ΔYb = tensor_split(ΔY)
    ΔXa, ΔXb = L.backward(ΔYa, ΔYb, Ya_, Yb_)[1:2]

    # Pass back gradients w.r.t. input X and from the residual block and 1x1 conv. layer
    return f, ΔXa, ΔXb, L.RB.W1.grad, L.RB.W2.grad, L.RB.W3.grad
end

# Invertible layers
RB0 = ResidualBlock(nx, ny, n_in, n_hidden, batchsize; fan=true)
L = reverse(CouplingLayerBasic(RB; logdet=true))
L01 = reverse(CouplingLayerBasic(RB; logdet=true))
L02 = reverse(CouplingLayerBasic(RB0; logdet=true))

# Gradient test w.r.t. input X0
Ya, Yb = L.forward(Xa, Xb)[1:2]
f0, ΔXa, ΔXb  = loss(L, Xa0, Xb0, Ya, Yb)[1:3]
h = 0.1f0
maxiter = 6
err5 = zeros(Float32, maxiter)
err6 = zeros(Float32, maxiter)

print("\nGradient test coupling layer\n")
for j=1:maxiter
    f = loss(L, Xa0 + h*dXa, Xb0 + h*dXb, Ya, Yb)[1]
    err5[j] = abs(f - f0)
    err6[j] = abs(f - f0 - h*dot(dXa, ΔXa) - h*dot(dXb, ΔXb))
    print(err5[j], "; ", err6[j], "\n")
    global h = h/2f0
end

@test isapprox(err5[end] / (err5[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err6[end] / (err6[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test w.r.t. weights of residual block
Ya, Yb = L.forward(Xa, Xb)[1:2]
Lini = deepcopy(L02)
dW1 = L.RB.W1.data - L02.RB.W1.data
dW2 = L.RB.W2.data - L02.RB.W2.data
dW3 = L.RB.W3.data - L02.RB.W3.data

f0, ΔXa, ΔXb, ΔW1, ΔW2, ΔW3 = loss(L02, Xa, Xb, Ya, Yb)
h = 0.1f0
maxiter = 4
err7 = zeros(Float32, maxiter)
err8 = zeros(Float32, maxiter)

print("\nGradient test coupling layer\n")
for j=1:maxiter
    L02.RB.W1.data = Lini.RB.W1.data + h*dW1
    L02.RB.W2.data = Lini.RB.W2.data + h*dW2
    L02.RB.W3.data = Lini.RB.W3.data + h*dW3
    f = loss(L02, Xa, Xb, Ya, Yb)[1]
    err7[j] = abs(f - f0)
    err8[j] = abs(f - f0 - h*dot(dW1, ΔW1) - h*dot(dW2, ΔW2) - h*dot(dW3, ΔW3))
    print(err7[j], "; ", err8[j], "\n")
    global h = h/2f0
end

@test isapprox(err7[end] / (err7[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err8[end] / (err8[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
module StiffKinetic
    using Revise

    include("problems.jl")
    include("space_derivatives.jl")
    include("time_schemes.jl")

    # export save_mesh
    export RelaxPb
    export posUpwind!, posWENO3!, posWENO5!
    export negUpwind!, negWENO3!, negWENO5!
    export posUpwind, posWENO3, posWENO5
    export negUpwind, negWENO3, negWENO5
    export fd_diff2_ord2!, fd_diff2_ord4!, fd_diff2_ord6!

    export swap_pbs!
    export imex_bdf1!, imex_bdf2!, imex_bdf3!, imex_bdf4!

end # module

max_stars_repo_name: tremelow/StiffKinetic

max_issues_repo_name: tremelow/StiffKinetic

max_forks_repo_name: tremelow/StiffKinetic


text:
using TuringModels

Turing.setadbackend(:reverse_diff)
#nbTuring.turnprogress(false);

d = CSV.read(joinpath(@__DIR__, "..", "..", "data", "Kline.csv"), delim=';');
size(d) # Should be 10x5

# New col log_pop, set log() for population data
d[!, :log_pop] = map((x) -> log(x), d[:, :population]);
d[!, :society] = 1:10;

# Turing model

@model m12_6(total_tools, log_pop, society) = begin

    # Total num of y
    N = length(total_tools)

    # priors
    α ~ Normal(0, 10)
    βp ~ Normal(0, 1)

    # Separate σ priors for each society
    σ_society ~ Truncated(Cauchy(0, 1), 0, Inf)

    # Number of unique societies in the data set
    N_society = length(unique(society)) #10

    # Vector of societies (1,..,10) which we'll set priors on
    α_society = Vector{Real}(undef, N_society)

    # For each society [1,..,10] set a prior N(0, σ_society)
    α_society ~ [Normal(0, σ_society)]

    for i ∈ 1:N
        λ = exp(α + α_society[society[i]] + βp*log_pop[i])
        total_tools[i] ~ Poisson(λ)
    end
end

# Sample

chns = sample(m12_6(d[:, :total_tools], d[:, :log_pop],
    d[:, :society]), NUTS(0.95), 1000);

# Results rethinking
m12_6rethinking = "
              Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
a              1.11   0.75      -0.05       2.24  1256    1
bp             0.26   0.08       0.13       0.38  1276    1
a_society[1]  -0.20   0.24      -0.57       0.16  2389    1
a_society[2]   0.04   0.21      -0.29       0.38  2220    1
a_society[3]  -0.05   0.19      -0.36       0.25  3018    1
a_society[4]   0.32   0.18       0.01       0.60  2153    1
a_society[5]   0.04   0.18      -0.22       0.33  3196    1
a_society[6]  -0.32   0.21      -0.62       0.02  2574    1
a_society[7]   0.14   0.17      -0.13       0.40  2751    1
a_society[8]  -0.18   0.19      -0.46       0.12  2952    1
a_society[9]   0.27   0.17      -0.02       0.52  2540    1
a_society[10] -0.10   0.30      -0.52       0.37  1433    1
sigma_society  0.31   0.13       0.11       0.47  1345    1
";

# Describe the posterior samples

describe(chns)

# End of `12/m12.6t.jl`


text:
using NonlinearSolve

""" Set the initial condition: rootfinding for finding the steady state solution """
function solve_ss(
    model,
    u0,
    p_ode,
    state_dim,
    batchsize;
    u0_guess = zeros(state_dim, batchsize),
    abstol = 1e-10,
)
    # backward pass with matrix valued steady state does not seem to work
    # ssprob = SteadyStateProblem((x,p,t) -> model(vcat(reshape(x, state_dim, batchsize), u0), p), u0_guess, p_ode)
    # x0 = reshape(Array(solve(ssprob, SSRootfind(nlsolve=(f,u0,abstol) -> (res=NLsolve.nlsolve(f,u0,ftol=abstol, method=method);res.zero)); abstol=abstol)), state_dim, batchsize)
    # return x0
    ssprob = NonlinearProblem{false}(
        (x, p) -> vec(model(vcat(reshape(x, state_dim, batchsize), u0), p)),
        vec(u0_guess),
        p_ode,
    )
    return reshape(
        Array(solve(ssprob, NewtonRaphson(), tol = abstol)),
        state_dim,
        batchsize,
    )
end

""" Set the initial condition: rootfinding for finding the steady state solution, but do not backprop through this """
function solve_ss_dropgrad(
    model,
    u0,
    p_ode,
    state_dim,
    batchsize;
    u0_guess = zeros(state_dim, batchsize),
    method = :newton,
    abstol = 1e-10,
)
    x0 = u0_guess
    Zygote.ignore() do
        # ssprob = SteadyStateProblem((x,p,t) -> model(vcat(reshape(x, state_dim, batchsize), u0), p), u0_guess, p_ode)
        # x0 = reshape(Array(solve(ssprob, SSRootfind(nlsolve=(f,u0,abstol) -> (res=NLsolve.nlsolve(f,u0,ftol=abstol, method=method);res.zero)); abstol=abstol)), state_dim, batchsize)
        ssprob = NonlinearProblem{false}(
            (x, p) -> vec(model(vcat(reshape(x, state_dim, batchsize), u0), p)),
            vec(u0_guess),
            p_ode,
        )
        x0 = reshape(
            Array(solve(ssprob, NewtonRaphson(), tol = abstol)),
            state_dim,
            batchsize,
        )
    end
    return x0
end

""" Set the initial condition: fixed value, default Float64 zeros """
function fixed_ic(
    model,
    u0,
    p_ode,
    state_dim,
    batchsize;
    u0_guess = zeros(state_dim, batchsize),
    method = :newton,
    abstol = 1e-8,
)
    return u0_guess
end

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
# -----------------------------------------------------------------------
# Direct Radiative Forcing From Methane (based on FUND model equations).
# -----------------------------------------------------------------------

# Create function to account for overlap in CH₄ and N₂O absoprtion bands.
function interact(M, N)
    d = 1.0 + (M * N)^0.75 * 2.01E-5 + (M * N)^1.52 * M * 5.31E-15
    return 0.47 * log(d)
end

@defcomp rf_ch4_direct_fund begin

    CH₄_0         = Parameter()             # Initial (pre-industrial) methane concentration (ppb).
    N₂O_0         = Parameter()             # Initial (pre-industrial) N₂O concentration (ppb).
    scale_CH₄     = Parameter()             # Scaling factor for uncertainty in methane radiative forcing.
    CH₄           = Parameter(index=[time]) # Atmospheric methane concetration for current period (ppb).

    rf_ch4_direct = Variable(index=[time])  # Direct forcing from atmospheric methane concentrations (Wm⁻²).

    function run_timestep(p, v, d, t)
        if is_first(t)
            # Set initial radiative forcing to 0.
            v.rf_ch4_direct[t] = 0.0
        else
            # Create temporary variable accounting for overlap between pre-industrial CH₄ and N₂O (for convenience).
            ch4n2o = interact(p.CH₄_0, p.N₂O_0)

            # Calculate direct CH₄ radiative forcing.
            v.rf_ch4_direct[t] = (0.036 * (sqrt(p.CH₄[t]) - sqrt(p.CH₄_0)) - interact(p.CH₄[t], p.N₂O_0) + ch4n2o) * p.scale_CH₄
        end
    end
end


text:
using HorizonSideRobots

function try_move!(r::Robot,side::HorizonSide)::Bool
    ort_side = left(side)
    n=0
    while isborder(r,side)
        if !isborder(r,ort_side)
            move!(r,ort_side)
            n += 1
        else
            break
        end
    end
    if isborder(r,side)
        movements!(r,inverse(ort_side),n)
        return false
    end
    move!(r,side)
    while isborder(r,inverse(ort_side))
        move!(r,side)
    end
    movements!(r, inverse(ort_side), n)
    return true
end

function movements!(r::Robot, side)
    n=0
    while !isborder(r,side)
        move!(r,side)
        n += 1
    end
end

function movements!(r::Robot, side, num_steps::Integer)
    for _ in 1:num_steps
        move!(r,side)
    end
end

inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2, 4))
left(side::HorizonSide) = HorizonSide(mod(Int(side)+1, 4))

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

abstract type AbstractRobot end

import HorizonSideRobots: move!, isborder, putmarker!, ismarker, temperature

move!(r::AbstractRobot, side) = move!(get(r), side)
isborder(r::AbstractRobot, side) = isborder(get(r), side)
putmarker!(r::AbstractRobot) = putmarker!(get(r))
ismarker(r::AbstractRobot) = ismarker(get(r))
temperature(r::AbstractRobot) = temperature(get(r))

# предполагается, что функция get() будет определена для каждого производного конкретного типа

# movements!(r::AbstractRobot, side) = movements!(get(r), side)
# movements!(r::AbstractRobot, side, num_steps::Integer) = movements!(get(r), side, num_steps)

function movements!(r::AbstractRobot, side)
    n=0
    while !isborder(r,side)
        move!(r,side)
        n += 1
    end
end

function movements!(r::AbstractRobot, side, num_steps::Integer)
    for _ in 1:num_steps
        move!(r,side)
    end
end
#-----------------------------------------------------

abstract type AbstractBorderRobot <: AbstractRobot end

function movements!(r::AbstractBorderRobot, side::HorizonSide)
    n=0
    while try_move!(r, side) # предполагается, что try_move! будет переопределяться для всех производных типов
        n+=1
    end
    return n
end

function movements!(r::AbstractBorderRobot, side::HorizonSide, n::Integer)
    for _ in 1:n
        try_move!(r, side) # как вариант - предполагается, что try_move! будет переопределяться для всех производных типов
    end
    return n
end

function try_move!(r::AbstractBorderRobot, side::HorizonSide)
    result = try_move!(get(r), side)
    if result == true
        action(r) # предполагается, что в каждом производном конткретном типе функция action определена
    end
    return result
end

#-------------------------------------------------------
struct PutmarkerRobot <: AbstractRobot
    robot::Robot
end

get(r::PutmarkerRobot) = r.robot

function move!(r::PutmarkerRobot, side)
    move!(r.robot, side)
    putmarker!(r)
end

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

mutable struct CountmarkersRobot <: AbstractRobot
    robot::Robot
    count::Int
end

get(r::CountmarkersRobot) = r.robot

function move!(r::CountmarkersRobot, side)
    move!(r.robot, side) # если бы тут было так: move!(r, side), то это был бы рекурсивный вызов, котрорый здесь не уместен !!!
    # можно было бы ещё и так: move!(get(r), side), разницы в этом никакой нет
    if ismarker(r)
        r.count += 1
    end
end
#----------------------------------

mutable struct CountmarkersBorderRobot <: AbstractBorderRobot
    robot::Robot
    count::Int
end

get(r::CountmarkersBorderRobot) = r.robot

function action(r::CountmarkersBorderRobot)
    if ismarker(r)
        r.count += 1
    end
end

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

struct PutmarkersBorderRobot <: AbstractBorderRobot
    robot::Robot
end

get(r::PutmarkersBorderRobot) = r.robot

action(r::PutmarkersBorderRobot) = putmarker!(r)

# ===================== 3 часть занятия =============================

mutable struct Coord
    x::Int
    y::Int
end
Coord()=Coord(0,0)

function move!(coord::Coord, side::HorizonSide)
    if side==Nord
        coord.y += 1
    elseif side==Sud
        coord.y -= 1
    elseif side==Ost
        coord.x += 1
    else #if side==West
        coord.x -= 1
    end
end

get_coord(coord::Coord) = (coord.x, coord.y)

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

mutable struct CoordRobot <: AbstractRobot
    robot::Robot
    coord::Coord
end

CoordRobot(r)=CoordRobot(r, Coord())

get(robot::CoordRobot) = robot.robot

function move!(robot::CoordRobot, side::HorizonSide)
    move!(get(robot), side)
    move!(robot.coord, side)
end

get_coord(r::CoordRobot)=get_coord(r.coord)
#--------------------------------------------------

# Проектирование праметрического абстрактного типа

abstract type AbstractBorderRobot2{TypeRobot} <: AbstractRobot end
# тип AbstractBorderRobot - уже есть, он не параметрический

function movements!(r::AbstractBorderRobot2, side::HorizonSide)
    n=0
    while try_move!(r, side) # предполагается, что try_move! будет переопределяться для всех производных типов
        n+=1
    end
    return n
end

function movements!(r::AbstractBorderRobot2, side::HorizonSide, n::Integer)
    for _ in 1:n
        try_move!(r, side) # как вариант - предполагается, что try_move! будет переопределяться для всех производных типов
    end
    return n
end

function try_move!(r::AbstractBorderRobot2, side::HorizonSide)
    result = try_move!(get(r), side)
    if result == true
        action(r) # предполагается, что в каждом производном конткретном типе функция action определена
    end
    return result
end

# Проектирование конкретного типа на базе абстрактного параметрического
struct CoordmarkersBorderRobot <: AbstractBorderRobot2{CoordRobot} # CoordRobot - ранее определенный тип
    robot::CoordRobot
    markers_coord::Vector{NTuple{2,Int}}
end

get(robot::CoordmarkersBorderRobot) = get(robot.robot)

function action(robot::CoordmarkersBorderRobot)
    if ismarker(robot)
        push!(robot.markers_coord, get_coord(robot))
    end
end

get_coord(r::CoordmarkersBorderRobot) = get_coord(r.robot)

#-------------------------------Дополнительно требуется определить:
get(r::Robot)=r


text:
mutable struct ConsoleManager <: GtkNotebook

    handle::Ptr{Gtk.GObject}
    main_window
    server::TCPServer
    port::UInt16
    top_module::Module
    watch_stdout_task::Task
    stdout
    stderr

    function ConsoleManager(main_window, top_module=GtkREPL)

        ntb = GtkNotebook()
        set_gtk_property!(ntb, :vexpand, true)
        port, server = RemoteGtkREPL.start_server()

        n = new(ntb.handle, main_window, server, port, top_module)
        Gtk.gobject_move_ref(n, ntb)
    end
end

function init!(console_mng::ConsoleManager)
    add_plus_button(console_mng)
    # add a task in Gtk's main main that writes stdout buffer in the console
    signal_connect(console_mng_button_press_cb,console_mng, "button-press-event",
    Cint, (Ptr{Gtk.GdkEvent},),false,console_mng.main_window)
    signal_connect(console_mng_switch_page_cb, console_mng, "switch-page", Nothing, (Ptr{Gtk.GtkWidget},Int32), false)
end

function init_stdout!(console_mng::ConsoleManager,watch_stdout_task,stdout,stderr)
    console_mng.watch_stdout_task = watch_stdout_task
    console_mng.stdout = stdout
    console_mng.stderr = stderr
end

function add_remote_console(main_window, mod=GtkREPL)
    port = console_manager(main_window).port
    id = length(console_manager(main_window)) + 1
    p = joinpath(HOMEDIR,"remote_console_startup.jl")
    #julia_path = Config.julia_path
    julia_path = ENV["_"] # ^_^
    s = "tell application \"Terminal\" to do script \"$(julia_path) -i --color=no \\\"$p\\\" $port $id $mod\""
    run(`osascript -e $s`)
end

function add_remote_console_cb(id, port)
    @info "GtkREPL: Starting console for port $port with id $id"

    c = try
        worker = connect(port)

        c = Console(id, main_window, worker)
        c.worker_port = port
        init!(c)

        #for some reason I need to warm-up things here, otherwise it bugs later on.
        @assert remotecall_fetch(identity, GtkREPL.worker(c),1) == 1

        showall(c.main_window)
        c
    catch err
        @warn err
    end

    RemoteGtkREPL.remotecall_fetch(println, worker(c),"Worker connected")
    "done"
end

@guarded (INTERRUPT) function console_mng_button_press_cb(widgetptr::Ptr, eventptr::Ptr, user_data)

    ntbook = convert(GtkNotebook, widgetptr)
    event = convert(Gtk.GdkEvent, eventptr)
    main_window = user_data #TODO ntbook should be a console manager with a MainWindow field?

    if rightclick(event)
        menu = buildmenu([
            MenuItem("Close Console", remove_console_cb),
            MenuItem("Add Console", add_console_cb)
            ],
            (ntbook, get_current_console(ntbook), main_window)
        )
        popup(menu,event)
        return INTERRUPT
    end

    return PROPAGATE
end

function console_mng_switch_page_cb(widgetptr::Ptr, pageptr::Ptr, pagenum::Int32, user_data)

    cm = convert(GtkNotebook, widgetptr)
    c  = convert(Gtk.GtkWidget, pageptr)

    on_console_mng_switch_page(cm, c)

    nothing
end

function current_console(m::ConsoleManager)
    c = m[index(m)]
    c isa Console && return c
    for c in m
        c isa Console && return c
    end
    c
end

function stop_console_redirect(main_window)

    t = main_window.console_manager.watch_stdout_task
    out = main_window.console_manager.stdout
    #err = main_window.console_manager.stderr

# The task end itself when is_running == false
#    try
#        Base.throwto(t, InterruptException())
#    end

    sleep(0.1)
    redirect_stdout(out)
    #redirect_stderr(err)
end


text:
global gh_auth = nothing

function github_auth(;token::AbstractString="")
   if !isempty(token)
      global gh_auth = GitHub.authenticate(token)
   elseif haskey(ENV,"GITHUB_TOKEN")
      global gh_auth = GitHub.authenticate(ENV["GITHUB_TOKEN"])
   else
      @info "Using anonymous github auth"
      global gh_auth = GitHub.AnonymousAuth()
   end
end

function github_repo(pkg::AbstractString; fork=nothing)
   isnothing(gh_auth) && github_auth()
   return GitHub.repo(pkg_url(pkg; full=false, fork=fork); auth=gh_auth)
end

function github_repo_exists(name::AbstractString, org::AbstractString)
   try
      github_repo(name; fork=org)
      return true
   catch
      return false
   end
end

parse_meta(file::AbstractString) = TOML.parsefile(file)

# we try to find a matching branch in the main repo for that pkg
# or in a given fork, or via a full 'url#branch'
function find_branch(pkg::AbstractString, branch::AbstractString; fork=nothing)
   isnothing(gh_auth) && github_auth()
   if startswith(branch, "https://")
      urlmatch = match(r"https://github\.com/([-\w]+)/\w+\.jl#(.*)", branch)
      if isnothing(urlmatch)
         @error "could not parse org and branch from $branch"
      end
      if urlmatch[1] == pkg_org(pkg)
         return find_branch(pkg, urlmatch[2])
      else
         return find_branch(pkg, urlmatch[2]; fork=urlmatch[1])
      end
   end
   @info "locating branch '$branch' for '$pkg'" * (isnothing(fork) ? "" : " (in fork '$fork')")
   if !isnothing(fork)
      try
         GitHub.reference(github_repo(pkg; fork=fork),
                                 "heads/$branch";
                                 auth=gh_auth)
         @info "  -> found in fork: '$fork'"
         return (pkg_url(pkg; full=true, fork=fork), branch, fork)
      catch
         @info "  -- not found in fork: '$fork'"
      end
   end
   if branch != "master"
      try
         GitHub.reference(github_repo(pkg), "heads/$branch"; auth=gh_auth)
         @info "  -> found in main repo"
         return (pkg_url(pkg; full=true), branch, nothing)
      catch
         @info "  -- not found in main repo"
      end
   end
   # fallback to master branch in default repo
   @warn "  ** branch $branch not found, using default 'master' branch"
   return (pkg_url(pkg; full=true), "master", nothing)
end

# generate a dict describing branches, includes, os and julia versions
# for use in github actions
function ci_matrix(meta::Dict{String,Any}; pr=0, fork=nothing, active_repo=nothing)
   isnothing(gh_auth) && github_auth()

   matrix = Dict{String,Any}(meta["env"])
   active_pkg = pkg_from_repo(active_repo)

   pr_branch = ""
   if pr > 0 && !isnothing(active_pkg) && isnothing(fork)
      @info "fetching $active_pkg PR #$pr."
      ghpr = GitHub.pull_request(github_repo(active_pkg), pr; auth=gh_auth)
      # check if this comes from a fork
      if ghpr.head.ref != "master"
         if ghpr.head.repo.full_name != pkg_url(active_pkg; full=false)
            fork = ghpr.head.user.login
         end
         pr_branch = ghpr.head.ref
      else
         @warn "PR branch name 'master' cannot be used for branch autodetection"
      end
      # TODO: we might even look into pr.body and parse the branch from there?
      # e.g.: look for such a line
      # SomePkg.jl: otherUser/SomePkg.jl#branchname
   end

   # for each package lookup branch with the same name in fork and main repo
   for (pkg,branches) in meta["pkgs"]
      # ignore currently active repo
      pkg == active_pkg && continue
      if !isempty(pr_branch)
         (url, branch, pkg_fork) = find_branch(pkg, pr_branch; fork=fork)
         # don't (re-)add 'master' (even from a fork)
         if branch != "master"
            push!(branches, isnothing(pkg_fork) ?  branch : "$url#$branch")
         end
      end
      if !isempty(branches)
         matrix[pkg] = [Dict("name" => "$pkg#$branch", "branch" => branch)
                           for branch in branches]
      end
   end

   # add includes for custom configurations
   matrix["include"] = []
   for (_,inc) in meta["include"]
      named_include = Dict()
      for (obj,val) in inc
         if obj in ("os","julia-version")
            named_include[obj] = val
         else
            named_include[obj] = Dict("name" => "$obj#$val", "branch" => val)
         end
      end
      push!(matrix["include"],named_include)
   end
   return matrix
end

function parse_job(job_json::AbstractString)
   job_dict = JSON.parse(job_json)
   pkgdict = Dict{String,Any}(pkg => val["branch"]
                                  for (pkg,val) in
                                  filter(p->(!in(first(p), ("os","julia-version"))),
                                         job_dict))
   return pkgdict
end

# this allows setting a github output variable 'matrix'
# which we can then use as input for the matrix-strategy
github_json(github_matrix::Dict{String,Any}) =
   "::set-output name=matrix::" * JSON.json(github_matrix)


text:
function calc_E_ED(iνₙ, ϵₖ, Vₖ, GImp, mP; full=false)
    full && @error "Not implemented for full GF."
    E_kin = 0.0
    E_pot = 0.0
    vk = sum(Vₖ .^ 2)
    Σ_hartree = mP.n * mP.U/2
    E_pot_tail = (mP.U^2)/2 * mP.n * (1-mP.n/2) - Σ_hartree*(Σ_hartree-mP.μ)
    E_kin_tail = vk

    for n in 1:length(GImp)
        Δ_n = sum((Vₖ .^ 2) ./ (iνₙ[n] .- ϵₖ))
        Σ_n = iνₙ[n] .- Δ_n .- 1.0 ./ GImp[n] .+ mP.μ
        E_kin += 2*real(GImp[n] * Δ_n - E_kin_tail/(iνₙ[n]^2))
        E_pot += 2*real(GImp[n] * Σ_n - E_pot_tail/(iνₙ[n]^2))
    end
    E_kin = E_kin .* (2/mP.β) - (mP.β/2) .* E_kin_tail
    E_pot = E_pot .* (1/mP.β) .+ 0.5*Σ_hartree .- (mP.β/4) .* E_pot_tail
    return E_kin, E_pot
end

function calc_E(Σ, qG, mP::ModelParameters, sP::SimulationParameters)
    #println("TODO: E_pot function has to be tested")
    #println("TODO: use GNew/GLoc/GImp instead of Sigma")
    #println("TODO: make frequency summation with sum_freq optional")
    νmax = size(Σ,1)
    νGrid = 0:(νmax-1)
    iν_n = iν_array(mP.β, νGrid)
    Σ_hartree = mP.n * mP.U/2
    Σ_corr = Σ .+ Σ_hartree

    E_kin_tail_c = [zeros(size(qG.ϵkGrid)), (qG.ϵkGrid .+ Σ_hartree .- mP.μ)]
    E_pot_tail_c = [zeros(size(qG.ϵkGrid)),
                    (mP.U^2 * 0.5 * mP.n * (1-0.5*mP.n) .+ Σ_hartree .* (qG.ϵkGrid .+ Σ_hartree .- mP.μ))]
    tail = [1 ./ (iν_n .^ n) for n in 1:length(E_kin_tail_c)]
    E_pot_tail = sum(E_pot_tail_c[i]' .* tail[i] for i in 1:length(tail))
    E_kin_tail = sum(E_kin_tail_c[i]' .* tail[i] for i in 1:length(tail))
    E_pot_tail_inv = sum((mP.β/2)  .* [Σ_hartree .* ones(size(qG.ϵkGrid)), (-mP.β/2) .* E_pot_tail_c[2]])
    E_kin_tail_inv = sum(map(x->x .* (mP.β/2) .* qG.ϵkGrid , [1, -(mP.β) .* E_kin_tail_c[2]]))

    G_corr = flatten_2D(G_from_Σ(Σ_corr, qG.ϵkGrid, νGrid, mP));
    E_pot = real.(G_corr .* Σ_corr .- E_pot_tail);
    E_kin = qG.ϵkGrid' .* real.(G_corr .- E_kin_tail);

    E_pot = [kintegrate(qG, 2 .* sum(E_pot[1:i,:], dims=[1])[1,:] .+ E_pot_tail_inv) for i in 1:νmax]
    E_kin = [kintegrate(qG, 4 .* sum(E_kin[1:i,:], dims=[1])[1,:] .+ E_kin_tail_inv) for i in 1:νmax]
    return E_kin, E_pot
end

"""

Specialized function for DGA potential energy. Better performance than calc_E.
"""
function calc_E_pot(qG, G, Σ, tail, tail_inv)
    E_pot = real.(G .* Σ .- tail);
    return kintegrate(qG, 2 .* sum(E_pot, dims=[1])[1,:] .+ tail_inv)[1]
end

function calc_E_pot_νn(qG, G, Σ, tail, tail_inv)
    E_pot = real.(G .* Σ .- tail);
    return [kintegrate(qG, 2 .* sum(E_pot[1:i,:], dims=[1])[1,:] .+ tail_inv) for i in 1:size(E_pot,1)]
end

function calc_E_kin(qG, G, ϵqGrid, tail, tail_inv)
    E_kin = ϵqGrid' .* real.(G .- tail)
    return kintegrate(qG, 4 .* sum(E_kin, dims=[1])[1,:] .+ tail_inv)
end


text:
module SolveBio
using Requests

include("utils.jl")
include("depository.jl")
include("dataset.jl")

API_HOST = get(ENV, "SOLVEBIO_API_HOST", "https://api.solvebio.com/")
API_KEY = get(ENV, "SOLVEBIO_API_KEY", "")
PKG_VERSION = v"0.0.2"

function login(api_key=API_KEY, api_host=API_HOST)
    # TODO: log in manually with username and password
    try
        res = request("GET", "v1/user")
        println("Logged in to $api_host as $(res["email"]).")
        return res
    catch e
        println("Login failed: $e")
    end
end

function request(method, path, query=Dict(), body=Dict(), args...)
    headers = {
        "Accept" => "application/json",
        "Content-Type" => "application/json",
    }

    if API_KEY != ""
        headers["Authorization"] = "Token $API_KEY"
    end

    # Slice of beginning slash
    if path[1] == '/'
        path = path[2:end]
    end

    useragent = "SolveBio Julia Client $PKG_VERSION [$VERSION $(Sys.MACHINE)]"
    headers["User-Agent"] = useragent

    uri = "$API_HOST$path"

    if method == "GET"
        res = get(uri; headers=headers, query=query, json=body)
    elseif method == "POST"
        res = post(uri; headers=headers, query=query, json=body)
    elseif method == "PUT"
        res = put(uri; headers=headers, query=query, json=body)
    elseif method == "PATCH"
        res = patch(uri; headers=headers, query=query, json=body)
    else
        error("Invalid request method: $method")
    end

    if (res.status < 200 || res.status >= 400)
        if res.status == 429
            error("API error: Too many requests, please retry in $(res.headers["retry-after"]) seconds")
        elseif res.status == 400
            content = formatSolveBioResponse(res, raw=false)
            error("API error: $(res.data)")
        else
            error("API error: $(res.status)")
        end
    end

    if (res.status == 204 || res.status == 301 || res.status == 302)
        return res
    end

    res = formatSolveBioResponse(res)
    return res
end

end