text:
# Random Feature/Fourier Expansion

"""
    RFE(D::Int, # number of basis functions
        n::Int, # number of input variables

        Ω::Matrix{T} where T <: AbstractFloat, # random feature coefficients
        b::Vector{T} where T <: AbstractFloat, # random feature phase

        c::Vector{T} where T <: AbstractFloat, # basis function coefficients
        P::Matrix{T} where T <: AbstractFloat, # P is the covariance matrix of the coefficients
        λ::AbstractFloat, # regularization parameter of the coefficients

        variable_offset::Bool,
        h::Vector{T} where T <: AbstractFloat, # history variable for computing offsets
        offset::AbstractFloat) # the currently used offset

A Random Feature/Fourier Expansion (RFE) is the approximation of a function y = f(x) using basis functions that are a random nonlinear mapping of the input: f(x) ≈ ∑ci gi(x).
The vector c contains the coefficients. The value of the basis functions are computed through g(x) = cos.(Ω*x + b), where Ω is of size D x n and Ω and b consist of randomly generated coefficients.
c (with inverse covariance matrix P) can be estimated from data through solving a simple linear least squares problem, typically estimated with regularization term λ||c||^2.

An offset is being calculated if variable_offset=true, so that f(x) + v = ∑(ci + hi) gi(x).
"""
mutable struct RFE
    D::Int # number of basis functions
    n::Int # number of input variables

    Ω::Matrix{T} where T <: AbstractFloat # random feature coefficients
    b::Vector{T} where T <: AbstractFloat # random feature phase

    c::Vector{T} where T <: AbstractFloat # basis function coefficients
    P::Matrix{T} where T <: AbstractFloat # P is the covariance matrix of the coefficients
    λ::AbstractFloat # regularization parameter of the coefficients

    variable_offset::Bool
    h::Vector{T} where T <: AbstractFloat # history variable for computing offsets
    offset::AbstractFloat

    function RFE(D,n,Ω,b,c,P,λ,variable_offset,h,offset)
        if size(Ω) != (D,n)
            DimensionMismatch("Dimensions of Ω ($(size(Ω,1)),$(size(Ω,2))) do not match D ($D),n ($n).") |> throw
        end
        if length(b) != D
            DimensionMismatch("Length of b does not match D.") |> throw
        end
        if length(c) != D
            DimensionMismatch("Length of c does not match D.") |> throw
        end
        if length(h) != D
            DimensionMismatch("Length of h does not match D.") |> throw
        end
        if size(P) != (D,D)
            DimensionMismatch("Size of P does not match (D,D).") |> throw
        end
        @assert λ >= 0
        @assert issymmetric(P)
        @assert all(eigvals(P) .> 0)

        new(D,n,Ω,b,c,P,λ,variable_offset,h,offset)
    end
end

"""
    RFE(D::Int, n::Int, distribution_Ω::Distributions.Distribution, distribution_b::Distributions.Distribution; λ::AbstractFloat=1E-6, variable_offset::Bool)

Pick the random coefficients of the matrix Ω and vector b from the distributions distribution_Ω and distribution_b.
"""
function RFE(D::Int, n::Int, distribution_Ω::Distributions.Distribution, distribution_b::Distributions.Distribution; λ::AbstractFloat=1E-6, variable_offset::Bool=true)
    Ω = rand(distribution_Ω, D, n)
    b = rand(distribution_b, D)
    c = zeros(Float64,D)
    P = Diagonal(ones(Float64,D)./λ) |> Matrix

    RFE(D,n,Ω,b,c,P,λ,variable_offset,zeros(Float64,D),0.)
end

"""
    RFE(D::Int,n::Int,σ_coef::AbstractFloat;λ::AbstractFloat=1E-6,variable_offset::Bool=true)

Pick the random coefficients of the matrix Ω from N(0,σI) and b from U(0,2π)
"""
function RFE(D::Int,n::Int,σ_coef::AbstractFloat;λ::AbstractFloat=1E-6,variable_offset::Bool=true)
    distribution_Ω = Distributions.Normal(0.,σ_coef)
    distribution_b = Distributions.Uniform(0.,2π)

    RFE(D,n,distribution_Ω, distribution_b, λ=λ, variable_offset=variable_offset)
end

"""
    featureexpand(rfe::RFE,x::Vector{T} where T <: AbstractFloat)

Compute g(x)::Vector
"""
function featureexpand(rfe::RFE,x::Vector{T} where T <: AbstractFloat)
    return cos.(rfe.Ω*x + rfe.b)
end

"""
    evaluateRFE(rfe::RFE,x::AbstractVector)

Compute ∑ci gi(x)
"""
function evaluateRFE(rfe::RFE,x::AbstractVector)
    return dot(rfe.c, featureexpand(rfe,x)) - rfe.offset
end

