text:
# BAM Record
# ==========

# NOTE: The order and size of fields are important; the data will be loaded
# using memory copy. See specs for the details.
type BAMRecord
    refid::Int32
    pos::Int32
    bin_mq_nl::UInt32
    flag_nc::UInt32
    l_seq::Int32
    next_refid::Int32
    next_pos::Int32
    tlen::Int32

    # variable length data
    data::Vector{UInt8}

    # filled bytes of data (≤ length(.data))
    datasize::Int

    # reference sequence names (shared)
    refseqnames::Vector{String}
end

function BAMRecord()
    return BAMRecord(-1, -1, 0, 0, 0, -1, -1, 0, UInt8[], 0, String[])
end

# NOTE: this does not copy `refseqnames`.
function Base.copy(rec::BAMRecord)
    return BAMRecord(
        rec.refid,
        rec.pos,
        rec.bin_mq_nl,
        rec.flag_nc,
        rec.l_seq,
        rec.next_refid,
        rec.next_pos,
        rec.tlen,
        copy(rec.data),
        rec.datasize,
        rec.refseqnames)
end

function Base.show(io::IO, rec::BAMRecord)
    println(summary(rec), ":")
    println(io, "reference name: ", refname(rec))
    println(io, "next reference name: ", nextrefname(rec))
    println(io, "position: ", position(rec))
    println(io, "next position: ", nextposition(rec))
    println(io, "mapping quality: ", mappingquality(rec))
    println(io, "flag: ", flag(rec))
    println(io, "template length: ", templatelength(rec))
    println(io, "sequence name: ", seqname(rec))
    println(io, "CIGAR string: ", cigar(rec))
    println(io, "sequence: ", sequence(rec))
    println(io, "base qualities: ", qualities(rec))
      print(io, "optional fields: ", optinal_fields(rec))
end

# the data size of fixed-length fields (.refid-.tlen)
const BAM_FIXED_FIELDS_BYTES = 32

"""
    ismapped(rec::BAMRecord)

Return `true` if and only if `rec` is mapped to a reference sequence.
"""
function ismapped(rec::BAMRecord)
    return rec.pos != -1
end

"""
    refid(rec::BAMRecord)

Return the index of a reference sequence that `rec` is mapped onto.

The index is 1-based and will be 0 for an alignment without mapping position.
"""
function refid(rec::BAMRecord)
    return rec.refid + 1
end

function nextrefid(rec::BAMRecord)
    return rec.next_refid + 1
end

"""
    refname(rec::BAMRecord)

Return the name of a reference sequence that `rec` is mapped onto.

If `rec` is unmapped, it returns `"*"` like SAM records.
"""
function refname(r::BAMRecord)
    id = refid(r)
    if id == 0
        return "*"
    else
        return r.refseqnames[id]
    end
end

function nextrefname(r::BAMRecord)
    id = nextrefid(r)
    if id == 0
        return "*"
    else
        return r.refseqnames[id]
    end
end

"""
    position(rec::BAMRecord)

Return the position of a mapped read.

The index is 1-based and will be 0 for an alignment without mapping position.
"""
function Base.position(rec::BAMRecord)
    return rec.pos + 1
end

function nextposition(rec::BAMRecord)
    return rec.next_pos + 1
end

"""
    mappingquality(rec::BAMRecord)

Return the mapping quality of the alignment `rec`.
"""
function mappingquality(rec::BAMRecord)
    return UInt8((rec.bin_mq_nl >> 8) & 0xff)
end

"""
    flag(rec::BAMRecord)

Return the flag of the alignment `rec`.
"""
function flag(rec::BAMRecord)
    return UInt16(rec.flag_nc >> 16)
end

"""
    templatelength(rec::BAMRecord)

Return the template length of the alignment `rec`.
"""
function templatelength(rec::BAMRecord)
    return rec.tlen
end

"""
    seqname(rec::BAMRecord)

Return the read name of the alignment `rec`.
"""
function seqname(rec::BAMRecord)
    # drop the last NUL character
    return unsafe_string(pointer(rec.data), max(seqname_length(rec) - 1, 0))
end

"""
    cigar_rle(rec::BAMRecord)

Return a run-length encoded tuple `(ops, lens)` of the CIGAR string.
See also `cigar`.
"""
function cigar_rle(rec::BAMRecord)
    offset = seqname_length(rec)
    ops = Bio.Align.Operation[]
    lens = Int[]
    for i in offset+1:4:offset+n_cigar_op(rec)*4
        x = unsafe_load(Ptr{UInt32}(pointer(rec.data, i)))
        op = Bio.Align.Operation(x & 0x0f)
        push!(ops, op)
        push!(lens, x >> 4)
    end
    return ops, lens
end

"""
    cigar(rec::BAMRecord)

Return a CIGAR string of the alignment `rec`. See also `cigar_rle`.
"""
function cigar(rec::BAMRecord)
    buf = IOBuffer()
    for (op, len) in zip(cigar_rle(rec)...)
        print(buf, len, Char(op))
    end
    return takebuf_string(buf)
end

"""
    sequence(rec::BAMRecord)

Return a DNA sequence of the alignment `rec`.
"""
function sequence(rec::BAMRecord)
    seqlen = sequence_length(rec)
    data = Vector{UInt64}(cld(seqlen, 16))
    src::Ptr{UInt64} = pointer(rec.data, seqname_length(rec) + n_cigar_op(rec) * 4 + 1)
    for i in 1:endof(data)
        # copy data flipping high and low nybble
        x = unsafe_load(src, i)
        data[i] = (x & 0x0f0f0f0f0f0f0f0f) << 4 | (x & 0xf0f0f0f0f0f0f0f0) >> 4
    end
    return Bio.Seq.DNASequence(data, 1:seqlen, false)
end

"""
    qualities(rec::BAMRecord)

Return base qualities of the alignment `rec`.
"""
function qualities(rec::BAMRecord)
    seqlen = sequence_length(rec)
    offset = seqname_length(rec) + n_cigar_op(rec) * 4 + cld(seqlen, 2)
    return [rec.data[i+offset] for i in 1:seqlen]
end

function Base.getindex(rec::BAMRecord, tag::AbstractString)
    checkkeytag(tag)
    return getvalue(rec.data, auxdata_position(rec), UInt8(tag[1]), UInt8(tag[2]))
end

function Base.setindex!(rec::BAMRecord, val, tag::AbstractString)
    checkkeytag(tag)
    setvalue!(rec.data, auxdata_position(rec), val, UInt8(tag[1]), UInt8(tag[2]))
    return rec
end

function Base.delete!(rec::BAMRecord, tag::AbstractString)
    checkkeytag(tag)
    deletevalue!(rec.data, auxdata_position(rec), UInt8(tag[1]), UInt8(tag[2]))
    return rec
end

function Base.haskey(rec::BAMRecord, tag::AbstractString)
    checkkeytag(tag)
    return findtag(rec.data, auxdata_position(rec), UInt8(tag[1]), UInt8(tag[2])) > 0
end

function optinal_fields(rec::BAMRecord)
    return AuxDataDict(rec.data[auxdata_position(rec):rec.datasize])
end

function auxdata_position(rec)
    seqlen = sequence_length(rec)
    return seqname_length(rec) + n_cigar_op(rec) * 4 + cld(seqlen, 2) + seqlen + 1
end

# Return the right-most position of alignment.
function rightmost_position(rec::BAMRecord)
    return Int32(position(rec) + alignment_length(rec) - 1)
end

# Return the length of alignment.
function alignment_length(rec::BAMRecord)
    offset = seqname_length(rec)
    length::Int = 0
    for i in offset+1:4:offset+n_cigar_op(rec)*4
        x = unsafe_load(Ptr{UInt32}(pointer(rec.data, i)))
        op = Bio.Align.Operation(x & 0x0f)
        if Bio.Align.ismatchop(op) || Bio.Align.isdeleteop(op)
            length += x >> 4
        end
    end
    return length
end

# Return the length of the read name.
function seqname_length(rec)
    return rec.bin_mq_nl & 0xff
end

# Return the number of CIGAR operations.
function n_cigar_op(rec)
    return rec.flag_nc & 0xffff
end

# Return the length of the DNA sequence.
function sequence_length(rec)
    return rec.l_seq
end


text:
using DataFrames
using CSV
using Random

@enum db_type output csv

mutable struct dbEngine
    id::String
    type:: db_type
    log:: Bool
    rows_df:: DataFrames.DataFrame
    results:: DataFrames.DataFrame
    meta :: Dict{String,String}
end

function db()
    id = randstring(12)
    return dbEngine(id, output,true, DataFrame(), DataFrame(), Dict())
end

function meta(db::dbEngine, key::String, value::String)
    push!(db.meta, (key => value))
end
# Note that constructing a DataFrame row by row is significantly less performant than constructing it all at once, or column by column. For many use-cases this will not matter, but for very large DataFrames this may be a consideration.
function setup(db::dbEngine)
    @debug "DB is setupped."
end

function finish(db::dbEngine)
    @debug "DB is finished."
    dbSaveData(db)
end

function store(db::dbEngine)
    date_format = Dates.format(now(), "yymmddHHMMSS_s")
    id = db.id
    CSV.write("results/rows.csv", db.rows_df, append=true)
    CSV.write("results/results.csv", db.results, append=true)
end

function info(db::dbEngine, text, data)
    @debug text data...
end

function record(db::dbEngine, data)
    info(db, "New row is inserted.", data )
    push!(db.rows_df, data)
end

function result(db::dbEngine, data)
    @debug "Result" data...
    push!(db.results, data)
end

ext: jl

max_stars_repo_name: ahrajabi/DiscO

max_issues_repo_name: ahrajabi/DiscO

max_forks_repo_name: ahrajabi/DiscO


text:
wdot(wx, x, wy, y) = dot(wx.*x, wy.*y)
wnorm(w, x) = norm(w.*x)

assess_convergence(f, ftol) = assess_convergence(NaN, NaN, f, NaN, ftol)
function assess_convergence(x,
                            x_previous,
                            f,
                            xtol,
                            ftol)
    x_converged, f_converged = false, false
    if norm(x-x_previous) <= xtol
        x_converged = true
    end
    if maximum(abs, f) <= ftol
        f_converged = true
    end

    return x_converged, f_converged
end

function check_isfinite(x::AbstractArray)
    if any(!isfinite, x)
        i = findall(!isfinite, x)
        throw(IsFiniteException(i))
    end
end

