text:
abstract type AbstractPrimaryEdgeCutModel end

@inline (mdl::AbstractPrimaryEdgeCutModel)(n_vertices, n_self_pins, n_cut_pins, k) = mdl(n_vertices, n_self_pins, n_cut_pins)

struct AffinePrimaryEdgeCutModel{Tv} <: AbstractPrimaryEdgeCutModel
    α::Tv
    β_vertex::Tv
    β_self_pin::Tv
    β_cut_pin::Tv
end

function AffinePrimaryEdgeCutModel(; α = false, β_vertex = false, β_self_pin = false, β_cut_pin = false)
    AffinePrimaryEdgeCutModel(promote(α, β_vertex, β_self_pin, β_cut_pin)...)
end

@inline cost_type(::Type{AffinePrimaryEdgeCutModel{Tv}}) where {Tv} = Tv

(mdl::AffinePrimaryEdgeCutModel)(n_vertices, n_self_pins, n_cut_pins, k) = mdl.α + n_vertices * mdl.β_vertex + n_self_pins * mdl.β_self_pin + n_cut_pins * mdl.β_cut_pin

struct PrimaryEdgeCutOracle{Ti, Lcp, Mdl} <: AbstractOracleCost{Mdl}
    pos::Vector{Ti}
    lcp::Lcp
    mdl::Mdl
end

oracle_model(ocl::PrimaryEdgeCutOracle) = ocl.mdl

function bound_stripe(A::SparseMatrixCSC, K, mdl::AffinePrimaryEdgeCutModel)
    @inbounds begin
        m, n = size(A)
        N = nnz(A)
        @assert mdl.β_vertex >= 0
        @assert mdl.β_self_pin >= 0
        @assert mdl.β_cut_pin >= 0
        c_hi = mdl.α + mdl.β_vertex * n + max(mdl.β_self_pin, mdl.β_cut_pin) * N
        c_lo = mdl.α + fld(mdl.β_vertex * n + min(mdl.β_self_pin, mdl.β_cut_pin) * N, K)
        return (c_lo, c_hi)
    end
end

function oracle_stripe(hint::AbstractHint, mdl::AbstractPrimaryEdgeCutModel, A::SparseMatrixCSC, Π; net=nothing, adj_A=nothing, kwargs...)
    @inbounds begin
        m, n = size(A)
        pos = A.colptr
        args, Ap = partwise(A, convert(MapPartition, Π))
        lcp = partwisecount!(hint, args..., pincount(hint, Ap, kwargs...); kwargs...)
        return PrimaryEdgeCutOracle(pos, lcp, mdl)
    end
end

@inline function (cst::PrimaryEdgeCutOracle{Ti, Mdl})(j::Ti, j′::Ti, k) where {Ti, Mdl}
    @inbounds begin
        w = cst.pos[j′] - cst.pos[j]
        l = cst.lcp(j, j′, k)
        return cst.mdl(j′ - j, l, w - l, k)
    end
end

@inline function (stp::Step{Ocl})(_j, _j′, _k) where {Ti, Mdl, Ocl <: PrimaryEdgeCutOracle{Ti, Mdl}}
    @inbounds begin
        cst = stp.ocl
        j = destep(_j)
        j′ = destep(_j′)
        k = destep(_k)
        w = cst.pos[j′] - cst.pos[j]
        l = Step(cst.lcp)(_j, _j′, _k)
        return cst.mdl(j′ - j, l, w - l, k)
    end
end


text:
"""
    getOverlap(MPSA, MPSB)

calculates the between to MPS encoded in MPSA and MPSB: < MPSB | MPSA >

# Arguments
- MPSA: vector containing the ket of the overlap
- MPSB: vector containing the bra of the overlap

return:

"""
function getOverlap(MPSA::Vector{<:Any}, MPSB::Vector{<:Any})
	Length = Int(size(MPSA)[1])

        @tensor T_multi[α_bra, α_ket, β_bra, β_ket] := (MPSA[1][α_ket, d, β_ket] *
							conj(MPSB[1][α_bra, d, β_bra]))

        for site = 2:Length
        	#println("add ",site)
        	@tensor begin

        		E_local[a_in',a_in,a_out',a_out] := (MPSA[site][a_in, d, a_out] *
        						     conj(MPSB[site][a_in', d, a_out']))

        		T_multi[α_bra, α_ket, β_bra, β_ket] := (T_multi[α_bra, α_ket, γ_bra, γ_ket] *
        							E_local[γ_bra, γ_ket, β_bra, β_ket])
        	end
        end

	return T_multi #(α_bra, α_ket, β_bra, β_ket)
end

"""
    getExpValue(MPSvec, MPOvec)

calculates the expectation value of an MPO : < MPSvec | MPOvec | MPSvec >

# Arguments
- MPSvec: vector of M's
- MPOvec: vector of MPO's

return:

"""
function getExpValue(MPSvec::Vector{<:Any}, MPOvec::Vector{MPO})
	OpString = stringMPO(MPOvec)
	expValue = 0
	for (i, x) in enumerate(OpString[1])
		test = applyTM_OP(MPSvec, MPSvec, x, MPOvec, Id(size(MPSvec[1])[1]))
		expValue += tr(test)
	end

	return expValue
end

"""
    singleSiteExpValue(MPSvec, C, Op, site; can_left = true)

calculates the expectation value of a single site Operator Op given
the local MPS-matrix MPSvec[site] and singular-values between the next(previous) bond
depending on the left(right) canonical condition of M specified by can_left.

# Arguments
- MPSvec: containing the local MPS-matrices
- C: singular values missing to bring M into mixed canonical form
- Op: local single site operator
- site: site on which Op should act
- can_left = true : is MPSvec[site] in left or right can. form.

return:

	expectation values of < MPSvec | Op_site | MPSvec >

"""
singleSiteExpValue(MPSvec::Vector{<:Any}, C::AbstractArray{<:Number},
		   Op::AbstractArray{<:Number}, site::Int; can_left=true) = singleSiteExpValue(MPSvec[site], C, Op, can_left = can_left)

function singleSiteExpValue(M::AbstractArray{<:T}, C::AbstractArray{<:T}, Op::AbstractArray{<:T}; can_left = true) where T<:Number
	if can_left
		@tensor tmp[β_ket,α_ket,d] := C[γ, β_ket] * M[α_ket, d, γ]
	else
		@tensor tmp[β_ket,α_ket,d] := C[α_ket, γ] * M[γ, d, β_ket]
	end

	@tensor tmp[d_ket,d_bra] := tmp[β,α,d_ket] * conj(tmp[β,α,d_bra])
	@tensor expValue[] := Op[d_bra,d_ket] * tmp[d_ket,d_bra]

	#expVlaue = 2.0
	return expValue
end

function getCoef(MPSvec)

	# to do get Coefficient vector
end

ext: jl

max_stars_repo_name: Gawatz/finiteMPS

max_issues_repo_name: Gawatz/finiteMPS

max_forks_repo_name: Gawatz/finiteMPS


text:
using JuLIPMaterials
using Test
using JuLIP

# use the same DataDep as JuLIP
test_pots = JuLIP.Deps.fetch_test_pots()

tests = [
    "test_greensfunctions.jl",
    "test_gfcorrectors3D.jl",
    "test_dislocations.jl",
    "test_elasticity.jl",
    "test_fracture.jl",
   "test_forceconstants.jl"
   # "test_cauchyborn1.jl",
]

println("≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡")
println("   JuLIPMaterials: Start Tests   ")
println("≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡")
@testset "JuLIPMaterials" begin
   for t in tests
      ts = @testset "$(t[6:end-3])" begin
         include(t)
      end
      println(ts)
   end
end


text:
nonparametricdistr_insts(ext::VecSpecExt) = [
    [(d, ext)
    for d in vcat(10, 500, 1000, 2500, 5000:5000:30000)] # includes compile run
    ]

insts = OrderedDict()
insts["nat"] = (nothing, vcat(nonparametricdistr_insts.([
    VecNegRtdet(),
    VecNegLog(),
    VecNegSqrt(),
    VecNegEntropy(),
    ])...))
insts["vecext"] = (nothing, vcat(nonparametricdistr_insts.([
    VecNegRtdetEFExp(),
    VecNegLogEF(),
    VecNegSqrtEF(),
    VecNegEntropyEF(),
    ])...))
return (NonparametricDistrJuMP, insts)


text:
using Test
using ArnoldiMethod: reflector!, restore_arnoldi!, Reflector, Arnoldi, reinitialize!, iterate_arnoldi!
using LinearAlgebra
using Random

@testset "Reflector $T" for T in (Float64, ComplexF64)
    n = 20
    x = rand(T, n)
    z = copy(x)
    τ = reflector!(z, n)
    z[n] = 1
    y = x .- τ * dot(z, x) .* z

    @inferred reflector!(copy(x), n)
    @test norm(y[1:n-1]) ≤ 10eps()
    @test abs(real(y[n])) ≈ norm(x)
    @test abs(imag(y[n])) ≤ eps()
    @test 1 ≤ real(τ) ≤ 2
    @test abs(τ - 1) ≤ 1
end

@testset "Trivial reflector $T" for T in (Float64, ComplexF64)
    x = [T(0), T(0), T(5)]
    z = copy(x)
    τ = reflector!(z, 3)
    @test iszero(τ)
end

@testset "rmul! $T" for T in (Float64, ComplexF64)
    A = rand(T, 4, 4)
    B = copy(A)

    # Implicit representation
    G = Reflector{T}(4)
    rand!(G.vec)
    τ = reflector!(G, 4)

    # Matrix representation
    z = [G.vec[1:3]; 1]
    H = I - τ * z * z'

    # Apply
    rmul!(A, G, 1, 4)

    @test A ≈ B * H'
end

@testset "lmul! $T" for T in (Float64, ComplexF64)
    A = rand(T, 4, 4)
    B = copy(A)

    # Implicit representation
    G = Reflector{T}(4)
    rand!(G.vec)
    τ = reflector!(G, 4)

    # Matrix representation
    z = [G.vec[1:3]; 1]
    H = I - τ * z * z'

    # Apply
    lmul!(G, A, 1, 4)

    @test A ≈ H * B
end

@testset "Restore Hessenberg matrix $T" for T in (Float64, ComplexF64)
    n, k = 10, 6
    Q = qr(randn(T, k, k)).Q * Matrix(1.0I, k, k)

    A = rand(T, n, n)
    arn = Arnoldi{T}(n, k)
    reinitialize!(arn)
    iterate_arnoldi!(A, arn, 1:k)

    # Do an orthogonal change of basis.

    H = copy(arn.H)
    H[1:k,1:k] .= Q' * H[1:k,1:k] * Q

    # Restore the relation
    restore_arnoldi!(H, 1, k, Q, Reflector{T}(k))

    # Only now do the change of basis in V with accumulated Q.
    W = [arn.V[:, 1:k] * Q arn.V[:, k+1]]
    @test norm(A * W[:, 1:k] - W[:, 1:k+1] * H) < 1e-14
end


text:
const HTML_minimal_template = mt"""
<!DOCTYPE html>

  <meta charset="utf-8">
  <meta name="viewport" content=\"{{{:viewport_content}}}\">

  <title>{{{:title}}}</title>

"""

HTML_minimal_dict = Dict(
  :lang => "en",
  :viewport_content => "width=device-width, initial-scale=1",
  :title => "Untitled",
)
# :bulma => "<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/bulma/0.7.5/css/bulma.min.css">"
# :fontawesome => ""

HTML_dict = DefaultDict(0, HTML_minimal_dict)

#=
res = Mustache.render(HTML_minimal_template, HTML_dict)
print(res)

<!DOCTYPE html>

  <meta charset="utf-8">
  <meta name="viewport" content="width=device-width, initial-scale=1">
  <title>Untitled</title>

using EzXML
root = root(parsehtml(res))
head = firstelement(root)
@assert head.name == "head"
body = lastelement(root)
@assert body.name == "body"

using Hyperscript, EzXML
julia> @tags div h1

julia> const entry = div.entry
WARNING: redefining constant entry

julia> entry(h1("An Important Announcement"))

An Important Announcement

julia> s = string(entry(h1("An Important Announcement")))
"

An Important Announcement

"

julia> p=parsehtml(s)
EzXML.Document(EzXML.Node(<HTML_DOCUMENT_NODE@0x000000000136f800>))

julia> b=lastelement(root(p))
EzXML.Node(<ELEMENT_NODE[body]@0x000000003c9b45f0>)

julia> prettyprint(b)

An Important Announcement

julia> string(b)
"

An Important Announcement

"

julia> for x in eachelement(root(parsehtml(res)))
          for y in eachelement(x)
            println(y.name)
           end
           end
meta
meta
title

julia> for x in eachelement(root(parsehtml(res)))
          for y in eachelement(x)
            for a in (attributes(y))
               println(a.name,": ",a.content)
           end
           end
           end

