text:
module TowardsNotebook

@time using Plots, DataFrames, Markdown
export notebook

function notebook()
    pv(fv::Float64, r::Float64, n::Float64) = fv / ((1 + r)^n)

    fv = 1.0
    years = [i for i in 0:20]
    qalyValue = DataFrame(
        :Year => years,
        :Value => [pv(fv, 0.035, Float64(n)) for n in years]
    )
    moneyValue = DataFrame(
        :Year => years,
        :Value => [pv(fv, 0.005, Float64(n)) for n in years]
    )

    @time p = plot(qalyValue.Year, qalyValue.Value, label = "QALY value loss at rate 3.5%/y")
    @time plot!(p, moneyValue.Year, moneyValue.Value, label = "Money value loss at rate 5%/y")

    @time savefig(p, "output.svg")

    md"""
    # QALY vs money discounting
    The below figure shows how not only money but also health outcomes (QALY)
    are expected to devaluate (due to incertanity and innovation)

    ![QALY vs money discounting](./output.svg)
    """
end

end # module

max_stars_repo_name: Jblew/towards-julia-notebook

max_issues_repo_name: Jblew/towards-julia-notebook

max_forks_repo_name: Jblew/towards-julia-notebook


text:
#!/usr/bin/env julia --project=@.
#
# TODO:
#   [x] basic features: --rs, -v, -i, -c
#   [x] additional features: --color, --mono with terminal detection
#   [ ] additional features multiple patterns, --and
#   [ ] multi-line record-separator
#   [ ] performance: use Boyre-moore or KMR for non-regex match
#   [ ] performance : mmap
#   [ ] performance: find-pattern-first algorithm for sparse matches
#
module JmlGrep

using ArgParse

const RS_DEFAULT = "^\$|^(=====*|-----*)\$"

re_dow="((Mon|Tue|Wed|Thu|Fri|Sat),?[ \t]+)?"
re_month="(Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Dec),?[ \t]*"
re_date="[0-9]{1,2},?"
re_time="[0-2][0-9]:[0-5][0-9](:[0-5][0-9])?"
re_year="(,?[ \t]20[0-9][0-9])?"
re_dty="$re_date[ \t]$re_time$re_year"
re_isodate="20[0-9][0-9]-(0[0-9]|11|12)-(0[1-9]|[12][0-9]|3[01])"
const REGEX_DATETIME = "^($re_dow$re_month$re_dty|$re_isodate)"

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--invert-match", "-v"
            help =  "Select non-matching lines (same as grep -v)."
            action = :store_true
        "--ignore-case",  "-i"
            help =  "Case-insensitive match (same as grep -i)."
            action = :store_true
        "--count", "-c"
            help =  "Print number of matches."
            action = :store_true
        "--rs", "-r"
            help = "Set input record separater. Default is: \"$RS_DEFAULT\""
            arg_type = String
            default = RS_DEFAULT
        # See https://github.com/JuliaLang/julia/pull/36671
        #"--color"
        #    help = "Highlight matches even when the output is not a terminal."
        #    action = :store_const
        #    constant = "always"
        #    default = "auto"
        #    dest_name = "highlight"
        "--mono"
            help = "Do not highlight matches even when the output is a terminal"
            arg_type = String
            action = :store_const
            constant = "never"
            default = "auto"
            dest_name = "highlight"
#        "--and", "-a"
#            help = "Print entries match with all, instead of any, of PATTERNS."
#            action = :store_true
        "--timestamp", "-t"
            help = "Synomym for --rs=REGEX_DATETIME, where the regex matches strings like '\n2014-12-31 12:34:56' or '\nDec 31 12:34:56'"
            action = :store_true
          "pattern"
              help = "An extended regex pattern."
              required = true
              nargs = '+'
              metavar = "PATTERN"
          "file"
              help = "Filename to search the pattern. Uses standard input if ommited."
              required = false
              nargs = '*'
              metavar = "FILE"
    end
    return parse_args(s)
end

struct RecLine
    line::String
    match::Union{RegexMatch,Nothing}
end

function print_record(rec, highlight::Bool)
    if highlight
        for recline in rec
            if isnothing(recline.match)
                println(recline.line)
            else
                pos = recline.match.offset
                sz = length(recline.match.match)
                str = recline.line
                print(str[1:pos-1])
                # See https://github.com/JuliaLang/julia/pull/36671
                printstyled(str[pos:(pos+sz)-1], bold=true, color=:red)
                println(str[pos+sz:end])
            end
        end
    else
        for recline in rec
            println(recline.line)
        end
    end
end

function mlgrep(input, rs::Regex, pat::Regex, highlight::Bool; invert=false::Bool, silent=false) :Int
    count = 0
    found = false
    rec = []
    for line in eachline(input)
        if !isnothing(match(rs, line))
            if found
                count += 1
                if !silent
                    print_record(rec, highlight)
                end
            end
            rec = []
            found = false
        end
        m = match(pat, line)
        push!(rec, RecLine(line, m))
        if xor(!isnothing(m), invert)
            found = true
        end
    end
    if found
        if !silent
            print_record(rec, highlight)
        end
        count += 1
    end

    return count
end

function pattern_filename_split(strings)
    patterns = [strings[1]]
    filenames = []
    is_pattern = true
    for s in strings[2:end]
        st = stat(s)
        if is_pattern && ( s == "-" || (st.mode != 0x0 && st.size > 0))
            is_pattern=false
        end
        if is_pattern
            push!(patterns, s)
        else
            push!(filenames, s)
        end
    end
    return (patterns, filenames)
end

function process_args(args)
    if length(args["pattern"]) == 0
        println(stderr, "At least one argument (a pattern) is required")
        exit(2)
    end

    patterns, filenames = pattern_filename_split(args["pattern"])

    if length(patterns) > 1
        println(stderr, "ERROR: only one pattern is supported at the moment.")
        exit(2)
    end

    if length(filenames) == 0
        filenames = ["-"]
    end

    if args["ignore-case"]
        reFlags = "i"
    else
        reFlags = ""
    end

    pat = Regex(patterns[1], reFlags)
    if args["timestamp"]
        rs = Regex(REGEX_DATETIME, reFlags)
    else
        rs = Regex(args["rs"], reFlags)
    end

    if args["highlight"] == "always"
        highlight = true
    elseif args["highlight"] == "never"
        highlight = false
    else
       highlight = isa(stdout, Base.TTY)
    end

    return (patterns, filenames, reFlags, rs, pat, highlight)
end

function main()
    args = parse_commandline()

    patterns, filenames, reFlags, rs, pat, highlight = process_args(args)

    count = 0
    for f in filenames
        if f == "-"
            input = stdin
        else
            input = open(f)
        end
        if args["count"]
            this_count = mlgrep(input, rs, pat, highlight, invert=args["invert-match"], silent=true)
            if length(filenames) > 1
                print("$f:")
            end
            println(this_count)
            count += this_count
        else
            count += mlgrep(input, rs, pat, highlight, invert=args["invert-match"])
        end
    end

    if count == 0
        exit(1)
    else
        exit(0)
    end
end

function julia_main()::Cint
    try
        ret = main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return ret
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

end #module

ext: jl


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    abstract type AbstractProposalDist

*BAT-internal, not part of stable public API.*

The following functions must be implemented for subtypes:

* `BAT.proposaldist_logpdf`
* `BAT.proposal_rand!`
* `ValueShapes.totalndof`, returning the number of DOF (i.e. dimensionality).
* `LinearAlgebra.issymmetric`, indicating whether p(a -> b) == p(b -> a) holds true.
"""
abstract type AbstractProposalDist end

"""
    proposaldist_logpdf(
        p::AbstractArray,
        pdist::AbstractProposalDist,
        v_proposed::AbstractVector,
        v_current:::AbstractVector
    )

*BAT-internal, not part of stable public API.*

Returns log(PDF) value of `pdist` for transitioning from current to proposed
variate/parameters.
"""
function proposaldist_logpdf end

# TODO: Implement proposaldist_logpdf for included proposal distributions

"""
    function proposal_rand!(
        rng::AbstractRNG,
        pdist::GenericProposalDist,
        v_proposed::Union{AbstractVector,VectorOfSimilarVectors},
        v_current::Union{AbstractVector,VectorOfSimilarVectors}
    )

*BAT-internal, not part of stable public API.*

Generate one or multiple proposed variate/parameter vectors, based on one or
multiple previous vectors.

Input:

* `rng`: Random number generator to use
* `pdist`: Proposal distribution to use
* `v_current`: Old values (vector or column vectors, if a matrix)

Output is stored in

* `v_proposed`: New values (vector or column vectors, if a matrix)

The caller must guarantee:

* `size(v_current, 1) == size(v_proposed, 1)`
* `size(v_current, 2) == size(v_proposed, 2)` or `size(v_current, 2) == 1`
* `v_proposed !== v_current` (no aliasing)