"""
    qrdelete!(Q, R, k)

Delete the left-most column of F = Q[:, 1:k] * R[1:k, 1:k] by updating Q and R.

Only Q[:, 1:(k-1)] and R[1:(k-1), 1:(k-1)] are valid on exit.
"""
function qrdelete!(Q::AbstractMatrix, R::AbstractMatrix, k::Int)
  n, m = size(Q)
  m == LinearAlgebra.checksquare(R) || throw(DimensionMismatch())
  1 ≤ k ≤ m || throw(ArgumentError())

  # apply Givens rotations
  for i in 2:k
      g = first(givens(R, i - 1, i, i))
      lmul!(g, R)
      rmul!(Q, g')
  end

  # move columns of R
  @inbounds for j in 1:(k-1)
    for i in 1:(k-1)
      R[i, j] = R[i, j + 1]
    end
  end

  Q, R
end

"""
    qradd!(Q, R, v, k)

Replace the right-most column of F = Q[:, 1:k] * R[1:k, 1:k] with v by updating Q and R.

This implementation modifies vector v as well. Only Q[:, 1:k] and R[1:k, 1:k] are valid on
exit.
"""
function qradd!(Q::AbstractMatrix, R::AbstractMatrix, v::AbstractVector, k::Int)
  n, m = size(Q)
  n == length(v) || throw(DimensionMismatch())
  m == LinearAlgebra.checksquare(R) || throw(DimensionMismatch())
  1 ≤ k ≤ m || throw(ArgumentError())

  @inbounds for i in 1:(k-1)
    q = view(Q, :, i)
    r = dot(q, v)

    R[i, k] = r
    axpy!(-r, q, v)
  end

  @inbounds begin
    d = norm(v)
    R[k, k] = d
    @. Q[:, k] = v / d
  end

  Q, R
end


text:
using LinearAlgebra
using MPI
using Random
using Rimu
using StaticArrays
using Statistics
using Test

using Rimu.RMPI
using Rimu.StatsTools
using Rimu.ConsistentRNG
using Rimu.RMPI: targetrank, mpi_synchronize!

const N_REPEATS = 5

"""
    rand_onr(N, M)

Generate random ONR with N particles in M sites.
"""
function rand_onr(N, M)
    result = zeros(MVector{M,Int})
    for _ in 1:N
        result[rand(1:M)] += 1
    end
    return result.data
end

"""
    correct_ranks(md)

Check if all entries in `md` are located on the correct rank.
"""
function correct_ranks(md)
    return mapreduce(|, pairs(md); init=true) do ((k, v))
        targetrank(k, mpi_size()) == mpi_rank()
    end
end

# Ignore all printing on ranks other than root.
if mpi_rank() != mpi_root
    redirect_stderr(devnull)
    redirect_stdout(devnull)
end

"""
    setup_dv(type, args...; kwargs...)

Create a local and distributed versions of dvec of type `type`. Ensure they both have the same
contents.
"""
function setup_dv(type, args...; md_kwargs=(;), kwargs...)
    v = type(args...; kwargs...)
    if mpi_rank() == mpi_root
        dv = MPIData(copy(v); md_kwargs...)
    else
        dv = MPIData(empty(v); md_kwargs...)
    end
    return v, dv
end

@testset "MPI tests" begin
    @testset "MPIData" begin
        for type in (InitiatorDVec, DVec)
            @testset "copy_to_local" begin
                dv = MPIData(type(mpi_rank() => 1, -1 => 10))
                loc = RMPI.copy_to_local(dv)
                @test length(loc) == mpi_size() + 1
                @test loc[-1] == 10 * mpi_size()
            end
            @testset "Single component $type" begin
                for i in 1:N_REPEATS
                    add = BoseFS((0,0,10,0,0))
                    H = HubbardMom1D(add)
                    Random.seed!(7350 * i)
                    v, dv = setup_dv(
                        type, [BoseFS(rand_onr(10, 5)) => 2 - 4rand() for _ in 1:100]
                    )
                    mpi_synchronize!(dv)

                    Random.seed!(1337 * i)
                    w, dw = setup_dv(
                        type, [BoseFS(rand_onr(10, 5)) => 2 - 4rand() for _ in 1:20]
                    )
                    mpi_synchronize!(dw)

                    @test correct_ranks(dv)
                    @test length(v) == length(dv)
                    @test correct_ranks(dw)
                    @test length(w) == length(dw)

                    @testset "Basics" begin
                        @test norm(v) ≈ norm(dv)
                        @test norm(v, 1) ≈ norm(dv, 1)
                        @test norm(v, 2) ≈ norm(dv, 2)
                        @test norm(v, Inf) ≈ norm(dv, Inf)
                        @test sum(values(v)) ≈ sum(values(dv))
                        f((k, v)) = (k == add) + v > 0
                        @test mapreduce(f, |, pairs(v); init=true) ==
                            mapreduce(f, |, pairs(dv); init=true)
                    end

                    @testset "Operations" begin
                        @test dot(v, w) ≈ dot(dv, dw)

                        @test dot(v, H, w) ≈ dot(v, H, dw)
                        @test dot(w, H, v) ≈ dot(w, H, dv)
                        @test dot(w, H, v) ≈ dot(dw, H, dv)
                        @test dot(w, H, v) ≈ dot(dw, H, v)

                        @test dot(freeze(dw), dv) ≈ dot(w, v)
                        @test dot(freeze(dv), dw) ≈ dot(v, w)

                        du = MPIData(H * dv)
                        u = H * v
                        @test correct_ranks(du)

                        @test length(u) == length(du)
                        @test norm(u, 1) ≈ norm(du, 1)
                        @test norm(u, 2) ≈ norm(du, 2)
                        @test norm(u, Inf) ≈ norm(du, Inf)
                    end
                end
            end
            @testset "Two-component $type" begin
                for i in 1:N_REPEATS
                    add = BoseFS2C((0,0,10,0,0), (0,0,2,0,0))
                    H = BoseHubbardMom1D2C(add)
                    Random.seed!(7350 * i)
                    v, dv = setup_dv(
                        type, [BoseFS2C(rand_onr(10, 5), rand_onr(2, 5)) => rand() for _ in 1:100]
                    )
                    mpi_synchronize!(dv)

                    Random.seed!(1337 * i)
                    w, dw = setup_dv(
                        type, [BoseFS2C(rand_onr(10, 5), rand_onr(2, 5)) => rand() for _ in 1:20]
                    )
                    mpi_synchronize!(dw)

                    @testset "Operations" begin
                        @test dot(v, w) ≈ dot(dv, dw)

                        @test dot(v, H, w) ≈ dot(v, H, dw)
                        @test dot(w, H, v) ≈ dot(w, H, dv)
                        G1 = G2Correlator(1)
                        G3 = G2Correlator(3)
                        @test dot(v, G1, w) ≈ dot(v, G1, dw)
                        @test dot(w, G3, v) ≈ dot(w, G3, dv)
                        @test dot(w, G3, v) ≈ dot(dw, G3, v)
                        @test dot(w, G3, v) ≈ dot(dw, G3, dv)

                        @test dot(freeze(dw), dv) ≈ dot(w, v)
                        @test dot(freeze(dv), dw) ≈ dot(v, w)

                        du = MPIData(H * dv)
                        u = H * v
                        @test correct_ranks(du)

                        @test length(u) == length(du)
                        @test norm(u, 1) ≈ norm(du, 1)
                        @test norm(u, 2) ≈ norm(du, 2)
                        @test norm(u, Inf) ≈ norm(du, Inf)

                        du = MPIData(G3 * dv)
                        u = G3 * v
                        @test correct_ranks(du)

                        @test length(u) == length(du)
                        @test norm(u, 1) ≈ norm(du, 1)
                        @test norm(u, 2) ≈ norm(du, 2)
                        @test norm(u, Inf) ≈ norm(du, Inf)
                    end
                end
            end
        end
        @testset "Communication strategies" begin
            # The idea here is to generate a big DVec and make sure all communication
            # strategies distribute it in the same manner.
            for i in 1:N_REPEATS
                sorted = map((
                    (; setup=RMPI.mpi_point_to_point),
                    (; setup=RMPI.mpi_one_sided, capacity=10_000),
                    (; setup=RMPI.mpi_all_to_all),
                )) do kw
                    Random.seed!(2021 * hash(mpi_rank()) * i)
                    source = DVec([BoseFS(rand_onr(10, 5)) => 2 - 4rand() for _ in 1:10_000])
                    target = MPIData(similar(source); kw...)
                    RMPI.mpi_combine_walkers!(target, source)
                    return target
                end

                for v in sorted
                    @test correct_ranks(v)
                end
                @test localpart(sorted[1]) == localpart(sorted[2])
                @test localpart(sorted[2]) == localpart(sorted[3])
            end
        end
    end

    @testset "Ground state energy estimates" begin
        # H1 = HubbardReal1D(BoseFS((1,1,1,1,1,1,1)); u=6.0)
        # E0 = eigsolve(H1, DVec(starting_address(H1) => 1.0), 1, :SR; issymmetric=true)[1][1]
        E0 = -4.6285244934941305
        mpi_seed_CRNGs!(1000_000_000)

        for (setup, kwargs) in (
            (RMPI.mpi_point_to_point, (;)),
            (RMPI.mpi_all_to_all, (;)),
            (RMPI.mpi_one_sided, (; capacity=1000)),
        )
            @testset "Regular with $setup and post-steps" begin
                H = HubbardReal1D(BoseFS((1,1,1,1,1,1,1)); u=6.0)
                dv = MPIData(
                    DVec(starting_address(H) => 3; style=IsDynamicSemistochastic());
                    setup,
                    kwargs...
                )

                post_step = (
                    ProjectedEnergy(H, dv),
                    SignCoherence(copy(localpart(dv))),
                    WalkerLoneliness(),
                    Projector(proj_1=Norm2Projector()),
                )
                df = lomc!(H, dv; post_step, laststep=5000).df

                # Shift estimate.
                Es, σs = mean_and_se(df.shift[2000:end])
                s_low, s_high = Es - 2σs, Es + 2σs
                # Projected estimate.
                r = ratio_of_means(df.hproj[2000:end], df.vproj[2000:end])
                p_low, p_high = pquantile(r, [0.0015, 0.9985])

                @test s_low < E0 < s_high
                @test p_low < E0 < p_high
                @test all(-1 .≤ df.coherence .≤ 1)
                @test all(0 .≤ df.loneliness .≤ 1)
            end
            @testset "Initiator with $setup" begin
                H = HubbardMom1D(BoseFS((0,0,0,7,0,0,0)); u=6.0)
                dv = MPIData(
                    InitiatorDVec(starting_address(H) => 3);
                    setup,
                    kwargs...
                )
                s_strat = DoubleLogUpdate(targetwalkers=100)
                df = lomc!(H, dv; laststep=5000, s_strat).df

                # Shift estimate.
                Es, _ = mean_and_se(df.shift[2000:end])
                @test E0 ≤ Es
            end
        end
    end

    @testset "Same seed gives same results" begin
        # Note: all_to_all will give different results here, as it sorts slightly differently.
        # The entries in the vectors are the same (tested above), but they appear in a
        # different order.
        add = BoseFS2C((1, 1, 1, 1, 1), (1, 0, 0, 0, 0))
        H = BoseHubbardReal1D2C(add; v=2)
        dv = DVec(add => 1)

        dv_ptp = MPIData(copy(dv); setup=RMPI.mpi_point_to_point)
        mpi_seed_CRNGs!(17)
        df_ptp = lomc!(H, dv_ptp).df

        dv_os = MPIData(copy(dv); setup=RMPI.mpi_one_sided, capacity=1000)
        mpi_seed_CRNGs!(17)
        df_os = lomc!(H, dv_os).df

        @test df_ptp == df_os
    end

    # Make sure all ranks came this far.
    @testset "Finish" begin
        @test MPI.Allreduce(true, &, mpi_comm())
    end

    @testset "example script" begin
        include("../scripts/BHM-example-mpi.jl")
        dfr = load_df("result.arrow")
        qmcdata = last(dfr, 5000)
        (qmcShift,qmcShiftErr) = mean_and_se(qmcdata.shift)
        @test qmcShift ≈ -6.5 atol=0.5

        # clean up
        rm("result.arrow", force=true)
    end
end

if mpi_rank() ≠ mpi_root
    redirect_stderr(stderr)
    redirect_stdout(stdout)
end


text:
# FIXME move to DFG
getPointDefault(V::InferenceVariable) = getPointIdentity(V)

function compare( c1::Channel,
                  c2::Channel;
                  skip::Vector{Symbol}=[] )
  #
  TP = true
  TP = TP && c1.state == c2.state
  TP = TP && c1.sz_max == c2.sz_max
  TP = TP && c1.data |> length == c2.data |> length
  # exit early if tests already failed
  !TP && (return false)
  # now check contents of data
  for i in 1:length(c1.data)
  TP = TP && c1.data[i] == c2.data[i]
  end
  return TP
end

compare(a::Int, b::Int) = a == b
compare(a::Bool, b::Bool) = a == b
compare(a::Dict, b::Dict) = a == b


text:
"""
    $(SIGNATURES)

Save a figure to a pdf file.

# Arguments

- p: Figure.
- io: If provided, display a short version of the file path there.
- `dataM`: If provided, write the data to a text file. A simple method of preserving underlying data.
"""
function figsave(p, filePath :: String;
    figNotes = nothing, io :: Union{Nothing, IO} = nothing,
    dataM = nothing)

    @assert !(p isa Tuple)  "Tuple input $p. Expecting a Figure";
    if !isempty(filePath)
		newPath = change_extension(filePath, ".pdf");
        fDir, fName = splitdir(newPath);
        if !isdir(fDir)
            mkpath(fDir);
        end
        # This prevents trouble with corrupt files
        isfile(newPath)  &&  rm(newPath);
        save(newPath, p);
        if !isnothing(io)
            pathStr = fpath_to_show(filePath);
            println(io, "Saved figure:  $pathStr");
        end

        save_fig_notes(figNotes, filePath);
        save_fig_data(dataM, filePath);
	end
end

figsave(p, fDir :: String, fName :: String; figNotes = nothing) =
    figsave(p, joinpath(fDir, fName); figNotes = figNotes);

## --------  These are indpendent of the plotting package

save_fig_notes(figNotes :: AbstractVector, fPath :: AbstractString) =
    save_text_file(fig_notes_path(fPath), figNotes);
save_fig_notes(figNotes :: AbstractString, fPath :: AbstractString) =
    save_fig_notes([figNotes], fPath);
function save_fig_notes(figNotes :: Nothing, fPath) end

function fig_notes_path(fPath :: AbstractString)
    newDir = fig_data_dir(fPath);
    fDir, fName = splitdir(fPath);
    newPath = joinpath(newDir, change_extension(fName, ".txt"));
    return newPath
end

# Saves the data in text format.
function save_fig_data(dataM, filePath :: AbstractString)
    newPath = fig_data_path(filePath);
    make_dir(newPath);
    if !isnothing(dataM)
        open(newPath, "w") do io
            println(io, dataM);
        end
    end
end

function fig_data_path(fPath :: AbstractString)
    newDir = fig_data_dir(fPath);
    fDir, fName = splitdir(fPath);
    newPath = joinpath(newDir, change_extension("data_" * fName, ".txt"));
    return newPath
end

function fig_data_dir(fPath :: AbstractString)
    fDir, fName = splitdir(fPath);
    newDir = joinpath(fDir, "plot_data");
    isdir(newDir)  ||  mkpath(newDir);
    return newDir
end

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


text:
# # 4-Potential

function EB(r, t, laser)
    R = geometry(laser).rotation_matrix
    r′ = rotate_coords(R, r)
    inv_c = immutable_cache(laser, :inv_c)
    if unit(eltype(r)) ≠ NoUnits
        λ = laser.λ
        r′ = uconvert.(unit(λ), r′)
    end

    E_B = real.(EB(r′, laser) .* g(r′[end], t, laser; inv_c))

    return inv_rotate.((R,), E_B)
end

function EB(r, laser)
    @assert length(r) == 3 "The laser is only defined in 3D"
    fill!(mutable_cache(laser), r)
    coords = required_coords(laser, r)

    E_x = Ex(laser, coords)
    update_cache!(laser, :Ex, E_x)
    E_y = Ey(laser, coords)
    update_cache!(laser, :Ey, E_y)
    E_z = Ez(laser, coords)
    update_cache!(laser, :Ez, E_z)

    ElectricField = SVector{3}(E_x, E_y, E_z)

    B_x = Bx(laser, coords)
    B_y = By(laser, coords)
    B_z = Bz(laser, coords)

    MagneticField = SVector{3}(B_x, B_y, B_z)

    return ElectricField, MagneticField
end

"""
    Fμν(x, laser)

Compute [the electromagnetic tensor](https://en.wikipedia.org/wiki/Electromagnetic_tensor) ``F`` in
covarinat matrix form at the spacetime 4-vector specified by `x`.
The ``x = (x^0 ≡ ct, x^1, x^2, x^3)`` convention is used.
```math
F^{μν} =
\\left[\\begin{array}{cccc}
    0   &  -E_x / c  & -E_y / c & -E_z / c \\\\
E_x / c &       0    &   -B_z   &    B_y   \\\\
E_y / c &      B_z   &     0    &   -B_x   \\\\
E_z / c &     -B_y   &    B_z   &     0
\\end{array}\\right]
```
"""
function Fμν(x, laser)
    c⁻¹ = immutable_cache(laser, :inv_c)
    r = x[begin+1:end]
    t = x[begin] * c⁻¹
    (Ex, Ey, Ez), (Bx, By, Bz) = EB(r, t, laser)

    𝟘 = zero(Bx)
    return @SMatrix [𝟘      -Ex*c⁻¹  -Ey*c⁻¹ -Ez*c⁻¹ ;
                     Ex*c⁻¹   𝟘      -Bz      By     ;
                     Ey*c⁻¹   Bz      𝟘      -Bx     ;
                     Ez*c⁻¹  -By      Bx      𝟘      ]
end


text:
function _transpose_front(t::AbstractTensorMap) # make TensorMap{S,N₁+N₂-1,1}
    I1 = TensorKit.codomainind(t)
    I2 = TensorKit.domainind(t)
    transpose(t, (I1..., reverse(Base.tail(I2))...), (I2[1],))
end
function _transpose_tail(t::AbstractTensorMap) # make TensorMap{S,1,N₁+N₂-1}
    I1 = TensorKit.codomainind(t)
    I2 = TensorKit.domainind(t)
    transpose(t, (I1[1],), (I2..., reverse(Base.tail(I1))...))
end
function _transpose_as(t1::AbstractTensorMap, t2::AbstractTensorMap{S,N1,N2}) where {S,N1,N2}
    I1 = (TensorKit.codomainind(t1)...,reverse(TensorKit.domainind(t1))...);

    A = ntuple(x->I1[x],N1);
    B = ntuple(x->I1[x+N1],N2);

    transpose(t1, A,B)
end

_firstspace(t::AbstractTensorMap) = space(t, 1)
_lastspace(t::AbstractTensorMap) = space(t, numind(t))

#given a hamiltonian with unit legs on the side, decompose it using svds to form a "localmpo"
function decompose_localmpo(inpmpo::AbstractTensorMap{PS,N,N},trunc = truncbelow(Defaults.tol)) where {PS,N}
    N == 2 && return [transpose(inpmpo,(1,2),(3,4))]

    leftind = (N+1,1,2)
    rightind = (ntuple(x->x+N+1,N-1)...,reverse(ntuple(x->x+2,N-2))...);
    (U,S,V) = tsvd(transpose(inpmpo,leftind,rightind),trunc = trunc)

    A = transpose(U,(2,3),(1,4));
    B = transpose(S*V,(1,reverse(ntuple(x->x+N,N-2))...),ntuple(x->x+1,N-1))
    return [A;decompose_localmpo(B)]
end

function add_util_leg(tensor::AbstractTensorMap{S,N1,N2}) where {S,N1,N2}
    ou = oneunit(_firstspace(tensor));

    util_front = isomorphism(storagetype(tensor),ou*codomain(tensor),codomain(tensor));
    util_back = isomorphism(storagetype(tensor),domain(tensor),domain(tensor)*ou);

    return util_front*tensor*util_back
end

"""
Take the L2 Tikhonov regularised inverse of a matrix `m`.

The regularisation parameter is the larger of `delta` (the optional argument that defaults
to zero) and square root of machine epsilon. The inverse is done using an SVD.
"""
function reginv(m, delta=zero(eltype(m)))
    delta = max(abs(delta), sqrt(eps(real(float(one(eltype(m)))))))
    U, S, Vdg = tsvd(m)
    Sinv = inv(real(sqrt(S^2 + delta^2*one(S))))
    minv = Vdg' * Sinv * U'
    return minv
end

function union_split(a::AbstractArray)
    T = reduce((a,b)->Union{a,b},typeof.(a))
    nA = similar(a,T);
    copy!(nA,a)
end
union_types(x::Union) = (x.a, union_types(x.b)...)
union_types(x::Type) = (x,)

function _embedders(spaces)
    totalspace = reduce(⊕,spaces);

    maps = [isometry(totalspace,first(spaces))];
    restmap = leftnull(first(maps));

    for sp in spaces[2:end]
        cm = isometry(domain(restmap),sp);

        push!(maps,restmap*cm);
        restmap = restmap*leftnull(cm);
    end

    maps
end

function _can_unambiguously_braid(sp::VectorSpace)
    s = sectortype(sp);

    #either the braidingstyle is bosonic
    BraidingStyle(s) isa Bosonic && return true

    #or there is only one irrep ocurring - the trivial one
    for sect in sectors(sp)
        sect == one(sect) || return false
    end
    return true
end


text:
import Term: tprintln
using Markdown

text = md"""
 # Markdown rendering in Term.jl
 ## two
 ### three
 #### four
 ##### five
 ###### six

 This is an example of markdown content rendered in Term.jl.
 You can use markdown syntax to make words **bold** and *italic* or insert `literals`.

 You markdown can include in-line latex ``\LaTeX  \frac{{1}}{{2}}`` and maths in a new line too:

 ```math
 f(a) = \frac{1}{2\pi}\int_{0}^{2\pi} (\alpha+R\cos(\theta))d\theta
 ```

 You can also have links: [Julia](http://www.julialang.org) and
 footnotes [^1] for your content [^named].

 And, of course, you can show some code too:

 ```julia
 function say_hi(x)
     print("Hellow World")
 end
 ```

 ---

 You can use "quotes" to highlight a section:

 > Multi-line quotes can be helpful to make a
 > paragram stand out, so that users won't miss it!
 > You can use **other inline syntax** in you `quotes` too.

 but if you really need to grab someone's attention, use admonitions:

 !!! note
     You can use different levels

 !!! warning
     to send different messages

 !!! danger
     to your reader

 !!! tip "Wow!"
     Turns out that admonitions can be pretty useful!
     What will you use them for?

 ---

 Of course you can have classic lists:
 * item one
 * item two
 * And a sublist:
     + sub-item one
     + sub-item two

 and ordered lists too:
 1. item one
 2. item two
 3. item three

 !!! note "Tables"
     You can use the [Markdown table syntax](https://www.markdownguide.org/extended-syntax/#tables)
     to insert tables - Term.jl will convert them to Table object!

 | Term | handles | tables|
 |:---------- | ---------- |:------------:|
 | Row `1`    | Column `2` |              |
 | *Row* 2    | **Row** 2  | Column ``3`` |

 ----

 This is where you print the content of your foot notes:

 [^1]: Numbered footnote text.

 [^note]:
     Named footnote text containing several toplevel elements.

 """

tprintln(text)

tprintln(md"""
Or you can use:

```
using Term: RenderableText
RenderableText(text)
```
""")


text:
# tangle.jl:
# =========
# Convert Julia AST:s to DAG:s in dag.jl format
#

load("utils/req.jl")
req("dag/dag.jl")

# == TangleContext ============================================================

typealias SymbolTable Dict{Symbol,Node}

type TangleContext
    symbols::SymbolTable  # current symbol bindings
    last_line::Node
    last_actions#::Nodes  # todo: add back once type inference bug is fixed

    TangleContext() = new(SymbolTable(), NoNode(), ActionNode[])
end

function emit_line(c::TangleContext, node::Node)
    c.last_line = isa(c.last_line, NoNode) ? node : KnotNode(c.last_line, node)
end
function emit_line(c::TangleContext, node::Node, name::Symbol)
    emitted = emit_line(c, node); emitted.name = name
    emitted
end

# == tangle ===================================================================

tangle(code) = raw_tangle(code)[1]
function raw_tangle(code)
    context = TangleContext()
    value = tangle(context, code)
    if is((value::Node).name, nothing);  value.name = :value;  end
#    sink = TupleNode(context.last_line, context.last_actions..., value)
    sink = TupleNode(context.last_actions..., value)
    #sink.name = :sink
    sink, context
end

tangle(c::TangleContext, exps::Vector) = Node[{tangle(c, ex) | ex in exps}...]

tangle(c::TangleContext, ex::Any) = LiteralNode(ex)   # literal
function tangle(c::TangleContext, name::Symbol)
    @setdefault c.symbols[name] SymNode(name, name==:(...) ? :symbol : :input)
end
function tangle(context::TangleContext, ex::Expr)
    if ex.head == :block    # exprs...
        value = LiteralNode(nothing)
        for subex in ex.args;  value = tangle(context, subex);  end

        return value
    elseif ex.head == :(=)  # assignment: lvalue = expr
        lhs = tangle_lhs(context, ex.args[1])
        rhs = tangle(context, ex.args[2])

        return entangle_assignment(context, lhs, rhs)
    elseif ex.head == :call
        fname = ex.args[1]
        op = @setdefault context.symbols[fname] SymNode(fname, :call)
        args = tangle(context, ex.args[2:end])

        return CallNode(op, args...)
    elseif ex.head == :ref;   return RefNode(tangle(context, ex.args)...)
    # todo: test!:
    elseif ex.head == :(...); return EllipsisNode(tangle(context, ex.args)...)
    # todo: test!:
    elseif ex.head == :tuple; return TupleNode(tangle(context, ex.args)...)
    elseif ex.head == :line;  return NoNode() # ignore line numbers
    end
    error("unexpected scalar rhs: ex = $ex")
end

# -- tangle_lhs --------------------------------------------------------------
# Tangle a scalar-valued lhs expr; return lone symbols unadorned

tangle_lhs(context::TangleContext, name::Symbol) = name
function tangle_lhs(context::TangleContext, ex::Expr)
    # lvalue for indexed assignment
    @expect is_expr(ex, :ref)
    oname = ex.args[1]
    output = @setdefault context.symbols[oname] SymNode(oname, :output)
    inds = tangle(context, ex.args[2:end])
    RefNode(output, inds...)
end

# -- entangle_assignment(context::TangleContext, lhs, rhs) --------------------
# Process assignment lhs = rhs.
# Returns value = rhs

function entangle_assignment(context::TangleContext, dest::Symbol, rhs::Node)
    # straight assignment: just store in symbol table
    context.symbols[dest] = emit_line(context, rhs, dest) # return rhs
end
function entangle_assignment(context::TangleContext, lhs::RefNode, rhs::Node)
    # Indexed assignment to output
    dest = get_A(lhs)::SymNode
    node = AssignNode(lhs, rhs, context.last_actions...)
    #     bind the AssignNode to the name of dest
    context.symbols[dest.val.name] = emit_line(context, node)
    context.last_actions = [node]
    rhs # and evaluate to the rhs
end

# == Some printing ============================================================

function print_list(list::Vector)
    for item in list; println("\t", item); end
end
print_symtable(st::SymbolTable) = (for (k, v) in st; println("\t$k = $v"); end)

function print_context(context::TangleContext)
    println("Symbols at end:")
    print_symtable(context.symbols)
end

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
module Foo

say_hello() = "hello"

end


text:
module Test_Collision_StarSetting

using Modia3D
using Modia3D.StaticArrays
using Modia3D.ModiaMath
using Modia3D.Unitful

# Materials
massProperties = "Aluminium"
vmat1   = Modia3D.Material(color="LightBlue", transparency=0.5, wireframe=false)
vblue   = Modia3D.Material(color="LightBlue" , transparency=0.5)
vred    = Modia3D.Material(color="Red", transparency=0.5)
vgreen  = Modia3D.Material(color="Green", transparency=0.5)
vblack  = Modia3D.Material(color="Black", transparency=0.5)
vyellow = Modia3D.Material(color="Yellow", transparency=0.5)
vpink   = Modia3D.Material(color="Pink", transparency=0.5)
cmat    = "Steel"
vMatAABB = Modia3D.Material(color="Grey", transparency=0.95)

# dimensions
Lx  = 0.7
Dx = 0.5
Ly    = 0.2
Lz   = 1.0

groundWidth  = 1.0
groundHeight = 0.01
cyl  = Modia3D.Cylinder(1.2*Ly, Ly/2)
Tcyl = ModiaMath.rot1(90u"°")
TGround = ModiaMath.rot1(-90u"°")

@assembly CollisionSolids(world) begin
world_f1 = Modia3D.Object3D(world , r=[0.5*Lx, 0.0, groundWidth/2])

blau1     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vblue;contactMaterial = cmat),    [ [0.0, 0.0, -Lz/2],[ 0.0, 0.0, Lz/2] ] )
blau2     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vblue; contactMaterial = cmat),   [ [0.0, 0.0, -Lz/2], [ 0.0, 0.0, Lz/2] ] )

