text:
# # Decaying barotropic QG turbulence over topography
#
#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/barotropicqg_decay_topography.ipynb).
#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/barotropicqg_decay_topography.ipynb).
#
# An example of decaying barotropic quasi-geostrophic turbulence over topography.

using FourierFlows, Plots, Printf, Random

using Statistics: mean
using FFTW: irfft

import GeophysicalFlows.BarotropicQG
import GeophysicalFlows.BarotropicQG: energy, enstrophy

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)
nothing # hide

# ## Numerical parameters and time-stepping parameters

      n = 128            # 2D resolution = n²
stepper = "FilteredRK4"  # timestepper
     dt = 0.05           # timestep
 nsteps = 2000           # total number of time-steps
 nsubs  = 10             # number of time-steps for intermediate logging/plotting (nsteps must be multiple of nsubs)
nothing # hide

# ## Physical parameters

L = 2π        # domain size
nothing # hide

# Define the topographic potential vorticity, ``\eta = f_0 h(x, y)/H``. The topography here
# is an elliptical mound at ``(1, 1)``, and an elliptical depression at ``(-1, -1)``.
σx, σy = 0.4, 0.8
topographicPV(x, y) = 3exp(-(x-1)^2/(2σx^2) -(y-1)^2/(2σy^2)) - 2exp(-(x+1)^2/(2σx^2) -(y+1)^2/(2σy^2))
nothing # hide

# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments. Not providing
# a viscosity coefficient `ν` leads to the module's default value: `ν=0`. In this
# example numerical instability due to accumulation of enstrophy in high wavenumbers
# is taken care with the `FilteredTimestepper` we picked. The topophic PV is prescribed
# via keyword argument `eta`.
prob = BarotropicQG.Problem(dev; nx=n, Lx=L, eta=topographicPV, dt=dt, stepper=stepper)
nothing # hide

# and define some shortcuts
sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
x, y = grid.x, grid.y
nothing # hide

# and let's plot the topographic PV:
contourf(grid.x, grid.y, params.eta',
          aspectratio = 1,
            linewidth = 0,
               levels = 10,
                    c = :balance,
                 clim = (-3, 3),
                xlims = (-grid.Lx/2, grid.Lx/2),
                ylims = (-grid.Ly/2, grid.Ly/2),
               xticks = -3:3,
               yticks = -3:3,
               xlabel = "x",
               ylabel = "y",
                title = "topographic PV η=f₀h/H")

# ## Setting initial conditions

# Our initial condition consist of a flow that has power only at wavenumbers with
# ``6 < \frac{L}{2\pi} \sqrt{k_x^2 + k_y^2} < 12`` and initial energy ``E_0``:

E₀ = 0.04 # energy of initial condition

K = @. sqrt(grid.Krsq)                          # a 2D array with the total wavenumber

Random.seed!(1234)
qih = randn(Complex{eltype(grid)}, size(sol))
@. qih = ifelse(K < 6  * 2π/L, 0, qih)
@. qih = ifelse(K > 12 * 2π/L, 0, qih)
qih *= sqrt(E₀ / energy(qih, grid))             # normalize qi to have energy E₀
qi = irfft(qih, grid.nx)

BarotropicQG.set_zeta!(prob, qi)
nothing # hide

# Let's plot the initial vorticity field:

p1 = heatmap(x, y, vars.q',
         aspectratio = 1,
              c = :balance,
           clim = (-8, 8),
          xlims = (-grid.Lx/2, grid.Lx/2),
          ylims = (-grid.Ly/2, grid.Ly/2),
         xticks = -3:3,
         yticks = -3:3,
         xlabel = "x",
         ylabel = "y",
          title = "initial vorticity ζ=∂v/∂x-∂u/∂y",
     framestyle = :box)

p2 = contourf(x, y, vars.psi',
        aspectratio = 1,
             c = :viridis,
        levels = range(-0.35, stop=0.35, length=10),
          clim = (-0.35, 0.35),
         xlims = (-grid.Lx/2, grid.Lx/2),
         ylims = (-grid.Ly/2, grid.Ly/2),
        xticks = -3:3,
        yticks = -3:3,
        xlabel = "x",
        ylabel = "y",
         title = "initial streamfunction ψ",
    framestyle = :box)

l = @layout Plots.grid(1, 2)
p = plot(p1, p2, layout=l, size=(900, 400))

# ## Diagnostics

# Create Diagnostics -- `energy` and `enstrophy` functions are imported at the top.
E = Diagnostic(energy, prob; nsteps=nsteps)
Z = Diagnostic(enstrophy, prob; nsteps=nsteps)
diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will  be updated every timestep.
nothing # hide

# ## Output

# We choose folder for outputing `.jld2` files.
filepath = "."
filename = joinpath(filepath, "decayingbetaturb.jld2")
nothing # hide

# Do some basic file management,
if isfile(filename); rm(filename); end
nothing # hide

# and then create Output.
get_sol(prob) = sol # extracts the Fourier-transformed solution
out = Output(prob, filename, (:sol, get_sol))
nothing # hide

# ## Visualizing the simulation

# We define a function that plots the vorticity and streamfunction and
# their corresponding zonal mean structure.

function plot_output(prob)
  ζ = prob.vars.zeta
  ψ = prob.vars.psi
  η = prob.params.eta

  pζ = heatmap(x, y, ζ',
       aspectratio = 1,
            legend = false,
                 c = :balance,
              clim = (-6, 6),
             xlims = (-grid.Lx/2, grid.Lx/2),
             ylims = (-grid.Ly/2, grid.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "vorticity ζ=∂v/∂x-∂u/∂y",
        framestyle = :box)

  contour!(pζ, x, y, η',
          levels=0.5:0.5:3,
          lw=2, c=:black, ls=:solid, alpha=0.7)

  contour!(pζ, x, y, η',
          levels=-2:0.5:-0.5,
          lw=2, c=:black, ls=:dash, alpha=0.7)

  pψ = contourf(x, y, ψ',
       aspectratio = 1,
            legend = false,
                 c = :viridis,
            levels = range(-0.7, stop=0.7, length=10),
              clim = (-0.7, 0.7),
             xlims = (-grid.Lx/2, grid.Lx/2),
             ylims = (-grid.Ly/2, grid.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "streamfunction ψ",
        framestyle = :box)

  l = @layout Plots.grid(1, 2)
  p = plot(pζ, pψ, layout = l, size = (900, 400))

  return p
end
nothing # hide

# ## Time-stepping the `Problem` forward

# We time-step the `Problem` forward in time.

startwalltime = time()

p = plot_output(prob)

anim = @animate for j = 0:round(Int, nsteps/nsubs)

  if j % (1000 / nsubs) == 0
    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])

    log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min",
      clock.step, clock.t, cfl, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)

    println(log)
  end

  p[1][1][:z] = vars.zeta
  p[1][:title] = "vorticity, t="*@sprintf("%.2f", clock.t)
  p[2][1][:z] = vars.psi

  stepforward!(prob, diags, nsubs)
  BarotropicQG.updatevars!(prob)
end

gif(anim, "barotropicqg_decay_topography.gif", fps=12)


text:
@discord_object struct Embed
    title::String
    type::String
    description::String
    url::String
    timestamp::Union{String, DateTime}
    color::Int
    footer::EmbedFooter
    image::EmbedImage
    thumbnail::EmbedThumbnail
    video::EmbedVideo
    provider::EmbedProvider
    author::EmbedAuthor
    fields::Vector{EmbedField}
end

ext: jl


text:
using Dates
using Printf

@enum GeneratorType begin
    LIMITABLE
    PILOTABLE
end

function Base.string(g::GeneratorType)
    if g==PILOTABLE
        return "Pilotable"
    elseif g==LIMITABLE
        return "Limitable"
    else
        throw( error("Undefined conversion of GeneratorType `", g, "` to a string") )
    end
end

function Base.parse(::Type{GeneratorType}, str::String)
    if lowercase(str) == "limitable"
        return LIMITABLE
    elseif  lowercase(str) == "pilotable"
        return PILOTABLE
    else
        throw( error("Unable to convert `", str, "` to a GeneratorType") )
    end
end

struct Generator
    id::String
    bus_id::String # pas directement un Bus, parce que ca fait des porblemes de references circulaires

    type::GeneratorType
    p_min::Float64
    p_max::Float64
    start_cost::Float64
    prop_cost::Float64
    dmo::Dates.Second
    dp::Dates.Second

    function Generator(id::String, bus_id::String,
                        type::GeneratorType,
                        p_min::Float64, p_max::Float64,
                        start_cost::Float64, prop_cost::Float64,
                        dmo::Dates.Second, dp::Dates.Second)
        return new(id, bus_id, type, p_min, p_max, start_cost, prop_cost, dmo, dp)
    end
end

function get_bus_id(gen::Generator)
    return gen.bus_id
end

function get_type(gen::Generator)
    return gen.type
end

function get_p_min(gen::Generator)
    return gen.p_min
end

function get_p_max(gen::Generator)
    return gen.p_max
end

function get_start_cost(gen::Generator)
    return gen.start_cost
end

function get_prop_cost(gen::Generator)
    return gen.prop_cost
end

function get_dmo(gen::Generator)
    return gen.dmo
end

function get_dp(gen::Generator)
    return gen.dp
end

################
##  HELPERS   ##
################

function is_limitable(generator::Generator)::Bool
    return get_type(generator) == LIMITABLE
end

function is_pilotable(generator::Generator)::Bool
    return get_type(generator) == PILOTABLE
end

function needs_commitment(generator::Generator)::Bool
    return (get_type(generator) == PILOTABLE) && (get_p_min(generator) > 0)
end

################
## INFO / LOG ##
################

function get_info(generator::Generator)::String
    info::String =
        generator.id
    return info
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using CUDA.CUSPARSE
using LinearAlgebra, SparseArrays

@testset "LinearAlgebra" begin
    @testset "$f(A)*b $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64],
                                 f in (identity, transpose, adjoint)
        n = 10
        alpha = rand()
        beta = rand()
        A = sprand(elty, n, n, rand())
        b = rand(elty, n)
        c = rand(elty, n)
        alpha = beta = 1.0
        c = zeros(elty, n)

        dA = CuSparseMatrixCSR(A)
        db = CuArray(b)
        dc = CuArray(c)

        # test with empty inputs
        @test Array(dA * CUDA.zeros(elty, n, 0)) == zeros(elty, n, 0)

        mul!(c, f(A), b, alpha, beta)
        mul!(dc, f(dA), db, alpha, beta)
        @test c ≈ collect(dc)

        A = A + transpose(A)
        dA = CuSparseMatrixCSR(A)

        @assert issymmetric(A)
        mul!(c, f(Symmetric(A)), b, alpha, beta)
        mul!(dc, f(Symmetric(dA)), db, alpha, beta)
        @test c ≈ collect(dc)
    end

    @testset "$f(A)*$h(B) $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64],
                                     f in (identity, transpose, adjoint),
                                     h in (identity, transpose, adjoint)
        if CUSPARSE.version() <= v"10.3.1" &&
            (h ∈ (transpose, adjoint) && f != identity) ||
            (h == adjoint && elty <: Complex)
            # These combinations are not implemented in CUDA10
            continue
        end

        n = 10
        alpha = rand()
        beta = rand()
        A = sprand(elty, n, n, rand())
        B = rand(elty, n, n)
        C = rand(elty, n, n)

        dA = CuSparseMatrixCSR(A)
        dB = CuArray(B)
        dC = CuArray(C)

        mul!(C, f(A), h(B), alpha, beta)
        mul!(dC, f(dA), h(dB), alpha, beta)
        @test C ≈ collect(dC)

        A = A + transpose(A)
        dA = CuSparseMatrixCSR(A)

        @assert issymmetric(A)
        mul!(C, f(Symmetric(A)), h(B), alpha, beta)
        mul!(dC, f(Symmetric(dA)), h(dB), alpha, beta)
        @test C ≈ collect(dC)
    end
end


text:
using LinearAlgebra
using SparseArrays
using Plots
using ForwardDiff

"This routine solves Poisson's equation -u''(x) = f(x) in weak form
(∫κ(x)*u'(x)*ϕi'(x) = ∫f(x)*ϕi(x) for i = 1,...,m)
using piecewise linear finite element methods. The matrix is assembled locally,
and boundary nodes are not removed."

m = 50 # number of elements

# α = 0 # Dirichlet BC
# β = 0 # Neumann BC
# f(x) = pi^2*sin(pi*x) #exp(2*sin(pi*x))
# κ(x) = 1

# Manufactured solution
uexact(x) = log(2+sin(pi*x))
κ(x) = 1 + .5*sin(pi*x)
dudx_exact(x) = ForwardDiff.derivative(uexact,x)
f(x) = -ForwardDiff.derivative(x->κ(x)*dudx_exact(x),x) # -d/dx (κ(x)*du/dx) = f
α = uexact(-1)
β = κ(1)*dudx_exact(1)

