text:
# AtmosRefStatePerturbations
#
# Computes perturbations from the reference state and outputs them
# on the specified interpolated grid.

import CUDA
using ..Atmos
using ..Mesh.Topologies
using ..Mesh.Grids
using ..Thermodynamics

"""
    setup_atmos_refstate_perturbations(
        ::ClimateMachineConfigType,
        interval::String,
        out_prefix::String;
        writer = NetCDFWriter(),
        interpol = nothing,
    )

Create and return a `DiagnosticsGroup` containing the
"AtmosRefStatePerturbations" diagnostics for both LES and GCM
configurations. All the diagnostics in the group will run at the
specified `interval`, be interpolated to the specified boundaries
and resolution, and written to files prefixed by `out_prefix`
using `writer`.
"""
function setup_atmos_refstate_perturbations(
    ::ClimateMachineConfigType,
    interval::String,
    out_prefix::String;
    writer = NetCDFWriter(),
    interpol = nothing,
)
    # TODO: remove this
    @assert !isnothing(interpol)

    return DiagnosticsGroup(
        "AtmosRefStatePerturbations",
        Diagnostics.atmos_refstate_perturbations_init,
        Diagnostics.atmos_refstate_perturbations_fini,
        Diagnostics.atmos_refstate_perturbations_collect,
        interval,
        out_prefix,
        writer,
        interpol,
    )
end

function vars_atmos_refstate_perturbations(m::AtmosModel, FT)
    @vars begin
        ref_state::vars_atmos_refstate_perturbations(m.ref_state, FT)
    end
end
vars_atmos_refstate_perturbations(::ReferenceState, FT) = @vars()
function vars_atmos_refstate_perturbations(rs::HydrostaticState, FT)
    @vars begin
        rho::FT
        pres::FT
        temp::FT
        et::FT
        qt::FT
    end
end
num_atmos_refstate_perturbation_vars(m, FT) =
    varsize(vars_atmos_refstate_perturbations(m, FT))
atmos_refstate_perturbation_vars(m, array) =
    Vars{vars_atmos_refstate_perturbations(m, eltype(array))}(array)

function atmos_refstate_perturbations!(
    atmos::AtmosModel,
    state,
    aux,
    thermo,
    vars,
)
    atmos_refstate_perturbations!(
        atmos.ref_state,
        atmos,
        state,
        aux,
        thermo,
        vars,
    )
    return nothing
end
function atmos_refstate_perturbations!(
    ::ReferenceState,
    ::AtmosModel,
    state,
    aux,
    thermo,
    vars,
)
    return nothing
end
function atmos_refstate_perturbations!(
    rs::HydrostaticState,
    atmos::AtmosModel,
    state,
    aux,
    thermo,
    vars,
)
    vars.ref_state.rho = state.ρ - aux.ref_state.ρ
    vars.ref_state.pres = thermo.pres - aux.ref_state.p
    vars.ref_state.temp = thermo.temp - aux.ref_state.T
    vars.ref_state.et =
        (state.ρe / state.ρ) - (aux.ref_state.ρe / aux.ref_state.ρ)
    # FIXME properly
    if atmos.moisture isa EquilMoist
        vars.ref_state.qt =
            (thermo.moisture.ρq_tot / state.ρ) -
            (aux.ref_state.ρq_tot / aux.ref_state.ρ)
    else
        vars.ref_state.qt = NaN
    end

    return nothing
end

function atmos_refstate_perturbations_init(dgngrp::DiagnosticsGroup, currtime)
    FT = eltype(Settings.Q)
    atmos = Settings.dg.balance_law
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)

    if mpirank == 0
        # get dimensions for the interpolated grid
        dims = dimensions(dgngrp.interpol)

        # adjust the level dimension for `planet_radius` on 'CubedSphere's
        if dgngrp.interpol isa InterpolationCubedSphere
            level_val = dims["level"]
            dims["level"] = (
                level_val[1] .- FT(planet_radius(Settings.param_set)),
                level_val[2],
            )
        end

        # set up the variables we're going to be writing
        vars = OrderedDict()
        varnames = map(
            s -> startswith(s, "ref_state.") ? s[11:end] : s,
            flattenednames(vars_atmos_refstate_perturbations(atmos, FT)),
        )
        for varname in varnames
            var = Variables[varname]
            vars[varname] = (tuple(collect(keys(dims))...), FT, var.attrib)
        end

        # create the output file
        dprefix = @sprintf(
            "%s_%s_%s_rank%04d",
            dgngrp.out_prefix,
            dgngrp.name,
            Settings.starttime,
            mpirank,
        )
        dfilename = joinpath(Settings.output_dir, dprefix)
        init_data(dgngrp.writer, dfilename, dims, vars)
    end

    return nothing
end