schwarz1  = Modia3D.Part(Modia3D.Solid(Modia3D.SolidCone(Dx,Lz; relativeTipDiameter=0.0), "Aluminium", vblack; contactMaterial = cmat),  [ [0.0, 0.0, 0.0]] )
rosa1     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vpink; contactMaterial = cmat),   [ [0.0, 0.0, -Lz/2],[ 0.0, 0.0, Lz/2] ] )
green1    = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), nothing, vgreen; contactMaterial = cmat),     [ [0.0, 0.0, -Lz/2],[ 0.0, 0.0, Lz/2] ] )
green2    = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), nothing, vgreen; contactMaterial = cmat),     [ [0.0, 0.0, -Lz/2]] )
gelb1     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vyellow; contactMaterial = cmat), [ [0.0, 0.0, -Lz/2],[ 0.0, 0.0, Lz/2] ] )
gelb2     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidEllipsoid(Lx,Ly,Lz), "Aluminium", vyellow; contactMaterial = cmat), [ [0.0, 0.0, -Lz/2] ] )
gelb3     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidSphere(Dx), "Aluminium", vyellow; contactMaterial = cmat),    [ [0.0, 0.0, -Lz/2] ] )
rotM1      = Modia3D.Part(Modia3D.Solid(Modia3D.SolidPipe(Dx,Lz), nothing, vred; contactMaterial = cmat),          [[0.0, 0.0, -Lz/2] , [ 0.0, 0.0,Lz/2]] )
rotM2      = Modia3D.Part(Modia3D.Solid(Modia3D.SolidCylinder(Dx, Lz), nothing, vred; contactMaterial = cmat),     [[0.0, 0.0, -Lz/2], [ 0.0, 0.0,Lz/2]])
neuBlau1     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vblue;contactMaterial = cmat),    [ [0.0, 0.0, -Lz/2],[ 0.0, 0.0, Lz/2] ] )
neuBlau2     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vblue; contactMaterial = cmat),   [ [0.0, 0.0, -Lz/2], [ 0.0, 0.0, Lz/2] ] )

