text:
"""
    module FMM3D

Wrappers for the Flatiron Institute's FMM3D
library.

# wrappers

All N-body codes return output in an
`FMMVals` structure. See documentation
of N-body codes for details.

N-body interactions with the Helmholtz kernel
- [`hfmm3d`](@ref): ``O(N)`` fast mutlipole code
- [`h3ddir`](@ref): ``O(N^2)`` direct code

N-body interactions with the Laplace kernel
- [`lfmm3d`](@ref): ``O(N)`` fast mutlipole code
- [`l3ddir`](@ref): ``O(N^2)`` direct code

N-body interactions with Stokes kernels
- [`stfmm3d`](@ref): ``O(N)`` fast mutlipole code
- [`st3ddir`](@ref): ``O(N^2)`` direct code

N-body interactions with an electromagnetic kernel
- [`emfmm3d`](@ref): ``O(N)`` fast mutlipole code
- [`em3ddir`](@ref): ``O(N^2)`` direct code

# lower level routines

For a list of lower level routines see the
documentation for [`lower_level_routs`](@ref)
"""
module FMM3D

using FMM3D_jll

export FMMVals, hfmm3d, lfmm3d, h3ddir, l3ddir, lowerlevel_routs
export stfmm3d, st3ddir, emfmm3d, em3ddir
export besseljs3d

"""

# Available lower level routines

* [`besseljs3d`](@ref) spherical Bessel function eval
"""
function lower_level_routs() end

# fortran input/return types

Fd = Ref{Float64}
Fi = Ref{Int32}
Fc = Ref{ComplexF64}

# common input types

TFN = Union{Array{Float64},Nothing}
TCN = Union{Array{ComplexF64},Nothing}

"""
    FMMVals()

Return type for FMM and direct computation function
calls. Fields are `nothing` on return unless requested.
See individual FMM/direct computation function
documentation for specifics.
"""
mutable struct FMMVals
    pot
    grad
    hess

    pottarg
    gradtarg
    hesstarg

    ier

    # stokes-specific return values

    pre
    pretarg

    # maxwell-specific return values
    E
    divE
    curlE

    Etarg
    divEtarg
    curlEtarg
end

function FMMVals()
    return FMMVals(nothing,nothing,nothing,
                   nothing,nothing,nothing,
                   nothing,
                   nothing,nothing,
                   nothing,nothing,nothing,
                   nothing,nothing,nothing)
end

"""
```julia
    anyfail, n, nt, ifcharge, ifdipole = (
        scalarfmm3dinputcheck(sources,charges,dipvecs,
                                targets,pg,pgt,nd) )
```

Input checking for Helmholtz and Laplace routines.
Checks the sizes of these arrays and the compatibility
of the various flags, nd, and provided arrays.
If something is off, a warning is issued and
`anyfail` is set to true. Non-fatal mistakes result in
a warning but `anyfail` remains false.

Output:
* `anyfail` - boolean is true if any checks fail, false otherwise
* n - integer, number of sources
* nt - integer, number of targets
* ifcharge - integer, is 1 if there are charges, 0 otherwise
* ifdipole - integer, is 1 if there are dipoles, 0 otherwise
"""
function scalarfmm3dinputcheck(sources,charges,dipvecs,targets,pg,pgt,nd)
    anyfail = false

    if (size(sources,1) != 3)
        @warn "sources array is wrong shape, computing nothing"
        anyfail = true
    end
    if (nd < 0)
        @warn "nd has invalid value, computing nothing"
        anyfail = true
    end
    if (pg > 3 || pg < 0)
        @warn "flag pg not in expected range, computing nothing"
        anyfail = true
    end
    if (pgt > 3 || pgt < 0)
        @warn "flag pgt not in expected range, computing nothing"
        anyfail = true
    end

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

    if targets != nothing
        if (size(targets,1) != 3)
            @warn "targets array is wrong shape, computing nothing"
            anyfail = true
        end
        nt = div(length(targets),3)
    end

    if charges != nothing
        if (div(length(charges),nd) != n)
            @warn "size of charges array incompatible with sources array and nd parameter, computing nothing"
            anyfail = true
        end
        ifcharge = 1
    else
        ifcharge = 0
    end

    if dipvecs != nothing
        if (div(length(dipvecs),nd) != n*3)
            @warn "size of dipvecs array incompatible with sources array and nd parameter, computing nothing"
            anyfail = true
        end
        ifdipole = 1
    else
        ifdipole = 0
    end

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

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

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

    return anyfail, n, nt, ifcharge, ifdipole
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}``)

# Optional Keyword Input

* `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.
    + Do not compute any quantities at sources if `pg == 0`
    + Potential (``u``) at sources evaluated if `pg == 1`.
    + Potential and gradient (``\\nabla u``) at sources evaluated if `pg == 2`
    + Potential, gradient, and Hessian (``\\nabla \\nabla u``) at sources evaluated if `pg == 3`
* `pgt::Integer` target eval flag.
    + Do not compute any quantities at targets if `pgt == 0`
    + Potential at targets evaluated if `pgt == 1`.
    + Potential and gradient at targets evaluated if `pgt == 2`
    + Potential, gradient, and Hessian at targets evaluated if `pgt == 3`
* `nd::Integer` number of densities

Note: if all default values are used for optional input, nothing
is computed.

# Output

`vals::FMMVals` with the fields

Note that the Hessian is returned as a length
6 vector at each point with the second derivatives
in the order: ``\\partial_{xx}``, ``\\partial_{yy}``,
``\\partial_{zz}``, ``\\partial_{xy}``,
``\\partial_{xz}``, ``\\partial_{yz}``.

* `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.hess::Array{ComplexF64}` size (nd,6,n) or (6,n) Hessian 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
* `vals.hesstarg::Array{ComplexF64}` size (nd,6,nt) or (6,nt) Hessian at target locations if requested
* `vals.ier` error flag as returned by FMM3D library. A value of 0 indicates a successful call.
Non-zero values may indicate insufficient memory available. See the documentation for the FMM3D library.
If not set (`nothing`), then FMM3D library was never called.
"""
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)

    zk = complex(zk)

    # default values

    vals = FMMVals()

    ifcharge = 0
    ifdipole = 0

    zero = ComplexF64(0)

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

    # check inputs

    if (pg == 3)
        @warn "Hessian not implemented for Helmholtz fmm, only computing potential and gradients at sources"
        pg = 2
    end
    if (pgt == 3)
        @warn "Hessian not implemented for Helmholtz fmm, only computing potential and gradients at targets"
        pgt = 2
    end

    anyfail, n, nt, ifcharge, ifdipole = (
        scalarfmm3dinputcheck(sources,charges,dipvecs,targets,pg,pgt,nd))

    if anyfail
        return vals
    end

    if (ifcharge == 0); charges = zero end
    if (ifdipole == 0); dipvecs = zero end
    if (nt == 0); targets = 0.0 end

    # allocate memory for return values

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

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

    if pg == 3
        if nd > 1
            hess = zeros(ComplexF64,nd,6,n)
        else
            hess = zeros(ComplexF64,6,n)
        end
    end

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

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

    if pgt == 3
        if nd > 1
            hesstarg = zeros(ComplexF64,nd,6,nt)
        else
            hesstarg = zeros(ComplexF64,6,nt)
        end
    end

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

    ier = Integer(0)
    iper = Integer(0)

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

    vals.ier = ier
    if (ier != 0); @warn "libfmm3d had an error, see vals.ier"; end

    # load requested values

    if pg == 1 || pg == 2 || pg == 3; vals.pot = pot end
    if pg == 2 || pg == 3; vals.grad = grad end
    if pg == 3; vals.hess = hess end
    if pgt == 1 || pgt == 2 || pgt == 3; vals.pottarg = pottarg end
    if pgt == 2 || pgt == 3; vals.gradtarg = gradtarg end
    if pgt == 3; vals.hesstarg = hesstarg 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. By convention this code only computes
the effect of sources on targets. If the value at sources is
also needed, the routine can be called again with targets equal
to the source locations.

```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``)

# Optional Keyword Input

* `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`
    + Potential, gradient, and Hessian at targets evaluated if `pgt == 3`
* `nd::Integer` number of densities
* `thresh::Float64` threshold for ignoring interactions when ``\\|x-x_{j}\\| \\leq thresh``

Note: if all default values are used for optional input, nothing
is computed.

# Output

`vals::FMMVals` with the fields

Note that the Hessian is returned as a length
6 vector at each point with the second derivatives
in the order: ``\\partial_{xx}``, ``\\partial_{yy}``,
``\\partial_{zz}``, ``\\partial_{xy}``,
``\\partial_{xz}``, ``\\partial_{yz}``.

* `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
* `vals.hesstarg::Array{ComplexF64}` size (nd,6,nt) or (6,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)

    zk = complex(zk)

    # default values

    vals = FMMVals()

    ifcharge = 0
    ifdipole = 0

    zero = ComplexF64(0)

    pottarg = zero
    gradtarg = zero
    hesstarg = zero

    # check inputs

    if (pgt == 3)
        @warn "Hessian not implemented for Helmholtz fmm, only computing potential and gradients at targets"
        pgt = 2
    end

    pg = 0
    anyfail, n, nt, ifcharge, ifdipole = (
        scalarfmm3dinputcheck(sources,charges,dipvecs,targets,pg,pgt,nd))

    if anyfail
        return vals
    end

    if (ifcharge == 0); charges = zero end
    if (ifdipole == 0); dipvecs = zero end
    if (nt == 0); targets = 0.0 end

    # allocate memory for return values

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

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

    if pgt == 3
        if nd > 1
            hesstarg = zeros(ComplexF64,nd,6,nt)
        else
            hesstarg = zeros(ComplexF64,6,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
    elseif pgt == 3
        if ifcharge == 1
            if ifdipole == 1
                h3ddirectcdh!(nd,zk,sources,charges,
                              dipvecs,n,targets,nt,
                              pottarg,gradtarg,hesstarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
                vals.hesstarg = hesstarg
            else
                h3ddirectch!(nd,zk,sources,charges,
                             n,targets,nt,pottarg,
                             gradtarg,hesstarg,thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
                vals.hesstarg = hesstarg
            end
        else
            if ifdipole == 1
                h3ddirectdh!(nd,zk,sources,
                              dipvecs,n,targets,nt,
                             pottarg,gradtarg,hesstarg,
                             thresh)
                vals.pottarg = pottarg
                vals.gradtarg = gradtarg
                vals.hesstarg = hesstarg
            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

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

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

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

"""
```julia
    vals = lfmm3d(eps,sources;charges=nothing,dipvecs=nothing,
                  targets=nothing,pg=0,pgt=0,nd=1)