# define spatial grid
x = LinRange(-1,1,m+2) # x_0, x_1, ..., x_m, x_{m+1} = x_0
# x = @. x + randn() / (2*m) # randomly perturb nodes

A = spzeros(m+2,m+2)
A_local = [1 -1;-1 1]
b = zeros(m+2)
for e = 1:m+1 # loop over the intervals
    ids = e:e+1
    h_e = x[e+1]-x[e]

    x_midpt = .5*(x[e+1]+x[e])
    A[ids,ids] .+= κ(x_midpt)*A_local * 1/h_e

    b[ids] .+= h_e / 2.0 * f(.5*(x[e]+x[e+1])) # f(.5*(x_i+x_i+1))
end

# impose Dirichlet BCs on both sides
b̃ = b - A[:,1]*α
b̃[1] = α
A[1,:] .= 0
A[:,1] .= 0
A[1,1] = 1

# b̃ = b̃ - A[:,m+2]*β
# b̃[m+2] = β
# A[m+2,:] .= 0
# A[:,m+2] .= 0
# A[m+2,m+2] = 1

b̃[m+2] += β # impose Neumann BC on the left

# convert b̃ from a SparseVector to a Vector
# (Julia doesn't currently allow "\" to be applied to SparseVectors)
u = A\Vector(b̃)

err = abs.(u - uexact.(x))
plot(x,u,mark=:dot,legend=false,title="Max error = $(maximum(err))")


text:
using Keras
using Flux
using Base.Test
using DataFlow:Call, constant, vertex, syntax

Keras.new_type(a, b) = Keras.new_type(a, nothing, b)

vcall(a...) = vertex(Call(), constant.(a)...)

a = rand(10,10)

@testset "Activation" begin
    #relu test
    temp = vcall(Keras.ops[:relu], a)
    @test vcall(temp, a) |> syntax |> eval == relu.(a)

    #tanh test
    temp = vcall(Keras.ops[:tanh], a)
    @test vcall(temp, a) |> syntax |> eval == tanh.(a)

    #sigmoid test
    temp = vcall(Keras.ops[:sigmoid], a)
    @test vcall(temp, a) |> syntax |> eval == sigmoid(a)

    #elu test
    temp = vcall(Keras.ops[:elu], a)
    @test vcall(:broadcast, temp, a) |> syntax |> eval == elu.(a)

    #softplus test
    temp = vcall(Keras.ops[:softplus], a)
    @test vcall(temp, a) |> syntax |> eval == softplus.(a)

    #softmax test
    a = rand(10)
    temp = vcall(Keras.ops[:softmax], a)
    #@test vcall(temp, a) |> syntax |> eval == softmax(a)

    #ZeroPadding2D test
    a = rand(5,5,1,1)
    d = Dict{Any, Any}()
    d["padding"] = [[2,2], [2,2]]
    b = Keras.new_type(:ZeroPadging2D, d)
    @test sum(meanpool(a, (1,1), pad=(2,2), stride=(1,1))) ==
                sum((Keras.ops[:ZeroPadding2D](b))(a))
    @test collect(size(meanpool(a, (1,1), pad=(2,2), stride=(1,1)))) - collect(size((Keras.ops[:ZeroPadding2D](b))(a))) ==
                [0,0,0,0]
end


text:
# Unit Test for by split vector type
using COSMO, Test, Random, LinearAlgebra

# This test is precision agnostic and will be run with T precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end

@testset "SplitVector" begin
for T in UnitTestFloats
    @testset "SplitVector (T = $(T))" begin
        rng = Random.MersenneTwister(13131)

        convex_set = COSMO.PsdConeTriangle{T}(3)
        composite_set = COSMO.CompositeConvexSet{T}([COSMO.ZeroSet{T}(3), COSMO.Nonnegatives{T}(1)])
        x1 = rand(T, 3)
        x2 = rand(T, 4)

        s1 = COSMO.SplitVector(x1, convex_set)
        s2 = COSMO.SplitVector(x2, composite_set)
        @test size(s1) == (3, )
        @test length(s2) == 4
        @test getindex(s1, 2) == x1[2]
        setindex!(s1, 5, 2)
        @test getindex(s1, 2) == 5.
        @test firstindex(s2) == 1
        @test lastindex(s2) == 4
        @test Base.showarg(IOBuffer(), s2, true ) == nothing
    end
end
end


text:
using Nextion

@testset "qrcode" begin
    nexSerial = NexSerial("/dev/ttyUSB0")

    nexQRcode = NexQRcode(nexSerial, Name("qr0"), cid=ComponentID(1))
    nexText = NexText(nexSerial, Name("t1"), cid=ComponentID(3))

    send(nexSerial, "page pg_qr")

    sleep(1)

    msg = "Hello"
    nexQRcode.text = msg
    nexText.text = msg

    @test nexQRcode.text = msg

    sleep(2)

    # text = "https://github.com/scls19fr/Nextion.jl"
    msg = "http://bit.ly/2vWS75r"
    nexText.text = msg
    #nexQRcode.textmaxlength = length(text)
    #nexQRcode.textmaxlength = 50
    #sleep(1)
    nexQRcode.text = msg

    close(nexSerial)
end


text:
using PastaQ
using ITensors
using Test
using Random

@testset "nqubits/nlayers/ngates" begin
  @test nqubits(("H", 2)) == 2
  @test nqubits(("CX", (2, 5))) == 5

  for i in 1:1000
    depth = 4
    N = rand(3:50)
    gates = randomcircuit(N; depth = depth, twoqubitgates="CX", onequbitgates="Rn")
    n = nqubits(gates)
    @test N == n
    @test PastaQ.nlayers(gates) == depth
    @test PastaQ.ngates(gates) == depth ÷ 2 * (2 * N ÷ 2 - 1 + 2 * N)
    gates = randomcircuit(N; depth = depth, twoqubitgates="CX", onequbitgates="Rn", layered=false)
    n = nqubits(gates)
    @test N == n
    @test PastaQ.nlayers(gates) == 1
    @test PastaQ.ngates(gates) == depth ÷ 2 * (2 * N ÷ 2 - 1 + 2 * N)
  end

  N = 3

  # MPS
  circuit = randomcircuit(N; depth = 4)
  ψ = runcircuit(circuit; full_representation = true)
  @test PastaQ.nsites(ψ) == N
  # MPO
  U = runcircuit(circuit; full_representation = true, process = true)
  @test PastaQ.nsites(U) == N
  # LPDO DM
  ρ = prod(randomstate(N; ξ = 2))
  @test PastaQ.nsites(ρ) == N

  # MPO Choi
  Λ = runcircuit(circuit; noise = ("DEP",(p=0.01,)), full_representation = true, process = true)
  @test PastaQ.nsites(Λ) == N
  Λ = prod(randomprocess(N; ξ = 2))
  @test PastaQ.nsites(Λ) == N
end

@testset "pre-defined ciruits" begin
  N = 10
  @test length(qft(N)) == sum(1:(N - 1)) + N
  @test length(qft(N; inverse=true)) == sum(1:(N - 1)) + N
  @test length(ghz(N)) == N
end

@testset "gatelayer" begin
  N = 10

  layer = gatelayer("X", N)
  @test length(layer) == N
  @test all(x -> x == "X", first.(layer))

  layer = gatelayer("X", 1:2:N)
  @test length(layer) == N ÷ 2

  layer = gatelayer("X", [1, 3, 5, 7])
  @test length(layer) == 4

  layer = gatelayer("Rx", N; θ=0.1)
  @test haskey.(last.(layer), :θ) == ones(N)
  @test values.(last.(layer)) == [(0.1,) for _ in 1:N]

  qarray = PastaQ.lineararray(N)
  layer = gatelayer("CX", qarray[1])
  @test length(layer) == length(qarray[1])
  @test all(x -> x == "CX", first.(layer))
  layer = gatelayer("CX", qarray[2])
  @test length(layer) == length(qarray[2])
  @test all(x -> x == "CX", first.(layer))
end

@testset "randommparams" begin
  N = 10
  pars = PastaQ.randomparams("Rx")
  @test haskey(pars, :θ)
  pars = PastaQ.randomparams("Ry")
  @test haskey(pars, :θ)
  pars = PastaQ.randomparams("Rz")
  @test haskey(pars, :ϕ)
  pars = PastaQ.randomparams("Rn")
  @test haskey(pars,:θ)
  @test haskey(pars,:ϕ)
  @test haskey(pars,:λ)

  pars = PastaQ.randomparams("RandomUnitary",2)
  @test haskey(pars,:random_matrix)
  @test size(pars[:random_matrix]) == (4,4)
  @test pars[:random_matrix] isa Matrix{ComplexF64}
end

@testset "randomlayer" begin
  N = 10
  layer = randomlayer("Rn", N)
  @test length(layer) == N
  @test all(x -> x == "Rn", first.(layer))

  qarray = PastaQ.lineararray(N)
  layer = randomlayer("RandomUnitary", qarray[1])
  @test length(layer) == length(qarray[1])
  @test all(x -> x == "RandomUnitary", first.(layer))
  layer = randomlayer("RandomUnitary", qarray[2])
  @test length(layer) == length(qarray[2])
  @test all(x -> x == "RandomUnitary", first.(layer))
end

@testset "random circuits" begin
  N = 30
  depth = 10
  circuit = randomcircuit(N; depth = depth, twoqubitgates="RandomUnitary")
  @test length(circuit) == depth
  for d in 1:depth
    @test all(x -> x == "RandomUnitary", first.(circuit[depth]))
  end

  circuit = randomcircuit(N;depth = depth, twoqubitgates="CX")
  @test PastaQ.nlayers(circuit) == depth
  for d in 1:depth
    @test all(x -> x == "CX", first.(circuit[depth]))
  end

  circuit = randomcircuit(N; depth = depth, twoqubitgates="CX", onequbitgates="Rn")
  @test size(circuit, 1) == depth

  circuit = randomcircuit(N; depth = depth, twoqubitgates="RandomUnitary", onequbitgates=["Rn", "X"]
  )
  @test size(circuit, 1) == depth

  Lx = 5
  Ly = 5
  circuit = randomcircuit((Lx,Ly); depth = depth, twoqubitgates="RandomUnitary", onequbitgates=["Rn", "X"])
  @test nqubits(circuit) == Lx*Ly

  Lx = 5
  Ly = 5
  circuit = randomcircuit((Lx,Ly); depth = depth,  twoqubitgates="RandomUnitary", onequbitgates=["Rn", "X"], rotated = true)
  @test nqubits(circuit) == Lx*Ly
end

@testset "dag circuit" begin

  N = 2
  depth = 4

  circuit = randomcircuit(N; depth = depth, twoqubitgates = "CX", onequbitgates = "Rn")
  U = PastaQ.array(runcircuit(circuit; process = true))
  dagcircuit = dag(circuit)
  V = PastaQ.array(runcircuit(dagcircuit; process = true))
  @test U ≈ V'

  circuit = randomcircuit(N; depth = depth, twoqubitgates = "RandomUnitary", onequbitgates = "Rn", layered = false)
  U = PastaQ.array(runcircuit(circuit; process = true))
  dagcircuit = dag(circuit)
  V = PastaQ.array(runcircuit(dagcircuit; process = true))
  @test U ≈ V'

  dagdagcircuit = dag(dagcircuit)
  Up = PastaQ.array(runcircuit(dagdagcircuit; process = true))
  @test U ≈ Up

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

error("Expected exception while running include")

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using DataFrames
using Compat

include("IndexedDF.jl")
typealias SparseMatrix SparseMatrixCSC{Float64, Int64}

export RelationData, addRelation!
export Relation, numData, numTest, assignToTest!, setTest!
export Entity, toStr, normalizeFeatures!, normalizeRows!
export EntityModel
export load_mf1c
export setPrecision!

type EntityModel
  sample::Matrix{Float64}  ## latent vectors (each column is one instance)

  mu    ::Vector{Float64}  ## mean
  Lambda::Matrix{Float64}  ## Precision matrix
  beta  ::Matrix{Float64}  ## parameter linking features to latent

  mu0   ::Vector{Float64}  ## Hyper-prior mean for NormalWishart
  b0    ::Float64          ## Hyper-prior for NormalWishart
  WI    ::Matrix{Float64}  ## Hyper-prior for NormalWishart (inverse of W)
  nu0   ::Float64          ## Hyper-prior for NormalWishart (degrees of freedom)

  uhat  ::Matrix{Float64}  ## temporary variable storing uhat (beta * feat)

  EntityModel() = new()
  EntityModel(num_latent::Int, num_instances::Int) = new(
    zeros(num_latent, num_instances), ## sample
    zeros(num_latent),    ## mu
    5*eye(num_latent),    ## Lambda
    zeros(0, num_latent), ## beta
    zeros(num_latent),    ## mu0
    2.0,                  ## b0
    eye(num_latent),      ## WI
    num_latent,           ## nu0
    zeros(0,0)            ## uhat
  )