Implementations of `proposal_rand!` must be thread-safe.
"""
function proposal_rand! end

struct GenericProposalDist{D<:Distribution{Multivariate},SamplerF,S<:Sampleable} <: AbstractProposalDist
    d::D
    sampler_f::SamplerF
    s::S

    function GenericProposalDist{D,SamplerF}(d::D, sampler_f::SamplerF) where {D<:Distribution{Multivariate},SamplerF}
        s = sampler_f(d)
        new{D,SamplerF, typeof(s)}(d, sampler_f, s)
    end

end

GenericProposalDist(d::D, sampler_f::SamplerF) where {D<:Distribution{Multivariate},SamplerF} =
    GenericProposalDist{D,SamplerF}(d, sampler_f)

GenericProposalDist(d::Distribution{Multivariate}) = GenericProposalDist(d, bat_sampler)

GenericProposalDist(D::Type{<:Distribution{Multivariate}}, varndof::Integer, args...) =
    GenericProposalDist(D, Float64, varndof, args...)

Base.similar(q::GenericProposalDist, d::Distribution{Multivariate}) =
    GenericProposalDist(d, q.sampler_f)

function Base.convert(::Type{AbstractProposalDist}, q::GenericProposalDist, T::Type{<:AbstractFloat}, varndof::Integer)
    varndof != totalndof(q) && throw(ArgumentError("q has wrong number of DOF"))
    q
end

get_cov(q::GenericProposalDist) = get_cov(q.d)
set_cov(q::GenericProposalDist, Σ::PosDefMatLike) = similar(q, set_cov(q.d, Σ))

function proposaldist_logpdf(
    pdist::GenericProposalDist,
    v_proposed::AbstractVector,
    v_current::AbstractVector
)
    params_diff = v_proposed .- v_current # TODO: Avoid memory allocation
    logpdf(pdist.d, params_diff)
end

function proposal_rand!(
    rng::AbstractRNG,
    pdist::GenericProposalDist,
    v_proposed::Union{AbstractVector,VectorOfSimilarVectors},
    v_current::Union{AbstractVector,VectorOfSimilarVectors}
)
    rand!(rng, pdist.s, flatview(v_proposed))
    params_new_flat = flatview(v_proposed)
    params_new_flat .+= flatview(v_current)
    v_proposed
end

ValueShapes.totalndof(pdist::GenericProposalDist) = length(pdist.d)

LinearAlgebra.issymmetric(pdist::GenericProposalDist) = issymmetric_around_origin(pdist.d)

struct GenericUvProposalDist{D<:Distribution{Univariate},T<:Real,SamplerF,S<:Sampleable} <: AbstractProposalDist
    d::D
    scale::Vector{T}
    sampler_f::SamplerF
    s::S
end

GenericUvProposalDist(d::Distribution{Univariate}, scale::Vector{<:AbstractFloat}, samplerF) =
    GenericUvProposalDist(d, scale, samplerF, samplerF(d))

GenericUvProposalDist(d::Distribution{Univariate}, scale::Vector{<:AbstractFloat}) =
    GenericUvProposalDist(d, scale, bat_sampler)

ValueShapes.totalndof(pdist::GenericUvProposalDist) = size(pdist.scale, 1)

LinearAlgebra.issymmetric(pdist::GenericUvProposalDist) = issymmetric_around_origin(pdist.d)

function BAT.proposaldist_logpdf(
    pdist::GenericUvProposalDist,
    v_proposed::Union{AbstractVector,VectorOfSimilarVectors},
    v_current::Union{AbstractVector,VectorOfSimilarVectors}
)
    params_diff = (flatview(v_proposed) .- flatview(v_current)) ./ pdist.scale  # TODO: Avoid memory allocation
    sum_first_dim(logpdf.(pdist.d, params_diff))  # TODO: Avoid memory allocation
end

function BAT.proposal_rand!(
    rng::AbstractRNG,
    pdist::GenericUvProposalDist,
    v_proposed::AbstractVector,
    v_current::AbstractVector
)
    v_proposed .= v_current
    dim = rand(rng, eachindex(pdist.scale))
    v_proposed[dim] += pdist.scale[dim] * rand(rng, pdist.s)
    v_proposed
end

abstract type ProposalDistSpec end

struct MvTDistProposal <: ProposalDistSpec
    df::Float64
end

MvTDistProposal() = MvTDistProposal(1.0)

(ps::MvTDistProposal)(T::Type{<:AbstractFloat}, varndof::Integer) =
    GenericProposalDist(MvTDist, T, varndof, convert(T, ps.df))

function GenericProposalDist(::Type{MvTDist}, T::Type{<:AbstractFloat}, varndof::Integer, df = one(T))
    Σ = PDMat(Matrix(ScalMat(varndof, one(T))))
    μ = Fill(zero(eltype(Σ)), varndof)
    M = typeof(Σ)
    d = Distributions.GenericMvTDist(convert(T, df), μ, Σ)
    GenericProposalDist(d)
end

struct UvTDistProposalSpec <: ProposalDistSpec
    df::Float64
end

(ps::UvTDistProposalSpec)(T::Type{<:AbstractFloat}, varndof::Integer) =
    GenericUvProposalDist(TDist(convert(T, ps.df)), fill(one(T), varndof))


text:
function tsort(nvertices, egress)
    # TODO could output list of independent nodes, such that parallel eval of
    # the lists are possible.
    color = zeros(Int, nvertices)
    order = Vector{Int}()
    stack = Vector{Int}()
    for v in 1:nvertices
        color[v] != 0 && continue
        push!(stack, v)
        while !isempty(stack)
            cur = stack[end]
            color[cur] = 1
            for child in egress(cur)
                assert(color[child] != 1)  # cycle (not a DAG)
                color[child] != 0 && continue
                push!(stack, child)
            end
            if cur == stack[end]
                color[cur] = 2
                push!(order, cur)
                pop!(stack)
            end
        end
    end
    reverse(order)
end

struct EvalTSort end
# TODO maybe pass Dag, assume static or handle Dag version number
# TODO this also assume Dag node are number 1 to N (start at 1)

# TODO fire interface should accept multiple nodes at once (event group)
function fire(eve::EvalTSort, d::Dag, src_nid, src_v)
    # TODO assume Dag is static (or at least no new node can fire during fire)
    #      SimpleTrait could help here
    # TODO caching of tsort result, reuse of `firing` memory etc
    t_order = tsort(length(d.nodes), nid -> [n for (n,_) in d.links[nid] if n != 0])

    # TODO maybe have array of nullable instead of dict, and iterate in zip
    firing = Dict{Int, Any}(src_nid=>src_v)
    for nid in t_order
        !haskey(firing, nid) && continue

        nid_v = firing[nid]
        # TODO could remove from firing (and stop when firing is empty)

        targets = Base.get(d.links, nid, [])
        for (dst_nid, f) in targets
            # TODO ugly hack to register external callback
            if dst_nid == 0
                f(nid_v)
                continue
            end

            # TODO store dst directly in targets?
            dst = d.nodes[dst_nid]
            # assert dst must be in rest of t_order

            tv = f(nid_v)  # feed builder edge buffer
            if !isnull(tv)
                # nodes are evaluated in sequence, eliding previous values
                # i.e. if a node produce multiple values (triggered by multiple
                # parents) only the latest produced value is kept.
                firing[dst.nid] = get(tv)
            end
        end
    end
end


text:
export tip, resolve, owner, name, canonical_name, ishead, move,
       remote_name, remote, isremote, upstream, set_upstream!

Base.delete!(b::GitBranch) = begin
    @check ccall((:git_branch_delete, libgit2), Cint, (Ptr{Void},), b)
    return nothing
end

Base.(:(==))(b1::GitBranch, b2::GitBranch) = canonical_name(b1) == canonical_name(b2)
Base.isequal(b1::GitBranch, b2::GitBranch) = b1 == b2

tip(b::GitBranch) = lookup(owner(b), target(resolve(b)))

function owner(b::GitBranch)
    repo_ptr = ccall((:git_reference_owner, libgit2), Ptr{Void}, (Ptr{Void},), b)
    # do not free repository handle, pointer to repo is not owned by julia runtime
    return GitRepo(repo_ptr, false)
end

function resolve(b::GitBranch)
    ref_ptr = Ptr{Void}[0]
    @check ccall((:git_reference_resolve, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}), ref_ptr, b)
    return GitReference(ref_ptr[1])
end

function name(b::GitBranch)
    name_ptr = Ptr{Uint8}[0]
    @check ccall((:git_branch_name, libgit2), Cint,
                 (Ptr{Ptr{Uint8}}, Ptr{Void}), name_ptr, b)
    return utf8(bytestring(name_ptr[1]))
end

ishead(b::GitBranch)  =
    bool(ccall((:git_branch_is_head, libgit2), Cint, (Ptr{Void},), b))

isremote(b::GitBranch) =
    bool(ccall((:git_reference_is_remote, libgit2), Cint, (Ptr{Void},), b))

canonical_name(b::GitBranch) =
    utf8(bytestring(ccall((:git_reference_name, libgit2), Ptr{Uint8}, (Ptr{Void},), b)))

function move(b::GitBranch, newname::String;
              force::Bool=false,
              sig::Union(Nothing, Signature)=nothing,
              logmsg::Union(Nothing, String)=nothing)
    branch_ptr = Ptr{Void}[0]
    if sig != nothing
        @check ccall((:git_branch_move, libgit2), Cint,
                      (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Uint8},
                       Cint, Ptr{SignatureStruct}, Ptr{Uint8}),
                       branch_ptr, b, newname, force, sig,
                       logmsg != nothing ? logmsg : C_NULL)
    else
        @check ccall((:git_branch_move, libgit2), Cint,
                      (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Uint8},
                       Cint, Ptr{SignatureStruct}, Ptr{Uint8}),
                       branch_ptr, b, newname, force, C_NULL,
                       logmsg != nothing ? logmsg : C_NULL)
    end
    return GitBranch(branch_ptr[1])
end

function remote_name(b::GitBranch)
    ref_ptr = Ptr{Void}[0]
    if isremote(b)
        ref_ptr[1] = b.ptr
    else
        err = ccall((:git_branch_upstream, libgit2), Cint,
                    (Ptr{Ptr{Void}}, Ptr{Void}), ref_ptr, b)
        if err == GitErrorConst.ENOTFOUND
            return nothing
        elseif err != GitErrorConst.GIT_OK
            if ref_ptr[1] != C_NULL
                ccall((:git_reference_free, libgit2), Void, (Ptr{Void},), ref_ptr[1])
            end
            throw(GitError(err))
        end
    end
    repo_ptr = ccall((:git_reference_owner, libgit2), Ptr{Void}, (Ptr{Void},), b)
    refname_ptr = ccall((:git_reference_name,  libgit2), Ptr{Uint8}, (Ptr{Void},), ref_ptr[1])
    buf_ptr = [BufferStruct()]
    err = ccall((:git_branch_remote_name, libgit2), Cint,
                (Ptr{BufferStruct}, Ptr{Void}, Ptr{Uint8}),
                buf_ptr, repo_ptr, refname_ptr)
    buf = buf_ptr[1]
    if err == GitErrorConst.GIT_OK
        str = utf8(bytestring(buf))
        ccall((:git_buf_free, libgit2), Void, (Ptr{BufferStruct},), buf_ptr)
    elseif err != GitErrorConst.GIT_OK
        if buf.ptr != C_NULL
            ccall((:git_buf_free, libgit2), Void, (Ptr{BufferStruct},), buf_ptr)
        end
        throw(GitError(err))
    end
    return str
end

function remote(b::GitBranch)
    if remote_name(b) != nothing
        return lookup_remote(owner(b), "origin")
    end
    return nothing
end

function upstream(b::GitBranch)
    isremote = ccall((:git_reference_is_remote, libgit2), Cint, (Ptr{Void},), b)
    if bool(isremote)
        return nothing
    end
    ubranch_ptr = Ptr{Void}[0]
    err = ccall((:git_branch_upstream, libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}), ubranch_ptr, b)
    if err == GitErrorConst.ENOTFOUND
        return nothing
    elseif err != GitErrorConst.GIT_OK
        throw(GitError(err))
    end
    return GitBranch(ubranch_ptr[1])
end

function set_upstream!(b::GitBranch, target::Nothing)
    @check ccall((:git_branch_set_upstream, libgit2), Cint,
                 (Ptr{Void}, Ptr{Uint8}), b, C_NULL)
    return b
end

function set_upstream!(b::GitBranch, target::Union(GitBranch, GitReference))
    name_ptr = Ptr{Uint8}[0]
    @check ccall((:git_branch_name, libgit2), Cint,
                 (Ptr{Ptr{Uint8}}, Ptr{Void}), name_ptr, target)
    @check ccall((:git_branch_set_upstream, libgit2), Cint,
                 (Ptr{Void}, Ptr{Uint8}), b, name_ptr[1])
    return b
end

function rename(b::GitBranch, new_name::String;
                force::Bool=false, sig=nothing, logmsg=nothing)
    return move(b, new_name, force=force, sig=sig, logmsg=logmsg)
end


text:
"""
    geoarray(filename; kw...) => AbstractGeoArray

Load a file path as `AbstractGeoArray`.

# Keywords

Passed to the constructor for the file type, and commmonly include:

- `crs`: crs to use instead of the detected crs.
- `mappedcrs`: CRS format like `EPSG(4326)` used in `Selectors` like `Between` and `At`, and
    for plotting. Can be any CRS `GeoFormat` from GeoFormatTypes.jl, like `WellKnownText`.
- `name`: `Symbol` name for the array.
- `dims`: `Tuple` of `Dimension`s for the array. Detected automatically, but can be passed in.
- `refdims`: `Tuple of` position `Dimension`s the array was sliced from.
- `missingval`: Value reprsenting missing values. Detected automatically when possible, but
    can be passed it.
- `metadata`: `Metadata` object for the array. Detected automatically but can be passed in.
"""
function geoarray end

geoarray(filename::AbstractString; kw...) = _constructor(geoarray, filename)(filename; kw...)

"""
    Base.write(filename::AbstractString, A::AbstractGeoArray; kw...)

