text:
abstract type NeuroPath{D,B} end

struct StudyPath{D} <: NeuroPath{D,String}
    dirname::D
    study::String

    StudyPath(d::Union{Nothing,String}, s::String) = new{typeof(d)}(d, s)
    StudyPath(d, s::AbstractString) = StudyPath(d, String(s))
    StudyPath(s) = StudyPath(nothing, s)
end

struct DataPath{D} <: NeuroPath{D,String}
    dirname::D
    basename::String

    function DataPath(d::Union{Nothing,StudyPath}, s::String)
        _check_data_name(s)
        new{typeof(d)}(d, s)
    end
    DataPath(d::Union{Nothing,StudyPath}, s) = DataPath(d, String(s))
    DataPath(s) = DataPath(nothing, s)
end

struct DerivativePath{D} <: NeuroPath{D,String}
    dirname::D
    pipeline::String

    DerivativePath(::Nothing, s::String) = new{Nothing}(nothing, s)
    function DerivativePath(d::DataPath, s::String)
        @assert basename(d) === "derivatives"
        return new{typeof(d)}(d, s)
    end
    DerivativePath(d, s::AbstractString) = DerivativePath(d, String(s))
    DerivativePath(s) = DerivativePath(nothing, s)
end

"""
    SubjectPath

Represents a path to data `parent` associated with subject `label`.
"""
struct SubjectPath{D} <: NeuroPath{D,String}
    dirname::D
    subject::String

    SubjectPath(d::Union{Nothing,StudyPath,DerivativePath}, s::String) = new{typeof(d)}(d, s)
    function SubjectPath(d::DataPath, s::String)
        if basename(d) === "derivatives"
            error("subject paths within the derivative directory must be within a pipeline.")
        end
        return new{DataPath}(d, s)
    end
    SubjectPath(d, s) = SubjectPath(d, String(s))
    SubjectPath(s) = SubjectPath(nothing, s)
end

"""
    SessionPath

Represents a path to data `parent` associated with session `label`.
"""
struct SessionPath{D} <: NeuroPath{D,String}
    dirname::D
    session::String

    SessionPath(d::Union{Nothing,SubjectPath}, s::String) = new{typeof(d)}(d, s)
    SessionPath(d::Union{Nothing,SubjectPath}, s) = SessionPath(d, String(s))
    SessionPath(s) = SessionPath(nothing, s)
end

"""
    ModalityPath

Represents a path to a modality directory `m` for a subject path `p`.
"""
struct ModalityPath{D} <: NeuroPath{D,String}
    dirname::D
    modality::String

    function ModalityPath(d::Union{Nothing,SubjectPath,SessionPath}, s::String)
        _check_modality_name(s)
        new{typeof(d)}(d, s)
    end
    ModalityPath(d::Union{Nothing,SubjectPath,SessionPath}, s) = ModalityPath(d, String(s))
    ModalityPath(s) = ModalityPath(nothing, s)
end

struct FilePath{D} <: NeuroPath{D,String}
    dirname::D
    basename::String

    FilePath(d::Union{Nothing,NeuroPath}, s::String) = new{typeof(d)}(d, s)
    FilePath(d::Union{Nothing,NeuroPath}, s) = FilePath(d, String(s))
    FilePath(s) = FilePath(nothing, s)
end

const File = FilePath{Nothing}

const Modality = ModalityPath{Nothing}

const Subject = SubjectPath{Nothing}

const Session = SessionPath{Nothing}

const Data = DataPath{Nothing}

const Study = StudyPath{Nothing}

const Derivative = DerivativePath{Nothing}

############
# Metadata #
############
struct Spreadsheet
    dict::OrderedDict{Symbol,Any}
    data::DataFrame

    Spreadsheet() = new(OrderedDict{Symbol,Any}(), DataFrame())
end

"""
    DataDescription

The file dataset_description.json is a JSON file describing the dataset. Every dataset MUST
include this file with the following fields:
"""
 struct DataDescription
    "`Name`: Name of the dataset (required)"
    Name::String
    "`BIDSVersion::String`: The version of the BIDS standard that was used. (required)"
    BIDSVersion::String
    "`HEDVersion::String`: If HED tags are used, the version of the HED schema used to validate HED tags for study. (recommended)"
    HEDVersion::String
    "`License::String`: What license is this dataset distributed under? (recommended)"
    License::String
    "`Authors::Vector{String}`: List of individuals who contributed to the creation/curation of the dataset (optional)."
    Authors::Vector{String}
    "`Acknowledgements`: Text acknowledging contributions of individuals or institutions beyond those listed in Authors or Funding (optional)."
    Acknowledgements::String
    "`HowToAcknowledge::String`: Instructions how researchers using this dataset should acknowledge the original authors. This field can also be used to define a publication that should be cited in publications that use the dataset (optional)."
    HowToAcknowledge::String
    "`Funding::String`: List of sources of funding (grant numbers) (optional)"
    Funding::String
    "`EthicsApprovals`::Vector{String}`: List of ethics committee approvals of the research protocols and/or protocol identifiers."
    EthicsApprovals::Vector{String}
    "`ReferencesAndLinks::Vector{String}`: List of references to publication that contain information on the dataset, or links (optional)."
    ReferencesAndLinks::Vector{String}
    "`DatasetDOI::String`:The Document Object Identifier of the dataset (not the corresponding paper) (optional)."
    DatasetDOI::String
end

struct StudyMetadata
    README::String
    LICENSE::String
    CHANGES::String
    dataset_description::DataDescription
    participants::Spreadsheet
end

struct SoftwareContainer
    Type::String
    Tag::String
    URI::String

    SoftwareContainer(Type="", Tag="", URI="") = new(Type, Tag, URI)
end

struct SourceInfo
    URL::String
    DOI::String
    Version::String

    SourceInfo(URL="", DOI="", Version="") = new(URL, DOI, Version)
end

struct GeneratedByInfo
    "Name of the pipeline or process that generated the outputs. Use \"Manual\" to indicate the derivatives were generated by hand, or adjusted manually after an initial run of an automated pipeline"
    Name::String
    "Version of the pipeline."
    Version::String
    "Plain-text description of the pipeline or process that generated the outputs. RECOMMENDED if Name is \"Manual\"."
    Description::String
    "URL where the code used to generate the derivatives may be found."
    CodeURL::String
    "Used to specify the location and relevant attributes of software container image used to produce the derivative. Valid keys in this object include Type, Tag and URI with string values."
    Container::SoftwareContainer

    function GeneratedByInfo(Name; Version="",Description="",CodeURL="",Container=SoftwareContainer())
        new(Name, Version, Description, CodeURL, Container)
    end
end

struct DerivedMetadata
    GeneratedBy::Vector{GeneratedByInfo}
    SourceDatasets::Vector{SourceInfo}

    DerivedDescription() = new(Vector{GeneratedByInfo}(), Vector{SourceInfo}())
end

Base.show(io::IO, ::MIME"text/plain", p::NeuroPath) = print(io, _print_path(p))
_print_path(p::Study) = "study($(basename(p)))"
_print_path(p::Session) = "session($(basename(p)))"
_print_path(p::Subject) = "subject($(basename(p)))"
_print_path(p::Modality) = "modality($(basename(p)))"
_print_path(p::Derivative) = "derivative($(basename(p)))"
_print_path(p::File) = "file($(basename(p)))"
_print_path(p::NeuroPath) = String(p)


text:
## Derivative related code

"""
    tangent(f::Function, c)

Returns a function describing the tangent line to the graph of `f` at `x=c`.

Example. Where does the tangent line intersect the y axis?
```
f(x) = sin(x)
tl(x) = tangent(f, pi/4)(x)  # or tl = tangent(f, pi/3) to use a non-generic function
tl(0)
```

Uses the automatic derivative of `f` to find the slope of the tangent line at `x=c`.

"""
tangent(f,c) = x -> f(c) + f'(c) * (x-c)

"""
    secant(f::Function, a, b)

Returns a function describing the secant line to the graph of `f` at `x=a` and `x=b`.

Example. Where does the secant line intersect the `y` axis?

```
f(x) = sin(x)
a, b = pi/4, pi/3
sl(x) = secant(f, a, b)(x)  # or sl = sl(f, a, b) to use a non-generic function
sl(0)
```

"""
secant(f, a, b) = x -> f(a) + (f(b) - f(a)) / (b-a) * (x - a)

"""
    D(f)

Function interface to `ForwardDiff.derivative`.

Also *overrides* `f'` to take take a derivative.
"""
function D(f, n::Int=1)
    n < 0 && throw(ArgumentError("n is a non-negative integer"))
    n == 0 && return f
    n == 1 && return t -> ForwardDiff.derivative(f, float(t))
    D(D(f), n-1)
end

## Create r' to mean the derivative for functions
## warning, this would be odd for [sin, cos]' as
## is it [sin', cos'] or the transpose...
Base.adjoint(r::Function) = D(r)

ext: jl


text:
# mnist_utils
#
# A set of utility functions for grabbing the MNIST dataset.

using MNIST

# This function does two things:
# 1) Filters the X and y data to only be those in labels
# 2) Does PCA to reduce the dimensionality of X to numdim.
#
# @returns - (X_train_reduced, y_train, X_test_reduced, y_test)
function mnist_reduce_dimension_pca(L::Int, labels::Array{Int})
    # setup training set
    Xtrain, ytrain = traindata()
    # subset out 7 and 9 only
    idx = find(x -> x in labels, ytrain)
    Xtrain = Xtrain[:,idx]'
    ytrain = ytrain[idx]
    # do SVD on Xtrain
    Xmean = mean(Xtrain, 1)
    Xtrain = Xtrain .- Xmean
    Utrain, Strain, Vtrain = svd(Xtrain)

    # prepare test data
    Xtest, ytest = testdata()
    idx = find(x -> x in labels, ytest)
    Xtest = Xtest[:,idx]'
    # make sure we use training data to do renormalization
    Xtest = Xtest .- Xmean
    Xtest = Xtest * Vtrain[:,1:L] * diagm(1./Strain[1:L])

    (Utrain[:,1:L], ytrain, Xtest, ytest)
end

# this function returns the eigenfaces from the PCA decomp
# on the training data
function mnist_eigenfaces(L::Int, labels::Array{Int})
    # setup training set
    Xtrain, ytrain = traindata()
    # subset out 7 and 9 only
    idx = find(x -> x in labels, ytrain)
    Xtrain = Xtrain[:,idx]'
    ytrain = ytrain[idx]
    # do SVD on Xtrain
    Xmean = mean(Xtrain, 1)
    Xtrain = Xtrain .- Xmean
    Utrain, Strain, Vtrain = svd(Xtrain)

    Vtrain[:,1:L]
end


text:
function Geoid(arg0::NormalizedSphericalHarmonicsProvider, arg1::ReferenceEllipsoid)
    return Geoid((NormalizedSphericalHarmonicsProvider, ReferenceEllipsoid), arg0, arg1)
end

function get_body_frame(obj::Geoid)
    return jcall(obj, "getBodyFrame", Frame, ())
end

function get_ellipsoid(obj::Geoid)
    return jcall(obj, "getEllipsoid", ReferenceEllipsoid, ())
end

