text:
module CLFFT

import OpenCL.cl

include("api.jl")
include("error.jl")

const SP_MAX_LEN = 1 << 24
const DP_MAX_LEN = 1 << 22

immutable CLFFTError
    code::Int
    desc::Symbol
    function CLFFTError(c::Integer)
        code = Int(c)
        new(code, _clfft_error_codes[code])
    end
end

Base.show(io::IO, err::CLFFTError) =
        Base.print(io, "CLFFTError(code=$(err.code), :$(err.desc))")

macro check(clfftFunc)
    quote
        local err::Cint
        err = $clfftFunc
        if err != api.CLFFT_SUCCESS
            throw(CLFFTError(err))
        end
        err
    end
end

version() = begin
    major = cl.CL_uint[0]
    minor = cl.CL_uint[0]
    patch = cl.CL_uint[0]
    api.clfftGetVersion(major, minor, patch)
    return VersionNumber(Int(major[1]), Int(minor[1]), Int(patch[1]))
end

function supported_radices()
  v = version()
  radices = [2,3,5]
  v ≥ v"2.8.0"  && push!(radices, 7)
  v ≥ v"2.12.0" && push!(radices, 11, 13)

  radices
end

# Module level library handle,
# when module is GC'd, the finalizer for SetupData is
# called to teardown the library state
const __handle = begin
    v = version()
    api.SetupData(v.major, v.minor, v.patch, 0)
end

# clFFT floating-point types:
typealias clfftNumber Union{Float64,Float32,Complex128,Complex64}
typealias clfftReal Union{Float64,Float32}
typealias clfftComplex Union{Complex128,Complex64}
typealias clfftDouble Union{Float64,Complex128}
typealias clfftSingle Union{Float32,Complex64}
typealias clfftTypeDouble Union{Type{Float64},Type{Complex128}}
typealias clfftTypeSingle Union{Type{Float32},Type{Complex64}}

typealias PlanHandle Csize_t

type Plan{T<:clfftNumber}
    # boxed handle (most api functions need address, setup/teardown need pointer)
    id::Array{PlanHandle,1}

    function Plan(plan::Array{PlanHandle,1})
        p = new(plan)
        finalizer(p, x -> begin
            if x.id[1] != C_NULL
                @check api.clfftDestroyPlan(x.id)
            end
            x.id[1] = C_NULL
        end)
        return p
    end
end

function Plan{T<:clfftNumber}(::Type{T}, ctx::cl.Context, sz::Dims)
    if length(sz) > 3
        throw(ArgumentError("Plans can have dimensions of 1,2 or 3"))
    end
    ndim = length(sz)
    lengths = Csize_t[0, 0, 0]
    total_length = 1
    radices = supported_radices()
    for i in 1:ndim
        s = sz[i]
        fs = keys(factor(s))
        if fs ⊆ radices
            lengths[i] = s
            total_length *= s
        else
            throw(ArgumentError("""Plans can only have dims that are
                                   powers of $(radices)"""))
        end
    end
    if T <: clfftSingle
        if total_length > SP_MAX_LEN
            throw(ArgumentError("""CLFFT supports single precision transform
                                   lengths up to $(SP_MAX_LEN)"""))
        end
    else
        if total_length > DP_MAX_LEN
            throw(ArgumentError("""clFFT supports double precision transform
                                 lengths up to $(DP_MAX_LEN)"""))
        end
    end
    ph = PlanHandle[0]
    err = api.clfftCreateDefaultPlan(ph, ctx.id,
                                     Int32(ndim), lengths)
    if err != api.CLFFT_SUCCESS
        if ph[1] != 0
            @check api.clfftDestroyPlan(ph)
        end
        throw(CLFFTError(err))
    end
    if T <: clfftSingle
        @check api.clfftSetPlanPrecision(ph[1], api.CLFFT_SINGLE)
    else
        @check api.clfftSetPlanPrecision(ph[1], api.CLFFT_DOUBLE)
    end
    @assert ph[1] != 0
    return Plan{T}(ph)
end

function Plan{T<:clfftNumber}(::Type{T}, ctx::cl.Context,
                              input::StridedArray{T}, region)
    ndim = length(region)
    if ndim > 3
        throw(ArgumentError("Plans can have dimensions of 1, 2, or 3"))
    end
    if ndims(input) != ndim
        throw(ArgumentError("input array and region must have the same dimensionality"))
    end
    insize = size(input)
    lengths = Csize_t[0,0,0]
    total_length = 1
    radices = supported_radices()
    for i in 1:ndim
        s = insize[i]
        fs = keys(factor(s))
        if fs ⊆ radices
            lengths[i] = s
            total_length *= s
        else
            throw(ArgumentError("""Plans can only have dims that are
                                   powers of $(radices)"""))
        end
    end
    if T <: clfftSingle
        if total_length > SP_MAX_LEN
            throw(ArgumentError("""CLFFT supports single precision transform
                                   lengths up to $(SP_MAX_LEN)"""))
        end
    else
        if total_length > DP_MAX_LEN
            throw(ArgumentError("""clFFT supports double precision transform
                                 lengths up to $(DP_MAX_LEN)"""))
        end
    end
    ph = PlanHandle[0]
    err = api.clfftCreateDefaultPlan(ph, ctx.id, Int32(ndim), lengths)
    if err != api.CLFFT_SUCCESS
        if ph[1] != 0
            @check api.clfftDestroyPlan(ph)
        end
        throw(CLFFTError(err))
    end
    if T <: clfftSingle
        @check api.clfftSetPlanPrecision(ph[1], api.CLFFT_SINGLE)
    else
        @check api.clfftSetPlanPrecision(ph[1], api.CLFFT_DOUBLE)
    end
    @assert ph[1] != 0
    plan = Plan{T}(ph)

    tdim = ndim
    tstrides = strides(input)
    # TODO : this works for dense arrays
    # need to test out for arbitrary regions
    tdistance = 0
    tbatchsize = 1

    set_layout!(plan, :interleaved, :interleaved)
    set_instride!(plan, tstrides)
    set_outstride!(plan, tstrides)
    set_distance!(plan, tdistance, tdistance)
    set_batchsize!(plan, tbatchsize)

    return plan
end

Plan{T<:clfftNumber}(::Type{T}, ctx::cl.Context, input::StridedArray{T}) =
        Plan(T, ctx, input, 1:ndims(input))

precision(p::Plan) = begin
    res = Cint[0]
    @check api.clfftGetPlanPrecision(p.id[1], res)
    if res[1] == 1
        return :single
    elseif res[1] == 2
        return :double
    elseif res[1] == 3
        return :single_fast
    elseif res[1] == 4
        return :double_fast
    else
        error("undefined")
    end
end

set_precision!(p::Plan, v::Symbol) = begin
    if v == :single
        @check api.clfftSetPlanPrecision(p.id[1], api.CLFFT_SINGLE)
    elseif v == :double
        @check api.clfftSetPlanPrecision(p.id[1], api.CLFFT_DOUBLE)
    elseif v == :single_fast
        @check api.clfftSetPlanPrecision(p.id[1], api.CLFFT_SINGLE_FAST)
    elseif v == :double_fast
        @check api.clfftSetPlanPrecision(p.id[1], api.CLFFT_DOUBLE_FAST)
    else
        error("unknown precision $v")
    end
    return p
end

layout(p::Plan) = begin
    i = Cint[0]
    o = Cint[0]
    @check api.clfftGetLayout(p.id[1], i, o)
    lout = x -> begin
        if x == 1
            return :interleaved
        elseif x == 2
            return :planar
        elseif x == 3
            return :hermitian_interleaved
        elseif x == 4
            return :hermitian_planar
        elseif x == 5
            return :real
        else
            error("undefined")
        end
    end
    (lout(i[1]), lout(o[1]))
end

set_layout!(p::Plan, in::Symbol, out::Symbol) = begin
    args = Cint[0, 0]
    if in == :interleaved
        args[1] = api.CLFFT_COMPLEX_INTERLEAVED
    elseif in == :planar
        args[1] = api.CLFFT_COMPLEX_PLANAR
    else
        throw(ArgumentError("in must be :interleaved or :planar"))
    end
    if in == :interleaved
        args[2] = api.CLFFT_COMPLEX_INTERLEAVED
    elseif in == :planar
        args[2] = api.CLFFT_COMPLEX_PLANAR
    else
        throw(ArgumentError("in must be :interleaved or :planar"))
    end
    @check api.clfftSetLayout(p.id[1], args[1], args[2])
    return p
end

result(p::Plan) = begin
    res = Cint[0]
    @check api.clfftGetResultLocation(p.id[1], res)
    if res[1] == 1
        return :inplace
    elseif res[1] == 2
        return :outofplace
    else
        error("undefined")
    end
end

set_result!(p::Plan, v::Symbol) = begin
    if v == :inplace
        @check api.clfftSetResultLocation(p.id[1], api.CLFFT_INPLACE)
    elseif v == :outofplace
        @check api.clfftSetResultLocation(p.id[1], api.CLFFT_OUTOFPLACE)
    else
        throw(ArgumentError("set_result must be :inplace or :outofplace"))
    end
    return p
end

scaling_factor(p::Plan, dir::Symbol) = begin
    res = Cint[0]
    d::Cint
    if dir == :forward
        d = Int32(-1)
    elseif dir == :backward
        d = Int32(1)
    else
        error("undefined")
    end
    scale = Float32[0]
    @check api.clfftGetPlanScale(p.id[1], d, scale)
    return scale[1]
end

set_scaling_factor!(p::Plan, dir::Symbol, f::AbstractFloat) = begin
    if dir == :forward
        d = Int32(-1)
    elseif d == :backward
        d = Int32(1)
    else
        error("undefined")
    end
    @check api.clfftsetPlanScale(p.id[1], d, float32(f))
    return p
end

set_batchsize!(p::Plan, n::Integer) = begin
    @assert n > 0
    @check api.clfftSetPlanBatchSize(p.id[1], convert(Csize_t, n))
    return p
end

batchsize(p::Plan) = begin
    res = Csize_t[0]
    @check api.clfftGetPlanBatchSize(p.id[1], res)
    return Int(res[1])
end

set_dim!(p::Plan, d::Integer) = begin
    @assert d > 0 && d <= 3
    @check api.clfftSetPlanDim(p.id[1], convert(Csize_t, d))
    return p
end

dim(p::Plan) = begin
    res  = Int32[0]
    size = cl.CL_uint[0]
    @check api.clfftGetPlanDim(p.id[1], res, size)
    return Int(res[1])
end

set_lengths!(p::Plan, dims::Dims) = begin
    ndim = length(dims)
    @assert ndim <= 3
    nd = Array(Csize_t, ndim)
    for (i, d) in enumerate(dims)
        nd[i] = d
    end
    @check api.clfftSetPlanLength(p.id[1], Int32(ndim), nd)
    return p
end

lengths(p::Plan) = begin
    d = dim(p)
    res = Array(Csize_t, d)
    @check api.clfftGetPlanLength(p.id[1], Int32(d), res)
    return map(Int, res)
end

instride(p::Plan) = begin
    d = dim(p)
    res = Array(Csize_t, d)
    @check api.clfftGetPlanInStride(p.id[1], Int32(d), res)
    return map(Int, res)
end

set_instride!(p::Plan, instrides) = begin
    d = length(instrides)
    @assert d == dim(p)
    strides = Csize_t[Int(s) for s in instrides]
    @check api.clfftSetPlanInStride(p.id[1], Int32(d), strides)
    return p
end

outstride(p::Plan) = begin
    d = dim(p)
    res = Array(Csize_t, d)
    @check api.clfftGetPlanOutStride(p.id[1], Int32(d), res)
    return map(Int, res)
end

set_outstride!(p::Plan, outstrides) = begin
    d = length(outstrides)
    @assert d == dim(p)
    strides = Csize_t[Int(s) for s in outstrides]
    @check api.clfftSetPlanInStride(p.id[1], Int32(d), strides)
    return p
end

distance(p::Plan) = begin
    indist = Csize_t[0]
    odist  = Csize_t[0]
    @check api.clfftGetPlanDistance(p.id[1], indist, odist)
    return (indist[1], odist[1])
end

set_distance!(p::Plan, indist::Integer, odist::Integer) = begin
    @assert indist >= 0 && odist >= 0
    i = UInt(indist)
    o = UInt(odist)
    @check api.clfftSetPlanDistance(p.id[1], i, o)
    return p
end

transpose_result(p::Plan) = begin
    res = Cint[0]
    @check api.clfftGetPlanTransposeResult(p.id[1], res)
    if res[1] == api.CLFFT_NOTRANSPOSE
        return false
    elseif res[1] == api.CLFFT_TRANSPOSED
        return true
    else
        error("undefined")
    end
end

set_transpose_result(p::Plan, transpose::Bool) = begin
    if transpose
        @check api.clfftSetTransposeResult(p.id[1], d,
                cl.CLFFT_TRANSPOSED)
    else
        @check api.clfftSetTransposeResult(p.id[1], d,
                cl.CLFFT_NOTRANSPOSE)
    end
end

tmp_buffer_size(p::Plan) = begin
    res = Csize_t[0]
    @check api.clfftGetTmpBufSize(p.id[1], res)
    return Int(res[1])
end

context(p::Plan) = begin
    res = Array(cl.CL_context)
    @check api.clfftGetPlanContext(p.id[1], res)
    return cl.Context(res[1])
end

bake!(p::Plan, qs::Vector{cl.CmdQueue}) = begin
    nqueues = length(qs)
    q_ids = [q.id for q in qs]
    @check api.clfftBakePlan(p.id[1], nqueues, q_ids, C_NULL, C_NULL)
    return p
end
bake!(p::Plan, q::cl.CmdQueue) = bake!(p, [q])