Write any [`AbstractGeoArray`](@ref) to file, guessing the backend from the file extension.

Keyword arguments are passed to the `write` method for the backend.
"""
function Base.write(
    filename::AbstractString, A::AbstractGeoArray; kw...
)
    write(filename, _constructor(geoarray, filename), A; kw...)
end

"""
    stack(filename; kw...) => AbstractGeoStack
    stack(filenames::NamedTuple; kw...) => AbstractGeoStack

Load a file path or a `NamedTuple` of paths as an `AbstractGeoStack`.

# Arguments

- `filename`: A `NamedTuple` of stack keys and `String` filenames.

# Keywords

Passed to the constructor for the file type, and commmonly include:

- `window`: A `Tuple` of `Dimension`/`Selector`/indices that will be applied to the
    contained arrays when they are accessed.
- `metadata`: `Metadata` as object.
- `child_kwargs`: A `NamedTuple` of keyword arguments to pass to the `childtype` constructor.
- `refdims`: `Tuple` of  position `Dimension` the array was sliced from.

```julia
files = (:temp="temp.tif", :pressure="pressure.tif", :relhum="relhum.tif")
stack = stack(files; child_kwargs=(mappedcrs=EPSG(4326),))
stack[:relhum][Lat(Contains(-37), Lon(Contains(144))
```
"""
function stack end

stack(filename::AbstractString; kw...) = _constructor(stack, filename)(filename; kw...)
stack(filenames::Union{Tuple,AbstractArray}; keys, kw...) =
    stack(NamedTuple{map(Symbol, keys)}(filenames); kw...)
function stack(filenames::NamedTuple; child_kwargs=(), kw...)
    if all(x -> splitext(x)[2] == splitext(first(filenames))[2], filenames)
        # All the same kind of file. We don't need to load them up front.
        DiskStack(filenames; childtype=_constructor(geoarray, first(filenames)), kw...)
    else
        # These files are different extensions, just load them all
        # as separate `AbstarctGeoArray` (which has some up front cost from
        # reading the dimensions). They are probably still disk-backed for the
        # actual array.
        arrays = map(filenames) do fn
            _constructor(geoarray, fn)(fn; child_kwargs...)
        end
        GeoStack(arrays; kw...)
    end
end

"""
    Base.write(filename::AbstractString, T::Type{<:AbstractGeoArray}, s::AbstractGeoStack)

Write any [`AbstractGeoStack`](@ref) to file, guessing the backend from the file extension.

Keyword arguments are passed to the `write` method for the backend.

If the source can't be saved as a stack-like object, individual array layers will be saved.
"""
function Base.write(filename::AbstractString, s::AbstractGeoStack; kw...)
    base, ext = splitext(filename)
    T = _constructor(stack, filename; throw=false)
    if T isa Nothing
        # Save as individual `geoarray`
        T = _constructor(geoarray, filename)
        for key in keys(s)
            fn = joinpath(string(base, "_", key, ext))
            write(fn, T, s[key])
        end
    else
        write(filename, _constructor(stack, filename), s; kw...)
    end
end

"""
    series(dirpath::AbstractString, dims; ext, child=geoarray, kw...) => AbstractGeoSeries
    series(filenames::AbstractVector{String}, dims; child=geoarray, kw...) => AbstractGeoSeries

Load a whole folder vector of filepaths as a `AbstractGeoSeries`.
`kw` are passed to the constructor.

# Arguments

- `dims`: A tuple of dimensions, possibly holding an index matching the files in the directory.
    By default it will be a numbered `Dim{:series}`.

# Keywords

- `child`: function to load the child object - may be `geoarray` or `stack`,
    `geoarray` by default.
- `ext`: an extension for the files in a directory. If empty, all files are loaded,
    and should be the same type.
"""
function series end

function series(dirpath::AbstractString, dims=(Dim{:series}(),); ext=nothing, child=geoarray, kw...)
    filepaths = filter_ext(dirpath, ext)
    series(filepaths, dims; child=child, kw...)
end
function series(filepaths::AbstractVector{<:AbstractString}, dims=(Dim{:series}(),); child=geoarray, kw...)
    childtype = _constructor(child, first(filepaths))
    GeoSeries(filepaths, dims; childtype=childtype, kw...)
end

# Support methods

const EXT = (GRD=(".grd", ".gri"), NCD=".nc", SMAP=".h5")

# The the constructor for a geoarray or stack, based on the
# filename extension. GDAL is the fallback for geoarray as it
# handles so many file types.
function _constructor(method::typeof(geoarray), filename; throw=true)
    _, extension = splitext(filename)
    return if extension in EXT[:GRD]
        GRDarray
    elseif extension == EXT[:NCD]
        _check_imported(:NCDatasets, :NCDarray, extension; throw=throw)
        NCDarray
    elseif extension == EXT[:SMAP]
        # In future we may need to examine the file and check if
        # it's a SMAP file or something else that uses .h5
        throw ? _no_gearray_error(extension) : nothing
    else # GDAL handles too many extensions to list, so just try it and see if it works
        _check_imported(:ArchGDAL, :GDALarray, extension; throw=throw)
        GDALarray
    end
end
function _constructor(method::typeof(stack), filename; throw=true)
    _, extension = splitext(filename)
    return if extension in EXT[:GRD]
        throw ? _no_stack_error(extension) : nothing
    elseif extension == EXT[:NCD]
        _check_imported(:NCDatasets, :NCDarray, extension; throw=throw)
        NCDstack
    elseif extension == EXT[:SMAP]
        _check_imported(:HDF5, :SMAPstack, extension; throw=throw)
        SMAPstack
    else
        throw ? _no_stack_error(extension) : nothing
    end
end

_no_stack_error(ext) =
    error("$ext files do not have named layers. Use `geoarray(filename)` to load, or `stack((key1=fn1, key2=fn2, ...)`")

_no_gearray_error(ext) =
    error("$ext files not have a single-layer implementation. Use `stack(filename)` to load")

function _check_imported(modulename, type, extension; throw=true)
    if type in names(GeoData)
        true
    else
        if throw
            error("Run `import $modulename` to enable loading $extension files.")
        else
            false
        end
    end
end


text:
using RandomMatrices
using Test

@testset "StochasticProcess" begin

seed!(1)
dx = 0.001

let
p = WhiteNoiseProcess(dx)
x, S = iterate(p)
@test typeof(x) == typeof(dx)
@test isa(S, Tuple{})
@test iterate(p, S) != nothing
end

let
p = BrownianProcess(dx)
x, S = iterate(p)
@test x == S
@test typeof(x) == typeof(dx)
@test iterate(p, S) != nothing
end

let
p = AiryProcess(dx, 2.0)
x, S = iterate(p)
@test S[1,1] == -2/dx^2
@test next!(p, S) != nothing
end

end # testset


text:
# General DEQ Utils
mutable struct DEQTrainingStats
    nfe::Int
end

"""
    get_and_clear_nfe!(model::AbstractDeepEquilibriumNetwork)

Return the number of function evaluations (NFE) and clear the counter.
"""
function get_and_clear_nfe!(model::AbstractDeepEquilibriumNetwork)
    nfe = model.stats.nfe
    model.stats.nfe = 0
    return nfe
end

"""
    SteadyStateAdjoint(reltol, abstol, maxiters; autojacvec=ZygoteVJP(),
                       linsolve=KrylovJL_GMRES(; rtol=reltol, atol=abstol, itmax=maxiters))

Creates SteadyStateAdjoint ([johnson2012notes](@cite)) with sensible defaults.

## Arguments

* `reltol`: Relative tolerance.
* `abstol`: Absolute tolerance.
* `maxiters`: Maximum number of iterations.
* `autojacvec`: Which backend to use for VJP.
* `linsolve`: Linear Solver from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl).
"""
function DiffEqSensitivity.SteadyStateAdjoint(reltol, abstol, maxiters; autojacvec=ZygoteVJP(),
                                              linsolve=KrylovJL_GMRES(; rtol=reltol, atol=abstol, itmax=maxiters))
    return SteadyStateAdjoint(; autodiff=true, autojacvec=autojacvec, linsolve=linsolve)
end

# Initialization
"""
    NormalInitializer(μ = 0.0f0, σ² = 0.01f0)

Initializes the weights of the network with a normal distribution. For DEQs the training is stable
if we use this as the Initialization
"""
function NormalInitializer(μ = 0.0f0, σ² = 0.01f0)
    return (dims...) -> randn(dims...) .* σ² .+ μ
end

# Wrapper for Discrete DEQs
"""
    DiscreteDEQSolver(solver=LimitedMemoryBroydenSolver; abstol=1e-8, reltol=1e-8, kwargs...)

Solver for Discrete DEQ Problem ([baideep2019](@cite)). A wrapper around `SSrootfind` to mimic the [`ContinuousDEQSolver`](@ref) API.

## Arguments

* `solver`: NonLinear Solver for the DEQ problem. (Default: [`LimitedMemoryBroydenSolver`](@ref))
* `abstol`: Absolute tolerance for termination. (Default: `1e-8`)
* `reltol`: Relative tolerance for termination. (Default: `1e-8`)
* `kwargs`: Additional keyword arguments passed to the solver.

!!! note
    There is no `mode` kwarg for [`DiscreteDEQSolver`](@ref). Instead solvers directly define their own termination condition.
    For [`BroydenSolver`](@ref) and [`LimitedMemoryBroydenSolver`](@ref), the termination conditions are `:abs_norm` &
    `:rel_deq_default` respectively.