visuElem     = Modia3D.Part(cyl,   [ [0.0, 0.0, -Lz/2], [ 0.0, 0.0, Lz/2] ] )

Modia3D.connect(world_f1, blau1.frames[1]; R=ModiaMath.rot2(160u"°"))
Modia3D.connect(world_f1, blau2.frames[1]; R=ModiaMath.rot2(45u"°"))

rev2 = Modia3D.Revolute(blau2.frames[2], rosa1.frames[1]; phi_start =-pi/2)
rev5 = Modia3D.Revolute(blau1.frames[2], schwarz1.frames[1]; phi_start =-pi/2)
rev3 = Modia3D.Revolute(rosa1.frames[2], green1.frames[1])

Modia3D.connect(neuBlau1.frames[2],neuBlau2.frames[1]; R=ModiaMath.rot1(130u"°"))
rev6 = Modia3D.Revolute(rosa1.frame0, neuBlau1.frames[1])

Modia3D.connect(green1.frame0, green2.frames[1]; R=ModiaMath.rot2(45u"°"))
rev4 = Modia3D.Revolute(green1.frames[2], gelb1.frames[1])

Modia3D.connect(gelb1.frames[2],gelb2.frames[1]; R=ModiaMath.rot2(130u"°"))
Modia3D.connect(gelb1.frames[2], gelb3.frames[1]; R=ModiaMath.rot2(90u"°"))

#Modia3D.connect(world, rotM1.frames[1]; fixed=false, R=ModiaMath.rot3(90u"°"))
Modia3D.connect(blau2.frame0, rotM1.frames[1]; fixed=false, R=ModiaMath.rot2(90u"°"))
Modia3D.connect(rotM1.frames[2], rotM2.frames[1]) #; R=ModiaMath.rot2(90u"°"))

Modia3D.connect(gelb1.frames[2], visuElem.frames[2])
end

@assembly CollisionTest() begin
  world = Modia3D.Object3D(Modia3D.CoordinateSystem(1.0))

  collisionSolids = CollisionSolids(world)
end

collisionTest = CollisionTest(sceneOptions=Modia3D.SceneOptions(visualizeFrames=true, defaultFrameLength=0.3,nz_max = 80))

Modia3D.initAnalysis!(collisionTest)

tStart=0.0
tEnd  =1.0

for time = range(tStart, stop=tEnd, length=101)
  s = Modia3D.linearMovement(2*Lx, tStart, tEnd, time)
  delta_phi = Modia3D.linearMovement(pi/3, tStart, tEnd, time)

  Modia3D.setAngle!(collisionTest.collisionSolids.rev2,-pi/2 + delta_phi)
  Modia3D.setAngle!(collisionTest.collisionSolids.rev3, delta_phi)
  Modia3D.setAngle!(collisionTest.collisionSolids.rev4,-delta_phi)
  Modia3D.setAngle!(collisionTest.collisionSolids.rev5, delta_phi)
#  Modia3D.setAngle!(collisionTest.collisionSolids.rev6, -2*delta_phi)

  Modia3D.set_r!(collisionTest.collisionSolids.rotM1.frames[1], [0, 0.0, -s])

  Modia3D.updatePosition!(collisionTest)
#  Modia3D.selectContactPairsWithEvent!(collisionTest)

#  Modia3D.setComputationFlag(collisionTest)
  Modia3D.visualize!(collisionTest,time)

#=
  time += 0.01
  s = Modia3D.linearMovement(2*Lx, tStart, tEnd, time)
  Modia3D.set_r!(rotM1, [0, 0, s])
  Modia3D.updatePosition!(world)
  Modia3D.selectContactPairsWithEvent!(world)
  Modia3D.setComputationFlag(world)
  Modia3D.visualize!(world,time)
=#
  time += 0.01
end

Modia3D.closeAnalysis!(collisionTest)
println("... success of Test_Collision_StarSetting.jl!")
end


text:
const CanonicalRank = [:superkingdom,
                       :kingdom,
                       :phylum,
                       :class,
                       :order,
                       :family,
                       :genus,
                       :species,
                       :subspecies,
                       :strain
                       ]


text:
module ThetaPinchNoncanonical

    import ElectromagneticFields.ThetaPinch

    export charged_particle_3d_ode, charged_particle_3d_sode, charged_particle_3d_iode,
           hamiltonian, toroidal_momentum, ϑ

    const qᵢ = [2.5, 0.0, 0.0, 0.0, 0.2, 0.1]

    ThetaPinch.@code() # inject magnetic field code

    include("charged_particle_3d_noncanonical.jl")

end


text:
using Test
using Statistics
import ArchGDAL;
const AG = ArchGDAL;