function enqueue_transform{T<:clfftNumber}(p::Plan,
                                           dir::Symbol,
                                           qs::Vector{cl.CmdQueue},
                                           input::cl.Buffer{T},
                                           output::Union{Void,cl.Buffer{T}};
                                           wait_for::Union{Void,Vector{cl.Event}}=nothing,
                                           tmp::Union{Void,cl.Buffer{T}}=nothing)
    if dir != :forward && dir != :backward
        throw(ArgumentError("Unknown direction $dir"))
    end
    q_ids = [q.id for q in qs]
    in_buff_ids  = [input.id]
    out_buff_ids = C_NULL
    if output != nothing
        out_buff_ids = [output.id]
    end
    nevts = 0
    evt_ids = C_NULL
    if wait_for != nothing
        nevts = length(wait_for)
        evt_ids = [evt.id for evt in wait_for]
    end
    nqueues = length(q_ids)
    out_evts = Array(cl.CL_event, nqueues)
    tmp_buffer = C_NULL
    if tmp != nothing
        tmp_buff_id = [tmp.id]
    end
    @check api.clfftEnqueueTransform(
                              p.id[1],
                              dir == :forward ? api.CLFFT_FORWARD : api.CLFFT_BACKWARD,
                              UInt32(nqueues),
                              q_ids,
                              UInt32(nevts),
                              evt_ids,
                              out_evts,
                              in_buff_ids,
                              out_buff_ids,
                              tmp_buffer)
    return [cl.Event(e_id) for e_id in out_evts]
end

end # module


text:
"""
    HelmholtzOutput
Return type for Helmholtz computations.
Fields are `nothing` on return unless requested.
See individual FMM/direct computation function
documentation for specifics.
"""
mutable struct HelmholtzOutput <: FMMVals
    pot
    grad

    pottarg
    gradtarg

    pgt
    pg
end
function HelmholtzOutput(pgt,pg=nothing)
    HelmholtzOutput(nothing,nothing,
                    nothing,nothing,
                    pgt,pg)
end

"""
```julia
    vals = hfmm3d(eps,zk,sources;charges=nothing,dipvecs=nothing,
                  targets=nothing,pg=0,pgt=0,nd=1)
```
This function computes the N-body Helmholtz interactions
in three dimensions where the interaction kernel is given
by ``e^{ikr}/r`` and its gradients. This is the
``O(N)`` fast multipole code which computes the interactions
to the requested precision.
```math
 u(x) = \\sum_{j=1}^{N} c_{j} \\frac{e^{ik \\|x-x_{j}\\|}}{\\|x-x_{j}\\|}
- v_{j} \\cdot \\nabla \\left( \\frac{e^{ik \\|x-x_{j}\\|}}{\\|x-x_{j}\\|}
\\right)  \\, ,
```

where ``c_{j}`` are the charge densities,
      ``v_{j}`` are the dipole orientation vectors, and
      ``x_{j}`` are the source locations.
      When ``x=x_{m}``, the term corresponding to
      ``x_{m}`` is dropped from the sum
# Input
* `eps::Float64` precision requested
* `zk::ComplexF64` Helmholtz parameter
* `sources::Array{Float64}` size (3,n) source locations (``x_{j}``)
* `charges::Array{ComplexF64}` size (nd,n) or (n) charge densities (c_{j})
* `dipvecs::Array{ComplexF64}` size (nd,3,n) or (3,n)) dipole orientation vectors (``v_{j}``)
* `targets::Array{Float64}` size (3,nt) target locations (``x``)
* `pg::Integer` source eval flag.
    + Potential (``u``) at sources evaluated if `pg == 1`.
    + Potential and gradient (``\\nabla u``) at sources evaluated if `pg == 2`
* `pgt::Integer` target eval flag.
    + Potential at targets evaluated if `pgt == 1`.
    + Potential and gradient at targets evaluated if `pgt == 2`
* `nd::Integer` number of densities
# Output

`vals::HelmholtzOutput` with the fields
* `vals.pot::Array{ComplexF64}` size (nd,n) or (n) potential at source locations if requested
* `vals.grad::Array{ComplexF64}` size (nd,3,n) or (3,n) gradient at source locations if requested
* `vals.pottarg::Array{ComplexF64}` size (nd,nt) or (nt) potential at target locations if requested
* `vals.gradtarg::Array{ComplexF64}` size (nd,3,nt) or (3,nt) gradient at target locations if requested
"""
function hfmm3d(eps::Float64,zk::Union{Float64,ComplexF64},
                sources::Array{Float64};
                charges::TCN=nothing,dipvecs::TCN=nothing,
                targets::TFN=nothing,pg::Integer=0,pgt::Integer=0,
                nd::Integer=1)

    # check inputs

    @assert size(sources,1) == 3
    @assert nd >= 0
    @assert (0 <= pg  && pg  <= 2) "flag pg not in expected range"
    @assert (0 <= pgt && pgt <= 2) "flag pgt not in expected range"

    # default values

    vals = HelmholtzOutput(pgt,pg)

    ifcharge = 0
    ifdipole = 0

    zero = ComplexF64(0)

    pot = zero
    grad = zero
    hess = zero
    pottarg = zero
    gradtarg = zero
    hesstarg = zero

    n = div(length(sources),3)
    nt = 0
    zk = complex(zk)

    if targets !== nothing
        @assert size(targets,1) == 3
        nt = div(length(targets),3)
    else
        targets = 0.0
    end

    if charges !== nothing
        @assert div(length(charges),nd) == n
        ifcharge = 1
    else
        charges = zero
        ifcharge = 0
    end

    if dipvecs !== nothing
        @assert div(length(dipvecs),nd) == n*3
        ifdipole = 1
    else
        dipvecs = zero
        ifdipole = 0
    end

    if ifcharge == 0 && ifdipole == 0
        @warn "no charges or dipoles provided, doing nothing"
        return vals
    end

    if (pg != 1 && pg != 2) && (pgt != 1 && pgt != 2)
        @warn "no output requested, doing nothing"
        return vals
    end

    if (pgt == 1 || pgt == 2) && targets == nothing
        @warn "target values requested but no targets provided"
    end

    # allocate memory for return values

    if pg == 1 || pg == 2
        if nd > 1
            pot = zeros(ComplexF64,nd,n)
        else
            pot = zeros(ComplexF64,n)
        end
    end

    if pg == 2
        if nd > 1
            grad = zeros(ComplexF64,nd,3,n)
        else
            grad = zeros(ComplexF64,3,n)
        end
    end

    if pgt == 1 || pgt == 2
        if nd > 1
            pottarg = zeros(ComplexF64,nd,nt)
        else
            pottarg = zeros(ComplexF64,nt)
        end
    end

    if pgt == 2
        if nd > 1
            gradtarg = zeros(ComplexF64,nd,3,nt)
        else
            gradtarg = zeros(ComplexF64,3,nt)
        end
    end

    # actually call the function
    #
    # fortran calling sequence:
    # subroutine hfmm3d(nd,eps,zk,nsource,source,ifcharge,
    #     $    charge,ifdipole,dipvec,ifpgh,pot,grad,hess,ntarg,
    #     $    targ,ifpghtarg,pottarg,gradtarg,hesstarg)
    #

    ccall((:hfmm3d_,libfmm3d),Cvoid,(Fi,Fd,Fc,Fi,Fd,Fi,Fc,
                                    Fi,Fc,Fi,Fc,Fc,Fc,Fi,
                                    Fd,Fi,Fc,Fc,Fc),
          nd,eps,zk,n,sources,ifcharge,charges,ifdipole,
          dipvecs,pg,pot,grad,hess,nt,targets,pgt,
          pottarg,gradtarg,hesstarg)

    # load requested values

    if pg == 1 || pg == 2; vals.pot = pot end
    if pg == 2; vals.grad = grad end
    if pgt == 1 || pgt == 2; vals.pottarg = pottarg end
    if pgt == 2; vals.gradtarg = gradtarg end

    return vals

end

