text:
# Types

export
    Interaction,
    GeneralInteraction,
    SpecificInteraction,
    AbstractCutoff,
    AbstractCoupler,
    NeighborFinder,
    Logger,
    InteractionList2Atoms,
    InteractionList3Atoms,
    InteractionList4Atoms,
    Atom,
    mass,
    AtomData,
    NeighborList,
    NeighborListVec,
    System

const DefaultFloat = Float64

"An interaction between atoms that contributes to forces on the atoms."
abstract type Interaction end

"""
A general interaction that will apply to all or most atom pairs.
Custom general interactions should sub-type this type.
"""
abstract type GeneralInteraction <: Interaction end

"""
A specific interaction between sets of specific atoms, e.g. a bond angle.
Custom specific interactions should sub-type this type.
"""
abstract type SpecificInteraction <: Interaction end

"""
A general type of cutoff encoding the approximation used for a potential.
Interactions can be parameterized by the cutoff behavior.
"""
abstract type AbstractCutoff end

"""
A way to keep properties of a system constant.
Custom temperature and pressure couplers should sub-type this type.
"""
abstract type AbstractCoupler end

"""
A way to find near atoms to save on simulation time.
Custom neighbor finders should sub-type this type.
"""
abstract type NeighborFinder end

"""
A way to record a property, e.g. the temperature, throughout a simulation.
Custom loggers should sub-type this type.
"""
abstract type Logger end

"""
    InteractionList2Atoms(is, js, inters)
    InteractionList2Atoms(inter_type)

A list of specific interactions between two atoms.
"""
struct InteractionList2Atoms{T}
    is::Vector{Int}
    js::Vector{Int}
    inters::T
end

"""
    InteractionList3Atoms(is, js, ks, inters)
    InteractionList3Atoms(inter_type)

A list of specific interactions between three atoms.
"""
struct InteractionList3Atoms{T}
    is::Vector{Int}
    js::Vector{Int}
    ks::Vector{Int}
    inters::T
end

"""
    InteractionList4Atoms(is, js, ks, ls, inters)
    InteractionList4Atoms(inter_type)

A list of specific interactions between four atoms.
"""
struct InteractionList4Atoms{T}
    is::Vector{Int}
    js::Vector{Int}
    ks::Vector{Int}
    ls::Vector{Int}
    inters::T
end

InteractionList2Atoms(T) = InteractionList2Atoms{Vector{T}}([], [], T[])
InteractionList3Atoms(T) = InteractionList3Atoms{Vector{T}}([], [], [], T[])
InteractionList4Atoms(T) = InteractionList4Atoms{Vector{T}}([], [], [], [], T[])

"""
    Atom(; <keyword arguments>)

An atom and its associated information.
Properties unused in the simulation or in analysis can be left with their
default values.
The types used should be bits types if the GPU is going to be used.

# Arguments
- `index::Int`: the index of the atom in the system.
- `charge::C=0.0`: the charge of the atom, used for electrostatic interactions.
- `mass::M=0.0u"u"`: the mass of the atom.
- `σ::S=0.0u"nm"`: the Lennard-Jones finite distance at which the inter-particle
    potential is zero.
- `ϵ::E=0.0u"kJ * mol^-1"`: the Lennard-Jones depth of the potential well.
"""
struct Atom{C, M, S, E}
    index::Int
    charge::C
    mass::M
    σ::S
    ϵ::E
end

function Atom(;
                index=1,
                charge=0.0,
                mass=0.0u"u",
                σ=0.0u"nm",
                ϵ=0.0u"kJ * mol^-1")
    return Atom(index, charge, mass, σ, ϵ)
end

"""
    mass(atom)

The mass of an atom.
"""
mass(atom::Atom) = atom.mass

AtomsBase.atomic_mass(atom::Atom) = mass(atom)

function Base.show(io::IO, a::Atom)
    print(io, "Atom with index ", a.index, ", charge=", a.charge,
            ", mass=", a.mass, ", σ=", a.σ, ", ϵ=", a.ϵ)
end

"""
    AtomData(atom_type, atom_name, res_number, res_name)

Data associated with an atom.
Storing this separately allows the atom types to be bits types and hence
work on the GPU.
"""
struct AtomData
    atom_type::String
    atom_name::String
    res_number::Int
    res_name::String
    element::String
end

function AtomData(;
                    atom_type="?",
                    atom_name="?",
                    res_number=1,
                    res_name="???",
                    element="?")
    return AtomData(atom_type, atom_name, res_number, res_name, element)
end

"""
    NeighborList()
    NeighborList(n, list)

Structure to contain pre-allocated neighbor lists.
"""
mutable struct NeighborList
    n::Int # Number of neighbors in list (n <= length(list))
    list::Vector{Tuple{Int, Int, Bool}}
end

NeighborList() = NeighborList(0, [])

function Base.empty!(nl::NeighborList)
    nl.n = 0
    return nl
end

function Base.push!(nl::NeighborList, element::Tuple{Int, Int, Bool})
    nl.n += 1
    if nl.n > length(nl.list)
        push!(nl.list, element)
    else
        nl.list[nl.n] = element
    end
    return nl
end

function Base.append!(nl::NeighborList, list::AbstractVector{Tuple{Int, Int, Bool}})
    for element in list
        push!(nl, element)
    end
    return nl
end

Base.append!(nl::NeighborList, nl_app::NeighborList) = append!(nl, @view(nl_app.list[1:nl_app.n]))

"""
    NeighborListVec(n, list)

Structure to contain neighbor lists for broadcasting.
"""
struct NeighborListVec{T}
    nbsi::Vector{Int} # Sorted ascending
    nbsj::Vector{Int}
    atom_bounds_i::Vector{Int}
    atom_bounds_j::Vector{Int}
    sortperm_j::Vector{Int}
    weights_14::T
end

"""
    System(; <keyword arguments>)

A physical system to be simulated.
Properties unused in the simulation or in analysis can be left with their
default values.
`atoms`, `atoms_data`, `coords` and `velocities` should have the same length.
This is a sub-type of `AbstractSystem` from AtomsBase.jl and implements the
interface described there.

# Arguments
- `atoms::A`: the atoms, or atom equivalents, in the system. Can be
    of any type but should be a bits type if the GPU is used.
- `atoms_data::AD`: other data associated with the atoms, allowing the atoms to
    be bits types and hence work on the GPU.
- `general_inters::GI=()`: the general interactions in the system, i.e.
    interactions between all or most atoms such as electrostatics. Typically a
    `Tuple`.
- `specific_inter_lists::SI=()`: the specific interactions in the system,
    i.e. interactions between specific atoms such as bonds or angles. Typically
    a `Tuple`.
- `coords::C`: the coordinates of the atoms in the system. Typically a
    vector of `SVector`s of 2 or 3 dimensions.
- `velocities::V=zero(coords)`: the velocities of the atoms in the system.
- `box_size::B`: the size of the box in which the simulation takes place.
    Typically a `SVector` of 2 or 3 dimensions.
- `neighbor_finder::NF=NoNeighborFinder()`: the neighbor finder used to find
    close atoms and save on computation.
- `loggers::L=Dict{String, Logger}()`: the loggers that record properties of
    interest during a simulation.
- `force_unit::F=u"kJ * mol^-1 * nm^-1"`: the unit of force of the system.
- `energy_unit::E=u"kJ * mol^-1"`: the unit of energy of the system.
- `gpu_diff_safe::Bool`: whether to use the GPU implementation. Defaults to
    `isa(coords, CuArray)`.
"""
mutable struct System{D, S, G, A, AD, GI, SI, C, V, B, NF, L, F, E} <: AbstractSystem{D, S}
    atoms::A
    atoms_data::AD
    general_inters::GI
    specific_inter_lists::SI
    coords::C
    velocities::V
    box_size::B
    neighbor_finder::NF
    loggers::L
    force_unit::F
    energy_unit::E
end

function System(;
                atoms,
                atoms_data=[],
                general_inters=(),
                specific_inter_lists=(),
                coords,
                velocities=zero(coords),
                box_size,
                neighbor_finder=NoNeighborFinder(),
                loggers=Dict{String, Logger}(),
                force_unit=u"kJ * mol^-1 * nm^-1",
                energy_unit=u"kJ * mol^-1",
                gpu_diff_safe=isa(coords, CuArray))
    D = length(box_size)
    S = eltype(atoms)
    A = typeof(atoms)
    AD = typeof(atoms_data)
    GI = typeof(general_inters)
    SI = typeof(specific_inter_lists)
    C = typeof(coords)
    V = typeof(velocities)
    B = typeof(box_size)
    NF = typeof(neighbor_finder)
    L = typeof(loggers)
    F = typeof(force_unit)
    E = typeof(energy_unit)
    return System{D, S, gpu_diff_safe, A, AD, GI, SI, C, V, B, NF, L, F, E}(
                    atoms, atoms_data, general_inters, specific_inter_lists,
                    coords, velocities, box_size, neighbor_finder, loggers,
                    force_unit, energy_unit)
end

Base.getindex(s::System, i::Integer) = s.atoms[i]
Base.length(s::System) = length(s.atoms)

AtomsBase.species(s::System) = s.atoms
AtomsBase.species(s::System, i::Integer) = s.atoms[i]

AtomsBase.position(s::System) = s.coords
AtomsBase.position(s::System, i::Integer) = s.coords[i]

AtomsBase.velocity(s::System) = s.velocities
AtomsBase.velocity(s::System, i::Integer) = s.velocities[i]

AtomsBase.boundary_conditions(::System{3}) = SVector(Periodic(), Periodic(), Periodic())
AtomsBase.boundary_conditions(::System{2}) = SVector(Periodic(), Periodic())

edges_to_box(bs::SVector{3}, z) = SVector{3}([SVector(bs[1], z, z), SVector(z, bs[2], z),
                                                SVector(z, z, bs[3])])
edges_to_box(bs::SVector{2}, z) = SVector{2}([SVector(bs[1], z), SVector(z, bs[2])])

function AtomsBase.bounding_box(s::System)
    bs = s.box_size
    z = zero(bs[1])
    bb = edges_to_box(bs, z)
    return unit(z) == NoUnits ? (bb)u"nm" : bb # Assume nm without other information
end

function Base.show(io::IO, ::MIME"text/plain", s::System)
    print(io, "System with ", length(s), " atoms, box size ", s.box_size)
end


text:
#=
# Parse json to dict
#TODO : fix broken data formats
function parse_json(filename,device_names)
    Devices =Dict{String,Any}()
    for name in device_names
        if isfile("$filename/x$name.json")
            temp = Dict()
            open("$filename/x$name.json", "r") do f
            global temp
            dicttxt = readstring(f)  # file information to string
            temp=JSON2.read(dicttxt,Dict{Any,Array{Dict}})  # parse and transform data
            Devices[name] = temp
            end
        end
    end
    return Devices
end
=#

ext: jl


text:
module LightGraphConvert

using CircuitViz
using LightGraphs

g = SimpleDiGraph()
LightGraphs.add_vertex!(g)
LightGraphs.add_vertex!(g)
LightGraphs.add_vertex!(g)
LightGraphs.add_edge!(g, 1, 2)
LightGraphs.add_edge!(g, 1, 3)

new = CircuitViz.convert_to_catlab_graph(g)
compile!(new)

end # module


text:
module TaylorModels

using Reexport
@reexport using TaylorSeries, IntervalArithmetic
# using IntervalRootFinding
using RecipesBase

using Markdown

import Base: setindex!, getindex, copy,
    ==, +, -, *, /, ^,
    zero, one, findfirst, #iszero,
    promote, show,
    inv, sqrt, exp, log, sin, cos, tan,
    asin, acos, atan, sinh, cosh, tanh

import TaylorSeries: integrate, get_order, evaluate, pretty_print

import IntervalArithmetic: showfull

# export Taylor1Model, bound, make_Taylor_model, TMcomposition,
#         taylor1_var, integrate, degree,
#         calculate_set, Taylor_step

export TaylorModel1, RTaylorModel1, TaylorModelN

export remainder, polynomial,
    rpa, fp_rpa, bound_absrem, bound_relrem

include("constructors.jl")
include("auxiliary.jl")
include("promotion.jl")
include("bounds.jl")
include("evaluate.jl")
include("rpa_functions.jl")
include("arithmetic.jl")
include("integration.jl")
include("recipe.jl")
include("show.jl")

# include("Taylor1/Taylor1.jl")
# include("TaylorN/TaylorN.jl")

end # module


text:
using NaiveNASflux, Flux, Test
using NaiveNASlib.Advanced, NaiveNASlib.Extend