function get_intersection_point(obj::Geoid, arg0::FieldLine, arg1::FieldVector3D, arg2::Frame, arg3::FieldAbsoluteDate)
    return jcall(obj, "getIntersectionPoint", FieldGeodeticPoint, (FieldLine, FieldVector3D, Frame, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function get_intersection_point(obj::Geoid, arg0::Line, arg1::Vector3D, arg2::Frame, arg3::AbsoluteDate)
    return jcall(obj, "getIntersectionPoint", GeodeticPoint, (Line, Vector3D, Frame, AbsoluteDate), arg0, arg1, arg2, arg3)
end

function get_undulation(obj::Geoid, arg0::jdouble, arg1::jdouble, arg2::AbsoluteDate)
    return jcall(obj, "getUndulation", jdouble, (jdouble, jdouble, AbsoluteDate), arg0, arg1, arg2)
end

function project_to_ground(obj::Geoid, arg0::TimeStampedPVCoordinates, arg1::Frame)
    return jcall(obj, "projectToGround", TimeStampedPVCoordinates, (TimeStampedPVCoordinates, Frame), arg0, arg1)
end

function project_to_ground(obj::Geoid, arg0::Vector3D, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "projectToGround", Vector3D, (Vector3D, AbsoluteDate, Frame), arg0, arg1, arg2)
end

function transform(obj::Geoid, arg0::FieldGeodeticPoint)
    return jcall(obj, "transform", FieldVector3D, (FieldGeodeticPoint,), arg0)
end

function transform(obj::Geoid, arg0::FieldVector3D, arg1::Frame, arg2::FieldAbsoluteDate)
    return jcall(obj, "transform", FieldGeodeticPoint, (FieldVector3D, Frame, FieldAbsoluteDate), arg0, arg1, arg2)
end

function transform(obj::Geoid, arg0::GeodeticPoint)
    return jcall(obj, "transform", Vector3D, (GeodeticPoint,), arg0)
end

function transform(obj::Geoid, arg0::Vector3D, arg1::Frame, arg2::AbsoluteDate)
    return jcall(obj, "transform", GeodeticPoint, (Vector3D, Frame, AbsoluteDate), arg0, arg1, arg2)
end


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

mutable struct IoK8sApiCoreV1ScopedResourceSelectorRequirement <: SwaggerModel
    operator::Any # spec type: Union{ Nothing, String } # spec name: operator
    scopeName::Any # spec type: Union{ Nothing, String } # spec name: scopeName
    values::Any # spec type: Union{ Nothing, Vector{String} } # spec name: values

    function IoK8sApiCoreV1ScopedResourceSelectorRequirement(;operator=nothing, scopeName=nothing, values=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ScopedResourceSelectorRequirement, Symbol("operator"), operator)
        setfield!(o, Symbol("operator"), operator)
        validate_property(IoK8sApiCoreV1ScopedResourceSelectorRequirement, Symbol("scopeName"), scopeName)
        setfield!(o, Symbol("scopeName"), scopeName)
        validate_property(IoK8sApiCoreV1ScopedResourceSelectorRequirement, Symbol("values"), values)
        setfield!(o, Symbol("values"), values)
        o
    end
end # type IoK8sApiCoreV1ScopedResourceSelectorRequirement

const _property_map_IoK8sApiCoreV1ScopedResourceSelectorRequirement = Dict{Symbol,Symbol}(Symbol("operator")=>Symbol("operator"), Symbol("scopeName")=>Symbol("scopeName"), Symbol("values")=>Symbol("values"))
const _property_types_IoK8sApiCoreV1ScopedResourceSelectorRequirement = Dict{Symbol,String}(Symbol("operator")=>"String", Symbol("scopeName")=>"String", Symbol("values")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiCoreV1ScopedResourceSelectorRequirement }) = collect(keys(_property_map_IoK8sApiCoreV1ScopedResourceSelectorRequirement))
Swagger.property_type(::Type{ IoK8sApiCoreV1ScopedResourceSelectorRequirement }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ScopedResourceSelectorRequirement[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ScopedResourceSelectorRequirement }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ScopedResourceSelectorRequirement[property_name]

function check_required(o::IoK8sApiCoreV1ScopedResourceSelectorRequirement)
    (getproperty(o, Symbol("operator")) === nothing) && (return false)
    (getproperty(o, Symbol("scopeName")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ScopedResourceSelectorRequirement }, name::Symbol, val)
end


text:
export MPIElasticPropagatorParams, MPIElasticSource, MPIElasticReceiver, MPIElasticPropagator, compute_PML_Params!,
        MPIElasticPropagatorSolver, MPISimulatedObservation!, MPIElasticReceiver, extract_local_patch
@with_kw mutable struct MPIElasticPropagatorParams
    # number of grids along x,y axis and time steps
    NX::Int64 = 0
    NY::Int64 = 0
    M::Int64 = 0
    N::Int64 = 0
    n::Int64 = 300
    NSTEP::Int64 = 2000 * 2

    # size of grid cell and time step
    DELTAX::Float64 = 10.
    DELTAY::Float64 = 10.
    DELTAT::Float64 = 2.e-3 / 2

    # PML boundary conditon
    f0::Float64 = 5. # source frequency
    vp_ref::Float64 = 2000.
    USE_PML_XMIN::Bool = true
    USE_PML_XMAX::Bool = true
    USE_PML_YMIN::Bool = true
    USE_PML_YMAX::Bool = true
    NPOINTS_PML::Int64 = 12
    NPOWER::Float64 = 2.
    K_MAX_PML::Float64 = 1.
    ALPHA_MAX_PML::Float64 = 2. * π * ( f0/2.)
    Rcoef::Float64 = 0.001 # Relative reflection coefficient

    # Auxilliary Data
    ax::Union{Missing, Array{Float64, 2}} = missing
    bx::Union{Missing, Array{Float64, 2}} = missing
    kx::Union{Missing, Array{Float64, 2}} = missing
    ay::Union{Missing, Array{Float64, 2}} = missing
    by::Union{Missing, Array{Float64, 2}} = missing
    ky::Union{Missing, Array{Float64, 2}} = missing

    # display params
    IT_DISPLAY::Int64 = 0

    # mpi
    II::Int64 = -1
    JJ::Int64 = -1

    # indices
    i_1j::Union{Missing, Array{Int64, 1}} = missing
    i_2j::Union{Missing, Array{Int64, 1}} = missing
    i1j::Union{Missing, Array{Int64, 1}} = missing
    i2j::Union{Missing, Array{Int64, 1}} = missing
    ij_1::Union{Missing, Array{Int64, 1}} = missing
    ij_2::Union{Missing, Array{Int64, 1}} = missing
    ij1::Union{Missing, Array{Int64, 1}} = missing
    ij2::Union{Missing, Array{Int64, 1}} = missing
    ij::Union{Missing, Array{Int64, 1}} = missing
    i1j1::Union{Missing, Array{Int64, 1}} = missing

end

"""
    MPIElasticSource

The input is the global index.

`srci`, `srcj`: nsrc
`srcv`: nt × nsrc, each row corresponds to a source time function
"""
mutable struct MPIElasticSource
    srci::Union{Missing, PyObject}
    srcj::Union{Missing, PyObject}
    srctype::Union{Missing,PyObject}
    srcv::Union{Missing, PyObject}
    function MPIElasticSource(param::MPIElasticPropagatorParams, srci::Array{Int64,1},
         srcj::Array{Int64,1}, srctype::Array{Int64, 1}, srcv::Union{Array{Float64,2}, PyObject})
        II, JJ = param.II, param.JJ
        nsrc = length(srci)
        srcv = constant(srcv)
        srcv = srcv'
        local_srci = Int64[]
        local_srcj = Int64[]
        local_srcv = PyObject[]
        local_srctype = Int64[]
        for i = 1:nsrc
            if ((II - 1) * param.n + 1 <= srci[i] <= II * param.n) &&
                    ((JJ - 1) * param.n + 1 <= srcj[i] <= JJ * param.n)
                    push!(local_srci, srci[i] - (II - 1) * param.n + 2)
                    push!(local_srcj, srcj[i] - (JJ - 1) * param.n + 2)
                    push!(local_srcv, srcv[i])
                    push!(local_srctype, srctype[i])
            end
        end
        local_srcv = hcat(local_srcv...)
        if length(local_srci)>0
            local_srci, local_srcj, local_srctype, local_srcv = convert_to_tensor([local_srci, local_srcj, local_srctype, local_srcv], [Int64, Int64, Int64, Float64])
            new(local_srci, local_srcj, local_srctype, local_srcv)
        else
            new(missing, missing, missing, missing)
        end
    end
end

mutable struct MPIElasticReceiver
    rcvi::Union{Missing,Array{Int64, 1}}
    rcvj::Union{Missing,Array{Int64, 1}}
    rcvtype::Union{Missing, Array{Int64, 1}}
    rcvv::Union{Missing, PyObject}
    function MPIElasticReceiver(param::MPIElasticPropagatorParams, rcvi::Array{Int64,1},
         rcvj::Array{Int64,1}, rcvtype::Array{Int64, 1})
        II, JJ = param.II, param.JJ
        nrcv = length(rcvi)
        local_rcvi = Int64[]
        local_rcvj = Int64[]
        local_rcvtype = Int64[]
        for i = 1:nrcv
            if ((II - 1) * param.n + 1 <= rcvi[i] <= II * param.n) &&
                    ((JJ - 1) * param.n + 1 <= rcvj[i] <= JJ * param.n)
                    push!(local_rcvi, rcvi[i] - (II - 1) * param.n)
                    push!(local_rcvj, rcvj[i] - (JJ - 1) * param.n)
                    push!(local_rcvtype, rcvtype[i])
            end
        end
        if length(local_rcvi)>0
            new(local_rcvi, local_rcvj, local_rcvtype, missing)
        else
            new(missing, missing, missing, missing)
        end
    end
end

"""
    ElasticPropagator

A structure holding all the information for acoustic simulation
"""
mutable struct MPIElasticPropagator
    param::MPIElasticPropagatorParams
    src::MPIElasticSource
    vx::PyObject
    vy::PyObject
    sigmaxx::PyObject
    sigmayy::PyObject
    sigmaxy::PyObject
    mem::PyObject
end

function get_mpi_id2(a, b, n)
    idx = Int64[]
    for i = 1:length(a)
        for j = 1:length(b)
            push!(idx, (n+4)*(a[i]-1)+b[j])
        end
    end
    idx
end

function extract_local_patch(param::MPIElasticPropagatorParams, v::Union{PyObject, Array{Float64,2}}; deps = missing, tag = 0)
    v = v[(param.II - 1) * param.n + 1:param.II * param.n, (param.JJ - 1) * param.n + 1:param.JJ * param.n]
    v = mpi_halo_exchange2(v, param.M, param.N, deps = deps, tag = tag)
end

function compute_PML_Params!(param::MPIElasticPropagatorParams)
    NX, NY = param.NX, param.NY
    n = param.n
    r = mpi_rank()
    param.M = div(NX,n)
    param.N = div(NY,n)
    param.II, param.JJ = div(r, param.N)+1, mod(r, param.N)+1
    if mpi_size() != param.M * param.N
        error("MPI Size Mismatch: mpi_size = $(mpi_size()), M = $(param.M), N = $(param.N)")
    end
    @assert NX>0 && NY>0
    @assert mod(NX, n)==0 && mod(NY, n)==0

    # ID
    k1 = 4:n+3
    k2 = 5:n+4
    k_1 = 2:n+1
    k_2 = 1:n
    kk = 3:n+2
    param.i_1j = get_mpi_id2(k_1, kk, n)
    param.i_2j = get_mpi_id2(k_2, kk, n)
    param.i1j = get_mpi_id2(k1, kk, n)
    param.i2j = get_mpi_id2(k2, kk, n)
    param.ij_1 = get_mpi_id2(kk, k_1, n)
    param.ij_2 = get_mpi_id2(kk, k_2, n)
    param.ij1 = get_mpi_id2(kk, k1, n)
    param.ij2 = get_mpi_id2(kk, k2, n)
    param.ij = get_mpi_id2(kk, kk, n)
    param.i1j1 = get_mpi_id2(k1, k1, n)

    ## define profile of absorption in PML region
    NPOINTS_PML = param.NPOINTS_PML
    NPOWER = param.NPOWER
    NX = param.NX
    NY = param.NY
    DELTAX = param.DELTAX
    DELTAY = param.DELTAY
    DELTAT = param.DELTAT
    K_MAX_PML = param.K_MAX_PML
    ALPHA_MAX_PML = param.ALPHA_MAX_PML
    USE_PML_XMIN = param.USE_PML_XMIN
    USE_PML_XMAX = param.USE_PML_XMAX
    USE_PML_YMIN = param.USE_PML_YMIN
    USE_PML_YMAX = param.USE_PML_YMAX
    Rcoef = param.Rcoef

    thickness_PML_x = NPOINTS_PML * DELTAX
    thickness_PML_y = NPOINTS_PML * DELTAY

    # reflection coefficient
    d0_x = - (NPOWER + 1) * param.vp_ref * log(Rcoef) / (2. * thickness_PML_x)
    d0_y = - (NPOWER + 1) * param.vp_ref * log(Rcoef) / (2. * thickness_PML_y)

    xoriginleft = thickness_PML_x
    xoriginright = (NX-0.5)*DELTAX - thickness_PML_x

    a_x = zeros(n)
    a_y = zeros(n)
    a_x_half = zeros(n)
    a_y_half = zeros(n)
    b_x = zeros(n)
    b_y = zeros(n)
    b_x_half = zeros(n)
    b_y_half = zeros(n)
    d_x = zeros(n)
    d_y = zeros(n)
    d_x_half = zeros(n)
    d_y_half = zeros(n)
    K_x = ones(n)
    K_y = ones(n)
    K_x_half = ones(n)
    K_y_half = ones(n)
    alpha_x = zeros(n)
    alpha_y = zeros(n)
    alpha_x_half = zeros(n)
    alpha_y_half = zeros(n)

    for i = 1:n
        xval = DELTAX * (i-1) + (param.II - 1) * param.n * DELTAX
        if (USE_PML_XMIN)

            abscissa_in_PML = xoriginleft - xval
            if (abscissa_in_PML >= 0.0)
                abscissa_normalized = abscissa_in_PML / thickness_PML_x
                d_x[i] = d0_x * abscissa_normalized^NPOWER
                K_x[i] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_x[i] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_x[i] = ALPHA_MAX_PML * abscissa_normalized
            end

            abscissa_in_PML = xoriginleft - (xval + DELTAX/2.)
            if (abscissa_in_PML >= 0.0)
                abscissa_normalized = abscissa_in_PML / thickness_PML_x
                d_x_half[i] = d0_x * abscissa_normalized^NPOWER
                K_x_half[i] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_x_half[i] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_x_half[i] = ALPHA_MAX_PML * abscissa_normalized
            end

        end

        if (USE_PML_XMAX)
            abscissa_in_PML = xval - xoriginright
            if (abscissa_in_PML >= 0)
                abscissa_normalized = abscissa_in_PML / thickness_PML_x
                d_x[i] = d0_x * abscissa_normalized^NPOWER
                K_x[i] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_x[i] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_x[i] = ALPHA_MAX_PML * abscissa_normalized
            end

            abscissa_in_PML = xval + DELTAX/2.0 - xoriginright
            if (abscissa_in_PML >= 0)
                abscissa_normalized = abscissa_in_PML / thickness_PML_x
                d_x_half[i] = d0_x * abscissa_normalized^NPOWER
                K_x_half[i] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_x_half[i] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_x_half[i] = ALPHA_MAX_PML * abscissa_normalized
            end
        end

        if (alpha_x[i] < 0)
            alpha_x[i] = 0
        end
        if (alpha_x_half[i] < 0)
            alpha_x_half[i] = 0
        end

        b_x[i] = exp(- (d_x[i] / K_x[i] + alpha_x[i]) * DELTAT)
        b_x_half[i] = exp(- (d_x_half[i] / K_x_half[i] + alpha_x_half[i]) * DELTAT)

        if (abs(d_x[i]) > 1e-6)
            a_x[i] = d_x[i] * (b_x[i] - 1.) / (K_x[i] * (d_x[i] + K_x[i] * alpha_x[i]))
        end
        if (abs(d_x_half[i]) > 1e-6)
            a_x_half[i] = d_x_half[i] * (b_x_half[i] - 1.) / (K_x_half[i] * (d_x_half[i] + K_x_half[i] * alpha_x_half[i]))
        end
    end

    yoriginbottom = thickness_PML_y
    yorigintop = (NY-0.5)*DELTAY - thickness_PML_y

    for j = 1:n

        yval = DELTAY * (j-1) + (param.JJ - 1) * param.n * DELTAX
        if param.USE_PML_YMIN
            abscissa_in_PML = yoriginbottom - yval
            if (abscissa_in_PML >= 0.0)
                abscissa_normalized = abscissa_in_PML / thickness_PML_y
                d_y[j] = d0_y * abscissa_normalized^NPOWER
                K_y[j] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_y[j] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_y[j] = ALPHA_MAX_PML * abscissa_normalized
            end
            abscissa_in_PML = yoriginbottom - (yval + DELTAY/2.)
            if abscissa_in_PML >= 0
                abscissa_normalized = abscissa_in_PML / thickness_PML_y
                d_y_half[j] = d0_y * abscissa_normalized^NPOWER
                K_y_half[j] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_y_half[j] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_y_half[j] = ALPHA_MAX_PML * abscissa_normalized
            end
        end

        if param.USE_PML_YMAX
            abscissa_in_PML = yval - yorigintop
            if abscissa_in_PML >= 0
                abscissa_normalized = abscissa_in_PML / thickness_PML_y
                d_y[j] = d0_y * abscissa_normalized^NPOWER
                K_y[j] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_y[j] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_y[j] = ALPHA_MAX_PML * abscissa_normalized
            end

            abscissa_in_PML = yval + DELTAY/2. - yorigintop
            if abscissa_in_PML >= 0
                abscissa_normalized = abscissa_in_PML / thickness_PML_y
                d_y_half[j] = d0_y * abscissa_normalized^NPOWER
                K_y_half[j] = 1. + (K_MAX_PML - 1.) * abscissa_normalized^NPOWER
                # alpha_y_half[j] = ALPHA_MAX_PML * (1. - abscissa_normalized)
                alpha_y_half[j] = ALPHA_MAX_PML * abscissa_normalized
            end
        end

        b_y[j] = exp(- (d_y[j] / K_y[j] + alpha_y[j]) * DELTAT)
        b_y_half[j] = exp(- (d_y_half[j] / K_y_half[j] + alpha_y_half[j]) * DELTAT)
        if abs(d_y[j]) > 1e-6
            a_y[j] = d_y[j] * (b_y[j] - 1.) / (K_y[j] * (d_y[j] + K_y[j] * alpha_y[j]))
        end
        if (abs(d_y_half[j]) > 1e-6)
            a_y_half[j] = d_y_half[j] * (b_y_half[j] - 1.) / (K_y_half[j] * (d_y_half[j] + K_y_half[j] * alpha_y_half[j]))
        end
    end
    # @show norm(d_x-d_y)

    ax = [a_x'; a_x_half']
    bx = [b_x'; b_x_half']
    kx = [K_x'; K_x_half']
    ay = [a_y'; a_y_half']
    by = [b_y'; b_y_half']
    ky = [K_y'; K_y_half']

    param.ax = ax
    param.bx = bx
    param.kx = kx
    param.ay = ay
    param.by = by
    param.ky = ky
    param
end

function MPIElasticPropagatorSolver(param::MPIElasticPropagatorParams, src::MPIElasticSource,
    ρ::Union{PyObject, Array{Float64, 2}}, λ::Union{PyObject, Array{Float64, 2}}, μ::Union{PyObject, Array{Float64, 2}};
    tag_offset::Int64 = 0, dep = nothing)
    ρ = reshape(convert_to_tensor(ρ),(-1,))
    λ = reshape(convert_to_tensor(λ),(-1,))
    μ = reshape(convert_to_tensor(μ),(-1,))
    n = param.n
    @assert size(ρ)==size(λ)==size(μ)==((n+4)^2,)
    @assert param.II>0
    if !ismissing(src.srcv)
        @assert size(src.srcv, 2)==length(src.srci)==length(src.srcj)==length(src.srctype)
        @assert size(src.srcv, 1)==param.NSTEP
    end
    NSTEP = param.NSTEP
    vx_arr = TensorArray(NSTEP+1);
    vx_arr = write(vx_arr, 1, constant(zeros(n^2)))
    vy_arr = TensorArray(NSTEP+1);
    vy_arr = write(vy_arr, 1, constant(zeros(n^2)))
    sigmaxx_arr = TensorArray(NSTEP+1);
    sigmaxx_arr = write(sigmaxx_arr, 1, constant(zeros(n^2)))
    sigmayy_arr = TensorArray(NSTEP+1);
    sigmayy_arr = write(sigmayy_arr, 1, constant(zeros(n^2)))
    sigmaxy_arr = TensorArray(NSTEP+1);
    sigmaxy_arr = write(sigmaxy_arr, 1, constant(zeros(n^2)))
    mem_arr = [TensorArray(NSTEP+1) for i = 1:8]
    for i = 1:8
        mem_arr[i] = write(mem_arr[i], 1, constant(zeros(n^2)))
    end
    ax,bx,kx,ay,by,ky = param.ax,param.bx,param.kx,param.ay,param.by,param.ky
    function condition(i, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arr...)
        i<=NSTEP+1
    end

    if isnothing(dep)
        dep = 0.0
    end

    function body(i, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arr...)
        mem = Array{PyObject}(undef, 8)
        σxx,σyy,σxy,vx,vy =
            read(sigmaxx_arr, i-1), read(sigmayy_arr, i-1), read(sigmaxy_arr, i-1),
            read(vx_arr, i-1), read(vy_arr, i-1)
        for k = 1:8
            m = read(mem_arr[k], i-1)
            mem[k] = _reshape_and_halo_exchange(m, param, 18*i+k-1+tag_offset, k>1 ? mem[k-1][1] + dep : λ[1] + ρ[1] + μ[1] + vy[1] + σxy[1] + dep)
        end

        mpi_idx = 18*i+8+tag_offset

        halo = (op, deps, mpi_idx)->begin
            op = _reshape_and_halo_exchange(op, param, mpi_idx, deps)
            mpi_idx = mpi_idx + 1
            op, mpi_idx
        end
        vx, mpi_idx = halo(vx, mem[8][1], mpi_idx)
        vy, mpi_idx = halo(vy, vx[1], mpi_idx)
        σxx, mpi_idx = halo(σxx, vy[1], mpi_idx)
        σyy, mpi_idx = halo(σyy, σxx[1], mpi_idx)
        σxy, mpi_idx = halo(σxy, σyy[1], mpi_idx)

        srci,srcj,srctype,srcv = ElasticSourceAtTimeT(src, i-1)
        σxx, σyy, σxy,vx, vy, mem, mpi_idx = one_step(param, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
            ax, bx, ay, by, kx, ky, srci, srcj, srcv, srctype, mpi_idx)

        sigmaxx_arr,sigmayy_arr,sigmaxy_arr,vx_arr,vy_arr =
            write(sigmaxx_arr, i, σxx), write(sigmayy_arr, i, σyy), write(sigmaxy_arr, i, σxy),
                    write(vx_arr, i, vx), write(vy_arr, i, vy)

        mem_arr_ = Array{PyObject}(undef, 8)
        for k = 1:8
            mem_arr_[k] = write(mem_arr[k], i, mem[k])
        end

        if param.IT_DISPLAY>0
            op = tf.cond(tf.equal(tf.math.floormod(i,param.IT_DISPLAY),1), ()->tf.print("time step = ", i-1, " out of $NSTEP"), ()->tf.no_op())
            i = bind(i, op)
        end

        i+1, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arr_...
    end
    i = constant(2, dtype=Int32)
    _, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arrs =
        while_loop(condition, body, [i, vx_arr, vy_arr, sigmaxx_arr, sigmayy_arr, sigmaxy_arr, mem_arr...])

    vx = reshape(stack(vx_arr), (param.NSTEP+1, n, n))
    vy = reshape(stack(vy_arr), (param.NSTEP+1, n, n))
    σxx = reshape(stack(sigmaxx_arr), (param.NSTEP+1, n, n))
    σyy = reshape(stack(sigmayy_arr), (param.NSTEP+1, n, n))
    σxy = reshape(stack(sigmaxy_arr), (param.NSTEP+1, n, n))
    MPIElasticPropagator(
        param, src, vx, vy, σxx, σyy, σxy, mem_arrs
    )
end

function _reshape_and_print(o, m, n, info)
    o = reshape(o, (m, n))
    o = print_tensor(o, info)
    reshape(o, (m*n,))
end

function _reshape_and_halo_exchange(o, param, tag, deps)
    @assert size(o)==(param.n^2,)
    o = reshape(o, (param.n, param.n))
    out = mpi_halo_exchange2(o, param.M, param.N, tag = tag, deps= deps)
    reshape(out, ((param.n+4)^2,))
end

function fw1(param::MPIElasticPropagatorParams, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
        ax, bx, ay, by, kx, ky, mpi_idx)
    NX, NY = param.NX, param.NY
    n = param.n
    DELTAX, DELTAY, DELTAT = param.DELTAX, param.DELTAY, param.DELTAT
    i_1j = param.i_1j
    i_2j = param.i_2j
    i1j = param.i1j
    i2j = param.i2j
    ij_1 = param.ij_1
    ij_2 = param.ij_2
    ij1 = param.ij1
    ij2 = param.ij2
    ij = param.ij

    # λ_ = 0.5*(λ[i1j]+λ[ij])
    # μ_ = 0.5*(μ[i1j]+μ[ij])
    λ_ = λ[ij]
    μ_ = μ[ij]
    λμ = λ_ + 2μ_
    value_dvx_dx = (27vx[i1j]-27vx[ij]-vx[i2j]+vx[i_1j])/(24*DELTAX)
    value_dvy_dy = (27*vy[ij]-27*vy[ij_1]-vy[ij1]+vy[ij_2]) / (24*DELTAY)
    # broadcast
    mem[1] =  makevector(mem[1], ij, adbroadcast(bx[2,:], mem[1][ij], 1) + adbroadcast(ax[2,:], value_dvx_dx, 1));
    mem[2] =  makevector(mem[2], ij, adbroadcast(by[1,:], mem[2][ij], 2) + adbroadcast(ay[1,:], value_dvy_dy, 2));
    value_dvx_dx = adbroadcast(value_dvx_dx, kx[2,:], 3) + mem[1][ij]
    value_dvy_dy = adbroadcast(value_dvy_dy, ky[1,:], 4) + mem[2][ij]
    σxx = scatter_add_op(σxx, ij, (λμ * value_dvx_dx + λ_ * value_dvy_dy) * DELTAT )
    σyy = scatter_add_op(σyy, ij, (λμ * value_dvy_dy + λ_ * value_dvx_dx) * DELTAT )

    σxx = _reshape_and_halo_exchange(σxx[ij], param, mpi_idx, σxy[1])
    σyy = _reshape_and_halo_exchange(σyy[ij], param, mpi_idx+1, σxx[1])

    return vx, vy, σxx, σyy,σxy, mem, mpi_idx+2
end

function fw2(param::MPIElasticPropagatorParams, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
    ax, bx, ay, by, kx, ky, mpi_idx)
    NX, NY = param.NX, param.NY
    n = param.n
    DELTAX, DELTAY, DELTAT = param.DELTAX, param.DELTAY, param.DELTAT
    i_1j = param.i_1j
    i_2j = param.i_2j
    i1j = param.i1j
    i2j = param.i2j
    ij_1 = param.ij_1
    ij_2 = param.ij_2
    ij1 = param.ij1
    ij2 = param.ij2
    ij = param.ij

    # μ_ = 0.5*(μ[ij]+μ[ij1])
    μ_ = μ[ij]
    value_dvy_dx = (27vy[ij]-27vy[i_1j]-vy[i1j]+vy[i_2j])/(24*DELTAX)
    value_dvx_dy = (27*vx[ij1]-27*vx[ij]-vx[ij2]+vx[ij_1]) / (24*DELTAY)

    # broadcast
    mem[3] = makevector(mem[3], ij,  adbroadcast(bx[1,:], mem[3][ij], 1) + adbroadcast(ax[1, :], value_dvy_dx, 1))
    mem[4] = makevector(mem[4], ij,  adbroadcast(by[2,:], mem[4][ij], 2) + adbroadcast(ay[2, :], value_dvx_dy, 2));

    value_dvy_dx = adbroadcast(value_dvy_dx, kx[1,:], 3) + mem[3][ij]
    value_dvx_dy = adbroadcast(value_dvx_dy, ky[1,:], 4) + mem[4][ij]

    σxy = scatter_add_op(σxy, ij, μ_*(value_dvy_dx + value_dvx_dy) * DELTAT )

    σxy = _reshape_and_halo_exchange(σxy[ij], param, mpi_idx, σyy[1])

    return vx, vy, σxx, σyy, σxy,mem, mpi_idx+1
end

function fw3(param::MPIElasticPropagatorParams, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
    ax, bx, ay, by, kx, ky, mpi_idx)
    NX, NY = param.NX, param.NY
    n = param.n
    DELTAX, DELTAY, DELTAT = param.DELTAX, param.DELTAY, param.DELTAT
    i_1j = param.i_1j
    i_2j = param.i_2j
    i1j = param.i1j
    i2j = param.i2j
    ij_1 = param.ij_1
    ij_2 = param.ij_2
    ij1 = param.ij1
    ij2 = param.ij2
    ij = param.ij

    value_dsigmaxx_dx = (27*σxx[ij]-27*σxx[i_1j]-σxx[i1j]+σxx[i_2j]) / (24*DELTAX);
    value_dsigmaxy_dy = (27*σxy[ij]-27*σxy[ij_1]-σxy[ij1]+σxy[ij_2]) / (24*DELTAY);

    mem[5] = makevector(mem[5],ij, adbroadcast(bx[1,:], mem[5][ij], 1) + adbroadcast(ax[1,:], value_dsigmaxx_dx, 1))
    mem[6] = makevector(mem[6],ij, adbroadcast(by[1,:], mem[6][ij], 2) + adbroadcast(ay[1,:], value_dsigmaxy_dy, 2))

    value_dsigmaxx_dx = adbroadcast(value_dsigmaxx_dx, kx[1,:], 3) + mem[5][ij]
    value_dsigmaxy_dy = adbroadcast(value_dsigmaxy_dy, ky[1,:], 4) + mem[6][ij]

    vx = scatter_add_op(vx, ij, (value_dsigmaxx_dx + value_dsigmaxy_dy) * DELTAT / ρ[ij])
    vx = _reshape_and_halo_exchange(vx[ij], param, mpi_idx, σxy[1])

    return vx, vy, σxx, σyy, σxy, mem, mpi_idx+1
end

function fw4(param::MPIElasticPropagatorParams, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
    ax, bx, ay, by, kx, ky, mpi_idx)
    NX, NY, n = param.NX, param.NY, param.n
    DELTAX, DELTAY, DELTAT = param.DELTAX, param.DELTAY, param.DELTAT
    i_1j = param.i_1j
    i_2j = param.i_2j
    i1j = param.i1j
    i2j = param.i2j
    ij_1 = param.ij_1
    ij_2 = param.ij_2
    ij1 = param.ij1
    ij2 = param.ij2
    ij = param.ij
    i1j1 = param.i1j1

    # ρ_ = 0.25 * (ρ[ij] + ρ[i1j] + ρ[i1j1] + ρ[ij1]);
    ρ_ = ρ[ij]

    value_dsigmaxy_dx = (27*σxy[i1j]-27*σxy[ij]-σxy[i2j]+σxy[i_1j]) / (24*DELTAX);
    value_dsigmayy_dy = (27*σyy[ij1]-27*σyy[ij]-σyy[ij2]+σyy[ij_1]) / (24*DELTAY);

    mem[7] = makevector(mem[7], ij, adbroadcast(bx[2,:], mem[7][ij], 1) + adbroadcast(ax[2,:] , value_dsigmaxy_dx, 1))
    mem[8] = makevector(mem[8], ij, adbroadcast(by[2,:], mem[8][ij], 2) + adbroadcast(ay[2,:] , value_dsigmayy_dy, 2))

    value_dsigmaxy_dx = adbroadcast(value_dsigmaxy_dx, kx[2,:], 3) + mem[7][ij]
    value_dsigmayy_dy = adbroadcast(value_dsigmayy_dy, ky[2,:], 4) + mem[8][ij]

    vy = scatter_add_op( vy, ij, (value_dsigmaxy_dx + value_dsigmayy_dy) * DELTAT / ρ_)
    vy = _reshape_and_halo_exchange(vy[ij], param, mpi_idx, vx[1])
    return vx, vy, σxx, σyy, σxy,mem, mpi_idx+1
end

function one_step(param::MPIElasticPropagatorParams, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
    ax, bx, ay, by, kx, ky, srci, srcj, srcv, srctype, mpi_idx)
    n = param.n
    @assert size(ρ)==size(λ)==size(μ)==size(vx)==size(vy)==size(σxx)==size(σyy)==size(σxy)==((n+4)^2,)

    add_source = load_op_and_grad("$(@__DIR__)/../deps/CustomOps/build/libADSeismic","add_source", multiple=true)
    vx, vy, σxx, σyy, σxy,mem, mpi_idx = fw1(param, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
                    ax, bx, ay, by, kx, ky, mpi_idx)
    vx, vy, σxx, σyy, σxy,mem, mpi_idx = fw2(param, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
                    ax, bx, ay, by, kx, ky, mpi_idx)
    vx, vy, σxx, σyy, σxy,mem, mpi_idx = fw3(param, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
                    ax, bx, ay, by, kx, ky, mpi_idx)
    vx, vy, σxx, σyy, σxy,mem, mpi_idx = fw4(param, ρ, λ, μ, vx, vy, σxx, σyy, σxy,mem,
                    ax, bx, ay, by, kx, ky, mpi_idx)
    if !ismissing(srci)
        σxx, σyy, σxy,vx, vy  =  add_source(σxx, σyy, σxy,vx, vy, constant(srci),constant(srcj),constant(srctype),constant(n+2),constant(n+2),srcv)
    end
    ij = param.ij
    σxx, σyy, σxy,vx, vy  = σxx[ij], σyy[ij], σxy[ij],vx[ij],vy[ij]
    for i = 1:8
        mem[i] = mem[i][ij]
    end
    return σxx, σyy, σxy,vx, vy, mem, mpi_idx
end

function ElasticSourceAtTimeT(src::MPIElasticSource,i::PyObject)
    if ismissing(src.srci)
        return missing, missing, missing, missing
    else
        return src.srci, src.srcj, src.srctype, src.srcv[i]
    end
end

"""
    SimulatedObservation!(ap::ElasticPropagator, rcv::ElasticReceiver)

Extract and save simulated displacement u into rcv::ElasticReceiver.
"""
function MPISimulatedObservation!(ep::MPIElasticPropagator, rcv::MPIElasticReceiver)
    if ismissing(rcv.rcvi)
        op = group([ep.vx, ep.vy, ep.sigmaxx, ep.sigmayy, ep.sigmaxy])
        rcv.rcvv = op
        return
    end

    N = (ep.param.n)*(ep.param.n)
    get_receive = load_op_and_grad("$(@__DIR__)/../deps/CustomOps/build/libADSeismic", "get_receive")
    u = get_receive(tf.reshape(ep.vx, (ep.param.NSTEP+1, N)),
                    tf.reshape(ep.vy, (ep.param.NSTEP+1, N)),
                    tf.reshape(ep.sigmaxx, (ep.param.NSTEP+1, N)),
                    tf.reshape(ep.sigmayy, (ep.param.NSTEP+1, N)),
                    tf.reshape(ep.sigmaxy, (ep.param.NSTEP+1, N)),
                    convert_to_tensor(rcv.rcvi), convert_to_tensor(rcv.rcvj), convert_to_tensor(rcv.rcvtype),
                        constant(ep.param.n-2), constant(ep.param.n-2))
    rcv.rcvv = u
end


text:
export ROC, auc

struct ROC{T <: Real} <: AbstractBinaryMetric
  labels::BitVector
  scores::Vector{T}
  AUC::T
  FPR::Vector{T}
  TPR::Vector{T}
  OP::Vector{T}
  function ROC(labels,scores)
    B = BinaryMetrics(labels, scores)
    OP = op(B)
    FPR = fpr(B, OP)
    TPR = tpr(B, OP)
    AUC = get_auc([0.0;FPR], [0.0;TPR])
    new{eltype(FPR)}(B.labels, B.scores, AUC, FPR, TPR, OP)
  end
end

function Base.show(io::IO, B::ROC)
  println(io,"Receiver Operating Characteristic (ROC) curve")
  AUC = round(auc(B),digits=2)
  print(io,"Samples:$(length(B)), P:$(positives(B)), N:$(negatives(B)), AUC:$(AUC)")
end

"""
  ROC(labels, scores)

Generates data for plotting a Receiver Operating Characteristic (ROC) curve.

* `labels` must be a boolean array e.g. either containing `true`,`false` or `0`,`1`.
* `scores` must be an array with the same length of `labels` containing the scores of the binary classifier.

You can generate a plot automatically using the following commands:
```julia
roc = ROC(labels, scores)

using Plots
plot(roc)
```

The Area Under Curve (AUC) FPR and TPR can be accessed by `auc(roc)`, `fpr(roc)` and `tpr(roc)` respectively.

"""
ROC(B::BinaryMetrics) = ROC(B.labels, B.scores)

auc(curve::ROC) = curve.AUC
fpr(curve::ROC) = curve.FPR
tpr(curve::ROC) = curve.TPR
op(curve::ROC) = curve.OP

function get_curve(curve::ROC)
  return [0.0;curve.FPR], [0.0;curve.TPR]
end


text:
# # tests for logistic regression: sum_i log(1+exp(-b_i<a_i,x>)) + λ \|x\|_1

@testset "logistic-l1" begin
    using LinearAlgebra
    using CIAOAlgorithms
    using ProximalOperators
    using ProximalAlgorithms: IterationTools
    using Base.Iterators: take

    T = Float64
    # create the two classes
    x_class1 = [
        5.1 3.5 1.4 0.2 1.0
        4.9 3.0 1.4 0.2 1.0
        4.7 3.2 1.3 0.2 1.0
        4.6 3.1 1.5 0.2 1.0
    ]
    x_class2 = [
        5.7 3.0 4.2 1.2 1.0
        5.7 2.9 4.2 1.3 1.0
        6.2 2.9 4.3 1.3 1.0
        5.1 2.5 3.0 1.1 1.0
    ]

    xs = [x_class1; x_class2]   # taking 500 sample from each
    ys = [fill(1.0, size(x_class1, 1)); fill(-1.0, size(x_class2, 1))]

    # preparations for the solver
    x_star = [0.0, 0.924160995722576, -1.1343956493097298, 0.0, 0.0]

    N, n = size(ys, 1), size(xs, 2)

    F = Vector{Precompose}(undef, 0)
    L = Vector{T}(undef, 0)
    for i = 1:N
        f = Precompose(LogisticLoss([ys[i]], 1.0), reshape(xs[i, :], 1, n), 1.0)
        push!(F, f)
        # compute L
        Lf = 0.25 * norm(xs[i, :])^2
        push!(L, Lf)
    end

    # compute L for the nonadaptive case
    lam = 1 / N
    g = NormL1(lam)
    x0 = ones(n)

    maxit = 9000
    tol = 1e-4

    @testset "Finito" begin
        # sweeping 1, 2, 3 for randomined, cyclic and shuffled sampling strategies, respectively.

        # basic finito
        @testset "nominal Finito" for sweeping in collect(1:3)
            solver = CIAOAlgorithms.Finito{T}(maxit = maxit, sweeping = sweeping)
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test norm(x_finito - x_star, Inf) < tol
        end

        # limited memory finito
        @testset "LFinito" for sweeping in collect(2:3)
            # @testset "cyclical" begin
            solver =
                CIAOAlgorithms.Finito{T}(maxit = maxit, sweeping = sweeping, LFinito = true)
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test norm(x_finito - x_star, Inf) < tol
        end

        # basic finito with minibatch
        vec_ref = [(1, 2), (2, 2), (3, 3)] # different samplings and batch sizes
        @testset "LFinito_minibatch" for (sweeping, batch) in vec_ref
            solver = CIAOAlgorithms.Finito{T}(
                maxit = maxit,
                sweeping = sweeping,
                minibatch = (true, batch),
            )
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test norm(x_finito - x_star, Inf) < tol
        end

        # limited memory finito with minibatch
        vec_ref = [(2, 1), (2, 2), (3, 3)] # different samplings and batch sizes
        @testset "LFinito_minibatch" for (sweeping, batch) in vec_ref
            solver = CIAOAlgorithms.Finito{T}(
                maxit = maxit,
                sweeping = sweeping,
                LFinito = true,
                minibatch = (true, batch),
            )
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test norm(x_finito - x_star, Inf) < tol
        end

        # test with user defined stepsizes
        @testset "γ and L as scalars" begin
            @testset "randomized" begin
                γ = N / maximum(L)
                solver = CIAOAlgorithms.Finito{T}(maxit = maxit, γ = γ)
                x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
                @test norm(x_finito - x_star, Inf) < tol
            end
            @testset "cyclic" begin
                solver = CIAOAlgorithms.Finito{T}(maxit = maxit)
                x_finito, it_finito = solver(x0, F = F, g = g, L = maximum(L), N = N)
                @test norm(x_finito - x_star, Inf) < tol
            end
        end

        # test the iterator
        @testset "the iterator" for LFinito in [true, false]
            solver = CIAOAlgorithms.Finito{T}(sweeping = 2, LFinito = LFinito, maxit = 10)
            iter = CIAOAlgorithms.iterator(solver, x0, F = F, g = g, L = L, N = N)
            @test iter.x0 === x0

            for state in take(iter, 2)
                @test solution(state) === state.z
                @test eltype(solution(state)) == T
            end
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test solution(IterationTools.loop(take(iter, 10))) == x_finito
        end
    end

    @testset "SVRG" begin
        γ = 1 / (10 * maximum(L))
        @testset "SVRG-Base" begin
            solver = CIAOAlgorithms.SVRG{T}(maxit = maxit, γ = γ)
            x_SVRG, it_SVRG = solver(x0, F = F, g = g, N = N)
            @test norm(x_SVRG - x_star) < tol
        end
        @testset "SVRG++" begin
            solver = CIAOAlgorithms.SVRG{T}(maxit = 16, γ = γ, m = N, plus = true)
            x_SVRG, it_SVRG = solver(x0, F = F, g = g, N = N)
            @test norm(x_SVRG - x_star) < tol
        end

        # test the iterator
        @testset "the iterator" begin
            solver = CIAOAlgorithms.SVRG{T}(γ = γ)
            iter = CIAOAlgorithms.iterator(solver, x0, F = F, g = g, N = N)
            @test iter.x0 === x0

            for state in take(iter, 2)
                @test solution(state) === state.z_full
                @test eltype(solution(state)) == T
            end
            next = iterate(iter) # next = (state, state)
            # one iteration with the solver
            solver = CIAOAlgorithms.SVRG{T}(γ = γ, maxit = 1)
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test solution(next[2]) == x_finito
        end
    end

    @testset "SAGA" begin

        ## test the solver
        @testset "SAGA-Base" begin
            solver = CIAOAlgorithms.SAGA{T}(maxit = maxit)
            x_SAGA, it_SAGA = solver(x0, F = F, g = g, N = N, L = L)
            norm(x_SAGA - x_star) < tol
        end
        @testset "SAGA-stepsize" begin
            γ = 1 / (3 * maximum(L))
            solver = CIAOAlgorithms.SAGA{T}(maxit = maxit, γ = γ)
            x_SAGA, it_SAGA = solver(x0, F = F, g = g, N = N)
            norm(x_SAGA - x_star) < tol
        end

        # test the iterator
        @testset "the iterator" begin
            γ = 1 / (3 * maximum(L))
            solver = CIAOAlgorithms.SAGA{T}(γ = γ)
            iter = CIAOAlgorithms.iterator(solver, x0, F = F, g = g, N = N)
            @test iter.x0 === x0

            for state in take(iter, 2)
                @test solution(state) === state.z
                @test eltype(solution(state)) == T
            end
            next = iterate(iter) # next = (state, state)
            # one iteration with the solver
            solver = CIAOAlgorithms.SAGA{T}(γ = γ, maxit = 1)
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test solution(next[2]) == x_finito
        end
    end

    @testset "SAG" begin
        # note that proximal SAG may not be theoretically convergent
        # maxit = 10000
        ## test the solver
        @testset "SAG-Base" begin
            solver = CIAOAlgorithms.SAG(T, maxit = maxit)
            x_SAG, it_SAG = solver(x0, F = F, g = g, N = N, L = L)
            norm(x_SAG - x_star) < tol
        end
        @testset "SAG" begin
            γ = 1 / (16 * maximum(L))
            solver = CIAOAlgorithms.SAG(T, maxit = maxit, γ = γ)
            x_SAG, it_SAG = solver(x0, F = F, g = g, N = N)
            norm(x_SAG - x_star) < tol
        end

        # test the iterator
        @testset "the iterator" begin
            γ = 1 / (16 * maximum(L))
            solver = CIAOAlgorithms.SAG(T, γ = γ)
            iter = CIAOAlgorithms.iterator(solver, x0, F = F, g = g, N = N)
            @test iter.x0 === x0

            for state in take(iter, 2)
                @test solution(state) === state.z
                @test eltype(solution(state)) == T
            end
            next = iterate(iter) # next = (state, state)
            # one iteration with the solver
            solver = CIAOAlgorithms.SAG(T, γ = γ, maxit = 1)
            x_finito, it_finito = solver(x0, F = F, g = g, L = L, N = N)
            @test solution(next[2]) == x_finito
        end
    end
end


text:
using Revise
using Basins
using Printf
using Plots
using DynamicalSystems

@inline @inbounds function duffing(u, p, t)
    d = p[1]; F = p[2]; omega = p[3]
    du1 = u[2]
    du2 = -d*u[2] + u[1] - u[1]^3 + F*sin(omega*t)
    return SVector{2}(du1, du2)
end

#F=0.3818791946308725; ω= 0.1966442953020134
F=0.2771812080536913; ω=0.1;  # smooth boundary
#ω=0.1617;F = 0.395
#ω=0.3;F = 0.1
ds = ContinuousDynamicalSystem(duffing, rand(2), [0.15, F, ω])
integ_df  = integrator(ds; reltol=1e-8, save_everystep=false)
xg = yg = range(-2.2,2.2,length=100)
@time bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω)

# Basin entropy
@show Sb,Sbb = basin_entropy(bsn; eps_x=20, eps_y=20)

# Wada merge Haussdorff distances
@time max_dist,min_dist = detect_wada_merge_method(xg, yg, bsn)
epsilon = xg[2]-xg[1]
@show dmax = max_dist/epsilon
@show dmin = min_dist/epsilon

# Wada grid
W = detect_wada_grid_method(integ_df, bsn; max_iter=8)
@show W[:,end]

# Uncertainty exponent for these parameter and grid
bd = box_counting_dim(xg, yg, bsn)
α = 2 - bd

D = uncertainty_exponent(bsn, integ_df)
@show 2-D

println("---------------")
println("---------------")
println("Basin Report: ")
println("---------------")
println("---------------")

@printf("Basin entropy %.2f \n", Sb)
@printf("Boundary Basin Entropy: %.2f\n", Sbb)
@printf("Uncertainty exponent: α= %.2f\n", α )
@printf("Box counting dim: bd= %.2f\n", bd)
@printf("Uncertainty dim estimator: d = %.2f\n", 2-D[1])
@printf("Number of basins: %d\n", bsn.Na)
@printf("Merge Method: Max fattening parameter: %.2f\n", dmax)
@printf("Wada Grid Method: W_Na = %.2f\n ", W[end,end] )

plot(xg, yg, bsn.basin', seriestype=:heatmap)


text:
using Test, PrivateExample4

@test hello("Julia") == "Hello, Julia"
@test domath(2.0) ≈ 7.0


text:
### A Pluto.jl notebook ###
# v0.17.2

# using Markdown
# using InteractiveUtils

# ╔═╡ 887893a3-0f17-47c4-ae40-3459fbf2a4f4
#=╠═╡ notebook_exclusive
md"""
# Exports
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 03bea6d7-c397-4a25-a5ff-af370e171f49
#=╠═╡ notebook_exclusive
md"""
# Functions
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 047bc2be-35bc-45e4-bca8-75beeda57057
function _cell_data(filesrc::String)
	names = split(filesrc,"#==#")
	if length(names) > 1
		return names
	else
		return names[1], ""
	end
end

# ╔═╡ 356a8eb0-4f70-11ec-3ca7-efdefb29fe4f
"""
	is_notebook_local(filesrc::AbstractString)::Bool
Returns `true` if Pluto is loaded and if the `filesrc` points to the currently running Pluto cell_id

To be fed the result of `@__FILE__` (or `__source__.file` used inside macros) as `filesrc`.\\
This works because the `@__FILE__` information contains the name of the Pluto notebook followed the string `#==#` and by the cell UUID in case this is called directly in a notebook (and not included from outside)
"""
function is_notebook_local(filesrc)::Bool
	# If this is ran from the notebook it is defined in Pluto, the output of __source__.file or @__FILE__ is split into filepath and cell_id with a #==# separator
	file, cell_id = _cell_data(filesrc)
	return isempty(cell_id) ? false : true
end

# ╔═╡ 35728756-c270-4455-9131-1636614f2466
"""
	only_in_nb(ex)
Executes the expression `ex` only if the macro is called from a running Pluto instance and ran directly from the source notebook file.

This is more strict than `PlutoHooks.@skip_as_script` as including a notebook with `@skip_as_script ex` from another notebook would still execute `ex`.\\
`@only_in_nb ex` instead only evaluates `ex` if the calling notebook is the original source notebook file.

See also: [`@only_out_nb`](@ref). [`PlutoHooks.@skip_as_script`](@ref).
"""
macro only_in_nb(ex)
	is_notebook_local(String(__source__.file::Symbol)) ? esc(ex) : nothing
end

# ╔═╡ fb6ce0fb-bae4-4120-9f33-4edcd36c3c2e
"""
	only_out_nb(ex)
Opposite of `@only_in_nb`

See also: [`@only_in_nb`](@ref). [`PlutoHooks.@only_as_script`](@ref).
"""
macro only_out_nb(ex)
	is_notebook_local(String(__source__.file::Symbol)) ? nothing : esc(ex)
end

# ╔═╡ fa0cf8bd-f783-44cb-9f32-115531adb7d4
"""
	plutodump(x::Union{Symbol, Expr}; maxdepth = 8)
Dumps a symbol or expression directly as text on the pluto cell output instead of in the stdout

See also: [`Meta.dump`](@ref)
"""
function plutodump(x::Union{Symbol, Expr}; maxdepth = 8)
	i = IOBuffer()
	Meta.dump(i, x; maxdepth)
	String(take!(i)) |> Text
end

# ╔═╡ 17f796de-e2db-41a8-9120-509ca54be657
"""
	@current_pluto_cell_id()
Returns the cell_id (as a string) of the cell where the macro is called. If not ran from within the pluto notebook containing the call, returns an empty string
"""
macro current_pluto_cell_id()
	_, cell_id = _cell_data(__source__.file::Symbol |> String)
	return isempty(cell_id) ? "" : cell_id
end

# ╔═╡ 4786625f-39d9-47f7-9f19-2410a91d52db
"""
	@current_pluto_notebook_file()
Returns the path of the notebook file of the cell where the macro is called. If not ran from within the pluto notebook containing the call, returns an empty string
"""
macro current_pluto_notebook_file()
	nbfile, cell_id = _cell_data(__source__.file::Symbol |> String)
	return isempty(cell_id) ? "" : nbfile
end

# ╔═╡ fd5f3b92-a1ac-449c-aafd-d29783152da3
export @only_in_nb, @only_out_nb, plutodump, @current_pluto_cell_id, @current_pluto_notebook_file

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0-rc2"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╟─887893a3-0f17-47c4-ae40-3459fbf2a4f4
# ╠═fd5f3b92-a1ac-449c-aafd-d29783152da3
# ╟─03bea6d7-c397-4a25-a5ff-af370e171f49
# ╠═047bc2be-35bc-45e4-bca8-75beeda57057
# ╠═356a8eb0-4f70-11ec-3ca7-efdefb29fe4f
# ╠═35728756-c270-4455-9131-1636614f2466
# ╠═fb6ce0fb-bae4-4120-9f33-4edcd36c3c2e
# ╠═fa0cf8bd-f783-44cb-9f32-115531adb7d4
# ╠═17f796de-e2db-41a8-9120-509ca54be657
# ╠═4786625f-39d9-47f7-9f19-2410a91d52db
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
sourcedir(res_dir) = Core.sourcedir(Core.ResDirTree(res_dir));
descfile(res_dir) = Core.descfile(Core.ResDirTree(res_dir));
logfile(res_dir) = Core.logfile(Core.ResDirTree(res_dir));
datafile(res_dir) = Core.datafile(Core.ResDirTree(res_dir));

max_stars_repo_name: josePereiro/MyResultsManager

max_issues_repo_name: josePereiro/MyResultsManager

max_forks_repo_name: josePereiro/MyResultsManager


text:
@doc raw"""
    add_variable(canonical::Canonical,
                      devices::D,
                      var_name::Symbol,
                      binary::Bool,
                      expression_name::Symbol,
                      sign::Float64)

Adds a variable to the optimization model and to the affine expressions contained
in the canonical model according to the specified sign. Based on the inputs, the variable can
be specified as binary.

# Bounds

``` lb_value_function <= varstart[name, t] <= ub_value_function ```

If binary = true:

``` varstart[name, t] in {0,1} ```

# LaTeX

``  lb \ge x^{device}_t \le ub \forall t ``

``  x^{device}_t \in {0,1} \forall t iff \text{binary = true}``

# Arguments
* canonical::Canonical : the canonical model built in PowerSimulations
* devices : Vector or Iterator with the devices
* var_name::Symbol : Base Name for the variable
* binary::Bool : Select if the variable is binary
* expression_name::Symbol : Expression_name name stored in canonical.expressions to add the variable
* sign::Float64 : sign of the addition of the variable to the expression_name. Default Value is 1.0

# Accepted Keyword Arguments
* ub_value_function : Provides the function over device to obtain the value for a upper_bound
* lb_value_function : Provides the function over device to obtain the value for a lower_bound. If the variable is meant to be positive define lb = x -> 0.0
* initial_value_function : Provides the function over device to obtain the warm start value

"""
function add_variable(canonical::Canonical,
                      devices::D,
                      var_name::Symbol,
                      binary::Bool,
                      expression_name::Union{Nothing,Symbol}=nothing,
                      sign::Float64=1.0; kwargs...) where {D<:Union{Vector{<:PSY.Device},
                                          IS.FlattenIteratorWrapper{<:PSY.Device}}}

    time_steps = model_time_steps(canonical)
    variable = add_var_container!(canonical, var_name, (PSY.get_name(d) for d in devices), time_steps)
    jvar_name = _remove_underscore(var_name)

    lb_f = get(kwargs, :lb_value, nothing)
    init_f= get(kwargs, :initial_value, nothing)
    ub_f = get(kwargs, :ub_value, nothing)

    for t in time_steps, d in devices
        name = PSY.get_name(d)
        variable[name, t] = JuMP.@variable(canonical.JuMPmodel,
                                        base_name="$(jvar_name)_{$(name), $(t)}",
                                        binary=binary
                                        )

        !isnothing(ub_f) && JuMP.set_upper_bound(variable[name, t], ub_f(d))
        !isnothing(lb_f) && !binary && JuMP.set_lower_bound(variable[name, t], lb_f(d))
        !isnothing(init_f) && JuMP.set_start_value(variable[name, t], init_f(d))

        if !(isnothing(expression_name))
        bus_number = PSY.get_number(PSY.get_bus(d))
        _add_to_expression!(get_expression(canonical, expression_name),
                            bus_number,
                            t,
                            variable[name, t],
                            sign)
        end
    end

    return

end

@doc raw"""
    set_variable_bounds(canonical::Canonical,
                            bounds::Vector{NamedMinMax},
                            var_name::Symbol)

Adds a bounds to a variable in the optimization model.

# Bounds

``` bounds.min <= varstart[name, t] <= bounds.max  ```

# LaTeX

``  x^{device}_t >= bound^{min} \forall t ``

``  x^{device}_t <= bound^{max} \forall t ``

# Arguments
* canonical::Canonical : the canonical model built in PowerSimulations
* bounds::Vector{NamedMinMax} : contains name of device (1) and its min/max (2)
* var_name::Symbol : Base Name for the variable

"""
function set_variable_bounds(canonical::Canonical,
                            bounds::Vector{NamedMinMax},
                            var_name::Symbol)

    for t in model_time_steps(canonical), (name, bound) in bounds
        var = canonical.variables[var_name][name, t]
        JuMP.set_upper_bound(var, bound.max)
        JuMP.set_lower_bound(var, bound.min)
    end

end


text:
export AbstractHLADataset, HLADataset

abstract type AbstractHLADataset end

struct HLADataset <: AbstractHLADataset
    name::String
    data::Vector{<: AbstractHLAData}
end

function Escape.HLADataset(::Val{:Test})
    sequence_dir = joinpath(@__DIR__, "..", "data", "test_sequences")
    fasta_path = joinpath(sequence_dir, "test_large.fasta")

    hla_types = rand(HLAType, 15)
    hla_types[5] = HLAType((parse_allele("A11"), parse_allele("A13"),
                            parse_allele("B53"), parse_allele("B56"),
                            missing, missing))
    hla_types[6] = HLAType((missing, missing, missing, missing, missing, missing))
    hla_data_1 = HLAData("protein_1", fasta_path, hla_types, missing, missing)
    hla_data_2 = HLAData("protein_2", fasta_path, rand(HLAType, 15), missing, missing)

    if !isfile(joinpath(sequence_dir, "test_input_1.jls"))
        input = Escape.compute_stan_input(hla_data_1)
        serialize(joinpath(sequence_dir, "test_input_1.jls"), input)
    end

    if !isfile(joinpath(sequence_dir, "test_input_2.jls"))
        input = Escape.compute_stan_input(hla_data_2)
        serialize(joinpath(sequence_dir, "test_input_2.jls"), input)
    end

    input_1 = deserialize(joinpath(sequence_dir, "test_input_1.jls"))
    input_2 = deserialize(joinpath(sequence_dir, "test_input_2.jls"))

    hla_data_1 = HLAData("protein_1", fasta_path, hla_types, missing, input_1)
    hla_data_2 = HLAData("protein_2", fasta_path, rand(HLAType, 15), missing, input_2)

    dataset = HLADataset("Test", [hla_data_1, hla_data_2])

    return dataset
end


text:
using FractionalCalculus
using Test

@testset "Test Fractional Integral" begin
    @test isapprox(fracint(x->x^4, 0.5, 0, 1.23, 1e-8, RL_Direct())[1], 1.163931646862977; atol = 0.00001)
    @test isapprox(fracint(x->x^4, 0.25, 0, 1.23, 1e-8, RL_Direct())[1], 1.64294134724409; atol = 0.00001)
    @test isapprox(fracint(x->x^4, 0.15, 0, 1.23, 1e-8, RL_Direct())[1], 1.878978485456182; atol = 0.00001)
    @test isapprox(fracint(x->sin(x), 0.5, 0, 0.97, 1e-8, RL_Direct())[1], 0.6442601778843754; atol = 0.00001)
    @test isapprox(fracint(x->sin(x), 0.25, 0, 0.97, 1e-8, RL_Direct())[1], 0.7448298999485293; atol = 0.00001)
    @test isapprox(fracint(x->sin(x), 0.15, 0, 0.97, 1e-8, RL_Direct())[1], 0.7802555196986999; atol = 0.00001)
end

@testset "Test Piecewise Fractional Integral" begin
    @test isapprox(fracint(x->x^4, 0.5, 1.23, 1e-4, RL_Piecewise()), 1.163931646862977; atol = 1e-5)
    @test isapprox(fracint(x->x^4, 0.25, 1.23, 1e-4, RL_Piecewise()), 1.64294134724409; atol = 1e-5)
    @test isapprox(fracint(sin, 0.5, 0.97, 1e-4, RL_Piecewise()), 0.6442601778843754; atol = 1e-5)
end


text:
using Revise, ApproxOperator

elements, nodes = importmsh("./msh/patchtest.msh")
nₚ = length(nodes[:x])
nₑ = length(elements["Ω"])

elements["Γ̃ᵍ"] = Element{:Seg2}(elements["Γᵍ"])
set𝓖!(elements["Ω"],:TriGI3)
set𝓖!(elements["Γᵍ"],:SegGI2)
set𝓖!(elements["Γ̃ᵍ"],:SegGI2)
elements["Γ̃ᵍ"] = Element{:Tri3}(elements["Ω"],elements["Γ̃ᵍ"])

prescribe!(elements["Ω"],:b,(x,y,z)->0.0)
# prescribe!(elements["Γᵗ₁"],:t₁,(x,y,z)->E/(1-ν))
# prescribe!(elements["Γᵗ₁"],:t₂,(x,y,z)->E/(1+ν))
# prescribe!(elements["Γᵗ₂"],:t₁,(x,y,z)->E/(1+ν))
# prescribe!(elements["Γᵗ₂"],:t₂,(x,y,z)->E/(1-ν))
prescribe!(elements["Γᵍ"],:g,(x,y,z)->1.0+x+y)
prescribe!(elements["Γ̃ᵍ"],:g,(x,y,z)->1.0+x+y)

coefficient = (:k=>1.0,:α=>0e0)
ops = [Operator(:∫∇v∇udΩ,coefficient...),
       Operator(:∫vbdΩ,coefficient...),
       Operator(:∫vtds,coefficient...),
       Operator(:∫∇𝑛vgdΓ,coefficient...),
       Operator(:∫vgdΓ,coefficient...),
       Operator(:H₁,coefficient...)]

k = zeros(nₚ,nₚ)
f = zeros(nₚ)

ops[1](elements["Ω"],k)
ops[4](elements["Γ̃ᵍ"],k,f)
ops[5](elements["Γᵍ"],k,f)

x = nodes[:x]
y = nodes[:y]
test = y'*f

d = k\f

push!(nodes,:d=>d)
prescribe!(elements["Ω"],:u,(x,y,z)->1.0+x+y)
prescribe!(elements["Ω"],:∂u∂x,(x,y,z)->1.0)
prescribe!(elements["Ω"],:∂u∂y,(x,y,z)->1.0)
h1,l2 = ops[6](elements["Ω"])

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
module BeamFEA

import LinearAlgebra: Symmetric, eigvals
import Interpolations: LinearInterpolation

export BeamProperties, Loads, NoLoads, Stiffness, Deflection
export fea, strain

const DOF = 6

"""
wrapper for linear interpolation
"""
function linterp(x, y, xpt)
    intp = LinearInterpolation(x, y)
    return intp(xpt)
end

"""
(private)
see doc: stiffness matrix
"""
function bendingstiffness(EI, L)

    e1 = EI[1]
    e2 = EI[2]

    k = [
        (6*e1 + 6*e2)/L^3   (4*e1 + 2*e2)/L^2   (-6*e1 - 6*e2)/L^3  (2*e1 + 4*e2)/L^2;
        0.0                 (3*e1 + e2)/L     (-4*e1 - 2*e2)/L^2  (e1 + e2)/L;
        0.0                 0.0                 (6*e1 + 6*e2)/L^3   (-2*e1 - 4*e2)/L^2;
        0.0                 0.0                 0.0                 (e1 + 3*e2)/L;
    ]

    return Symmetric(k)
end

"""
(private)
see doc: stiffness matrix
"""
function axialstiffness(EA, L)

    k = (EA[1] + EA[2])/(2*L)*[
        1.0 -1.0;
        -1.0  1.0
    ]

    return Symmetric(k)
end

"""
(private)
see doc: Distributed Loads
"""
function bendingloads(P, L)

    p1 = P[1]
    p2 = P[2]

    f = L*[
        7.0/20*p1 + 3.0/20*p2;
        L*(1.0/20*p1 + 1.0/30*p2);
        3.0/20*p1 + 7.0/20*p2;
        L*(-1.0/30*p1 - 1.0/20*p2);
    ]

    return f
end

"""
(private)
see doc: Distributed Loads
"""
function axialloads(P, L)

    p1 = P[1]
    p2 = P[2]

    f = L*[
        p1/3.0 + p2/6.0;
        p1/6.0 + p2/3.0;
    ]

    return f
end

"""
(private)
see doc: Inertial Matrix
"""
function bendinginertial(rhoA, L)

    a1 = rhoA[1]
    a2 = rhoA[2]

    m = 1.0/420*[
        (120*a1 + 36*a2)*L  (15*a1 + 7*a2)*L^2      (27*a1 + 27*a2)*L   (-7*a1 - 6*a2)*L^2;
        0.0                 (2.5*a1 + 1.5*a2)*L^3   (6*a1 + 7*a2)*L^2   (-1.5*a1 - 1.5*a2)*L^3;
        0.0                 0.0                     (36*a1 + 120*a2)*L  (-7*a1 - 15*a2)*L^2;
        0.0                 0.0                     0.0                 (1.5*a1 + 2.5*a2)*L^3;
    ]

    return Symmetric(m)
end

"""
(private)
see doc: Inertial Matrix
"""
function axialinertial(rhoA, L)

    r1 = rhoA[1]
    r2 = rhoA[2]

    m = L*[
        (1.0/4*r1 + 1.0/12*r2)   (1.0/12*r1 + 1.0/12*r2);
        (1.0/12*r1 + 1.0/12*r2)  (1.0/12*r1 + 1.0/4*r2)
    ]

    return m
end

"""
(private)
computes FEM matrices for one 12-dof beam element

q = [x1, thetax1, y1, thetay1, z1, thetaz1, ... repeat for 2]
x is axial
"""
function beam_matrix(L, EIy, EIz, EA, GJ, rhoA, rhoJ, Px, Py, Pz)

    # initialize
    type = typeof(EIy[1])
    K = zeros(type, 2*DOF, 2*DOF)  # stiffness matrix
    M = zeros(type, 2*DOF, 2*DOF)  # inertial matrix
    F = zeros(type, 2*DOF)  # force vector

    # --- axial ----
    idx = [1, 7]

    kx = axialstiffness(EA, L)
    K[idx, idx] = kx

    mx = axialinertial(rhoA, L)
    M[idx, idx] = mx

    fx = axialloads(Px, L)
    F[idx] = fx

    # --- torsion ---
    idx = [2, 8]

    kt = axialstiffness(GJ, L)
    K[idx, idx] = kt

    mt = axialinertial(rhoJ, L)
    M[idx, idx] = mt

    # no distributed torsional loads

    # --- bending in y ---
    idx = [3, 4, 9, 10]

    ky = bendingstiffness(EIy, L)
    K[idx, idx] = ky

    my = bendinginertial(rhoA, L)
    M[idx, idx] = my

    fy = bendingloads(Py, L)
    F[idx] = fy

    # --- bending in z ---
    idx = [5, 6, 11, 12]

    kz = bendingstiffness(EIz, L)
    K[idx, idx] = kz

    M[idx, idx] = my  # inertial is same in y and z

    fz = bendingloads(Pz, L)
    F[idx] = fz

    return K, M, F
end

struct BeamProperties{TF, TA<:AbstractVector{TF}}

    x::TA  # sometimes this is a range
    EIy::Vector{TF}
    EIz::Vector{TF}
    EA::Vector{TF}
    GJ::Vector{TF}
    rhoA::Vector{TF}
    rhoJ::Vector{TF}

    # EIx, EA etc. are arrays of length nodes
end

struct Loads{TF, TA<:AbstractVector{TF}}

    x::TA
    Px::Vector{TF}
    Py::Vector{TF}
    Pz::Vector{TF}
    Fx::Vector{TF}
    Fy::Vector{TF}
    Fz::Vector{TF}
    Mx::Vector{TF}
    My::Vector{TF}
    Mz::Vector{TF}
end

NoLoads(x) = Loads(x, zeros(length(x)), zeros(length(x)), zeros(length(x)), zeros(length(x)),
    zeros(length(x)), zeros(length(x)), zeros(length(x)), zeros(length(x)), zeros(length(x)))

# addedK: additional stiffness (or infinite stiffness for a rigid connection)
struct Stiffness{TF}

    kx::Vector{TF}
    ky::Vector{TF}
    kz::Vector{TF}
    kthetax::Vector{TF}
    kthetay::Vector{TF}
    kthetaz::Vector{TF}
end

struct Deflection{TF}
    x::Vector{TF}
    y::Vector{TF}
    z::Vector{TF}
    thetax::Vector{TF}
    thetay::Vector{TF}
    thetaz::Vector{TF}
end

struct ShearBending{TF}
    Nx::Vector{TF}
    Vy::Vector{TF}
    Vz::Vector{TF}
    Tx::Vector{TF}
    My::Vector{TF}
    Mz::Vector{TF}
end

"""
    fea_analysis(beam, loads, stiffness)

Beam finite element analysis.

**Inputs**
- beam::BeamProperties
- loads::Loads
- stiffness::Stiffness

**Outputs**
- defl::Deflection
- freq: natural frequencies
"""
function fea(beam, loads, stiffness)

    nodes = length(beam.x)
    elements = nodes - 1

    # ---- interpolate loads onto beam ----
    Px = linterp(loads.x, loads.Px, beam.x)
    Py = linterp(loads.x, loads.Py, beam.x)
    Pz = linterp(loads.x, loads.Pz, beam.x)
    Fx = linterp(loads.x, loads.Fx, beam.x)
    Fy = linterp(loads.x, loads.Fy, beam.x)
    Fz = linterp(loads.x, loads.Fz, beam.x)
    Mx = linterp(loads.x, loads.Mx, beam.x)
    My = linterp(loads.x, loads.My, beam.x)
    Mz = linterp(loads.x, loads.Mz, beam.x)
    # ----------------------------------

    # ---- global to fea convention ----
    My, Mz = -Mz, My
    # EIy/z and kthetay/z changes done in place
    # ---------------------------------

    # --- assemble global matrices -----
    type = typeof(beam.EIy[1])
    K = zeros(type, DOF*nodes, DOF*nodes)
    M = zeros(type, DOF*nodes, DOF*nodes)
    F = zeros(type, DOF*nodes)

    # pull out submatrix for each element
    for i = 1:elements

        Ksub, Msub, Fsub = beam_matrix(beam.x[i+1] - beam.x[i],
            beam.EIz[i:i+1], beam.EIy[i:i+1], # swapped y and z (global -> fea convention)
            beam.EA[i:i+1], beam.GJ[i:i+1], beam.rhoA[i:i+1], beam.rhoJ[i:i+1],
            Px[i:i+1], Py[i:i+1], Pz[i:i+1])

        start = (i-1)*DOF
        idx = start+1:start+2*DOF
        K[idx, idx] += Ksub
        M[idx, idx] += Msub
        F[idx] += Fsub
    end

    # ---- add point loads ----
    for i = 1:nodes
        start = (i-1)*DOF

        F[start + 1] += Fx[i]
        F[start + 2] += Mx[i]
        F[start + 3] += Fy[i]
        F[start + 4] += My[i]
        F[start + 5] += Fz[i]
        F[start + 6] += Mz[i]
    end

    # ---- add point inertial elements -----

    # --- apply boundary conditions ----

    # add additional stiffness (and save nodes to remove)
    save = trues(DOF*nodes)  # nodes to keep
    for i = 1:nodes

        start = (i-1)*DOF

        kvec = [stiffness.kx[i], stiffness.kthetax[i], stiffness.ky[i], stiffness.kthetaz[i],  # swapped y and z (global -> fea convention)
            stiffness.kz[i], stiffness.kthetay[i]]  # order for convenience

        for j = 1:DOF  # iterate through each DOF

            if kvec[j] == Inf  # if rigid, remove this index
                save[start+j] = false
            else  # otherwise add directly to stiffness matrix
                K[start+j, start+j] += kvec[j]
            end
        end
    end

    K = K[save, save]
    M = M[save, save]
    F = F[save]

    # denote symmeric more efficient linear solve and eigenanalysis
    K = Symmetric(K)
    M = Symmetric(M)

    # ---- compute deflections -----
    # initialize to zero so we can keep rigid deflections
    delta = zeros(DOF*nodes)

    try
        deltasub = K\F
        delta[save] = deltasub  # insert nonzero deflections
    catch err
        @warn "structure improperly constrained if deflections are desired."
    end

    defl = Deflection(delta[1:DOF:end], delta[3:DOF:end], delta[5:DOF:end],
        delta[2:DOF:end], delta[6:DOF:end], -delta[4:DOF:end])  # swapped y and z (fea -> global convention)

    # ----- compute eigenvalues -----
    lambda = eigvals(K, M)  # eigenvalues are omega^2 (currently not using eigenvectors)

    # don't save rigid modes
    lambda = lambda[lambda .> 1e-6]

    # convert to freq
    freq = sqrt.(lambda) / (2*pi)

    return defl, freq
end

"""
Currently it assumes that the free end is at index 1.  TODO: allow either.
"""
function strain(y, z, beam, loads)

    # initialize
    n = length(beam.x)
    type = typeof(beam.EIy[1])
    Nx = Vector{type}(undef, n)
    Vy = Vector{type}(undef, n)
    Vz = Vector{type}(undef, n)
    Tx = Vector{type}(undef, n)
    My = Vector{type}(undef, n)
    Mz = Vector{type}(undef, n)

    # unpack
    x = beam.x
    Px = loads.Px
    Py = loads.Py
    Pz = loads.Pz
    Fxpt = loads.Fx
    Fypt = loads.Fy
    Fzpt = loads.Fz
    Mxpt = loads.Mx
    Mypt = loads.My
    Mzpt = loads.Mz

    # integrate
    for i = 1:n-1
        Nx[i+1] = Nx[i] + Fxpt[i] + (x[i+1] - x[i])*(Px[i] + Px[i+1])/2.0
        Vy[i+1] = Vy[i] + Fypt[i] + (x[i+1] - x[i])*(Py[i] + Py[i+1])/2.0
        Vz[i+1] = Vz[i] + Fzpt[i] + (x[i+1] - x[i])*(Pz[i] + Pz[i+1])/2.0

        Tx[i+1] = Tx[i] + Mxpt[i]
        My[i+1] = My[i] + Mypt[i] - (x[i+1] - x[i])*(Vz[i] + Fzpt[i]) - (x[i+1] - x[i])^2*(2*Pz[i] + Pz[i+1])/6.0
        Mz[i+1] = Mz[i] + Mzpt[i] + (x[i+1] - x[i])*(Vy[i] + Fypt[i]) + (x[i+1] - x[i])^2*(2*Py[i] + Py[i+1])/6.0
    end

    # strain
    epsilon_axial = Tx./beam.EA - Mz./beam.EIz.*y + My./beam.EIy.*z

    return epsilon_axial, ShearBending(Nx, Vy, Vz, Tx, My, Mz)
end

end # module


text:
using VectorizationBase, LoopVectorization
using VectorizationBase: data

# @generated to use VectorizationBase's API for supporting 1.5 and 1.6+
@generated function readcyclecounter()
    decl = "declare i64 @llvm.readcyclecounter()"
    instr = "%res = call i64 @llvm.readcyclecounter()\nret i64 %res"
    VectorizationBase.llvmcall_expr(decl, instr, :Int64, :(Tuple{}), "i64", String[], Symbol[])
end

@generated function volatile(x::Vec{W,T}) where {W,T}
    typ = VectorizationBase.LLVM_TYPES[T]
    vtyp = "<$W x $typ>"

    suffix = T == Float32 ? "ps" : "pd"
    sideeffect_str = """%res = call <$W x $(typ)> asm sideeffect "", "=v,v"(<$W x $(typ)> %0)
                               ret <$W x $(typ)> %res"""
    quote
        $(Expr(:meta, :inline))
        Vec(Base.llvmcall($sideeffect_str, NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}}}, VectorizationBase.data(x)))
    end