"""
```julia
    vals = h3ddir(zk,sources,targets;charges=nothing,
                    dipvecs=nothing,pgt=0,nd=1,
                    thresh=1e-16)
```
This function computes the N-body Helmholtz interactions
in three dimensions where the interaction kernel is given
by ``e^{ikr}/r`` and its gradients. This is the
``O(N^2)`` direct evaluation code.
```math
 u(x) = \\sum_{j=1}^{N} c_{j} \\frac{e^{ik \\|x-x_{j}\\|}}{\\|x-x_{j}\\|}
- v_{j} \\cdot \\nabla \\left( \\frac{e^{ik \\|x-x_{j}\\|}}{\\|x-x_{j}\\|}
\\right)  \\, ,
```

where ``c_{j}`` are the charge densities,
      ``v_{j}`` are the dipole orientation vectors, and
      ``x_{j}`` are the source locations.
      When ``x=x_{m}``, the term corresponding to
      ``x_{m}`` is dropped from the sum
# Input
* `zk::ComplexF64` Helmholtz parameter
* `sources::Array{Float64}` size (3,n) source locations (``x_{j}``)
* `targets::Array{Float64}` size (3,nt) target locations (``x``)
* `charges::Array{ComplexF64}` size (nd,n) or (n) charge densities (c_{j})
* `dipvecs::Array{ComplexF64}` size (nd,3,n) or (3,n)) dipole orientation vectors (``v_{j}``)
* `pgt::Integer` target eval flag.
    + Potential at targets evaluated if `pgt == 1`.
    + Potential and gradient at targets evaluated if `pgt == 2`
* `nd::Integer` number of densities
* `thresh::Float64` threshold for ignoring interactions when ``\\|x-x_{j}\\| \\leq thresh``
# Output

`vals::HelmholtzOutput` with the fields
* `vals.pottarg::Array{ComplexF64}` size (nd,nt) or (nt) potential at target locations if requested
* `vals.gradtarg::Array{ComplexF64}` size (nd,3,nt) or (3,nt) gradient at target locations if requested
"""
function h3ddir(zk::Union{ComplexF64,Float64},sources::Array{Float64},
                targets::Array{Float64};
                charges::TCN=nothing,dipvecs::TCN=nothing,
                pgt::Integer=0,nd::Integer=1,
                thresh::Float64=1e-16)

    # check inputs

    @assert size(sources,1) == 3
    @assert size(targets,1) == 3
    @assert nd >= 0
    @assert (0 <= pgt || pgt <= 2) "flag pgt not in expected range"

    # default values

    vals = HelmholtzOutput(pgt)

    ifcharge = 0
    ifdipole = 0

    zero = ComplexF64(0)

    pottarg = zero
    gradtarg = zero
    hesstarg = zero

    n = div(length(sources),3)
    nt = div(length(targets),3)
    zk = complex(zk)

    if charges !== nothing
        @assert div(length(charges),nd) == n
        ifcharge = 1
    else
        charges = zero
        ifcharge = 0
    end

    if dipvecs !== nothing
        @assert div(length(dipvecs),nd) == n*3
        ifdipole = 1
    else
        dipvecs = zero
        ifdipole = 0
    end

    if ifcharge == 0 && ifdipole == 0
        @warn "no charges or dipoles provided, doing nothing"
        return vals
    end

    if (pgt != 1 && pgt != 2)
        @warn "no output requested, doing nothing"
        return vals
    end

    # allocate memory for return values

    if pgt == 1 || pgt == 2
        if nd > 1
            pottarg = zeros(ComplexF64,nd,nt)
        else
            pottarg = zeros(ComplexF64,nt)
        end
    end

    if pgt == 2
        if nd > 1
            gradtarg = zeros(ComplexF64,nd,3,nt)
        else
            gradtarg = zeros(ComplexF64,3,nt)
        end
    end

    # dispatch to appropriate wrapper

    if pgt == 1
        if ifcharge == 1
            if ifdipole == 1
                h3ddirectcdp!(nd,zk,sources,charges,
                              dipvecs,n,targets,nt,
                              pottarg,thresh)
                vals.pottarg = pottarg
            else
                h3ddirectcp!(nd,zk,sources,charges,
                             n,targets,nt,pottarg,
                             thresh)
                vals.pottarg = pottarg
            end
        else
            if ifdipole == 1
                h3ddirectdp!(nd,zk,sources,
                             dipvecs,n,targets,nt,
                             pottarg,thresh)
                vals.pottarg = pottarg
            end
        end
    elseif pgt == 2
        if ifcharge == 1
            if ifdipole == 1
                h3ddirectcdg!(nd,zk,sources,charges,
                              dipvecs,n,targets,nt,
                              pottarg,gradtarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
            else
                h3ddirectcg!(nd,zk,sources,charges,
                             n,targets,nt,pottarg,
                             gradtarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
            end
        else
            if ifdipole == 1
                h3ddirectdg!(nd,zk,sources,
                              dipvecs,n,targets,nt,
                             pottarg,gradtarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
            end
        end
    end

    return vals

end

function h3ddirectcp!(nd,zk,sources,charges,n,targets,nt,
                      pot,thresh)
    ccall((:h3ddirectcp_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fi,
                                          Fd,Fi,Fc,Fd),
          nd,zk,sources,charges,n,targets,nt,pot,thresh)
    return
end

function h3ddirectdp!(nd,zk,sources,dipvecs,n,targets,nt,
                      pot,thresh)
    ccall((:h3ddirectdp_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fi,
                                          Fd,Fi,Fc,Fd),
          nd,zk,sources,dipvecs,n,targets,nt,pot,thresh)
    return
end

function h3ddirectcdp!(nd,zk,sources,charges,dipvecs,n,targets,
                       nt,pot,thresh)
    ccall((:h3ddirectcdp_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fc,Fi,
                                          Fd,Fi,Fc,Fd),
          nd,zk,sources,charges,dipvecs,n,targets,nt,pot,thresh)
    return
end

function h3ddirectcg!(nd,zk,sources,charges,n,targets,nt,
                      pot,grad,thresh)
    ccall((:h3ddirectcg_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fi,
                                          Fd,Fi,Fc,Fc,Fd),
          nd,zk,sources,charges,n,targets,nt,pot,grad,thresh)
    return
end

function h3ddirectdg!(nd,zk,sources,dipvecs,n,targets,nt,
                      pot,grad,thresh)
    ccall((:h3ddirectdg_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fi,
                                          Fd,Fi,Fc,Fc,Fd),
          nd,zk,sources,dipvecs,n,targets,nt,pot,grad,thresh)
    return
end

function h3ddirectcdg!(nd,zk,sources,charges,dipvecs,n,targets,
                       nt,pot,grad,thresh)
    ccall((:h3ddirectcdg_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fc,Fi,
                                          Fd,Fi,Fc,Fc,Fd),
          nd,zk,sources,charges,dipvecs,n,targets,nt,pot,grad,thresh)
    return
end

ext: jl


text:
fit(E::NFXP, U, S, D) = fit_nfxp(U, S, D, E)
fit_nfxp(U, S, D) = fit_nfxp(U, S, D, NFXP(U, S, D))
function fit_nfxp{Tu<:AbstractUtility, Ts<:Stocas.AbstractState}(U::Tu,
                  S::Ts,
                  D::Data,
                  E::NFXP;
                  ## Keyword arguments
                  # Convergence parameters
                  ε_outer = 1e-8,
                  ε_inner = 1e-11,
                  method = NewtonTrustRegion(),
                  options = nothing)
    options = options == nothing ? Optim.Options(show_trace=E.verbose, extended_trace=E.verbose, g_tol = ε_outer) : options

    conv_info = ConvergenceNFXP()

    # Cache variables
    c = NFXPCacheVars(U.nvar, D.nobs)
    c.Θᵏ .= U.Θ
    c.last_Θᵏ .= copy(c.Θᵏ)+1.

      llᵈ(x)           =   ll(x,                     c.last_Θᵏ, U, S, D, E, conv_info)
     ∇llᵈ(gradient, x) =  ∇ll!(gradient, c.score, x, c.last_Θᵏ, U, S, D, E, conv_info)
    ∇²llᵈ(hessian, x)  = ∇²ll!(hessian,  c.score, x, c.last_Θᵏ, U, S, D, E, conv_info)

    likelihood = TwiceDifferentiable(llᵈ, ∇llᵈ, ∇²llᵈ, c.Θᵏ)
    result = optimize(likelihood, c.Θᵏ, method, options)
    # Return gradient and hessian at estimated parameters
    gradient = zeros(U.nvar)
    hessian = zeros(U.nvar, U.nvar)
     ∇llᵈ(gradient, c.Θᵏ)
    ∇²llᵈ(hessian, c.Θᵏ)
    EstimationResults(E,
                      -Optim.minimum(result)*D.nobs,
                      -gradient*D.nobs,
                      -hessian*D.nobs,
                      likelihood,
                      Optim.minimizer(result),
                      conv_info,
                      result,
                      D.nobs,
                      likelihood)
end

""" Updates the choice specific utilities based on most recent parameters. """
update_utilities!(U) = error("Need to implement update_utilities! for $(typeof(U))")
function update_utilities!(U::LinearUtility)
    for i in 1:length(U.Z)
        U.U[i] .= U.Z[i]*U.Θ
    end
end


text:
using Judycon, Test

#
# 1   2---3   4---5
# |   |   |   |   |
# 6---7   8   9   10
#

G1 = QuickFind(10);
G2 = QuickUnion(10, false, false);
G3 = QuickUnion(10, true, false);
G4 = QuickUnion(10, false, true);
G5 = QuickUnion(10, true, true);

function test(G)

    connect!(G, 1, 6)
    connect!(G, 6, 7)
    connect!(G, 7, 2)
    connect!(G, 2, 3)
    connect!(G, 3, 8)

    connect!(G, 9, 4)
    connect!(G, 4, 5)
    connect!(G, 5, 10)

    # `pts1` should form one set of connected components and `pts2` another, respectively.
    pts1 = [1, 2, 3, 6, 7, 8]
    pts2 = [4, 5, 9, 10]
    for p in pts1
        for q in pts1
            @test isconnected(G, p, q)
            @test isconnected(G, q, p)
        end
        for q in pts2
            @test !isconnected(G, p, q)
            @test !isconnected(G, q, p)
        end
    end
end

@testset "Test Judycon.jl" begin
    @testset "Test QuickFind" begin test(G1) end
    @testset "Test QuickUnion without weighting and path compression" begin test(G2) end
    @testset "Test QuickUnion with weighting and without path compression" begin test(G3) end
    @testset "Test QuickUnion without weighting and with path compression" begin test(G4) end
    @testset "Test QuickUnion with weighting and path compression" begin test(G5) end
end


text:
# using Pkg
# pkg"activate .."
# push!(LOAD_PATH,"../src/")
using ProjectFlow, Documenter

makedocs(;
    modules=[ProjectFlow],
    authors="Data Psycho",
    repo="https://github.com/DataPsycho/ProjectFlow.jl/blob/{commit}{path}#L{line}",
    sitename="ProjectFlow.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://DataPsycho.github.io/ProjectFlow.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/DataPsycho/ProjectFlow.jl",
)


text:
    """
    Cartesian 3×3 matrix

        Matrix3()
        Matrix3(a11,a12,a13,a21,a23,a23,a31,a32,a33)*
        Matrix3(matrix)

    (*) definition is row major, overriding the column major definition of `SMatrix{}`

    """
    const Matrix3 = SMatrix{3, 3, Float64}

    Matrix3() = Matrix3(fill(0.0, 3,3))
    # Define matrix elements in row major order
    Matrix3(a11,a12,a13,a21,a22,a23,a31,a32,a33) = Matrix3([a11 a12 a13; a21 a22 a23; a31 a32 a33])

    export Ô, Î

    """
    Zero 3×3 matrix ``\\hat{O} = \\pmatrix{0 & 0 & 0 \\\\ 0 & 0 & 0 \\\\ 0 & 0 & 0}``

        const Ô = Matrix3(
            0,0,0,
            0,0,0,
            0,0,0)

    """
    const Ô = Matrix3()
    """
    Identity 3×3 matrix ``\\hat{I} = \\pmatrix{1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & 0 & 1}``

        const Î = Matrix3(
            1,0,0,
            0,1,0,
            0,0,1)

    """
    const Î = Matrix3(I)

    function Base.inv(m::Matrix3)::Matrix3
        t2 = m[1,1]*m[2,2]*m[3,3]
        t3 = m[1,2]*m[2,3]*m[3,1]
        t4 = m[1,3]*m[2,1]*m[3,2]
        t7 = m[1,1]*m[2,3]*m[3,2]
        t8 = m[1,2]*m[2,1]*m[3,3]
        t9 = m[1,3]*m[2,2]*m[3,1]
        t6 = 1.0/(t2+t3+t4-t7-t8-t9)
        return Matrix3(
             t6*(m[2,2]*m[3,3]-m[2,3]*m[3,2]),
            -t6*(m[1,2]*m[3,3]-m[1,3]*m[3,2]),
             t6*(m[1,2]*m[2,3]-m[1,3]*m[2,2]),
            -t6*(m[2,1]*m[3,3]-m[2,3]*m[3,1]),
             t6*(m[1,1]*m[3,3]-m[1,3]*m[3,1]),
            -t6*(m[1,1]*m[2,3]-m[1,3]*m[2,1]),
             t6*(m[2,1]*m[3,2]-m[2,2]*m[3,1]),
            -t6*(m[1,1]*m[3,2]-m[1,2]*m[3,1]),
             t6*(m[1,1]*m[2,2]-m[1,2]*m[2,1]))
    end

    function solve(m::Matrix3, v::Vector3)::Vector3
        t = (-m[1,2]*m[3,3]+m[1,3]*m[3,2])*m[2,1]+(-m[2,2]*m[1,3]+m[2,3]*m[1,2])*m[3,1]+(m[2,2]*m[3,3]-m[2,3]*m[3,2])*m[1,1]
        return Vector3(
            (v[1]*( m[2,2]*m[3,3] - m[2,3]*m[3,2]) + v[2]*(-m[1,2]*m[3,3] + m[1,3]*m[3,2]) + v[3]*(-m[2,2]*m[1,3] + m[2,3]*m[1,2]))/t,
            (v[1]*(-m[2,1]*m[3,3] + m[2,3]*m[3,1]) + v[2]*( m[1,1]*m[3,3] - m[1,3]*m[3,1]) + v[3]*( m[2,1]*m[1,3] - m[2,3]*m[1,1]))/t,
            (v[1]*( m[2,1]*m[3,2] - m[2,2]*m[3,1]) + v[2]*(-m[1,1]*m[3,2] + m[1,2]*m[3,1]) + v[3]*(-m[2,1]*m[1,2] + m[2,2]*m[1,1]))/t)
    end

    function Base.:\(A::Matrix3, b::Vector3)::Vector3
        return solve(A,b)
    end

    function Base.getproperty(m::Matrix3, name::Symbol)
        if name===:A11; return m[1,1]; end
        if name===:A12; return m[1,2]; end
        if name===:A13; return m[1,3]; end
        if name===:A21; return m[2,1]; end
        if name===:A22; return m[2,2]; end
        if name===:A23; return m[2,3]; end
        if name===:A31; return m[3,1]; end
        if name===:A32; return m[3,2]; end
        if name===:A33; return m[3,3]; end
        if name===:Trace; return m[1,1]+m[2,2]+m[3,3]; end
        if name===:Determinat; return m[1,1]*(m[2,2]*m[3,3]-m[2,3]*m[3,2])+m[1,2]*(m[2,3]*m[3,1]-m[2,1]*m[3,3])+m[1,3]*(m[2,1]*m[3,2]-m[2,2]*m[3,1]); end
        if name===:Row1; return Vector3(m[1,1], m[1,2], m[1,3]); end
        if name===:Row2; return Vector3(m[2,1], m[2,2], m[2,3]); end
        if name===:Row3; return Vector3(m[3,1], m[3,2], m[3,3]); end
        if name===:Column1; return Vector3(m[1,1], m[2,1], m[3,1]); end
        if name===:Column2; return Vector3(m[1,2], m[2,2], m[3,2]); end
        if name===:Column3; return Vector3(m[1,3], m[2,3], m[3,3]); end
        return getfield(m, name)
    end

    function Base.propertynames(m::Matrix3)
        return (:A11, :A12, :A13, :A21, :A22, :A23, :A31, :A32, :A33,
        :Row1, :Row2, :Row3, :Column1, :Column2, :Column3, :Trace, :Determinat)
    end


text:
type GeneralBlackScholesType <: BlackScholesType end
type BlackScholesMertonType <: BlackScholesType end

type GeneralizedBlackScholesProcess{Y1 <: YieldTermStructure, Y2 <: YieldTermStructure, B <: BlackVolTermStructure, D <: AbstractDiscretization, BST <: BlackScholesType} <: AbstractBlackScholesProcess
  x0::Quote
  riskFreeRate::Y1
  dividendYield::Y2
  blackVolatility::B
  localVolatility::LocalConstantVol
  disc::D
  isStrikeDependent::Bool
  blackScholesType::BST
end

function GeneralizedBlackScholesProcess(x0::Quote, riskFreeRate::YieldTermStructure, dividendYield::YieldTermStructure,
          blackVolatility::BlackConstantVol, disc::AbstractDiscretization = EulerDiscretization())
  localVolatility = LocalConstantVol(blackVolatility.referenceDate, black_vol(blackVolatility, 0.0, x0.value), blackVolatility.dc)

  return GeneralizedBlackScholesProcess(x0, riskFreeRate, dividendYield, blackVolatility, localVolatility, disc, true, GeneralBlackScholesType())
end

# type BlackScholesMertonProcess{Y1 <: YieldTermStructure, Y2 <: YieldTermStructure, B <: BlackVolTermStructure, D <: AbstractDiscretization} <: AbstractBlackScholesProcess
#   x0::Quote
#   riskFreeRate::Y1
#   dividendYield::Y2
#   blackVolatility::B
#   localVolatility::LocalConstantVol
#   disc::D
#   isStrikeDependent::Bool
# end

function BlackScholesMertonProcess(x0::Quote, riskFreeRate::YieldTermStructure, dividendYield::YieldTermStructure,
          blackVolatility::BlackConstantVol, disc::AbstractDiscretization = EulerDiscretization())
  localVolatility = LocalConstantVol(blackVolatility.referenceDate, black_vol(blackVolatility, 0.0, x0.value), blackVolatility.dc)

  return GeneralizedBlackScholesProcess(x0, riskFreeRate, dividendYield, blackVolatility, localVolatility, disc, true, BlackScholesMertonType())
end

function drift(process::GeneralizedBlackScholesProcess, t::Float64, x::Float64)
  sigma = diffusion(process, t, x)
  t1 = t + 0.0001

  return forward_rate(process.riskFreeRate, t, t1, ContinuousCompounding(), NoFrequency()).rate -
          forward_rate(process.dividendYield, t, t1, ContinuousCompounding(), NoFrequency()).rate -
          0.5 * sigma * sigma
end

diffusion(process::GeneralizedBlackScholesProcess, t::Float64, x::Float64) = local_vol(process.localVolatility, t, x)

apply(process::GeneralizedBlackScholesProcess, x0::Float64, dx::Float64) = x0 * exp(dx)

expectation(process::GeneralizedBlackScholesProcess, ::Float64, ::Float64, ::Float64) = error("not implemented")

get_x0(process::AbstractBlackScholesProcess) = process.x0.value

function evolve(process::GeneralizedBlackScholesProcess, t0::Float64, x0::Float64, dt::Float64, dw::Float64)
  if process.isStrikeDependent
    var = black_variance(process.blackVolatility, t0 + dt, 0.01) - black_variance(process.blackVolatility, t0, 0.01)
    drift_ = (forward_rate(process.riskFreeRate, t0, t0 + dt, ContinuousCompounding(), NoFrequency()).rate -
              forward_rate(process.dividendYield, t0, t0 + dt, ContinuousCompounding(), NoFrequency()).rate) *
              dt - 0.5 * var

    return x0 * exp(sqrt(var) * dw + drift_)
  else
    return apply(process, x0, drift(process.disc, process, t0, x0, dt) + std_deviation(process, t0, x0, dt) * dw)
  end
end

state_variable(process::AbstractBlackScholesProcess) = process.x0

get_time(process::AbstractBlackScholesProcess, d::Date) = year_fraction(process.riskFreeRate.dc, reference_date(process.riskFreeRate), d)


text:
#=
    hilbert
    Copyright © 2015 jaidev <jaidev@newton>

    Distributed under terms of the MIT license.
=#

module hht

function hilbert(x)
    N = length(x)
    Xf = fft(x)
    h = zeros(N)
    h[1] = 1
    if N % 2 == 0
        h[floor(N / 2)] = 1
        h[1:floor(N / 2)] = 2
    else
        h[1:floor((N + 1) / 2)] = 2
    end
    x = ifft(Xf .* h)
    return x
end

end


text:
# KronigPenny

"""
    mutable struct KronigPennyPotential <: Potential

    KronigPennyPotential(v0, ρ) # Default constructor

holds parameters of Kronig-Penny potential.

Fields
* `v0` : potential height in units of ``E_{1}^{(0)}``
* `ρ` : barrier width in units of period ``a``,  ``\\rho = \\dfrac{b}{a}``, ``0 < \\rho < 1``
  * Note that a position ``x`` is expressed in units of ``a`` throughout this package.

The constructor `KronigPennyPotential(v0, ρ)` confirms that `0 ≤ ρ ≤ 1`, otherwise throws an error.
"""
mutable struct KronigPennyPotential <: Potential
   v0:: Real
   ρ:: Real # ρ = b/a

   KronigPennyPotential(v0, ρ) = 0 ≤ ρ ≤ 1 ? new(v0,ρ) : error("ρ < 0 or 1 < ρ")
end

"""
    get_potential(::KronigPennyPotential)

returns a function ``v`` to evaluate potential ``v(x)`` as a position ``x``,
such that:
```math
\\begin{aligned}
v(x) & = \\begin{cases}
v_{0} &
\\text{inside well, i.e.,}
 \\dfrac{1-\\rho}{2} \\le \\dfrac{x}{a} \\le \\dfrac{1+\\rho}{2},
 \\\\
0 & \\text{outside well}\\end{cases}
\\\\
v(x+a) &= v(x)
\\end{aligned}
```
* Note that a position ``x`` is expressed in units of ``a`` throughout this package.
"""
function get_potential(pot::KronigPennyPotential)
   return x -> (1-pot.ρ)/2 <= mod(x,1)< (1+pot.ρ)/2 ? 0 : pot.v0
end

"""
    mutable struct KronigPennyModel <: Model

holds parameters of Kronig-Penny model.

* Fields given to the constructor:
  * `potential`
  * `nmax`
  * `Ka` : wavenumber multiplied by `a`, period
* Fields calculated by constructor:
  * `mmax` : size of Hamiltonian matrix
  * `qnum` : iterator of quantum numbers
* Fields calculated by constructor and `update!` methods:
  * `hnm` : hamiltonian matrix
"""
mutable struct KronigPennyModel <: Model
   potential::KronigPennyPotential
   nmax::Int64  # maximum of quantum numbers
   mmax::Int64  # size of Hamiltonian
   Ka::Float64
   qnum::Alternates
   hnm::Matrix
end

"""
    KronigPennyModel(::KronigPennyPotential, nmax, Ka)

is a constructor of Kronig-Penny model, and defines other fields: `qnum`, `nmax`, and `hnm`
* Mandantory parameters:
  * `pot` : potential
  * `nmax` : maximum of quantum numbers
  * `Ka` : wavenumber multiplied by `a`, period
"""
function KronigPennyModel(pot::KronigPennyPotential, nmax, Ka)
   qnum=Alternates(nmax)
   mmax=length(qnum)
   hnm=zeros(Float64,(mmax,mmax))
   model_lem=KronigPennyModel(pot,nmax,mmax,Ka,qnum,hnm)
   constuctMatrix(model_lem)
   model_lem
end

"""
    update!(model_lem::KronigPennyModel; Ka=nothing)

updates Kronnig-Penny model.

* Optional parameters
  * `Ka` : wavenumber multiplied by `a`, period
"""
function update!(model_lem::KronigPennyModel; Ka=nothing)
   if !isnothing(Ka)
      model_lem.Ka=Ka
   end
   constuctMatrix(model_lem)
   model_lem
end

"""
    constuctMatrix(model_lem::KronigPennyModel)

computes and fills Hamiltonian matrix fields `hnm` in `model_lem`.

```math
h_{nm} = \\begin{cases}
\\left(2n + \\dfrac{Ka}{\\pi}\\right)^{2} + v_{0} (1-\\rho) &
\\text{for}\\; n = m\\;\\text{(diagonal elements)} \\\\
v_{0}
\\dfrac{(-1)^{m-n+1}}{\\pi} \\dfrac{\\sin \\pi(m-n)\\rho}{m-n}
 &
 \\text{for}\\; n \\neq m\\;\\text{(off-diagonal elements)}\\end{cases}
```
"""
function constuctMatrix(model_lem::KronigPennyModel)
   qnum=model_lem.qnum
   hnm=model_lem.hnm
   Ka=model_lem.Ka
   potential=model_lem.potential
   v0=potential.v0
   ρ=potential.ρ

   for (i,m) in enumerate(qnum)
      for (j,n) in Iterators.take(enumerate(qnum), i-1)
         hnm[i,j] = 0
      end
   end

   for (i,m) in enumerate(qnum)
      hnm[i,i] = (2m+Ka/π)^2 + v0*(1-ρ)
   end

   for (i,m) in enumerate(qnum)
      for (j,n) in Iterators.take(enumerate(qnum), i-1)
         hnm[i,j] = v0*(-1)^(m-n+1)/π * sinpi((m-n)ρ)/(m-n)
      end
      for (j,n) in Iterators.drop(enumerate(qnum), i)
         hnm[i,j] = v0*(-1)^(m-n+1)/π * sinpi((m-n)ρ)/(m-n)
      end
   end
end


text:
using Random

function get_mask(region::AbstractArray, scramble_rate)
    mask = zeros(1, 1, 1, size(region)[2:end]...)
    num_masked = Integer(floor(scramble_rate * size(region)[2] * size(region)[3]))
    idx = [(i,j) for i in 1:size(region)[2], j in 1:size(region)[3]]
    idx = shuffle(idx)[1:num_masked]

    for n in idx
        mask[1,1,1,n...] = -Inf
    end

    return mask
end


text:
# QuEST.jl/deps/build.jl
# MIT License

#
# This script builds the package.
#

using Clang
using Clang.LibClang.Clang_jll

"""
True, if expert installation has been selected
"""
const EXPERT = haskey(ENV, "QUEST_EXPERT") # && ENV["QUEST_EXPERT"] == "1" ? true : false

"""
Execute commands to build QuEST
"""
function _auxBuild(makePrecision ::Int, precision ::String ; isWindows ::Bool) ::Nothing
    mkdir("build"*precision)
    cd("build"*precision)

    # `run()` with `wait=true` throws an error if anything goes wrong,
    # e.g., non-zero exit status

    if isWindows
        run(`cmake -DPRECISION=$makePrecision .. -G "MinGW Makefiles"`; wait=true)
        libclang_include = joinpath("..","QuEST","include") |> normpath
        libclang_headers = [ joinpath(libclang_include,header)   for header in readdir(libclang_include) if endswith(header, ".h")]
        wc = init(; headers = libclang_headers,
                  output_file = joinpath(@__DIR__, "libclang_api.jl"),
                  common_file = joinpath(@__DIR__, "libclang_common.jl"),
                  clang_includes = vcat(libclang_include, CLANG_INCLUDE),
                  clang_args = ["-I", joinpath(libclang_include, ".."), "-DQuEST_PREC=$(makePrecision)"],
                  header_wrapped = (root, current)->root == current,
                  header_library = x->"libclang",
                  clang_diagnostics = true,
                  )
        run(wc)

    else
        run(`cmake -DPRECISION=$makePrecision ..` ; wait=true)
        libclang_include = joinpath("..","QuEST","include") |> normpath
        libclang_headers = [ joinpath(libclang_include,header)   for header in readdir(libclang_include) if endswith(header, ".h")]
        wc = init(; headers = libclang_headers,
                  output_file = joinpath(@__DIR__, "libclang_api.jl"),
                  common_file = joinpath(@__DIR__, "libclang_common.jl"),
                  clang_includes = vcat(libclang_include, CLANG_INCLUDE),
                  clang_args = ["-I", joinpath(libclang_include, ".."), "-DQuEST_PREC=$(makePrecision)"],
                  header_wrapped = (root, current)->root == current,
                  header_library = x->"libclang",
                  clang_diagnostics = true,
                  )
        run(wc)
    end

    run(`make` ; wait=true)

    cd("..")
end

"""
Clone repository, make, Clang
"""
function build(;isWindows::Bool) ::Nothing
    @info "Cloning QuEST..."
    run(`git clone https://github.com/Ketita/QuEST.git` ; wait=true)
    cd("QuEST")
    @info "Build with 32-bit floats..."
    _auxBuild(1,"32" ; isWindows)
    @info "Build with 64-bit floats..."
    _auxBuild(2,"64" ; isWindows)
    @info "Builds successful."
end

#
# File "main" section
#

if EXPERT
    @warn "Expert installation: Make sure the libraries `libQuEST_32` and `libQuEST_64` are loadable."
else
    @warn "Non-expert installation: Downloading QuEST and building it with default settings"

    if ispath("QuEST") && !isempty(readdir("QuEST"))
        rm("QuEST" ; force=true,recursive=true)
    end

    if Sys.isunix()
        build(isWindows=false)
    elseif Sys.iswindows()
        build(isWindows=true)
    else
        error("OS not supported.")
    end
end #^ else

#EOF


text:
using Test

# Symmetric monoidal category
#############################

A, B = Ob(FreeSymmetricMonoidalCategory, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)

# Domains and codomains
@test dom(otimes(f,g)) == otimes(dom(f),dom(g))
@test codom(otimes(f,g)) == otimes(codom(f),codom(g))
@test dom(braid(A,B)) == otimes(A,B)
@test codom(braid(A,B)) == otimes(B,A)
@test σ(A, B) == braid(A,B)

# Associativity and unit
I = munit(FreeSymmetricMonoidalCategory.Ob)
@test otimes(A,I) == A
@test otimes(I,A) == A
@test otimes(otimes(A,B),A) == otimes(A,otimes(B,A))
@test otimes(otimes(f,g),f) == otimes(f,otimes(g,f))

# Extra syntax
@test otimes(A,B,A) == otimes(otimes(A,B),A)
@test otimes([A,B,A]) == otimes(otimes(A,B),A)
@test otimes(f,f,f) == otimes(otimes(f,f),f)
@test otimes([f,f,f]) == otimes(otimes(f,f),f)
@test otimes(FreeSymmetricMonoidalCategory.Ob[]) == I
@test_throws MethodError otimes([])
@test A⊗B == otimes(A,B)
@test f⊗g == otimes(f,g)
@test A⊗B⊗A == otimes(A,B,A)

# Extra functions
@test collect(A) == [A]
@test collect(otimes(A,B)) == [A,B]
@test collect(I) == []
@test typeof(collect(I)) == Vector{FreeSymmetricMonoidalCategory.Ob}
@test ndims(A) == 1
@test ndims(otimes(A,B)) == 2
@test ndims(I) == 0

# String format
@test string(otimes(f,g)) == "otimes(f,g)"
@test string(compose(otimes(f,f),otimes(g,g))) == "compose(otimes(f,f),otimes(g,g))"

# S-expressions
@test sexpr(I) == "(munit)"
@test sexpr(otimes(A,B)) == "(otimes :A :B)"
@test sexpr(otimes(f,g)) == "(otimes :f :g)"
@test sexpr(compose(otimes(f,f),otimes(g,g))) == "(compose (otimes :f :f) (otimes :g :g))"

# Infix notation (Unicode)
@test unicode(I) == "I"
@test unicode(otimes(A,B)) == "A⊗B"
@test unicode(otimes(f,g)) == "f⊗g"
@test unicode(compose(otimes(f,f),otimes(g,g))) == "(f⊗f)⋅(g⊗g)"
@test unicode(otimes(compose(f,g),compose(g,f))) == "(f⋅g)⊗(g⋅f)"

# Infix notation (LaTeX)
@test latex(I) == "I"
@test latex(otimes(A,B)) == "A \\otimes B"
@test latex(otimes(f,g)) == "f \\otimes g"
@test latex(compose(otimes(f,f),otimes(g,g))) ==
  "\\left(f \\otimes f\\right) \\cdot \\left(g \\otimes g\\right)"
@test latex(otimes(compose(f,g),compose(g,f))) ==
  "\\left(f \\cdot g\\right) \\otimes \\left(g \\cdot f\\right)"
@test latex(braid(A,B)) == "\\sigma_{A,B}"
@test latex(σ(A,B)) == "\\sigma_{A,B}"

# Cartesian category
####################

A, B = Ob(FreeCartesianCategory, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)

# Domains and codomains
@test dom(mcopy(A)) == A
@test codom(mcopy(A)) == otimes(A,A)
@test dom(delete(A)) == A
@test codom(delete(A)) == I

# Derived syntax
@test pair(f,f) == compose(mcopy(A), otimes(f,f))
@test proj1(A,B) == otimes(id(A), delete(B))
@test proj2(A,B) == otimes(delete(A), id(B))

# LaTeX notation
@test latex(mcopy(A)) == "\\Delta_{A}"
@test latex(delete(A)) == "\\lozenge_{A}"

# Biproduct category
####################

A, B = Ob(FreeBiproductCategory, :A, :B)

# LaTeX notation
@test latex(mmerge(A)) == "\\nabla_{A}"
@test latex(create(A)) == "\\square_{A}"

# Cartesian closed category
###########################

A, B, C = Ob(FreeCartesianClosedCategory, :A, :B, :C)
f = Hom(:f, otimes(A,B), C)

# Domains and codomains
@test dom(ev(A,B)) == otimes(hom(A,B),A)
@test codom(ev(A,B)) == B
@test dom(curry(A,B,f)) == A
@test codom(curry(A,B,f)) == hom(B,C)

# LaTeX notation
@test latex(hom(A,B)) == "{B}^{A}"
@test latex(hom(otimes(A,B),C)) == "{C}^{A \\otimes B}"
@test latex(hom(A,otimes(B,C))) == "{\\left(B \\otimes C\\right)}^{A}"
@test latex(ev(A,B)) == "\\mathrm{eval}_{A,B}"
@test latex(curry(A,B,f)) == "\\lambda f"

# Compact closed category
#########################

A, B, C = Ob(FreeCompactClosedCategory, :A, :B, :C)
I = munit(FreeCompactClosedCategory.Ob)
f, g = Hom(:f, A, B), Hom(:g, B, A)

# Duals
@test dual(otimes(A,B)) == otimes(dual(B),dual(A))
@test dual(I) == I
@test dual(dual(A)) == A
@test dual(otimes(dual(A),dual(B))) == otimes(B,A)
@test hom(A,B) == otimes(B,dual(A))

# Mates
@test mate(mate(f)) == f
@test mate(compose(f,g)) == compose(mate(g),mate(f))
@test mate(otimes(f,g)) == otimes(mate(g),mate(f))

# Domains and codomains
@test dom(dunit(A)) == I
@test codom(dunit(A)) == otimes(dual(A), A)
@test dom(dcounit(A)) == otimes(A, dual(A))
@test codom(dcounit(A)) == I
@test dom(mate(f)) == dual(B)
@test codom(mate(f)) == dual(A)

f = Hom(:f, otimes(A,B), C)
@test dom(ev(A,B)) == otimes(B,dual(A),A)
@test codom(ev(A,B)) == B
@test dom(curry(A,B,f)) == A
@test codom(curry(A,B,f)) == otimes(C,dual(B))

# LaTeX notation
@test latex(dual(A)) == "{A}^*"
@test latex(dunit(A)) == "\\eta_{A}"
@test latex(dcounit(A)) == "\\varepsilon_{A}"
@test latex(mate(f)) == "{f}^*"

# Dagger category
#################

A, B = Ob(FreeDaggerCategory, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)

# Domains and codomains
@test dom(dagger(f)) == B
@test codom(dagger(f)) == A

# Dagger
@test dagger(compose(f,g)) == compose(dagger(g),dagger(f))
@test dagger(id(A)) == id(A)
@test dagger(dagger(f)) == f

# LaTeX notation
@test latex(dagger(f)) == "{f}^\\dagger"
#@test latex(dagger(compose(f,g))) == "{\\left(f \\cdot g\\right)}^\\dagger"

# Dagger compact category
#########################

A, B = Ob(FreeDaggerCompactCategory, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)

# Dagger
@test dagger(compose(f,g)) == compose(dagger(g),dagger(f))
@test dagger(id(A)) == id(A)
@test dagger(dagger(f)) == f
@test dagger(otimes(f,g)) == otimes(dagger(f),dagger(g))
@test (dagger(otimes(compose(f,g),compose(g,f))) ==
       otimes(compose(dagger(g),dagger(f)),compose(dagger(f),dagger(g))))

# Traced monoidal category
##########################

A, B, X, Y = Ob(FreeTracedMonoidalCategory, :A, :B, :X, :Y)
f = Hom(:f, X⊗A, X⊗B)

# Domains and codomains
@test dom(trace(X, A, B, f)) == A
@test codom(trace(X, A, B, f)) == B
#@test_throws SyntaxDomainError trace(X, A, B, Hom(:g, A⊗X, B⊗X))

# LaTeX notation
@test latex(trace(X, A, B, f)) == "\\operatorname{Tr}_{A,B}^{X} \\left(f\\right)"


text:
module CLIMAParameters

export AbstractParameterSet
export AbstractEarthParameterSet

include("file_parsing.jl")
include("types.jl")
include("UniversalConstants.jl")

Base.broadcastable(aps::AbstractParameterSet) = Ref(aps)

# Function stubs
include(joinpath("Planet", "Planet.jl"))
include(joinpath("Atmos", "Atmos.jl"))
include(joinpath("SurfaceFluxes", "SurfaceFluxes.jl"))
include(joinpath("SubgridScale", "SubgridScale.jl"))

# Define values
include(joinpath("Planet", "planet_parameters.jl"))
include(joinpath("Atmos", "atmos_parameters.jl"))
include(joinpath("SurfaceFluxes", "surface_fluxes_parameters.jl"))
include(joinpath("SubgridScale", "subgrid_scale_parameters.jl"))

end # module

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
# Extend broadcasting for Sequence

struct SequenceStyle <: Base.Broadcast.BroadcastStyle end

Base.Broadcast.broadcastable(a::Sequence) = a

Base.Broadcast.BroadcastStyle(::Type{<:Sequence}) = SequenceStyle()
Base.Broadcast.BroadcastStyle(::Base.Broadcast.AbstractArrayStyle, ::SequenceStyle) = SequenceStyle()
Base.Broadcast.BroadcastStyle(::SequenceStyle, ::Base.Broadcast.AbstractArrayStyle) = SequenceStyle()

function Base.similar(bc::Base.Broadcast.Broadcasted{SequenceStyle}, ::Type{ElType}) where {ElType}
    space = _find_space(bc)
    return Sequence(space, similar(Vector{ElType}, dimension(space)))
end

_find_space(bc::Base.Broadcast.Broadcasted) = _find_space(bc.args)
_find_space(args::Tuple) = _find_space(_find_space(args[1]), _find_space(Base.tail(args)))
_find_space(a::Sequence) = a.space
_find_space(::Tuple{}) = nothing
_find_space(::Any) = nothing
_find_space(space::VectorSpace, ::Nothing) = space
_find_space(::Nothing, space::VectorSpace) = space
_find_space(::Nothing, ::Nothing) = nothing
function _find_space(s₁::T, s₂::S) where {T<:VectorSpace,S<:VectorSpace}
    s₁ == s₂ || return throw(DomainError)
    return convert(promote_type(T, S), s₁)
end

function Base.copyto!(dest::Sequence, bc::Base.Broadcast.Broadcasted)
    axes(dest) == axes(bc) || return throw(DimensionMismatch)
    bc′ = Base.Broadcast.preprocess(dest, bc)
    @inbounds @simd for i ∈ eachindex(bc′)
        dest.coefficients[i] = bc′[i]
    end
    return dest
end

Base.@propagate_inbounds Base.Broadcast._broadcast_getindex(a::Sequence, i::Int) =
    Base.Broadcast._broadcast_getindex(a.coefficients, i)

Base.@propagate_inbounds Base.Broadcast._broadcast_getindex(a::Sequence, I::CartesianIndex) =
    Base.Broadcast._broadcast_getindex(a.coefficients, I)

# to allow a[...] .= f.(...)
Base.@propagate_inbounds Base.Broadcast.dotview(a::Sequence, α) = view(a, α)

# Extend broadcasting for LinearOperator

struct LinearOperatorStyle <: Base.Broadcast.BroadcastStyle end

Base.Broadcast.broadcastable(A::LinearOperator) = A

Base.Broadcast.BroadcastStyle(::Type{<:LinearOperator}) = LinearOperatorStyle()
Base.Broadcast.BroadcastStyle(::Base.Broadcast.AbstractArrayStyle, ::LinearOperatorStyle) = LinearOperatorStyle()
Base.Broadcast.BroadcastStyle(::LinearOperatorStyle, ::Base.Broadcast.AbstractArrayStyle) = LinearOperatorStyle()

function Base.similar(bc::Base.Broadcast.Broadcasted{LinearOperatorStyle}, ::Type{ElType}) where {ElType}
    domain, codomain = _find_domain_codomain(bc)
    return LinearOperator(domain, codomain, similar(Matrix{ElType}, dimension(codomain), dimension(domain)))
end

_find_domain_codomain(bc::Base.Broadcast.Broadcasted) = _find_domain_codomain(bc.args)
_find_domain_codomain(args::Tuple) = _find_domain_codomain(_find_domain_codomain(args[1]), _find_domain_codomain(Base.tail(args)))
_find_domain_codomain(A::LinearOperator) = (A.domain, A.codomain)
_find_domain_codomain(::Tuple{}) = nothing
_find_domain_codomain(::Any) = nothing
_find_domain_codomain(domain_codomain::NTuple{2,VectorSpace}, ::Nothing) = domain_codomain
_find_domain_codomain(::Nothing, domain_codomain::NTuple{2,VectorSpace}) = domain_codomain
_find_domain_codomain(::Nothing, ::Nothing) = nothing
function _find_domain_codomain(s₁::Tuple{T₁,S₁}, s₂::Tuple{T₂,S₂}) where {T₁<:VectorSpace,S₁<:VectorSpace,T₂<:VectorSpace,S₂<:VectorSpace}
    s₁ == s₂ || return throw(DomainError)
    return @inbounds (convert(promote_type(T₁, T₂), s₁[1]), convert(promote_type(S₁, S₂), s₁[2]))
end

function Base.copyto!(dest::LinearOperator, bc::Base.Broadcast.Broadcasted)
    axes(dest) == axes(bc) || return throw(DimensionMismatch)
    bc′ = Base.Broadcast.preprocess(dest, bc)
    @inbounds @simd for i ∈ eachindex(bc′)
        dest.coefficients[i] = bc′[i]
    end
    return dest
end

Base.@propagate_inbounds Base.Broadcast._broadcast_getindex(A::LinearOperator, i::Int) =
    Base.Broadcast._broadcast_getindex(A.coefficients, i)

Base.@propagate_inbounds Base.Broadcast._broadcast_getindex(A::LinearOperator, I::CartesianIndex) =
    Base.Broadcast._broadcast_getindex(A.coefficients, I)

# to allow A[...] .= f.(...)
Base.@propagate_inbounds Base.Broadcast.dotview(A::LinearOperator, α, β) = view(A, α, β)


text:
apiUrl = get(ENV,"API_URL","https://api.navability.io")
userId = get(ENV,"USER_ID","guest@navability.io")
robotId = get(ENV,"ROBOT_ID","IntegrationRobot")
sessionId = get(ENV,"SESSION_ID","TestSession_"*string(uuid4())[1:8])

@testset "nva-sdk-standard-api-testset" begin

    client = NavAbilityHttpsClient(apiUrl)
    context = Client(userId,robotId,sessionId)

    addVariable(client, context,  "x0", :Pose2)
    addFactor(client, context, ["x0"], NVA.PriorPose2(Z=FullNormal([0.,1.,0], diagm([1.,1,1]))))

    addVariable(client, context,  :x1, :Pose2)
    addFactor(client, context, [:x0,:x1], NVA.Pose2Pose2(Z=FullNormal([1.,0.,0], diagm([1.,1,1]))))

    addFactor(client, context, [:x1], NVA.PriorPose2(Z=FullNormal([1.,0.,0], diagm([1.,1,1]))); nullhypo=0.1)

    addFactor(client, context, ["x0"], NVA.PriorPose2(Z=FullNormal([0.5,0.5,0], diagm([1.,1,1]))))

    addVariable(client, context,  :x2_a, :Pose2)
    addVariable(client, context,  :x2_b, :Pose2)
    addFactor(client, context, ["x1", "x2_a", "x2_b"], NVA.PriorPose2(Z=FullNormal([0.5,0.5,0], diagm([1.,1,1]))); multihypo=[1,0.1,0.9])

    addVariable(client, context,  "y0", :Position1)
    addFactor(client, context, ["y0"], NVA.Prior(Z=Normal(0.0, 1.0)))

    addVariable(client, context,  "z0", :Point2)
    addFactor(client, context, ["z0"], NVA.PriorPoint2(Z=FullNormal([1.,0.], diagm([1.,1]))))

    addVariable(client, context,  "z1", :Point2)
    addFactor(client, context, ["z0", "z1"], NVA.Point2Point2(Z=FullNormal([1.,0.], diagm([1.,1]))))

end


text:
include("clonetask.jl")
include("brokentask.jl")
include("tarray.jl")
include("tarray2.jl")


text:
export findw
"""
function findw(param,M;maxIterQP,out,pwr)

Find a frequency vector w in Omega_W that solves the constrained linear least
squares problem

min_w | M*w - d |_{Gamma,pwr}

Input:

	1. param     - StrainReconParam
	2. M         - binary matrix
    3. maxIterQP - maximum number of iterations for asqp
	4. out       - show intermediate results for asqp/iplp if > 1
	5. pwr       - exponent in the likelihood function

Output:

	1. wsol      - solution vector for constrained least squares problem

"""
function findw(param::StrainReconParam,M::BitArray{2};maxIterQP::Int=5,out=-1,pwr::Float64=param.pwr)

	(m,n) = size(M)
	R     = param.R
	d  = param.d

	# equality constraint
	(Ci,ci,Ce,ce) = getConstraints(param)

	if pwr==2.0
		# solve equality constrained QP problem and see if constraints are satisfied
		Q   = M'*((R'*R)*M) # diagm(collect(0:n-1));
		c   = M'*(R'*(R*d))
		x0  = 1./collect(1:n)
		x0  ./= sum(x0)
		wsol, iter, Wc = asqp(Q,-c,Ce,ce,Ci,ci,x0,maxIter=maxIterQP,out=out)
	elseif pwr==1.0
        # solve equality constrained linear programming with interior point method
        # return wsol, Lagrange multipliers, slack variables, history of iterations
        w0  = 1./collect(1:n)
		w0  ./= sum(w0)
		# set up linear system to be solved by IPLP
	    ni = size(Ci,1)
	    ne = size(Ce,1)
	    Ae = [M -eye(m) eye(m) zeros(m,ni);
	          Ce zeros(ne,2*m) zeros(ne,ni);
	          Ci zeros(ni,2*m) -eye(ni)]
	    ae = [vec(d); vec(ce); vec(ci); ]

	    c  = [zeros(n); diag(R); diag(R); zeros(ni)]
	    # how to initialize Lagrange mutlipliers and slack variables?
		x0 = rand(n+2*m+ni)
		y0 = rand(size(Ae,1))
	    s0 = rand(size(x0))

        wsol, lm, s, his = iplpPathFollowing(Ae,ae,c,x0,y0,s0,maxIter=100,out=-1,ftol=1e-12,mutol=1e-12)
		wsol = wsol[1:n]
	end

	return wsol