function assertlayer(l, Wexp, bexp)
    @test size(Wexp) == size(weights(l))
    if bexp isa Flux.Zeros
        @test bias(l) isa Flux.Zeros
    else
        @test size(bexp) == size(bias(l))
    end
    @test Wexp == weights(l)
    @test bexp == bias(l)
end

@testset "NaiveNASflux.jl" begin

    @info "Testing util"
    include("util.jl")

    @info "Testing select"
    include("select.jl")

    @info "Testing mutable"
    include("mutable.jl")

    @info "Testing vertex"
    include("vertex.jl")

    @info "Testing neuronutility"
    include("neuronutility.jl")

    @info "Testing examples"
    include("examples.jl")

    if Int !== Int32
        import Documenter
        Documenter.doctest(NaiveNASflux)
    end
end


text:
using GLM, DataFrames
include("misc.jl")

abstract type FitGoodness end
abstract type InfoMeasure <: FitGoodness end

abstract type AIC <: InfoMeasure end
abstract type AICc <: InfoMeasure end
abstract type BIC <: InfoMeasure end

function AIC(dfrm::RegressionModel)
    n = size(dfrm.mf.df)[1] # number of cases
    if dfrm.mf.terms.intercept
        k = length(dfrm.mf.terms.terms) + 1
    else
        k = length(dfrm.mf.terms.terms)
    end
    # number of parameters, including intercept
    return n*log(deviance(dfrm.model.rr)/n) + 2*k
end

function AICc(dfrm::RegressionModel)
    n = size(dfrm.mf.df)[1]
    if dfrm.mf.terms.intercept
        k = length(dfrm.mf.terms.terms) + 1
    else
        k = length(dfrm.mf.terms.terms)
    end
    return AIC(dfrm) + (2*k*(k+1)) / (n-k-1)
end

function BIC(dfrm::RegressionModel)
    n = size(dfrm.mf.df)[1]
    if dfrm.mf.terms.intercept
        k = length(dfrm.mf.terms.terms) + 1
    else
        k = length(dfrm.mf.terms.terms)
    end
    return n*log(deviance(dfrm.model.rr)/n) + log(n)*k
end

struct extractrhs
    rhsarray::Array
    rhsstring::AbstractString
end

# Extract the right hand side of a RegressionModel formula
function extractrhs(dfrm::RegressionModel)
    regressors = dfrm.mf.terms.terms
    rhsarray = similar(regressors, AbstractString)
    for i in 1:length(regressors)
        rhsarray[i] = string(regressors[i])
        rhsarray[i] = replace(rhsarray[i], " ", "")
    end
    rhsstring = ""
    for var in rhsarray
        rhsstring = rhsstring*"+"*var
    end
    if !dfrm.mf.terms.intercept
        rhsarray = [0, rhsarray]
        rhsstring = "0+"*rhsstring
    end
    return extractrhs(rhsarray, rhsstring)
end

abstract type VarSel end

struct add1 <: VarSel
    aic::Float64
    add::AbstractString
    model::RegressionModel
end

struct drop1 <: VarSel
    aic::Float64
    drop::AbstractString
    model::RegressionModel
end

function add1(dfrm::RegressionModel, scope::AbstractString, data::DataFrame)
    # e.g. scope = "X1+X2+X1&X2"
    replace(scope, " ", "") # remove whitespace in scope
    scope = split(scope, "+")
    index = find(notnull, scope)
    scope = scope[index]

    lhs = dfrm.mf.terms.eterms[1] # response
    rhs = extractrhs(dfrm)

    add = ""
    aic = AIC(dfrm)
    model = dfrm
    for var in scope
        if var in rhs.rhsarray
            continue # if var already in model, try next var
        end
        rhsnew = rhs.rhsstring * "+" * var
        rhsnew = rhsnew[2:end]
        rhsnew = parse(rhsnew)
        fnew = Formula(lhs, rhsnew)
        newfit = fit(LinearModel, fnew, data)
        newaic = AIC(newfit)
        if newaic < aic
            aic = newaic
            add = var
            f = fnew
            model = newfit
        end
    end
    if add == ""
        println("No term added")
        return add1(aic, "N/A", dfrm)
    else
        println("Add $add with AIC = $aic")
        return add1(aic, add, model)
    end
end

function drop1(dfrm::RegressionModel, scope::AbstractString)
    # e.g. scope = "X1+X2+X1&X2"
    replace(scope, " ", "") # remove whitespace in scope
    scope = split(scope, "+")
    index = find(notnull, scope)
    scope = scope[index]

    lhs = dfrm.mf.terms.eterms[1] # response
    rhs = extractrhs(dfrm) # extract rhs

    drop = ""
    aic = AIC(dfrm)
    model = dfrm
    for var in scope
        if var in rhs.rhsarray
            var = "+"*var
            rhsnew = replace(rhs.rhsstring, var, "")
        else
            error("$var is not in the original model, please modify scope")
        end
        rhsnew = rhsnew[2:end]
        rhsnew = parse(rhsnew)
        fnew = Formula(lhs, rhsnew)
        newfit = fit(LinearModel, fnew, dfrm.mf.df)
        newaic = AIC(newfit)
        if newaic < aic
            aic = newaic
            drop = var[2:end]
            f = fnew
            model = newfit
        end
    end
    if drop == ""
        println("No term dropped")
        return drop1(aic, "N/A", dfrm)
    else
        println("Drop $drop with AIC = $aic")
    return drop1(aic, drop, model)
    end
end

drop1(dfrm::RegressionModel) = drop1(dfrm, extractrhs(dfrm).rhsstring)

struct step <: VarSel
    aic::Float64
    model::RegressionModel
end

function step(dfrm::RegressionModel, scope::AbstractString, data::DataFrame,
              direction::AbstractString, trace::Bool=false)
    directions = ["backward", "forward", "both"]
    if ~ (direction in directions)
        error("direction must be one of 'backward', 'forward', or 'both'")
    end
    # turn scope to a string array
    scope = split(scope, "+")
    index = find(notnull, scope)
    scope = scope[index]
    # initialize variables
    todrop = ""
    aic = AIC(dfrm)
    model = dfrm
    if direction == "backward"
        drop = drop1(dfrm, scope)
        while drop.drop != ""
            scope = replace(scope, drop.drop, "")
            drop = drop1(drop.model, scope)
        end
        if drop == ""
            println("No term dropped")
            return step(drop.aic, dfrm)
        else
            return step(drop.aic, model)
        end

    elseif direction == "forward"
        add = add1(dfrm, scope, data)
        while add.add != ""
            scope = replace(scope, add.add, "")
            add = add1(add.model, scope, data)
        end
        if add == ""
            println("No term dropped")
            return step(add.aic, dfrm)
        else
            return step(add.aic, model)
        end

    else # both
        drop = drop1(dfrm, scope)
        add = add1(dfrm, scope)
        while drop.drop != "" | add.add != ""
            if drop.aic < add.aic
                model = drop.model
                aic = drop.aic
                scopenew = replace(scope, drop.drop, "")
            else
                model = add.model
                aic = add.aic
                scopenew = replace(scope, add.add, "")
            end
            drop = drop1(model, scopenew)
            add = add1(model, scopenew)
        end
        if scopenew == scope
            println("No term changed")
            return step(AIC(dfrm), dfrm)
        else
            return step(aic, model)
        end
    end
end

function step(dfrm::RegressionModel, direction::AbstractString; trace::Bool=false)
    directions = ["backward", "both"]
    if ~ (direction in directions)
        error("If no scope specified, direction must be 'backward' or 'both'")
    end
    todrop = ""
    aic = AIC(dfrm)
    model = dfrm
    if direction == "backward"
        drop = drop1(dfrm)
        while drop.aic < aic
            aic = drop.aic
            todrop = drop.drop
            model = drop.model
            drop = drop1(model)
        end
        if todrop == ""
            println("No term dropped")
            return step(aic, dfrm)
        else
            return step(aic, model)
        end

    else # both
        scope = extractrhs(dfrm).rhsstring
        data = dfrm.mf.df
        toadd = ""
        drop = drop1(dfrm)
        add = add1(dfrm, scope, data)
        while drop.aic < aic || add.aic < aic
            if drop.aic < add.aic
                aic = drop.aic
                todrop = drop.drop
                model = drop.model
            else
                aic = add.aic
                toadd = add.add
                model = add.model
            end
            drop = drop1(model)
            add = add1(model, scope, data)
        end
        if todrop == "" && toadd == ""
            println("No term changed")
            return step(AIC(dfrm), dfrm)
        else
            return step(aic, model)
        end
    end
end

max_stars_repo_name: wakakusa/RegTools

max_issues_repo_name: wakakusa/RegTools

max_forks_repo_name: wakakusa/RegTools


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Sequential(transforms)

A transform where `transforms` are applied in sequence.

# Examples

```julia
Scale(low=0.2, high=0.8) → EigenAnalysis(:VDV)
Select(:x, :z) → ZScore() → EigenAnalysis(:V) → Scale(low=0, high=1)
```
"""
struct Sequential <: Transform
  transforms::Vector{Transform}
end

isrevertible(s::Sequential) = all(isrevertible, s.transforms)

function apply(s::Sequential, table)
  allcache = []
  current  = table
  for transform in s.transforms
    current, cache = apply(transform, current)
    push!(allcache, cache)
  end
  current, allcache
end

function revert(s::Sequential, newtable, cache)
  allcache = deepcopy(cache)
  current  = newtable
  for transform in reverse(s.transforms)
    current = revert(transform, current, pop!(allcache))
  end
  current
end

function reapply(s::Sequential, table, cache)
  # basic checks
  ntrans = length(s.transforms)
  ncache = length(cache)
  @assert ntrans == ncache "invalid cache for transform"

  current = table
  for (ctransform, ccache) in zip(s.transforms, cache)
    current = reapply(ctransform, current, ccache)
  end
  current
end

"""
    transform₁ → transform₂ → ⋯ → transformₙ