end

@compat type Entity{FT,R}
  F::FT
  FF
  use_FF::Bool
  Frefs::Vector{Future}
  relations::Vector{R}
  count::Int64
  name::AbstractString

  modes::Vector{Int}
  modes_other::Vector{Vector{Int}}

  lambda_beta::Float64
  lambda_beta_sample::Bool
  mu::Float64   ## Hyper-prior for lambda_beta
  nu::Float64   ## Hyper-prior for lambda_beta

  model::EntityModel
  @compat Entity(F, relations::Vector{R}, count::Int64, name::AbstractString, lb::Float64=1.0, lb_sample::Bool=true, mu=1.0, nu=1e-3) = new(F, zeros(0,0), false, Future[], relations, count, name, Int[], Vector{Int}[], lb, lb_sample, mu, nu)
end

Entity(name::AbstractString; F=zeros(0,0), lambda_beta=1.0) = Entity{Any,Relation}(F::Any, Relation[], 0, name, lambda_beta)

## initializes the model parameters
function initModel!(entity::Entity, num_latent::Int64; lambda_beta::Float64 = NaN)
  m = EntityModel()
  entity.model = m

  m.sample = zeros(num_latent, entity.count)
  m.mu     = zeros(num_latent)
  m.Lambda = 5 * eye(num_latent)
  if hasFeatures(entity)
    m.beta = zeros( size(entity.F, 2), num_latent )
    m.uhat = zeros(num_latent, entity.count)
  else
    m.beta = zeros( 0, num_latent )
  end

  m.mu0    = zeros(num_latent)
  m.b0     = 2.0
  m.WI     = eye(num_latent)
  m.nu0    = num_latent

  if ! isnan(lambda_beta)
    entity.lambda_beta = lambda_beta
  end

  return nothing
end

hasFeatures(entity::Entity) = ! isempty(entity.F)

function toStr(en::Entity)
  if ! isdefined(en, :model)
    return string(en.name[1:min(3,end)], "[]")
  end
  return string(
    en.name[1:min(3,end)],
    "[",
       @sprintf("U:%6.2f", vecnorm(en.model.sample)),
       hasFeatures(en) ? @sprintf(" β:%3.2f", vecnorm(en.model.beta)) :"",
       hasFeatures(en) && en.lambda_beta_sample ? @sprintf(" λ=%1.1f", en.lambda_beta) :"",
    "]")
end

type RelationModel
  alpha_sample::Bool
  alpha_nu0::Float64
  alpha_lambda0::Float64

  lambda_beta::Float64

  alpha::Float64
  beta::Vector{Float64}
  mean_value::Float64
end

RelationModel(alpha::Float64=1.0, lambda_beta::Float64=1.0) = RelationModel(false, 2.0, 1.0, lambda_beta, alpha, zeros(0), 0.0)

type RelationTemp
  linear_values::Vector{Float64} ## mean_value + F * beta
  FF::Matrix{Float64}

  RelationTemp() = new()
end

type Relation
  data::IndexedDF
  F
  entities::Vector{Entity}
  name::AbstractString

  test_vec::DataFrame
  test_F
  test_label::Vector{Bool}
  class_cut::Float64

  model::RelationModel
  temp::RelationTemp

  Relation(data::IndexedDF, name::AbstractString, class_cut=0.0, alpha=1.0) = new(data, (), Entity[], name, data.df[Int[],:], (), Bool[], class_cut, RelationModel(alpha))
  ## Relation with already setup entities
  function Relation(data::DataFrame, name::AbstractString, entities=Entity[]; class_cut=0.0,
                    dims = Int64[maximum(data[:,i]) for i in 1 : size(data,2)-1])
    if ! isempty(entities)
      size(data, 2) - 1 == length(entities) || throw(ArgumentError("data has $(size(data,2)) columns but needs to have $(length(entities)) which is number of entities + 1"))
      for i in 1:length(entities)
        en = entities[i]
        ## update the counts of the entity
        if en.count == 0
          en.count = dims[i]
        elseif en.count > dims[i]
          dims[i] = en.count
        elseif en.count < dims[i]
          throw(ArgumentError("Entity $(en.name) has smaller count $(en.count) than the largest id in the data $(dims[i]). Set entity.count manually before creating the relation."))
        end
      end
    end
    return new(IndexedDF(data, dims), (), entities, name, data[Int[],:], (), Bool[], class_cut, RelationModel())
  end
end

## creating Relation from sparse matrix
function Relation(data::SparseMatrixCSC, name::AbstractString, entities=Entity[]; class_cut=0.0)
  length(entities) == 2 || throw(ArgumentError("For matrix relation the number of entities has to be 2."))
  U, V, X = findnz(data)
  dims = Int[size(data,1), size(data,2)]
  df = DataFrame(E1=U, E2=V, values=X)
  return Relation(df, name, entities, class_cut = class_cut, dims = dims)
end

import Base.size
size(r::Relation) = tuple([length(x) for x in r.data.index]...)
size(r::Relation, d::Int) = length(r.data.index[d])
numData(r::Relation) = nnz(r.data)
numTest(r::Relation) = size(r.test_vec, 1)
hasFeatures(r::Relation) = ! isempty(r.F)
function toStr(r::Relation)
  if ! isdefined(r, :model)
    return string(r.name[1:min(4,end)], "[]")
  end
  return string(
    r.name[1:min(4,end)],
    "[",
       @sprintf("α=%2.1f", r.model.alpha),
       hasFeatures(r) ? @sprintf(" β:%2.1f", vecnorm(r.model.beta)) :"",
    "]")
end

function assignToTest!(r::Relation, ntest::Int64)
  test_id  = sample(1:size(r.data.df,1), ntest; replace=false)
  assignToTest!(r, test_id)
end

function setPrecision!(r::Relation, precision)
  r.model.alpha = precision
end

function assignToTest!(r::Relation, test_id::Vector{Int64})
  test_vec = r.data.df[test_id, :]
  r.data   = removeSamples(r.data, test_id)
  r.test_vec    = test_vec
  r.test_label  = r.test_vec[:,end] .< r.class_cut
  if hasFeatures(r)
    r.test_F = r.F[test_id,:]
    train    = ones(Bool, size(r.F, 1) )
    train[test_id] = false
    r.F      = r.F[train,:]
  end
  nothing
end

function setTest!(r::Relation, test_df::DataFrame, test_feat = false)
  if hasFeatures(r) && test_feat == false
    throw(ArgumentError("Relation has features, please supply features with test data:\nsetTest(rel, test_df, test_features"))
  end
  if hasFeatures(r) && size(r.F, 2) != size(test_feat, 2)
    throw(ArgumentError("The test_feat must have the same number of columns as relation.F."))
  end
  if hasFeatures(r) && size(test_f, 1) != size(test_feat, 1)
    throw(ArgumentError("The test_feat must have the same number of rows as test_df."))
  end
  if size(test_df, 2) != size(r.data.df, 2)
    throw(ArgumentError("The number of columns in test_df must be the same as in relation.data.df."))
  end
  r.test_vec   = test_df
  r.test_label = r.test_vec[:,end] .< r.class_cut
  if hasFeatures(r)
    r.test_F = test_feat
  end
  nothing
end

function setTest!(r::Relation, test_mat::SparseMatrixCSC)
  if hasFeatures(r)
    throw(ArgumentError("Cannot add test set using SparseMatrixCSC when relation has features. Use DataFrame instead."))
  end
  if size(r.data.df, 2) != 3
    throw(ArgumentError("Relation must have 2 entities if using SparseMatrixCSC for test set."))
  end
  test_df = DataFrame()
  tnz = findnz(test_mat)

  test_df[ names(r.data.df)[1] ] = tnz[1]
  test_df[ names(r.data.df)[2] ] = tnz[2]
  test_df[ names(r.data.df)[3] ] = tnz[3]

  r.test_vec   = test_df
  r.test_label = r.test_vec[:,end] .< r.class_cut
  nothing
end

type RelationData
  entities::Vector{Entity}
  relations::Vector{Relation}

  RelationData() = new( Entity[], Relation[] )

  function RelationData(Am::IndexedDF; feat1=zeros(0,0), feat2=zeros(0,0), entity1="E1", entity2="E2", relation="Rel", ntest=0, class_cut=log10(200), alpha=5.0, alpha_sample=false, lambda_beta=1.0)
    r  = alpha_sample ?Relation(Am, relation, class_cut) :Relation(Am, relation, class_cut, alpha)
    e1 = Entity{isempty(feat1) ? Any :typeof(feat1), Relation}( feat1, [r], size(r,1), entity1, lambda_beta )
    e2 = Entity{isempty(feat2) ? Any :typeof(feat2), Relation}( feat2, [r], size(r,2), entity2, lambda_beta )
    if ! isempty(feat1) && size(feat1,1) != size(r,1)
      throw(ArgumentError("Number of rows in feat1 $(size(feat1,1)) must equal number of rows in the relation $(size(Am,1))"))
    end
    if ! isempty(feat2) && size(feat2,1) != size(Am,2)
      throw(ArgumentError("Number of rows in feat2 $(size(feat2,1)) must equal number of columns in the relation $(size(Am,2))"))
    end
    push!(r.entities, e1)
    push!(r.entities, e2)
    return new( Entity[e1, e2], Relation[r] )
  end
end

function RelationData(Am::DataFrame; rname="R1", class_cut=log10(200), alpha=5.0)
  dims = Int64[maximum(Am[:,i]) for i in 1 : size(Am,2)-1]
  idf  = IndexedDF(Am, dims)
  rd   = RelationData()
  push!(rd.relations, Relation(idf, rname, class_cut, alpha))
  for d in 1:length(dims)
    en = Entity{Any, Relation}( zeros(0,0), [rd.relations[1]], size(idf,d), string(names(Am)[d]))
    push!(rd.entities, en)
    push!(rd.relations[1].entities, en)
  end
  return rd
end

function rep_int(x, times)
  out = zeros(eltype(x), sum(times))
  idx = 1
  for i in 1:size(x,1)
    out[idx : idx+times[i]-1] = x[i]
    idx += times[i]
  end
  return out
end

function RelationData(M::SparseMatrixCSC{Float64,Int64}; kw...)
  dims = size(M)
  cols = rep_int(1:size(M,2), M.colptr[2:end] - M.colptr[1:end-1])
  df   = DataFrame( row=M.rowval, col=cols, value=nonzeros(M) )
  idf  = IndexedDF(df, dims)
  return RelationData(idf; kw...)
end

function RelationData(r::Relation)
  rd = RelationData()
  addRelation!(rd, r)
  return rd
end

## computes F * beta
function F_mul_beta{F}(en::Entity{F,Relation})
  if ! isempty(en.Frefs)
    return Frefs_mul_B(en.Frefs, en.model.beta)
  else
    return en.F * en.model.beta
  end
end

## computes F' * beta
function Ft_mul_B{F}(en::Entity{F,Relation}, B::Matrix{Float64})
  if ! isempty(en.Frefs)
    return Frefs_t_mul_B(en.Frefs, B)
  else
    return At_mul_B(en.F, B)
  end
end