```

This function computes the N-body Laplace interactions
in three dimensions where the interaction kernel is given
by ``\$1/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} / \\|x-x_{j}\\| + v_{j}
\\cdot \\nabla( 1/\\|x-x_{j}\\|)  \\, ,
u(x) = \\sum_{j=1}^{N} c_{j} \\frac{1}{\\|x-x_{j}\\|}
- v_{j} \\cdot \\nabla \\left( \\frac{1}{\\|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
* `sources::Array{Float64}` size (3,n) source locations (``x_{j}``)

# Optional Keyword Input

* `charges::Array{Float64}` size (nd,n) or (n) charge densities (c_{j})
* `dipvecs::Array{Float64}` 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

Note: if all default values are used for optional input, nothing
is computed.

# Output

`vals::FMMVals` with the fields

Note that the Hessian is returned as a length
6 vector at each point with the second derivatives
in the order: ``\\partial_{xx}``, ``\\partial_{yy}``,
``\\partial_{zz}``, ``\\partial_{xy}``,
``\\partial_{xz}``, ``\\partial_{yz}``.

* `vals.pot::Array{Float64}` size (nd,n) potential at source locations if requested
* `vals.grad::Array{Float64}` size (nd,3,n) gradient at source locations if requested
* `vals.hess::Array{Float64}` size (nd,6,n) Hessian at source locations if requested
* `vals.pottarg::Array{Float64}` size (nd,nt) potential at target locations if requested
* `vals.gradtarg::Array{Float64}` size (nd,3,nt) gradient at target locations if requested
* `vals.hesstarg::Array{Float64}` size (nd,6,nt) Hessian at target locations if requested
* `vals.ier` error flag as returned by FMM3D library. A value of 0 indicates a successful call.
Non-zero values may indicate insufficient memory available. See the documentation for the FMM3D library.
If not set (`nothing`), then FMM3D library was never called.
"""
function lfmm3d(eps::Float64,sources::Array{Float64};charges::TFN=nothing,
                dipvecs::TFN=nothing,targets::TFN=nothing,
                pg::Integer=0,pgt::Integer=0,nd::Integer=1)

    # default values

    vals = FMMVals()

    ifcharge = 0
    ifdipole = 0

    zero = Float64(0)

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

    # check inputs

    anyfail, n, nt, ifcharge, ifdipole = (
        scalarfmm3dinputcheck(sources,charges,dipvecs,targets,pg,pgt,nd))

    if anyfail
        return vals
    end

    if (ifcharge == 0); charges = zero end
    if (ifdipole == 0); dipvecs = zero end
    if (nt == 0); targets = 0.0 end

    # allocate memory for return values

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

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

    if pg == 3
        if nd > 1
            hess = zeros(Float64,nd,6,n)
        else
            hess = zeros(Float64,6,n)
        end
    end

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

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

    if pgt == 3
        if nd > 1
            hesstarg = zeros(Float64,nd,6,nt)
        else
            hesstarg = zeros(Float64,6,nt)
        end
    end

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

    ier = Integer(0)
    iper = Integer(0)

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

    if (ier != 0); @warn "libfmm3d had an error, see vals.ier"; end
    vals.ier = ier

    # load requested values

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

    return vals

end

"""
```julia
    vals = l3ddir(sources,targets;charges=nothing,
                dipvecs=nothing,pgt=0,nd=1,
                thresh=1e-16)
```

This function computes the N-body Laplace interactions
in three dimensions where the interaction kernel is given
by ``1/r`` and its gradients.  This is the
``O(N^2)`` direct evaluation code. By convention this code only computes
the effect of sources on targets. If the value at sources is
also needed, the routine can be called again with targets equal
to the source locations.

```math
 u(x) = \\sum_{j=1}^{N} c_{j} \\frac{1}{\\|x-x_{j}\\|}
