text:
#
# This file is part of the Simulate.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

"a state machine is undefined (after creation)"
struct Undefined <: ClockState end

"a state machine is idle"
struct Idle <: ClockState end

"a state machine is empty"
struct Empty <: ClockState end

"a state machine is busy"
struct Busy <: ClockState end

"a state machine is halted"
struct Halted <: ClockState end

"`Init(info)`: Init event with some info."
struct Init <: ClockEvent
    info::Any
end

"`Setup(vars::Array{Symbol,1}, scope::Module)`: setup a logger with some info."
struct Setup <: ClockEvent
    vars::Array{Symbol,1}
    scope::Module
end

"`Step()`: command"
struct Step <: ClockEvent end

"`Run()`: command"
struct Run <: ClockEvent
    duration::Float64
end

"`Register()`: command from master to an active clock."
struct Register{T <: Union{DiscreteEvent,DiscreteCond,Sample}} <: ClockEvent
    x::T
end

"`Forward()`: forward an event or sample from an active clock to master"
struct Forward{T <: Union{DiscreteEvent,DiscreteCond,Sample}} <: ClockEvent
    x::T
    id::Int
end

"`Sync()`: event for syncing of parallel clocks"
struct Sync <: ClockEvent end

"`Start()`: event used for transferring start data to an active clock"
struct Start <: ClockEvent
    x::Any
end

"`Stop()`: command"
struct Stop <: ClockEvent end

"`Resume()`: command"
struct Resume <: ClockEvent end

"`Clear()`: command"
struct Clear <: ClockEvent end

"`Query()`: command, causes an active clock to send its clock data."
struct Query <: ClockEvent end

"`Diag()`: command, causes an active clock to send the last stacktrace."
struct Diag <: ClockEvent end

"`Reset()`: command"
struct Reset <: ClockEvent end

"`Response()`: response from an active clock"
struct Response <: ClockEvent
    x::Any
end

"`Done()`: event from an active clock containing the elapsed time in ns"
struct Done <: ClockEvent
    t::UInt
end

"""
    step!(A::AbstractClock, q::ClockState, σ::ClockEvent)

Default transition for clock state machines.

This is called if no otherwise defined transition occurs.

# Arguments
- `A::AbstractClock`: state machine for which a transition is called
- `q::ClockState`:  state of the state machine
- `σ::ClockEvent`:  event, triggering the transition
"""
function step!(A::AbstractClock, q::ClockState, σ::ClockEvent)
    println(stderr, "Warning: undefined transition ",
            "$(typeof(A)), ::$(typeof(q)), ::$(typeof(σ)))")
end


text:
@doc """
    simple_estimator(model::Clock, T::Real, Js::AbstractArray)
    simple_estimator(model::XY, T::Real, Js::AbstractArray)

Returns the following observables as `Dict{String, Any}`

# Observables
- `"Energy"`
    - Energy per spin (site)
- `"Energy^2"`
- `"|Magnetization|"`
    - Absolute value of total magnetization per spin (order paremeter)
- `"|Magnetization|^2"`
- `"|Magnetization|^4"`
- `"Magnetization x"`
    - x component of total magnetization per spin (order paremeter)
- `"|Magnetization x|"`
- `"Magnetization x^2"`
- `"Magnetization x^4"`
- `"Magnetization y"`
    - y component of total magnetization per spin (order paremeter)
- `"|Magnetization y|"`
- `"Magnetization y^2"`
- `"Magnetization y^4"`
- `"Helicity Modulus x"`
- `"Helicity Modulus y"`
"""
function simple_estimator(model::Clock, T::Real, Js::AbstractArray, _=nothing)
    nsites = numsites(model)
    nbonds = numbonds(model)
    D = dim(model)
    invN = 1.0/nsites
    beta = 1.0/T

    M = zeros(2)
    E = 0.0
    U1 = zeros(2)
    U2 = zeros(2)

    @inbounds for s in 1:nsites
        M[1] += model.cosines[model.spins[s]]
        M[2] += model.sines[model.spins[s]]
    end

    @inbounds for b in bonds(model)
        i, j = source(b), target(b)
        dir = bonddirection(b)
        dt = mod1(model.spins[j] - model.spins[i], model.Q)
        E -= model.cosines[dt] * Js[bondtype(b)]

        for d in 1:D
            U1[d] += model.cosines[dt] * dir[d]^2
            U2[d] += model.sines[dt] * dir[d]
        end
    end
    for d in 1:2
        M[d] *= invN
        U1[d] -= beta * U2[d]^2
        U1[d] *= invN
    end
    E /= nsites

    res = Measurement()
    for (i,c) in enumerate(["x","y"])
        res["Magnetization $c"] = M[i]
        res["|Magnetization $c|"] = abs(M[i])
        res["Magnetization $c^2"] = M[i]^2
        res["Magnetization $c^4"] = M[i]^4
        res["Helicity Modulus $c"] = U1[i]
    end
    M2 = sum(abs2, M)
    res["|Magnetization|"] = sqrt(M2)
    res["|Magnetization|^2"] = M2
    res["|Magnetization|^4"] = M2^2
    res["Energy"] = E
    res["Energy^2"] = E^2

    return res
end

ext: jl


text:
using Base: Char
# one-body matrix, nearest neighbour, PBC
function kinetic_matrix_generator(NsX::Int64, NsY::Int64, t::Float64)

    KMatrix = zeros(NsX * NsY, NsX * NsY)

    for i = 1 : NsX
        for j = 1 : NsY
            # indices of nearest neighbours (nn) of (i, j)
            nn_up = mod(i - 2, NsX) + 1
            nn_dn = mod(i, NsX) + 1
            nn_lf = mod(j - 2, NsY) + 1
            nn_rg = mod(j, NsY) + 1

            # parse indices into the kinetic matrix
            KMatrix[(i - 1) * NsY + j, (nn_up - 1) * NsY + j] = -t
            KMatrix[(i - 1) * NsY + j, (nn_dn - 1) * NsY + j] = -t
            KMatrix[(i - 1) * NsY + j, (i - 1) * NsY + nn_lf] = -t
            KMatrix[(i - 1) * NsY + j, (i - 1) * NsY + nn_rg] = -t
        end
    end

    return KMatrix

end

# auxiliary field matrix generator
@inline function auxfield_matrix_generator(σ::Array{Int64,1}, γ::Float64, Δτ::Float64, U::Float64)

    AFMatrixUp = diagm(0 => exp.(2 * γ * σ .- Δτ * U / 2))
    AFMatrixDn = diagm(0 => exp.(-2 * γ * σ .- Δτ * U / 2))
    return AFMatrixUp, AFMatrixDn

end

function singlestep_matrix(σ::Array{Int64,1}, lattice::Lattice, qmc::QMC)

    AFMatrixUp, AFMatrixDn = auxfield_matrix_generator(σ, qmc.γ, qmc.Δτ, lattice.U)
    Bl_Up = qmc.KMatrixExp * AFMatrixUp * qmc.KMatrixExp
    Bl_Dn = qmc.KMatrixExp * AFMatrixDn * qmc.KMatrixExp

    return Bl_Up, Bl_Dn

end