See also: [`ContinuousDEQSolver`](@ref)
"""
function DiscreteDEQSolver(solver=LimitedMemoryBroydenSolver; abstol=1e-8, reltol=1e-8, kwargs...)
    solver = solver(; kwargs..., reltol=reltol, abstol=abstol)
    return SSRootfind(; nlsolve=(f, u0, abstol) -> solver(f, u0))
end

# For MultiScale DEQs
function split_array_by_indices(x::AbstractVector, idxs)
    return collect((x[(i + 1):j] for (i, j) in zip(idxs[1:(end - 1)], idxs[2:end])))
end

function split_array_by_indices(x::AbstractMatrix, idxs)
    return collect((x[(i + 1):j, :] for (i, j) in zip(idxs[1:(end - 1)], idxs[2:end])))
end

Zygote.@adjoint function split_array_by_indices(x, idxs)
    res = split_array_by_indices(x, idxs)
    function split_array_by_indices_sensitivity(Δ)
        is_nothings = Δ .=== nothing
        if any(is_nothings)
            Δ[is_nothings] .= zero.(res[is_nothings])
        end
        return (vcat(Δ...), nothing)
    end
    return res, split_array_by_indices_sensitivity
end

# Zygote Fix
function Zygote.accum(x::NTuple{N,T}, y::AbstractVector{T}) where {N,T<:AbstractArray}
    return Zygote.accum.(x, y)
end

function Zygote.accum(x::AbstractVector{T}, y::NTuple{N,T}) where {N,T<:AbstractArray}
    return Zygote.accum.(x, y)
end

function Zygote.accum(x::AbstractVector{T}, y::NTuple{N,Nothing}) where {N,T<:AbstractArray}
    return Zygote.accum.(x, y)
end

# General Utils
@inline function _init_identity_matrix(x::AbstractArray{T}, scale::T=T(1)) where {T}
    x_ = vec(x)
    return _init_identity_matrix!(x_ .* x_', scale)
end

@inline function _init_identity_matrix!(x::AbstractMatrix{T}, scale::T=T(1)) where {T}
    x .= zero(T)
    idxs = diagind(x)
    @. @view(x[idxs]) = scale * true
    return x
end

@inline function _norm(x; dims=Colon())
    return sqrt.(sum(abs2, x; dims=dims))
end

# Compute norm over all dimensions except `except_dim`
@inline function _norm(x::AbstractArray{T,N}, except_dim) where {T,N}
    dims = filter(i -> i != except_dim, 1:N)
    return _norm(x; dims=dims)
end

flatten_merge(x, y) = (x..., y...)
flatten_merge(x::T, y::T) where {T<:AbstractArray} = (x, y)
flatten_merge(x::NTuple{N,T}, y::T) where {N,T<:AbstractArray} = (x..., y)
flatten_merge(x::T, y::NTuple{N,T}) where {N,T<:AbstractArray} = (x, y...)
flatten_merge(x::NTuple{N,T}, y) where {N,T<:AbstractArray} = (x, y...)
flatten_merge(x, y::NTuple{N,T}) where {N,T<:AbstractArray} = (x..., y)
function flatten_merge(x::NTuple{N,T}, y::NTuple{N,T}) where {N,T<:AbstractArray}
    return (x, y)
end


text:
## index.jl  methods for NamedArray that keep the names (some checking may be done)

## (c) 2013, 2014 David A. van Leeuwen

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

# Keep names for consistently named vectors, or drop them
function Base.hcat{T}(V::NamedVector{T}...)
    keepnames=true
    V1=V[1]
    firstnames = names(V1,1)
    for i=2:length(V)
        keepnames &= names(V[i],1)==firstnames
    end
    a = hcat(map(a -> a.array, V)...)
    if keepnames
        colnames = [string(i) for i=1:size(a,2)]
        NamedArray(a, (firstnames, colnames), (V1.dimnames[1], "hcat"))
    else
        NamedArray(a)
    end
end

## helper function for broadcast
function verify_names(a::NamedArray...)
    nargs = length(a)
    @assert nargs>1
    bigi = indmax(map(length, a)) # find biggest dimension
    big = a[bigi]
    for i in setdiff(bigi, 1:nargs)
        println("i ", i)
        for d=1:ndims(a[i])
            println("d ", d)
            if size(a[i],d) > 1
                @assert names(big,d) == names(a(i),d)
            end
        end
    end
    return bigi::Int, big
end

## broadcast
import Base.broadcast, Base.broadcast!
function broadcast(f::Function, a::NamedArray...)
    ## verify that the names are consistent
    bigi, big = verify_names(a...)
    arrays = map(x->x.array, a)
    NamedArray(broadcast(f, arrays...), big.dicts, big.dimnames)
end

function broadcast!(f::Function, dest::NamedArray, a::NamedArray, b::NamedArray...)
    ab = tuple(a, b...)
    ## verify that the names are consistent, we assume dest is the right size
    bigi, big = verify_names(ab...)
    arrays = map(x->x.array, ab)
    broadcast!(f, dest.array, arrays...)
    dest
end

## keep names intact
for f in (:sin, :cos, :tan, :sind, :cosd, :tand, :sinpi, :cospi, :sinh, :cosh, :tanh, :asin, :acos, :atan, :asind, :acosd, :sec, :csc, :cot, :secd, :cscd, :cotd, :asec, :acsc, :asecd, :acscd, :acotd, :sech, :csch, :coth, :asinh, :acosh, :atanh, :asech, :acsch, :acoth, :sinc, :cosc, :deg2rad, :log, :log2, :log10, :log1p, :exp, :exp2, :exp10, :expm1, :ceil, :floor, :trunc, :round, :abs, :abs2, :sign, :signbit, :sqrt, :isqrt, :cbrt, :erf, :erfc, :erfcx, :erfi, :dawson, :erfinv, :erfcinv, :real, :imag, :conj, :angle, :cis, :gamma, :lgamma, :digamma, :invdigamma, :trigamma, :airyai, :airyprime, :airyaiprime, :airybi, :airybiprime, :besselj0, :besselj1, :bessely0, :bessely1, :eta, :zeta)
    eval(Expr(:import, :Base, f))
    @eval ($f)(a::NamedArray) = NamedArray(($f)(a.array), a.dicts, a.dimnames)
end

## reorder names
import Base: sort, sort!
function sort!(a::NamedVector; kws...)
    i = sortperm(a.array; kws...)
    newnames = names(a, 1)[i]
    empty!(a.dicts[1])
    for (ind, k) in enumerate(newnames)
        a.dicts[1][k] = ind
    end
    a.array = a.array[i]
    a
end

function sort(a::NamedVector; kws...)
    return sort!(copy(a); kws...)
end

## Note: I can't think of a sensible way to define sort!(a::NamedArray, dim>1)

## drop name of sorted dimension, as each index along that dimension is sorted individually
function sort(a::NamedArray, dim::Integer; kws...)
    if ndims(a)==1 && dim==1
        return sort(a; kws...)
    else
        n = allnames(a)
        n[dim] = [string(i) for i in 1:size(a,dim)]
        return NamedArray(sort(a.array, dim; kws...), tuple(n...), a.dimnames)
    end
end


text:
Vector{u_T}(undef, n, d) where u_T = Vector{u_T}[Vector{u_T}(undef, d) for i = 1:n]
Vector{u_T}(undef, n₂, n₁, d) where u_T = Vector{Vector{u_T}}[Vector{u_T}(undef, n₁, d) for i = 1:n₂]

zero!(v::AbstractVector) = fill!(v, zero(eltype(v)))
function zero!(v::AbstractVector{<:AbstractVector})
    for i in eachindex(v)
        zero!(v[i])
    end
    return v
end

function norm!(v::AbstractVector{<:AbstractVector})
    r = zero(eltype(v))
    for i in eachindex(v)
        r += norm(v[i])
    end
    return r
end


text:
module TestShowDispatch
using PatternDispatch

# @pattern f(x::Int) = x^2
# @pattern f(3)      = 3
# @pattern f(x~::AbstractString) = 5

@pattern f(x,y) = 1
@pattern f(1,y) = 2
@pattern f(1,2) = 3

println("==== Full dispatch: ====")
show_dispatch(f)

println("\n==== Dispatch for (Int,Any): ====")
show_dispatch(f, (Int,Any))

end


text:
include("preamble.jl")

# ----------------------------
# Plotting a 3D polyhedron
# ----------------------------

# Note: in this figure we use non-latex fonts for the ticks labels because
# the Makie plot ignores L"x_3"

using LazySets, Symbolics
import Polyhedra, CDDLib, Plots
using Colors

# projected polyhedron from symbolic constraints
vars = @variables x1, x2, x3  # create symbols
P = HPolyhedron([2*x1>=0, 3*x2+1.7*x3>=0], vars)
Q = project(P, [2, 3])  # 2D projection (x2 and x3)

fig = p(xlim=(-4, 4), ylim=(-4, 4), ratio=1, xlab="x₂", ylab="x₃",
           legendfontsize=20,
           tickfont=font(20), # , "Times"),
           guidefontsize=20,
           xguidefont=font(20), # , "Times"),
           yguidefont=font(20), # , "Times"),
           xtick = ([-4, 0, 4], ["-4", "0", "4"]),
           ytick = ([-4, 0, 4], ["-4", "0", "4"]),
           bottom_margin=0mm, left_margin=0mm, right_margin=0mm, top_margin=2mm, size=(500, 380))

plot!(fig, Q, color=:deepskyblue)  # plot the 2D projection

savefig("polyhedron2D.pdf")

# plot the 3D set directly
using GLMakie

# intersection with a bounding box
B = BallInf(zeros(3), 5.0)
R = intersection(P, B)

# obtain the mesh
backend = LazySets.default_polyhedra_backend(R);
P_poly_mesh = LazySets._plot3d_helper(R, backend);

fig3d = Figure()
ax = Axis3(fig3d[1, 1], xlabel="x₁", xlabelsize=40, xticklabelsize=40,
                    ylabel="x₂", ylabelsize=40, yticklabelsize=40,
                    zlabel="x₃", zlabelsize=40, zticklabelsize=40,
                    aspect=:data,
                    height=400)

mesh!(P_poly_mesh, alpha=1.0, color=:deepskyblue, colormap=:viridis, colorrange=Makie.Automatic(), interpolate=false,
                   linewidth=1, overdraw=false, shading=true, transparency=true, visible=true)

wireframe!(P_poly_mesh, color = :red, linewidth=1.5, linestyle=:dash, axis=(type=Axis3,))

Makie.save("polyhedron3D.png", fig3d)

# f = Figure()
# fig = plot3d(R)  # plot in three dimensions
# current_figure()
# display(fig) # display in a pop-up window


text:
# Logically, wild bootstrap tests perform estimation at two stages, once as part of the bootstrap DGP, once in each bootstrap replication
# The StrEstimator "class" and its three "children" hold the estimation logic for the OLS, Anderson-Rubin, and IV/GMM cases

function perp(A::AbstractMatrix)
  F = eigen(Symmetric(A*invsym(A'A)*A'))
  F.vectors[:, abs.(F.values) .< 1000eps(eltype(A))]
end

# R₁ is constraints. R is attack surface for null; only needed when using FWL for WRE
# for DGP regression, R₁ is maintained constraints + null if imposed while R should have 0 nrows
# for replication regressions R₁ is maintained constraints, R is null
function setR!(o::StrEstimator{T}, parent::StrBootTest{T}, R₁::AbstractMatrix{T}, R::Union{UniformScaling{Bool},AbstractMatrix{T}}=Matrix{T}(undef,0,0)) where {T,E}
  o.restricted = nrows(R₁) > 0
	if o.restricted > 0
		singular, invR₁R₁ = invsymsingcheck(R₁ * R₁')
		singular && throw(ErrorException("Null hypothesis or model constraints are inconsistent or redundant."))
	  o.R₁invR₁R₁ = R₁'invR₁R₁
	  F = eigen(Symmetric(o.R₁invR₁R₁ * R₁))
	  o.R₁perp = F.vectors[:, abs.(F.values) .< 1000*eps(T)]  # eigenvectors orthogonal to span of R₁; foundation for parameterizing subspace compatible with constraints
  else
	  o.R₁invR₁R₁ = Matrix{T}(undef, parent.kZ, 0)  # and R₁perp = I
  end

  if !iszero(o.κ)
	  RR₁perp = Matrix{T}([R ; zeros(T, parent.kY₂, parent.kX₁) I])  # rows to prevent partialling out of endogenous regressors; convert sparse matrix produced by constructor to dense

	  o.restricted && (RR₁perp *= o.R₁perp)

	  F = eigen(Symmetric(RR₁perp'pinv(RR₁perp * RR₁perp')*RR₁perp)); val = abs.(F.values) .> 1000*eps(T)
	  o.Rpar   = F.vectors[:, val]  # defines attack surface for null
	  o.restricted && (o.Rpar = o.R₁perp * o.Rpar)  # fold model constraint factors into Rpar, RperpX
	  o.RRpar = R * o.Rpar
		o.RparX = o.Rpar[1:parent.kX₁,:]  # part of Rpar that refers to Y₂
		o.RparY = o.Rpar[parent.kX₁+1:end,:]  # part of Rpar that refers to Y₂
	  o.R₁invR₁R₁Y = o.R₁invR₁R₁[parent.kX₁+1:end,:]
	  o.RR₁invR₁R₁ = R * o.R₁invR₁R₁

		if !(o.isDGP && parent.WREnonARubin)  # DGP regressions will just copy the replication-regression Zperp, X2par for speed
			o.RperpX = F.vectors[:, abs.(F.values) .< 1000*eps(T)]
			o.restricted && (o.RperpX = o.R₁perp * o.RperpX)  # fold model constraint factors into Rpar, RperpX
			o.RperpX = o.RperpX[1:parent.kX₁,:]  # Zperp=Z*RperpX; though formally a multiplier on Z, it will only extract exogenous components, in X₁, since all endogenous ones will be retained
			o.RperpXperp = perp(o.RperpX)
		end
  end
	nothing
end

# stuff that can be done before r set, and depends only on exogenous variables, which are fixed throughout all bootstrap methods
function InitVarsOLS!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::AbstractMatrix{T}) where T # Rperp is for replication regression--no null imposed
  o.y₁par = parent.y₁
	o.ü₁ = Vector{T}(undef, parent.Nobs)
  H = Symmetric(parent.X₁'parent.X₁)
  o.invH = Symmetric(pinv(H))
  R₁AR₁ = iszero(nrows(o.R₁perp)) ? o.invH : Symmetric(o.R₁perp * invsym(o.R₁perp'H*o.R₁perp) * o.R₁perp')  # for DGP regression
  o.β̈₀ = R₁AR₁ * (parent.X₁'o.y₁par)
  o.∂β̈∂r = R₁AR₁ * H * o.R₁invR₁R₁ - o.R₁invR₁R₁

  o.A = iszero(nrows(Rperp)) ? o.invH : Symmetric(Rperp * invsym(Rperp'H*Rperp) * Rperp')  # for replication regression
  o.AR = o.A * parent.R'
  (parent.scorebs || parent.robust) && (o.XAR = parent.X₁ * o.AR)
	nothing
end

function InitVarsARubin!(o::StrEstimator{T}, parent::StrBootTest{T}) where T
	o.y₁par = Vector{T}(undef, parent.Nobs)
	o.ü₁    = Vector{T}(undef, parent.Nobs)

	X₂X₁ = parent.X₂'parent.X₁
  H = Symmetric([parent.X₁'parent.X₁ X₂X₁' ; X₂X₁ parent.X₂'parent.X₂])
  o.A = invsym(H)
  o.AR = o.A * parent.R'
  (parent.scorebs || parent.robust) && (o.XAR = X₁₂B(parent, parent.X₁, parent.X₂, o.AR))

  R₁AR₁ = iszero(nrows(o.R₁perp)) ? o.A : Symmetric(o.R₁perp * invsym(o.R₁perp'H*o.R₁perp) * o.R₁perp')
  o.β̈₀   = R₁AR₁ * [parent.X₁'parent.y₁ ; parent.X₂'parent.y₁]
  o.∂β̈∂r = R₁AR₁ * [parent.X₁'parent.Y₂ ; parent.X₂'parent.Y₂]
	nothing
end

function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::AbstractMatrix{T}...) where T
  !isempty(Rperp) && (o.Rperp = Rperp[1])

	if o.isDGP && parent.WREnonARubin
		o.Zperp = parent.Repl.Zperp
		o.Xpar₁ = parent.Repl.Xpar₁
		isdefined(parent.Repl, :X₁) && (o.X₁ = parent.Repl.X₁)
		isdefined(parent.Repl, :X₂) && (o.X₂ = parent.Repl.X₂)
		isdefined(parent.Repl, :Y₂) && (o.Y₂ = parent.Repl.Y₂)
		isdefined(parent.Repl, :y₁) && (o.y₁ = parent.Repl.y₁)
		o.invZperpZperp = parent.Repl.invZperpZperp
		o.XX = parent.Repl.XX
		o.kX = parent.Repl.kX
		o.invXX = parent.Repl.invXX
		o.XY₂ = parent.Repl.XY₂
		o.Y₂Y₂ = parent.Repl.Y₂Y₂
		o.S✻Y₂Y₂ = parent.Repl.S✻Y₂Y₂
		o.Y₂y₁ = parent.Repl.Y₂y₁
		o.X₂y₁ = parent.Repl.X₂y₁
		o.X₁y₁ = parent.Repl.X₁y₁
		o.y₁y₁ = parent.Repl.y₁y₁
		o.invZperpZperpZperpX₁ = parent.Repl.invZperpZperpZperpX₁
		o.invZperpZperpZperpX₂ = parent.Repl.invZperpZperpZperpX₂
		o.invZperpZperpZperpy₁ = parent.Repl.invZperpZperpZperpy₁
		o.invZperpZperpZperpY₂ = parent.Repl.invZperpZperpZperpY₂
		o.S✻⋂X₁Y₂ = parent.Repl.S✻⋂X₁Y₂
		o.S✻⋂X₂Y₂ = parent.Repl.S✻⋂X₂Y₂
		o.S✻⋂ZperpY₂ = parent.Repl.S✻⋂ZperpY₂
		o.S✻Y₂y₁ = parent.Repl.S✻Y₂y₁
	else
		o.kZperp = ncols(o.RperpX)
		o.Zperp = parent.X₁ * o.RperpX
		o.S✻⋂ZperpZperp = panelcross(o.Zperp, o.Zperp, parent.info✻⋂)
		o.invZperpZperp = iszero(ncols(o.RperpX)) ? Symmetric(Matrix{T}(undef,0,0)) : invsym(sumpanelcross(o.S✻⋂ZperpZperp))

		o.Xpar₁ = parent.X₁ * o.RperpXperp  # X∥ := [Xpar₁ X₂]
		S✻⋂X₁Zperp = panelcross(o.Xpar₁, o.Zperp, parent.info✻⋂)
		ZperpX₁ = sumpanelcross(S✻⋂X₁Zperp)'
		o.invZperpZperpZperpX₁ = o.invZperpZperp * ZperpX₁

		S✻⋂X₂Zperp = panelcross(parent.X₂, o.Zperp, parent.info✻⋂)
		ZperpX₂ = sumpanelcross(S✻⋂X₂Zperp)'
		o.ZperpX = [ZperpX₁ ZperpX₂]
		o.S✻⋂XZperp = [S✻⋂X₁Zperp; S✻⋂X₂Zperp]

		o.invZperpZperpZperpX₂ = o.invZperpZperp * ZperpX₂
		o.invZperpZperpZperpX = [o.invZperpZperpZperpX₁ o.invZperpZperpZperpX₂]

		if parent.NFE>0 && (parent.liml || !isone(parent.κ) || parent.bootstrapt)  ||  parent.granular && parent.robust && parent.bootstrapt  ||  !o.liml && !isempty(Rperp)
			o.X₁ = o.Xpar₁ - o.Zperp * o.invZperpZperpZperpX₁  # shrink and FWL-process X₁; do it as an O(N) operation because it can be so size-reducing
			o.X₂ = o.Zperp * o.invZperpZperpZperpX₂; o.X₂ .= parent.X₂ .- o.X₂  # FWL-process X₂
		end

		S✻⋂X₁X₁ = panelcross(o.Xpar₁, o.Xpar₁, parent.info✻⋂)  # S⋂(X₂:*M_(Z⟂)X₁∥)
		S✻⋂X₂X₁ = panelcross(parent.X₂, o.Xpar₁, parent.info✻⋂)
		S✻⋂X₂X₂ = panelcross(parent.X₂, parent.X₂, parent.info✻⋂)
		X₂X₁ = sumpanelcross(S✻⋂X₂X₁) - ZperpX₂'o.invZperpZperp * ZperpX₁
		X₁X₁ = Symmetric(sumpanelcross(S✻⋂X₁X₁)) - Symmetric(ZperpX₁'o.invZperpZperp * ZperpX₁)
		X₂X₂ = Symmetric(sumpanelcross(S✻⋂X₂X₂)) - Symmetric(ZperpX₂'o.invZperpZperp * ZperpX₂)
		o.S✻⋂XX = [[S✻⋂X₁X₁ S✻⋂X₂X₁'] ; [S✻⋂X₂X₁ S✻⋂X₂X₂]]  # [a b; c d] syntax would call hvcat() to concatenate horizontally along dim 2 rather than 3
		o.XX = Symmetric([X₁X₁ X₂X₁' ; X₂X₁ X₂X₂])
		o.kX = ncols(o.XX)
		o.invXX = invsym(o.XX)

		o.S✻⋂ZperpY₂ = panelcross(o.Zperp, parent.Y₂, parent.info✻⋂)
		ZperpY₂ = sumpanelcross(o.S✻⋂ZperpY₂)
		o.invZperpZperpZperpY₂ = o.invZperpZperp * ZperpY₂
		((parent.NFE>0 && (parent.liml || !isone(parent.κ) || parent.bootstrapt)) || (parent.robust && parent.bootstrapt && parent.granular)) &&
			(o.Y₂ = parent.Y₂ - o.Zperp * o.invZperpZperpZperpY₂)
		o.S✻⋂Zperpy₁ = panelcross(o.Zperp, parent.y₁, parent.info✻⋂)
		Zperpy₁ = sumpanelcross(o.S✻⋂Zperpy₁)
		o.invZperpZperpZperpy₁ = o.invZperpZperp * Zperpy₁
		((parent.NFE>0 && (parent.liml || !isone(parent.κ) || parent.bootstrapt)) || (parent.scorebs || parent.robust && parent.bootstrapt && parent.granular)) &&
			(o.y₁ = parent.y₁ - o.Zperp * o.invZperpZperpZperpy₁)

	  o.S✻⋂X₁Y₂ = panelcross(o.Xpar₁, parent.Y₂, parent.info✻⋂)
	  o.S✻⋂X₂Y₂ = panelcross(parent.X₂, parent.Y₂, parent.info✻⋂)
		o.S✻⋂XY₂ = [o.S✻⋂X₁Y₂; o.S✻⋂X₂Y₂]
		o.XY₂ = sumpanelcross(o.S✻⋂XY₂) - o.invZperpZperpZperpX'ZperpY₂
		o.S✻Y₂y₁ = panelcross(parent.Y₂, parent.y₁, parent.info✻)
	  o.Y₂y₁ = sumpanelcross(o.S✻Y₂y₁) - ZperpY₂'o.invZperpZperpZperpy₁
		o.S✻Y₂Y₂ = panelcross(parent.Y₂, parent.Y₂, parent.info✻)
		o.Y₂Y₂ = Symmetric(sumpanelcross(o.S✻Y₂Y₂)) - Symmetric(ZperpY₂'o.invZperpZperpZperpY₂)
		S✻⋂X₂y₁ = panelcross(parent.X₂, parent.y₁, parent.info✻⋂)
		o.X₂y₁ = reshape(sumpanelcross(S✻⋂X₂y₁), :) - ZperpX₂'o.invZperpZperpZperpy₁
		S✻⋂X₁y₁ = panelcross(o.Xpar₁, parent.y₁, parent.info✻⋂)
		o.S✻⋂Xy₁ = [S✻⋂X₁y₁; S✻⋂X₂y₁]
		o.X₁y₁ = reshape(sumpanelcross(S✻⋂X₁y₁), :) - ZperpX₁'o.invZperpZperpZperpy₁
		o.S✻y₁y₁ = reshape(panelcross(parent.y₁, parent.y₁, parent.info✻), :)
		o.y₁y₁ = sum(o.S✻y₁y₁) - Zperpy₁'o.invZperpZperpZperpy₁
	end

  o.Z = X₁₂B(parent, parent.X₁, parent.Y₂, o.Rpar)  # Z∥

	X₁par = parent.X₁ * o.RparX  # XXX expressible as a linear combination of Xpar₁??
	S✻⋂X₁Zpar = panelcross(o.Xpar₁  , X₁par, parent.info✻⋂) + o.S✻⋂X₁Y₂ * o.RparY
	S✻⋂X₂Zpar = panelcross(parent.X₂, X₁par, parent.info✻⋂) + o.S✻⋂X₂Y₂ * o.RparY
	o.S✻⋂XZpar = [S✻⋂X₁Zpar; S✻⋂X₂Zpar]
	X₁Zpar = sumpanelcross(S✻⋂X₁Zpar)
	X₂Zpar = sumpanelcross(S✻⋂X₂Zpar)
	S✻⋂ZperpX₁par = panelcross(o.Zperp, X₁par, parent.info✻⋂)
	ZperpX₁par = sumpanelcross(S✻⋂ZperpX₁par)
	o.S✻⋂ZperpZpar = S✻⋂ZperpX₁par + o.S✻⋂ZperpY₂ * o.RparY
  ZperpZpar = ZperpX₁par + sumpanelcross(o.S✻⋂ZperpY₂) * o.RparY
	o.invZperpZperpZperpZpar = o.invZperpZperp * ZperpZpar

	S✻X₁pary₁ = panelcross(X₁par, parent.y₁, parent.info✻)
	o.S✻Zpary₁ = S✻X₁pary₁ + o.RparY' * o.S✻Y₂y₁
	o.Zy₁ = sumpanelcross(S✻X₁pary₁) + o.RparY' * o.Y₂y₁ - ZperpX₁par'o.invZperpZperpZperpy₁

	S✻X₁parY₂ = panelcross(X₁par, parent.Y₂, parent.info✻)
	o.XZ = [X₁Zpar - o.invZperpZperpZperpX₁'ZperpZpar ; X₂Zpar - o.invZperpZperpZperpX₂'ZperpZpar]
  o.S✻ZparY₂ = S✻X₁parY₂ + o.RparY' * o.S✻Y₂Y₂
  tmp = S✻X₁parY₂ * o.RparY; o.S✻ZparZpar = panelcross(X₁par, X₁par, parent.info✻) + tmp + tmp' + o.RparY' * o.S✻Y₂Y₂ * o.RparY
  o.ZZ = Symmetric(sumpanelcross(o.S✻ZparZpar)) - Symmetric(ZperpZpar'o.invZperpZperpZperpZpar)

  o.invXXXZ = o.invXX * o.XZ
  o.ZXinvXXXZ = o.XZ'o.invXXXZ  # symmetric but converting to Symmetric() hampers type inference in the one place it's used

  if o.restricted
		_ZR₁ = X₁₂B(parent, parent.X₁, parent.Y₂, o.R₁invR₁R₁)
		S✻⋂X₁ZR₁ = panelcross(o.Xpar₁, _ZR₁, parent.info✻⋂)
		S✻⋂X₂ZR₁ = panelcross(parent.X₂, _ZR₁, parent.info✻⋂)
		o.S✻⋂XZR₁ = [S✻⋂X₁ZR₁ ; S✻⋂X₂ZR₁]
		o.S✻⋂ZperpZR₁ = panelcross(o.Zperp, _ZR₁, parent.info✻⋂)
		o.ZperpZR₁ = sumpanelcross(o.S✻⋂ZperpZR₁)
		o.invZperpZperpZperpZR₁ = o.invZperpZperp * o.ZperpZR₁
		o.ZR₁ = _ZR₁ - o.Zperp * o.invZperpZperpZperpZR₁
	  o.X₁ZR₁    = sumpanelcross(S✻⋂X₁ZR₁) - o.invZperpZperpZperpX₁'o.ZperpZR₁
	  o.X₂ZR₁    = sumpanelcross(S✻⋂X₂ZR₁) - o.invZperpZperpZperpX₂'o.ZperpZR₁
		o.S✻ZR₁Z   = panelcross(_ZR₁, o.Z, parent.info✻)
	  o.ZR₁Z     = sumpanelcross(o.S✻ZR₁Z) - o.ZperpZR₁'o.invZperpZperp * ZperpZpar
		o.S✻ZR₁Y₂  = panelcross(_ZR₁, parent.Y₂, parent.info✻)
	  o.ZR₁Y₂    = sumpanelcross(o.S✻ZR₁Y₂) - o.ZperpZR₁'o.invZperpZperpZperpY₂
		o.S✻ZR₁y₁  = panelcross(_ZR₁, parent.y₁, parent.info✻)
	  o.twoZR₁y₁ = 2 * (sumpanelcross(o.S✻ZR₁y₁) - o.ZperpZR₁'o.invZperpZperpZperpy₁)
		o.S✻ZR₁ZR₁ = panelcross(_ZR₁, _ZR₁, parent.info✻)
	  o.ZR₁ZR₁   = Symmetric(sumpanelcross(o.S✻ZR₁ZR₁)) - Symmetric(o.ZperpZR₁'o.invZperpZperp * o.ZperpZR₁)
  else
	  o.Y₂y₁par    = o.Y₂y₁
	  o.X₂y₁par    = o.X₂y₁
	  o.X₁y₁par    = o.X₁y₁
	  o.Zy₁par     = o.Zy₁
	  o.y₁pary₁par = o.y₁y₁
	  o.Xy₁par     = [o.X₁y₁ ; o.X₂y₁]
  end
	(parent.scorebs || parent.robust && parent.bootstrapt && parent.granular) &&
		(o.y₁par = copy(o.y₁))

	if o.isDGP
		if parent.scorebs
			o.ü₁ = Vector{T}(undef, parent.Nobs)
		elseif parent.robust && parent.bootstrapt && parent.granular
			o.Ü₂ = Matrix{T}(undef, parent.Nobs, parent.kY₂)
			o.u⃛₁ = Vector{T}(undef, parent.Nobs)
		end
	end

	o.Z .-= o.Zperp * o.invZperpZperpZperpZpar

  o.V =  o.invXX * o.XZ  # in 2SLS case, estimator is (V' XZ)^-1 * (V'Xy₁). Also used in k-class and liml robust VCV by Stata convention

  if o.isDGP
		parent.WREnonARubin && (o.ZY₂ = sumpanelcross(o.S✻ZparY₂) - ZperpZpar'o.invZperpZperpZperpY₂)
		o.H_2SLS = Symmetric(o.V'o.XZ)  # Hessian
		(o.liml || !isone(o.κ)) && (o.H_2SLSmZZ = o.H_2SLS - o.ZZ)
		!o.liml && MakeH!(o, parent, !isempty(Rperp)) # DGP is liml except possibly when getting confidence peak for A-R plot; but liml=0 when exactly id'd, for then κ=1 always and Hessian doesn't depend on r₁ and can be computed now
	else
	  o.kZ = ncols(o.Rpar)
	  o.Yendog = [true colsum(o.RparY .!= zero(T)).!=0]  # columns of Y = [y₁par Zpar] that are endogenous (normally all)
	end
	nothing
end

# do most of estimation; for liml r₁ must be passed now in order to solve eigenvalue problem involving it
# inconsistency: for replication regression of Anderson-Rubin, r₁ refers to the *null*, not the maintained constraints, because that's what affects the endogenous variables
# For OLS, compute β̈₀ (β̈ when r=0) and ∂β̈∂r without knowing r₁, for efficiency
# For WRE, should only be called once for the replication regressions, since for them r₁ is the unchanging model constraints
function EstimateOLS!(o::StrEstimator, r₁::AbstractVector)
  o.β̈ = o.β̈₀ - o.∂β̈∂r * r₁
	nothing
end

function EstimateARubin!(o::StrEstimator{T}, parent::StrBootTest{T}, r₁::AbstractVector) where T
  o.β̈ = o.β̈₀ - o.∂β̈∂r * r₁
  o.y₁par .= parent.y₁ .- parent.Y₂ * r₁
	nothing
end

function MakeH!(o::StrEstimator{T}, parent::StrBootTest{T}, makeXAR::Bool=false) where T
  H = isone(o.κ) ? o.H_2SLS : o.ZZ + o.κ * o.H_2SLSmZZ
  o.invH = invsym(H)
  if makeXAR  # for replication regression in score bootstrap of IV/GMM
	  o.A = ncols(o.Rperp)>0 ? Symmetric(o.Rperp * invsym(o.Rperp'H*o.Rperp) * o.Rperp') : o.invH
	  o.AR = o.A * (o.Rpar'parent.R')
	  o.XAR = X₁₂B(parent, o.X₁, o.X₂, o.V * o.AR)
  end
	nothing
end

function EstimateIV!(o::StrEstimator{T}, parent::StrBootTest{T}, r₁::AbstractVector) where T
  if o.restricted
    o.y₁pary₁par = o.y₁y₁ - (o.twoZR₁y₁'r₁)[1] + r₁'o.ZR₁ZR₁ * r₁
	  o.Y₂y₁par = o.Y₂y₁  - o.ZR₁Y₂'r₁
	  o.X₂y₁par = o.X₂y₁ - o.X₂ZR₁ * r₁
	  o.X₁y₁par = o.X₁y₁ - o.X₁ZR₁ * r₁
	  o.Zy₁par  = o.Zy₁ -  o.ZR₁Z'r₁
	  o.Xy₁par  = [o.X₁y₁par ; o.X₂y₁par]
	  (parent.scorebs || parent.robust && parent.bootstrapt && parent.granular) &&
			(o.y₁par .= o.y₁ .- o.ZR₁ * r₁)
  end

  o.invXXXy₁par = o.invXX * o.Xy₁par
  o.ZXinvXXXy₁par = o.XZ'o.invXXXy₁par
  o.YY   = Symmetric([[o.y₁pary₁par          ] o.Zy₁par'        ; o.Zy₁par        Matrix(o.ZZ)])
  o.YPXY = Symmetric([[o.invXXXy₁par'o.Xy₁par] o.ZXinvXXXy₁par' ; o.ZXinvXXXy₁par  o.ZXinvXXXZ])

  if o.isDGP
	  if o.liml
	    o.κ = 1/(1 - real(eigvals(invsym(o.YY) * o.YPXY)[1]))  # like Fast & Wild (81), but more stable, at least in Mata
	    !iszero(o.fuller) && (o.κ -= o.fuller / (parent._Nobs - parent.kX))
	    MakeH!(o, parent)
	  end

	  o.β̈ = o.invH * (isone(o.κ) ? o.ZXinvXXXy₁par : o.κ * (o.ZXinvXXXy₁par - o.Zy₁par) + o.Zy₁par)
	  o.t₁Y = o.R₁invR₁R₁Y * r₁
  elseif parent.WREnonARubin  # if not score bootstrap of IV/GMM
	  o.Rt₁ = o.RR₁invR₁R₁ * r₁
  end
	nothing
end

@inline function MakeResidualsOLSARubin!(o::StrEstimator{T}, parent::StrBootTest{T}) where T
  o.ü₁ .= o.y₁par .- X₁₂B(parent, parent.X₁, parent.X₂, o.β̈)
	nothing
end

function MakeResidualsIV!(o::StrEstimator{T}, parent::StrBootTest{T}) where T
  if parent.scorebs
		o.ü₁ .= o.y₁par .- o.Z * o.β̈
	else
		_β = [1 ; -o.β̈]
	  uu = _β'o.YY * _β

	  Xu = o.Xy₁par - o.XZ * o.β̈  # after DGP regression, compute Y₂ residuals by regressing Y₂ on X while controlling for y₁ residuals, done through FWL
	  negXuinvuu = Xu / -uu
	  o.Π̂ = invsym(o.XX + negXuinvuu * Xu') * (negXuinvuu * (o.Y₂y₁par - o.ZY₂'o.β̈)' + o.XY₂)
		o.γ̈ = o.RparY * o.β̈ + o.t₁Y
	  if parent.robust && parent.bootstrapt && parent.granular
			o.Ü₂ .= o.Y₂ .- X₁₂B(parent, o.X₁, o.X₂, o.Π̂)
			o.u⃛₁ .= o.y₁par .- o.Z * o.β̈ .+ o.Ü₂ * o.γ̈
		end
  end
	nothing
end

# non-WRE stuff that only depends on r in A-R case, for test stat denominators in replication regressions
# since the non-AR OLS code never creates an object for replication regresssions, in that case this is called on the DGP regression object
# depends on results of Estimate() only when doing OLS-style bootstrap on an overidentified IV/GMM regression--score bootstrap or A-R. Then κ from DGP liml affects Hessian, H.
function InitTestDenoms!(o::StrEstimator{T}, parent::StrBootTest{T}) where T
  if parent.bootstrapt && (parent.scorebs || parent.robust)
	  (parent.granular || parent.purerobust) && (o.WXAR = o.XAR)  # XXX simplify

	  if parent.robust && parent.NFE>0 && !(parent.FEboot || parent.scorebs) && parent.granular < parent.NErrClustCombs  # make first factor of second term of (64) for c=⋂ (c=1)
	    !isdefined(o, :WXAR) && (o.WXAR = o.XAR)  # XXX simplify
	    o.CT_XAR = [crosstabFEt(parent, view(o.WXAR,:,d), parent.info⋂) for d ∈ 1:parent.dof]
	  end
  end
	nothing
end


text:
######################################################################
# range.jl
#
# range only sensing
######################################################################

type RangeOnly <: Sensor
	noise_sigma::Float64		# noise std deviation

	# for discretized measurements
	num_bins::Int
	bin_range::UnitRange{Int64}

	function RangeOnly(noise_sigma::Real, num_bins::Int)
		new( float(noise_sigma), num_bins, 0:(num_bins-1) )
	end
	RangeOnly(noise_sigma::Real) = RangeOnly(noise_sigma, 1000)
end

function observe(m::SearchDomain, ro::RangeOnly, p::Pose)
	dx = (p[1] - m.theta[1])
	dy = (p[2] - m.theta[2])
	return sqrt(dx*dx + dy*dy) + ro.noise_sigma*randn()
end

# just returns a density. used in particle filtering
function O(ro::RangeOnly, theta::LocTuple, p::Pose, o::Float64)
	dx = (p[1] - theta[1])
	dy = (p[2] - theta[2])
	o_diff = sqrt(dx*dx + dy*dy) - o

	d = Normal(0, ro.noise_sigma)
	p = pdf(d, o_diff)
	return p
end

# just returns a density. used in particle filtering
# TODO: make this be correct
function O(ro::RangeOnly, theta::LocTuple, p::Pose, o::Int)
	dx = (p[1] - theta[1])
	dy = (p[2] - theta[2])
	o_diff = sqrt(dx*dx + dy*dy) - o

	d = Normal(0, ro.noise_sigma)
	p = pdf(d, o_diff)
	return p
end

# here, num_bins isn't important; we just bin to nearest integer
obs2bin(o::Float64, s::RangeOnly) = round(Int, o, RoundDown)


text:
using InteractionNetworkModels, BenchmarkTools, Plots, Distributions, StatsPlots
using Distances, Distances, ProgressMeter

d = Hamming()
V = 20
mode = rand(Bool,V,V)
model = CER(
    mode, 0.05
)
n = 10
data = draw_sample(model, n)

plot(map(x->d(x,mode), data))

G_prior = CER(mode, 0.5)
α_prior = ScaledBeta(1,10,0,0.5)
plot(α_prior)

posterior = CerPosterior(data, G_prior, α_prior)

mcmc = CerPosteriorMcmc(τ=0.001, ε=0.01)

@time out = mcmc(
    posterior,
    desired_samples=10000, burn_in=0, lag=1,
    G_init=rand_perturb(mode, 0.001),
    α_init=0.02
)
acceptance_prob(mcmc)
plot(out, mode)

plot(log_posterior_prob(out))

hline!([model.α], subplot=2)

# Verify alpha is unbiased
function sim(
    true_model::CER,
    mcmc::CerPosteriorMcmc,
    n_train::Int,
    n_reps::Int
    )

    # α_true = true_model.α
    α_ests = Float64[]
    @showprogress for i in 1:n_reps
        data = draw_sample(true_model, n_train)
        G_prior = CER(true_model.mode, 0.5)
        α_prior = ScaledBeta(1,10,0,0.5)
        posterior = CerPosterior(data, G_prior, α_prior)
        out = mcmc(
            posterior,
            G_init=rand_perturb(true_model.mode, 0.001),
            α_init=true_model.α+0.1,
            loading_bar=false
        )
        push!(α_ests, mean(out.α_sample))
    end
    return α_ests
end

mcmc = CerPosteriorMcmc(
    τ=0.001, ε=0.01,
    desired_samples=250,
    burn_in=5_000,
)
ests = sim(model, mcmc, 20, 500)
boxplot(ests)
hline!([model.α])

max_stars_repo_path: Examples/Graph Models/CER_test.jl

max_issues_repo_path: Examples/Graph Models/CER_test.jl

max_forks_repo_path: Examples/Graph Models/CER_test.jl


text:
## Exercise 3-3
## A function object is a value you can assign to a variable or pass as an argument. For example, dotwice is a function that takes a function object as an argument and calls it twice:

## function dotwice(f)
##    f()
##    f()
## end

## Here’s an example that uses dotwice to call a function named printspam twice.

## function printspam()
##    println("spam")
## end

## dotwice(printspam)

## 1. Type this example into a script and test it.
println("Ans 1:")
function dotwice(func)
    func()
    func()
end

function printspam()
    println("Spam")
end

dotwice(printspam)

## 2. Modify dotwice so that it takes two arguments, a function object and a value, and calls the function twice, passing the value as an argument.
println("Ans 2:")
function dotwice(func, arg)
    func(arg)
    func(arg)
end

# extra
function dotwice(func, args...)
    func(args...)
    func(args...)
end

dotwice(println, "dotwice")

## 3. Copy the definition of printtwice from earlier in this chapter to your script.
println("Ans 3:")
function printtwice(bruce)
    println(bruce)
    println(bruce)
end

## 4. Use the modified version of dotwice to call printtwice twice, passing "spam" as an argument.
println("Ans 4:")
dotwice(printtwice, "spam")

## 5. Define a new function called dofour that takes a function object and a value and calls the function four times, passing the value as a parameter. There should be only two statements in the body of this function, not four.
println("Ans 5:")
function dofour(func, arg)
    dotwice(func, arg)
    dotwice(func, arg)
end

# extra
function dofour(func, args...)
    dotwice(func, args...)
    dotwice(func, args...)
end

dofour(printtwice, "spam")

println("End.")


text:
"""
    iht!(v::IHTVariables, x::BEDFile, y, k)