end
@inline volatile(x::VecUnroll) = VecUnroll(VectorizationBase.fmap(volatile, data(x)))
@inline volatile(x::Tuple) = map(volatile, x)
# @generated function volatile(x::Vec{W,T}, x::Vec{W,T}) where {W,T}
#     typ = VectorizationBase.LLVM_TYPES[T]
#     vtyp = "<$W x $typ>"

#     suffix = T == Float32 ? "ps" : "pd"
#     sideeffect_str = """%res = call <$W x $(typ)> asm sideeffect "", "=v,v"(<$W x $(typ)> %0, <$W x $(typ)> %1)
#                                ret <$W x $(typ)> %res"""
#     quote
#         $(Expr(:meta, :inline))
#         Vec(Base.llvmcall($sideeffect_str, NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}}}, VectorizationBase.data(x)))
#     end
# end

num_vectors(::VecUnroll{N}) where {N} = N+1
num_vectors(::Vec) = 1
function unrolltest(f::F, vs::Vararg{Any,K}) where {F,K}
    cc = readcyclecounter()
    # num_iter = 1_048_576
    num_iter = 4_194_304
    for i ∈ 1:num_iter
        volatile(f(map(volatile, vs)...))
    end
    cycles = readcyclecounter() - cc
    cycles / (num_vectors(first(vs)) * num_iter)