@testset "test_ospy_examples.jl" begin
    """
    function to copy fields (not the data) from one layer to another
    parameters:
      fromLayer: layer object that contains the fields to copy
      toLayer: layer object to copy the fields into
    """
    function copyfields(fromlayer, tolayer)
        featuredefn = AG.layerdefn(fromlayer)
        for i in 0:(AG.nfield(featuredefn)-1)
            fd = AG.getfielddefn(featuredefn, i)
            if AG.gettype(fd) == OFTReal
                # to deal with errors like
                # ERROR: GDALError (Warning, code 1):
                # Value 18740682.1600000001 of field SHAPE_AREA of
                # feature 1 not successfully written. Possibly due
                # to too larger number with respect to field width
                fwidth = AG.getwidth(fd)
                if fwidth != 0
                    AG.setwidth!(fd, fwidth + 1)
                end
            end
            AG.createfield!(tolayer, fd)
        end
    end

    """
    function to copy attributes from one feature to another
    (this assumes the features have the same attribute fields!)
    parameters:
      fromFeature: feature object that contains the data to copy
      toFeature: feature object that the data is to be copied into
    """
    function copyattributes(fromfeature, tofeature)
        for i in 0:(AG.nfield(fromfeature)-1)
            if AG.isfieldset(fromfeature, i)
                try
                    AG.setfield!(tofeature, i, AG.getfield(fromfeature, i))
                catch
                    println(fromfeature)
                    println(tofeature)
                    println("$i: $(AG.getfield(fromfeature, i))")
                end
            end
        end
    end

    # function reproject(inFN, inEPSG, outEPSG)
    #     AG.fromEPSG(inEPSG) do inspatialref
    #     AG.fromEPSG(outEPSG) do outspatialref
    #     AG.createcoordtrans(inspatialref, outspatialref) do coordtrans
    #     AG.read(inFN) do inDS
    #     AG.create("", "MEMORY") do outDS
    #         inlayer = AG.getlayer(inDS, 0)
    #         outlayer = AG.createlayer(
    #                         name = "outlayer",
    #                         dataset = outDS,
    #                         geom = AG.getgeomtype(AG.layerdefn(inlayer)))
    #         copyfields(inlayer, outlayer)
    #         featuredefn = AG.layerdefn(outlayer)
    #         for infeature in inlayer
    #             geom = AG.getgeom(infeature)
    #             AG.createfeature(featuredefn) do outfeature
    #                 AG.setgeom!(outfeature, AG.transform!(geom, coordtrans))
    #                 copyattributes(infeature, outfeature)
    #                 AG.createfeature(outlayer, outfeature)
    #         end end
    #         println(outlayer)
    #     end
    #     end
    #     end
    #     println(AG.toWKT(AG.morphtoESRI!(outspatialref)))
    #     end
    #     end
    # end

    @testset "Homework 1" begin
        AG.read("ospy/data1/sites.shp") do input
            #reference: http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw1a.py
            for (i, feature) in enumerate(AG.getlayer(input, 0))
                id = AG.getfield(feature, 0)
                cover = AG.getfield(feature, 1)
                (x, y) = AG.getpoint(AG.getgeom(feature, 0), 0)
                @test id == i
                @test 4e5 <= x <= 5e5
                @test 4.5e6 <= y <= 5e6
                @test cover in
                      ("shrubs", "trees", "rocks", "grass", "bare", "water")
            end

            #reference: http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw1b.py
            # version 1
            AG.create(AG.getdriver("MEMORY")) do output
                inlayer = AG.getlayer(input, 0)
                outlayer = AG.createlayer(
                    name = "hw1b",
                    dataset = output,
                    geom = AG.wkbPoint,
                )
                inlayerdefn = AG.layerdefn(inlayer)
                AG.addfielddefn!(outlayer, AG.getfielddefn(inlayerdefn, 0))
                AG.addfielddefn!(outlayer, AG.getfielddefn(inlayerdefn, 1))
                for infeature in inlayer
                    id = AG.getfield(infeature, 0)
                    @test AG.asint64(infeature, 0) == id
                    cover = AG.getfield(infeature, 1)
                    if cover == "trees"
                        AG.createfeature(outlayer) do outfeature
                            AG.setgeom!(outfeature, AG.getgeom(infeature))
                            AG.setfield!(outfeature, 0, id)
                            return AG.setfield!(outfeature, 1, cover)
                        end
                    end
                end
                @test sprint(print, output) == """
                GDAL Dataset (Driver: Memory/Memory)
                File(s):

                Number of feature layers: 1
                  Layer 0: hw1b (wkbPoint)
                """
            end

            # version 2
            AG.create(AG.getdriver("MEMORY")) do output
                AG.executesql(
                    input,
                    """SELECT id, cover FROM sites
                       WHERE cover = 'trees' """,
                ) do results
                    @test sprint(print, results) == """
                    Layer: sites
                      Geometry 0 (): [wkbPoint], POINT (449959.840851...), ...
                         Field 0 (id): [OFTInteger], 2, 6, 9, 14, 19, 20, 22, 26, 34, 36, 41
                         Field 1 (cover): [OFTString], trees, trees, trees, trees, trees, trees, ...
                    """
                    return AG.copy(results, name = "hw1b", dataset = output)
                end
                @test sprint(print, output) == """
                GDAL Dataset (Driver: Memory/Memory)
                File(s):

                Number of feature layers: 1
                  Layer 0: hw1b (wkbPoint)
                """
            end
        end
    end

    @testset "Homework 2" begin
        # http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw2a.py
        open("ospy/data2/ut_counties.txt", "r") do file
            AG.create(AG.getdriver("MEMORY")) do output
                layer = AG.createlayer(
                    name = "hw2a",
                    dataset = output,
                    geom = AG.wkbPolygon,
                )
                @test sprint(print, layer) == """
                Layer: hw2a
                  Geometry 0 (): [wkbPolygon]
                """
                AG.createfielddefn("name", AG.OFTString) do fielddefn
                    AG.setwidth!(fielddefn, 30)
                    return AG.addfielddefn!(layer, fielddefn)
                end
                @test sprint(print, layer) == """
                Layer: hw2a
                  Geometry 0 (): [wkbPolygon]
                     Field 0 (name): [OFTString]
                """
                for line in readlines(file)
                    (name, coords) = split(line, ":")
                    coordlist = split(coords, ",")
                    AG.createfeature(layer) do feature
                        AG.setfield!(feature, 0, name)
                        AG.createpolygon() do poly
                            ring = AG.createlinearring()
                            for xy in map(split, coordlist)
                                AG.addpoint!(
                                    ring,
                                    parse(Float64, xy[1]),
                                    parse(Float64, xy[2]),
                                )
                            end
                            AG.addgeom!(poly, ring)
                            return AG.setgeom!(feature, poly)
                        end
                    end
                end
                @test sprint(print, layer) == """
                Layer: hw2a
                  Geometry 0 (): [wkbPolygon], POLYGON ((-111.50278...), ...
                     Field 0 (name): [OFTString], Cache, Box Elder, Rich, Weber, Morgan, ...
                """

                # input = output
                # # http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw2b.py
                # AG.fromEPSG(4269) do inspatialref
                # AG.fromEPSG(26912) do outspatialref
                # AG.createcoordtrans(inspatialref, outspatialref) do coordtrans
                # AG.create("", "MEMORY") do output
                #     inlayer = AG.getlayer(input, 0)
                #     outlayer = AG.createlayer(
                #         name = "hw2b",
                #         dataset = output,
                #         geom = AG.wkbPolygon
                #     )
                #     infeaturedefn = AG.layerdefn(inlayer)
                #     nameindex = AG.findfieldindex(infeaturedefn, "name")
                #     fielddefn = AG.getfielddefn(infeaturedefn, nameindex)
                #     AG.createfield!(outlayer, fielddefn)
                #     for infeature in inlayer
                #         AG.createfeature(outlayer) do outfeature
                #             geom = AG.getgeom(infeature)
                #             AG.setgeom!(outfeature, AG.transform!(geom, coordtrans))
                #             AG.setfield!(outfeature,0,AG.getfield(infeature, nameindex))
                #             println(outfeature)
                #     end end
                #     println(layer)
                # end
                # end
                # println(AG.toWKT(AG.morphtoESRI!(outspatialref)))
                # end
                # end
            end
        end
    end

    @testset "Homework 3" begin
        #reference: http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw3a.py
        AG.read("ospy/data1/sites.shp") do sitesDS
            AG.read("ospy/data3/cache_towns.shp") do townsDS
                siteslayer = AG.getlayer(sitesDS, 0)
                townslayer = AG.getlayer(townsDS, 0)
                AG.setattributefilter!(townslayer, "NAME = 'Nibley'")
                AG.getfeature(townslayer, 0) do nibleyFeature
                    AG.buffer(AG.getgeom(nibleyFeature), 1500) do bufferGeom
                        AG.setspatialfilter!(siteslayer, bufferGeom)
                        @test [AG.getfield(f, "ID") for f in siteslayer] == [26]
                    end
                end
            end
        end

        #reference: http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw3b.py
        # commented out until https://github.com/visr/GDAL.jl/issues/30 is resolved
        # for inFN in readdir("./ospy/data3/")
        #     if endswith(inFN, ".shp")
        #         reproject("./ospy/data3/$(inFN)", 26912, 4269)
        # end end
    end

    AG.read("ospy/data4/aster.img") do ds
        #reference: http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw4a.py
        @testset "Homework 4a" begin
            AG.read("ospy/data1/sites.shp") do shp
                shplayer = AG.getlayer(shp, 0)
                id = AG.findfieldindex(AG.layerdefn(shplayer), "ID")

                transform = AG.getgeotransform(ds)
                xOrigin = transform[1]
                yOrigin = transform[4]
                pixelWidth = transform[2]
                pixelHeight = transform[6]

                results = fill(0, 42, 3)
                for (i, feature) in enumerate(shplayer)
                    geom = AG.getgeom(feature)
                    x = AG.getx(geom, 0)
                    y = AG.gety(geom, 0)
                    # compute pixel offset
                    xOffset = round(Int, (x - xOrigin) / pixelWidth)
                    yOffset = round(Int, (y - yOrigin) / pixelHeight)
                    # create a string to print out
                    @test AG.getfield(feature, id) == i
                    for j in 1:AG.nraster(ds)
                        data = AG.read(ds, j, xOffset, yOffset, 1, 1)
                        results[i, j] = data[1, 1]
                    end
                end
                @test maximum(results) == 100
                @test minimum(results) == 0
                @test mean(results) ≈ 64.98412698412699
                @test std(results) ≈ 22.327734905980627
            end
        end

        #reference: http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw4b.py
        @testset "Homework 4b" begin
            @testset "version 1" begin
                count = 0
                total = 0
                data = AG.read(ds, 1)
                for (cols, rows) in AG.windows(AG.getband(ds, 1))
                    window = data[cols, rows]
                    count = count + sum(window .> 0)
                    total = total + sum(window)
                end
                @test total / count ≈ 76.33891347095299
                @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
            end

            @testset "version 2" begin
                band = AG.getband(ds, 1)
                count = 0
                total = 0
                buffer = Array{AG.pixeltype(band)}(undef, AG.blocksize(band)...)
                for (cols, rows) in AG.windows(band)
                    AG.rasterio!(band, buffer, rows, cols)
                    data = buffer[1:length(cols), 1:length(rows)]
                    count += sum(data .> 0)
                    total += sum(data)
                end
                @test total / count ≈ 76.33891347095299
                @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
            end

            @testset "version 3" begin
                count = 0
                total = 0
                # BufferIterator uses a single buffer, so this loop cannot be parallelized
                for data in AG.bufferwindows(AG.getband(ds, 1))
                    count += sum(data .> 0)
                    total += sum(data)
                end
                @test total / count ≈ 76.33891347095299
                @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
            end
        end

        #reference: http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw5a.py
        @testset "Homework 5a" begin
            @time begin
                rows = AG.height(ds)
                cols = AG.width(ds)
                bands = AG.nraster(ds)

                # get the band and block sizes
                inband2 = AG.getband(ds, 2)
                inband3 = AG.getband(ds, 3)
                (xbsize, ybsize) = AG.blocksize(inband2)

                buffer2 = Array{Float32}(undef, ybsize, xbsize)
                buffer3 = Array{Float32}(undef, ybsize, xbsize)
                ndvi = Array{Float32}(undef, ybsize, xbsize)
                AG.create(
                    AG.getdriver("MEM"),
                    width = cols,
                    height = rows,
                    nbands = 1,
                    dtype = Float32,
                ) do outDS
                    for ((i, j), (nrows, ncols)) in AG.blocks(inband2)
                        AG.rasterio!(inband2, buffer2, j, i, ncols, nrows)
                        AG.rasterio!(inband3, buffer3, j, i, ncols, nrows)
                        data2 = buffer2[1:nrows, 1:ncols]
                        data3 = buffer3[1:nrows, 1:ncols]
                        for row in 1:nrows, col in 1:ncols
                            denominator = data2[row, col] + data3[row, col]
                            if denominator > 0
                                numerator = data3[row, col] - data2[row, col]
                                ndvi[row, col] = numerator / denominator
                            else
                                ndvi[row, col] = -99
                            end
                        end
                        # write the data
                        AG.write!(outDS, ndvi, 1, j, i, ncols, nrows)
                    end
                    @test sprint(print, outDS) == """
                    GDAL Dataset (Driver: MEM/In Memory Raster)
                    File(s):

                    Dataset (width x height): 5665 x 5033 (pixels)
                    Number of raster bands: 1
                      [GA_Update] Band 1 (Undefined): 5665 x 5033 (Float32)
                    """
                    # flush data to disk, set the NoData value and calculate stats
                    outband = AG.getband(outDS, 1)
                    @test sprint(print, outband) == """
                    [GA_Update] Band 1 (Undefined): 5665 x 5033 (Float32)
                        blocksize: 5665×1, nodata: nothing, units: 1.0px + 0.0
                        overviews: """
                    AG.setnodatavalue!(outband, -99)
                    # georeference the image and set the projection
                    AG.setgeotransform!(outDS, AG.getgeotransform(ds))
                    return AG.setproj!(outDS, AG.getproj(ds))

                    # build pyramids
                    # gdal.SetConfigOption('HFA_USE_RRD', 'YES')
                    # AG.buildoverviews!(outDS,
                    #                   Cint[2,4,8,16,32,64,128], # overview list
                    #                   # bandlist (omit to include all bands)
                    #                   resampling="NEAREST")     # resampling method
                end
            end
        end
    end

    #reference: http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw5b.py
    @testset "Homework 5b" begin
        AG.read("ospy/data5/doq1.img") do ds1
            AG.read("ospy/data5/doq2.img") do ds2
                # read in doq1 and get info about it
                band1 = AG.getband(ds1, 1)
                rows1 = AG.height(ds1)
                cols1 = AG.width(ds1)

                # get the corner coordinates for doq1
                transform1 = AG.getgeotransform(ds1)
                minX1 = transform1[1]
                maxY1 = transform1[4]
                pixelWidth1 = transform1[2]
                pixelHeight1 = transform1[6]
                maxX1 = minX1 + (cols1 * pixelWidth1)
                minY1 = maxY1 + (rows1 * pixelHeight1)

                # read in doq2 and get info about it
                band2 = AG.getband(ds2, 1)
                rows2 = AG.height(ds2)
                cols2 = AG.width(ds2)

                # get the corner coordinates for doq1
                transform2 = AG.getgeotransform(ds2)
                minX2 = transform1[1]
                maxY2 = transform1[4]
                pixelWidth2 = transform1[2]
                pixelHeight2 = transform1[6]
                maxX2 = minX2 + (cols2 * pixelWidth2)
                minY2 = maxY2 + (rows2 * pixelHeight2)

                # get the corner coordinates for the output
                minX = min(minX1, minX2)
                maxX = max(maxX1, maxX2)
                minY = min(minY1, minY2)
                maxY = max(maxY1, maxY2)

                # get the number of rows and columns for the output
                cols = round(Int, (maxX - minX) / pixelWidth1)
                rows = round(Int, (maxY - minY) / abs(pixelHeight1))

                # compute the origin (upper left) offset for doq1
                xOffset1 = round(Int, (minX1 - minX) / pixelWidth1)
                yOffset1 = round(Int, (maxY1 - maxY) / pixelHeight1)

                # compute the origin (upper left) offset for doq2
                xOffset2 = round(Int, (minX2 - minX) / pixelWidth1)
                yOffset2 = round(Int, (maxY2 - maxY) / pixelHeight1)

                dtype = AG.pixeltype(band1)
                data1 = Array{dtype}(undef, rows, cols)
                data2 = Array{dtype}(undef, rows, cols)
                # create the output image
                AG.create(
                    AG.getdriver("MEM"),
                    width = cols,
                    height = rows,
                    nbands = 1,
                    dtype = AG.pixeltype(band1),
                ) do dsout
                    # read in doq1 and write it to the output
                    AG.rasterio!(band1, data1, 0, 0, cols1, rows1)
                    AG.write!(dsout, data1, 1, xOffset1, yOffset1, cols, rows)

                    # read in doq2 and write it to the output
                    AG.rasterio!(band2, data2, 0, 0, cols2, rows2)
                    AG.write!(dsout, data2, 1, xOffset2, yOffset2, cols, rows)

                    @test sprint(print, dsout) == """
                    GDAL Dataset (Driver: MEM/In Memory Raster)
                    File(s):

                    Dataset (width x height): 4500 x 3000 (pixels)
                    Number of raster bands: 1
                      [GA_Update] Band 1 (Undefined): 4500 x 3000 (UInt8)
                    """
                    # compute statistics for the output
                    bandout = AG.getband(dsout, 1)
                    @test sprint(print, bandout) == """
                    [GA_Update] Band 1 (Undefined): 4500 x 3000 (UInt8)
                        blocksize: 4500×1, nodata: nothing, units: 1.0px + 0.0
                        overviews: """

                    # set the geotransform and projection on the output
                    geotransform = [minX, pixelWidth1, 0, maxY, 0, pixelHeight1]
                    AG.setgeotransform!(dsout, geotransform)
                    return AG.setproj!(dsout, AG.getproj(ds1))

                    # build pyramids for the output
                    # gdal.SetConfigOption('HFA_USE_RRD', 'YES')
                    # buildoverviews not supported for in-memory rasters
                    # AG.buildoverviews!(dsout,
                    #                    Cint[2,4,8,16],       # overview list
                    #                   # bandlist (omit to include all bands)
                    #                    resampling="NEAREST") # resampling method
                end
            end
        end
    end
