text:
# TODO Add VPRKpSecondary tableau.

"Parameters for right-hand side function of variational partitioned Runge-Kutta methods."
mutable struct ParametersVPRKpSecondary{DT, TT, ET <: VODE{DT,TT}, D, S} <: AbstractParametersVPRK{DT,TT,ET,D,S}
    equ::ET
    tab::TableauVPRK{TT}
    Δt::TT

    ω::Matrix{TT}
    R::Vector{TT}

    t̅::TT
    q̅::Vector{DT}
    p̅::Vector{DT}
end

function ParametersVPRKpSecondary(equ::ET, tab::TableauVPRK{TT}, Δt::TT) where {DT, TT, ET <: VODE{DT,TT}}
    D = equ.d
    S = tab.s

    # compute reduction matrix
    ω = zeros(TT, S-1, S)

    for i in 1:(S-1)
        for j in 1:S
            ω[i,j] = tab.q.b[j] * tab.q.c[j]^(i-1)
        end
    end

    R = convert(Vector{TT}, [1, tab.R∞])
    q = zeros(DT,D)
    p = zeros(DT,D)

    ParametersVPRKpSecondary{DT,TT,ET,D,S}(equ, tab, Δt, ω, R, 0, q, p)
end

@doc raw"""
Variational partitioned Runge-Kutta integrator with projection on secondary constraint.

The VPRK integrator solves the following system of equations for the internal stages,
```math
\begin{align*}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, \big( V_{n,j} + \Lambda_{n,j} \big) , \\
P_{n,i} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{a}_{ij} \, \big( F_{n,j} + \nabla \vartheta (Q_{n,j})^{T} \cdot \Lambda_{n,j} \big) - d_i \lambda , \\
0 &= \sum \limits_{i=1}^{s} d_i V_i , \\
0 &= \sum \limits_{j=1}^{s} \omega_{ij} \Psi_{n,j} ,
\end{align*}
```
with definitions
```math
\begin{align*}
P_{n,i} &= \dfrac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , \\
F_{n,i} &= \dfrac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , \\
\Psi_{n,i} &= \psi(Q_{n,i}, V_{n,i}, P_{n,i}, F_{n,i}) ,
\end{align*}
```
and update rule
```math
\begin{align*}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, \big( V_{n,i} + \Lambda_{n,i} \big) , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i} \, \big( F_{n,i} + \nabla \vartheta (Q_{n,j})^{T} \cdot \Lambda_{n,j} \big) , \\
0 &= \phi (q_{n+1}, p_{n+1}) ,
\end{align*}
```
satisfying the symplecticity conditions
```math
\begin{align*}
b_{i} \bar{a}_{ij} + \bar{b}_{j} a_{ji} &= b_{i} \bar{b}_{j} , &
\bar{b}_i &= b_i ,
\end{align*}
```
the primary constraint,
```math
\begin{align*}
\phi(q,p) = p - \vartheta (q) = 0 ,
\end{align*}
```
at the final solution ``(q_{n+1}, p_{n+1})``,
and super positions of the secondary constraints,
```math
\begin{align*}
\psi(q,\dot{q},p,\dot{p})
= \dot{p} - \dot{q} \cdot \nabla \vartheta (q)
= \big( \nabla \vartheta (q) - \nabla \vartheta^{T} (q) \big) \cdot \dot{q} - \nabla H (q)
= 0,
\end{align*}
```
which, evaluated at the internal stages, read
```math
\begin{align*}
\Psi_{n,j} = \big( \nabla \vartheta (Q_{n,j}) - \nabla \vartheta^{T} (Q_{n,j}) \big) \cdot V_{n,j} - \nabla H (Q_{n,j}) .
\end{align*}
```
Here, ``\omega`` is a ``(s-1) \times s`` matrix, chosen such that the resulting
method has optimal order.
The vector ``d`` is zero for Gauss-Legendre methods and needs to be chosen
appropriately for Gauss-Lobatto methods (for details see documentation of
VPRK methods).
"""
struct IntegratorVPRKpSecondary{DT, TT, PT <: ParametersVPRKpSecondary{DT,TT},
                                        ST <: NonlinearSolver{DT},
                                        IT <: InitialGuessPODE{DT,TT}, D, S} <: AbstractIntegratorVPRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    cache::IntegratorCacheVPRK{DT,D,S}
end

function IntegratorVPRKpSecondary(equation::ET, tableau::TableauVPRK{TT}, Δt::TT) where {DT, TT, ET <: VODE{DT,TT}}
    D = equation.d
    M = equation.n
    S = tableau.s

    # create params
    params = ParametersVPRKpSecondary(equation, tableau, Δt)

    # create solver
    solver = create_nonlinear_solver(DT, 2*D*S, params)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create cache for internal stage vectors and update vectors
    cache = IntegratorCacheVPRK{DT,D,S}(true)

    # create integrator
    IntegratorVPRKpSecondary{DT, TT, typeof(params), typeof(solver), typeof(iguess), D, S}(
                params, solver, iguess, cache)
end

function compute_stages_vprk!(x, q, p, Q, V, Λ, P, F, R, Φ, params)
    # copy x to V
    compute_stages_v_vprk!(x, V, params)

    # copy x to Λ
    compute_stages_λ_vprk!(x, Λ, params)

    # compute Q
    compute_stages_q_vprk!(q, Q, V, Λ, params)

    # compute p̅, R and Ψ
    compute_projection_vprk!(q, p, Q, V, Λ, R, Φ, params)

    # compute P and F
    compute_stages_p_vprk!(Q, V, P, F, params)
end

function compute_stages_q_vprk!(q::Vector{ST}, Q::Vector{Vector{ST}},
                V::Vector{Vector{ST}}, Λ::Vector{Vector{ST}},
                params::ParametersVPRKpSecondary{DT,TT,ET,D,S}) where {ST,DT,TT,ET,D,S}

    # @assert D == size(Q,1) == size(V,1) == size(Λ,1) == length(q̅)
    @assert S == length(Q) == length(V) == length(Λ)

    local y1::ST
    local y2::ST
    local y3::ST
    local y4::ST

    # compute Q
    for i in 1:S
        for k in 1:D
            y1 = y2 = y3 = y4 = 0
            for j in 1:S
                y1 += params.tab.q.a[i,j] * V[j][k]
                y2 += params.tab.q.â[i,j] * V[j][k]
                y3 += params.tab.q.a[i,j] * Λ[j][k]
                y4 += params.tab.q.â[i,j] * Λ[j][k]
            end
            Q[i][k] = params.q̅[k] + params.Δt * (y1 + y2 + y3 + y4)
        end
    end

    # compute q
    for k in 1:D
        y1 = y2 = y3 = y4 = 0
        for j in 1:S
            y1 += params.tab.q.b[j] * V[j][k]
            y2 += params.tab.q.b̂[j] * V[j][k]
            y3 += params.tab.q.b[j] * Λ[j][k]
            y4 += params.tab.q.b̂[j] * Λ[j][k]
        end
        q[k] = params.q̅[k] + params.Δt * (y1 + y2 + y3 + y4)
    end
end