charset: utf-8
name: viewport
content: width=device-width, initial-scale=1

=#


text:
using Test, YaoBlocks, YaoArrayRegister
using YaoBase, BitBasis

struct MockedTag{BT,N} <: TagBlock{BT,N}
    content::BT

    MockedTag(x::BT) where {N,BT<:AbstractBlock{N}} = new{BT,N}(x)
end

@test nqubits(MockedTag(X)) == nqubits(X)
@test nqubits(MockedTag(kron(X, Y))) == nqubits(kron(X, Y))

@test getiparams(MockedTag(phase(0.1))) == ()
@test getiparams(MockedTag(cache(phase(0.1)))) == ()
@test getiparams(MockedTag(Rx(0.1))) == ()

@test parameters(MockedTag(phase(0.1))) == [0.1]
@test parameters(MockedTag(cache(phase(0.1)))) == [0.1]
@test parameters(MockedTag(Rx(0.1))) == [0.1]

@test occupied_locs(MockedTag(chain(3, put(1 => X), put(3 => X)))) == occupied_locs(chain(
    3,
    put(1 => X),
    put(3 => X),
))

@testset "scale apply" begin
    # factor
    @test factor(-X) == -1
    @test factor(-2 * X) == -2

    # type
    @test -X isa Scale{<:Val}
    @test -im * X isa Scale{<:Val}
    @test -2 * X isa Scale{<:Number}

    # apply!
    reg = rand_state(1)
    @test apply!(copy(reg), -X) ≈ -apply!(copy(reg), X)
    @test apply!(copy(reg), (-im * Y)') ≈ apply!(copy(reg), im * Y)
    @test apply!(copy(reg), (-2im * Y)') ≈ apply!(copy(reg), 2im * Y)

    # mat
    @test mat(-X) ≈ -mat(X)
    @test mat((-im * Y)') ≈ (-im * mat(Y))'
    @test mat((-2im * Y)') ≈ (-2im * mat(Y))'

    # other
    @test chsubblocks(-X, Y) == -Y
    @test cache_key(Scale(Val(-1), X)) == cache_key(Scale(-1, X))
    @test copy(-X) == -X
    @test copy(-X) isa Scale{<:Val}
end

@testset "properties" begin
    xg = put(100, 3 => X)
    yg = put(100, 3 => Y)
    yg = put(100, 3 => Y)
    @test_throws Exception !ishermitian(im * xg)
    @test ishermitian(3 * xg)
    @test_throws Exception !isreflexive(im * xg)
    @test isreflexive(-1 * xg)
    @test isunitary(im * xg)
    @test_throws Exception !isunitary(2 * xg)
    @test !iscommute(2 * xg, 2 * yg)
    @test iscommute(2 * xg, 2 * xg)
end

@testset "daggered" begin
    dg = Daggered(ConstGate.T)
    @test dg isa Daggered
    @test apply!(product_state(bit"1"), dg) ≈ apply!(product_state(bit"1"), ConstGate.Tdag)
    @test chsubblocks(dg, X) == Daggered(X)
end


text:
# Allocator for persistent memory regions.
#
# Each allocation initializes a new file for reading and writing.

# Global constant mapping pointers to the allocated size and backing file.
const PtrMap = Dict{Ptr{Nothing}, Tuple{Int, String}}()

# The directory where pmem backed files are created.
const FILE_DIR = Ref(get(ENV, "JULIA_PA_DIR", "/mnt/public"))

# Constants copied from libpmem.h
const PMEM_FILE_CREATE = (1 << 0)

"""
    alloc(bytes::Integer, [file]) -> (Ptr{Nothing}, String)

Allocate a persistent memory region of size `bytes`.
Returns a tuple of the pointer to the newly allocated region and the backing file name.

If `file` is not provided, a temporary filename will be used.
"""
function alloc(bytes::Integer, file = basename(tempname()))
    # Create the whole file path.
    filepath = joinpath(FILE_DIR[], file)

    # Call into libpmem
    ptr = ccall(
        (:pmem_map_file, libpmem),
        Ptr{Nothing},
        (Cstring, Csize_t, Cint, Base.Cmode_t, Csize_t, Cint),
        filepath,
        bytes,
        PMEM_FILE_CREATE,
        0o666,
        Ptr{Cvoid}(),
        Ptr{Cvoid}(),
    )

    # Save the pointer mapping and return.
    PtrMap[ptr] = (bytes, filepath)
    return (ptr, filepath)
end

"""
    free(ptr::Ptr{Nothing}, delete_file=false)

Free the persistent memory backing `ptr`.
If `delete_file = true`, also delete the backing file form the file system.
"""
function free(ptr::Ptr{Nothing}, delete_file = false)
    # Get the mapped size and file
    bytes, filepath = PtrMap[ptr]

    ccall(
        (:pmem_unmap, libpmem),
        Cvoid,
        (Ptr{Cvoid}, Csize_t),
        ptr,
        bytes,
    )

    # Clean up this pointer mapping.
    delete!(PtrMap, ptr)

    # If requested, clean up the backing file
    delete_file && rm(filepath)
    return nothing
end


text:
module Ping

using Discord

function on_message_create(c::Client, event::MessageCreate)
    if event.message.content == "ping"
       reply(c, event.message, "Pong!")
    elseif event.message.content == "pong"
       reply(c, event.message, "Ping!")
    end
end

function main()
    c = Client(ENV["DISCORD_TOKEN"])
    open(c)
    add_handler!(c, MessageCreate, on_message_create; tag=:ping)
    return c
end

end

if abspath(PROGRAM_FILE) == @__FILE__
    c = Ping.main()
    wait(c)
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((5, 1, 2, true, true, true),),
    ((5, 1, 2, false, true, true),),
    ((5, 2, 1, false, true, true),),
    ((5, 1, 2, true, false, true), (default_tol_relax = 100,)),
    ((5, 1, 2, true, true, false),),
    ((:iris, 2, true, true, true),),
    ]
insts["fast"] = [
    ((50, 1, 4, true, true, true),),
    ((50, 1, 10, true, true, true),),
    ((50, 1, 50, true, true, true),),
    ((100, 1, 250, true, true, true),),
    ((50, 2, 2, true, true, true),),
    ((200, 2, 20, true, true, true),),
    ((50, 2, 2, false, true, true),),
    ((50, 2, 2, true, false, true),),
    ((50, 2, 2, true, true, false),),
    ((500, 3, 14, true, true, true),),
    ((20, 4, 3, false, true, false),),
    ((20, 4, 3, true, true, true),),
    ((100, 8, 2, true, true, true),),
    ((100, 8, 2, false, true, true),),
    ((100, 8, 2, true, false, true),),
    ((100, 8, 2, true, true, false),),
    ((250, 4, 6, true, true, true),),
    ((250, 4, 6, false, true, true),),
    ((250, 4, 6, true, false, true),),
    ((200, 32, 2, true, true, true),),
    ((:iris, 4, true, true, true),),
    ((:iris, 5, true, true, true),),
    ((:iris, 6, true, true, true),),
    ((:iris, 4, false, true, true),),
    ((:iris, 4, true, false, true),),
    ((:iris, 4, true, true, false),),
    ((:cancer, 4, true, true, true),),
    ((:cancer, 4, false, true, true),),
    ((:cancer, 4, true, false, true),),
    ]
return (DensityEstNative, insts)


text:
using RecipesBase
using Random: GLOBAL_RNG
export SamplingPlot

# function xs_ys_barline(walls, ys)
#     @argcheck length(walls) == length(ys) + 1
#     xs_ret = Float64[]
#     ys_ret = Float64[]
#
#     push!(xs_ret, walls[1])
#     push!(ys_ret, 0)
#
#     for i in 1:length(walls)-1
#         x_left = walls[i]
#         x_right = walls[i+1]
#         y = ys[i]
#         push!(xs_ret, x_left)
#         push!(ys_ret, y)
#         push!(xs_ret, x_right)
#         push!(ys_ret, y)
#     end
#     push!(xs_ret, walls[end])
#     push!(ys_ret, 0)
#
#     xs_ret, ys_ret
# end
function xs_ys_barline(walls, ys)
    @argcheck length(walls) == length(ys) + 1
    xs_ret = Float64[]
    ys_ret = Float64[]

    push!(xs_ret, walls[1])
    push!(ys_ret, 0)

    for i in 1:length(walls)-1
        x_left = walls[i]
        x_right = walls[i+1]
        y = ys[i]
        push!(xs_ret, x_left)
        push!(ys_ret, y)
        push!(xs_ret, x_right)
        push!(ys_ret, y)
        push!(xs_ret, x_right)
        push!(ys_ret, 0)
    end

    xs_ret, ys_ret
end

function spacings(xs)
    xs[2:end] - xs[1:end-1]
end

@recipe function plot(grid::VegasGrid{1})
    # seriestype --> :vline
    # walls = first(grid.boundaries)
    xs = first(grid.boundaries)
    if get(plotattributes, :seriestype, nothing) == :vline
        xs
    else
        ys = 1 ./ spacings(xs)
        xs_ys_barline(xs, ys)
    end
end

@qstruct SamplingPlot(f, domain; neval::Int=100, axis::Int=1)

@recipe function plot(p::SamplingPlot)
    seriestype --> :scatter
    axis = p.axis
    f = p.f; dom = p.domain
    s = draw(GLOBAL_RNG, dom)
    pt = s.position
    x = pt[p.axis]
    y = f(pt) * s.weight
    xs = typeof(x)[]
    ys = typeof(y)[]
    for _ in 1:p.neval
        s = draw(GLOBAL_RNG, dom)
        pt =  s.position
        y = f(pt) * s.weight
        x = pt[p.axis]
        push!(xs, x)
        push!(ys, y)
    end
    xs, ys
end


text:
# From David Goldberg's blog post
# https://tech.ebayinc.com/engineering/fast-approximate-logarithms-part-iii-the-formulas/
@inline function fastlog2(x::Float32)::Float32
    # (x-1)*(a*(x-1) + b)/((x-1) + c) (line 8 of table 2)
    a =   0.338953f0
    b =   2.198599f0
    c =   1.523692f0
    #
    # Assume IEEE representation, which is sgn(1):exp(8):frac(23)
    # representing (1+frac)*2^(exp-127)  Call 1+frac the significand
    #

    # get exponent
    ux1i = reinterpret(UInt32, x)
    exp = (ux1i & 0x7F800000) >> 23
    # actual exponent is exp-127, will subtract 127 later

    greater = ux1i & 0x00400000  # true if signif > 1.5
    if greater !== 0x00000000
        ux2i = (ux1i & 0x007FFFFF) | 0x3f000000
        signif = reinterpret(Float32, ux2i)
        fexp = exp - 126f0    # 126 instead of 127 compensates for division by 2
        signif = signif - 1.0f0
    else
        ux2i = (ux1i & 0x007FFFFF) | 0x3f800000
        signif = reinterpret(Float32, ux2i)
        fexp = exp - 127f0
        signif = signif - 1.0f0
    end
    lg2 = fexp + signif*(a*signif + b)/(signif + c)
    return lg2
end

# Translated from OpenLibm but less accurate because I forced the tableau to be
# Float32, whereas OpenLibm uses Float64
#
# https://github.com/JuliaMath/openlibm/blob/cca41bc1abd01804afa4862bbd2c79cc9803171a/src/s_exp2f.c
const EXP2FT = (Float32(0x1.6a09e667f3bcdp-1),
                Float32(0x1.7a11473eb0187p-1),
                Float32(0x1.8ace5422aa0dbp-1),
                Float32(0x1.9c49182a3f090p-1),
                Float32(0x1.ae89f995ad3adp-1),
                Float32(0x1.c199bdd85529cp-1),
                Float32(0x1.d5818dcfba487p-1),
                Float32(0x1.ea4afa2a490dap-1),
                Float32(0x1.0000000000000p+0),
                Float32(0x1.0b5586cf9890fp+0),
                Float32(0x1.172b83c7d517bp+0),
                Float32(0x1.2387a6e756238p+0),
                Float32(0x1.306fe0a31b715p+0),
                Float32(0x1.3dea64c123422p+0),
                Float32(0x1.4bfdad5362a27p+0),
                Float32(0x1.5ab07dd485429p+0))
@inline function _exp2(x::Float32)
    TBLBITS = UInt32(4)
    TBLSIZE = UInt32(1 << TBLBITS)

    redux = Float32(0x1.8p23f) / TBLSIZE
    P1    = Float32(0x1.62e430p-1f)
    P2    = Float32(0x1.ebfbe0p-3f)
    P3    = Float32(0x1.c6b348p-5f)
    P4    = Float32(0x1.3b2c9cp-7f)

    # Reduce x, computing z, i0, and k.
    t::Float32 = x + redux
    i0 = reinterpret(UInt32, t)
    i0 += TBLSIZE ÷ UInt32(2)
    k::UInt32 = unsafe_trunc(UInt32, (i0 >> TBLBITS) << 20)
    i0 &= TBLSIZE - UInt32(1)
    t -= redux
    z = x - t
    twopk = Float32(reinterpret(Float64, UInt64(0x3ff00000 + k) << 32))

    # Compute r = exp2(y) = exp2ft[i0] * p(z).
    tv = EXP2FT[i0+UInt32(1)]
    u = tv * z
    tv = tv + u * (P1 + z * P2) + u * (z * z) * (P3 + z * P4)

    # Scale by 2**(k>>20)
    return tv * twopk
end

"""
    fastpow(x::Real, y::Real) -> Float32
"""
@inline function fastpow(x::Real, y::Real)
    if iszero(x)
        return 0f0
    elseif isinf(x) && isinf(y)
        return Float32(Inf)
    else
        return _exp2(convert(Float32,y) * fastlog2(convert(Float32, x)))
    end
end
@inline fastpow(x, y) = x^y


text:
using Test
using MPI
using Random

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
sz = MPI.Comm_size(comm)
filename = MPI.bcast(tempname(), 0, comm)

MPI.Barrier(comm)

# Assume `data` is a subarray of a (3sz, 2sz, 5) global array that should be
# written as a contiguous block to a single file. Note that `data` is a
# discontiguous view of the global array.
data = ArrayType(rand(Float64, 3, 2, 5))
offset = (3rank, 2rank, 0)
etype = MPI.Datatype(eltype(data))
filetype = MPI.Types.create_subarray((3sz, 2sz, 5), size(data), offset, etype)
MPI.Types.commit!(filetype)

fh = MPI.File.open(comm, filename, read=true, write=true, create=true)

# Collective write
disp = 0
MPI.File.set_view!(fh, disp, etype, filetype)
MPI.File.write_all(fh, data)

# Noncollective write
disp = sizeof(data) * sz
MPI.File.set_view!(fh, disp, etype, filetype)
MPI.File.write(fh, data)

MPI.File.sync(fh)

# Noncollective read
data_read = similar(data)
disp = 0
MPI.File.set_view!(fh, disp, etype, filetype)
MPI.File.read!(fh, data_read)
@test data_read == data

# Collective read
disp = sizeof(data) * sz
fill!(data_read, 0)
MPI.File.set_view!(fh, disp, etype, filetype)
MPI.File.read_all!(fh, data_read)
@test data_read == data

close(fh)


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
	  invR₁R₁ = invsym(R₁ * R₁')
	  all(iszero.(diag(invR₁R₁))) && 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([R
		                  zeros(T, parent.kY₂, parent.kX₁) I])  # rows to prevent partialling out of endogenous regressors
	  nrows(R₁)>0 && (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]  # perp and par of RR₁perp
    o.RperpX = F.vectors[:, abs.(F.values) .< 1000*eps(T)]

	  if nrows(R₁) > 0  # fold model constraint factors into Rpar, RperpX
	    o.Rpar   = o.R₁perp * o.Rpar
	    o.RperpX = o.R₁perp * o.RperpX
	  end

	  o.RRpar = R * o.Rpar
	  o.RperpX = o.RperpX[1:parent.kX₁,:]  # Zperp=X₁*RperpX
		o.RperpXperp = perp(o.RperpX)
		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₁
  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₁
  H = symcross(parent.X₁, parent.wt)
  o.invH = Symmetric(inv(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₁ * crossvec(parent.X₁, parent.wt, 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
  X₂X₁ = cross(parent.X₂, parent.wt, parent.X₁)
  H = Symmetric([symcross(parent.X₁, parent.wt) X₂X₁' ; X₂X₁ symcross(parent.X₂, parent.wt)])
  o.A = inv(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₁ * [crossvec(parent.X₁, parent.wt, parent.y₁) ; crossvec(parent.X₂, parent.wt, parent.y₁)]
  o.∂β̂∂r = R₁AR₁ * [cross(parent.X₁, parent.wt, parent.Y₂) ; cross(parent.X₂, parent.wt, parent.Y₂)]
	nothing
end

@inline sumpanelsum(X::Array{T,3} where T) = dropdims(sum(X, dims=2), dims=2)

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

  o.Zperp = parent.X₁ * o.RperpX
	o.S⋂ZperpZperp = panelcross11(parent, o.Zperp, o.Zperp, parent.info⋂)  #  XXX exploit symmetry?  XXX missing weights
  o.invZperpZperp = iszero(length(o.Zperp)) ? Symmetric(Matrix{T}(undef,0,0)) : inv(Symmetric(sumpanelsum(o.S⋂ZperpZperp)))
	#=o.X₁=# _X₁ = parent.X₁ * o.RperpXperp
	S⋂X₁Zperp = panelcross11(parent, _X₁, o.Zperp, parent.info⋂)
	ZperpX₁ = sumpanelsum(S⋂X₁Zperp)'

	invZperpZperpZperpX₁ = o.invZperpZperp * ZperpX₁
	o.X₁ = _X₁ - o.Zperp * invZperpZperpZperpX₁  # shrink and FWL-process X₁; do it as an O(N) operation because it can be so size-reducing

	S⋂X₂X₁ = panelcross11(parent, parent.X₂, _X₁, parent.info⋂)  # S⋂(X₂:*M_(Z⟂)X₁∥)
	S⋂X₂Zperp = panelcross11(parent, parent.X₂, o.Zperp, parent.info⋂)
	ZperpX₂ = sumpanelsum(S⋂X₂Zperp)'
	o.ZperpX = [ZperpX₁ ZperpX₂]

	invZperpZperpZperpX₂ = o.invZperpZperp * ZperpX₂
	o.X₂ = o.Zperp * invZperpZperpZperpX₂; o.X₂ .= parent.X₂ .- o.X₂  # FWL-process X₂
	o.invZperpZperpZperpX = [invZperpZperpZperpX₁ invZperpZperpZperpX₂]

	X₂X₁ = sumpanelsum(S⋂X₂X₁) - ZperpX₂'o.invZperpZperp * ZperpX₁  # cross(o.X₂, parent.wt, o.X₁)
  X₁X₁ = Symmetric(cross(_X₁, parent.wt, _X₁) - ZperpX₁'o.invZperpZperp * ZperpX₁)
	X₂X₂ = symcross(parent.X₂, parent.wt) - Symmetric(ZperpX₂'o.invZperpZperp * ZperpX₂)
	o.XX = Symmetric([X₁X₁ X₂X₁' ; X₂X₁ X₂X₂])
  o.kX = ncols(o.XX)
  o.invXX = invsym(o.XX)

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

	S⋂X₁Z = panelcross11(parent,       _X₁, _Z, parent.info⋂)
	S⋂X₂Z = panelcross11(parent, parent.X₂, _Z, parent.info⋂)
	S⋂XZ  = [S⋂X₁Z; S⋂X₂Z]
	X₁Z = sumpanelsum(S⋂X₁Z)
	X₂Z = sumpanelsum(S⋂X₂Z)
	S⋂X₂Zperp = panelcross11(parent, parent.X₂, o.Zperp, parent.info⋂)
	ZperpX₂ = sumpanelsum(S⋂X₂Zperp)'
	o.S⋂XZperp = [S⋂X₁Zperp; S⋂X₂Zperp]
	S⋂ZperpZ = panelcross11(parent, o.Zperp, _Z, parent.info⋂)
  o.ZperpZ = sumpanelsum(S⋂ZperpZ)
	o.invZperpZperpZperpZ = o.invZperpZperp * o.ZperpZ

#	!o.isDGP && ((o.LIML || !parent.robust || !isone(o.κ)) || parent.granular) &&
		(o.Z = _Z - o.Zperp * o.invZperpZperpZperpZ)
	o.ZperpY₂ = cross(o.Zperp, parent.wt, parent.Y₂)
	o.invZperpZperpZperpY₂ = o.invZperpZperp * o.ZperpY₂
  o.Y₂ = parent.Y₂ - o.Zperp * o.invZperpZperpZperpY₂
	o.S⋂Zperpy₁ = panelcross11(parent, o.Zperp, parent.y₁, parent.info⋂)
  Zperpy₁ = reshape(sumpanelsum(o.S⋂Zperpy₁), :)
	o.invZperpZperpZperpy₁ = o.invZperpZperp * Zperpy₁
	o.y₁ = parent.y₁ - o.Zperp * o.invZperpZperpZperpy₁

  o.X₁Y₂ = cross(_X₁, parent.wt, parent.Y₂) - ZperpX₁'o.invZperpZperp * o.ZperpY₂  # cross(o.X₁, parent.wt, o.Y₂)
  o.X₂Y₂ = cross(parent.X₂, parent.wt, parent.Y₂) .- ZperpX₂'o.invZperpZperp * o.ZperpY₂  # cross(o.X₂, parent.wt, o.Y₂)
  o.XY₂ = [o.X₁Y₂ ; o.X₂Y₂]
  o.Y₂y₁ = crossvec(parent.Y₂, parent.wt, parent.y₁) - o.ZperpY₂'o.invZperpZperpZperpy₁  # crossvec(o.Y₂, parent.wt, o.y₁)
	S⋂X₂y₁ = panelcross11(parent, parent.X₂, parent.y₁, parent.info⋂)
  o.X₂y₁ = reshape(sumpanelsum(S⋂X₂y₁), :) - ZperpX₂'o.invZperpZperpZperpy₁  # crossvec(o.X₂, parent.wt, o.y₁)
	S⋂X₁y₁ = panelcross11(parent, _X₁, parent.y₁, parent.info⋂)
  o.X₁y₁ = reshape(sumpanelsum(S⋂X₁y₁), :) - ZperpX₁'o.invZperpZperpZperpy₁  # crossvec(o.X₁, parent.wt, o.y₁)
  o.y₁y₁ = cross(parent.y₁, parent.wt, parent.y₁)[1] -  Zperpy₁'o.invZperpZperpZperpy₁  # cross(o.y₁, parent.wt, o.y₁)[1]
  o.Zy₁  = crossvec(_Z, parent.wt, parent.y₁) - o.ZperpZ'o.invZperpZperpZperpy₁  # crossvec(o.Z, parent.wt, o.y₁)
	o.XZ   = [X₁Z - ZperpX₁'o.invZperpZperpZperpZ ; X₂Z - ZperpX₂'o.invZperpZperpZperpZ]
  o.ZY₂ =  cross(_Z, parent.wt, parent.Y₂) - o.ZperpZ'o.invZperpZperp * o.ZperpY₂  # cross(o.Z, parent.wt, o.Y₂)
  o.ZZ  =  symcross(_Z, parent.wt) - Symmetric(o.ZperpZ'o.invZperpZperpZperpZ)

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

  if o.restricted
		#=o.ZR₁=# _ZR₁ = X₁₂B(parent, parent.X₁, parent.Y₂, o.R₁invR₁R₁)
		S⋂X₁_ZR₁ = panelcross11(parent, _X₁, _ZR₁, parent.info⋂)
		S⋂X₂_ZR₁ = panelcross11(parent, parent.X₂, _ZR₁, parent.info⋂)
		o.S⋂X_ZR₁ = [S⋂X₁_ZR₁; S⋂X₂_ZR₁]
		o.S⋂Zperp_ZR₁ = panelcross11(parent, o.Zperp, _ZR₁, parent.info⋂)
		o.Zperp_ZR₁ = sumpanelsum(o.S⋂Zperp_ZR₁)
		o.invZperpZperpZperpZR₁ = o.invZperpZperp * o.Zperp_ZR₁
		#=o.ZR₁ .-= =# o.ZR₁ = _ZR₁ - o.Zperp * o.invZperpZperpZperpZR₁
	  o.X₁ZR₁    = sumpanelsum(S⋂X₁_ZR₁) - ZperpX₁'o.invZperpZperp * o.Zperp_ZR₁  # cross(o.X₁, parent.wt, o.ZR₁)
	  o.X₂ZR₁    = sumpanelsum(S⋂X₂_ZR₁) - ZperpX₂'o.invZperpZperp * o.Zperp_ZR₁  # cross(o.X₂, parent.wt, o.ZR₁)
	  o.ZZR₁     = cross(       _Z, parent.wt, _ZR₁) - o.ZperpZ'o.invZperpZperp * o.Zperp_ZR₁  # cross(o.Z , parent.wt, o.ZR₁)
	  o.ZR₁Y₂    = cross(_ZR₁, parent.wt, parent.y₁) - o.Zperp_ZR₁'o.invZperpZperpZperpy₁  # cross(o.ZR₁, parent.wt, o.Y₂)
	  o.twoR₁Zy₁ = 2 * (crossvec(_ZR₁, parent.wt, parent.y₁) - o.Zperp_ZR₁'o.invZperpZperpZperpy₁)  # 2crossvec(o.ZR₁, parent.wt, o.y₁)
	  o.ZR₁ZR₁   = symcross(_ZR₁, parent.wt) - Symmetric(o.Zperp_ZR₁'o.invZperpZperp * o.Zperp_ZR₁)  # symcross(o.ZR₁, parent.wt)
o._ZR₁=_ZR₁
  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₁]
	  o.y₁par      = o.y₁
  end
o._X₁=_X₁
o._Z=_Z
  o.V =  o.invXX * o.XZ  # in 2SLS case, StrEstimator is (V' XZ)^-1 * (V'Xy₁). Also used in k-class and LIML robust VCV by Stata convention
  o.H_2SLS = Symmetric(o.V'o.XZ)  # Hessian
  (o.LIML || o.κ ≠ 1) && (o.H_2SLSmZZ = o.H_2SLS - o.ZZ)

  if o.isDGP
	  !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; #=o.RparY.type==identity ? fill(true, o.kZ) :=# mapreduce(v->v.≠0, .|, eachcol(o.RparY))]  # columns of Y = [y₁par Zpar] that are endogenous (normally all)

	  if parent.robust && parent.bootstrapt  # for WRE replication regression, prepare for CRVE
			S⋂X₁y₁ = panelcross11(parent, _X₁, parent.y₁, parent.info⋂)  # S⋂(M_(Z⟂)X₁∥ :* y₁)  XXX missing weights
			S⋂X₂y₁ = panelcross11(parent, parent.X₂, parent.y₁, parent.info⋂)  # S⋂(X₂∥:*y₁)
			o.S⋂Xy₁ = [S⋂X₁y₁; S⋂X₂y₁]

			o.S⋂YX       = Vector{Matrix{T}}(undef, o.kZ+1)
			o.S⋂PXYZperp = Vector{Matrix{T}}(undef, o.kZ+1)

			o.XinvXX = X₁₂B(parent, o.X₁, o.X₂, o.invXX); o.PXZ = X₁₂B(parent, o.X₁, o.X₂, o.invXXXZ)
			if parent.haswt
				o.PXZ .*= parent.wt
				o.XinvXX .*= parent.wt
			end

			if parent.NFE>0
				o.CT_FE⋂PY = Vector{Matrix{T}}(undef, o.kZ+1)  # XXX would array comprehension break type stability?
				@inbounds for i ∈ 1:o.kZ
					o.CT_FE⋂PY[i+1] = crosstabFEt(parent, view(o.PXZ,:,i), parent.info⋂) .* parent.invFEwt
				end
			end

			invZperpZperpZperpX = o.invZperpZperp * [ZperpX₁ ZperpX₂]
			S⋂YX = S⋂XZ - o.S⋂XZperp * o.invZperpZperpZperpZ - invZperpZperpZperpX' * (S⋂ZperpZ - o.S⋂ZperpZperp * o.invZperpZperpZperpZ)  # S⋂(M_Zperp[Z or y₁] .* P_(MZperpX)])

			FillingT₀ = o.invXXXZ' * S⋂YX
			o.FillingT₀ = Matrix{Matrix{T}}(undef, o.kZ+1, o.kZ+1)  # fixed component of groupwise term in sandwich filling  XXX would array comprehension break type stability?
			@inbounds for i ∈ 1:o.kZ, j ∈ 1:o.kZ
				o.FillingT₀[i+1,j+1] = reshape(view(FillingT₀,i,:,j),:,1)
			end

			S⋂PXYZperp = o.invXXXZ' * (o.S⋂XZperp - invZperpZperpZperpX' * o.S⋂ZperpZperp)
			@inbounds for i ∈ 1:o.kZ  # precompute various clusterwise sums
				o.S⋂PXYZperp[i+1] = S⋂PXYZperp[i,:,:]  # S⋂(P_(MZperpX) * Z .* Zperp)
				if !parent.granular
					o.S⋂YX[i+1] = S⋂YX[:,:,i]
				end
	    end

			o.FillingT₀₀ = o.invXXXZ' * (o.S⋂Xy₁ - o.ZperpX' * o.invZperpZperp * o.S⋂Zperpy₁ - (o.S⋂XZperp - o.ZperpX' * o.invZperpZperp * o.S⋂ZperpZperp) * o.invZperpZperpZperpy₁)
			o.restricted &&
				(o.∂FillingT₀∂r = o.invXXXZ' * (o.S⋂X_ZR₁ - o.S⋂XZperp * o.invZperpZperp * o.Zperp_ZR₁ - o.ZperpX' * o.invZperpZperp * o.S⋂Zperp_ZR₁))
			if !parent.granular
				S⋂Xy₁ = [S⋂X₁y₁; S⋂X₂y₁]
				o.S⋂y₁X₀ = dropdims(S⋂Xy₁ - o.S⋂XZperp * o.invZperpZperpZperpy₁ - o.ZperpX' * o.invZperpZperp * (o.S⋂Zperpy₁ - o.S⋂ZperpZperp * o.invZperpZperpZperpy₁); dims=3)
				o.restricted &&
					(o.∂S⋂y₁X∂r = o.S⋂X_ZR₁ - o.S⋂XZperp * o.invZperpZperp * o.Zperp_ZR₁ - o.ZperpX' * o.invZperpZperp * (o.S⋂Zperp_ZR₁ - o.S⋂ZperpZperp * o.invZperpZperp * o.Zperp_ZR₁))
			end
	  end
  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{T} where T, 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.twoR₁Zy₁'r₁)[1] + r₁'o.ZR₁ZR₁ * r₁
	  o.y₁par   = o.y₁ - o.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.ZZR₁ * r₁
	  o.Xy₁par  = [o.X₁y₁par ; o.X₂y₁par]
  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₁

	  if parent.robust && parent.bootstrapt  # prepare WRE replication regressions
			tmp = o.restricted ? o.FillingT₀₀ - o.∂FillingT₀∂r * r₁ : o.FillingT₀₀ # @panelsum(parent, o.PXZ, o.y₁par, parent.info⋂)

	    for i ∈ 1:o.kZ
	  	  o.FillingT₀[i+1,1] = reshape(view(tmp,i,:,1),:,1)
	    end
	    if !parent.granular
	  		o.S⋂YX[1] = o.S⋂y₁X₀
				o.restricted && (o.S⋂YX[1] .-= dropdims(o.∂S⋂y₁X∂r * r₁, dims=3))   # S⋂(M_Zperp*y₁ .* P_(MZperpX)])
			end
	  end
  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
  o.ü₁ = o.y₁par - o.Z * o.β̈

  if !parent.scorebs
	  _β = [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.Ü₂ = X₁₂B(parent, o.X₁, o.X₂, o.Π̂); o.Ü₂ .= o.Y₂ .- o.Ü₂

	  o.γ̈ = o.RparY * o.β̈ + o.t₁Y
		o.u⃛₁ = o.Ü₂ * o.γ̈; o.u⃛₁ .+= o.ü₁
  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 = vHadw(o.XAR, parent.wt))

	  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 = vHadw(o.XAR, parent.wt))
	    o.CT_XAR = [crosstabFEt(parent, view(o.WXAR,:,d), parent.info⋂) for d ∈ 1:parent.dof]
	  end
  end
	nothing
end


text:
#!/usr/bin/env julia

#
#  This is the main file for data assimilation runs.
#
#
# 1. loads up a configuration file,
# 2. obtains data from a WRF model,
# 3. construct covariate vectors
# 3. reads in observations and metadata for a list of stations,
# 4. runs the moisture model and the assimilation mechanism.
#
#

using Calendar
import Calendar.CalendarTime

using Storage
import Storage.setup_tag, Storage.spush, Storage.next_frame, Storage.flush_frame

using Stations
import Stations.Station, Stations.Observation, Stations.load_station_info,
       Stations.load_station_data, Stations.build_observation_data, Stations.register_to_grid, Stations.nearest_grid_point, Stations.obs_value, Stations.obs_station_id

using Kriging
import Kriging.trend_surface_model_kriging

using NetCDF
using WRF

using FM
import FM.FMModel, FM.advance_model, FM.kalman_update

function create_fm_variables(file_name)

    # first extract the dimensions
    nc = NetCDF.open(file_name)
    we_dim = nc.dim["west_east"]
    sn_dim = nc.dim["south_north"]
    t_dim = nc.dim["Time"]
    NetCDF.close(nc)

    var_names = ASCIIString[ "FM1", "FM10", "FM100" ]
    long_names = ["1hr fuel moisture", "10hr fuel moisture", "100hr fuel moisture"]
    for i in 1:3
        if !haskey(nc.vars, var_names[i])
            atts = {"coordinates"=>"XLONG XLAT","units"=>"kg/kg","stagger"=>"","MemoryOrder"=>"XY ","description"=>long_names[i]}
            NetCDF.nccreate(file_name, var_names[i], atts, we_dim, sn_dim, t_dim)
            println("INFO: created variable $(var_names[i])")
        else
            println("INFO: variable $(var_names[i]) already existed")
        end
    end
end

function main(args)

    # the arguments passed to the julia program do not include program name
    if length(args) != 1
        println("Usage: julia run_data_assimilation.jl cfg_file")
        exit(1)
    end

    t1 = Calendar.now()
    println("INFO: run_data_assimilation.jl started on $t1")

    ### Read configuration file and setup the system
    println("INFO: reading in config file $(args[1])")
    cfg = evalfile(args[1])

    # create the output directory if it does not exist
    println("INFO: output directory is ", cfg["output_dir"])
    !ispath(cfg["output_dir"]) && mkdir(cfg["output_dir"])

    # configure Storage mechanism
    Storage.sopen(cfg["output_dir"], "moisture_model_v2_diagnostics.txt", "frame")

    # setup Storage & output policies for interesting quantities
    setup_tag("mt", true, true)

    setup_tag("fm1_model_state_assim", false, false)
    setup_tag("fm1_model_var", false, false)
    setup_tag("fm10_model_state", false, false)
    setup_tag("fm10_model_state_assim", false, false)
    setup_tag("fm10_model_var", false, false)
    setup_tag("fm100_model_state_assim", false, false)
    setup_tag("fm100_model_var", false, false)

    setup_tag("fm10_model_na_state", false, false)
    setup_tag("fm10_model_deltas", false, false)

    setup_tag("kriging_beta", true, true)
    setup_tag("kriging_xtx_cond", false, true)
    setup_tag("kriging_field", false, false)
    setup_tag("kriging_variance", false, false)
    setup_tag("kriging_sigma2_eta", true, true)
    setup_tag("kriging_iters", false, true)
    setup_tag("kriging_subzero_s2_estimates", false, true)
    setup_tag("model_raws_mae", false, true)
    setup_tag("model_raws_mae_assim", false, true)
    setup_tag("model_na_raws_mae", false, true)

    # co-located model/model_na/kriging field/observation
    setup_tag("kriging_obs", false, false)
    setup_tag("kriging_obs_station_ids", false, false)
    setup_tag("kriging_obs_ngp", false, false)
    setup_tag("kriging_errors", true, true)

    setup_tag("kalman_gain_fm10", false, false)

    ### Load WRF output data
    t1 = Calendar.now()
    println("INFO: configuration complete, loading WRF data.")

    # read in data from the WRF output file pointed to by cfg
    w = WRF.load_wrf_data(cfg["wrf_output"], ["HGT"])

    # the terrain height need not be stored for all time points
    WRF.slice_field(w, "HGT")

    # extract WRF fields
    lat, lon = WRF.lat(w), WRF.lon(w)
    dsize = size(lat)
    wtm = WRF.times(w)
    println("INFO: WRF grid size is $(dsize[1]) x $(dsize[2]) and found $(length(wtm)) timepoints.")

    # retrieve equilibria and rain (these are already precomputed)
    Ed, Ew = WRF.field(w, "Ed"), WRF.field(w, "Ew")
    rain = WRF.field(w, "RAIN")
    hgt = WRF.field(w, "HGT")
    T = WRF.interpolated_field(w, "T2")
    P = WRF.interpolated_field(w, "PSFC")

    t2 = Calendar.now()
    println("INFO: WRF output loaded, sliced and diced [$(t2-t1)].")

    ### Load observation data from stations
    io = open(join([cfg["station_info_dir"], cfg["station_info"]], "/"), "r")
    station_ids = filter(x -> x[1] != '#', map(x -> strip(x), readlines(io)))
    close(io)

    # load each station from its info and observation files
    stations = Station[]
    for sid in station_ids
        s = load_station_info(join([cfg["station_info_dir"], string(sid, ".info")], "/"))
        load_station_data(s, join([cfg["station_data_dir"], string(sid, ".obs")], "/"))
	register_to_grid(s, lat, lon)
#        println("STATION: $(s.id), $(s.loc), ngp is $(s.ngp) with lat $(lat[s.ngp[1], s.ngp[2]]) and lon $(lon[s.ngp[1], s.ngp[2]]).")
        push!(stations, s)
    end

    # build the observation data from stations
    obs_fm10 = build_observation_data(stations, "FM")
    obs_times = keys(obs_fm10)

    t3 = Calendar.now()
    println("INFO: Station data loaded and preprocessed [$(t3 - t2)].")

    ### Initialize model

    # number of simulated fuel components
    Nf = 3

    # construct initial conditions (FIXME: can we do better here?)
    E = squeeze(0.5 * (Ed[:,:,2] + Ew[:,:,2]), 3)

    # set up parameters
    Q = diagm(cfg["Q"])
    P0 = diagm(cfg["P0"])
    mV = zeros(Float64, dsize)
    pred = zeros(Float64, dsize)
    mresV = zeros(Float64, dsize)
    mid = zeros(Int32, dsize)
    Kg = zeros(Float64, (dsize[1], dsize[2], 9))
    K = zeros(Float64, dsize)
    V = zeros(Float64, dsize)

    # prepare static & time-varying covariates
    cov_ids = cfg["covariates"]
    st_covar_map = [:lon => lon,
                    :lat => lat,
                    :elevation => hgt,
                    :constant => ones(Float64, dsize) ]
    dyn_covar_map = [:temperature => T, :pressure => P, :rain => rain]
    Xd3 = length(cov_ids) + 1
    X = zeros(Float64, (dsize[1], dsize[2], Xd3))
    Xr = zeros(Float64, (dsize[1], dsize[2], Xd3))
    for i in 2:Xd3
        cov_id = cov_ids[i-1]
        if haskey(st_covar_map, cov_id)
            println("INFO: processing static covariate $cov_id.")
            v = st_covar_map[cov_id]
            Xr[:,:,i] = v
        elseif haskey(dyn_covar_map, cov_id)
            println("INFO: found dynamic covariate $(cov_id).")
        else
            error("ERROR: unknown covariate $(cov_id) encountered, fatal.")
        end
    end
    println("INFO: there are $Xd3 covariates (including model state).")

    t1 = Calendar.now()
    println("INFO: starting simulation at $t1 ...")
    dt = (wtm[2] - wtm[1]).millis / 1000
    assim_time_win = cfg["assimilation_time_window"]
    println("INFO: time step from WRF is $dt s, assimilation time window is $assim_time_win s.")

    # construct model grid from fuel parameters
    Tk = [ 1.0, 10.0, 100.0 ]
    models = Array(FMModel, dsize)
    models_na = Array(FMModel, dsize)
    for i in 1:dsize[1]
    	for j in 1:dsize[2]
            geo_loc = (lat[i,j], lon[i,j])
	    models[i,j] = FMModel(geo_loc, Nf, E[i,j], P0, Tk)
	    models_na[i,j] = FMModel(geo_loc, Nf, E[i,j], P0, Tk)
	end
    end

    # create the new variables in the wrfout file for storing fuel moisture
    println("INFO: creating FM variables in wrfout file")
    create_fm_variables(cfg["wrf_output"])
    nc = NetCDF.open(cfg["wrf_output"], NetCDF.NC_WRITE)

    ###  Run the model and data assimilation
    for t in 2:length(wtm)
    	mt = wtm[t]
        spush("mt", mt)

        # run the model update (in parallel if possible)
	for i in 1:dsize[1]
	    for j in 1:dsize[2]
	        advance_model(models[i,j], Ed[i, j, t-1], Ew[i, j, t-1], rain[i, j, t-1], dt, Q)
		advance_model(models_na[i,j], Ed[i, j, t-1], Ew[i, j, t-1], rain[i, j, t-1], dt, Q)
	    end
	end

        # store the model state in an array (and store in output frame)
        fm10_model_state = [ models[i,j].m_ext[2] for i=1:dsize[1], j=1:dsize[2] ]
        fm10_model_na_state = [ models_na[i,j].m_ext[2] for i=1:dsize[1], j=1:dsize[2] ]
        fm10_model_var = [ models[i,j].P[2,2] for i=1:dsize[1], j=1:dsize[2] ]

        spush("fm10_model_state", fm10_model_state)
        spush("fm10_model_na_state", fm10_model_na_state)
        spush("fm10_model_var", fm10_model_var)

        # if observation data for this timepoint is available
        obs_i = Observation[]
        tm_valid_now = filter(x -> abs((mt - x).millis) / 1000.0 <= assim_time_win/2, obs_times)

        # gather all observations
        for tvn in tm_valid_now append!(obs_i, obs_fm10[tvn]) end

        # exclude zero observations - must be sensor failure
        obs_i = filter(x -> obs_value(x) > 0, obs_i)

        # if there are no valid observations, continue with next time step, else run kriging
        if length(obs_i) > 0

            # set the current fm10 model state as the covariate
            X[:,:,1] = fm10_model_state

#            fm10_norm = sum(fm10_model_state.^2)^0.5
            println("INFO: assimilating $(length(obs_i)) observations.")

            # loop over additional covariates
            for i in 2:Xd3
                cov_id = cov_ids[i-1]
                if haskey(st_covar_map, cov_id)
                    # just copy and rescale corresponding static covariate
                    X[:,:,i] = Xr[:,:,i]
                elseif haskey(dyn_covar_map, cov_id)
                    # retrieve the field pointed to by the dynamic covariate id
                    F = dyn_covar_map[cov_id]
                    X[:,:,i] = squeeze(F[:,:,t], 3)
                else
                    error("FATAL: found unknown covariate.")
                end
            end

            # store diagnostic information
            ngp_list = map(x -> nearest_grid_point(x), obs_i)
            stat_ids = map(x -> obs_station_id(x), obs_i)
            m_at_obs = Float64[X[i, j, 1] for (i,j) in  ngp_list]
            m_na_at_obs = Float64[models_na[i,j].m_ext[2] for (i,j) in ngp_list]
            raws = Float64[obs_value(o) for o in obs_i]

            spush("model_raws_mae", mean(abs(m_at_obs - raws)))
            spush("model_na_raws_mae", mean(abs(m_na_at_obs - raws)))

            spush("kriging_obs", raws)
            spush("kriging_obs_station_ids", stat_ids)
            spush("kriging_obs_ngp", ngp_list)

            # compute the kriging estimates and fill in pre-allocated arrays
            trend_surface_model_kriging(obs_i, X, K, V)

            # push diagnostic outputs
            spush("kriging_field", K)
            spush("kriging_variance", V)

            # execute the Kalman update at each grid point
            Kp = zeros(1)
            Vp = zeros(1,1)
            Kg = zeros(Float64, dsize)
            fuel_types = [2]
            for i in 1:dsize[1]
                for j in 1:dsize[2]
                    Kp[1] = K[i,j]
                    Vp[1,1] = V[i,j]
                    Kg[i,j] = kalman_update(models[i,j], Kp, Vp, fuel_types)[1,1]
                end
            end

            # push the fm10 model state after the assimilation
            fm10_model_state = [ models[i,j].m_ext[2] for i=1:dsize[1], j=1:dsize[2] ]
            spush("fm10_model_state_assim", fm10_model_state)

            # retrieve adjustments to time constants and to equilibria
            fm10_adj = zeros(Float64, (6, dsize[1], dsize[2]))
            for i in 1:dsize[1]
                for j in 1:dsize[2]
                    fm10_adj[:,i,j] = models[i,j].m_ext[Nf+1:2*Nf+3]
                end
            end
            fm10_adj_max = [ max(abs(fm10_adj[i,:,:])) for i in 1:6 ]
            spush("fm10_model_deltas", fm10_adj_max)

            # gather model values at ngp points after assimilation
            m_at_obs = Float64[fm10_model_state[i, j] for (i,j) in  ngp_list]
            spush("model_raws_mae_assim", mean(abs(m_at_obs - raws)))

            spush("kalman_gain_fm10", Kg)

            # move to the next storage frame
            next_frame()

        end # if there is anything to assimilate

        # store current assimilated state in the wrfout file
        fm_stor = zeros(Float64, (dsize[1], dsize[2], 1))
        fm_stor[:,:,1] = [ models[i,j].m_ext[1] for i=1:dsize[1], j=1:dsize[2] ]
        NetCDF.putvar(nc, "FM1", [1, 1, t], fm_stor)
        fm_stor[:,:,1] = [ models[i,j].m_ext[2] for i=1:dsize[1], j=1:dsize[2] ]
        NetCDF.putvar(nc, "FM10", [1, 1, t], fm_stor)
        fm_stor[:,:,1] = [ models[i,j].m_ext[3] for i=1:dsize[1], j=1:dsize[2] ]
        NetCDF.putvar(nc, "FM100", [1, 1, t], fm_stor)

    end # for each time point

    # Close down the storage system
    Storage.sclose()

    NetCDF.close(nc)

    t2 = Calendar.now()
    println("INFO: simulation completed at $t2 after $(t2-t1).")

end

main(ARGS)


text:
# /Users/kartikeygupta/Workspace/Julia/Gym/NES/deps/roms/smb.nes
import Base: parse
MAGIC = (0x4e, 0x45, 0x53, 0x1a)

struct ROM
    raw_data::Array{UInt8, 1}
    header::Array{UInt8, 1}
end

function ROM(rom_path::AbstractString)
    io = open(rom_path, "r")
    buffer_size = 10
    buffer = zeros(UInt8, buffer_size)
    raw_data = UInt8[]
    while !eof(io)
        readbytes!(io, buffer)
        push!(raw_data, buffer...)
    end
    header = raw_data[1:16]
    !all(MAGIC .== header[1:4]) && error("ROM missing magic number in the header.")
    sum(header[12:end]) != 0 && error("ROM Header zero fill bytes are not zero.")

    ROM(raw_data, header)
end

parse(chr::Char) = Base.parse(UInt8, chr)
parse(str::AbstractString) = Base.parse(UInt8, str)

function Base.getproperty(rom::ROM, sym::Symbol)

    """Return the header of the ROM file as bytes."""
    if sym == :header
        return rom.raw_data[1:16]

    """Return the size of the PRG ROM in KB"""
    elseif sym == :prg_rom_size
        return 16 * rom.header[5]

    """Return the size of the CHR ROM in KB"""
    elseif sym == :chr_rom_size
        return 8 * rom.header[6]

    """Return the flags at the 6th byte of the header."""
    elseif sym == :flags_6
        return string(rom.header[7], base=2, pad=8)

    """Return the flags at the 7th byte of the header."""
    elseif sym == :flags_7
        return string(rom.header[8], base=2, pad=8)

    """Return the size of the PRG RAM in KB."""
    elseif sym == :prg_ram_size
        return 8 * (rom.header[9] == 0 ? 1 : rom.header[9])

    """Return the flags at the 9th byte of the header."""
    elseif sym == :flags_9
        return string(rom.header[10], base=2, pad=8)

    """
    Return the flags at the 10th byte of the header.

    Notes:
        - these flags are not part of official specification.
        - ignored in this emulator

    """
    elseif sym == :flags_10
        return string(rom.header[11], base=2, pad=8)

#-------------------------HEADER FLAGS-----------------------------------------#

    """Return the mapper number this ROM uses."""
    elseif sym == :mapper
        return parse(rom.flags_7[1:4] * rom.flags_6[1:4])

    """Return a boolean determining if the ROM ignores mirroring."""
    elseif sym == :is_ignore_mirroring
        return rom.flags_6[5] |> parse |> Bool

    """Return a boolean determining if the ROM has a trainer block."""
    elseif sym == :has_trainer
        return rom.flags_6[6] |> parse |> Bool

    """Return a boolean determining if the ROM has a battery-backed RAM."""
    elseif sym == :has_battery_backed_ram
        return rom.flags_6[7] |> parse |> Bool

    """Return the mirroring mode this ROM uses."""
    elseif sym == :is_vertical_mirroring
        return rom.flags_6[8] |> parse |> Bool

    """
    Return whether this cartridge uses PlayChoice-10.

    Note:
        - Play-Choice 10 uses different color palettes for a different PPU
        - ignored in this emulator

    """
    elseif sym == :has_play_choice_10
        return rom.flags_7[7] |> parse |> Bool

    """
    Return whether this cartridge has VS Uni-system.

    Note:
        VS Uni-system is for ROMs that have a coin slot (Arcades).
        - ignored in this emulator

    """
    elseif sym == :has_vs_unisystem
        return rom.flags_7[8] |> parse |> Bool

    """Returns the TV System this ROM supports"""
    elseif sym == :is_pal
        return rom.flags_9[8] |> parse |> Bool

#----------------------------------ROM------------------------------------------

    """The inclusive starting index of the trainer ROM."""
    elseif sym == :trainer_rom_start
        return 17

    """The inclusive stopping index of the trainer ROM."""
    elseif sym == :trainer_rom_stop
        return rom.has_trainer ? rom.trainer_rom_start + 512 : rom.trainer_rom_start - 1

    """Return the trainer ROM of the ROM file"""
    elseif sym == :trainer_rom
        return rom.raw_data[rom.trainer_rom_start:rom.trainer_rom_stop]

    """The inclusive starting index of the PRG ROM."""
    elseif sym == :prg_rom_start
        return rom.trainer_rom_stop + 1

    """The inclusive stopping index of the PRG ROM."""
    elseif sym == :prg_rom_stop
        return rom.prg_rom_start + rom.prg_rom_size * 2^10 - 1

    """Return the PRG ROM of the ROM file."""
    elseif sym == :prg_rom
        try
            return rom.raw_data[rom.prg_rom_start : rom.prg_rom_stop]
        catch BoundsError
            error("Failed to read PRG-ROM on ROM")
        end

    """The inclusive starting index of the CHR ROM."""
    elseif sym == :chr_rom_start
        return rom.prg_rom_stop + 1

    """The inclusive stopping index of the CHR ROM."""
    elseif sym == :chr_rom_stop
        return rom.chr_rom_start + rom.chr_rom_size * 2^10 - 1

    """Return the CHR ROM of the ROM file."""
    elseif sym == :chr_rom
        try
            return rom.raw_data[rom.chr_rom_start:rom.chr_rom_stop]
        catch BoundsError
            error("Failed to read CHR-ROM on ROM")
        end
    else
        return getfield(rom, sym)
    end
end


text:
export parse_value,
get_pdfcontentops

import Base: peek

_error(msg, io) = error(msg*" at $(position(io)) found $(_peekb(io))")
_error(msg, io, cs...) = error(msg*" at $(position(io)) next: $(_peekb(io)) curr: $cs")

#This function is for testing only
function parse_data(filename)
    ps=util_open(filename,"r")
    try
        while(!eof(ps))
            println(parse_value(ps, x->(length(x), String(x))))
            chomp_space!(ps)
        end
    finally
        util_close(ps)
    end
end

#=
Given a `IOStream`, after possibly any amount of whitespace, return the next
parseable value.
=#
function parse_value(ps::IO, fparse_more=x->(length(x), nothing))
    chomp_space!(ps)
    byte = UInt8(peek(ps))
    byte == LEFT_PAREN ? parse_string(ps) :
    byte == LESS_THAN  ? parse_xstring(ps) :
    byte == PERCENT    ? parse_comment(ps) :
    byte == SOLIDUS    ? parse_name(ps) :
    byte == MINUS_SIGN ||
        byte == PLUS_SIGN || byte == PERIOD ? parse_number(ps) :
    ispdfdigit(byte)   ? try_parse_indirect_reference(ps) :
    byte == LEFT_SB    ? parse_array(ps) :
    parse_pdfOpsOrConst(ps, fparse_more)
end

function parse_comment(ps::IO)
    b = Vector{UInt8}()
    skip(ps,1)  # skip comment marker
    while true
        c = advance!(ps)
        if is_crorlf(c)
            break
        end
        push!(b, c)
    end
    chomp_space!(ps)
    return CosComment(b)
end

function parse_name(ps::IO)
    b = UInt8[]
    skipv(ps,SOLIDUS)  #skip solidus and ensure it
    while true
        c = ps |> _peekb
        if ispdfspace(c) || ispdfdelimiter(c)
            break
        elseif (c == NUMBER_SIGN)
            skip(ps,1)
            #Now look for 2 hex numbers
            c1 = ps |> _peekb
            skip(ps,1)
            c2 = ps |> _peekb
            if ispdfxdigit(c1) && ispdfxdigit(c2)
                c = UInt8(gethexval(c1)*16+gethexval(c2))
            else
                _error(E_UNEXPECTED_CHAR, ps, c1, c2)
            end
        end
        skip(ps,1)
        push!(b, c)
    end
    chomp_space!(ps)
    return CosName(String(b))
end

function parse_pdfOpsOrConst(ps::IO, fparse_more::Function)
    b = UInt8[]
    mark(ps)
    try
        while !eof(ps)
            c = ps |> _peekb
            (ispdfspace(c) || ispdfdelimiter(c)) && break
            skip(ps, 1)
            push!(b, c)
        end
        ns = chomp_space!(ps)
        # It's ok to skip all spaces till eof
        # but if eof is not reached read only one space
        # If no space is found then do not advance
        reset_marker = ns > 0 && !eof(ps)
        nused, ret = length(b)+1, get_pdfconstant(b)
        if ret === nothing
            nused, ret = fparse_more(b)
            if nused < length(b)
                reset_marker = true
            else
                nused += 1
            end
        end
        if reset_marker
            reset(ps)
            skip(ps, nused)
        end
        return ret
    finally
        unmark(ps)
    end
end

function get_pdfconstant(b::Vector{UInt8})
    b == [LATIN_T,LATIN_R,LATIN_U,LATIN_E]              && return CosTrue
    b == [LATIN_F,LATIN_A, LATIN_L, LATIN_S, LATIN_E]   && return CosFalse
    b == [LATIN_N, LATIN_U, LATIN_L, LATIN_L]           && return CosNull
    return nothing
end

function parse_array(ps::IO)
    result = CosArray()
    @inbounds skip(ps,1)  # Skip over opening '['
    chomp_space!(ps)
    if ps |> _peekb != RIGHT_SB  # special case for empty array
        while true
            push!(result.val, parse_value(ps))
            chomp_space!(ps)
            b = ps |> _peekb
            b == RIGHT_SB && break
        end
    end

    @inbounds skip(ps,1)
    chomp_space!(ps)
    result
end

function read_octal_escape!(c, ps)
    local n::UInt8 = getnumval(c)
    for _ in 1:2
        b = ps |> _peekb
        !ispdfodigit(b) && return n
        n = (n << 3) + getnumval(b)
        skip(ps,1)
    end
    return n
end

function parse_string(ps::IO)
  b = UInt8[]
  skip(ps,1)  # skip opening quote
  local paren_cnt = 0
  while true
    c = advance!(ps)

    if c == BACKSLASH
      c = advance!(ps)
      if ispdfodigit(c) #Read octal digits
        append!(b, read_octal_escape!(c,ps))
      elseif is_crorlf(c) #ignore the solidus, EOLs and move on
        chomp_space!(ps)
      else
        c = get(ESCAPES, c, 0x00)
        c == 0x00 && error(E_BAD_ESCAPE)
        push!(b, c)
      end
      continue
    elseif c == LEFT_PAREN
      paren_cnt+=1
    elseif c == RIGHT_PAREN
      if (paren_cnt > 0)
        paren_cnt-=1
      else
        chomp_space!(ps)
        return CosLiteralString(b)
      end
    end
    push!(b, c)
  end
end

function parse_xstring(ps::IO)
    b = UInt8[]
    skip(ps, 1)  # skip open LT

    while true
        c = advance!(ps)
        c == LESS_THAN && return parse_dict(ps)
        # As per section 7.3.4.3 space must be ignored
        ispdfspace(c) && continue
        if ispdfxdigit(c)
            push!(b, c)
        elseif c == GREATER_THAN
            length(b) % 2 != 0 && push!(b, NULL)
            chomp_space!(ps)
            return CosXString(b)
        else
            _error(E_UNEXPECTED_CHAR, ps, c)
        end
    end
end

function parse_dict(ps::IO)
    #Move the cursor beyond < char
    chomp_space!(ps)

    dict=CosDict()

    while(true)
        # Empty dict File 431 stillhq
        if ps |> _peekb == SOLIDUS
            key = parse_name(ps)
            chomp_space!(ps)

            val = parse_value(ps)
            if (val !== CosNull)
                dict.val[key] = val
            end
        end

        chomp_space!(ps)

        c = ps |> _peekb
        (c == SOLIDUS) && continue
        skip(ps, 1)
        if c == GREATER_THAN
            skipv(ps, GREATER_THAN)
            break
        end
    end
    chomp_space!(ps)
    return dict
end

function ensure_line_feed_eol(ps::IO)
  c = advance!(ps)
  if (c == RETURN)
      skipv(ps,LINE_FEED)
  elseif (c == LINE_FEED)
      return c
  else
      _error(E_UNEXPECTED_CHAR, ps, c)
  end
end

#=
Read the internal stream data and externalize to a temp file.
If it's already an externalized stream then false is returned.
The value can be stored in the stream object attribute so that the reverse
process will be carried out for serialization.
=#
function read_internal_stream_data(ps::IO, extent::CosDict, len::Int)
    get(extent, cn"F") != CosNull && return false

    (path, io) = get_tempfilepath()
    try
        data = read(ps, len)
        write(io, data)
    finally
        util_close(io)
    end

    #Ensuring all the data is written to a file
    set!(extent, cn"F", CosLiteralString(path))

    filter = get(extent, cn"Filter")
    if (filter != CosNull)
        set!(extent, cn"FFilter", filter)
        set!(extent, cn"Filter", CosNull)
    end

    parms = get(extent, cn"DecodeParms")
    if (parms != CosNull)
        set!(extent, cn"FDecodeParms",  parms)
        set!(extent, cn"DecodeParms", CosNull)
    end

    return true
end

mutable struct CosObjectLoc
    loc::Int
    stm::IDDNRef{CosObjectStream}
    obj::CosObject
    CosObjectLoc(l, s=CosNull, o=CosNull) = new(l, s, o)
end

process_stream_length(stmlen::CosInt,
                      ps::IO,
                      hoffset::Int,
                      xref::Dict{CosIndirectObjectRef, CosObjectLoc})=stmlen

function process_stream_length(stmlen::CosIndirectObjectRef,
                               ps::IO,
                               hoffset::Int,
                               xref::Dict{CosIndirectObjectRef, CosObjectLoc})
    cosObjectLoc = xref[stmlen]
    if (cosObjectLoc.obj === CosNull)
        seek(ps, cosObjectLoc.loc + hoffset)
        cosObjectLoc.obj = parse_indirect_obj(ps, hoffset, xref)
    end
    return cosObjectLoc.obj
end

function postprocess_indirect_object(ps::IO, hoffset::Int, obj::CosDict,
                                     xref::Dict{CosIndirectObjectRef,
                                                CosObjectLoc})
    if locate_keyword!(ps, STREAM) == 0
        ensure_line_feed_eol(ps)
        pos = position(ps)

        stmlen = get(obj, cn"Length")

        lenobj = process_stream_length(stmlen, ps, hoffset, xref)

        len = get(lenobj)

        if (lenobj != stmlen)
            set!(obj, cn"Length", lenobj)
        end

        seek(ps, pos)

        # Here you can make sure file data is decoded into a file
        # later it can be made into a memory based on size etc.
        # Since, these are temporary files the spec is system file only
        isInternal = read_internal_stream_data(ps, obj, len)

        obj = CosStream(obj, isInternal)

        #Now eat away the ENDSTREAM token
        chomp_space!(ps)
        skipv(ps,ENDSTREAM)
        obj = createIfObjectStream(obj)
    end
    return obj
end

postprocess_indirect_object(ps::IO,
                            hoffset::Int,
                            obj::CosObject,
                            xref::Dict{CosIndirectObjectRef, CosObjectLoc}) = obj

function parse_indirect_obj(ps::IO,
                            hoffset::Int,
                            xref::Dict{CosIndirectObjectRef, CosObjectLoc})
    chomp_space!(ps)
    objn = parse_unsignednumber(ps).val
    chomp_space!(ps)
    genn = parse_unsignednumber(ps).val
    chomp_space!(ps)
    skipv(ps, OBJ)
    obj = parse_value(ps)
    chomp_space!(ps)
    obj = postprocess_indirect_object(ps, hoffset, obj, xref)
    chomp_space!(ps)
    skipv(ps,ENDOBJ)
    return CosIndirectObject(objn, genn, obj)
end

function parse_indirect_ref(ps::IO)
    chomp_space!(ps)
    objn = parse_unsignednumber(ps).val
    chomp_space!(ps)
    genn = parse_unsignednumber(ps).val
    chomp_space!(ps)
    skipv(ps, LATIN_UPPER_R)
    # There has to be a PDF delimiter or space after this.
    # This may fail the content parser due to the RG operator
    b = ps |> _peekb
    ispdfdelimiter(b) || ispdfspace(b) ||
        error(E_INVALID_DELIMITER)
    chomp_space!(ps)
    return CosIndirectObjectRef(objn, genn)
end

function try_parse_indirect_reference(ps::IO)
    nobj = parse_number(ps)
    nobj isa CosFloat && return nobj
    chomp_space!(ps)
    mark(ps)
    if ispdfdigit(ps |> _peekb)
        objn = nobj.val
        genn = parse_unsignednumber(ps).val
        chomp_space!(ps)
        if (ps |> _peekb == LATIN_UPPER_R)
            skip(ps, 1)
            # There has to be a PDF delimiter or space after this.
            # This may fail the content parser due to the RG operator
            b = ps |> _peekb
            if ispdfdelimiter(b) || ispdfspace(b)
                unmark(ps)
                chomp_space!(ps)
                return CosIndirectObjectRef(objn, genn)
            else
                reset(ps)
                chomp_space!(ps)
                return nobj
            end
        else
            reset(ps)
            chomp_space!(ps)
            return nobj
        end
    else
        unmark(ps)
        chomp_space!(ps)
        return nobj
    end
end

#=
Return `true` if the given bytes vector, starting at `from` and ending at `to`,
has a leading zero.
=#
function hasleadingzero(bytes::Vector{UInt8}, from::Int, to::Int)
    c = bytes[from]
    from + 1 < to && c == UInt8('-') &&
            bytes[from + 1] == DIGIT_ZERO && ispdfdigit(bytes[from + 2]) ||
    from < to && to > from + 1 && c == DIGIT_ZERO &&
            ispdfdigit(bytes[from + 1])
end

#=
Parse a float from the given bytes vector, starting at `from` and ending at the
byte before `to`. Bytes enclosed should all be ASCII characters.
=#
float_from_bytes(bytes::Vector{UInt8}) = tryparse(Float64, String(bytes))

#=
Parse an integer from the given bytes vector, starting at `from` and ending at
the byte before `to`. Bytes enclosed should all be ASCII characters.
=#
function int_from_bytes(bytes::Vector{UInt8})
    from = 1
    to = length(bytes)
    @inbounds isnegative = bytes[from] == MINUS_SIGN ? (from += 1; true) : false
    num = Int(0)
    @inbounds for i in from:to
        num = Int(10) * num + Int(bytes[i] - DIGIT_ZERO)
    end
    return ifelse(isnegative, -num, num)
end

function number_from_bytes(ps::IO, isint::Bool, bytes::Vector{UInt8})
    from = 1
    to = length(bytes)
    if isint
        @inbounds if to == from && bytes[from] == MINUS_SIGN
            _error(E_BAD_NUMBER, ps)
        end
        num = int_from_bytes(bytes)
        return CosInt(num)
    else
        res = float_from_bytes(bytes)
        res === nothing && _error(E_BAD_NUMBER, ps)
        return CosFloat(res)
    end
end

function parse_unsignednumber(ps::IO)
    number = UInt8[]
    isint = true
    while true
        c = ps |> _peekb
        if ispdfdigit(c)
            push!(number, UInt8(c))
        else
            break
        end
        skip(ps,1)
    end
    chomp_space!(ps)
    return number_from_bytes(ps, isint, number)
end

function parse_number(ps::IO)
    number = UInt8[]
    isint = true

    while true
        c = ps |> _peekb

        if ispdfdigit(c) || c == MINUS_SIGN
            push!(number, UInt8(c))
        elseif c == PLUS_SIGN
        elseif c==DECIMAL_POINT
            push!(number, UInt8(c))
            isint = false
        else
            break
        end

        skip(ps, 1)
    end
    chomp_space!(ps)
    return number_from_bytes(ps, isint, number)
end


text:
using FileIO
using Images
using Logging

export
    run_tesseract

# Tesseract settings
command = "tesseract"
psm_valid_range = 0:14
oem_valid_range = 0:4

"""Util to check and inform whether Tesseract is installed or not"""
function check_tesseract_installed()
    try
        read(`$command --version`, String);
    catch IOError
        @error "Tesseract is not properly installed. Command $command not recognized"
    end
end

check_tesseract_installed()

"""Util to get version of Tesseract installed"""
function get_tesseract_version()
    info = read(`$command --version`, String)
    version = split(info, "\n")[1]
    return version
end

@info "Using Tesseract version: $(get_tesseract_version())"

"""
    run_tesseract(args...; kwargs...) -> Bool

Wrapper function to run Tesseract over a image stored in disk,
    and write the results in a given path.
Errors / Warnings are reported through `Logging`, so no exceptions are thrown.

# Arguments
- `input_path::String`: Path to the image to be processed
- `output_path::String`: Path to the text result to be written

# Keywords
- `lang::Union{String, Nothing}` Language to be configured in Tesseract (optional).
- `psm::Integer`: Page segmentation modes (PSM):
    - `psm=0`:   Orientation and script detection (OSD) only.
    - `psm=1`:   Automatic page segmentation with OSD.
    - `psm=2`:   Automatic page segmentation, but no OSD, or OCR.
    - `psm=3`:   Fully automatic page segmentation, but no OSD. (Default)
    - `psm=4`:   Assume a single column of text of variable sizes.
    - `psm=5`:   Assume a single uniform block of vertically aligned text.
    - `psm=6`:   Assume a single uniform block of text.
    - `psm=7`:   Treat the image as a single text line.
    - `psm=8`:   Treat the image as a single word.
    - `psm=9`:   Treat the image as a single word in a circle.
    - `psm=10`:  Treat the image as a single character.
    - `psm=11`:  Sparse text. Find as much text as possible in no particular order.
    - `psm=12`:  Sparse text with OSD.
    - `psm=13`:  Raw line. Treat the image as a single text line,
        bypassing hacks that are Tesseract-specific.
- `oem::Integer`: OCR Engine modes (OEM):
    - `oem=0`:   Legacy engine only.
    - `oem=1`:   Neural nets LSTM engine only. (Default)
    - `oem=2`:   Legacy + LSTM engines.
    - `oem=3`:   Default, based on what is available.
- `kwargs`: Other key-value pairs to be sent to Tesseract command as "-c" config variables.
    You can check the options with `tesseract --print-parameters`.

# Returns
- `Bool`: indicating whether execution was successful or not

# Examples
```julia-repl
julia> using OCReract;
julia> img_path = "/path/to/img.png";
julia> out_path = "/tmp/tesseract_result.txt";
julia> run_tesseract(img_path, out_path, psm=3, oem=1)
```

"""
function run_tesseract(
    input_path::String,
    output_path::String;
    lang::Union{String, Nothing}=nothing,
    psm::Integer=3,
    oem::Integer=1,
    kwargs...
)
    # Check arguments
    if !isfile(input_path)
        @error "Input path '$input_path' doesn't exist!"
        return false
    end

    if isfile(output_path)
        @warn "Output path '$output_path' already exists!"
    end

    if !(psm in psm_valid_range)
        psm = 3
        @warn "PSM parameter not in valid range: [$(string(psm_valid_range))]. "*
              "Changing to default value: psm=$psm"
    end

    if !(oem in oem_valid_range)
        oem = 1
        @warn "OEM parameter not in valid range: [$(string(oem_valid_range))]. "*
              "Changing to default value: oem=$oem"
    end

    # Build command to execute
    cmd = "tesseract $input_path $output_path"

    # Add language if specified
    if lang !== nothing
        cmd = join([cmd, "-l $lang"], " ")
    end

    # Add config variables from kwargs
    config_vars = ""
    for (k, v) in kwargs
        config_vars *= "-c $k=$v "
    end

    # Join arguments to final command
    cmd = join([cmd, "--oem $oem", "--psm $psm", config_vars], " ")
    @debug "Running command '$cmd' ..."

    # Run command
    try
        run(`$(split(cmd))`)
    catch e
        @error "Error ocurred while running Tesseract! $e"
        return false
    end

    # NOTE: tesseract ALWAYS add a suffix ".txt" to the result
    # so, rename the resulting file to call it as the user defined
    output_file = output_path*".txt"
    try
        mv(output_file, output_path, force=true)
    catch e
        # At least, inform that to the user
        @warn "Result was stored in '$output_file'"
    end

    return true
end

"""
    run_tesseract(args...; kwargs...) -> String

Function to run Tesseract over an image in memory, and get the results in a `String`.
Errors / Warnings are reported through `Logging`, so no exceptions are thrown.

# Arguments
- `image`: Image to be processed, in a format compatible with `Images` module.

# Keywords
- `lang::Union{String, Nothing}` Language to be configured in Tesseract (optional)
- `psm::Integer`: Page segmentation modes (PSM):
    - `psm=0`:   Orientation and script detection (OSD) only.
    - `psm=1`:   Automatic page segmentation with OSD.
    - `psm=2`:   Automatic page segmentation, but no OSD, or OCR.
    - `psm=3`:   Fully automatic page segmentation, but no OSD. (Default)
    - `psm=4`:   Assume a single column of text of variable sizes.
    - `psm=5`:   Assume a single uniform block of vertically aligned text.
    - `psm=6`:   Assume a single uniform block of text.
    - `psm=7`:   Treat the image as a single text line.
    - `psm=8`:   Treat the image as a single word.
    - `psm=9`:   Treat the image as a single word in a circle.
    - `psm=10`:  Treat the image as a single character.
    - `psm=11`:  Sparse text. Find as much text as possible in no particular order.
    - `psm=12`:  Sparse text with OSD.
    - `psm=13`:  Raw line. Treat the image as a single text line,
        bypassing hacks that are Tesseract-specific.
- `oem::Integer`: OCR Engine modes (OEM):
    - `oem=0`:   Legacy engine only.
    - `oem=1`:   Neural nets LSTM engine only. (Default)
    - `oem=2`:   Legacy + LSTM engines.
    - `oem=3`:   Default, based on what is available.
- `kwargs`: Other key-value pairs to be sent to Tesseract command as "-c" config variables.
    You can check the options with `tesseract --print-parameters`.

# Returns
- `String`: text extracted, or empty string in case an error occurs

# Examples
```julia-repl
julia> using Images;
julia> using OCReract;
julia> img_path = "/path/to/img.png";
julia> img = Images.load(img_path);
julia> res_text = run_tesseract(img, psm=3, oem=1);
julia> println(strip(res_text));
```
"""
function run_tesseract(
    image;
    lang::Union{String, Nothing}=nothing,
    psm::Integer=3,
    oem::Integer=1,
    kwargs...
)
    input_path = tempname() * ".png"   # PNG image resulting
    output_path = tempname() * ".txt"  # TXT file resulting

    # Save image into disk to be handled by Tesseract
    FileIO.save(input_path, image)

    # Run Tesseract!
    run_tesseract(input_path, output_path, lang=lang, psm=psm, oem=oem; kwargs...)

    # Read output into a string
    txt = ""
    open(output_path, "r") do f
        txt = read(f, String)
    end

    # Now remove these temporary files generated
    try
        @debug "Removing temporary files generated..."
        rm(input_path)
        rm(output_path)
    catch e
        @warn "Error ocurred while removing temporary files! $e"
    end

    return txt
end


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

name = "NEO"
version = v"21.44.21506"

# Collection of sources required to build this package
sources = [
    GitSource("https://github.com/intel/compute-runtime.git",
              "586eaeaf2ebceeb9c0095ecd2f774a272038afbc"),
]

# Bash recipe for building across all platforms
script = raw"""
cd compute-runtime
install_license LICENSE.md

# work around compilation failures
## already defined in gmmlib
sed -i '/__stdcall/d' shared/source/gmm_helper/gmm_lib.h
## extend LD_LIBRARY_PATH, don't overwrite it
find . \( -name CMakeLists.txt -or -name '*.cmake' \) -exec \
    sed -i 's/LD_LIBRARY_PATH=/LD_LIBRARY_PATH=$ENV{LD_LIBRARY_PATH}:/g' '{}' \;

CMAKE_FLAGS=()

# Release build for best performance
CMAKE_FLAGS+=(-DCMAKE_BUILD_TYPE=Release)

# Install things into $prefix
CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})

# NOTE: NEO currently can't cross compile because of its IGC dependency
CMAKE_FLAGS+=(-DCMAKE_CROSSCOMPILING:BOOL=OFF)

# Explicitly use our cmake toolchain file
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})

# Don't run tests
CMAKE_FLAGS+=(-DSKIP_UNIT_TESTS:Bool=true)

# we don't care about cl_intel_va_api_media_sharing
CMAKE_FLAGS+=(-DDISABLE_LIBVA:Bool=true)

# enable support for the DG1
CMAKE_FLAGS+=(-DSUPPORT_DG1:Bool=true)

# libigc installs libraries and pkgconfig rules in lib64, so look for them there.
# FIXME: shouldn't BinaryBuilder do this?
export PKG_CONFIG_PATH=${prefix}/lib64/pkgconfig:${prefix}/lib/pkgconfig

cmake -B build -S . -GNinja ${CMAKE_FLAGS[@]}
ninja -C build -j ${nproc} install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    # NEO is 64-bit only: https://github.com/intel/compute-runtime/issues/179
    # and does not support musl: https://github.com/intel/compute-runtime/issues/265
    Platform("x86_64", "linux", libc="glibc"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("ocloc", :ocloc),
    LibraryProduct("libigdrcl", :libigdrcl, ["lib/intel-opencl", "lib64/intel-opencl"]),
    LibraryProduct("libze_intel_gpu", :libze_intel_gpu),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("gmmlib_jll"; compat="=21.2.1"),
    Dependency("libigc_jll"; compat="=1.0.8744"),
    Dependency("oneAPI_Level_Zero_Headers_jll", v"1.2.43"; compat="~1.2"),  # XXX: don't specify patch version
]

# GCC 4 has constexpr incompatibilities
# GCC 7 triggers: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=79929
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"8")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using CUDA
using CUDA.CUSPARSE
using ExaPF
using KernelAbstractions
using LinearAlgebra
using Random
using SparseArrays
using Test
using TimerOutputs

import ExaPF: Parse, PowerSystem

Random.seed!(2713)

case = "case14.raw"
# case = "ACTIVSg70K.raw"

@testset "Powerflow residuals and Jacobian" begin
    # read data
    to = TimerOutputs.TimerOutput()
    datafile = joinpath(dirname(@__FILE__), case)
    data = Parse.parse_raw(datafile)
    BUS_B, BUS_AREA, BUS_VM, BUS_VA, BUS_NVHI, BUS_NVLO, BUS_EVHI,
    BUS_EVLO, BUS_TYPE = Parse.idx_bus()
    bus = data["BUS"]
    nbus = size(bus, 1)

    # obtain V0 from raw data
    V = Array{Complex{Float64}}(undef, nbus)
    T = Vector
    for i in 1:nbus
        V[i] = bus[i, BUS_VM]*exp(1im * pi/180 * bus[i, BUS_VA])
    end

    # form Y matrix
    Ybus, Yf_br, Yt_br, Yf_tr, Yt_tr = PowerSystem.makeYbus(data);

    Vm = abs.(V)
    Va = angle.(V)
    bus = data["BUS"]
    gen = data["GENERATOR"]
    load = data["LOAD"]
    nbus = size(bus, 1)
    ngen = size(gen, 1)
    nload = size(load, 1)

    ybus_re, ybus_im = ExaPF.Spmat{T}(Ybus)
    SBASE = data["CASE IDENTIFICATION"][1]
    Sbus = PowerSystem.assembleSbus(gen, load, SBASE, nbus)
    pbus = real(Sbus)
    qbus = imag(Sbus)

    ref, pv, pq = PowerSystem.bustypeindex(bus, gen)
    npv = size(pv, 1);
    npq = size(pq, 1);

    @testset "Computing residuals" begin
        F = zeros(Float64, npv + 2*npq)
        # First compute a reference value for resisual computed at V
        F♯ = ExaPF.residualFunction(V, Ybus, Sbus, pv, pq)
        # residual_polar! uses only binary types as this function is meant
        # to be deported on the GPU
        ExaPF.residualFunction_polar!(F, Vm, Va,
            ybus_re.nzval, ybus_re.colptr, ybus_re.rowval,
            ybus_im.nzval, ybus_im.colptr, ybus_im.rowval,
            pbus, qbus, pv, pq, nbus)
        @test F ≈ F♯
    end
    @testset "Computing Jacobian of residuals" begin
        F = zeros(Float64, npv + 2*npq)
        # Compute Jacobian at point V manually and use it as reference
        J = ExaPF.residualJacobian(V, Ybus, pv, pq)
        J♯ = copy(J)

        # Then, create a JacobianAD object
        coloring = ExaPF.matrix_colors(J)
        jacobianAD = ExaPF.AD.JacobianAD(J, coloring, F, Vm, Va, pv, pq)
        # and compute Jacobian with ForwardDiff
        ExaPF.AD.residualJacobianAD!(
            jacobianAD, ExaPF.residualFunction_polar!, Vm, Va,
            ybus_re, ybus_im, pbus, qbus, pv, pq, nbus, to)
        @test jacobianAD.J ≈ J♯
    end
end

@testset "Wrapping of iterative solvers" begin
    nblocks = 2
    n, m = 32, 32
    to = TimerOutputs.TimerOutput()

    # Add a diagonal term for conditionning
    A = randn(n, m) + 15I
    x♯ = randn(m)
    b = A * x♯

    # Be careful: all algorithms work with sparse matrix
    A = sparse(A)
    # Init iterators for KernelAbstractions
    if has_cuda_gpu()
        iterators = zip([CPU(), CUDADevice()],
                        [Array{Float64}, CuVector{Float64}],
                        [SparseMatrixCSC, CuSparseMatrixCSR])
    else
        iterators = zip([CPU()], [Array{Float64}], [SparseMatrixCSC])
    end
    @testset "Run iterative solvers on device $(device)" for (device, V, SM) in iterators
        As = SM(A)
        bs = convert(V, b)
        x0 = convert(V, zeros(m))
        xs♯ = convert(V, x♯)

        # TODO: currently Preconditioner takes as input only sparse matrix
        # defined on the main memory.
        precond = ExaPF.Precondition.Preconditioner(A, nblocks, device)

        # First test the custom implementation of BICGSTAB
        @testset "BICGSTAB" begin
            # Need to update preconditioner before resolution
            ExaPF.Precondition.update(As, precond, to)
            P = precond.P
            fill!(x0, 0.0)
            x_sol, n_iters = ExaPF.bicgstab(As, bs, P, x0, to)
            @test n_iters <= m
            @test x_sol ≈ xs♯ atol=1e-6
        end
        @testset "Interface for iterative algorithm ($algo)" for algo in ExaPF.list_solvers(device)
            fill!(x0, 0.0)
            n_iters = ExaPF.Iterative.ldiv!(x0, As, bs, algo, precond, to)
            @test n_iters <= m
            @test x0 ≈ xs♯ atol=1e-6
        end
    end
end

@testset "Powerflow solver" begin
    datafile = joinpath(dirname(@__FILE__), case)
    nblocks = 8
    # Create a network object:
    pf = ExaPF.PowerSystem.PowerNetwork(datafile)
    x = ExaPF.PowerSystem.get_x(pf)
    u = ExaPF.PowerSystem.get_u(pf)
    p = ExaPF.PowerSystem.get_p(pf)
    target = CPU()
    @testset "[CPU] Powerflow solver $precond" for precond in ExaPF.list_solvers(target)
        sol, has_conv, res = solve(pf, x, u, p, nblocks, precond, device=target)
        @test has_conv
        @test res < 1e-6
    end

    if has_cuda_gpu()
        target = CUDADevice()
        @testset "[GPU] Powerflow solver $precond" for precond in ExaPF.list_solvers(target)
            sol, conv, res = solve(pf, x, u, p, nblocks, precond, device=target)
            @test conv
            @test res < 1e-6
        end
    end
end

# # Not working yet. Will check whether Ipopt and reduced method match in objective
# @testset "rgm_3bus" begin
#    include("../scripts/rgm_3bus.jl")
#    @show red_cost = cfun(xk, uk, p)
#    include("../scripts/ipopt.jl")
#    @show ipopt_cost = cfun(xk, uk, p)
#    gap = abs(red_cost - ipopt_cost)
#    println("gap = abs(red_cost - ipopt_cost): $gap = abs($red_cost - $ipopt_cost)")
#    @test gap ≈ 0.0
# end

# @testset "rgm_3bus_ref" begin
#    include("../scripts/rgm_3bus_ref.jl")
#    @show red_cost = cfun(xk, uk, p)
#    include("../scripts/ipopt_ref.jl")
#    @show ipopt_cost = cfun(xk, uk, p)
#    gap = abs(red_cost - ipopt_cost)
#    println("gap = abs(red_cost - ipopt_cost): $gap = abs($red_cost - $ipopt_cost)")
#    @test gap ≈ 0.0
# end


text:
using DemoCards: democard, walkpage, DemoPage

@testset "regexes" begin
    root = joinpath("assets", "regexes")

    @testset "title" begin
        config_1_files = ("example_1.md",
                          "example_3.md",
                          "example_4.md",
                          "example_6.md",
                          "example_7.md",
                          "example_8.md",
                          "example_1.jl",
                          "example_4.jl",
                          "example_5.jl",
                          "example_6.jl",
                          "example_8.jl",
                          "example_9.jl",
                          "example_10.jl",
                          "example_11.jl",
                          "example_12.jl",
                          "example_13.jl")
        foreach(config_1_files) do filename
            card = democard(joinpath(root, filename))
            @test card.title == "This is a title"
            @test card.id == "This-is-a-title"
            @test card.description == "This is a description"
        end

        config_2_files =  ("example_2.md", "example_2.jl", "example_3.jl")
        foreach(config_2_files) do filename
            card = democard(joinpath(root, filename))
            @test card.title == "This is a title"
            @test card.id == "custom_id"
            @test card.description == "This is a description"
        end

        config_3_files = ("example_5.md", "example_7.jl")
        foreach(config_3_files) do filename
            card = democard(joinpath(root, filename))
            @test card.title == "This is a title"
            @test card.id == "This-is-a-title"
            @test card.description == "This is a description that spreads along multiple lines"
        end
    end

    # only lines before body contents can be potential titles
    config = DemoCards.parse(Val(:Julia), joinpath(root, "title_before_body.jl"))
    @test haskey(config, "description")
    @test config["description"] == "This is parsed as a description line --- [1]"
    @test !haskey(config, "title")
    @test !haskey(config, "id")
    @test !haskey(config, "cover")

    config = DemoCards.parse(Val(:Markdown), joinpath(root, "title_before_body.md"))
    @test haskey(config, "description")
    @test config["description"] == "This is parsed as a description line --- [1]"
    @test !haskey(config, "title")
    @test !haskey(config, "id")
    @test !haskey(config, "cover")

    # check if upstream changes on "Edit On GitHub" button breaks the redirect_url logic
    try
        example_url = "https://juliadocs.github.io/Documenter.jl/stable/"
        regex_url = r"http[s]?://.*index\.md"

        contents = String(HTTP.get(example_url; readtimeout=3, retry=false).body)

        m = match(DemoCards.regex_edit_on_github, contents)
        @test isa(m, RegexMatch)
        if m isa RegexMatch
            @test isa(match(regex_url, m.captures[1]), RegexMatch)
        end
    catch err
        # if network is down...
        @warn err
    end
end

@testset "url_redirect" begin
    source = "quickstart"
    build_url = "https://github.com/johnnychen94/DemoCards.jl/blob/master/docs/src/quickstart/usage_example/julia_demos/2.cover_on_the_fly.md"
    src_url = "https://github.com/johnnychen94/DemoCards.jl/blob/master/docs/quickstart/usage_example/julia_demos/2.cover_on_the_fly.jl"
    @test src_url == DemoCards.get_source_url(build_url, source, "2.cover_on_the_fly.jl", "src")

    source = joinpath("theme_gallery", "grid")
    build_url = "https://github.com/johnnychen94/DemoCards.jl/blob/master/docs/src/grid/grid_section_1/grid_subsection_1/grid_card_1.md"
    src_url = "https://github.com/johnnychen94/DemoCards.jl/blob/master/docs/theme_gallery/grid/grid_section_1/grid_subsection_1/grid_card_1.md"
    @test src_url == DemoCards.get_source_url(build_url, source, "grid_card_1.md", "src")
end

@testset "walkpage" begin
    page = DemoPage(joinpath("assets", "page", "hidden"))
    reference = "Hidden" => ["hidden1.jl", "hidden2.md", "normal.md"]
    @test reference == walkpage(page) do item
        basename(item.path)
    end
    reference = "Hidden" => ["Sec" => ["hidden1.jl", "hidden2.md", "normal.md"]]
    @test reference == walkpage(page; flatten=false) do item
        basename(item.path)
    end

    page = DemoPage(joinpath("assets", "page", "one_card"))
    reference = "One card" => ["card.md"]
    @test reference == walkpage(page) do item
        basename(item.path)
    end
    reference = "One card" => ["Section" => ["Subsection" => ["card.md"]]]
    @test reference == walkpage(page; flatten=false) do item
        basename(item.path)
    end
end


text:
using MarketData

rng = 1:50:500

apple1  = Stock(:AAPL, cl.values[rng[1]], 100)
apple2  = Stock(:AAPL, cl.values[rng[2]], -100)
apple3  = Stock(:AAPL, cl.values[rng[3]], 100)
apple4  = Stock(:AAPL, cl.values[rng[4]], -100)
apple5  = Stock(:AAPL, cl.values[rng[5]], 100)
apple6  = Stock(:AAPL, cl.values[rng[6]], -100)
apple7  = Stock(:AAPL, cl.values[rng[7]], 100)
apple8  = Stock(:AAPL, cl.values[rng[8]], -100)
apple9  = Stock(:AAPL, cl.values[rng[9]], 100)
apple10 = Stock(:AAPL, cl.values[rng[10]], -100)

port = [apple1,
        apple2,
        apple3,
        apple4,
        apple5,
        apple6,
        apple7,
        apple8,
        apple9,
        apple10]

dts = [cl.timestamp[rng[1]],
       cl.timestamp[rng[2]],
       cl.timestamp[rng[3]],
       cl.timestamp[rng[4]],
       cl.timestamp[rng[5]],
       cl.timestamp[rng[6]],
       cl.timestamp[rng[7]],
       cl.timestamp[rng[8]],
       cl.timestamp[rng[9]],
       cl.timestamp[rng[10]]]

res = Timestamp[]

for d in 1:10
push!(res, Timestamp(dts[d], port[d]))
end

vals = zeros(10);

for r in 1:10
    #vals[r] = -1*(res[r].value.quantity * res[r].value.basis)
    vals[r] = res[r].value.quantity * res[r].value.basis
end