Create a [`Sequential`](@ref) transform with
`[transform₁, transform₂, …, transformₙ]`.
"""
→(t1::Transform, t2::Transform)   = Sequential([t1, t2])
→(t1::Transform, t2::Sequential)  = Sequential([t1; t2.transforms])
→(t1::Sequential, t2::Transform)  = Sequential([t1.transforms; t2])
→(t1::Sequential, t2::Sequential) = Sequential([t1.transforms; t2.transforms])


text:
# Plot benchmark results

using CSV
using DataFrames
using Gadfly
using Cairo
using Fontconfig

df = CSV.read("benchmarks.csv")

benchmarks = [
    "Parse PDB 1CRN",
    "Parse PDB 1HTQ",
    "Parse mmCIF 1CRN",
    "Parse mmCIF 1HTQ",
    "Parse MMTF 1CRN",
    "Parse MMTF 1HTQ",
    "Count",
    "Distance",
    "Ramachandran",
]
benchind(b) = findfirst(x -> x == b, benchmarks)
df_sorted = sort(df, (order(:Benchmark, by=benchind), :Runtime))

theme = Theme(
    background_color="white",
    panel_stroke="black",
    major_label_color="black",
    minor_label_color="black",
    highlight_width=0mm,
)

p = Gadfly.with_theme(theme) do
    plot(df_sorted,
        x=:Package,
        y=:Runtime,
        color=:Benchmark,
        Scale.y_log10,
        Guide.xlabel(nothing),
        Guide.ylabel("Runtime / s"),
        Geom.point,
        Geom.line,
    )
end

draw(PNG("plot/plot.png", dpi=300), p)

max_stars_repo_name: wojdyr/pdb-benchmarks

max_issues_repo_name: wojdyr/pdb-benchmarks

max_forks_repo_name: wojdyr/pdb-benchmarks


text:
import ..Cos: CosXString

export pdFontIsBold,
    pdFontIsItalic,
    pdFontIsFixedW,
    pdFontIsAllCap,
    pdFontIsSmallCap

using Rectangle

#=
Sample CMaps are available now as 8.cmap and 16.cmap in the test/files directory
for 8 and 16-bit toUnicode CMaps.
CMaps can have both 8 and 16 bit ranges in the same CMap file as well.
=#

const beginbfchar = b"beginbfchar"
const endbfchar   = b"endbfchar"
const beginbfrange = b"beginbfrange"
const endbfrange = b"endbfrange"
const begincodespacerange = b"begincodespacerange"
const endcodespacerange = b"endcodespacerange"

mutable struct CMap
    code_space::IntervalTree{UInt8,
                             Union{CosNullType, IntervalTree{UInt8, CosNullType}}}
    range_map::IntervalTree{UInt8,
                            Union{CosObject, IntervalTree{UInt8, CosObject}}}
    function CMap()
        cs = IntervalTree{UInt8,
                          Union{CosNullType, IntervalTree{UInt8, CosNullType}}}()
        rm = IntervalTree{UInt8,
                          Union{CosObject, IntervalTree{UInt8, CosObject}}}()
        new(cs, rm)
    end
end

function show(io::IO, cmap::CMap)
    show(io, "Code Space:\n")
    show(io, cmap.code_space)
    show(io, "Range Map:\n")
    show(io, cmap.range_map)
end

mutable struct FontUnicodeMapping
    encoding::Dict
    cmap::CMap
    hasCMap::Bool
    FontUnicodeMapping() = new(Dict(), CMap(), false)
end

function merge_encoding!(fum::FontUnicodeMapping, encoding::CosName,
                         doc::CosDoc, font::CosObject)
    encoding_mapping =
        encoding == cn"WinAnsiEncoding"   ? WINEncoding_to_Unicode :
        encoding == cn"MacRomanEncoding"  ? MACEncoding_to_Unicode :
        encoding == cn"MacExpertEncoding" ? MEXEncoding_to_Unicode :
                                            STDEncoding_to_Unicode
    merge!(fum.encoding, encoding_mapping)
    return fum
end

# for type 0 use cmap.
# for symbol and zapfdingbats - use font encoding
# for others use STD Encoding
# Reading encoding from the font files in case of Symbolic fonts are not
# supported.
# Font subset is addressed with font name identification.
function merge_encoding!(fum::FontUnicodeMapping, encoding::CosNullType,
                        doc::CosDoc, font::CosObject)
    subtype  = cosDocGetObject(doc, font, cn"Subtype")
    subtype != cn"Type1" && subtype != cn"MMType1" && return fum
    basefont = cosDocGetObject(doc, font, cn"BaseFont")
    basefont_with_subset = CDTextString(basefont)
    basefont_str = rsplit(basefont_with_subset, '+';limit=2)[end]
    enc = basefont_str == "Symbol"       ? SYMEncoding_to_Unicode :
          basefont_str == "ZapfDingbats" ? ZAPEncoding_to_Unicode :
                                           STDEncoding_to_Unicode
    merge!(fum.encoding, enc)
    return fum
end

function merge_encoding!(fum::FontUnicodeMapping,
                         encoding::Union{CosDict, CosIndirectObject{CosDict}},
                         doc::CosDoc, font::CosObject)
    baseenc = cosDocGetObject(doc, encoding, cn"BaseEncoding")
    merge_encoding!(fum, baseenc, doc, font)
    # Add the Differences
    diff = cosDocGetObject(doc, encoding, cn"Differences")
    diff === CosNull && return fum
    values = get(diff)
    d = Dict()
    cid = -1
    for v in values
        if v isa CosInt
            cid = get(v)
        else
            @assert cid != -1
            d[cid] = v
            cid += 1
        end
    end
    dict_to_unicode = dict_remap(d, AGL_Glyph_to_Unicode)
    merge!(fum.encoding, dict_to_unicode)
    return fum
end

function merge_encoding!(fum::FontUnicodeMapping, doc::CosDoc, font::CosObject)
    encoding = cosDocGetObject(doc, font, cn"Encoding")
    merge_encoding!(fum, encoding, doc, font)
    toUnicode = cosDocGetObject(doc, font, cn"ToUnicode")
    toUnicode == CosNull && return fum
    merge_encoding!(fum, toUnicode, doc, font)
end

function merge_encoding!(fum::FontUnicodeMapping,
                         cmap::CosIndirectObject{CosStream},
                         doc::CosDoc, font::CosObject)
    stm_cmap = get(cmap)
    try
        fum.cmap = read_cmap(stm_cmap)
        fum.hasCMap = true
    finally
        util_close(stm_cmap)
    end
    return fum
end

function get_glyph_id_mapping(cosdoc::CosDoc, cosfont::CosObject)
    glyph_name_id = Dict{CosName, UInt8}()
    (cosfont === CosNull) && return glyph_name_id
    subtype = get(cosfont, cn"Subtype")
    (subtype === cn"Type0") && return glyph_name_id
    baseenc = cosDocGetObject(cosdoc, cosfont, cn"BaseEncoding")
    encoding_mapping =
        baseenc == cn"WinAnsiEncoding"   ? GlyphName_to_WINEncoding :
        baseenc == cn"MacRomanEncoding"  ? GlyphName_to_MACEncoding :
        baseenc == cn"MacExpertEncoding" ? Glyphname_to_MEXEncoding :
                                           GlyphName_to_STDEncoding
    merge!(glyph_name_id, encoding_mapping)

    diff = cosDocGetObject(cosdoc, cosfont, cn"Differences")
    diff === CosNull && return glyph_name_id
    values = get(diff)
    d = Dict()
    cid = -1
    for v in values
        if v isa CosInt
            cid = get(v)
        else
            @assert cid != -1
            glyph_name_id[v] = cid
            cid += 1
        end
    end
    return glyph_name_id
end

get_encoded_string(s::CosString, fum::Nothing) = CDTextString(s)

get_encoded_string(s::CosString, fum::FontUnicodeMapping) =
    get_encoded_string(Vector{UInt8}(s), fum)

@inline function get_encoded_string(v::Vector{UInt8}, fum::FontUnicodeMapping)
    length(v) == 0 && return ""
    fum.hasCMap && return get_encoded_string(v, fum.cmap)
    return String(NativeEncodingToUnicode(v, fum.encoding))
end

function get_unicode_chars(b::UInt8, i::Interval, v::CosObject)
    f = i.lo
    l = i.hi
    if v isa CosXString
        bytes = Vector{UInt8}(v)
        carr = get_unicode_chars(bytes)
        carr[1] += (b - f)  # Only one char should be generated here
    elseif v isa CosArray
        @assert v isa CosArray
        arr = get(v)
        xstr = arr[b - f + 1]
        @assert xstr isa CosXString
        bytes = Vector{UInt8}(xstr)
        carr = get_unicode_chars(bytes)
    else
        @assert 1 == 0
    end
    return carr
end

function get_unicode_chars(barr::Vector{UInt8})
    l = length(barr)
    nb = 0
    retarr = Vector{Char}()
    while nb < l
        b1 = barr[1]
        b2 = barr[2]
        nb += 2
        c::UInt32 = 0
        if 0xD8  <= b1 <= 0xDB
            # UTF-16 Supplementary plane = 4 bytes
            b1 -= 0xD8
            c = b1
            c = (c << 8) + b2
            b3 = barr[3]
            b4 = barr[4]
            nb += 2
            if 0xDC <= b3 <= 0xDF
                b3 -= 0xDC
                c1 = b3
                c1 = (c1 << 8) + b4
                c = (c << 10) + c1
                c += 0x10000
            end
        else
            c = b1
            c = (c << 8) + b2
        end
        push!(retarr, Char(c))
    end
    return retarr
end

get_encoded_string(s::CosString, cmap::CMap) =
    get_encoded_string(Vector{UInt8}(s), cmap::CMap)

function get_encoded_string(barr::Vector{UInt8}, cmap::CMap)
    cs = cmap.code_space
    rm = cmap.range_map
    l = length(barr)
    b1 = b2 = 0x0
    carr = Vector{Char}()
    retarr = Vector{Char}()
    i = 0
    while i < l
        b1 = barr[i+=1]
        xs = intersect(cs, Interval(b1, b1))
        # When byte range is not in code space we should not return NUL.
        if length(xs) == 0
            push!(carr, Char(0))
            continue
        end
        # Some cmaps do not call out single byte ranges explicitly in the
        # code space. So may need to decipher the existence of a single
        # byte vs 2-byte code from the range map. See `else` below.
        itree = xs[1][2]
        # This case is very clearly a single byte range
        if itree === CosNull
            itv = intersect(rm, Interval(b1, b1))
            if length(itv) > 0
                carr = get_unicode_chars(b1, itv[1][1], itv[1][2])
            else
                push!(carr, Char(0))
            end
        else
            itree1 = intersect(rm, Interval(b1, b1))
            if length(itree1) == 0
                push!(carr, Char(0))
                continue
            end
            # This is a single byte range case
            if itree1[1][2] isa CosObject
                carr = get_unicode_chars(b1, itree1[1][1], itree1[1][2])
            else
                b2 = barr[i+=1]
                itv = intersect(itree1[1][2], Interval(b2, b2))
                if length(itv) > 0
                    carr = get_unicode_chars(b2, itv[1][1], itv[1][2])
                else
                    push!(carr, Char(0))
                end
            end
        end
        append!(retarr, carr)
    end
    return retarr
end

cmap_command(b::Vector{UInt8}) =
    b != beginbfchar && b != beginbfrange && b != begincodespacerange ?
        nothing : Symbol(String(b))

function on_cmap_command!(stm::IO, command::Symbol,
                         params::Vector{CosInt}, cmap::CMap)
    n = get(pop!(params))
    o1, o2, o3 = CosNull, CosNull, CosNull
    for i = 1:n
        o1 = parse_value(stm)
        @assert isa(o1, CosXString)
        d1 = Vector{UInt8}(o1)
        o2 = (command == :beginbfchar) ? o1 : parse_value(stm)
        @assert isa(o2, CosXString)
        d2 = Vector{UInt8}(o2)
        if (command != :begincodespacerange)
            o3 = parse_value(stm)
            @assert isa(o3, CosXString) || isa(o3, CosArray)
            l = length(d1)
            if l == 1
                cmap.range_map[Interval(d1[1], d2[1])] = o3
            else
                imap = get!(cmap.range_map, Interval(d1[1], d2[1]),
                            IntervalTree{UInt8, CosObject}())
                imap[Interval(d1[2], d2[2])] = o3
            end
        else
            l = length(d1)
            if l == 1
                cmap.code_space[Interval(d1[1], d2[1])] = CosNull
            else
                imap = IntervalTree{UInt8, CosNullType}()
                imap[Interval(d1[2], d2[2])] = CosNull
                cmap.code_space[Interval(d1[1], d2[1])] = imap
            end
        end
    end
    return cmap
end

on_cmap_command!(stm::IO, command::CosObject,
                 params::Vector{CosInt}, cmap::CMap) = nothing

function read_cmap(stm::IO)
    tcmap = CMap()
    params = Vector{CosInt}()
    while !eof(stm)
        obj = parse_value(stm, cmap_command)
        if isa(obj, CosInt)
            push!(params, obj)
        end
        (obj == :beginbfchar || obj == :beginbfrange ||
         obj == :begincodespacerange) &&
            on_cmap_command!(stm, obj, params, tcmap)
    end
    return tcmap
end

struct CIDWidth
    imap::IntervalTree{UInt16, Int}
    dw::Int
    CIDWidth(m::IntervalTree{UInt16, Int}, tdw::Int) = new(m, tdw)
end

CIDWidth(m::IntervalTree{UInt16, Int}) = CIDWidth(m, 1000)
CIDWidth(tdw::Int) = CIDWidth(IntervalTree{UInt16, Int}(), tdw)
CIDWidth() = CIDWidth(1000)

mutable struct PDFont
    doc::PDDoc
    obj::CosObject
    widths::Union{AdobeFontMetrics, Vector{Int}, CIDWidth}
    fum::FontUnicodeMapping
    glyph_name_id::Dict{CosName, UInt8}
    flags::UInt32
    fontname::CosName
    @inline function PDFont(doc::PDDoc, cosfont::CosObject)
        fum = FontUnicodeMapping()
        merge_encoding!(fum, doc.cosDoc, cosfont)
        widths = get_font_widths(doc.cosDoc, cosfont)
        glyph_name_id = get_glyph_id_mapping(doc.cosDoc, cosfont)
        flags = get_font_flags(doc, cosfont, widths)
        fontname = get_font_name(doc, cosfont, widths)
        return new(doc, cosfont, widths, fum, glyph_name_id, flags, fontname)
    end
end

INIT_CODE(::CIDWidth) = 0x0000
SPACE_CODE(w::CIDWidth) = get_character_code(cn"space", w)
INIT_CODE(x) = 0x00
SPACE_CODE(x) = get_character_code(cn"space", x)

pdFontIsBold(pdfont::PDFont)     = (pdfont.flags & 0x80000000) > 0
pdFontIsItalic(pdfont::PDFont)   = (pdfont.flags & 0x00000040) > 0
pdFontIsFixedW(pdfont::PDFont)   = (pdfont.flags & 0x00000001) > 0
pdFontIsAllCap(pdfont::PDFont)   = (pdfont.flags & 0x00010000) > 0
pdFontIsSmallCap(pdfont::PDFont) = (pdfont.flags & 0x00020000) > 0

# Not supported FD attribute in CIDFonts
@inline function get_font_flags(doc::PDDoc, cosfont::CosObject, widths)
    flags = 0x00000000
    refdesc = get(cosfont, cn"FontDescriptor")
    refdesc === CosNull && return get_font_flags(widths)
    cosflags = cosDocGetObject(doc.cosDoc, refdesc, cn"Flags")
    cfweight = cosDocGetObject(doc.cosDoc, refdesc, cn"FontWeight")
    cfname   = cosDocGetObject(doc.cosDoc, refdesc, cn"FontName")
    cfweight !== CosNull && get(cfweight) >= 700 && (flags |= 0x80000000)
    cfname   !== CosNull &&
        (occursin("Bold", string(cfname)) ||
         occursin("bold", string(cfname))) &&
        (flags |= 0x80000000)
    cosflags !== CosNull && (flags += UInt32(get(cosflags)))
    return flags
end
get_font_flags(x) = 0x00000000

@inline function get_font_name(doc::PDDoc, cosfont::CosObject, widths)
    refdesc = get(cosfont, cn"FontDescriptor")
    refdesc === CosNull && return get_font_name(cosfont, widths)
    return cosDocGetObject(doc.cosDoc, refdesc, cn"FontName")
end
#Not implemented for CIDFonts
get_font_name(cosfont::CosObject, ::CIDWidth) = cn""
function get_font_name(cosfont::CosObject, x)
    basefname = get(cosfont, cn"BaseFont")
    basefname === CosNull &&
        error("Non-standard 14 fonts having no font descriptor")
    return basefname
end

function get_character_code(name::CosName, pdfont::PDFont)
    length(pdfont.glyph_name_id) > 0 &&
        return get(pdfont.glyph_name_id, name, INIT_CODE(pdfont.widths))
    return get_character_code(name, pdfont.widths)
end

get_character_code(name::CosName, w::CIDWidth) =
    UInt16(get(AGL_Glyph_to_Unicode, name, INIT_CODE(w)))

get_character_code(name::CosName, w) =
    get(GlyphName_to_STDEncoding, name, INIT_CODE(w))

get_encoded_string(s, pdfont::PDFont) = get_encoded_string(s, pdfont.fum)

get_char(barr, w) = iterate(barr)
function get_char(barr, w::CIDWidth)
    next = iterate(barr)
    next === nothing && return nothing
    (b1, state) = next
    next = iterate(barr, state)
    @assert next !== nothing "Error in obtaining character data"
    (b2, state) = next
    return (UInt16(b1*0x0100 + b2), state)
end

function get_char(barr, state, w::CIDWidth)
    next = iterate(barr, state)
    next === nothing && return nothing
    (b1, state) = next
    next = iterate(barr, state)
    @assert next !== nothing "Error in obtaining character data"
    (b2, state) = next
    return (UInt16(b1*0x0100 + b2), state)
end
get_char(barr, state, w) = iterate(barr, state)

function get_string_width(barr::Vector{UInt8}, widths, pc, tfs, tj, tc, tw)
    totalw = 0f0
    next = get_char(barr, widths)
    while next !== nothing
        c, st = next
        w = get_character_width(c, widths)
        kw = get_kern_width(pc, c, widths)
        w = (w - tj)*tfs / 1000f0 + ((c == SPACE_CODE(widths)) ? tw : tc)
        w += kw
        pc = c
        tj = 0f0
        totalw += w
        next = get_char(barr, st, widths)
    end
    return totalw
end

function get_TextBox(ss::Vector{Union{CosString,CosNumeric}},
                     pdfont::PDFont,
                     tfs, tc, tw, th)
    totalw = 0f0
    tj = 0f0
    text = ""
    for s in ss
        if s isa CosString
            prev_char = INIT_CODE(pdfont.widths)
            t = String(get_encoded_string(s, pdfont))
            if (-tj) > 180 && length(t) > 0 && t[1] != ' ' &&
                length(text) > 0 && text[end] != ' '
                text *= " "
            end
            text *= t
            barr = Vector{UInt8}(s)
            totalw += get_string_width(barr, pdfont.widths, prev_char,
                                       tfs, tj, tc, tw)
            tj = 0f0
        elseif s isa CosNumeric
            tj = s |> get |> Float32
        end
    end
    totalw *= th
    return text, totalw, tfs
end

function get_character_width(cid::UInt16, w::CIDWidth)
    itv = intersect(w.imap, Interval(cid, cid))
    length(itv) == 0 && return w.dw
    return itv[1][2]
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module ModiaMath.DAE

Interface between the [`ModiaMath.SimulationEngine`](@ref) and the index 1 DAE model.
A DAE model is a `struct` that has a required field
`simulationState::`[`ModiaMath.SimulationState`](@ref)
in which the main properties of the DAE model are reported to the simulation engine:

```julia
# DAE model ModelName
mutable struct ModelName <: ModiaMath.AbstractSimulationModel
    simulationState::ModiaMath.SimulationState

    # other definitions (e.g. parameters of model)
end
```

The following functions can be called in the DAE model to inquire
**information about the simulation state**:

- [`ModiaMath.getTime`](@ref)
- [`ModiaMath.getStartTime`](@ref)
- [`ModiaMath.getStopTime`](@ref)
- [`ModiaMath.getTolerance`](@ref)
- [`ModiaMath.get_is_constraint`](@ref)
- [`ModiaMath.compute_der_fc`](@ref)
- [`ModiaMath.isInitial`](@ref)
- [`ModiaMath.isTerminal`](@ref)
- [`ModiaMath.isEvent`](@ref)
- [`ModiaMath.isZeroCrossing`](@ref)
- [`ModiaMath.isAfterSimulationStart`](@ref)
- [`ModiaMath.isStoreResult`](@ref)
- [`ModiaMath.isLogInfos`](@ref)
- [`ModiaMath.isLogWarnings`](@ref)
- [`ModiaMath.isLogEvents`](@ref)

The following functions can be called in the DAE model to set
**properties in the simulation engine**:

- [`ModiaMath.setNominal!`](@ref)
- [`ModiaMath.setNextEvent!`](@ref)
- [`ModiaMath.positive!`](@ref)
- [`ModiaMath.negative!`](@ref)
- [`ModiaMath.change!`](@ref)
- [`ModiaMath.edge!`](@ref)

The following functions can be either called in the DAE model
or they can be called on a simulation model
(before or after `ModiaMath.simulate!(simulationModel, ...)` is called).

- [`ModiaMath.logOn!`](@ref)
- [`ModiaMath.logOff!`](@ref)
- [`ModiaMath.setLogCategories!`](@ref)

# Main developer

[Martin Otter](https://rmc.dlr.de/sr/de/staff/martin.otter/),
[DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module DAE

export EventRestart, NoRestart, Restart, FullRestart, Terminate
export getVariableName, getResultNames, defaultVariableName

export getTime, getStartTime, getStopTime, getTolerance, get_is_constraint, compute_der_fc
export isInitial, isTerminal, isEvent, isZeroCrossing, isAfterSimulationStart, isStoreResult, setNominal!
export setNextEvent!, positive!, negative!, change!, edge!
export getSimulationResult

export SimulationState
export StructureOfDAE, DAE_LinearDerivativesAndConstraints, DAE_ExplicitDerivatives, DAE_NoSpecialStructure

# export InitInfo, EventInfo, reset!
# export EventHandler
# export getEventResidues!
# export reinitialize!, eventIteration!, initialize!, getResidues!
# export computeAndStoreResult!, terminate!, processEvent!, getEventIndicators

import ..ModiaMath
using  LinearAlgebra
import DataFrames

"""
    @enum EventRestart NoRestart Restart FullRestart Terminate