end


text:
export demean!

using Statistics: mean

function demean!(A::AbstractArray)
    μ = mean(A)
    A .-= μ
    return nothing
end


text:
# See **[`ProjectiveTransforms`](#).**

"""
    ProjectiveTransforms(size, [augmentations])

Pipeline step that resizes images and keypoints to `size`.

In context [`Training`](#), applies `augmentations`.
"""
@with_kw_noshow struct ProjectiveTransforms <: PipelineStep
    traintfm
    validtfm
    inferencetfm
end

function ProjectiveTransforms(
        size;
        augmentations = Identity(),
        inferencefactor = 1,
        buffered = true)
    tfms = (
        augmentations |> RandomResizeCrop(size),
        CenterResizeCrop(size),
        ResizePadDivisible(size, inferencefactor),
    )

    if buffered
        tfms = (
            BufferedThreadsafe(tfms[1]),
            BufferedThreadsafe(tfms[2]),
            # Inference transform is not buffered since it can have
            # varying sizes
            tfms[3]
        )
    end

    return ProjectiveTransforms(tfms...)
end

function Base.show(io::IO, spatial::ProjectiveTransforms)
    outsize = _parenttfm(spatial.validtfm).transforms[1].crop.size
    print(io, "ProjectiveTransforms($(outsize))")
end

function run(spatial::ProjectiveTransforms, context, datas::Tuple)
    items = makespatialitems(datas)
    tfm = _gettfm(spatial, context)
    tdatas = itemdata.(DataAugmentation.apply(tfm, items))
    return tdatas
end

function run!(bufs, spatial::ProjectiveTransforms, context, datas::Tuple)
    items = makespatialitems(datas)
    tfm = _gettfm(spatial, context)
    tdatas = DataAugmentation.apply(tfm, items) |> itemdata
    _copyrec!(bufs, tdatas)
    return bufs
end

run(spatial::ProjectiveTransforms, context, data) = run(spatial, context, (data,)) |> only
run!(buf, spatial::ProjectiveTransforms, context, data) = run!((buf,), spatial, context, (data,)) |> only

## Utils

_gettfm(spatial::ProjectiveTransforms, context::Training) = spatial.traintfm
_gettfm(spatial::ProjectiveTransforms, context::Validation) = spatial.validtfm
_gettfm(spatial::ProjectiveTransforms, context::Inference) = spatial.inferencetfm

function makespatialitems(datas::Tuple)
    if datas[begin] isa Item
        return makespatialitems(datas, getbounds(datas[begin]))
    else
        return makespatialitems(datas, makebounds(size(datas[begin])))
    end
end
function makespatialitems(datas::Tuple, bounds)
    return Tuple(makeitem(data, bounds) for data in datas)
end

"""
    makeitem(data, args...)

Tries to assign a `DataAugmentation.Item` from `data` based on its type.
`args` are passed to the chosen  `Item` constructor.

- `AbstractMatrix{<:Colorant}` -> `Image`
- `Vector{<:Union{Nothing, SVector}}` -> `Keypoints`
"""
makeitem(data, args...) = itemtype(data)(data, args...)
makeitem(item::Item, args...) = item
makeitem(datas::Tuple, args...) = Tuple(makeitem(data, args...) for data in datas)

itemtype(data::AbstractMatrix{<:Colorant}) = DataAugmentation.Image
itemtype(data::Vector{<:Union{Nothing, SVector}}) = DataAugmentation.Keypoints


text:
using StatsBase, StaticArrays, BitBasis, Random
export measure, measure!, measure_remove!, measure_collapseto!, select, select!

function _measure(rng::AbstractRNG, pl::AbstractVector, nshots::Int)
    N = log2i(length(pl))
    sample(rng, basis(BitStr64{N}), Weights(pl), nshots)
end

function _measure(rng::AbstractRNG, pl::AbstractMatrix, nshots::Int)
    B = size(pl, 2)
    res = Matrix{BitStr64{log2i(length(pl))}}(undef, nshots, B)
    for ib in 1:B
        @inbounds res[:, ib] = _measure(rng, view(pl, :, ib), nshots)
    end
    return res
end

YaoBase.measure(
    rng::AbstractRNG,
    ::ComputationalBasis,
    reg::ArrayReg{1},
    ::AllLocs;
    nshots::Int = 1,
) = _measure(rng, reg |> probs, nshots)

function YaoBase.measure(
    rng::AbstractRNG,
    ::ComputationalBasis,
    reg::ArrayReg{B},
    ::AllLocs;
    nshots::Int = 1,
) where {B}
    pl = dropdims(sum(reg |> rank3 .|> abs2, dims = 2), dims = 2)
    return _measure(rng, pl, nshots)
end

function YaoBase.measure_remove!(
    rng::AbstractRNG,
    ::ComputationalBasis,
    reg::ArrayReg{B},
    ::AllLocs,
) where {B}
    state = reg |> rank3
    nstate = similar(reg.state, 1 << nremain(reg), B)
    pl = dropdims(sum(state .|> abs2, dims = 2), dims = 2)
    res = Vector{BitStr64{nactive(reg)}}(undef, B)
    @inbounds for ib in 1:B
        ires = _measure(rng, view(pl, :, ib), 1)[]
        # notice ires is `BitStr` type, can be use as indices directly.
        nstate[:, ib] = view(state, Int64(ires) + 1, :, ib) ./ sqrt(pl[Int64(ires)+1, ib])
        res[ib] = ires
    end
    reg.state = reshape(nstate, 1, :)
    return res
end

function YaoBase.measure!(
    rng::AbstractRNG,
    ::ComputationalBasis,
    reg::ArrayReg{B},
    ::AllLocs,
) where {B}
    state = reg |> rank3
    nstate = zero(state)
    res = measure_remove!(rng, reg)
    _nstate = reshape(reg.state, :, B)
    for ib in 1:B
        @inbounds nstate[Int64(res[ib])+1, :, ib] .= view(_nstate, :, ib)
    end
    reg.state = reshape(nstate, size(state, 1), :)
    return res
end

function YaoBase.measure_collapseto!(
    rng::AbstractRNG,
    ::ComputationalBasis,
    reg::ArrayReg{B},
    ::AllLocs;
    config::Integer = 0,
) where {B}
    state = rank3(reg)
    M, N, B1 = size(state)
    nstate = zero(state)
    res = measure_remove!(rng, reg)
    nstate[Int(config)+1, :, :] = reshape(reg.state, :, B)
    reg.state = reshape(nstate, M, N * B)
    return res
end

import YaoBase: select, select!
select(r::ArrayReg{B}, bits::AbstractVector{T}) where {B,T<:Integer} =
    ArrayReg{B}(r.state[Int64.(bits).+1, :])
select(r::ArrayReg{B}, bit::Integer) where {B} = select(r, [bit])

function select!(r::ArrayReg, bits::AbstractVector{T}) where {T<:Integer}
    r.state = r.state[Int64.(bits).+1, :]
    return r
end

select!(r::ArrayReg, bit::Integer) = select!(r, [bit])


text:
Random.seed!(0)
@testset "SparseMatrixCSC: Constructors B=$B T=$T" for B in [1, 2, 4, 8],
                                                       T in [Int32, Float32,
                                                        Complex{Float32},
                                                        SMatrix{2,2,Float32,
                                                                2 * 2}]

    leni = 9 * B
    lenj = 11 * B
    N = bspzeros(Val(B), T, leni, lenj)
    @test size(N) == (leni, lenj)
    @test eltype(N) == T

    A = bsprand(Val(B), T, leni, lenj, 0.3)
    @test size(A) == (leni, lenj)
    @test eltype(A) == T
    if !(T <: Integer)
        C = bsprandn(Val(B), T, leni, lenj, 0.3)
        @test size(C) == (leni, lenj)
        @test eltype(C) == T
    end

    Af = Matrix(A)
    @test eltype(Af) == T
    @test size(Af) == size(A)
    A′ = blocksparse(Val(B), Af)
    @test size(A′) == size(Af)
    @test eltype(A′) == T
    @test A′ == A

    As = SparseMatrixCSC(A)
    @test eltype(As) == T
    @test size(As) == (leni, lenj)
    A′ = blocksparse(Val(B), As)
    @test eltype(A′) == T
    @test A′ == A

    zf = rand(T, leni, lenj)::Matrix{T}
    z = blocksparse(Val(B), zf)
    @test size(z) == size(zf)
    @test eltype(z) == T
    zf′ = Array(z)
    @test eltype(zf′) == T
    @test zf′ == zf

    zs = sprand(T, leni, lenj, 0.3)::SparseMatrixCSC{T}
    z = blocksparse(Val(B), zs)
    @test size(z) == size(zs)
    @test Matrix(zs) == Matrix(z)
    @test eltype(z) == T
    zs′ = SparseMatrixCSC(z)
    @test eltype(zs′) == T
    @test zs′ == zs
end

Random.seed!(0)
@testset "SparseMatrixCSC: Collection B=$B T=$T" for B in [1, 2, 4, 8],
                                                     T in [Int32, Float32,
                                                      Complex{Float32},
                                                      SMatrix{2,2,Float32,
                                                              2 * 2}]

    leni = 9 * B
    lenj = 11 * B
    A = bsprand(Val(B), T, leni, lenj, 0.3)
    As = SparseMatrixCSC(A)
    # @test A == As
    for i in eachindex(A)
        @test A[i] == As[i]
    end
end