end

# @generated function vapply!(f::F, y, x, ::Val{U}) where {F,U}
#     quote
#         @avx unroll=$U for j ∈ 1:1024
#             y[j] = f(x[j])
#         end
#     end
# end

# vector_init(::Val{N}, ::Type{T}) where {N,T} = VectorizationBase.zero_vecunroll(StaticInt(N), pick_vector_width(T), T, VectorizationBase.register_size())
# vector_init(::Val{1}, ::Type{T}) where {T} = VectorizationBase.vzero(pick_vector_width(T), T)

# @generated function unrolltest(f::F, x::AbstractVector{T}, ::Val{U}) where {F,U,T}
#     quote
#         cc = readcyclecounter()
#         for i ∈ 1:8192
#             s = vector_init(Val{$U}(), $T)
#             @avx unroll=$U for j ∈ 1:512
#                 s += f(x[j])
#             end
#             volatile(s)
#         end
#         cycles = readcyclecounter() - cc
#         pick_vector_width(T) * cycles / (512 * 8192)
#     end
# end

# @generated function unrolltest!(f::F, y::AbstractVector{T}, x::AbstractVector{T}, ::Val{U}) where {F,U,T}
#     quote
#         cc = readcyclecounter()
#         for i ∈ 1:8192
#             @avx unroll=$U for j ∈ 1:512
#                 y[j] = f(x[j])
#             end
#         end
#         cycles = readcyclecounter() - cc
#         pick_vector_width(T) * cycles / (512 * 8192)
#     end
# end