end


text:
function rang(A; isStep = false)::Int
    if isStep == false
        convert!(A)
    end
    rows, _ = size(A)
    for i in 1:rows
        if A[i, i] == 0
            return i-1
        end
    end
    return rows
end


text:
"""
A Julia suite for chaos and nonlinear dynamics
"""
module DynamicalSystems

using StaticArrays

export SVector, SMatrix, @SVector, @SMatrix, Size

using Reexport

@reexport using DelayEmbeddings
@reexport using DynamicalSystemsBase
@reexport using ChaosTools
@reexport using RecurrenceAnalysis

display_update = true
update_name = "update_v1.1.0"

if display_update
if !isfile(joinpath(@__DIR__, update_name))
printstyled(stdout,
"""
\nUpdate message: DynamicalSystems v1.1

A new package has joined DynamicalSystems: RecurrenceAnalysis !

This new package offers tools to compute and analyze
recurrences in your timeseries!

Check out the documentation for more!\n
"""; color = :light_magenta)
touch(joinpath(@__DIR__, update_name))
end
end

end # module


text:
using GraphUtilities, Graphs, GenericTensorNetworks
using GraphUtilities: save_property, load_property, SmallGraphConfig
using Test

@testset "file io" begin
    # creating file
    graph = SmallGraphConfig("petersen")
    config = GraphProblemConfig(problem="IndependentSet", graph=graph)
    folder = GraphUtilities.foldername("data", config; create=true)
    @test isdir(folder)
    @test isfile(joinpath(folder, "info.toml"))
    graph = SmallGraphConfig("petersen")
    config = GraphProblemConfig(problem="IndependentSet", graph=graph)
    folder2 = GraphUtilities.foldername("data", config; create=false)
    @test folder2 == folder

    # save load problem code
    gp = problem_instance(config)
    GraphUtilities.save_code(folder, gp)
    gp2 = GraphUtilities.load_code(config, folder)
    for field in fieldnames(typeof(gp))
        @test getfield(gp, field) == getfield(gp2, field)
    end

    # save load configs (table)
    prop = ConfigsMax(; tree_storage=false)
    res = solve(gp, prop)[]
    fc = joinpath(folder, "table")
    !isdir(fc) && mkdir(fc)
    @show fc
    GraphUtilities.saveconfigs(fc, [res.n], [res.c])
    fname = joinpath(fc, "size_4.0.dat")
    @test isfile(fname)
    configs = GraphUtilities.loadconfigs(fc, [res.n]; tree_storage=false, bitlength=10)
    @test configs == [res.c]

    # TREE STORAGE
    # creating file
    graph = smallgraph(:petersen)
    gp = IndependentSet(graph)
    prop = ConfigsMax(; tree_storage=true)

    # save load configs (tree)
    res = solve(gp, prop)[]
    fc = joinpath(folder, "tree")
    !isdir(fc) && mkdir(fc)
    GraphUtilities.saveconfigs(fc, [res.n], [res.c])
    fname = joinpath(fc, "size_4.0.dat")
    @test isfile(fname)
    configs = GraphUtilities.loadconfigs(fc, [res.n]; tree_storage=true, bitlength=10)
    @test configs == [res.c]