# Calculate Q', D', T' for Q'D'T' = QDT * B or B * QDT
# Q => orthogonal matrix
# D => diagonal matrix with descending entries
# arguments Q, D, T would be overwritten
function QRCP_update!(
    Q::Array{Float64,2}, D::Array{Float64,2}, T::Array{Float64,2},
    B::Array{Float64,2}, direction::Char
    )

    # Q'D'T' = B * QDT
    if direction == 'L'
        BQD = B * Q * D
        # column-pivoted QR (QRCP) decomposition
        QRCP_BQD = LinearAlgebra.qr(BQD, Val(true))
        Q .= QRCP_BQD.Q * I
        D .= diagm(0 => diag(QRCP_BQD.R))
        T .= (inv(D) * QRCP_BQD.R) * (QRCP_BQD.P' * T)

        return Q, D, T

    # Q'D'T' = QDT * B
    elseif direction == 'R'
        DTB = D * T * B
        # column-pivoted QR (QRCP) decomposition
        QRCP_DTB = LinearAlgebra.qr(DTB, Val(true))
        Q .= Q * QRCP_DTB.Q
        D .= diagm(0 => diag(QRCP_DTB.R))
        T .= (inv(D) * QRCP_DTB.R) * QRCP_DTB.P'

        return Q, D, T

    else
        ArgumentError("direction can only be 'L' or 'R'")
    end

end

# full propagation using column-pivoted QR (QRCP) decomposition
function matrix_product_stablized(
    σlist::Array{Int64,2}, lattice::Lattice, qmc::QMC
    )

    MatProdUp_Q = Matrix{Float64}(I, lattice.Ns, lattice.Ns)
    MatProdUp_D = Matrix{Float64}(I, lattice.Ns, lattice.Ns)
    MatProdUp_T = Matrix{Float64}(I, lattice.Ns, lattice.Ns)

    MatProdDn_Q = Matrix{Float64}(I, lattice.Ns, lattice.Ns)
    MatProdDn_D = Matrix{Float64}(I, lattice.Ns, lattice.Ns)
    MatProdDn_T = Matrix{Float64}(I, lattice.Ns, lattice.Ns)

    for i = 1 : div(qmc.L, qmc.cluster_size)

        MatProdUp = Matrix{Float64}(I, lattice.Ns, lattice.Ns)
        MatProdDn = Matrix{Float64}(I, lattice.Ns, lattice.Ns)

        for j = 1 : qmc.cluster_size
            σ = σlist[:, (i - 1) * qmc.cluster_size + j]
            #AFMatrixUp, AFMatrixDn = auxfield_matrix_generator(σ, qmc.γ, qmc.Δτ, lattice.U)
            #MatProdUp = Symmetric(qmc.KMatrixExp * AFMatrixUp * qmc.KMatrixExp) * MatProdUp
            #MatProdDn = Symmetric(qmc.KMatrixExp * AFMatrixDn * qmc.KMatrixExp) * MatProdDn
            Bj_Up, Bj_Dn = singlestep_matrix(σ, lattice, qmc)
            MatProdUp = Bj_Up * MatProdUp
            MatProdDn = Bj_Dn * MatProdDn
        end

        QRCP_update!(MatProdUp_Q, MatProdUp_D, MatProdUp_T,
                    MatProdUp, 'L')
        QRCP_update!(MatProdDn_Q, MatProdDn_D, MatProdDn_T,
                    MatProdDn, 'L')

    end

    return MatProdUp_Q, MatProdUp_D, MatProdUp_T, MatProdDn_Q, MatProdDn_D, MatProdDn_T

end

# generate overlap "tensor" / column major
function overlap_coefficient(Ns::Int64,
    UUp::Array{T1,2}, UDn::Array{T2,2}
) where {T1<:FloatType, T2<:FloatType}

    UUp_Inv, UDn_Inv = inv(UUp), inv(UDn)
    UUp_Array = complex(zeros(Ns, Ns, Ns))
    UDn_Array = complex(zeros(Ns, Ns, Ns))

    for i = 1 : Ns
        for j = 1 : Ns
            for k = 1 : Ns
                UUp_Array[k ,i, j] = UUp_Inv[k, i] * UUp[j, k]
                UDn_Array[k, i, j] = UDn_Inv[k, i] * UDn[j, k]
            end
        end
    end

    return UUp_Array, UDn_Array

end

max_stars_repo_name: TongSericus/CanonicalAFQMC

max_issues_repo_name: TongSericus/CanonicalAFQMC

max_forks_repo_name: TongSericus/CanonicalAFQMC


text:
Base.isequal(model1::CoreModel, model2::CoreModel) =
    isequal(model1.S, model2.S) &&
    isequal(model1.b, model2.b) &&
    isequal(model1.c, model2.c) &&
    isequal(model1.xl, model2.xl) &&
    isequal(model1.xu, model2.xu) &&
    isequal(model1.rxns, model2.rxns) &&
    isequal(model1.mets, model2.mets)

Base.copy(model::CoreModel) =
    CoreModel(model.S, model.b, model.c, model.xl, model.xu, model.rxns, model.mets)

Base.isequal(model1::CoreModelCoupled, model2::CoreModelCoupled) =
    isequal(model1.lm, model2.lm) &&
    isequal(model1.C, model2.C) &&
    isequal(model1.cl, model2.cl) &&
    isequal(model1.cu, model2.cu)

Base.copy(model::CoreModelCoupled) = CoreModelCoupled(model.lm, model.C, model.cl, model.cu)


text:
struct MatrixAngularCentralGaussian{N,T,n,k} <: ParameterizedMeasure{N}
    par::NamedTuple{N,T}
    manifold::Stiefel{n,k,ℝ}
end

function MatrixAngularCentralGaussian(n, k; kwargs...)
    par = NamedTuple(kwargs)
    return MatrixAngularCentralGaussian(par, Stiefel(n, k))
end

Manifolds.base_manifold(d::MatrixAngularCentralGaussian) = getfield(d, :manifold)

function MeasureTheory.basemeasure(μ::MatrixAngularCentralGaussian{(:Σ⁻¹,)})
    return normalize(Hausdorff(base_manifold(μ), μ.Σ⁻¹))
end
function MeasureTheory.basemeasure(μ::MatrixAngularCentralGaussian{(:L,)})
    return normalize(Hausdorff(base_manifold(μ), μ.L))
end

# Chikuse, 2003 Eq. 2.4.3
function MeasureTheory.logdensity(
    d::MatrixAngularCentralGaussian{(:Σ⁻¹,)}, x::AbstractMatrix
)
    n, k = representation_size(base_manifold(d))
    Σ⁻¹ = d.Σ⁻¹
    return -n//2 * log(dot(x, Σ⁻¹, x)) + k//2 * logdet(Σ⁻¹)
end
function MeasureTheory.logdensity(d::MatrixAngularCentralGaussian{(:L,)}, x::AbstractMatrix)
    n, k = representation_size(base_manifold(d))
    L = d.L
    z = LowerTriangular(L) \ x
    return -n//2 * log(dot(z, z)) - k * logdet(L)
end

function Base.rand(rng::AbstractRNG, ::Type, d::MatrixAngularCentralGaussian{(:L,)})
    L = d.L
    z = randn!(rng, similar(d.L, representation_size(base_manifold(d))))
    y = lmul!(LowerTriangular(L), z)
    return y * (Symmetric(y'y))^(-1//2)
end


text:
# source: https://github.com/umpirsky/country-list/blob/0007d7ce1736f15ed691d1d35c7c32801a769753/data/en_US/country.csv
const countries = [
    "Afghanistan",
    "Åland Islands",
    "Albania",
    "Algeria",
    "American Samoa",
    "Andorra",
    "Angola",
    "Anguilla",
    "Antarctica",
    "Antigua & Barbuda",
    "Argentina",
    "Armenia",
    "Aruba",
    "Australia",
    "Austria",
    "Azerbaijan",
    "Bahamas",
    "Bahrain",
    "Bangladesh",
    "Barbados",
    "Belarus",
    "Belgium",
    "Belize",
    "Benin",
    "Bermuda",
    "Bhutan",
    "Bolivia",
    "Bosnia & Herzegovina",
    "Botswana",
    "Bouvet Island",
    "Brazil",
    "British Indian Ocean Territory",
    "British Virgin Islands",
    "Brunei",
    "Bulgaria",
    "Burkina Faso",
    "Burundi",
    "Cambodia",
    "Cameroon",
    "Canada",
    "Cape Verde",
    "Caribbean Netherlands",
    "Cayman Islands",
    "Central African Republic",
    "Chad",
    "Chile",
    "China",
    "Christmas Island",
    "Cocos (Keeling) Islands",
    "Colombia",
    "Comoros",
    "Congo - Brazzaville",
    "Congo - Kinshasa",
    "Cook Islands",
    "Costa Rica",
    "Côte d’Ivoire",
    "Croatia",
    "Cuba",
    "Curaçao",
    "Cyprus",
    "Czechia",
    "Denmark",
    "Djibouti",
    "Dominica",
    "Dominican Republic",
    "Ecuador",
    "Egypt",
    "El Salvador",
    "Equatorial Guinea",
    "Eritrea",
    "Estonia",
    "Eswatini",
    "Ethiopia",
    "Falkland Islands",
    "Faroe Islands",
    "Fiji",
    "Finland",
    "France",
    "French Guiana",
    "French Polynesia",
    "French Southern Territories",
    "Gabon",
    "Gambia",
    "Georgia",
    "Germany",
    "Ghana",
    "Gibraltar",
    "Greece",
    "Greenland",
    "Grenada",
    "Guadeloupe",
    "Guam",
    "Guatemala",
    "Guernsey",
    "Guinea",
    "Guinea-Bissau",
    "Guyana",
    "Haiti",
    "Heard & McDonald Islands",
    "Honduras",
    "Hong Kong SAR China",
    "Hungary",
    "Iceland",
    "India",
    "Indonesia",
    "Iran",
    "Iraq",
    "Ireland",
    "Isle of Man",
    "Israel",
    "Italy",
    "Jamaica",
    "Japan",
    "Jersey",
    "Jordan",
    "Kazakhstan",
    "Kenya",
    "Kiribati",
    "Kuwait",
    "Kyrgyzstan",
    "Laos",
    "Latvia",
    "Lebanon",
    "Lesotho",
    "Liberia",
    "Libya",
    "Liechtenstein",
    "Lithuania",
    "Luxembourg",
    "Macao SAR China",
    "Madagascar",
    "Malawi",
    "Malaysia",
    "Maldives",
    "Mali",
    "Malta",
    "Marshall Islands",
    "Martinique",
    "Mauritania",
    "Mauritius",
    "Mayotte",
    "Mexico",
    "Micronesia",
    "Moldova",
    "Monaco",
    "Mongolia",
    "Montenegro",
    "Montserrat",
    "Morocco",
    "Mozambique",
    "Myanmar (Burma)",
    "Namibia",
    "Nauru",
    "Nepal",
    "Netherlands",
    "New Caledonia",
    "New Zealand",
    "Nicaragua",
    "Niger",
    "Nigeria",
    "Niue",
    "Norfolk Island",
    "North Korea",
    "North Macedonia",
    "Northern Mariana Islands",
    "Norway",
    "Oman",
    "Pakistan",
    "Palau",
    "Palestinian Territories",
    "Panama",
    "Papua New Guinea",
    "Paraguay",
    "Peru",
    "Philippines",
    "Pitcairn Islands",
    "Poland",
    "Portugal",
    "Puerto Rico",
    "Qatar",
    "Réunion",
    "Romania",
    "Russia",
    "Rwanda",
    "Samoa",
    "San Marino",
    "São Tomé & Príncipe",
    "Saudi Arabia",
    "Senegal",
    "Serbia",
    "Seychelles",
    "Sierra Leone",
    "Singapore",
    "Sint Maarten",
    "Slovakia",
    "Slovenia",
    "Solomon Islands",
    "Somalia",
    "South Africa",
    "South Georgia & South Sandwich Islands",
    "South Korea",
    "South Sudan",
    "Spain",
    "Sri Lanka",
    "St. Barthélemy",
    "St. Helena",
    "St. Kitts & Nevis",
    "St. Lucia",
    "St. Martin",
    "St. Pierre & Miquelon",
    "St. Vincent & Grenadines",
    "Sudan",
    "Suriname",
    "Svalbard & Jan Mayen",
    "Sweden",
    "Switzerland",
    "Syria",
    "Taiwan",
    "Tajikistan",
    "Tanzania",
    "Thailand",
    "Timor-Leste",
    "Togo",
    "Tokelau",
    "Tonga",
    "Trinidad & Tobago",
    "Tunisia",
    "Turkey",
    "Turkmenistan",
    "Turks & Caicos Islands",
    "Tuvalu",
    "U.S. Outlying Islands",
    "U.S. Virgin Islands",
    "Uganda",
    "Ukraine",
    "United Arab Emirates",
    "United Kingdom",
    "United States",
    "Uruguay",
    "Uzbekistan",
    "Vanuatu",
    "Vatican City",
    "Venezuela",
    "Vietnam",
    "Wallis & Futuna",
    "Western Sahara",
    "Yemen",
    "Zambia",
    "Zimbabwe",
]


text:
## Example script for testing linear-strength vortex and source panel methods
using LinearAlgebra
using Base.Iterators
using Seaborn
using AeroMDAO

## Tests
wedge_pts    = @. Point2D([0.5, 0.0, -0.5, 0.0, 0.5],
                          [-0.0, -0.5, 0.0, 0.5, 0.0])
wedge_panels = @. Panel2D(wedge_pts[1:end-1], wedge_pts[2:end])
uniform      = Uniform2D(1.0, 0.0)

## Source
A = source_influence_matrix(wedge_panels)
b = neumann_boundary_condition(wedge_panels, velocity(uniform))
σs = A \ b

##
sum(σs)

##
x_domain, y_domain = (-1, 1), (-1, 1)
grid_size = 50
x_dom, y_dom = range(x_domain..., length = grid_size), range(y_domain..., length = grid_size)
grid = product(x_dom, y_dom)
wedge = panel_points(wedge_panels)

##
#
source_vels = [ velocity(uniform) .+ sum(source_velocity.(σs[1:end-1], σs[2:end], wedge_panels, x, y)) for (x, y) in grid ]
# source_vels = [ sum(source_velocity.(1., 1., wedge_panels, x, y)) for (x, y) in grid ]
speeds  = @. sqrt(first(source_vels)^2 + last(source_vels)^2)

cock = plot(first.(wedge), last.(wedge))
streamplot(first.(grid)', last.(grid)', first.(source_vels)', last.(source_vels)', color = speeds', density = 2.5)

show()

## Vortex
A = vortex_influence_matrix(wedge_panels)
b = neumann_boundary_condition(wedge_panels, velocity(uniform))
γs = A \ b

##
sum(γs)

##
vortex_vels = [ sum(vortex_velocity(1., 1., panel, x, y) for panel in wedge_panels) for (x, y) in grid ]

vortex_vels = [ velocity(uniform) .+ sum(vortex_velocity.(γs[1:end-1], γs[2:end], wedge_panels, x, y)) for (x, y) in grid ]

speeds  = @. sqrt(first(vortex_vels)^2 .+ last(vortex_vels)^2)

plot(first.(wedge), last.(wedge))
streamplot(first.(grid)', last.(grid)', first.(vortex_vels)', last.(vortex_vels)', color = speeds', density = 2.5)

show()

##
vortex_source_vels = [ sum(total_velocity.(σs[1:end-1], σs[2:end], γs[1:end-1], γs[2:end], wedge_panels, x, y)) for (x, y) in grid ]
speeds     = @. sqrt(first(vortex_source_vels)^2 + last(vortex_source_vels)^2)

plot(first.(wedge), last.(wedge))
quiver(first.(grid), last.(grid), first.(vortex_source_vels), last.(vortex_source_vels), speeds)

show()


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03xp(datfile, io=stdout)
    ZERO = 0.0e0
    ONE = 1.0e0
    NIN = 5
    NOUT = 6
    NMAX = 200
    LDA = NMAX
    LDB = NMAX
    LDQ = NMAX
    LDRES = NMAX
    LDWORK = NMAX
    LDZ = NMAX
    A = Array{Float64,2}(undef, LDA,NMAX)
    B = Array{Float64,2}(undef, LDA,NMAX)
    Q = Array{Float64,2}(undef, LDQ,NMAX)
    Z = Array{Float64,2}(undef, LDZ,NMAX)
    ALPHAR = Array{Float64,1}(undef, NMAX)
    ALPHAI = Array{Float64,1}(undef, NMAX)
    BETA = Array{Float64,1}(undef, NMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    RES = Array{Float64,2}(undef, LDRES,3*NMAX)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    ILO = parse(BlasInt, vs[2])
    IHI = parse(BlasInt, vs[3])
    if ( N<=0 || N>NMAX )
        @error "illegal N=$N"
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 1

    #FOREIGN.dlacpy!( 'All', N, N, A, LDA, RES(1,N+1), LDRES )
    RES[1:N,N+1:2N] .= A[1:N,1:N]

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       B[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    close(f)
# interp call 2

    #FOREIGN.dlacpy!( 'All', N, N, B, LDB, RES(1,2*N+1), LDRES )
    RES[1:N,2*N+1:3N] .= B[1:N,1:N]

# interp call 3

    INFO = SLICOT.mb03xp!('S', 'I', 'I', N, ILO, IHI, A, B, Q, Z, ALPHAR, ALPHAI, BETA, LDWORK)
    @test INFO == 0
    INFO == 0 || return

# interp output 1
    println(io, "A:")
    _nc = N
    _nr = N
    show(io, "text/plain", A[1:_nr,1:_nc])
    println(io)

# interp call 4

    #FOREIGN.dgemm!( 'No Transpose', 'No Transpose', N, N, N, ONE, RES(1,N+1), LDRES, Z, LDZ, ZERO, RES, LDRES )
    BLAS.gemm!('N','N',ONE,RES[1:N,N+1:2N],Z[1:N,1:N],ZERO,view(RES,1:N,1:N))

# interp call 5

    #FOREIGN.dgemm!( 'No Transpose', 'No Transpose', N, N, N, -ONE, Q, LDQ, A, LDA, ONE, RES, LDRES )
    BLAS.gemm!('N','N',-ONE,Q[1:N,1:N],A[1:N,1:N],ONE,view(RES,1:N,1:N))
    resid = norm(RES[1:N,1:N])
    println(io, "Frobenius norm of A Z - Q S: $resid")
    @test resid < 1e-12
    println(io, "B:")
    _nc = N
    _nr = N
    show(io, "text/plain", B[1:_nr,1:_nc])
    println(io)

# interp call 6

    #FOREIGN.dgemm!( 'No Transpose', 'No Transpose', N, N, N, ONE, RES(1,2*N+1), LDRES, Q, LDQ, ZERO, RES, LDRES )
    BLAS.gemm!('N','N',ONE,RES[1:N,2*N+1:3N],Q[1:N,1:N],ZERO,view(RES,1:N,1:N))

# interp call 7

    #FOREIGN.dgemm!( 'No Transpose', 'No Transpose', N, N, N, -ONE, Z, LDZ, B, LDB, ONE, RES, LDRES )
    BLAS.gemm!('N','N',-ONE,Z[1:N,1:N],B[1:N,1:N],ONE,view(RES,1:N,1:N))
    resid = norm(RES[1:N,1:N])
    println(io, "Frobenius norm of B Q - Z T:  $resid")
    @test resid < 1e-12
    println(io, "Q:")
    _nc = N
    _nr = N
    show(io, "text/plain", Q[1:_nr,1:_nc])
    println(io)

# interp call 8

    #FOREIGN.dgemm!( 'Transpose', 'No Transpose', N, N, N, ONE, Q, LDQ, Q, LDQ, ONE, RES, LDRES )
    BLAS.gemm!('T','N',ONE,Q[1:N,1:N],Q[1:N,1:N],ZERO,view(RES,1:N,1:N))
    resid =  norm(RES[1:N,1:N] - I)
    println(io, "orth. residual norm of Q: $resid")# interp output 4
    @test resid < 1e-12
    println(io, "Z:")
    _nc = N
    _nr = N
    show(io, "text/plain", Z[1:_nr,1:_nc])
    println(io)

# interp call 9

    #FOREIGN.dgemm!( 'Transpose', 'No Transpose', N, N, N, ONE, Z, LDZ, Z, LDZ, ONE, RES, LDRES )
    BLAS.gemm!('T','N',ONE,Z[1:N,1:N],Z[1:N,1:N],ZERO,view(RES,1:N,1:N))
    resid = norm(RES[1:N,1:N] - I)
    println(io, "orth. residual norm of Z: $resid")# unable to translate write statement:
    @test resid < 1e-12
#   in do block [('70', 'I', 'N')]
    #   write ALPHAR(I), ALPHAI(I), BETA(I)
    println(io, "ALPHA:")
    show(io, "text/plain", ALPHAR[1:N]+im*ALPHAI[1:N])
    println(io)
    println(io, "BETA:")
    show(io, "text/plain", BETA[1:N])
    println(io)

end # run_mb03xp()


text:
using ArbNumerics, Plots
prec = setprecision(BigFloat, 256)

#a = BigFloat(0.506073569036822351319599371053047956980141736828203749380990114218225638827)
#b = BigFloat(0.02328852830307032054478158044023918735669943648088852646123182739831022528)
#c = BigFloat(0.121205692738975111744666848150620569782497212127938371936404761693002104361)

#function T(x)
#    if x <= (0.3)
#        return (a+cbrt((x-(1/8))))*exp(-x) + b
#    end
#    return c*(10*x*exp((-10*x)/3))^19 + b
#end
##
#epsilonc = BigFloat(0.7)
#tauc = BigFloat(0.8)

#function T(x)
#    t = x-(epsilonc /(2*pi))*sin(2*pi*x)+tauc
#    if t < 1
#        return t
#    else
#        return t-1
#    end
#end

a = BigFloat(3.93)
function T(x)
   t = a*x*(1-x)
    if t < 1
        return t
    else
        return t -1
    end
end

T(0.2)
#typeof(T(0.2))
##

Size = 256
noise_size = 0.05

P, gamma = NoiseFourier.Fourier1D.assemble_matrix(T, Size; x_0 = 0, x_1 = 1, T = BigFloat)
##
gamma

##
D = NoiseFourier.Fourier1D.noise_matrix(noise_size, Size)
M = D*P
typeof(M)
##
maxco = 0
for i in 1:Size
    if D[i,i]> maxco
        maxco = D[i,i]
    end
end

function Gamma(sigma,n)
    setprecision(256)
    return BigFloat((exp((-n^2*sigma^2)/2))/(sigma*n*sqrt(2*pi)))
end

#Approximation Error
gamma1 = gamma * maxco + Size*Gamma(noise_size, 2048)
typeof(gamma1)
##

#Discretization Error
prec = setprecision(BigFloat, 64)
P1, asd = NoiseFourier.Fourier1D.assemble_matrix(T, Size; x_0 = 0, x_1 = 1, T = Float64)
##
D1 = NoiseFourier.Fourier1D.noise_matrix(noise_size, Size)
M1 = D1*P1
v2, vinf = NoiseFourier.Fourier1D.rigorous_norm(M1, k=500)
##
#print(asd)
print(v2)
##

function discr_error(v)
    setprecision(256)
    Ci = BigFloat(0.0)
    t = size(v)[1]
    for i in 1:t-1
        Ci += v[i]
    end
    return (1/(1-v[t])) * Ci * (1 + 1/(noise_size*sqrt(2*pi)) + Gamma(noise_size,Size)) * Gamma(noise_size, Size)
end

##
discr_error(v2)
##

#Power method
prec = setprecision(BigFloat, 256)
v = zeros(2*Size+1)
v[1] = BigFloat(1.0)
for i in 1:1000
    @info(i)
    v = M * v
    v /= norm(v,2)
end
lambda = abs(adjoint(v)*M*v)
##
norm(M*v -v,2)
##
#Eigen error
epsilon = BigFloat(norm(M*v -v,2))
##
eigenerr = epsilon + gamma1/2

##
#Plot
function idx_to_freq(j,k)
    # k is the size of the matrix
    if j>=0 && j <=k ÷ 2+1
        return j-1
    end
    if j >k ÷ 2+1
        return j-k-1
    end
end

function measure(x,v)
    t = real(v[1])
    for i in 2:Size
        a = 2*real(v[i])
        b = 2*imag(v[i])
        k = idx_to_freq(i,Size)
        t += a*cos(2*pi*x*k) - b*sin(2*pi*x*k)
    end
    return t
end

p = 50000
x = range(0,1, length = p)
y = zeros(p)
for i in 1:p
    y[i] = measure(x[i],v)
end
plot(x,y,ylims=(0,5),fmt = :png)#,xticks=25:5:75)

##
#Total error
#err = discr_error(v2) + gamma1 + eigenerr
##
print(gamma1)

##
print(discr_error(v2))
##
print(eigenerr)

max_stars_repo_path: examples/Belousov- Zhabotinsky.jl

max_issues_repo_path: examples/Belousov- Zhabotinsky.jl

max_forks_repo_path: examples/Belousov- Zhabotinsky.jl


text:
#==
!!! note

    This code does not need to be run by normal users.

This script converts the FIRI-2018 xlsx from https://figshare.com/s/357cb03b3e5bed649bbc
to CSV format.

This new file is then put into a .tar.gz and uploaded to my figshare account to be used
as a source URL for artifact. The original data is under a CC BY 4.0 license which gives
freedom to share (copy and redistribute in any medium or format) and adapt as long as
appropriate credit is provided.

Finally, ArtifactUtils.jl is used to generate the Artifact.toml file for FIRITools.jl.
==#

# These packages are not all in FIRITools Project.toml
using Artifacts, ArtifactUtils
using CSV, DataFrames, XLSX
using Tar, CodecZlib

using FIRITools

# User friendly URL is: https://figshare.com/s/357cb03b3e5bed649bbc
firi2018_url = "https://ndownloader.figshare.com/files/11823206?private_link=357cb03b3e5bed649bbc"

firi_xlsx = "firi2018.xlsx"
firi_targz = "firi.tar.gz"
firi_dir = "firi"
isdir(firi_dir) || mkdir(firi_dir)

download("$(firi2018_url)", firi_xlsx)

# Reading the XLSX is very slow
df = DataFrame(XLSX.readtable(firi_xlsx, "FIRI-new_vs_altitude")...)

CSV.write(joinpath(firi_dir, "firi2018.csv"), df)

# Tarball the csv
targz = open(firi_targz, "w")
tar = GzipCompressorStream(targz)
Tar.create(firi_dir, tar)
close(tar)

# Cleanup
rm(firi_xlsx)
rm(firi_dir, recursive=true)

# Now post firi.tar.gz to figshare.
# Friendly URL: https://figshare.com/s/33a146841a3f74a74590
# Then (don't skip posting to figshare, because it downloads the tarball!)...

# Write the Artifacts.toml file
add_artifact!(
    "Artifacts.toml",
    "firi",
    "https://ndownloader.figshare.com/files/26640602?private_link=33a146841a3f74a74590",
    force=true,
)

rm(firi_targz)


text:
requiredTypes = [ "CareerStatus" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( dirname( Base.source_path() ), reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end

export CareerPath

type CareerPath

    # Status table
    Path::Vector{CareerStatus}

    # the activated levels with : academic level being 0
    # academic level    0       Affiliation         3
    # Rank              1       Job                 4
    # Assignment        2       SubJob              5
    ActivatedLevels::Vector{Int}

    Length::Int

    RecruitmentAge::Int

    CareerPath() = new(Vector{CareerStatus}(), Vector{Int}(), 0, 19)
    function CareerPath(careerpath::CareerPath)
        cp = new()
        cp.Path = copy(careerpath.Path)
        cp.ActivatedLevels = copy(careerpath.ActivatedLevels)
        cp.Length = careerpath.Length
        cp.RecruitmentAge = careerpath.RecruitmentAge
        return cp
    end
end

LevelsCodes = Dict{DataType, Int}(AcademicLevel => 0,
                                  Rank => 1,
                                  Assignment => 2,
                                  Affiliation => 3,
                                  Job => 4,
                                  SubJob => 5)

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
using Libfive
using Test

@testset "Libfive.jl" begin
    # All of these trees are manually managed and must be freed
    # In higher-level languages, freeing can be attached to the object's
    # destructor and run during garbage collection.
    x = Libfive.libfive_tree_x()
    y = Libfive.libfive_tree_y()

    x2 = Libfive.libfive_tree_unary( Libfive.libfive_opcode_enum("square"), x)
    y2 = Libfive.libfive_tree_unary( Libfive.libfive_opcode_enum("square"), y)
    s  = Libfive.libfive_tree_binary(Libfive.libfive_opcode_enum("add"), x2, y2)

    one = Libfive.libfive_tree_const(1)
    out = Libfive.libfive_tree_binary(Libfive.libfive_opcode_enum("sub"), s, one)

    # Select a 2D region to export, then write an SVG
    R = Libfive.libfive_region2(Libfive.libfive_interval(-2,2), Libfive.libfive_interval(-2,2))

    Libfive.libfive_tree_save_slice(out, R, 0, 10, "circle.svg");

    # Finally, clean up all of the intermediate trees
    #Libfive.libfive_tree_free(x);
    #Libfive.libfive_tree_free(y);
    #Libfive.libfive_tree_free(x2);
    #Libfive.libfive_tree_free(y2);
    #Libfive.libfive_tree_free(s);
    #Libfive.libfive_tree_free(one);
    #Libfive.libfive_tree_free(out);
end


text:
"""

Solve an expression for any zeros or a system of expressions passed a vector.

Examples:

```
x,y, a,b,c = symbols("x, y, a, b, c", real=true)
solve(x^2 - x + 1)
solve(a*x^2  + b*x + c, x)
solve([x - y - 1, x + y - 2], [x,y])
solve([x - y - 1, x + y - 2])   # same as above, as it will assume all free variables.
```

Can solve a single expression for a lone free variable, or if a variable is passed as a second argument, solve for that variable in terms of the others.

When passed a vector of expressions, `solve` looks for solutions to
the system of m equations with n unknowns, where the equations and
unknowns are passed as `Vector{Sym}`. If the unknowns are not
specified, all the free symbols will be used.

The `SymPy` docs say this about `solve`:

> Note If solve returns [] or raises NotImplementedError, it
> doesn’t mean that the equation has no solutions. It just means
> that it couldn’t find any. Often this means that the solutions
> cannot be represented symbolically. For example, the equation
> `x=cos(x)` has a solution, but it cannot be represented
> symbolically using standard functions.
>
> In fact, solve makes no guarantees whatsoever about the completeness
> of the solutions it finds. Much of solve is heuristics, which may find
> some solutions to an equation or system of equations, but not all of
> them.

The return value depends on the inputs:

* If there is one equation with one specified variable (either explicit, or because `free_symbols` returns only one variable), the return value is an array of solutions.

* Otherwise, if there is a unique solution found a dictionary is returned.

* Otherwise, if there is 0 or more than one solution found, an array of dictionaries is returned.

[Note: this may change. The current arrangement is perhaps more convenient, but having a consistent return type has other advantages.]

Note: The individual components of the array display more nicely than the array.

Reference: [SymPy Docs](http://docs.sympy.org/0.7.5/modules/solvers/solvers.html#algebraic-equations)

"""
function solve(ex::Union{T,Vector{T}};  kwargs...) where {T<:Sym}
    ## No symbols specified? Find them
    xs = free_symbols(ex)
    if length(xs) ==0
        error("The expression has no free variables")
    elseif length(xs) == 1
        xs = xs[1]
    end
    solve(ex, xs; kwargs...)
end

## solve for a single variable, Return Sym[]
function solve(ex::Sym, x::Sym; kwargs...)
    a = sympy_meth(:solve, ex, x;  kwargs...)

    ## Way too much work here to finesse into a nice enough output
    ## (Issue comes from solving without specifying variable when 2 or more variables in the expression
    isa(a, Dict) && return(a)
    length(a) == 0 && return(a)
    if (length(a) > 0 && isa(a[1], Dict))
        d = Dict()
        for kv in a
            for (k,v) in kv
                d[k] = v
            end
        end
        d
    else
        Sym[v for v in a]
    end
end

function _mapdict(out, xs) ## can be a tuple if m=n
    d = Dict{Sym, Sym}()
    [d[xs[i]] = out[i] for i in 1:length(out)]
    d
end
function _mapdict(out::Dict,xs=nothing)
    d = Dict{Sym,Sym}()
    [d[k]=v for (k,v) in out]
    d
end

## Solve for a single variable from equations. Return Dict{Sym,Sym}
function solve(exs::Union{T,Vector{T}}, x::Sym; kwargs...) where {T<:Sym}

    solve(exs, [x;]; kwargs...)
end

function solve(exs::Union{T,Vector{T}}, xs::Vector{S}; kwargs...) where {T<:Sym,S<:Sym}
    a = sympy_meth(:solve, exs, xs;  kwargs...)
    ## nicer output
    if isa(a, Dict)
        _mapdict(a)              # XXX type unstable! should be array... Should we change? XXX
    else
        [_mapdict(_underscore,xs) for _underscore in a]
    end
end

## Numeric solutions

"""
Numerically solve for a zero of an expression.

Examples:
```
solve(x^5 - x -1) # inconclusive, though `N.(solve(x^5 - x - 1))` works
nsolve(x^5 - x - 1, 1)
```

Returns symbolic values. Use `N`, or some other means, to convert to floating point.

Reference: [SymPy Docs](http://docs.sympy.org/0.7.5/modules/solvers/solvers.html#algebraic-equations)
"""
nsolve(ex::Sym, x::Sym, x0::Number, args...; kwargs...) = sympy_meth(:nsolve, ex, x, x0, args...; kwargs...)
nsolve(ex::Sym, x0::Number,  args...; kwargs...) =  sympy_meth(:nsolve, ex, x0,  args...; kwargs...)
function nsolve(ex::Vector{Sym}, x::Vector{Sym}, x0::Vector{T}; kwargs...) where {T <: Number}
    out = sympy_meth(:nsolve, tuple(ex...), tuple(x...), tuple(x0...); kwargs...)
    ## ans is matrix object -- convert
    convert(Array{Sym}, sympy["Matrix"](out))
end
export nsolve

"""
`solveset(ex, sym, ...)`

    SymPy is moving to solveset to replace solve (with 1.0)
    http://docs.sympy.org/latest/modules/solvers/solveset.html

The `solveset` function has the advantage of a consistent return type, a set.

The main interface is `solveset(eqn(s), variable(s); domain=...)`

Example
```
out = solveset(x^5 - x - 1)    # set of `CRootof`
as = elements(out)             # An array of same
N.(as)                         # numeric estimates, as would be `N.(solve(x^5 - x - 1))`
solveset(x^2 + 1, x)           # +/- i
solveset(x^2 + 1, x, domain=S.Reals)  # empty set
```
"""
function solveset(ex::Union{T,Vector{T}};  kwargs...) where {T<:Sym}
    ## No symbols specified? Find them
    xs = free_symbols(ex)
    if length(xs) ==0
        error("The expression has non free variables")
    elseif length(xs) == 1
        xs = xs[1]
    end
    solveset(ex, xs; kwargs...)
end

solveset_sympy_methods = (:solveset,
                          :solveset_real, :solveset_complex,
                          :domain_check)

## solveset does not currently allow simulataneous equations
## Solve for a single variable from equations. Return Dict{Sym,Sym}
# function solveset{T<:Sym}(exs::(Union{T,Vector{T}}), x::Sym; kwargs...)
#     solveset(exs, [x;]; kwargs...)
# end

# function solveset{T<:Sym,S<:Sym}(exs::(Union{T,Vector{T}}), xs::Vector{S}; kwargs...)
#     a = sympy_meth(:solveset, convert(SymMatrix, exs), xs;  kwargs...)
#     ## finesse output
#     a
# end

"""

`linsolve`: solve linear system of equations, return a set

Linsolve can be used to solve:

* An augmented matrix. Input the matrix:

Example:
```
@vars x y a b
M = [2a 1 3; 3b 2 1]   # 2ax + y = 3; 3bx+2y=1
linsolve(M, x, y)
```

The matrix must be of `Matrix{Sym}`, so if used on a numeric problem, the matrix must be coerced, as in
`Sym[1 2 3; 3 2 1]`.

* A system of equations. The system is specified as a tuple.

Example
```
linsolve((2a*x+y-3, 3b*x+2y-1), x, y)
linsolve([2a*x+y-3, 3b*x+2y-1], x, y)
```

To get arrays from the sets, the `elements` function is useful:

```
eqs = (2a*x+y-3, 3b*x+2y-1)
as = linsolve(eqs, x, y) # a FiniteSet
bs = elements(as)        # Array{Any}, but really Array of FiniteSets
cs = elements(bs[1])     # Array of symbolic objects
map(simplify, subs(eqs, x => cs[1], y => cs[2]))  # [0,0]. If

* A system `Ax=b`. A tuple is used to pass in `M` and `b`.

Example
```
M = [2a 1; 3b 2]
B = Sym[3;1]
linsolve((M, B), x, y)
```

[cf.](http://docs.sympy.org/dev/modules/solvers/solveset.html#sympy.solvers.solveset.linsolve)
"""
linsolve(exs::Matrix{T}, args...; kwargs...) where {T<:Sym} = sympy_meth(:linsolve, exs, args...; kwargs...)
linsolve(exs::NTuple{N,T}, args...; kwargs...) where {T<:Sym, N} = sympy_meth(:linsolve, exs, args...; kwargs...)
linsolve(exs::Vector{T}, args...; kwargs...) where {T<:Sym} = sympy_meth(:linsolve, tuple(exs...), args...; kwargs...)
linsolve(exs::Tuple{Array{T,M}, N}, args...; kwargs...) where {T<:Sym, M, N} = sympy_meth(:linsolve, exs, args...; kwargs...)

## may not be there
"""
`nonlinsolve`: solve non-linear system of equations, return a set (may not be defined)

[cf.](http://docs.sympy.org/dev/modules/solvers/solvers.html)
"""
nonlinsolve(exs::Vector{T}, args...; kwargs...) where {T<:Sym} = sympy_meth(:nonlinsolve, exs, args...; kwargs...)

export linsolve, nonlinsolve


text:
import ROCAnalysis

mid_of_vec(v::AbstractVector) = v[max(1, round(Int, length(v)/2))]

@recipe function plot(roc::ROCAnalysis.Roc, auc::Bool = true)
    :legend     --> false
    :xlabel     --> "False Positive Rate"
    :ylabel     --> "True Positive Rate"
    :linestyle  --> [:solid :dash]
    :linealpha  --> [1.0 0.5]
    :label      --> ["ROC" ""]

    x = roc.pfa
    y = 1 .- roc.pmiss

    # add fill to the diagonal
    :fillrange  --> Any[0 nothing]
    :fillalpha  --> 0.3

    if auc
        auc_ann = @sprintf("AUC %1.3f", 1 - ROCAnalysis.auc(roc))
        :annotation --> (mid_of_vec(x), mid_of_vec(y), text(auc_ann, :left, :top, 9, :monospace))
    end

    # xargs, yargs
    Any[x, 0:1], Any[y, 0:1]
end


text:
module VEDA
    include("api.jl")

    #using .API

    include("enum.jl")
    include("error.jl")
    include("devices.jl")

    export VEModule, VEFunction, VEContext

    # Hardcoded linker path: from binutils-ve
    const nld_path = "/opt/nec/ve/bin/nld"

    # TODO: Find these at build time.

    mutable struct VEContext
        handle::VEDAcontext

        function VEContext(dev)
            r_handle = Ref{VEDAcontext}()
            @check vedaDevicePrimaryCtxRetain(r_handle, dev)
            ctx = new(r_handle[])
            finalizer(ctx) do ctx
                vedaDevicePrimaryCtxRelease(ctx.handle)
            end
        end
    end

    mutable struct VEModule
        handle::VEDAmodule

        function VEModule(vso::String)
            r_handle = Ref{VEDAmodule}()
            @check vedaModuleLoad(r_handle, vso)
            mod = new(r_handle[])
            finalizer(mod) do mod
                vedaModuleUnload(mod.handle)
            end
            return mod
        end
    end

    function VEModule(obj::Base.CodeUnits{UInt8, String})
        vso = mktemp() do path_o, io_o
            write(io_o, obj)
            flush(io_o)
            vso = path_o*".vso"
            run(`$nld_path -shared -o $vso $path_o`)
            vso
        end
        mod = VEModule(vso)
        rm(vso)
        mod
    end

    mutable struct VEFunction
        handle::VEDAfunction
        mod::VEModule

        function VEFunction(mod, fname)
            r_handle = Ref{VEDAfunction}()
            @check vedaModuleGetFunction(r_handle, mod.handle, fname)
            new(r_handle[], mod)
        end
    end

    include("execution.jl")

    const pctx = Ref{VEContext}()

    ## VEDA device library
    const libcache = Base.WeakKeyDict{VEContext, VEModule}()

    function __init__()
        isempty(libveda) && return
        # limiting to one core for now, as we're not using OpenMP
        ENV["VE_OMP_NUM_THREADS"] = 1
        # TODO: Do a lazy init
        vedaInit(0)
        ctx = VEContext(0)
        pctx[] = ctx
        vedaCtxSetCurrent(ctx.handle)
        #atexit(vedaExit)
    end
end # module


text:
include(joinpath(pwd(), "src/models/box.jl"))
include(joinpath(pwd(), "src/constraints/contact.jl"))

# Horizon
T = 26

# Time step
tf = 2.5
h = tf / (T-1)

# Bounds

# ul <= u <= uu
_uu = Inf * ones(model.m)
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
mrp_init = MRP(UnitQuaternion(RotY(0.0) * RotX(0.0)))
mrp_corner = MRP(UnitQuaternion(RotY(-1.0 * atan(1.0 / sqrt(2.0))) * RotX(pi / 4.0)))

q1 = [model.r, model.r, model.r, mrp_init.x, mrp_init.y, mrp_init.z]
qT = [0.0, 0.0, model.r * sqrt(3.0), mrp_corner.x, mrp_corner.y, mrp_corner.z]

x1 = [q1; q1]

xl, xu = state_bounds(model, T, x1 = x1)

# Objective
obj_penalty = PenaltyObjective(100.0, model.m)

Qq = Diagonal(ones(model.nq))
Q = cat(0.5 * Qq, 0.5 * Qq, dims = (1, 2))
QT = cat(0.5 * Qq, 100.0 * Diagonal(ones(model.nq)), dims = (1, 2))
R = Diagonal([1.0e-1 * ones(model.nu)..., zeros(model.m - model.nu)...])

obj_tracking = quadratic_tracking_objective(
    [t < T ? Q : QT for t = 1:T],
    [R for t = 1:T-1],
    [[zeros(model.nq); qT] for t = 1:T],
    [zeros(model.m) for t = 1:T]
    )

obj = MultiObjective([obj_tracking, obj_penalty])

# Constraints
con_contact = contact_constraints(model, T)

# Problem
prob = problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con_contact
               )

# Trajectory initialization
X0 = linear_interp(x1, [qT; qT], T) # linear interpolation on state
U0 = [0.001 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
Z0 = pack(X0, U0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem

@time Z̄ = solve(prob, copy(Z0), tol = 1.0e-3, c_tol = 1.0e-3)

check_slack(Z̄, prob)

X̄, Ū = unpack(Z̄, prob)

# Visualize
include(joinpath(pwd(), "src/models/visualize.jl"))
vis = Visualizer()
open(vis)
visualize!(vis, model, state_to_configuration(X̄), Δt = h)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using Constraints
using CompositionalNetworks
using Test

@testset "Constraints.jl" begin

# Empty constraint
c = Constraint()
@test concept(c, []) == true
@test error_f(c, []) == 0.0

# Test on usual_constraints
include("usual_constraints.jl")

end


text:
"""
    _powershell() -> AbstractMatrix{<:Colorant}

Paste an image from clipboard using `powershell`
"""
function _powershell()
    mktempdir() do dir
        # Define path
        path_png = joinpath(dir, "clipboard.png")

        # Compose command & run
        addtype = `Add-Type -AssemblyName System.Windows.Forms\;`
        getimg = `\$img=\[Windows.Forms.Clipboard\]::GetImage\(\)\;`
        saveimg = `if \(\$img -ne \$null\)\{\$img.Save\(\"$(path_png)\"\)\}`
        cmd = `powershell $addtype $getimg $saveimg`
        run(cmd)

        # Paste from clipboard
        if isfile(path_png)
            img = load(path_png)
            return img
        else
            error("No image in clipboard")
        end
    end
end

"""
    _powershell(img::AbstractMatrix{<:Colorant})

Copy an image to clipboard using `powershell`
"""
function _powershell(img::AbstractMatrix{<:Colorant})
    mktempdir() do dir
        # Define path
        path_png = joinpath(dir, "clipboard.png")

        # Save image
        save(path_png, img)

        # Compose command & run
        addtype = `Add-Type -AssemblyName System.Windows.Forms\;`
        adddrawing = `\[Reflection.Assembly\]::LoadWithPartialName\(\'System.Drawing\'\)\;`
        getfile = `\$file = get-item\(\"$(path_png)\"\)\;`
        getimg = `\$img = \[System.Drawing.Image\]::Fromfile\(\$file\)\;`
        copyimg = `\[System.Windows.Forms.Clipboard\]::SetImage\(\$img\)\;`
        cmd = `powershell $addtype $adddrawing $getfile $getimg $copyimg`
        run(cmd)
    end
end


text:
export BrickletNFCReaderGetTagIDLowLevel
struct BrickletNFCReaderGetTagIDLowLevel
    tag_type::Integer
    tag_id_length::Integer
    tag_id_data::Vector{Integer}
end

export BrickletNFCReaderGetState
struct BrickletNFCReaderGetState
    state::Integer
    idle::Bool
end

export BrickletNFCReaderReadNDEFLowLevel
struct BrickletNFCReaderReadNDEFLowLevel
    ndef_length::Integer
    ndef_chunk_offset::Integer
    ndef_chunk_data::Vector{Integer}
end

export BrickletNFCReaderReadPageLowLevel
struct BrickletNFCReaderReadPageLowLevel
    data_length::Integer
    data_chunk_offset::Integer
    data_chunk_data::Vector{Integer}
end

export BrickletNFCCardemuGetState
struct BrickletNFCCardemuGetState
    state::Integer
    idle::Bool
end

export BrickletNFCP2PGetState
struct BrickletNFCP2PGetState
    state::Integer
    idle::Bool
end

export BrickletNFCP2PReadNDEFLowLevel
struct BrickletNFCP2PReadNDEFLowLevel
    ndef_length::Integer
    ndef_chunk_offset::Integer
    ndef_chunk_data::Vector{Integer}
end

export BrickletNFCSimpleGetTagIDLowLevel
struct BrickletNFCSimpleGetTagIDLowLevel
    tag_type::Integer
    tag_id_length::Integer
    tag_id_data::Vector{Integer}
    last_seen::Integer
end

export BrickletNFCSPITFPErrorCount
struct BrickletNFCSPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletNFCIdentity
struct BrickletNFCIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletNFCReaderGetTagID
struct BrickletNFCReaderGetTagID
    tag_type::Integer
    tag_id::Vector{Integer}
end

export BrickletNFCSimpleGetTagID
struct BrickletNFCSimpleGetTagID
    tag_type::Integer
    tag_id::Vector{Integer}
    last_seen::Integer
end

export BrickletNFC
"""
NFC tag read/write, NFC P2P and Card Emulation
"""
mutable struct BrickletNFC <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletNFC(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_nfc")
        deviceInternal = package.BrickletNFC(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_mode
"""
    $(SIGNATURES)

Sets the mode. The NFC Bricklet supports four modes:

* Off
* Card Emulation (Cardemu): Emulates a tag for other readers
* Peer to Peer (P2P): Exchange data with other readers
* Reader: Reads and writes tags

If you change a mode, the Bricklet will reconfigure the hardware for this mode.
Therefore, you can only use functions corresponding to the current mode. For
example, in Reader mode you can only use Reader functions.
"""
function set_mode(device::BrickletNFC, mode)
    device.deviceInternal.set_mode(mode)
end

export get_mode
"""
    $(SIGNATURES)

Returns the mode as set by :func:`Set Mode`.
"""
function get_mode(device::BrickletNFC)
    return device.deviceInternal.get_mode()
end

export reader_request_tag_id
"""
    $(SIGNATURES)

After you call :func:`Reader Request Tag ID` the NFC Bricklet will try to read
the tag ID from the tag. After this process is done the state will change.
You can either register the :cb:`Reader State Changed` callback or you can poll
:func:`Reader Get State` to find out about the state change.

If the state changes to *ReaderRequestTagIDError* it means that either there was
no tag present or that the tag has an incompatible type. If the state
changes to *ReaderRequestTagIDReady* it means that a compatible tag was found
and that the tag ID has been saved. You can now read out the tag ID by
calling :func:`Reader Get Tag ID`.

If two tags are in the proximity of the NFC Bricklet, this
function will cycle through the tags. To select a specific tag you have
to call :func:`Reader Request Tag ID` until the correct tag ID is found.

In case of any *ReaderError* state the selection is lost and you have to
start again by calling :func:`Reader Request Tag ID`.
"""
function reader_request_tag_id(device::BrickletNFC)
    device.deviceInternal.reader_request_tag_id()
end

export reader_get_tag_id_low_level
"""
    $(SIGNATURES)

Returns the tag type and the tag ID. This function can only be called if the
NFC Bricklet is currently in one of the *ReaderReady* states. The returned tag ID
is the tag ID that was saved through the last call of :func:`Reader Request Tag ID`.

To get the tag ID of a tag the approach is as follows:

1. Call :func:`Reader Request Tag ID`
2. Wait for state to change to *ReaderRequestTagIDReady* (see :func:`Reader Get State` or
   :cb:`Reader State Changed` callback)
3. Call :func:`Reader Get Tag ID`
"""
function reader_get_tag_id_low_level(device::BrickletNFC)
    return device.deviceInternal.reader_get_tag_id_low_level()
end

export reader_get_state
"""
    $(SIGNATURES)

Returns the current reader state of the NFC Bricklet.

On startup the Bricklet will be in the *ReaderInitialization* state. The
initialization will only take about 20ms. After that it changes to *ReaderIdle*.

The Bricklet is also reinitialized if the mode is changed, see :func:`Set Mode`.

The functions of this Bricklet can be called in the *ReaderIdle* state and all of
the *ReaderReady* and *ReaderError* states.

Example: If you call :func:`Reader Request Page`, the state will change to
*ReaderRequestPage* until the reading of the page is finished. Then it will change
to either *ReaderRequestPageReady* if it worked or to *ReaderRequestPageError* if it
didn't. If the request worked you can get the page by calling :func:`Reader Read Page`.

The same approach is used analogously for the other API functions.
"""
function reader_get_state(device::BrickletNFC)
    return device.deviceInternal.reader_get_state()
end

export reader_write_ndef_low_level
"""
    $(SIGNATURES)

Writes NDEF formated data.

This function currently supports NFC Forum Type 2 and 4.

The general approach for writing a NDEF message is as follows:

1. Call :func:`Reader Request Tag ID`
2. Wait for state to change to *ReaderRequestTagIDReady* (see
   :func:`Reader Get State` or :cb:`Reader State Changed` callback)
3. If looking for a specific tag then call :func:`Reader Get Tag ID` and check
   if the expected tag was found, if it was not found got back to step 1
4. Call :func:`Reader Write NDEF` with the NDEF message that you want to write
5. Wait for state to change to *ReaderWriteNDEFReady* (see :func:`Reader Get State`
   or :cb:`Reader State Changed` callback)
"""
function reader_write_ndef_low_level(device::BrickletNFC, ndef_length, ndef_chunk_offset, ndef_chunk_data)
    device.deviceInternal.reader_write_ndef_low_level(ndef_length, ndef_chunk_offset, ndef_chunk_data)
end

export reader_request_ndef
"""
    $(SIGNATURES)

Reads NDEF formated data from a tag.

This function currently supports NFC Forum Type 1, 2, 3 and 4.

The general approach for reading a NDEF message is as follows:

1. Call :func:`Reader Request Tag ID`
2. Wait for state to change to *RequestTagIDReady* (see :func:`Reader Get State`
   or :cb:`Reader State Changed` callback)
3. If looking for a specific tag then call :func:`Reader Get Tag ID` and check if the
   expected tag was found, if it was not found got back to step 1
4. Call :func:`Reader Request NDEF`
5. Wait for state to change to *ReaderRequestNDEFReady* (see :func:`Reader Get State`
   or :cb:`Reader State Changed` callback)
6. Call :func:`Reader Read NDEF` to retrieve the NDEF message from the buffer
"""
function reader_request_ndef(device::BrickletNFC)
    device.deviceInternal.reader_request_ndef()
end

export reader_read_ndef_low_level
"""
    $(SIGNATURES)

Returns the NDEF data from an internal buffer. To fill the buffer
with a NDEF message you have to call :func:`Reader Request NDEF` beforehand.
"""
function reader_read_ndef_low_level(device::BrickletNFC)
    return device.deviceInternal.reader_read_ndef_low_level()
end

export reader_authenticate_mifare_classic_page
"""
    $(SIGNATURES)

Mifare Classic tags use authentication. If you want to read from or write to
a Mifare Classic page you have to authenticate it beforehand.
Each page can be authenticated with two keys: A (``key_number`` = 0) and B
(``key_number`` = 1). A new Mifare Classic
tag that has not yet been written to can be accessed with key A
and the default key ``[0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF]``.

The approach to read or write a Mifare Classic page is as follows:

1. Call :func:`Reader Request Tag ID`
2. Wait for state to change to *ReaderRequestTagIDReady* (see :func:`Reader Get State`
   or :cb:`Reader State Changed` callback)
3. If looking for a specific tag then call :func:`Reader Get Tag ID` and check if the
   expected tag was found, if it was not found got back to step 1
4. Call :func:`Reader Authenticate Mifare Classic Page` with page and key for the page
5. Wait for state to change to *ReaderAuthenticatingMifareClassicPageReady* (see
   :func:`Reader Get State` or :cb:`Reader State Changed` callback)
6. Call :func:`Reader Request Page` or :func:`Reader Write Page` to read/write page

The authentication will always work for one whole sector (4 pages).
"""
function reader_authenticate_mifare_classic_page(device::BrickletNFC, page, key_number, key)
    device.deviceInternal.reader_authenticate_mifare_classic_page(page, key_number, key)
end

export reader_write_page_low_level
"""
    $(SIGNATURES)

Writes a maximum of 8192 bytes starting from the given page. How many pages are written
depends on the tag type. The page sizes are as follows:

* Mifare Classic page size: 16 byte
* NFC Forum Type 1 page size: 8 byte
* NFC Forum Type 2 page size: 4 byte
* NFC Forum Type 3 page size: 16 byte
* NFC Forum Type 4: No pages, page = file selection (CC or NDEF, see below)

The general approach for writing to a tag is as follows:

1. Call :func:`Reader Request Tag ID`
2. Wait for state to change to *ReaderRequestTagIDReady* (see :func:`Reader Get State` or
   :cb:`Reader State Changed` callback)
3. If looking for a specific tag then call :func:`Reader Get Tag ID` and check if the
   expected tag was found, if it was not found got back to step 1
4. Call :func:`Reader Write Page` with page number and data
5. Wait for state to change to *ReaderWritePageReady* (see :func:`Reader Get State` or
   :cb:`Reader State Changed` callback)

If you use a Mifare Classic tag you have to authenticate a page before you
can write to it. See :func:`Reader Authenticate Mifare Classic Page`.

NFC Forum Type 4 tags are not organized into pages but different files. We currently
support two files: Capability Container file (CC) and NDEF file.

Choose CC by setting page to 3 or NDEF by setting page to 4.
"""
function reader_write_page_low_level(device::BrickletNFC, page, data_length, data_chunk_offset, data_chunk_data)
    device.deviceInternal.reader_write_page_low_level(page, data_length, data_chunk_offset, data_chunk_data)
end

export reader_request_page
"""
    $(SIGNATURES)

Reads a maximum of 8192 bytes starting from the given page and stores them into a buffer.
The buffer can then be read out with :func:`Reader Read Page`.
How many pages are read depends on the tag type. The page sizes are
as follows:

* Mifare Classic page size: 16 byte
* NFC Forum Type 1 page size: 8 byte
* NFC Forum Type 2 page size: 4 byte
* NFC Forum Type 3 page size: 16 byte
* NFC Forum Type 4: No pages, page = file selection (CC or NDEF, see below)

The general approach for reading a tag is as follows:

1. Call :func:`Reader Request Tag ID`
2. Wait for state to change to *RequestTagIDReady* (see :func:`Reader Get State`
   or :cb:`Reader State Changed` callback)
3. If looking for a specific tag then call :func:`Reader Get Tag ID` and check if the
   expected tag was found, if it was not found got back to step 1
4. Call :func:`Reader Request Page` with page number
5. Wait for state to change to *ReaderRequestPageReady* (see :func:`Reader Get State`
   or :cb:`Reader State Changed` callback)
6. Call :func:`Reader Read Page` to retrieve the page from the buffer

If you use a Mifare Classic tag you have to authenticate a page before you
can read it. See :func:`Reader Authenticate Mifare Classic Page`.

NFC Forum Type 4 tags are not organized into pages but different files. We currently
support two files: Capability Container file (CC) and NDEF file.

Choose CC by setting page to 3 or NDEF by setting page to 4.
"""
function reader_request_page(device::BrickletNFC, page, length)
    device.deviceInternal.reader_request_page(page, length)
end

export reader_read_page_low_level
"""
    $(SIGNATURES)

Returns the page data from an internal buffer. To fill the buffer
with specific pages you have to call :func:`Reader Request Page` beforehand.
"""
function reader_read_page_low_level(device::BrickletNFC)
    return device.deviceInternal.reader_read_page_low_level()
end

export cardemu_get_state
"""
    $(SIGNATURES)

Returns the current cardemu state of the NFC Bricklet.

On startup the Bricklet will be in the *CardemuInitialization* state. The
initialization will only take about 20ms. After that it changes to *CardemuIdle*.

The Bricklet is also reinitialized if the mode is changed, see :func:`Set Mode`.

The functions of this Bricklet can be called in the *CardemuIdle* state and all of
the *CardemuReady* and *CardemuError* states.

Example: If you call :func:`Cardemu Start Discovery`, the state will change to
*CardemuDiscover* until the discovery is finished. Then it will change
to either *CardemuDiscoverReady* if it worked or to *CardemuDiscoverError* if it
didn't.

The same approach is used analogously for the other API functions.
"""
function cardemu_get_state(device::BrickletNFC)
    return device.deviceInternal.cardemu_get_state()
end

export cardemu_start_discovery
"""
    $(SIGNATURES)

Starts the discovery process. If you call this function while a NFC
reader device is near to the NFC Bricklet the state will change from
*CardemuDiscovery* to *CardemuDiscoveryReady*.

If no NFC reader device can be found or if there is an error during
discovery the cardemu state will change to *CardemuDiscoveryError*. In this case you
have to restart the discovery process.

If the cardemu state changes to *CardemuDiscoveryReady* you can start the NDEF message
transfer with :func:`Cardemu Write NDEF` and :func:`Cardemu Start Transfer`.
"""
function cardemu_start_discovery(device::BrickletNFC)
    device.deviceInternal.cardemu_start_discovery()
end

export cardemu_write_ndef_low_level
"""
    $(SIGNATURES)

Writes the NDEF message that is to be transferred to the NFC peer.

The maximum supported NDEF message size in Cardemu mode is 255 byte.

You can call this function at any time in Cardemu mode. The internal buffer
will not be overwritten until you call this function again or change the
mode.
"""
function cardemu_write_ndef_low_level(device::BrickletNFC, ndef_length, ndef_chunk_offset, ndef_chunk_data)
    device.deviceInternal.cardemu_write_ndef_low_level(ndef_length, ndef_chunk_offset, ndef_chunk_data)
end

export cardemu_start_transfer
"""
    $(SIGNATURES)

You can start the transfer of a NDEF message if the cardemu state is *CardemuDiscoveryReady*.

Before you call this function to start a write transfer, the NDEF message that
is to be transferred has to be written via :func:`Cardemu Write NDEF` first.

After you call this function the state will change to *CardemuTransferNDEF*. It will
change to *CardemuTransferNDEFReady* if the transfer was successful or
*CardemuTransferNDEFError* if it wasn't.
"""
function cardemu_start_transfer(device::BrickletNFC, transfer)
    device.deviceInternal.cardemu_start_transfer(transfer)
end

export p2p_get_state
"""
    $(SIGNATURES)

Returns the current P2P state of the NFC Bricklet.

On startup the Bricklet will be in the *P2PInitialization* state. The
initialization will only take about 20ms. After that it changes to *P2PIdle*.

The Bricklet is also reinitialized if the mode is changed, see :func:`Set Mode`.

The functions of this Bricklet can be called in the *P2PIdle* state and all of
the *P2PReady* and *P2PError* states.

Example: If you call :func:`P2P Start Discovery`, the state will change to
*P2PDiscover* until the discovery is finished. Then it will change
to either P2PDiscoverReady* if it worked or to *P2PDiscoverError* if it
didn't.

The same approach is used analogously for the other API functions.
"""
function p2p_get_state(device::BrickletNFC)
    return device.deviceInternal.p2p_get_state()
end

export p2p_start_discovery
"""
    $(SIGNATURES)

Starts the discovery process. If you call this function while another NFC
P2P enabled device is near to the NFC Bricklet the state will change from
*P2PDiscovery* to *P2PDiscoveryReady*.

If no NFC P2P enabled device can be found or if there is an error during
discovery the P2P state will change to *P2PDiscoveryError*. In this case you
have to restart the discovery process.

If the P2P state changes to *P2PDiscoveryReady* you can start the NDEF message
transfer with :func:`P2P Start Transfer`.
"""
function p2p_start_discovery(device::BrickletNFC)
    device.deviceInternal.p2p_start_discovery()
end

export p2p_write_ndef_low_level
"""
    $(SIGNATURES)

Writes the NDEF message that is to be transferred to the NFC peer.

The maximum supported NDEF message size for P2P transfer is 255 byte.

You can call this function at any time in P2P mode. The internal buffer
will not be overwritten until you call this function again, change the
mode or use P2P to read an NDEF messages.
"""
function p2p_write_ndef_low_level(device::BrickletNFC, ndef_length, ndef_chunk_offset, ndef_chunk_data)
    device.deviceInternal.p2p_write_ndef_low_level(ndef_length, ndef_chunk_offset, ndef_chunk_data)
end

export p2p_start_transfer
"""
    $(SIGNATURES)

You can start the transfer of a NDEF message if the P2P state is *P2PDiscoveryReady*.

Before you call this function to start a write transfer, the NDEF message that
is to be transferred has to be written via :func:`P2P Write NDEF` first.

After you call this function the P2P state will change to *P2PTransferNDEF*. It will
change to *P2PTransferNDEFReady* if the transfer was successfull or
*P2PTransferNDEFError* if it wasn't.

If you started a write transfer you are now done. If you started a read transfer
you can now use :func:`P2P Read NDEF` to read the NDEF message that was written
by the NFC peer.
"""
function p2p_start_transfer(device::BrickletNFC, transfer)
    device.deviceInternal.p2p_start_transfer(transfer)
end

export p2p_read_ndef_low_level
"""
    $(SIGNATURES)

Returns the NDEF message that was written by a NFC peer in NFC P2P mode.

The NDEF message is ready if you called :func:`P2P Start Transfer` with a
read transfer and the P2P state changed to *P2PTransferNDEFReady*.
"""
function p2p_read_ndef_low_level(device::BrickletNFC)
    return device.deviceInternal.p2p_read_ndef_low_level()
end

export set_detection_led_config
"""
    $(SIGNATURES)

Sets the detection LED configuration. By default the LED shows
if a card/reader is detected.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is off.
"""
function set_detection_led_config(device::BrickletNFC, config)
    device.deviceInternal.set_detection_led_config(config)
end

export get_detection_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Detection LED Config`
"""
function get_detection_led_config(device::BrickletNFC)
    return device.deviceInternal.get_detection_led_config()
end

export set_maximum_timeout
"""
    $(SIGNATURES)

Sets the maximum timeout.

This is a global maximum used for all internal state timeouts. The timeouts depend heavily
on the used tags etc. For example: If you use a Type 2 tag and you want to detect if
it is present, you have to use :func:`Reader Request Tag ID` and wait for the state
to change to either the error state or the ready state.

With the default configuration this takes 2-3 seconds. By setting the maximum timeout to
100ms you can reduce this time to ~150-200ms. For Type 2 this would also still work
with a 20ms timeout (a Type 2 tag answers usually within 10ms). A type 4 tag can take
up to 500ms in our tests.

If you need a fast response time to discover if a tag is present or not you can find
a good timeout value by trial and error for your specific tag.

By default we use a very conservative timeout, to be sure that any tag can always
answer in time.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function set_maximum_timeout(device::BrickletNFC, timeout)
    device.deviceInternal.set_maximum_timeout(timeout)
end

export get_maximum_timeout
"""
    $(SIGNATURES)

Returns the timeout as set by :func:`Set Maximum Timeout`

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function get_maximum_timeout(device::BrickletNFC)
    return device.deviceInternal.get_maximum_timeout()
end

export simple_get_tag_id_low_level
"""
    $(SIGNATURES)

.. versionadded:: 2.0.6\$nbsp;(Plugin)
"""
function simple_get_tag_id_low_level(device::BrickletNFC, index)
    return device.deviceInternal.simple_get_tag_id_low_level(index)
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletNFC)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletNFC, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletNFC)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletNFC, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletNFC, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletNFC, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletNFC)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletNFC)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletNFC)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletNFC, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletNFC)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletNFC)
    return device.deviceInternal.get_identity()
end

export reader_get_tag_id
"""
Returns the tag type and the tag ID. This function can only be called if the
NFC Bricklet is currently in one of the *ReaderReady* states. The returned tag ID
is the tag ID that was saved through the last call of :func:`Reader Request Tag ID`.

To get the tag ID of a tag the approach is as follows:

1. Call :func:`Reader Request Tag ID`
2. Wait for state to change to *ReaderRequestTagIDReady* (see :func:`Reader Get State` or
   :cb:`Reader State Changed` callback)
3. Call :func:`Reader Get Tag ID`
"""
function reader_get_tag_id(device::BrickletNFC)

    ret = reader_get_tag_id_low_level(device, )

        return ReaderGetTagID(ret.tag_type, ret.tag_id_data[:ret.tag_id_length])
end

export reader_write_ndef
"""
Writes NDEF formated data.

This function currently supports NFC Forum Type 2 and 4.

The general approach for writing a NDEF message is as follows:

1. Call :func:`Reader Request Tag ID`
2. Wait for state to change to *ReaderRequestTagIDReady* (see
   :func:`Reader Get State` or :cb:`Reader State Changed` callback)
3. If looking for a specific tag then call :func:`Reader Get Tag ID` and check
   if the expected tag was found, if it was not found got back to step 1
4. Call :func:`Reader Write NDEF` with the NDEF message that you want to write
5. Wait for state to change to *ReaderWriteNDEFReady* (see :func:`Reader Get State`
   or :cb:`Reader State Changed` callback)
"""
function reader_write_ndef(device::BrickletNFC, ndef)

        ndef = convert(Vector{UInt8}, ndef)

    if length(ndef) > 65535
        throw(TinkerforgeInvalidParameterError("NDEF can be at most 65535 items long"))
    end

    ndef_length = length(ndef)
    ndef_chunk_offset = 0

    if ndef_length == 0
        ndef_chunk_data = [0] * 60
        ret = reader_write_ndef_low_level(device, ndef_length, ndef_chunk_offset, ndef_chunk_data)
    else
        lock(device.stream_lock) do
            while ndef_chunk_offset < ndef_length
                ndef_chunk_data = create_chunk_data(ndef, ndef_chunk_offset, 60, 0)
                ret = reader_write_ndef_low_level(device, ndef_length, ndef_chunk_offset, ndef_chunk_data)
                ndef_chunk_offset += 60
            end
        end
    end

        return ret
end

export reader_read_ndef
"""
Returns the NDEF data from an internal buffer. To fill the buffer
with a NDEF message you have to call :func:`Reader Request NDEF` beforehand.
"""
function reader_read_ndef(device::BrickletNFC)

    lock(device.stream_lock) do
        ret = reader_read_ndef_low_level(device, )
            ndef_length = ret.ndef_length
        ndef_out_of_sync = ret.ndef_chunk_offset != 0

        ndef_data = ret.ndef_chunk_data

        while !ndef_out_of_sync && length(ndef_data) < ndef_length
            ret = reader_read_ndef_low_level(device, )
                ndef_length = ret.ndef_length
            ndef_out_of_sync = ret.ndef_chunk_offset != length(ndef_data)
            ndef_data += ret.ndef_chunk_data
        end

        if ndef_out_of_sync # discard remaining stream to bring it back in-sync
            while ret.ndef_chunk_offset + 60 < ndef_length
                ret = reader_read_ndef_low_level(device, )
                    ndef_length = ret.ndef_length
            end

            throw(TinkerforgeStreamOutOfSyncError("NDEF stream is out-of-sync"))
        end
    end

        return ndef_data[:ndef_length]
end

export reader_write_page
"""
Writes a maximum of 8192 bytes starting from the given page. How many pages are written
depends on the tag type. The page sizes are as follows:

* Mifare Classic page size: 16 byte
* NFC Forum Type 1 page size: 8 byte
* NFC Forum Type 2 page size: 4 byte
* NFC Forum Type 3 page size: 16 byte
* NFC Forum Type 4: No pages, page = file selection (CC or NDEF, see below)

The general approach for writing to a tag is as follows:

1. Call :func:`Reader Request Tag ID`
2. Wait for state to change to *ReaderRequestTagIDReady* (see :func:`Reader Get State` or
   :cb:`Reader State Changed` callback)
3. If looking for a specific tag then call :func:`Reader Get Tag ID` and check if the
   expected tag was found, if it was not found got back to step 1
4. Call :func:`Reader Write Page` with page number and data
5. Wait for state to change to *ReaderWritePageReady* (see :func:`Reader Get State` or
   :cb:`Reader State Changed` callback)

If you use a Mifare Classic tag you have to authenticate a page before you
can write to it. See :func:`Reader Authenticate Mifare Classic Page`.

NFC Forum Type 4 tags are not organized into pages but different files. We currently
support two files: Capability Container file (CC) and NDEF file.

Choose CC by setting page to 3 or NDEF by setting page to 4.
"""
function reader_write_page(device::BrickletNFC, page, data)

        page = UInt16(page)
    data = convert(Vector{UInt8}, data)

    if length(data) > 65535
        throw(TinkerforgeInvalidParameterError("Data can be at most 65535 items long"))
    end

    data_length = length(data)
    data_chunk_offset = 0

    if data_length == 0
        data_chunk_data = [0] * 58
        ret = reader_write_page_low_level(device, page, data_length, data_chunk_offset, data_chunk_data)
    else
        lock(device.stream_lock) do
            while data_chunk_offset < data_length
                data_chunk_data = create_chunk_data(data, data_chunk_offset, 58, 0)
                ret = reader_write_page_low_level(device, page, data_length, data_chunk_offset, data_chunk_data)
                data_chunk_offset += 58
            end
        end
    end

        return ret
end

export reader_read_page
"""
Returns the page data from an internal buffer. To fill the buffer
with specific pages you have to call :func:`Reader Request Page` beforehand.
"""
function reader_read_page(device::BrickletNFC)

    lock(device.stream_lock) do
        ret = reader_read_page_low_level(device, )
            data_length = ret.data_length
        data_out_of_sync = ret.data_chunk_offset != 0

        data_data = ret.data_chunk_data

        while !data_out_of_sync && length(data_data) < data_length
            ret = reader_read_page_low_level(device, )
                data_length = ret.data_length
            data_out_of_sync = ret.data_chunk_offset != length(data_data)
            data_data += ret.data_chunk_data
        end

        if data_out_of_sync # discard remaining stream to bring it back in-sync
            while ret.data_chunk_offset + 60 < data_length
                ret = reader_read_page_low_level(device, )
                    data_length = ret.data_length
            end

            throw(TinkerforgeStreamOutOfSyncError("Data stream is out-of-sync"))
        end
    end

        return data_data[:data_length]
end

export cardemu_write_ndef
"""
Writes the NDEF message that is to be transferred to the NFC peer.

The maximum supported NDEF message size in Cardemu mode is 255 byte.

You can call this function at any time in Cardemu mode. The internal buffer
will not be overwritten until you call this function again or change the
mode.
"""
function cardemu_write_ndef(device::BrickletNFC, ndef)

        ndef = convert(Vector{UInt8}, ndef)

    if length(ndef) > 65535
        throw(TinkerforgeInvalidParameterError("NDEF can be at most 65535 items long"))
    end

    ndef_length = length(ndef)
    ndef_chunk_offset = 0

    if ndef_length == 0
        ndef_chunk_data = [0] * 60
        ret = cardemu_write_ndef_low_level(device, ndef_length, ndef_chunk_offset, ndef_chunk_data)
    else
        lock(device.stream_lock) do
            while ndef_chunk_offset < ndef_length
                ndef_chunk_data = create_chunk_data(ndef, ndef_chunk_offset, 60, 0)
                ret = cardemu_write_ndef_low_level(device, ndef_length, ndef_chunk_offset, ndef_chunk_data)
                ndef_chunk_offset += 60
            end
        end
    end

        return ret
end

export p2p_write_ndef
"""
Writes the NDEF message that is to be transferred to the NFC peer.

The maximum supported NDEF message size for P2P transfer is 255 byte.

You can call this function at any time in P2P mode. The internal buffer
will not be overwritten until you call this function again, change the
mode or use P2P to read an NDEF messages.
"""
function p2p_write_ndef(device::BrickletNFC, ndef)

        ndef = convert(Vector{UInt8}, ndef)

    if length(ndef) > 65535
        throw(TinkerforgeInvalidParameterError("NDEF can be at most 65535 items long"))
    end

    ndef_length = length(ndef)
    ndef_chunk_offset = 0

    if ndef_length == 0
        ndef_chunk_data = [0] * 60
        ret = p2p_write_ndef_low_level(device, ndef_length, ndef_chunk_offset, ndef_chunk_data)
    else
        lock(device.stream_lock) do
            while ndef_chunk_offset < ndef_length
                ndef_chunk_data = create_chunk_data(ndef, ndef_chunk_offset, 60, 0)
                ret = p2p_write_ndef_low_level(device, ndef_length, ndef_chunk_offset, ndef_chunk_data)
                ndef_chunk_offset += 60
            end
        end
    end

        return ret
end

export p2p_read_ndef
"""
Returns the NDEF message that was written by a NFC peer in NFC P2P mode.

The NDEF message is ready if you called :func:`P2P Start Transfer` with a
read transfer and the P2P state changed to *P2PTransferNDEFReady*.
"""
function p2p_read_ndef(device::BrickletNFC)

    lock(device.stream_lock) do
        ret = p2p_read_ndef_low_level(device, )
            ndef_length = ret.ndef_length
        ndef_out_of_sync = ret.ndef_chunk_offset != 0

        ndef_data = ret.ndef_chunk_data

        while !ndef_out_of_sync && length(ndef_data) < ndef_length
            ret = p2p_read_ndef_low_level(device, )
                ndef_length = ret.ndef_length
            ndef_out_of_sync = ret.ndef_chunk_offset != length(ndef_data)
            ndef_data += ret.ndef_chunk_data
        end

        if ndef_out_of_sync # discard remaining stream to bring it back in-sync
            while ret.ndef_chunk_offset + 60 < ndef_length
                ret = p2p_read_ndef_low_level(device, )
                    ndef_length = ret.ndef_length
            end

            throw(TinkerforgeStreamOutOfSyncError("NDEF stream is out-of-sync"))
        end
    end

        return ndef_data[:ndef_length]
end

export simple_get_tag_id
"""
.. versionadded:: 2.0.6\$nbsp;(Plugin)
"""
function simple_get_tag_id(device::BrickletNFC, index)

        index = UInt8(index)

    ret = simple_get_tag_id_low_level(device, index)

        return SimpleGetTagID(ret.tag_type, ret.tag_id_data[:ret.tag_id_length], ret.last_seen)
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletNFC, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using Distributed
using Test
addprocs(2)

@everywhere begin

using Distributed, Dagger, DaggerGPU
using CUDA, AMDGPU

function myfunc(X)
    @assert !(X isa Array)
    X
end

end

function generate_thunks()
    as = [delayed(x->x+1)(1) for i in 1:10]
    delayed((xs...)->[sum(xs)])(as...)
end

@testset "CUDA" begin
    if !DaggerGPU.cancompute(:CUDA)
        @warn "No CUDA devices available, skipping tests"
    else
        cuproc = DaggerGPU.processor(:CUDA)
        b = generate_thunks()
        opts = Dagger.Sch.ThunkOptions(;proctypes=[cuproc])
        c_pre = delayed(myfunc; options=opts)(b)
        c = delayed(sum; options=opts)(b)

        opts = Dagger.Sch.ThunkOptions(;proctypes=[Dagger.ThreadProc])
        d = delayed(identity; options=opts)(c)
        @test collect(d) == 20
    end
end

@testset "ROCm" begin
    if !DaggerGPU.cancompute(:ROC)
        @warn "No ROCm devices available, skipping tests"
    else
        rocproc = DaggerGPU.processor(:ROC)
        b = generate_thunks()
        opts = Dagger.Sch.ThunkOptions(;proctypes=[rocproc])
        c_pre = delayed(myfunc; options=opts)(b)
        c = delayed(sum; options=opts)(b)

        opts = Dagger.Sch.ThunkOptions(;proctypes=[Dagger.ThreadProc])
        d = delayed(identity; options=opts)(c)
        @test collect(d) == 20
    end
end


text:
using HDF5
using MLDatasets
using MultivariateStats
# include("DPMMSubClusters.jl")
using MAT
using Random
using CSV
using NPZ
# using DPMMSubCluste~rs

DS_PATH = "/vildata/vlad/splitnet/experiments_data/"
mkpath(DS_PATH)

##Mnist
function create_mnist(PATH, max_dim=128)
    mnist_x, mnist_y = MNIST.traindata()
    reshaped_x = copy(Float32.(reshape(mnist_x,:,size(mnist_x)[end])))
    M = MultivariateStats.fit(PCA, reshaped_x; maxoutdim=max_dim) #64 provides 0.86 var presereved
    # M = MultivariateStats.fit(PCA, reshaped_x; maxoutdim=72) #64 provides 0.86 var presereved
	tranformed_x = transform(M, reshaped_x)

    MNIST_PATH = PATH * "mnist_dim_"*string(max_dim)*".h5"

    ispath(MNIST_PATH) && rm(MNIST_PATH)
    h5write(MNIST_PATH, "data", tranformed_x)
    h5write(MNIST_PATH, "labels", mnist_y)

    # MAT_PATH = PATH*"mnist.mat"
    # matwrite(MAT_PATH, Dict(
	# "data" => tranformed_x,
	# "labels" => mnist_y))

	npzwrite(PATH*"mnist_x.npy", tranformed_x')
	npzwrite(PATH*"mnist_y.npy", mnist_y)
end

##Mnist
function create_fashion_mnist(PATH, max_dim=128)
    fmnist_x, fmnist_y = FashionMNIST.traindata()
    reshaped_x = copy(Float32.(reshape(fmnist_x,:,size(fmnist_x)[end])))
    M = MultivariateStats.fit(PCA, reshaped_x; maxoutdim=max_dim) #64 provides 0.86 var presereved
    tranformed_x = transform(M, reshaped_x)

    MNIST_PATH = PATH * "fmnist_dim_"*string(max_dim)*".h5"

    ispath(MNIST_PATH) && rm(MNIST_PATH)
    h5write(MNIST_PATH, "data", tranformed_x)
    h5write(MNIST_PATH, "labels", fmnist_y)

    MAT_PATH = PATH*"fmnist.mat"
    matwrite(MAT_PATH, Dict(
	"data" => tranformed_x,
	"labels" => fmnist_y))

	npzwrite(PATH*"fmnist_x.npy", tranformed_x')
	npzwrite(PATH*"fmnist_y.npy", fmnist_y)
end

##cifar 10
function create_cifar10(PATH, max_dim=88)
    cifar_x, cifar_y = CIFAR10.traindata()
    reshaped_x = copy(Float32.(reshape(cifar_x,:,size(cifar_x)[end])))
    M = MultivariateStats.fit(PCA, reshaped_x; maxoutdim=max_dim) #128 preserves 0.918 of the variance
	# M = MultivariateStats.fit(PCA, reshaped_x; maxoutdim=16) #128 preserves 0.918 of the variance
    tranformed_x = transform(M, reshaped_x)

    CIFAR_PATH = PATH * "cifar10_dim_"*string(max_dim)*".h5"

    ispath(CIFAR_PATH) && rm(CIFAR_PATH)
    h5write(CIFAR_PATH, "data", tranformed_x)
    h5write(CIFAR_PATH, "labels", cifar_y)

	MAT_PATH = PATH*"cifar10.mat"
	matwrite(MAT_PATH, Dict(
	"data" => tranformed_x,
	"labels" => cifar_y))

	npzwrite(PATH*"cifar_x.npy", tranformed_x')
	npzwrite(PATH*"cifar_y.npy", cifar_y)
end

##cifar 100
function create_cifar100(PATH, max_dim=96)
    cifar_x, _, cifar_y = CIFAR100.traindata()
    reshaped_x = copy(Float32.(reshape(cifar_x,:,size(cifar_x)[end])))
    M = MultivariateStats.fit(PCA, reshaped_x; maxoutdim=max_dim) #128 preserves 0.918 of the variance
    tranformed_x = transform(M, reshaped_x)

    CIFAR_PATH = PATH * "cifar100_dim_"*string(max_dim)*".h5"

    ispath(CIFAR_PATH) && rm(CIFAR_PATH)
    h5write(CIFAR_PATH, "data", tranformed_x)
    h5write(CIFAR_PATH, "labels", cifar_y)

	MAT_PATH = PATH*"cifar100.mat"
	matwrite(MAT_PATH, Dict(
	"data" => tranformed_x,
	"labels" => cifar_y))

	npzwrite(PATH*"cifar100_x.npy", tranformed_x')
	npzwrite(PATH*"cifar100_y.npy", cifar_y)
end

function create_svhn(PATH, max_dim=64)
    train_x, train_y = SVHN2.traindata()
    reshaped_x = copy(Float32.(reshape(train_x,:,size(train_x)[end])))
    M = MultivariateStats.fit(PCA, reshaped_x; maxoutdim=max_dim) #128 preserves 0.918 of the variance
    tranformed_x = transform(M, reshaped_x)

    CIFAR_PATH = PATH * "svhn_dim_"*string(max_dim)*".h5"

    ispath(CIFAR_PATH) && rm(CIFAR_PATH)
    h5write(CIFAR_PATH, "data", tranformed_x)
    h5write(CIFAR_PATH, "labels", train_y)

	MAT_PATH = PATH*"svhn.mat"
	matwrite(MAT_PATH, Dict(
	"data" => tranformed_x,
	"labels" => train_y))

	npzwrite(PATH*"svhn_x.npy", tranformed_x')
	npzwrite(PATH*"svhn_y.npy", train_y)
end

function create_2d_gaussian(PATH)
	Random.seed!(1) #generate same data
    x,y,clusters =  DPMMSubClusters.generate_gaussian_data(100000,2,10,Float32(20.0))

    SAVE_PATH = PATH * "gaussian2d.h5"

    ispath(SAVE_PATH) && rm(SAVE_PATH)
    h5write(SAVE_PATH, "data", x)
    h5write(SAVE_PATH, "labels", y)

	MAT_PATH = PATH*"gaussian2d.mat"
	matwrite(MAT_PATH, Dict(
	"data" => x,
	"labels" => y))

	npzwrite(PATH*"gaussian2d_x.npy", x')
	npzwrite(PATH*"gaussian2d_y.npy", y)

end

function create_3d_gaussian(PATH)
	Random.seed!(1234) #generate same data
    x,y,clusters =  DPMMSubClusters.generate_gaussian_data(1000000,3,6,Float32(3.0))

    SAVE_PATH = PATH * "gaussian3d.h5"

    ispath(SAVE_PATH) && rm(SAVE_PATH)
    h5write(SAVE_PATH, "data", x)
    h5write(SAVE_PATH, "labels", y)

	MAT_PATH = PATH*"gaussian3d.mat"
	matwrite(MAT_PATH, Dict(
	"data" => x,
	"labels" => y))

	npzwrite(PATH*"gaussian3d_x.npy", x')
	npzwrite(PATH*"gaussian3d_y.npy", y)

end

function create_5d_gaussian(PATH)
	Random.seed!(1234) #generate same data
    x,y,clusters = DPMMSubClusters.generate_gaussian_data(1000000,5,40,Float32(2.0))

    SAVE_PATH = PATH * "gaussian5d.h5"

    ispath(SAVE_PATH) && rm(SAVE_PATH)
    h5write(SAVE_PATH, "data", x)
    h5write(SAVE_PATH, "labels", y)

	MAT_PATH = PATH*"gaussian5d.mat"
	matwrite(MAT_PATH, Dict(
	"data" => x,
	"labels" => y))

	npzwrite(PATH*"gaussian5d_x.npy", x')
	npzwrite(PATH*"gaussian5d_y.npy", y)
end

function create_10d_gaussian(PATH)
	Random.seed!(1234) #generate same data
    x,y,clusters = DPMMSubClusters.generate_gaussian_data(400000,10,100,Float32(0.1))

    SAVE_PATH = PATH * "gaussian10d.h5"

    ispath(SAVE_PATH) && rm(SAVE_PATH)
    h5write(SAVE_PATH, "data", x)
    h5write(SAVE_PATH, "labels", y)

	MAT_PATH = PATH*"gaussian10d.mat"
	matwrite(MAT_PATH, Dict(
	"data" => x,
	"labels" => y))

	npzwrite(PATH*"gaussian10d_x.npy", x')
	npzwrite(PATH*"gaussian10d_y.npy", y)
end

function create_30d_gaussian(PATH)
	Random.seed!(1234) #generate same data
    x,y,clusters = DPMMSubClusters.generate_gaussian_data(1000000,30,80,Float32(0.00000000001))

    SAVE_PATH = PATH * "gaussian30d.h5"

    ispath(SAVE_PATH) && rm(SAVE_PATH)
    h5write(SAVE_PATH, "data", x)
    h5write(SAVE_PATH, "labels", y)

	MAT_PATH = PATH*"gaussian30d.mat"
	matwrite(MAT_PATH, Dict(
	"data" => x,
	"labels" => y))

	npzwrite(PATH*"gaussian30d_x.npy", x')
	npzwrite(PATH*"gaussian30d_y.npy", y)
end

function create_100d_gaussian(PATH)
	Random.seed!(1234) #generate same data
    x,y,clusters = DPMMSubClusters.generate_gaussian_data(100000,100,50,Float32(0.001))

    SAVE_PATH = PATH * "gaussian100d.h5"

    ispath(SAVE_PATH) && rm(SAVE_PATH)
    h5write(SAVE_PATH, "data", x)
    h5write(SAVE_PATH, "labels", y)

	MAT_PATH = PATH*"gaussian100d.mat"
	matwrite(MAT_PATH, Dict(
	"data" => x,
	"labels" => y))

	npzwrite(PATH*"gaussian100d_x.npy", x')
	npzwrite(PATH*"gaussian100d_y.npy", y)
end

function create_otto(PATH, PCA = false)
	ds = CSV.read("/home/dinari/smart_splits_experiments/dpgmm_smartsplits/otto-group-product-classification-challenge/train.csv")
	gt = ds[:,end]
	factorize_dict = Dict([v=>k for (k,v) in enumerate(unique(gt))])
	gt = (x -> factorize_dict[x]).(gt)

	data = Float32.(convert(Array,ds[:,2:end-1])')

	M = MultivariateStats.fit(MultivariateStats.PCA, data; maxoutdim=128)
	x_reshaped = transform(M, data)
	if PCA
		data = x_reshaped
	end
	SAVE_PATH = PATH * "otto.h5"

	ispath(SAVE_PATH) && rm(SAVE_PATH)
	h5write(SAVE_PATH, "data", data)
	h5write(SAVE_PATH, "labels", gt)

	MAT_PATH = PATH*"otto.mat"
	matwrite(MAT_PATH, Dict(
	"data" => data,
	"labels" => gt))

	npzwrite(PATH*"otto_x.npy", data')
	npzwrite(PATH*"otto_y.npy", gt)
end

function create_frogs(PATH)
	ds = CSV.read("/home/dinari/smart_splits_experiments/raw_datasets/Anuran Calls (MFCCs)/Frogs_MFCCs.csv")
	gt = ds[:,end-1]
	factorize_dict = Dict([v=>k for (k,v) in enumerate(unique(gt))])
	gt = (x -> factorize_dict[x]).(gt)

	data = Float32.(convert(Array,ds[:,1:end-4])')

	M = MultivariateStats.fit(PCA, data; maxoutdim=4) #64 provides 0.86 var presereved
	tranformed_x = transform(M, data)

	# data = tranformed_x

	SAVE_PATH = PATH * "frogs.h5"

	ispath(SAVE_PATH) && rm(SAVE_PATH)
	h5write(SAVE_PATH, "data", data)
	h5write(SAVE_PATH, "labels", gt)

	MAT_PATH = PATH*"frogs.mat"
	matwrite(MAT_PATH, Dict(
	"data" => data,
	"labels" => gt))

	npzwrite(PATH*"frogs_x.npy", data')
	npzwrite(PATH*"frogs_y.npy", gt)
end

function create_wine(PATH)
	ds = CSV.read("/home/dinari/smart_splits_experiments/raw_datasets/winequality-white.csv")
	gt = ds[:,end]
	factorize_dict = Dict([v=>k for (k,v) in enumerate(unique(gt))])
	gt = (x -> factorize_dict[x]).(gt)

	data = Float32.(convert(Array,ds[:,1:end-1])')

	M = MultivariateStats.fit(PCA, data; maxoutdim=10) #64 provides 0.86 var presereved
	tranformed_x = transform(M, data)

	data = tranformed_x

	SAVE_PATH = PATH * "wine.h5"

	ispath(SAVE_PATH) && rm(SAVE_PATH)
	h5write(SAVE_PATH, "data", data)
	h5write(SAVE_PATH, "labels", gt)

	MAT_PATH = PATH*"wine.mat"
	matwrite(MAT_PATH, Dict(
	"data" => data,
	"labels" => gt))

	npzwrite(PATH*"wine_x.npy", data')
	npzwrite(PATH*"wine_y.npy", gt)
end

function create_yeast(PATH)
	ds = CSV.read("/home/dinari/smart_splits_experiments/raw_datasets/yeast.data"; delim=',')
	gt = ds[:,end]
	factorize_dict = Dict([v=>k for (k,v) in enumerate(unique(gt))])
	gt = (x -> factorize_dict[x]).(gt)

	data = Float32.(convert(Array,ds[:,3:end-1])')

	M = MultivariateStats.fit(PCA, data; maxoutdim=3) #64 provides 0.86 var presereved
	tranformed_x = transform(M, data)

	data = tranformed_x

	SAVE_PATH = PATH * "yeast.h5"

	ispath(SAVE_PATH) && rm(SAVE_PATH)
	h5write(SAVE_PATH, "data", data)
	h5write(SAVE_PATH, "labels", gt)

	MAT_PATH = PATH*"yeast.mat"
	matwrite(MAT_PATH, Dict(
	"data" => data,
	"labels" => gt))

	npzwrite(PATH*"yeast_x.npy", data')
	npzwrite(PATH*"yeast_y.npy", gt)
end

function create_bees(PATH)
	ds = CSV.read("/home/dinari/smart_splits_experiments/raw_datasets/Tetragonula.csv")
	gt = ds[:,end]
	factorize_dict = Dict([v=>k for (k,v) in enumerate(unique(gt))])
	gt = (x -> factorize_dict[x]).(gt)

	data = Float32.(convert(Array,ds[:,1:end-3])') ./ 10000

	# M = MultivariateStats.fit(PCA, data; maxoutdim=3) #64 provides 0.86 var presereved
	# tranformed_x = transform(M, data)
	#
	# data = tranformed_x

	SAVE_PATH = PATH * "bees.h5"

	ispath(SAVE_PATH) && rm(SAVE_PATH)
	h5write(SAVE_PATH, "data", data)
	h5write(SAVE_PATH, "labels", gt)

	MAT_PATH = PATH*"bees.mat"
	matwrite(MAT_PATH, Dict(
	"data" => data,
	"labels" => gt))

	npzwrite(PATH*"bees_x.npy", data')
	npzwrite(PATH*"bees_y.npy", gt)
end

# Synthetic Data:
# create_2d_gaussian(DS_PATH)
# create_3d_gaussian(DS_PATH)
# create_5d_gaussian(DS_PATH)
# create_10d_gaussian(DS_PATH)
# create_30d_gaussian(DS_PATH)
# create_100d_gaussian(DS_PATH)

# Real Data:
# create_mnist(DS_PATH)
# create_fashion_mnist(DS_PATH, 128)
# create_svhn(DS_PATH)
create_cifar10(DS_PATH, 128)
create_cifar100(DS_PATH, 128)
# create_otto(DS_PATH)
# create_frogs(DS_PATH)
# create_wine(DS_PATH)
# create_yeast(DS_PATH)
# create_bees(DS_PATH)


text:
str1 = "string"
println(str1)
str2 = """String"""
println(str2)
str3 = """ この記号  " " を使いたいときは3つの方を使う"""
println(str3)

println("'a'はchr型として扱われる")
println(typeof('a'))

num = 10
println("この値は$num")
println("$num + $num = $(num+num)")

println("文字列結合")
println(string("私","は","犬です"))

println(string("私","は","犬です",num))

println("文字結合は*を"*"使えばできるらしい")

chr1 = 'a'
chr2 = 'b'
println(string(chr1,chr2))
println(typeof(string(chr1,chr2)))
println(chr1*chr2)
println(typeof(chr1*chr2))


text:
# Implement MotifEnumeration

const NUCLEOTIDES = ('A', 'C', 'G', 'T')

convert_to_intlist(line) = map(x -> parse(Int, x), split(line))

calc_hamming_distance(s1, s2) = count(((c1, c2),) -> c1 != c2, zip(s1, s2))

function calc_immediate_neighbours(pattern)
    neighbours = Set{String}()
    for (ix, letter) in enumerate(pattern)
        for nucleotide in NUCLEOTIDES
            if nucleotide != letter
                pattern_array = collect(pattern)
                pattern_array[ix] = nucleotide
                push!(neighbours, join(pattern_array))
            end
        end
    end
    neighbours
end

function generate_neighbourhood(text, d)
    neighbourhood = Set{String}([text])
    for _ in 1:d
        neighbours_of_neighbours = Set{String}()
        for pattern in neighbourhood
            immediate_neighbours = calc_immediate_neighbours(pattern)
            union!(neighbours_of_neighbours, immediate_neighbours)
        end
        union!(neighbourhood, neighbours_of_neighbours)
    end
    neighbourhood
end

function does_k_mer_occurs_approximately_in_text(text, kmer, k, d)
    for ix in 1:(length(text)-k+1)
        substring = text[ix:ix+k-1]
        if calc_hamming_distance(substring, kmer) <= d
            return true
        end
    end
    false
end

function motif_enumeration(texts, k, d)
    all_texts_approximately_contain_pattern(kmer) =
        all(does_k_mer_occurs_approximately_in_text(text, kmer, k, d) for text in texts)

    motifs = Set{String}()
    for text in texts
        for ix in 1:(length(text)-k+1)
            pattern = text[ix:ix+k-1]
            neighbours = generate_neighbourhood(pattern, d)
            for neighbour in neighbours
                if all_texts_approximately_contain_pattern(neighbour)
                    push!(motifs, neighbour)
                end
            end
        end
    end
    motifs
end

function main()
    k, d = convert_to_intlist(readline())
    texts = readlines()
    result = motif_enumeration(texts, k, d)
    println(join(result, " "))
end

main()


text:
#=

Code for representing sparse matrices in IJV format.
It is often more efficient to manipulate these.
Our graph generators with _ijv at the end of their name produce these.

=#

import Base: ==, hash, +, *, transpose, adjoint
import SparseArrays: nnz

mutable struct IJV{Tv,Ti}
    n::Ti
    nnz::Ti
    i::Array{Ti,1}
    j::Array{Ti,1}
    v::Array{Tv,1}
end

==(a::IJV, b::IJV) =
    a.n == b.n &&
    a.nnz == b.nnz &&
    a.i == b.i &&
    a.j == b.j &&
    a.v == b.v

nnz(a::IJV) = a.nnz

hash(a::IJV) =
    hash((a.n, a.nnz, a.i, a.j, a.v), hash(IJV))

hash(a::IJV, h::UInt) = hash(hash(a), h)

function +(a::IJV, b::IJV)
    @assert a.n == b.n
    IJV(a.n, a.nnz+b.nnz,
        [a.i; b.i],
        [a.j; b.j],
        [a.v; b.v])
end

function *(a::IJV, x::Number)
    ijv = deepcopy(a)
    ijv.v .*= x

    return ijv
end

*(x::Number, a::IJV) = *(a, x)

transpose(ijv::IJV) = IJV(ijv.n, ijv.nnz, ijv.j, ijv.i, ijv.v)
adjoint(ijv::IJV) = IJV(ijv.n, ijv.nnz, ijv.j, ijv.i, adjoint.(ijv.v))

"""
    ijv = IJV(A::SparseMatrixCSC)
Convert a sparse matrix to an IJV.
"""
function IJV(A::SparseMatrixCSC{Tv, Ti}) where {Tv, Ti}
    (ai,aj,av) = findnz(A)
    IJV{Tv,Ti}(A.n, nnz(A), ai, aj, av)
end

import SparseArrays.sparse

sparse(ijv::IJV) = sparse(ijv.i, ijv.j, ijv.v, ijv.n, ijv.n)

compress(ijv::IJV) = IJV(sparse(ijv))