If used with a `BEDFile` object `x`, then the temporary arrays `b0`, `Xb`, `Xb0`, and `r` housed in the `IHTVariables` object `v` are all initialized as `SharedArray`s of the proper dimensions.
The additional optional arguments are:

- `pids`, a vector of process IDs. Defaults to `procs(x)`.
"""
function iht!{T <: Float}(
    v     :: IHTVariables{T},
    x     :: BEDFile{T},
    y     :: DenseVector{T},
    k     :: Int;
    pids  :: Vector{Int} = procs(x),
    iter  :: Int = 1,
    nstep :: Int = 50,
)
    # compute indices of nonzeroes in beta
    _iht_indices(v, k)

    # if support has not changed between iterations,
    # then xk and gk are the same as well
    # avoid extracting and computing them if they have not changed
    # one exception: we should always extract columns on first iteration
    if !isequal(v.idx, v.idx0) || iter < 2
        decompress_genotypes!(v.xk, x, v.idx)
    end

    # store relevant components of gradient
    fill_perm!(v.gk, v.df, v.idx)  # gk = g[v.idx]

    # now compute subset of x*g
    BLAS.gemv!('N', one(T), v.xk, v.gk, zero(T), v.xgk)

    # warn if xgk only contains zeros
    all(v.xgk .== zero(T)) && warn("Entire active set has values equal to 0")

    # compute step size
    mu = (sumabs2(v.gk) / sumabs2(v.xgk)) :: T
#    mu = _iht_stepsize(v, k) :: T

    # notify problems with step size
    isfinite(mu) || throw(error("Step size is not finite, is active set all zero?"))
    mu <= eps(typeof(mu))  && warn("Step size $(mu) is below machine precision, algorithm may not converge correctly")

    # compute gradient step
    _iht_gradstep(v, mu, k)

    # update xb
    #PLINK.A_mul_B!(v.xb, x, v.b, v.idx, k, pids=pids)
    PLINK.A_mul_B!(v.xb, x, v.b, v.idx, k)

    # calculate omega
    omega_top, omega_bot = _iht_omega(v)

    # backtrack until mu sits below omega and support stabilizes
    mu_step = 0
    while _iht_backtrack(v, omega_top, omega_bot, mu, mu_step, nstep)

        # stephalving
        mu /= 2

        # warn if mu falls below machine epsilon
        mu <= eps(typeof(mu)) && warn("Step size equals zero, algorithm may not converge correctly")

        # recompute gradient step
        copy!(v.b, v.b0)
        _iht_gradstep(v, mu, k)

        # recompute xb
        #PLINK.A_mul_B!(v.xb, x, v.b, v.idx, k, pids=pids)
        PLINK.A_mul_B!(v.xb, x, v.b, v.idx, k)

        # calculate omega
        omega_top, omega_bot = _iht_omega(v)

        # increment the counter
        mu_step += 1
    end

    return mu::T, mu_step::Int
end

"""
    L0_reg(x::BEDFile, y, k)