@generated function compute_projection_vprk!(q::Vector{ST}, p::Vector{ST},
                Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, Λ::Vector{Vector{ST}},
                R::Vector{Vector{ST}}, Ψ::Vector{Vector{ST}},
                params::ParametersVPRKpSecondary{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    # create temporary vectors
    dH = zeros(ST,D)
    Ω = zeros(ST,D,D)

    quote
        local t₀::TT = params.t̅
        local t₁::TT = params.t̅ + params.Δt
        local tᵢ::TT

        # compute p=ϑ(q)
        params.equ.ϑ(t₁, q, p)

        for i in 1:S
            tᵢ = t₀ + params.Δt * params.tab.p.c[i]

            params.equ.g(tᵢ, Q[i], Λ[i], R[i])
            params.equ.Ω(tᵢ, Q[i], $Ω)
            params.equ.∇H(tᵢ, Q[i], $dH)

            # TODO Check if ω() returns Ω or Ω^T -> this decides the sign on dH below
            #      Seems to be Ω^T -> then dH should be subtracted

            Ψ[i] .= $Ω * V[i] .- $dH
        end
    end
end

function compute_rhs_vprk!(b::Vector{ST}, P::Vector{Vector{ST}}, F::Vector{Vector{ST}}, R::Vector{Vector{ST}},
                params::ParametersVPRKpSecondary{DT,TT,ET,D,S}) where {ST,DT,TT,ET,D,S}

    local z1::ST
    local z2::ST
    local z3::ST
    local z4::ST

    # compute b = - [(P-AF-AR)]
    for i in 1:S
        for k in 1:D
            z1 = z2 = z3 = z4 = 0
            for j in 1:S
                z1 += params.tab.p.a[i,j] * F[j][k]
                z2 += params.tab.p.â[i,j] * F[j][k]
                z3 += params.tab.p.a[i,j] * R[j][k]
                z4 += params.tab.p.â[i,j] * R[j][k]
            end
            b[D*(i-1)+k] = (P[i][k] - params.p̅[k]) - params.Δt * (z1 + z2 + z3 + z4)
        end
    end
end

function compute_rhs_vprk_projection!(b::Vector{ST}, p::Vector{ST},
                F::Vector{Vector{ST}}, R::Vector{Vector{ST}}, Ψ::Vector{Vector{ST}}, offset::Int,
                params::ParametersVPRKpSecondary{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    local z1::ST
    local z2::ST
    local z3::ST
    local z4::ST

    for i in 1:S
        for k in 1:D
            b[offset+D*(i-1)+k] = Ψ[i][k]
        end
    end

    # TODO # reactivate averaging of secondary constraints and primary constraint of solution

    # for i in 1:(S-1)
    #     for k in 1:D
    #         b[offset+D*(i-1)+k] = 0
    #         for j in 1:S
    #             b[offset+D*(i-1)+k] += params.ω[i,j] * Ψ[j][k]
    #         end
    #     end
    # end
    #
    # for k in 1:D
    #     z1 = z2 = z3 = z4 = 0
    #     for j in 1:S
    #         z1 += params.tab.p.b[j] * F[j][k]
    #         z2 += params.tab.p.b̂[j] * F[j][k]
    #         z3 += params.tab.p.b[j] * R[j][k]
    #         z4 += params.tab.p.b̂[j] * R[j][k]
    #     end
    #     b[offset+D*(S-1)+k] = (p[k] - params.p̅[k]) - params.Δt * (z1 + z2 + z3 + z4)
    # end
end

"Compute stages of variational partitioned Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST},
                params::ParametersVPRKpSecondary{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    cache = IntegratorCacheVPRK{ST, D, S}(true)

    function_stages = quote
        compute_stages_vprk!(x, $cache.q̃, $cache.p̃,
                                $cache.Q, $cache.V, $cache.Λ,
                                $cache.P, $cache.F, $cache.R,
                                $cache.Φ, params)

        # compute b = [P-AF-AR]
        compute_rhs_vprk!(b, $cache.P, $cache.F, $cache.R, params)

        # compute b = Φ
        compute_rhs_vprk_projection!(b, $cache.p̃, $cache.F, $cache.R, $cache.Φ, D*S, params)

        compute_rhs_vprk_correction!(b, $cache.V, params)
    end

    return function_stages
end

function initial_guess!(int::IntegratorVPRKpSecondary, sol::AtomicSolutionPODE)
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.ṽ,
                              tableau(int).q.c[i])

        for k in eachdim(int)
            int.solver.x[ndims(int)*(0*nstages(int)+i-1)+k] = int.cache.ṽ[k]
            int.solver.x[ndims(int)*(1*nstages(int)+i-1)+k] = 0
        end
    end
end

"Integrate ODE with variational partitioned Runge-Kutta integrator."
function integrate_step!(int::IntegratorVPRKpSecondary{DT,TT}, sol::AtomicSolutionPODE{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, sol)

    # compute initial guess
    initial_guess!(int, sol)

    # reset solution
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute final update
    compute_stages_vprk!(int.solver.x, int.cache.q̃, int.cache.p̃,
                          int.cache.Q, int.cache.V, int.cache.Λ,
                          int.cache.P, int.cache.F, int.cache.R,
                          int.cache.Φ, int.params)

    # compute unprojected solution
    update_solution!(sol.q, sol.q̃, int.cache.V, tableau(int).q.b, tableau(int).q.b̂, timestep(int))
    update_solution!(sol.p, sol.p̃, int.cache.F, tableau(int).p.b, tableau(int).p.b̂, timestep(int))

    # add projection to solution
    update_solution!(sol.q, sol.q̃, int.cache.Λ, tableau(int).q.b, tableau(int).q.b̂, timestep(int))
    update_solution!(sol.p, sol.p̃, int.cache.R, tableau(int).p.b, tableau(int).p.b̂, timestep(int))

    # copy solution to initial guess
    update!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)
end


text:
@testset "The Symmetric Matrices" begin
    #
    # Despite checking both cases of Eucliean Manifolds, this can serve as a
    # prototype testcase suite for your manifold
    #
    M = Symmetric(2)
    v1 = [1. 0.;0. 1.]
    v2 = [2. 1.; 1. 2.]
    d1 = [1. 1.;1. 1.]
    d2 = [1. 0.; 0. 0.]

    x1 = SymPoint(v1)
    x2 = SymPoint(v2)
    ξ1 = SymTVector(d1)
    ξ2 = SymTVector(d2)

    @test getValue(x1) == v1
    @test getValue(ξ1) == d1

    @test distance(M,x1,x2) == norm(v1-v2)

    @test dot(M,x1,ξ1,ξ2) == dot(d1,d2)

    @test exp(M,x1, ξ1) == SymPoint(v1 + d1)
    @test exp(M,x2, ξ2) == SymPoint(v2 + d2)

    @test log(M,x1,x2) == SymTVector(v2-v1)

    @test manifoldDimension(M) == 3
    @test manifoldDimension(M) == manifoldDimension(x1)

    @test norm(M,x1,ξ1) == norm(d1)

    @test parallelTransport(M,x1,x2,ξ1) == ξ1

    @test validateMPoint(M,x1)
    @test_throws DomainError !validateMPoint(M, SymPoint([1. 0.; 1. 1.]))
    @test_throws DomainError !validateMPoint(M, SymPoint([1. 0. 0.; 0. 1. 0.; 0. 0. 1.]))
    @test validateTVector(M,x1,ξ1)
    @test_throws DomainError validateTVector(M,x1,SymTVector([1. 0.; 1. 1.]))
    @test_throws DomainError validateTVector(M,x1,SymTVector([1. 0. 0.; 0. 1. 0.; 0. 0. 1.]))

    @test typeofTVector(typeof(x1)) == SymTVector{Float64}
    @test typeofMPoint(typeof(ξ1)) == SymPoint{Float64}
    @test typeofTVector(x1) == SymTVector{Float64}
    @test typeofMPoint(ξ1) == SymPoint{Float64}

    @test typicalDistance(M) == sqrt(2)

    @test zeroTVector(M,x1) == SymTVector(zeros(2,2))

    @test "$M" == "The Manifold of 2-by-2 symmetric matrices."
    @test "$(x1)" == "Sym([1.0 0.0; 0.0 1.0])"
    @test "$(ξ1)" == "SymT([1.0 1.0; 1.0 1.0])"
end

ext: jl


text:
#=== Absorbing Boundary ===#
function ModExpand(para::AbstractArray, ext::Int64, iflag::Union{Void,Int64})
    if iflag == 1 #free surface
        temp_vertical2 = repmat(para[end:end,:],ext,1)
        para = vcat(para, temp_vertical2)
        temp_horizon1 = repmat(para[:,1],1,ext)
        temp_horizon2 = repmat(para[:,end],1,ext)
        para = hcat(temp_horizon2,para,temp_horizon1)
    elseif iflag == 2 # unlimited medium
        temp_vertical1 = repmat(para[1:1,:],ext,1)
        temp_vertical2 = repmat(para[end:end,:],ext,1)
        para = vcat(temp_vertical1,para,temp_vertical2)
        temp_horizon1 = repmat(para[:,1],1,ext)
        temp_horizon2 = repmat(para[:,end],1,ext)
        para = hcat(temp_horizon1,para,temp_horizon2)
    elseif iflag == nothing
        temp_horizon1 = repmat(para[:,1],1,ext)
        temp_horizon2 = repmat(para[:,end],1,ext)
        para = hcat(temp_horizon1,para,temp_horizon2)
    end
    return para
end

############ elastic medium 2d ############
function model{T1,T2,T3<:Real}(
    pvel::Union{String,<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64}}},
    svel::Union{Void,String,<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64}}},
    rho::Union{<:Real,String},
    DZ::Union{<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    HX::T1,
    pkf::T2,
    T::T3,
    ext::Int64,
    iflag::Int64,
    dx::Union{Void,<:Real},
    dz::Union{Void,<:Real},
    dt::Union{Void,<:Real},
    nT::Union{Void,Int64})

    if typeof(pvel) == String
        if !(typeof(DZ) == Int64 && typeof(HX) == Int64)
            error("Depth and Horzion have to be Integer grid number")
        elseif dx == nothing && dz == nothing
            error("Set up grid size dx and dz")
        elseif !(typeof(DZ) == Int64 && typeof(HX) == Int64)&&dx == nothing && dz == nothing
            error("1. Depth and Horzion have to be Integer grid number. 2.Set up grid size dx and dz.")
        end
        if pvel[end-2:end] == "bin"
            fid = open(pvel,"r")
            pvel = read(fid,Float64, DZ*HX,1)
            close(fid)
        elseif pvel[end-2:end] == "dat"
            fid = open(pvel,"r")
            pvel = readdlm(fid)
            close(fid)
        else error("P velocity file is not valid")
        end
        if svel == nothing
            svel = pvel./sqrt(3);
        elseif svel[end-2:end] == "bin"
            fid = open(svel,"r")
            svel = read(fid,Float64, DZ*HX,1)
            close(fid)
        elseif svel[end-2:end] == "dat"
            fid = open(svel,"r")
            svel = readdlm(fid)
            close(fid)
        else error("S velocity file is not valid")
        end
        max_pvel = maximum(pvel)
        min_pvel = minimum(pvel)
        max_svel = maximum(svel)
        min_svel = minimum(svel)

        if !(pkf < 2*min_svel/3/dx && pkf < 2*min_svel/3/dz)
            error("Peak frequency is too high for selected grid size")
        end
        if dt == nothing
            dt = 0.406 * (min_svel / 10 / pkf) / max_pvel
            nT = Int64(round(T/dt))
        elseif dt !=nothing
            if dt > 0.6*(min_svel / 10 / pkf)/max_pvel
                error("Simulation does not satisfy stablity criterion:
                maximum(pvel)*dt/dz < 0.6.")
            end
            if nT != nothing
                T= nT*dt
            elseif nT == nothing
                nT = Int64(round(T/dt))
            end
        end
        pvel = Float64.(reshape(pvel, DZ,HX))
        svel = Float64.(reshape(svel, DZ,HX))
        if typeof(rho) == String
            if rho[end-2:end] == "bin"
                fid = open(rho,"r")
                rho = read(fid,Float64, DZ*HX,1)
                close(fid)
            elseif rho[end-2:end] == "dat"
                fid = open(rho,"r")
                rho = readdlm(rho)
                close(fid)
            else error("Density file is not valid")
            end
            Rho = Float64.(reshape(rho, DZ,HX))
        elseif typeof(rho)<:Real
            Rho = Float64.(zeros(DZ,HX) .+ rho)
        else error("Give a density value or file path")
        end
        nDZ = DZ
        nHX = HX
        DZ = DZ*dz
        HX = HX*dx

    elseif typeof(pvel) != String
        if svel == nothing
            svel = pvel./sqrt(3);
        end
        if !(length(pvel)==length(svel)==length(DZ))
            error("Layer numbers of pvel, svel and depth should be the same. Please give a depth for each layer.")
        end
        max_pvel = maximum(pvel)
        min_pvel = minimum(pvel)
        max_svel = maximum(svel)
        min_svel = minimum(svel)
        # Spatial interval and time interval
        if !(max_svel > 0 && max_pvel > 0 && min_svel > 0 && max_pvel > 0)
            error("Medium velocities have to be larger than 0")
        else
            dz = min_svel / 10 / pkf
            dx = min_svel / 10 / pkf
        end
        if dt == nothing
            dt = 0.406 * (min_svel / 10 / pkf) / max_pvel
            nT = Int64(round(T/dt))
        elseif dt !=nothing
            if dt > 0.6*(min_svel / 10 / pkf)/max_pvel
                error("Simulation does not satisfy stablity criterion:
                maximum(P velocity)*dt/dz < 0.6.")
            end
            if nT != nothing
                T= nT*dt
            elseif nT == nothing
                nT = Int64(round(T/dt))
            end
        end
        nDZ = 0
        nlayerdep = zeros(Int64,length(DZ))
        if typeof(DZ) <: Array{<:Real,1} || typeof(DZ) <: Real
            for i = 1:length(DZ)
                nlayerdep[i] = Int64(round(DZ[i]/dz))
                nDZ = nDZ + nlayerdep[i]
            end
            DZ = sum(DZ)
        else
            nlayerdep[1] = Int64(round(DZ[1]/dz))
            nDZ = Int64(round(DZ[1]/dz))
            for i = 2:length(DZ)
                nlayerdep[i] = Int64(round((DZ[i]-DZ[i-1])/dz))
                nDZ = nDZ + nlayerdep[i]
            end
            DZ = DZ[end]
        end

        nHX = Int64(round(HX/dx))
        Pvel = zeros(Float64,nDZ,nHX)
        Svel = zeros(Float64,nDZ,nHX)
        if typeof(rho) == String
            if rho[end-2:end] == "bin"
                fid = open(rho,"r")
                rho = read(fid,Float64, nDZ*nHX,1)
                close(fid)
            elseif rho[end-2:end] == "dat"
                fid = open(rho,"r")
                rho = readdlm(rho)
                close(fid)
            else error("Density file is not valid")
            end
            Rho = Float64.(reshape(rho, nDZ,nHX))
        elseif typeof(rho)<:Real
            Rho = Float64.(zeros(nDZ,nHX) .+ rho)
        else error("Give a density value or file path")
        end

        # velocity interval and depth interval

        Pvel[1:nlayerdep[1],:] = Float64.(pvel[1])
        Svel[1:nlayerdep[1],:] = Float64.(svel[1])
        SumDep = 0
        for iL = 2 : length(pvel)
            SumDep = SumDep + nlayerdep[iL-1]
            Pvel[SumDep+1:SumDep+nlayerdep[iL],:] = Float64.(pvel[iL])
            Svel[SumDep+1:SumDep+nlayerdep[iL],:] = Float64.(svel[iL])
        end
        pvel = Float64.(Pvel)
        svel = Float64.(Svel)
    end

    pvel = ModExpand(pvel, ext, iflag) # 1 : free surface; 2: unlimited surface
    svel = ModExpand(svel, ext, iflag)
    rho = ModExpand(Rho, ext, iflag)
    BDnDZ,BDnHX = size(pvel)
    lambda = Array{Float64,2}(BDnDZ,BDnHX)
    mu = Array{Float64,2}(BDnDZ,BDnHX)
    Rho = Array{Float64,2}(BDnDZ,BDnHX)
    for i in 1 : BDnDZ
        for j in 1 : BDnHX
            lambda[i,j] = pvel[i,j]^2*rho[i,j]  - 2*svel[i,j]^2*rho[i,j]
            mu[i,j] = svel[i,j]^2*rho[i,j]
            Rho[i,j] = rho[i,j]
        end
    end

    return elastic2d(pvel,svel,Rho,lambda,mu,dx,dz,dt,DZ,HX,nDZ,nHX,BDnDZ,BDnHX,T,nT,pkf,ext,iflag)
end

############ acoustic medium 2d ############
  function model{T1,T2,T3<:Real}(
    pvel::Union{String,<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    rho::Union{<:Real,String},
    DZ::Union{<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    HX::T1,
    pkf::T2,
    T::T3,
    ext::Int64,
    iflag::Int64,
    dx::Union{Void,<:Real},
    dz::Union{Void,<:Real},
    dt::Union{Void,<:Real},
    nT::Union{Void,Int64})

    if typeof(pvel) == String
        if !(typeof(DZ) == Int64 && typeof(HX) == Int64)
            error("Depth and Horzion have to be Integer grid number")
        elseif dx == nothing && dz == nothing
            error("Set up grid size dx and dz")
        elseif !(typeof(DZ) == Int64 && typeof(HX) == Int64)&&dx == nothing && dz == nothing
            error("1. Depth and Horzion have to be Integer grid number. 2.Set up grid size dx and dz.")
        end
        if pvel[end-2:end] == "bin"
            fid = open(pvel,"r")
            pvel = read(fid,Float64, DZ*HX,1)
            close(fid)
        elseif pvel[end-2:end] == "dat"
            fid = open(pvel,"r")
            pvel = readdlm(fid)
            close(fid)
        else error("P velocity file is not valid")
        end

        max_pvel = maximum(pvel)
        min_pvel = minimum(pvel)

        if !(pkf < 2*min_pvel/3/dx && pkf < 2*min_pvel/3/dz)
            error("Peak frequency is too high for selected grid size")
        end
        if dt == nothing
            dt = 0.406 * (min_pvel / 10 / pkf) / max_pvel
            nT = Int64(round(T/dt))
        elseif dt !=nothing
            if dt > 0.6*(min_pvel / 10 / pkf)/max_pvel
                error("Simulation does not satisfy stablity criterion:
                maximum(pvel)*dt/dz < 0.6.")
            end
            if nT != nothing
                T= nT*dt
            elseif nT == nothing
                nT = Int64(round(T/dt))
            end
        end
        pvel = Float64.(reshape(pvel, DZ,HX))
        if typeof(rho) == String
            if rho[end-2:end] == "bin"
                fid = open(rho,"r")
                rho = read(fid,Float64, DZ*HX,1)
                close(fid)
            elseif rho[end-2:end] == "dat"
                fid = open(rho,"r")
                rho = readdlm(rho)
                close(fid)
            else error("Density file is not valid")
            end
            Rho = Float64.(reshape(rho, DZ,HX))
        elseif typeof(rho)<:Real
            Rho = Float64.(zeros(DZ,HX) .+ rho)
        else error("Give a density value or file path")
        end

        nDZ = DZ
        nHX = HX
        DZ = DZ*dz
        HX = HX*dx

    elseif typeof(pvel) != String

        if !(length(pvel)==length(DZ))
            error("Layer numbers of pvel, svel and depth should be the same. Please give a depth for each layer.")
        end
        max_pvel = maximum(pvel)
        min_pvel = minimum(pvel)
        # Spatial interval and time interval
        if !(max_pvel > 0 && max_pvel > 0)
            error("Medium velocities have to be larger than 0")
        else
            dz = min_pvel / 10 / pkf
            dx = min_pvel / 10 / pkf
        end
        if dt == nothing
            dt = 0.406 * (min_pvel / 10 / pkf) / max_pvel
            nT = Int64(round(T/dt))
        elseif dt !=nothing
            if dt > 0.6*(min_svel / 10 / pkf)/max_pvel
                error("Simulation does not satisfy stablity criterion:
                maximum(P velocity)*dt/dz < 0.6.")
            end
            if nT != nothing
                T= nT*dt
            elseif nT == nothing
                nT = Int64(round(T/dt))
            end
        end
        nDZ = 0
        nlayerdep = zeros(Int64,length(DZ))
        if typeof(DZ) <: Array{<:Real,1} || typeof(DZ) <: Real
            for i = 1:length(DZ)
                nlayerdep[i] = Int64(round(DZ[i]/dz))
                nDZ = nDZ + nlayerdep[i]
            end
            DZ = sum(DZ)
        else
            nlayerdep[1] = Int64(round(DZ[1]/dz))
            nDZ = Int64(round(DZ[1]/dz))
            for i = 2:length(DZ)
                nlayerdep[i] = Int64(round((DZ[i]-DZ[i-1])/dz))
                nDZ = nDZ + nlayerdep[i]
            end
            DZ = DZ[end]
        end

        nHX = Int64(round(HX/dx))
        Pvel = zeros(Float64,nDZ,nHX)

        # velocity interval and depth interval

        Pvel[1:nlayerdep[1],:] = Float64.(pvel[1])
        SumDep = 0
        for iL = 2 : length(pvel)
            SumDep = SumDep + nlayerdep[iL-1]
            Pvel[SumDep+1:SumDep+nlayerdep[iL],:] = Float64.(pvel[iL])
        end
        pvel = Float64.(Pvel)
    end

    if typeof(rho) == String
        if rho[end-2:end] == "bin"
            fid = open(rho,"r")
            rho = read(fid,Float64, nDZ*nHX,1)
            close(fid)
        elseif rho[end-2:end] == "dat"
            fid = open(rho,"r")
            rho = readdlm(rho)
            close(fid)
        else error("Density file is not valid")
        end
        Rho = Float64.(reshape(rho, nDZ,nHX))
    elseif typeof(rho)<:Real
        Rho = Float64.(zeros(nDZ,nHX) .+ rho)
    else error("Give a density value or file path")
    end
    pvel = ModExpand(pvel, ext, iflag) # 1 : free surface; 2: unlimited surface
    rho = ModExpand(Rho, ext, iflag)

    BDnDZ, BDnHX = size(pvel)
    lambda = Array{Float64,2}(BDnDZ,BDnHX)
    Rho = Array{Float64,2}(BDnDZ,BDnHX)
    for i in 1 : BDnDZ
        for j in 1 : BDnHX
            lambda[i,j] = pvel[i,j]^2*rho[i,j]
            Rho[i,j] = rho[i,j]
        end
    end
    return acoustic2d(pvel,Rho,lambda,dx,dz,dt,DZ,HX,nDZ,nHX,BDnDZ,BDnHX,T,nT,pkf,ext,iflag)
end

############ elastic medium 3d ############
 function model{T1,T2,T3,T4<:Real}(
    pvel::Union{String,<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    svel::Union{Void,String,<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    rho::Union{<:Real,String},
    DZ::Union{<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    HX::T1,
    HY::T2,
    pkf::T3,
    T::T4,
    ext::Int64,
    iflag::Int64,
    dx::Union{Void,<:Real},
    dy::Union{Void,<:Real},
    dz::Union{Void,<:Real},
    dt::Union{Void,<:Real},
    nT::Union{Void,Int64})

    if typeof(pvel) == String
        if !(typeof(DZ) == Int64 && typeof(HX) == Int64 && typeof(HY) == Int64)
            error("Depth and Horzion have to be Integer grid number")
        elseif dx == nothing && dz == nothing && dy == nothing
            error("Set up grid size dx, dy and dz")
        elseif !(typeof(DZ) == Int64 && typeof(HX) == Int64 && typeof(HY) == Int64) && dx == nothing && dz == nothing && dy == nothing
            error("1. Depth and Horzion have to be Integer grid number. 2.Set up grid size dx, dy and dz.")
        end
        if pvel[end-2:end] == "bin"
            fid = open(pvel,"r")
            pvel = read(fid,Float64, DZ*HX*HY,1)
            close(fid)
        elseif pvel[end-2:end] == "dat"
            fid = open(pvel,"r")
            pvel = readdlm(fid)
            close(fid)
        else error("P velocity file is not valid")
        end
        if svel == nothing
            svel = pvel./sqrt(3);
        elseif svel[end-2:end] == "bin"
            fid = open(svel,"r")
            svel = read(fid,Float64, DZ*HX*HY,1)
            close(fid)
        elseif svel[end-2:end] == "dat"
            fid = open(svel,"r")
            svel = readdlm(fid)
            close(fid)
        else error("S velocity file is not valid")
        end
        max_pvel = maximum(pvel)
        min_pvel = minimum(pvel)
        max_svel = maximum(svel)
        min_svel = minimum(svel)

        if !(pkf < 2*min_svel/3/dx && pkf < 2*min_svel/3/dz && pkf < 2*min_svel/3/dy)
            error("Peak frequency is too high for selected grid size")
        end
        if dt == nothing
            dt = 0.406 * (min_svel / 10 / pkf) / max_pvel
            nT = Int64(round(T/dt))
        elseif dt != nothing
            if dt > 0.6 * (min_svel / 10 / pkf) / max_pvel
                error("Simulation does not satisfy stablity criterion:
                maximum(pvel)*dt/dz < 0.6.")
            end
            if nT != nothing
                T= nT*dt
            elseif nT == nothing
                nT = Int64(round(T/dt))
            end
        end
        Pvel = reshape(pvel, DZ,HX,HY)
        Svel = reshape(svel, DZ,HX,HY)
        if typeof(rho) == String
            if rho[end-2:end] == "bin"
                fid = open(rho,"r")
                rho = read(fid,Float64, DZ*HX*HY,1)
                close(fid)
            elseif rho[end-2:end] == "dat"
                fid = open(rho,"r")
                rho = readdlm(rho)
                close(fid)
            else error("Density file is not valid")
            end
            Rho = Float64.(reshape(rho, DZ,HX,HY))
        elseif typeof(rho)<:Real
            Rho = Float64.(zeros(DZ,HX,HY) .+ rho)
        else error("Give a density value or file path")
        end
        nDZ = DZ
        nHX = HX
        nHY = HY
        DZ = DZ*dz
        HX = HX*dx
        HY = HY*dy

    elseif typeof(pvel) != String
        if svel == nothing
            svel = pvel/sqrt(3);
        end
        if !(length(pvel)==length(svel)==length(DZ))
            error("Layer numbers of pvel, svel and depth should be the same. Please give a depth for each layer.")
        end
        max_pvel = maximum(pvel)
        min_pvel = minimum(pvel)
        max_svel = maximum(svel)
        min_svel = minimum(svel)
        # Spatial interval and time interval
        if !(max_svel > 0 && max_pvel > 0 && min_svel > 0 && max_pvel > 0)
            error("Medium velocities have to be larger than 0")
        else
            dz = min_svel / 10 / pkf
            dx = min_svel / 10 / pkf
            dy = min_svel / 10 / pkf
        end
        if dt == nothing
            dt = 0.406 * (min_svel / 10 / pkf) / max_pvel
            nT = Int64(round(T/dt))
        elseif dt != nothing
            if dt > 0.6*(min_svel / 10 / pkf)/max_pvel
                error("Simulation does not satisfy stablity criterion:
                maximum(P velocity)*dt/dz < 0.6.")
            end
            if nT != nothing
                T= nT*dt
            elseif nT == nothing
                nT = Int64(round(T/dt))
            end
        end
        nDZ = 0
        nlayerdep = zeros(Int64,length(DZ))
        if typeof(DZ) <: Array{<:Real,1} || typeof(DZ) <: Real
            for i = 1:length(DZ)
                nlayerdep[i] = Int64(round(DZ[i]/dz))
                nDZ = nDZ + nlayerdep[i]
            end
            DZ = sum(DZ)
        else
            nlayerdep[1] = Int64(round(DZ[1]/dz))
            nDZ = Int64(round(DZ[1]/dz))
            for i = 2:length(DZ)
                nlayerdep[i] = Int64(round((DZ[i]-DZ[i-1])/dz))
                nDZ = nDZ + nlayerdep[i]
            end
            DZ = DZ[end]
        end

        nHX = Int64(round(HX/dx))
        nHY = Int64(round(HY/dy))
        Pvel = zeros(nDZ,nHX,nHY)
        Svel = zeros(nDZ,nHX,nHY)
        if typeof(rho) == String
            if rho[end-2:end] == "bin"
                fid = open(rho,"r")
                rho = read(fid,Float64, nDZ*nHX*nHY,1)
                close(fid)
            elseif rho[end-2:end] == "dat"
                fid = open(rho,"r")
                rho = readdlm(rho)
                close(fid)
            else error("Density file is not valid")
            end
            Rho = Float64.(reshape(rho, nDZ,nHX,nHY))
        elseif typeof(rho)<:Real
            Rho = Float64.(zeros(nDZ,nHX,nHY) .+ rho)
        else error("Give a density value or file path")
        end
        # velocity interval and depth interval

        Pvel[1:nlayerdep[1],:,:] = pvel[1]
        Svel[1:nlayerdep[1],:,:] = svel[1]
        SumDep = 0
        for iL = 2 : length(pvel)
            SumDep = SumDep + nlayerdep[iL-1]
            Pvel[SumDep+1:SumDep+nlayerdep[iL],:,:] = pvel[iL]
            Svel[SumDep+1:SumDep+nlayerdep[iL],:,:] = svel[iL]
        end
    end
    if iflag == 1 # free surface
        BDnDZ = nDZ + ext
    elseif iflag == 2
        BDnDZ = nDZ + 2*ext
    end
    BDnHX = nHX + 2*ext
    BDnHY = nHY + 2*ext
    pvel = zeros(BDnDZ,BDnHX,BDnHY)
    svel = zeros(BDnDZ,BDnHX,BDnHY)
     rho = zeros(BDnDZ,BDnHX,BDnHY)
    # 1 : free surface; 2: unlimited surface
    pvel[:,:,1] = ModExpand(Pvel[:,:,1], ext, iflag)
    svel[:,:,1] = ModExpand(Svel[:,:,1], ext, iflag)
    rho[:,:,1] = ModExpand( Rho[:,:,1], ext, iflag)
    for i in 2:BDnHY
        pvel[:,:,i] = pvel[:,:,1]
        svel[:,:,i] = svel[:,:,1]
         rho[:,:,i] = rho[:,:,1]
    end
    for i in 1:BDnDZ
        pvel[i,:,:] = ModExpand(pvel[i,:,ext+1:ext+nHY], ext, nothing)
        svel[i,:,:] = ModExpand(svel[i,:,ext+1:ext+nHY], ext, nothing)
         rho[i,:,:] = ModExpand( rho[i,:,ext+1:ext+nHY], ext, nothing)
    end
    # dz= Float64(dz)
    # dx= Float64(dx)
    # dy= Float64(dy)
    # dt= Float64(dt)
    lambda = Array{Float64,3}(BDnDZ,BDnHX,BDnHY)
    mu = Array{Float64,3}(BDnDZ,BDnHX,BDnHY)
    Rho = Array{Float64,3}(BDnDZ,BDnHX,BDnHY)
    for i in 1 : BDnDZ
        for j in 1 : BDnHX
            for k in 1 : BDnHY
                lambda[i,j,k] = pvel[i,j,k]^2*rho[i,j,k] - 2*svel[i,j,k]^2*rho[i,j,k]
                mu[i,j,k] = svel[i,j,k]^2*rho[i,j,k]
                Rho[i,j,k] = rho[i,j,k]
            end
        end
    end
    return elastic3d(pvel,svel,Rho,lambda,mu,dx,dy,dz,dt,
    DZ,HX,HY,nDZ,nHX,nHY,BDnDZ,BDnHX,BDnHY,T,nT,pkf,ext,iflag)
end

############ acoustic medium 3d ############
 function model{T1,T2,T3,T4<:Real}(
    pvel::Union{String,<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    rho::Union{<:Real,String},
    DZ::Union{<:Real,Array{<:Real,1},StepRange{Int64,Int64},UnitRange{Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},
    HX::T1,
    HY::T2,
    pkf::T3,
    T::T4,
    ext::Int64,
    iflag::Int64,
    dx::Union{Void,<:Real},
    dy::Union{Void,<:Real},
    dz::Union{Void,<:Real},
    dt::Union{Void,<:Real},
    nT::Union{Void,Int64})

    if typeof(pvel) == String
        if !(typeof(DZ) == Int64 && typeof(HX) == Int64 && typeof(HY) == Int64)
            error("Depth and Horzion have to be Integer grid number")
        elseif dx == nothing && dz == nothing && dy == nothing
            error("Set up grid size dx and dz")
        elseif !(typeof(DZ) == Int64 && typeof(HX) == Int64 && typeof(HY) == Int64) && dx == nothing && dz == nothing && dx == nothing
            error("1. Depth and Horzion have to be Integer grid number. 2.Set up grid size dx and dz.")
        end
        if pvel[end-2:end] == "bin"
            fid = open(pvel,"r")
            pvel = read(fid,Float64, DZ*HX*HY,1)
            close(fid)
        elseif pvel[end-2:end] == "dat"
            fid = open(pvel,"r")
            pvel = readdlm(fid)
            close(fid)
        else error("P velocity file is not valid")
        end

        max_pvel = maximum(pvel)
        min_pvel = minimum(pvel)

        if !(pkf < 2*min_pvel/3/dx && pkf < 2*min_pvel/3/dz && pkf < 2*min_pvel/3/dy)
            error("Peak frequency is too high for selected grid size")
        end
        if dt == nothing
            dt = 0.406 * (min_pvel / 10 / pkf) / max_pvel
            nT = Int64(round(T/dt))
        elseif dt !=nothing
            if dt > 0.6*(min_pvel / 10 / pkf)/max_pvel
                error("Simulation does not satisfy stablity criterion:
                maximum(pvel)*dt/dz < 0.6.")
            end
            if nT != nothing
                T= nT*dt
            elseif nT == nothing
                nT = Int64(round(T/dt))
            end
        end
        Pvel = reshape(pvel, DZ,HX,HY)
        if typeof(rho) == String
            if rho[end-2:end] == "bin"
                fid = open(rho,"r")
                rho = read(fid,Float64, DZ*HX*HY,1)
                close(fid)
            elseif rho[end-2:end] == "dat"
                fid = open(rho,"r")
                rho = readdlm(rho)
                close(fid)
            else error("Density file is not valid")
            end
            Rho = Float64.(reshape(rho, DZ,HX,HY))
        elseif typeof(rho)<:Real
            Rho = Float64.(zeros(DZ,HX,HY) .+ rho)
        else error("Give a density value or file path")
        end
        nDZ = DZ
        nHX = HX
        nHY = HY
        DZ = DZ*dz
        HX = HX*dx
        HY = HY*dy

    elseif typeof(pvel) != String

        if !(length(pvel)==length(DZ))
            error("Layer numbers of pvel, svel and depth should be the same. Please give a depth for each layer.")
        end
        max_pvel = maximum(pvel)
        min_pvel = minimum(pvel)
        # Spatial interval and time interval
        if !(max_pvel > 0 && max_pvel > 0)
            error("Medium velocities have to be larger than 0")
        else
            dz = min_pvel / 10 / pkf
            dx = min_pvel / 10 / pkf
            dy = min_pvel / 10 / pkf
        end
        if dt == nothing
            dt = 0.406 * (min_pvel / 10 / pkf) / max_pvel
            nT = Int64(round(T/dt))
        elseif dt !=nothing
            if dt > 0.6*(min_pvel / 10 / pkf)/max_pvel
                error("Simulation does not satisfy stablity criterion:
                maximum(P velocity)*dt/dz < 0.6.")
            end
            if nT != nothing
                T= nT*dt
            elseif nT == nothing
                nT = Int64(round(T/dt))
            end
        end
        nDZ = 0
        nlayerdep = zeros(Int64,length(DZ))
        if typeof(DZ) <: Array{<:Real,1} || typeof(DZ) <: Real
            for i = 1:length(DZ)
                nlayerdep[i] = Int64(round(DZ[i]/dz))
                nDZ = nDZ + nlayerdep[i]
            end
            DZ = sum(DZ)
        else
            nlayerdep[1] = Int64(round(DZ[1]/dz))
            nDZ = Int64(round(DZ[1]/dz))
            for i = 2:length(DZ)
                nlayerdep[i] = Int64(round((DZ[i]-DZ[i-1])/dz))
                nDZ = nDZ + nlayerdep[i]
            end
            DZ = DZ[end]
        end

        nHX = Int64(round(HX/dx))
        nHY = Int64(round(HY/dy))
        Pvel = zeros(nDZ,nHX,nHY)
        if typeof(rho) == String
            if rho[end-2:end] == "bin"
                fid = open(rho,"r")
                rho = read(fid,Float64, nDZ*nHX*nHY,1)
                close(fid)
            elseif rho[end-2:end] == "dat"
                fid = open(rho,"r")
                rho = readdlm(rho)
                close(fid)
            else error("Density file is not valid")
            end
            Rho = Float64.(reshape(rho, nDZ,nHX,nHY))
        elseif typeof(rho)<:Real
            Rho = Float64.(zeros(nDZ,nHX,nHY) .+ rho)
        else error("Give a density value or file path")
        end
        # velocity interval and depth interval
        for i in 1:nlayerdep[1]
            Pvel[i,:,:] = Pvel[i,:,:] .+ pvel[1]
        end
        SumDep = 0
        for iL = 2 : length(pvel)
            SumDep = SumDep + nlayerdep[iL-1]
            Pvel[SumDep+1:SumDep+nlayerdep[iL],:,:] = Pvel[SumDep+1:SumDep+nlayerdep[iL],:,:] .+ pvel[iL]
        end
    end

    # 1 : free surface; 2: unlimited surface
    if iflag == 1
        BDnDZ = nDZ + ext
    elseif iflag == 2
        BDnDZ = nDZ + 2*ext
    end
    BDnHX = nHX + 2*ext
    BDnHY = nHY + 2*ext
    pvel = zeros(BDnDZ,BDnHX,BDnHY)
     rho = zeros(BDnDZ,BDnHX,BDnHY)
    # 1 : free surface; 2: unlimited surface
    pvel[:,:,1] = ModExpand(Pvel[:,:,1], ext, iflag)
    rho[:,:,1] = ModExpand( Rho[:,:,1], ext, iflag)
    for i in 2:BDnHY
        pvel[:,:,i] = pvel[:,:,1]
         rho[:,:,i] = rho[:,:,1]
    end
    for i in 1:BDnDZ
        pvel[i,:,:] = ModExpand(pvel[i,:,ext+1:ext+nHY], ext, nothing)
         rho[i,:,:] = ModExpand( rho[i,:,ext+1:ext+nHY], ext, nothing)
    end

    lambda = Array{Float64,3}(BDnDZ,BDnHX,BDnHY)
    Rho = Array{Float64,3}(BDnDZ,BDnHX,BDnHY)
    for i in 1 : BDnDZ
        for j in 1 : BDnHX
            for k in 1 : BDnHY
                lambda[i,j,k] = pvel[i,j,k]^2*rho[i,j,k]
                Rho[i,j,k] = rho[i,j,k]
            end
        end
    end
    return acoustic3d(pvel,Rho,lambda,dx,dy,dz,dt,
    DZ,HX,HY,nDZ,nHX,nHY,BDnDZ,BDnHX,BDnHY,T,nT,pkf,ext,iflag)
end


text:
using CSV, DataFrames, CategoricalArrays, Statistics, Chain, Dates
using Plots, StatsPlots
gr(dpi=300)

date_format = DateFormat("y-m-d H:M:S")

df =  CSV.read("data/flu_pt_raw_tweets.csv", DataFrame;
               header=["index", "id", "created_at", "text", "user",
                       "place", "user_place", "country", "coordinates",
                       "undefined_col", "undefined_col2", "undefined_col3"],
               types=[Int, Int, DateTime, String, Int,
                      String, String, String, String,
                      String, Int, Int],
                dateformat=date_format)

println(describe(df))

# filter!(:created_at => <(Dates.DateTime(2021), df)

# Total Tweets per day
gdf = @chain df begin
    transform(:created_at => ByRow(Dates.Date) => :created_at)
    groupby(:created_at)
end

tweet_plot = @df combine(gdf, nrow) plot(:created_at,
                            :nrow,
                            label=false,
                            xlabel="Day",
                            ylabel="Total Tweets Scraped",
                            yformatter=y -> string(round(Int64, y ÷ 1_000)) * "K",
                            xrotation=45)

savefig(tweet_plot, "images/tweets_scraped.png")

# Twitter Predictions
tweets_pred = CSV.read("data/tweets_pred.csv", DataFrame)
tweets_pred_plot = @df tweets_pred plot(:date,
                            :predicted,
                            label=false,
                            xlabel="Day",
                            ylabel="Total Tweets Predicted",
                            yformatter=y -> string(round(Int64, y ÷ 1_000)) * "K",
                            xrotation=45)
savefig(tweet_plot, "images/tweets_predicted.png")

# caso full here: https://data.brasil.io/dataset/covid19/caso_full.csv.gz
cases_br = @chain CSV.read("data/caso_full.csv", DataFrame) begin
    filter([:date, :city] => (date, city) -> date < Dates.Date("2021-01-01") && date > Dates.Date("2020-04-01") && ismissing(city), _)
    groupby(:date)
    combine(
        [:estimated_population_2019,
         :last_available_confirmed_per_100k_inhabitants,
         :last_available_deaths,
         :new_confirmed,
         :new_deaths] .=> sum .=>
         [:estimated_population_2019,
         :last_available_confirmed_per_100k_inhabitants,
         :last_available_deaths,
         :new_confirmed,
         :new_deaths]
    )
end

brazil_covid = @df cases_br plot(:date,
                            :last_available_deaths,
                            label="Total Deaths",
                            xlabel="Day",
                            legend=:topleft,
                            # ylabel="Total Tweets Predicted",
                            yformatter=y -> string(round(Int64, y ÷ 1_000)) * "K",
                            xrotation=45)
@df cases_br plot!(brazil_covid, :date, :new_confirmed, label="Daily Confirmed Cases")
savefig(brazil_covid, "images/brazil_2020.png")

max_stars_repo_name: codatmo/Brazil-Tweet-Classifier

max_issues_repo_name: codatmo/Brazil-Tweet-Classifier

max_forks_repo_name: codatmo/Brazil-Tweet-Classifier


text:
export sphplot

@recipe f(::Type{Vector{SVectorF{1}}}, array::Vector{SVectorF{1}}) =
    [x[] for x in array]

@userplot SPHPlot
@recipe function f(p::SPHPlot)
    x, y = p.args
    seriestype := :scatter
    legend := false
    x, y
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
function _end(obj::RealMatrixPreservingVisitor)
    return jcall(obj, "end", jdouble, ())
end

function start(obj::RealMatrixPreservingVisitor, arg0::jint, arg1::jint, arg2::jint, arg3::jint, arg4::jint, arg5::jint)
    return jcall(obj, "start", void, (jint, jint, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4, arg5)
end

function visit(obj::RealMatrixPreservingVisitor, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "visit", void, (jint, jint, jdouble), arg0, arg1, arg2)
end


text:
using Matching
using Test

agent_a = Agent("a")
agent_b = Agent("b")
agent_c = Agent("c")
agent_1 = Agent("1")
agent_2 = Agent("2")
agent_3 = Agent("3")

rankings = RankedPreferences()
add_rank!(rankings, agent_a, [agent_2, agent_1, agent_3])
add_rank!(rankings, agent_b, [agent_3, agent_2, agent_1])
add_rank!(rankings, agent_c, [agent_1, agent_3, agent_2])

add_rank!(rankings, agent_1, [agent_b, agent_a, agent_c])
add_rank!(rankings, agent_2, [agent_c, agent_b, agent_a])
add_rank!(rankings, agent_3, [agent_a, agent_c, agent_b])

matching_1 = gale_shapley([agent_1, agent_2, agent_3], [agent_a, agent_b, agent_c], rankings)
matching_1 = matching_1.results
possible_stable_matching_1 = matching_1[agent_a] == agent_1 && matching_1[agent_b] == agent_2 && matching_1[agent_c] == agent_3
possible_stable_matching_2 = matching_1[agent_a] == agent_3 && matching_1[agent_b] == agent_1 && matching_1[agent_c] == agent_2
@test possible_stable_matching_2 || possible_stable_matching_1

# all from the proposing pool should get their first choice
matching_2 = gale_shapley([agent_a, agent_b, agent_c], [agent_1, agent_2, agent_3], rankings)
print(matching_2)
matching_2 = matching_2.results
@test matching_2[agent_2] == agent_a
@test matching_2[agent_3] == agent_b
@test matching_2[agent_1] == agent_c

# test for non-trivial cycles (i.e. > 2 people) and multiple rounds
rankings= RankedPreferences()
add_rank!(rankings, agent_a, [agent_1, agent_2, agent_3])
add_rank!(rankings, agent_b, [agent_1, agent_3, agent_2])
add_rank!(rankings, agent_c, [agent_1, agent_3, agent_2])

add_rank!(rankings, agent_1, [agent_c, agent_b, agent_a])
add_rank!(rankings, agent_2, [agent_c, agent_b, agent_a])
add_rank!(rankings, agent_3, [agent_c, agent_a, agent_b])
matching_3 = ttc_with_counter([agent_a, agent_b, agent_c], [agent_1, agent_2, agent_3], rankings)
print(matching_3)
matching_3 = matching_3.results
@test matching_3[agent_a] == agent_2
@test matching_3[agent_b] == agent_3
@test matching_3[agent_c] == agent_1

# test capacity > 1
agent_a = Agent("a")
agent_b = Agent("b")
agent_c = Agent("c")
agent_1 = Agent("1", 2, 1)
agent_2 = Agent("2")
agent_3 = Agent("3")
rankings= RankedPreferences()
add_rank!(rankings, agent_a, [agent_1, agent_2, agent_3])
add_rank!(rankings, agent_b, [agent_1, agent_3, agent_2])
add_rank!(rankings, agent_c, [agent_1, agent_3, agent_2])

add_rank!(rankings, agent_1, [agent_c, agent_b, agent_a])
add_rank!(rankings, agent_2, [agent_c, agent_b, agent_a])
add_rank!(rankings, agent_3, [agent_c, agent_a, agent_b])
matching_3 = ttc_with_counter([agent_a, agent_b, agent_c], [agent_1, agent_2, agent_3], rankings)
print(matching_3)
matching_3 = matching_3.results
@test matching_3[agent_a] == agent_2
@test matching_3[agent_b] == agent_1
@test matching_3[agent_c] == agent_1

# test weighted capacity deduction

agent_a = Agent("a")
agent_b = Agent("b")
agent_c = Agent("c", 1, 2)
agent_1 = Agent("1", 2, 1)
agent_2 = Agent("2")
agent_3 = Agent("3")
rankings= RankedPreferences()
add_rank!(rankings, agent_a, [agent_1, agent_2, agent_3])
add_rank!(rankings, agent_b, [agent_1, agent_3, agent_2])
add_rank!(rankings, agent_c, [agent_1, agent_3, agent_2])

add_rank!(rankings, agent_1, [agent_c, agent_b, agent_a])
add_rank!(rankings, agent_2, [agent_c, agent_b, agent_a])
add_rank!(rankings, agent_3, [agent_c, agent_a, agent_b])
matching_3 = ttc_with_counter([agent_a, agent_b, agent_c], [agent_1, agent_2, agent_3], rankings)
print(matching_3)
matching_3 = matching_3.results
@test matching_3[agent_a] == agent_2
@test matching_3[agent_b] == agent_3
@test matching_3[agent_c] == agent_1


text:
using LAR

V = Verts([[1,2] [3,4] [5,6]])
m,n = size(V)
mymodel = Lar()
fieldnames(Lar)
fieldnames(Verts)

V = Float64[[1,2] [3,4] [5,6]]
vs = Verts(V)

mymodel = Lar()
fieldnames(Lar)
fieldnames(Verts)

V = readcsv("/Users/paoluzzi/Documents/dev/mtj/testfile.V",Float64)
EV = readcsv("/Users/paoluzzi/Documents/dev/mtj/testfile.EV",Int32)
EV1 = readcsv("/Users/paoluzzi/Documents/dev/mtj/testfile1.EV")

lar = JSON.parsefile("/Users/paoluzzi/Documents/dev/mtj/larfile.json";
		dicttype=Dict, use_mmap=true)

lar = JSON.parse("""
	{"V":[[0.,0.,0.],[1.,0.,0.],[1.,1.,0.],[0.,1.,0.]],
	"EV":[[1,2],[2,3],[3,4],[4,1]],
	"FV":[[4,1,3,2]]}
""")

larstring = """
	{"V" : [[5,0],[7,1],[9,0],[13,2],[15,4],[17,8],[14,9],[13,10],
	[11,11],[9,10],[5,9],[7,9],[3,8],[0,6],[2,3],[2,1],[8,3],
	[10,2],[13,4],[14,6],[13,7],[12,10],[11,9],[9,7],[7,7],[4,7],
	[2,6],[3,5],[4,2],[6,3],[11,4],[12,6],[12,7],[10,5],[8,5],
	[7,6],[5,5]],

	"FV" : [[0,1,16,28,29],[0,15,28],[1,2,17],[1,16,17,33],[2,3,17],
	[3,4,18,19],[3,17,18,30],[4,5,19],[5,6,19],[6,7,20,21,22,32],
	[6,19,20],[7,8,21],[8,9,21,22],[9,11,23,24],[9,22,23],
	[10,11,24,25],[10,12,25],[12,13,25,26],[13,14,27],[13,26,27],
	[14,15,28],[14,27,28,29,36],[16,29,34],[16,33,34],[17,30,33],
	[18,19,31],[18,30,31],[19,20,31,32],[22,23,32,33],[23,24,34,35],
	[23,33,34],[24,25,27,36],[24,35,36],[25,26,27],[29,34,35],
	[29,35,36],[30,31,32,33],[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]]}
""";

lar = json2larmodel(larstring)
fieldnames(lar)
typeof(lar.Lar)
lar.Lar.VV, lar.Lar.EV, lar.Lar.CV
lar.Lar.FV
map(Array{Int,1}, lar.Lar.FV)

EV1 = readcsv("/Users/paoluzzi/Documents/dev/mtj/testfile1.EV")
EV1 = listOfList(EV1)

model1 = json2larmodel("""
	{"V":[[0.,0.,0.],[1.,0.,0.],[1.,1.,0.],[0.,1.,0.]],
	"EV":[[1,2],[2,3],[3,4],[4,1]],
	"FV":[[4,1,3,2]]}
""");

mod2 = json2larmodel(larstring);

M_1 = cellComplex(model1.Lar.EV);
M_2 = cellComplex(model1.Lar.FV);
print(full(M_1))
print(full(M_2))

cube = p.COLOR(p.RED)(p.CUBE(0.5))

p.VIEW(cube)

V = Float64[[1,2] [3,4] [5,6]]
a = PyObject(V)

verts = PyObject(a[:tolist]())

println("V =\n",V)
println("EV =\n",EV)
println("EV1 =\n",EV1)

methods(larview)

larmodel = p.larCuboids((2,2,2),true);
verts,cells,chainbases = importmodel(larmodel);
larview(verts,cells.EV);
larview(verts,cells.FV);
larview(verts,cells.CV);
larview(verts,cells.VV);

v,cv = p.larCuboids((1,2,1))
viewexploded(v',cv'+1) # input are REQUIRED column-major AND 1-based

v,cv = p.larCuboids((20,20,20));
viewexploded(v',cv'+1)

larmodel = p.larCuboids((2,2,2),full=true);
(v,cc)=larmodel;
(vv,ev,fv,cv) = cc;
cc = p.CAT([PyObject(vv),PyObject(ev),PyObject(fv),PyObject(cv)]);
scaleargs=2;
viewexploded(v',rebase(cc),scaleargs);

mod1 = json2larmodel("""
	{"V":[[0.,0.,0.],[1.,0.,0.],[1.,1.,0.],[0.,1.,0.]],
	"EV":[[1,2],[2,3],[3,4],[4,1]],
	"FV":[[4,1,3,2]]}
""")

larview(mod1.Verts,mod1.Lar.FV)
larview(mod2.Verts,rebase(mod2.Lar.FV))
larview(        embed(mod2.Verts),rebase(mod2.Lar.FV))
viewexploded(mod2.Verts,rebase(mod2.Lar.FV))

larmodel = p.larCuboids((1,1));
(v,fv)=larmodel;
larview(translate(v',[0.5; 0.5]),rebase(fv))

larmodel = p.larCuboids((2,2,2));
(v,cv)=larmodel;
viewexploded(translate(v',[-2; 0; 0]),rebase(cv),1.2)

A = reshape(1:30,3,10)
s = [1,1/2,1/3]
scale( A, s )

larmodel = p.larCuboids((2,2,2));
(v,cv)=larmodel;
viewexploded(scale(v',[-1; -1; -1]),rebase(cv),1.2)

larmodel = p.larCuboids((1,1,1));
(v,cv)=larmodel;
viewexploded(rotate([1,1,1],v'),rebase(cv),1.2);


text:
using Test, DiscreteExteriorCalculus, DiscretePDEs
using LinearAlgebra: norm
using UniqueVectors: UniqueVector

function disks!(r1::Real, r2::Real, factor::Real, N::Int)
    @assert (N == 2) || (N == 3)
    file_name = joinpath(@__DIR__, "disks.geo")
    if N == 2
        geo_write!(file_name, characteristic_length_factor=factor,
            footer="""
            Circle(1) = {0, 0, 0, $r1, 0, 2*Pi};
            Circle(2) = {0, 0, 0, $r2, 0, 2*Pi};
            Curve Loop(1) = {2};
            Curve Loop(2) = {1};
            Plane Surface(1) = {1, 2};
            Physical Curve("Conductor") = {1};
            Physical Surface("Dielectric") = {1};
            """)
    elseif N == 3
        geo_write!(file_name, characteristic_length_factor=factor,
            footer="""
            Sphere(1) = {0, 0, 0, $r1, -Pi/2, Pi/2, 2*Pi};
            Sphere(2) = {0, 0, 0, $r2, -Pi/2, Pi/2, 2*Pi};
            BooleanDifference{ Volume{2}; Delete; }{ Volume{1}; Delete; }
            Physical Surface("Conductor") = {1};
            Physical Volume("Dielectric") = {2};
            """)
    end
    return file_name
end

function voltage(R, outer_radius, ϵ, N)
    @assert (N == 2) || (N == 3)
    if N == 2
        return (log(outer_radius) - log(norm(R)))/(2π * ϵ)
    elseif N == 3
        return (1/norm(R) - 1/outer_radius)/(4π * ϵ)
    end
end

initialize!()
for (N, factor) in [(2, .3), (3, .45)]
    r1, r2 = 1.5, 3
    file_name = disks!(r1, r2, factor, N)
    gmsh_open!(file_name)
    K = N+1
    mesh!(K)

    node_tags, points, tcomp = get_triangulated_complex(N, K)
    group_dict = get_physical_groups(node_tags, points)
    @test typeof(tcomp) <: TriangulatedComplex{N, K}
    comp = tcomp.complex
    orient!(comp)
    sources = [get_charge_source(comp, group_dict["Conductor"])]
    _, exterior = boundary_components_connected(comp)
    m = Metric(N)
    mesh = Mesh(tcomp, circumcenter(m))

    ϵ = 2
    ϵ_form = get_material(comp, group_dict["Dielectric"], ϵ, 2)
    bbox, null_basis = electrostatics_blackbox(m, mesh, sources, exterior, ϵ_form)

    elastance = impedance_matrix(bbox, null_basis)[1,1]
    charge = 2
    φ, q = solve_statics(bbox, null_basis, [charge])
    @testset "concentric $N-disks" begin
        @test isapprox(elastance[1,1], voltage(r1, r2, ϵ, N), rtol= N == 2 ? 1e-4 : 4e-2)
        @test bbox.Y[1] * φ == q
        inner_charges = q[[findfirst(isequal(c), comp.cells[1]) for c in sources[1]]]
        @test sum(inner_charges) ≈ charge
        @test all(abs.(diff(inner_charges)) .< charge * (N == 2 ? .0003 : .005))
        outer_charges = q[[findfirst(isequal(c), comp.cells[1]) for c in exterior.cells[1]]]
        @test sum(outer_charges) ≈ -charge
        @test all(abs.(diff(outer_charges)) .< charge * (N == 2 ? .0003 : .005))
        @test isapprox(φ, charge * voltage.([c.points[1].coords
            for c in comp.cells[1]], r2, ϵ, N), rtol= N == 2 ? 7e-4 : 4e-2)
    end
    if false # set to true to display plots with Gmsh
        ρ = source_density(m, mesh, q, 1)
        vec_E = [N == 2 ? [p..., 0] : p for p in electric_field(m, mesh, φ)]
        comp_points = UniqueVector([c.points[1] for c in comp.cells[1]])
        ordering = [findfirst(isequal(p), comp_points) for p in points]
        add_field!("Electric field", node_tags, vec_E[ordering])
        inds = N == 2 ? (1:length(points)) : findall([p.coords[1] > 0 for p in points])
        add_field!("Electric field magnitude", node_tags[inds], norm.(vec_E)[ordering][inds])
        dd_field!("Voltage", node_tags[inds], φ[ordering][inds])
        add_field!("Charge density", node_tags[inds], ρ[ordering][inds])
        gui!()
    end
    rm(file_name)
end


text:
module Borg

include("BorgC.jl")

export createProblem, solve

function eval_wrapper(vars_raw::Ptr{Float64}, objs_raw::Ptr{Float64}, consts_raw::Ptr{Float64}, userParams::Ptr{Void})
    problem = unsafe_pointer_to_objref(userParams)::BorgProblem

    vars = pointer_to_array(vars_raw, problem.numberOfVariables)
    objs = pointer_to_array(objs_raw, problem.numberOfObjectives)
    consts = pointer_to_array(consts_raw, problem.numberOfConstraints)

    problem.eval_fn(vars, objs, consts)

    return
end

type BorgProblem
	ref::BorgC.BORG_Problem

	numberOfVariables::Int64
	numberOfObjectives::Int64
	numberOfConstraints::Int64

	eval_fn::Function

	function BorgProblem(
            numberOfVariables::Int,
            numberOfObjectives::Int,
            numberOfConstraints::Int,
            eval_fn::Function)

		problem = new()

        problem.numberOfVariables = numberOfVariables
        problem.numberOfObjectives = numberOfObjectives
        problem.numberOfConstraints = numberOfConstraints
        problem.eval_fn = eval_fn
		problem.ref = BorgC.BORG_Problem_create(numberOfVariables, numberOfObjectives, numberOfConstraints, eval_wrapper, problem)

        return problem
	end
end

function createProblem(
    varLowerBounds::Vector{Float64},
    varUpperBounds::Vector{Float64},
    objEpsilons::Vector{Float64},
    numberOfConstraints::Int,
    eval_fn::Function)

    if length(varLowerBounds) != length(varUpperBounds)
        error("varLowerBounds and varUpperBounds must be of same length.")
    end

    numberOfVariables = length(varLowerBounds)
    numberOfObjectives = length(objEpsilons)

    problem = BorgProblem(numberOfVariables, numberOfObjectives, numberOfConstraints, eval_fn)

    for i in 1:numberOfVariables
        BorgC.BORG_Problem_set_bounds(problem.ref, i, varLowerBounds[i], varUpperBounds[i])
    end

    for i in 1:numberOfObjectives
        BorgC.BORG_Problem_set_epsilon(problem.ref, i, objEpsilons[i])
    end

    return problem
end

function solve(problem::BorgProblem, maxIt::Int)
    archive = BorgC.BORG_Algorithm_run(problem.ref, maxIt)
    return archive
end

end


text:
@testset "Simplicial Complex" begin
    cplx = SimplicialComplex()
    @test eltype(cplx) == Simplex
    @test length(cells(cplx)) == 0
    @test length(cells(cplx,2)) == 0
    @test size(cplx) == ()

    cplx = SimplicialComplex(Simplex(1,2,3), Simplex(2,4), Simplex(4))

    @test eltype(cplx) == Simplex
    @test length(cplx) == 9
    @test length(cells(cplx)) == 3
    @test length(cells(cplx,2)) == 1
    @test length(cells(cplx,3)) == 0
    @test size(cplx) == (4,4,1)
    @test size(cplx,-1) == 0
    @test size(cplx,0) == 4
    @test size(cplx,1) == 4
    @test size(cplx,2) == 1
    @test size(cplx,3) == 0
    @test Simplex(4) ∈ cplx
    @test ComputationalHomology.position(cplx, Simplex(4)) == 4
    @test Simplex(10) ∉ cplx

    @test first(Iterators.drop(cplx, 5)) == Simplex(1,3)

    cplx = SimplicialComplex(Simplex(1,2,3),
                       Simplex(2,4),
                       Simplex(3,4),
                       Simplex(5,4),
                       Simplex(6))

    # (co)faces
    fs = [Simplex(4), Simplex(3)]
    @testset "Faces" for fi in faces(cplx, hash(Simplex(3,4)))
        @test cplx[fi,0] ∈ fs
        @test fi ∈ map(hash, fs)
    end
    @test length(faces(cplx, hash(Simplex(1)))) == 0
    fs = [Simplex(2,4), Simplex(3,4), Simplex(5,4)]
    @testset "Cofaces" for fi in cofaces(cplx, hash(Simplex(4)))
        @test cplx[fi,1] ∈ fs
        @test fi ∈ map(hash, fs)
    end
    @test length(cofaces(cplx, hash(Simplex(1,2,3)))) == 0

    # boundary matrix
    @test convert(Matrix, boundary(cplx, 2, Int)) == [1 -1 1 0 0 0]'

    # simplex values
    cplx = SimplicialComplex()

    splx = push!(cplx, Simplex(1))
    @test ComputationalHomology.position(cplx, splx[1]) == 1
    @test dim(splx[1]) == 0
    @test sum(size(cplx)) == 1
    @test indices(cells(cplx, 0)[1]) == [1]

    splx = push!(cplx, Simplex(1,2))
    @test splx[1] == Simplex(1,2)
    @test dim(splx[1]) == 1
    @test sum(size(cplx)) == 2
    @test sort!(indices(cells(cplx, 1)[1])) == [1,2]

    cplx = SimplicialComplex()
    splxs = push!(cplx, Simplex(1, 2), recursive=true)
    @test splxs[1] == Simplex(1,2)
    @test sort!(indices(splxs[1])) == [1, 2]
    @test splxs[3] == Simplex(2)
    @test indices(splxs[3]) == [2]
    @test splxs[2] == Simplex(1)
    @test indices(splxs[2]) == [1]
    @test sum(size(cplx)) == 3
    @test size(cplx, 0) == 2
    @test size(cplx, 1) == 1

    splxs = push!(cplx, Simplex(1,3), recursive=true)
    @test splxs[1] == Simplex(1,3)
    @test sort!(indices(splxs[1])) == [1,3]
    @test splxs[2] == Simplex(3)
    @test indices(splxs[2]) == [3]
    @test sum(size(cplx)) == 5
    @test size(cplx, 0) == 3
    @test size(cplx, 1) == 2

    @test Simplex(1,3) in cplx
    @test first(cells(cplx,1)) in cplx
    @test !(Simplex(2,3) in cplx)
    @test !(Simplex(2,3,3) in cplx)

    io = IOBuffer()
    write(io, cplx)
    seekstart(io)
    cplx2 = read(io, SimplicialComplex(), Simplex{0,Int})
    @test size(cplx2, 0) == 3
    @test size(cplx2, 1) == 2
    @test all(cells(cplx2,1) .== cells(cplx,1))

    A = ComputationalHomology.adjacency_matrix(ComputationalHomology.sphere(1), UInt8)
    @test collect(A) == UInt8[0x00 0x01; 0x01 0x00]
    @test eltype(A) == UInt8
end


text:
@assert supertype(Int32) == Signed
@assert supertype(Signed) == Integer
@assert supertype(Integer) == Real
@assert supertype(Symbol) == Any
@assert Set(subtypes(Type)) == Set([DataType, TypeConstructor, Union])

@assert Float64 <: Real
@assert isa(subtypes(Type), Array)
@assert isa(Array, Type)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct RouteFilterRulePropertiesFormat <: SwaggerModel
    access::Any # spec type: Union{ Nothing, String } # spec name: access
    routeFilterRuleType::Any # spec type: Union{ Nothing, String } # spec name: routeFilterRuleType
    communities::Any # spec type: Union{ Nothing, Vector{String} } # spec name: communities
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function RouteFilterRulePropertiesFormat(;access=nothing, routeFilterRuleType=nothing, communities=nothing, provisioningState=nothing)
        o = new()
        validate_property(RouteFilterRulePropertiesFormat, Symbol("access"), access)
        setfield!(o, Symbol("access"), access)
        validate_property(RouteFilterRulePropertiesFormat, Symbol("routeFilterRuleType"), routeFilterRuleType)
        setfield!(o, Symbol("routeFilterRuleType"), routeFilterRuleType)
        validate_property(RouteFilterRulePropertiesFormat, Symbol("communities"), communities)
        setfield!(o, Symbol("communities"), communities)
        validate_property(RouteFilterRulePropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type RouteFilterRulePropertiesFormat

const _property_map_RouteFilterRulePropertiesFormat = Dict{Symbol,Symbol}(Symbol("access")=>Symbol("access"), Symbol("routeFilterRuleType")=>Symbol("routeFilterRuleType"), Symbol("communities")=>Symbol("communities"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_RouteFilterRulePropertiesFormat = Dict{Symbol,String}(Symbol("access")=>"String", Symbol("routeFilterRuleType")=>"String", Symbol("communities")=>"Vector{String}", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ RouteFilterRulePropertiesFormat }) = collect(keys(_property_map_RouteFilterRulePropertiesFormat))
Swagger.property_type(::Type{ RouteFilterRulePropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_RouteFilterRulePropertiesFormat[name]))}
Swagger.field_name(::Type{ RouteFilterRulePropertiesFormat }, property_name::Symbol) =  _property_map_RouteFilterRulePropertiesFormat[property_name]

const _allowed_RouteFilterRulePropertiesFormat_access = ["Allow", "Deny"]

const _allowed_RouteFilterRulePropertiesFormat_routeFilterRuleType = ["Community"]

function check_required(o::RouteFilterRulePropertiesFormat)
    (getproperty(o, Symbol("communities")) === nothing) && (return false)
    true
end

function validate_property(::Type{ RouteFilterRulePropertiesFormat }, name::Symbol, val)
    if name === Symbol("access")
        Swagger.validate_param(name, "RouteFilterRulePropertiesFormat", :enum, val, _allowed_RouteFilterRulePropertiesFormat_access)
    end
    if name === Symbol("routeFilterRuleType")
        Swagger.validate_param(name, "RouteFilterRulePropertiesFormat", :enum, val, _allowed_RouteFilterRulePropertiesFormat_routeFilterRuleType)
    end
end


text:
using GraphVisualize
using LightGraphs
using Base.Test

g = WheelGraph(10)
plot(g, observe=true)

max_stars_repo_name: JuliaGraphs/GraphVisualization

max_issues_repo_name: JuliaGraphs/GraphVisualization

max_forks_repo_name: JuliaGraphs/GraphVisualization


text:
module NewtonsMethod

using Distributions, LinearAlgebra, Statistics

function newtonroot(f, f_prime; x_0=0.5, tolerance=1E-7, maxiter=1000)
    normdiff=Inf
    iter=1
    x_old=x_0

    while normdiff > tolerance && iter <= maxiter
        x_new=x_old-f(x_old)/f_prime(x_old)
        normdiff = norm(x_new - x_old)
        x_old=x_new
        iter+=1
    end

    if iter <= maxiter
        return (value = x_old, normdiff=normdiff, iter=iter)
    else
        return (value =nothing, normdiff=normdiff, iter=iter)
    end
end

export(newtonroot)

using ForwardDiff

function newtonroot1(f; x_0=0.5, tolerance=1E-7, maxiter=1000)

    D(f) = x -> ForwardDiff.derivative(f, x)
    f_autoprime = D(f)

    normdiff=Inf
    iter=1
    x_old=x_0

    while normdiff > tolerance && iter <= maxiter
        x_new=x_old-f(x_old)/f_prime(x_old)
        normdiff = norm(x_new - x_old)
        x_old=x_new
        iter+=1
    end

    if iter <= maxiter
        return (value = x_old, normdiff=normdiff, iter=iter)
    else
        return (value =nothing, normdiff=normdiff, iter=iter)
    end

end

export(newtonroot1)

end # module


text:
function EstimatedMeasurement(arg0::ObservedMeasurement, arg1::jint, arg2::jint, arg3::Vector{SpacecraftState}, arg4::Vector{TimeStampedPVCoordinates})
    return EstimatedMeasurement((ObservedMeasurement, jint, jint, Vector{SpacecraftState}, Vector{TimeStampedPVCoordinates}), arg0, arg1, arg2, arg3, arg4)
end

function get_count(obj::EstimatedMeasurement)
    return jcall(obj, "getCount", jint, ())
end

function get_date(obj::EstimatedMeasurement)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_derivatives_drivers(obj::EstimatedMeasurement)
    return jcall(obj, "getDerivativesDrivers", Stream, ())
end

function get_estimated_value(obj::EstimatedMeasurement)
    return jcall(obj, "getEstimatedValue", Vector{jdouble}, ())
end

function get_iteration(obj::EstimatedMeasurement)
    return jcall(obj, "getIteration", jint, ())
end

function get_observed_measurement(obj::EstimatedMeasurement)
    return jcall(obj, "getObservedMeasurement", ObservedMeasurement, ())
end

function get_observed_value(obj::EstimatedMeasurement)
    return jcall(obj, "getObservedValue", Vector{jdouble}, ())
end

function get_parameter_derivatives(obj::EstimatedMeasurement, arg0::ParameterDriver)
    return jcall(obj, "getParameterDerivatives", Vector{jdouble}, (ParameterDriver,), arg0)
end

function get_participants(obj::EstimatedMeasurement)
    return jcall(obj, "getParticipants", Vector{TimeStampedPVCoordinates}, ())
end

function get_state_derivatives(obj::EstimatedMeasurement, arg0::jint)
    return jcall(obj, "getStateDerivatives", Vector{Vector{jdouble}}, (jint,), arg0)
end

function get_state_size(obj::EstimatedMeasurement)
    return jcall(obj, "getStateSize", jint, ())
end

function get_states(obj::EstimatedMeasurement)
    return jcall(obj, "getStates", Vector{SpacecraftState}, ())
end

function get_status(obj::EstimatedMeasurement)
    return jcall(obj, "getStatus", EstimatedMeasurement_Status, ())
end

function get_time_offset(obj::EstimatedMeasurement)
    return jcall(obj, "getTimeOffset", jdouble, ())
end

function set_estimated_value(obj::EstimatedMeasurement, arg0::Vector{jdouble})
    return jcall(obj, "setEstimatedValue", void, (Vector{jdouble},), arg0)
end

function set_parameter_derivatives(obj::EstimatedMeasurement, arg0::ParameterDriver, arg1::Vector{jdouble})
    return jcall(obj, "setParameterDerivatives", void, (ParameterDriver, Vector{jdouble}), arg0, arg1)
end

function set_state_derivatives(obj::EstimatedMeasurement, arg0::jint, arg1::Vector{Vector{jdouble}})
    return jcall(obj, "setStateDerivatives", void, (jint, Vector{Vector{jdouble}}), arg0, arg1)
end

function set_status(obj::EstimatedMeasurement, arg0::EstimatedMeasurement_Status)
    return jcall(obj, "setStatus", void, (EstimatedMeasurement_Status,), arg0)
end


text:
function find_std_headers()
    headers = String[]
    @static if Sys.isapple()
        xcode_path = strip(read(`xcode-select --print-path`, String))
        sdk_path = strip(read(`xcrun --show-sdk-path`, String))
        occursin("Xcode", xcode_path) &&
            (xcode_path *= "/Toolchains/XcodeDefault.xctoolchain/")
        didfind = false
        lib = joinpath(xcode_path, "usr", "lib", "c++", "v1")
        inc = joinpath(xcode_path, "usr", "include", "c++", "v1")
        isdir(lib) && (push!(headers, lib); didfind = true)
        isdir(inc) && (push!(headers, inc); didfind = true)
        if isdir("/usr/include")
            push!(headers, "/usr/include")
        else
            isdir(joinpath(sdk_path, "usr", "include"))
            push!(headers, joinpath(sdk_path, "usr", "include"))
        end
        didfind ||
            error("Could not find C++ standard library. Is XCode or CommandLineTools installed?")
    end
    return headers
end

"""
    find_dependent_headers(headers::Vector{T}, args::Vector, general_ops::Dict) where {T<:AbstractString}
Return a vector of headers to which those missing dependent headers are added.
"""
function find_dependent_headers(headers::Vector{T}, args::Vector, general_ops::Dict) where {T<:AbstractString}
    blacklist = get(general_ops, "header_blacklist", [])
    flags = CXTranslationUnit_DetailedPreprocessingRecord
    flags |= CXTranslationUnit_SkipFunctionBodies
    all_headers = copy(headers)
    dependent_headers = T[]
    new_headers = T[]
    blocked_headers = T[]
    idx = Index(false)
    while true
        empty!(new_headers)
        for header in all_headers
            tu = parse_header(idx, header, args, flags)
            GC.@preserve tu begin
                tu_cursor = getTranslationUnitCursor(tu)
                header_name = spelling(tu_cursor)
                header_dir = dirname(header_name)
                for cursor in children(tu_cursor)
                    is_inclusion_directive(cursor) || continue
                    file = getIncludedFile(cursor)
                    file_name = get_filename(file) |> normpath
                    (isempty(file_name) || !isfile(file_name)) && continue

                    dir = dirname(file_name)
                    file_name ∈ all_headers && continue
                    if startswith(header_dir, dir) || startswith(dir, header_dir)
                        isempty(header_dir) && continue
                        file_name ∈ new_headers && continue
                        if any(x->!isempty(x) && endswith(file_name, x), blacklist)
                            file_name ∉ blocked_headers && push!(blocked_headers, file_name)
                            continue
                        end
                        push!(new_headers, file_name)
                    end
                end
            end
        end
        isempty(new_headers) && break
        foreach(new_headers) do h
            @info "found dependent header: $h"
        end
        append!(dependent_headers, new_headers)
        append!(all_headers, new_headers)
    end

    for file in blocked_headers
        @info "skipped a dependent file: $file because it's in the blacklist."
    end

    return normpath.(dependent_headers)
end

"""
    create_context(headers::Vector, args::Vector=String[], options::Dict=Dict())
Create a context from a vector of paths of headers, a vector of compiler flags and
a option dict.
"""
function create_context(headers::Vector, args::Vector=String[], options::Dict=Dict())
    ctx = Context(options)

    general_options = get(options, "general", Dict())
    if get(general_options, "auto_detect_system_headers", true)
        sys_headers = "-I" .* find_std_headers()
        args = vcat(sys_headers, args)
    end

    if get(general_options, "use_clang_headers", true)
        clang_incs = ["-I"*CLANG_INCLUDE]
        args = vcat(clang_incs, args)
    end

    dependent_headers = find_dependent_headers(headers, args, general_options)

    parse_headers!(ctx, headers, args)

    push!(ctx.passes, CollectTopLevelNode(ctx.trans_units, dependent_headers))
    push!(ctx.passes, LinkTypedefToAnonymousTagType())
    push!(ctx.passes, IndexDefinition())
    push!(ctx.passes, CollectNestedRecord())
    push!(ctx.passes, FindOpaques())
    push!(ctx.passes, ResolveDependency())
    push!(ctx.passes, RemoveCircularReference())
    push!(ctx.passes, TopologicalSort())
    push!(ctx.passes, IndexDefinition())
    push!(ctx.passes, ResolveDependency())
    push!(ctx.passes, CatchDuplicatedAnonymousTags())
    push!(ctx.passes, CodegenPreprocessing())
    if get(general_options, "smart_de_anonymize", true)
        push!(ctx.passes, DeAnonymize())
    end
    push!(ctx.passes, Audit())
    push!(ctx.passes, Codegen())
    push!(ctx.passes, CodegenMacro())
    # push!(ctx.passes, CodegenPostprocessing())
    if get(general_options, "auto_mutability", false)
        push!(ctx.passes, TweakMutability())
    end

    # support old behavior
    api_file = get(general_options, "output_api_file_path", "")
    common_file = get(general_options, "output_common_file_path", "")

    output_file_path = get(general_options, "output_file_path", "")
    if isempty(api_file) && isempty(common_file)
        if !isempty(output_file_path)
            push!(ctx.passes, ProloguePrinter(output_file_path))
            push!(ctx.passes, GeneralPrinter(output_file_path))
            push!(ctx.passes, EpiloguePrinter(output_file_path))
        else
            # print to stdout if there is no `output_file_path`
            # this is handy when playing in REPL
            push!(ctx.passes, StdPrinter())
        end
    else
        # TODO: impl
    end

    return ctx
end


text:
#!usr/bin/env julia

include("chemical_parse.jl")

str = "CH4 K3[Fe(CN)6] C(NH2)3CH(CH3)2"
str = "CH4 O2 CO2 H2O"
str = "HCl Na2CO3 NaCl H2O CO2"

function lindep{T}(args::Vector{T}...)
    length(args) > length(args[1]) && false
    rank(hcat(args...)) == length(args)
end

## function chemical_solve(str::AbstractString)
M, row, column = chemical_parse(str)
rk = rank(M)
if rk == length(row)
    return "No solution."
else

end
## end

max_stars_repo_name: sclereid/collection

max_issues_repo_name: sclereid/collection

max_forks_repo_name: sclereid/collection


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IpsecPolicy <: SwaggerModel
    saLifeTimeSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: saLifeTimeSeconds
    saDataSizeKilobytes::Any # spec type: Union{ Nothing, Int32 } # spec name: saDataSizeKilobytes
    ipsecEncryption::Any # spec type: Union{ Nothing, String } # spec name: ipsecEncryption
    ipsecIntegrity::Any # spec type: Union{ Nothing, String } # spec name: ipsecIntegrity
    ikeEncryption::Any # spec type: Union{ Nothing, String } # spec name: ikeEncryption
    ikeIntegrity::Any # spec type: Union{ Nothing, String } # spec name: ikeIntegrity
    dhGroup::Any # spec type: Union{ Nothing, String } # spec name: dhGroup
    pfsGroup::Any # spec type: Union{ Nothing, String } # spec name: pfsGroup

    function IpsecPolicy(;saLifeTimeSeconds=nothing, saDataSizeKilobytes=nothing, ipsecEncryption=nothing, ipsecIntegrity=nothing, ikeEncryption=nothing, ikeIntegrity=nothing, dhGroup=nothing, pfsGroup=nothing)
        o = new()
        validate_property(IpsecPolicy, Symbol("saLifeTimeSeconds"), saLifeTimeSeconds)
        setfield!(o, Symbol("saLifeTimeSeconds"), saLifeTimeSeconds)
        validate_property(IpsecPolicy, Symbol("saDataSizeKilobytes"), saDataSizeKilobytes)
        setfield!(o, Symbol("saDataSizeKilobytes"), saDataSizeKilobytes)
        validate_property(IpsecPolicy, Symbol("ipsecEncryption"), ipsecEncryption)
        setfield!(o, Symbol("ipsecEncryption"), ipsecEncryption)
        validate_property(IpsecPolicy, Symbol("ipsecIntegrity"), ipsecIntegrity)
        setfield!(o, Symbol("ipsecIntegrity"), ipsecIntegrity)
        validate_property(IpsecPolicy, Symbol("ikeEncryption"), ikeEncryption)
        setfield!(o, Symbol("ikeEncryption"), ikeEncryption)
        validate_property(IpsecPolicy, Symbol("ikeIntegrity"), ikeIntegrity)
        setfield!(o, Symbol("ikeIntegrity"), ikeIntegrity)
        validate_property(IpsecPolicy, Symbol("dhGroup"), dhGroup)
        setfield!(o, Symbol("dhGroup"), dhGroup)
        validate_property(IpsecPolicy, Symbol("pfsGroup"), pfsGroup)
        setfield!(o, Symbol("pfsGroup"), pfsGroup)
        o
    end
end # type IpsecPolicy

const _property_map_IpsecPolicy = Dict{Symbol,Symbol}(Symbol("saLifeTimeSeconds")=>Symbol("saLifeTimeSeconds"), Symbol("saDataSizeKilobytes")=>Symbol("saDataSizeKilobytes"), Symbol("ipsecEncryption")=>Symbol("ipsecEncryption"), Symbol("ipsecIntegrity")=>Symbol("ipsecIntegrity"), Symbol("ikeEncryption")=>Symbol("ikeEncryption"), Symbol("ikeIntegrity")=>Symbol("ikeIntegrity"), Symbol("dhGroup")=>Symbol("dhGroup"), Symbol("pfsGroup")=>Symbol("pfsGroup"))
const _property_types_IpsecPolicy = Dict{Symbol,String}(Symbol("saLifeTimeSeconds")=>"Int32", Symbol("saDataSizeKilobytes")=>"Int32", Symbol("ipsecEncryption")=>"String", Symbol("ipsecIntegrity")=>"String", Symbol("ikeEncryption")=>"String", Symbol("ikeIntegrity")=>"String", Symbol("dhGroup")=>"String", Symbol("pfsGroup")=>"String")
Base.propertynames(::Type{ IpsecPolicy }) = collect(keys(_property_map_IpsecPolicy))
Swagger.property_type(::Type{ IpsecPolicy }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IpsecPolicy[name]))}
Swagger.field_name(::Type{ IpsecPolicy }, property_name::Symbol) =  _property_map_IpsecPolicy[property_name]

const _allowed_IpsecPolicy_ipsecEncryption = ["None", "DES", "DES3", "AES128", "AES192", "AES256", "GCMAES128", "GCMAES192", "GCMAES256"]

const _allowed_IpsecPolicy_ipsecIntegrity = ["MD5", "SHA1", "SHA256", "GCMAES128", "GCMAES192", "GCMAES256"]

const _allowed_IpsecPolicy_ikeEncryption = ["DES", "DES3", "AES128", "AES192", "AES256"]

const _allowed_IpsecPolicy_ikeIntegrity = ["MD5", "SHA1", "SHA256", "SHA384"]

const _allowed_IpsecPolicy_dhGroup = ["None", "DHGroup1", "DHGroup2", "DHGroup14", "DHGroup2048", "ECP256", "ECP384", "DHGroup24"]

const _allowed_IpsecPolicy_pfsGroup = ["None", "PFS1", "PFS2", "PFS2048", "ECP256", "ECP384", "PFS24"]

function check_required(o::IpsecPolicy)
    (getproperty(o, Symbol("saLifeTimeSeconds")) === nothing) && (return false)
    (getproperty(o, Symbol("saDataSizeKilobytes")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IpsecPolicy }, name::Symbol, val)
    if name === Symbol("ipsecEncryption")
        Swagger.validate_param(name, "IpsecPolicy", :enum, val, _allowed_IpsecPolicy_ipsecEncryption)
    end
    if name === Symbol("ipsecIntegrity")
        Swagger.validate_param(name, "IpsecPolicy", :enum, val, _allowed_IpsecPolicy_ipsecIntegrity)
    end
    if name === Symbol("ikeEncryption")
        Swagger.validate_param(name, "IpsecPolicy", :enum, val, _allowed_IpsecPolicy_ikeEncryption)
    end
    if name === Symbol("ikeIntegrity")
        Swagger.validate_param(name, "IpsecPolicy", :enum, val, _allowed_IpsecPolicy_ikeIntegrity)
    end
    if name === Symbol("dhGroup")
        Swagger.validate_param(name, "IpsecPolicy", :enum, val, _allowed_IpsecPolicy_dhGroup)
    end
    if name === Symbol("pfsGroup")
        Swagger.validate_param(name, "IpsecPolicy", :enum, val, _allowed_IpsecPolicy_pfsGroup)
    end
end


text:
###################################
### Julia OpenStreetMap Package ###
### MIT License                 ###
### Copyright 2014              ###
###################################

module OpenStreetMap

using Reexport
@reexport using Geodesy
using LightXML
using LibExpat
using Winston
using Graphs
using Compat

export parseMapXML, getOSMData, getBounds
export plotMap, cropMap!
export findIntersections, nearestNode, segmentHighways, highwaySegments
export roadways, walkways, cycleways, classify
export createGraph, shortestRoute, fastestRoute, routeEdges
export nodesWithinRange, nodesWithinDrivingDistance, nodesWithinDrivingTime
export findHighwaySets, findIntersectionClusters, replaceHighwayNodes!
export simCityGrid

include("types.jl")
include("classes.jl")
include("layers.jl")
include("speeds.jl")

include("parseMap.jl")
include("nodes.jl")
include("bounds.jl")
include("highways.jl")
include("features.jl")
include("buildings.jl")

include("crop.jl")
include("plot.jl")
include("intersections.jl")
include("transforms.jl")
include("routing.jl")

include("simulate.jl")

include("deprecated.jl")

end # module OpenStreetMap


text:
## Multiply 2 matrices
# checks for invalid dimensions
# uses 1st method (element dot product method)

function outer(A,B)
  return A * B'
end

function matmul(A,B)

# convert any 1D column vectors into 1*x 2D matrices
  if ndims(A) = 1
    A = A[:,:]
  end

  if ndims(B) = 1
    B = B[:,:]
  end

# define row and column sizes
# for each matrix
  rowA, colA = size(A)
  rowB, colB = size(B)

# check that dimensions are valid for multiplication, otherwise throw error
  if colA != rowB
    throw(DimensionMismatch("Matrices of dimensions $(size(A)) and $(size(B)) cannot be multiplied together"))
  end

# initialise zero matrix to hold result
# matrix size is the outer dimensions of the two matrices
  result = zeros(rowA, colB)

# calculate A*B, populate zero matrix with results
# uses the outer products of each row with its corresponding column
# and adds the results together
  for i in 1:rowB # could also use colA, as colA = rowB
      result += outer(A[:, i], B[i, :])
  end

# return the result
  return result
end

max_stars_repo_name: jawuku/julia

max_issues_repo_name: jawuku/julia

max_forks_repo_name: jawuku/julia


text:
## ------------------------------------------------------------------
genlabel(prefix::String, n::Int = 8) = string(prefix, randstring(n))
genlabel(n::Int = 8) = randstring(n)

macro genlabel!()

    macroreg = _macro_call_regex("genlabel!")
    label = genlabel()
    rep = string("\"", genlabel(), "\"")
    _replace_line(__source__, macroreg, rep)
    return :($label)
end


text:
Random.seed!(123)

rowvec(y) = y
rowvec(y::Vector) = reshape(y, 1, length(y))

@test MLJFlux.MLJModelInterface.istransparent(Flux.ADAM(0.1))

@testset "nrows" begin
    Xmatrix = rand(10, 3)
    X = MLJBase.table(Xmatrix)
    @test MLJFlux.nrows(X) == 10
    @test MLJFlux.nrows(Tables.columntable(X)) == 10
end

@testset "collate" begin
    # NeuralNetworRegressor:
    Xmatrix = broadcast(x->round(x, sigdigits=2), rand(10, 3))
    # convert to a column table:
    X = MLJBase.table(Xmatrix)

    y = rand(10)
    model = MLJFlux.NeuralNetworkRegressor()
    model.batch_size= 3
    @test MLJFlux.collate(model, X, y) ==
        ([Xmatrix'[:,1:3], Xmatrix'[:,4:6], Xmatrix'[:,7:9], Xmatrix'[:,10:10]],
         rowvec.([y[1:3], y[4:6], y[7:9], y[10:10]]))

    # NeuralNetworClassifier:
    y = categorical(['a', 'b', 'a', 'a', 'b', 'a', 'a', 'a', 'b', 'a'])
    model = MLJFlux.NeuralNetworkClassifier()
    model.batch_size = 3
    data = MLJFlux.collate(model, X, y)

    @test data == ([Xmatrix'[:,1:3], Xmatrix'[:,4:6],
                    Xmatrix'[:,7:9], Xmatrix'[:,10:10]],
                   rowvec.([[1 0 1; 0 1 0], [1 0 1; 0 1 0], [1 1 0; 0 0 1],
                            reshape([1; 0], (2,1))]))

    # MultitargetNeuralNetworRegressor:
    ymatrix = rand(10, 2)
    y = MLJBase.table(ymatrix) # a rowaccess table
    model = MLJFlux.NeuralNetworkRegressor()
    model.batch_size= 3
    @test MLJFlux.collate(model, X, y) ==
        ([Xmatrix'[:,1:3], Xmatrix'[:,4:6], Xmatrix'[:,7:9], Xmatrix'[:,10:10]],
         rowvec.([ymatrix'[:,1:3], ymatrix'[:,4:6], ymatrix'[:,7:9],
                  ymatrix'[:,10:10]]))

    y = Tables.columntable(y) # try a columnaccess table
    @test MLJFlux.collate(model, X, y) ==
        ([Xmatrix'[:,1:3], Xmatrix'[:,4:6], Xmatrix'[:,7:9], Xmatrix'[:,10:10]],
         rowvec.([ymatrix'[:,1:3], ymatrix'[:,4:6],
                  ymatrix'[:,7:9], ymatrix'[:,10:10]]))

    # ImageClassifier
    Xmatrix = coerce(rand(6, 6, 1, 10), GrayImage)
    y = categorical(['a', 'b', 'a', 'a', 'b', 'a', 'a', 'a', 'b', 'a'])
    model = MLJFlux.ImageClassifier(batch_size=2)

    data = MLJFlux.collate(model, Xmatrix, y)
    @test  first.(data) == (Float32.(cat(Xmatrix[1], Xmatrix[2], dims=4)),
                            rowvec.([1 0;0 1]))

    expected_y = [[1 0;0 1], [1 1;0 0], [0 1; 1 0], [1 1;0 0], [0 1; 1 0]]
    for i=1:5
        @test Int.(last(data)[i]) == expected_y[i]
    end

end

Xmatrix = rand(100, 5)
X = MLJBase.table(Xmatrix)
y = Xmatrix[:, 1] + Xmatrix[:, 2] + Xmatrix[:, 3] +
    Xmatrix[:, 4] + Xmatrix[:, 5]

data = [(Xmatrix'[:,1:20], y[1:20]),
        (Xmatrix'[:,21:40], y[21:40]),
        (Xmatrix'[:,41:60], y[41:60]),
        (Xmatrix'[:,61:80], y[61:80]),
        (Xmatrix'[:, 81:100], y[81:100])]

data = ([Xmatrix'[:,1:20], Xmatrix'[:,21:40], Xmatrix'[:,41:60],
         Xmatrix'[:,61:80], Xmatrix'[:,81:100]],
            rowvec.([y[1:20], y[21:40], y[41:60], y[61:80], y[81:100]]))

# construct two chains with identical state, except one has
# dropout and the other does not:
chain_yes_drop = Flux.Chain(Flux.Dense(5, 15),
                            Flux.Dropout(0.2),
                            Flux.Dense(15, 8),
                            Flux.Dense(8, 1))

chain_no_drop = deepcopy(chain_yes_drop)
chain_no_drop.layers[2].p = 1.0

test_input = rand(Float32, 5, 1)

# check both chains have same behaviour before training:
@test chain_yes_drop(test_input) == chain_no_drop(test_input)

epochs = 10

@testset_accelerated "fit! and dropout" accel begin

    move = MLJFlux.Mover(accel)

    Random.seed!(123)

    _chain_yes_drop, history = MLJFlux.fit!(chain_yes_drop,
                                            Flux.Optimise.ADAM(0.001),
                                            Flux.mse,
                                            epochs,
                                            0,
                                            0,
                                            0,
                                            accel,
                                            data[1],
                                            data[2])

    println()

    Random.seed!(123)

    _chain_no_drop, history = MLJFlux.fit!(chain_no_drop,
                                           Flux.Optimise.ADAM(0.001),
                                           Flux.mse,
                                           epochs,
                                           0,
                                           0,
                                           0,
                                           accel,
                                           data[1],
                                           data[2])

    # check chains have different behaviour after training:
    @test !(_chain_yes_drop(test_input) ≈ _chain_no_drop(test_input))

    # check chain with dropout is deterministic outside of training
    # (if we do not differentiate):
    @test all(_chain_yes_drop(test_input) ==
              _chain_yes_drop(test_input) for i in 1:1000)

    @test length(history) == epochs + 1

end

# to control chain initialization:
myinit(n, m) = reshape(float(1:n*m), n , m)

mutable struct TESTBuilder <: MLJFlux.Builder end
MLJFlux.build(builder::TESTBuilder, n_in, n_out) =
    Flux.Chain(Flux.Dense(n_in, n_out, init=myinit))

@testset_accelerated "issue #152" accel begin

    # data:
    n = 100
    d = 5
    Xmat = rand(Float64, n, d)
    X = MLJBase.table(Xmat);
    y = X.x1 .^2 + X.x2 .* X.x3 - 4 * X.x4

    # train a model on all the data using batch size > 1:
    model = MLJFlux.NeuralNetworkRegressor(builder = TESTBuilder(),
                                   batch_size=25,
                                   epochs=1,
                                   loss=Flux.mse)
    mach = machine(model, X, y)
    fit!(mach, verbosity=0)

    # extract the pre-training loss computed in the `fit!(chain, ...)` method:
    pretraining_loss = report(mach).training_losses[1]

    # compute by hand:
    chain0 = myinit(1, d)
    pretraining_yhat = Xmat*chain0' |> vec
    @test y isa Vector && pretraining_yhat isa Vector
    pretraining_loss_by_hand =  MLJBase.l2(pretraining_yhat, y) |> mean

    # compare:
    @test pretraining_loss ≈ pretraining_loss_by_hand

end