- v_{j} \\cdot \\nabla \\left( \\frac{1}{\\|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

* `sources::Array{Float64}` size (3,n) source locations (``x_{j}``)
* `targets::Array{Float64}` size (3,nt) target locations (``x``)

# Optional Keyword Input

* `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``

Note: if all default values are used for optional input, nothing
is computed.

# Output

`vals::FMMVals` 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 l3ddir(sources::Array{Float64},targets::Array{Float64};
                charges::TFN=nothing,dipvecs::TFN=nothing,
                pgt::Integer=0,nd::Integer=1,
                thresh::Float64=1e-16)

    # default values

    vals = FMMVals()

    ifcharge = 0
    ifdipole = 0

    zero = Float64(0)

    pottarg = zero
    gradtarg = zero
    hesstarg = zero

    # check inputs

    pg = 0
    anyfail, n, nt, ifcharge, ifdipole = (
        scalarfmm3dinputcheck(sources,charges,dipvecs,targets,pg,pgt,nd))

    if anyfail
        return vals
    end

    if (ifcharge == 0); charges = zero end
    if (ifdipole == 0); dipvecs = zero end
    if (nt == 0); targets = 0.0 end

    # allocate memory for return values

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

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

    # dispatch to appropriate wrapper

    if pgt == 1
        if ifcharge == 1
            if ifdipole == 1
                l3ddirectcdp!(nd,sources,charges,
                              dipvecs,n,targets,nt,
                              pottarg,thresh)
            else
                l3ddirectcp!(nd,sources,charges,
                             n,targets,nt,pottarg,
                             thresh)
            end
        else
            if ifdipole == 1
                l3ddirectdp!(nd,sources,
                             dipvecs,n,targets,nt,
                             pottarg,thresh)
            end
        end
    elseif pgt == 2
        if ifcharge == 1
            if ifdipole == 1
                l3ddirectcdg!(nd,sources,charges,
                              dipvecs,n,targets,nt,
                              pottarg,gradtarg,thresh)
            else
                l3ddirectcg!(nd,sources,charges,
                             n,targets,nt,pottarg,
                             gradtarg,thresh)
            end
        else
            if ifdipole == 1
                l3ddirectdg!(nd,sources,
                              dipvecs,n,targets,nt,
                             pottarg,gradtarg,thresh)
            end
        end
    end

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

    return vals

end

function l3ddirectcp!(nd,sources,charges,n,targets,nt,
                      pot,thresh)
    ccall((:l3ddirectcp_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd),
          nd,sources,charges,n,targets,nt,pot,thresh)
    return
end

function l3ddirectdp!(nd,sources,dipvecs,n,targets,nt,
                      pot,thresh)
    ccall((:l3ddirectdp_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd),
          nd,sources,dipvecs,n,targets,nt,pot,thresh)
    return
end

function l3ddirectcdp!(nd,sources,charges,dipvecs,n,targets,
                       nt,pot,thresh)
    ccall((:l3ddirectcdp_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd),
          nd,sources,charges,dipvecs,n,targets,nt,pot,thresh)
    return
end

function l3ddirectcg!(nd,sources,charges,n,targets,nt,
                      pot,grad,thresh)
    ccall((:l3ddirectcg_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd,Fd),
          nd,sources,charges,n,targets,nt,pot,grad,thresh)
    return
end

function l3ddirectdg!(nd,sources,dipvecs,n,targets,nt,
                      pot,grad,thresh)
    ccall((:l3ddirectdg_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd,Fd),
          nd,sources,dipvecs,n,targets,nt,pot,grad,thresh)
    return
end

function l3ddirectcdg!(nd,sources,charges,dipvecs,n,targets,
                       nt,pot,grad,thresh)
    ccall((:l3ddirectcdg_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fd,Fi,
                                          Fd,Fi,Fd,Fd,Fd),
          nd,sources,charges,dipvecs,n,targets,nt,pot,grad,thresh)
    return
end

"""
    function fj, fjder = besseljs3d(nterms,z;scale=1.0,ifder=0)

This subroutine evaluates the first `nterms` spherical Bessel
functions, and if requested, their derivatives.
It incorporates a scaling parameter `scale` so that

      	fjs_n(z)=j_n(z)/SCALE^n
      	fjder_n(z)=\\frac{\\partial fjs_n(z)}{\\partial z}

# Input

* `nterms::Integer` order of expansion of output array `fjs`
* `z::ComplexF64` argument of the spherical Bessel functions
* `scale::Float64` scaling factor
* `ifder::Integer1` flag indicating whether to calculate `fjder`
      	          0	NO
      	          1	YES
OUTPUT:

* `fjs::Array{ComplexF64}` array of length `nterms+1` of scaled Bessel functions.
* `fjder::Array{ComplexF64}` array of derivatives of scaled Bessel functions, if requested.
"""
function besseljs3d(nterms::Integer,z::ComplexF64;scale::Float64=1.0,
                    ifder::Integer=0)

    if (nterms < 0)
        @warn "nterms has invalid value, computing nothing"
        fjs = nothing; fjder = nothing;
        return fjs, fjder
    end
    if (ifder !=0 && ifder != 1)
        @warn "unexpected value in ifder, no ders computed"
    end

    fjs = Array{ComplexF64}(undef,nterms+1)
    fjder = ComplexF64(0)

    if ifder == 1
        fjder = Array{ComplexF64}(undef,nterms+1)
    end

    # fortran interface
    # subroutine besseljs3d(nterms,z,scale,fjs,ifder,fjder)

    ccall((:besseljs3d_,libfmm3d),Cvoid,(Fi,Fc,Fd,Fc,Fi,Fc),
          nterms,z,scale,fjs,ifder,fjder)

    if ifder != 1
        fjder = nothing
    end

    return fjs, fjder

end

"""
```julia
    anyfail, n, nt, ifstoklet, ifstrslet = (
        stfmm3dinputcheck(sources,stoklet,strslet,strsvec,
                                targets,ppreg,ppregt,nd) )
```

Input checking for Stokes routines.
Checks the sizes of these arrays and the compatibility
of the various flags, nd, and provided arrays.
If something is off, a warning is issued and
`anyfail` is set to true. Non-fatal mistakes result in
a warning but `anyfail` remains false.

Output:
* `anyfail` - boolean is true if any checks fail, false otherwise
* n - integer, number of sources
* nt - integer, number of targets
* ifstoklet - integer, is 1 if there are Stokeslets, 0 otherwise
* ifstrslet - integer, is 1 if there are (type I) stresslets, 0 otherwise
"""
function stfmm3dinputcheck(sources,stoklet,strslet,strsvec,targets,ppreg,ppregt,nd)
    anyfail = false

    if (size(sources,1) != 3)
        @warn "sources array is wrong shape, computing nothing"
        anyfail = true
    end
    if (nd < 0)
        @warn "nd has invalid value, computing nothing"
        anyfail = true
    end
    if (ppreg > 3 || ppreg < 0)
        @warn "flag ppreg not in expected range, computing nothing"
        anyfail = true
    end
    if (ppregt > 3 || ppregt < 0)
        @warn "flag ppregt not in expected range, computing nothing"
        anyfail = true
    end

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

    if targets != nothing
        if (size(targets,1) != 3)
            @warn "targets array is wrong shape, computing nothing"
            anyfail = true
        end
        nt = div(length(targets),3)
    end

    if stoklet != nothing
        if (div(length(stoklet),nd) != 3*n)
            @warn "size of stoklet array incompatible with sources array and nd parameter, computing nothing"
            anyfail = true
        end
        ifstoklet = 1
    else
        ifstoklet = 0
    end

    if (strslet != nothing || strsvec != nothing)
        if (strsvec == nothing || strslet == nothing)
            @warn "if stresslets in calculation, need both strslet and strsvec arrays, computing nothing"
            anyfail = true
        end
        if (div(length(strslet),nd) != n*3)
            @warn "size of strslet array is incompatible with sources array and nd paramter, computing nothing"
            anyfail = true
        end
        if (div(length(strsvec),nd) != n*3)
            @warn "size of strsvec array is incompatible with sources array and nd paramter, computing nothing"
            anyfail = true
        end
        ifstrslet = 1
    else
        ifstrslet = 0
    end

    if ifstoklet == 0 && ifstrslet == 0
        @warn "no Stokeslets or stresslets provided, doing nothing"
        anyfail = true
    end

    if (ppreg != 1 && ppreg != 2 && ppreg != 3) && (ppregt != 1 && ppregt != 2 && ppregt != 3)
        @warn "no output requested, doing nothing"
        anyfail = true
    end

    if (ppregt == 1 || ppregt == 2 || ppregt == 3) && targets == nothing
        if !anyfail
            @warn "target values requested but no targets provided, proceeding anyway"
        else
            @warn "target values requested but no targets provided. There are other points of failure. Computing nothing."
        end
    end

    return anyfail, n, nt, ifstoklet, ifstrslet
end

"""
```julia
    vals = stfmm3d(eps,sources;stoklet=nothing,strslet=nothing,
                   strsvec=nothing,targets=nothing,ppreg=0,
                   ppregt=0,nd=1)
```

This function computes the N-body Stokes interactions
in three dimensions where the interaction kernels are
the Stokeslet and stresslet (see below). This is the
``O(N)`` fast multipole code which computes the interactions
to the requested precision.

We take the following conventions for the Stokes kernels:

For a source ``y`` and target ``x``, let ``r_i = x_i-y_i``
 and let ``r = \\sqrt{r_1^2 + r_2^2 + r_3^2}``

The Stokeslet, ``G_{ij}``, and its associated pressure tensor,
``P_j``, (without the ``1/4\\pi`` scaling) are

```math
G_{ij}(x,y) = (r_i r_j)/(2r^3) + \\delta_{ij}/(2r) \\; ,
\\quad
P_j(x,y) = r_j/r^3
```
The (Type I) stresslet, ``T_{ijk}``, and its associated
pressure tensor, ``\\Pi_{jk}``, (without the ``1/4\\pi``
 scaling) are

```math
 T_{ijk}(x,y) = -3 r_i r_j r_k/ r^5 \\; , \\quad
 PI_{jk} = -2 \\delta_{jk} + 6 r_j r_k/r^5
```

The output of this routine gives the velocity
```math
 u_i(x) = \\sum_{m=1}^{N} G_{ij}(x,x_m) \\sigma^{(m)}_j
    + T_{ijk}(x,x_m) \\mu^{(m)}_j \\nu^{(m)}_k \\, ,
```
and, optionally, the pressure

```math
 p(x) = \\sum_{m=1}^{N} P_{j}(x,x_m) \\sigma^{(m)}_j
    + \\Pi{jk}(x,x_m) \\mu^{(m)}_j \\nu^{(m)}_k \\, ,
```
where ``\\sigma^{(m)}`` is a Stokeslet density 3-vector,
``\\mu^{(m)}`` and ``\\nu^{(m)}`` are the stresslet
density and orientation 3-vectors, and
``x_{m}`` are the source locations.
When ``x=x_{m}``, the term corresponding to
``x_{m}`` is dropped from the sum.
The gradient of the velocity can also be computed by request.

# Input

* `eps::Float64` precision requested
* `sources::Array{Float64}` size (3,n) source locations (``x_{m}``)

# Optional Keyword Input

* `stoklet::Array{Float64}` size (nd,3,n) or (3,n) Stokeslet densities (``\\sigma``)
* `strslet::Array{Float64}` size (nd,3,n) or (3,n) Stresslet densities (``\\mu'')  if either strsvec or strslet is provided, the other must be as well
* `strsvec::Array{Float64}` size (nd,3,n) or (3,n) Stresslet orientations (``\\nu``) if either strsvec or strslet is provided, the other must be as well
* `targets::Array{Float64}` size (3,nt) target locations (``x``)
* `ppreg::Integer` source eval flag.
    + Velocity (``u``) at sources evaluated if `ppreg == 1`.
    + Velocity and pressure (``p``) at sources evaluated if `ppreg == 2`
    + Velocity, pressure, and velocity gradient (``\\nabla u``)
    at sources evaluated if `ppreg == 3`
* `ppregt::Integer` source eval flag.
    + Velocity (``u``) at targets evaluated if `ppregt == 1`.
    + Velocity and pressure (``p``) at targets evaluated if `ppregt == 2`
    + Velocity, pressure, and velocity gradient (``\\nabla u``)
    at targets evaluated if `ppregt == 3`
* `nd::Integer` number of densities of each type

Note: if all default values are used for optional input, nothing
is computed.

# Output

`vals::FMMVals` with the fields

* `vals.pot::Array{Float64}` size (nd,3,n) or (3,n) velocity at source locations if requested
* `vals.grad::Array{Float64}` size (nd,3,3,n) or (3,n) gradient of velocity at source locations if requested.
`grad[l,i,j,k]` is the i-th component of the gradient of the j-th component of the velocity
for the l-th density at the k-th source location
* `vals.pre::Array{Float64}` size (nd,n) or (n) pressure at source locations if requested
* `vals.pottarg::Array{Float64}` size (nd,3,nt) or (3,nt) velocity at target locations if requested
* `vals.gradtarg::Array{Float64}` size (nd,3,3,nt) or (3,nt) gradient of velocity at target locations if requested.
`gradtarg[l,i,j,k]` is the i-th component of the gradient of the j-th component of the velocity
for the l-th density at the k-th target location
* `vals.pretarg::Array{Float64}` size (nd,nt) or (nt) pressure at target locations if requested
* `vals.ier` error flag as returned by FMM3D library. A value of 0 indicates a successful call.
Non-zero values may indicate insufficient memory available. See the documentation for the FMM3D library.
If not set (`nothing`), then FMM3D library was never called.
"""
function stfmm3d(eps::Float64,
                sources::Array{Float64};
                 stoklet::TFN=nothing,strslet::TFN=nothing,
                 strsvec::TFN=nothing,
                 targets::TFN=nothing,ppreg::Integer=0,
                 ppregt::Integer=0,nd::Integer=1)

    # default values

    vals = FMMVals()

    ifstoklet = 0
    ifstrslet = 0

    zero = Float64(0)
    zero3 = zeros(Float64,3)

    pot = zero
    grad = zero
    pre = zero
    pottarg = zero
    gradtarg = zero
    pretarg = zero

    ier = Integer(0)

    # check inputs

    anyfail, n, nt, ifstoklet, ifstrslet = (
        stfmm3dinputcheck(sources,stoklet,strslet,strsvec,targets,ppreg,ppregt,nd))

    if anyfail
        return vals
    end

    if (ifstoklet == 0); stoklet = zero3 end
    if (ifstrslet == 0); strslet = zero3; strsvec=zero3 end
    if (nt == 0); targets = zero3 end

    # allocate memory for return values

    if ppreg == 1 || ppreg == 2 || ppreg == 3
        if nd > 1
            pot = zeros(Float64,nd,3,n)
        else
            pot = zeros(Float64,3,n)
        end
    end

    if ppreg == 2 || ppreg == 3
        if nd > 1
            pre = zeros(Float64,nd,n)
        else
            pre = zeros(Float64,n)
        end
    end

    if ppreg == 3
        if nd > 1
            grad = zeros(Float64,nd,3,3,n)
        else
            grad = zeros(Float64,3,3,n)
        end
    end

    if ppregt == 1 || ppregt == 2 || ppregt == 3
        if nd > 1
            pottarg = zeros(Float64,nd,3,nt)
        else
            pottarg = zeros(Float64,3,nt)
        end
    end

    if ppregt == 2 || ppregt == 3
        if nd > 1
            pretarg = zeros(Float64,nd,nt)
        else
            pretarg = zeros(Float64,nt)
        end
    end

    if ppregt == 3
        if nd > 1
            gradtarg = zeros(Float64,nd,3,3,nt)
        else
            gradtarg = zeros(Float64,3,3,nt)
        end
    end

    # actually call the function
    #
    # fortran calling sequence
    #  subroutine stfmm3d(nd, eps,
    # $                 nsource, source,
    # $                 ifstoklet, stoklet, ifstrslet, strslet, strsvec,
    # $                 ifppreg, pot, pre, grad, ntarg, targ,
    # $                 ifppregtarg, pottarg, pretarg, gradtarg,ier)

    ccall((:stfmm3d_,libfmm3d),Cvoid,(Fi,Fd,Fi,Fd,Fi,Fd,
                                    Fi,Fd,Fd,Fi,Fd,Fd,Fd,Fi,
                                    Fd,Fi,Fd,Fd,Fd,Fi),
          nd,eps,n,sources,ifstoklet,stoklet,ifstrslet,
          strslet,strsvec,ppreg,pot,pre,grad,nt,targets,ppregt,
          pottarg,pretarg,gradtarg,ier)

    # copy over error message

    if (ier != 0); @warn "libfmm3d had an error, see vals.ier" end
    vals.ier = ier

    # load requested values

    if ppreg == 1 || ppreg == 2 || ppreg == 3; vals.pot = pot end
    if ppreg == 2 || ppreg == 3; vals.pre = pre end
    if ppreg == 3; vals.grad = grad end
    if ppregt == 1 || ppregt == 2 || ppregt == 3; vals.pottarg = pottarg end
    if ppregt == 2 || ppregt == 3; vals.pretarg = pretarg end
    if ppregt == 3; vals.gradtarg = gradtarg end

    return vals

end

"""
```julia
    vals = st3ddir(sources,targets;stoklet=nothing,strslet=nothing,
                   strsvec=nothing,ppregt=0,nd=1,thresh=1e-16)
```

This function computes the N-body Stokes interactions
in three dimensions where the interaction kernels are
the Stokeslet and stresslet (see below). This is the
``O(N^2)`` direct evaluation code. By convention this code only computes
the effect of sources on targets. If the value at sources is
also needed, the routine can be called again with targets equal
to the source locations.

We take the following conventions for the Stokes kernels:

For a source ``y`` and target ``x``, let ``r_i = x_i-y_i``
 and let ``r = \\sqrt{r_1^2 + r_2^2 + r_3^2}``

The Stokeslet, ``G_{ij}``, and its associated pressure tensor,
``P_j``, (without the ``1/4\\pi`` scaling) are

```math
G_{ij}(x,y) = (r_i r_j)/(2r^3) + \\delta_{ij}/(2r) \\; ,
\\quad
P_j(x,y) = r_j/r^3
```
The (Type I) stresslet, ``T_{ijk}``, and its associated
pressure tensor, ``\\Pi_{jk}``, (without the ``1/4\\pi``
 scaling) are

```math
 T_{ijk}(x,y) = -3 r_i r_j r_k/ r^5 \\; , \\quad
 PI_{jk} = -2 \\delta_{jk} + 6 r_j r_k/r^5
```

The output of this routine gives the velocity
```math
 u_i(x) = \\sum_{m=1}^{N} G_{ij}(x,x_m) \\sigma^{(m)}_j
    + T_{ijk}(x,x_m) \\mu^{(m)}_j \\nu^{(m)}_k \\, ,
```
and, optionally, the pressure

```math
 p(x) = \\sum_{m=1}^{N} P_{j}(x,x_m) \\sigma^{(m)}_j
    + \\Pi{jk}(x,x_m) \\mu^{(m)}_j \\nu^{(m)}_k \\, ,
```
where ``\\sigma^{(m)}`` is a Stokeslet density 3-vector,
``\\mu^{(m)}`` and ``\\nu^{(m)}`` are the stresslet
density and orientation 3-vectors, and
``x_{m}`` are the source locations.
When ``x=x_{m}``, the term corresponding to
``x_{m}`` is dropped from the sum.
The gradient of the velocity can also be computed by request.

# Input

* `sources::Array{Float64}` size (3,n) source locations (``x_{m}``)
* `targets::Array{Float64}` size (3,nt) target locations (``x``)

# Optional Keyword Input

* `stoklet::Array{Float64}` size (nd,3,n) or (3,n) Stokeslet densities (``\\sigma``)
* `strslet::Array{Float64}` size (nd,3,n) or (3,n) Stresslet densities (``\\mu'')  if either strsvec or strslet is provided, the other must be as well
* `strsvec::Array{Float64}` size (nd,3,n) or (3,n) Stresslet orientations (``\\nu``) if either strsvec or strslet is provided, the other must be as well
* `ppregt::Integer` source eval flag.
    + Velocity (``u``) at targets evaluated if `ppregt == 1`.
    + Velocity and pressure (``p``) at targets evaluated if `ppregt == 2`
    + Velocity, pressure, and velocity gradient (``\\nabla u``)
    at targets evaluated if `ppregt == 3`
* `nd::Integer` number of densities of each type

Note: if all default values are used for optional input, nothing
is computed.

# Output

`vals::FMMVals` with the fields

* `vals.pottarg::Array{Float64}` size (nd,3,nt) or (3,nt) velocity at target locations if requested
* `vals.gradtarg::Array{Float64}` size (nd,3,3,nt) or (3,nt) gradient of velocity at target locations if requested.
`gradtarg[l,i,j,k]` is the i-th component of the gradient of the j-th component of the velocity
for the l-th density at the k-th target location
* `vals.pretarg::Array{Float64}` size (nd,nt) or (nt) pressure at target locations if requested
"""
function st3ddir(sources::Array{Float64},targets::Array{Float64};
                 stoklet::TFN=nothing,strslet::TFN=nothing,
                 strsvec::TFN=nothing,
                 ppregt::Integer=0,nd::Integer=1,thresh=1e-16)

    # default values

    vals = FMMVals()

    ifstoklet = 0
    ifstrslet = 0

    zero = Float64(0)
    zero3 = zeros(Float64,3)

    pottarg = zero
    gradtarg = zero
    pretarg = zero

    ier = Integer(0)

    # check inputs

    ppreg=0
    anyfail, n, nt, ifstoklet, ifstrslet = (
        stfmm3dinputcheck(sources,stoklet,strslet,strsvec,targets,ppreg,ppregt,nd))

    if anyfail
        return vals
    end

    if (ifstoklet == 0); stoklet = zero3 end
    if (ifstrslet == 0); strslet = zero3 end
    if (nt == 0); targets = zero3 end

    # allocate memory for return values

    # for now we always allocate all quantities
    # and return requested

    if nd > 1
        pottarg = zeros(Float64,nd,3,nt)
    else
        pottarg = zeros(Float64,3,nt)
    end

    if nd > 1
        pretarg = zeros(Float64,nd,nt)
    else
        pretarg = zeros(Float64,nt)
    end

    if nd > 1
        gradtarg = zeros(Float64,nd,3,3,nt)
    else
        gradtarg = zeros(Float64,3,3,nt)
    end

    # dispatch to appropriate wrapper

    if ifstoklet == 1
        if ifstrslet == 1
            istress = 1
            st3ddirectstokstrsg!(nd,sources,stoklet,
                                 istress,strslet,strsvec,
                                 n,targets,nt,pottarg,
                                 pretarg,gradtarg,thresh)
        else
            st3ddirectstokg!(nd,sources,stoklet,
                             n,targets,nt,pottarg,
                             pretarg,gradtarg,thresh)
        end
    else
        if ifstrslet == 1
            istress = 1
            stoklet = zeros(Float64,nd,3,n)
            st3ddirectstokstrsg!(nd,sources,stoklet,
                             istress,strslet,strsvec,
                             n,targets,nt,pottarg,
                             pretarg,gradtarg,thresh)
        end
    end

    if ppregt == 1 || ppregt == 2 || ppregt == 3; vals.pottarg = pottarg end
    if ppregt == 2 || ppregt == 3; vals.pretarg = pretarg end
    if ppregt == 3; vals.gradtarg = gradtarg end

    return vals

end

function st3ddirectstokstrsg!(nd,sources,stoklet,istress,
                              strslet,strsvec,n,targets,nt,pottarg,
                              pretarg,gradtarg,thresh)
    ccall((:st3ddirectstokstrsg_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,Fd,Fd,Fi,
                                                  Fd,Fi,Fd,Fd,Fd,Fd),
          nd,sources,stoklet,istress,
                              strslet,strsvec,n,targets,nt,pottarg,
          pretarg,gradtarg,thresh)
    return
end

function st3ddirectstokg!(nd,sources,stoklet,n,targets,nt,pottarg,
                              pretarg,gradtarg,thresh)
    ccall((:st3ddirectstokg_,libfmm3d),Cvoid,(Fi,Fd,Fd,Fi,Fd,Fi,Fd,
                                              Fd,Fd,Fd),
          nd,sources,stoklet,n,targets,nt,pottarg,
          pretarg,gradtarg,thresh)
    return
end

"""
```julia
    anyfail, n, nt, ifA, ifB, iflam = (
        emfmm3dinputcheck(sources,A,B,lambda,
                        targets,ifE,ifdivE,ifcurlE,
                        ifEtarg,ifdivEtarg,
                        ifcurlEtarg,nd) )
```

Input checking for electromagnetics routines.
Checks the sizes of these arrays and the compatibility
of the various flags, nd, and provided arrays.
If something is off, a warning is issued and
`anyfail` is set to true. Non-fatal mistakes result in
a warning but `anyfail` remains false.

Output:
* `anyfail` - boolean is true if any checks fail, false otherwise
* `n` - integer, number of sources
* `nt` - integer, number of targets
* `ifA` - integer, is 1 if there is an A vector, 0 otherwise
* `ifB` - integer, is 1 if there is a B vector, 0 otherwise
* `iflam` - integer, is 1 if there is a lambda density, 0 otherwise
"""
function emfmm3dinputcheck(sources,A,B,lambda,targets,
                           ifE,ifdivE,ifcurlE,ifEtarg,ifdivEtarg,
                           ifcurlEtarg,nd)
    anyfail = false

    if (size(sources,1) != 3)
        @warn "sources array is wrong shape, computing nothing"
        anyfail = true
    end
    if (nd < 0)
        @warn "nd has invalid value, computing nothing"
        anyfail = true
    end

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

    if targets != nothing
        if (size(targets,1) != 3)
            @warn "targets array is wrong shape, computing nothing"
            anyfail = true
        end
        nt = div(length(targets),3)
    end

    if A != nothing
        if (div(length(A),nd) != 3*n)
            @warn "size of A vector array incompatible with sources array and nd parameter, computing nothing"
            anyfail = true
        end
        ifA = 1
    else
        ifA = 0
    end

    if B != nothing
        if (div(length(B),nd) != 3*n)
            @warn "size of B vector array incompatible with sources array and nd parameter, computing nothing"
            anyfail = true
        end
        ifB = 1
    else
        ifB = 0
    end

    if lambda != nothing
        if (div(length(lambda),nd) != n)
            @warn "size of lambda array incompatible with sources array and nd parameter, computing nothing"
            anyfail = true
        end
        iflam = 1
    else
        iflam = 0
    end

    if ifA == 0 && ifB == 0 && iflam == 0
        @warn "none of A vectors, B vectors, or lambdas provided at sources, doing nothing"
        anyfail = true
    end

    if !ifEtarg && !ifdivEtarg && !ifcurlEtarg
        @warn "no output requested, doing nothing"
        anyfail = true
    end

    if (ifEtarg || ifdivEtarg || ifcurlEtarg) && targets == nothing
        @warn "target values requested but no targets provided, doing nothing for these"
    end

    return anyfail, n, nt, ifA, ifB, iflam
end

"""
```julia
    vals = emfmm3d(eps,zk,sources;A=nothing,B=nothing,lambda=nothing,
                ifE=false,ifdivE=false,ifcurlE=false,
                ifEtarg=false,ifdivEtarg=false,ifcurlEtarg=false,
                nd=1,targets=nothing)
```

This function computes N-body electromagnetic interactions
in three dimensions where the interaction kernels are
the Helmholtz kernel and its curl and gradient (see below). This is the
``O(N)`` fast multipole code which computes the interactions
to the requested precision.

The Helmholtz Green's function (without the ``1/4 \\pi`` scaling)
is
```math
 G_k(x,y) = \\frac{e^{ik \\|x-y\\|}}{\\|x-y\\|} \\, ,
```

This routine computes the sum

```math
E(x) = \\sum_{j=1}^m \\nabla \\times (G_k(x,x_{j}) A_j) +
G_k(x,x_{j}) B_j + \\nabla G_k(x,x_{j}) \\lambda_j
```

where ``A_{j}`` and ``B_{j}`` are vector densities
      ``\\lambda_{j}`` is a scalar density, 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_{m}``)

# Optional Keyword Input

* `targets::Array{Float64}` size (3,nt) target locations (``x``)
* `A::Array{ComplexF64}` size (nd,3,n) or (3,n) vector densities
* `B::Array{ComplexF64}` size (nd,3,n) or (3,n) vector densities
* `lambda::Array{ComplexF64}` size (nd,n) or (n) scalar densities
* `ifE::Bool` E field evaluation flag, if true evaluate E at source points.
* `ifdivE::Bool` divergence of E field evaluation flag, if true evaluate divergence of E at source points.
* `ifcurlE::Bool` curl of E field evaluation flag, if true evaluate curl of E at source points.
* `ifEtarg::Bool` E field evaluation flag, if true evaluate E at target points.
* `ifdivEtarg::Bool` divergence of E field evaluation flag, if true evaluate divergence of E at target points.
* `ifcurlEtarg::Bool` curl of E field evaluation flag, if true evaluate curl of E at target points.
* `nd::Integer` number of densities of each type

Note: if all default values are used for optional input, nothing
is computed.

# Output

`vals::FMMVals` with the fields

* `vals.E::Array{ComplexF64}` size (nd,3,n) or (3,n) E field at source locations if requested
* `vals.divE::Array{ComplexF64}` size (nd,n) or (n) divergence of E field at source locations if requested
* `vals.curlE::Array{ComplexF64}` size (nd,3,n) or (3,n) curl of E field at source locations if requested
* `vals.Etarg::Array{ComplexF64}` size (nd,3,nt) or (3,nt) E field at target locations if requested
* `vals.divEtarg::Array{ComplexF64}` size (nd,nt) or (nt) divergence of E field at target locations if requested
* `vals.curlEtarg::Array{ComplexF64}` size (nd,3,nt) or (3,nt) curl of E field at target locations if requested
* `vals.ier` error flag as returned by FMM3D library. A value of 0 indicates a successful call.
Non-zero values may indicate insufficient memory available. See the documentation for the FMM3D library.
If not set (`nothing`), then FMM3D library was never called.
"""
function emfmm3d(eps::Float64,zk::Union{Float64,ComplexF64},sources::Array{Float64};
                 A::TCN=nothing,B::TCN=nothing,lambda::TCN=nothing,
                 ifE::Bool=false,ifdivE::Bool=false,
                 ifcurlE::Bool=false,ifEtarg::Bool=false,
                 ifdivEtarg::Bool=false,ifcurlEtarg::Bool=false,
                 nd::Integer=1,targets::TFN=nothing)

    zk = complex(zk)

    # default values

    vals = FMMVals()

    zero = ComplexF64(0)
    zero3 = zeros(ComplexF64,3)

    E1 = zero3
    divE1 = zero
    curlE1 = zero3

    ier = Integer(0)

    # check inputs

    anyfail, n, nt, ifA, ifB, iflam = (
        emfmm3dinputcheck(sources,A,B,lambda,
                          targets,ifE,ifdivE,ifcurlE,
                          ifEtarg,ifdivEtarg,
                          ifcurlEtarg,nd) )

    if anyfail
        return vals
    end

    if (ifA == 0); A = zero3 end
    if (ifB == 0); B = zero3 end
    if (iflam == 0); lambda = zero end
    if (nt == 0); targets = zero3 end

    # allocate memory for return values

    anysrc = (ifE || ifdivE || ifcurlE)
    anytarg = (ifEtarg || ifdivEtarg || ifcurlEtarg)

    ntot = 0
    if anysrc; ntot = ntot+n; end
    if anytarg; ntot = ntot+nt; end

    targ1 = zeros(Float64,3,ntot)
    itstart = 1
    itend = nt
    if anysrc
        targ1[:,1:n] = sources
        itstart = itstart+n
        itend = itend+n
    end
    if anytarg
        targ1[:,itstart:itend] = targets
    end

    anyE = ifE || ifEtarg
    anydivE = ifdivE || ifdivEtarg
    anycurlE = ifcurlE || ifcurlEtarg

    ifE1 = 0; ifdivE1 = 0; ifcurlE1 = 0

    if anyE
        if nd > 1
            E1 = zeros(ComplexF64,nd,3,ntot)
        else
            E1 = zeros(ComplexF64,3,ntot)
        end
        ifE1 = 1
    end

    if anydivE
        if nd > 1
            divE1 = zeros(ComplexF64,nd,ntot)
        else
            divE1 = zeros(ComplexF64,ntot)
        end
        ifdivE1 = 1
    end

    if anyE
        if nd > 1
            curlE1 = zeros(ComplexF64,nd,3,ntot)
        else
            curlE1 = zeros(ComplexF64,3,ntot)
        end
        ifcurlE1 = 1
    end

# actually call the function
#
# fortran calling sequence
#  subroutine emfmm3d(nd,eps,zk,ns,source,ifa_vect,a_vect,&
#   ifb_vect,b_vect,iflambda,lambda,nt,targets,ifE,E,ifcurlE,curlE,&
#   ifdivE,divE,ier)

    ccall((:emfmm3d_,libfmm3d),Cvoid,(Fi,Fd,Fc,Fi,Fd,Fi,Fc,
                                      Fi,Fc,Fi,Fc,Fi,Fd,Fi,Fc,
                                      Fi,Fc,Fi,Fc,Fi),
          nd,eps,zk,n,sources,ifA,A,ifB,B,iflam,lambda,
          ntot,targ1,ifE1,E1,ifcurlE1,curlE1,ifdivE1,divE1,ier)

    # copy over error message

    if (ier != 0); @warn "libfmm3d had an error, see vals.ier" end
    vals.ier = ier

    # load requested values

    if (nd == 1)
        if ifE; vals.E = E1[:,1:n] end
        if ifdivE; vals.divE = divE1[1:n] end
        if ifcurlE; vals.curlE = curlE1[:,1:n] end

        if ifEtarg; vals.Etarg = E1[:,itstart:itend] end
        if ifdivEtarg; vals.divEtarg = divE1[itstart:itend] end
        if ifcurlEtarg; vals.curlEtarg = curlE1[:,itstart:itend] end
    else
        if ifE; vals.E = E1[:,:,1:n] end
        if ifdivE; vals.divE = divE1[:,1:n] end
        if ifcurlE; vals.curlE = curlE1[:,:,1:n] end

        if ifEtarg; vals.Etarg = E1[:,:,itstart:itend] end
        if ifdivEtarg; vals.divEtarg = divE1[:,itstart:itend] end
        if ifcurlEtarg; vals.curlEtarg = curlE1[:,:,itstart:itend] end
    end
    return vals

end

"""
```julia
    vals = em3ddir(zk,sources,targets;A=nothing,B=nothing,lambda=nothing,
                ifEtarg=false,ifdivEtarg=false,ifcurlEtarg=false,
                nd=1,thresh=1e-16)
```

This function computes N-body electromagnetic interactions
in three dimensions where the interaction kernels are
the Helmholtz kernel and its curl and gradient (see below). This is the
``O(N^2)`` direct code. By convention this code only computes
the effect of sources on targets. If the value at sources is
also needed, the routine can be called again with targets equal
to the source locations.

The Helmholtz Green's function (without the ``1/4 \\pi`` scaling)
is
```math
 G_k(x,y) = \\frac{e^{ik \\|x-y\\|}}{\\|x-y\\|} \\, ,
```

This routine computes the sum

```math
E(x) = \\sum_{j=1}^m \\nabla \\times (G_k(x,x_{j}) A_j) +
G_k(x,x_{j}) B_j + \\nabla G_k(x,x_{j}) \\lambda_j
```

where ``A_{j}`` and ``B_{j}`` are vector densities
      ``\\lambda_{j}`` is a scalar density, 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_{m}``)
* `targets::Array{Float64}` size (3,nt) target locations (``x``)

# Optional Keyword Input

* `A::Array{ComplexF64}` size (nd,3,n) or (3,n) vector densities
* `B::Array{ComplexF64}` size (nd,3,n) or (3,n) vector densities
* `lambda::Array{ComplexF64}` size (nd,n) or (n) scalar densities
* `ifEtarg::Bool` E field evaluation flag, if true evaluate E at target points.
* `ifdivEtarg::Bool` divergence of E field evaluation flag, if true evaluate divergence of E at target points.
* `ifcurlEtarg::Bool` curl of E field evaluation flag, if true evaluate curl of E at target points.
* `nd::Integer` number of densities of each type

Note: if all default values are used for optional input, nothing
is computed.

# Output

`vals::FMMVals` with the fields

* `vals.Etarg::Array{ComplexF64}` size (nd,3,nt) or (3,nt) E field at target locations if requested
* `vals.divEtarg::Array{ComplexF64}` size (nd,nt) or (nt) divergence of E field at target locations if requested
* `vals.curlEtarg::Array{ComplexF64}` size (nd,3,nt) or (3,nt) curl of E field at target locations if requested
"""
function em3ddir(zk::Union{Float64,ComplexF64},sources::Array{Float64},
                 targets::Array{Float64};
                 A::TCN=nothing,B::TCN=nothing,lambda::TCN=nothing,
                 ifEtarg::Bool=false,
                 ifdivEtarg::Bool=false,ifcurlEtarg::Bool=false,
                 nd::Integer=1,thresh::Float64=1e-16)

    zk = complex(zk)

    # default values

    vals = FMMVals()

    zero = ComplexF64(0)
    zero3 = zeros(ComplexF64,3)

    Etarg = zero3
    divEtarg = zero
    curlEtarg = zero3

    ier = Integer(0)

    # check inputs

    ifE = false; ifdivE = false; ifcurlE = false

    anyfail, n, nt, ifA, ifB, iflam = (
        emfmm3dinputcheck(sources,A,B,lambda,
                          targets,ifE,ifdivE,ifcurlE,
                          ifEtarg,ifdivEtarg,
                          ifcurlEtarg,nd) )
    if anyfail
        return vals
    end

    if (ifA == 0); A = zero3 end
    if (ifB == 0); B = zero3 end
    if (iflam == 0); lambda = zero end
    if (nt == 0); targets = zero3 end

    if (nt == 0)
        return vals
    end

    # allocate memory for return values

    ifE1 = Integer(0); ifcurlE1 = Integer(0)
    ifdivE1 = Integer(0)

    if ifEtarg
        if nd > 1
            Etarg = zeros(ComplexF64,nd,3,nt)
        else
            Etarg = zeros(ComplexF64,3,nt)
        end
        ifE1 = 1
    end

    if ifcurlEtarg
        if nd > 1
            curlEtarg = zeros(ComplexF64,nd,3,nt)
        else
            curlEtarg = zeros(ComplexF64,3,nt)
        end
        ifcurlE1 = 1
    end

    if ifdivEtarg
        if nd > 1
            divEtarg = zeros(ComplexF64,nd,nt)
        else
            divEtarg = zeros(ComplexF64,nt)
        end
        ifdivE1 = 1
    end

# actually call routine
#
# subroutine em3ddirect(nd,zk,ns,source,ifa_vect,a_vect,&
#  ifb_vect,b_vect,iflambda,lambda,nt,targets,ifE,E,ifcurlE,curlE,&
#  ifdivE,divE,thresh)

    ccall((:em3ddirect_,libfmm3d),Cvoid,(Fi,Fc,Fi,Fd,Fi,Fc,
                                         Fi,Fc,Fi,Fc,Fi,Fd,
                                         Fi,Fc,Fi,Fc,Fi,Fc,
                                         Fd),
          nd,zk,n,sources,ifA,A,ifB,B,iflam,lambda,nt,targets,
          ifE1,Etarg,ifcurlE1,curlEtarg,ifdivE1,divEtarg,thresh)

    if ifEtarg; vals.Etarg = Etarg end
    if ifdivEtarg; vals.divEtarg = divEtarg end
    if ifcurlEtarg; vals.curlEtarg = curlEtarg end

    return vals

end

end # module


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

"""
Traverse the entries in a tree and its subtrees in post or pre order.

Function parameter should have following signature:

    (Cstring, Ptr{Void}, Ptr{Void}) -> Cint
"""
function treewalk(f::Function, tree::GitTree, payload=Any[], post::Bool = false)
    cbf = cfunction(f, Cint, (Cstring, Ptr{Void}, Ptr{Void}))
    cbf_payload = Ref{typeof(payload)}(payload)
    @check ccall((:git_tree_walk, :libgit2), Cint,
                  (Ptr{Void}, Cint, Ptr{Void}, Ptr{Void}),
                   tree.ptr, post, cbf, cbf_payload)
    return cbf_payload
end

function repository(tree::GitTree)
    repo_ptr = ccall((:git_tree_owner, :libgit2), Ptr{Void},
                 (Ptr{Void},), tree.ptr)
    return GitRepo(repo_ptr)
end

function filename(te::GitTreeEntry)
    str = ccall((:git_tree_entry_name, :libgit2), Cstring, (Ptr{Void},), te.ptr)
    str != C_NULL && return unsafe_string(str)
    return nothing
end

function filemode(te::GitTreeEntry)
    return ccall((:git_tree_entry_filemode, :libgit2), Cint, (Ptr{Void},), te.ptr)
end

function entrytype(te::GitTreeEntry)
    otype = ccall((:git_tree_entry_type, :libgit2), Cint, (Ptr{Void},), te.ptr)
    return getobjecttype(otype)
end

function entryid(te::GitTreeEntry)
    oid_ptr = ccall((:git_tree_entry_id, :libgit2), Cint, (Ptr{Void},), te.ptr)
    return GitHash(oid_ptr[])
end

function (::Type{T}){T<:GitObject}(repo::GitRepo, te::GitTreeEntry)
    obj_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_tree_entry_to_object, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ref{Void}),
                   obj_ptr_ptr, repo.ptr, te.ptr)
    return T(repo, obj_ptr_ptr[])
end

function Base.show(io::IO, te::GitTreeEntry)
    println(io, "GitTreeEntry:")
    println(io, "Entry name: ", filename(te))
    println(io, "Entry type: ", entrytype(te))
    println(io, "Entry OID: ", entryid(te))
end

ext: jl


text:
"""
Determines if of Image Type
"""
function isa_image(x::Type{T}) where T<:Matrix
    eltype(T) <: Union{Colorant, Colors.Fractional}
end
isa_image(x::Matrix) = isa_image(typeof(x))
# if !isdefined(Images, :ImageAxes)
#     include_string("""
#     isa_image(x::Images.Image) = true
#     """)
# end
isa_image(x) = false

# Splits a dictionary in two dicts, via a condition
function Base.split(condition::Function, associative::Associative)
    A = similar(associative)
    B = similar(associative)
    for (key, value) in associative
        if condition(key, value)
            A[key] = value
        else
            B[key] = value
        end
    end
    A, B
end

function assemble_robj(data, program, bb, primitive, pre_fun, post_fun)
    pre = if pre_fun != nothing
        () -> (GLAbstraction.StandardPrerender(); pre_fun())
    else
        GLAbstraction.StandardPrerender()
    end
    robj = RenderObject(data, program, pre, nothing, bb, nothing)
    post = if haskey(data, :instances)
        GLAbstraction.StandardPostrenderInstanced(data[:instances], robj.vertexarray, primitive)
    else
        GLAbstraction.StandardPostrender(robj.vertexarray, primitive)
    end
    robj.postrenderfunction = if post_fun != nothing
        () -> begin
            post()
            post_fun()
        end
    else
        post
    end
    robj
end

function assemble_shader(data)
    shader = data[:shader]
    delete!(data, :shader)
    default_bb = Signal(centered(AABB))
    bb  = get(data, :boundingbox, default_bb)
    if bb == nothing || isa(bb, Signal{Void})
        bb = default_bb
    end
    glp = get(data, :gl_primitive, GL_TRIANGLES)
    robj = assemble_robj(
        data, shader, bb, glp,
        get(data, :prerender, nothing),
        get(data, :postrender, nothing)
    )
    Context(robj)
end

function y_partition_abs(area, amount)
    a = round(Int, amount)
    p = const_lift(area) do r
        (
            SimpleRectangle{Int}(0, 0, r.w, a),
            SimpleRectangle{Int}(0, a, r.w, r.h - a)
        )
    end
    return map(first, p), map(last, p)
end
function x_partition_abs(area, amount)
    a = round(Int, amount)
    p = const_lift(area) do r
        (
            SimpleRectangle{Int}(0, 0, a, r.h),
            SimpleRectangle{Int}(a, 0, r.w - a, r.h)
        )
    end
    return map(first, p), map(last, p)
end

function y_partition(area, percent)
    amount = percent / 100.0
    p = const_lift(area) do r
        (
            SimpleRectangle{Int}(0, 0, r.w, round(Int, r.h*amount)),
            SimpleRectangle{Int}(0, round(Int, r.h*amount), r.w, round(Int, r.h*(1-amount)))
        )
    end
    return map(first, p), map(last, p)
end
function x_partition(area, percent)
    amount = percent / 100.0
    p = const_lift(area) do r
        (
            SimpleRectangle{Int}(0, 0, round(Int, r.w*amount), r.h ),
            SimpleRectangle{Int}(round(Int, r.w*amount), 0, round(Int, r.w*(1-amount)), r.h)
        )
    end
    return map(first, p), map(last, p)
end

glboundingbox(mini, maxi) = AABB{Float32}(Vec3f0(mini), Vec3f0(maxi)-Vec3f0(mini))
function default_boundingbox(main, model)
    main == nothing && return Signal(AABB{Float32}(Vec3f0(0), Vec3f0(1)))
    const_lift(*, model, AABB{Float32}(main))
end
AABB(a::GPUArray) = AABB{Float32}(gpu_data(a))
AABB{T}(a::GPUArray) where {T} = AABB{T}(gpu_data(a))

"""
Returns two signals, one boolean signal if clicked over `robj` and another
one that consists of the object clicked on and another argument indicating that it's the first click
"""
function clicked(robj::RenderObject, button::MouseButton, window::Screen)
    @materialize mouse_hover, mouse_buttons_pressed = window.inputs
    leftclicked = const_lift(mouse_hover, mouse_buttons_pressed) do mh, mbp
        mh[1] == robj.id && mbp == Int[button]
    end
    clicked_on_obj = keepwhen(leftclicked, false, leftclicked)
    clicked_on_obj = const_lift((mh, x)->(x,robj,mh), mouse_hover, leftclicked)
    leftclicked, clicked_on_obj
end

"""
Returns a boolean signal indicating if the mouse hovers over `robj`
"""
is_hovering(robj::RenderObject, window::Screen) =
    droprepeats(const_lift(is_same_id, mouse2id(window), robj))

function dragon_tmp(past, mh, mbp, mpos, robj, button, start_value)
    diff, dragstart_index, was_clicked, dragstart_pos = past
    over_obj = mh[1] == robj.id
    is_clicked = mbp == Int[button]
    if is_clicked && was_clicked # is draggin'
        return (dragstart_pos-mpos, dragstart_index, true, dragstart_pos)
    elseif over_obj && is_clicked && !was_clicked # drag started
        return (Vec2f0(0), mh[2], true, mpos)
    end
    return start_value
end

"""
Returns a signal with the difference from dragstart and current mouse position,
and the index from the current ROBJ id.
"""
function dragged_on(robj::RenderObject, button::MouseButton, window::Screen)
    @materialize mouse_buttons_pressed, mouseposition = window.inputs
    mousehover = mouse2id(window)
    mousedown = const_lift(GLAbstraction.singlepressed, mouse_buttons_pressed, GLFW.MOUSE_BUTTON_LEFT)
    condition = const_lift(is_same_id, mousehover, robj)
    dragg = GLAbstraction.dragged(mouseposition, mousedown, condition)
    filterwhen(mousedown, (value(dragg), 0), map(dragg) do d
        d, value(mousehover).index
    end)
end

points2f0(positions::Vector{T}, range::Range) where {T} = Point2f0[Point2f0(range[i], positions[i]) for i=1:length(range)]

extrema2f0(x::Array{T,N}) where {T<:Intensity,N} = Vec2f0(extrema(reinterpret(Float32,x)))
extrema2f0(x::Array{T,N}) where {T,N} = Vec2f0(extrema(x))
extrema2f0(x::GPUArray) = extrema2f0(gpu_data(x))
function extrema2f0(x::Array{T,N}) where {T<:Vec,N}
    _norm = map(norm, x)
    Vec2f0(minimum(_norm), maximum(_norm))
end

function mix_linearly(a::C, b::C, s) where C<:Colorant
    RGBA{Float32}((1-s)*comp1(a)+s*comp1(b), (1-s)*comp2(a)+s*comp2(b), (1-s)*comp3(a)+s*comp3(b), (1-s)*alpha(a)+s*alpha(b))
end

color_lookup(cmap, value, mi, ma) = color_lookup(cmap, value, (mi, ma))
function color_lookup(cmap, value, color_norm)
    mi,ma = color_norm
    scaled = clamp((value-mi)/(ma-mi), 0, 1)
    index = scaled * (length(cmap)-1)
    i_a, i_b = floor(Int, index)+1, ceil(Int, index)+1
    mix_linearly(cmap[i_a], cmap[i_b], scaled)
end

"""
Converts index arrays to the OpenGL equivalent.
"""
to_index_buffer(x::GLBuffer) = x
to_index_buffer(x::TOrSignal{Int}) = x
to_index_buffer(x::VecOrSignal{UnitRange{Int}}) = x
to_index_buffer(x::TOrSignal{UnitRange{Int}}) = x
"""
For integers, we transform it to 0 based indices
"""
to_index_buffer(x::Vector{I}) where {I<:Integer} = indexbuffer(map(i-> Cuint(i-1), x))
function to_index_buffer(x::Signal{Vector{I}}) where I<:Integer
    x = map(x-> Cuint[i-1 for i=x], x)
    gpu_mem = GLBuffer(value(x), buffertype = GL_ELEMENT_ARRAY_BUFFER)
    preserve(const_lift(update!, gpu_mem, x))
    gpu_mem
end
"""
If already GLuint, we assume its 0 based (bad heuristic, should better be solved with some Index type)
"""
to_index_buffer(x::Vector{I}) where {I<:GLuint} = indexbuffer(x)
function to_index_buffer(x::Signal{Vector{I}}) where I<:GLuint
    gpu_mem = GLBuffer(value(x), buffertype = GL_ELEMENT_ARRAY_BUFFER)
    preserve(const_lift(update!, gpu_mem, x))
    gpu_mem
end
function to_index_buffer(x::Signal{Vector{I}}) where I <: Face{2, GLIndex}
    gpu_mem = GLBuffer(value(x), buffertype = GL_ELEMENT_ARRAY_BUFFER)
    preserve(const_lift(update!, gpu_mem, x))
    gpu_mem
end
to_index_buffer(x) = error(
    "Not a valid index type: $(typeof(x)).
    Please choose from Int, Vector{UnitRange{Int}}, Vector{Int} or a signal of either of them"
)

"""
Creates a moving average and discards values to close together.
If discarded return (false, p), if smoothed, (true, smoothed_p).
"""
function moving_average(p, cutoff,  history, n = 5)
    if length(history) > 0
        if norm(p - history[end]) < cutoff
            return false, p # don't keep point
        end
    end
    if length(history) == 5
        # maybe better to just keep a current index
        history[1:5] = circshift(view(history, 1:5), -1)
        history[end] = p
    else
        push!(history, p)
    end
    true, sum(history) ./ length(history)# smooth
end

function layoutlinspace(n::Integer)
    if n == 1
        1:1
    else
        linspace(1/n, 1, n)
    end
end
xlayout(x::Int) = zip(layoutlinspace(x), Iterators.repeated(""))
function xlayout(x::AbstractVector{T}) where T <: AbstractFloat
    zip(x, Iterators.repeated(""))
end

function xlayout(x::AbstractVector)
    zip(layoutlinspace(length(x)), x)
end
function ylayout(x::AbstractVector)
    zip(layoutlinspace(length(x)), x)
end
function ylayout(x::AbstractVector{T}) where T <: Tuple
    sizes = map(first, x)
    values = map(last, x)
    zip(sizes, values)
end
function IRect(x, y , w, h)
    SimpleRectangle(
        round(Int, x),
        round(Int, y),
        round(Int, w),
        round(Int, h),
    )
end

function layout_rect(area, lastw, lasth, w, h)
    wp = widths(area)
    xmin = wp[1] * lastw
    ymin = wp[2] * lasth
    xmax = wp[1] * w
    ymax = wp[2] * h
    xmax = max(xmin, xmax)
    xmin = min(xmin, xmax)
    ymax = max(ymin, ymax)
    ymin = min(ymin, ymax)
    IRect(xmin, ymin, xmax - xmin, ymax - ymin)
end

function layoutscreens(parent, layout;
        title_height = 6mm, text_color = default(RGBA),
        background_color = RGBA(1f0, 1f0, 1f0), stroke_color = RGBA(0f0, 0f0, 0f0, 0.4f0),
        kw_args...
    )
    layout = reverse(layout) # we start from bottom to top, while lists are written top to bottom
    lastw, lasth = 0, 0
    result = Vector{Screen}[]
    for (h, xlist) in ylayout(layout)
        result_x = Screen[]
        for (w, title) in xlayout(xlist)
            area = const_lift(layout_rect, parent.area, lastw, lasth, w, h)

            screen = Screen(parent; area = area, kw_args...)
            push!(result_x, screen)
            if !isempty(value(title))
                tarea = map(area) do area
                    IRect(0, area.h - title_height, area.w, title_height)
                end
                title_screen = Screen(screen, area = tarea, color = background_color)
                robj = visualize(
                    title, relative_scale = title_height * 0.7,
                    direction = 1, gap = Vec3f0(1mm, 0, 0),# in case it's a list!
                    color = text_color
                )
                gap = title_height * 0.15
                GLAbstraction.transform!(robj, translationmatrix(Vec3f0(gap, gap, 0)))
                _view(robj, title_screen, camera = :fixed_pixel)
                if stroke_color != nothing
                    _view(visualize(
                        map(a-> Point2f0[(0, 0), (a.w, 0)], area), :linesegment,
                        thickness = 1f0, color = stroke_color
                    ), title_screen, camera = :fixed_pixel)
                end
            end
            lastw = w
        end
        lastw = 0; lasth = h
        push!(result, result_x)
    end
    result
end


text:
@testset "planck_spectra" begin

    ν_grid = collect((1e7/778):0.015:(1e7/755))
    T = 3777

    println("Verifying planck function correctness...")

    Ts = [290, 1000, 3777]
    peak_ν = [568.693, 1961.01, 7406.74]
    peak_λ = [9.99225, 2.89775, 0.76721]
    peak_L_ν = [138.636, 5684.38, 306284]
    peak_L_λ = [8.40098, 4095.81, 3.14829e6]

    for i in 1:length(Ts)

        T_curr = Ts[i]

        peak_ν_true = peak_ν[i]
        peak_L_true = peak_L_ν[i]

        peak_calc = findmax(planck_spectrum_wn(T_curr)[:,2])

        @test isapprox(peak_calc[1], peak_L_true, rtol=0.0001, atol=1)
        @test isapprox(peak_calc[2], peak_ν_true, rtol=0.0001, atol=1)
        @test isapprox(peak_L_λ[i], planck_spectrum_wl(T_curr, [peak_λ[i]])[1], rtol=0.0001, atol=1)

    end

end


text:
_tall(_, ::Type{Tuple{}}) = true
_tall(f, ::Type{Tuple{T}}) where {T} = f(T)
_tall(f, ::Type{T}) where {T<:Tuple} =
    f(Base.tuple_type_head(T)) && _tall(f, Base.tuple_type_tail(T))

_map(f, a::Tuple) = map(f, a)
_map(f, a::NamedTuple{names}) where {names} = NamedTuple{names}(_map(f, Tuple(a)))

_map(f, (x,)::Tuple{Any}, (y,)::Tuple{Any}) = (f(x, y),)
_map(f, a::Tuple, b::Tuple) = (f(a[1], b[1]), _map(f, Base.tail(a), Base.tail(b))...)
_map(f, a::NamedTuple{names}, b::NamedTuple{names}) where {names} =
    NamedTuple{names}(_map(f, Tuple(a), Tuple(b)))

struct DisallowedOperation <: Exception
    f
    x
end

function Base.showerror(io::IO, err::DisallowedOperation)
    print(io, "$(err.f) on ")
    Base.showarg(io, err.x, false)
    print(io, " not allowed")
end


text:
"""
rollout.jl
"""

show_nice(x) = show(IOContext(stdout), "text/plain", x)

function rollout!(solver::iLQRSolver{IR}, α) where {IR}
    N = solver.N
    X = solver.X
    U = solver.U
    X_tmp = solver.X_tmp
    U_tmp = solver.U_tmp
    δx = solver.X_tmp[N + 1]
    δu = solver.U_tmp[N]
    ts = solver.ts
    K = solver.K
    d = solver.d
    J = 0.

    X_tmp[1] .= X[1]
    for k = 1:N-1
        # handle feedforward
        dt = ts[k + 1] - ts[k]
        state_diff!(δx, solver.model, X_tmp[k], X[k])
	δu .= d[k]
        δu .*= α
	mul!(δu, K[k], δx, 1., 1.)
        U_tmp[k] .= U[k]
        U_tmp[k] .+= δu
        discrete_dynamics!(X_tmp[k + 1], IR, solver.model, X_tmp[k],
                              U_tmp[k], ts[k], dt)
        # compute cost
        J += cost(solver.obj, X_tmp, U_tmp, k)
        max_x = norm(X_tmp[k + 1], Inf)
        if max_x > solver.opts.max_state_value || isnan(max_x)
            solver.stats.status = STATE_LIMIT
            return 0., true
        end
        max_u = norm(U_tmp[k], Inf)
        if max_u > solver.opts.max_control_value || isnan(max_u)
            solver.stats.status = CONTROL_LIMIT
            return 0., true
        end
    end
    J += cost(solver.obj, X_tmp, U_tmp, N)
    solver.stats.status = UNSOLVED
    return J, false
end


text:
# if no seed, each task will have different result

using MPI
using Random

function runc()
    MPI.Init()

    comm = MPI.COMM_WORLD
    my_rank = MPI.Comm_rank(comm)
    Random.seed!(1234)

    s = randn(10_000)
    @show s[1:10]

    MPI.Finalize()
end

runc()

# In uneven  randn(),1.why rank0 different: because alpha=randn() already use seed.
# 2.why rank_last is different: because rand(my_rank_size) use more seed than others


text:
########################################################################################
# Post-Processing: Smoothing and uncertainty calibration
########################################################################################
function calibrate!(integ::ODEFilterIntegrator)

    @unpack state_estimates, sigmas = integ
    for s in state_estimates
        s.Σ .*= sigmas[end]
    end
end

function smooth_all!(integ::ODEFilterIntegrator)

    @unpack state_estimates, times, sigmas = integ
    @unpack A!, Q!, Precond, InvPrecond = integ.constants
    @unpack Ah, Qh = integ.cache
    # x_pred is just used as a cache here

    for i in length(state_estimates)-1:-1:2
        dt = times[i+1] - times[i]

        P = Precond(dt)
        PI = InvPrecond(dt)

        A!(Ah, dt)
        Q!(Qh, dt)
        Qh .*= sigmas[i]

        state_estimates[i] = P * state_estimates[i]
        smooth!(state_estimates[i], P*state_estimates[i+1], Ah, Qh, integ, PI)
        any(isnan.(state_estimates[i].μ)) && error("NaN mean after smoothing")
        any(isnan.(state_estimates[i].Σ)) && error("NaN cov after smoothing")
        state_estimates[i] = PI * state_estimates[i]
    end
end

function smooth!(x_curr, x_next, Ah, Qh, integ, PI=I)

    @unpack d, q = integ.constants
    @unpack x_tmp = integ.cache

    if all((Qh) .< eps(eltype(Qh)))
        @warn "smooth: Qh is really small! The system is basically deterministic, so we just \"predict backwards\"."
        return inv(Ah) * x_next
    end

    # Prediction: t -> t+1
    mul!(x_tmp.μ, Ah, x_curr.μ)
    x_tmp.Σ .= Ah * x_curr.Σ * Ah' .+ Qh

    # Smoothing
    cov_before = copy(x_curr.Σ)
    cov_pred = copy(x_tmp.Σ)
    P_p = Symmetric(cov_pred)
    P_p_inv = inv(P_p)
    G = x_curr.Σ * Ah' * P_p_inv
    x_curr.μ .+= G * (x_next.μ .- x_tmp.μ)

    # Vanilla:
    cov_diff = x_next.Σ .- x_tmp.Σ
    GDG = G * cov_diff * G'
    x_tmp.Σ .= x_curr.Σ .+ GDG

    assert_nonnegative_diagonal(x_curr.Σ)

    return nothing
end

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
# Column-based operations: Some operations on rows are faster when considering columns

# In `map`, the output shouldn't alias inputs, so copies are made
Base.map(::typeof(identity), t::Table) = copy(t)

Base.map(::typeof(merge), t::Table) = copy(t)

function Base.map(::typeof(merge), t1::Table, t2::Table)
    return copy(Table(merge(columns(t1), columns(t2))))
end

function Base.map(f::GetProperty, t::Table)
    return copy(f(t))
end

@inline function Base.map(f::GetProperties, t::Table)
    return copy(f(t))
end

@inline function Base.map(f::Compute{names}, t::Table) where {names}
    # minimize number of columns before iterating over the rows
    map(f, GetProperties(names)(t))
end

@inline function Base.map(f::Compute{names}, t::Table{<:NamedTuple{names}}) where {names}
    # efficient to iterate over rows with a minimal number of columns
    if length(names) == 1 # unwrap in the simple cases
        return map(f.f, getproperty(names[1])(t))
    elseif length(names) == 2
        return map(f.f, getproperty(names[1])(t), getproperty(names[2])(t))
    end

    invoke(map, Tuple{Function, typeof(t)}, f, t)
end

@generated function Base.map(s::Select{names}, t::Table) where {names}
    exprs = [:($(names[i]) = map(s.fs[$i], t)) for i in 1:length(names)]

    return :(Table($(Expr(:tuple, exprs...))))
end

# In `mapview`, the output should alias the inputs
SplitApplyCombine.mapview(::typeof(merge), t::Table) = t

function SplitApplyCombine.mapview(::typeof(merge), t1::Table, t2::Table)
    return Table(merge(columns(t1), columns(t2)))
end

@inline function SplitApplyCombine.mapview(f::GetProperty, t::Table)
    return f(t)
end

@inline function SplitApplyCombine.mapview(f::GetProperties, t::Table)
    return f(t)
end

@inline function SplitApplyCombine.mapview(f::Compute{names}, t::Table) where {names}
    # minimize number of columns before iterating over the rows
    mapview(f, GetProperties(names)(t))
end

@inline function SplitApplyCombine.mapview(f::Compute{names}, t::Table{<:NamedTuple{names}}) where {names}
    # efficient to iterate over rows with a minimal number of columns
    if length(names) == 1 # unwrap in the simple cases (consider 2-argument version)
        return mapview(f.f, getproperty(names[1])(t))
    end

    invoke(mapview, Tuple{Function, typeof(t)}, f, t)
end

@generated function SplitApplyCombine.mapview(s::Select{names}, t::Table) where {names}
    exprs = [:($(names[i]) = mapview(s.fs[$i], t)) for i in 1:length(names)]

    return :(Table($(Expr(:tuple, exprs...))))
end

# broadcast

@inline function Broadcast.broadcasted(::Broadcast.DefaultArrayStyle{N}, ::typeof(merge), ts::Table{<:Any, N}...) where {N}
	Table(merge(map(columns, ts)...))
end

@inline function Broadcast.broadcasted(::Broadcast.DefaultArrayStyle{N}, f::GetProperty, t::Table{<:Any, N}) where {N}
	return f(t)
end

@inline function Broadcast.broadcasted(::Broadcast.DefaultArrayStyle{N}, f::GetProperties, t::Table{<:Any, N}) where {N}
    return f(t)
end

@inline function Broadcast.broadcasted(style::Broadcast.DefaultArrayStyle{N}, f::Compute{names}, t::Table{<:NamedTuple, N}) where {N, names}
    # minimize number of columns before iterating over the rows
    return Broadcast.broadcasted(style, f, GetProperties(names)(t))
end

@inline function Broadcast.broadcasted(style::Broadcast.DefaultArrayStyle{N}, f::Compute{names}, t::Table{<:NamedTuple{names}, N}) where {N, names}
    # efficient to iterate over rows with a minimal number of columns
    if length(names) == 1 # unwrap in the simple cases
        return Broadcast.broadcasted(f.f, getproperty(names[1])(t))
    elseif length(names) == 2
        return Broadcast.broadcasted(f.f, getproperty(names[1])(t), getproperty(names[2])(t))
    end

    invoke(Broadcast.broadcasted, Tuple{typeof(style), Function, typeof(t)}, style, f, t)
end

@inline function Broadcast.broadcasted(::Broadcast.DefaultArrayStyle{N}, f::Select, t::Table{<:Any, N}) where {N}
    return mapview(f, t)
end

# I'm not 100% sure how wise this pattern is...
Broadcast.materialize(t::Table) = Table(map(_materialize, columns(t)))
_materialize(x) = Broadcast.materialize(x)
_materialize(x::MappedArray) = copy(x)

# mapreduce

function Base.mapreduce(f::GetProperty, op, t::Table; kwargs...)
    return mapreduce(identity, op, f(t); kwargs...)
end

function Base.mapreduce(f::GetProperties, op, t::Table; kwargs...)
    return mapreduce(identity, op, f(t); kwargs...)
end

function Base.mapreduce(f::Compute{names}, op, t::Table; kwargs...) where {names}
    # minimize number of columns before iterating over the rows
    t2 = GetProperties(names)(t)
    return mapreduce(f, op, t2; kwargs...)
end

function Base.mapreduce(f::Compute{names}, op, t::Table{<:NamedTuple{names}}; kwargs...) where {names}
    # efficient to iterate over rows with a minimal number of columns
    if length(names) == 1 # unwrap in the simple cases
        return mapreduce(f.f, op, getproperty(names[1])(t))
    elseif length(names) == 2
        return mapreduce(f.f, op, getproperty(names[1])(t), getproperty(names[2])(t))
    end

    invoke(mapreduce, Tuple{Function, typeof(op), typeof(t)}, f, op, t; kwargs...)
end

# `filter(f, t)` defaults to `t[map(f, t)]`

function Base.filter(f::GetProperty, t::Table)
    return @inbounds t[f(t)::AbstractArray{Bool}]
end

# findall

function Base.findall(f::GetProperty, t::Table)
    return findall(identity, f(t))
end

function Base.findall(f::Compute{names}, t::Table) where {names}
    # minimize number of columns before iterating over the rows
    return findall(f, GetProperties(names)(t))
end

function Base.findall(f::Compute{names}, t::Table{<:NamedTuple{names}}) where {names}
    # efficient to iterate over rows with a minimal number of columns
    if length(names) == 1 # unwrap in the simple cases
        return findall(f.f, getproperty(names[1])(t))
    end

    invoke(findall, Tuple{Function, typeof(t)}, f, t)
end


text:
module Ecobee

export generate_pin, fetch_oauth_tokens, refresh_tokens,
    fetch_thermostat_ids, fetch_data,
    column_title, DEFAULT_COLUMNS, ALL_COLUMNS

using HTTP
using JSON
using Dates

include("auth.jl")
include("columns.jl")

function generate_pin(api_key, scope=WRITE_SCOPE)
    query = "response_type=ecobeePin&client_id=$api_key&scope=$scope"
    r = make_request("/authorize", query)
    return (pin=r["ecobeePin"], auth_code=r["code"], expires=r["expires_in"])
end

function fetch_oauth_tokens(api_key, auth_code)
    query = "grant_type=ecobeePin&code=$auth_code&client_id=$api_key"
    r = make_request("/token", query, method="POST")
    return (refresh_token=r["refresh_token"], access_token=r["access_token"], expires=r["expires_in"])
end

function refresh_tokens(refresh_token, api_key)
    query = "grant_type=refresh_token&code=$refresh_token&client_id=$api_key"
    r = make_request("/token", query, method="POST")
    return (refresh_token=r["refresh_token"], access_token=r["access_token"])
end

function fetch_thermostat_ids(access_token)
    body = JSON.json(Dict(
        "selection" => Dict(
            "selectionType" => "registered",
            "selectionMatch" => ""
        )
    ))
    r = make_request("/1/thermostat", "format=json&body=$body", access_token)
    return map(t -> t["identifier"], r["thermostatList"])
end

function fetch_data(access_token, thermostat_ids; columns=DEFAULT_COLUMNS, days_ago_start=1, days_ago_end=0)
    days_ago_date = d -> string(Dates.today() - Dates.Day(d))
    body = JSON.json(Dict(
        "startDate" => days_ago_date(days_ago_start),
        "endDate" => days_ago_date(days_ago_end),
        "columns" => join(columns, ","),
        "selection" => Dict(
            "selectionType" => "thermostats",
            "selectionMatch" => join(thermostat_ids, ","),
            "includeElectricity" => true
        )
    ))
    r = make_request("/1/runtimeReport", "format=json&body=$body", access_token)
    data = Dict()
    for report in r["reportList"]
        thermostat_id = report["thermostatIdentifier"]
        data[thermostat_id] = report["rowList"]
    end
    return data
end

function column_title(column)
    return COLUMN_TITLE_LOOKUP[column]
end

function make_request(path, query, access_token=nothing; method="GET")
    url = "https://api.ecobee.com$path"
    headers = access_token != nothing ? Dict("Authorization" => "Bearer $access_token") : []
    r = HTTP.request(method, url, headers; query=query)
    return JSON.parse(String(r.body))
end

end


text:
__precompile__()
module HierarchicalMatrices
    using LinearAlgebra

    BLOCKRANK(T::Type{R}) where {R<:Real} = 2round(Int, half(T)*log(3+sqrt(T(8)), inv(eps(T))))
    BLOCKRANK(T::Type{C}) where {C<:Complex} = BLOCKRANK(real(T))
    BLOCKSIZE(T) = 4BLOCKRANK(T)

    import Base: convert, view, size
    import Base: copy, getindex, setindex!, show, one, zero, inv, isless
    import Base: div, rem
    import Base: broadcast, Matrix, promote_op
    import Base: +, -, *, /, \, ==, !=
    import LinearAlgebra: Factorization, rank, norm, cond, istriu, istril, issymmetric, ishermitian,
                                   transpose, adjoint
    import LinearAlgebra.BLAS: @blasfunc, libblas, BlasInt, BlasFloat, BlasReal, BlasComplex

    export BLOCKSIZE, BLOCKRANK, Block
    export AbstractLowRankMatrix, AbstractBarycentricMatrix
    export LowRankMatrix, BarycentricMatrix2D, EvenBarycentricMatrix
    export @hierarchical, barycentricmatrix, blocksize, blockgetindex
    export add_col, add_col!, lrzeros, indsplit

    """
    Compute a typed 0.5.
    """
    half(x::Number) = oftype(x, 0.5)
    half(x::Integer) = half(float(x))
    half(::Type{T}) where {T<:Number} = convert(T, 0.5)
    half(::Type{T}) where {T<:Integer} = half(AbstractFloat)

    """
    Compute a typed 2.
    """
    two(x::Number) = oftype(x,2)
    two(::Type{T}) where {T<:Number} = convert(T, 2)

    mul!(args...) = LinearAlgebra.mul!(args...)

    include("LowRankMatrix.jl")
    include("BarycentricMatrix.jl")
    include("block.jl")
    include("hierarchical.jl")
    include("HierarchicalMatrix.jl")
    include("KernelMatrix.jl")
    include("algebra.jl")

end # module


text:
using .DataFrames

function meltdata(U::AbstractDataFrame, colgroups_::Vector{Col.GroupedColumn})
    um, un = size(U)

    colgroups = Set(colgroups_)

    # Figure out the size of the new melted matrix
    allcolumns = Set{Symbol}(names(U))

    vm = um
    grouped_columns = Set{Symbol}()
    for colgroup in colgroups
        if colgroup.columns===nothing # null => group all columns
            vm *= un
            grouped_columns = copy(allcolumns)
        else
            for j in colgroup.columns
                if !isa(j, Symbol)
                    error("DataFrame columns can only be grouped by (Symbol) names")
                end
                push!(grouped_columns, j)
            end
            vm *= length(colgroup.columns)
        end
    end

    ungrouped_columns = setdiff(allcolumns, grouped_columns)
    vn = length(colgroups) + length(ungrouped_columns)

    V = AbstractArray[]
    vnames = Symbol[]
    colmap = Dict{Any, Int}()

    # allocate vectors for grouped columns
    for (j, colgroup) in enumerate(colgroups)
        cols = colgroup.columns===nothing ? allcolumns : colgroup.columns

        # figure the grouped common column type
        firstcol = U[first(cols)]
        eltyp = eltype(firstcol)
        vectyp = isa(firstcol, Vector) ? Vector : DataVector
        for col in cols
            eltyp = promote_type(eltyp, typeof(U[col]))
            if !isa(U[col], Vector)
                vectyp = DataVector
            end
        end

        push!(V, eltyp == Vector ? Array{eltyp}(undef, vm) : Array{Union{Nothing,eltyp}}(undef, vm))
        name = gensym()
        push!(vnames, name)
        colmap[colgroup] = j
    end

    # allocate vectors for ungrouped columns
    for (j, col) in enumerate(ungrouped_columns)
        push!(V, similar(U[col], vm))
        colmap[col] = j + length(colgroups)
        push!(vnames, col)
    end

    # Indicator columns for each colgroup
    col_indicators = Array{Symbol}(undef, vm, length(colgroups))
    row_indicators = Array{Int}(undef, vm, length(colgroups))

    colidxs = [colgroup.columns===nothing ? collect(allcolumns) : colgroup.columns
               for colgroup in colgroups]

    vi = 1
    for ui in 1:um
        for colidx in Iterators.product(colidxs...)
            # copy grouped columns
            for (vj, uj) in enumerate(colidx)
                V[vj][vi] = U[ui, uj]
                col_indicators[vi, vj] = uj
                row_indicators[vi, vj] = ui
            end

            # copy ungrouped columns
            for (vj, uj) in enumerate(ungrouped_columns)
                V[vj + length(colgroups)][vi] = U[ui, uj]
            end

            vi += 1
        end
    end

    df = DataFrame(; collect(zip(vnames, V))...)
    return MeltedData(U, df, row_indicators, col_indicators, colmap)
end

evalmapping(source::MeltedData{T}, arg::Col.GroupedColumnValue) where T<:AbstractDataFrame =
    source.melted_data[source.colmap[Col.GroupedColumn(arg.columns)]]

evalmapping(source::AbstractDataFrame, arg::Symbol) = source[arg]
evalmapping(source::AbstractDataFrame, arg::AbstractString) = evalmapping(source, Symbol(arg))
evalmapping(source::AbstractDataFrame, arg::Integer) = source[arg]
evalmapping(source::AbstractDataFrame, arg::Expr) = with(source, arg)


text:
struct V3{T}
    x::T
    y::T
    z::T
end
@inline Base.zero(::Type{V3{T}}) where T = V3(zero(T), zero(T), zero(T))
@inline Base.zero(::V3{T}) where T = zero(V3{T})
@inline Base.isapprox(a::V3, b::V3; kwargs...) = isapprox(a.x, b.x; kwargs...) && isapprox(a.y, b.y; kwargs...) && isapprox(a.z, b.z; kwargs...)
function Base.getindex(v::V3, i::Integer)
    @boundscheck i==1 || i==2 || i==3
    if i==1
        v.x
    elseif i==2
        v.y
    else
        v.z
    end
end

# +=*/
for OP in [:+, :-]
    @inline @eval function Base.$OP(v::V3, w::V3)
        V3($OP(v.x, w.x), $OP(v.y, w.y), $OP(v.z, w.z))
    end
    @eval @i @inline @eval function :(+=)($OP)(y!::V3, v::V3, w::V3)
        y!.x += $OP(v.x, w.x)
        y!.y += $OP(v.y, w.y)
        y!.z += $OP(v.z, w.z)
    end
end

@i @inline @eval function :(+=)(identity)(y!::V3, w::V3)
    y!.x += w.x
    y!.y += w.y
    y!.z += w.z
end

@inline function Base.:(*)(v::V3, x::Real)
    V3(v.x*x, v.y*x, v.z*x)
end
@inline function Base.:(*)(x::Real, v::V3)
    v*x
end

@i @inline function :(+=)(*)(y!::V3, v::V3, x::Real)
    y!.x += v.x * x
    y!.y += v.y * x
    y!.z += v.z * x
end
@i @inline function :(+=)(*)(y!::V3, x::Real, v::V3)
    y! += v*x
end

# norms and distances
@inline function LinearAlgebra.norm2(v::V3)
    sqrt(v.x^2 + v.y^2 + v.z^2)
end
@inline function sqnorm2(v::V3)
    v.x^2 + v.y^2 + v.z^2
end
@inline function sqdistance(r1::V3, r2::V3)
    (r1.x - r2.x)^2 + (r1.y-r2.y)^2 + (r1.z-r2.z)^2
end
@inline function distance(r1::V3, r2::V3)
    sqrt((r1.x - r2.x)^2 + (r1.y-r2.y)^2 + (r1.z-r2.z)^2)
end

@i @inline function :(+=)(sqnorm2)(y!::Real, v::V3)
    y! += v.x^2
    y! += v.y^2
    y! += v.z^2
end

using NiLang.AD
for TP in [:Real, :GVar]
    @eval @i @inline function :(+=)(sqdistance)(y!::$TP, r1::V3{T}, r2::V3{T}) where T
        @routine @invcheckoff begin
            @zeros T dx dy dz
            dx += r1.x - r2.x
            dy += r1.y - r2.y
            dz += r1.z - r2.z
        end
        y! += dx^2
        y! += dy^2
        y! += dz^2
        ~@routine
    end
end