end

@testset "save property" begin
    # create a folder
    graph = SmallGraphConfig("petersen")
    config = GraphProblemConfig(; problem="IndependentSet", graph=graph)
    folder = GraphUtilities.foldername("data", config; create=true)

    # instantiate a graph problem
    gp = problem_instance(config)
    GraphUtilities.save_code(folder, gp)
    gp2 = GraphUtilities.load_code(config, folder)

    for property in [SizeMax(), SizeMin(), SizeMax(3), SizeMin(3), CountingAll(),
        CountingMax(), CountingMin(), CountingMax(3), CountingMin(3), GraphPolynomial(),
        SingleConfigMax(; bounded=false), SingleConfigMin(; bounded=false), SingleConfigMax(3; bounded=false), SingleConfigMin(3; bounded=false),
        SingleConfigMax(; bounded=true), SingleConfigMin(; bounded=true), SingleConfigMax(3; bounded=true), SingleConfigMin(3; bounded=true),
        ConfigsMax(; bounded=true, tree_storage=true), ConfigsMin(; bounded=true, tree_storage=false),
        ConfigsMax(; bounded=false, tree_storage=true), ConfigsMin(; bounded=false, tree_storage=false),
        ConfigsMax(3; bounded=true, tree_storage=true), ConfigsMin(3; bounded=true, tree_storage=false),
        ConfigsMax(3; bounded=false, tree_storage=true), ConfigsMin(3; bounded=false, tree_storage=false),
        ConfigsAll{true}(), ConfigsAll{false}()
        ]
        @show property
        res = solve(gp2, property)[]
        save_property(folder, property, res)
        @test load_property(folder, property) == res
    end