function reset!(data::RelationData, num_latent; lambda_beta=NaN, compute_ff_size = 6500, cg_pids=Int[myid()])
  for en in data.entities
    initModel!(en, num_latent, lambda_beta = lambda_beta)
    en.modes = Int64[ find(en2 -> en2 == en, r.entities)[1] for r in en.relations ]
    en.modes_other = Vector{Int64}[ find(en2 -> en2 != en, r.entities) for r in en.relations ]
    if hasFeatures(en)
      if size(en.F, 2) <= compute_ff_size
        en.FF = full(At_mul_B(en.F, en.F))
        en.use_FF = true
      else
        ## setup CG on julia threads
        init_Frefs!(en, num_latent, cg_pids)
        en.use_FF = false
      end
    end
  end
  for r in data.relations
    r.model.mean_value    = valueMean(r.data)
    r.temp = RelationTemp()
    if hasFeatures(r) && size(r.F,2) <= compute_ff_size
      r.temp.linear_values = r.model.mean_value * ones(numData(r))
      r.temp.FF = full(r.F' * r.F)
    end
  end
end

function init_Frefs!(en::Entity, num_latent::Int, pids::Vector{Int})
  Fns = nonshared(en.F)
  empty!(en.Frefs)
  for i = 1:min(length(pids), num_latent)
    pid = pids[i]
    push!(en.Frefs, @spawnat pid copyto(Fns, Int[pid]))
  end
end

function init_Frefs!(en::Entity, num_latent::Int, pids::Vector{Vector{Int}})
  Fns = nonshared(en.F)
  empty!(en.Frefs)
  for i = 1:min(length(pids), num_latent)
    pid   = pids[i][1]
    mulpids = pids[i][2:end]
    isempty(mulpids) && push!(mulpids, pid)
    push!(en.Frefs, @spawnat pid copyto(Fns, mulpids))
  end
end

function init_Frefs!(en::Entity, num_latent::Int, pids::Dict{Int, Vector{Int}})
  Fns = nonshared(en.F)
  empty!(en.Frefs)
  for e in pids
    pid     = e[1]
    mulpids = e[2]
    push!(en.Frefs, @spawnat pid copyto(Fns, mulpids))
  end
end

function addRelation!(rd::RelationData, r::Relation)
  if length(size(r)) != length(r.entities)
    throw(ArgumentError("Relation has $(length(r.entities)) entities but its data implies $(size(r))."))
  end
  push!(rd.relations, r)
  ## adding entities
  for i in 1:length(r.entities)
    en = r.entities[i]
    if en.count == 0
      ## updating entity count
      en.count = size(r, i)
    elseif en.count != size(r, i)
      throw(ArgumentError("Entity $(en.name) has $(en.count) instances, relation $(r.name) has data for $(size(r.entities, i))."))
    end
    if ! any(rd.entities .== en)
      push!(rd.entities, en)
    end
    if ! any(en.relations .!= r)
      push!(en.relations, r)
    end
  end
  return nothing
end

import Base.show
function show(io::IO, rd::RelationData)
  println(io, "[Relations]")
  for r in rd.relations
    @printf(io, "%10s: %s, #known = %d, #test = %d, α = %s", r.name, join([e.name for e in r.entities], "--"), numData(r), numTest(r), r.model.alpha_sample ?"sample" :@sprintf("%.2f", r.model.alpha))
    hasFeatures(r) && @printf(io, ", #feat = %d", size(r.F,2))
    @printf(io, "\n")
  end

  println(io, "[Entities]")
  for en in rd.entities
    @printf(io, "%10s: %6d ", en.name, en.count)
    if hasFeatures(en)
      print(io, "with ", size(en.F,2), " features (λ = ",
                en.lambda_beta_sample ? "sample" :@sprintf("%1.1f", en.lambda_beta),")")
    else
      print(io, "with no features")
    end
    println(io)
  end
end

function show(io::IO, r::Relation)
  print(io, "[Relation]")
  @printf(io, " %s: %s, #known = %d, #test = %d, α = %s", r.name, join([e.name for e in r.entities], "--"), numData(r), numTest(r), r.model.alpha_sample ?"sample" :@sprintf("%.2f", r.model.alpha))
  hasFeatures(r) && @printf(io, ", #feat = %d", size(r.F,2))
end

function show(io::IO, en::Entity)
  print(io, "[Entity]")
  @printf(io, " %s: %6d ", en.name, en.count)
  if hasFeatures(en)
    print(io, "with ", size(en.F,2), " features (λ = ",
              en.lambda_beta_sample ? "sample" :@sprintf("%1.1f", en.lambda_beta),")")
  else
    print(io, "with no features")
  end
end

function normalizeFeatures!(entity::Entity)
  diagsq  = sqrt(vec( sum(entity.F .^ 2,1) ))
  entity.F  = entity.F * spdiagm(1.0 ./ diagsq)
  return
end

function normalizeRows!(entity::Entity)
  diagf    = sqrt(vec( sum(entity.F.^2, 2) ))
  entity.F = spdiagm( 1.0 ./ diagf ) * entity.F
end

function load_mf1c(;ic50_file     = "chembl_19_mf1c/chembl-IC50-346targets.csv",
                   cmp_feat_file  = "chembl_19_mf1c/chembl-IC50-compound-feat.csv",
                   normalize_feat = false,
                   alpha_sample   = false)
  ## reading IC50 matrix
  X = readtable(ic50_file, header=true)
  rename!(X, [:row, :col], [:compound, :target])

  dims = [maximum(X[:compound]), maximum(X[:target])]

  X[:, :value] = log10(X[:, :value]) + 1e-5
  idx          = sample(1:size(X,1), round(Int, floor(20/100 * size(X,1))); replace=false)
  probe_vec    = X[idx,:]
  X            = X[setdiff(1:size(X,1), idx), :]

  ## reading feature matrix
  feat = readtable(cmp_feat_file, header=true)
  F    = sparse(feat[:compound], feat[:feature], 1.0)

  #Am = sparse( X[:compound], X[:target], X[:value])
  Xi = IndexedDF(X, dims)

  ## creating data object
  data = RelationData(Xi, feat1 = F, alpha_sample = alpha_sample)
  data.relations[1].test_vec    = probe_vec
  data.relations[1].test_label  = data.relations[1].test_vec[:,3] .< log10(200)

  if normalize_feat != false
    if normalize_feat == "rows"
      normalizeRows!(data.entities[1])
    else
      normalizeFeatures!(data.entities[1])
    end
  end

  return data
end


text:
using Base: @propagate_inbounds

using OffsetArrays
using JLD2

using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.Fields
using Oceananigans.Fields: show_location

import Oceananigans: short_show

struct FieldTimeSeries{X, Y, Z, K, A, T, N, D, G, B, χ} <: AbstractDataField{X, Y, Z, A, G, T, N}
                   data :: D
           architecture :: A
                   grid :: G
    boundary_conditions :: B
                  times :: χ
                   name :: String
               filepath :: String

    function FieldTimeSeries{X, Y, Z}(backend::K, data::D, arch::A, grid::G, bcs::B, times::χ, name, filepath, N) where {X, Y, Z, K, D, A, G, B, χ}
        T = eltype(grid)
        return new{X, Y, Z, K, A, T, N, D, G, B, χ}(data, arch, grid, bcs, times, name, filepath)
    end
end

# Include the time dimension.
@inline Base.size(fts::FieldTimeSeries) = (size(location(fts), fts.grid)..., length(fts.times))

@propagate_inbounds Base.getindex(f::FieldTimeSeries{X, Y, Z, InMemory}, i, j, k, n) where {X, Y, Z} = f.data[i, j, k, n]

"""
    FieldTimeSeries(filepath, name; architecture=CPU(), backend=InMemory())

Returns a `FieldTimeSeries` for the field `name` describing a field's time history from a JLD2 file
located at `filepath`. Note that model output must have been saved with halos. The `InMemory` backend
will store the data fully in memory as a 4D multi-dimensional array while the `OnDisk` backend will
lazily load field time snapshots when the `FieldTimeSeries` is indexed linearly.
"""
FieldTimeSeries(filepath, name; architecture=CPU(), backend=InMemory()) =
    FieldTimeSeries(filepath, name, architecture, backend)

function FieldTimeSeries(filepath, name, architecture, backend::InMemory)
    file = jldopen(filepath)

    grid = file["serialized/grid"]
    Hx, Hy, Hz = halo_size(grid)

    iterations = parse.(Int, keys(file["timeseries/t"]))
    times = [file["timeseries/t/$i"] for i in iterations]

    LX, LY, LZ = location = file["timeseries/$name/serialized/location"]

    Nt = length(times)
    data_size = size(file["timeseries/$name/0"])

    ArrayType = array_type(architecture)
    raw_data = zeros(data_size..., Nt) |> ArrayType
    data = offset_data(raw_data, grid, location)

    for (n, iter) in enumerate(iterations)
        data.parent[:, :, :, n] .= file["timeseries/$name/$iter"] |> ArrayType
    end

    bcs = file["timeseries/$name/serialized/boundary_conditions"]

    close(file)

    return FieldTimeSeries{LX, LY, LZ}(backend, data, architecture, grid, bcs, times, name, abspath(filepath), ndims(data))
end

function FieldTimeSeries(filepath, name, architecture, backend::OnDisk)
    file = jldopen(filepath)

    grid = file["serialized/grid"]
    iterations = parse.(Int, keys(file["timeseries/t"]))
    times = [file["timeseries/t/$i"] for i in iterations]

    data = nothing
    LX, LY, LZ = file["timeseries/$name/serialized/location"]
    bcs = file["timeseries/$name/serialized/boundary_conditions"]

    close(file)

    return FieldTimeSeries{LX, LY, LZ}(backend, data, architecture, grid, bcs, times, name, abspath(filepath), 4)
end

Base.getindex(fts::FieldTimeSeries{X, Y, Z, InMemory}, n::Int) where {X, Y, Z} =
    Field((X, Y, Z), fts.architecture, fts.grid, fts.boundary_conditions, fts.data[:, :, :, n])

function Base.getindex(fts::FieldTimeSeries{X, Y, Z, OnDisk}, n::Int) where {X, Y, Z}
    file = jldopen(fts.filepath)
    iter = keys(file["timeseries/t"])[n]
    raw_data = file["timeseries/$(fts.name)/$iter"] |> array_type(fts.architecture)
    close(file)

    loc = (X, Y, Z)
    field_data = offset_data(raw_data, fts.grid, loc)
    return Field(loc, fts.architecture, fts.grid, fts.boundary_conditions, field_data)
end

backend_str(::InMemory) = "InMemory"
backend_str(::OnDisk) = "OnDisk"

short_show(fts::FieldTimeSeries{X, Y, Z, K}) where {X, Y, Z, K} =
    string("$(join(size(fts), "×")) FieldTimeSeries{$(backend_str(K()))} located at $(show_location(fts))")

Base.show(io::IO, fts::FieldTimeSeries{X, Y, Z, K, A}) where {X, Y, Z, K, A} =
    print(io, "$(short_show(fts))\n",
          "├── filepath: $(fts.filepath)\n",
          "├── architecture: $A\n",
          "└── grid: $(short_show(fts.grid))")


text:
# CUDA extensions to the C language

# TODO: "CUDA C programming guide" > "C language extensions" lists mathematical functions,
#       without mentioning libdevice. Is this implied, by NVCC always using libdevice,
#       or are there some natively-supported math functions as well?

# yes: https://docs.nvidia.com/cuda/cuda-math-api/group__CUDA__MATH__INTRINSIC__DOUBLE.html
# see /home/tbesard/CUDA/toolkit/current/include/sm_20_intrinsics.h

include(joinpath("cuda_intrinsics", "memory_shared.jl"))
include(joinpath("cuda_intrinsics", "indexing.jl"))
include(joinpath("cuda_intrinsics", "synchronization.jl"))
include(joinpath("cuda_intrinsics", "warp_vote.jl"))
include(joinpath("cuda_intrinsics", "warp_shuffle.jl"))
include(joinpath("cuda_intrinsics", "output.jl"))
include(joinpath("cuda_intrinsics", "assertion.jl"))
include(joinpath("cuda_intrinsics", "memory_dynamic.jl"))

max_stars_repo_name: Cyvadra/CUDAnative

max_issues_repo_name: Cyvadra/CUDAnative

max_forks_repo_name: Cyvadra/CUDAnative


text:
# Structs and constants related to the codes_handle type

export Message, getbytes, writemessage, missingvalue, clone, data, maskedvalues

mutable struct codes_handle
end

mutable struct Message
    ptr::Ptr{codes_handle}
end

"""
    Message(index::Index)
    Message(file::GribFile)

Retrieve the next message.

Returns `nothing` if there are no more messages.
"""
function Message(index::Index)::Union{Message, Nothing}
    eref = Ref(Int32(0))
    ptr = ccall((:codes_handle_new_from_index, eccodes), Ptr{codes_handle},
                (Ptr{codes_index}, Ref{Cint}), index.ptr, eref)
    if eref[] != 0 && eref[] != CODES_END_OF_INDEX
        throw(ErrorException("Error: $(errors[-eref[]+1])"))
    end
    if ptr == C_NULL || eref[] == CODES_END_OF_INDEX
        return nothing
    else
        m = Message(ptr)
        finalizer(destroy, m)
        return m
    end
end

function Message(file::GribFile)::Union{Message, Nothing}
    eref = Ref(Int32(0))
    ptr = ccall((:codes_grib_handle_new_from_file, eccodes), Ptr{codes_handle},
                (Ptr{codes_context}, Ptr{File}, Ref{Cint}),
                C_NULL, file.ptr, eref)
    file.pos += 1
    errorcheck(eref[])
    if ptr == C_NULL
        return nothing
    else
        m = Message(ptr)
        finalizer(destroy, m)
        return m
    end
end

"""
    Message(raw::Vector{UInt8})
    Message(messages::Vector{Vector{UInt8}})

Create a Message from a bytes array or array of bytes arrays.

The Message will share its underlying data with the inputed data.
"""
function Message(raw::Vector{UInt8})
    ptr = ccall((:codes_handle_new_from_message, eccodes), Ptr{codes_handle},
                (Ptr{codes_context}, Ptr{UInt8}, Csize_t),
                C_NULL, raw, length(raw))
    if ptr == C_NULL
        throw(ErrorException("Invalid message, could not create Message"))
    end
    m = Message(ptr)
    finalizer(destroy, m)
    return m
end

function Message(messages::Vector{Vector{UInt8}})
    sizes = map(length, messages)
    eref = Ref(Int32(0))
    ptr = ccall((:codes_grib_handle_new_from_multi_message, eccodes), Ptr{codes_handle},
                (Ptr{codes_context}, Ptr{Ptr{UInt8}}, Ptr{Csize_t}, Ptr{Cint}),
                C_NULL, messages, sizes, eref)
    errorcheck(eref[])
    m = Message(ptr)
    finalizer(destroy, m)
    return m
end

"""
    getbytes(message::Message)

Get a coded representation of the Message.
"""
function getbytes(message::Message)
    lenref = Ref(Csize_t(0))
    vp = Ref{Ptr{UInt8}}(C_NULL)
    err = ccall((:codes_get_message, eccodes), Cint,
                (Ptr{codes_handle}, Ptr{Ptr{UInt8}}, Ref{Csize_t}),
                message.ptr, vp, lenref)
    errorcheck(err)
    return unsafe_wrap(Vector{UInt8}, vp[], lenref[])
end

"""
    writemessage(handle::Message, filename::AbstractString; mode="w")

Write the message respresented by `handle` to the file at `filename`.

`mode` is a mode as described by `Base.open`.
"""
function writemessage(handle::Message, filename::AbstractString; mode="w")
    err = ccall((:codes_write_message, eccodes), Cint, (Ptr{codes_handle}, Cstring, Cstring),
                handle.ptr, filename, mode)
    errorcheck(err)
end

function Base.haskey(msg::Message, key::AbstractString)
    try
        msg[key]
        return true
    catch e
        return false
    end
end

"""
    missingvalue(msg::Message)

Return the missing value of the message. If one isn't included returns 1e30.
"""
function missingvalue(msg::Message)
    if haskey(msg, "missingValue")
        return msg["missingValue"]
    else
        return 1e30
    end
end

"""
    maskedvalues(msg::Message)

Return the values of the message masked so the missing value is `missing`.
"""
function maskedvalues(msg::Message)
    missingval = missingvalue(msg)
    vals = msg["values"]
    masked = Array{Union{Float64, Missing}, 2}(undef, size(vals)...)
    for i in eachindex(vals)
        @inbounds masked[i] = vals[i] == missingval ? missing : vals[i]
    end
    return masked
end

""" Get the native type of a key. """
function getnativetype(handle::Message, key::AbstractString)::DataType
    typeref = Ref(Int32(0))
    err = ccall((:codes_get_native_type, eccodes), Cint, (Ptr{codes_handle}, Cstring, Ref{Cint}),
                handle.ptr, key, typeref)
    errorcheck(err)
    typeint = typeref[]

    if typeint == CODES_TYPE_LONG
        return Clong
    elseif typeint == CODES_TYPE_STRING
        return String
    elseif typeint == CODES_TYPE_DOUBLE
        return Float64
    elseif typeint == CODES_TYPE_BYTES
        return Vector{UInt8}
    elseif typeint == CODES_TYPE_MISSING
        return Missing
    else
        return Tuple{}
    end
end

""" Get the number of values in a key. """
function getsize(message::Message, key::AbstractString)
    lenref = Ref(Csize_t(0))
    err = ccall((:codes_get_size, eccodes), Cint, (Ptr{codes_handle}, Cstring, Ref{Csize_t}),
                 message.ptr, key, lenref)
    errorcheck(err)
    return lenref[]
end

""" Get the value of a key. """
function Base.getindex(handle::Message, key::AbstractString)
    typ = getnativetype(handle, key)
    len = getsize(handle, key)
    if typ == String
        bufr = UInt(1024)
        bufref = Ref(bufr)
        mesg = zeros(Cuchar, bufr)
        err = ccall((:codes_get_string, eccodes), Cint,
                    (Ptr{codes_handle}, Cstring, Ptr{Cuchar}, Ref{Csize_t}),
                    handle.ptr, key, mesg, bufref)

        mesgstr = rstrip(transcode(String, mesg), '\0')
        errorcheck(err)
        return mesgstr
    elseif typ == Float64
        if len == 1
            valref = Ref(Float64(0))
            err = ccall((:codes_get_double, eccodes), Cint,
                        (Ptr{codes_handle}, Cstring, Ref{Cdouble}),
                        handle.ptr, key, valref)
            errorcheck(err)
            return valref[]
        else
            lenref = Ref(len)
            vals = Vector{Float64}(undef, len)
            err = ccall((:codes_get_double_array, eccodes), Cint,
                       (Ptr{codes_handle}, Cstring, Ref{Cdouble}, Ref{Csize_t}),
                        handle.ptr, key, vals, lenref)
            errorcheck(err)
            if key == "values"
                nj = handle["Nj"]
                ni = handle["Ni"]
                columnmajor = handle["jPointsAreConsecutive"] != 0
                if columnmajor
                    return reshape(vals, nj, ni)
                else
                    return reshape(vals, ni, nj)
                end
            else
                return vals
            end
        end
    elseif typ == Clong
        if len == 1
            valref = Ref(Clong(0))
            err = ccall((:codes_get_long, eccodes), Cint, (Ptr{codes_handle}, Cstring, Ref{Clong}),
                        handle.ptr, key, valref)
            errorcheck(err)
            return valref[]
        else
            lenref = Ref(len)
            vals = Vector{Clong}(undef, len)
            err = ccall((:codes_get_long_array, eccodes), Cint,
                       (Ptr{codes_handle}, Cstring, Ref{Clong}, Ref{Csize_t}),
                        handle.ptr, key, vals, lenref)
            errorcheck(err)
            return vals
        end
    elseif typ == Vector{UInt8}
        bufr = UInt(1024)
        bufref = Ref(bufr)
        mesg = Vector{UInt8}(undef, bufr)
        err = ccall((:codes_get_bytes, eccodes), Cint,
                    (Ptr{codes_handle}, Cstring, Ptr{Cuchar}, Ref{Csize_t}),
                    handle.ptr, key, mesg, bufref)
        errorcheck(err)
        return mesg[1:bufref[]]
    else
        throw(ErrorException("Could not get key: unsupported type"))
    end
end

"""
    setindex!(message::Message, value, key::AbstractString)

Set the value of a key.
"""
function Base.setindex!(handle::Message, value::Integer, key::AbstractString)
    err = ccall((:codes_set_long, eccodes), Cint, (Ptr{codes_handle}, Cstring, Clong),
                 handle.ptr, key, value)
    errorcheck(err)
end

function Base.setindex!(handle::Message, value::AbstractFloat, key::AbstractString)
    err = ccall((:codes_set_double, eccodes), Cint, (Ptr{codes_handle}, Cstring, Cdouble),
                 handle.ptr, key, value)
    errorcheck(err)
end

function Base.setindex!(handle::Message, value::AbstractString, key::AbstractString)
    lenref = Ref(Csize_t(length(value)))
    err = ccall((:codes_set_string, eccodes), Cint,
                (Ptr{codes_handle}, Cstring, Cstring, Ref{Csize_t}),
                 handle.ptr, key, value, lenref)
    errorcheck(err)
end

function Base.setindex!(handle::Message, value::Vector{UInt8}, key::AbstractString)
    lenref = Ref(Csize_t(length(value)))
    err = ccall((:codes_set_bytes, eccodes), Cint,
                (Ptr{codes_handle}, Cstring, Ref{UInt8}, Ref{Csize_t}),
                 handle.ptr, key, value, lenref)
    errorcheck(err)
end

function Base.setindex!(handle::Message, value::Array{Float64}, key::AbstractString)
    len = Csize_t(length(value))
    err = ccall((:codes_set_double_array, eccodes), Cint,
                (Ptr{codes_handle}, Cstring, Ref{Cdouble}, Csize_t),
                handle.ptr, key, value, len)
    errorcheck(err)
end

function Base.setindex!(handle::Message, value::Array{Clong}, key::AbstractString)
    len = Csize_t(length(value))
    err = ccall((:codes_set_long_array, eccodes), Cint,
                (Ptr{codes_handle}, Cstring, Ref{Clong}, Csize_t),
                 handle.ptr, key, value, len)
    errorcheck(err)
end

function Base.setindex!(handle::Message, value::Array{String}, key::AbstractString)
    len = Csize_t(length(value))
    err = ccall((:codes_set_string_array, eccodes), Cint,
                (Ptr{codes_handle}, Cstring, Ref{Cstring}, Csize_t),
                 handle.ptr, key, value, len)
    errorcheck(err)
end

""" Duplicate a message. """
function clone(handle::Message)::Union{Message, Nothing}
    ptr = ccall((:codes_handle_clone, eccodes), Ptr{codes_handle}, (Ptr{codes_handle},), handle.ptr)
    if ptr == C_NULL
        return nothing
    else
        m = Message(ptr)
        finalizer(destroy, m)
        return m
    end
end

# Aliases for clone for more Julian style
Message(handle::Message) = clone(handle)
Base.deepcopy(handle::Message) = clone(handle)

"""
    data(handle::Message)

Retrieve the longitudes, latitudes, and values from the message.

# Example
```
GribFile(filename) do f
    msg = Message(f)
    lons, lats, values = data(msg)
end
```
"""
function data(handle::Message)::Tuple{Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}}
    # Allocate arrays
    npoints = handle["numberOfPoints"]
    lons = Vector{Float64}(undef, npoints)
    lats = Vector{Float64}(undef, npoints)
    values = Vector{Float64}(undef, npoints)
    err = ccall((:codes_grib_get_data, eccodes), Cint,
                (Ptr{codes_handle}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
                handle.ptr, lats, lons, values)
    errorcheck(err)

    nj = handle["Nj"]
    ni = handle["Ni"]
    columnmajor = handle["jPointsAreConsecutive"] != 0
    if columnmajor
        return reshape(lons, nj, ni), reshape(lats, nj, ni), reshape(values, nj, ni)
    else
        return reshape(lons, ni, nj), reshape(lats, ni, nj), reshape(values, ni, nj)
    end
end

""" Safely destroy a message handle. """
function destroy(handle::Message)
    if handle.ptr != C_NULL
        err = ccall((:codes_handle_delete, eccodes), Cint, (Ptr{codes_handle},), handle.ptr)
        errorcheck(err)
        handle.ptr == C_NULL
    end
end

# Functions related to printing
function Base.show(io::IO, mime::MIME"text/plain", message::Message)
    dispkeys = ["date", "gridType", "stepRange", "typeOfLevel", "level",
                 "shortName", "name"]
    dispkeys = [k for k in dispkeys if haskey(message, k)]
    offsets = [9, 15, 10, 18, 6, 10, 4]
    line1 = strip(join([rpad(k, offsets[i]) for (i, k) in enumerate(dispkeys)])) * "\n"
    line2 = strip(join([rpad(message[k], offsets[i]) for (i, k) in enumerate(dispkeys)]))
    print(io, line1 * line2)
end

function Base.show(io::IO, message::Message)
    bytes = getbytes(message)
    str = bytes2hex(bytes)
    print(io, str)
end


text:
using NumericalMethodsforEngineers
using Test

@sym begin
  ClearAll(xi, yi, Y, Ydot, Ydotdot, ytilde1, ytilde2, R)
  xi = [0, 1//3, 2//3, 1]
  yi = [0, a, b, 1]
  Y(x_) := LagrangePolynomial(xi, yi)
  #
  # Can be re-formulated as Y(x) = F(x) + C1(a) * Ψ1(x) + C2(a) * Ψ2(x)
  #
  # F(x_) := (1/2)*(9x^3 - 9x^2 + 2x)
  # C1(a_) := (27/2)*(a-b)
  # Ψ1(x_) := x^3 - x^2
  # C2(a_) := -(9/2)*(2a-b)
  # Ψ2(x_) := x^2 - x
  #
  # Y(x_) := F(x) + C1(a)*Ψ1(x) + C2(a)*Ψ2(x)
  #
  R(x_) = Simplify(D(Y(x), x, 2) - 3*x - 4*Y(x))
  sol = Solve([R(xi[2]), R(xi[3])], [a, b])
  ytilde(x_) = Y(x) ./ Flatten(sol)
  SetJ(r, ToString(Simplify(R(x))))
  SetJ(s, "[$(sol[1][1][2]), $(sol[1][2][2])]")
  SetJ(t, ToString(Simplify(ytilde(x))))
end

println("\n\nExample 7.15b: y'' = 3x + 4y, y(0)=0, y(1)=1")
println("by 2-point collocation Weighted Residual Method")
@sym Println("\nY(x): ", Y(x), "\n")
@sym Println("R(x) = ", R(x), "\n")
println("(a, b) = $(s)\n")
@sym Println("ytilde_2pt_collocation(x) = ", Simplify(ytilde(x)), "\n")
println()

@eval rf_2pt_collocation_2(x, a, b) = $(Meta.parse(r))
@eval ytilde_2pt_collocation(x) = $(Meta.parse(t))
@eval (a, b) = $(Meta.parse(s))

@test t == "(1/806)*x*(239 - 252x + 819*x^2)"
@test rf_2pt_collocation_2(1//3, a, b) <= eps()
@test rf_2pt_collocation_2(2//3, a, b) <= eps()


text:
using CoilFunctions
using CoilFunctions: ideal_fill, AWG_Chart
using Unitful
using Test

@testset "CoilGeometry" begin
    @test CoilGeometry(10u"mm",20u"mm",30u"mm") == CoilGeometry(10.0u"mm",20.0u"mm",30.0u"mm")
    @test CoilGeometry(10u"mm",20u"mm",30u"mm", 1.5u"mm") == CoilGeometry(13u"mm",17.0u"mm",27u"mm")
    @test CoilGeometry(10u"mm",20u"mm",30u"mm") ≈ CoilGeometry(prevfloat(10.0)u"mm",nextfloat(20.0)u"mm",30.0u"mm")
end

@testset "Ideal Fill" begin
    let
        fillfactor=0.9
        awg=AWG_Chart[30]
        a=CoilGeometry(10u"mm",10u"mm"+10u"mm",30.0u"mm")
        b=CoilGeometry(10u"mm",10u"mm"+10u"mm"*sqrt(fillfactor),30.0u"mm"*sqrt(fillfactor))
        @test ideal_fill(a,awg) == CoilFunctions.WindingLayers(19, 101, 101, 1010, 909, 1919)
        @test ideal_fill(a,awg,fillfactor) == CoilFunctions.WindingLayers(18, 96, 96, 864, 864, 1728)
        @test ideal_fill(a,awg,fillfactor) == ideal_fill(b,awg)

        b=CoilGeometry(10u"mm",10u"mm"+2.01awg.d_insul,10.01awg.d_insul)
        @test ideal_fill(b,awg) == WindingLayers(1, 10, 9, 10, 0, 10)
    end
end

@testset "Winding Associativity" begin
    let
        g = (layers, odd, id, gauge) -> begin
            f = (layers, odd, even, id, gauge) -> begin
                coil = enclosewinding(layers, odd, even, id, gauge)
                winding = ideal_fill(coil, gauge)
                winding.layers == layers && winding.odd_layer == odd && winding.even_layer == even
            end
            all(f(layers,odd,odd,id,gauge) for layers in layers for odd in odd for id in id) &&
            all(f(layers,odd,odd-1,id,gauge) for layers in layers for odd in odd for id in id)
        end
        @test all(x->g(1:100,1:200,range(10u"mm";stop=120u"mm",length=17),x),values(AWG_Chart))
    end
end

@testset "Layer Associativity" begin
    let
        f = (n,gauge) -> begin
            d=gauge.d_insul
            w=CoilFunctions.oddwidth(n,d) |> u"mm"
            _n =CoilFunctions.oddturns(w,d)
            # @info (n,_n,d,w)
            n == _n
        end
        @test all(f(n,gauge) for n in 1:1000 for (key,gauge) in CoilFunctions.AWG_Chart)
    end
    let
        f = (n,gauge) -> begin
            d=gauge.d_insul
            w=CoilFunctions.evenwidth(n,d) |> u"mm"
            _n =CoilFunctions.eventurns(w,d)
            # @info (n,_n,d,w)
            n == _n
        end
        @test all(f(n,gauge) for n in 1:1000 for (key,gauge) in CoilFunctions.AWG_Chart)
    end
    let
        f = (n,gauge) -> begin
            d=gauge.d_insul
            h=CoilFunctions.layerheight(n,d) |> u"mm"
            _n =CoilFunctions.ideal_layers(h,d)
            # @info (n=n,n2=_n,d=d,h=h)
            n == _n
        end
        @test all(f(n,gauge) for n in 1:100 for (key,gauge) in CoilFunctions.AWG_Chart)
    end
end

@testset "Basic Calculations" begin
    @test CoilFunctions.resistance(1000u"ft"/π,CoilFunctions.AWG_Chart[30]) ≈ 102.7u"Ω"
end

@testset "Assertion Errors" begin
    @test_throws AssertionError enclosewinding(10,10,8,10u"mm",0.2u"mm")
    #TODO enclosedwinding from coil
end


text:
module TestIO
    using Base.Test
    using DataFrames, Compat

    #test_group("We can read various file types.")

    data = joinpath(dirname(@__FILE__), "data")

    filenames = ["$data/blanklines/blanklines.csv",
                 "$data/compressed/movies.csv.gz",
                 "$data/headeronly/headeronly.csv",
                 "$data/newlines/embedded_os9.csv",
                 "$data/newlines/embedded_osx.csv",
                 "$data/newlines/embedded_windows.csv",
                 "$data/padding/space_after_delimiter.csv",
                 "$data/padding/space_around_delimiter.csv",
                 "$data/padding/space_before_delimiter.csv",
                 "$data/quoting/empty.csv",
                 "$data/quoting/escaping.csv",
                 "$data/quoting/quotedcommas.csv",
                 "$data/scaling/10000rows.csv",
                 "$data/utf8/corrupt_utf8.csv",
                 "$data/utf8/short_corrupt_utf8.csv",
                 "$data/utf8/utf8.csv"]

    for filename in filenames
        try
            df = readtable(filename)
        catch
            error(@sprintf "Failed to read %s\n" filename)
        end
    end

    #test_group("We get the right size, types, values for a basic csv.")

    filename = "$data/scaling/movies.csv"
    df = readtable(filename)

    @test size(df) == (58788, 25)

    @test df[1, 1] === 1
    @test df[1, 2] == "\$"
    @test df[1, 3] === 1971
    @test df[1, 4] === 121
    @test df[1, 5] === NA
    @test df[1, 6] === 6.4
    @test df[1, 7] === 348
    @test df[1, 8] === 4.5
    @test df[1, 9] === 4.5
    @test df[1, 10] === 4.5
    @test df[1, 11] === 4.5
    @test df[1, 12] === 14.5
    @test df[1, 13] === 24.5
    @test df[1, 14] === 24.5
    @test df[1, 15] === 14.5
    @test df[1, 16] === 4.5
    @test df[1, 17] === 4.5
    @test df[1, 18] == ""
    @test df[1, 19] === 0
    @test df[1, 20] === 0
    @test df[1, 21] === 1
    @test df[1, 22] === 1
    @test df[1, 23] === 0
    @test df[1, 24] === 0
    @test df[1, 25] === 0

    @test df[end, 1] === 58788
    @test df[end, 2] == "xXx: State of the Union"
    @test df[end, 3] === 2005
    @test df[end, 4] === 101
    @test df[end, 5] === 87000000
    @test df[end, 6] === 3.9
    @test df[end, 7] === 1584
    @test df[end, 8] === 24.5
    @test df[end, 9] === 4.5
    @test df[end, 10] === 4.5
    @test df[end, 11] === 4.5
    @test df[end, 12] === 4.5
    @test df[end, 13] === 14.5
    @test df[end, 14] === 4.5
    @test df[end, 15] === 4.5
    @test df[end, 16] === 4.5
    @test df[end, 17] === 14.5
    @test df[end, 18] == "PG-13"
    @test df[end, 19] === 1
    @test df[end, 20] === 0
    @test df[end, 21] === 0
    @test df[end, 22] === 0
    @test df[end, 23] === 0
    @test df[end, 24] === 0
    @test df[end, 25] === 0

    #test_group("readtable handles common separators and infers them from extensions.")

    df1 = readtable("$data/separators/sample_data.csv")
    df2 = readtable("$data/separators/sample_data.tsv")
    df3 = readtable("$data/separators/sample_data.wsv")
    df4 = readtable("$data/separators/sample_data_white.txt", separator = ' ')

    @test df1 == df2
    @test df2 == df3
    @test df3 == df4

    readtable("$data/quoting/quotedwhitespace.txt", separator = ' ')

    #test_group("readtable handles common newlines.")

    df = readtable("$data/newlines/os9.csv")
    @test isequal(readtable("$data/newlines/osx.csv"), df)
    @test isequal(readtable("$data/newlines/windows.csv"), df)

    @test isequal(df, readtable("$data/newlines/os9.csv", skipblanks = false))
    @test isequal(df, readtable("$data/newlines/osx.csv", skipblanks = false))
    @test isequal(df, readtable("$data/newlines/windows.csv", skipblanks = false))

    #test_group("readtable treats rows as specified.")

    df1 = readtable("$data/comments/before_after_data.csv", allowcomments = true)
    df2 = readtable("$data/comments/middata.csv", allowcomments = true)
    df3 = readtable("$data/skiplines/skipfront.csv", skipstart = 3)
    df4 = readtable("$data/skiplines/skipfront.csv", skipstart = 4, header = false)
    names!(df4, names(df1))
    df5 = readtable("$data/comments/before_after_data_windows.csv", allowcomments = true)
    df6 = readtable("$data/comments/middata_windows.csv", allowcomments = true)
    df7 = readtable("$data/skiplines/skipfront_windows.csv", skipstart = 3)
    df8 = readtable("$data/skiplines/skipfront_windows.csv", skipstart = 4, header = false)
    names!(df8, names(df1))
    # df9 = readtable("$data/skiplines/skipfront.csv", skipstart = 3, skiprows = 5:6)
    # df10 = readtable("$data/skiplines/skipfront.csv", skipstart = 3, header = false, skiprows = [4, 6])
    # names!(df10, names(df1))

    @test df2 == df1
    @test df3 == df1
    @test df4 == df1

    # Windows EOLS
    @test df5 == df1
    @test df6 == df1
    @test df7 == df1
    @test df8 == df1

    # @test df9 == df1[3:end]
    # @test df10 == df1[[1, 3:end]]

    function normalize_eol!(df)
        for (name, col) in eachcol(df)
            if eltype(col) <: AbstractString
                df[name] = map(s -> replace(s, "\r\n", "\n"), col)
            end
        end
        df
    end

    osxpath = "$data/skiplines/complex_osx.csv"
    winpath = "$data/skiplines/complex_windows.csv"

    opts1 = @compat Dict{Any,Any}(:allowcomments => true)
    opts2 = @compat Dict{Any,Any}(:skipstart => 4, :skiprows => [6, 7, 12, 14, 17], :skipblanks => false)

    df1 = readtable(osxpath; opts1...)
    # df2 = readtable(osxpath; opts2...)
    df1w = readtable(winpath; opts1...)
    # df2w = readtable(winpath; opts2...)

    # @test df2 == df1
    @test normalize_eol!(df1w) == df1
    # @test normalize_eol!(df2w) == df1

    opts1[:nrows] = 3
    opts2[:nrows] = 3

    @test readtable(osxpath; opts1...) == df1[1:3, :]
    # @test readtable(osxpath; opts2...) == df1[1:3, :]
    @test normalize_eol!(readtable(winpath; opts1...)) == df1[1:3, :]
    # @test readtable(winpath; opts2...) == df1[1:3, :]

    # opts2[:header] = false
    # opts2[:skipstart] = 5

    # df2b = readtable(path; opts2...)
    # names!(df2b, names(df1))

    # @test df2b == df1[1:3]

    #test_group("readtable handles custom delimiters.")

    readtable("$data/skiplines/skipfront.csv", allowcomments = true, commentmark = '%')

    readtable("$data/separators/sample_data.csv", quotemark = Char[])
    @test_throws ErrorException readtable("$data/newlines/embedded_osx.csv", quotemark = Char[])
    df = readtable("$data/quoting/single.csv", quotemark = ['\''])
    @test df == readtable("$data/quoting/mixed.csv", quotemark = ['\'', '"'])

    # df = readtable("$data/decimal/period.csv")
    # @test df[2, :A] == 0.3
    # @test df[2, :B] == 4.0

    # @test df == readtable("$data/decimal/comma.tsv", decimal = ',')

    #test_group("readtable column names.")

    ns = [:Var1, :Var2, :Var3, :Var4, :Var5]
    df = readtable("$data/typeinference/mixedtypes.csv")
    names!(df, ns)
    @test df == readtable("$data/typeinference/mixedtypes.csv", names = ns)

    df = readtable("$data/separators/sample_data.csv", header = false, names = ns[1:3])
    @test df[1, :Var1] == 0
    df = readtable("$data/separators/sample_data.csv", names = ns[1:3])
    @test df[1, :Var1] == 1

    #test_group("Properties of data frames returned by readtable method.")

    # Readtable ignorepadding
    io = IOBuffer("A , \tB  , C\n1 , \t2, 3\n")
    @test readtable(io, ignorepadding = true) == DataFrame(A = 1, B = 2, C = 3)

    # Readtable c-style escape options

    df = readtable("$data/escapes/escapes.csv", allowescapes = true)
    @test df[1, :V] == "\t\r\n"
    @test df[2, :V] == "\\\\t"
    @test df[3, :V] == "\\\""

    df = readtable("$data/escapes/escapes.csv")
    @test df[1, :V] == "\\t\\r\\n"
    @test df[2, :V] == "\\\\t"
    @test df[3, :V] == "\\\""

    # df = readtable("$data/escapes/escapes.csv", escapechars = ['"'], nrows = 2)
    # @test df[1, :V] == "\\t\\r\\n"
    # @test df[2, :V] == "\\\\\\\\t"

    # Readtable with makefactors active should only make factors from columns
    # of strings.
    filename = "$data/factors/mixedvartypes.csv"
    df = readtable(filename, makefactors = true)

    @test typeof(df[:factorvar]) == PooledDataArray{Compat.UTF8String,UInt32,1}
    @test typeof(df[:floatvar]) == DataArray{Float64,1}

    # Readtable shouldn't silently drop data when reading highly compressed gz.
    df = readtable("$data/compressed/1000x2.csv.gz")
    @test size(df) == (1000, 2)

    # Readtable type inference
    filename = "$data/typeinference/bool.csv"
    df = readtable(filename)
    @test typeof(df[:Name]) == DataArray{Compat.UTF8String,1}
    @test typeof(df[:IsMale]) == DataArray{Bool,1}
    @test df[:IsMale][1] == true
    @test df[:IsMale][4] == false

    filename = "$data/typeinference/standardtypes.csv"
    df = readtable(filename)
    @test typeof(df[:IntColumn]) == DataArray{Int,1}
    @test typeof(df[:IntlikeColumn]) == DataArray{Float64,1}
    @test typeof(df[:FloatColumn]) == DataArray{Float64,1}
    @test typeof(df[:BoolColumn]) == DataArray{Bool,1}
    @test typeof(df[:StringColumn]) == DataArray{Compat.UTF8String,1}

    filename = "$data/typeinference/mixedtypes.csv"
    df = readtable(filename)
    @test typeof(df[:c1]) == DataArray{Compat.UTF8String,1}
    @test df[:c1][1] == "1"
    @test df[:c1][2] == "2.0"
    @test df[:c1][3] == "true"
    @test typeof(df[:c2]) == DataArray{Float64,1}
    @test df[:c2][1] == 1.0
    @test df[:c2][2] == 3.0
    @test df[:c2][3] == 4.5
    @test typeof(df[:c3]) == DataArray{Compat.UTF8String,1}
    @test df[:c3][1] == "0"
    @test df[:c3][2] == "1"
    @test df[:c3][3] == "f"
    @test typeof(df[:c4]) == DataArray{Bool,1}
    @test df[:c4][1] == true
    @test df[:c4][2] == false
    @test df[:c4][3] == true
    @test typeof(df[:c5]) == DataArray{Compat.UTF8String,1}
    @test df[:c5][1] == "False"
    @test df[:c5][2] == "true"
    @test df[:c5][3] == "true"

    # Readtable defining column types
    filename = "$data/definedtypes/mixedvartypes.csv"

    df = readtable(filename)
    @test typeof(df[:n]) == DataArray{Int,1}
    @test df[:n][1] == 1
    @test typeof(df[:s]) == DataArray{Compat.UTF8String,1}
    @test df[:s][1] == "text"
    @test typeof(df[:f]) == DataArray{Float64,1}
    @test df[:f][1] == 2.3
    @test typeof(df[:b]) == DataArray{Bool,1}
    @test df[:b][1] == true

    df = readtable(filename, eltypes = [Int64, Compat.UTF8String, Float64, Bool])
    @test typeof(df[:n]) == DataArray{Int64,1}
    @test df[:n][1] == 1
    @test typeof(df[:s]) == DataArray{Compat.UTF8String,1}
    @test df[:s][1] == "text"
    @test df[:s][4] == "text ole"
    @test typeof(df[:f]) == DataArray{Float64,1}
    @test df[:f][1] == 2.3
    @test typeof(df[:b]) == DataArray{Bool,1}
    @test df[:b][1] == true
    @test df[:b][2] == false

    df = readtable(filename, eltypes = [Int64, Compat.UTF8String, Float64, Compat.UTF8String])
    @test typeof(df[:n]) == DataArray{Int64,1}
    @test df[:n][1] == 1.0
    @test isna(df[:s][3])
    @test typeof(df[:f]) == DataArray{Float64,1}
    # Float are not converted to int
    @test df[:f][1] == 2.3
    @test df[:f][2] == 0.2
    @test df[:f][3] == 5.7
    @test typeof(df[:b]) == DataArray{Compat.UTF8String,1}
    @test df[:b][1] == "T"
    @test df[:b][2] == "FALSE"

    # Readtable name normalization
    abnormal = "\u212b"
    ns = [:Å, :_B_C_, :_end]
    @test !in(Symbol(abnormal), ns)

    io = IOBuffer(abnormal*",%_B*\tC*,end\n1,2,3\n")
    @test names(readtable(io)) == ns

    # With normalization disabled
    io = IOBuffer(abnormal*",%_B*\tC*,end\n1,2,3\n")
    @test names(readtable(io, normalizenames=false)) == [Symbol(abnormal),Symbol("%_B*\tC*"),:end]

    # Test writetable with NA and compare to the results
    tf = tempname()
    isfile(tf) && rm(tf)
    df = DataFrame(A = @data([1,NA]), B = @data(["b", NA]))
    writetable(tf, df)
    @test readcsv(tf) == ["A" "B"; 1 "b"; "NA" "NA"]

    # Test writetable with nastring set and compare to the results
    isfile(tf) && rm(tf)
    writetable(tf, df, nastring="none")
    @test readcsv(tf) == ["A" "B"; 1 "b"; "none" "none"]

    # Test writetable with append
    df1 = DataFrame(a = @data([1, 2, 3]), b = @data([4, 5, 6]))
    df2 = DataFrame(a = @data([1, 2, 3]), b = @data([4, 5, 6]))
    df3 = DataFrame(a = @data([1, 2, 3]), c = @data([4, 5, 6])) # 2nd column mismatch
    df3b = DataFrame(a = @data([1, 2, 3]), b = @data([4, 5, 6]), c = @data([4, 5, 6])) # number of columns mismatch

    # Would use joinpath(tempdir(), randstring()) to get around tempname
    # creating a file on Windows, but Julia 0.3 has no srand() to unset the
    # seed set in test/data.jl -- annoying for local testing.
    tf = tempname()
    isfile(tf) && rm(tf)

    # Written as normal if file doesn't exist
    writetable(tf, df1, append = true)
    @test readtable(tf) == df1

    # Written as normal if file is empty
    open(io -> print(io, ""), tf, "w")
    writetable(tf, df1, append = true)
    @test readtable(tf) == df1

    # Appends to existing file if append == true
    writetable(tf, df1)
    writetable(tf, df2, header = false, append = true)
    @test readtable(tf) == vcat(df1, df2)

    # Overwrites file if append == false
    writetable(tf, df1)
    writetable(tf, df2)
    @test readtable(tf) == df2

    # Enforces matching column names iff append == true && header == true
    writetable(tf, df2)
    @test_throws KeyError writetable(tf, df3, append = true)
    writetable(tf, df3, header = false, append = true)

    # Enforces matching column count if append == true
    writetable(tf, df3)
    @test_throws DimensionMismatch writetable(tf, df3b, header = false, append = true)

    # Quotemarks are escaped
    tf = tempname()
    isfile(tf) && rm(tf)

    df = DataFrame(a = ["who's"]) # We have a ' in our string

    # Make sure the ' doesn't get escaped for no reason
    writetable(tf, df)
    @test readtable(tf) == df

    # Make sure the ' does get escaped when needed
    writetable(tf, df, quotemark='\'')
    @test readstring(tf) == "'a'\n'who\\'s'\n"

    ### Tests for nonstandard string literals
    # Test basic @csv_str usage
    df1 = csv"""
        name,  age, squidPerWeek
        Alice,  36,         3.14
        Bob,    24,         0
        Carol,  58,         2.71
        Eve,    49,         7.77
        """
    @test size(df1) == (4, 3)
    @test names(df1) == [:name, :age, :squidPerWeek]
    @test df1[1] == ["Alice","Bob","Carol","Eve"]
    @test df1[2] == [36,24,58,49]
    @test df1[3] == [3.14,0,2.71,7.77]
    @test typeof(df1[1]) <: DataArray

    # Test @wsv_str
    df2 = wsv"""
        name  age squidPerWeek
        Alice  36         3.14
        Bob    24         0
        Carol  58         2.71
        Eve    49         7.77
        """
    @test df2 == df1

    # Test @tsv_str
    df3 = tsv"""
        name	age	squidPerWeek
        Alice	36	3.14
        Bob	24	0
        Carol	58	2.71
        Eve	49	7.77
        """
    @test df3 == df1

    # csv2 can't be tested until non-'.' decimals are implemented
    #df4 = csv2"""
    #    name;  age; squidPerWeek
    #    Alice;  36;         3,14
    #    Bob;    24;         0
    #    Carol;  58;         2,71
    #    Eve;    49;         7,77
    #    """
    #@test df4 == df1

    # Test 'f' flag
    df5 = csv"""
        name,  age, squidPerWeek
        Alice,  36,         3.14
        Bob,    24,         0
        Carol,  58,         2.71
        Eve,    49,         7.77
        """f
    @test typeof(df5[1]) <: PooledDataArray

    # Test 'c' flag
    df6 = csv"""
        name,  age, squidPerWeek
        Alice,  36,         3.14
        Bob,    24,         0
        #Carol,  58,         2.71
        Eve,    49,         7.77
        """c
    @test df6 == df1[[1,2,4],:]

    # Test 'H' flag
    df7 = csv"""
        Alice,  36,         3.14
        Bob,    24,         0
        Carol,  58,         2.71
        Eve,    49,         7.77
        """H
    @test names(df7) == [:x1,:x2,:x3]
    @test Array(df7) == Array(df1)

    # Test multiple flags at once
    df8 = csv"""
        Alice,  36,         3.14
        Bob,    24,         0
        #Carol,  58,         2.71
        Eve,    49,         7.77
        """fcH
    @test typeof(df8[1]) <: PooledDataArray
    @test names(df8) == [:x1,:x2,:x3]
    @test Array(df8) == Array(df1[[1,2,4],:])

    # Test invalid flag
    # Need to wrap macro call inside eval to prevent the error from being
    # thrown prematurely
    @test_throws ArgumentError eval(:(csv"foo,bar"a))
end


text:
import Pkg

using DifferentialEquations
using LinearAlgebra
using Random
using Plots

using Statistics
using Interpolations

using Revise
include("SimulationBasedStability.jl")
include("structures.jl")
#include("functinos.jl")
xs = 1:0.2:5
f(x) = sin(x)
A = [[f(x),x] for x in xs]
interp_linear = LinearInterpolation(xs, A,Gridded(Linear()))
interp_linear(3.1) # approximately log(3.1)
StateSmaplingTime=-0.2:0.1:1

interp_linear = LinearInterpolation(StateSmaplingTime,  [u0 for x in StateSmaplingTime])
LinearInterpolation(StateSmaplingTime,[CompRand(u0) for t in StateSmaplingTime])
collect(StateSmaplingTime)
function delay_mathieu_model(du,u,h,p,t)
    δ, ϵ, b, κ, τ = p

    du[1] = - κ * u[1] - (δ + ϵ * cos(t)) * u[2] + b * h(p, t-τ)[2]
    du[2] = u[1]
end

#h(p, t; idxs) = (rand()-0.5) .+ 1im*(rand()-0.5)*10.0
h(p, t) = [0.0,0.0]# .* (1.0+1.0im)
δ=3.0
ϵ=2.0
b=-0.15
κ=0.1
τ=2pi
T=2pi
p=( δ, ϵ, b, κ, τ )
lags = [τ]
taumax=maximum(lags)
tspan = (0.0, T) # The end of the integration time considert to be the timeperiod of the system.

#u0=[1.0+1.0im,1.0+1.0im]
u0=[0.0,0.0];
prob = DDEProblem(delay_mathieu_model, u0, h, tspan, p; constant_lags = lags);
alg = MethodOfSteps(Tsit5());

#abc=solve(prob,alg,reltol=1e-3,abstol=1e-3)

DelayMathieu_egi_prblem=dynamic_problem(prob,alg,taumax,Historyresolution=40,eigN=2);

compute_eig!(DelayMathieu_egi_prblem);

for k=1:10
iterate!(DelayMathieu_egi_prblem);
ei=compute_eig!(DelayMathieu_egi_prblem)
@show maximum(abs.(ei.values))
end

#alg=Rosenbrock23(autodiff=false)
#alg=Rodas4(autodiff=false)

#big_u0=big.(u0)

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

DelayMathieu_egi_prblem=dynamic_problem(prob,alg,taumax,Historyresolution=500,eigN=4);

Nshift=0;
Si,Vi =SVi1real(DelayMathieu_egi_prblem,1)
plot(DelayMathieu_egi_prblem.StateSmaplingTime.+ T*(Nshift-1),real.(Si),linewidth=2)
plot!(DelayMathieu_egi_prblem.StateSmaplingTime.+ T*(Nshift),real.(Vi),linewidth=2)

ei=compute_eig!(DelayMathieu_egi_prblem)
@show (abs.(ei))
#DelayMathieu_egi_prblem.StateCombinations[:] .= [1.0,-10.0,0.0,1.0]

A =DelayMathieu_egi_prblem.StateCombinations
plot!(DelayMathieu_egi_prblem.StateSmaplingTime.+ T*(Nshift),real.(Vi*A),linewidth=5)

Nshift += 1;
iterate!(DelayMathieu_egi_prblem);
Si,Vi =SVi1real(DelayMathieu_egi_prblem,1) #Si2,Vi2
plot!(DelayMathieu_egi_prblem.StateSmaplingTime.+ T*(Nshift-1),real.(Si),linewidth=2)
plot!(DelayMathieu_egi_prblem.StateSmaplingTime.+ T*(Nshift),real.(Vi),linewidth=2)

plot!(legend = :bottomleft)
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#Vi*DelayMathieu_egi_prblem.StateCombinations .- Si2

SV=Si'*Vi;
SS=Si'*Si;

Hi=SS\SV;

using GenericSchur
S = schur(Hi .+ 0im)
VL = eigvecs(S, left=true)
VR = eigvecs(S, left=false)

A' * diagm(μs) * A
-Hi

Snorm=[norm(Si[:,k]) for k in 1:DelayMathieu_egi_prblem.eigN ];

solset=DelayMathieu_egi_prblem.SolutionSet;
plot(solset[1])
plot!.(solset[2:end-1])
plot!(solset[end])

#<<<<<<<<<<<<<<<<<<<<<<<<<<<
plot()
for kk=1:6
s=DelayMathieu_egi_prblem.SolutionSet[kk]
ti=collect(-taumax:0.01:T)
SS=getvalues(s,ti)
Nshift = 0;
plot!(ti .+ Nshift*T,real.(getvalues(s,ti)))
end

iterate!(DelayMathieu_egi_prblem);
ei=compute_eig!(DelayMathieu_egi_prblem)
#@show ei.values
@show maximum(abs.(ei.values))
Nshift += 1;
for kk=1:6
s=DelayMathieu_egi_prblem.SolutionSet[kk]
SS=getvalues(s,ti)
plot!(ti .+ Nshift*T,real.(getvalues(s,ti)))
end
plot!()

#<<<<<<<<<<<<<<<<<<<<<<<<<<<

SV=Si'*Vi;
SS=Si'*Si;

Snorm=[norm(Si[:,k]) for k in 1:dp.eigN ];

solset=DelayMathieu_egi_prblem.SolutionSet;
plot(solset[1])
plot!.(solset[2:end-1])
plot!(solset[end])

#using StaticArrays
#A = @SMatrix rand(3,2)
# f(u,p,t)=A*u

Ai=diagm(0 => ones(Float64,eigvecnum))

u0=1.0+.5im

rand(Complex{typeof(u0[1])}, 10)

Ntime=10
taumax=maximum(prob.constant_lags)
timelocations=-taumax  :taumax/(Ntime-1) :0

#Vi=[solset[solind](t+T) for t in timelocations, solind in 1:eigvecnum];#initialization of the Starting Vector
#------------------------------

for kiter in 1:25
#Si=Vi*Ai;
solsetnew=deepcopy(solset);
for k in 1:eigvecnum
    #h(p, t; idxs) = histopryremapp(t+T,Ai,solset,k)[idxs]#+(rand(size(idxs)[1]) .+ 1im*rand(size(idxs)[1]) .- 0.5 .- 0.5im)*10.0 # overwrites the original function even inside the prob.variable
    h(p, t;) = histopryremapp(t+T,Ai,solset,k)
    prob.u0 .= h(p, 0)
   solsetnew[k] = solve(prob,alg,reltol=1e-3,abstol=1e-3);
end
solset=solsetnew;
Si=[solset[solind](t) for t in timelocations, solind in 1:eigvecnum];

Vi=[solset[solind](t+T) for t in timelocations, solind in 1:eigvecnum];#initialization of the Starting Vector

SV=Si'*Vi;
SS=Si'*Si;

#VV=Vi'*Vi;
#SV=Si'*Vi;

Snorm=[norm(Si[:,k]) for k in 1:eigvecnum ];

Hi=SS\SV
#Hi=SV\VV
#Si\Vi
Eigs=eigen(Hi);
@show μs=Eigs.values
Ai=Eigs.vectors  ./ Snorm ./ Eigs.values;
end
plot(solset[1])
plot!.(solset[2:end-1])
plot!(solset[end])

#plot(1:eigvecnum,1:eigvecnum,real.(Ai))

#------------------------------

Ntime=100
taumax=maximum(prob.constant_lags)
timelocations=-taumax  :taumax/(Ntime-1) :0

Vi=[solset[solind](t) for t in timelocations, solind in 1:eigvecnum];

plot(solset[1].t)

    # figure(9)
    # clf()
    # surf(abs.(Hi))

    Si=Vi*Eigs.vectors
    for iloc in 1:eigvecnum
        sinorminv=1/norm(Si[:,iloc])
        Si[:,iloc] *= sinorminv
    end
    push!(eigsol_inter,Eigs.values)
println("---------------------------------")

solset = [solve(prob,alg) for k in 1:eigvecnum]

plot(solset[1])
plot!.(solset[2:end-1])
plot!(solset[end])

size(solset)

@show Si=[sol(t+T) for t in timelocations];

solset[1](2)

@show Si=[sol(t+T) for t in timelocations];

sum([Si'[i]*Si[i] for i in 1:Ntime])
Si'*Si
plot(timelocations,Si)
Si'[end]*Si[end]
# initialization of the Eigen vector
N = 50
eigvecnum = 6
# δ, ϵ, b, κ, τ = [2.0,1.0,0.02,0.0,2*pi] # [1.0,2.0,1,0.1]
# p=[δ, ϵ, b, κ, τ]
p=[δ,1.0,b,0.0,2*pi]
#
Si = rand(N, eigvecnum)
eigsol_inter=Array{Array{ComplexF64}}(undef,0)
#@time begin
for k in 1:6
    #global(Si)
    Vi=mappinggeneration(Si,p)
    Hi=Si\Vi
    # figure(9)
    # clf()
    # surf(abs.(Hi))
    Eigs=eigen(Hi)
    Si=Vi*Eigs.vectors
    for iloc in 1:eigvecnum
        sinorminv=1/norm(Si[:,iloc])
        Si[:,iloc] *= sinorminv
    end
    push!(eigsol_inter,Eigs.values)

    # figure(14)
    # plot(k*ones(size(Eigs.values)),abs.(Eigs.values), linestyle="none", marker="o",color="green", markerfacecolor="blue", markersize=4)

    #figure(4)
    # clf()
    #plot(real(Eigs.values),imag(Eigs.values), linestyle="none", marker="o",color="green", markerfacecolor="blue", markersize=4)
    # fi=0:0.01:2*pi
    # plot(sin.(fi),cos.(fi))
    #pause(0.0001)
end

function mappinggeneration(Si,p)
    # figure(3)
    # clf()
    Vi = zero(Si)
    τ=p[5]
    N=size(Si,1)
    for i in 1:size(Si,2)

        hist = historymaker(Si[:,i],τ,N)

        # p = [1.0,2.0,1,0.1]

        T = 2 * pi
        lags = [τ]
        tspan = (0.0, T)
        Δt = 1e-6
        u0 = [(hist(p, 0, idxs = 2) - hist(p, -Δt, idxs = 2)) / Δt, hist(p, 0, idxs = 2)]
        prob = DDEProblem(delay_mathieu_model, u0, hist, tspan, p; constant_lags = lags)

        alg = MethodOfSteps(Tsit5())
        sol = solve(prob, alg,reltol=1e-3,abstol=1e-3)#,reltol=1e-8,abstol=1e-8)

        tplot = LinRange(T - τ, T, N)
        solplot = [sol(ttt)[2] for ttt in tplot]
        Vi[:,i] .= solplot

        # figure(3)
        # # clf()
        # tplot = LinRange(-τ, 0, N)
        # solplot = [hist(p, ttt, idxs = 1) for ttt in tplot]
        # plot(tplot, real.(solplot), marker="o",color="green", markerfacecolor="blue", markersize=4)
        #
        # # tplot = LinRange(0, T - τ, N)
        # # solplot = [sol(ttt)[1] for ttt in tplot]
        # # plot(tplot, real.(solplot), marker="o",color="green", markerfacecolor="blue", markersize=4)
        #
        # tplot = LinRange(T - τ, T, N)
        # solplot = [sol(ttt)[2] for ttt in tplot]
        # plot(tplot, real.(solplot), marker="o",color="green", markerfacecolor="blue", markersize=4)

    end

    return Vi

end

function delayoscillmapp(δ, b)

# initialization of the Eigen vector
N = 50
eigvecnum = 6
# δ, ϵ, b, κ, τ = [2.0,1.0,0.02,0.0,2*pi] # [1.0,2.0,1,0.1]
# p=[δ, ϵ, b, κ, τ]
p=[δ,1.0,b,0.0,2*pi]
#
Si = rand(N, eigvecnum)
eigsol_inter=Array{Array{ComplexF64}}(undef,0)
#@time begin
for k in 1:6
    #global(Si)
    Vi=mappinggeneration(Si,p)
    Hi=Si\Vi
    # figure(9)
    # clf()
    # surf(abs.(Hi))
    Eigs=eigen(Hi)
    Si=Vi*Eigs.vectors
    for iloc in 1:eigvecnum
        sinorminv=1/norm(Si[:,iloc])
        Si[:,iloc] *= sinorminv
    end
    push!(eigsol_inter,Eigs.values)

    # figure(14)
    # plot(k*ones(size(Eigs.values)),abs.(Eigs.values), linestyle="none", marker="o",color="green", markerfacecolor="blue", markersize=4)

    #figure(4)
    # clf()
    #plot(real(Eigs.values),imag(Eigs.values), linestyle="none", marker="o",color="green", markerfacecolor="blue", markersize=4)
    # fi=0:0.01:2*pi
    # plot(sin.(fi),cos.(fi))
    #pause(0.0001)
end

#end

# figure(16)
# clf()
# for keig in 1:min(10,eigvecnum)#1:2
#     scatter(1:length(eigsol_inter),log.(abs.([abs(e[keig])-abs(eigsol_inter[end][keig]) for e in eigsol_inter])) ./ log(10))
#     #scatter(1:length(eigsol_inter),log.(abs.([e[keig] for e in eigsol_inter])))
# end
#
#
# figure(4)
# clf()
# plot(real(eigsol_inter[end]),imag(eigsol_inter[end]), linestyle="none", marker="o",color="green", markerfacecolor="blue", markersize=4)
# fi=0:0.01:2*pi
# plot(sin.(fi),cos.(fi))
# pause(0.1)

return abs(eigsol_inter[end][1])-1

end

using MDBM

using PyPlot;
pygui(true);

mymdbm=MDBM_Problem(delayoscillmapp,[-1:0.25:9,-1.5:0.25:1.5])

for kinter in 1:4
MDBM.solve!(mymdbm,1)

x_eval,y_eval=getevaluatedpoints(mymdbm)

x_sol,y_sol=getinterpolatedsolution(mymdbm)

fig = figure(112);clf()
scatter(x_eval,y_eval,s=5)
scatter(x_sol,y_sol,s=5)
pause(0.01)
end


text: const UInt64Digit = Ref(true)


text:
export models

const allmodels = [
  PekerisRayModel,
  RaySolver,
  Bellhop
]

"""
    models()
    models(env::UnderwaterEnvironment)

List available models, optionally for a given environment.
"""
function models(env=missing)
  mlist = []
  for m ∈ allmodels
    try
      check(m, env)
      push!(mlist, m)
    catch ex
      # don't add to list
    end
  end
  mlist
end

"""
$(SIGNATURES)
Register model in the list of available models.
"""
addmodel!(mtype) = push!(allmodels, mtype)