If used with a `BEDFile` object `x`, then the temporary floating point arrays are all initialized as `SharedArray`s of the proper dimensions.
The additional optional arguments are:

- `pids`, a vector of process IDs. Defaults to `procs()`.
- `mask_n`, an `Int` vector used as a bitmask for crossvalidation purposes. Defaults to a vector of ones.
"""
function L0_reg{T <: Float, V <: DenseVector}(
    x        :: BEDFile{T},
    y        :: V,
    k        :: Int;
    pids     :: Vector{Int} = procs(x),
    v        :: IHTVariables{T, V} = IHTVariables(x, y, k),
    mask_n   :: Vector{Int} = ones(Int, size(y)),
    tol      :: T     = convert(T, 1e-4),
    max_iter :: Int   = 100,
    max_step :: Int   = 50,
    quiet    :: Bool  = true
)

    # start timer
    tic()

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

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

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

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

    # initialize booleans
    converged = false             # scaled_norm < tol?

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

    # calculate the gradient
    At_mul_B!(v.df, x, v.r, mask_n, pids=pids)

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

    # main loop
    for mm_iter = 1:max_iter

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

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

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

            # stop timer
            mm_time = toq()

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

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

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

        # the IHT kernel gives us an updated x*b
        # use it to recompute residuals and gradient
#        update_r_grad!(v, x, y, pids=pids)
        difference!(v.r, y, v.xb)
        mask!(v.r, mask_n, 0, zero(T), n=n)

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

        # update loss, objective, and gradient
        next_loss = sumabs2(sdata(v.r)) / 2

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

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

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

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

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

            # stop time
            mm_time = toq()

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

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

        # algorithm is unconverged at this point.
        # if algorithm is in feasible set, then rho should not be changing
        # check descent property in that case
        # if rho is not changing but objective increases, then abort
        if next_loss > loss + tol
            !quiet && print_descent_error(mm_iter, loss, next_loss)
            throw(ErrorException("Descent failure!"))
        end
    end # end main loop
end # end function

"""
    iht_path(x::BEDFile, y, path)