Random.seed!(0)
@testset "SparseMatrixCSC: Vector space B=$B T=$T" for B in [1, 2, 4, 8],
                                                       T in [Int32, Float32,
                                                        Complex{Float32},
                                                        SMatrix{2,2,Float32,
                                                                2 * 2}]

    leni = 9 * B
    lenj = 11 * B

    N = bspzeros(Val(B), T, leni, lenj)
    A = bsprand(Val(B), T, leni, lenj, 0.3)
    C = bsprand(Val(B), T, leni, lenj, 0.3)
    D = bsprand(Val(B), T, leni, lenj, 0.3)
    a = rand(eltype(T))
    b = rand(eltype(T))

    @test A == A
    @test C ≠ A

    @test +A == A
    @test -(-A) == A
    @test N + A == A
    @test A + N == A
    @test N - A == -A
    @test A - N == A
    @test A + C == C + A
    @test (A + C) + D == A + (C + D)

    @test zero(a) * A == N
    @test one(a) * A == A
    @test (-one(a)) * A == -A
    @test a * A == A * a
    @test a * (b * A) ≈ (a * b) * A

    if T <: Union{AbstractFloat,Complex{<:AbstractFloat}}
        @test inv(a) * (a * A) ≈ A
        @test A / a ≈ A * inv(a)
        @test a \ A ≈ inv(a) * A
    end

    @test (a + b) * A ≈ a * A + b * A
    @test a * (A + C) ≈ a * A + a * C
end

Random.seed!(0)
@testset "SparseMatrixCSC: Category B=$B T=$T" for B in [1, 2, 4, 8],
                                                   T in [Int32, Float32,
                                                    Complex{Float32},
                                                    SMatrix{2,2,Float32,2 * 2}]

    len = 9 * B

    N = bspzeros(Val(B), T, len, len)
    E = one(N)
    A = bsprand(Val(B), T, len, len, 0.3)
    C = bsprand(Val(B), T, len, len, 0.3)
    D = bsprand(Val(B), T, len, len, 0.3)
    a = rand(eltype(T))
    b = rand(eltype(T))

    @test N * A == N
    @test A * N == N

    @test E * A == A
    @test A * E == A

    @test (A * C) * D ≈ A * (C * D)

    @test a * (A * C) ≈ (a * A) * C

    @test (A + C) * D ≈ A * D + C * D
    @test A * (C + D) ≈ A * C + A * D
end

Random.seed!(0)
@testset "SparseMatrixCSC: Operator B=$B T=$T" for B in [1, 2, 4, 8],
                                                   T in [Int32, Float32,
                                                    Complex{Float32},
                                                    SMatrix{2,2,Float32,2 * 2}]

    len = 9 * B

    n = bspzeros(Val(B), T, len)
    x = bsprand(Val(B), T, len, 0.3)
    y = bsprand(Val(B), T, len, 0.3)
    N = bspzeros(Val(B), T, len, len)
    E = one(N)
    A = bsprand(Val(B), T, len, len, 0.3)
    C = bsprand(Val(B), T, len, len, 0.3)
    a = rand(eltype(T))
    b = rand(eltype(T))

    @test N * x == n
    @test E * x == x

    @test A * n == n

    @test A * (x + y) ≈ A * x + A * y
    @test (A + C) * x ≈ A * x + C * x

    @test (a * A) * x ≈ a * (A * x)
    @test (A * a) * x ≈ A * (a * x)

    if !(T <: Integer)          # doesn't make sense otherwise
        if T <: Union{AbstractFloat,Complex} # UMFPACK only supports some types
            # We can't do LU decompositions! We lost.
            # r = A \ x
            # @test A * r ≈ x
        end
    end
end


text:
include("iPredict.jl/src/iPredict.jl")
using .iPredict
using DataFrames
using Plots, Plots.Measures
using JLD2
using DelimitedFiles
using StatsBase

function make_hist(path, titleArg)
    Arr = readdlm(path)

    if (titleArg == "3 Isoforms" &&
        path == "figures/figure2_data/H9_96_scnorm_counts.txt_3_isoforms_array.txt")
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 10),
        ylabel = "Probability Density",
        normalize=:pdf,
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)

    elseif titleArg == "Real" #|| titleArg == "alpha=0.45, beta=1.45"
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 10),
        normalize=:pdf,
        xlabel = "Mean No. Isoforms\nper Gene per Cell",
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)

    else
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 10),
        normalize=:pdf,
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)
    end
    vline!([mean(Arr)], color = :black, linewidth = 8)
    return p1
end