let
    vx = Vec(ntuple(_ -> 10randn(), pick_vector_width(Float64))...);
    vu2 = VectorizationBase.VecUnroll(ntuple(_ -> Vec(ntuple(_ -> 10randn(), pick_vector_width(Float64))...), Val(2)));
    vu4 = VectorizationBase.VecUnroll(ntuple(_ -> Vec(ntuple(_ -> 10randn(), pick_vector_width(Float64))...), Val(4)));
    vu8 = VectorizationBase.VecUnroll(ntuple(_ -> Vec(ntuple(_ -> 10randn(), pick_vector_width(Float64))...), Val(8)));
    for unaryf ∈ [log, log2, log10, log1p, exp, exp2, exp10, expm1, sin, cos]
        rt1 = unrolltest(f, vx)
        rt2 = unrolltest(f, vu2)
        rt4 = unrolltest(f, vu4)
        rt8 = unrolltest(f, vu8)
    end
    for binaryf ∈ [+, *, ^]
        rt1 = unrolltest(f, vx, vx)
        rt2 = unrolltest(f, vu2, vu2)
        rt4 = unrolltest(f, vu4, vu4)
        rt8 = unrolltest(f, vu8, vu8)
    end
end

let
    f, io = mktemp()
    W = Int(VectorizationBase.pick_vector_width(Float64))
    code_native(io, exp, (VecUnroll{1,W,Float64,Vec{W,Float64}},); debuginfo=:none)
    close(io)
    run(`llvm-mca -mcpu=$(Sys.CPU_NAME) -output-asm-variant=1 -bottleneck-analysis $f`)