If used with a `BEDFile` object `x`, then the temporary arrays are all initialized as `SharedArray`s of the proper dimensions.
The additional optional arguments are:

- `pids`, a vector of process IDs. Defaults to `procs(x)`.
- `mask_n`, an `Int` vector used as a bitmask for crossvalidation purposes. Defaults to a vector of ones.
"""
function iht_path{T <: Float}(
    x        :: BEDFile{T},
    y        :: DenseVector{T},
    path     :: DenseVector{Int};
    pids     :: Vector{Int} = procs(x),
    mask_n   :: Vector{Int} = ones(Int, size(y)),
    tol      :: T    = convert(T, 1e-4),
    max_iter :: Int  = 100,
    max_step :: Int  = 50,
    quiet    :: Bool = true
)

    # size of problem?
    n,p = size(x)

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

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

    # preallocate temporary arrays
    v = IHTVariables(x, y, 1)

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

        # model size?
        q = path[i]

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

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

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

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

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

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

    # return a sparsified copy of the models
    return betas
end

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

If used with a `BEDFile` object `x`, then the additional optional arguments are:

- `pids`, a vector of process IDs. Defaults to `procs(x)`.
"""
function one_fold{T <: Float}(
    x        :: BEDFile{T},
    y        :: DenseVector{T},
    path     :: DenseVector{Int},
    folds    :: DenseVector{Int},
    fold     :: Int;
    pids     :: Vector{Int} = procs(x),
    tol      :: T    = convert(T, 1e-4),
    max_iter :: Int  = 100,
    max_step :: Int  = 50,
    quiet    :: Bool = true
)
    # dimensions of problem
    n,p = size(x)

    # make vector of indices for folds
    test_idx = folds .== fold
    test_size = sum(test_idx)

    # train_idx is the vector that indexes the TRAINING set
    train_idx = !test_idx
    mask_n    = convert(Vector{Int}, train_idx)
    mask_test = convert(Vector{Int}, test_idx)

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

    # tidy up
    #gc()

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

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

    # allocate the arrays for the test set
    xb = SharedArray(T, (n,), init = S -> S[localindexes(S)] = zero(T), pids=pids) :: SharedVector{T}
    b  = SharedArray(T, (p,), init = S -> S[localindexes(S)] = zero(T), pids=pids) :: SharedVector{T}
    r  = SharedArray(T, (n,), init = S -> S[localindexes(S)] = zero(T), pids=pids) :: SharedVector{T}

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

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

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

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

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

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

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

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

    return myerrors :: Vector{T}