H9_24_1 = make_hist("figures/figure2_data/H9_24_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
H9_24_2 = make_hist("figures/figure2_data/H9_24_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
H9_24_3 = make_hist("figures/figure2_data/H9_24_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
H9_24_4 = make_hist("figures/figure2_data/H9_24_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
H9_24_real = make_hist("figures/figure2_data/H9_24_scnorm_counts.txt_real_isoforms_array.txt", "Real")

H9_96_1 = make_hist("figures/figure2_data/H9_96_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
H9_96_2 = make_hist("figures/figure2_data/H9_96_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
H9_96_3 = make_hist("figures/figure2_data/H9_96_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
H9_96_4 = make_hist("figures/figure2_data/H9_96_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
H9_96_real = make_hist("figures/figure2_data/H9_96_scnorm_counts.txt_real_isoforms_array.txt", "Real")

pA = plot(H9_96_1,H9_24_1,
H9_96_2,H9_24_2,
H9_96_3,H9_24_3,
H9_96_4,H9_24_4,
H9_96_real,H9_24_real,
layout = (5,2), legend = false, size = (500*8, 600*8),
margin = 45mm)

savefig(pA, "figures/supplementary_figs/figure1A.png")

################################################################################
# Panel B
################################################################################

include("iPredict.jl/src/TableProcessing.jl")
include("iPredict.jl/src/FindKm.jl")
include("iPredict.jl/src/Simulate.jl")
include("iPredict.jl/src/Plots.jl")
include("iPredict.jl/src/pValues.jl")
import Cairo
import Fontconfig
import DataFrames
using Distributions
using DelimitedFiles

#Set probability(FP|N) and probability(FN|P)
pFN = 0.04
pFP = 0.01

#Function that performs global level predictions
function pDropoutInfo(countsMatrixPath::String, GeneIsoTablePath::String,
      output::String, filteringThreshold::Int64)

    #Read in counts matrix
    println("Data wrangling commencing!")
    println("Reading in counts matrix.") #Quick
    countsMatrix = ReadCountsMatrix(countsMatrixPath)

    #Read in gene-isoform table
    println("Reading in gene-isoform relationships")
    GeneIsoformRelationships = ReadGeneIsoformRelationships(GeneIsoTablePath)

    #Make sure gene-isoform table matches counts matrix
    if CheckGenesMatchIsoforms(countsMatrix, GeneIsoformRelationships) == false
        throw(ArgumentError("Your counts matrix doesn't match your list of genes
        and isoforms - there is something wrong with your input files."))
    end

    #Add gene information to isoforms table
    println("Adding Gene information to countsMatrix")
    countsMatrixWithGenes = AddGeneInfo(countsMatrix, GeneIsoformRelationships)

    #Clean data for finding Km
    println("Cleaning data for finding Km")
    MMdf = CleanDataForMM(countsMatrix)

    #Find Km
    println("Calculating Km....")
    MMParams = FindKm(MMdf)
    Km = MMParams[1]
    sigma = MMParams[2]

    println("Converting counts to dict")
    GenesWithNIsoforms = KeepGenesWithNIsoforms(countsMatrixWithGenes,
    filteringThreshold)
    workingDict = GenesWithNIsoforms[1]
    unexpr_dict = GenesWithNIsoforms[2]
    MaxNumIsos = GenesWithNIsoforms[3]

    #Get ranking and dropout info for each gene
    pDropoutArr::Array{Float64} = []
    for (key, value) in workingDict
        append!(pDropoutArr, FindPDropout(value, Km))
    end

    pDropoutDistributionPlot = pDropoutDistribution(pDropoutArr, output)
    params = fit(Beta, pDropoutArr)
    return (params.α, params.β, pDropoutDistributionPlot[1], pDropoutDistributionPlot[2])
end

H9_96 = pDropoutInfo("data/H9_96_scnorm_counts.txt", "data/gencode_human_gene_isoform_relationships.txt", "figures/pdrop_scnorm_H9_96.png", 4)
H9_24 = pDropoutInfo("data/H9_24_scnorm_counts.txt", "data/gencode_human_gene_isoform_relationships.txt", "figures/pdrop_scnorm_H9_24.png", 4)

pB = plot(H9_96[3], H9_24[3],
       H9_96[4], H9_24[4],
        size = (500 * 8, 240 * 8), layout = (2,2), legend = false, margin = 35mm)

        @show H9_96[1]
        @show H9_96[2]
        @show H9_24[1]
        @show H9_24[2]

savefig(pB, "figures/supplementary_figs/figure1B.png")

################################################################################
# Panel C
################################################################################

v0 = pdf(Beta(0.45, 1.45), 0:0.01:1)
v1 = pdf(Beta(0.7, 1.2), 0:0.01:1)
v2 = pdf(Beta(0.95, 0.95), 0:0.01:1)
v3 = pdf(Beta(1.2, 0.7), 0:0.01:1)
v4 = pdf(Beta(1.45, 0.45), 0:0.01:1)
pC = plot(0:0.01:1, [v0,v1,v2,v3,v4], label = ["alpha=0.45, beta=1.45",
       "alpha=0.7, beta=1.2",
       "alpha=0.95, beta=0.95",
       "alpha=1.2, beta=0.7",
       "alpha=1.45, beta=0.45"], size = (250 * 8, 240 * 8), axis = (font(60)), margin = 35mm,
       titlefontsize=80, legendfontsize = 60, linestyle=:solid, linealpha=0.5, linewidth=2*8, ylabel = "Normalised Frequency", xlabel = "Probability", title = "Beta Distributions")
savefig(pC, "figures/supplementary_figs/figure1C.png")

################################################################################
# Panel D
################################################################################

function make_hist_MM(path, titleArg)
    Arr = readdlm(path)

    if (titleArg == "3 Isoforms" &&
        path == "figures/figure2_data/H9_96_scnorm_counts.txt_3_isoforms_array.txt")
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 1.2),
        ylims = (0, 5),
        ylabel = "Probability Density",
        normalize=:pdf,
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)

    elseif titleArg == "Real" #|| titleArg == "alpha=0.45, beta=1.45"
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 1.2),
        ylims = (0, 5),
        normalize=:pdf,
        xlabel = "Mean No. Isoforms\nper Gene per Cell",
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)

    else
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 1.2),
        ylims = (0, 5),
        normalize=:pdf,
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)
    end
    vline!([mean(Arr)], color = :black, linewidth = 8)
    return p1
end

v0 =  make_hist_MM("figures/figure2_data/0.45_1.45_H9_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=0.45, beta=1.45")
v1 =  make_hist_MM("figures/figure2_data/0.7_1.2_H9_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=0.7, beta=1.2")
v2 =  make_hist_MM("figures/figure2_data/0.95_0.95_H9_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=0.95, beta=0.95")
v3 =  make_hist_MM("figures/figure2_data/1.2_0.7_H9_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=1.2, beta=0.7")
v4 =  make_hist_MM("figures/figure2_data/1.45_0.45_H9_24_scnorm_counts.txt_1_isoforms_array.txt", "alpha=1.45, beta=0.45")

pD = plot(v4, v3, v2, v1, v0,
       layout = (5,1), legend = false, size = (250*8, 600*8),
       axis = (font(180)), titlefontsize=160, margin = 25mm)

savefig(pD, "figures/supplementary_figs/figure1D.png")

####################################
# Overlap plots
####################################

function make_overlap_hist(path, titleArg)

    Arr = readdlm(path)

    if titleArg == "4 Isoforms"

        p1 = histogram(Arr,
        title = titleArg,
        color = palette(:default)[2],
        xlims = (0, 1),
        ylims = (0, 10),
        ylabel = "Probability Density",
        xlabel = "% Overlap with Ground Truth",
        normalize=:pdf,
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)
        vline!([mean(Arr)], color = :black, linewidth = 8)
        return p1

    else

        p1 = histogram(Arr,
        title = titleArg,
        color = palette(:default)[2],
        xlims = (0, 1),
        ylims = (0, 10),
        ylabel = "Probability Density",
        normalize=:pdf,
        bins = 100,
        axis = (font(100)), titlefontsize=100, yticks = 0:5:10,
        margin = 25mm)
        vline!([mean(Arr)], color = :black, linewidth = 8)
        return p1
    end

end

H9_24_1 = make_overlap_hist("figures/figure2_data/H9_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
H9_24_2 = make_overlap_hist("figures/figure2_data/H9_24_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
H9_24_3 = make_overlap_hist("figures/figure2_data/H9_24_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
H9_24_4 = make_overlap_hist("figures/figure2_data/H9_24_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

H9_96_1 = make_overlap_hist("figures/figure2_data/H9_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
H9_96_2 = make_overlap_hist("figures/figure2_data/H9_96_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
H9_96_3 = make_overlap_hist("figures/figure2_data/H9_96_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
H9_96_4 = make_overlap_hist("figures/figure2_data/H9_96_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

pA = plot(H9_96_1,H9_24_1,
H9_96_2,H9_24_2,
H9_96_3,H9_24_3,
H9_96_4,H9_24_4,
layout = (4,2), legend = false, size = (400*8, 600*8))

savefig(pA, "figures/supplementary_figs/figure2_overlap.png")

########################
# Old 1D overlap hists
########################

v0 =  make_overlap_hist("figures/figure2_data/0.45_1.45_H9_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=0.45, beta=1.45")
v1 =  make_overlap_hist("figures/figure2_data/0.7_1.2_H9_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=0.7, beta=1.2")
v2 =  make_overlap_hist("figures/figure2_data/0.95_0.95_H9_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=0.95, beta=0.95")
v3 =  make_overlap_hist("figures/figure2_data/1.2_0.7_H9_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=1.2, beta=0.7")
v4 =  make_overlap_hist("figures/figure2_data/1.45_0.45_H9_24_scnorm_counts.txt_1_isoforms_overlap_array.txt", "alpha=1.45, beta=0.45")

pD = plot(v4, v3, v2, v1, v0,
       layout = (5,1), legend = false, size = (250*8, 600*8),
       axis = (font(180)), titlefontsize=160, margin = 30mm)

savefig(pD, "figures/supplementary_figs/figure2D_overlap.png")

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
right_url = "https://github.com/jpjones76/SeisIO-TestData/trunk/SVN_test"
wrong_url = "https://github.com/jpjones76/SeisIO-TestData/trunk/DOESNT_EXIST"
svn_dir = "SVN_test"
svn_file = joinpath(svn_dir, "test")

try
  get_svn(right_url, svn_dir)
  @test isfile(svn_file)
  str = readline(svn_file)
  @test str == "."
  rm(svn_file)
  rm(svn_dir)
  @test_throws ErrorException get_svn(wrong_url, "SVN_test")
catch err
  @warn(err)
end


text:
module BinaryBuilderBase

using Pkg, Pkg.Artifacts, Random, Libdl, InteractiveUtils
using Base.BinaryPlatforms
using Downloads
using JSON, OutputCollectors, Scratch

# Re-export useful stuff from Base.BinaryPlatforms:
export HostPlatform, platform_dlext, valid_dl_path, arch, libc,
       libgfortran_version, libstdcxx_version, cxxstring_abi,
       detect_libgfortran_version, detect_libstdcxx_version, detect_cxxstring_abi,
       call_abi, wordsize, triplet, select_platform, platforms_match,
       Platform, AnyPlatform

export AbstractSource, AbstractDependency, SetupSource, PatchSource,
    resolve_jlls, coerce_dependency, coerce_source, Runner,
    generate_compiler_wrappers!, preferred_runner, CompilerShard, UserNSRunner,
    DockerRunner, choose_shards, exeext, preferred_libgfortran_version,
    preferred_cxxstring_abi, gcc_version, available_gcc_builds, getversion,
    getpkg, aatriplet,
    nbits, proc_family, storage_dir, extract_kwargs, extract_fields,
    download_source, setup_workspace, setup_dependencies, update_registry,
    getname, cleanup_dependencies, compress_dir, prepare_for_deletion,
    run_interactive, sourcify, dependencify, with_logfile, get_concrete_platform

include("compat.jl")

include("ArchiveUtils.jl")
include("Sources.jl")
include("Dependencies.jl")
include("Prefix.jl")
include("Products.jl")
include("Platforms.jl")

include("Runner.jl")
include("BuildToolchains.jl")
include("Rootfs.jl")
include("squashfs_utils.jl")
include("UserNSRunner.jl")
include("DockerRunner.jl")

include("utils.jl")

# This is the location that all binary builder-related files are stored under.
# downloads, unpacked .tar.gz shards, mounted shards, ccache cache, etc....
function storage_dir(args::AbstractString...)
    global storage_cache
    dir = joinpath(storage_cache, args...)
    mkpath(dirname(dir))
    return dir
end
ccache_dir() = storage_dir("ccache")

"""
    enable_apple_file()

Return the path to file that, if exists, indicates that the user accepts to
download macOS SDK.  The file is automatically created when the package is
loaded if the environment variable `BINARYBUILDER_AUTOMATIC_APPLE` is set to
`"true"`.
"""
enable_apple_file() = storage_dir("enable_apple")

# These globals store important information such as where we're downloading
# the rootfs to, and where we're unpacking it.  These constants are initialized
# by `__init__()` to allow for environment variable overrides from the user.
storage_cache = ""
use_squashfs = false
allow_ecryptfs = false
use_ccache = false
bootstrap_list = Symbol[]

function get_bbb_version(dir=@__DIR__, uuid="7f725544-6523-48cd-82d1-3fa08ff4056e")
    # Get BinaryBuilder.jl's version and git sha
    version = Pkg.TOML.parsefile(joinpath(dir, "..", "Project.toml"))["version"]
    try
        # get the gitsha if we can
        repo = LibGit2.GitRepo(dirname(@__DIR__))
        gitsha = string(LibGit2.GitHash(LibGit2.GitCommit(repo, "HEAD")))
        return VersionNumber("$(version)-git-$(gitsha[1:10])")
    catch
        try
            # Settle for the treehash otherwise
            env = Pkg.Types.Context().env
            bb_uuid = Pkg.Types.UUID(uuid)
            treehash = bytes2hex(env.manifest[bb_uuid].tree_hash.bytes)
            return VersionNumber("$(version)-tree-$(treehash[1:10])")
        catch
            # Something went so wrong, we can't get any of that.
            return VersionNumber(version)
        end
    end
end

"""
    versioninfo()

Helper function to print out some debugging information
"""
function versioninfo(; name=@__MODULE__, version=get_bbb_version())
    @info("Julia versioninfo(): ")
    InteractiveUtils.versioninfo()

    @info("$(name).jl version: $(version)")

    @static if Sys.isunix()
        @info("Kernel version: $(readchomp(`uname -r`))")
    end

    # Dump if some important directories are encrypted:
    @static if Sys.islinux()
        print_enc(n, path) = begin
            is_encrypted, mountpoint = is_ecryptfs(path)
            if is_encrypted
                @info("$n is encrypted on mountpoint $mountpoint")
            else
                @info("$n is NOT encrypted on mountpoint $mountpoint")
            end
        end

        print_enc("pkg dir", dirname(@__FILE__))
        print_enc("storage dir", storage_dir())
    end

    # Dump any relevant environment variables:
    @info("Relevant environment variables:")
    env_var_suffixes = [
        "AUTOMATIC_APPLE",
        "USE_SQUASHFS",
        "STORAGE_DIR",
        "RUNNER",
        "ALLOW_ECRYPTFS",
        "USE_CCACHE",
    ]
    for e in env_var_suffixes
        envvar = "BINARYBUILDER_$(e)"
        if haskey(ENV, envvar)
            @info("  $(envvar): \"$(ENV[envvar])\"")
        end
    end

    # Print out the preferred runner stuff here:
    @info("Preferred runner: $(preferred_runner())")

    # Try to run 'echo julia' in Linux x86_64 environment
    @info("Trying to run `echo hello julia` within a Linux x86_64 environment...")

    runner = preferred_runner()(
        pwd();
        cwd="/workspace/",
        platform=Platform("x86_64", "linux"),
        verbose=true
    )
    run_interactive(runner, `/bin/bash -c "echo hello julia"`)

    # If we use ccache, dump the ccache stats
    if use_ccache
        @info("ccache stats:")
        runner = preferred_runner()(
            pwd();
            cwd="/workspace/",
            platform=Platform("x86_64", "linux"),
        )
        run_interactive(runner, `/usr/bin/ccache -s`)
    end
    return nothing
end

function __init__()
    global runner_override, use_squashfs, allow_ecryptfs
    global use_ccache, storage_cache

    # Allow the user to override the default value for `storage_dir`
    storage_cache = get(ENV, "BINARYBUILDER_STORAGE_DIR",
                        abspath(joinpath(@__DIR__, "..", "deps")))

    # If the user has signalled that they really want us to automatically
    # accept apple EULAs, do that.
    if get(ENV, "BINARYBUILDER_AUTOMATIC_APPLE", "") == "true"
        touch(enable_apple_file())
    end

    # If the user has overridden our runner selection algorithms, honor that
    runner_override = lowercase(get(ENV, "BINARYBUILDER_RUNNER", ""))
    if runner_override == "unprivileged"
        runner_override = "userns"
    end
    if !(runner_override in ["", "userns", "privileged", "docker"])
        @warn("Invalid runner value $runner_override, ignoring...")
        runner_override = ""
    end

    # If the user has asked for squashfs mounting instead of tarball mounting,
    # use that here.  Note that on Travis, we default to using squashfs, unless
    # BINARYBUILDER_USE_SQUASHFS is set to "false", which overrides this
    # default. If we are not on Travis, we default to using tarballs and not
    # squashfs images as using them requires `sudo` access.
    if get(ENV, "BINARYBUILDER_USE_SQUASHFS", "") == "false"
        use_squashfs = false
    elseif get(ENV, "BINARYBUILDER_USE_SQUASHFS", "") == "true"
        use_squashfs = true
    else
        # If it hasn't been specified, but we're on Travis, default to "on"
        if get(ENV, "TRAVIS", "") == "true"
            use_squashfs = true
        end

        # If it hasn't been specified but we're going to use the docker runner,
        # then set `use_squashfs` to `true` here.
        if preferred_runner() == DockerRunner
            # Conversely, if we're dock'ing it up, don't use it.
            use_squashfs = false
        elseif runner_override == "privileged"
            # If we're forcing a privileged runner, go ahead and default to squashfs
            use_squashfs = true
        end
    end

    # If the user has signified that they want to allow mounting of ecryptfs
    # paths, then let them do so at their own peril.
    if get(ENV, "BINARYBUILDER_ALLOW_ECRYPTFS", "") == "true"
        allow_ecryptfs = true
    end

    # If the user has enabled `ccache` support, use it!
    if get(ENV, "BINARYBUILDER_USE_CCACHE", "false") == "true"
        use_ccache = true
    end
end

end # module BinaryBuilderBase


text:
"""
    add_sub_heatmap!(f, m)

Internal function to add a heatmap subroutine to the GLE script.
This subroutine only generates one column (resp. one row) of boxes, not the
full grid (that is done by calling the routine repeatedly).
"""
function add_sub_heatmap!(f::Figure, hm::Heatmap, hashid::UInt)
    #=
    zij = 0
    cij = ""
    for i = 1 to ndata(ds)
        zij = datayvalue(ds, i)
        if zij <= 1 then
            cij = firstcolor
        else if zij <=2 then
            cij = secondcolor
        ...
        else
            cij = missingcolor
        end

        move and  draw box(zij, cij)
    end
    =#
    ifpart = """
        if zij <= 1 then
        \t\t    cij\$ = \"$(col2str(hm.cmap[1]))\"
        """
    for k ∈ 2:length(hm.cmap)-1
        ifpart *= """
            \n\t\telse if zij <= $k then
            cij\$ = \"$(col2str(hm.cmap[k]))\"
            """
    end
    ifpart *= """
        \n\t\telse
            cij\$ = \"$(col2str(hm.cmiss))\"
        end if"""

    box_default = """
        amove xg((j-1)*bw) yg(1-i*bh)
        \t\tbox xg(bw)-xg(0) yg(bh)-yg(0) nobox fill cij\$
        """
    box_transpose = """
        amove xg((i-1)*bh) yg(1-j*bw)
        \t\tbox xg(bh)-xg(0) yg(bw)-yg(0) nobox fill cij\$
        """
    boxpart = ifelse(hm.transpose, box_transpose, box_default)

    f.subroutines["hm_$hashid"] = """
        sub hm_$hashid j ds\$ bw bh
            local zij = 0
            local cij = \"\"
            for i = 1 to ndata(ds\$)
                zij = datayvalue(ds\$,i)
                $ifpart
                $boxpart
            next i
        end sub
        """
    return
end