end


text:
using Libdl

function parsehtml(input::AbstractString; strict=false, preserve_whitespace=false)
    result_ptr = ccall((:gumbo_parse,libgumbo),
                       Ptr{CGumbo.Output},
                       (Cstring,),
                       input)
    goutput::CGumbo.Output = unsafe_load(result_ptr)

    if strict && goutput.errors.length > 0
        throw(InvalidHTMLException("input html was invalid"))
    end
    doc = document_from_gumbo(goutput, preserve_whitespace)
    default_options = Libdl.dlsym(Gumbo_jll.libgumbo_handle, :kGumboDefaultOptions)

    ccall((:gumbo_destroy_output,libgumbo),
          Cvoid,
          (Ptr{Cvoid}, Ptr{CGumbo.Output}),
          default_options, result_ptr)
    return doc
end

# turn a gumbo vector into a Julia vector
# of Ptr{T} where T is the struct contained
# by the gumbo vector
gvector_to_jl(T,gv::CGumbo.Vector) = unsafe_wrap(Array, convert(Ptr{Ptr{T}},gv.data),
                                                     (Int(gv.length),))

# convert a vector of pointers to GumboAttributes to
# a Dict AbstractString => AbstractString
function attributes(av::Vector{Ptr{CGumbo.Attribute}})
    result = Dict{AbstractString,AbstractString}()
    for ptr in av
        ga::CGumbo.Attribute = unsafe_load(ptr)
        result[unsafe_string(ga.name)] = unsafe_string(ga.value)
    end
    return result