end


text:
# ### m10.10t.jl

using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking
Turing.setprogress!(false)

delim=';'
df = CSV.read(sr_datadir("Kline.csv"), DataFrame; delim);

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

# New col contact_high, set binary values 1/0 if high/low contact
df[!, :contact_high] = map(x -> ifelse(x=="high", 1, 0), df[:, :contact]);

# This is supposed to be a "bad" model since we take non-centered data for the
# predictor log_pop
@model ppl10_10(total_tools, log_pop, contact_high) = begin
    α ~ Normal(0, 100)
    βp ~ Normal(0, 1)
    βc ~ Normal(0, 1)
    βpc ~ Normal(0, 1)

    for i ∈ 1:length(total_tools)
        λ = exp(α + βp*log_pop[i] + βc*contact_high[i] +
            βpc*contact_high[i]*log_pop[i])
        total_tools[i] ~ Poisson(λ)
    end
end;

m10_10t = ppl10_10(df.total_tools, df.log_pop, df.contact_high)
nchains = 4; sampler = NUTS(0.65); nsamples=2000
chns10_10t = mapreduce(c -> sample(m10_10t, sampler, nsamples), chainscat, 1:nchains)

# Rethinking result

m_10_10s_result = "
     mean   sd  5.5% 94.5% n_eff Rhat
 a    0.94 0.37  0.36  1.53  3379    1
 bp   0.26 0.04  0.21  0.32  3347    1
 bc  -0.08 0.84 -1.41  1.23  2950    1
 bpc  0.04 0.09 -0.10  0.19  2907    1
";

# End of m10.10t.jl


text:
using TerminalUI
import TerminalUI: ListWidget

listw = Border(ListWidget([1:100;]),"List")
tty = Base.Terminals.TTYTerminal("xterm",STDIN,STDOUT,STDERR)
wait(FullScreenDialog(listw,tty))


text:
#!/bin/env julia

data = readstring("../data/saved_txt/chr1_train.txt")
chars = unique(data)
data_size, vocab_size = length(data), length(vocab_size)
print("data has $data_size characters, %vocab_size unique.")
char_to_ix = Dict{Char,Int64}(ch=>i for (i,ch) in enumerate(chars))
ix_to_char = map(reverse, collect(char_to_ix)) |> Dict

# hyperparameters
hidden_size = 100
seq_length = 25
learning_rate = 1e-1

# model parameters
Wxh = randn(hidden_size, vocab_size)*0.01
Whh = randn(hidden_size, hidden_size)*0.01
Why = randn(vocab_size,  hidden_size)*0.01
bh  = zeros(hidden_size, 1)
by  = zeros(vocab_size, 1)

@doc """
  inputs,targets are both list of integers.
  hprev is Hx1 array of initial hidden state
  returns the loss, gradients on model parameters, and last hidden state
""" ->
function lossFun(inputs::Array{Int32,1}, targets::Array{Int32,1}, hprev::Array{Int32,1})
    xs,hs,ys,ps = Dict(),Dict(),Dict(),Dict()
    hs[0] = copy(hprev)
    loss = 0.0
    # Forward pass
    for t in 1:length(inputs)
        xs[t] = zeros(vocab_size, 1) # opt
        xs[t][inputs[t]] = 1 # input vector is xs[t] with one-hot encoding
        hs[t] = tanh(Wxh*xs[t] + Whh*hs[t-1] + bh)
        ys[t] = Why*hs[t] + by
        ps[t] = exp(ys[t]) / sum(ys[t])
        loss += -log(ps[t][targets[t]])
    end
    # Backward pass
    dWxh, dWhh, dWhy = zeros(size(Wxh)),zeros(size(Whh)),zeros(size(Why))
    dbh, dby = zeros(size(bh)), zeros(size(by))
    dhnext = zeros(size(hs[1]))
    for t in length(inputs):-1:1
        dy = copy(ps[t])
        dy[targets[t]] -= 1
        dWhy += dy*hs[t]'
        dby += dy
        dh = Why'*dy + dhnext
        dhraw = (1 - hs[t]*hs[t]) * dh
        dbh  += dhraw
        dWxh += dhraw * xs[t]'
        dWhh += dhraw * hs[t-1]'
        dhnext = Whh' * dhraw
    end
    for dparam in [dWxh, dWhh, dWhy, dbh, dby]
        dparam = clip(dparam)
    end
    loss, dWxh, dWhh, dWhy, dbh, dby, hs[length(inputs)]
end

function clip(d::Array{Float32,1}, l::Float32, h::Float32)
    map(x->clip(x, l, h), d)
end
function clip(d::Float32, l::Float32, h::Float32)
    if d < l
        return l
    elseif d > h
        return h
    else
        return d
    end
end

function sample{T<:AbstractFloat}(p::Array{T,1})
    θ = rand()
    c = 0.0
    for (i, v) in enumerate(p)
        c += v
        c >= θ && break
    end
    i
end

@doc """
  sample a sequence of integers from the model
  h is memory state, seed_ix is seed letter for first time step
""" ->
function sample(h::Array{Float32,1}, seed_ix, n)
    ixes = Int32[]
    x = zeros(vocab_size,1)
    x[seed_ix] = 1
    for t in 1:n
        h = tanh(Wxh*x + Whh*h + bh)
        y = Why*y + by
        p = exp(y) ./ sum(exp(y))
        ix = sample(p)
        x = zeros(vocab_size,1) # OPT
        x[ix] = 1
        push!(ixes, ix)
    end
    ixes
end

n, p = 0, 0
mWxh, mWhh, mWhy = zeros(size(Wxh)), zeros(size(Whh)), zeros(size(Why)) # ???
mbh, mby = zeros(size(bh)), zeros(size(by))
smooth_loss = - log(1.0/vocab_size) * seq_length # a vector or a interger ???

while true
    if p+seq_length+1 >= length(data) || n == 0
        hprev = zeros(hidden_size,1)
        p = 0
    end
    inputs = map(char->char_to_ix[char], data[p:p+seq_length+1])
    targets = string(inputs[2:end], data[p+seq_length+2])

    if n % 100 == 0
        sample_ix = sample(hprev, inputs[0], 200)
        txt = join(map(ix->ix_to_char[ix], sample_ix))
        println("----\n $txt \n----")
    end

    loss, dWxh, dWhy, dbh, dby, hprev = lossFun(inputs, targets, hprev)
    if n % 100 == 0
        println("iter: %n, loss: $smooth_loss")
    end

    for (param, dparam, mem) in zip([Wxh, Whh, Why, bh, by],
                                    [dWxh, dWhy, dbh, dby],
                                    [mWxh, mWhh, mWhy, mby])
        mem += dparam * dparam
        param += -learning_rate * dparam / sqrt(mem + 1e-8)
    end

    p += seq_length
    n += 1
end


text:
"""
    Filter

A wrapper around a function that takes a log `Record` and returns
a bool whether to skip logging it.

# Fields
`f::Function`: a function that should return a bool given a `Record`
"""
struct Filter
    f::Function
end

function (filter::Filter)(rec::Record)::Bool
    return filter.f(rec)
end


text:
module Layout

using Bukdu

function layout(title, script, style, body)
    render(HTML, """
<html>
<head>
  <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
  <title>Bukdu $title</title>
  $script
  $style
</head>
<body>
$body
</body>
</html>
""")
end

end # module Layout

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
module RunTests

include("constraint_tests.jl")
include("dominant_set_tests.jl")
include("automaton_tests.jl")
include("sequence_tests.jl")
include("printing_tests.jl")

end # module


text:
#!/usr/bin/env julia

using DataFrames

import Celeste.AccuracyBenchmark
import Celeste.ArgumentParse

parser = ArgumentParse.ArgumentParser()
ArgumentParse.add_argument(
    parser,
    "--write-error-csv",
    help="Write raw errors to the given CSV file",
)
ArgumentParse.add_argument(
    parser,
    "ground_truth_csv",
    help="Ground truth CSV catalog",
)
ArgumentParse.add_argument(
    parser,
    "celeste_prediction_csv",
    help="Celeste predictions CSV catalog",
)
parsed_args = ArgumentParse.parse_args(parser, ARGS)

truth = AccuracyBenchmark.read_catalog(parsed_args["ground_truth_csv"])
predictions = AccuracyBenchmark.read_catalog(parsed_args["celeste_prediction_csv"])
uncertainty_df = AccuracyBenchmark.get_uncertainty_df(truth, predictions)
scores = AccuracyBenchmark.score_uncertainty(uncertainty_df)
println(repr(scores))

if haskey(parsed_args, "write-error-csv")
    writetable(parsed_args["write-error-csv"], uncertainty_df)
end


text:
function matmul!(A, B, R)
    for i = 1:size(A)[1], j = 1:size(B)[2], k = 1:size(B)[1]
        R[i, j] += A[i, k] * B[k, j]
    end
end

max_stars_repo_name: mjaglarz/BA-thesis

max_issues_repo_name: mjaglarz/BA-thesis

max_forks_repo_name: mjaglarz/BA-thesis


text:
#!/usr/bin/env julia

using Distributed
addprocs(${ncpus})

@everywhere using FlashWeave

using CSV
using DataFrames
metadata_path = "fmt-control-Family_sample_metadata.tsv"

# data
data_path = "${otu_file}"
new_data_path = "${otu_file.baseName}.tsv"
if data_path != new_data_path
    cp(data_path, new_data_path, force=true, follow_symlinks=true)
end
metadata_path = "${sample_metadata}"
# NOTE: Transposing directly here does not work if the column names are long
metadata = CSV.read(metadata_path, DataFrame)
metadata_transpose = DataFrame(
    [[names(metadata)]; collect.(eachrow(metadata))], [:column; Symbol.(axes(metadata, 1))]
)
new_metadata_path = "${sample_metadata.baseName}_transposed.tsv"
CSV.write(new_metadata_path, metadata_transpose; writeheader=false, delim="\\t")

# options
sensitive = ${sensitive}
heterogeneous = ${heterogeneous}
FDR = ${fdr_correction}

netw_results = learn_network(
new_data_path,
new_metadata_path,
sensitive=sensitive,
heterogeneous=heterogeneous,
FDR=FDR,
transposed=true,
)

# You need to use networkx to export the .gml to table
output_file = "${meta.id}_network.gml"
save_network(output_file, netw_results)

max_stars_repo_name: segrelab/MiCoNE-synthetic-data

max_issues_repo_name: segrelab/MiCoNE-synthetic-data

max_forks_repo_name: dileep-kishore/MiCoNE


text:
# Tuple equality tests
@test (0.0,1.0) ≈ (0.00000000001,0.99999999)
@test !((0.0,1.0) ≈ (0.00001,0.9999))
@test [(0.0,1.0),(2.00000001,2.49999999)] ≈ [(0.00000001,0.999999999),(2.0,2.5)]
@test !([(0.0,1.0),(2.0,3.0)] ≈ [(0.0000001,0.999999)])


