text:
# Continuous Moving Bed Reconstruction
#
# All data and results are in single-precision (32-bit) float, to
# reduce memory requirements.

using NFFT
using MAT
include("correct.jl")
require("Julia Recon Scripts/coords.jl")
require("Julia Recon Scripts/pnfft.jl")

function coilcompress{T}(x::Array{T,3}; nvirt=6)
  # Geometric coil compression with alignment to middle
  println("compressing data to $nvirt virtual coils")
  nchan, nro, npe = size(x)
  y = Array(T, nvirt, nro, npe)
  n = div(npe,2)  # middle profile, the starting point
  U0 = Array(T, nchan, nvirt)
  Uprev = Array(T, nchan, nvirt)
  for p in n
    t = squeeze(x[:,:,p], 3)
    C = t*t'
    U, s, V = svd(C)
    U0[:] = U[:,1:nvirt]
    y[:,:,p] = U0'*t
  end
  Uprev[:] = U0
  for p in n+1:npe
    t = squeeze(x[:,:,p], 3)
    C = t*t'
    U, s, V = svd(C)
    U = U[:,1:nvirt]
    W, s, V = svd(U'*Uprev)
    U = U*W*V'  # rotate U to align with previous slice
    y[:,:,p] = U'*t
    #Uprev[:] = U  # took this out bc seems unnecessary
  end
  Uprev[:] = U0
  for p in n-1:-1:1
    t = squeeze(x[:,:,p], 3)
    C = t*t'
    U, s, V = svd(C)
    U = U[:,1:nvirt]
    W, s, V = svd(U'*Uprev)
    U = U*W*V' # rotate U to align with previous slice
    y[:,:,p] = U'*t
    #Uprev[:] = U # took this out bc seems unnecessary
  end
  y
end

function test_coilcompress()
  # ADJUSTABLE PARAMETERS
  golden = true # this should correspond to the sequence
  npe_per_z = 128
  sliceΔz = 64
  fcc_phase = 1.4580  # from SIN file
  fcc_scale = 1.1784
  infile = "1305_QBC_GA_12mm_15_3.7.mat"
  #infile = "/Users/dss/Dropbox/1305_QBC_GA_12mm_15_3.7.mat"
  #infile = "/Users/dss/Data/combi/CMT_XTEND_20141216/1219_QBC_GA_12mm_15_3.7_min.mat"
  infile = "/home/dss/Data/combi/CMT_XTEND_20141216/1231_QBC_GA_12mm_15_3.7_min.mat"

  # load some COMBI data
  println("reading data from $infile")
  mat = matread(infile)
  data = mat["data"]
  info = mat["info"]  # in case we want to pull params from INFO
  data = coilcompress(data)
end

function walshcombineworker{T,N}(x::Array{T,N})
  # Performs Walsh adaptive coil combination on arbitrary N-d arrays
  # Assumes:
  #   first dimension is channel dimension
  #   number of dimensions N is at least 2
  @assert N > 1 "input must have two or more dimensions"
  insize = size(x)
  nchan = insize[1]
  nvox = prod(insize[2:end])
  W = Array(T, nchan, nvox)
  tstart = time()
  y = Array(T, nvox)
  for i in 1:nvox
    #if mod(i,div(nvox,20)) == 0
    #  pct = 100i / float(nvox)
    #  telapsed = time() - tstart
    #  trem = 100telapsed/pct - telapsed
    #  @printf "channel combination %.0f%% complete, %.1f s elapsed, %.1f s remaining\n" pct telapsed trem
    #end
    C = x[:,i]*x[:,i]'
    U, s, V = svd(C)
    y[i] = (U[:,1]'*x[:,i])[1]
  end
  return y
end

walshcombine{T}(x::Array{T,2}) = walshcombineworker(x) # single threaded
walshcombine{T}(x::Array{T,3}) = walshcombineworker(x) # single threaded

function walshcombine{T}(x::Array{T,4})  # multi-threaded
  # Wraps Walsh combination code for volume data and manages parallel execution
  # Assumes:
  #   First dimension is the channel dimension.
  #   All workers are available for work.
  nchan, nx, ny, nz = size(x)
  if nchan == 1  # skip fancy combination if only one channel
    return squeeze(x[1,:,:,:], 1)
  end
  y = Array(T, nx, ny, nz)
  nw = nworkers()  # see how many workers available
  if nw > 1
    blas_set_num_threads(1) # don't compete with BLAS threading
    reflist = Any[]
    workerids = workers()
    slices = [[w:nw:nz] for w in 1:nw] # divvy up data by slices
    for w in 1:nw
      println("worker $w gets slices $(slices[w])")
      r = @spawnat workerids[w] walshcombineworker(x[:,:,:,slices[w]])
      push!(reflist, r)   # save the ref for later 'fetch'ing
    end
    y = Array(T, nx, ny, nz)
    for w in 1:nw  # for each worker
      y[:,:,slices[w]] = fetch(reflist[w])  # get result and store
    end
  else  # run single threaded
    y = walshcombine(x)
  end
  y
end

max_stars_repo_path: code/Julia Recon Scripts/coilcomb.jl

max_issues_repo_path: code/Julia Recon Scripts/coilcomb.jl

max_forks_repo_path: code/Julia Recon Scripts/coilcomb.jl


text:
# ------------------------------------
# muscl
# ------------------------------------
function muscl(Qbase, QbaseU, QbaseD, QbaseL, QbaseR, QconU, QconD, QconL, QconR,
                cellxmax, cellymax, vecAx, vecAy, nval, icell, specific_heat_ratio, volume, nodes)

    #=
    Q system
    VecA system is sama
    Q is defined on(near) boundary

            nodes[2,3]  QU[2,3]
                |--------------------|
                |       QD[2,3]      |
                |                    |
        QL[2,2] | QR[2,2]    QL[3,2] | QR[3,2]
                |      cell(2,2)     |
                |                    |
                |       QU[2,2]      |
                |--------------------|
            nodes[2,2]  QD[2,2]   nodes[3,2]

    =#
    for l in 1:nval
        for j in 1+icell-1:cellymax-icell+1
            for i in 1+icell:cellxmax-icell+1

                #vecAx_xav = 0.25*( vecAx[i,j,1] + vecAx[i-1,j,1] + vecAx[i,j+1,1] + vecAx[i-1,j+1,1] )
                #vecAx_yav = 0.25*( vecAx[i,j,2] + vecAx[i-1,j,2] + vecAx[i,j+1,2] + vecAx[i-1,j+1,2] )
                #volume_av = 0.5*( volume[i-1,j] + volume[i-1,j+1])

                # ------------------------------
                # L
                # ------------------------------
                #=
                vim = 0.5*( volume[i-2,j] + volume[i-2,j-1])
                vi  = 0.5*( volume[i-1,j] + volume[i-1,j-1])
                vip = 0.5*( volume[i,  j] + volume[i,  j-1])

                vim_p = 0.5*( volume[i-2,j+1] + volume[i-2,j])
                vi_p  = 0.5*( volume[i-1,j+1] + volume[i-1,j])
                vip_p = 0.5*( volume[i,  j+1] + volume[i,  j])

                lim = 0.5 * ((1/vecAy[i-2,  j,1]^2 + 1/vecAy[i-2,  j,2]^2)^0.5 * vim_p
                            +(1/vecAy[i-2,j-1,1]^2 + 1/vecAy[i-2,j-1,2]^2)^0.5 * vim)
                li  = 0.5 * ((1/vecAy[i-1,  j,1]^2 + 1/vecAy[i-1,  j,2]^2)^0.5 * vi_p
                            +(1/vecAy[i-1,j-1,1]^2 + 1/vecAy[i-1,j-1,2]^2)^0.5 * vi)
                lip = 0.5 * ((1/vecAy[  i,  j,1]^2 + 1/vecAy[  i,  j,2]^2)^0.5 * vip_p
                            +(1/vecAy[  i,j-1,1]^2 + 1/vecAy[  i,j-1,2]^2)^0.5 * vip)
                =#

                lim = 0.5 * (((nodes[i-2,j+1,1] - nodes[i-1,j+1,1])^2 + (nodes[i-2,j+1,2] - nodes[i-1,j+1,2])^2)^0.5
                            +((nodes[i-2,  j,1] - nodes[i-1,  j,1])^2 + (nodes[i-2,  j,2] - nodes[i-1,  j,2])^2)^0.5)
                li  = 0.5 * (((nodes[i-1,j+1,1] - nodes[  i,j+1,1])^2 + (nodes[i-1,j+1,2] - nodes[  i,j+1,2])^2)^0.5
                            +((nodes[i-1,  j,1] - nodes[  i,  j,1])^2 + (nodes[i-1,  j,2] - nodes[  i,  j,2])^2)^0.5)
                lip = 0.5 * (((nodes[  i,j+1,1] - nodes[i+1,j+1,1])^2 + (nodes[  i,j+1,2] - nodes[i+1,j+1,2])^2)^0.5
                            +((nodes[  i,  j,1] - nodes[i+1,  j,1])^2 + (nodes[  i,  j,2] - nodes[i+1,  j,2])^2)^0.5)

                dip = 2*(Qbase[i,j,l]   - Qbase[i-1,j,l]) / (li  + lip)
                dim = 2*(Qbase[i-1,j,l] - Qbase[i-2,j,l]) / (lim + li)
                QbaseL[i,j,l] = Qbase[i-1,j,l] + 0.5*li*minmod(dim, dip)

                # ------------------------------
                # R
                # ------------------------------
                #=
                vip   = 0.5*( volume[i+1,  j] + volume[i+1,j-1])
                vip_p = 0.5*( volume[i+1,j+2] + volume[i+1,  j])

                lim = li
                li  = lip
                lip = 0.5 * ((1/vecAy[i,j+2,1]^2 + 1/vecAy[i,j+2,2]^2)^0.5 * vip_p
                            +(1/vecAy[i,j+1,1]^2 + 1/vecAy[i,j+1,2]^2)^0.5 * vip)
                =#
                lim = li
                li  = lip
                lip = 0.5 * (((nodes[i+1,j+1,1] - nodes[i+2,j+1,1])^2 + (nodes[i+1,j+1,2] - nodes[i+2,j+1,2])^2)^0.5
                            +((nodes[i+1,  j,1] - nodes[i+2,  j,1])^2 + (nodes[i+1,  j,2] - nodes[i+2,  j,2])^2)^0.5)

                dip = 2*(Qbase[i+1,j,l] - Qbase[i,j,l]) / (li  + lip)
                dim = 2*(Qbase[i,j,l]   - Qbase[i-1,j,l]) / (lim + li)
                QbaseR[i,j,l] = Qbase[i,j,l]  - 0.5*li*minmod(dip, dim)

                #QbaseL[i,j,l] = Qbase[i-1,j,l]
                #QbaseR[i,j,l] = Qbase[i,j,l]
            end
        end
    end

    for l in 1:nval
        for j in 1+icell:cellymax-icell + 1
            for i in 1+icell-1:cellxmax-icell+1
                # ------------------------------
                # D
                # ------------------------------
                #=
                vim = 0.5*( volume[i,j-2] + volume[i-1,j-2])
                vi  = 0.5*( volume[i,j-1] + volume[i-1,j-1])
                vip = 0.5*( volume[i,  j] + volume[i-1,  j])

                vim_p = 0.5*( volume[i+1,j-2] + volume[i,j-2])
                vi_p  = 0.5*( volume[i+1,j-1] + volume[i,j-1])
                vip_p = 0.5*( volume[i+1,  j] + volume[i,  j])

                lim = 0.5 * ((1/vecAx[  i,j-2,1]^2 + 1/vecAx[  i,j-2,2]^2)^0.5 * vim_p
                            +(1/vecAx[i-1,j-2,1]^2 + 1/vecAx[i-1,j-2,2]^2)^0.5 * vim)
                li  = 0.5 * ((1/vecAx[  i,j-1,1]^2 + 1/vecAx[  i,j-1,2]^2)^0.5 * vi_p
                            +(1/vecAx[i-1,j-1,1]^2 + 1/vecAx[i-1,j-1,2]^2)^0.5 * vi)
                lip = 0.5 * ((1/vecAx[  i,  j,1]^2 + 1/vecAx[  i,  j,2]^2)^0.5 * vip_p
                            +(1/vecAx[i-1,  j,1]^2 + 1/vecAx[i-1,  j,2]^2)^0.5 * vip)
                =#

                lim = 0.5 * (((nodes[i+1,j-2,1] - nodes[i+1,j-1,1])^2 + (nodes[i+1,j-2,2] - nodes[i+1,j-1,2])^2)^0.5
                            +((nodes[  i,j-2,1] - nodes[  i,j-1,1])^2 + (nodes[  i,j-2,2] - nodes[  i,j-1,2])^2)^0.5)
                li  = 0.5 * (((nodes[i+1,j-1,1] - nodes[i+1,  j,1])^2 + (nodes[i+1,j-1,2] - nodes[i+1,  j,2])^2)^0.5
                            +((nodes[  i,j-1,1] - nodes[  i,  j,1])^2 + (nodes[  i,j-1,2] - nodes[  i,  j,2])^2)^0.5)
                lip = 0.5 * (((nodes[i+1,  j,1] - nodes[i+1,j+1,1])^2 + (nodes[i+1,  j,2] - nodes[i+1,j+1,2])^2)^0.5
                            +((nodes[  i,  j,1] - nodes[  i,j+1,1])^2 + (nodes[  i,  j,2] - nodes[  i,j+1,2])^2)^0.5)

                dip = 2*(Qbase[i,j,l]   - Qbase[i,j-1,l]) / (li  + lip)
                dim = 2*(Qbase[i,j-1,l] - Qbase[i,j-2,l]) / (lim + li)
                QbaseD[i,j,l] = Qbase[i,j-1,l] + 0.5*li*minmod(dim, dip)

                # ------------------------------
                # U
                # ------------------------------
                #=
                vip   = 0.5*( volume[  i,j+1] + volume[i-1,j+1])
                vip_p = 0.5*( volume[i+1,j+1] + volume[  i,j+1])

                lim = li
                li  = lip
                lip = 0.5 * ((1/vecAx[i+2,j,1]^2 + 1/vecAx[i+2,j,2]^2)^0.5 * vip_p
                            +(1/vecAx[i+1,j,1]^2 + 1/vecAx[  i,j,2]^2)^0.5 * vip)
                =#
                lim = li
                li  = lip
                lip = 0.5 * (((nodes[i+1,j+1,1] - nodes[i+1,j+2,1])^2 + (nodes[i+1,j+1,2] - nodes[i+1,j+2,2])^2)^0.5
                            +((nodes[  i,j+1,1] - nodes[  i,j+2,1])^2 + (nodes[  i,j+1,2] - nodes[  i,j+2,2])^2)^0.5)

                dip = 2*(Qbase[i,j+1,l] - Qbase[i,j,l])   / (li  + lip)
                dim = 2*(Qbase[i,j,l]   - Qbase[i,j-1,l]) / (lim + li)
                QbaseU[i,j,l] = Qbase[i,j,l]  - 0.5*li*minmod(dip, dim)

                #QbaseD[i,j,l] = Qbase[i,j-1,l]
                #QbaseU[i,j,l] = Qbase[i,j,l]
            end
        end
    end

    QconU = base_to_conservative(QbaseU, QconU, cellxmax, cellymax, specific_heat_ratio)
    QconD = base_to_conservative(QbaseD, QconD, cellxmax, cellymax, specific_heat_ratio)
    QconL = base_to_conservative(QbaseL, QconL, cellxmax, cellymax, specific_heat_ratio)
    QconR = base_to_conservative(QbaseR, QconR, cellxmax, cellymax, specific_heat_ratio)

    #=
    println("")
    println(Qbase[49,50,2])
    println(QbaseL[50,50,2])
    println(Qbase[50,50,2])
    =#

    return QbaseU, QbaseD, QbaseL, QbaseR, QconU, QconD, QconL, QconR
end

function minmod(a,b)
    ans = sign(a) * max(0, min(a*sign(b), b*sign(a)))
    return ans
end

ext: jl


text:
###### This file is part of the MendelImpute.jl package.
###### It contains relevant code to solve the least squares criteria
###### `x ≈ hᵢ + hⱼ`.

"""
    compute_optimal_haplotypes!(haplotype1, haplotype2, compressed_Hunique, ...)

Computes the best haplotype pair `(hᵢ, hⱼ)` for each genotype vector `x` and
stores result in `haplotype1` and `haplotype2`.

# Arguments
- `haplotype1`: Person `i` strand1 haplotype in window `w` is `haplotype1[i][w]`
- `haplotype2`: Person `i` strand2 haplotype in window `w` is `haplotype2[i][w]`
- `compressed_Hunique`: A `CompressedHaplotypes` object
- `X`: the full genotype matrix possibly with missings. Each column is an
    individual.
- `X_pos`: Position of each SNP in `X`
- `typed_snps_error`: `typed_snps_error[i]` is the % of SNPs correctly imputed.
    `1` is best, `0` is worse.
- `stepscreen`: Boolean indicating whether to use a stepwise heuristic to screen
    for top haplotypes instead of performing global search
- `tf`: This option solves the least squares objective on only
    `tf` unique haplotypes.
- `scale_allelefreq` Boolean indicating whether to give rare SNPs more weight
    scaled by `wᵢ = 1 / √2p(1-p)` where max weight is 2.
- `max_haplotypes` Maximum number of haplotypes for using to global search.
    This number should be specified along with `stepscreen` or `tf`.
- `rescreen` Boolean for a more rigorous global search.
- `timers`

# Timers:
- `t1` = screening for top haplotypes (0 unless `tf` option is on)
- `t2` = BLAS3 mul! to get M and N
- `t3` = haplopair search
- `t4` = rescreen time (0 unless `rescreen = true`)
- `t5` = initializing missing
- `t6` = allocating and viewing
- `t7` = index conversion
"""
function compute_optimal_haplotypes!(
    haplotype1::AbstractVector,
    haplotype2::AbstractVector,
    compressed_Hunique::CompressedHaplotypes,
    X::AbstractMatrix,
    X_pos::AbstractVector,
    stepscreen::Union{Nothing, Int},
    tf::Union{Nothing, Int}, # thinning factor
    scale_allelefreq::Bool,
    max_haplotypes::Int,
    rescreen::Bool
    )
    # constants
    people = size(X, 2)
    ref_snps = length(compressed_Hunique.pos)
    max_width, max_d = MendelImpute.max_dim(compressed_Hunique)
    windows = length(haplotype1[1])
    threads = Threads.nthreads()
    inv_sqrt_allele_var = nothing
    winranges = compressed_Hunique.X_window_range
    overlap = compressed_Hunique.overlap

    # allocate working arrays
    timers = [zeros(7*8) for _ in 1:threads] # 8 for spacing
    timers[1][48] += @elapsed begin # time for allocating
        happair1 = [ones(Int32, people)               for _ in 1:threads]
        happair2 = [ones(Int32, people)               for _ in 1:threads]
        hapscore = [zeros(Float32, people)            for _ in 1:threads]
        Xwork    = [zeros(Float32, max_width, people) for _ in 1:threads]
        Hwork    = [zeros(Float32, max_width, max_d)  for _ in 1:threads]
        if !isnothing(tf)
            maxindx = [zeros(Int32, tf)         for _ in 1:threads]
            maxgrad = [zeros(Float32, tf)       for _ in 1:threads]
            Hk = [zeros(Float32, max_width, tf) for _ in 1:threads]
            Xi = [zeros(Float32, max_width)     for _ in 1:threads]
            M  = [zeros(Float32, tf, tf)        for _ in 1:threads]
            N  = [zeros(Float32, tf)            for _ in 1:threads]
            R  = [zeros(Float32, max_d, people) for _ in 1:threads]
        elseif !isnothing(stepscreen)
            maxindx = [zeros(Int32,   stepscreen) for _ in 1:threads]
            maxgrad = [zeros(Float32, stepscreen) for _ in 1:threads]
            M  = [zeros(Float32, max_d, max_d)    for _ in 1:threads]
            Nt = [zeros(Float32, max_d, people)   for _ in 1:threads]
        else
            M = [zeros(Float32, max_d,  max_d) for _ in 1:threads]
            N = [zeros(Float32, people, max_d) for _ in 1:threads]
        end
    end
    pmeter = Progress(windows, 5, "Computing optimal haplotypes...")

    # for w in 1:windows
    Threads.@threads for w in 1:windows
        id = Threads.threadid()
        t6 = @elapsed begin
            Hw_aligned = compressed_Hunique.CW_typed[w].uniqueH
            Xrange = extend_to_overlap_range(compressed_Hunique, w, overlap)
            Xw_aligned = view(X, Xrange, :)
            d = size(Hw_aligned, 2)
        end

        # weight snp by inverse allele variance if requested
        t6 += @elapsed if scale_allelefreq
            Hw_range = compressed_Hunique.Hstart[w]:(w ==
                windows ? ref_snps : compressed_Hunique.Hstart[w + 1] - 1)
            Hw_snp_pos = indexin(X_pos[winranges[w]],
                compressed_Hunique.pos[Hw_range])
            inv_sqrt_allele_var = compressed_Hunique.altfreq[Hw_snp_pos]
            map!(x -> x < 0.15 ? 1.98 : 1 / sqrt(2*x*(1-x)),
                inv_sqrt_allele_var, inv_sqrt_allele_var) # set min pᵢ = 0.15
        end

        # compute top haplotype pairs for each sample in current window
        if !isnothing(stepscreen) && d > max_haplotypes
            # find hᵢ via stepwise regression, then find hⱼ via global search
            t1, t2, t3, t4, t5, t6 = haplopair_stepscreen!(Xw_aligned,
                Hw_aligned, r=stepscreen,
                inv_sqrt_allele_var=inv_sqrt_allele_var, happair1=happair1[id],
                happair2=happair2[id], hapscore=hapscore[id],
                maxindx=maxindx[id], maxgrad=maxgrad[id], Xwork=Xwork[id],
                Hwork=Hwork[id], Nt=Nt[id], M=M[id])
        elseif !isnothing(tf) && d > max_haplotypes
            # haplotype thinning: search all (hᵢ, hⱼ) pairs where hᵢ ≈ x ≈ hⱼ
            t1, t2, t3, t4, t5, t6 = haplopair_thin_BLAS2!(Xw_aligned,
                Hw_aligned, allele_freq=inv_sqrt_allele_var, keep=tf,
                happair1=happair1[id], happair2=happair2[id],
                hapscore=hapscore[id], maxindx=maxindx[id], maxgrad=maxgrad[id],
                Xi=Xi[id], N=N[id], Hk=Hk[id], M=M[id], Xwork=Xwork[id],
                Hwork=Hwork[id], R=R[id])
        elseif rescreen
            # global search + searching ||x - hᵢ - hⱼ|| on observed entries
            t1, t2, t3, t4, t5, t6 = haplopair_rescreen!(Xw_aligned,
                Hw_aligned, happair1=happair1[id], happair2=happair2[id],
                hapscore=hapscore[id], Xwork=Xwork[id])
        else
            # global search
            t1, t2, t3, t4, t5, t6 = haplopair!(Xw_aligned, Hw_aligned,
                inv_sqrt_allele_var=inv_sqrt_allele_var, happair1=happair1[id],
                happair2=happair2[id], hapscore=hapscore[id], Xwork=Xwork[id],
                Hwork=Hwork[id], M=M[id], N = N[id])
        end

        # save result
        t7 = @elapsed save_haplotypes!(haplotype1, haplotype2, happair1[id],
            happair2[id], compressed_Hunique, w)

        # record timings and haplotypes (× 8 to avoid false sharing)
        timers[id][8]  += t1
        timers[id][16] += t2
        timers[id][24] += t3
        timers[id][32] += t4
        timers[id][40] += t5
        timers[id][48] += t6
        timers[id][56] += t7

        # update progress
        next!(pmeter)
    end

    return sum(timers) ./ threads
end

"""
    save_haplotypes!(haplotype1, haplotype2, happair1, happair2, ...)

Helper function to convert `happair`s (which index off unique haplotypes) to
indices of full haplotype pool, and store them in `haplotype`s.

# Arguments
- `haplotype1`: Person `i` strand1 haplotype in window `w` is `haplotype1[i][w]`
- `haplotype2`: Person `i` strand2 haplotype in window `w` is `haplotype2[i][w]`
- `happair1`: Optimal haplotype pair in strand1 of current window, indexes off
    of unique haplotypes.
- `happair2`: Optimal haplotype pair in strand2 of current window, indexes off
    of unique haplotypes.
- `compressed_Hunique`: A `CompressedHaplotypes` object
- `window` current window.
"""
function save_haplotypes!(
    haplotype1::Vector{Vector{Int32}},
    haplotype2::Vector{Vector{Int32}},
    happair1::Vector{Int32},
    happair2::Vector{Int32},
    compressed_Hunique::CompressedHaplotypes,
    window::Int,
    )

    people = length(haplotype1)
    @inbounds for i in 1:people
        haplotype1[i][window] = unique_idx_to_complete_idx(
            happair1[i], window, compressed_Hunique)
        haplotype2[i][window] = unique_idx_to_complete_idx(
            happair2[i], window, compressed_Hunique)
    end
    return nothing
end

# """
# Records optimal-redundant haplotypes for each window.

# Warning: This function is called in a multithreaded loop. If you modify this
# function you must check whether imputation accuracy is affected (when run with
# >1 threads).

# # Arguments:
# - `window_idx`: window in current chunk
# - `window_overall`: window index in terms of every windows
# """
# function compute_redundant_haplotypes!(
#     redundant_haplotypes::Vector{OptimalHaplotypeSet},
#     Hunique::CompressedHaplotypes,
#     happair1::AbstractVector,
#     happair2::AbstractVector,
#     window_idx::Int,
#     window_overall::Int,
#     storage1 = falses(nhaplotypes(Hunique)),
#     storage2 = falses(nhaplotypes(Hunique))
#     )

#     people = length(redundant_haplotypes)

#     @inbounds for k in 1:people
#         # convert happairs from unique idx to complete idx
#         Hi_idx = unique_idx_to_complete_idx(happair1[k], window_overall,
#             Hunique)
#         Hj_idx = unique_idx_to_complete_idx(happair2[k], window_overall,
#             Hunique)

#         # strand1
#         storage1 .= false
#         if haskey(Hunique.CW_typed[window_overall].hapmap, Hi_idx)
#             h1_set = Hunique.CW_typed[window_overall].hapmap[Hi_idx]
#             for i in h1_set
#                 storage1[i] = true
#             end
#         else
#             storage1[Hi_idx] = true # Hi_idx is singleton (i.e. unique)
#         end

#         # strand2
#         storage2 .= false
#         if haskey(Hunique.CW_typed[window_overall].hapmap, Hj_idx)
#             h2_set = Hunique.CW_typed[window_overall].hapmap[Hj_idx]
#             for i in h2_set
#                 storage2[i] = true
#             end
#         else
#             storage2[Hj_idx] = true # Hj_idx is singleton (i.e. unique)
#         end

#         # redundant_haplotypes[k].strand1[window_idx] = copy(storage1)
#         # redundant_haplotypes[k].strand2[window_idx] = copy(storage2)
#         if isassigned(redundant_haplotypes[k].strand1, window_idx)
#             redundant_haplotypes[k].strand1[window_idx] .= storage1
#             redundant_haplotypes[k].strand2[window_idx] .= storage2
#         else
#             redundant_haplotypes[k].strand1[window_idx] = copy(storage1)
#             redundant_haplotypes[k].strand2[window_idx] = copy(storage2)
#         end
#     end

#     return nothing
# end

"""
    screen_flanking_windows!(haplotype1, haplotype2, compressed_Hunique, X)

For each window's haplotype pair, tests whether adjacent window's haplotype
pairs produce better error on the observed entries.

# Arguments
- `haplotype1`: Person `i` strand1 haplotype in window `w` is `haplotype1[i][w]`
- `haplotype2`: Person `i` strand2 haplotype in window `w` is `haplotype2[i][w]`
- `compressed_Hunique`: A `CompressedHaplotypes` object
- `X`: the full genotype matrix possibly with missings. Each column is an
    individual.
"""
function screen_flanking_windows!(
    haplotype1::AbstractVector,
    haplotype2::AbstractVector,
    compressed_Hunique::CompressedHaplotypes,
    X::AbstractMatrix
    )

    people = length(haplotype1)
    haplotypes = nhaplotypes(compressed_Hunique)
    width = compressed_Hunique.width
    windows = length(haplotype1[1])
    winranges = compressed_Hunique.X_window_range

    for w in 1:windows
        Hw_aligned = compressed_Hunique.CW_typed[w].uniqueH
        Xw_aligned = view(X, winranges[w], :)

        for i in 1:people
            # calculate observed error for current pair
            h1_curr_complete = haplotype1[i][w] # complete index
            h2_curr_complete = haplotype2[i][w] # complete index
            h1_curr = complete_idx_to_unique_typed_idx(h1_curr_complete, w,
                compressed_Hunique) # unique index
            h2_curr = complete_idx_to_unique_typed_idx(h2_curr_complete, w,
                compressed_Hunique) # unique index
            curr_err = observed_error(Xw_aligned, i, Hw_aligned, h1_curr,
                h2_curr)

            # consider previous pair
            if w != 1
                h1_prev_complete = haplotype1[i][w - 1]
                h2_prev_complete = haplotype2[i][w - 1]
                h1_prev = complete_idx_to_unique_typed_idx(h1_prev_complete,
                    w, compressed_Hunique)
                h2_prev = complete_idx_to_unique_typed_idx(h2_prev_complete,
                    w, compressed_Hunique)
                prev_err = observed_error(Xw_aligned, i, Hw_aligned, h1_prev,
                    h2_prev)
                if prev_err < curr_err
                    h1_curr, h2_curr, curr_err = h1_prev, h2_prev, prev_err
                end
            end

            # consider next pair
            if w != windows
                h1_next_complete = haplotype1[i][w + 1]
                h2_next_complete = haplotype2[i][w + 1]
                h1_next = complete_idx_to_unique_typed_idx(h1_next_complete,
                    w, compressed_Hunique) # unique index
                h2_next = complete_idx_to_unique_typed_idx(h2_next_complete,
                    w, compressed_Hunique) # unique index
                next_err = observed_error(Xw_aligned, i, Hw_aligned, h1_next,
                    h2_next)
                if next_err < curr_err
                    h1_curr, h2_curr, curr_err = h1_next, h2_next, next_err
                end
            end

            # convert from unique idx back to complete idx
            H1_idx = unique_idx_to_complete_idx(h1_curr, w, compressed_Hunique)
            H2_idx = unique_idx_to_complete_idx(h2_curr, w, compressed_Hunique)
            haplotype1[i][w] = H1_idx
            haplotype2[i][w] = H2_idx
        end
    end

    return nothing
end

# uses dynamic programming. Only the first 1000 haplotype pairs will be saved.
# function compute_redundant_haplotypes!(
#     redundant_haplotypes::Vector{Vector{Vector{T}}},
#     Hunique::CompressedHaplotypes,
#     happair1::AbstractVector,
#     happair2::AbstractVector,
#     window_idx::Int,
#     window_overall::Int,
#     storage1 = falses(nhaplotypes(Hunique)),
#     storage2 = falses(nhaplotypes(Hunique))
#     ) where T <: Tuple{Int32, Int32}

#     people = length(redundant_haplotypes)

#     @inbounds for k in 1:people
#         # convert happairs from unique idx to complete idx
#         Hi_idx = unique_idx_to_complete_idx(happair1[k], window_overall, Hunique)
#         Hj_idx = unique_idx_to_complete_idx(happair2[k], window_overall, Hunique)

#         # find haplotypes that match Hi_idx and Hj_idx on typed snps
#         h1_set = get(Hunique.CW_typed[window_overall].hapmap, Hi_idx, Hi_idx)
#         h2_set = get(Hunique.CW_typed[window_overall].hapmap, Hj_idx, Hj_idx)

#         # save first 1000 haplotype pairs
#         for h1 in h1_set, h2 in h2_set
#             if length(redundant_haplotypes[k][window_idx]) < 1000
#                 push!(redundant_haplotypes[k][window_idx], (h1, h2))
#             else
#                 break
#             end
#         end
#     end

#     return nothing
# end

"""
    haplopair(X, H)

Calculate the best pair of haplotypes in `H` for each individual in `X`.
Missing data is initialized as 2x alternate allele freq.

# Input
* `X`: `p × n` genotype matrix possibly with missings. Each column is an individual.
* `H`: `p × d` haplotype matrix. Each column is a haplotype.

# Optional arguments
* `happair1`: optimal haplotype pairs for strand 1.
* `happair2`: optimal haplotype pairs for strand 2.
* `hapscore`: haplotyping score. 0 means best. Larger means worse.
* `inv_sqrt_allele_var`: If not nothing, SNP `i` will be scaled according to
    `inv_sqrt_allele_var[i]`.
* `Xwork`: A `Float32` version of `X`, missing entries are filled with `2maf`
* `Hwork`: A `Float32` version of `H`.
* `M`: `d x d` matrix with entries `M[i, j] = 2dot(Hwork[:, i], Hwork[:, j]) +
    sumabs2(Hwork[:, i]) + sumabs2(Hwork[:, j])`. Only the upper triangular
    part of `M` is used.
* `N`: `n x d` matrix `2Xwork'Hwork`.
"""
function haplopair!(
    X::AbstractMatrix, # p × n
    H::AbstractMatrix; # p × d
    # preallocated vectors
    happair1::AbstractVector = ones(Int32, size(X, 2)), # length n
    happair2::AbstractVector = ones(Int32, size(X, 2)), # length n
    hapscore::AbstractVector = Vector{Float32}(undef, size(X, 2)), # length n
    inv_sqrt_allele_var::Union{Nothing, AbstractVector} = nothing, # length p
    # preallocated matrices
    Xwork::AbstractMatrix{Float32}=Matrix{Float32}(undef,size(X,1),size(X,2)), # p × n
    Hwork::AbstractMatrix{Float32}=convert(Matrix{Float32}, H),                # p × d
    M::AbstractMatrix{Float32}=Matrix{Float32}(undef, size(H, 2), size(H, 2)), # d × d
    N::AbstractMatrix{Float32}=Matrix{Float32}(undef, size(X, 2), size(H, 2)), # n × d
    )
    p,  n = size(X)
    pp, d = size(H)
    p == pp || error("size of H doesn't match X! Expected $p got $pp")

    # create views for matrices
    t6 = @elapsed begin
        Mwork_view = view(M, 1:d, 1:d)
        Nwork_view = view(N, :, 1:d)
        Hwork_view = view(Hwork, 1:p, 1:d)
        Xwork_view = view(Xwork, 1:p, :)
        copyto!(Hwork_view, H) # update Hwork_view, others are updated later
    end

    # initializes missing
    t5 = @elapsed initXfloat!(Xwork_view, X)

    t2, t3 = haplopair!(Xwork_view, Hwork_view, Mwork_view, Nwork_view,
        happair1, happair2, hapscore, inv_sqrt_allele_var)
    t1 = t4 = 0.0 # no time spent on haplotype thinning or rescreening

    return t1, t2, t3, t4, t5, t6
end

"""
    haplopair!(X, H, M, N, happair, hapscore)

Calculate the best pair of haplotypes in `H` for each individual in `X`.
Overwite `M` by `M[i, j] = 2dot(H[:, i], H[:, j]) + sumabs2(H[:, i]) +
sumabs2(H[:, j])`, `N` by `2X'H`, `happair` by optimal haplotype pair, and
`hapscore` by objective value from the optimal haplotype pair.

# Input
* `X`: `p x n` genotype matrix. Each column is an individual.
* `H`: `p x d` haplotype matrix. Each column is a haplotype.
* `M`: overwritten by `M[i, j] = 2dot(H[:, i], H[:, j]) + sumabs2(H[:, i]) +
    sumabs2(H[:, j])`.
* `N`: overwritten by `n x d` matrix `2X'H`.
* `happair`: optimal haplotype pair. `X[:, k] ≈ H[:, happair[k, 1]] +
    H[:, happair[k, 2]]`.
* `hapscore`: haplotyping score. 0 means best. Larger means worse.
"""
function haplopair!(
    X::AbstractMatrix{Float32},
    H::AbstractMatrix{Float32},
    M::AbstractMatrix{Float32},
    N::AbstractMatrix{Float32},
    happair1::AbstractVector{Int32},
    happair2::AbstractVector{Int32},
    hapscore::AbstractVector{Float32},
    inv_sqrt_allele_var::Union{Nothing, AbstractVector}
    )

    p, n, d = size(X, 1), size(X, 2), size(H, 2)

    # assemble M (upper triangular only)
    t2 = @elapsed begin
        if !isnothing(inv_sqrt_allele_var)
            H .*= inv_sqrt_allele_var # wᵢ = 1/√2p(1-p)
        end
        mul!(M, Transpose(H), H)
        for j in 1:d, i in 1:(j - 1) # off-diagonal
            @inbounds M[i, j] = 2M[i, j] + M[i, i] + M[j, j]
        end
        for j in 1:d # diagonal
            @inbounds M[j, j] *= 4
        end

        # assemble N
        if !isnothing(inv_sqrt_allele_var)
            H .*= inv_sqrt_allele_var # wᵢ = 1/2p(1-p)
        end
        mul!(N, Transpose(X), H)
        @simd for I in eachindex(N)
            N[I] *= 2
        end
    end

    # computational routine
    t3 = @elapsed haplopair!(happair1, happair2, hapscore, M, N)

    # supplement the constant terms in objective
    t3 += @elapsed begin
        @inbounds for j in 1:n
            @simd for i in 1:p
                hapscore[j] += abs2(X[i, j])
            end
        end
    end

    return t2, t3
end

"""
    haplopair!(happair, hapscore, M, N)

Calculate the best pair of haplotypes pairs in the filtered haplotype panel
for each individual in `X` using sufficient statistics `M` and `N`.

# Note
The best haplotype pairs are column indices of the filtered haplotype panels.

# Input
* `happair1`: optimal haplotype in strand 1 for each individual.
* `happair2`: optimal haplotype in strand 2 for each individual.
* `hapmin`: minimum offered by the optimal haplotype pair.
* `M`: `d x d` matrix with entries `M[i, j] = 2dot(H[:, i], H[:, j]) +
    sumabs2(H[:, i]) + sumabs2(H[:, j])`, where `H` is the haplotype matrix
    with haplotypes in columns. Only the upper triangular part of `M` is used.
* `N`: `n x d` matrix `2X'H`, where `X` is the genotype matrix with individuals
    in columns.
"""
function haplopair!(
    happair1::AbstractVector{Int32},
    happair2::AbstractVector{Int32},
    hapmin::AbstractVector{Float32},
    M::AbstractMatrix{Float32},
    N::AbstractMatrix{Float32},
    )

    n, d = size(N)
    chunks = div(n, 4)
    fill!(hapmin, Inf32)
    @inbounds for k in 1:d, j in 1:k
        Mjk = M[j, k]
        i = 1
        # loop over individuals
        for chunk in 1:chunks # unrolled loop, equivalent to for i in 1:n
            i2, i3, i4 = i + 1, i + 2, i + 3
            score = Mjk - N[i, j] - N[i, k]
            if score < hapmin[i]
                hapmin[i], happair1[i], happair2[i] = score, j, k
            end
            score = Mjk - N[i2, j] - N[i2, k]
            if score < hapmin[i2]
                hapmin[i2], happair1[i2], happair2[i2] = score, j, k
            end
            score = Mjk - N[i3, j] - N[i3, k]
            if score < hapmin[i3]
                hapmin[i3], happair1[i3], happair2[i3] = score, j, k
            end
            score = Mjk - N[i4, j] - N[i4, k]
            if score < hapmin[i4]
                hapmin[i4], happair1[i4], happair2[i4] = score, j, k
            end
            i += 4
        end
        # handle remaining terms
        while i ≤ n
            score = Mjk - N[i, j] - N[i, k]
            if score < hapmin[i]
                hapmin[i], happair1[i], happair2[i] = score, j, k
            end
            i += 1
        end
    end
end

"""
    fillmissing!(Xm, Xwork, H, haplopairs)

Fill in missing genotypes in `X` according to haplotypes. Non-missing genotypes
remain same.

# Input
* `Xm`: `p x n` genotype matrix with missing values. Each column is an
    individual.
* `Xwork`: `p x n` genotype matrix where missing values are filled with sum of 2
    haplotypes.
* `H`: `p x d` haplotype matrix. Each column is a haplotype.
* `happair`: pair of haplotypes. `X[:, k] = H[:, happair[1][k]] +
    H[:, happair[2][k]]`.
"""
function fillmissing!(
    Xm::AbstractMatrix{Union{U, Missing}},
    Xwork::AbstractMatrix{T},
    H::AbstractMatrix{T},
    happairs::Vector{Vector{Tuple{Int32, Int32}}},
    ) where {T, U}

    p, n = size(Xm)
    best_discrepancy = typemax(eltype(Xwork))

    for j in 1:n, happair in happairs[j]
        discrepancy = zero(T)
        for i in 1:p
            if ismissing(Xm[i, j])
                tmp = H[i, happair[1]] + H[i, happair[2]]
                discrepancy += abs2(Xwork[i, j] - tmp)
                Xwork[i, j] = tmp
            end
        end
        if discrepancy < best_discrepancy
            best_discrepancy = discrepancy
        end
    end
    return best_discrepancy
end

"""
    initXfloat!(Xfloat, X)

Initializes the matrix `Xfloat` where missing values of matrix `X` by `2 x`
allele frequency and nonmissing entries of `X` are converted to type `Float32`
for subsequent BLAS routines.

# Input
* `X` is a `p x n` genotype matrix. Each column is an individual.
* `Xfloat` is the `p x n` matrix of X where missing values are filled by 2x
    allele frequency.
"""
function initXfloat!(
    Xfloat::AbstractMatrix,
    X::AbstractMatrix
    )

    T = Float32
    p, n = size(X)

    @inbounds for i in 1:p
        # allele frequency
        cnnz = zero(T)
        csum = zero(T)
        for j in 1:n
            if !ismissing(X[i, j])
                cnnz += one(T)
                csum += convert(T, X[i, j])
            end
        end
        # set missing values to 2freq, unless cnnz is 0
        imp = (cnnz == 0 ? zero(T) : csum / cnnz)
        for j in 1:n
            if ismissing(X[i, j])
                Xfloat[i, j] = imp
            else
                Xfloat[i, j] = convert(T, X[i, j])
            end
        end
    end

    any(isnan, Xfloat) && error("Xfloat has NaN during initialization!")
    any(isinf, Xfloat) && error("Xfloat has Inf during initialization!")
    any(ismissing, Xfloat) && error("Xfloat has Missing during initialization!")

    return nothing
end

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
## lightweight IntervalArithemetic over BigFloat following paper
## (https://arxiv.org/pdf/1104.1362.pdf) on root refinement
## https://github.com/JuliaIntervals/IntervalArithmetic.jl is the
## proper source for such calculations
struct 𝐈{T}
    lo::T
    hi::T
end

## An interval can be iterated over to return the end points
Base.length(I::𝐈) = 2
function Base.iterate(I::𝐈, state=nothing)
    isnothing(state) && return (I.lo, 1)
    isone(state) && return (I.hi, 2)
    return nothing
end

function Base.show(io::IO, I::𝐈)
    a, b = I
    L = maximum(precision, (a,b))
    w = b-a

    w <= 0 && (print(io, "{$a}"); return nothing)
    L′ = max(0, ceil(Int, -log2(b-a))) + 5
    ds = ceil(Int, -log10(b-a)) + 5

    sbs = ("₋","","","₀","₁","₂","₃","₄","₅","₆","₇","₈","₉")
    iob = IOBuffer()
    for i in string(L)
        print(iob, sbs[Int(i)-44])
    end
    ind = String(take!(iob))

    setprecision(L′) do
        print(io, "[$(round(a, digits=ds))…, $(round(b, digits=ds))…]$ind")
    end
    nothing
end

Base.diff(a::𝐈) = a.hi-a.lo
midpoint(a::𝐈) = a.lo + (a.hi - a.lo)/2
Base.sign(a::𝐈) = a.hi < 0 ? -1 : a.lo > 0 ? 1 : 0
minabs(a::𝐈) = a.hi < 0 ? -a.hi : (a.lo > 0) ? a.lo : zero(a.lo)
maxabs(a::𝐈) = max(abs(a.lo), abs(a.hi))
Base.isinf(a::𝐈) = isinf(a.lo) || isinf(a.hi)
Base.isnan(a::𝐈) = isnan(a.lo) || isnan(a.hi)
function Base.zero(a::𝐈)
    L = maximum(precision, (a.lo, a.hi))
    ball(0,L)
end
zero!(a::𝐈) = (zero!(a.lo); zero!(a.hi); nothing)
precision(I::𝐈) = maximum(precision, I)

function ball!(𝐱::𝐈, x, ϵ::BigFloat)
    x′ = BigFloat(x, precision=precision(ϵ)) #deepcopy(x)
    set!(𝐱.lo, x′, RoundDown)
    rsub!(𝐱.lo, ϵ, RoundDown)
    set!(𝐱.hi, x′, RoundUp)
    radd!(𝐱.hi, ϵ, RoundUp)
    𝐱
end

function ball!(𝐱, x, L::Int)
    ϵ = epsᴸ(L)
    ball!(𝐱, x, ϵ)
end

function ball(x, ϵ::BigFloat)
    𝐱 = ball(precision(ϵ))
    ball!(𝐱, x, ϵ)
end

function ball(x,L::Int)
    𝐱 = ball(L)
    ball!(𝐱, x, L)
end

function ball(x::Int,L::Int)
    T = BigFloat
    return 𝐈(T(x, precision=L), T(x, precision=L))
    𝐱 = ball(L)
    x′ = T(x, precision=L) #deepcopy(x)
    ball!(𝐱, x′, L)
    𝐱
end

ball(L::Int) = 𝐈(BigFloat(;precision=L), BigFloat(;precision=L))
ball(ϵ::BigFloat) = I(BigFloat(;precision=precision(ϵ)), BigFloat(;precision=precision(ϵ)))

Base.one(::Type{T}) where {T <: 𝐈} = 𝐈(BigFloat("1"), BigFloat("1"))

function Base.copy!(a::𝐈, b::𝐈)
    set!(a.lo, b.lo)
    set!(a.hi, b.hi)
end

function swap!(a::T, b::T) where {T <: 𝐈}
    swap!(a.lo, b.lo)
    swap!(a.hi, b.hi)
end

function Base.:+(a::𝐈, b::𝐈)
    c = deepcopy(a)
    radd!(c, b)
    c
end

function radd!(a::𝐈, b::𝐈)
    d₁,u₁ = a.lo, a.hi
    d₂,u₂ = b.lo, b.hi
    T = BigFloat
    radd!(d₁, d₂, RoundDown)
    radd!(u₁, u₂, RoundUp)
end

function radd!(a::𝐈, b::Number)
    T = BigFloat
    radd!(a.lo, b, RoundDown)
    radd!(a.hi, b, RoundUp)
end

function Base.:-(a::𝐈)
    d₁,u₁ = a.lo, a.hi
    𝐈(-a.hi, -a.lo)
end

Base.:-(a::𝐈, b::𝐈) = a + (-b)

function p(a::𝐈, b::𝐈)
    down = d1d2, u1u2, u1d2, d1u2
    up   = d1d2, u1u2, u1d2, d1u2

    d2d1
end

# right multiplication a <- a*b;
# needs three BigFloat numbers to be non-allocating
# working, d, u make non-allocating
function LinearAlgebra.rmul!(a::𝐈, b::𝐈, working′=nothing, d′=nothing, u′=nothing) # working a bigplot
    d₁,u₁ = a.lo, a.hi
    d₂,u₂ = b.lo, b.hi

    L = maximum(precision, (d₁,u₁,d₂,u₂))
    working = working′ != nothing ? working′ : BigFloat(;precision=L)
    d = d′ == nothing ? BigFloat(;precision=L) : d′
    u = u′ == nothing ? BigFloat(;precision=L) : u′

    set!(d, d₁); set!(u, d₁)
    rmul!(d, d₂, RoundDown); rmul!(u, d₂, RoundUp)

    set!(working, u₁)
    rmul!(working, u₂, RoundDown)
    min!(d, working)

    set!(working, u₁)
    rmul!(working, d₂, RoundDown)
    min!(d, working)

    set!(working, d₁)
    rmul!(working, u₂, RoundDown)
    min!(d, working)

    set!(working, d₁)
    rmul!(working, u₂, RoundUp)
    max!(u, working)

    set!(working, u₁)
    rmul!(working, d₂, RoundUp)
    max!(u, working)

    set!(working, u₁)
    rmul!(working, u₂, RoundUp)
    max!(u, working)

    set!(d₁,d)
    set!(u₁,u)

    a
end

function LinearAlgebra.rmul!(a::𝐈, b::Number)
    T = BigFloat
    rmul!(a.lo, b, RoundDown)
    rmul!(a.hi, b, RoundUp)

    if b < 0
        tmp = a.lo
        set!(a.lo, a.hi)
        set!(a.hi, tmp)
    end

    a
end

function inv!(b::𝐈)
    T = BigFloat
    c = deepcopy(b.lo)
    L = precision(c)
    set!(b.lo, rdiv!(T(1,precision=L), b.hi, RoundDown))
    set!(b.hi, rdiv!(T(1,precision=L), c, RoundUp))
    b
end

function LinearAlgebra.rdiv!(a::𝐈, b::𝐈)
    c = deepcopy(b)
    inv!(c)
    rmul!(a, c)
end

function abs!(a::𝐈)
    if a.hi < 0
        c = deepcopy(a.lo)
        swap!(a.lo, -a.hi)
        swap!(a.hi, -c)
    end
end

# muladd with extra storage to pass to `rmul!`
function muladd!(a::𝐈, b::𝐈, c::𝐈, w=BigFloat(), d=BigFloat(), u=BigFloat())
    rmul!(a, b, w, d, u)
    radd!(a, c)
    nothing
end

# 0 is uncertain
maybe0(x::𝐈) = x.lo < zero(x.lo) < x.hi
function maybe0(x::𝐈, L::Int)
    setprecision(L) do
        ϵ = eps(BigFloat)
        return !(x.hi < -ϵ || x.lo > ϵ)
        x.lo < -ϵ && x.hi > ϵ
    end
end


text:
"""
Package: SQLdf

    sqldf(query::String)::DataFrame

Execute R sqldf and return a julia DataFrame.

Columns in the DataFrame must have a type other than Any. In order to work with dates expressions like
\"""select strftime("%Y", datetime_column, "unixepoch") as year from T\""" may be used.

# Arguments
`query`: SQL query handled by R sqldf

# Returns
Julia DataFrame with the results from R sqldf

# Examples
```julia-repl
julia> T = DataFrame(C1 = [1,2], C2 = ["a","b"])

julia> query = \"""
               select *
               from T
               where C2 = "a"
               \""";

julia> sqldf(query)
1×2 DataFrame
 Row │ C1     C2
     │ Int64  String
─────┼───────────────
   1 │     1  a

julia> @sqldf query
1×2 DataFrame
 Row │ C1     C2
     │ Int64  String
─────┼───────────────
   1 │     1  a

```
"""
function sqldf(query::String)::DataFrame

    # Normalize Query
    query = replace(query,"\n"=>" ")
    nq = split(query," ")
    nq = nq[nq .!= ""]

    # Extract tables in query
    tables = String[]
    nt = false
    for w in nq
        nt = nt ? push!(tables,w)==0 : lowercase(w) in ["from","join"]
    end

    # Prepare R
    for t in tables
        Main.eval(Main.Meta.parse("@rput "*t))
    end

    # Retrieve and Return
    rquery = "c7e881694dfe63f2 = sqldf('"*query*"')"
    rquery = replace(rquery,"\"" => "\\\"")

    eval(Meta.parse("R\""*rquery*"\""))
    R"colnames(c7e881694dfe63f2) <- make.unique(colnames(c7e881694dfe63f2),sep='_')"
    res = rcopy(R"c7e881694dfe63f2")

    # Clean Up
    R"rm(c7e881694dfe63f2)"
    for t in tables
        eval(Meta.parse("R"*"\"rm("*t*")\""))
    end
    R"gc()"

    return res
end

macro sqldf(query)
    return :( sqldf($(esc(query))) )
end


text:
using PQEcolaPoint

# setup DTE -
T = 72.0
DTE = range(T, stop=0.0,step=-1.0) |> collect

# setup underlying range -
Sₒ = 79.74
IV = 0.4456
σ = Sₒ*(IV)*(√(T/365))
L = round(Sₒ - 2.96*σ)
U = round(Sₒ + 2.96*σ)
S = range(L,stop=U, step=0.5) |> collect

# setup IV array -
IV_array = range(0.1, stop = 1.2, step = 0.05) |> collect

# setup strike coordinates -
K = range(60.0, stop = 125.0, step = 5.0) |> collect

# build JIT lattice model -
JIT_lattice = build(CRRJITContractPremiumLatticeModel, PutContractModel, S, K, IV_array, DTE; number_of_levels = 80)

# create a point -
test_point = build(PQContractPremiumLatticePoint; i = 10, j = 5, s = 20, d = 10)

# decode -
result_tuple = decode(JIT_lattice,test_point)

# encode a point -
encoded_point = encode(JIT_lattice; S = 42.5, IV = 0.30, K = 105.0, DTE = 63.0)
results_tuple_2 = decode(JIT_lattice, encoded_point)


text:
gw = GridWorld(sx=2, sy=2, rs=[GridWorldState(1,1)], rv=[10.0])
s = initial_state(gw, MersenneTwister(4))

policy = ValuePolicy(gw)

for s in iterator(states(gw))
    @test action(policy, s) isa action_type(gw)
end

@get_requirements ValuePolicy(gw)
@get_requirements action(policy, s)


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

# What is included in this file:
# - the `fmi2ComponentState`--enum which mirrors the internal FMU state (state-machine, not the system state)
# - the `FMU2ComponentEnvironment`- and `FMU2Component`-struct
# - the `FMU2`-struct
# - string/enum-converters for FMI-attribute-structs (e.g. `fmi2StatusToString`, ...)

# this is a custom type to catch the internal mode of the component
@enum fmi2ComponentState begin
    fmi2ComponentStateInstantiated       # after instantiation
    fmi2ComponentStateInitializationMode # after finishing initialization
    fmi2ComponentStateEventMode
    fmi2ComponentStateContinuousTimeMode
    fmi2ComponentStateTerminated
    fmi2ComponentStateError
    fmi2ComponentStateFatal
end

"""
ToDo.
"""
mutable struct FMU2ComponentEnvironment
    logStatusOK::Bool
    logStatusWarning::Bool
    logStatusDiscard::Bool
    logStatusError::Bool
    logStatusFatal::Bool
    logStatusPending::Bool

    function FMU2ComponentEnvironment()
        inst = new()
        inst.logStatusOK = true
        inst.logStatusWarning = true
        inst.logStatusDiscard = true
        inst.logStatusError = true
        inst.logStatusFatal = true
        inst.logStatusPending = true
        return inst
    end
end

"""
The mutable struct represents an instantiated instance of an FMU in the FMI 2.0.2 Standard.
"""
mutable struct FMU2Component
    compAddr::Ptr{Nothing}
    fmu # ::FMU2
    state::fmi2ComponentState
    componentEnvironment::FMU2ComponentEnvironment

    loggingOn::Bool
    callbackFunctions::fmi2CallbackFunctions
    instanceName::String
    continuousStatesChanged::fmi2Boolean
    eventInfo::fmi2EventInfo

    t::fmi2Real             # the system time
    t_offset::fmi2Real      # time offset between simulation environment and FMU
    x::Union{Array{fmi2Real, 1}, Nothing}   # the system states (or sometimes u)
    ẋ::Union{Array{fmi2Real, 1}, Nothing}   # the system state derivative (or sometimes u̇)
    ẍ::Union{Array{fmi2Real, 1}, Nothing}   # the system state second derivative
    #u::Array{fmi2Real, 1}   # the system inputs
    #y::Array{fmi2Real, 1}   # the system outputs
    #p::Array{fmi2Real, 1}   # the system parameters
    z::Array{fmi2Real, 1}   # the system event indicators
    z_prev::Union{Nothing, Array{fmi2Real, 1}}   # the last system event indicators

    realValues::Dict

    x_vrs::Array{fmi2ValueReference, 1}   # the system state value references
    ẋ_vrs::Array{fmi2ValueReference, 1}   # the system state derivative value references
    u_vrs::Array{fmi2ValueReference, 1}   # the system input value references
    y_vrs::Array{fmi2ValueReference, 1}   # the system output value references
    p_vrs::Array{fmi2ValueReference, 1}   # the system parameter value references

    # deprecated
    jac_ẋy_x::Matrix{fmi2Real}
    jac_ẋy_u::Matrix{fmi2Real}
    jac_x::Array{fmi2Real}
    jac_u::Union{Array{fmi2Real}, Nothing}
    jac_t::fmi2Real

    # linearization jacobians
    A::Union{Matrix{fmi2Real}, Nothing}
    B::Union{Matrix{fmi2Real}, Nothing}
    C::Union{Matrix{fmi2Real}, Nothing}
    D::Union{Matrix{fmi2Real}, Nothing}

    jacobianUpdate!             # function for a custom jacobian constructor (optimization)
    skipNextDoStep::Bool    # allows skipping the next `fmi2DoStep` like it is not called
    senseFunc::Symbol       # :auto, :full, :sample, :directionalDerivatives, :adjointDerivatives

    # constructor

    function FMU2Component()
        inst = new()
        inst.state = fmi2ComponentStateInstantiated
        inst.t = -Inf
        inst.t_offset = 0.0
        inst.eventInfo = fmi2EventInfo()

        inst.senseFunc = :auto

        inst.x = nothing
        inst.ẋ = nothing
        inst.ẍ = nothing
        inst.z_prev = nothing

        inst.realValues = Dict()
        inst.x_vrs = Array{fmi2ValueReference, 1}()
        inst.ẋ_vrs = Array{fmi2ValueReference, 1}()
        inst.u_vrs = Array{fmi2ValueReference, 1}()
        inst.y_vrs = Array{fmi2ValueReference, 1}()
        inst.p_vrs = Array{fmi2ValueReference, 1}()

        # initialize further variables
        inst.jac_x = Array{fmi2Real, 1}()
        inst.jac_u = nothing
        inst.jac_t = -1.0
        inst.jac_ẋy_x = zeros(fmi2Real, 0, 0)
        inst.jac_ẋy_u = zeros(fmi2Real, 0, 0)
        inst.skipNextDoStep = false

        inst.A = nothing
        inst.B = nothing
        inst.C = nothing
        inst.D = nothing

        return inst
    end

    function FMU2Component(compAddr, fmu)
        inst = FMU2Component()
        inst.compAddr = compAddr
        inst.fmu = fmu
        return inst
    end
end

"""
Overload the Base.show() function for custom printing of the FMU2Component.
"""
Base.show(io::IO, c::FMU2Component) = print(io,
"FMU:            $(c.fmu.modelDescription.modelName)
InstanceName:   $(isdefined(c, :instanceName) ? c.instanceName : "[not defined]")
Address:        $(c.compAddr)
State:          $(c.state)
Logging:        $(c.loggingOn)
FMU time:       $(c.t)
FMU states:     $(c.x)"
)

"""
A mutable struct representing the excution configuration of a FMU.
For FMUs that have issues with calls like `fmi2Reset` or `fmi2FreeInstance`, this is pretty useful.
"""
mutable struct FMU2ExecutionConfiguration
    terminate::Bool     # call fmi2Terminate before every training step / simulation
    reset::Bool         # call fmi2Reset before every training step / simulation
    setup::Bool         # call setup functions before every training step / simulation
    instantiate::Bool   # call fmi2Instantiate before every training step / simulation
    freeInstance::Bool  # call fmi2FreeInstance after every training step / simulation

    loggingOn::Bool
    externalCallbacks::Bool

    handleStateEvents::Bool                 # handle state events during simulation/training
    handleTimeEvents::Bool                  # handle time events during simulation/training

    assertOnError::Bool                     # wheter an exception is thrown if a fmi2XXX-command fails (>= fmi2StatusError)
    assertOnWarning::Bool                   # wheter an exception is thrown if a fmi2XXX-command warns (>= fmi2StatusWarning)

    autoTimeShift::Bool                     # wheter to shift all time-related functions for simulation intervals not starting at 0.0

    sensealg                                # algorithm for sensitivity estimation over solve call
    useComponentShadow::Bool                # whether FMU outputs/derivatives/jacobians should be cached for frule/rrule (useful for ForwardDiff)
    rootSearchInterpolationPoints::UInt     # number of root search interpolation points
    inPlace::Bool                           # whether faster in-place-fx should be used
    useVectorCallbacks::Bool                # whether to vector (faster) or scalar (slower) callbacks

    maxNewDiscreteStateCalls::UInt          # max calls for fmi2NewDiscreteStates before throwing an exception

    function FMU2ExecutionConfiguration()
        inst = new()

        inst.terminate = true
        inst.reset = true
        inst.setup = true
        inst.instantiate = false
        inst.freeInstance = false

        inst.loggingOn = false
        inst.externalCallbacks = false

        inst.handleStateEvents = true
        inst.handleTimeEvents = true

        inst.assertOnError = false
        inst.assertOnWarning = false

        inst.autoTimeShift = true

        inst.sensealg = nothing # auto
        inst.useComponentShadow = false
        inst.rootSearchInterpolationPoints = 100
        inst.inPlace = true
        inst.useVectorCallbacks = true

        inst.maxNewDiscreteStateCalls = 100

        return inst
    end
end

# default for a "healthy" FMU - this is the fastetst
FMU_EXECUTION_CONFIGURATION_RESET = FMU2ExecutionConfiguration()
FMU_EXECUTION_CONFIGURATION_RESET.terminate = true
FMU_EXECUTION_CONFIGURATION_RESET.reset = true
FMU_EXECUTION_CONFIGURATION_RESET.instantiate = false
FMU_EXECUTION_CONFIGURATION_RESET.freeInstance = false

# if your FMU has a problem with "fmi2Reset" - this is default
FMU_EXECUTION_CONFIGURATION_NO_RESET = FMU2ExecutionConfiguration()
FMU_EXECUTION_CONFIGURATION_NO_RESET.terminate = false
FMU_EXECUTION_CONFIGURATION_NO_RESET.reset = false
FMU_EXECUTION_CONFIGURATION_NO_RESET.instantiate = true
FMU_EXECUTION_CONFIGURATION_NO_RESET.freeInstance = true

# if your FMU has a problem with "fmi2Reset" and "fmi2FreeInstance" - this is for weak FMUs (but slower)
FMU_EXECUTION_CONFIGURATION_NO_FREEING = FMU2ExecutionConfiguration()
FMU_EXECUTION_CONFIGURATION_NO_FREEING.terminate = false
FMU_EXECUTION_CONFIGURATION_NO_FREEING.reset = false
FMU_EXECUTION_CONFIGURATION_NO_FREEING.instantiate = true
FMU_EXECUTION_CONFIGURATION_NO_FREEING.freeInstance = false

"""
ToDo
"""
struct FMU2Event
    t::Union{Float32, Float64}
    indicator::UInt

    x_left::Union{Array{Float64, 1}, Array{Float32, 1}, Nothing}
    x_right::Union{Array{Float64, 1}, Array{Float32, 1}, Nothing}

    function FMU2Event(t::Union{Float32, Float64}, indicator::UInt = 0, x_left::Union{Array{Float64, 1}, Array{Float32, 1}, Nothing} = nothing, x_right::Union{Array{Float64, 1}, Array{Float32, 1}, Nothing} = nothing)
        inst = new(t, indicator, x_left, x_right)
        return inst
    end
end

"""
Overload the Base.show() function for custom printing of the FMU2.
"""
Base.show(io::IO, e::FMU2Event) = print(io, e.indicator == 0 ? "Time-Event @ $(round(e.t; digits=4))s" : "State-Event #$(e.indicator) @ $(round(e.t; digits=4))s")

"""
ToDo
"""
mutable struct FMU2Solution <: FMUSolution
    fmu                                             # FMU2
    success::Bool

    states                                          # ODESolution

    values
    valueReferences::Union{Array, Nothing}          # Array{fmi2ValueReference}

    events::Array{FMU2Event, 1}

    function FMU2Solution(fmu)
        inst = new()

        inst.fmu = fmu
        inst.success = false
        inst.states = nothing
        inst.values = nothing
        inst.valueReferences = nothing

        inst.events = []

        return inst
    end
end

"""
Overload the Base.show() function for custom printing of the FMU2.
"""
function Base.show(io::IO, sol::FMU2Solution)
    print(io, "Model name:\n\t$(sol.fmu.modelDescription.modelName)\nSuccess:\n\t$(sol.success)\n")

    if sol.states != nothing
        print(io, "States [$(length(sol.states))]:\n")
        if length(sol.states.u) > 10
            for i in 1:9
                print(io, "\t$(sol.states.t[i])\t$(sol.states.u[i])\n")
            end
            print(io, "\t...\n\t$(sol.states.t[end])\t$(sol.states.u[end])\n")
        else
            for i in 1:length(sol.states)
                print(io, "\t$(sol.states.t[i])\t$(sol.states.u[i])\n")
            end
        end
    end

    if sol.values != nothing
        print(io, "Values [$(length(sol.values.saveval))]:\n")
        if length(sol.values.saveval) > 10
            for i in 1:9
                print(io, "\t$(sol.values.t[i])\t$(sol.values.saveval[i])\n")
            end
            print(io, "\t...\n\t$(sol.values.t[end])\t$(sol.values.saveval[end])\n")
        else
            for i in 1:length(sol.values.saveval)
                print(io, "\t$(sol.values.t[i])\t$(sol.values.saveval[i])\n")
            end
        end
    end

    if sol.events != nothing
        print(io, "Events [$(length(sol.events))]:\n")
        if length(sol.events) > 10
            for i in 1:9
                print(io, "\t$(sol.events[i])\n")
            end
            print(io, "\t...\n\t$(sol.events[end])\n")
        else
            for i in 1:length(sol.events)
                print(io, "\t$(sol.events[i])\n")
            end
        end
    end

end

"""
The mutable struct representing a FMU and all it instantiated instances in the FMI 2.0.2 Standard.
Also contains the paths to the FMU and ZIP folder as well als all the FMI 2.0.2 function pointers.
"""
mutable struct FMU2 <: FMU
    modelName::String
    fmuResourceLocation::String

    modelDescription::fmi2ModelDescription

    type::fmi2Type
    components::Array # {fmi2Component}

    # c-functions
    cInstantiate::Ptr{Cvoid}
    cGetTypesPlatform::Ptr{Cvoid}
    cGetVersion::Ptr{Cvoid}
    cFreeInstance::Ptr{Cvoid}
    cSetDebugLogging::Ptr{Cvoid}
    cSetupExperiment::Ptr{Cvoid}
    cEnterInitializationMode::Ptr{Cvoid}
    cExitInitializationMode::Ptr{Cvoid}
    cTerminate::Ptr{Cvoid}
    cReset::Ptr{Cvoid}
    cGetReal::Ptr{Cvoid}
    cSetReal::Ptr{Cvoid}
    cGetInteger::Ptr{Cvoid}
    cSetInteger::Ptr{Cvoid}
    cGetBoolean::Ptr{Cvoid}
    cSetBoolean::Ptr{Cvoid}
    cGetString::Ptr{Cvoid}
    cSetString::Ptr{Cvoid}
    cGetFMUstate::Ptr{Cvoid}
    cSetFMUstate::Ptr{Cvoid}
    cFreeFMUstate::Ptr{Cvoid}
    cSerializedFMUstateSize::Ptr{Cvoid}
    cSerializeFMUstate::Ptr{Cvoid}
    cDeSerializeFMUstate::Ptr{Cvoid}
    cGetDirectionalDerivative::Ptr{Cvoid}

    # Co Simulation function calls
    cSetRealInputDerivatives::Ptr{Cvoid}
    cGetRealOutputDerivatives::Ptr{Cvoid}
    cDoStep::Ptr{Cvoid}
    cCancelStep::Ptr{Cvoid}
    cGetStatus::Ptr{Cvoid}
    cGetRealStatus::Ptr{Cvoid}
    cGetIntegerStatus::Ptr{Cvoid}
    cGetBooleanStatus::Ptr{Cvoid}
    cGetStringStatus::Ptr{Cvoid}

    # Model Exchange function calls
    cEnterContinuousTimeMode::Ptr{Cvoid}
    cGetContinuousStates::Ptr{Cvoid}
    cGetDerivatives::Ptr{Cvoid}
    cSetTime::Ptr{Cvoid}
    cSetContinuousStates::Ptr{Cvoid}
    cCompletedIntegratorStep::Ptr{Cvoid}
    cEnterEventMode::Ptr{Cvoid}
    cNewDiscreteStates::Ptr{Cvoid}
    cGetEventIndicators::Ptr{Cvoid}
    cGetNominalsOfContinuousStates::Ptr{Cvoid}

    # paths of ziped and unziped FMU folders
    path::String
    binaryPath::String
    zipPath::String

    # execution configuration
    executionConfig::FMU2ExecutionConfiguration
    hasStateEvents::Union{Bool, Nothing}
    hasTimeEvents::Union{Bool, Nothing}

    # c-libraries
    libHandle::Ptr{Nothing}
    callbackLibHandle::Ptr{Nothing}

    # START: experimental section (to FMIFlux.jl)
    dependencies::Matrix{Union{fmi2DependencyKind, Nothing}}

    # END: experimental section

    # DEPRECATED
    instanceName::String

    # Constructor
    function FMU2()
        inst = new()
        inst.components = []
        inst.callbackLibHandle = C_NULL

        inst.hasStateEvents = nothing
        inst.hasTimeEvents = nothing

        inst.executionConfig = FMU_EXECUTION_CONFIGURATION_NO_RESET

        return inst
    end
end

"""
Overload the Base.show() function for custom printing of the FMU2.
"""
Base.show(io::IO, fmu::FMU2) = print(io,
"Model name:        $(fmu.modelDescription.modelName)
Type:              $(fmu.type)"
)

"""
Formats a fmi2Status to String.
"""
function fmi2StatusToString(status::fmi2Status)
    if status == fmi2StatusOK
        return "OK"
    elseif status == fmi2StatusWarning
        return "Warning"
    elseif status == fmi2StatusDiscard
        return "Discard"
    elseif status == fmi2StatusError
        return "Error"
    elseif status == fmi2StatusFatal
        return "Fatal"
    elseif status == fmi2StatusPending
        return "Pending"
    else
        return "Unknown"
    end
end

"""
Formats a fmi2Status to Integer.
"""
function fmi2StatusToString(status::Integer)
    if status == fmi2StatusOK
        return "OK"
    elseif status == fmi2StatusWarning
        return "Warning"
    elseif status == fmi2StatusDiscard
        return "Discard"
    elseif status == fmi2StatusError
        return "Error"
    elseif status == fmi2StatusFatal
        return "Fatal"
    elseif status == fmi2StatusPending
        return "Pending"
    else
        return "Unknown"
    end
end

function fmi2CausalityToString(c::fmi2Causality)
    if c == fmi2CausalityParameter
        return "parameter"
    elseif c == fmi2CausalityCalculatedParameter
        return "calculatedParameter"
    elseif c == fmi2CausalityInput
        return "input"
    elseif c == fmi2CausalityOutput
        return "output"
    elseif c == fmi2CausalityLocal
        return "local"
    elseif c == fmi2CausalityIndependent
        return "independent"
    else
        @assert false "fmi2CausalityToString(...): Unknown causality."
    end
end

function fmi2StringToCausality(s::AbstractString)
    if s == "parameter"
        return fmi2CausalityParameter
    elseif s == "calculatedParameter"
        return fmi2CausalityCalculatedParameter
    elseif s == "input"
        return fmi2CausalityInput
    elseif s == "output"
        return fmi2CausalityOutput
    elseif s == "local"
        return fmi2CausalityLocal
    elseif s == "independent"
        return fmi2CausalityIndependent
    else
        @assert false "fmi2StringToCausality($(s)): Unknown causality."
    end
end

function fmi2VariabilityToString(c::fmi2Variability)
    if c == fmi2VariabilityConstant
        return "constant"
    elseif c == fmi2VariabilityFixed
        return "fixed"
    elseif c == fmi2VariabilityTunable
        return "tunable"
    elseif c == fmi2VariabilityDiscrete
        return "discrete"
    elseif c == fmi2VariabilityContinuous
        return "continuous"
    else
        @assert false "fmi2VariabilityToString(...): Unknown variability."
    end
end

function fmi2StringToVariability(s::AbstractString)
    if s == "constant"
        return fmi2VariabilityConstant
    elseif s == "fixed"
        return fmi2VariabilityFixed
    elseif s == "tunable"
        return fmi2VariabilityTunable
    elseif s == "discrete"
        return fmi2VariabilityDiscrete
    elseif s == "continuous"
        return fmi2VariabilityContinuous
    else
        @assert false "fmi2StringToVariability($(s)): Unknown variability."
    end
end

function fmi2InitialToString(c::fmi2Initial)
    if c == fmi2InitialApprox
        return "approx"
    elseif c == fmi2InitialExact
        return "exact"
    elseif c == fmi2InitialCalculated
        return "calculated"
    else
        @assert false "fmi2InitialToString(...): Unknown initial."
    end
end

function fmi2StringToInitial(s::AbstractString)
    if s == "approx"
        return fmi2InitialApprox
    elseif s == "exact"
        return fmi2InitialExact
    elseif s == "calculated"
        return fmi2InitialCalculated
    else
        @assert false "fmi2StringToInitial($(s)): Unknown initial."
    end
end

function fmi2DependencyKindToString(c::fmi2DependencyKind)
    if c == fmi2DependencyKindDependent
        return "dependent"
    elseif c == fmi2DependencyKindConstant
        return "constant"
    elseif c == fmi2DependencyKindFixed
        return "fixed"
    elseif c == fmi2DependencyKindTunable
        return "tunable"
    elseif c == fmi2DependencyKindDiscrete
        return "discrete"
    else
        @assert false "fmi2DependencyKindToString(...): Unknown dependency kind."
    end
end

function fmi2StringToDependencyKind(s::AbstractString)
    if s == "dependent"
        return fmi2DependencyKindDependent
    elseif s == "exact"
        return fmi2DependencyKindConstant
    elseif s == "fixed"
        return fmi2DependencyKindFixed
    elseif s == "tunable"
        return fmi2DependencyKindTunable
    elseif s == "discrete"
        return fmi2DependencyKindDiscrete
    else
        @assert false "fmi2StringToDependencyKind($(s)): Unknown dependency kind."
    end
end


text:
# Acrobot
T = Float64
model = TrajectoryOptimization.Dynamics.acrobot_model
model_d = rk3(model)
n = model.n; m = model.m

# costs
Qf = 100.0*Diagonal(I,n)
Q = (1.0e-2)*Diagonal(I,n)
R = (1e-2)*Diagonal(I,m)

x0 = [0.; 0.; 0.; 0.]
xf = [π; 0.; 0.; 0.];
N = 151
dt = 0.01
U0 = [ones(m) for k = 1:N-1]
obj = TrajectoryOptimization.LQRObjective(Q,R,Qf,xf,N)
goal = goal_constraint(xf)

acrobot = TrajectoryOptimization.Problem(model_d, obj, x0=x0, xf=xf, N=N, dt=dt)
acrobot.constraints[N] += goal
initial_controls!(acrobot, U0)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
module GitLab

@info """
    Loading src/gitlab.jl via Requires.jl. JSON2 wasn't added to Project.toml to avoid
    compat issues in projects where GitLab isn't used.
    """

import JSON2

using Codex
using HTTP
using HTTP: Response

export enforce_schedules, enforce_variables

struct Config
    token::AbstractString
    url::AbstractString
end

struct Project
    config::Config
    # Not using urlencoded NAMESPACE/PROJECT_NAME for simplicity.
    project_id::Int
end
project_url(p::Project) = "$(p.config.url)/api/v4/projects/$(p.project_id)"

struct Schedule
    project::Project
    schedule_id::Int
end
schedule_url(s::Schedule) =
    "$(project_url(s.project))/pipeline_schedules/$(s.schedule_id)"

auth_header(config::Config) = Dict("PRIVATE-TOKEN" => config.token)
auth_header(p::Project) = auth_header(p.config)
auth_header(s::Schedule) = auth_header(s.project)

# All values are coverted to string to avoid errors when using `active = false`.
form(param::NamedTuple) = HTTP.Form(nt2dict(Codex.apply(string, param)))

resp2json(r::Response) = JSON2.read(String(r.body))

"""
    list_schedules(p::Project) -> Vector{NamedTuple}

Returns an array of named tuples `(id = ..., description = ..., ...)`.
For details, see
<https://docs.gitlab.com/ee/api/pipeline_schedules.html#get-all-pipeline-schedules>.
"""
function list_schedules(p::Project)::Vector{NamedTuple}
    endpoint = "$(project_url(p))/pipeline_schedules"
    r = HTTP.request("GET", endpoint, auth_header(p))
    list = resp2json(r)
    return list
end
n_schedules(p::Project) = length(GitLab.list_schedules(p))

"""
    create_schedule(p::Project, param::NamedTuple) -> NamedTuple

For `param`, see
<https://docs.gitlab.com/ee/api/pipeline_schedules.html#create-a-new-pipeline-schedule>.
"""
function create_schedule(p::Project, param::NamedTuple)::NamedTuple
    endpoint = "$(project_url(p))/pipeline_schedules"
    r = HTTP.request("POST", endpoint, auth_header(p), form(param))
    nt = resp2json(r)
    for key in keys(param)
        @assert param[key] == nt[key] "$key: $(param[key]) != $(nt[key])"
    end
    return nt
end

"""
    delete_schedule(s::Schedule) -> NamedTuple

Deletes `schedule_id` for `project_id`.
"""
function delete_schedule(s::Schedule)::NamedTuple
    r = HTTP.request("DELETE", schedule_url(s), auth_header(s))
    return resp2json(r)
end

"""
    create_schedule_variable(v::Variable, param::NamedTuple) -> NamedTuple

For param, see
<https://docs.gitlab.com/ee/api/pipeline_schedules.html#create-a-new-pipeline-schedule-variable>.
This codebase does not define edit or delete methods because the variables
are deleted via `enforce_schedules`.
"""
function create_schedule_variable(s::Schedule, param::NamedTuple)::NamedTuple
    endpoint = "$(schedule_url(s))/variables"
    r = HTTP.request("POST", endpoint, auth_header(s), form(param))
    nt = resp2json(r)
    @assert param[:value] == nt[:value]
    nt
end

"""
    create_schedule_variable(s::Schedule, params::Vector) -> Vector

Set multiple variables via `params`.
For example,
```
params = [
    (key = "key1", value = "value1"),
    (key = "key2", value = "value2")
]
```
"""
function create_schedule_variable(s::Schedule, params::Vector)::Vector
    [create_schedule_variable(s::Schedule, param) for param in params]
end

"""
    enforce_schedules(p::Project, params::Vector; variables=[]) -> Vector{NamedTuple}

Enforce schedules for `p` as defined by `params`.
Variables per schedule can be set via `variables`.
For example,
```
params = [
    (description = "A", ...),
    (description = "B", ...)
]
variables = [
    [(key = "key1", value = "value1")], # for schedule A.
    [(key = "key1", value = "value2")] # for schedule B.
]
```
"""
function enforce_schedules(p::Project, params::Vector; variables=[])::Vector{NamedTuple}
    [delete_schedule(Schedule(p, nt.id)) for nt in list_schedules(p)]
    @assert n_schedules(p) == 0
    list = [create_schedule(p, param) for param in params]
    if 0 < length(variables)
        [create_schedule_variable(Schedule(p, t[2].id), variables[t[1]]) for t in enumerate(list)]
    end
    list
end

function list_variables(p::Project)::Vector
    endpoint = "$(project_url(p))/variables"
    r = HTTP.request("GET", endpoint, auth_header(p))
    resp2json(r)
end

function delete_variable(p::Project, key::String)::NamedTuple
    endpoint = "$(project_url(p))/variables/$key"
    r = HTTP.request("DELETE", endpoint, auth_header(p))
    resp2json(r)
end

function create_variable(p::Project, param::NamedTuple)::NamedTuple
    endpoint = "$(project_url(p))/variables"
    r = HTTP.request("POST", endpoint, auth_header(p), form(param))
    resp2json(r)
end

"""
    enforce_variables(p::Project, variables::Vector) -> Vector

Enforces variables for project `p`.
For example,
```
variables = [
    (key = "key1", value = "value1", protected = true),
    (key = "key2", value = "value2", protected = false)
]
```
This method can fail with `key` exists.
Possibly because the create request is too quick.
"""
function enforce_variables(p::Project, variables::Vector)::Vector
    list = list_variables(p)
    [delete_variable(p, x.key) for x in list]
    @assert length(list_variables(p)) == 0
    [create_variable(p, param) for param in variables]
end

end # module


text:
@testset "Test the 2-by-2 symmetric positive definite matrices" begin
using Random: seed!
# explicitly compute an easy exp
M = SymmetricPositiveDefinite(2)
x = SPDPoint([1.0 0.0; 0.0 1.0])
ξ = SPDTVector([1.0 0.0; 0.0 0.0])
η = SPDTVector([1.0 0.0; 0.0 1.0])
y = exp(M,x,ξ)
@test norm( getValue(y) - [ℯ 0;0 1.0]) ≈ 0 atol=10.0^(-16)
# check that with base the base must mach.
z = SPDPoint([1.0 0.0; 0.0 0.1])
ξE = TVectorE(ξ,x);
@test_throws DomainError exp(M,z,ξE)
# check that log is the inverse of exp and stores the base point correctly
# if that's activated
xT = MPointE(x);
ξT2 = log(M,x,y)
@test norm( getValue(ξ) - getValue(ξT2) ) ≈ 0 atol=10.0^(-16)
@test distance(M, getBasePoint(ξE), x) ≈ 0 atol=10.0^(-16)
# test norm
@test norm(M,x,ξ) ≈ distance(M,x,y) atol=10.0^(-16)
# test parallel transport
@test norm(getValue( parallelTransport(M,y,x,log(M,y,x)) ) + getValue(ξ) ) ≈ 0 atol=3*10.0^(-16)
@test parallelTransport(M,x,x,ξ) == ξ
# random
seed!(42)
@test validateMPoint(M,randomMPoint(M))
@test validateTVector(M,x, randomTVector(M,x))
@test validateTVector(M,x, randomTVector(M,x, :Gaussian))
@test validateTVector(M,x, randomTVector(M,x, :Rician))
#
# Test tangent ONB
#
n = 4
M2 = SymmetricPositiveDefinite(n)
x2 = SPDPoint(one(zeros(n,n)))
ξ2 = SPDTVector(getValue(x2))
y2 = exp(M2,x2,ξ2)
Ξ,κ = tangentONB(M2,x2,ξ2)
# test orthogonality
@test all([ dot(M2,x2,Ξ[i],Ξ[j]) for i=1:length(Ξ) for j=i+1:length(Ξ) ] .== 0)
# test normality
@test all( [1-dot(M2,x2,Ξ[i],Ξ[i]) for i=1:length(Ξ)] .< 3*10^(-16) )
Ξ2,κ2 = tangentONB(M2,x2,y2)
@test sum( norm.(Ref(M2),Ref(x2), Ξ .- Ξ2) ) ≈ 0
@test sum(κ - κ2) ≈ 0
#
# test types
@test typeofTVector(x) == SPDTVector{Float64}
@test typeofTVector(SPDPoint{Float64}) == SPDTVector{Float64}
@test typeofMPoint(ξ) == SPDPoint{Float64}
@test typeofMPoint(SPDTVector{Float64}) == SPDPoint{Float64}
#
@test typicalDistance(M) ≈ sqrt(3)
#
@test_throws DomainError validateMPoint(M2,x)
@test_throws DomainError validateMPoint(M, SPDPoint([0. 1.;-1. 0.]))
@test_throws DomainError validateMPoint(M, SPDPoint([0. 0.;0. 0.]))
@test_throws DomainError validateTVector(M,x,ξ2)
@test_throws DomainError validateTVector(M,x,SPDTVector([0. 1.; -2. 0.]) )
#
@test norm(getValue(zeroTVector(M,x))) == 0
@test validateTVector(M,x,zeroTVector(M,x))
# show
@test "$M" == "The Manifold of 2-by-2 symmetric positive definite matrices"
@test "$x" == "SPD($(getValue(x)))"
@test "$ξ" == "SPDT($(getValue(ξ)))"
# Test Matrix trait
@test getValue(x+y) == getValue(x)+getValue(y)
@test getValue(x-y) == getValue(x)-getValue(y)
@test getValue(x*y) == getValue(x)*getValue(y)
@test getValue(ξ+η) == getValue(ξ)+getValue(η)
@test getValue(ξ-η) == getValue(ξ)-getValue(η)
@test getValue(ξ*η) == getValue(ξ)*getValue(η)
@test transpose(x) == SPDPoint(Matrix(transpose(getValue(x))))
@test transpose(ξ) == SPDTVector(Matrix(transpose(getValue(ξ))))
end


text:
@assert length(ARGS) == 2
file_loc = joinpath(dirname(@__FILE__), "innerjoin_performance.jl")
llen = ARGS[1]
rlen = ARGS[2]

for a3 in ["str", "int", "pool", "cat"],
    a4 in ["uniq", "dup", "manydup"],
    a5 in ["sort", "rand"],
    a6 in ["1", "2"]
    a4 == "uniq" && a3 in ["pool", "cat"] && continue
    run(`julia $file_loc $llen $rlen $a3 $a4 $a5 $a6`)
end


text:
module Jordan
using Hecke

function rref_with_trans(M::T) where {T <: MatElem{S} where {S <: FieldElem}}
  #does row ops
  n = hcat(M, identity_matrix(base_ring(M), nrows(M)))
  rref!(n)
  s = nrows(n)
  while s>0 && iszero(sub(n, s:s, 1:ncols(M)))
    s -= 1
  end
  return s, sub(n, 1:nrows(M), 1:ncols(M)), sub(n, 1:nrows(M), ncols(M)+1:ncols(n))
end

function quo_gens(A, B)
  n, B = rref(B')
  B = sub(B, 1:n, 1:ncols(B))'
  fl, T = can_solve_with_solution(A, B)  # A T = B
  d, r = rref(T')
  t = sub(r, 1:d, 1:ncols(r))
  t = Hecke.reduce_mod(identity_matrix(base_ring(t), ncols(t)), t)
  d = rref!(t)
  t = t[1:d, :]
  return A*t'
end

function rational_normal_form(M::T) where {T <: MatElem{S} where {S <: FieldElem}}
  f = minpoly(M)
  lf = factor(f)
  all_b = zero_matrix(base_ring(M), nrows(M), 0)
  for (g, v) = lf.fac
    m = (g^v)(M)
    this_g = zero_matrix(base_ring(M), nrows(M), 0)
    while true
      d, km = nullspace(m)
      km = quo_gens(km, this_g)
      if ncols(km) == 0
        break
      end
      ex = v
      while true
        b = quo_gens(km, g(M)*km)
        if ncols(b) > 0
          break
        end
        @show ex -= 1
        km = g(M)*km
        d, km = rref(km')
        km = km[1:d, :]'
      end
      b = b[:, 1]
    #any 0 ne b in km/g(M)km will generate a sigma cyclic subspace
    # b, Mb M^2b ... g(M)b, M g(M) b, ... g^2(M)b, ....
      ve = deepcopy(b)
      b = zero_matrix(base_ring(M), nrows(M), 0)
      for k = 1:ex
        for j = 1:degree(g)
          @show k, j
          b = hcat(b, (M^(j-1))*((g^(k-1))(M))*ve)
        end
      end
      @show this_g = hcat(this_g, b)
    end
    all_b = hcat(all_b, this_g)
  end
  return all_b
end

end


text:
export BoundedLattice, boundary, isperiodic

"""
    BoundedLattice{N, BC} <: AbstractLattice{N}

A lattice with boundary conditions `BC` on `N` dimension.
"""
abstract type BoundedLattice{N, BC <: BoundaryCondition} <: AbstractLattice{N} end

# Bounded Lattice Traits
"""
    boundary(lattice) -> BoundaryCondition

Returns lattice's boundary condition.
"""
function boundary end

"""
    isperiodic(lattice) -> Bool

Check if this lattice's boundary condition is periodic.
"""
function isperiodic end

boundary(::Type{<:BoundedLattice{N, BC}}) where {N, BC} = BC

isperiodic(::Type{<:BoundedLattice{N, Periodic}}) where N = true
isperiodic(::Type{<:BoundedLattice}) = false

# forward instance traits
isperiodic(::T) where {T <: BoundedLattice} = isperiodic(T)
boundary(::T) where {T <: BoundedLattice} = boundary(T)

nameof(::BoundedLattice) = "Bounded Lattice"

function show(io::IO, lattice::BoundedLattice)
    indent = get(io, :indent, 0)
    println(io, " " ^ indent, nameof(lattice))
    println(io, " " ^ indent, "boundary: ", boundary(lattice))
    print(io, " " ^ indent, "size:     ", size(lattice))
end


text:
using OmegaCore
using OmegaRJMCMC
using Distributions

"From Green et. al"
function poissonmodel()
  # Bayesian model over two models
  # 1. Poisson omodel
  # 2. Negative binomeial model

  function f(ω)
    k = 1 ~ DiscreteUniform(1, 2) # model indicator, selects model
    N = 10
    αλ, βλ = 25, 10
    λ = 2 ~ Gamma(αλ, βλ)
    λ_ = λ(ω)
    if k(ω) == 1
      Y = 5 ~ Poisson(λ_)
    else
      ακ, βκ = 1, 10
      κ = (3 ~ Gamma(ακ, βκ))
      κ_ = κ(ω)
      # Distributions has different parameterization
      # https://wiki.analytica.com/index.php?title=Negative_binomial_distribution#Alternate_parameterizations
      Y = 5 ~ NegativeBinomial(κ_, λ_ / (λ_ + κ_))
    end
    Y(ω)
    # Mv(Y, N)(ω)
  end

  # TODO: Mv needs to work
  # TODO: get data for Y
  function h(ω::T) where
    μ = 0.5
    x = 0.2
    σ = 0.1 # Parameter
    u = Normal(0, σ)
    if ω[κ] == 1
      (k => 2, λ => ω[λ], κ => μ * exp(u))
    else
      (k => 1, )
    end
  end

  randsample((k, λ), Y == Ydata, Constant(h),j, alg = RJMCMC)

end

## FIXME / TODO
# Should ω be RandVar -> τ
# or RandVarId -> τ
#

# randsample(poissonmodel())

#   # This is literal pointwise semanics
#   q(ω) -> Poisson(λ(ω))

#   So that makes sense, but
#   we want kind of parity with parameters
#   Poission(0.4) produces a ciid class
#   lift(Poisson)(λ) produces a distribution over ciid classes

#   q(id, ω) = Poisson(λ(ω))(id, ω)

#   into

#   (id, ω) -> Poisson(_)(id, ω)

#   The problem with this reduction is that we would need to ensure hte \oemgas are consistent

#   What's to stop me form violating that.
#   In this snse i cant see it as a valid independent procedure.

#   u might say the general mechanism is just partial application.

#   I'm not sure what it should produce.
#   On the one hand the lifted notation makes sense because its the literal
#   extention of lifting.

#   On the other hand, the version that produces a `ciid` makes sense too.
#   Practically, if we have Mv(i)

#   Option 1: Partial application-
#   Option 2: specialise to this classes
#   Option 3: just use special types

#   @inline Space.recurse(mv::Mv{<:Distribution}, id, ω) =
#     recurse(Mv(mv.dist(ω)), id, ω)

#   # Model 1

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
module Util
  hasarg(s::Symbol, kwargs) = length(filter((x -> x[1] == s), kwargs)) > 0
  function getarg(s::Symbol, kwargs, defaults)
    for (k, v) in kwargs
      if k == s; return v end
    end
    if !haskey(defaults, s); error("Unknown keyword argument $s"); end
    deepcopy(defaults[s])
  end

  function heretempfile(prefix::String="")
    newname() = abspath(joinpath(prefix, basename(tempname())))
    filename = newname()
    while ispath(filename)
      filename = newname()
    end
    filename
  end
end
using .Util


text:
using DDCModelsExample
using Test

@testset "DDCbasic.jl" begin

    # Define number of periods and number of firms
    nperiod = 100
    nfirm = 10000

    # Define support of X
    supportX = collect(1:5)
    nsupportX = length(supportX)

    # Define Π
    Π = 1 ./ (1 .+ abs.(repeat(supportX',nsupportX,1) .- repeat(supportX,1,nsupportX)))
    Π ./= sum(Π;dims=2)

    # Define parameters
    β = [-0.1*nsupportX;0.2]
    δ = [0.0;1.0]
    ρ = 0.95

    # Construct model
    m = DDCbasic(;supportX,Π,β,δ,ρ)

    # Simulate dataset
    X,choices = gen_data(m,nperiod,nfirm)

    # Estimate Π
    Πestim = estimateΠ(X,m.supportX)

    # Test whether estimated Π is approximately the same as
    # true Π with which data was simulated
    @test maximum(abs.(Πestim .- m.Π)) < 0.01

    # Estimate model
    mestim = estimate_model(m,X,choices; startvals = [-1;-0.1;0.5])
    # Calculate standard errors
    stderr = calc_stderr(mestim,X,choices)

    # Test whether coefficient estimates within 2 std
    # errors of true value
    @test maximum(abs.(m.β .- mestim.β) .< 2 .* stderr[1:2])
    @test abs(m.δ[2] - m.δ[2])  < 2*stderr[3]

end

max_stars_repo_name: rgreminger/DDCModelsExample

max_issues_repo_name: rgreminger/DDCModelsExample

max_forks_repo_name: rgreminger/DDCModelsExample


text:
# module Border

using OffsetArrays, CatIndices

abstract AbstractBorder

immutable NoPad{T} <: AbstractBorder
    border::T
end
NoPad() = NoPad(nothing)

"""
    NoPad()
    NoPad(border)

Indicates that no padding should be applied to the input array, or that you have already pre-padded the input image. Passing a `border` object allows you to preserve "memory" of a border choice; it can be retrieved by indexing with `[]`.

# Example

    np = NoPad(Pad(:replicate))
    imfilter!(out, img, kernel, np)

runs filtering directly, skipping any padding steps.  Every entry of
`out` must be computable using in-bounds operations on `img` and
`kernel`.
"""
NoPad

Base.getindex(np::NoPad) = np.border

"""
`Pad` is a type that stores choices about padding. Instances must set `style`, a
Symbol specifying the boundary conditions of the image, one of:

- `:replicate` (repeat edge values to infinity)
- `:circular` (image edges "wrap around")
- `:symmetric` (the image reflects relative to a position between pixels)
- `:reflect` (the image reflects relative to the edge itself)

The default value is `:replicate`.

It's worth emphasizing that padding is most straightforwardly specified as a string,

    imfilter(img, kernel, "replicate")

rather than

    imfilter(img, kernel, Pad(:replicate))
"""
immutable Pad{N} <: AbstractBorder
    style::Symbol
    lo::Dims{N}    # number to extend by on the lower edge for each dimension
    hi::Dims{N}    # number to extend by on the upper edge for each dimension
end

const valid_borders = ("replicate", "circular", "reflect", "symmetric")

"""
    Pad(style::Symbol, m, n, ...)

Pad the input image symmetrically, `m` pixels at the lower and upper edge of dimension 1, `n` pixels for dimension 2, and so forth.
"""
Pad(style::Symbol, both::Int...) = Pad(style, both, both)
Pad(both::Int...) = Pad(:replicate, both, both)

"""
    Pad(style::Symbol, (m,n))

Pad the input image symmetrically, `m` pixels at the lower and upper edge of dimension 1, `n` pixels for dimension 2.
"""
Pad(style::Symbol, both::Dims) = Pad(style, both, both)
Pad(both::Dims) = Pad(:replicate, both, both)

"""
    Pad(style::Symbol, lo::Dims, hi::Dims)

Pad the input image by `lo` pixels at the lower edge, and `hi` pixels at the upper edge.
"""
Pad(lo::Dims, hi::Dims) = Pad(:replicate, lo, hi)
Pad(style::Symbol, lo::Tuple{}, hi::Tuple{}) = Pad{0}(style, lo, hi)
Pad{N}(style::Symbol, lo::Dims{N}, hi::Tuple{}) = Pad(style, lo, ntuple(d->0,Val{N}))
Pad{N}(style::Symbol, lo::Tuple{}, hi::Dims{N}) = Pad(style, ntuple(d->0,Val{N}), hi)
Pad(style::Symbol, lo::AbstractVector{Int}, hi::AbstractVector{Int}) = Pad(style, (lo...,), (hi...,))

Pad(style::Symbol, inds::Indices) = Pad(style, map(lo,inds), map(hi,inds))

"""
    Pad(style, kernel)
    Pad(style)(kernel)

Given a filter array `kernel`, determine the amount of padding from the `indices` of `kernel`.
"""
(p::Pad{0})(kernel) = Pad(p.style, calculate_padding(kernel))
(p::Pad{0})(kernel, img, ::Alg) = p(kernel)

# Padding for FFT: round up to next size expressible as 2^m*3^n
function (p::Pad{0})(kernel, img, ::FFT)
    inds = calculate_padding(kernel)
    newinds = map(padfft, inds, map(length, indices(img)))
    Pad(p.style, newinds)
end
function padfft(indk::AbstractUnitRange, l::Integer)
    lk = length(indk)
    range(first(indk), nextprod([2,3], l+lk)-l+1)
end

function padindices{_,N}(img::AbstractArray{_,N}, border::Pad)
    throw(ArgumentError("$border lacks the proper padding sizes for an array with $(ndims(img)) dimensions"))
end
function padindices{_,N}(img::AbstractArray{_,N}, border::Pad{N})
    _padindices(border, border.lo, indices(img), border.hi)
end
function padindices{P<:Pad}(img::AbstractArray, ::Type{P})
    throw(ArgumentError("must supply padding sizes to $P"))
end

# The 3-argument map is not inferrable, so do it manually
@inline _padindices(border, lo, inds, hi) =
    (padindex(border, lo[1], inds[1], hi[1]),
     _padindices(border, tail(lo), tail(inds), tail(hi))...)
_padindices(border, ::Tuple{}, ::Tuple{}, ::Tuple{}) = ()

"""
    padarray([T], img, border) --> imgpadded

Generate a padded image from an array `img` and a specification
`border` of the boundary conditions and amount of padding to
add. `border` can be a `Pad`, `Fill`, or `Inner` object.

Optionally provide the element type `T` of `imgpadded`.
"""
padarray(img::AbstractArray, border::Pad)  = padarray(eltype(img), img, border)
function padarray{T}(::Type{T}, img::AbstractArray, border::Pad)
    inds = padindices(img, border)
    # like img[inds...] except that we can control the element type
    newinds = map(Base.indices1, inds)
    dest = similar(img, T, newinds)
    @unsafe for I in CartesianRange(newinds)
        J = CartesianIndex(map((i,x)->x[i], I.I, inds))
        dest[I] = img[J]
    end
    dest
end
padarray{P}(img, ::Type{P}) = img[padindices(img, P)...]      # just to throw the nice error

Base.ndims{N}(::Pad{N}) = N

# Make these separate types because the dispatch almost surely needs to be different
"""
    Inner()
    Inner(lo, hi)

Indicate that edges are to be discarded in filtering, only the interior of the result it to be returned.

# Example:

    imfilter(img, kernel, Inner())
"""
immutable Inner{N} <: AbstractBorder
    lo::Dims{N}
    hi::Dims{N}
end

"""
    NA()
    NA(lo, hi)

Choose filtering using "NA" (Not Available) boundary conditions. This
is most appropriate for filters that have only positive weights, such
as blurring filters. Effectively, the output pixel value is normalized
in the following way:

              filtered img with Fill(0) boundary conditions
    output =  ---------------------------------------------
              filtered 1   with Fill(0) boundary conditions

As a consequence, filtering has the same behavior as
`nanmean`. Indeed, invalid pixels in `img` can be marked as `NaN` and
then they are effectively omitted from the filtered result.
"""
immutable NA{N} <: AbstractBorder
    lo::Dims{N}
    hi::Dims{N}
end

for T in (:Inner, :NA)
    @eval begin
        (::Type{$T})(both::Int...) = $T(both, both)
        (::Type{$T}){N}(both::Dims{N}) = $T(both, both)
        (::Type{$T})(lo::Tuple{}, hi::Tuple{}) = $T{0}(lo, hi)
        (::Type{$T}){N}(lo::Dims{N}, hi::Tuple{}) = $T{N}(lo, ntuple(d->0,Val{N}))
        (::Type{$T}){N}(lo::Tuple{}, hi::Dims{N}) = $T{N}(ntuple(d->0,Val{N}), hi)
        (::Type{$T}){N}(inds::Indices{N}) = $T{N}(map(lo,inds), map(hi,inds))
        (::Type{$T{N}}){N}(lo::AbstractVector, hi::AbstractVector) = $T{N}((lo...,), (hi...,))
        (::Type{$T})(lo::AbstractVector, hi::AbstractVector) = $T((lo...,), (hi...,)) # not inferrable

        (p::$T{0})(kernel, img, ::Alg) = p(kernel)
        (p::$T{0})(kernel) = $T(calculate_padding(kernel))
    end
end

padarray(img, border::Inner) = padarray(eltype(img), img, border)
padarray{T}(::Type{T}, img::AbstractArray{T}, border::Inner) = copy(img)
padarray{T}(::Type{T}, img::AbstractArray, border::Inner) = copy!(similar(Array{T}, indices(img)), img)

"""
    Fill(val)
    Fill(val, lo, hi)

Pad the edges of the image with a constant value, `val`.

Optionally supply the extent of the padding, see `Pad`.

# Example:

    imfilter(img, kernel, Fill(zero(eltype(img))))
"""
immutable Fill{T,N} <: AbstractBorder
    value::T
    lo::Dims{N}
    hi::Dims{N}

    Fill(value::T) = new(value)
    Fill(value::T, lo::Dims{N}, hi::Dims{N}) = new(value, lo, hi)
end

Fill{T}(value::T) = Fill{T,0}(value)
Fill{T,N}(value::T, lo::Dims{N}, hi::Dims{N}) = Fill{T,N}(value, lo, hi)
Fill(value, lo::AbstractVector, hi::AbstractVector) = Fill(value, (lo...,), (hi...,))
Fill{T,N}(value::T, inds::Base.Indices{N}) = Fill{T,N}(value, map(lo,inds), map(hi,inds))
Fill(value, kernel) = Fill(value, calculate_padding(kernel))

(p::Fill)(kernel) = Fill(p.value, kernel)
(p::Fill)(kernel, img, ::Alg) = Fill(p.value, kernel)
function (p::Fill)(kernel, img, ::FFT)
    inds = calculate_padding(kernel)
    newinds = map(padfft, inds, map(length, indices(img)))
    Fill(p.value, newinds)
end

function padarray{T}(::Type{T}, img::AbstractArray, border::Fill)
    throw(ArgumentError("$border lacks the proper padding sizes for an array with $(ndims(img)) dimensions"))
end
function padarray{T,S,_,N}(::Type{T}, img::AbstractArray{S,N}, f::Fill{_,N})
    A = similar(arraytype(img, T), map((l,r,h)->first(r)-l:last(r)+h, f.lo, indices(img), f.hi))
    try
        fill!(A, f.value)
    catch
        error("Unable to fill! an array of element type $(eltype(A)) with the value $(f.value). Supply an appropriate value to `Fill`, such as `zero(eltype(A))`.")
    end
    A[indices(img)...] = img
    A
end
padarray(img::AbstractArray, f::Fill) = padarray(eltype(img), img, f)

# There are other ways to define these, but using `mod` makes it safe
# for cases where the padding is bigger than length(inds)
"""
    padindex(border::Pad, lo::Integer, inds::AbstractUnitRange, hi::Integer)

Generate an index-vector to be used for padding. `inds` specifies the image indices along a particular axis; `lo` and `hi` are the amount to pad on the lower and upper, respectively, sides of this axis. `border` specifying the style of padding.

You can specialize this for custom `Pad{:name}` types to generate
arbitrary (cartesian) padding types.
"""
function padindex(border::Pad, lo::Integer, inds::AbstractUnitRange, hi::Integer)
    if border.style == :replicate
        return vcat(fill(first(inds), lo), PinIndices(inds), fill(last(inds), hi))
    elseif border.style == :circular
        return modrange(extend(lo, inds, hi), inds)
    elseif border.style == :symmetric
        I = [inds; reverse(inds)]
        r = modrange(extend(lo, inds, hi), 1:2*length(inds))
        return I[r]
    elseif border.style == :reflect
        I = [inds; last(inds)-1:-1:first(inds)+1]
        return I[modrange(extend(lo, inds, hi), 1:2*length(inds)-2)]
    else
        error("border style $(border.style) unrecognized")
    end
end

function inner(lo::Integer, inds::AbstractUnitRange, hi::Integer)
    first(inds)+lo:last(inds)-hi
end

lo(o::Integer) = max(-o, zero(o))
lo(r::AbstractUnitRange) = lo(first(r))

hi(o::Integer) = max(o, zero(o))
hi(r::AbstractUnitRange) = hi(last(r))

# extend(lo::Integer, inds::AbstractUnitRange, hi::Integer) = CatIndices.URange(first(inds)-lo, last(inds)+hi)
function extend(lo::Integer, inds::AbstractUnitRange, hi::Integer)
    newind = first(inds)-lo:last(inds)+hi
    OffsetArray(newind, newind)
end

calculate_padding(kernel) = indices(kernel)
@inline function calculate_padding(kernel::Tuple{Any, Vararg{Any}})
    inds = accumulate_padding(indices(kernel[1]), tail(kernel)...)
    if hasiir(kernel) && hasfir(kernel)
        inds = map(doublepadding, inds)
    end
    inds
end

hasiir(kernel) = _hasiir(false, kernel...)
_hasiir(ret) = ret
_hasiir(ret, kern, kernel...) = _hasiir(ret, kernel...)
_hasiir(ret, kern::AnyIIR, kernel...) = true

hasfir(kernel) = _hasfir(false, kernel...)
_hasfir(ret) = ret
_hasfir(ret, kern, kernel...) = true
_hasfir(ret, kern::AnyIIR, kernel...) = _hasfir(ret, kernel...)

function doublepadding(ind::AbstractUnitRange)
    f, l = first(ind), last(ind)
    f = f < 0 ? 2f : f
    l = l > 0 ? 2l : l
    f:l
end

accumulate_padding(inds::Indices, kernel1, kernels...) =
    accumulate_padding(expand(inds, indices(kernel1)), kernels...)
accumulate_padding(inds::Indices) = inds

modrange(x, r::AbstractUnitRange) = mod(x-first(r), length(r))+first(r)
modrange(A::AbstractArray, r::AbstractUnitRange) = map(x->modrange(x, r), A)

arraytype{T}(A::AbstractArray, ::Type{T}) = Array{T}  # fallback
arraytype(A::BitArray, ::Type{Bool}) = BitArray

interior(A, kernel) = _interior(indices(A), indices(kernel))
interior(A, factkernel::Tuple) = _interior(indices(A), accumulate_padding(indices(factkernel[1]), tail(factkernel)...))
function _interior{N}(indsA::NTuple{N}, indsk)
    indskN = fill_to_length(indsk, 0:0, Val{N})
    map(intersect, indsA, shrink(indsA, indsk))
end

next_shrink(inds::Indices, ::Tuple{}) = inds
function next_shrink(inds::Indices, kernel::Tuple)
    kern = first(kernel)
    iscopy(kern) && return next_shrink(inds, tail(kernel))
    shrink(inds, kern)
end

"""
    expand(inds::Indices, kernel)
    expand(inds::Indices, indskernel::Indices)

Expand an image region `inds` to account for necessary padding by `kernel`.
"""
expand(inds::Indices, kernel) = expand(inds, calculate_padding(kernel))
expand(inds::Indices, pad::Indices) = firsttype(map_copytail(expand, inds, pad))
expand(ind::AbstractUnitRange, pad::AbstractUnitRange) = oftype(ind, first(ind)+first(pad):last(ind)+last(pad))
expand(ind::Base.OneTo, pad::AbstractUnitRange) = expand(UnitRange(ind), pad)

"""
    shrink(inds::Indices, kernel)
    shrink(inds::Indices, indskernel)

Remove edges from an image region `inds` that correspond to padding needed for `kernel`.
"""
shrink(inds::Indices, kernel) = shrink(inds, calculate_padding(kernel))
shrink(inds::Indices, pad::Indices) = firsttype(map_copytail(shrink, inds, pad))
shrink(ind::AbstractUnitRange, pad::AbstractUnitRange) = oftype(ind, first(ind)-first(pad):last(ind)-last(pad))
shrink(ind::Base.OneTo, pad::AbstractUnitRange) = shrink(UnitRange(ind), pad)

allocate_output{T}(::Type{T}, img, kernel, border) = similar(img, T)
function allocate_output{T}(::Type{T}, img, kernel, ::Inner{0})
    inds = interior(img, kernel)
    similar(img, T, inds)
end
allocate_output(img, kernel, border) = allocate_output(filter_type(img, kernel), img, kernel, border)

"""
    map_copytail(f, a::Tuple, b::Tuple)

Apply `f` to paired elements of `a` and `b`, copying any tail elements
when the lengths of `a` and `b` are not equal.
"""
@inline map_copytail(f, a::Tuple, b::Tuple) = (f(a[1], b[1]), map_copytail(f, tail(a), tail(b))...)
map_copytail(f, a::Tuple{}, b::Tuple{}) = ()
map_copytail(f, a::Tuple, b::Tuple{}) = a
map_copytail(f, a::Tuple{}, b::Tuple) = b

function firsttype(t::Tuple)
    T = typeof(t[1])
    map(T, t)
end
firsttype(::Tuple{}) = ()

# end