"""
    atmos_refstate_perturbations_collect(dgngrp, currtime)

Perform a global grid traversal to compute various diagnostics.
"""
function atmos_refstate_perturbations_collect(
    dgngrp::DiagnosticsGroup,
    currtime,
)
    dg = Settings.dg
    atmos = dg.balance_law

    # FIXME properly
    if !isa(atmos.ref_state, HydrostaticState)
        @warn """
            Diagnostics $(dgngrp.name): has useful output only for `HydrostaticState`
            """
    end

    Q = Settings.Q
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)
    grid = dg.grid
    topology = grid.topology
    N = polynomialorder(grid)
    Nq = N + 1
    Nqk = dimensionality(grid) == 2 ? 1 : Nq
    npoints = Nq * Nq * Nqk
    nrealelem = length(topology.realelems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    # get needed arrays onto the CPU
    if array_device(Q) isa CPU
        ArrayType = Array
        state_data = Q.realdata
        aux_data = dg.state_auxiliary.realdata
    else
        ArrayType = CUDA.CuArray
        state_data = Array(Q.realdata)
        aux_data = Array(dg.state_auxiliary.realdata)
    end
    FT = eltype(state_data)

    # Compute thermo variables
    thermo_array = Array{FT}(undef, npoints, num_thermo(atmos, FT), nrealelem)
    @visitQ nhorzelem nvertelem Nqk Nq begin
        state = extract_state(dg, state_data, ijk, e, Prognostic())
        aux = extract_state(dg, aux_data, ijk, e, Auxiliary())

        thermo = thermo_vars(atmos, view(thermo_array, ijk, :, e))
        compute_thermo!(atmos, state, aux, thermo)
    end

    # Interpolate the state and thermo variables.
    interpol = dgngrp.interpol
    istate =
        ArrayType{FT}(undef, interpol.Npl, number_states(atmos, Prognostic()))
    interpolate_local!(interpol, Q.realdata, istate)

    if interpol isa InterpolationCubedSphere
        # TODO: get indices here without hard-coding them
        _ρu, _ρv, _ρw = 2, 3, 4
        project_cubed_sphere!(interpol, istate, (_ρu, _ρv, _ρw))
    end

    iaux = ArrayType{FT}(undef, interpol.Npl, number_states(atmos, Auxiliary()))
    interpolate_local!(interpol, dg.state_auxiliary.realdata, iaux)

    ithermo = ArrayType{FT}(undef, interpol.Npl, num_thermo(atmos, FT))
    interpolate_local!(interpol, ArrayType(thermo_array), ithermo)

    # FIXME: accumulating to rank 0 is not scalable
    all_state_data = accumulate_interpolated_data(mpicomm, interpol, istate)
    all_aux_data = accumulate_interpolated_data(mpicomm, interpol, iaux)
    all_thermo_data = accumulate_interpolated_data(mpicomm, interpol, ithermo)

    if mpirank == 0
        # get dimensions for the interpolated grid
        dims = dimensions(dgngrp.interpol)

        # set up the array for the diagnostic variables based on the interpolated grid
        nlong = length(dims["long"][1])
        nlat = length(dims["lat"][1])
        nlevel = length(dims["level"][1])

        perturbations_array = Array{FT}(
            undef,
            nlong,
            nlat,
            nlevel,
            num_atmos_refstate_perturbation_vars(atmos, FT),
        )

        @visitI nlong nlat nlevel begin
            statei = Vars{vars_state(atmos, Prognostic(), FT)}(view(
                all_state_data,
                lo,
                la,
                le,
                :,
            ))
            auxi = Vars{vars_state(atmos, Auxiliary(), FT)}(view(
                all_aux_data,
                lo,
                la,
                le,
                :,
            ))
            thermoi = thermo_vars(atmos, view(all_thermo_data, lo, la, le, :))

            perturbations = atmos_refstate_perturbation_vars(
                atmos,
                view(perturbations_array, lo, la, le, :),
            )
            atmos_refstate_perturbations!(
                atmos,
                statei,
                auxi,
                thermoi,
                perturbations,
            )
        end

        varvals = OrderedDict()
        varnames = map(
            s -> startswith(s, "ref_state.") ? s[11:end] : s,
            flattenednames(vars_atmos_refstate_perturbations(atmos, FT)),
        )
        for (vari, varname) in enumerate(varnames)
            varvals[varname] = perturbations_array[:, :, :, vari]
        end

        # write output
        append_data(dgngrp.writer, varvals, currtime)
    end

    MPI.Barrier(mpicomm)
    return nothing
end # function collect

function atmos_refstate_perturbations_fini(dgngrp::DiagnosticsGroup, currtime) end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using LinearAlgebra

function so3ToVec(so3mat)
    if( size(so3mat)!=(3,3) )
        println("!Warning! : Input matrix should be 3x3.")
        return
    end
    if( !((so3mat[1,1]==so3mat[2,2]==so3mat[3,3]==0)) )
        println("!Warning! : Diagonals should be zero.")
        return
    end
    if( so3mat[2,1]!=-so3mat[1,2] ||
        so3mat[3,1]!=-so3mat[1,3] ||
        so3mat[3,2]!=-so3mat[2,3] )
        println("!Warning! : Entered matrix isn't a skew symmetric matrix.")
        return
    end
    ret_vec = [so3mat[3,2];
               so3mat[1,3];
               so3mat[2,1] ]
        return ret_vec
end

function vecToSo3(omg)
    if( size(omg)!=(3,) )
        println("!Warning! : Input matrix should be 3x1.")
        return
    end

    ret_ss_matrix = [0      -omg[3]   omg[2];
                     omg[3]    0     -omg[1];
                    -omg[2]  omg[1]      0];
    return ret_ss_matrix
end

using LinearAlgebra

function so3ToVec(so3mat)
    if( size(so3mat)!=(3,3) )
        println("!Warning! : Input matrix should be 3x3.")
        return
    end
    if( !((so3mat[1,1]==so3mat[2,2]==so3mat[3,3]==0)) )
        println("!Warning! : Diagonals should be zero.")
        return
    end
    if( so3mat[2,1]!=-so3mat[1,2] ||
        so3mat[3,1]!=-so3mat[1,3] ||
        so3mat[3,2]!=-so3mat[2,3] )
        println("!Warning! : Entered matrix isn't a skew symmetric matrix.")
        return
    end
    ret_vec = [so3mat[3,2];
               so3mat[1,3];
               so3mat[2,1] ]
        return ret_vec
end

function VecToso3(omg)
    if( size(omg)!=(3,) )
        println("!Warning! : Input matrix should be 3x1.")
        return
    end

    ret_ss_matrix = [0      -omg[3]   omg[2];
                     omg[3]    0     -omg[1];
                    -omg[2]  omg[1]      0];
    return ret_ss_matrix
end
#****************************************************************#
# VeysiADN 17 Apr 2021
# Advanced Robotics Homework 4 Codes
#****************************************************************#
using LinearAlgebra
const lAlgebra = LinearAlgebra
#****************************************************************#
# Returns normalized version of given matrix/array
function Normalize(N)
    return N/lAlgebra.norm(N)
end
#****************************************************************#
# Returns 4x4 twist matrix ‘se3mat’ (an element of Lie Algebra se3)
# from given twist vector (= 6x1 vector) ‘xi
function VecTose3( xi )
    w_skew = [ 0   -xi[3]   xi[2];
              xi[3]   0    -xi[1];
             -xi[2]  xi[1]    0 ]
    return vcat(hcat(w_skew,xi[4:6]), zeros(1, 4))
end
#****************************************************************#
# Returns twist vector ‘xi’ from given twist matrix ‘se3mat’
function se3ToVec( se3mat )
    omg = [se3mat[3, 2], se3mat[1, 3],se3mat[2, 1]]
    v = se3mat[1:3, 4]
    return vcat(omg,v)
end
#****************************************************************#
# Calculates and return rotation matrix ‘R’ by
# evaluating matrix exponential of 3x3 skew symmetric matrix ‘so3mat’

function MatrixExp3(so3mat)
        omg = [so3mat[3, 2], so3mat[1, 3],so3mat[2, 1]]
        theta = lAlgebra.norm(omg)
        omg_n = so3mat / theta
        return lAlgebra.I + sin(theta) * omg_n + (1 - cos(theta)) * omg_n * omg_n
end

#****************************************************************#
# Calculates and return transformation matrix ‘g’ by
# evaluating matrix exponential of 4x4 twist matrix ‘se3mat’
function MatrixExp6( se3mat )
    omg = [se3mat[3, 2], se3mat[1, 3],se3mat[2, 1]]
    theta = lAlgebra.norm(omg)
    omg_n = se3mat[1:3, 1:3] / theta
    temp  = hcat(MatrixExp3(se3mat[1:3, 1:3]),
                         (lAlgebra.I * theta +
                          (1 - cos(theta)) * omg_n +
                          (theta - sin(theta)) * omg_n * omg_n) *
                         se3mat[1:3, 4] / theta)
    return    vcat(temp,[0 0 0 1])
end

#****************************************************************#
# VeysiADN 12 May 2021
# Advanced Robotics Homework 5 Codes
#****************************************************************#
using LinearAlgebra
const lAlgebra = LinearAlgebra
#****************************************************************#
# Computes exponential coordinate ‘so3mat’
# (skew symmetric matrix) of the rotation matrix R

function MatrixLog3(R)
    acosinput = (lAlgebra.tr(R) - 1) / 2
    if acosinput >= 1
        return zeros(3, 3)
    else
        theta = acos(acosinput)
        return theta / 2 / sin(theta) * (R - R')
    end
end
#****************************************************************#
# Computes exponential coordinate ‘se3mat’
# (twist matrix) of the rigid body transformation matrix T
function MatrixLog6(T)
    if( size(T)!=(4,4) )
        println("!Warning! : Transformation matrix should be 4x4.")
        return
    end
    R = T[1:3,1:3]
    t = T[1:3,4]
    omg_h = MatrixLog3(R)
    w = [omg_h[3,2];omg_h[1,3];omg_h[2,1]]
    if omg_h == zeros(3, 3)
        return vcat(hcat(zeros(3, 3), T[1:3, 4]), [0 0 0 0])
    else
        theta = acos((lAlgebra.tr(R) - 1) / 2)
        return vcat(hcat(omg_h,
                         (I - omg_h / 2 +
                          (1 / theta - 1 / tan(theta / 2) / 2) *
                          omg_h * omg_h / theta) * T[1:3, 4]),
                    [0 0 0 0])
    end
end
#****************************************************************#

#****************************************************************#
# VeysiADN 12 May 2021
# Advanced Robotics Homework 5 Codes
#****************************************************************#
using LinearAlgebra
#****************************************************************#
# Returns Adjoint matrix with given rigid body
# transformation matrix g
function Adjoint(g)

    R=g[1:3,1:3];
    p=g[1:3,4];
    p_hat = [ 0       -p[3]     p[2];
             p[3]       0      -p[1];
            -p[2]      p[1]       0 ]

    return [R zeros(3,3);p_hat*R R]
end
#****************************************************************#
# VeysiADN 4 Apr 2021
# Advanced Robotics Homework 3 Codes
#****************************************************************#
using LinearAlgebra
#****************************************************************#
# Returns 3 x 3 rotation matrix representing
# rotation about unit vector w by angle theta.
# theta should be in radian
function Rodrigues(w,theta)
    if( length(w) != 3)
        println("!Warning! : W matrix should have 3 elements.")
        return
    end
    if(size(w)==(1,3))
        w=w';
    end
    w_hat = [0      -w[3]    w[2];
             w[3]     0     -w[1];
            -w[2]    w[1]      0];
    rotation_matrix = cos(theta)*Matrix(1I,3,3) + sin(theta)*w_hat+(1-cos(theta))*w*w';
    return rotation_matrix;
end
#****************************************************************#
# Builds rigid transformation matrix T from given
#rotation matrix R and translation t
function RpToTrans(R,t)
    if( size(R)!=(3,3) )
        println("!Warning! : Rotation matrix should be 3x3.")
        return
    end
    if( length(t) != 3 )
        println("!Warning! : Translation matrix should have 3 elements.")
        return
    end
    if(size(t)==(1,3))
        t=t';
    end
    R_B_T = hcat(R,t) ;
    R_B_T = vcat(R_B_T,[0 0 0 1]) ;
    return R_B_T
end
#****************************************************************#
#Extracts rotation matrix R and translation t from given /
# rigid transformation matrix T
function TransToRp(T)
    if( size(T)!=(4,4) )
        println("!Warning! : Transformation matrix should be 4x4.")
        return
    end
    R = T[1:3,1:3]
    t = T[1:3,4]
    return R,t;
end
#****************************************************************#
# Inverts rigid transformation matrix T
function TransInv(T)
    if( size(T)!=(4,4) )
        println("!Warning! : W matrix should be 4x4.")
        return
    end
    if( det(T)==0 )
        println("!Warning! : Inverse doesn't exist (Determinant is equal to zero).")
        return
    end
    R,t = TransToRp(T);
    r_inv =hcat(R',-(R'*t))
    r_inv = vcat(r_inv,[0 0 0 1]);
    return r_inv;
end
#****************************************************************#

ext: jl

max_stars_repo_name: veysiadn/AdvancedRoboticsJulia

max_issues_repo_name: veysiadn/AdvancedRoboticsJulia

max_forks_repo_name: veysiadn/AdvancedRoboticsJulia


text:
"""
    trackpath(H::AbstractHomotopy, startvalue, algorithm::PredictorCorrector, start, finish; kwargs...)

Tracks the path ``z(t)`` which is implicitly defined via ``H(x,t)`` from
`z(start)=H(startvalue, start)` to `z(finish)= H(result, finish)` with the given
PredictorCorrector `algorithm`.

We reformulate the path ``z(t)`` as a path z(γ(s)) where γ: [1, 0] -> line from start to finish

## Optional arguments
* `maxiterations=10000`
* `tolerance=1e-4`: The tolerance used for the correction step
* `refinement_tolerance=1e-8`: The result will be refined to the given accuracy
* `initial_steplength=0.01`: The initial_steplength (w.r.t s)
* `correction_step_maxiterations=3`
* `successfull_steps_until_steplength_increase=3`
* `steplength_increase_factor=2.0`
* `steplength_decrease_factor=0.5`
* `record_trace=false`
* `record_steps=false`
* `debug=false`: Print *(quite a lot)* different values to debug behaviour.
"""
function trackpath(
    H::AbstractHomotopy{T},
    startvalue::Vector{T},
    algorithm::AbstractPredictorCorrectorAlgorithm,
    start::S,
    finish::S;kwargs...
) where {S<:Number,T<:Number}
    J_H = differentiate(H)
    ∂H∂t = ∂t(H)
    trackpath(H, J_H, ∂H∂t, startvalue, algorithm, start, finish; kwargs...)
end

function trackpath(
    H::AbstractHomotopy{T},
    J_H,
    ∂H∂t,
    startvalue::Vector{T},
    algorithm::AbstractPredictorCorrectorAlgorithm,
    start::S,
    finish::S;
    maxiterations=10000,
    tolerance=1e-4,
    refinement_tolerance=1e-12,
    initial_steplength=0.05,
    successfull_steps_until_steplength_increase=3,
    steplength_increase_factor=2.0,
    steplength_decrease_factor=0.5,
    record_trace=false,
    record_steps=false,
    correction_step_maxiterations=3,
    debug=false
) where {S<:Number,T<:Number}
    # We want to track a path z(t) from z(start) to z(finish).
    # We reformulate it as a problem z(γ(s)) where γ: [1, 0] -> line from start to finish
    γlength = finish - start
    γ(s) = start + (1 - s) * γlength

    # some setup before we start the main loop
    x = copy(startvalue)
    u = similar(x)

    steps::Vector{S} = [];
    trace::Vector{Vector{T}} = [];

    s = 1.0
    steplength = initial_steplength
    sucessive_successes = 0
    k = 0

    while s > 0 && k < maxiterations
        if record_steps
            push!(steps, s)
        end
        if record_trace
            push!(trace, x)
        end

        Δs= min(steplength, s)

        u .= predict(algorithm, H, J_H, ∂H∂t, x, γ(s), Δs * γlength)

        if debug
            println("")
            println("t: $t Δt: $Δt, iteration: $k")
            println("x: $x, norm: $(norm(evaluate(H,x,t)))")
            println("predicted: $x, norm: $(norm(evaluate(H,x,t)))")
        end

        converged = correct!(u, algorithm, H, J_H, x, γ(s - Δs), tolerance, correction_step_maxiterations)

        if converged
            x .= u
            s -= Δs

            sucessive_successes += 1
            if sucessive_successes == successfull_steps_until_steplength_increase
                steplength *= steplength_increase_factor
                sucessive_successes = 0
            end
        else
            if debug
                println("correction failed!")
            end
            sucessive_successes = 0
            steplength *= steplength_decrease_factor
        end

        k += 1
    end

    if s ≈ 0 && norm(evaluate(H, x, γ(0.0))) < tolerance
        refinement_converged = correct!(u, algorithm, H, J_H, x, γ(0.0), refinement_tolerance, 10)
        if refinement_converged
            x .= u
        else

        end
    end

    if k >= maxiterations
        retcode = :MaxIterations
    elseif norm(evaluate(H, x, γ(0.0))) > refinement_tolerance
        retcode = :Diverged
    else
        retcode = :Success
    end

    PathResult(x, retcode, startvalue, k, γ(s), steps, trace)
end


text:
export ft_getstatus

function ft_getstatus(ft_handle::Culong)
  amountinrxqueue = Ref{Cuint}()
  amountintxqueue = Ref{Cuint}()
  eventstatus = Ref{Cuint}()
  ft_status = ccall((:FT_GetStatus, d2xx),
                     Cuint,
                     (Culong,Ref{Cuint},Ref{Cuint},Ref{Cuint}),
                     ft_handle, amountinrxqueue, amountintxqueue, eventstatus)
  checkstatus(ft_status)
  return (convert(Int32,amountinrxqueue[]),
          convert(Int32,amountintxqueue[]),convert(Int32,eventstatus[]))
end


text:
using LinearAlgebra
using NLopt
using Printf

mutable struct ExpFitBuffer
    p::Int
    Nt::Int
    Neval::Int

    r::Vector{Float64}
    t::Vector{Float64}
    w::Vector{Float64}
    wt::Vector{Float64}
    wrt::Vector{Float64}

    A::Vector{ComplexF64}
    B::Vector{ComplexF64}
    C::Vector{Float64}
    E::Vector{ComplexF64}
    G::Matrix{ComplexF64}
    M::Matrix{ComplexF64}
    N::Matrix{ComplexF64}
    Q::Matrix{ComplexF64}
    H::AbstractMatrix{ComplexF64}
    Bt::Matrix{ComplexF64}
    Θ::Matrix{ComplexF64}

    function ExpFitBuffer(r::Vector, t::AbstractVector, p::Integer, w::Vector)
        Nt = length(r)
        @assert Nt == length(t)

        r = copy(r)
        t = copy(t)
        wt = w.*t
        wrt = wt.*r

        A = fill(0.0im, p)
        C = fill(0.0, p)
        B = fill(0.0im, p)
        E = fill(0.0im, p)
        G = fill(0.0im, p, p)
        M = fill(0.0im, p, p)
        N = fill(0.0im, p, p)
        Q = fill(0.0im, p, p)
        if p == 1
            H = fill(0.0im, p, p)
        else
            H = Tridiagonal(M)
        end
        Bt = fill(0.0im, Nt, p)
        Θ = fill(0.0im, Nt, p)

        new(p, Nt, 0, r, t, w, wt, wrt, A, B, C, E, G, M, N, Q, H, Bt, Θ)
    end
end

function ARMA_gradient(grad::Vector, buffer::ExpFitBuffer, A::Vector{T},
    B::Vector{T}, fmodel::Vector) where {T<:Number}
    t = buffer.t
    p = buffer.p

    for i = 1:p
        buffer.Bt[:,i] = B[i].^t
        buffer.Θ[:,i] = A[i]/B[i] * t .* buffer.Bt[:,i]
    end

    buffer.M[:,:] = buffer.Bt'*(buffer.w.*buffer.Bt)
    buffer.N[:,:] = buffer.Bt'*(buffer.wt.*buffer.Bt)
    for i=1:p; buffer.N[i,:] /= B[i]; end

    for i=1:p
        buffer.E[i] = dot(buffer.wt, buffer.Bt[:,i])/B[i]
    end

    for i=1:p
        for j=1:p
            buffer.Q[i,j] = A[i] * buffer.N[i,j]
        end
    end

    buffer.G[:,:] = pinv(buffer.M) * (Diagonal(buffer.E-buffer.N*A)-buffer.Q)

    H = buffer.H
    for i=1:2:p-1
        ID = 1.0/(B[i]-B[i+1])
        H[i,i] = -B[i]*ID
        H[i+1,i] = B[i+1]*ID
        H[i,i+1] = -ID
        H[i+1,i+1] = ID
    end
    if p%2 == 1
        H[end,end] = B[end]^2
    end

    Jacobian = real((buffer.Θ + buffer.Bt*buffer.G)*buffer.H)
    grad[:] = 2*(Jacobian' * (buffer.w.*(fmodel-buffer.r)))
end

function ARMA_objective(C::Vector, grad::Vector, buffer::ExpFitBuffer)
    p = length(C)
    B = C2B(C)
    A = findA(buffer.t, buffer.r, B, w=buffer.w)
    fmodel = exponential_model(buffer.t, A, B)

    if length(grad) > 0
        ARMA_gradient(grad, buffer, A, B, fmodel)
    end
    buffer.Neval += 1
    sum(buffer.w .* (fmodel - buffer.r).^2)
end

function optimize_exponentials(data::Vector, w::Vector, guessC::Vector)
    p = length(guessC)
    if p == 0
        return sum(data[2:end].^2 .* w[2:end]), Float64[]
    end

    opt = Opt(:LD_MMA, p)
    opt = Opt(:LN_BOBYQA, p)
    xtol_abs!(opt, 1e-6)
    ftol_rel!(opt, 1e-5)
    maxeval!(opt, 100p)

    # Constraint bounds on the quadratic representation of the bases.
    # As all bases are in the unit circle, the linear coefficients are in (-2,2)
    # and the quadratic ones are in (-1,1).
    lb = fill(-2.0, size(guessC))
    ub = fill(2.0, size(guessC))
    lb[2:2:end] .= -1
    ub[2:2:end] .= +1
    if p%2 == 1
        lb[end] = -Inf
        ub[end] = +Inf
    end
    lower_bounds!(opt, lb)
    upper_bounds!(opt, ub)
    step = 0.1*(ub-lb)
    if step[end] > 1
        step[end] = guessC[end]*0.5
    end
    initial_step!(opt, step)

    t = 0:length(data)-1
    buffer = ExpFitBuffer(data, t, p, w)
    min_objective!(opt, (x,grad)->ARMA_objective(x, grad, buffer))
    minf1 = ARMA_objective(guessC, [], buffer)
    # println("Before opt: $minf1 at $guessC")
    try
        (minf,minx,ret) = optimize(opt, guessC)
        # println("After opt: $minf at $minx after $(buffer.Neval) iterations (returned $ret)")
        # println("Moved by $(minx-guessC), improved by $(minf1-minf)\n")
        return minf, minx
    catch e
        @printf("Error in %d-order fit: %s\n", p, e)
        return minf1, guessC
    end
end


text:
# This file is part of the Steganography package
# (http://github.com/davidssmith/Steganography.jl).
#
# The MIT License (MIT)
#
# Copyright (c) 2017 David Smith
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

module Steganography

using Compat

export embed, extract, setlastbits, getlastbits, setlast7, setlast8,
    getlast7, getlast8

const version = v"0.0.2"

function flattencomplex{T<:Complex,N}(x::Array{T,N})
    if T == Complex128
        T2 = Float64
    elseif T == Complex64
        T2 = Float32
    elseif T == Complex32
        T2 = Float16
    else
        error("type $T unhandled in Steganography.flattencomplex()")
    end
    y = Array{T2}(2length(x))
    for k in 1:length(x)
        y[2k-1] = real(x[k])
        y[2k] = imag(x[k])
    end
    return y
end
function unflattencomplex{T<:Real,N}(y::Array{T,N})
    if T == Float64
        T2 = Complex128
    elseif T == Float32
        T2 = Complex64
    elseif T == Float16
        T2 = Complex32
    else
        error("type $T unhandled in Steganography.unflattencomplex()")
    end
    x = Array{T2}(div(length(y),2))
    for k in 1:length(x)
        x[k] = complex(y[2k-1],y[2k])
    end
    return x
end

@compat function setlastbits{T<:Integer}(i::T, n::UInt8, nbits::UInt8)
    S = typeof(i)
    j = (i >> nbits) << nbits
    j | (n & ((S(1) << nbits) - S(1)))
end
@compat setlastbits{T<:AbstractFloat}(x::T, n::UInt8, nbits::UInt8) = reinterpret(T, setlastbits(reinterpret(Unsigned, x), n, nbits))
@compat setlast8{T}(x::T, n::UInt8) = setlastbits(x, n, UInt8(8))
@compat setlast7{T}(x::T, n::UInt8) = setlastbits(x, n, UInt8(7))

@compat function getlastbits{T}(i::T, nbits::UInt8)
    S = typeof(i)
    return UInt8(i & ((S(1) << nbits) - S(1)))
end
@compat getlastbits{T<:AbstractFloat}(x::T, nbits::UInt8) = getlastbits(reinterpret(Unsigned, x), nbits)

@compat getlast8{T<:AbstractFloat}(x::T) = UInt8(reinterpret(Unsigned, x) & 0xff)
@compat getlast8{T}(x::T) = UInt8(x & 0xff)   # catch all

@compat getlast7{T<:AbstractFloat}(x::T) = UInt8(reinterpret(Unsigned, x) & 0x7f)
@compat getlast7{T}(x::T) = UInt8(x & 0x7f)   # catch all

@compat function embed{T,N}(data::Array{T,N}, text::Array{UInt8,1}; ignorenonascii::Bool=true)
    if T <: Complex
        data = flattencomplex(data)
    end
    @assert length(text) <= length(data)
    y = copy(data)   # make sure we have enough space
    for j in 1:length(text)
        @assert text[j] != 0x04
        if !ignorenonascii
            @assert text[j] <= 0x7f
        end
        if text[j] > 0x7f
            println(text[j], " ", Char(text[j]), " ", hex(text[j]))
            y[j] = setlast7(data[j], UInt8(0))
        else
            y[j] = setlast7(data[j], text[j])
        end
    end
    if length(text) < length(data)
        y[length(text)+1] = setlast7(data[length(text)+1], 0x04) # ASCII 0x04 means 'end of transmission'
    end
    if T <: Complex
        return unflattencomplex(y)
    else
        return y
    end
end

@compat function extract{T,N}(s::Array{T,N})
    if T <: Complex
        s = flattencomplex(s)
   end
    t = Array{UInt8}(length(s))
    k = 1
    while true
        t[k] = getlast7(s[k])
        if t[k] == 0x04
            break
        end
        k += 1
    end
    t[1:k-1]
end

end


text:
cond(c, t, f) = c ? t() : f()

allsuccs(ir, (id,chs)) =
  union(map(b -> b.id, successors(block(ir, id))),
        allsuccs.((ir,), chs)...)

function captures(ir, (id, chs), cs = Dict())
  bl = block(ir, id)
  captures.((ir,), chs, (cs,))
  cs[id] = setdiff(union([cs[i] for (i, _) in chs]..., usages(bl)), definitions(bl))
  return cs
end

function return_thunk(x)
  ir = IR()
  return!(ir, xcall(:getindex, argument!(ir), 1))
  Expr(:lambda, ir, x)
end

function br_thunk(args...)
  ir = IR()
  self = argument!(ir)
  return!(ir, xcall([xcall(:getindex, self, i) for i = 1:length(args)]...))
  Expr(:lambda, ir, args...)
end

function functionalbranches!(bl, pr, labels)
  if length(branches(bl)) == 1
    br = branches(bl)[1]
    if !isreturn(br)
      r = push!(pr, Expr(:call, labels[br.block], br.args...))
      empty!(branches(pr.to))
      return!(pr, r)
    end
  else
    @assert length(branches(bl)) == 2
    f, t = branches(bl)
    function brfunc(br)
      isreturn(br) && return return_thunk(returnvalue(br))
      isempty(arguments(br)) && return labels[br.block]
      br_thunk(labels[br.block], arguments(br)...)
    end
    r = push!(pr, xcall(IRTools, :cond, f.condition, brfunc(t), brfunc(f)))
    empty!(branches(pr.to))
    return!(pr, r)
  end
end

function _functional(ir, tree, vars = [], cs = captures(ir, tree))
  id = tree[1]
  bl = IR(block(ir, id))
  labels = Dict()
  labels[id] = self = id == 1 ? arguments(bl)[1] : argument!(bl, at = 1)
  pr = Pipe(bl)
  for (i, v) in enumerate(vars)
    v′ = push!(pr, xcall(:getindex, self, i))
    v isa Integer ? (labels[v] = v′) : substitute!(pr, v, substitute(pr, v′))
  end
  for _ in pr end
  for t in reverse(tree[2])
    bs = filter(id -> haskey(labels, id), allsuccs(ir, t))
    λ = Expr(:lambda, _functional(ir, t, [bs..., cs[t[1]]...]),
             map(b -> labels[b], bs)..., cs[t[1]]...)
    labels[t[1]] = push!(pr, λ)
  end
  functionalbranches!(bl, pr, labels)
  return finish(pr)
end

functional(ir) = _functional(explicitbranch!(ir), domtree(ir))


text:
#==============================================================================#
# DynamoDBStreams.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/streams.dynamodb-2012-08-10.normal.json
#==============================================================================#

__precompile__()

module DynamoDBStreams

using AWSCore

"""
    using AWSSDK.DynamoDBStreams.describe_stream
    describe_stream([::AWSConfig], arguments::Dict)
    describe_stream([::AWSConfig]; StreamArn=, <keyword arguments>)

    using AWSCore.Services.streams_dynamodb
    streams_dynamodb([::AWSConfig], "DescribeStream", arguments::Dict)
    streams_dynamodb([::AWSConfig], "DescribeStream", StreamArn=, <keyword arguments>)

# DescribeStream Operation

Returns information about a stream, including the current status of the stream, its Amazon Resource Name (ARN), the composition of its shards, and its corresponding DynamoDB table.

**Note**
> You can call `DescribeStream` at a maximum rate of 10 times per second.

Each shard in the stream has a `SequenceNumberRange` associated with it. If the `SequenceNumberRange` has a `StartingSequenceNumber` but no `EndingSequenceNumber`, then the shard is still open (able to receive more stream records). If both `StartingSequenceNumber` and `EndingSequenceNumber` are present, then that shard is closed and can no longer receive more data.

# Arguments

## `StreamArn = ::String` -- *Required*
The Amazon Resource Name (ARN) for the stream.

## `Limit = ::Int`
The maximum number of shard objects to return. The upper limit is 100.

## `ExclusiveStartShardId = ::String`
The shard ID of the first item that this operation will evaluate. Use the value that was returned for `LastEvaluatedShardId` in the previous operation.

# Returns

`DescribeStreamOutput`

# Exceptions

`ResourceNotFoundException` or `InternalServerError`.

# Example: To describe a stream with a given stream ARN

The following example describes a stream with a given stream ARN.

Input:
```
[
    "StreamArn" => "arn:aws:dynamodb:us-west-2:111122223333:table/Forum/stream/2015-05-20T20:51:10.252"
]
```

Output:
```
Dict(
    "StreamDescription" => Dict(
        "CreationRequestDateTime" => "Wed May 20 13:51:10 PDT 2015",
        "KeySchema" => [
            Dict(
                "AttributeName" => "ForumName",
                "KeyType" => "HASH"
            ),
            Dict(
                "AttributeName" => "Subject",
                "KeyType" => "RANGE"
            )
        ],
        "Shards" => [
            Dict(
                "SequenceNumberRange" => Dict(
                    "EndingSequenceNumber" => "20500000000000000910398",
                    "StartingSequenceNumber" => "20500000000000000910398"
                ),
                "ShardId" => "shardId-00000001414562045508-2bac9cd2"
            ),
            Dict(
                "ParentShardId" => "shardId-00000001414562045508-2bac9cd2",
                "SequenceNumberRange" => Dict(
                    "EndingSequenceNumber" => "820400000000000001192334",
                    "StartingSequenceNumber" => "820400000000000001192334"
                ),
                "ShardId" => "shardId-00000001414576573621-f55eea83"
            ),
            Dict(
                "ParentShardId" => "shardId-00000001414576573621-f55eea83",
                "SequenceNumberRange" => Dict(
                    "EndingSequenceNumber" => "1683700000000000001135967",
                    "StartingSequenceNumber" => "1683700000000000001135967"
                ),
                "ShardId" => "shardId-00000001414592258131-674fd923"
            ),
            Dict(
                "ParentShardId" => "shardId-00000001414592258131-674fd923",
                "SequenceNumberRange" => Dict(
                    "StartingSequenceNumber" => "2574600000000000000935255"
                ),
                "ShardId" => "shardId-00000001414608446368-3a1afbaf"
            )
        ],
        "StreamArn" => "arn:aws:dynamodb:us-west-2:111122223333:table/Forum/stream/2015-05-20T20:51:10.252",
        "StreamLabel" => "2015-05-20T20:51:10.252",
        "StreamStatus" => "ENABLED",
        "StreamViewType" => "NEW_AND_OLD_IMAGES",
        "TableName" => "Forum"
    )
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/streams-dynamodb-2012-08-10/DescribeStream)
"""
@inline describe_stream(aws::AWSConfig=default_aws_config(); args...) = describe_stream(aws, args)

@inline describe_stream(aws::AWSConfig, args) = AWSCore.Services.streams_dynamodb(aws, "DescribeStream", args)

@inline describe_stream(args) = describe_stream(default_aws_config(), args)

"""
    using AWSSDK.DynamoDBStreams.get_records
    get_records([::AWSConfig], arguments::Dict)
    get_records([::AWSConfig]; ShardIterator=, <keyword arguments>)

    using AWSCore.Services.streams_dynamodb
    streams_dynamodb([::AWSConfig], "GetRecords", arguments::Dict)
    streams_dynamodb([::AWSConfig], "GetRecords", ShardIterator=, <keyword arguments>)

# GetRecords Operation

Retrieves the stream records from a given shard.

Specify a shard iterator using the `ShardIterator` parameter. The shard iterator specifies the position in the shard from which you want to start reading stream records sequentially. If there are no stream records available in the portion of the shard that the iterator points to, `GetRecords` returns an empty list. Note that it might take multiple calls to get to a portion of the shard that contains stream records.

**Note**
> `GetRecords` can retrieve a maximum of 1 MB of data or 1000 stream records, whichever comes first.

# Arguments

## `ShardIterator = ::String` -- *Required*
A shard iterator that was retrieved from a previous GetShardIterator operation. This iterator can be used to access the stream records in this shard.

## `Limit = ::Int`
The maximum number of records to return from the shard. The upper limit is 1000.

# Returns

`GetRecordsOutput`

# Exceptions

`ResourceNotFoundException`, `LimitExceededException`, `InternalServerError`, `ExpiredIteratorException` or `TrimmedDataAccessException`.

# Example: To retrieve all the stream records from a shard

The following example retrieves all the stream records from a shard.

Input:
```
[
    "ShardIterator" => "arn:aws:dynamodb:us-west-2:111122223333:table/Forum/stream/2015-05-20T20:51:10.252|1|AAAAAAAAAAEvJp6D+zaQ...  <remaining characters omitted> ..."
]
```

Output:
```
Dict(
    "NextShardIterator" => "arn:aws:dynamodb:us-west-2:111122223333:table/Forum/stream/2015-05-20T20:51:10.252|1|AAAAAAAAAAGQBYshYDEe ... <remaining characters omitted> ...",
    "Records" => [
        Dict(
            "awsRegion" => "us-west-2",
            "dynamodb" => Dict(
                "ApproximateCreationDateTime" => "1.46480646E9",
                "Keys" => Dict(
                    "ForumName" => Dict(
                        "S" => "DynamoDB"
                    ),
                    "Subject" => Dict(
                        "S" => "DynamoDB Thread 3"
                    )
                ),
                "SequenceNumber" => "300000000000000499659",
                "SizeBytes" => 41,
                "StreamViewType" => "KEYS_ONLY"
            ),
            "eventID" => "e2fd9c34eff2d779b297b26f5fef4206",
            "eventName" => "INSERT",
            "eventSource" => "aws:dynamodb",
            "eventVersion" => "1.0"
        ),
        Dict(
            "awsRegion" => "us-west-2",
            "dynamodb" => Dict(
                "ApproximateCreationDateTime" => "1.46480527E9",
                "Keys" => Dict(
                    "ForumName" => Dict(
                        "S" => "DynamoDB"
                    ),
                    "Subject" => Dict(
                        "S" => "DynamoDB Thread 1"
                    )
                ),
                "SequenceNumber" => "400000000000000499660",
                "SizeBytes" => 41,
                "StreamViewType" => "KEYS_ONLY"
            ),
            "eventID" => "4b25bd0da9a181a155114127e4837252",
            "eventName" => "MODIFY",
            "eventSource" => "aws:dynamodb",
            "eventVersion" => "1.0"
        ),
        Dict(
            "awsRegion" => "us-west-2",
            "dynamodb" => Dict(
                "ApproximateCreationDateTime" => "1.46480646E9",
                "Keys" => Dict(
                    "ForumName" => Dict(
                        "S" => "DynamoDB"
                    ),
                    "Subject" => Dict(
                        "S" => "DynamoDB Thread 2"
                    )
                ),
                "SequenceNumber" => "500000000000000499661",
                "SizeBytes" => 41,
                "StreamViewType" => "KEYS_ONLY"
            ),
            "eventID" => "740280c73a3df7842edab3548a1b08ad",
            "eventName" => "REMOVE",
            "eventSource" => "aws:dynamodb",
            "eventVersion" => "1.0"
        )
    ]
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/streams-dynamodb-2012-08-10/GetRecords)
"""
@inline get_records(aws::AWSConfig=default_aws_config(); args...) = get_records(aws, args)

@inline get_records(aws::AWSConfig, args) = AWSCore.Services.streams_dynamodb(aws, "GetRecords", args)

@inline get_records(args) = get_records(default_aws_config(), args)

"""
    using AWSSDK.DynamoDBStreams.get_shard_iterator
    get_shard_iterator([::AWSConfig], arguments::Dict)
    get_shard_iterator([::AWSConfig]; StreamArn=, ShardId=, ShardIteratorType=, <keyword arguments>)

    using AWSCore.Services.streams_dynamodb
    streams_dynamodb([::AWSConfig], "GetShardIterator", arguments::Dict)
    streams_dynamodb([::AWSConfig], "GetShardIterator", StreamArn=, ShardId=, ShardIteratorType=, <keyword arguments>)

# GetShardIterator Operation

Returns a shard iterator. A shard iterator provides information about how to retrieve the stream records from within a shard. Use the shard iterator in a subsequent `GetRecords` request to read the stream records from the shard.

**Note**
> A shard iterator expires 15 minutes after it is returned to the requester.

# Arguments

## `StreamArn = ::String` -- *Required*
The Amazon Resource Name (ARN) for the stream.

## `ShardId = ::String` -- *Required*
The identifier of the shard. The iterator will be returned for this shard ID.

## `ShardIteratorType = "TRIM_HORIZON", "LATEST", "AT_SEQUENCE_NUMBER" or "AFTER_SEQUENCE_NUMBER"` -- *Required*
Determines how the shard iterator is used to start reading stream records from the shard:

*   `AT_SEQUENCE_NUMBER` - Start reading exactly from the position denoted by a specific sequence number.

*   `AFTER_SEQUENCE_NUMBER` - Start reading right after the position denoted by a specific sequence number.

*   `TRIM_HORIZON` - Start reading at the last (untrimmed) stream record, which is the oldest record in the shard. In DynamoDB Streams, there is a 24 hour limit on data retention. Stream records whose age exceeds this limit are subject to removal (trimming) from the stream.

*   `LATEST` - Start reading just after the most recent stream record in the shard, so that you always read the most recent data in the shard.

## `SequenceNumber = ::String`
The sequence number of a stream record in the shard from which to start reading.

# Returns

`GetShardIteratorOutput`

# Exceptions

`ResourceNotFoundException`, `InternalServerError` or `TrimmedDataAccessException`.

# Example: To obtain a shard iterator for the provided stream ARN and shard ID

The following example returns a shard iterator for the provided stream ARN and shard ID.

Input:
```
[
    "ShardId" => "00000001414576573621-f55eea83",
    "ShardIteratorType" => "TRIM_HORIZON",
    "StreamArn" => "arn:aws:dynamodb:us-west-2:111122223333:table/Forum/stream/2015-05-20T20:51:10.252"
]
```

Output:
```
Dict(
    "ShardIterator" => "arn:aws:dynamodb:us-west-2:111122223333:table/Forum/stream/2015-05-20T20:51:10.252|1|AAAAAAAAAAEvJp6D+zaQ...  <remaining characters omitted> ..."
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/streams-dynamodb-2012-08-10/GetShardIterator)
"""
@inline get_shard_iterator(aws::AWSConfig=default_aws_config(); args...) = get_shard_iterator(aws, args)

@inline get_shard_iterator(aws::AWSConfig, args) = AWSCore.Services.streams_dynamodb(aws, "GetShardIterator", args)

@inline get_shard_iterator(args) = get_shard_iterator(default_aws_config(), args)

"""
    using AWSSDK.DynamoDBStreams.list_streams
    list_streams([::AWSConfig], arguments::Dict)
    list_streams([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.streams_dynamodb
    streams_dynamodb([::AWSConfig], "ListStreams", arguments::Dict)
    streams_dynamodb([::AWSConfig], "ListStreams", <keyword arguments>)

# ListStreams Operation

Returns an array of stream ARNs associated with the current account and endpoint. If the `TableName` parameter is present, then `ListStreams` will return only the streams ARNs for that table.

**Note**
> You can call `ListStreams` at a maximum rate of 5 times per second.

# Arguments

## `TableName = ::String`
If this parameter is provided, then only the streams associated with this table name are returned.

## `Limit = ::Int`
The maximum number of streams to return. The upper limit is 100.

## `ExclusiveStartStreamArn = ::String`
The ARN (Amazon Resource Name) of the first item that this operation will evaluate. Use the value that was returned for `LastEvaluatedStreamArn` in the previous operation.

# Returns

`ListStreamsOutput`

# Exceptions

`ResourceNotFoundException` or `InternalServerError`.

# Example: To list all of the stream ARNs

The following example lists all of the stream ARNs.

Input:
```
[

]
```

Output:
```
Dict(
    "Streams" => [
        Dict(
            "StreamArn" => "arn:aws:dynamodb:us-wesst-2:111122223333:table/Forum/stream/2015-05-20T20:51:10.252",
            "StreamLabel" => "2015-05-20T20:51:10.252",
            "TableName" => "Forum"
        ),
        Dict(
            "StreamArn" => "arn:aws:dynamodb:us-west-2:111122223333:table/Forum/stream/2015-05-20T20:50:02.714",
            "StreamLabel" => "2015-05-20T20:50:02.714",
            "TableName" => "Forum"
        ),
        Dict(
            "StreamArn" => "arn:aws:dynamodb:us-west-2:111122223333:table/Forum/stream/2015-05-19T23:03:50.641",
            "StreamLabel" => "2015-05-19T23:03:50.641",
            "TableName" => "Forum"
        )
    ]
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/streams-dynamodb-2012-08-10/ListStreams)
"""
@inline list_streams(aws::AWSConfig=default_aws_config(); args...) = list_streams(aws, args)

@inline list_streams(aws::AWSConfig, args) = AWSCore.Services.streams_dynamodb(aws, "ListStreams", args)

@inline list_streams(args) = list_streams(default_aws_config(), args)

end # module DynamoDBStreams

#==============================================================================#
# End of file
#==============================================================================#


text:
function alternate_case(c)
   global upper
   if upper
      c = uppercase(c)
      upper = false
   else
      c = lowercase(c)
      upper = true
   end
   return c
end

while true
   global upper

   print("> ")

   in = readline()

   upper = false

   println(join(map(alternate_case, split(in, "")), ""))
end

max_stars_repo_name: Coalpha/lang-experimentation

max_issues_repo_name: Coalpha/lang-experimentation

max_forks_repo_name: Coalpha/lang-experimentation


text:
using CUDA, Base.Test

@test devcount() > 0

dev = CuDevice(0)
if capability(dev) < v"2.0"
    warn("native execution not supported on SM < 2.0")
else
    include("native.jl")
end


text:
# The eigensystem when no symmetry is present
function eigensystem(ω::T, medium::PhysicalMedium{T,3}, species::Species{T,3}, ::WithoutSymmetry{3};
        basis_order::Int = 2,
        basis_field_order::Int = 2*basis_order,
        numberofparticles::Number = Inf,
        kws...) where {T<:AbstractFloat}

    scale_number_density = one(T) - one(T) / numberofparticles

    k = ω/medium.c
    sps = species

    S = length(sps)
    L = basis_order
    L1 = basis_field_order
    len = (L1+1)^2 * (L+1)^2 * S
    MM_mat = Matrix{Complex{T}}(undef,len,len)

    t_matrices = get_t_matrices(medium, sps, ω, L)
    t_diags = diag.(t_matrices)
    len(order::Int) = basisorder_to_basislength(Acoustic{T,3},order)

    as = [
        s1.exclusion_distance * outer_radius(s1) + s2.exclusion_distance * outer_radius(s2)
    for s1 in sps, s2 in sps]

    function M_component(keff,Ns,l,m,l2,m2,s1,dl,dm,l1,m1,s2)::Complex{T}
        minl3 = max(abs(m1-m2),abs(dl-l),abs(l1-l2))
        maxl3 = min(abs(dl+l),abs(l1+l2))

        (m == dm && l == dl && m1 == m2 && l1 == l2 && s1 == s2 ? 1.0 : 0.0) +
        if minl3 <= maxl3
            as[s1,s2] * scale_number_density * number_density(sps[s2]) * t_diags[s1][len(l)] *
            sum(l3 ->
                gaunt_coefficient(l,m,dl,dm,l3,m1-m2) *
                gaunt_coefficient(l1,m1,l2,m2,l3,m1-m2) * Ns[l3+1,s1,s2]
            , minl3:maxl3)
        else
            zero(Complex{T})
        end
    end

    function MM(keff::Complex{T})::Matrix{Complex{T}}
        Ns = [
            kernelN3D(l3,k*as[s1,s2],keff*as[s1,s2])
        for l3 = 0:min(2L1,2L), s1 = 1:S, s2 = 1:S] ./  (keff^2.0 - k^2.0)

        # The order of the indices below is important
        ind2 = 1
        for s2 = 1:S for dl = 0:L for dm = -dl:dl for l1 = 0:L1 for m1 = -l1:l1
            ind1 = 1
            for s1 = 1:S for l = 0:L for m = -l:l for l2 = 0:L1 for m2 = -l2:l2
                MM_mat[ind1, ind2] = M_component(keff,Ns,l,m,l2,m2,s1,dl,dm,l1,m1,s2)
                ind1 += 1
            end end end end end
            ind2 += 1
        end end end end end
        return MM_mat
    end

    return MM
end

function eigensystem(ω::T, medium::PhysicalMedium{T,3}, species::Species{T,3}, ::AbstractAzimuthalSymmetry;
        basis_order::Int = 2,
        basis_field_order::Int = 2*basis_order,
        numberofparticles::Number = Inf,
        kws...) where {T<:AbstractFloat}

    scale_number_density = one(T) - one(T) / numberofparticles

    k = ω/medium.c
    sps = species

    S = length(sps)
    L = basis_order
    L1 = basis_field_order

    len = Int(1 - L*(2 + L)*(L - 3*L1 - 2)/3 + L1) * S
    MM_mat = Matrix{Complex{T}}(undef,len,len)

    t_matrices = get_t_matrices(medium, sps, ω, L)
    t_diags = diag.(t_matrices)
    len(order::Int) = basisorder_to_basislength(Acoustic{T,3},order)

    as = [
        s1.exclusion_distance * outer_radius(s1) + s2.exclusion_distance * outer_radius(s2)
    for s1 in sps, s2 in sps]

    # the index for the T-matrix below needs to be changed when seperating correctly the 2D and 3D case.
    function M_component(keff,Ns,l,m,l2,s1,dl,dm,l1,s2)::Complex{T}
        minl3 = max(abs(m-dm),abs(dl-l),abs(l1-l2))
        maxl3 = min(dl+l,l1+l2)

        (m == dm && l == dl && l1 == l2 && s1 == s2 ? 1.0 : 0.0) +
        if minl3 <= maxl3
            as[s1,s2] * scale_number_density * number_density(species[s2]) * t_diags[s1][len(l)] *
            sum(l3 ->
                gaunt_coefficient(l,m,dl,dm,l3,m-dm) *
                gaunt_coefficient(l1,-dm,l2,-m,l3,m-dm) * Ns[l3+1,s1,s2]
            , minl3:maxl3)
        else
            zero(Complex{T})
        end
    end

    # The order of the indices below is important
    function MM(keff::Complex{T})::Matrix{Complex{T}}
        Ns = [
            kernelN3D(l3,k*as[s1,s2],keff*as[s1,s2])
        for l3 = 0:min(2L1,2L), s1 = 1:S, s2 = 1:S] ./ (keff^2.0 - k^2.0)

        ind2 = 1
        for s2 = 1:S for dl = 0:L for dm = -dl:dl for l1 = abs(dm):L1
            ind1 = 1
            for s1 = 1:S for l = 0:L for m = -l:l for l2 = abs(m):L1
                MM_mat[ind1, ind2] = M_component(keff,Ns,l,m,l2,s1,dl,dm,l1,s2)
                ind1 += 1
            end end end end
            ind2 += 1
        end end end end
        return MM_mat
    end

    return MM
end

function eigensystem(ω::T, medium::PhysicalMedium{T,3}, species::Species{T,3}, ::RadialSymmetry{3};
        basis_order::Int = 2,
        kws...) where {T<:AbstractFloat}

        MM = eigensystem(ω, medium, species, PlanarAzimuthalSymmetry{3}(); basis_order = basis_order, kws...)

        factors = [Complex{T}(im)^l * sqrt(T(2l + 1)) for l = 0:basis_order, s = 1:length(species)][:]

        MMR(keff::Complex{T}) = MM(keff) * diagm(factors)

        return MMR
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -42.99754755332367 <= q <= 40.38858433100043)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.8925450194703983 + 0.6110890427800215*1/(1 + exp(-(-0.16626426792027038 + 0.44845618413236066*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + 0.6799320945085388*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + -0.7514544713734312*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + 0.795720190818876*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + 0.1740395626916258*1/(1 + exp(-(0.8063063231762757 + 0.4905142959916642*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + 0.7092497630992938*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + 0.8897979792487347*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + -0.11224890387090136*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + -0.8017072131372296*1/(1 + exp(-(0.9210510999614945 + -0.4935968211971775*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + -0.33547745483400737*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + 0.1703969445706579*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + -0.2742027660464488*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + -0.7966067754659663*1/(1 + exp(-(-0.7328670881809516 + 0.6102429954796835*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + -0.6353392761755483*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + -0.7830650903096812*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + 0.1627621376749886*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))))))))) + 1/(1 + exp(-(-0.7008872510227278 + -0.258173144464545*1/(1 + exp(-(-0.16626426792027038 + 0.44845618413236066*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + 0.6799320945085388*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + -0.7514544713734312*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + 0.795720190818876*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + 0.49550787111793193*1/(1 + exp(-(0.8063063231762757 + 0.4905142959916642*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + 0.7092497630992938*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + 0.8897979792487347*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + -0.11224890387090136*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + -0.8462943069272191*1/(1 + exp(-(0.9210510999614945 + -0.4935968211971775*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + -0.33547745483400737*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + 0.1703969445706579*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + -0.2742027660464488*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + -0.05869131996649157*1/(1 + exp(-(-0.7328670881809516 + 0.6102429954796835*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + -0.6353392761755483*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + -0.7830650903096812*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + 0.1627621376749886*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))))))))) + 1/(1 + exp(-(-0.2643911691309153 + -0.31660511055293883*1/(1 + exp(-(-0.16626426792027038 + 0.44845618413236066*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + 0.6799320945085388*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + -0.7514544713734312*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + 0.795720190818876*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + -0.09136517532701127*1/(1 + exp(-(0.8063063231762757 + 0.4905142959916642*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + 0.7092497630992938*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + 0.8897979792487347*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + -0.11224890387090136*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + -0.36400896662330684*1/(1 + exp(-(0.9210510999614945 + -0.4935968211971775*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + -0.33547745483400737*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + 0.1703969445706579*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + -0.2742027660464488*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + 0.4876638787119827*1/(1 + exp(-(-0.7328670881809516 + 0.6102429954796835*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + -0.6353392761755483*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + -0.7830650903096812*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + 0.1627621376749886*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))))))))) + 1/(1 + exp(-(0.6679332747384503 + -0.6061375406574494*1/(1 + exp(-(-0.16626426792027038 + 0.44845618413236066*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + 0.6799320945085388*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + -0.7514544713734312*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + 0.795720190818876*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + -0.7532640510772808*1/(1 + exp(-(0.8063063231762757 + 0.4905142959916642*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + 0.7092497630992938*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + 0.8897979792487347*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + -0.11224890387090136*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + 0.9959397069462161*1/(1 + exp(-(0.9210510999614945 + -0.4935968211971775*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + -0.33547745483400737*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + 0.1703969445706579*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + -0.2742027660464488*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + -0.8544693505755525*1/(1 + exp(-(-0.7328670881809516 + 0.6102429954796835*1/(1 + exp(-(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5])))) + -0.6353392761755483*1/(1 + exp(-(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5])))) + -0.7830650903096812*1/(1 + exp(-(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5])))) + 0.1627621376749886*1/(1 + exp(-(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5])))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function swap!(i, j, x)
  tmp = deepcopy(x[i])
  x[i] = deepcopy(x[j])
  x[j] = tmp
  return x
end

function println_mu_sig2(states, i)
  println("    mu*0 for state $(i): $(-cumsum(states[i].theta.delta[false]))")
  println("    mu*1 for state $(i): $(cumsum(states[i].theta.delta[true]))")
  println("    sig2 for state $(i): $(states[i].theta.sig2)")
end

function swapchains!(states, loglike, temperatures;
                     paircounts, swapcounts, verbose=0, randpair=0.0)
  """
  randpair: The proportion of time to propose swapping random pairs.
            If set to 0.0, then pairs are swapped from hottest to coolest
            or coolest to hottest. If set to 1.0, then random pairs are
            formed and proposed to be swapped.

  See: https://academic.oup.com/gji/article/196/1/357/585739
  """
  nchains = length(states)
  @assert (nchains == length(temperatures)) && (mod(nchains, 2) == 0)

  pairs = if randpair > rand()
    Iterators.partition(Random.shuffle(1:nchains), 2)
  else
    p = zip(1:(nchains - 1), 2:nchains)
    rand(Bool) ? p : Iterators.reverse(p)
  end

  for (i, j) in pairs
    # Increment pair counts
    paircounts[i, j] += 1
    paircounts[j, i] += 1

    si, sj = states[i].theta, states[j].theta
    ti, tj = temperatures[i], temperatures[j]
    log_accept_ratio = MCMC.WSPT.compute_log_accept_ratio(loglike,
                                                          (si, sj),
                                                          (ti, tj),
                                                          verbose=verbose-2)
    should_swap_chains = log_accept_ratio > log(rand())

    if verbose > 1
      println("swap log accept ratio ($(i), $(j)): $(log_accept_ratio)")
      if verbose > 2
        println_mu_sig2(states, i)
        println_mu_sig2(states, j)
      end
    end

    if should_swap_chains
      swap!(i, j, states)
      # Increment swap-counts matrix
      swapcounts[i, j] += 1
      swapcounts[j, i] += 1
      if verbose > 0
        println("Swapped chains $(i) and $(j)")
      end
    end
  end

  return
end


text:
############################################################
# joAbstractFosterLinearOperator subtypes ##################
############################################################

############################################################
# joLooseMatrix #################################################

export joLooseMatrix, joLooseMatrixException

# type definition
"""
joLooseMatrix type

# TYPE PARAMETERS
- DDT::DataType : domain DataType
- RDT::DataType : range DataType

# FIELDS
- name::String : given name
- m::Integer : # of rows
- n::Integer : # of columns
- fop::Function : forward matrix
- fop_T::Function : transpose matrix
- fop_A::Function : adjoint matrix
- fop_C::Function : conj matrix
- iop::Nullable{Function} : inverse for fop
- iop_T::Nullable{Function} : inverse for fop_T
- iop_A::Nullable{Function} : inverse for fop_A
- iop_C::Nullable{Function} : inverse for fop_C

"""
struct joLooseMatrix{DDT<:Number,RDT<:Number} <: joAbstractFosterLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    fop::Function    # forward
    fop_T::Function  # transpose
    fop_A::Function  # adjoint
    fop_C::Function  # conj
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
end

# type exception
struct joLooseMatrixException <: Exception
    msg :: String
end

############################################################
# joLooseLinearFunction ####################################

export joLooseLinearFunction, joLooseLinearFunctionException

# type definition
"""
joLooseLinearFunction type

# TYPE PARAMETERS
- DDT::DataType : domain DataType
- RDT::DataType : range DataType

# FIELDS
- name::String : given name
- m::Integer : # of rows
- n::Integer : # of columns
- fop::Function : forward function
- fop_T::Nullable{Function} : transpose function
- fop_A::Nullable{Function} : adjoint function
- fop_C::Nullable{Function} : conj function
- fMVok : whether fops are rady to handle mvec
- iop::Nullable{Function} : inverse for fop
- iop_T::Nullable{Function} : inverse for fop_T
- iop_A::Nullable{Function} : inverse for fop_A
- iop_C::Nullable{Function} : inverse for fop_C
- iMVok::Bool : whether iops are rady to handle mvec

"""
struct joLooseLinearFunction{DDT<:Number,RDT<:Number} <: joAbstractFosterLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    fop::Function              # forward
    fop_T::Nullable{Function}  # transpose
    fop_A::Nullable{Function}  # adjoint
    fop_C::Nullable{Function}  # conj
    fMVok::Bool                # forward can do mvec
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
    iMVok::Bool
end

# type exceptions
struct joLooseLinearFunctionException <: Exception
    msg :: String
end


text:
# Utilities for working with HSA agents

mutable struct HSAAgent
    agent::HSA.Agent
end

const DEFAULT_AGENT = Ref{HSAAgent}()

### @cfunction callbacks ###

function agent_iterate_isas_cb(isa::HSA.ISA, ptr::Ptr{HSA.ISA})
    unsafe_store!(ptr, isa)
    return HSA.STATUS_INFO_BREAK
end

function iterate_agents_cb(agent::HSA.Agent, agents)
    push!(agents, HSAAgent(agent))
    return HSA.STATUS_SUCCESS
end

"Determines if a memory region can be used for kernarg allocations."
function get_kernarg_memory_region_cb(region::HSA.Region, data::Ptr{HSA.Region})
    segment = Ref{HSA.RegionSegment}()
    HSA.region_get_info(region, HSA.REGION_INFO_SEGMENT, segment)
    if (segment[] != HSA.REGION_SEGMENT_GLOBAL)
        return HSA.STATUS_SUCCESS
    end

    flags = Ref{HSA.RegionGlobalFlag}()
    HSA.region_get_info(region, HSA.REGION_INFO_GLOBAL_FLAGS, flags)
    if (flags[] & HSA.REGION_GLOBAL_FLAG_KERNARG > 0)
        unsafe_store!(data, region)
        return HSA.STATUS_INFO_BREAK
    end

    return HSA.STATUS_SUCCESS
end

"Determines if a memory region can be used for fine grained allocations."
function get_fine_grained_memory_region_cb(region::HSA.Region, data::Ptr{HSA.Region})
    segment = Ref{HSA.RegionSegment}()
    HSA.region_get_info(region, HSA.REGION_INFO_SEGMENT, segment)
    if (segment[] != HSA.REGION_SEGMENT_GLOBAL)
        return HSA.STATUS_SUCCESS
    end

    flags = Ref{HSA.RegionGlobalFlag}()
    HSA.region_get_info(region, HSA.REGION_INFO_GLOBAL_FLAGS, flags)
    if (flags[] & HSA.REGION_GLOBAL_FLAG_FINE_GRAINED > 0)
        unsafe_store!(data, region)
        return HSA.STATUS_INFO_BREAK
    end

    return HSA.STATUS_SUCCESS
end

function Base.show(io::IO, agent::HSAAgent)
    print(io, "HSAAgent($(agent.agent)): Name=$(get_name(agent)), Type=$(device_type(agent))")
end

function get_agents()
    agents = Ref(Vector{HSAAgent}())
    GC.@preserve agents begin
        func = @cfunction(iterate_agents_cb, HSA.Status,
                (HSA.Agent, Ref{Vector{HSAAgent}}))
        HSA.iterate_agents(func, agents) |> check
        _agents = agents[]
    end
    return _agents
end
get_agents(kind::Symbol) =
    filter(agent->device_type(agent)==kind, get_agents())

get_default_agent() = DEFAULT_AGENT[]
function set_default_agent!(kind::Symbol)
    DEFAULT_AGENT[] = first(get_agents(kind))
end
set_default_agent!() = set_default_agent!(:gpu)

function get_name(agent::HSAAgent)
    #len = Ref(0)
    #hsa_agent_get_info(agent.agent, HSA.AGENT_INFO_NAME_LENGTH, len) |> check
    #FIXME: this could be better?
    name = Vector{UInt8}(undef, 64)
    getinfo(agent.agent, HSA.AGENT_INFO_NAME, name) |> check
    return rstrip(String(name), '\0')
end
function profile(agent::HSAAgent)
    profile = Ref{HSA.Profile}()
    getinfo(agent.agent, HSA.AGENT_INFO_PROFILE, profile) |> check
    return profile[]
end

function device_type(agent::HSAAgent)
    devtype = Ref{UInt32}()
    getinfo(agent.agent, HSA.AGENT_INFO_DEVICE, devtype) |> check

    if HSA.DeviceType(devtype[]) == HSA.DEVICE_TYPE_CPU
        return :cpu
    elseif HSA.DeviceType(devtype[]) == HSA.DEVICE_TYPE_GPU
        return :gpu
    elseif HSA.DeviceType(devtype[]) == HSA.DEVICE_TYPE_DSP
        return :dsp
    end

end

# TODO: Get all ISAs
function get_first_isa(agent::HSAAgent)
    isa = Ref{HSA.ISA}()
    func = @cfunction(agent_iterate_isas_cb, HSA.Status,
        (HSA.ISA, Ptr{HSA.ISA}))
    ret = HSA.agent_iterate_isas(agent.agent, func, isa)
    @assert ret == HSA.STATUS_INFO_BREAK "Failed to find an agent ISA"

    len = Ref{Cuint}(0)
    getinfo(isa[], HSA.ISA_INFO_NAME_LENGTH, len) |> check
    name = repeat(" ", len[])
    getinfo(isa[], HSA.ISA_INFO_NAME, name) |> check
    isa_name = string(rstrip(last(split(name, "-")), '\0'))
    return isa_name
end

function get_region(agent::HSAAgent, kind::Symbol)
    region = newref!(Ref{HSA.Region}, typemax(UInt64))
    if kind == :finegrained
        func = @cfunction(get_fine_grained_memory_region_cb, HSA.Status,
                (HSA.Region, Ptr{HSA.Region}))
    elseif kind == :kernarg
        func = @cfunction(get_kernarg_memory_region_cb, HSA.Status,
                (HSA.Region, Ptr{HSA.Region}))
    else
        error("Region kind $kind not supported")
    end
    HSA.agent_iterate_regions(agent.agent, func, region)
    check((region[].handle == typemax(UInt64)) ? HSA.STATUS_ERROR : HSA.STATUS_SUCCESS)
    return region
end


text:
@testset "DQMC" begin
    m = HubbardModelAttractive(dims=1, L=8);

    # constructors
    dqmc = DQMC(m; beta=5.0)
    @test m.L == 8 && m.dims == 1
    @test typeof(m) == MonteCarlo.HubbardModelAttractive{MonteCarlo.Chain}
    d = Dict(:dims=>2,:L=>3)
    m = HubbardModelAttractive(d)
    @test typeof(m) == MonteCarlo.HubbardModelAttractive{MonteCarlo.SquareLattice}
    @test m.L == 3 && m.dims == 2

    # generic checkerboard
    sq = MonteCarlo.SquareLattice(4);
    @test MonteCarlo.build_checkerboard(sq) == ([1.0 3.0 5.0 7.0 9.0 11.0 13.0 15.0 1.0 2.0 4.0 6.0 9.0 10.0 12.0 14.0 2.0 3.0 4.0 5.0 8.0 10.0 11.0 16.0 6.0 7.0 8.0 12.0 13.0 14.0 15.0 16.0; 2.0 4.0 6.0 8.0 10.0 12.0 14.0 16.0 5.0 3.0 8.0 7.0 13.0 11.0 16.0 15.0 6.0 7.0 1.0 9.0 12.0 14.0 15.0 13.0 10.0 11.0 5.0 9.0 1.0 2.0 3.0 4.0; 1.0 5.0 9.0 13.0 17.0 21.0 25.0 29.0 2.0 3.0 8.0 11.0 18.0 19.0 24.0 27.0 4.0 6.0 7.0 10.0 16.0 20.0 22.0 31.0 12.0 14.0 15.0 23.0 26.0 28.0 30.0 32.0], UnitRange[1:8, 9:16, 17:24, 25:32], 4)

    m = HubbardModelAttractive(dims=2, L=8, mu=0.5)
    mc1 = DQMC(m, beta=5.0)
    mc2 = DQMC(m, beta=5.0, checkerboard=false)
    mc2.conf = deepcopy(mc1.conf)
    MonteCarlo.init_hopping_matrices(mc1, m)
    MonteCarlo.init_hopping_matrices(mc2, m)
    MonteCarlo.build_stack(mc1)
    MonteCarlo.build_stack(mc2)
    @test MonteCarlo.slice_matrix(mc1, m, 1, 1.) == MonteCarlo.slice_matrix(mc2, m, 1, 1.)

    mc = DQMC(m, beta=5.0, checkerboard=true, delta_tau=0.1)
    MonteCarlo.init_hopping_matrices(mc, m)
    hop_mat_exp_chkr = foldl(*,mc.s.chkr_hop_half) * sqrt.(mc.s.chkr_mu)
    r = MonteCarlo.effreldiff(mc.s.hopping_matrix_exp,hop_mat_exp_chkr)
    r[findall(x -> x==zero(x), hop_mat_exp_chkr)] .= 0.
    @test maximum(MonteCarlo.absdiff(mc.s.hopping_matrix_exp,hop_mat_exp_chkr)) <= mc.p.delta_tau

    # initial greens test
    mc = DQMC(m, beta=5.0, safe_mult=1)
    MonteCarlo.build_stack(mc)
    MonteCarlo.propagate(mc)
    # With this we effectively test calculate_greens without wrap_greens
    greens, = MonteCarlo.calculate_greens_and_logdet(mc, mc.s.current_slice+1)
    MonteCarlo.wrap_greens!(mc, greens, mc.s.current_slice+1, -1)
    @test greens ≈ mc.s.greens
    # here with a single implied wrap
    greens, = MonteCarlo.calculate_greens_and_logdet(mc, mc.s.current_slice)
    @test maximum(MonteCarlo.absdiff(greens, mc.s.greens)) < 1e-13

    # wrap greens test
    for k in 0:9
        MonteCarlo.wrap_greens!(mc, mc.s.greens, mc.s.current_slice - k, -1)
    end
    greens, = MonteCarlo.calculate_greens_and_logdet(mc, mc.s.current_slice-10)
    @test maximum(MonteCarlo.absdiff(greens, mc.s.greens)) < 1e-9

    include("slice_matrices.jl")
end


text:
# Design optimized routines to compute pressure and its Jacobian if the vortices are placed according to the method of images, and we evaluate the pressure field about the x axis.

export measure_state_symmetric, symmetric_pressure#, analytical_symmetric_jacobian_strength

function measure_state_symmetric(state, t, config::VortexConfig; withfreestream::Bool = false)
	if withfreestream == false
    	return symmetric_pressure(real(config.ss), state_to_lagrange(state, config)[1], t)
	else
		freestream = Freestream(config.U)
		return symmetric_pressure(real(config.ss), state_to_lagrange(state, config)[1], freestream, t)
	end
end

# Symmetric pressure calculation for point vortices
function symmetric_pressure(target::Vector{Float64}, source::T, t) where T <: Vector{PotentialFlow.Points.Point{Float64, Float64}}

    Nv = size(source, 1)
    Ny = size(target, 1)

    press = zeros(Ny)

    # Quadratic term

    @inbounds for (J, bJ) in enumerate(source)
        tmpJ = bJ.S*imag(bJ.z)
        for (i, xi) in enumerate(target)
            diJ = inv(abs2(xi-bJ.z))
            press[i] += tmpJ*diJ
        end
    end
    press .= deepcopy((-0.5/π^2)*press.^2)

    sourcevels = zeros(ComplexF64, Nv)
    @inbounds for (J, bJ) in enumerate(source)
        zJ = bJ.z
        ΓJ = bJ.S
        for (K, bK) in enumerate(source)
            if K != J
                zK = bK.z
                ΓK = bK.S
                sourcevels[J] += (ΓK*imag(zK))/(π*(zJ - zK)*(zJ - conj(zK)))
            end
        end
        # Contribution of the mirrored vortex
        sourcevels[J] += ΓJ/(4*π*imag(zJ))
    end

    # Unsteady term
    @inbounds for (J, bJ) in enumerate(source)
        zJ = bJ.z
        ΓJ = bJ.S
        for (i, xi) in enumerate(target)
            press[i] += 2*real((-im*ΓJ)*inv(2*π*(xi - zJ))*conj(sourcevels[J]))
        end
    end

    return press
end

# Symmetric pressure calculation for point vortices with freestream
function symmetric_pressure(target::Vector{Float64}, source::T, freestream, t) where T <: Vector{PotentialFlow.Points.Point{Float64, Float64}}

    Nv = size(source, 1)
    Ny = size(target, 1)

	U = freestream.U

    press = zeros(Ny)

    # Quadratic term
    @inbounds for (J, bJ) in enumerate(source)
        tmpJ = bJ.S*imag(bJ.z)
        for (i, xi) in enumerate(target)
            diJ = inv(abs2(xi-bJ.z))
            press[i] += tmpJ*diJ
        end
    end

	press .=  -0.5*deepcopy(abs2.(1/π*press .+ conj(U)))

	sourcevels = zeros(ComplexF64, Nv)
    @inbounds for (J, bJ) in enumerate(source)
        zJ = bJ.z
        ΓJ = bJ.S
        for (K, bK) in enumerate(source)
            if K != J
                zK = bK.z
                ΓK = bK.S
                sourcevels[J] += (ΓK*imag(zK))/(π*(zJ - zK)*(zJ - conj(zK)))
            end
        end
        # Contribution of the mirrored vortex
        sourcevels[J] += ΓJ/(4*π*imag(zJ)) + conj(U)
    end

    # Unsteady term
    @inbounds for (J, bJ) in enumerate(source)
        zJ = bJ.z
        ΓJ = bJ.S
        for (i, xi) in enumerate(target)
            press[i] += 2*real((-im*ΓJ)*inv(2*π*(xi - zJ))*conj(sourcevels[J]))
        end
    end

    return press
end

# Symmetric pressure calculation for regularized vortices
function symmetric_pressure(target::Vector{Float64}, source::T, t) where T <: Vector{PotentialFlow.Blobs.Blob{Float64, Float64}}

	source = deepcopy(source)

    Nv = size(source, 1)
    Ny = size(target, 1)

	δ = source[1].δ

    press = zeros(Ny)

    # Quadratic term
    @inbounds for (J, bJ) in enumerate(source)
		zJ = bJ.z
        tmpJ = bJ.S*imag(zJ)
        for (i, xi) in enumerate(target)
            diJ = abs2(xi - zJ) + δ^2
            press[i] += tmpJ/diJ
		end
	end

	press .= deepcopy((-0.5/π^2)*press.^2)

    sourcevels = zeros(ComplexF64, Nv)
    @inbounds for (J, bJ) in enumerate(source)
        zJ = bJ.z
		yJ = imag(zJ)
        ΓJ = bJ.S
        for (K, bK) in enumerate(source)
            if K != J
                zK = bK.z
                ΓK = bK.S
                sourcevels[J] += ΓK*(conj(zJ - zK))/(abs2(zJ - zK) + δ^2)
				sourcevels[J] += -ΓK*(conj(zJ) - zK)/(abs2(zJ - conj(zK)) + δ^2)
            end
        end
		sourcevels[J] *= -im/(2*π)
        # Contribution of the mirrored vortex
		sourcevels[J] += (ΓJ*yJ/π)/(4*(yJ)^2 + δ^2)
    end

    # Unsteady term
    @inbounds for (J, bJ) in enumerate(source)
        zJ = bJ.z
        ΓJ = bJ.S
        for (i, xi) in enumerate(target)
            press[i] -= ΓJ/π*imag(sourcevels[J]/(xi - conj(zJ)))
        end
    end
    return press
end

# Symmetric pressure calculation for regularized vortices with freestream
function symmetric_pressure(target::Vector{Float64}, source::T, freestream, t) where T <: Vector{PotentialFlow.Blobs.Blob{Float64, Float64}}

	source = deepcopy(source)

    Nv = size(source, 1)
    Ny = size(target, 1)

	δ = source[1].δ
	U = freestream.U

    press = zeros(Ny)

    # Quadratic term
    @inbounds for (J, bJ) in enumerate(source)
		zJ = bJ.z
        tmpJ = bJ.S*imag(zJ)
        for (i, xi) in enumerate(target)
            diJ = abs2(xi - zJ) + δ^2
            press[i] += tmpJ/diJ
		end
	end

	press .=  -0.5*deepcopy(abs2.(1/π*press .+ conj(U)))

    sourcevels = zeros(ComplexF64, Nv)
    @inbounds for (J, bJ) in enumerate(source)
        zJ = bJ.z
		yJ = imag(zJ)
        ΓJ = bJ.S
        for (K, bK) in enumerate(source)
            if K != J
                zK = bK.z
                ΓK = bK.S
                sourcevels[J] += ΓK*(conj(zJ - zK))/(abs2(zJ - zK) + δ^2)
				sourcevels[J] += -ΓK*(conj(zJ) - zK)/(abs2(zJ - conj(zK)) + δ^2)
            end
        end
		sourcevels[J] *= -im/(2*π)
        # Contribution of the mirrored vortex
		sourcevels[J] += (ΓJ*yJ/π)/(4*(yJ)^2 + δ^2) + conj(U)
    end

    # Unsteady term
    @inbounds for (J, bJ) in enumerate(source)
        zJ = bJ.z
        ΓJ = bJ.S
        for (i, xi) in enumerate(target)
            press[i] -= ΓJ/π*imag(sourcevels[J]/(xi - conj(zJ)))
        end
    end
    return press
end


text:
const VIO_PIX_FMT_DEF_ELTYPE_LU = Dict{Cint,DataType}(
    AV_PIX_FMT_GRAY8 => Gray{N0f8},
    AV_PIX_FMT_GRAY10LE => Gray{N6f10},
    AV_PIX_FMT_RGB24 => RGB{N0f8},
    AV_PIX_FMT_GBRP10LE => RGB{N6f10},
)

const VIO_SUPPORTED_PIX_FMTS = collect(keys(VIO_PIX_FMT_DEF_ELTYPE_LU))
# This array is terminated by AV_PIX_FMT_NONE so it can be used by libav functions
const VIO_SUPPORTED_PIX_FMTS_AVARRAY = push!(copy(VIO_SUPPORTED_PIX_FMTS), AV_PIX_FMT_NONE)

strip_interpretation(::Type{X}) where {T,X<:Gray{T}} = strip_interpretation(T)
strip_interpretation(::Type{X}) where {T,X<:Normed{T}} = strip_interpretation(T)
strip_interpretation(::Type{T}) where {T} = T

const VIO_DEF_ELTYPE_PIX_FMT_LU =
    Dict{DataType,Int32}((strip_interpretation(v) => k for (k, v) in pairs(VIO_PIX_FMT_DEF_ELTYPE_LU)))

struct VioColorspaceDetails
    color_primaries::Cint
    color_trc::Cint
    colorspace::Cint
    color_range::Cint
end

VioColorspaceDetails() =
    VioColorspaceDetails(AVCOL_PRI_UNSPECIFIED, AVCOL_TRC_UNSPECIFIED, AVCOL_SPC_UNSPECIFIED, AVCOL_RANGE_UNSPECIFIED)

const VIO_DEFAULT_COLORSPACE_DETAILS =
    VioColorspaceDetails(AVCOL_PRI_BT709, AVCOL_TRC_IEC61966_2_1, AVCOL_SPC_RGB, AVCOL_RANGE_JPEG)

const VIO_DEFAULT_TRANSFER_COLORSPACE_DETAILS = Dict{Cint,VioColorspaceDetails}(
    AV_PIX_FMT_GRAY8 => VIO_DEFAULT_COLORSPACE_DETAILS,
    AV_PIX_FMT_GRAY10LE => VIO_DEFAULT_COLORSPACE_DETAILS,
    AV_PIX_FMT_RGB24 => VIO_DEFAULT_COLORSPACE_DETAILS,
    AV_PIX_FMT_GBRP10LE => VIO_DEFAULT_COLORSPACE_DETAILS,
)

# avarray_dst_pix_fmts MUST be terminated properly, see `avcodec_find_best_pix_fmt_of_list`
function _vio_determine_best_pix_fmt(src_pix_fmt, avarray_dst_pix_fmts = VIO_SUPPORTED_PIX_FMTS_AVARRAY; loss_flags = 0)
    loss_ptr = Ref{Cint}(loss_flags)
    dst_pix_fmt = avcodec_find_best_pix_fmt_of_list(avarray_dst_pix_fmts, src_pix_fmt, 0, loss_ptr)
    return dst_pix_fmt, loss_ptr[]
end

is_pixel_type_supported(pxfmt) = haskey(VIO_PIX_FMT_DEF_ELTYPE_LU, pxfmt)

is_eltype_transfer_supported(::Type{T}) where {T} = haskey(VIO_DEF_ELTYPE_PIX_FMT_LU, strip_interpretation(T))

get_transfer_pix_fmt(::Type{T}) where {T} = VIO_DEF_ELTYPE_PIX_FMT_LU[strip_interpretation(T)]

@inline function bytes_of_uint16(x::UInt16)
    msb = convert(UInt8, (x & 0xFF00) >> 8)
    lsb = convert(UInt8, x & 0x00FF)
    return msb, lsb
end

@inline function unsafe_store_uint16_in_le!(data, x, offset)
    msb, lsb = bytes_of_uint16(x)
    unsafe_store!(data, lsb, offset)
    return unsafe_store!(data, msb, offset + 1)
end

@inline uint16_from_bytes(msb, lsb) = (convert(UInt16, msb) << 8) | lsb

function load_uint16_from_le_bytes(data, pos)
    lsb = unsafe_load(data, pos)
    msb = unsafe_load(data, pos + 1)
    return uint16_from_bytes(msb, lsb)
end

load_n6f10_from_le_bytes(data, pos) = reinterpret(N6f10, load_uint16_from_le_bytes(data, pos))

function set_basic_frame_properties!(frame, width, height, format)
    frame.format = format
    frame.width = width
    frame.height = height
    ret = av_frame_get_buffer(frame, 0) # Allocate picture buffers
    return ret < 0 && error("Could not allocate the video frame data")
end

#############################
## AVFrame -> Julia Buffer ##
#############################

# Transfer bytes from AVFrame to buffer
function _transfer_frame_bytes_to_img_buf!(buf::AbstractArray{UInt8}, frame, bytes_per_pixel)
    width = frame.width
    height = frame.height
    ip = frame.data[1]
    op = pointer(buf)
    out_linesize = width * bytes_per_pixel
    for h in 1:height
        out_line_p = op + (h - 1) * out_linesize
        in_line_p = ip + (h - 1) * frame.linesize[1]
        unsafe_copyto!(out_line_p, in_line_p, out_linesize)
    end
    return buf
end

# Read a 8 bit monochrome frame
function transfer_frame_to_img_buf!(buf::AbstractArray{UInt8}, frame, bytes_per_pixel)
    target_format = frame.format
    if target_format == AV_PIX_FMT_GRAY8
        _transfer_frame_bytes_to_img_buf!(buf, frame, bytes_per_pixel)
    else
        unsupported_retrieval_format(target_format)
    end
    return buf
end

# Read a 10 bit monochrome frame
function transfer_frame_to_img_buf!(buf::AbstractArray{UInt16}, frame, bytes_per_pixel)
    target_format = frame.format
    if target_format == AV_PIX_FMT_GRAY10LE
        if ENDIAN_BOM != 0x04030201
            error("Reading AV_PIX_FMT_GRAY10LE on big-endian machines not yet supported")
        end
        _transfer_frame_bytes_to_img_buf!(reinterpret(UInt8, buf), frame, bytes_per_pixel)
    else
        unsupported_retrieval_format(target_format)
    end
    return buf
end

transfer_frame_to_img_buf!(buf::AbstractArray{X}, args...) where {T,X<:Normed{T}} =
    transfer_frame_to_img_buf!(reinterpret(T, buf), args...)

transfer_frame_to_img_buf!(buf::AbstractArray{<:Gray}, args...) = transfer_frame_to_img_buf!(channelview(buf), args...)

# Read a 8 bit RGB frame
function transfer_frame_to_img_buf!(buf::AbstractArray{RGB{N0f8}}, frame, bytes_per_pixel)
    target_format = frame.format
    if target_format == AV_PIX_FMT_RGB24
        _transfer_frame_bytes_to_img_buf!(reinterpret(UInt8, buf), frame, bytes_per_pixel)
    else
        unsupported_retrieval_format(target_format)
    end
    return buf
end

# Read a 10 bit RGB frame
function transfer_frame_to_img_buf!(buf::AbstractArray{RGB{N6f10}}, frame, ::Any)
    target_format = frame.format
    if target_format == AV_PIX_FMT_GBRP10LE
        bytes_per_sample = 2
        width = frame.width
        height = frame.height
        size(buf) == (width, height) || error("buffer wrong size")
        ls = frame.linesize
        @inbounds linesizes = ntuple(i -> ls[i], 3)
        @inbounds for r in 1:height
            line_offsets = (r - 1) .* linesizes
            @simd ivdep for c in 1:width
                col_pos = (c - 1) * bytes_per_sample + 1
                line_poss = line_offsets .+ col_pos
                rg = load_n6f10_from_le_bytes(frame.data[1], line_poss[1])
                rb = load_n6f10_from_le_bytes(frame.data[2], line_poss[2])
                rr = load_n6f10_from_le_bytes(frame.data[3], line_poss[3])
                buf[c, r] = RGB{N6f10}(rr, rg, rb) # scanline-major!
            end
        end
    else
        unsupported_retrieval_format(target_format)
    end
    return buf
end

transfer_frame_to_img_buf!(buf::PermutedDimsArray, frame, bytes_per_sample) =
    transfer_frame_to_img_buf!(parent(buf), frame, bytes_per_sample)

#############################
## Julia Buffer -> AVFrame ##
#############################

_unsupported_append_encode_type() = error("Array element type not supported")

# bytes_per_sample is the number of bytes per pixel sample, not the size of the
# element type of img
function transfer_img_bytes_to_frame_plane!(
    data_ptr,
    img::AbstractArray{UInt8},
    px_width,
    px_height,
    data_linesize,
    bytes_per_sample = 1,
)
    stride(img, 1) == 1 || error("stride(img, 1) must be equal to one")
    img_line_nbytes = px_width * bytes_per_sample
    if data_linesize == img_line_nbytes
        # When both sizes match the buffers are both contiguous, so can be transferred in one go
        unsafe_copyto!(data_ptr, pointer(img), length(img))
    else
        @inbounds for r in 1:px_height
            data_line_ptr = data_ptr + (r - 1) * data_linesize
            img_line_ptr = pointer(img, img_line_nbytes * (r - 1) + 1)
            unsafe_copyto!(data_line_ptr, img_line_ptr, img_line_nbytes)
        end
    end
end

function make_into_sl_col_mat(img::AbstractVector{<:Union{RGB,Unsigned}}, width, height)
    img_mat = reshape(img, (height, width))
    return PermutedDimsArray(img_mat, (2, 1))
end
make_into_sl_col_mat(img::AbstractMatrix{<:Union{RGB,Unsigned}}, args...) = PermutedDimsArray(img, (2, 1))

function transfer_sl_col_img_buf_to_frame!(frame, img::AbstractArray{UInt8})
    frame.format == AV_PIX_FMT_RGB24 || _unsupported_append_encode_type()
    width = frame.width
    height = frame.height
    ls = frame.linesize[1]
    data_p = frame.data[1]
    @inbounds for r in 1:height
        line_offset = (r - 1) * ls
        @simd for c in 1:width
            val = img[c, r]
            row_offset = line_offset + c
            for s in 0:2
                @preserve frame unsafe_store!(data_p, val, row_offset + s)
            end
        end
    end
end

function transfer_img_buf_to_frame!(frame, img::AbstractArray{UInt8}, scanline_major)
    pix_fmt = frame.format
    if pix_fmt == AV_PIX_FMT_GRAY8
        width = frame.width
        height = frame.height
        ls = frame.linesize[1]
        data_p = frame.data[1]
        if scanline_major
            @preserve frame transfer_img_bytes_to_frame_plane!(data_p, img, width, height, ls)
        else
            @inbounds for r in 1:height
                line_offset = (r - 1) * ls
                @simd for c in 1:width
                    @preserve frame unsafe_store!(data_p, img[r, c], line_offset + c)
                end
            end
        end
    elseif pix_fmt == AV_PIX_FMT_RGB24
        # This is specifying each RGB triple with one UInt8 "luma" value
        if scanline_major
            transfer_sl_col_img_buf_to_frame!(frame, img)
        else
            img_sl_col_mat = make_into_sl_col_mat(img)
            transfer_sl_col_img_buf_to_frame!(frame, img_sl_col_mat)
        end
    else
        _unsupported_append_encode_type()
    end
end

function transfer_img_buf_to_frame!(frame, img::AbstractArray{UInt16}, scanline_major)
    bytes_per_sample = 2
    pix_fmt = frame.format
    if pix_fmt == AV_PIX_FMT_GRAY10LE
        width = frame.width
        height = frame.height
        ls = frame.linesize[1]
        datap = frame.data[1]
        if scanline_major
            if ENDIAN_BOM != 0x04030201
                error("""
    Writing scanline_major AV_PIX_FMT_GRAY10LE on big-endian machines not yet
    supported, use scanline_major = false
    """)
            end
            img_raw = reinterpret(UInt8, img)
            @preserve frame transfer_img_bytes_to_frame_plane!(datap, img_raw, width, height, ls, bytes_per_sample)
        else
            @inbounds for r in 1:height
                line_offset = (r - 1) * ls
                @simd for c in 1:width
                    px_offset = line_offset + bytes_per_sample * (c - 1) + 1
                    @preserve frame unsafe_store_uint16_in_le!(datap, img[r, c], px_offset)
                end
            end
        end
    else
        _unsupported_append_encode_type()
    end
end

simple_rawview(a::AbstractArray{X}) where {T,X<:Normed{T}} = reinterpret(T, a)

transfer_img_buf_to_frame!(frame, img::AbstractArray{<:Normed}, args...) =
    transfer_img_buf_to_frame!(frame, simple_rawview(img), args...)

transfer_img_buf_to_frame!(frame, img::AbstractArray{<:Gray}, args...) =
    transfer_img_buf_to_frame!(frame, channelview(img), args...)

function transfer_img_buf_to_frame!(frame, img::AbstractMatrix{RGB{N0f8}}, scanline_major)
    pix_fmt = frame.format
    if pix_fmt == AV_PIX_FMT_RGB24
        width = frame.width
        height = frame.height
        lss = frame.linesize
        fdata = frame.data
        bytes_per_pixel = 3
        if scanline_major
            data_p = fdata[1]
            @preserve frame transfer_img_bytes_to_frame_plane!(
                data_p,
                reinterpret(UInt8, img),
                width,
                height,
                lss[1],
                bytes_per_pixel,
            )
        else
            @inbounds for h in 1:height
                line_offset = (h - 1) * lss[1]
                @simd for w in 1:width
                    px = img[h, w]
                    px_offset = line_offset + (w - 1) * bytes_per_pixel + 1
                    unsafe_store!(fdata[1], reinterpret(px.r), px_offset)
                    unsafe_store!(fdata[1], reinterpret(px.g), px_offset + 1)
                    unsafe_store!(fdata[1], reinterpret(px.b), px_offset + 2)
                end
            end
        end
    else
        _unsupported_append_encode_type()
    end
end

function transfer_sl_col_img_buf_to_frame!(frame, img::AbstractArray{RGB{N6f10}})
    pix_fmt = frame.format
    pix_fmt == AV_PIX_FMT_GBRP10LE || _unsupported_append_encode_type()
    bytes_per_sample = 2
    # Luma
    fdata = frame.data
    linesize_y = frame.linesize[1]
    width = frame.width
    height = frame.height
    all(size(img) .>= (width, height)) || error("img is not large enough")
    ls = frame.linesize
    @inbounds linesizes = ntuple(i -> ls[i], 3)
    @inbounds for h in 1:height
        line_offsets = (h - 1) .* linesizes
        @simd ivdep for w in 1:width
            col_pos = (w - 1) * bytes_per_sample + 1
            line_poss = line_offsets .+ col_pos
            input_px = img[w, h]
            unsafe_store_uint16_in_le!(fdata[1], reinterpret(green(input_px)), line_poss[1])
            unsafe_store_uint16_in_le!(fdata[2], reinterpret(blue(input_px)), line_poss[2])
            unsafe_store_uint16_in_le!(fdata[3], reinterpret(red(input_px)), line_poss[3])
        end
    end
end

function transfer_img_buf_to_frame!(frame, img::AbstractArray{RGB{N6f10}}, scanline_major)
    if scanline_major
        transfer_img = img
    else
        transfer_img = img_sl_col_mat = make_into_sl_col_mat(img)
    end
    return transfer_sl_col_img_buf_to_frame!(frame, transfer_img)
end

# Fallback
transfer_img_buf_to_frame!(frame, img, scanline_major) = _unsuported_append_encode_type()


text:
include("./src/AssociationRules.jl")

v = [sample(1:20, 5, replace = false) for x = 1:10_000];

@time res = par_apriori_gen(v[1:100], 4);
# Profile.clear_malloc_data()

@time res = par_apriori_gen(v, 4);

@time res = par_apriori_gen3(v[1:100], 4, 100);
# Profile.clear_malloc_data()

@time res = par_apriori_gen3(v, 4, 10_000);