end

function pfold(
    T          :: Type,
    xfile      :: String,
    xtfile     :: String,
    x2file     :: String,
    yfile      :: String,
    meanfile   :: String,
    precfile   :: String,
    path       :: DenseVector{Int},
    folds      :: DenseVector{Int},
    pids       :: Vector{Int},
    q          :: Int;
    max_iter   :: Int  = 100,
    max_step   :: Int  = 50,
    quiet      :: Bool = true,
    header     :: Bool = false
)

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

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

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

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

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

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

        # loop over all workers
        for worker in pids

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

                # asynchronously distribute tasks
                @async begin
                    while true

                        # grab next fold
                        current_fold = nextidx()

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

                        # report distribution of fold to worker and device
                        quiet || print_with_color(:blue, "Computing fold $current_fold on worker $worker.\n\n")
                        r = remotecall_fetch(worker) do
                            processes = [worker]
                            x = BEDFile(T, xfile, x2file, meanfile, precfile, pids=processes, header=header)
                            y = SharedArray(abspath(yfile), T, (x.geno.n,), pids=processes) :: SharedVector{T}

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

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

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

function pfold(
    T        :: Type,
    xfile    :: String,
    x2file   :: String,
    yfile    :: String,
    path     :: DenseVector{Int},
    folds    :: DenseVector{Int},
    pids     :: Vector{Int},
    q        :: Int;
    max_iter :: Int  = 100,
    max_step :: Int  = 50,
    quiet    :: Bool = true,
    header   :: Bool = false
)

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

    # do not allow crossvalidation with fewer than 3 folds
    q > 2 || throw(ArgumentError("Number of folds q = $q must be at least 3."))

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

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

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

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

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

        # loop over all workers
        for worker in pids

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

                # asynchronously distribute tasks
                @async begin
                    while true

                        # grab next fold
                        current_fold = nextidx()

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

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

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

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

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

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

"""
    cv_iht(xfile, xtfile, x2file, yfile, meanfile, precfile, [q=max(3, min(CPU_CORES,5)), path=collect(1:min(p,20)), folds=cv_get_folds(n,q), pids=procs()])

This variant of `cv_iht()` performs `q`-fold crossvalidation with a `BEDFile` object loaded by `xfile`, `xtfile`, and `x2file`,
with column means stored in `meanfile` and column precisions stored in `precfile`.
The continuous response is stored in `yfile` with data particioned by the `Int` vector `folds`.
The folds are distributed across the processes given by `pids`.
The dimensions `n` and `p` are inferred from BIM and FAM files corresponding to the BED file path `xpath`.
"""
function cv_iht(
    T        :: Type,
    xfile    :: String,
    xtfile   :: String,
    x2file   :: String,
    yfile    :: String,
    meanfile :: String,
    precfile :: String;
    q        :: Int = cv_get_num_folds(3,5),
    path     :: DenseVector{Int} = begin
           # find p from the corresponding BIM file, then make path
            bimfile = xfile[1:(endof(xfile)-3)] * "bim"
            p       = countlines(bimfile)
            collect(1:min(20,p))
            end,
    folds    :: DenseVector{Int} = begin
           # find n from the corresponding FAM file, then make folds
            famfile = xfile[1:(endof(xfile)-3)] * "fam"
            n       = countlines(famfile)
            cv_get_folds(n, q)
            end,
    pids     :: Vector{Int} = procs(),
    tol      :: Float = convert(T, 1e-4),
    max_iter :: Int   = 100,
    max_step :: Int   = 50,
    quiet    :: Bool  = true,
    header   :: Bool  = false
)

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

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

    # compute folds in parallel
    mses = pfold(T, xfile, xtfile, x2file, yfile, meanfile, precfile, path, folds, pids, q, max_iter=max_iter, max_step=max_step, quiet=quiet, header=header)

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

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

    # recompute ideal model
    # first load data on *all* processes
    x = BEDFile(T, xfile, xtfile, x2file, meanfile, precfile, header=header, pids=pids) :: BEDFile{T}
    y = SharedArray(abspath(yfile), T, (x.geno.n,), pids=pids) :: SharedVector{T}

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

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

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

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

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

# encodes default type FLoat64 for previous function
### 22 Sep 2016: Julia v0.5 warns that this conflicts with cv_iht for GPUs
### since this is no longer the default interface for cv_iht with CPUs,
### then it is commented out here
#cv_iht(xfile::String, xtfile::String, x2file::String, yfile::String, meanfile::String, precfile::String; q::Int = max(3, min(CPU_CORES, 5)), path::DenseVector{Int} = begin bimfile=xfile[1:(endof(xfile)-3)] * "bim"; p=countlines(bimfile); collect(1:min(20,p)) end, folds::DenseVector{Int} = begin famfile=xfile[1:(endof(xfile)-3)] * "fam"; n=countlines(famfile); cv_get_folds(n, q) end, pids::DenseVector{Int}=procs(), tol::Float64=1e-4, max_iter::Int=100, max_step::Int=50, quiet::Bool=true, header::Bool=false) = cv_iht(Float64, xfile, xtfile, x2file, yfile, meanfile, precfile, path=path, folds=folds, q=q, pids=pids, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet, header=header)

"""
    cv_iht(T::Type, xfile, x2file, yfile, [q=max(3, min(CPU_CORES,5)), path=collect(1:min(p,20)), folds=cv_get_folds(n,q), pids=procs()])

An abbreviated call to `cv_iht` that calculates means, precs, and transpose on the fly.
"""
function cv_iht(
    T        :: Type,
    xfile    :: String,
    x2file   :: String,
    yfile    :: String;
    q        :: Int = cv_get_num_folds(3,5),
    path     :: DenseVector{Int} = begin
           # find p from the corresponding BIM file, then make path
            bimfile = xfile[1:(endof(xfile)-3)] * "bim"
            p       = countlines(bimfile)
            collect(1:min(20,p))
            end,
    folds    :: DenseVector{Int} = begin
           # find n from the corresponding FAM file, then make folds
            famfile = xfile[1:(endof(xfile)-3)] * "fam"
            n       = countlines(famfile)
            cv_get_folds(n, q)
            end,
    pids     :: Vector{Int} = procs(),
    tol      :: Float = convert(T, 1e-4),
    max_iter :: Int   = 100,
    max_step :: Int   = 50,
    quiet    :: Bool  = true,
    header   :: Bool  = false
)

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

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

    # compute folds in parallel
    mses = pfold(T, xfile, x2file, yfile, path, folds, pids, q, max_iter=max_iter, max_step=max_step, quiet=quiet, header=header)

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

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

    # recompute ideal model
    # first load data on *all* processes
    x = BEDFile(T, xfile, x2file, header=header, pids=pids) :: BEDFile{T}
    y = SharedArray(abspath(yfile), T, (x.geno.n,), pids=pids) :: SharedVector{T}

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

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

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

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

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

"""
    cv_iht(xfile, x2file, yfile)

The default call to `cv_iht`. Here `xfile` points to the PLINK BED file stored on disk, `x2file` points to the nongenetic covariates stored in a delimited file, and `yfile` points to the response variable stored in a **binary** file.

Important optional arguments and defaults include:

- `q`, the number of crossvalidation folds. Defaults to `max(3, min(CPU_CORES,5))`
- `path`, an `Int` vector that contains the model sizes to test. Defaults to `collect(1:min(p,20))`, where `p` is the number of genetic predictors read from the PLINK BIM file.
- `folds`, an `Int` vector that specifies the fold structure. Defaults to `cv_get_folds(n,q)`, where `n` is the number of cases read from the PLINK FAM file.
- `pids`, an `Int` vector of process IDs. Defaults to `procs()`.
"""
cv_iht(xfile::String, x2file::String, yfile::String; q::Int = cv_get_num_folds(3,5), path::DenseVector{Int} = begin bimfile=xfile[1:(endof(xfile)-3)] * "bim"; p=countlines(bimfile); collect(1:min(20,p)) end, folds::DenseVector{Int} = begin famfile=xfile[1:(endof(xfile)-3)] * "fam"; n=countlines(famfile); cv_get_folds(n, q) end, pids::Vector{Int}=procs(), tol::Float64=1e-4, max_iter::Int=100, max_step::Int=50, quiet::Bool=true, header::Bool=false) = cv_iht(Float64, xfile, x2file, yfile, path=path, folds=folds, q=q, pids=pids, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet, header=header)