Define how to continue or restart integration after an event. Usually, `Restart`
should be used. Only in special cases, the other flags are useful.

- `NoRestart`, continue integration without restarting the integrator
- `Restart`, restart integrator
- `FullRestart`, restart integrator and optionally exchange simulationState (so dimensions may change)
- `Terminate`, terminate integration
"""
@enum EventRestart NoRestart Restart FullRestart Terminate

@enum VariableCategory Category_X Category_DERX Category_W

# include code
include("events.jl")
include("returnStructs.jl")
include("simulationState.jl")
include("functionsForUserModels.jl")

end


text:
export crotinvKernel, Kx
export dKernel, dcrotinvKernel

# TODO: combine with other group action (currently f(phi,x) not rot(phi)*x)
function Krot(phi)
    c = cos(phi)
    s = sin(phi)
    [c -s 0 0; s c 0 0; 0 0 c -s; 0 0 s c]
end

### invariant Kernel
struct crotinvKernel{K <: KernelFunctions.Kernel, GKP }  <: MatrixKernel
    kernel::K
    gkparams::GKP
end

# constructor
function crotinvKernel(ker::K; N::Int) where K
    gkparams = NamedTuple{(:x, :weights)}(gauss(N, 0, 2pi))
#     gkparams = zip(gkparams.x, gkparams.weights)
#     gkparams = Iterators.product(gkparams, gkparams)
    crotinvKernel{K, typeof(gkparams)}(ker, gkparams)
end

function (crker::crotinvKernel)(x1::Array{T,1},x2::Array{T,1}) where T <: Real
    phis = crker.gkparams.x

    rx1 = Krot.(phis) .* [x1]
    rx2 = Krot.(phis) .* [x2]

    return crker(rx1, rx2)
end

function (crker::crotinvKernel)(x1::Array{T,1},rx2::Array{Array{T,1},1}) where T <: Real
    phis = crker.gkparams.x

    rx1 = Krot.(phis) .* [x1]

    return crker(rx1, rx2)
end

function (crker::crotinvKernel)(rx1::Array{Array{T,1},1}, rx2::Array{Array{T,1},1}) where T <: Real
    w = crker.gkparams.weights
    ker = crker.kernel
    iterIdx = Iterators.product(1:length(w), 1:length(w))

    tmp = mapreduce(iI -> w[iI[1]] * ker(rx1[iI[1]], rx2[iI[2]]) * w[iI[2]], + , iterIdx)
    return tmp
end

struct Kx{K, rZ}
    rker::K
    rotZ::rZ
end

function Kx(ker::K, Z) where K <: MatrixKernel
    phis = ker.gkparams.x
    tmpRots = Krot.(phis)
    rotZ = [tmpRots .* [z] for z in Z]

    Kx{K, typeof(rotZ)}(ker, rotZ)
end

function (Kx::Kx)()
    rZ = Kx.rotZ
#     w = Kx.rker.gkparams.weights
#     ker = Kx.rker.kernel
#     iterIdx = Iterators.product(1:length(w), 1:length(w))

    N = length(Kx.rotZ)
    K = zeros(N, N)

    Threads.@threads for i in 1:N
        for j in 1:i
            tmp = Kx.rker(rZ[i], rZ[j])
#             tmp = mapreduce(iI -> w[iI[1]] * ker(rZ[i][iI[1]], rZ[j][iI[2]]) * w[iI[2]], + , iterIdx)
            K[i,j] = K[j,i] = tmp
        end
    end
    return K
end

function (Kx::Kx)(x)
    rZ = Kx.rotZ
#     w = Kx.rker.gkparams.weights
    ker = Kx.rker.kernel
#     iterIdx = Iterators.product(1:length(w), 1:length(w))

    # phis = Kx.rker.gkparams.x
    # rx = Krot.(phis) .* [x]

    D = length(ker(x,x))

    N = length(Kx.rotZ)
    K = Zygote.Buffer([0.], D,N)
#     zeros(1, N)

#     Threads.@threads for i in 1:N
    for i in 1:N
#         tmp = mapreduce(iI -> w[iI[1]] * ker(rx[iI[1]], rZ[i][iI[2]]) * w[iI[2]], + , iterIdx)
        tmp = Kx.rker(x, rZ[i])
        for (j, val) in enumerate(tmp)
            K[j,i] = val
        end
#         display(K[:, i])
    end
    return copy(K)
end

# some extra stuff
function (mker::crotinvKernel)(w::Array{T, 1}) where {T <: Real}
    typeof(mker)(mker.kernel(w), mker.gkparams)
end

function (mker::crotinvKernel)(x1, x2)
    typeof(mker)(mker.kernel(w), mker.gkparams)
end

import KernelFunctions: kernelmatrix
function kernelmatrix(ker::crotinvKernel, X)
    tmpKx = Kx(ker, X)
    return tmpKx()
end

####################################################
### The derivative kernels
## The base kernel
struct dKernel{dK} <: KernelFunctions.Kernel
    dkernel::dK
end

# constructor
function dKernel(ker::K) where K <: KernelFunctions.Kernel
    dker(x1,x2) = gradient(x1->ker(x1,x2), x1)
    dKernel{typeof(dker)}(dker)
end

function (dker::dKernel)(x1,x2)
    dker.dkernel(x1,x2)[1]
end

## the invariant derivative kernel
struct dcrotinvKernel{K <: KernelFunctions.Kernel, GKP, DI }  <: MatrixKernel
    kernel::K
    gkparams::GKP
    derividx::DI
end
# TODO: possibly make things a little more efficient by only computing part of the kernel

function dcrotinvKernel(ker::K; N::Int, di=:) where K <: dKernel
    gkparams = NamedTuple{(:x, :weights)}(gauss(N, 0, 2pi))
    dcrotinvKernel{K, typeof(gkparams), typeof(di)}(ker, gkparams, di)
end

function (crker::dcrotinvKernel)(x1::Array{T,1},x2::Array{T,1}) where T <: Real
    phis = crker.gkparams.x
    R = GaussianProcessODEs.Krot.(phis)
    rx1 = R .* [x1]
    rx2 = R .* [x2]
    return crker(R, rx1, rx2)
end

function (crker::dcrotinvKernel)(x1::Array{T,1},rx2::Array{Array{T,1},1}) where T <: Real
    phis = crker.gkparams.x
    R = GaussianProcessODEs.Krot.(phis)
    rx1 = R .* [x1]
    return crker(R, rx1, rx2)
end

function _rh(R, di)
    transpose(R)[di, :]
end

function (crker::dcrotinvKernel)(R, rx1::Array{Array{T,1},1}, rx2::Array{Array{T,1},1}) where T <: Real
    w = crker.gkparams.weights
    dker = crker.kernel
    di = crker.derividx
    iterIdx = Iterators.product(1:length(w), 1:length(w))
#
    tmp = mapreduce(iI -> w[iI[1]] * _rh(R[iI[1]], di) * dker(rx1[iI[1]], rx2[iI[2]]) * w[iI[2]], + , iterIdx)
    return tmp
end


text:
struct Prop3DAcoIsoDenQ_DEO2_FDTD
    p::Ptr{Cvoid}
end

function Prop3DAcoIsoDenQ_DEO2_FDTD(;
        nz=0,
        ny=0,
        nx=0,
        nsponge=60,
        nbz=512,
        nby=8,
        nbx=8,
        dz=5.0,
        dy=8.0,
        dx=5.0,
        dt=0.001,
        nthreads=Sys.CPU_THREADS,
        freesurface=true,
        freqQ=5.0,
        qMin=0.1,
        qInterior=100.0)
    nz,ny,nx,nsponge,nbz,nby,nbx,nthreads = map(x->round(Int,x), (nz,ny,nx,nsponge,nbz,nby,nbx,nthreads))
    dz,dy,dx,dt = map(x->Float32(x), (dz,dy,dx,dt))

    @assert nx > 0
    @assert ny > 0
    @assert nz > 0
    @assert nsponge > 0
    @assert nthreads > 0
    @assert nbx > 0
    @assert nby > 0
    @assert nbz > 0

    fs = freesurface ? 1 : 0

    p = ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_alloc, libprop3DAcoIsoDenQ_DEO2_FDTD),
        Ptr{Cvoid},
        (Clong, Clong,    Clong, Clong, Clong, Clong,   Cfloat, Cfloat, Cfloat, Cfloat, Clong, Clong, Clong),
         fs,    nthreads, nx,    ny,    nz,    nsponge, dx,     dy,     dz,     dt,     nbx,   nby,   nbz)

    ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_SetupDtOmegaInvQ, libprop3DAcoIsoDenQ_DEO2_FDTD),
        Cvoid,
        (Ptr{Cvoid}, Cfloat, Cfloat, Cfloat),
         p,          freqQ,  qMin,   qInterior)

    Prop3DAcoIsoDenQ_DEO2_FDTD(p)
end

free(prop::Prop3DAcoIsoDenQ_DEO2_FDTD) = ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_free, libprop3DAcoIsoDenQ_DEO2_FDTD), Cvoid, (Ptr{Cvoid},), prop.p)

function size(prop::Prop3DAcoIsoDenQ_DEO2_FDTD)
    nz = ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_getNz, libprop3DAcoIsoDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    ny = ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_getNy, libprop3DAcoIsoDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    nx = ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_getNx, libprop3DAcoIsoDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    (nz,ny,nx)
end

for _f in (:V, :B, :PSpace, :PCur, :POld, :DtOmegaInvQ)
    symf = "Prop3DAcoIsoDenQ_DEO2_FDTD_get" * string(_f)
    @eval $(_f)(prop::Prop3DAcoIsoDenQ_DEO2_FDTD) = unsafe_wrap(Array, ccall(($symf, libprop3DAcoIsoDenQ_DEO2_FDTD), Ptr{Float32}, (Ptr{Cvoid},), prop.p), size(prop), own=false)
end

abstract type Prop3DAcoIsoDenQ_DEO2_FDTD_Model end
struct Prop3DAcoIsoDenQ_DEO2_FDTD_Model_V <: Prop3DAcoIsoDenQ_DEO2_FDTD_Model end
struct Prop3DAcoIsoDenQ_DEO2_FDTD_Model_B <: Prop3DAcoIsoDenQ_DEO2_FDTD_Model end
struct Prop3DAcoIsoDenQ_DEO2_FDTD_Model_VB <: Prop3DAcoIsoDenQ_DEO2_FDTD_Model end

propagateforward!(prop::Prop3DAcoIsoDenQ_DEO2_FDTD) = ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_TimeStep, libprop3DAcoIsoDenQ_DEO2_FDTD), Cvoid, (Ptr{Cvoid},), prop.p)
propagateadjoint!(prop::Prop3DAcoIsoDenQ_DEO2_FDTD) = propagateforward!(prop) # self-adjoint

scale_spatial_derivatives!(prop::Prop3DAcoIsoDenQ_DEO2_FDTD) =
    ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_ScaleSpatialDerivatives, libprop3DAcoIsoDenQ_DEO2_FDTD), Cvoid, (Ptr{Cvoid},), prop.p)

# v in model-space
function forwardBornInjection!(prop::Prop3DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop3DAcoIsoDenQ_DEO2_FDTD_Model_V, dmodel, wavefields)
    ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_ForwardBornInjection_V, libprop3DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat}),
         prop.p,     dmodel["v"], wavefields["pspace"])
end

# v,b in model-space
function forwardBornInjection!(prop::Prop3DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop3DAcoIsoDenQ_DEO2_FDTD_Model_VB, dmodel, wavefields)
    ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_ForwardBornInjection_VB, libprop3DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat},        Ptr{Cfloat}),
         prop.p,     dmodel["v"], dmodel["b"], wavefields["pold"], wavefields["pspace"])
end

# b in model-space
function forwardBornInjection!(prop::Prop3DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop3DAcoIsoDenQ_DEO2_FDTD_Model_B, dmodel, wavefields)
    ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_ForwardBornInjection_B, libprop3DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},        Ptr{Cfloat}),
         prop.p,     dmodel["b"], wavefields["pold"], wavefields["pspace"])
end

# v in model-space
function adjointBornAccumulation!(prop::Prop3DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop3DAcoIsoDenQ_DEO2_FDTD_Model_V, imagingcondition::ImagingConditionStandard, dmodel, wavefields)
    ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_AdjointBornAccumulation_V, libprop3DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat}),
         prop.p,     dmodel["v"], wavefields["pspace"])
end

# v,b in model-space
function adjointBornAccumulation!(prop::Prop3DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop3DAcoIsoDenQ_DEO2_FDTD_Model_VB, imagingcondition::ImagingConditionStandard, dmodel, wavefields)
    ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_AdjointBornAccumulation_VB, libprop3DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat},        Ptr{Cfloat}),
         prop.p,     dmodel["v"], dmodel["b"], wavefields["pold"], wavefields["pspace"])
end

# b in model-space
function adjointBornAccumulation!(prop::Prop3DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop3DAcoIsoDenQ_DEO2_FDTD_Model_B, imagingcondition::ImagingConditionStandard, dmodel, wavefields)
    ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_AdjointBornAccumulation_B, libprop3DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},        Ptr{Cfloat}),
         prop.p,     dmodel["b"], wavefields["pold"], wavefields["pspace"])
end

# v in model-space with FWI wavefield separation
function adjointBornAccumulation!(prop::Prop3DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop3DAcoIsoDenQ_DEO2_FDTD_Model_V, imagingcondition::ImagingConditionWaveFieldSeparationFWI, dmodel, wavefields)
    ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_AdjointBornAccumulation_wavefieldsep, libprop3DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},          Clong),
         prop.p,     dmodel["v"], wavefields["pspace"], 1)
end

# v in model-space with RTM wavefield separation
function adjointBornAccumulation!(prop::Prop3DAcoIsoDenQ_DEO2_FDTD, modeltype::Prop3DAcoIsoDenQ_DEO2_FDTD_Model_V, imagingcondition::ImagingConditionWaveFieldSeparationRTM, dmodel, wavefields)
    ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_AdjointBornAccumulation_wavefieldsep, libprop3DAcoIsoDenQ_DEO2_FDTD), Cvoid,
        (Ptr{Cvoid}, Ptr{Cfloat}, Ptr{Cfloat},         Clong),
        prop.p,     dmodel["v"], wavefields["pspace"], 0)
end

function show(io::IO, prop::Prop3DAcoIsoDenQ_DEO2_FDTD)
    nx = ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_getNx, libprop3DAcoIsoDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    ny = ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_getNy, libprop3DAcoIsoDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    nz = ccall((:Prop3DAcoIsoDenQ_DEO2_FDTD_getNz, libprop3DAcoIsoDenQ_DEO2_FDTD), (Clong), (Ptr{Cvoid},), prop.p)
    write(io, "Prop3DAcoIsoDenQ_DEO2_FDTD -- nx,ny,nz; $nx,$ny,$nz")
end

function tilesize(::Type{Prop3DAcoIsoDenQ_DEO2_FDTD}; nz, ny, nx, rngy=4:2:64, rngx=4:2:64, N=10, nthreads=Sys.CPU_THREADS-1)
    speeds = Matrix{Float64}(undef, length(rngy), length(rngx))
    for (iby,by) in enumerate(rngy), (ibx,bx) in enumerate(rngx)
        @info "by=$by, bx=$bx"
        prop = Prop3DAcoIsoDenQ_DEO2_FDTD(nz=nz, ny=ny, nx=nx, nbz=nz, nby=by, nbx=bx, nthreads=nthreads)

        v = WaveFD.V(prop)
        b = WaveFD.B(prop)
        pcur = WaveFD.PCur(prop)
        pold = WaveFD.POld(prop)

        v .= 1500
        b .= 1
        rand!(pcur)
        rand!(pold)

        t = 0.0
        set_zero_subnormals(true)
        for i = 1:N
            t += @elapsed propagateforward!(prop)
        end
        set_zero_subnormals(false)
        speeds[iby,ibx] = N*prod(size(prop))/1000/1000/t
        free(prop)
    end
    rngx,rngy,speeds
end

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
############################################################################################
# Constraints

############################################################################################
@testset "Constraints - Bijector" begin
    for output in outputtypes
        for flattentype in flattentypes
            flatdefault = FlattenDefault(; output = output, flattentype = flattentype)
            val = 2.
            constraint = Bijectors.bijector(Gamma(2,2))
            ReConstructor(constraint, val)
            reconstruct = ReConstructor(flatdefault, constraint, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test x_unflatAD isa eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test x_unflatAD2 isa eltype(x_flatAD)
################################################################################
# Transforms
            _transform, _inversetransform = construct_transform(constraint, val)
            x_transformed = _transform(val)
            x_inversetransformed = _inversetransform(x_transformed)

            transformer = TransformConstructor(constraint, val)
            val_uncon = unconstrain(transformer, val)
            val_con = constrain(transformer, val_uncon)

            @test transformer.unconstrain == _transform
            @test transformer.constrain == _inversetransform
            @test x_transformed == val_uncon
            @test x_inversetransformed == val_con

            @test x_inversetransformed isa typeof(val)
            val_unconstrained = unconstrain(_transform, val)
            val_constrained = constrain(_inversetransform, val_unconstrained)
            @test val_constrained isa typeof(val)
            @test val_constrained ≈ val atol = _TOL
            logabs = log_abs_det_jac(_transform, val)
            @test logabs isa AbstractFloat
            @test _check(_RNG, constraint, val)
        end
    end
end

############################################################################################
@testset "Constraints - Distribution" begin
    for output in outputtypes
        for flattentype in flattentypes
            flatdefault = FlattenDefault(; output = output, flattentype = flattentype)
            val = 2.
            constraint = Gamma(2,2)
            ReConstructor(constraint, val)
            reconstruct = ReConstructor(flatdefault, constraint, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test x_unflatAD isa eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test x_unflatAD2 isa eltype(x_flatAD)
################################################################################
# Transforms
            _transform, _inversetransform = construct_transform(constraint, val)
            x_transformed = _transform(val)
            x_inversetransformed = _inversetransform(x_transformed)

            transformer = TransformConstructor(constraint, val)
            val_uncon = unconstrain(transformer, val)
            val_con = constrain(transformer, val_uncon)

            @test transformer.unconstrain == _transform
            @test transformer.constrain == _inversetransform
            @test x_transformed == val_uncon
            @test x_inversetransformed == val_con

            @test x_inversetransformed isa typeof(val)
            val_unconstrained = unconstrain(_transform, val)
            val_constrained = constrain(_inversetransform, val_unconstrained)
            @test val_constrained isa typeof(val)
            @test val_constrained ≈ val atol = _TOL
            logabs = log_abs_det_jac(_transform, val)
            @test logabs isa AbstractFloat
            @test _check(_RNG, constraint, val)
        end
    end
end

############################################################################################
@testset "Constraints - Constrained" begin
    for output in outputtypes
        for flattentype in flattentypes
            flatdefault = FlattenDefault(; output = output, flattentype = flattentype)
            val = 2.
            constraint = Constrained(1.,3.)
            ReConstructor(constraint, val)
            reconstruct = ReConstructor(flatdefault, constraint, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test x_unflatAD isa eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test x_unflatAD2 isa eltype(x_flatAD)
################################################################################
# Transforms
            _transform, _inversetransform = construct_transform(constraint, val)
            x_transformed = _transform(val)
            x_inversetransformed = _inversetransform(x_transformed)

            transformer = TransformConstructor(constraint, val)
            val_uncon = unconstrain(transformer, val)
            val_con = constrain(transformer, val_uncon)

            @test transformer.unconstrain == _transform
            @test transformer.constrain == _inversetransform
            @test x_transformed == val_uncon
            @test x_inversetransformed == val_con

            @test x_inversetransformed isa typeof(val)
            val_unconstrained = unconstrain(_transform, val)
            val_constrained = constrain(_inversetransform, val_unconstrained)
            @test val_constrained isa typeof(val)
            @test val_constrained ≈ val atol = _TOL
            logabs = log_abs_det_jac(_transform, val)
            @test logabs isa AbstractFloat
            @test _check(_RNG, constraint, val)
        end
    end
end

############################################################################################
@testset "Constraints - Unconstrained" begin
    for output in outputtypes
        for flattentype in flattentypes
            flatdefault = FlattenDefault(; output = output, flattentype = flattentype)
            val = 2.
            constraint = Unconstrained()
            ReConstructor(constraint, val)
            reconstruct = ReConstructor(flatdefault, constraint, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test x_unflatAD isa eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test x_unflatAD2 isa eltype(x_flatAD)
################################################################################
# Transforms
            _transform, _inversetransform = construct_transform(constraint, val)
            x_transformed = _transform(val)
            x_inversetransformed = _inversetransform(x_transformed)

            transformer = TransformConstructor(constraint, val)
            val_uncon = unconstrain(transformer, val)
            val_con = constrain(transformer, val_uncon)

            @test transformer.unconstrain == _transform
            @test transformer.constrain == _inversetransform
            @test x_transformed == val_uncon
            @test x_inversetransformed == val_con

            @test x_inversetransformed isa typeof(val)
            val_unconstrained = unconstrain(_transform, val)
            val_constrained = constrain(_inversetransform, val_unconstrained)
            @test val_constrained isa typeof(val)
            @test val_constrained ≈ val atol = _TOL
            logabs = log_abs_det_jac(_transform, val)
            @test logabs isa AbstractFloat
            @test _check(_RNG, constraint, val)
        end
    end
end

############################################################################################
@testset "Constraints - Fixed" begin
    for output in outputtypes
        for flattentype in flattentypes
            flatdefault = FlattenDefault(; output = output, flattentype = flattentype)
            val = zeros(2,3,4)
            constraint = Fixed()
            ReConstructor(constraint, val)
            reconstruct = ReConstructor(flatdefault, constraint, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
#            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
#            @test x_unflatAD isa eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
#            @test x_unflatAD2 isa eltype(x_flatAD)
################################################################################
# Transforms
            _transform, _inversetransform = construct_transform(constraint, val)
            x_transformed = _transform(val)
            x_inversetransformed = _inversetransform(x_transformed)

            transformer = TransformConstructor(constraint, val)
            val_uncon = unconstrain(transformer, val)
            val_con = constrain(transformer, val_uncon)

            @test transformer.unconstrain == _transform
            @test transformer.constrain == _inversetransform
            @test x_transformed == val_uncon
            @test x_inversetransformed == val_con

            @test x_inversetransformed isa typeof(val)
            val_unconstrained = unconstrain(_transform, val)
            val_constrained = constrain(_inversetransform, val_unconstrained)
            @test val_constrained isa typeof(val)
            @test val_constrained ≈ val atol = _TOL
            logabs = log_abs_det_jac(_transform, val)
            @test logabs isa AbstractFloat
            @test _check(_RNG, constraint, val)
        end
    end
end

############################################################################################
@testset "Constraints - CorrelationMatrix" begin
    for output in outputtypes
        for flattentype in flattentypes
            flatdefault = FlattenDefault(; output = output, flattentype = flattentype)
            val = zeros(3,3)
            val[1,1] = val[2,2] = val[3,3] = 1.0
            val
            val[2,1] = val[1, 2] = 0.12
            val[3,1] = val[1, 3] = 0.13
            val[3,2] = val[2, 3] = 0.14
            val
            constraint = CorrelationMatrix()
            ReConstructor(constraint, val)
            reconstruct = ReConstructor(flatdefault, constraint, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
            @test length(x_flat) == 3
            @test x_flat == output.([0.12, 0.13, 0.14])
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test eltype(x_unflatAD) == eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test eltype(x_unflatAD2) == eltype(x_flatAD)

################################################################################
# Check if errors if flattened parameter not correct size
            @test_throws ArgumentError flatten(reconstruct, zeros(4,4))
            @test_throws ArgumentError flatten(reconstruct, zeros(2,2))
            @test_throws ArgumentError flattenAD(reconstruct, zeros(4,4))
            @test_throws ArgumentError flattenAD(reconstruct,  zeros(2,3))
            @test_throws ArgumentError unflatten(reconstruct, zeros(length(x_flat)+1))
            @test_throws ArgumentError unflatten(reconstruct, zeros(length(x_flat)-1))
            @test_throws ArgumentError unflattenAD(reconstruct, zeros(length(x_flat)+1))
            @test_throws ArgumentError unflattenAD(reconstruct, zeros(length(x_flat)-1))

################################################################################
# Transforms
            _transform, _inversetransform = construct_transform(constraint, val)
            x_transformed = _transform(val)
            x_inversetransformed = _inversetransform(x_transformed)

            transformer = TransformConstructor(constraint, val)
            val_uncon = unconstrain(transformer, val)
            val_con = constrain(transformer, val_uncon)

            @test transformer.unconstrain == _transform
            @test transformer.constrain == _inversetransform
            @test x_transformed == val_uncon
            @test x_inversetransformed == val_con

            val_unconstrained = unconstrain(_transform, val)
            #Upper triangular
            @test val_unconstrained[3,1] == 0.0
            val_constrained = constrain(_inversetransform, val_unconstrained)
            @test val_constrained isa typeof(val)
            @test val_constrained ≈ val atol = _TOL
            logabs = log_abs_det_jac(_transform, val)
            @test logabs isa AbstractFloat
            @test _check(_RNG, constraint, val)
        end
    end
end

############################################################################################
@testset "Constraints - CovarianceMatrix" begin
    for output in outputtypes
        for flattentype in flattentypes
            flatdefault = FlattenDefault(; output = output, flattentype = flattentype)
            val = zeros(3,3)
            val[1,1] = 1.5
            val[2,2] = 1.6
            val[3,3] = 1.7
            val
            val[2,1] = val[1, 2] = 0.12
            val[3,1] = val[1, 3] = 0.13
            val[3,2] = val[2, 3] = 0.14
            val
            constraint = CovarianceMatrix()
            ReConstructor(constraint, val)
            reconstruct = ReConstructor(flatdefault, constraint, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
            @test length(x_flat) == 6
            @test x_flat == output.([1.50, 0.12, 0.13, 1.60, 0.14, 1.70])
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test eltype(x_unflatAD) == eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test eltype(x_unflatAD2) == eltype(x_flatAD)

################################################################################
# Check if errors if flattened parameter not correct size
            @test_throws ArgumentError flatten(reconstruct, zeros(4,4))
            @test_throws ArgumentError flatten(reconstruct, zeros(2,2))
            @test_throws ArgumentError flattenAD(reconstruct, zeros(4,4))
            @test_throws ArgumentError flattenAD(reconstruct,  zeros(2,3))
            @test_throws ArgumentError unflatten(reconstruct, zeros(length(x_flat)+1))
            @test_throws ArgumentError unflatten(reconstruct, zeros(length(x_flat)-1))
            @test_throws ArgumentError unflattenAD(reconstruct, zeros(length(x_flat)+1))
            @test_throws ArgumentError unflattenAD(reconstruct, zeros(length(x_flat)-1))

################################################################################
# Transforms
            _transform, _inversetransform = construct_transform(constraint, val)
            x_transformed = _transform(val)
            x_inversetransformed = _inversetransform(x_transformed)

            transformer = TransformConstructor(constraint, val)
            val_uncon = unconstrain(transformer, val)
            val_con = constrain(transformer, val_uncon)

            @test transformer.unconstrain == _transform
            @test transformer.constrain == _inversetransform
            @test x_transformed == val_uncon
            @test x_inversetransformed == val_con

            val_unconstrained = unconstrain(_transform, val)
            #Lower triangular
            @test val_unconstrained[1,3] == 0.0
            val_constrained = constrain(_inversetransform, val_unconstrained)
            @test val_constrained isa typeof(val)
            @test val_constrained ≈ val atol = _TOL
            logabs = log_abs_det_jac(_transform, val)
            @test logabs isa AbstractFloat
            @test _check(_RNG, constraint, val)
################################################################################
# Custom Matrix flattening
            _tag = tag(val, false, true)
            _vals4 = flatten_Symmetric(val, _tag)
            @test length(_vals4) == 6
            @test all(Symmetric_from_flatten(_vals4, _tag) .≈ val)
        end
    end
end

############################################################################################
@testset "Constraints - Simplex" begin
    for output in outputtypes
        for flattentype in flattentypes
            flatdefault = FlattenDefault(; output = output, flattentype = flattentype)
            val = [.1, .2, .7]
            constraint = Simplex(val)
            ReConstructor(constraint, val)
            reconstruct = ReConstructor(flatdefault, constraint, val)
# Flatten
            x_flat = flatten(reconstruct, val)
            @test x_flat isa AbstractVector
            @test eltype(x_flat) == output
            @test length(x_flat) == 2
            @test x_flat == output.([.1, .2])
# Flatten AD
            x_flatAD = flattenAD(reconstruct, val)
            @test x_flatAD isa AbstractVector
            @test eltype(x_flatAD) == eltype(val)
# Unflatten
            x_unflat = unflatten(reconstruct, x_flat)
            @test x_unflat isa typeof(val)
# Unflatten AD
            x_unflatAD = unflattenAD(reconstruct, x_flat)
            @test eltype(x_unflatAD) == eltype(x_flat)
            x_unflatAD2 = unflattenAD(reconstruct, x_flatAD)
            @test eltype(x_unflatAD2) == eltype(x_flatAD)

################################################################################
# Check if errors if flattened parameter not correct size
            @test_throws ArgumentError flatten(reconstruct, zeros(4))
            @test_throws ArgumentError flatten(reconstruct, zeros(2))
            @test_throws ArgumentError flattenAD(reconstruct, zeros(4))
            @test_throws ArgumentError flattenAD(reconstruct,  zeros(2))
            @test_throws ArgumentError unflatten(reconstruct, zeros(length(x_flat)+1))
            @test_throws ArgumentError unflatten(reconstruct, zeros(length(x_flat)-1))
            @test_throws ArgumentError unflattenAD(reconstruct, zeros(length(x_flat)+1))
            @test_throws ArgumentError unflattenAD(reconstruct, zeros(length(x_flat)-1))

################################################################################
# Transforms
            _transform, _inversetransform = construct_transform(constraint, val)
            x_transformed = _transform(val)
            x_inversetransformed = _inversetransform(x_transformed)

            transformer = TransformConstructor(constraint, val)
            val_uncon = unconstrain(transformer, val)
            val_con = constrain(transformer, val_uncon)

            @test transformer.unconstrain == _transform
            @test transformer.constrain == _inversetransform
            @test x_transformed == val_uncon
            @test x_inversetransformed == val_con

            @test x_inversetransformed isa typeof(val)
            val_unconstrained = unconstrain(_transform, val)
            val_constrained = constrain(_inversetransform, val_unconstrained)
            @test val_constrained isa typeof(val)
            @test sum(val_constrained .≈ val) == 3
            @test sum(val_constrained) ≈ 1.0
            logabs = log_abs_det_jac(_transform, val)
            @test logabs isa AbstractFloat
            @test _check(_RNG, constraint, val)
################################################################################
# Custom Matrix flattening
            _vals4 = flatten_Simplex(val)
            @test length(_vals4) == length(val)-1
            _vals5 = Simplex_from_flatten!(zeros(length(val)), _vals4)
            _vals6 = Simplex_from_flatten(_vals4)
            @test sum(_vals5) ≈ 1.0
            @test sum(_vals6) ≈ 1.0
            @test all(output.(val) .≈ _vals5)
            @test all(output.(val) .≈ _vals6)
        end
    end
end


text:
###################################################################
#  Mark-related definitions & functions
###################################################################

########### mark-level configuration ##############

config_mark_spec = Dict(
  :filled           => Bool,
  :color            => AbstractString,
  :fill             => AbstractString,
  :stroke           => AbstractString,
  :opacity          => Real,
  :fillOpacity      => Real,
  :strokeOpacity    => Real,
  :strokeWidth      => Real,
  :strokeDash       => Vector,
  :strokeDashOffset => Vector,
  :stacked          => AbstractString,
  :interpolate      => AbstractString,
  :tension          => Real,
  :orient           => AbstractString,
  :barSize          => Real,
  :shape            => AbstractString,
  :size             => Real,
  :tickSize         => Real,
  :tickThickness    => Real,
  :angle            => Real,
  :align            => AbstractString,
  :baseline         => AbstractString,
  :dx               => Real,
  :dy               => Real,
  :radius           => Real,
  :theta            => Real,
  :font             => AbstractString,
  :fontSize         => Real,
  :fontStyle        => AbstractString,
  :fontWeight       => AbstractString,
  :text             => AbstractString,
  :format           => AbstractString,
  :shortTimeLabels  => Bool
)

"""
Sets mark properties