"""
    updateRFE!(rfe::RFE,x::Vector{T} where T <: AbstractFloat,y::AbstractFloat)

The RFE is updated in a recursive manner so that adding new measurements is easy.
The vector x is the input, and the vector y is the measured output.

See: page 37 of the dissertation of Laurens Bliek - Automatic Tuning of Photonic Beamformers
"""
function updateRFE!(rfe::RFE,x::Vector{T} where T <: AbstractFloat,y::AbstractFloat)
    an = featureexpand(rfe, x)
    γn = 1.0 / (1.0 + dot(an,rfe.P*an))
    gn = γn * rfe.P * an
    rfe.c[:] = rfe.c + gn*(y - dot(an,rfe.c)) # update coeffients
    rfe.P[:] = rfe.P - (gn * gn') ./ γn
    rfe.P[:] = (rfe.P + rfe.P')/2
    return an, gn
end

"""
    downdateRFE!(rfe::RFE,x::Vector{T} where T <: AbstractFloat,y::AbstractFloat)

The RFE is downdated in a recursive manner so we can remove the effect of past measurements if we use a sliding window.
The vector x is the input to be removed, and the vector y is the measured output to be removed.

See: page 111 of the dissertation of Laurens Bliek - Automatic Tuning of Photonic Beamformers
"""
function downdateRFE!(rfe::RFE,x::AbstractVector{Float64},y::AbstractFloat)
    an = featureexpand(rfe, x)
    γn = 1.0 / (1.0 - dot(an,rfe.P*an))
    gn = γn * rfe.P * an
    rfe.c[:] = rfe.c - gn*(y - dot(an,rfe.c)) # update coeffients
    rfe.P[:] = rfe.P + (gn * gn') ./ γn
    rfe.P[:] = (rfe.P + rfe.P')/2
    return an, gn
end


text:
# Julia wrapper over typeGRAMPCparam, which contains all problem-specific parameters
mutable struct CGRAMPCparam
    # Problem size variables
    Nx::typeInt
    Nu::typeInt
    Np::typeInt
    Ng::typeInt
    Nh::typeInt
    NgT::typeInt
    NhT::typeInt
    Nc::typeInt

    x0::Ptr{typeRNum}
    xdes::Ptr{typeRNum}

    u0::Ptr{typeRNum}
    udes::Ptr{typeRNum}
    umax::Ptr{typeRNum}
    umin::Ptr{typeRNum}

    p0::Ptr{typeRNum}
    pmax::Ptr{typeRNum}
    pmin::Ptr{typeRNum}

    Thor::typeRNum
    Tmax::typeRNum
    Tmin::typeRNum

    dt::typeRNum
    t0::typeRNum
end

# Julia wrapper typeGRAMPCopt, which contains all algorithm-related options
mutable struct CGRAMPCopt
    Nhor::typeInt
    MaxGradIter::typeInt
    MaxMultIter::typeInt
    ShiftControl::typeInt

    TimeDiscretization::typeInt

    IntegralCost::typeInt
    TerminalCost::typeInt
    IntegratorCost::typeInt

    # ODE Integrator options
    Integrator::typeInt
    IntegratorRelTol::typeRNum
    IntegratorAbsTol::typeRNum
    IntegratorMinStepSize::typeRNum
    IntegratorMaxSteps::typeInt
    FlagsRodas::Ptr{typeInt} # 8-element vector

    # Line search options
    LineSearchType::typeInt
    LineSearchExpAutoFallback::typeInt
    LineSearchMax::typeRNum
    LineSearchMin::typeRNum
    LineSearchInit::typeRNum
    LineSearchAdaptAbsTol::typeRNum
    LineSearchAdaptFactor::typeRNum
    LineSearchIntervalTol::typeRNum
    LineSearchIntervalFactor::typeRNum

    OptimControl::typeInt
    OptimParam::typeInt
    OptimParamLineSearchFactor::typeRNum
    OptimTime::typeInt
    OptimTimeLineSearchFactor::typeRNum

    ScaleProblem::typeInt
    xScale::Ptr{typeRNum}
    xOffset::Ptr{typeRNum}
    uScale::Ptr{typeRNum}
    uOffset::Ptr{typeRNum}
    pScale::Ptr{typeRNum}
    pOffset::Ptr{typeRNum}

    TScale::typeRNum
    TOffset::typeRNum
    JScale::typeRNum

    cScale::Ptr{typeRNum};

    EqualityConstraints::typeInt
    InequalityConstraints::typeInt
    TerminalEqualityConstraints::typeInt
    TerminalInequalityConstraints::typeInt
    ConstraintsHandling::typeInt

    ConstraintsAbsTol::Ptr{typeRNum};

    MultiplierMax::typeRNum
    MultiplierDampingFactor::typeRNum
    PenaltyMax::typeRNum
    PenaltyMin::typeRNum
    PenaltyIncreaseFactor::typeRNum
    PenaltyDecreaseFactor::typeRNum
    PenaltyIncreaseThreshold::typeRNum
    AugLagUpdateGradientRelTol::typeRNum

    ConvergenceCheck::typeInt
    ConvergenceGradientRelTol::typeRNum
end

# Julia wrapper for typeGRAMPCsol, which contains the (public) solution data
mutable struct CGRAMPCsol
    xnext::Ptr{typeRNum}
    unext::Ptr{typeRNum}
    pnext::Ptr{typeRNum}

    Tnext::typeRNum
    J::Ptr{typeRNum} # A 2-element array
    cfct::typeRNum
    pen::typeRNum

    iter::Ptr{typeInt}

    status::typeInt
end

# Julia wrapper for typeGRAMPCrws, which contains the (private) workspace data
mutable struct CGRAMPCrws
    t::Ptr{typeRNum}
    tls::Ptr{typeRNum}

    x::Ptr{typeRNum}
    adj::Ptr{typeRNum}
    dcdx::Ptr{typeRNum}

    u::Ptr{typeRNum}
    uls::Ptr{typeRNum}
    uprev::Ptr{typeRNum}
    gradu::Ptr{typeRNum}
    graduprev::Ptr{typeRNum}
    dcdu::Ptr{typeRNum}

    p::Ptr{typeRNum}
    pls::Ptr{typeRNum}
    pprev::Ptr{typeRNum}
    gradp::Ptr{typeRNum}
    gradpprev::Ptr{typeRNum}
    dcdp::Ptr{typeRNum}

    T::typeRNum
    Tprev::typeRNum
    gradT::typeRNum
    gradTprev::typeRNum
    dcdt::typeRNum

    mult::Ptr{typeRNum}
    pen::Ptr{typeRNum}
    cfct::Ptr{typeRNum}
    cfctprev::Ptr{typeRNum}
    cfctAbsTol::Ptr{typeRNum}

    lsAdapt::Ptr{typeRNum}
    lsExplicit::Ptr{typeRNum}
    rwsScale::Ptr{typeRNum}

    lrwsGeneral::typeInt
    rwsGeneral::Ptr{typeRNum}

    lworkRodas::typeInt
    liworkRodas::typeInt

    rparRodas::Ptr{typeRNum}
    iparRodas::Ptr{typeInt}
    workRodas::Ptr{typeRNum}
    iworkRodas::Ptr{typeInt}
end

# Julia wrapper for typeGRAMPC, which is the main solver type
mutable struct CGRAMPC
    param::Ptr{GRAMPC.CGRAMPCparam}
    opt::Ptr{GRAMPC.CGRAMPCopt}
    sol::Ptr{GRAMPC.CGRAMPCsol}
    rws::Ptr{GRAMPC.CGRAMPCrws}
    userparam::Ptr{typeUSERPARAM};
end

ext: jl


text:
@parameters t
const TIME = t
D = Differential(t)

@variables E[1:2](t) F[1:2](t)

# Linear resistance (:R)
@parameters R
R_dict = Dict(
    :description => "Generalised Linear Resistor",
    :numports => 1,
    :parameters => OrderedDict(
        R => "Resistance"
    ),
    :equations => [0 ~ E[1]/R - F[1]]
)

# Linear capacitor (:C)
@parameters C
@variables q(t)
C_dict = Dict(
    :description => "Generalised Linear Resistor",
    :numports => 1,
    :parameters => OrderedDict(
        C => "Capacitance"
    ),
    :state_vars => OrderedDict(
        q => "Generalised Position"
    ),
    :equations => [
        0 ~ q/C - E[1],
        D(q) ~ F[1]
    ]
)

# Linear inductance (:I)
@parameters L
@variables p(t)
I_dict = Dict(
    :description => "Generalised Linear Inductor",
    :numports => 1,
    :parameters => OrderedDict(
        L => "Inductance"
    ),
    :state_vars => OrderedDict(
        p => "Generalised Momentum"
    ),
    :equations => [
        0 ~ p/L - F[1],
        D(p) ~ E[1]
    ]
)

# Source of effort (:Se)
@parameters e
Se_dict = Dict(
    :description => "Effort source",
    :numports => 1,
    :parameters => OrderedDict(
        e => "Effort"
    ),
    :equations => [0 ~ e - E[1]]
)

# Source of flow (:Sf)
@parameters f
Sf_dict = Dict(
    :description => "Flow source",
    :numports => 1,
    :parameters => OrderedDict(
        e => "Flow"
    ),
    :equations => [0 ~ f - F[1]]
)

# Transformer (:TF)
@parameters r
TF_dict = Dict(
    :description => "Linear Transformer",
    :numports => 2,
    :parameters => OrderedDict(
        r => "Winding ratio"
    ),
    :equations => [
        0 ~ E[2] - r*E[1],
        0 ~ F[1] - r*F[2]
    ]
)

# Chemical species (:Ce)
@parameters k R T
@variables q_1(t)
Ce_dict = Dict(
    :description => "Chemical species",
    :numports => 1,
    :parameters => OrderedDict(
        k => "Biochemical Constant; exp(mu_0/RT)/V",
        R => "Universal Gas Constant",
        T => "Temperature"
    ),
    :state_vars => OrderedDict(
        q => "Molar Quantity"
    ),
    :equations => [
        0 ~ R*T*log(k*q) - E[1],
        D(q) ~ F[1]
    ]
)

# Normalised chemical species (:ce)
@parameters k
@variables q(t)
ce_dict = Dict(
    :description => "Chemical species (normalised)",
    :numports => 1,
    :parameters => OrderedDict(
        k => "Biochemical Constant; exp(mu_0/RT)/V"
    ),
    :state_vars => OrderedDict(
        q => "Molar Quantity"
    ),
    :equations => [
        0 ~ log(k*q) - E[1],
        D(q) ~ F[1]
    ]
)

# Chemical reaction (:Re)
@parameters r R T
Re_dict = Dict(
    :description => "Biochemical reaction",
    :numports => 2,
    :parameters => OrderedDict(
        r => "Reaction Rate",
        R => "Universal Gas Constant",
        T => "Temperature"
    ),
    :equations => [
        0 ~ F[1] + F[2],
        0 ~ F[1] - r*(exp(E[1]/R/T) - exp(E[2]/R/T))
    ]
)

# Normalised chemical reaction (:re)
@parameters r R T
re_dict = Dict(
    :description => "Biochemical reaction (normalised)",
    :numports => 2,
    :parameters => OrderedDict(
        r => "Reaction Rate"
    ),
    :equations => [
        0 ~ F[1] + F[2],
        0 ~ F[1] - r*(exp(E[1]) - exp(E[2]))
    ]
)

standard_library = Dict(
    :R => R_dict,
    :C => C_dict,
    :I => I_dict,
    :Se => Se_dict,
    :Sf => Sf_dict,
    :TF => TF_dict,
    :Ce => Ce_dict,
    :ce => ce_dict,
    :Re => Re_dict,
    :re => re_dict
)

max_stars_repo_name: jedforrest/BondGraphs

max_issues_repo_name: jedforrest/BondGraphs

max_forks_repo_name: jedforrest/BondGraphs


text:
function pvcam_get_ver(pvcam_version::Ref{UInt16})::Bool
    ccall((:pl_pvcam_get_ver, :libpvcam), UInt16, (Ref{UInt16},), pvcam_version)
end

function pvcam_init()::Bool
    ccall((:pl_pvcam_init, :libpvcam), UInt16, ())
end

function pvcam_uninit()::Bool
    ccall((:pl_pvcam_uninit, :libpvcam), UInt16, ())
end

function cam_check(hcam::Integer)::Bool
    ccall((:pl_cam_check, :libpvcam), UInt16, (Int16,), hcam)
end

function cam_close(hcam::Integer)::Bool
    ccall((:pl_cam_close, :libpvcam), UInt16, (Int16,), hcam)
end

function cam_get_name(cam_num::Integer, camera_name::Vector{UInt8})::Bool
    ccall((:pl_cam_get_name, :libpvcam), UInt16, (Int16, Ref{UInt8}), cam_num, camera_name)
end

function cam_get_total(totl_cams::Ref{Int16})::Bool
    ccall((:pl_cam_get_total, :libpvcam), UInt16, (Ref{Int16},), totl_cams)
end

function cam_open(camera_name, hcam::Ref{Int16}, o_mode::Integer)::Bool
    ccall((:pl_cam_open, :libpvcam), UInt16, (Cstring, Ref{Int16}, Int16), camera_name, hcam, o_mode)
end

function cam_register_callback(hcam::Integer, callback_event::Integer, callback::Ref{Cvoid})::Bool
    ccall((:pl_cam_register_callback, :libpvcam), UInt16, (Int16, Cint, Ref{Cvoid}), hcam, callback_event, callback)
end

function cam_register_callback_ex(hcam::Integer, callback_event::Integer, callback::Ref{Cvoid}, context::Ref{Cvoid})::Bool
    ccall((:pl_cam_register_callback_ex, :libpvcam), UInt16, (Int16, Cint, Ref{Cvoid}, Ref{Cvoid}), hcam, callback_event, callback, context)
end

function cam_register_callback_ex2(hcam::Integer, callback_event::Integer, callback::Ref{Cvoid})::Bool
    ccall((:pl_cam_register_callback_ex2, :libpvcam), UInt16, (Int16, Cint, Ref{Cvoid}), hcam, callback_event, callback)
end

function cam_register_callback_ex3(hcam::Integer, callback_event::Integer, callback::Ref{Cvoid}, context::Ref{Cvoid})::Bool
    ccall((:pl_cam_register_callback_ex3, :libpvcam), UInt16, (Int16, Cint, Ref{Cvoid}, Ref{Cvoid}), hcam, callback_event, callback, context)
end

function cam_deregister_callback(hcam::Integer, callback_event::Integer)::Bool
    ccall((:pl_cam_deregister_callback, :libpvcam), UInt16, (Int16, Cint), hcam, callback_event)
end

function error_code()
    ccall((:pl_error_code, :libpvcam), Int16, ())
end

function error_message(err_code::Integer, msg)::Bool
    ccall((:pl_error_message, :libpvcam), UInt16, (Int16, Ref{UInt8}), err_code, msg)
end

function get_param(hcam::Integer, param_id::Integer, param_attribute::Integer, param_value::Ref{Cvoid})::Bool
    ccall((:pl_get_param, :libpvcam), UInt16, (Int16, UInt32, Int16, Ref{Cvoid}), hcam, param_id, param_attribute, param_value)
end

function set_param(hcam::Integer, param_id::Integer, param_value::Ref{Cvoid})::Bool
    ccall((:pl_set_param, :libpvcam), UInt16, (Int16, UInt32, Ref{Cvoid}), hcam, param_id, param_value)
end

function get_enum_param(hcam::Integer, param_id::Integer, index::Integer, value::Ref{Cint}, desc, length::Integer)::Bool
    ccall((:pl_get_enum_param, :libpvcam), UInt16, (Int16, UInt32, UInt32, Ref{Cint}, Cstring, UInt32), hcam, param_id, index, value, desc, length)
end

function enum_str_length(hcam::Integer, param_id::Integer, index::Integer, length::Ref{UInt32})::Bool
    ccall((:pl_enum_str_length, :libpvcam), UInt16, (Int16, UInt32, UInt32, Ref{UInt32}), hcam, param_id, index, length)
end

function pp_reset(hcam::Integer)::Bool
    ccall((:pl_pp_reset, :libpvcam), UInt16, (Int16,), hcam)
end

function create_smart_stream_struct(array::Ref{Ptr{smart_stream_type}}, entries::Integer)::Bool
    ccall((:pl_create_smart_stream_struct, :libpvcam), UInt16, (Ref{Ptr{smart_stream_type}}, UInt16), array, entries)
end

function release_smart_stream_struct(array::Ref{Ptr{smart_stream_type}})::Bool
    ccall((:pl_release_smart_stream_struct, :libpvcam), UInt16, (Ref{Ptr{smart_stream_type}},), array)
end

function create_frame_info_struct(new_frame::Ref{Ptr{FRAME_INFO}})::Bool
    ccall((:pl_create_frame_info_struct, :libpvcam), UInt16, (Ref{Ptr{FRAME_INFO}},), new_frame)
end

function release_frame_info_struct(frame_to_delete::Ref{FRAME_INFO})::Bool
    ccall((:pl_release_frame_info_struct, :libpvcam), UInt16, (Ref{FRAME_INFO},), frame_to_delete)
end

function exp_setup_seq(hcam::Integer, exp_total::Integer, rgn_total::Integer, rgn_array::Ref{rgn_type}, exp_mode::Integer, exposure_time::Integer, exp_bytes::Ref{UInt32})::Bool
    ccall((:pl_exp_setup_seq, :libpvcam), UInt16, (Int16, UInt16, UInt16, Ref{rgn_type}, Int16, UInt32, Ref{UInt32}), hcam, exp_total, rgn_total, rgn_array, exp_mode, exposure_time, exp_bytes)
end

function exp_start_seq(hcam::Integer, pixel_stream::Vector{UInt16})::Bool
    ccall((:pl_exp_start_seq, :libpvcam), UInt16, (Int16, Ref{UInt16}), hcam, pixel_stream)
end

function exp_setup_cont(hcam::Integer, rgn_total::Integer, rgn_array::Vector{rgn_type}, exp_mode::Integer, exposure_time::Integer, exp_bytes::Ref{UInt32}, buffer_mode::PL_CIRC_MODES)::Bool
    ccall((:pl_exp_setup_cont, :libpvcam), UInt16, (Int16, UInt16, Ref{rgn_type}, Int16, UInt32, Ref{UInt32}, Int16), hcam, rgn_total, rgn_array, exp_mode, exposure_time, exp_bytes, buffer_mode)
end

function exp_start_cont(hcam::Integer, pixel_stream::Vector{UInt16}, size::Integer)::Bool
    ccall((:pl_exp_start_cont, :libpvcam), UInt16, (Int16, Ref{UInt16}, UInt32), hcam, pixel_stream, size)
end

function exp_check_status(hcam::Integer, status::Ref{Int16}, bytes_arrived::Ref{UInt32})::Bool
    ccall((:pl_exp_check_status, :libpvcam), UInt16, (Int16, Ref{Int16}, Ref{UInt32}), hcam, status, bytes_arrived)
end

function exp_check_cont_status(hcam::Integer, status::Ref{Int16}, bytes_arrived::Ref{UInt32}, buffer_cnt::Ref{UInt32})::Bool
    ccall((:pl_exp_check_cont_status, :libpvcam), UInt16, (Int16, Ref{Int16}, Ref{UInt32}, Ref{UInt32}), hcam, status, bytes_arrived, buffer_cnt)
end

function exp_check_cont_status_ex(hcam::Integer, status::Ref{Int16}, byte_cnt::Ref{UInt32}, buffer_cnt::Ref{UInt32}, frame_info::Ref{FRAME_INFO})::Bool
    ccall((:pl_exp_check_cont_status_ex, :libpvcam), UInt16, (Int16, Ref{Int16}, Ref{UInt32}, Ref{UInt32}, Ref{FRAME_INFO}), hcam, status, byte_cnt, buffer_cnt, frame_info)
end

function exp_get_latest_frame(hcam::Integer, frame::Ref{Ptr{UInt16}})::Bool
    ccall((:pl_exp_get_latest_frame, :libpvcam), UInt16, (Int16, Ref{Ptr{UInt16}}), hcam, frame)
end

function exp_get_latest_frame_ex(hcam::Integer, frame::Ref{Ptr{UInt16}}, frame_info::Ref{FRAME_INFO})::Bool
    ccall((:pl_exp_get_latest_frame_ex, :libpvcam), UInt16, (Int16, Ref{Ptr{UInt16}}, Ref{FRAME_INFO}), hcam, frame, frame_info)
end

function exp_get_oldest_frame(hcam::Integer, frame::Ref{Ptr{Cvoid}})::Bool
    ccall((:pl_exp_get_oldest_frame, :libpvcam), UInt16, (Int16, Ref{Ptr{Cvoid}}), hcam, frame)
end

function exp_get_oldest_frame_ex(hcam::Integer, frame::Ref{Ptr{Cvoid}}, frame_info::Ref{FRAME_INFO})::Bool
    ccall((:pl_exp_get_oldest_frame_ex, :libpvcam), UInt16, (Int16, Ref{Ptr{Cvoid}}, Ref{FRAME_INFO}), hcam, frame, frame_info)
end

function exp_unlock_oldest_frame(hcam::Integer)::Bool
    ccall((:pl_exp_unlock_oldest_frame, :libpvcam), UInt16, (Int16,), hcam)
end

function exp_stop_cont(hcam::Integer, cam_state::PL_CCS_ABORT_MODES)::Bool
    ccall((:pl_exp_stop_cont, :libpvcam), UInt16, (Int16, Int16), hcam, cam_state)
end

function exp_abort(hcam::Integer, cam_state::Integer)::Bool
    ccall((:pl_exp_abort, :libpvcam), UInt16, (Int16, Int16), hcam, cam_state)
end

function exp_finish_seq(hcam::Integer, pixel_stream::Ref{Cvoid}, hbuf::Integer)::Bool
    ccall((:pl_exp_finish_seq, :libpvcam), UInt16, (Int16, Ref{Cvoid}, Int16), hcam, pixel_stream, hbuf)
end

function io_script_control(hcam::Integer, addr::Integer, state::Cdouble, location::Integer)::Bool
    ccall((:pl_io_script_control, :libpvcam), UInt16, (Int16, UInt16, Cdouble, UInt32), hcam, addr, state, location)
end

function io_clear_script_control(hcam::Integer)::Bool
    ccall((:pl_io_clear_script_control, :libpvcam), UInt16, (Int16,), hcam)
end

function md_frame_decode(pDstFrame::Ref{md_frame}, pSrcBuf::Ref{Cvoid}, srcBufSize::Integer)::Bool
    ccall((:pl_md_frame_decode, :libpvcam), UInt16, (Ref{md_frame}, Ref{Cvoid}, UInt32), pDstFrame, pSrcBuf, srcBufSize)
end

function md_frame_recompose(pDstBuf::Ref{Cvoid}, offX::Integer, offY::Integer, dstWidth::Integer, dstHeight::Integer, pSrcFrame::Ref{md_frame})::Bool
    ccall((:pl_md_frame_recompose, :libpvcam), UInt16, (Ref{Cvoid}, UInt16, UInt16, UInt16, UInt16, Ref{md_frame}), pDstBuf, offX, offY, dstWidth, dstHeight, pSrcFrame)
end

function md_create_frame_struct_cont(pFrame::Ref{Ptr{md_frame}}, roiCount::Integer)::Bool
    ccall((:pl_md_create_frame_struct_cont, :libpvcam), UInt16, (Ref{Ptr{md_frame}}, UInt16), pFrame, roiCount)
end

function md_create_frame_struct(pFrame::Ref{Ptr{md_frame}}, pSrcBuf::Ref{Cvoid}, srcBufSize::Integer)::Bool
    ccall((:pl_md_create_frame_struct, :libpvcam), UInt16, (Ref{Ptr{md_frame}}, Ref{Cvoid}, UInt32), pFrame, pSrcBuf, srcBufSize)
end

function md_release_frame_struct(pFrame::Ref{md_frame})::Bool
    ccall((:pl_md_release_frame_struct, :libpvcam), UInt16, (Ref{md_frame},), pFrame)
end

function md_read_extended(pOutput::Ref{md_ext_item_collection}, pExtMdPtr::Ref{Cvoid}, extMdSize::Integer)::Bool
    ccall((:pl_md_read_extended, :libpvcam), UInt16, (Ref{md_ext_item_collection}, Ref{Cvoid}, UInt32), pOutput, pExtMdPtr, extMdSize)
end


text:
function NonNegativeConstraint(arg0::jboolean)
    return NonNegativeConstraint((jboolean,), arg0)
end

function is_restricted_to_non_negative(obj::NonNegativeConstraint)
    return jcall(obj, "isRestrictedToNonNegative", jboolean, ())
end


text:
using Pkg
Pkg.activate(".")
using DataFrames

global ULLR_required_columns = [ :upper_left_x, :lower_right_x, :upper_left_y, :lower_right_y ]

"""
    has_required_columns( df::DataFrame, required_columns::Vector )::Bool

Returns `true` if a `DataFrame` instance has all of the `required_columns` present.

"""
function has_required_columns( df::DataFrame, required_columns::Vector )::Bool
    df_names = names( df )
    result = true
    for required_column in required_columns
        if !( required_column in df_names )
            result = false
            break
        end
    end
    return result
end

"""
    highest_vote!( df::DataFrame, class_columns::Vector )::DataFrame

Finds highest vote amongst all specified class columns and creates a new column with that variable name.

"""
function highest_vote!( df::DataFrame, class_columns::Vector )::DataFrame
    df_classes = df[ !, class_columns ]
    df[ !, :Highest_Vote ] = map( row -> argmax( row ), eachrow( df_classes ) )
    df_classes = df[ !, vcat( class_columns, :Highest_Vote ) ]
    df[ !, :Highest_Prediction ] = map( row -> row[row.Highest_Vote], eachrow( df_classes ) )
    return df
end

highest_vote!( dfs::Vector{DataFrame}, class_columns::Vector )::Vector{DataFrame} = [ highest_vote!(df, class_columns) for df in dfs ]

"""
    barycenter_boxes!( df::DataFrame )::DataFrame

Will convert a DataFrame containing the following variables: center_x, center_y, width, height to be formatted in upper left lower right coordinates.

"""
function barycenter_boxes!( df::DataFrame )::DataFrame
    req_cols = [ :center_x, :center_y, :width, :height ]
    @assert( has_required_columns( df, req_cols ), "Input dataframe does not have the required columns: $req_cols. Please check spelling and capitalization of column names." )
    df[ !, [ :width, :height] ] ./= 2
    df[ !, :upper_left_x ] = df.center_x - df.width
    df[ !, :lower_right_x] = df.center_x + df.width
    df[ !, :upper_left_y ] = df.center_y - df.height
    df[ !, :lower_right_y] = df.center_y + df.height
    select!( df, Not( req_cols ) )
    return df
end

"""
    LLUR_boxes!( df::DataFrame )::DataFrame
"""
function LLUR_boxes!( df::DataFrame )::DataFrame
    req_cols = [ :lower_left_x, :upper_right_x, :lower_left_y, :upper_right_y ]
    new_cols = ULLR_required_columns
    @assert( has_required_columns( df, req_cols ), "Input dataframe does not have the required columns: $req_cols. Please check spelling and capitalization of column names." )
    rename!( df, Dict( req_cols .=> new_cols ) )
    return df
end

"""
    area( df::DataFrame )::DataFrame

Calculates the areas of all Boxes in a DataFrame as a new column.

"""
function area!( df::DataFrame )::DataFrame
    @assert( has_required_columns( df, ULLR_required_columns ), "Input dataframe does not have the required columns: $ULLR_required_columns. Please check spelling and capitalization of column names and ensure the format is ULLR." )
    df[ !, :Area ] = ( df.lower_right_x - df.upper_left_x ) .* ( df.lower_right_y - df.upper_left_y )
    return df
end

area!( dfs::Vector{DataFrame} )::Vector{DataFrame} = [ area!(df) for df in dfs ]

"""
    intersection_area( a::DataFrame, b::DataFrame )::DataFrame

Calculates the area of intersection between 2 DataFrames containing ULLR boxes.
This returns a new DataFrame

"""
function intersection_over_union( a::DataFrame, b::DataFrame )::DataFrame
    @assert( has_required_columns( a, ULLR_required_columns ), "Input dataframe does not have the required columns: $ULLR_required_columns. Please check spelling and capitalization of column names and ensure the format is ULLR." )
    @assert( has_required_columns( b, ULLR_required_columns ), "Input dataframe does not have the required columns: $ULLR_required_columns. Please check spelling and capitalization of column names and ensure the format is ULLR." )
    @assert( has_required_columns( a, [:Area] ), "Input dataframe does not have the required columns: Area. Please run `area()` on the individual dataframes." )
    @assert( has_required_columns( b, [:Area] ), "Input dataframe does not have the required columns: Area. Please run `area()` on the individual dataframes." )
    len_a, len_b = first.(size.( [ a, b ] ) )
    total_len = len_a * len_b

    output = DataFrame( (   :ID_A => Vector{Int}( undef, total_len ),
                            :ID_B => Vector{Int}( undef, total_len ),
                            :class_labels => Vector{Symbol}( undef, total_len ),
                            :Intersection_Area => Vector{Float64}( undef, total_len ),
                            :Intersection_Over_Union => Vector{Float64}( undef, total_len ),
                            :Vote_Matches => Vector{Bool}( undef, total_len ),
                            :Highest_Prediction => Vector{Float64}( undef, total_len),
                            :Highest_Vote => Vector{Symbol}( undef, total_len)
                        ) )
    current_row = 1
    for ( row_a_ID, row_a ) in enumerate( eachrow( a ) ), ( row_b_ID, row_b ) in enumerate( eachrow( b ) )
        #current_row = row_a_ID + ( len_a * ( row_b_ID - 1 ) )
        output[ current_row, :ID_A ] = row_a_ID
        output[ current_row, :ID_B ] = row_b_ID
        output[ current_row, :class_labels ] = row_b.class_labels
        output[ current_row, :Vote_Matches ] = row_a.Highest_Vote == row_b.class_labels
        output[ current_row, :Highest_Prediction ] = row_a.Highest_Prediction #row_b[!,row_a.Highest_Vote]
        output[ current_row, :Highest_Vote ] = row_a.Highest_Vote
        output[ current_row, :Intersection_Area ] = max( 0, min( row_a.lower_right_x, row_b.lower_right_x ) - max( row_a.upper_left_x, row_b.upper_left_x ) ) .*
                                                    max( 0, min( row_a.lower_right_y, row_b.lower_right_y ) - max( row_a.upper_left_y, row_b.upper_left_y ) )
        output[ current_row, :Intersection_Over_Union ] = output[ current_row, :Intersection_Area ] / ( row_a.Area + row_b.Area - output[ current_row, :Intersection_Area ])
        current_row += 1
    end
    #remove intermediate calculation column Intersection_Area
    select!(output, Not(:Intersection_Area))
    #Split data by each unique prediction
    #   -> sort by intersection over union (for each prediction all others are constant)
    #   -> grab the highest IoU and toss the rest.
    return combine( map( x -> sort(x, :Intersection_Over_Union, rev = true)[1,:], groupby( output, :ID_A ) ) )
end

mutable struct ObjectDetectionStats
    predictions::Vector{DataFrame}
    groundtruth::Vector{DataFrame}
    IoU_map::Vector{DataFrame}
    classnames::Vector{Symbol}
    len::Int
end

"""
    potential_scores( ods::ObjectDetectionStats; IoU_threshold = 0.5 )

"""
function potential_scores( cur_IoU_df::DataFrame; IoU_threshold = 0.5 )
    cur_IoU_df[ :iou_threshold ] = cur_IoU_df.Intersection_Over_Union .> IoU_threshold
    #is TP or FN?
    cur_IoU_df[ :TP_or_FN ] = cur_IoU_df.Vote_Matches .& cur_IoU_df.iou_threshold
    #Ensure we do not double dip from a GT prediction. Only one prediction can map to a GT here...
    #Group by ground truth class label
    #   -> sort by the highest predicted class score
    #       _-> find first prediced class matches above an IoU threshold
    for gdf in groupby( cur_IoU_df, :class_labels )
        sort!( cur_IoU_df, :Highest_Prediction, rev = true)
        kept_first = false
        for row in eachrow( gdf )
            if row.TP_or_FN
                if kept_first
                    row.TP_or_FN = false
                end
                kept_first = true
            end
        end
    end
    #Bifurcate DF on the possibility of true positives
    TP_or_FN, FP_or_FN = DataFrame.( collect( groupby( cur_IoU_df, :TP_or_FN ) ) )
    return TP_or_FN, FP_or_FN
end

"""
    ObjectDetectionStats(   predictions::Vector{DataFrame},
                            groundtruth::Vector{DataFrame},
                            classnames::Vector{Symbol} )::ObjectDetectionStats

Instantiates an `ObjectDetectionStats` object.

"""
function ObjectDetectionStats(  predictions::Vector{DataFrame},
                                groundtruth::Vector{DataFrame},
                                classnames::Vector{Symbol} )::ObjectDetectionStats
    len = length(predictions)
    @assert len == length(groundtruth) "Prediction and groundtruth Vector's must have the same length."
    area!( predictions )
    highest_vote!( predictions, classnames )
    area!( groundtruth )
    IoU_map = [ intersection_over_union( p, gt ) for (p, gt) in zip( predictions, groundtruth) ]
    return ObjectDetectionStats( predictions, groundtruth, IoU_map, classnames, len )
end

# each prediction maps to every GT
# predictions sorted by their confidence scores (highest to lowest)
# Each prediction can ONLY map to one GT
#  for each fixed prediction -> [GT] the prediction score is static
#  Need to find "most likely hit" via IoU.
#So sort by confidence and IoU
# Now for each group of predictions, remove everything but the top ranked ones
#Now we have the most representative prediction for each prediction to [GT]
#Find TP-or-FN mappings
#Find FP-or-FN mappings

using DataFrames
a = DataFrame(Dict([ :a => [1,1,1,2,2,2,3,3,3], :b => randn(9)]))
x,z,t = DataFrame.( collect( groupby(a, :a) ) )

a

KeyDF = DataFrame( Dict( [    :a               => [0.8, 0.1, 0.1],
                              :b               => [0.1, 0.7, 0.1],
                              :c               => [0.1, 0.2, 0.8],
                              :upper_left_x     => [12.5, 10.0 ,20.0],
                              :lower_right_x    => [17.5, 20.0, 30.0],
                              :upper_left_y     => [7.50, 5.00, 30.0],
                              :lower_right_y    => [12.5, 15.0, 35.0]
                         ] ) )

GT = DataFrame( Dict( [     :class_labels     => [:a, :b, :c],
                                :upper_left_x     => [12.5, 10.0, 16.0],
                                :lower_right_x    => [17.5, 20.0, 60.0],
                                :upper_left_y     => [7.50, 5.00, 16.0],
                                :lower_right_y    => [12.5, 15.0, 20.0]
                      ] ) )

ods = ObjectDetectionStats(  [KeyDF], [GT], [:a,:b,:c] )

score( ods, detection_threshold = 0.5, IoU_threshold = 0.5 )

@show ods.IoU_map[1]

#@benchmark join( KeyDF, KeyDFB, kind = :cross; makeunique = true )
#new_df = intersection_over_union( area( KeyDF ), area( KeyDFB ) )
@show( new_df )

# Pkg.add("DataFrames")
# using DataFrames
#
# b = DataFrame(  (:A => [1,2,1,2,2,2,3,3,3], :B => collect(1:9) )  )
#
# for gdf in groupby( b, :A )
#     gdf.B = map(x -> (x.B % 2) == 1, eachrow(gdf))
#     println(gdf)
# end
#
# b


text:
using FreqTables
using Base.Test

x = repeat(["a", "b", "c", "d"], outer=[100]);
y = repeat(["A", "B", "C", "D"], inner=[10], outer=[10]);

@test freqtable(x).array == [100, 100, 100, 100]
@test freqtable(y).array == [100, 100, 100, 100]
@test freqtable(x, y).array == [30 20 30 20;
                                30 20 30 20;
                                20 30 20 30;
                                20 30 20 30]

@test freqtable(x, y,
                subset=1:20,
                weights=repeat([1, .5], outer=[10])).array == [3.0 2.0
                                                               1.5 1.0
                                                               2.0 3.0
                                                               1.0 1.5]

using DataArrays
xpda = PooledDataArray(x)
ypda = PooledDataArray(y)

@test freqtable(xpda).array == [100, 100, 100, 100]
@test freqtable(ypda).array == [100, 100, 100, 100]
@test freqtable(xpda, ypda).array == [30 20 30 20;
                                      30 20 30 20;
                                      20 30 20 30;
                                      20 30 20 30]

xpda[1] = NA
ypda[[1, 10, 20, 400]] = NA

@test freqtable(xpda).array == [99, 100, 100, 100]
@test freqtable(ypda).array == [98, 99, 100, 99]
@test freqtable(xpda, ypda).array == [29 20 30 20;
                                      29 20 30 20;
                                      20 30 20 30;
                                      20 29 20 29]

@test freqtable(xpda, usena=true).array == [99, 100, 100, 100, 1]
@test freqtable(ypda, usena=true).array == [98, 99, 100, 99, 4]
@test freqtable(xpda, ypda, usena=true).array == [29 20 30 20 0;
                                                  29 20 30 20 1;
                                                  20 30 20 30 0;
                                                  20 29 20 29 2;
                                                  0 0 0 0 1]

using RDatasets
iris = dataset("datasets", "iris")
iris[:LongSepal] = iris[:SepalLength] .> 5.0
@test freqtable(iris, :Species, :LongSepal).array == [28 22
                                                       3 47
                                                       1 49]
@test freqtable(iris, :Species, :LongSepal,
                subset=iris[:PetalLength] .< 4.0).array ==[28 22
                                                            3  8]

# Issue #5
@test freqtable([Set(1), Set(2)]).array == [1, 1]
@test freqtable([Set(1), Set(2)], [Set(1), Set(2)]).array == eye(2)


text:
using Gaugefields
using MPI

const pes = Tuple(parse.(Int64,ARGS[1:4]))
const mpi = parse(Bool,ARGS[5])

function test()
    NX = 8*2
    NY = 8*2
    NZ = 8*2
    NT = 8*2
    Nwing = 0
    NC = 3

    Dim = 4

    if mpi
        PEs = pes#(1,1,1,2)
        U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot",mpi=true,PEs = PEs,mpiinit = false)

    else
        U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot")

    end

    temp1 = similar(U[1])
    temp2 = similar(U[1])
    temp3 = similar(U[1])

    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    β = 5.7

    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println(" plaq_t = $plaq_t")

    g = Gradientflow(U,eps = 0.01,mpi=mpi)

    for itrj=1:100
        @time flow!(U,g)
        if itrj % 10 == 0
            plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
            println("$itrj plaq_t = $plaq_t")
            #poly = calculate_Polyakov_loop(U,temp1,temp2)
            #println("$itrj polyakov loop = $(real(poly)) $(imag(poly))")
        end
    end

end
test()


text:
abstract type RDD end
abstract type SingleRDD <: RDD end
abstract type PairRDD <: RDD end

"Pure wrapper around JavaRDD"
mutable struct JavaRDD <: SingleRDD
    jrdd::JJavaRDD
end

"Pure wrapper around JavaPairRDD"
mutable struct JavaPairRDD <: PairRDD
    jrdd::JJavaPairRDD
end

"""
Julia type to handle RDDs. Can handle pipelining of operations to reduce interprocess IO.
"""
mutable struct PipelinedRDD <: SingleRDD
    parentrdd::RDD
    func::Function
    jrdd::JJuliaRDD
end

"""
Julia type to handle Pair RDDs. Can handle pipelining of operations to reduce interprocess IO.
"""
mutable struct PipelinedPairRDD <: PairRDD
    parentrdd::RDD
    func::Function
    jrdd::JJuliaPairRDD
end

mutable struct RDDIterator
    itr::JavaObject{Symbol("java.util.Iterator")}
    l::jint
end

function Base.iterate(x::RDDIterator, state=iterate(x.itr) )
    if state == nothing
        return nothing
    end
    value, itr_state = state
    b = jcall(Spark.JJuliaRDD, "writeToByteArray", Vector{jbyte}, (JObject,), value)
    b=reinterpret(Vector{UInt8}, b)
    return readobj(IOBuffer(b))[2], iterate(x.itr, itr_state)
end

Base.IteratorSize(::Type{RDDIterator}) = Base.HasLength()
Base.length(x::RDDIterator) = x.l
Base.IteratorEltype(::Type{RDDIterator}) = Base.EltypeUnknown()

"""
Params:
 * parentrdd - parent RDD
 * func - function of type `(index, iterator) -> iterator` to apply to each partition
"""
function PipelinedRDD(parentrdd::RDD, func::Function)
    pipelined_func = create_pipeline_command(parentrdd, func)
    command_ser = reinterpret(Vector{jbyte}, serialized(pipelined_func))
    jrdd = jcall(JJuliaRDD, "fromRDD", JJuliaRDD,
                 (JRDD, Vector{jbyte}),
                 get_root_rdd(parentrdd), command_ser)
    PipelinedRDD(parentrdd, pipelined_func, jrdd)
end

"""
Params:
 * parentrdd - parent RDD
 * func - function of type `(index, iterator) -> iterator` to apply to each partition
"""
function PipelinedPairRDD(parentrdd::RDD, func::Function)
    pipelined_func = create_pipeline_command(parentrdd, func)
    command_ser = reinterpret(Vector{jbyte}, serialized(pipelined_func))
    jrdd = jcall(JJuliaPairRDD, "fromRDD", JJuliaPairRDD,
                 (JRDD, Vector{jbyte}),
                 get_root_rdd(parentrdd), command_ser)
    PipelinedPairRDD(parentrdd, pipelined_func, jrdd)
end

# as_java_rdd returns a JJavaRDD or JJavaPairRDD class
as_java_rdd(rdd::RDD) = rdd.jrdd
as_java_rdd(rdd::PipelinedRDD) = jcall(rdd.jrdd, "asJavaRDD", JJavaRDD, ())
as_java_rdd(rdd::PipelinedPairRDD) = jcall(rdd.jrdd, "asJavaPairRDD", JJavaPairRDD, ())

# as_java_rdd returns a JRDD class
as_rdd(rdd::RDD) = jcall(rdd.jrdd, "rdd", JRDD, ())
as_rdd(rdd::PipelinedRDD) = rdd.jrdd
as_rdd(rdd::PipelinedPairRDD) = rdd.jrdd

# as_root_rdd returns the root RDD suitable for contructing JJuliaRDD or JJuliaPairRDD objects
get_root_rdd(rdd::RDD) = as_rdd(rdd)
get_root_rdd(rdd::PipelinedRDD) = get_root_rdd(rdd.parentrdd)
get_root_rdd(rdd::PipelinedPairRDD) = get_root_rdd(rdd.parentrdd)

Base.show(io::IO, rdd::JavaRDD) = print(io, "JavaRDD()")
Base.show(io::IO, rdd::JavaPairRDD) = print(io, "JavaPairRDD()")
Base.show(io::IO, rdd::PipelinedRDD) =  print(io, "PipelinedRDD($(rdd.parentrdd))")
Base.show(io::IO, rdd::PipelinedPairRDD) = print(io, "PipelinedPairRDD($(rdd.parentrdd))")

" chain 2 partion functions together "
function chain_function(parent_func, child_func)
    function pipelined_func(split, iterator)
        return child_func(split, parent_func(split, iterator))
    end
    pipelined_func
end

create_pipeline_command(rdd::RDD, func) = func
create_pipeline_command(rdd::PipelinedRDD, func) = chain_function(rdd.func, func)
create_pipeline_command(rdd::PipelinedPairRDD, func) = chain_function(rdd.func, func)

Base.reinterpret(::Type{Array{jbyte,1}}, bytes::Array{UInt8,1}) =
    jbyte[reinterpret(jbyte, b) for b in bytes]

Base.reinterpret(::Type{Array{UInt8,1}}, bytes::Array{jbyte,1}) =
    UInt8[reinterpret(UInt8, b) for b in bytes]

"""
Apply function `f` to each partition of `rdd`. `f` should be of type
`(index, iterator) -> iterator`
"""
function map_partitions_with_index(rdd::RDD, f::Function)
    return PipelinedRDD(rdd, f)
end

function add_index_param(f::Function)
    function func(idx, it)
        f(it)
    end
    func
end

"""
Apply function `f` to each partition of `rdd`. `f` should be of type
`(iterator) -> iterator`
"""
function map_partitions(rdd::RDD, f::Function)
    return PipelinedRDD(rdd, add_index_param(f))
end

"""
Apply function `f` to each partition of `rdd`. `f` should be of type
`(iterator) -> iterator`
"""
function map_partitions_pair(rdd::RDD, f::Function)
    return PipelinedPairRDD(rdd, add_index_param(f))
end

"""
creates a function that operates on a partition from an
element by element map function
"""
function create_map_function(f::Function)
    function func(idx, it)
        (f(i) for i in it)
    end
    return func
end

"Apply function `f` to each element of `rdd`"
function map(rdd::RDD, f::Function)
    return PipelinedRDD(rdd, create_map_function(f))
end

"Apply function `f` to each element of `rdd`"
function map_pair(rdd::RDD, f::Function)
    return PipelinedPairRDD(rdd, create_map_function(f))
end

"""
Similar to `map`, but each input item can be mapped to 0 or more
output items (so `f` should return an iterator rather than a single item)
"""
function flat_map(rdd::RDD, f::Function)
    # return PipelinedRDD(rdd, create_flat_map_function(f))
    return PipelinedRDD(rdd, (idx, it) -> Iterators.flatten(Iterators.map(f, it)))
end

"""
Similar to `map`, but each input item can be mapped to 0 or more
output items (so `f` should return an iterator of pairs rather than a single item)
"""
function flat_map_pair(rdd::RDD, f::Function)
    return PipelinedPairRDD(rdd, create_flat_map_function(f))
end

function create_filter_function(f::Function)
    function func(idx, it)
        Iterators.filter(f, it)
    end
    return func
end

Base.filter(rdd::SingleRDD, f::Function) = PipelinedRDD(rdd, create_filter_function(f))
Base.filter(rdd::PairRDD, f::Function) = PipelinedPairRDD(rdd, create_filter_function(f))

"Reduce elements of `rdd` using specified function `f`"
function reduce(rdd::RDD, f::Function)
    process_attachments(context(rdd))
    locally_reduced = map_partitions(rdd, it->reduction_function(f, it))
    subresults = collect(locally_reduced)
    subresults = collect_itr(locally_reduced)
    return reduce(f, subresults)
end

function reduction_function(f, it)
    try
        return [reduce(f, it)]
    catch
        return []
    end
end

"Get SparkContext of this RDD"
function context(rdd::RDD)
    ssc = jcall(rdd.jrdd, "context", JSparkContext, ())
    jsc = jcall(JJavaSparkContext, "fromSparkContext",
                JJavaSparkContext, (JSparkContext,), ssc)
    return SparkContext(jsc)
end

"""Collects the RDD to the Julia process, by serialising all values
via a byte array"""
function collect_internal(rdd::RDD, static_java_class, result_class)
    process_attachments(context(rdd))
    jbyte_arr = jcall(static_java_class, "collectToJulia", Vector{jbyte},
                 (result_class,),
                 as_java_rdd(rdd))

    byte_arrs = reinterpret(Vector{UInt8}, jbyte_arr)
    val = readobj(IOBuffer(byte_arrs))[2]
    return val
end

"""Collects the RDD to the Julia process, via an Julia iterator
that fetches each row at a time. This prevents creation of a byte
array containing all rows at a time."""
function collect_internal_itr(rdd::RDD, static_java_class, result_class)
    process_attachments(context(rdd))
    lst = jcall(static_java_class, "collectToJuliaItr", @jimport(java.util.List),
                 (result_class,),
                 as_java_rdd(rdd))
    itr = jcall(lst, "iterator", @jimport(java.util.Iterator), ())
    l = jcall(lst, "size", jint, ())
    return RDDIterator(itr, l)
end

"""
Collect all elements of `rdd` on a driver machine
"""
function collect(rdd::SingleRDD)
    collect_internal(rdd, JJuliaRDD, JJavaRDD)
end

"""
Collect all elements of `rdd` on a driver machine
"""
function collect_itr(rdd::PairRDD)
    collect_internal_itr(rdd, JJuliaPairRDD, JJavaPairRDD)
end

"""
Collect all elements of `rdd` on a driver machine
"""
function collect_itr(rdd::SingleRDD)
    collect_internal_itr(rdd, JJuliaRDD, JJavaRDD)
end

"""
Collect all elements of `rdd` on a driver machine
"""
function collect(rdd::PairRDD)
    collect_internal(rdd, JJuliaPairRDD, JJavaPairRDD)
end

"Count number of elements in this RDD"
function count(rdd::RDD)
    process_attachments(context(rdd))
    return jcall(rdd.jrdd, "count", jlong, ())
end

"Persist this RDD with the default storage level (MEMORY_ONLY)"
function cache(rdd::SingleRDD)
    JavaRDD(jcall(as_java_rdd(rdd), "cache", JJavaRDD, ()))
end

"Persist this RDD with the default storage level (MEMORY_ONLY)"
function cache(rdd::PairRDD)
    JavaPairRDD(jcall(as_java_rdd(rdd), "cache", JJavaPairRDD, ()))
end

"Create a pair RDD with every combination of the values of rdd1 and rdd2"
function cartesian(rdd1::SingleRDD, rdd2::SingleRDD)
    jprdd = jcall(JJuliaRDD, "cartesianSS", JJavaPairRDD,
                 (JJavaRDD, JJavaRDD),
                 as_java_rdd(rdd1),  as_java_rdd(rdd2))

    return JavaPairRDD(jprdd)
end

"When called on a dataset of (K, V) pairs, returns a dataset of (K, [V]) pairs."
function group_by_key(rdd::PairRDD)
    jprdd = jcall(as_java_rdd(rdd), "groupByKey", JJavaPairRDD, ())
    return JavaPairRDD(jprdd)
end

"""Return a new RDD that has exactly num_partitions partitions."""
function repartition(rdd::T, num_partitions::Integer) where {T<:RDD}
    (Tjr, Tr) = (T <: PairRDD) ? (JJavaPairRDD, JavaPairRDD) : (JJavaRDD, JavaRDD)
    jrdd = jcall(as_java_rdd(rdd), "repartition", Tjr, (jint,), num_partitions)
    return Tr(jrdd)
end

"""Return a new RDD that is reduced into num_partitions partitions."""
function coalesce(rdd::T, num_partitions::Integer; shuffle::Union{Nothing,Bool}=nothing) where {T<:RDD}
    (Tjr, Tr) = (T <: PairRDD) ? (JJavaPairRDD, JavaPairRDD) : (JJavaRDD, JavaRDD)
    if shuffle === nothing
        jrdd = jcall(as_java_rdd(rdd), "coalesce", Tjr, (jint,), num_partitions)
    else
        jrdd = jcall(as_java_rdd(rdd), "coalesce", Tjr, (jint,jboolean), num_partitions, shuffle)
    end
    return Tr(jrdd)
end

"""
When called on a dataset of (K, V) pairs, returns a dataset of (K, V) pairs where the
values for each key are aggregated using the given reduce function func,
which must be of type (V,V) => V.
"""
function reduce_by_key(rdd::PairRDD, f::Function)
    grouped = group_by_key(rdd)
    function func(it)
        (it[1], reduce(f, it[2]))
    end
    return map_pair(grouped, func)
end

"""Returns the number of partitions of this RDD."""
num_partitions(rdd::Union{PipelinedRDD,PipelinedPairRDD}) = jcall(rdd.jrdd, "getNumPartitions", jint, ())
num_partitions(rdd::JavaRDD) = jcall(JRDDUtils, "getNumPartitions", jint, (JJavaRDD,), as_java_rdd(rdd))
num_partitions(rdd::JavaPairRDD) = jcall(JRDDUtils, "getNumPartitions", jint, (JJavaPairRDD,), as_java_rdd(rdd))

"Return the id of the rdd"
function id(rdd::RDD)
    jcall(rdd.jrdd, "id", jint, ())
end

"""Return an RDD created by piping elements to a forked external process."""
function pipe(rdd::RDD, command::String)
    jrdd = jcall(as_java_rdd(rdd), "pipe", JJavaRDD, (JString,), command)
    return JavaRDD(jrdd)
end

function pipe(rdd::RDD, command::Vector{String})
    jrdd = jcall(as_java_rdd(rdd), "pipe", JJavaRDD, (JList,), convert(JArrayList, command, JString))
    return JavaRDD(jrdd)
end

function pipe(rdd::RDD, command::Vector{String}, env::Dict{String,String})
    jrdd = jcall(as_java_rdd(rdd), "pipe", JJavaRDD, (JList,JMap), convert(JArrayList, command, JString), convert(JHashMap, JString, JString, env))
    return JavaRDD(jrdd)
end


text:
struct Unit{S, T} # S = suffix symbol
    value::T
end
Unit{S}(value) where {S} = Unit{S, typeof(value)}(value)

Base.show(io::IO, x::Unit{S}) where {S} = print(io, x.value, S)

# diagonal dispatch for unit + unit, unit - unit
Base.:+(x::Unit{S}, y::Unit{S}) where {S} = Unit{S}(x.value + y.value)
Base.:-(x::Unit{S}, y::Unit{S}) where {S} = Unit{S}(x.value - y.value)

# scalar * unit, unit / scalar
Base.:*(x::Number, y::Unit{S}) where {S} = Unit{S}(x * y.value)
Base.:/(x::Unit{S}, y::Number) where {S} = Unit{S}(x.value / y)

# calc() expressions
struct Calc
    expr::String
    Calc(expr) = new("($expr)")
end
Base.show(io::IO, x::Calc) = print(io, "calc", x.expr)

# default to calc() for mismatched units
Base.:+(x::Unit,   y::Unit) = Calc("$x + $y")
Base.:-(x::Unit,   y::Unit) = Calc("$x - $y")

# unit + calc(), calc() + unit
Base.:+(x::Unit,   y::Calc) = Calc("$x + $(y.expr)")
Base.:+(x::Calc,   y::Unit) = Calc("$(x.expr) + $y")

# unit - calc(), calc() - unit,
Base.:-(x::Unit,   y::Calc) = Calc("$x - $(y.expr)")
Base.:-(x::Calc,   y::Unit) = Calc("$(x.expr) - $y")

# scalar * calc(), calc() / scalar
Base.:*(x::Number, y::Calc) = Calc("$x * $(y.expr)")
Base.:/(x::Calc, y::Number) = Calc("$(x.expr) / $y")

# concise scalar * unit construction, e.g. 5px
struct BareUnit{T} end
Base.:*(x::T, ::BareUnit{U}) where {U <: Unit, T<:Number} = U{T}(x)

# common css units (ex, ch excluded)
const px = BareUnit{Unit{:px}}()
const pt = BareUnit{Unit{:pt}}()
const em = BareUnit{Unit{:em}}()
const rem = BareUnit{Unit{:rem}}()
const vh = BareUnit{Unit{:vh}}()
const vw = BareUnit{Unit{:vw}}()
const vmin = BareUnit{Unit{:vmin}}()
const vmax = BareUnit{Unit{:vmax}}()
const pc = BareUnit{Unit{Symbol("%")}}()


text:
using Documenter, SurrogateModelOptim

makedocs(;
        modules = [SurrogateModelOptim],
        format = Documenter.HTML(   assets = ["assets/favicon.ico"],
                                prettyurls = get(ENV, "CI", nothing) == "true"),
        sitename = "SurrogateModelOptim.jl",
        strict = true,
        clean = true,
        checkdocs = :none,
        pages = [
                "Home" => "index.md",
                "Manual" => [
                        "smoptimize.md",
                        "surrogate_model.md",
                        "model_infill.md",
                        "options.md",
                        "examples.md",
                        "benchmark.md",
                ]
        ]
)

deploydocs(
    repo = "github.com/MrUrq/SurrogateModelOptim.jl.git"
)


text:
function collapse_data(x::Vector{DataFrame})
    base = copy(x[1])
    refnames = Set(names(base))

    # The first instance of each key is always taken. For the sake of speed, we
    # don't bother checking for equality or any of that crap.
    for i in 2:length(x)
        current = x[i]
        for n in names(current)
            if !(n in refnames)
                insertcols!(base, n => current[!,n])
            end
        end
    end

    return base
end

function collapse_data(x::Vector{Dict{String,Any}})
    base = copy(x[1])

    # The first instance of each key is always taken. For the sake of speed, we
    # don't bother checking for equality or any of that crap.
    for i in 2:length(x)
        current = x[i]
        for (key, val) in current
            if !haskey(base, key)
                base[key] = val
            end
        end
    end

    return base
end

"""
    vcat(A::Vararg{SummarizedExperiment})

Vertically concatenate one or more `SummarizedExperiment` objects.
The input objects must satisfy the following constraints:

- All objects must have the same number of columns, which are assumed to be in the same order.
- All objects must have `DataFrame`s in their `rowdata` with the same type and names of all columns (though they may be ordered differently).
- All objects must have the same names and types of assays;
  for a given assay name, the dimensions of the corresponding arrays across all `A` should be the same except for the first dimension.

This function returns a single `SummarizedExperiment` instance where the number of rows is equal to the sum of the number of rows across all objects in `A`.
The number of columns in the output object is the same as the number of columns in any object in `A`.
The order of columns in the output `rowdata` is the same as that of the first object.
The output `coldata` is created by combining columns horizontally across `coldata` of all objects in `A`;
if columns have duplicate names, only the first instance of each column is retained.

# Examples
```jldoctest
julia> using SummarizedExperiments

julia> x = exampleobject(20, 10);

julia> y = exampleobject(30, 10);

julia> z = vcat(x, y)
50x10 SummarizedExperiment
  assays(3): foo bar whee
  rownames: Gene1 Gene2 ... Gene29 Gene30
  rowdata(2): name Type
  colnames: Patient1 Patient2 ... Patient9 Patient10
  coldata(3): name Treatment Response
  metadata(1): version
```
"""
function Base.vcat(A::Vararg{SummarizedExperiment})
    if length(A) < 2
        return A[1] # no-op in this case.
    end

    collected_rowdata = [rowdata(A[1])];
    collected_coldata = [coldata(A[1])];
    collected_metadata = [metadata(A[1])];
    collected_assays = OrderedDict{String,Vector{AbstractArray}}();
    for (key, val) in assays(A[1])
        collected_assays[key] = AbstractArray[val]
    end

    for i in 2:length(A)
        push!(collected_rowdata, rowdata(A[i]))
        push!(collected_coldata, coldata(A[i]))

        for (key, val) in assays(A[i])
            if !haskey(collected_assays, key)
                throw(KeyError("could not find '" * key * "' in assays of object " * string(i)))
            end

            ref = collected_assays[key][1]
            if length(size(ref)) != length(size(val)) || size(ref)[2] != size(val)[2]
                throw(DimensionMismatch("inconsistent dimensions for assay '" * key * "' in object " * string(i)))
            end

            push!(collected_assays[key], val)
        end
    end

    full_rowdata = vcat(collected_rowdata...)
    full_assays = OrderedDict{String,AbstractArray}()
    for (key, val) in collected_assays
        full_assays[key] = vcat(val...)
    end

    full_coldata = collapse_data(collected_coldata)
    full_metadata = collapse_data(collected_metadata)

    return SummarizedExperiment(full_assays, full_rowdata, full_coldata, full_metadata)
end

"""
    hcat(A::Vararg{SummarizedExperiment})

Horizontally concatenate one or more `SummarizedExperiment` objects.
The input objects must satisfy the following constraints:

- All objects must have the same number of rows, which are assumed to be in the same order.
- All objects must have `DataFrame`s in their `coldata` with the same type and names of all columns (though they may be ordered differently).
- All objects must have the same names and types of assays;
  for a given assay name, the dimensions of the corresponding arrays across all `A` should be the same except for the second dimension.

This function returns a single `SummarizedExperiment` instance where the number of columns is equal to the sum of the number of columns across all objects in `A`.
The number of rows in the output object is the same as the number of rows in any object in `A`.
The order of columns in the output `coldata` is the same as that of the first object.
The output `rowdata` is created by combining columns horizontally across `rowdata` of all objects in `A`;
if columns have duplicate names, only the first instance of each column is retained.

# Examples
```jldoctest
julia> using SummarizedExperiments

julia> x = exampleobject(20, 20);

julia> y = exampleobject(20, 30);

julia> z = hcat(x, y)
20x50 SummarizedExperiment
  assays(3): foo bar whee
  rownames: Gene1 Gene2 ... Gene19 Gene20
  rowdata(2): name Type
  colnames: Patient1 Patient2 ... Patient29 Patient30
  coldata(3): name Treatment Response
  metadata(1): version
```
"""
function Base.hcat(A::Vararg{SummarizedExperiment})
    if length(A) < 2
        return A[1] # no-op in this case.
    end

    collected_rowdata = [rowdata(A[1])];
    collected_coldata = [coldata(A[1])];
    collected_metadata = [metadata(A[1])];
    collected_assays = OrderedDict{String,Vector{AbstractArray}}();
    for (key, val) in assays(A[1])
        collected_assays[key] = AbstractArray[val]
    end

    for i in 2:length(A)
        push!(collected_rowdata, rowdata(A[i]))
        push!(collected_coldata, coldata(A[i]))

        for (key, val) in assays(A[i])
            if !haskey(collected_assays, key)
                throw(KeyError("could not find '" * key * "' in assays of object " * string(i)))
            end

            ref = collected_assays[key][1]
            if length(size(ref)) != length(size(val)) || size(ref)[1] != size(val)[1]
                throw(DimensionMismatch("inconsistent dimensions for assay '" * key * "' in object " * string(i)))
            end

            push!(collected_assays[key], val)
        end
    end

    full_coldata = vcat(collected_coldata...)
    full_assays = OrderedDict{String,AbstractArray}()
    for (key, val) in collected_assays
        full_assays[key] = hcat(val...)
    end

    full_rowdata = collapse_data(collected_rowdata)
    full_metadata = collapse_data(collected_metadata)

    return SummarizedExperiment(full_assays, full_rowdata, full_coldata, full_metadata)
end


text:
using CompScienceMeshes
using BEAST
using LinearAlgebra
using Test

# Γ = meshcuboid(1.0, 1.0, 1.0, 0.25)
# Γ = meshsphere(1.0, 0.3)
Γ = meshrectangle(1.0, 1.0, 1/3, 3)
Γ1 = CompScienceMeshes.rotate(Γ, pi/2*point(1,0,0))
Γ2 = CompScienceMeshes.weld(Γ,Γ1)

sol = 1.0
SL0 = TDMaxwell3D.singlelayer(speedoflight=sol, numdiffs=0)
SL1 = TDMaxwell3D.singlelayer(speedoflight=sol, numdiffs=1)
DL0 = TDMaxwell3D.doublelayer(speedoflight=sol, numdiffs=0)

Δt, Nt = 0.6, 80
duration = 20 * Δt
delay = 1.5 * duration
amplitude = 1.0
gaussian = creategaussian(duration, delay, amplitude)
direction, polarisation = ẑ, x̂
E = BEAST.planewave(polarisation, direction, gaussian, 1.0)
E1 = BEAST.planewave(polarisation, direction, derive(gaussian), 1.0)

δ = timebasisdelta(Δt, Nt)
T0 = timebasisshiftedlagrange(Δt, Nt, 0)
T1 = timebasisshiftedlagrange(Δt, Nt, 1)
T2 = timebasisshiftedlagrange(Δt, Nt, 2)
T3 = timebasisshiftedlagrange(Δt, Nt, 3)

iT1 = integrate(T1)
iT2 = integrate(T2)

X = raviartthomas(Γ)
Y = buffachristiansen(Γ, ibscaled=false)

fr1, store1 = BEAST.allocatestorage(SL0, X⊗δ, X⊗T1,
    BEAST.Val{:bandedstorage},
    BEAST.LongDelays{:compress})
fr2, store2 = BEAST.allocatestorage(SL0, X⊗δ, X⊗T1,
    BEAST.Val{:bandedstorage},
    BEAST.LongDelays{:ignore})

BEAST.assemble!(SL0, X⊗δ, X⊗T1, store1); Z1 = fr1()
BEAST.assemble!(SL0, X⊗δ, X⊗T1, store2); Z2 = fr2()

@test norm(Z1-Z2,Inf) < 1e-12

fr3, store7 = BEAST.allocatestorage(SL1, X⊗δ, X⊗T2,
    BEAST.Val{:bandedstorage},
    BEAST.LongDelays{:compress})
fr4, store8 = BEAST.allocatestorage(SL1, X⊗δ, X⊗T2,
    BEAST.Val{:bandedstorage},
    BEAST.LongDelays{:ignore})

BEAST.assemble!(SL1, X⊗δ, X⊗T2, store7); Z3 = fr3()
BEAST.assemble!(SL1, X⊗δ, X⊗T2, store8); Z4 = fr4()

# x = rand(size(Z1)[1:2]...)
# csx = cumsum(x, dims = 2)
#
# j, k_start = 10, 2
# y1 = zeros(eltype(Z1), size(Z1,1))
# y2 = zeros(eltype(Z2), size(Z2,1))
# BEAST.convolve!(y1, Z1, x, csx, j, k_start)
# y2 = convolve(Z2, x, j, k_start)
# @show norm(y2 - y1)

b = assemble(E, X⊗δ)
b1 = assemble(E1, X⊗δ)

W1 = inv(Z1[:,:,1])
W2 = inv(Z2[:,:,1])
@test norm(W1-W2,Inf) < 1e-12

W3 = inv(Z3[:,:,1])
W4 = inv(Z4[:,:,1])
@test norm(W3-W4) < 1e-12

x1 = marchonintime(W1, Z1, b, Nt)
x2 = marchonintime(W1, Z2, b, Nt)
@test norm(x1-x2,Inf) < 1e-12

x3 = marchonintime(W3, Z3, b1, Nt)
x4 = marchonintime(W4, Z4, b1, Nt)
@test norm(x3-x4,Inf) < 1e-12

@test norm(vec(x1)-vec(x3)) / norm(vec(x3)) < 0.1
@test norm(vec(x1)-vec(x4)) / norm(vec(x4)) < 0.1

X2 = raviartthomas(Γ2)
Y2 = raviartthomas(Γ2)
DLh = (DL0, X2⊗δ, X2⊗iT2)

fr9, store9 = BEAST.allocatestorage(DLh...,
    BEAST.Val{:bandedstorage},
    BEAST.LongDelays{:compress})
fr10, store10 = BEAST.allocatestorage(DLh...,
    BEAST.Val{:bandedstorage},
    BEAST.LongDelays{:ignore})

@time BEAST.assemble!(DLh..., store9)
@time BEAST.assemble!(DLh..., store10)

Z9 = fr9()
Z10 = fr10()

@test norm(Z9-Z10,Inf) < 1e-12

iDLh = (integrate(DL0), Y2⊗δ, X2⊗T1)
fr11, store11 = BEAST.allocatestorage(iDLh..., BEAST.Val{:densestorage}, BEAST.LongDelays{:ignore})
fr12, store12 = BEAST.allocatestorage(iDLh..., BEAST.Val{:bandedstorage}, BEAST.LongDelays{:compress})
fr13, store13 = BEAST.allocatestorage(iDLh..., BEAST.Val{:bandedstorage}, BEAST.LongDelays{:ignore})

BEAST.assemble!(iDLh..., store11); Z11 = fr11()
BEAST.assemble!(iDLh..., store12); Z12 = fr12()
BEAST.assemble!(iDLh..., store13); Z13 = fr13()

@test norm(Z11-Z12,Inf) < 1e-8
@test norm(Z11-Z13,Inf) < 1e-8


text:
# Reexport SinkingParticles as they are useful outside too
@reexport module SinkingParticles

using Unitful
using LinearAlgebra, SparseArrays
using ..GridTools

"""
    buildPFD(w, DIV, Iabove)

Builds the particle flux divergence operator `PFD`
for a particle sinking speed w.

Schematic of a grid cell:
```
     top  ┌─────────────────────────────┐   ┬
          │      ↓ w            ↓ Φ     │   │
          │ (sinking speed)    (flux)   │   │
          │                             │   │
          │            x                │   dz
          │      (particle conc.)       │   │
          │                             │   │
          │                             │   │
  bottom  └─────────────────────────────┘   ┴
```

- dz is the height of grid cell [m]
- w is the particle sinking speed at the top of the grid cell [m s⁻¹]
- Φ is the flux at the top of the grid cell [mmol m⁻² s⁻¹]
- x is the particle concentration of the cell [mmol m⁻³]

The PFD operator is defined by

    PFD * x ≃ dΦ/dz,

i.e., so that applied to x it approximates dΦ/dz.
It is calculated as the matrix product of the DIV and FLUX operators:

    PFD = DIV * FLUX.

DIV, the divergence operator (d/dz only here), is defined by

    DIV * Φ = 1/dz * (Φ - Φ[below]) ≃ dΦ/dz.

And FLUX, the sinking particle flux operator, is defined by

    FLUX * x = w * x[above] ≃ Φ.

# Example

```julia-repl
julia> PFD = buildPFD(w)
```
"""
buildPFD(w, DIV, Iabove) = DIV * buildFLUX(w, Iabove)

function buildPFD(grid, wet3D; sinking_speed=1.0)
    w = sinking_speed
    iwet = findall(vec(wet3D))
    DIV = buildDIV(wet3D, iwet, grid)
    Iabove = buildIabove(wet3D, iwet)
    return DIV * buildFLUX(w, Iabove)
end

"""
    buildDIV(wet3D, iwet, grid)

Build the DIV operator such that

    DIV * Φ = 1/dz * (Φ - Φ[below]) ≃ dΦ/dz.
"""
function buildDIV(wet3D, iwet, grid)
    Ibelow = buildIbelow(wet3D, iwet)
    dz = ustrip.(grid.δz_3D[iwet])
    return sparse(Diagonal(dz.^(-1))) * (Ibelow - I) # divergence with positive downwards
end

"""
    buildFLUX(w, Iabove)

Build the `FLUX` operator for a particle sinking speed ``w(r)``

Note `w` is the sinking speed at the top of each grid cell.
(`w` can also be a constant (Float64 type).)

The `FLUX` operator is defined by

    FLUX * x = w * x(above) ≃ Φ.
"""
buildFLUX(w::Vector, Iabove) = sparse(Diagonal(w)) * Iabove
buildFLUX(w::Number, Iabove) = w * Iabove

export buildDIV, buildPFD, buildFLUX

end


text:
using JuMP
using ParameterJuMP
using Polyhedra
mutable struct BalancedComplexPolytope{T, VT <: AbstractVector{T}, D<:Polyhedra.FullDim}
    d::D
    points::Vector{VT}
    factory::Union{Nothing, JuMP.OptimizerFactory}
    model::Union{Nothing, JuMP.Model}
    real_z::Union{Nothing, Vector{ParameterJuMP.ParameterRef}}
    imag_z::Union{Nothing, Vector{ParameterJuMP.ParameterRef}}
    t_0::Union{Nothing, JuMP.VariableRef}
    function BalancedComplexPolytope{T, VT, D}(
        d::Polyhedra.FullDim, points::Polyhedra.PointIt,
        factory::Union{Nothing, JuMP.OptimizerFactory}=nothing) where {T, VT, D}
        new{T, VT, D}(Polyhedra.FullDim_convert(D, d),
                      Polyhedra.lazy_collect(points), factory, nothing, nothing, nothing, nothing)
    end
end
function BalancedComplexPolytope(d::Polyhedra.FullDim, points::Polyhedra.PointIt, args...)
    return BalancedComplexPolytope{Polyhedra.coefficient_type(eltype(points)), eltype(points),
                                   typeof(d)}(d, points, args...)
end

Polyhedra.FullDim(v::BalancedComplexPolytope) = v.d

function _reset_model(p::BalancedComplexPolytope)
    p.model = nothing
    p.real_z = nothing
    p.imag_z = nothing
    p.t_0 = nothing
end

function _build_model(brp::BalancedComplexPolytope)
    n = length(brp.points)
    brp.model = ParameterJuMP.ModelWithParams(brp.factory)
    # Real part of `t`
    α = @variable(brp.model, [1:n])
    # Imaginary part of `t`
    β = @variable(brp.model, [1:n])
    q = @variable(brp.model, [1:n])
    # |t| ≤ q
    for i in 1:n
        @constraint(brp.model,  [q[i], α[i], β[i]] in SecondOrderCone())
    end
    # sum |t| ≤ sum q ≤ 1
    @constraint(brp.model, sum(q) ≤ 1)
    return sum(α[i] .* real.(brp.points[i]) for i in 1:n) - sum(β[i] .* imag.(brp.points[i]) for i in 1:n),
           sum(α[i] .* imag.(brp.points[i]) for i in 1:n) + sum(β[i] .* real.(brp.points[i]) for i in 1:n)
end

function _fix(p::BalancedComplexPolytope, v)
    JuMP.fix.(brp.real_z, real.(v))
    JuMP.fix.(brp.imag_z, imag.(v))
end

function _parametrized_model(brp::BalancedComplexPolytope, v)
    real_hull, imag_hull = _build_model(brp)
    brp.real_z = ParameterJuMP.add_parameters(brp.model, collect(real.(v)))
    brp.imag_z = ParameterJuMP.add_parameters(brp.model, collect(imag.(v)))
    @constraint(brp.model, brp.real_z .== real_hull)
    @constraint(brp.model, brp.imag_z .== imag_hull)
end

function _ratio_model(brp::BalancedComplexPolytope, v)
    real_hull, imag_hull = _build_model(brp)
    brp.t_0 = @variable(brp.model)
    @constraint(brp.model, brp.t_0 .* real.(v) .== real_hull)
    @constraint(brp.model, brp.t_0 .* imag.(v) .== imag_hull)
end


text:
function get_unitcell_honeycomb() :: Unitcell

    # define basis sites
    basis    = Vector{SVector{3, Float64}}(undef, 2)
    basis[1] = SVector{3, Float64}(0.0, 0.0, 0.0)
    basis[2] = SVector{3, Float64}(1.0, 0.0, 0.0)

    # define Bravais vectors
    vectors    = Vector{SVector{3, Float64}}(undef, 3)
    vectors[1] = SVector{3, Float64}(3.0 / 2.0,  sqrt(3.0) / 2.0, 0.0)
    vectors[2] = SVector{3, Float64}(3.0 / 2.0, -sqrt(3.0) / 2.0, 0.0)
    vectors[3] = SVector{3, Float64}(      0.0,              0.0, 0.0)

    # define bonds for basis sites
    bonds     = Vector{Vector{SVector{4, Int64}}}(undef, length(basis))
    bonds1    = Vector{SVector{4, Int64}}(undef, 3)
    bonds1[1] = SVector{4, Int64}( 0, -1,  0,  1)
    bonds1[2] = SVector{4, Int64}(-1,  0,  0,  1)
    bonds1[3] = SVector{4, Int64}( 0,  0,  0,  1)
    bonds[1]  = bonds1
    bonds2    = Vector{SVector{4, Int64}}(undef, 3)
    bonds2[1] = SVector{4, Int64}( 0,  1,  0, -1)
    bonds2[2] = SVector{4, Int64}( 1,  0,  0, -1)
    bonds2[3] = SVector{4, Int64}( 0,  0,  0, -1)
    bonds[2]  = bonds2

    # build unitcell
    uc = Unitcell(basis, vectors, bonds)

    return uc
end


text:
"""
    atom_exchange(flux_dict::Dict{String, Float64}, model::StandardModel)

Return a dictionary mapping the flux of atoms across the boundary of the model
given `flux_dict` (the solution of a constraint based analysis) of reactions in `model`.
"""
function atom_exchange(flux_dict::Dict{String,Float64}, model::StandardModel)
    atom_flux = Dict{String,Float64}()
    for (rxn_id, flux) in flux_dict
        if is_boundary(model.reactions[rxn_id])
            for (met, stoich) in model.reactions[rxn_id].metabolites
                adict = get_atoms(model.metabolites[met])
                for (atom, stoich) in adict
                    atom_flux[atom] = get(atom_flux, atom, 0.0) + flux * stoich
                end
            end
        end
    end
    return atom_flux
end

"""
    get_exchanges(flux_dict::Dict{String, Float64}; top_n=Inf, ignorebound=_constants.default_reaction_bound, verbose=true)

Display the `top_n` producing and consuming exchange fluxes.
If `top_n` is not specified (by an integer), then all are displayed.
Ignores infinite (problem upper/lower bound) fluxes (set with ignorebound).
When `verbose` is false, the output is not printed out.
Return these reactions (id => ) in two dictionaries: `consuming`, `producing`
"""
function exchange_reactions(
    flux_dict::Dict{String,Float64},
    model::StandardModel;
    top_n = Inf,
    ignorebound = _constants.default_reaction_bound,
    verbose = true,
)
    consuming = Dict{String,Float64}()
    producing = Dict{String,Float64}()

    for (k, v) in flux_dict
        if is_boundary(model.reactions[k])
            if v < 0 # consuming
                consuming[k] = v
            elseif v > 0 # producing
                producing[k] = v
            else # no flux
                continue
            end
        end
    end

    if verbose
        # Do consuming
        ks = collect(keys(consuming))
        vs = [consuming[k] for k in ks]
        inds = sortperm(vs)
        n_max = length(ks)
        println("Consuming fluxes: ")
        ii = 0 # counter
        for i in inds
            if vs[i] > -ignorebound
                println(ks[i], " = ", round(vs[i], digits = 6))
                ii += 1
            end
            if ii > top_n
                break
            end
        end
        # Do producing
        ks = collect(keys(producing))
        vs = [producing[k] for k in ks]
        inds = sortperm(vs, rev = true)
        n_max = length(ks)
        println("Producing fluxes: ")
        ii = 0 # counter
        for i in inds
            if vs[i] < ignorebound
                println(ks[i], " = ", round(vs[i], digits = 6))
                ii += 1
            end
            if ii > top_n
                break
            end
        end
    end

    return consuming, producing
end

"""
    metabolite_fluxes(flux_dict::Dict{String, Float64}, model::StandardModel)

Return two dictionaries of metabolite `id`s mapped to reactions that consume or
produce them given the flux distribution supplied in `fluxdict`.
"""
function metabolite_fluxes(flux_dict::Dict{String,Float64}, model::StandardModel)
    S = stoichiometry(model)
    met_flux = Dict{String,Float64}()
    rxnids = reactions(model)
    metids = metabolites(model)

    producing = Dict{String,Dict{String,Float64}}()
    consuming = Dict{String,Dict{String,Float64}}()
    for (row, metid) in enumerate(metids)
        for (col, rxnid) in enumerate(rxnids)
            mf = flux_dict[rxnid] * S[row, col]
            # ignore zero flux
            if mf < -_constants.tolerance # consuming rxn
                if haskey(consuming, metid)
                    consuming[metid][rxnid] = mf
                else
                    consuming[metid] = Dict(rxnid => mf)
                end
            elseif mf > _constants.tolerance
                if haskey(producing, metid)
                    producing[metid][rxnid] = mf
                else
                    producing[metid] = Dict(rxnid => mf)
                end
            end
        end
    end
    return consuming, producing
end

"""
    set_bound(index, optimization_model;
        ub=_constants.default_reaction_rate,
        lb=-_constants.default_reaction_rate)

Helper function to set the bounds of variables.
The JuMP `set_normalized_rhs` function is a little confusing,
so this function simplifies setting constraints. In short, JuMP
uses a normalized right hand side representation of constraints,
which means that lower bounds have their sign flipped. This function
does this for you, so you don't have to remember to do this whenever you
change the constraints.

Just supply the constraint `index` and the JuMP model (`opt_model`) that
will be solved, and the variable's bounds will be set to `ub` and `lb`.
"""
function set_bound(
    vind,
    opt_model;
    ub = _constants.default_reaction_rate,
    lb = -_constants.default_reaction_rate,
)
    set_normalized_rhs(opt_model[:lbs][vind], -lb)
    set_normalized_rhs(opt_model[:ubs][vind], ub)
end

"""
    get_bound_vectors(opt_model)

Returns vectors of the lower and upper bounds of `opt_model` constraints, where
`opt_model` is a JuMP model constructed by e.g.
[`make_optimization_model`](@ref) or [`flux_balance_analysis`](@ref).

"""
function get_bound_vectors(opt_model)
    lbconref = opt_model[:lbs]
    ubconref = opt_model[:ubs]
    lbs = zeros(length(lbconref))
    for i in eachindex(lbs)
        lbs[i] = -normalized_rhs(lbconref[i])
    end
    ubs = [normalized_rhs(ubconref[i]) for i in eachindex(ubconref)]

    return lbs, ubs
end


text:
"""
    em_coefficients(prim, E, B, mr, ld, rL, dt)

Calculate electromagnetic coeffcients in hyperbolic Maxwell's equations
"""
function em_coefficients(
    prim::AbstractMatrix{X},
    E::AbstractVector{Y},
    B::AbstractVector{Z},
    mr,
    lD,
    rL,
    dt,
) where {X<:Real,Y<:Real,Z<:Real}

    if eltype(prim) <: Int
        A = zeros(9, 9)
        b = zeros(9)
    else
        A = similar(prim, 9, 9)
        A .= 0.0
        b = similar(prim, 9)
        b .= 0.0
    end

    A[1, 1] = -1.0 / (2.0 * rL)
    A[2, 2] = -1.0 / (2.0 * rL)
    A[3, 3] = -1.0 / (2.0 * rL)
    A[4, 1] = mr / (2.0 * rL)
    A[5, 2] = mr / (2.0 * rL)
    A[6, 3] = mr / (2.0 * rL)
    A[7, 1] = 1.0 / (dt)
    A[8, 2] = 1.0 / (dt)
    A[9, 3] = 1.0 / (dt)

    A[1, 4] = 1.0 / (dt)
    A[1, 5] = -B[3] / (2.0 * rL)
    A[1, 6] = B[2] / (2.0 * rL)
    A[2, 4] = B[3] / (2.0 * rL)
    A[2, 5] = 1.0 / (dt)
    A[2, 6] = -B[1] / (2.0 * rL)
    A[3, 4] = -B[2] / (2.0 * rL)
    A[3, 5] = B[1] / (2.0 * rL)
    A[3, 6] = 1.0 / (dt)

    A[4, 7] = 1.0 / (dt)
    A[4, 8] = mr * B[3] / (2.0 * rL)
    A[4, 9] = -mr * B[2] / (2.0 * rL)
    A[5, 7] = -mr * B[3] / (2.0 * rL)
    A[5, 8] = 1.0 / (dt)
    A[5, 9] = mr * B[1] / (2.0 * rL)
    A[6, 7] = mr * B[2] / (2.0 * rL)
    A[6, 8] = -mr * B[1] / (2.0 * rL)
    A[6, 9] = 1.0 / (dt)

    A[7, 4] = prim[1, 1] / (2.0 * rL * lD^2)
    A[8, 5] = prim[1, 1] / (2.0 * rL * lD^2)
    A[9, 6] = prim[1, 1] / (2.0 * rL * lD^2)
    A[7, 7] = -(prim[1, 2] * mr) / (2.0 * rL * lD^2)
    A[8, 8] = -(prim[1, 2] * mr) / (2.0 * rL * lD^2)
    A[9, 9] = -(prim[1, 2] * mr) / (2.0 * rL * lD^2)

    b[1] =
        prim[2, 1] / (dt) + E[1] / (2.0 * rL) - B[2] * prim[4, 1] / (2.0 * rL) +
        B[3] * prim[3, 1] / (2.0 * rL)
    b[2] =
        prim[3, 1] / (dt) + E[2] / (2.0 * rL) - B[3] * prim[2, 1] / (2.0 * rL) +
        B[1] * prim[4, 1] / (2.0 * rL)
    b[3] =
        prim[4, 1] / (dt) + E[3] / (2.0 * rL) - B[1] * prim[3, 1] / (2.0 * rL) +
        B[2] * prim[2, 1] / (2.0 * rL)
    b[4] =
        prim[2, 2] / (dt) - mr * E[1] / (2.0 * rL) + mr * B[2] * prim[4, 2] / (2.0 * rL) -
        mr * B[3] * prim[3, 2] / (2.0 * rL)
    b[5] =
        prim[3, 2] / (dt) - mr * E[2] / (2.0 * rL) + mr * B[3] * prim[2, 2] / (2.0 * rL) -
        mr * B[1] * prim[4, 2] / (2.0 * rL)
    b[6] =
        prim[4, 2] / (dt) - mr * E[3] / (2.0 * rL) + mr * B[1] * prim[3, 2] / (2.0 * rL) -
        mr * B[2] * prim[2, 2] / (2.0 * rL)
    b[7] =
        E[1] / (dt) - prim[1, 1] * prim[2, 1] / (2.0 * rL * lD^2) +
        prim[1, 2] * prim[2, 2] * mr / (2.0 * rL * lD^2)
    b[8] =
        E[2] / (dt) - prim[1, 1] * prim[3, 1] / (2.0 * rL * lD^2) +
        prim[1, 2] * prim[3, 2] * mr / (2.0 * rL * lD^2)
    b[9] =
        E[3] / (dt) - prim[1, 1] * prim[4, 1] / (2.0 * rL * lD^2) +
        prim[1, 2] * prim[4, 2] * mr / (2.0 * rL * lD^2)

    return A, b

end


text:
using GridapMHD: hunt

hunt(
  nc=(4,4),
  L=1.0,
  B=(0.,50.,0.),
  debug=false,
  vtk=true,
  nsums=50,
  path=".",
  nruns=1,
  verbose=true,
  title="warmup_hunt",
  solver=:petsc,
  kmap=1,
  petsc_options="-snes_monitor -ksp_error_if_not_converged true -ksp_converged_reason -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps"
)

hunt(
  nc=(4,4),
  np=(1,1),
  backend=:mpi,
  L=1.0,
  B=(0.,50.,0.),
  debug=false,
  vtk=true,
  nsums=50,
  path=".",
  nruns=1,
  verbose=true,
  title="warmup_hunt_mpi",
  solver=:petsc,
  kmap=1,
  petsc_options="-snes_monitor -ksp_error_if_not_converged true -ksp_converged_reason -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps"
)


text:
module FESolversTests

using Test
using Gridap.Arrays
using Gridap.Algebra
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Integration
using Gridap.Fields
using Gridap.FESpaces
using Gridap.CellData

domain =(0,1,0,1,0,1)
partition = (3,3,3)
model = CartesianDiscreteModel(domain,partition)
trian = get_triangulation(model)
degree = 4
quad = CellQuadrature(trian,degree)

order = 2

order = 2
reffe = ReferenceFE(lagrangian,Float64,order)
V = FESpace(model,reffe,dirichlet_tags=[1,10])
U = V
f(x) = x[2]

v = get_cell_shapefuns(V)
u = get_cell_shapefuns_trial(U)

cellmat = integrate(∇(v)⊙∇(u),quad)
cellvec = integrate(v⊙f,quad)
cellids = collect(1:num_cells(trian))

assem = SparseMatrixAssembler(U,V)
matdata = ([cellmat],[cellids],[cellids])
vecdata = ([cellvec],[cellids])
A =  assemble_matrix(assem,matdata)
b =  assemble_vector(assem,vecdata)
x = A \ b
x0 = zeros(length(x))

op = AffineFEOperator(U,V,A,b)
solver = LinearFESolver()
test_fe_solver(solver,op,x0,x)
uh = solve(solver,op)
@test get_free_values(uh) ≈ x
uh = solve(op)
@test get_free_values(uh) ≈ x

solver = NonlinearFESolver()
test_fe_solver(solver,op,x0,x)
uh = solve(solver,op)
@test get_free_values(uh) ≈ x

end # module


text:
module twister
using math
using math.pi
include("src/compiler/index.jl")

# Twister is a pi computations-software
# Copyright 2020: developed by Timo Sarkar

function main
    macro = solve
        pi.compute()
    end
    # println "a"
    @solve = (1/3 + 1/5 + 1/7 + 1/9)
end

export twister.solve
export twister
export module.twister

max_stars_repo_name: timo-cmd/DANCE

max_issues_repo_name: timo-cmd/Twister

max_forks_repo_name: timo-cmd/DANCE


text:
# NOTE: PlanarAzimuthalSymmetry() does not included all possible wavenumbers
function wavenumbers_bisection(ω::T, medium::PhysicalMedium{Dim}, species::Species{T,Dim};
        symmetry::AbstractSymmetry{Dim} = PlanarAzimuthalSymmetry{Dim}(),
        tol::T = 1e-5,
        num_wavenumbers = 3,
        box_k::Vector{Vector{T}} = box_keff(ω, medium, species; tol = tol),
        bisection_mesh_points::Int = Int(round(2 - log(tol))) + 25 * num_wavenumbers,
        fixedparameters::Optim.FixedParameters = NelderMeadparameters(),
        optimoptions::Optim.Options{T} = Optim.Options(
            g_tol = tol^T(3), x_tol=tol^T(2),
            iterations = 800 - Int(round(log(tol)*40))
        ),
        kws...) where {T,Dim}

    ko = real(ω / medium.c)

    disp = dispersion_complex(ω, medium, species, symmetry; kws...)
    # disp = dispersion_complex(ω, medium, species, symmetry; basis_order = basis_order)

    function f(x, y)
        z = disp(x + y*im)
        [real(z),imag(z)]
    end

    x = LinRange(box_k[1][1],box_k[1][2],bisection_mesh_points);
    y = LinRange(box_k[2][1],box_k[2][2], Int(round(bisection_mesh_points/2)));

    mymdbm = MDBM_Problem(f,[Axis(x,"x"),Axis(y,"y")]);

    iteration = 1 #number of refinements (resolution doubling)
    solve!(mymdbm,iteration; interpolationorder = 0);

    x_sol,y_sol=getinterpolatedsolution(mymdbm);
    sols = map(x_sol,y_sol) do x,y
        [x,y]
    end
    sols = reduce_kvecs(sols, ko * sqrt(tol) * T(10))

    # using Plots
    # x_eval,y_eval=getevaluatedpoints(mymdbm)
    # scatter(x_eval[1:7:end],y_eval[1:7:end], label="eval", markersize = 0.0, markerstrokewidth = 0.0)
    # scatter!(x_sol,y_sol, label="sol", markersize = 1.0, markerstrokewidth = 2.0)

    # keep only those with the smallest attenuation
    sols = sort(sols, by = vec -> vec[2])
    sols = sols[1:min(length(sols),10*num_wavenumbers)]

    # Add some guess from asymptotics
    kφ = wavenumber_low_volumefraction(ω, medium, species)

    k_real = max(abs(real(kφ)), real(ko), sqrt(eps(T)))
    k_imag = [imag(kφ), sqrt(eps(T))]

    k_vecs = [[x,y] for x in LinRange(-k_real,k_real,10) for y in k_imag]
    sols = [sols; k_vecs]

    function f_vec(vec::Vector{T})
        abs(disp(vec[1] + vec[2]*im))
    end

    dx = T(0.5) * maximum(abs.(x_sol)) / sqrt(length(sols))
    dy = T(0.5) * maximum(abs.(y_sol)) / sqrt(length(sols))

    # Here we refine the roots
    sols = map(sols) do vec
        inner_optimizer = NelderMead(
            initial_simplex =  MySimplexer(dx,dy),
            parameters  = fixedparameters
        )

        res = optimize(f_vec, vec, inner_optimizer, optimoptions)
        if res.minimum < T(10)*tol || (Optim.converged(res) && res.minimum < T(100)*tol)
            res.minimizer
        else
            [zero(T),-one(T)]
        end
    end
    sols = deleteat!(sols, findall(v-> v == [zero(T),-one(T)], sols))
    sols = reduce_kvecs(sols, T(10) * tol * ko)

    k_effs = [sol[1] + sol[2]*im for sol in sols]
    k_effs = sort(k_effs, by=imag)

    return k_effs
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
using ValueHistories

function msg(args...; newline = true)
    print("   --> ", args...)
    newline && println()
end

n = 100 # increase to match benchmark

msg("Baseline: $n loops that accumulates a Float64")

function f(n)
    tmp = 0.
    for i=1:n
        tmp += i * 3.
    end
    tmp
end

@time f(n)
@time f(n)

#-----------------------------------------------------------

for T in [History, QHistory]

    msg("$(T.name.name): $n loops tracking accumulator as Float64")

    function g(_history,n)
        tmp = 0.
        for i=1:n
            tmp += i * 3.
            push!(_history, i, tmp)
        end
        tmp
    end

    _history = T(Float64)
    @time g(_history,n)

    _history = T(Float64)
    @time g(_history,n)

    msg("$(T.name.name): Converting result into arrays")

    @time x,y = get(_history)
    @time x,y = get(_history)
end

#-----------------------------------------------------------

msg("MVHistory: $n loops tracking accumulator as Float64 and String")

function g(_history,n)
    tmp = 0.
    for i=1:n
        tmp += i * 3.
        push!(_history, :myint, i, tmp)
        push!(_history, :mystr, i, string(tmp))
    end
    tmp
end

_history = MVHistory(History)
@time g(_history,n)

_history = MVHistory(History)
@time g(_history,n)

msg("MVHistory: Converting result into arrays")

@time x,y = get(_history, :mystr)
@time x,y = get(_history, :mystr)

@assert length(y) == n
@assert typeof(y) <: Vector{String}


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

#using Compat

const OBJC_NEW_PROPERTIES = 1

const laszip_BOOL = Cint
const laszip_U8   = Cuchar
const laszip_U16  = UInt16
const laszip_U32  = UInt32
const laszip_U64  = Culonglong
const laszip_I8   = UInt8
const laszip_I16  = Int16
const laszip_I32  = Cint
const laszip_I64  = Clonglong
const laszip_CHAR = UInt8
const laszip_F32  = Cfloat
const laszip_F64  = Cdouble
const laszip_POINTER = Ptr{Cvoid}

mutable struct laszip_geokey
	key_id::UInt16
	tiff_tag_location::UInt16
	count::UInt16
	value_offset::UInt16
end

struct Array_16_UInt8
	d1::UInt8
	d2::UInt8
	d3::UInt8
	d4::UInt8
	d5::UInt8
	d6::UInt8
	d7::UInt8
	d8::UInt8
	d9::UInt8
	d10::UInt8
	d11::UInt8
	d12::UInt8
	d13::UInt8
	d14::UInt8
	d15::UInt8
	d16::UInt8
end

zero(::Type{Array_16_UInt8}) = begin  # /Users/j/.julia/v0.4/Clang/src/wrap_c.jl, line 266:
		Array_16_UInt8(fill(zero(UInt8),16)...)
	end

struct Array_32_UInt8
	d1::UInt8
	d2::UInt8
	d3::UInt8
	d4::UInt8
	d5::UInt8
	d6::UInt8
	d7::UInt8
	d8::UInt8
	d9::UInt8
	d10::UInt8
	d11::UInt8
	d12::UInt8
	d13::UInt8
	d14::UInt8
	d15::UInt8
	d16::UInt8
	d17::UInt8
	d18::UInt8
	d19::UInt8
	d20::UInt8
	d21::UInt8
	d22::UInt8
	d23::UInt8
	d24::UInt8
	d25::UInt8
	d26::UInt8
	d27::UInt8
	d28::UInt8
	d29::UInt8
	d30::UInt8
	d31::UInt8
	d32::UInt8
end

zero(::Type{Array_32_UInt8}) = begin  # /Users/j/.julia/v0.4/Clang/src/wrap_c.jl, line 266:
		Array_32_UInt8(fill(zero(UInt8),32)...)
	end

mutable struct laszip_vlr
	reserved::UInt16
	user_id::Array_16_UInt8
	record_id::UInt16
	record_length_after_header::UInt16
	description::Array_32_UInt8
	data::Ptr{Cuchar}
end

struct Array_8_UInt8
	d1::UInt8
	d2::UInt8
	d3::UInt8
	d4::UInt8
	d5::UInt8
	d6::UInt8
	d7::UInt8
	d8::UInt8
end

zero(::Type{Array_8_UInt8}) = begin  # /Users/j/.julia/v0.4/Clang/src/wrap_c.jl, line 266:
		Array_8_UInt8(fill(zero(UInt8),8)...)
	end

struct Array_5_UInt32
	d1::UInt32
	d2::UInt32
	d3::UInt32
	d4::UInt32
	d5::UInt32
end

zero(::Type{Array_5_UInt32}) = begin  # /Users/j/.julia/v0.4/Clang/src/wrap_c.jl, line 266:
		Array_5_UInt32(fill(zero(UInt32),5)...)
	end

struct Array_15_Culonglong
	d1::Culonglong
	d2::Culonglong
	d3::Culonglong
	d4::Culonglong
	d5::Culonglong
	d6::Culonglong
	d7::Culonglong
	d8::Culonglong
	d9::Culonglong
	d10::Culonglong
	d11::Culonglong
	d12::Culonglong
	d13::Culonglong
	d14::Culonglong
	d15::Culonglong
end

zero(::Type{Array_15_Culonglong}) = begin  # /Users/j/.julia/v0.4/Clang/src/wrap_c.jl, line 266:
		Array_15_Culonglong(fill(zero(Culonglong),15)...)
	end

mutable struct laszip_header
	file_source_ID::UInt16
	global_encoding::UInt16
	project_ID_GUID_data_1::UInt32
	project_ID_GUID_data_2::UInt16
	project_ID_GUID_data_3::UInt16
	project_ID_GUID_data_4::Array_8_UInt8
	version_major::Cuchar
	version_minor::Cuchar
	system_identifier::Array_32_UInt8
	generating_software::Array_32_UInt8
	file_creation_day::UInt16
	file_creation_year::UInt16
	header_size::UInt16
	offset_to_point_data::UInt32
	number_of_variable_length_records::UInt32
	point_data_format::Cuchar
	point_data_record_length::UInt16
	number_of_point_records::UInt32
	number_of_points_by_return::Array_5_UInt32
	x_scale_factor::Cdouble
	y_scale_factor::Cdouble
	z_scale_factor::Cdouble
	x_offset::Cdouble
	y_offset::Cdouble
	z_offset::Cdouble
	max_x::Cdouble
	min_x::Cdouble
	max_y::Cdouble
	min_y::Cdouble
	max_z::Cdouble
	min_z::Cdouble
# LAS 1.3 and higher only
	start_of_waveform_data_packet_record::Culonglong

# LAS 1.4 and higher only
	start_of_first_extended_variable_length_record::Culonglong
	number_of_extended_variable_length_records::UInt32
	extended_number_of_point_records::Culonglong
	extended_number_of_points_by_return::Array_15_Culonglong

# optional
	user_data_in_header_size::UInt32
	user_data_in_header::Ptr{Cuchar}

# optional VLRs
	vlrs::Ptr{laszip_vlr}

# optional
	user_data_after_header_size::UInt32
	user_data_after_header::Ptr{Cuchar}
end

struct Array_4_UInt16
	d1::UInt16
	d2::UInt16
	d3::UInt16
	d4::UInt16
end

zero(::Type{Array_4_UInt16}) = begin  # /Users/j/.julia/v0.4/Clang/src/wrap_c.jl, line 266:
		Array_4_UInt16(fill(zero(UInt16),4)...)
	end

struct Array_7_Cuchar
	d1::Cuchar
	d2::Cuchar
	d3::Cuchar
	d4::Cuchar
	d5::Cuchar
	d6::Cuchar
	d7::Cuchar
end

zero(::Type{Array_7_Cuchar}) = begin  # /Users/j/.julia/v0.4/Clang/src/wrap_c.jl, line 266:
		Array_7_Cuchar(fill(zero(Cuchar),7)...)
	end

struct Array_29_Cuchar
	d1::Cuchar
	d2::Cuchar
	d3::Cuchar
	d4::Cuchar
	d5::Cuchar
	d6::Cuchar
	d7::Cuchar
	d8::Cuchar
	d9::Cuchar
	d10::Cuchar
	d11::Cuchar
	d12::Cuchar
	d13::Cuchar
	d14::Cuchar
	d15::Cuchar
	d16::Cuchar
	d17::Cuchar
	d18::Cuchar
	d19::Cuchar
	d20::Cuchar
	d21::Cuchar
	d22::Cuchar
	d23::Cuchar
	d24::Cuchar
	d25::Cuchar
	d26::Cuchar
	d27::Cuchar
	d28::Cuchar
	d29::Cuchar
end

zero(::Type{Array_29_Cuchar}) = begin  # /Users/j/.julia/v0.4/Clang/src/wrap_c.jl, line 266:
		Array_29_Cuchar(fill(zero(Cuchar),29)...)
	end

mutable struct laszip_point
	X::Cint
	Y::Cint
	Z::Cint
	intensity::UInt16
	return_number::Cuchar
	number_of_returns::Cuchar
	scan_direction_flag::Cuchar
	edge_of_flight_line::Cuchar
	classification::Cuchar
	synthetic_flag::Cuchar
	keypoint_flag::Cuchar
	withheld_flag::Cuchar
	scan_angle_rank::UInt8
	user_data::Cuchar
	point_source_ID::UInt16
# LAS 1.4 only
	extended_scan_angle::Int16
	extended_point_type::Cuchar
	extended_scanner_channel::Cuchar
	extended_classification_flags::Cuchar
	extended_classification::Cuchar
	extended_return_number::Cuchar
	extended_number_of_returns::Cuchar
# for 8 byte alignment of the GPS time
	dummy::Array_7_Cuchar
	gps_time::Cdouble
	rgb::Array_4_UInt16
	wave_packet::Array_29_Cuchar
	num_extra_bytes::Cint
	extra_bytes::Ptr{Cuchar}
end

#=
mutable struct laszip_dll
	header::laszip_header
	p_count::Clonglong
	npoints::Clonglong
	point::laszip_point
	point_items::Ptr{Ptr{UInt8}}
	FILE* file;
	ByteStreamIn* streamin;
	LASreadPoint* reader;
	ByteStreamOut* streamout;
	LASwritePoint* writer;
	CHAR error[1024];
	CHAR warning[1024];
end
=#

max_stars_repo_name: joa-quim/Laszip

max_issues_repo_name: joa-quim/Laszip

max_forks_repo_name: joa-quim/Laszip


text:
@testset "K" begin
    @testset "ktotal" begin
        try
            furnsh(path(CORE, :pck))
            @test  ktotal("ALL") == 1
        finally
            kclear()
        end
        @test ktotal("ALL") == 0
    end
    @testset "kdata" begin
        try
            furnsh(path(CORE, :pck))
            file, ftype, source, handle = kdata(1, "ALL")
            @test ftype == "TEXT"
        finally
            kclear()
        end
        @test kdata(1, "ALL") === nothing
    end
    @testset "kinfo" begin
        try
            furnsh(path(CORE, :pck))
            filetype, source, handle = kinfo(path(CORE, :pck))
            @test filetype == "TEXT"
        finally
            kclear()
        end
    end
    @testset "kplfrm" begin
        kclear()
        try
            furnsh(path(CASSINI, :fk))
            cell = kplfrm(-1)
            @test length(cell) == 57
        finally
            kclear()
        end
    end
    @testset "kxtrct" begin
        terms = ["FROM", "TO", "BEGINNING", "ENDING"]
        str = "FROM 1 October 1984 12:00:00 TO 1 January 1987"

        exp_string = " TO 1 January 1987"
        exp_substr = "1 October 1984 12:00:00"
        act_string, act_substr = kxtrct("FROM", terms, str)
        @test act_string == exp_string
        @test act_substr == exp_substr

        exp_string = "FROM 1 October 1984 12:00:00"
        exp_substr = "1 January 1987"
        act_string, act_substr = kxtrct("TO", terms, str)
        @test act_string == exp_string
        @test act_substr == exp_substr

        terms = ["ADDRESS:", "PHONE:", "NAME:"]
        str = "ADDRESS: 4800 OAK GROVE DRIVE PHONE: 354-4321"

        exp_string = " PHONE: 354-4321"
        exp_substr = "4800 OAK GROVE DRIVE"
        act_string, act_substr = kxtrct("ADDRESS:", terms, str)
        @test act_string == exp_string
        @test act_substr == exp_substr

        exp_string = "ADDRESS: 4800 OAK GROVE DRIVE"
        exp_substr = "354-4321"
        act_string, act_substr = kxtrct("PHONE:", terms, str)
        @test act_string == exp_string
        @test act_substr == exp_substr

        @test kxtrct("NAME:", terms, str) === nothing
    end
end


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

using BAT
using Test

using LinearAlgebra
using ValueShapes, ArraysOfArrays, Distributions, ForwardDiff

import Cuba

@testset "transformed_density" begin
    @testset "distribution transforms" begin
        function test_uv_transformed(target_type::Type{<:Distribution}, source_dist::Distribution)
            trafo = BAT.DistributionTransform(target_type, source_dist)
            @testset "$(typeof(trafo.source_dist))) to $(typeof(trafo.target_dist))" begin
                if target_type == Uniform
                    @test trafo.target_dist isa BAT.StandardUvUniform
                elseif target_type == Normal
                    @test trafo.target_dist isa BAT.StandardUvNormal
                end

                target_dist = trafo.target_dist
                source_dist = trafo.source_dist

                source_x = mean(target_dist) + std(target_dist) / 2
                @test @inferred(trafo(source_x)) isa Real
                target_x = trafo(source_x)
                @test @inferred(trafo(fill(source_x))) == fill(target_x)

                source_X = rand(source_dist, 10^5)
                target_X = @inferred broadcast(trafo, source_X)
                @test isapprox(@inferred(broadcast(inv(trafo), (target_X))), source_X, atol = 10^-8)
                @test isapprox(mean(target_X), mean(target_dist), atol = 0.05)
                @test isapprox(var(target_X), var(target_dist), atol = 0.1)

                @test @inferred(trafo(convert(AbstractDensity, source_dist))) isa BAT.TransformedDensity
                density = trafo(convert(AbstractDensity, source_dist))

                @test isapprox(@inferred(inv(density.trafo)(target_x)), source_x, atol = 10^-5)
                @test isapprox(@inferred(inv(density.trafo)(fill(target_x))), fill(source_x), atol = 10^-5)

                @test minimum(target_dist) <= @inferred(bat_initval(density, InitFromTarget())).result <= maximum(target_dist)
                @test all(minimum(target_dist) .<= @inferred(bat_initval(density, 100, InitFromTarget())).result .<= maximum(target_dist))

                fix_nni(x::T) where {T<:Real} = x <= BAT.near_neg_inf(T) ? T(-Inf) : x

                tX = if target_type == Uniform
                    [0+eps(Float64), 0.25, 0.75, 1-eps(Float64)]
                elseif target_type == Normal
                    [-Inf, -2.1, -1.2, 0.0, 1.2, 2.1, Inf]
                else
                    @assert false
                end

                @test isapprox(fix_nni.(logdensityof(density).(tX)), logpdf.(target_dist, tX), atol = 1e-10)
                @test @inferred(logdensityof(density)(target_x)) isa Real
                @test @inferred(logdensityof(unshaped(density))(unshaped(target_x))) isa Real
                !any(isnan, @inferred(broadcast(ForwardDiff.derivative, logdensityof(density), tX)))

                tX_finite = tX[findall(isfinite, fix_nni.(logdensityof(density).(tX)))]
                @test isapprox(@inferred(broadcast(ForwardDiff.derivative, logdensityof(density), tX_finite)), broadcast(ForwardDiff.derivative, x -> logpdf(target_dist, x), tX_finite), atol = 10^-7)

                @test minimum(target_dist) <= bat_findmode(density, MaxDensityLBFGS(trafo = NoDensityTransform())).result <= maximum(target_dist)

                if trafo.target_dist isa Union{BAT.StandardUvUniform,BAT.StandardMvUniform}
                    @test isapprox(bat_integrate(density, VEGASIntegration(trafo = NoDensityTransform())).result, 1, rtol = 10^-7)
                end
            end
        end

        test_uv_transformed(Uniform, Weibull())
        test_uv_transformed(Uniform, Normal(2, 4))
        test_uv_transformed(Uniform, Uniform(-2, 3))

        test_uv_transformed(Normal, Weibull())
        test_uv_transformed(Normal, Normal(2, 4))
        test_uv_transformed(Normal, Uniform(-2, 3))

        src_d = NamedTupleDist(a = Exponential(), b = [4.2, 3.3], c = Weibull(), d = [Normal(1, 3), Normal(3, 2)], e = Uniform(-2, 3), f = MvNormal([0.3, -2.9], [1.7 0.5; 0.5 2.3]))
        trafo = @inferred(BAT.DistributionTransform(Normal, src_d))
        density = @inferred(trafo(convert(AbstractDensity, trafo.source_dist)))
        @test isfinite(@inferred logdensityof(density)(@inferred(bat_initval(density)).result))
        @test isapprox(cov(@inferred(bat_initval(density, 10^4)).result), I(totalndof(density)), rtol = 0.1)

        samples_is = bat_sample(density, MCMCSampling(mcalg = HamiltonianMC(), trafo = NoDensityTransform(), nsteps = 10^4)).result
        @test isapprox(cov(samples_is), I(totalndof(density)), rtol = 0.1)
        samples_os = inv(trafo).(samples_is)
        @test all(isfinite, logpdf.(Ref(src_d), samples_os.v))
        @test isapprox(cov(unshaped.(samples_os)), cov(unshaped(src_d)), rtol = 0.1)
        @test isapprox(mean(unshaped.(samples_os)), mean(rand(unshaped(src_d), 10^5), dims = 2), rtol = 0.1)

        primary_dist = NamedTupleDist(a = Normal(), b = Weibull(), c = 5)
        f_secondary = x -> NamedTupleDist(y = Normal(x.a, x.b), z = MvNormal([1.3 0.5; 0.5 2.2]))
        prior = HierarchicalDistribution(f_secondary, primary_dist)
        likelihood = varshape(prior)(MvNormal(Diagonal(fill(1.0, totalndof(varshape(prior))))))
        density = PosteriorDensity(likelihood, prior)
        hmc_samples = bat_sample(density, MCMCSampling(mcalg = HamiltonianMC(), trafo = PriorToGaussian(), nsteps = 10^4)).result
        is_samples = bat_sample(density, PriorImportanceSampler(nsamples = 10^4)).result
        @test isapprox(mean(unshaped.(hmc_samples)), mean(unshaped.(is_samples)), rtol = 0.1)
        @test isapprox(cov(unshaped.(hmc_samples)), cov(unshaped.(is_samples)), rtol = 0.2)
    end
end