end

function elem_tag(ge::CGumbo.Element)
    tag = CGumbo.TAGS[ge.tag+1]  # +1 is for 1-based julia indexing
    if tag == :unknown
        ot = ge.original_tag
        tag = split(unsafe_string(ot.data, ot.length)[2:end-1])[1] |> Symbol
    end
    tag
end

function gumbo_to_jl(parent::HTMLNode, ge::CGumbo.Element, preserve_whitespace)
    tag = elem_tag(ge)
    attrs = attributes(gvector_to_jl(CGumbo.Attribute,ge.attributes))
    children = HTMLNode[]
    res = HTMLElement{tag}(children, parent, attrs)
    for childptr in gvector_to_jl(CGumbo.Node{Int},ge.children)
        node = load_node(childptr, preserve_whitespace)
        if in(typeof(node).parameters[1], [CGumbo.Element, CGumbo.Text])
            push!(children, gumbo_to_jl(res, node.v, preserve_whitespace))
        end
    end
    res
end

function gumbo_to_jl(parent::HTMLNode, gt::CGumbo.Text, preserve_whitespace)
    HTMLText(parent, unsafe_string(gt.text))
end

# this is a fallback method that should only be called to construct
# the root of a tree
gumbo_to_jl(ge::CGumbo.Element, preserve_whitespace) = gumbo_to_jl(NullNode(), ge, preserve_whitespace)

# load a GumboNode struct into memory as the appropriate Julia type
# this involves loading it once as a CGumbo.Node{Int} in order to
# figure out what the correct type actually is, and then reloading it as
# that type
function load_node(nodeptr::Ptr, preserve_whitespace=false)
    precursor = unsafe_load(reinterpret(Ptr{CGumbo.Node{Int}},nodeptr))
    # TODO clean this up with a Dict in the CGumbo module
    correctptr = if precursor.gntype == CGumbo.ELEMENT
        reinterpret(Ptr{CGumbo.Node{CGumbo.Element}},nodeptr)
    elseif precursor.gntype == CGumbo.TEXT
        reinterpret(Ptr{CGumbo.Node{CGumbo.Text}},nodeptr)
    elseif precursor.gntype == CGumbo.DOCUMENT
        reinterpret(Ptr{CGumbo.Node{CGumbo.Document}},nodeptr)
    elseif preserve_whitespace && precursor.gntype == CGumbo.WHITESPACE
        reinterpret(Ptr{CGumbo.Node{CGumbo.Text}},nodeptr)
    else
        # TODO this is super sketchy and should realistically be an
        # error
        nodeptr
    end
    unsafe_load(correctptr)
end

# transform gumbo output into Julia data
function document_from_gumbo(goutput::CGumbo.Output, preserve_whitespace)
    # TODO convert some of these typeasserts to better error messages?
    gnode::CGumbo.Node{CGumbo.Document} = load_node(goutput.document, preserve_whitespace)
    gdoc = gnode.v
    doctype = unsafe_string(gdoc.name)
    groot::CGumbo.Node{CGumbo.Element} = load_node(goutput.root, preserve_whitespace)
    root = gumbo_to_jl(groot.v, preserve_whitespace)  # already an element
    HTMLDocument(doctype, root)
end


text:
@testset "propagate.jl" begin
    println("             > propagate.jl")
    @testset "GPU" begin
        m_out = 1024
        m_in  = 1024

        state = BinomialState(128, m_out, m_in)
        model = BinomialModel(128, m_out)

        for i in 1:100
            propagate!(state, model, 0.1f0)
            @test minimum(state.n) >= 0
            @test minimum(state.k) >= 0
            @test minimum(state.n .- state.k) >= 0
            @test minimum(state.n .- state.n) >= 0
        end

        if RUN_BENCHMARKS
            println("")
            println("Benchmarking function propagate!: should take about 4ms")
            display(@benchmark CUDA.@sync propagate!($state, $model, 0.1f0))
            println("")
        end

        dt    = 0.1f0
        m_out = 4
        m_in  = 4
        N     = 10 .* CUDA.ones(Int, m_out)
        p     = cu([0f0, 0f0, 1f0, 1f0])
        q     = CUDA.rand(m_out)
        σ     = CUDA.rand(m_out)
        τ     = cu([1f3, 1f-4, 1f-4, 1f3])
        model = BinomialModel(N, p, q, σ, τ);
        n     = 5 .* CUDA.ones(Int, m_out, m_in)
        k     = 4 .* CUDA.ones(Int, m_out, m_in)

        state = BinomialState(n, k)

        propagate!(state, model, dt)
        n = Array(state.n)
        k = Array(state.k)

        # When τ is much larger than dt, no vesicle is refilled
        @test all(n[1,:] .== 1)
        @test all(n[4,:] .== 1)
        # When τ is much shorter than dt, all vesicles are refilled
        @test all(n[2,:] .== 10)
        @test all(n[3,:] .== 10)

        # When the release probability is 0, no vesicle is released
        @test all(k[1,:] .== 0)
        @test all(k[2,:] .== 0)
        # When the release probability is 1, all vesicles are released
        @test all(k[3,:] .== 10)
        @test all(k[4,:] .== 1)
    end#GPU

    @testset "CPU" begin
        m_out = 16
        m_in  = 16

        state = BinomialState(128, m_out, m_in, :cpu)
        model = BinomialModel(128, m_out, :cpu)

        for i in 1:10
            propagate!(state, model, 0.1)
            @test minimum(state.n) >= 0
            @test minimum(state.k) >= 0
            @test minimum(state.n .- state.k) >= 0
            @test minimum(state.n .- state.n) >= 0
        end

        dt    = 0.1
        m_out = 4
        m_in  = 4
        N     = 10 .* ones(Int, m_out)
        p     = [0., 0., 1., 1.]
        q     = rand(m_out)
        σ     = rand(m_out)
        τ     = [1e3, 1e-4, 1e-4, 1e3]
        model = BinomialModel(N, p, q, σ, τ);
        n     = 5 .* ones(Int, m_out, m_in)
        k     = 4 .* ones(Int, m_out, m_in)

        state = BinomialState(n, k)

        propagate!(state, model, dt)
        n = Array(state.n)
        k = Array(state.k)

        # When τ is much larger than dt, no vesicle is refilled
        @test all(n[1,:] .== 1)
        @test all(n[4,:] .== 1)
        # When τ is much shorter than dt, all vesicles are refilled
        @test all(n[2,:] .== 10)
        @test all(n[3,:] .== 10)

        # When the release probability is 0, no vesicle is released
        @test all(k[1,:] .== 0)
        @test all(k[2,:] .== 0)
        # When the release probability is 1, all vesicles are released
        @test all(k[3,:] .== 10)
        @test all(k[4,:] .== 1)
    end#CPU
end


text:
@testset "Linear disjoint" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  _K, _ = NumberField([x^2 - 2, x^2 - 3], "a", cached = false)
  K, _ = simple_extension(_K)
  L, b = NumberField(x^2 - 2, "b", cached = false)
  @test !islinearly_disjoint(K, L)

  M, c = NumberField(x^2 - 3, "c", cached = false)
  @test islinearly_disjoint(L, M)
end

@testset "Random" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  K, a = NumberField(x^32 + 2, "a")

  b = @inferred rand([a], -10:10)
  @test b isa nf_elem
  @test_throws ErrorException rand(nf_elem[], -10:10)

  b = @inferred rand(basis(K), 1:100, 10)
  @test count(!iszero, (coeff(b, i) for i in 0:31)) <= 10
  @test_throws ErrorException rand(nf_elem[], -10:10, 5)
  @test_throws ErrorException rand([a, a^2], -10:10, -10)
  @test_throws ErrorException rand(basis(K), -10:10, 100)

  @inferred rand!(b, basis(K), 1:100, 20)
  @test count(!iszero, (coeff(b, i) for i in 0:31)) <= 20
  @test_throws ErrorException rand!(b, basis(K), 1:100, 120)
  @test_throws ErrorException rand!(b, basis(K), 1:100, -100)

  @inferred rand!(b, basis(K), 1:100)
  @test_throws ErrorException rand!(b, nf_elem[], 1:100)
end