config_mark(;properties...)

- `filled` (`Bool`) :	Whether the shape's color should be used as fill color instead of stroke color. See mark for a usage example. Default value: true for all marks except point and false for point. Applicable for: bar, point, circle, square, and area marks.
- `color` (`AbstractString`) :	The color of the mark – either fill or stroke color based on the filled mark config. Default value: ■ blue (""#4682b4")
- `fill` (`AbstractString`) :	The fill color. This config will be overridden by color channel’s specified or mapped values if filled is true. Default value: (None)
- `stroke` (`AbstractString`) :	The stroke color. This config will be overridden by color channel’s specified or mapped values if filled is false. Default value: (None)

Opacity

- `opacity` (`Real`) : The overall opacity (value between [0,1]). Default value: 0.7 for non-aggregate plots with point, tick, circle, or square marks and 1 otherwise.
- `fillOpacity` (`Real`) :	The fill opacity (value between [0,1]). Default value: 1
- `strokeOpacity` (`Real`) :	The stroke opacity (value between [0,1]). Default value: 1

Stroke Style

- `strokeWidth` (`Real`) : The stroke width, in pixels.
- `strokeDash` (`Vector{Real}`) : An array of alternating stroke, space lengths for creating dashed or dotted lines.
- `strokeDashOffset` (`Vector{Real}`) : The offset (in pixels) into which to begin drawing with the stroke dash array.

Stacking (for Bar and Area)

- `stacked` (`AbstractString`) :	Stacking modes for bar and area marks.

  • zero - stacking with baseline offset at zero value of the scale (for creating typical stacked bar and area chart).
  • normalize - stacking with normalized domain (for creating normalized stacked bar and area chart).
  • center - stacking with center baseline (for streamgraph).
  • none - No-stacking. This will produces layered bar and area chart. Default value: zero if applicable.

Interpolation (for Line and Area Marks)

- `interpolate` (`AbstractString`) :	The line interpolation method to use. One of "linear", "step-before", "step-after", "basis", "basis-open", "basis-closed", "bundle", "cardinal", "cardinal-open", "cardinal-closed", "monotone". For more information about each interpolation method, please see D3’s line interpolation.
- `tension` (`Real`) :	Depending on the interpolation type, sets the tension parameter. (See D3’s line interpolation.)

Orientation (for Bar, Tick, Line, and Area Marks)

- `orient` (`AbstractString`) :	The orientation of a non-stacked bar, area, and line charts. The value is either "horizontal", or "vertical" (default). For bar and tick, this determines whether the size of the bar and tick should be applied to x or y dimension. For area, this property determines the orient property of the Vega output. For line, this property determines the path order of the points in the line if path channel is not specified. For stacked charts, this is always determined by the orientation of the stack; therefore explicitly specified value will be ignored.

Bar Config

- `barSize` (`Real`) :	The size of the bars (width for vertical bar charts and height for horizontal bar chart). Default value: bandSize-1 if the bar’s x or y axis is an ordinal scale. (This provides 1 pixel offset between bars.) and 2 for if both x and y scales have linear scales.

Point Config

- `shape` (`AbstractString`) : The symbol shape to use. One of "circle", "square", "cross", "diamond", "triangle-up", or "triangle-down". Default value: "circle"

Point Size Config (for Point, Circle, and Square Marks)

- `size` (`Real`) :	The pixel area each the point. For example: in the case of circles, the radius is determined in part by the square root of the size value. Default value: 30

Tick Config

- `tickSize` (`Real`) : The size of the ticks (height of the ticks for horizontal dot plots and strip plots and width of the ticks for vertical dot plots and strip plots). Default value: 2/3*bandSize (This will provide offset between band equals to the width of the tick.)
- `tickThickness` (`Real`) : Thickness of the tick mark. Default value: 1

Text Config

Text Position

- `angle` (`Real`) : The rotation angle of the text, in degrees.
- `align` (`AbstractString`) :	The horizontal alignment of the text. One of left, right, center.
- `baseline` (`AbstractString`) :	The vertical alignment of the text. One of top, middle, bottom.
- `dx` (`Real`) :	The horizontal offset, in pixels, between the text label and its anchor point. The offset is applied after rotation by the angle property.
- `dy` (`Real`) :	The vertical offset, in pixels, between the text label and its anchor point. The offset is applied after rotation by the angle property.
- `radius` (`Real`) :	Polar coordinate radial offset, in pixels, of the text label from the origin determined by the x and y properties.
- `theta` (`Real`) :	Polar coordinate angle, in radians, of the text label from the origin determined by the x and y properties. Values for theta follow the same convention of arc mark startAngle and endAngle properties: angles are measured in radians, with 0 indicating “north”.

Font Style

- `font` (`AbstractString`) :	The typeface to set the text in (e.g., Helvetica Neue).
- `fontSize` (`Real`) : : The font size, in pixels. The default value is 10.
- `fontStyle` (`AbstractString`) :	The font style (e.g., italic).
- `fontWeight` (`AbstractString`) :	The font weight (e.g., bold).

Text Value and Format

- `text` (`AbstractString`) :	Placeholder text if the text channel is not specified ("Abc" by default).
- `format` (`AbstractString`) :	The formatting pattern for text value. If not defined, this will be determined automatically
- `shortTimeLabels` (`Bool`) : Whether month names and weekday names should be abbreviated.

"""
config_mark(;properties...) = _mkvis((:config, :mark), config_mark_spec, properties)

######## Mark type choice ################

mark_spec = Dict(:mark => AbstractString)

mark_bar()    = _mkvis(tuple(), mark_spec, [(:mark, "bar")])
mark_circle() = _mkvis(tuple(), mark_spec, [(:mark, "circle")])
mark_square() = _mkvis(tuple(), mark_spec, [(:mark, "square")])
mark_tick()   = _mkvis(tuple(), mark_spec, [(:mark, "tick")])
mark_line()   = _mkvis(tuple(), mark_spec, [(:mark, "line")])
mark_area()   = _mkvis(tuple(), mark_spec, [(:mark, "area")])
mark_point()  = _mkvis(tuple(), mark_spec, [(:mark, "point")])
mark_text()   = _mkvis(tuple(), mark_spec, [(:mark, "text")])


text:
using PyCall
using DataFrames
using Pandas: query, read_csv, read_excel, read_pickle
using RCall
using Pipe
#using Winston

using TSne
using Clustering
using StatsBase

@pyimport re
@pyimport pandas as pd
@pyimport numpy as np
@pyimport pyspark.sql as sparksql

@pyimport scipy.stats as sps

@pyimport sklearn.preprocessing as sklpp

unshift!(PyVector(pyimport("sys")["path"]),
            "/home/benmo/OneDrive/GitHub/PyPackages/Algos")
@pyimport ModelsML

include("./myPandasIO.jl")

function fw_theorem(RᵢRₘ,i="AAPL",m="SP500";q=0.005)

    for stock in RᵢRₘ[:columns][:tolist]()
        qᵤ = py"$RᵢRₘ[$stock].quantile(1-$q)"
        qₗ = py"$RᵢRₘ[$stock].quantile($q)"
        RᵢRₘ[stock] = py"$RᵢRₘ[$stock].mask($RᵢRₘ[$stock] > $qᵤ, $qᵤ)"
        RᵢRₘ[stock] = py"$RᵢRₘ[$stock].mask($RᵢRₘ[$stock] < $qₗ, $qₗ)"
    end

    X = [ones(RᵢRₘ[m][:shape][1]) RᵢRₘ[m][:values]]
    Mᵢ = I - X*inv(X'*X)*X'

    stocks = py"list(filter(lambda x: x != $m, $RᵢRₘ.columns))"

    R̃ᵢ = RᵢRₘ[:copy]()

    for stock in stocks
        R̃ᵢ[stock] = Mᵢ * py"$RᵢRₘ[$stock].values"
    end
        R̃ᵢ = py"$R̃ᵢ[$stocks]"
end

function laglead(df; lag=15, lead=5)
    df = df[:dropna]()
    newDict = Dict()
    for col in df[:columns]
        newDict[col] = pd.DataFrame([])
        for i in -lag:lead
            newDict[col] = pd.concat((newDict[col],
                                        df[col][:shift](-i)),axis=1)
        end
        newDict[col] = newDict[col][:dropna]()
    end
    return newDict
end

function classifyRⱼ(βⱼ; perplx=50, nclust=16)
    dims = tsne(Matrix(βⱼ), 2, perplx)
    km1 = kmeans(transpose(dims), nclust; maxiter=100000, display=:iter)

    dimsT = dims'
    coeffs = [dimsT[:,i] for i in 1:length(dimsT[1,1:end])]
    tsneDF = DataFrame(Dict(zip(names(βⱼ), coeffs)))
    classesRow  = DataFrame(Dict(zip(names(βⱼ), km1.assignments')))

    classes = [tsneDF; classesRow]

    stockCats = hcat(
            DataFrames.columns(
            classes[3,:]) .|> (x) -> x[1],
            string.(names(classes))
            )
    return stockCats, dims, km1
end

function clean_Vol(βⱼ,rtns, stocksVol; startY=2013)
    dates = rtns[:index][:tolist]()

    #histogram(sample(techRtnJ[!isnan(techRtnJ[2]),2],500), opacity=0.45)
    #histogram!(sample(techRtnJ[!isnan(techRtnJ[1]),1],500), opacity=0.45)
    startY = string(startY)
    stocksVol = py"""$stocksVol[$startY:][$βⱼ.columns]"""
    isNa = stocksVol[:isna]()[:sum]()

    stocksVol = py"$stocksVol[$isNa[($isNa < 30)].index]"
    stocksVol = stocksVol[:dropna]()

    stocksVolμ = hcat(stocksVol[:mean]()[:tolist](),
                    stocksVol[:mean]()[:index][:tolist]())

    #histogram(techVolμ[:,1])

    sVol = stocksVolμ[stocksVolμ[:,1] .< 2.5e5,:]
    mVol = stocksVolμ[2.5e5 .< stocksVolμ[:,1] .< 1e7,:]
    lVol =  stocksVolμ[stocksVolμ[:,1] .> 1e7,:]

    nacols = py"$βⱼ.dropna(axis=1,how='all').columns.tolist()"
    βⱼ = py"$βⱼ[$nacols].loc[$nacols]"

    sVol_pd = py"$βⱼ[$(sVol[:,2])].loc[$(sVol[:,2])].dropna()"
    mVol_pd = py"$βⱼ[$(mVol[:,2])].loc[$(mVol[:,2])].dropna()"
    lVol_pd =  py"$βⱼ[$(lVol[:,2])].loc[$(lVol[:,2])].dropna()"

    dupcols = sVol_pd[:drop_duplicates]()[:index][:tolist]()
    sVol_pd = py"$sVol_pd[$dupcols].loc[$dupcols]"

    dupcols = mVol_pd[:drop_duplicates]()[:index][:tolist]()
    mVol_pd = py"$mVol_pd[$dupcols].loc[$dupcols]"

    dupcols = lVol_pd[:drop_duplicates]()[:index][:tolist]()
    lVol_pd = py"$lVol_pd[$dupcols].loc[$dupcols]"

    allstocks = []
    for x in [sVol_pd, mVol_pd, lVol_pd]
        allstocks = vcat(allstocks, x[:columns][:tolist]())
    end

    stocksRtn = py"$rtns[$allstocks]"
    βⱼ = py"$βⱼ[$allstocks]"

    return stocksRtn, stocksVol, sVol_pd, mVol_pd, lVol_pd

end

function sampleStocks(stockNames, sVol_pd, mVol_pd, lVol_pd, ωₛ, ωₘ, ωₗ; N=10)
    namesₛ, namesₘ, namesₗ =  (sVol_pd, mVol_pd, lVol_pd) .|> (x) ->
                        x[:columns][:tolist]()
    namesₛ, namesₘ, namesₗ = (namesₛ, namesₘ, namesₗ) .|>
                                        (x -> filter((z) -> z ∈ stockNames, x))

    ωₛ, ωₘ, ωₗ = (namesₛ, namesₘ, namesₗ) |> (x) ->  (ωₛ/length(x[1]), ωₘ/length(x[2]), ωₗ/length(x[3]))

    ωₛ = fill(ωₛ, length(namesₛ))
    ωₘ = fill(ωₘ, length(namesₘ))
    ωₗ = fill(ωₗ, length(namesₗ))

    samplex = [sample([namesₛ; namesₘ; namesₗ],Weights([ωₛ; ωₘ; ωₗ])) for i in 1:N]

end

max_stars_repo_name: benzomo/PyPackages

max_issues_repo_name: benzomo/PyPackages

max_forks_repo_name: benzomo/PyPackages


text:
function _prod_diagonals!(Y, A, B)
    @inbounds for j ∈ 1:size(A, 1)
        Y[j] = 0
        @inbounds for i ∈ 1:size(A, 2)
            Y[j] += A[j, i] * B[i, j]
        end
    end
    Y
end

function random_rotation(p)
    mat = randn(p, p)
    qr(mat).Q
end


text:
"""
    boundary_conditions_x!(density,bcond_x)

Apply the boundary conditions `bcond_x` on the x-axis of the matrix `density`.
Currently, only the boundary conditions `"periodic"`, `"Neumann"` and `"reflecting"` are implemented.

"""
function boundary_conditions_x!(density,bcond_x)
    if bcond_x == "periodic"
        density[1,:] = density[end-1,:]
        density[end,:] = density[2,:]
    elseif bcond_x == "Neumann"
        density[1,:] = density[2,:]
        density[end,:] = density[end-1,:]
    elseif bcond_x == "reflecting"
        density[1,:] = -density[2,:]
        density[end,:] = -density[end-1,:]
    else
        ArgumentError("Boundary condition not defined!")
    end
end

"""
    boundary_conditions_y!(density,bcond_y)

Apply the boundary conditions `bcond_y` on the y-axis of the matrix `density`.
Currently, only the boundary conditions `"periodic"`, `"Neumann"` and `"reflecting"` are implemented.
"""
function boundary_conditions_y!(density,bcond_y)
    if bcond_y == "periodic"
        density[:,1] = density[:,end-1]
        density[:,end] = density[:,2]
    elseif bcond_y == "Neumann"
        density[:,1] = density[:,2]
        density[:,end] = density[:,end-1]
    elseif bcond_y == "reflecting"
        density[:,1] = -density[:,2]
        density[:,end] = -density[:,end-1]
    else
        ArgumentError("Boundary condition not defined!")
    end
end

"""
    boundary_conditions_xy!(density,bcond_x,bcond_y)

Apply successively the functions `boundary_conditions_x!` and `boundary_conditions_x!` to
the matrix `density`.
"""
function boundary_conditions_xy!(density,bcond_x,bcond_y)
    boundary_conditions_x!(density,bcond_x)
    boundary_conditions_y!(density,bcond_y)
end

"""
    boundary_conditions_ρuv!(ρ,u,v,bcond_x,bcond_y)

Apply the boundary conditions to the data `(ρ,u,v)`.
"""
function boundary_conditions_ρuv!(ρ,u,v,bcond_x,bcond_y)
    if bcond_x == "periodic" || bcond_x == "Neumann"
        boundary_conditions_x!(ρ,bcond_x)
        boundary_conditions_x!(u,bcond_x)
        boundary_conditions_x!(v,bcond_x)
    elseif bcond_x == "reflecting"
        boundary_conditions_x!(ρ,"Neumann")
        boundary_conditions_x!(u,"reflecting")
        boundary_conditions_x!(v,"Neumann")
    else
        ArgumentError("Boundary condition not defined!")
    end

    if bcond_y == "periodic" || bcond_y == "Neumann"
        boundary_conditions_y!(ρ,bcond_y)
        boundary_conditions_y!(u,bcond_y)
        boundary_conditions_y!(v,bcond_y)
    elseif bcond_y == "reflecting"
        boundary_conditions_y!(ρ,"Neumann")
        boundary_conditions_y!(u,"Neumann")
        boundary_conditions_y!(v,"reflecting")
    else
        ArgumentError("Boundary condition not defined!")
    end
end


text:
"""
    HMCDA(n_iters::Int, n_adapt::Int, delta::Float64, lambda::Float64)

Hamiltonian Monte Carlo sampler wiht Dual Averaging algorithm.

Usage:

```julia
HMCDA(1000, 200, 0.65, 0.3)
```

Example:

```julia
# Define a simple Normal model with unknown mean and variance.
@model gdemo(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0, sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end

sample(gdemo([1.5, 2]), HMCDA(1000, 200, 0.65, 0.3))
```
"""
mutable struct HMCDA{T} <: Hamiltonian
  n_iters   ::  Int       # number of samples
  n_adapt   ::  Int       # number of samples with adaption for epsilon
  delta     ::  Float64   # target accept rate
  lambda    ::  Float64   # target leapfrog length
  space     ::  Set{T}    # sampling space, emtpy means all
  gid       ::  Int       # group ID
end
function HMCDA(n_adapt::Int, delta::Float64, lambda::Float64, space...)
    _space = isa(space, Symbol) ? Set([space]) : Set(space)
    return HMCDA(1, n_adapt, delta, lambda, _space, 0)
end
function HMCDA(n_iters::Int, delta::Float64, lambda::Float64)
    n_adapt_default = Int(round(n_iters / 2))
    n_adapt = n_adapt_default > 1000 ? 1000 : n_adapt_default
    return HMCDA(n_iters, n_adapt, delta, lambda, Set(), 0)
end
function HMCDA(alg::HMCDA, new_gid::Int)
    return HMCDA(alg.n_iters, alg.n_adapt, alg.delta, alg.lambda, alg.space, new_gid)
end
HMCDA{T}(alg::HMCDA, new_gid::Int) where {T} = HMCDA(alg, new_gid)
function HMCDA(n_iters::Int, n_adapt::Int, delta::Float64, lambda::Float64)
    return HMCDA(n_iters, n_adapt, delta, lambda, Set(), 0)
end
function HMCDA(n_iters::Int, n_adapt::Int, delta::Float64, lambda::Float64, space...)
    _space = isa(space, Symbol) ? Set([space]) : Set(space)
    return HMCDA(n_iters, n_adapt, delta, lambda, _space, 0)
end
function HMCDA(n_iters::Int, n_adapt::Int, δ::Float64, λ::Float64, space::Set, gid::Int)
    return HMCDA{eltype(space)}(n_iters, n_adapt, δ, λ, space, gid)
end

function step(model, spl::Sampler{<:HMCDA}, vi::VarInfo, is_first::Bool)
    if is_first
        if ~haskey(spl.info, :wum)
            spl.alg.gid != 0 && link!(vi, spl)

            init_warm_up_params(vi, spl)

            θ = vi[spl]
            ϵ = spl.alg.delta > 0 ?
                find_good_eps(model, vi, spl) :       # heuristically find optimal ϵ
                spl.info[:pre_set_ϵ]
            vi[spl] = θ

            spl.alg.gid != 0 && invlink!(vi, spl)
            push!(spl.info[:wum][:ϵ], ϵ)
            update_da_μ(spl.info[:wum], ϵ)
        end

        push!(spl.info[:accept_his], true)
    else
        # Set parameters
        λ = spl.alg.lambda
        ϵ = spl.info[:wum][:ϵ][end]
        @debug "current ϵ: $ϵ"

        spl.info[:lf_num] = 0   # reset current lf num counter

        @debug "X-> R..."
        if spl.alg.gid != 0
            link!(vi, spl)
            runmodel!(model, vi, spl)
        end

        grad_func = gen_grad_func(vi, spl, model)
        lj_func = gen_lj_func(vi, spl, model)
        rev_func = gen_rev_func(vi, spl)
        log_func = gen_log_func(spl)

        θ, lj, stds = vi[spl], vi.logp, spl.info[:wum][:stds]

        θ_new, lj_new, is_accept, τ_valid, α = _hmc_step(
            θ, lj, lj_func, grad_func, ϵ, λ, stds; rev_func=rev_func, log_func=log_func)

        if PROGRESS[] && spl.alg.gid == 0
            stds_str = string(spl.info[:wum][:stds])
            stds_str = length(stds_str) >= 32 ? stds_str[1:30]*"..." : stds_str
            haskey(spl.info, :progress) && ProgressMeter.update!(
                spl.info[:progress],
                spl.info[:progress].counter;
                showvalues = [(:ϵ, ϵ), (:α, α), (:pre_cond, stds_str)],
            )
        end

        @debug "decide whether to accept..."
        if is_accept
            push!(spl.info[:accept_his], true)
            vi[spl] = θ_new
            setlogp!(vi, lj_new)
        else
            push!(spl.info[:accept_his], false)
            vi[spl] = θ
            setlogp!(vi, lj)
        end

        # QUES: why do we need the 2nd condition here (Kai)
        if spl.alg.delta > 0
            # TODO: figure out whether or not the condition below is needed
            # if spl.alg.delta > 0 && τ_valid > 0    # only do adaption for HMCDA
            adapt!(spl.info[:wum], α, vi[spl], adapt_ϵ = true)
        end

        @debug "R -> X..."
        spl.alg.gid != 0 && invlink!(vi, spl)
    end
    return vi
end

function _hmc_step(
    θ::AbstractVector{<:Real},
    lj::Real,
    lj_func,
    grad_func,
    ϵ::Real,
    λ::Real,
    stds::AbstractVector{<:Real};
    rev_func=nothing,
    log_func=nothing,
)

    θ_dim = length(θ)

    @debug "sampling momentums..."
    p = _sample_momentum(θ_dim, stds)

    @debug "recording old values..."
    H = _find_H(θ, p, lj, stds)

    τ = max(1, round(Int, λ / ϵ))
    @debug "leapfrog for $τ steps with step size $ϵ"
    θ_new, p_new, τ_valid = _leapfrog(θ, p, τ, ϵ, grad_func; rev_func=rev_func, log_func=log_func)

    @debug "computing new H..."
    lj_new = lj_func(θ_new)
    H_new = (τ_valid == 0) ? Inf : _find_H(θ_new, p_new, lj_new, stds)

    @debug "deciding wether to accept and computing accept rate α..."
    is_accept, logα = mh_accept(H, H_new)

    if is_accept
        θ = θ_new
        lj = lj_new
    end

    return θ, lj, is_accept, τ_valid, exp(logα)
end

function _hmc_step(
    θ::AbstractVector{<:Real},
    lj::Real,
    lj_func,
    grad_func,
    τ::Int,
    ϵ::Real,
    stdsstds::AbstractVector{<:Real},
)
    return _hmc_step(θ, lj, lj_func, grad_func, ϵ, τ * ϵ, stds)
end