text:
inpt = """turn off 660,55 through 986,197
turn off 341,304 through 638,850
turn off 199,133 through 461,193
toggle 322,558 through 977,958
toggle 537,781 through 687,941
turn on 226,196 through 599,390
turn on 240,129 through 703,297
turn on 317,329 through 451,798
turn on 957,736 through 977,890
turn on 263,530 through 559,664
turn on 158,270 through 243,802
toggle 223,39 through 454,511
toggle 544,218 through 979,872
turn on 313,306 through 363,621
toggle 173,401 through 496,407
toggle 333,60 through 748,159
turn off 87,577 through 484,608
turn on 809,648 through 826,999
toggle 352,432 through 628,550
turn off 197,408 through 579,569
turn off 1,629 through 802,633
turn off 61,44 through 567,111
toggle 880,25 through 903,973
turn on 347,123 through 864,746
toggle 728,877 through 996,975
turn on 121,895 through 349,906
turn on 888,547 through 931,628
toggle 398,782 through 834,882
turn on 966,850 through 989,953
turn off 891,543 through 914,991
toggle 908,77 through 916,117
turn on 576,900 through 943,934
turn off 580,170 through 963,206
turn on 184,638 through 192,944
toggle 940,147 through 978,730
turn off 854,56 through 965,591
toggle 717,172 through 947,995
toggle 426,987 through 705,998
turn on 987,157 through 992,278
toggle 995,774 through 997,784
turn off 796,96 through 845,182
turn off 451,87 through 711,655
turn off 380,93 through 968,676
turn on 263,468 through 343,534
turn on 917,936 through 928,959
toggle 478,7 through 573,148
turn off 428,339 through 603,624
turn off 400,880 through 914,953
toggle 679,428 through 752,779
turn off 697,981 through 709,986
toggle 482,566 through 505,725
turn off 956,368 through 993,516
toggle 735,823 through 783,883
turn off 48,487 through 892,496
turn off 116,680 through 564,819
turn on 633,865 through 729,930
turn off 314,618 through 571,922
toggle 138,166 through 936,266
turn on 444,732 through 664,960
turn off 109,337 through 972,497
turn off 51,432 through 77,996
turn off 259,297 through 366,744
toggle 801,130 through 917,544
toggle 767,982 through 847,996
turn on 216,507 through 863,885
turn off 61,441 through 465,731
turn on 849,970 through 944,987
toggle 845,76 through 852,951
toggle 732,615 through 851,936
toggle 251,128 through 454,778
turn on 324,429 through 352,539
toggle 52,450 through 932,863
turn off 449,379 through 789,490
turn on 317,319 through 936,449
toggle 887,670 through 957,838
toggle 671,613 through 856,664
turn off 186,648 through 985,991
turn off 471,689 through 731,717
toggle 91,331 through 750,758
toggle 201,73 through 956,524
toggle 82,614 through 520,686
toggle 84,287 through 467,734
turn off 132,367 through 208,838
toggle 558,684 through 663,920
turn on 237,952 through 265,997
turn on 694,713 through 714,754
turn on 632,523 through 862,827
turn on 918,780 through 948,916
turn on 349,586 through 663,976
toggle 231,29 through 257,589
toggle 886,428 through 902,993
turn on 106,353 through 236,374
turn on 734,577 through 759,684
turn off 347,843 through 696,912
turn on 286,699 through 964,883
turn on 605,875 through 960,987
turn off 328,286 through 869,461
turn off 472,569 through 980,848
toggle 673,573 through 702,884
turn off 398,284 through 738,332
turn on 158,50 through 284,411
turn off 390,284 through 585,663
turn on 156,579 through 646,581
turn on 875,493 through 989,980
toggle 486,391 through 924,539
turn on 236,722 through 272,964
toggle 228,282 through 470,581
toggle 584,389 through 750,761
turn off 899,516 through 900,925
turn on 105,229 through 822,846
turn off 253,77 through 371,877
turn on 826,987 through 906,992
turn off 13,152 through 615,931
turn on 835,320 through 942,399
turn on 463,504 through 536,720
toggle 746,942 through 786,998
turn off 867,333 through 965,403
turn on 591,477 through 743,692
turn off 403,437 through 508,908
turn on 26,723 through 368,814
turn on 409,485 through 799,809
turn on 115,630 through 704,705
turn off 228,183 through 317,220
toggle 300,649 through 382,842
turn off 495,365 through 745,562
turn on 698,346 through 744,873
turn on 822,932 through 951,934
toggle 805,30 through 925,421
toggle 441,152 through 653,274
toggle 160,81 through 257,587
turn off 350,781 through 532,917
toggle 40,583 through 348,636
turn on 280,306 through 483,395
toggle 392,936 through 880,955
toggle 496,591 through 851,934
turn off 780,887 through 946,994
turn off 205,735 through 281,863
toggle 100,876 through 937,915
turn on 392,393 through 702,878
turn on 956,374 through 976,636
toggle 478,262 through 894,775
turn off 279,65 through 451,677
turn on 397,541 through 809,847
turn on 444,291 through 451,586
toggle 721,408 through 861,598
turn on 275,365 through 609,382
turn on 736,24 through 839,72
turn off 86,492 through 582,712
turn on 676,676 through 709,703
turn off 105,710 through 374,817
toggle 328,748 through 845,757
toggle 335,79 through 394,326
toggle 193,157 through 633,885
turn on 227,48 through 769,743
toggle 148,333 through 614,568
toggle 22,30 through 436,263
toggle 547,447 through 688,969
toggle 576,621 through 987,740
turn on 711,334 through 799,515
turn on 541,448 through 654,951
toggle 792,199 through 798,990
turn on 89,956 through 609,960
toggle 724,433 through 929,630
toggle 144,895 through 201,916
toggle 226,730 through 632,871
turn off 760,819 through 828,974
toggle 887,180 through 940,310
toggle 222,327 through 805,590
turn off 630,824 through 885,963
turn on 940,740 through 954,946
turn on 193,373 through 779,515
toggle 304,955 through 469,975
turn off 405,480 through 546,960
turn on 662,123 through 690,669
turn off 615,238 through 750,714
turn on 423,220 through 930,353
turn on 329,769 through 358,970
toggle 590,151 through 704,722
turn off 884,539 through 894,671
toggle 449,241 through 984,549
toggle 449,260 through 496,464
turn off 306,448 through 602,924
turn on 286,805 through 555,901
toggle 722,177 through 922,298
toggle 491,554 through 723,753
turn on 80,849 through 174,996
turn off 296,561 through 530,856
toggle 653,10 through 972,284
toggle 529,236 through 672,614
toggle 791,598 through 989,695
turn on 19,45 through 575,757
toggle 111,55 through 880,871
turn off 197,897 through 943,982
turn on 912,336 through 977,605
toggle 101,221 through 537,450
turn on 101,104 through 969,447
toggle 71,527 through 587,717
toggle 336,445 through 593,889
toggle 214,179 through 575,699
turn on 86,313 through 96,674
toggle 566,427 through 906,888
turn off 641,597 through 850,845
turn on 606,524 through 883,704
turn on 835,775 through 867,887
toggle 547,301 through 897,515
toggle 289,930 through 413,979
turn on 361,122 through 457,226
turn on 162,187 through 374,746
turn on 348,461 through 454,675
turn off 966,532 through 985,537
turn on 172,354 through 630,606
turn off 501,880 through 680,993
turn off 8,70 through 566,592
toggle 433,73 through 690,651
toggle 840,798 through 902,971
toggle 822,204 through 893,760
turn off 453,496 through 649,795
turn off 969,549 through 990,942
turn off 789,28 through 930,267
toggle 880,98 through 932,434
toggle 568,674 through 669,753
turn on 686,228 through 903,271
turn on 263,995 through 478,999
toggle 534,675 through 687,955
turn off 342,434 through 592,986
toggle 404,768 through 677,867
toggle 126,723 through 978,987
toggle 749,675 through 978,959
turn off 445,330 through 446,885
turn off 463,205 through 924,815
turn off 417,430 through 915,472
turn on 544,990 through 912,999
turn off 201,255 through 834,789
turn off 261,142 through 537,862
turn off 562,934 through 832,984
turn off 459,978 through 691,980
turn off 73,911 through 971,972
turn on 560,448 through 723,810
turn on 204,630 through 217,854
turn off 91,259 through 611,607
turn on 877,32 through 978,815
turn off 950,438 through 974,746
toggle 426,30 through 609,917
toggle 696,37 through 859,201
toggle 242,417 through 682,572
turn off 388,401 through 979,528
turn off 79,345 through 848,685
turn off 98,91 through 800,434
toggle 650,700 through 972,843
turn off 530,450 through 538,926
turn on 428,559 through 962,909
turn on 78,138 through 92,940
toggle 194,117 through 867,157
toggle 785,355 through 860,617
turn off 379,441 through 935,708
turn off 605,133 through 644,911
toggle 10,963 through 484,975
turn off 359,988 through 525,991
turn off 509,138 through 787,411
toggle 556,467 through 562,773
turn on 119,486 through 246,900
turn on 445,561 through 794,673
turn off 598,681 through 978,921
turn off 974,230 through 995,641
turn off 760,75 through 800,275
toggle 441,215 through 528,680
turn off 701,636 through 928,877
turn on 165,753 through 202,780
toggle 501,412 through 998,516
toggle 161,105 through 657,395
turn on 113,340 through 472,972
toggle 384,994 through 663,999
turn on 969,994 through 983,997
turn on 519,600 through 750,615
turn off 363,899 through 948,935
turn on 271,845 through 454,882
turn off 376,528 through 779,640
toggle 767,98 through 854,853
toggle 107,322 through 378,688
turn off 235,899 through 818,932
turn on 445,611 through 532,705
toggle 629,387 through 814,577
toggle 112,414 through 387,421
toggle 319,184 through 382,203
turn on 627,796 through 973,940
toggle 602,45 through 763,151
turn off 441,375 through 974,545
toggle 871,952 through 989,998
turn on 717,272 through 850,817
toggle 475,711 through 921,882
toggle 66,191 through 757,481
turn off 50,197 through 733,656
toggle 83,575 through 915,728
turn on 777,812 through 837,912
turn on 20,984 through 571,994
turn off 446,432 through 458,648
turn on 715,871 through 722,890
toggle 424,675 through 740,862
toggle 580,592 through 671,900
toggle 296,687 through 906,775"""

function parseInput(line, grid)
	tokens = split(line)
	endpoint = split(tokens[end], ',')
	startpoint = split(tokens[end-2], ',')
	# @printf("%s\t%d,%d\t%d,%d\n", line, parse(Int, startpoint[1]),parse(Int, startpoint[2]), parse(Int, endpoint[1]), parse(Int, endpoint[2]))
	for x in parse(Int, startpoint[1])+1:parse(Int, endpoint[1])+1, y in parse(Int, startpoint[2])+1:parse(Int, endpoint[2])+1
		# println(x, " ", y)
		if tokens[end-3] == "off"
			if grid[x,y]>0 grid[x,y] -= 1 end
		elseif tokens[end-3] == "on"
			grid[x,y] += 1
		elseif tokens[end-3] == "toggle"
			grid[x,y] += 2
		else
			println("Could not identify ", tokens[end-3])
		end
	end
end

grid = [ 0 for x in 1:1000, y in 1:1000 ]
# inpt = "turn on 499,499 through 500,500"
for line in split(inpt, "\n")
	parseInput(line, grid)
end
println(sum(grid))

max_stars_repo_name: jstrube/adventofcode

max_issues_repo_name: jstrube/adventofcode

max_forks_repo_name: jstrube/adventofcode


text:
# validation of properties and code

export InvalidIRError

function check_method(@nospecialize(job::CompilerJob))
    isa(job.source.f, Core.Builtin) && throw(KernelError(job, "function is not a generic function"))

    # get the method
    ms = Base.methods(job.source.f, job.source.tt)
    isempty(ms)   && throw(KernelError(job, "no method found"))
    length(ms)!=1 && throw(KernelError(job, "no unique matching method"))
    m = first(ms)

    # kernels can't return values
    if job.source.kernel
        cache = ci_cache(job)
        mt = method_table(job)
        interp = GPUInterpreter(cache, mt, job.source.world)
        rt = Base.return_types(job.source.f, job.source.tt, interp)[1]
        if rt != Nothing
            throw(KernelError(job, "kernel returns a value of type `$rt`",
                """Make sure your kernel function ends in `return`, `return nothing` or `nothing`.
                   If the returned value is of type `Union{}`, your Julia code probably throws an exception.
                   Inspect the code with `@device_code_warntype` for more details."""))
        end
    end

    return
end

# The actual check is rather complicated
# and might change from version to version...
function hasfieldcount(@nospecialize(dt))
    try
        fieldcount(dt)
    catch
        return false
    end
    return true
end

function explain_nonisbits(@nospecialize(dt), depth=1; maxdepth=10)
    dt===Module && return ""    # work around JuliaLang/julia#33347
    depth > maxdepth && return ""
    hasfieldcount(dt) || return ""
    msg = ""
    for (ft, fn) in zip(fieldtypes(dt), fieldnames(dt))
        if !isbitstype(ft)
            msg *= "  "^depth * ".$fn is of type $ft which is not isbits.\n"
            msg *= explain_nonisbits(ft, depth+1)
        end
    end
    return msg
end

function check_invocation(@nospecialize(job::CompilerJob), entry::LLVM.Function)
    # make sure any non-isbits arguments are unused
    real_arg_i = 0
    sig = Base.signature_type(job.source.f, job.source.tt)::Type
    for (arg_i,dt) in enumerate(sig.parameters)
        isghosttype(dt) && continue
        Core.Compiler.isconstType(dt) && continue
        real_arg_i += 1

        if !isbitstype(dt)
            throw(KernelError(job, "passing and using non-bitstype argument",
                """Argument $arg_i to your kernel function is of type $dt, which is not isbits:
                    $(explain_nonisbits(dt))"""))
        end
    end

    return
end

## IR validation

const IRError = Tuple{String, StackTraces.StackTrace, Any} # kind, bt, meta

struct InvalidIRError <: Exception
    job::CompilerJob
    errors::Vector{IRError}
end

const RUNTIME_FUNCTION = "call to the Julia runtime"
const UNKNOWN_FUNCTION = "call to an unknown function"
const POINTER_FUNCTION = "call through a literal pointer"
const DELAYED_BINDING  = "use of an undefined name"
const DYNAMIC_CALL     = "dynamic function invocation"

function Base.showerror(io::IO, err::InvalidIRError)
    print(io, "InvalidIRError: compiling ", err.job.source, " resulted in invalid LLVM IR")
    for (kind, bt, meta) in err.errors
        print(io, "\nReason: unsupported $kind")
        if meta !== nothing
            if kind == RUNTIME_FUNCTION || kind == UNKNOWN_FUNCTION || kind == POINTER_FUNCTION || kind == DYNAMIC_CALL
                print(io, " (call to ", meta, ")")
            elseif kind == DELAYED_BINDING
                print(io, " (use of '", meta, "')")
            end
        end
        Base.show_backtrace(io, bt)
    end
    return