@testset "Polynomial" begin
  Qx, x = QQ["x"]
  K, _a = NumberField(x^3 - 3*x - 1, "a")
  Kt, t = K["t"]
  f = t^4+(-28*_a^2 + 26*_a + 124)*t^2+(81*_a^2 + 504*_a + 936)
  @test @inferred isirreducible(f)

  f = x^3-7*x^2+6*x-1
  K, a = number_field(f, "a")
  Ky, y = K["y"]
  h = y^3+(15037//140*a^2 - 109//40*a - 915//14)*y^2+(16375724527//78400*a^2 - 993527643//4900*a + 393774209//11200)*y+(107296943419277//878080*a^2 - 2594040461688323//21952000*a + 17784340885567//878080)
  @assert isirreducible(h)
end

@testset "Is integral" begin
  Qx, x = FlintQQ["x"]
  f = x^2 + 1
  K, a = number_field(f, "a")

  @test Hecke.isintegral(a) == true
  @test Hecke.isintegral(fmpq(1, 2)*a) == false

  g = x^3 + 3
  L, b = number_field([f, g], "b")

  @test Hecke.isintegral(b[1]) == true
  @test Hecke.isintegral(fmpq(1, 2)*b[1]) == false
end

@testset "Compositum" begin
  Qx, x = FlintQQ["x"]
  f = x^2 + 1
  K, a = number_field(f, "a")
  L, b = number_field(x^2-3, "b")
  C, mK, mL = compositum(K, L)

  @test iszero(K.pol(mK(gen(K))))
  @test iszero(L.pol(mL(gen(L))))
end

@testset "PolyFactor" begin
  Zx, x = FlintZZ["x"]
  k, a = number_field(swinnerton_dyer(3, x))
  kt, t = k["t"]

  g = swinnerton_dyer(8, x)
  @test length(factor((t^2-a)*(t^3-a-1))) == 2 #Trager
  @test length(factor((t^2-a)*(t^3-a-1)*(t+a^2+1)*(t^5+t+a))) == 4 #Zassenhaus
  @test length(factor(change_base_ring(k, g))) == 8 # van Hoeij

  K, a = NumberField(x - 1, "a")
  Kt, t = K["t"]
  f = t^5 -3 * t^4 - 104 * t^3 + 312 * t^2 + 400t -1200
  @test length(factor(f)) == 5
  @test length(factor(f*t)) == 6

  #Tommys
  K, a = number_field(x^2 - x - 4)
  Ky, y = K["y"]
  f = y^16+(39)*y^14+(449)*y^12+(1794)*y^10+(2830)*y^8+(1794)*y^6+(449)*y^4+(39)*y^2+(1)
  @test length(factor(f)) == 2
end


text:
# Start Test Script
using ILUZero, LinearAlgebra, SparseArrays, Test

function cg(A, b; M=I)
    n = size(A, 1)
    x = zeros(n)
    r = copy(b)
    z = zeros(n)
    ldiv!(z, M, r)
    p = copy(z)
    γ = dot(r, z)
    k = 0
    tired = false
    solved = false
    while !(solved || tired)
        Ap = A * p
        α = γ / dot(p, Ap)
        x = x + α * p
        r = r - α * Ap
        ldiv!(z, M, r)
        γ_new = dot(r, z)
        β = γ_new / γ
        γ = γ_new
        p = z + β * p
        k = k + 1
        tired = k > n
        solved = norm(r) ≤ 1e-8 * norm(b)
    end
    return x, k
end

function test_solve()
    n = 100
    tA = sprandn(n, n, .1) + 10.0 * I
    A = tA' * tA
    ilu_prec = ilu0(A)
    b = rand(n)
    x1, k1 = cg(A, b)
    println("No preconditioning: ", k1, " iterations")
    x2, k2 = cg(A, b, M=ilu_prec)
    println("Preconditioned: ", k2, " iterations")
    return k1 > k2
end

function test_substitutions()
    n = 100
    tA = sprandn(n, n, .1) + 10.0 * I
    A = tA' * tA
    ilu_prec = ilu0(A)
    b = rand(n)
    x = ilu_prec \ b
    x1 = zeros(n)
    x2 = zeros(n)
    forward_substitution!(x1, ilu_prec, b)
    backward_substitution!(x2, ilu_prec, x1)
    return (x2 == x)
end

@test test_solve()
@test test_substitutions()


text:
import Statistics
import Statistics: require_one_based_indexing
import NaNStatistics: _nanquantile!
# quantile, quantile!, _quantile,

# const q_missing = missing;
const q_missing = NaN; #-9999.0;

function quantile!(q::AbstractArray, v::AbstractVector, p::AbstractArray;
  sorted::Bool=false, alpha::Real=1.0, beta::Real=alpha)
  require_one_based_indexing(q, v, p)
  if size(p) != size(q)
    throw(DimensionMismatch("size of p, $(size(p)), must equal size of q, $(size(q))"))
  end
  isempty(q) && return q

  minp, maxp = extrema(p)
  _quantilesort!(v, sorted, minp, maxp)

  for (i, j) in zip(eachindex(p), eachindex(q))
    @inbounds q[j] = _quantile(v, p[i], alpha=alpha, beta=beta)
  end
  return q
end

function quantile!(v::AbstractVector, p::Union{AbstractArray,Tuple{Vararg{Real}}};
  sorted::Bool=false, alpha::Real=1.0, beta::Real=alpha)
  if !isempty(p)
    minp, maxp = extrema(p)
    _quantilesort!(v, sorted, minp, maxp)
  end
  return map(x -> _quantile(v, x, alpha=alpha, beta=beta), p)
end

quantile!(v::AbstractVector, p::Real; sorted::Bool=false, alpha::Real=1.0, beta::Real=alpha) =
  _quantile(_quantilesort!(v, sorted, p, p), p, alpha=alpha, beta=beta)

# Function to perform partial sort of v for quantiles in given range
function _quantilesort!(v::AbstractArray, sorted::Bool, minp::Real, maxp::Real)
  isempty(v) && return #throw(ArgumentError("empty data vector"))
  require_one_based_indexing(v)

  if !sorted
    lv = length(v)
    lo = floor(Int, minp * (lv))
    hi = ceil(Int, 1 + maxp * (lv))

    # only need to perform partial sort
    sort!(v, 1, lv, Base.Sort.PartialQuickSort(lo:hi), Base.Sort.Forward)
  end
  if (sorted && (ismissing(v[end]) || (v[end] isa Number && isnan(v[end])))) ||
     any(x -> ismissing(x) || (x isa Number && isnan(x)), v)
    throw(ArgumentError("quantiles are undefined in presence of NaNs or missing values"))
  end
  return v
end

# Core quantile lookup function: assumes `v` sorted
@inline function _quantile(v::AbstractVector, p::Real; alpha::Real=1.0, beta::Real=alpha)
  0 <= p <= 1 || throw(ArgumentError("input probability out of [0,1] range"))
  0 <= alpha <= 1 || throw(ArgumentError("alpha parameter out of [0,1] range"))
  0 <= beta <= 1 || throw(ArgumentError("beta parameter out of [0,1] range"))
  require_one_based_indexing(v)

  n = length(v)

  # @assert n > 0 # this case should never happen here
  n == 0 && return q_missing

  m = alpha + p * (one(alpha) - alpha - beta)
  aleph = n * p + oftype(p, m)
  j = clamp(trunc(Int, aleph), 1, n - 1)
  γ = clamp(aleph - j, 0, 1)

  if n == 1
    a = v[1]
    b = v[1]
  else
    a = v[j]
    b = v[j+1]
  end

  if isfinite(a) && isfinite(b)
    return a + γ * (b - a)
  else
    return (1 - γ) * a + γ * b
  end
end

_quantile(v::Nothing, p::Real; alpha::Real=1.0, beta::Real=alpha) = q_missing;
_quantile(v::Vector{Union{}}, p::Real; alpha::Real=1.0, beta::Real=alpha) = q_missing;

quantile!(v::Vector{Union{}}, p::Float64; kw...) = q_missing;

quantile2(itr, p; sorted::Bool=false, alpha::Real=1.0, beta::Real=alpha) =
  quantile!(collect(itr), p, sorted=sorted, alpha=alpha, beta=beta)

quantile2(v::AbstractVector, p; sorted::Bool=false, alpha::Real=1.0, beta::Real=alpha) =
  quantile!(sorted ? v : Base.copymutable(v), p; sorted=sorted, alpha=alpha, beta=beta)

# fix quantile missing
quantile2(x::Vector{Missing}, p; kw...) = repeat([q_missing], length(p))

function Quantile(array::AbstractArray{<:Real}, probs=[0, 0.25, 0.5, 0.75, 1]; dims=1, missval=nothing)
  if missval === nothing
    mapslices(x -> quantile2(x, probs), array, dims=dims)
  else
    array = to_missing(array, missval)
    Quantile(array, probs, dims=dims)
  end
end

function Quantile(array::AbstractNanArray, probs=[0, 0.25, 0.5, 0.75, 1]; dims=1)
  mapslices(x -> begin
      ans = quantile2(skipmissing(x), probs)
    end, array, dims=dims)
end

# `dims`: symbol
"""
  Quantile(array::AbstractArray{<:Real}, probs=[0, 0.25, 0.5, 0.75, 1]; dims=1, missval = nothing)
  Quantile(array::AbstractNanArray, probs=[0, 0.25, 0.5, 0.75, 1]; dims=1)

# Examples
```julia
arr = rand(200, 200, 365);
d = DimArray(arr, ["lon", "lat", "time"]);
probs = [0.5, 0.9];
Quantile(d, probs; dims = :time)
```
"""
function Quantile(da::AbstractDimArray, probs=[0, 0.25, 0.5, 0.75, 1]; dims=1)
  # dims = dimnum2(array, dims)
  if eltype(dims) <: Integer
    dims = name(da.dims)[dims]
  end
  if dims isa String
    dims = Symbol(dims)
  end

  bands = collect(name(da.dims))
  r = mapslices(x -> quantile(x, probs), da, dims=dims)
  r = DimArray(r.data, bands) # dimension error, rebuild it
  set(r, dims => :prob)
end

"""
  nanquantile(A::AbstractArray{T,N}, probs::Vector{<:Real};
    dims::Integer=1, type = Float64) where {T,N}

# Examples
```julia
x = rand(4, 4, 201);
probs = [0.9, 0.99, 0.9999]

r1 = Quantile(x, probs, dims=3);
r2 = nanquantile(x, probs, dims=3);
```
"""
function nanquantile(A::AbstractArray{T,N}, probs::Vector{<:Real};
  dims::Integer=1, type=Float64) where {T,N}

  Aₜ = copyto!(Array{T,N}(undef, size(A)), A)

  Size = size(A) |> collect
  Size[dims] = length(probs)

  res = zeros(type, Size...)
  for k = 1:length(probs)
    q = probs[k]
    selectdim(res, dims, k) .= _nanquantile!(Aₜ, q, dims)
  end
  res
end

export quantile2, nanquantile, Quantile


text:
function gen_β(β, η, D, R)
    α = 2.0 .- β .^ (-  η - 1.0 )
    mask = R .<= 1.0 ./ α
    s = 1.0 / (η + 1.0)
    βq = [ mask[i] ?  (R[i] * α[i])^s : (1.0 / (2.0 - R[i]*α[i]))^s for i in 1:D]
    βq
end

"""
    SBX_crossover(vector1, vector2, bounds, η=15, p_variable = 0.9)

Simulated binomial crossover for given two `Vectors{Real}`.
"""
function SBX_crossover(vector1, vector2, bounds, η=15, p_variable = 0.9)
    xu = view(bounds, 2,:)
    xl = view(bounds, 1,:)
    D = length(vector1)

    do_crossover = ones(Bool, D)
    do_crossover[rand(D) .> p_variable] .= false
    do_crossover[ abs.( vector2 - vector1 ) .<= eps() ] .= false

    y1 = min.( vector1, vector2 )
    y2 = max.( vector1, vector2 )
    Δ = max.(eps(), y2 - y1)

    R = rand(D)

    β = @. 1.0 + (2.0 * (y1 - xl) / Δ)
    βq = gen_β(β, η, D, R)
    c1 = @. 0.5*(y1 + y2 -  βq*Δ)

    β = @. 1.0 + (2.0 * (y1 - xl) / Δ)
    βq = gen_β(β, η, D, R)
    c2 = @. 0.5*(y1 + y2 +  βq*Δ)

    # swap
    mask = rand(Bool, D)
    cc = copy(c1)
    c1[mask] = c2[mask]
    c2[mask] = cc[mask]

    cc1 = copy(vector1)
    cc1[do_crossover] = c1[do_crossover]
    cc2 = copy(vector2)
    cc2[do_crossover] = c2[do_crossover]

    reset_to_violated_bounds!(cc1, bounds)
    reset_to_violated_bounds!(cc2, bounds)

    return cc1, cc2
end

"""
    DE_crossover(x, u, CR)

Binomial crossover between x and u for Differential Evolution with probability CR, i.e.,
`v[j] = u[j]` if `rand() < CR`, otherwise `v[j] = x[j]`. Return `v`.
"""
function DE_crossover(x, u, CR)
    D = length(x)
    # binomial crossover
    v = zeros(D)
    j_rand = rand(1:D)

    # binomial crossover
    for j = 1:D
        if rand() < CR || j == j_rand
            v[j] = u[j]
        else
            v[j] = x[j]
        end
    end

    return v
end


text:
function readTSP(path::AbstractString)
  raw = read(path, String)
  checkEOF(raw)
  return _generateTSP(raw)
end

readTSPLIB(instance::Symbol) = readTSP(joinpath(TSPLIB95_path,string(instance)*".tsp"))

function _generateTSP(raw::AbstractString)
  _dict = keyextract(raw, tsp_keys)
  name = _dict["NAME"]
  dimension = parse(Int,_dict["DIMENSION"])
  weight_type = _dict["EDGE_WEIGHT_TYPE"]
  dxp = false

  if weight_type == "EXPLICIT" && haskey(_dict,"EDGE_WEIGHT_SECTION")
    explicits = parse.(Float64, split(_dict["EDGE_WEIGHT_SECTION"]))
    weights = explicit_weights(_dict["EDGE_WEIGHT_FORMAT"],explicits)
    #Push display data to nodes if possible
    if haskey(_dict,"DISPLAY_DATA_SECTION")
      coords = parse.(Float64, split(_dict["DISPLAY_DATA_SECTION"]))
      n_r = convert(Integer,length(coords)/dimension)
      nodes = reshape(coords,(n_r,dimension))'[:,2:end]
      dxp = true
    else
      nodes = zeros(dimension,2)
    end
  elseif haskey(_dict,"NODE_COORD_SECTION")
    coords = parse.(Float64, split(_dict["NODE_COORD_SECTION"]))
    n_r = convert(Integer,length(coords)/dimension)
    nodes = reshape(coords,(n_r,dimension))'[:,2:end]
    weights = calc_weights(_dict["EDGE_WEIGHT_TYPE"], nodes)
  end

  fFX = fullFit(weights)
  pFX = partFit(weights)

  if haskey(Optimals, Symbol(name))
    optimal = Optimals[Symbol(name)]
  else
    optimal = -1
  end

  TSP(name,dimension,weight_type,weights,nodes,dxp,fFX,pFX,optimal)
end

function keyextract(raw::T,ks::Array{T}) where T<:AbstractString
  pq = PriorityQueue{T,Tuple{Integer,Integer}}()
  vals = Dict{T,T}()
  for k in ks
    idx = findfirst(k,raw)
    idx != nothing && enqueue!(pq,k,extrema(idx))
  end
  while length(pq) > 1
    s_key, s_pts = peek(pq)
    dequeue!(pq)
    f_key, f_pts = peek(pq)
    rng = (s_pts[2]+1):(f_pts[1]-1)
    vals[s_key] = strip(replace(raw[rng],":"=>""))
  end
  return vals
end

function explicit_weights(key::AbstractString,data::Vector{Float64})
  w = @match key begin
    "UPPER_DIAG_ROW"  => vec2UDTbyRow(data)
    "LOWER_DIAG_ROW"  => vec2LDTbyRow(data)
    "UPPER_DIAG_COL"  => vec2UDTbyCol(data)
    "LOWER_DIAG_COL"  => vec2LDTbyCol(data)
    "UPPER_ROW"       => vec2UTbyRow(data)
    "LOWER_ROW"       => vec2LTbyRow(data)
    "FULL_MATRIX"     => vec2FMbyRow(data)
  end
  if !in(key,["FULL_MATRIX"])
    w.+=w'
  end
  return w
end

function calc_weights(key::AbstractString, data::Matrix)
  w = @match key begin
    "EUC_2D"  => euclidian(data[:,1], data[:,2])
    "MAN_2D"  => manhattan(data[:,1], data[:,2])
    "MAX_2D"  => max_norm(data[:,1], data[:,2])
    "GEO"     => geo(data[:,1], data[:,2])
    "ATT"     => att_euclidian(data[:,1], data[:,2])
    "CEIL_2D" => ceil_euclidian(data[:,1], data[:,2])
    _         => error("Distance function type $key is not supported.")
  end

  return w
end

function checkEOF(raw::AbstractString)
  n = findlast("EOF",raw)
  if n == nothing
    throw("EOF not found")
  end
  return
end