end

function check_ir(job, args...)
    errors = check_ir!(job, IRError[], args...)
    unique!(errors)
    if !isempty(errors)
        throw(InvalidIRError(job, errors))
    end

    return
end

function check_ir!(job, errors::Vector{IRError}, mod::LLVM.Module)
    for f in functions(mod)
        check_ir!(job, errors, f)
    end

    return errors
end

function check_ir!(job, errors::Vector{IRError}, f::LLVM.Function)
    for bb in blocks(f), inst in instructions(bb)
        if isa(inst, LLVM.CallInst)
            check_ir!(job, errors, inst)
        end
    end

    return errors
end

const libjulia = Ref{Ptr{Cvoid}}(C_NULL)

function check_ir!(job, errors::Vector{IRError}, inst::LLVM.CallInst)
    bt = backtrace(inst)
    dest = called_value(inst)
    if isa(dest, LLVM.Function)
        fn = LLVM.name(dest)

        # some special handling for runtime functions that we don't implement
        if fn == "jl_get_binding_or_error"
            try
                m, sym, _ = operands(inst)
                sym = first(operands(sym::ConstantExpr))::ConstantInt
                sym = convert(Int, sym)
                sym = Ptr{Cvoid}(sym)
                sym = Base.unsafe_pointer_to_objref(sym)
                push!(errors, (DELAYED_BINDING, bt, sym))
            catch e
                @debug "Decoding arguments to jl_get_binding_or_error failed" inst bb=LLVM.parent(inst)
                push!(errors, (DELAYED_BINDING, bt, nothing))
            end
        elseif fn == "jl_invoke"
            try
                f, args, nargs, meth = operands(inst)
                meth = first(operands(meth::ConstantExpr))::ConstantInt
                meth = convert(Int, meth)
                meth = Ptr{Cvoid}(meth)
                meth = Base.unsafe_pointer_to_objref(meth)::Core.MethodInstance
                push!(errors, (DYNAMIC_CALL, bt, meth.def))
            catch e
                @debug "Decoding arguments to jl_invoke failed" inst bb=LLVM.parent(inst)
                push!(errors, (DYNAMIC_CALL, bt, nothing))
            end
        elseif fn == "jl_apply_generic"
            try
                f, args, nargs, _ = operands(inst)
                f = first(operands(f))::ConstantInt # get rid of inttoptr
                f = convert(Int, f)
                f = Ptr{Cvoid}(f)
                f = Base.unsafe_pointer_to_objref(f)
                push!(errors, (DYNAMIC_CALL, bt, f))
            catch e
                @debug "Decoding arguments to jl_apply_generic failed" inst bb=LLVM.parent(inst)
                push!(errors, (DYNAMIC_CALL, bt, nothing))
            end

        # detect calls to undefined functions
        elseif isdeclaration(dest) && !LLVM.isintrinsic(dest) && !isintrinsic(job, fn)
            # figure out if the function lives in the Julia runtime library
            if libjulia[] == C_NULL
                paths = filter(Libdl.dllist()) do path
                    name = splitdir(path)[2]
                    startswith(name, "libjulia")
                end
                libjulia[] = Libdl.dlopen(first(paths))
            end

            if Libdl.dlsym_e(libjulia[], fn) != C_NULL
                push!(errors, (RUNTIME_FUNCTION, bt, LLVM.name(dest)))
            else
                push!(errors, (UNKNOWN_FUNCTION, bt, LLVM.name(dest)))
            end
        end

    elseif isa(dest, InlineAsm)
        # let's assume it's valid ASM

    elseif isa(dest, ConstantExpr)
        # detect calls to literal pointers
        if occursin("inttoptr", string(dest))
            # extract the literal pointer
            ptr_arg = first(operands(dest))
            @compiler_assert isa(ptr_arg, ConstantInt) job
            ptr_val = convert(Int, ptr_arg)
            ptr = Ptr{Cvoid}(ptr_val)

            if !valid_function_pointer(job, ptr)
                # look it up in the Julia JIT cache
                frames = ccall(:jl_lookup_code_address, Any, (Ptr{Cvoid}, Cint,), ptr, 0)
                if length(frames) >= 1
                    @compiler_assert length(frames) == 1 job frames=frames
                    fn, file, line, linfo, fromC, inlined = last(frames)
                    push!(errors, (POINTER_FUNCTION, bt, fn))
                else
                    push!(errors, (POINTER_FUNCTION, bt, nothing))
                end
            end
        end
    end

    return errors
end


text:
function mm_spec_kj6_dimspec()::MPCCDimSpec
    dimspec = MPCCDimSpec(2, 1, 2, 0, 2, 1, 0)
    return dimspec
end

function mm_spec_kj6_defn(x, pr, ps)
    f = exp(-x[1]+x[2])

    ce = Vector{Num}()

    ci = Vector{Num}(undef, 2)
    ci[1] = ((x[1]-2)^2 + (x[2]+1)^2 -2pr[1] -6)
    ci[2] = (1 - x[1])

    F = Matrix{Num}(undef, 2, 1)
    F[1,1] = x[1]
    F[2,1] = x[2]

    return MPCCDefinition(f, ce, ci, F)
end

function mm_spec_kj6_nzmask()::MPCCModelNZMask
    nz_mask = MPCCModelNZMask(
        # ce
        [],
        # ci
        [   Set{Int64}([1, 2]),
            Set{Int64}([1]) ],
        # F
        [
            # F1
            [   Set{Int64}([1]) ],
            # F2
            [   Set{Int64}([2]) ]
        ] )
    return nz_mask
end

function mm_spec_kj6_testvectors()
    return Vector{MPCCModelTestVector}(undef, 0)
end

function mm_spec_kj6_knownsols()

    function local_kj6_knownsol(pr::Vector{R}, ps::Vector{Int64}) where {R <: Real}
        @assert ( 1 == length(pr) && 0 == length(ps) )
        t = pr[1]
        if ( t <= -2 )
            x_sol = [1.0, 0.0]
            f_sol = exp(-1)
        elseif ( t <= -(1//2) && t > -2 )
            x_sol = [ (-sqrt(5 + 2t) + 2), 0]
            f_sol = exp(-x_sol[1])
        else
            x_sol = [ 0, (sqrt(6 + 2*t - (0-2)^2) - 1) ]
            f_sol = exp(x_sol[2])
        end

        return (x_sol, f_sol)
    end

    return [ local_kj6_knownsol ]
end

function mm_spec_kj6_parameterisations(t)
    defn1 = MPCCParameterisationDefn(
            Vector{Num}([ t ]),
            (-5.0, 2.0),
            "Standard"
        )

    return  Vector{MPCCParameterisationDefn}( [ defn1 ] )
end


text:
module Satisfactory

using JuMP, Cbc, DataFrames, JSON3, LightGraphs, LayeredLayouts, Plots
using InteractiveUtils: subtypes

include("products.jl") ; export Product #, subtypes(Product)
include("recipes.jl") ; export recipes, dependantRecipes, allRecipes
include("models.jl") ; export maxModelLinear, maxModelMIP
include("utils.jl")

export baseRecipes, baseResources, maxResources
export maximize!, maximizeDiscrete!

const allRecipes = Recipe[]
const baseRecipes = Recipe[]
const harvestedProducts = (AlienCarapace, AlienOrgans, FICSITCoupon, FlowerPetals, GreenPowerSlug, Leaves,
    Mycelia, PurplePowerSlug, SAMOre, Wood, YellowPowerSlug)
const dictProductRecipes = Dict(p => Set{Tuple{Recipe, Float64}}() for p in subtypes(Product))
const dictProductDependantRecipes = Dict(p => Set{Tuple{Recipe, Float64}}() for p in subtypes(Product))
const baseResources = (Limestone, IronOre, CopperOre, CateriumOre, Coal, RawQuartz, Sulfur, Bauxite, Uranium, Water, CrudeOil)
const maxResources = Dict(zip(baseResources, (52860, 70380, 28860, 11040, 31680, 10500, 3300, 7800, 1200, 100_000, 7500)))

readData()

function maximize!(::Type{T} ; resources::Dict, alternates = String[],
    allowMultiRecipes = true, nodeLabelSize = 4, edgeLabelSize = 4) where T <: Product
    m, x, y = maxModelLinear(T ; resources, alternates, allowMultiRecipes)
    JuMP.set_silent(m)
    optimize!(m);
    print_results(m, 1., nodeLabelSize, edgeLabelSize, T)
end

function maximizeDiscrete!(::Type{T}, frac = 1/4 ; resources::Dict, alternates = String[],
    allowMultiRecipes = true, nodeLabelSize = 4, edgeLabelSize = 4, minProductionThreshold = 0.) where T <: Product

    m, x, y = maxModelMIP(T, frac ; resources, alternates, allowMultiRecipes, minProductionThreshold)
    JuMP.set_silent(m)
    optimize!(m);
    print_results(m, frac, nodeLabelSize, edgeLabelSize, T)
end

end


text:
# #####################################################################
# ####################### Mixed Potentials (GaN) ######################
# #####################################################################
# function virial(r::Vector{<:Real}, p::GaN, Type1::Symbol, Type2::Symbol)
#     f = force(r, p, Type1, Type2)
#     return f[1] * r[1] + f[2] * r[2] + f[3] * r[3]
# end

# function virial(r::Vector{Atom}, p::MixedPotential)
#     n = length(r)
#     v = 0.0
#     for i = 1:(n-1)
#         ri = r[i].Position
#         for j = (i+1):n
#             rj = r[j].Position
#             rtemp = ri - rj
#             v +=  virial(rtemp, p, r[i].Type, r[j].Type)
#         end
#     end
#     return v
# end

# function virial(c::Configuration, p::MixedPotential)
#     return virial(c.Atoms, p)
# end

# function virial(r::Vector{Configuration}, p::MixedPotential)
#     n = length(r)
#     v = zeros(n)
#     for i = 1:n
#         v[i] = virial(r[i].Atoms, p)
#     end
#     return v
# end

# function virial_stress(r::Vector{<:Real}, p::GaN, Type1::Symbol, Type2::Symbol)
#     f = force(r, p, Type1, Type2)
#     vi = r * f'
#     return [vi[1, 1], vi[2, 2], vi[3, 3], vi[2, 3], vi[1, 3], vi[1,2]]
# end

# function virial_stress(r::Vector{Atom}, p::MixedPotential)
#     n = length(r)
#     v = zeros(6)
#     for i = 1:(n-1)
#         ri = r[i].Position
#         for j = (i+1):n
#             rj = r[j].Position
#             rtemp = ri - rj
#             v +=  virial_stress(rtemp, p, r[i].Type, r[j].Type)
#         end
#     end
#     return v
# end

# function virial_stress(c::Configuration, p::MixedPotential)
#     return virial_stress(c.Atoms, p)
# end

# function virial_stress(r::Vector{Configuration}, p::MixedPotential)
#     n = length(r)
#     v = [zeros(6) for i = 1:n]
#     for i = 1:n
#         v[i] = virial_stress(r[i].Atoms, p)
#     end
#     return v
# end


text:
ProjDir = joinpath(dirname(@__FILE__), "..", "Examples", "Line3")
cd(ProjDir) do

  println("Moving to directory: $(ProjDir)")

  isdir("tmp") &&
    rm("tmp", recursive=true);

  include(joinpath(ProjDir, "jline3.jl"))

  isdir("tmp") &&
    rm("tmp", recursive=true);

end


text:
import MacroTools
using Base

const TRACE_SYMBOL = :trace
const STACK_SYMBOL = :handlers_stack
const RETURN_KEY = :(:_return)

"""
    @jyro

Takes a function and transforms it into a `MinijyroModel`. It replaces all the
`~` expressions in the function body and replaces them with `sample!` statements.
Furthermore, the new generated function will return a `Dict` which is called the
`trace` of the program. If the function has any return values they will now be
stored in `trace[:_return]`.
"""
macro jyro(expr)
    fn_dict = MacroTools.splitdef(expr)
    pushfirst!(fn_dict[:args], STACK_SYMBOL)

    loop_var = gensym()
    body = translate_tilde(fn_dict[:body])
    body = translate_return(body, loop_var)

    fn_dict[:body] = quote
        $TRACE_SYMBOL = Dict{Any,Any}($(RETURN_KEY) => nothing)
        for $(loop_var) in $(STACK_SYMBOL)[end:-1:1]
            enter!($TRACE_SYMBOL, $(loop_var))
        end
        $(body.args...)
        # NOTE: If the last line of (body.args...) is already a return statement
        #       then the compiler is smart enough to remove this code.
        for $(loop_var) in $STACK_SYMBOL
            exit!($TRACE_SYMBOL, $(loop_var))
        end
        return $TRACE_SYMBOL
    end

    model_name = fn_dict[:name]
    fn_name = gensym(model_name)
    fn_dict[:name] = fn_name
    fn_dict[:rtype] = Dict{Any,Any} # Change return type of generated function to a Dict.
    fn_expr = esc(MacroTools.combinedef(fn_dict))

    return quote
        $fn_expr

        $(esc(model_name)) = MinijyroModel(AbstractHandler[], $(esc(fn_name)))
    end
end

function translate_tilde(expr)
    # NOTE: Code adapted from https://github.com/probcomp/Gen/blob/master/src/dsl/dsl.jl
    MacroTools.postwalk(expr) do e
        if MacroTools.@capture(e, {addr_} ~ rhs_)
            # NOTE: This does not sample
            :(Minijyro.sample!($TRACE_SYMBOL, $STACK_SYMBOL, $(addr), $rhs))
        elseif MacroTools.@capture(e, lhs_ ~ rhs_)
            name_symbol = QuoteNode(lhs)
            :($lhs = Minijyro.sample!($TRACE_SYMBOL, $STACK_SYMBOL, $name_symbol, $rhs))
        else
            e
        end
    end
end

function translate_return(expr, loop_var)
    MacroTools.postwalk(expr) do e
        if MacroTools.@capture(e, return r_)
            quote
                for $(loop_var) in $STACK_SYMBOL
                    exit!($TRACE_SYMBOL, $(loop_var))
                end
                $(TRACE_SYMBOL)[$(RETURN_KEY)] = $r
                return $TRACE_SYMBOL
            end
        else
            e
        end
    end
end

# TODO: Macro for handling code segments