text:
using CUDAdrv

using Compat

## pointer

# construction
DevicePtr{Void}()
@test_throws InexactError DevicePtr{Void}(C_NULL)

# conversion
Base.unsafe_convert(DevicePtr{Void}, C_NULL)
Base.unsafe_convert(Ptr{Void}, DevicePtr{Void}())
@test_throws InexactError convert(Ptr{Void}, DevicePtr{Void}())
@test_throws InexactError convert(DevicePtr{Void}, C_NULL)

let
    nullptr = DevicePtr{Void}()

    @test eltype(DevicePtr{Void}) == Void
    @test eltype(nullptr) == Void
    @test isnull(nullptr)

    @test_throws InexactError convert(Ptr{Void}, nullptr)
    @test_throws InexactError convert(DevicePtr{Void}, C_NULL)
end

## errors

let
    ex = CuError(0)
    @test CUDAdrv.name(ex) == :SUCCESS
    @test CUDAdrv.description(ex) == "Success"

    io = IOBuffer()
    showerror(io, ex)
    str = takebuf_string(io)

    @test contains(str, "0")
    @test contains(str, "Success")
end

let
    ex = CuError(0, "foobar")

    io = IOBuffer()
    showerror(io, ex)
    str = takebuf_string(io)

    @test contains(str, "foobar")
end

## base

CUDAdrv.@apicall(:cuDriverGetVersion, (Ptr{Cint},), Ref{Cint}())

@test_throws ErrorException CUDAdrv.@apicall(:cuNonexisting, ())

@test_throws ErrorException CUDAdrv.@apicall(:cuDummyAvailable, ())
@test_throws CUDAdrv.CuVersionError CUDAdrv.@apicall(:cuDummyUnavailable, ())

CUDAdrv.trace(prefix=" ")

@test_throws ErrorException eval(
    quote
        foo = :bar
        CUDAdrv.@apicall(foo, ())
    end
)

typealias CuDevice_t Cint
try
    CUDAdrv.@apicall(:cuDeviceGet, (Ptr{CuDevice_t}, Cint), Ref{CuDevice_t}(), devcount())
catch e
    e == CUDAdrv.ERROR_INVALID_DEVICE || rethrow(e)
end

CUDAdrv.vendor()

dev = CuDevice(0)
ctx = CuContext(dev)

## version

CUDAdrv.version()

## devices

name(dev)
totalmem(dev)
attribute(dev, CUDAdrv.MAX_THREADS_PER_BLOCK)
capability(dev)
list_devices()

## context

@test device(ctx) == dev
synchronize(ctx)

let
    ctx2 = CuContext(dev)       # implicitly pushes
    @test current_context() == ctx2
    @test_throws ArgumentError device(ctx)

    push(ctx)
    @test current_context() == ctx

    pop()
    @test current_context() == ctx2

    pop()
    @test current_context() == ctx

    destroy(ctx2)
end

## module

let
    md = CuModuleFile(joinpath(dirname(@__FILE__), "ptx/vadd.ptx"))

    vadd = CuFunction(md, "vadd")

    unload(md)
end

CuModuleFile(joinpath(dirname(@__FILE__), "ptx/vadd.ptx")) do md
    vadd = CuFunction(md, "vadd")
end

let
    f = open(joinpath(dirname(@__FILE__), "ptx/vadd.ptx"))
    ptx = readstring(f)
    close(f)

    md = CuModule(ptx)
    vadd = CuFunction(md, "vadd")
    unload(md)
end

try
    CuModule("foobar")
catch ex
    ex == CUDAdrv.ERROR_INVALID_IMAGE  || rethrow(ex)
end

let
    md = CuModuleFile(joinpath(dirname(@__FILE__), "ptx/global.ptx"))

    var = CuGlobal{Int32}(md, "foobar")
    @test eltype(var) == Int32
    @test eltype(typeof(var)) == Int32

    @test_throws ArgumentError CuGlobal{Int64}(md, "foobar")

    set(var, Int32(42))
    @test get(var) == Int32(42)

    unload(md)
end

let
    link = CuLink()

    # regular string
    open(joinpath(dirname(@__FILE__), "ptx/empty.ptx")) do f
        addData(link, "vadd_parent", readstring(f), CUDAdrv.PTX)
    end

    # string as vector of bytes
    open(joinpath(dirname(@__FILE__), "ptx/empty.ptx")) do f
        addData(link, "vadd_parent", convert(Vector{UInt8}, readstring(f)), CUDAdrv.PTX)
    end

    # string containing \0
    @test_throws ArgumentError addData(link, "vadd_parent", "\0", CUDAdrv.PTX)
    @test_throws CuError addData(link, "vadd_parent", "\0", CUDAdrv.OBJECT)

    # vector of bytes containing \0
    @test_throws ArgumentError addData(link, "vadd_parent", convert(Vector{UInt8}, "\0"), CUDAdrv.PTX)
    @test_throws CuError addData(link, "vadd_parent", convert(Vector{UInt8}, "\0"), CUDAdrv.OBJECT)

    destroy(link)
end

let
    link = CuLink()
    addFile(link, joinpath(dirname(@__FILE__), "ptx/vadd_child.ptx"), CUDAdrv.PTX)
    open(joinpath(dirname(@__FILE__), "ptx/vadd_parent.ptx")) do f
        addData(link, "vadd_parent", readstring(f), CUDAdrv.PTX)
    end
    obj = complete(link)

    md = CuModule(obj)
    destroy(link)

    vadd = CuFunction(md, "vadd")
    unload(md)
end

## memory

let
    ptr = cualloc(Int, 1)
    free(ptr)

    @test_throws ArgumentError cualloc(Function, 1)   # abstract
    @test_throws ArgumentError cualloc(Array{Int}, 1) # non-leaftype
    @test_throws ArgumentError cualloc(Int, 0)
end

let
    type MutablePtrFree
        foo::Int
        bar::Int
    end
    ptr = cualloc(MutablePtrFree, 1)
    copy!(ptr, MutablePtrFree(0,0))
    free(ptr)
end

let
    type MutableNonPtrFree
        foo::Int
        bar::String
    end
    ptr = cualloc(MutableNonPtrFree, 1)
    @test_throws ArgumentError copy!(ptr, MutableNonPtrFree(0,""))
    free(ptr)
end

let
    dev_array = CuArray(Int32, 10)
    cumemset(dev_array.ptr, UInt32(0), 10)
    host_array = Array(dev_array)

    @test all(x -> x==0, host_array)

    free(dev_array)
end

## stream

let
    s = CuStream()
    synchronize(s)
    destroy(s)

    synchronize(default_stream())
end

## execution

let
    @test CUDAdrv.CuDim3((3,2,1)) == CUDAdrv.CuDim3(3,2,1)
    @test CUDAdrv.CuDim3((3,2))   == CUDAdrv.CuDim3(3,2,1)
    @test CUDAdrv.CuDim3(3)       == CUDAdrv.CuDim3(3,1,1)
end

let
    md = CuModuleFile(joinpath(dirname(@__FILE__), "ptx/vectorops.ptx"))
    vadd = CuFunction(md, "vadd")
    vmul = CuFunction(md, "vmul")
    vsub = CuFunction(md, "vsub")
    vdiv = CuFunction(md, "vdiv")

    a = rand(Float32, 10)
    b = rand(Float32, 10)
    ad = CuArray(a)
    bd = CuArray(b)

    # Addition
    let
        c = zeros(Float32, 10)
        cd = CuArray(c)
        cudacall(vadd, 10, 1, (DevicePtr{Cfloat},DevicePtr{Cfloat},DevicePtr{Cfloat}), ad, bd, cd)
        c = Array(cd)
        @test c ≈ a+b
        free(cd)
    end

    # Subtraction
    let
        c = zeros(Float32, 10)
        cd = CuArray(c)
        cudacall(vsub, 10, 1, (DevicePtr{Cfloat},DevicePtr{Cfloat},DevicePtr{Cfloat}), ad, bd, cd)
        c = Array(cd)
        @test c ≈ a-b
        free(cd)
    end

    # Multiplication
    let
        c = zeros(Float32, 10)
        cd = CuArray(c)
        cudacall(vmul, 10, 1, (DevicePtr{Cfloat},DevicePtr{Cfloat},DevicePtr{Cfloat}), ad, bd, cd)
        c = Array(cd)
        @test c ≈ a.*b
        free(cd)
    end

    # Division
    let
        c = zeros(Float32, 10)
        cd = CuArray(c)
        cudacall(vdiv, 10, 1, (DevicePtr{Cfloat},DevicePtr{Cfloat},DevicePtr{Cfloat}), ad, bd, cd)
        c = Array(cd)
        @test c ≈ a./b
        free(cd)
    end

    free(ad)
    free(bd)
    unload(md)
end

## events

let
    start = CuEvent()
    stop = CuEvent()
    record(start)
    record(stop)
    synchronize(stop)
    @test elapsed(start, stop) > 0
    destroy(start)
    destroy(stop)
end

## array

let
    # Inner constructors
    CuArray{Int,1}((1,))
    CuArray{Int,1}((1,), DevicePtr{Int}())

    # Outer constructors
    CuArray{Int}(1)
    CuArray{Int}((1,2))
    CuArray{Int}(1, DevicePtr{Int}())
    CuArray{Int}((1,2), DevicePtr{Int}())
    CuArray(Int, 1)
    CuArray(Int, (1,2))

    # Conversions
    @test Base.unsafe_convert(DevicePtr{Int}, CuArray{Int,1}((1,), DevicePtr{Int}())) == DevicePtr{Int}()
    @test pointer(CuArray{Int,1}((1,), DevicePtr{Int}())) == DevicePtr{Int}()

    # Negative test cases
    a = rand(Float32, 10)
    ad = CuArray(Float32, 5)
    @test_throws ArgumentError copy!(ad, a)
    @test_throws ArgumentError copy!(a, ad)

    # Utility
    @test ndims(ad) == 1
    @test size(ad, 1) == 5
    @test size(ad, 2) == 1
    @test eltype(ad) == Float32
    @test eltype(typeof(ad)) == Float32

    free(ad)
end

let
    # ghost type
    @test_throws ArgumentError CuArray([x->x*x for i=1:10])

    # non-isbits elements
    @test_throws ArgumentError CuArray(["foobar" for i=1:10])
    @test_throws ArgumentError CuArray(Function, 10)
    @test_throws ArgumentError CuArray(Function, (10, 10))
end

## profile

@cuprofile begin end

destroy(ctx)


text:
Index = AbstractVector{I} where I <: Integer

struct IndexedKroneckerProduct{T<:Any,TK<:AbstractKroneckerProduct} <: GeneralizedKroneckerProduct{T}
    K::TK
    p
    q
    r
    t
    function IndexedKroneckerProduct(K::AbstractKroneckerProduct, p::Index, q::Index, r::Index,
                            t::Index)
        if order(K) != 2
            throw(DimensionMismatch(
                "Indexed Kronecker only implemented for second order Kronecker systems"
            ))
        elseif !(length(p) == length(q) && length(r) == length(t))
            throw(DimensionMismatch("Indices should have matching lengths"))
        elseif minimum((minimum.((p, q, r, t)))) <= 0
            throw(BoundsError("Negative indices not allowed"))
        end
        A, B = getmatrices(K)
        (m, n) = size(B)
        (k, l) = size(A)
        if !(maximum(p) ≤ m && maximum(q) ≤ k)
            throw(BoundsError("Indices exeed matrix bounds"))
        end
        if !(maximum(r) ≤ n && maximum(t) ≤ l)
            throw(BoundsError("Indices exeed matrix bounds"))
        end
        return new{eltype(K),typeof(K)}(K, p, q, r, t)
    end
end

# Best way to do multi-line strings?

_dm(str) = throw(DimensionMismatch(str))
_ibe(str) = throw(BoundsError(str))

getindices(K::IndexedKroneckerProduct) = K.p, K.q, K.r, K.t

getmatrices(K::IndexedKroneckerProduct) = getmatrices(K.K)

Base.size(K::IndexedKroneckerProduct) = (length(K.p), length(K.t))

function Base.getindex(K::IndexedKroneckerProduct, i::Int, j::Int)
    A, B = getmatrices(K)
    p, q, r, t = getindices(K)
    return B[p[i],r[j]] * A[q[i],t[j]]
end

function Base.getindex(K::AbstractKroneckerProduct, p::Index, q::Index, r::Index, t::Index)
    N, M = getmatrices(K)
    a, b = size(M)
    c, d = size(N)
    return IndexedKroneckerProduct(K, p, q, r, t)
end

Base.eltype(K::IndexedKroneckerProduct) = eltype(K.K)

function Base.collect(K::IndexedKroneckerProduct)
    f, e = size(K)
    Kexplicit = Array{eltype(K)}(undef, f, e)
    N, M = getmatrices(K.K)
    p = K.p
    q = K.q
    r = K.r
    t = K.t
    for h in 1:f, g in 1:e
        Kexplicit[h, g] = N[q[h], t[g]] * M[p[h], r[g]]
    end
    return Kexplicit
end

function genvectrick!(M, N, v, u, p, q, r, t)
    # computes N ⊗ M

    a, b = size(M)
    c, d = size(N)
    e = length(v)
    f = length(u)

    @assert maximum(p) ≤ a && maximum(q) ≤ c
    @assert maximum(r) ≤ b && maximum(t) ≤ d
    u .= 0  # reset for inplace
    if a * e + d * f < c * e + b *f
        # compute T = VM'
        T = zeros(eltype(v), d, a)
        @simd for h in 1:e
            i, j = r[h], t[h]
            @simd for k in 1:a
                @inbounds T[j,k] += v[h] * M[k,i]
            end
        end
        @simd for h in 1:f
            i, j = p[h], q[h]
            @simd for k in 1:d
                @inbounds u[h] += N[j,k] * T[k,i]
            end
        end
    else
        # compute S = NV
        S = zeros(eltype(v), d, a)
        @simd for h in 1:e
            i, j = r[h], t[h]
            @simd for k in 1:c
                @inbounds S[k,j] += v[h] * N[k,j]
            end
        end
        @simd for h in 1:f
            i, j = p[h], q[h]
            @simd for k in 1:b
                @inbounds u[h] += S[j,k] * M[i,k]
            end
        end
    end
    return u
end

function genvectrick(M, N, v, p, q, r, t)
    # computes N ⊗ M
    f = length(p)
    u = zeros(eltype(v), f)
    return genvectrick!(M, N, v, u, p, q, r, t)
end

function Base.:*(K::IndexedKroneckerProduct, v::AbstractVector)
    if length(v) != size(K, 2)
        throw(DimensionMismatch(string(
            "Size indexed Kronecker system ($(size(K))) does not ",
            "match vector size ($(size(v)))",
        )))
    end
    N, M = getmatrices(K.K)
    p, q, r, t = getindices(K)
    return genvectrick(M, N, v, p, q, r, t)
end

ext: jl


text:
printstyled("  read_meta equivalencies\n", color=:light_green)
printstyled("    full (XML, RESP, dataless)\n", color=:light_green)

fname         = "JRO.sacpz"
dataless_name = "CC.dataless"
sxml_file     = path*"/SampleFiles/XML/fdsnws-station_2019-09-11T06_26_58Z.xml"
resp_file     = path*"/SampleFiles/SEED/JRO.resp"
dataless_file = path*"/SampleFiles/SEED/"*dataless_name
dataless_wc   = path*"/SampleFiles/SEED/CC.*"
sacpz_file    = path*"/SampleFiles/SAC/"*fname
sacpz_wc      = path*"/SampleFiles/SAC/JRO.sacp*"

S1 = read_meta("sxml", sxml_file, s="2016-01-01T00:00:00", memmap=true, msr=true)
S2 = read_meta("resp", resp_file, memmap=true, units=true)
S3 = read_meta("dataless", dataless_file, memmap=true, s="2016-01-01T00:00:00", units=true)[56:58]
S4 = read_meta("sacpz", sacpz_file, memmap=true)

@test_throws ErrorException read_meta("dataless", sxml_file)
@test_throws ErrorException read_meta("deez", "nutz.sac")
@test_throws ErrorException read_meta("deez", sxml_file)

C = Array{Char,2}(undef, 10, 3)
fill!(C, ' ')
SA = Array{String, 2}(undef, 10, 3)
for k = 1:3
  fstr = ""
  n = 0
  # println("    ===== ", S1.id[k], " =====")
  R1 = S1.resp[k]
  R2 = S2.resp[k]
  R3 = S3.resp[k]
  for f in fieldnames(MultiStageResp)
    fstr = uppercase(String(f))
    n += 1
    f1 = getfield(R1, f)
    f2 = getfield(R2, f)
    f3 = getfield(R3, f)
    t = min(isequal(f1, f2), isequal(f1, f3))
    if t
      C[n,k] = '='
    else

      # Check for same size, type
      L = length(f1)
      if L != length(f2) || L != length(f3)
        C[n,k] = 'f'
        continue
      elseif typeof(f1) != typeof(f2) || typeof(f1) != typeof(f3)
        C[n,k] = 'f'
        continue
      end

      # Check for approximately equal fields
      T = falses(L)
      for i = 1:L
        i1 = getindex(f1, i)
        i2 = getindex(f2, i)
        i3 = getindex(f3, i)
        t2 = min(isequal(i1, i2), isequal(i1, i3))
        if t2 == true
          T[i] = true
        else
          T1 = typeof(i1)
          if T1 != typeof(i2) || T1 != typeof(i3)
            C[n,k] = 'f'
            break

          # Easyfor a bitstype
          elseif isbitstype(T1)
            if isapprox(i1, i2) && isapprox(i1, i3)
              C[n,k] = '≈'
              T[i] = true
              continue
            end
          else
            FF = fieldnames(T1)

            # Only possible for a String in these Types
            if isempty(FF)
              C[n,k] = 'f'

            # Check for approximately equal subfields
            else
              TT = falses(length(FF))
              for (j,g) in enumerate(FF)
                j1 = getfield(i1, g)
                j2 = getfield(i2, g)
                j3 = getfield(i3, g)

                # Dimension mismatch
                if !(length(j1) == length(j2) == length(j3))
                  C[n,k] = 'f'

                # True mismatch
                elseif min(isapprox(j1, j2), isapprox(j1, j3)) == false
                  C[n,k] = 'f'

                # Approx. equality
                else
                  TT[j] = true
                end
              end

              # If they're all approximate, set T[i] to true
              if minimum(TT) == true
                T[i] = true
              end
            end
          end
        end
      end

      # If they're all approximate, set C[n,k] to ≈
      if minimum(T) == true
        C[n,k] = '≈'
      end
    end
    SA[n,k] = (k == 1 ? lpad(fstr * ": ", 12) : "") * lpad(C[n,k], 5)
    @test C[n,k] in ('≈', '=')
  end
end
println("")
println(" "^12,
        lpad(S1.id[1][end-3:end], 6),
        lpad(S1.id[2][end-3:end], 6),
        lpad(S1.id[3][end-3:end], 6))
println(" "^12, "|", "="^5, "|", "="^5, "|", "="^5)
for i = 1:size(SA,1)
  println(join(SA[i,:], " "))
end
println("")

printstyled("    one-stage (SACPZ, XML)\n", color=:light_green)
S1 = read_meta("sxml", sxml_file, s="2016-01-01T00:00:00", msr=false)
S4 = read_meta("sacpz", sacpz_file)
K = Array{Char,2}(undef, 3, 3)
fill!(K, ' ')
SA = Array{String, 2}(undef, 3, 3)
for k = 1:3
  n = 0
  R1 = S1.resp[k]
  R2 = S4.resp[k]
  for f in fieldnames(PZResp)
    fstr = uppercase(String(f))
    (f == :f0) && continue # SACPZ lacks f0 for some reason
    n += 1
    f1 = getfield(R1, f)
    f2 = getfield(R2, f)
    t = isequal(f1, f2)
    if t == true
      K[n,k] = '='
    else

      # Check for same size, type
      L = length(f1)
      if L != length(f2)
        K[n,k] = 'f'
        continue
      elseif typeof(f1) != typeof(f2)
        K[n,k] = 'f'
        continue
      elseif isapprox(f1, f2)
        K[n,k] = '≈'
      end
    end
    @test K[n,k] in ('≈', '=')
    SA[n,k] = (k == 1 ? lpad(fstr * ": ", 12) : "") * lpad(K[n,k], 5)
  end
end
println("")
println(" "^12,
        lpad(S1.id[1][end-3:end], 6),
        lpad(S1.id[2][end-3:end], 6),
        lpad(S1.id[3][end-3:end], 6))
println(" "^12, "|", "="^5, "|", "="^5, "|", "="^5)
for i = 1:size(SA,1)
  println(join(SA[i,:], " "))
end
println("")

S1 = read_meta("sacpz", sacpz_file)
S2 = read_meta("sacpz", sacpz_wc)
S3 = read_meta("dataless", dataless_wc, s="2016-01-01T00:00:00")
for f in SeisIO.datafields
  if (f in (:src, :notes)) == false
    @test isequal(getfield(S1,f), getfield(S2,f))
  end
end

# test here to track +meta logging
printstyled("    logging to :notes\n", color=:light_green)
for i in 1:S1.n
  @test any([occursin(fname, n) for n in S1.notes[i]])
end
for i in 1:S2.n
  @test any([occursin(fname, n) for n in S2.notes[i]])
end
for i in 1:S3.n
  @test any([occursin(dataless_name, n) for n in S3.notes[i]])
end

printstyled("    array of inputs\n", color=:light_green)
S = read_meta("sacpz", [sacpz_file, sacpz_wc])


text:
using FastTransforms, Makie, Random

include("SLPDE.jl")
include("etdrk4.jl")
include("evaluate_lambda.jl")
include("sphcesaro.jl")

Random.seed!(0)

f = 0.8
ε = 0.075
τ = 7.8125
E = 4.0

ue = ε^2*E/(1-f)*sqrt(4π)
ve = (1-f)/(E*ε^2)*sqrt(4π)

n = 64
U0 = zeros(4n, 8n-1), zeros(4n, 8n-1)
U0[1][1:n, 1:2n-1] = sphrandn(Float64, n, 2n-1)/(π*n)*ue/100
U0[2][1:n, 1:2n-1] = sphrandn(Float64, n, 2n-1)/(π*n)*ve/100
U0[1][1] += ue
U0[2][1] += ve

n = size(U0[1], 1)

θ = [0;(0.5:n-0.5)/n;1]
φ = [(0:2n-2)*2/(2n-1);2]
x = Float32[cospi(φ)*sinpi(θ) for θ in θ, φ in φ]
y = Float32[sinpi(φ)*sinpi(θ) for θ in θ, φ in φ]
z = Float32[cospi(θ) for θ in θ, φ in φ]

P = plan_sph2fourier(U0[1])
PS = plan_sph_synthesis(U0[1])

F = deepcopy(U0[1])
Ft = zeros(Float32, size(F, 1)+2, size(F, 2)+1)
lmul!(P, F)
lmul!(PS, F)
fillF!(Ft, F)
G = deepcopy(U0[2])
Gt = zeros(Float32, size(G, 1)+2, size(G, 2)+1)
lmul!(P, G)
lmul!(PS, G)
fillF!(Gt, G)

println("These are the extrema of the first component: ", extrema(F))
println(minimum(F)/(ue/sqrt(4π)),"  ",maximum(F)/(ue/sqrt(4π)))
scene = Scene(resolution = (1200, 1200));
surf = surface!(scene, x, y, z, color = Ft, colormap = :viridis, colorrange = extrema(F))[end];
update_cam!(scene, Vec3f0(2.5), Vec3f0(0), Vec3f0(0, 0, 1))
scene.center = false
scene
Makie.save("plots/LBRu0.jpeg", scene)
Makie.save("plots/NBRu0.jpeg", scene)

println("These are the extrema of the second component: ", extrema(G))
println(minimum(G)/(ve/sqrt(4π)),"  ",maximum(G)/(ve/sqrt(4π)))
scene = Scene(resolution = (1200, 1200));
surf = surface!(scene, x, y, z, color = Gt, colormap = :viridis, colorrange = extrema(G))[end];
update_cam!(scene, Vec3f0(2.5), Vec3f0(0), Vec3f0(0, 0, 1))
scene.center = false
scene
Makie.save("plots/LBRv0.jpeg", scene)
Makie.save("plots/NBRv0.jpeg", scene)

T = 16.0
n = round(Int, 10*T)

SLPDE = SemiLinearPDE((LaplaceBeltrami(ε),LaplaceBeltrami(inv(sqrt(τ)))), (NonlinearOperator((u,v)->ε^2*E-u+f*u^2*v),NonlinearOperator((u,v)->(u-u^2*v)/(τ*ε^2))), U0)
io = VideoStream(scene);
ETDRK4(SLPDE, T, n, surf, io; colorrange=:notfixed)
Makie.save("plots/LBR.gif", io)

SLPDE = SemiLinearPDE((NonlocalLaplaceBeltrami(ε, 0.0, 1.0),NonlocalLaplaceBeltrami(inv(sqrt(τ)), 0.0, 1.0)), (NonlinearOperator((u,v)->ε^2*E-u+f*u^2*v),NonlinearOperator((u,v)->(u-u^2*v)/(τ*ε^2))), U0)
io = VideoStream(scene);
ETDRK4(SLPDE, T, n, surf, io; colorrange=:notfixed)
Makie.save("plots/NBR.gif", io)


text:
#=
deepCondense : remove empty mapping of row, column, and value, and return the condensed of the input Assoc.
=#
function deepCondense(A::Assoc)
    Anew = condense(A)
    row,col,val = findnz(Anew.A)
    uniVal = sort!(unique(val))
    val = Array{Int64,1}(pmap(x -> searchsortedfirst(uniVal,x), val))
    #At this point val is the mapping to uniVal
    Anew.A = sparse(row,col,val)

    if A.val == [1.0] #Checking if the A.val mapping needs to be done.
        Anew.val = Array{Union{AbstractString,Number},1}(uniVal)
    else
        uniVal = pmap(x -> A.newval[x],uniVal)
        Anew.val = Array{Union{AbstractString,Number},1}(uniVal)
    end
    return Anew

end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
module ChainedFixes

@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end ChainedFixes

using Base: Fix1, Fix2, tail

export @nfix, and, ⩓, or, pipe_chain, ⩔

include("chained_core.jl")
using .ChainedCore

end # module


text:
using DistStat, Random, Test, LinearAlgebra

type=[Float32,Float64]

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    using CUDA
    ArrayType = CuArray
else
    ArrayType = Array
end

for T in type
    A=ArrayType{T}(reshape(collect(1:36),4,9))
    B=ArrayType{T}(reshape(collect(-7:28),4,9))

    A_dist=distribute(A); B_dist=distribute(B)

    @test isapprox(LinearAlgebra.dot(A_dist,B_dist),LinearAlgebra.dot(A,B))

    A_vec=vec(A); B_vec=vec(B)

    @test isapprox(LinearAlgebra.dot(A_dist,B_dist),LinearAlgebra.dot(A_vec,B_vec))

end


text:
# cfitsio.jl - C-style interface to CFITSIO functions:
#
# - Function names closely mirror the C interface (e.g., `fits_open_file()`).
# - Functions operate on `FITSFile`, a thin wrapper for `fitsfile` C struct
#   (`FITSFile` has concept of "current HDU", as in CFITSIO).
# - Note that the wrapper functions *do* check the return status from CFITSIO
#   and throw an error with the appropriate message.
#
#
# Syntax compatibility notes:
#
# - `convert(Vector{T}, x)` can be changed to `Vector{T}(x)` once
#   v0.3 is no longer supported, or can be changed to
#   `@compat Vector{T}(x)` once syntax support is in Compat.
# - `convert(Cint, x)` can be changed to `Cint(x)` once v0.3 is not supported
#   or `@compat Cint(x)` once syntax support is in Compat.
#
# The following table gives the correspondances between CFITSIO "types",
# the BITPIX keyword and Julia types.
#
#     -------------------------------------------------
#     CODE  CFISTIO         Julia     Comments
#     -------------------------------------------------
#           int             Cint
#           long            Clong
#           LONGLONG        Int64     64-bit integer
#     -------------------------------------------------
#     -------- FITS BITPIX ----------------------------
#        8  BYTE_IMG        Uint8
#       16  SHORT_IMG       Int16
#       32  LONG_IMG        Int32
#       64  LONGLONG_IMG    Int64
#      -32  FLOAT_IMG       Float32
#      -64  DOUBLE_IMG      Float64
#      -------- cfitsio "aliases" ---------------------
#       10  SBYTE_IMG       Int8     written as: BITPIX = 8, BSCALE = 1,
#                                                BZERO = -128
#       20  USHORT_IMG      Uint16   written as: BITPIX = 16, BSCALE = 1,
#                                                BZERO = 32768
#       40  LONG_IMG        Uint32   written as: BITPIX = 32, BSCALE = 1,
#                                                BZERO = 2147483648
#     -------------------------------------------------
#     -------- FITS TABLE DATA TYPES ------------------
#        1  TBIT
#       11  TBYTE           Cuchar = Uint8
#       12  TSBYTE          Cchar = Int8
#       14  TLOGICAL        Bool
#       16  TSTRING         Compat.ASCIIString
#       20  TUSHORT         Cushort
#       21  TSHORT          Cshort
#       30  TUINT           Cuint
#       31  TINT            Cint
#       40  TULONG          Culong
#       41  TLONG           Clong
#       42  TFLOAT          Cfloat
#       81  TLONGLONG       Int64
#       82  TDOUBLE         Cdouble
#       83  TCOMPLEX        Complex{Cfloat}
#      163  TDBLCOMPLEX     Complex{Cdouble}
#     -------------------------------------------------
#

isdefined(Base, :__precompile__) && __precompile__()

module Libcfitsio

using Compat

export FITSFile,
       fits_assert_open,
       fits_clobber_file,
       fits_close_file,
       fits_copy_image_section,
       fits_create_ascii_tbl,
       fits_create_binary_tbl,
       fits_create_file,
       fits_create_img,
       fits_delete_file,
       fits_delete_key,
       fits_delete_record,
       fits_delete_rows,
       fits_file_mode,
       fits_file_name,
       fits_get_hdrspace,
       fits_get_hdu_num,
       fits_get_hdu_type,
       fits_get_img_dim,
       fits_get_img_equivtype,
       fits_get_img_size,
       fits_get_img_type,
       fits_get_num_cols,
       fits_get_num_hdus,
       fits_get_num_rows,
       fits_get_rowsize,
       fits_get_colnum,
       fits_get_coltype,
       fits_get_eqcoltype,
       fits_get_version,
       fits_read_tdim,
       fits_hdr2str,
       fits_insert_rows,
       fits_movabs_hdu,
       fits_movrel_hdu,
       fits_movnam_hdu,
       fits_open_data,
       fits_open_file,
       fits_open_image,
       fits_open_table,
       fits_read_col,
       fits_read_descript,
       fits_read_keyn,
       fits_read_keyword,
       fits_read_pix,
       fits_read_record,
       fits_read_subset,
       fits_update_key,
       fits_write_col,
       fits_write_comment,
       fits_write_history,
       fits_write_key,
       fits_write_pix,
       fits_write_record,
       fits_write_tdim

if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("FITSIO not properly installed. Please run Pkg.build(\"FITSIO\")")
end

const TYPE_FROM_BITPIX = Dict{Cint, DataType}()
for (T, code) in ((UInt8,     8), # BYTE_IMG
                  (Int16,    16), # SHORT_IMG
                  (Int32,    32), # LONG_IMG
                  (Int64,    64), # LONGLONG_IMG
                  (Float32, -32), # FLOAT_IMG
                  (Float64, -64), # DOUBLE_IMG
                  (Int8,     10), # SBYTE_IMG
                  (UInt16,   20), # USHORT_IMG
                  (UInt32,   40)) # ULONG_IMG
    local value = convert(Cint, code)
    @eval begin
        TYPE_FROM_BITPIX[$value] = $T
        bitpix_from_type(::Type{$T}) = $value
    end
end

for (T, code) in ((UInt8,       11),
                  (Int8,        12),
                  (Bool,        14),
                  (Compat.ASCIIString, 16),
                  (Cushort,     20),
                  (Cshort,      21),
                  (Cuint,       30),
                  (Cint,        31),
                  (Int64,       81),
                  (Float32,     42),
                  (Float64,     82),
                  (Complex64,   83),
                  (Complex128, 163))
    @eval cfitsio_typecode(::Type{$T}) = convert(Cint, $code)
end

# Above, we don't define a method for Clong because it is either Cint (Int32)
# or Int64 depending on the platform, and those methods are already defined.
# Culong is either UInt64 or Cuint depending on platform. Only define it if
# not already defined.
if Culong !== Cuint
    cfitsio_typecode(::Type{Culong}) = convert(Cint, 40)
end

# -----------------------------------------------------------------------------
# FITSFile type

type FITSFile
    ptr::Ptr{Void}

    function FITSFile(ptr::Ptr{Void})
        f = new(ptr)
        finalizer(f, fits_close_file)
        f
    end
end

# -----------------------------------------------------------------------------
# error messaging

function fits_assert_open(f::FITSFile)
    if f.ptr == C_NULL
        error("attempt to access closed FITS file")
    end
end

function fits_get_errstatus(status::Cint)
    msg = Vector{UInt8}(31)
    ccall((:ffgerr,libcfitsio), Void, (Cint,Ptr{UInt8}), status, msg)
    unsafe_string(pointer(msg))
end

function fits_assert_ok(status::Cint)
    if status != 0
        error(fits_get_errstatus(status))
    end
end

fits_get_version() = ccall((:ffvers, libcfitsio), Cfloat, (Ptr{Cfloat},), &0.)

# -----------------------------------------------------------------------------
# file access & info functions

function fits_create_file(filename::AbstractString)
    ptr = Ref{Ptr{Void}}()
    status = Ref{Cint}(0)
    ccall((:ffinit,libcfitsio), Cint, (Ref{Ptr{Void}},Ptr{UInt8},Ref{Cint}),
          ptr, filename, status)
    fits_assert_ok(status[])
    FITSFile(ptr[])
end

fits_clobber_file(filename::AbstractString) = fits_create_file("!"*filename)

for (a,b) in ((:fits_open_data, "ffdopn"),
              (:fits_open_file, "ffopen"),
              (:fits_open_image,"ffiopn"),
              (:fits_open_table,"fftopn"))
    @eval begin
        function ($a)(filename::AbstractString, mode::Integer=0)
            ptr = Ref{Ptr{Void}}()
            status = Ref{Cint}(0)
            ccall(($b,libcfitsio), Cint,
                  (Ref{Ptr{Void}},Ptr{UInt8},Cint,Ref{Cint}),
                  ptr, filename, mode, status)
            fits_assert_ok(status[])
            FITSFile(ptr[])
        end
    end
end

for (a,b) in ((:fits_close_file, "ffclos"),
              (:fits_delete_file,"ffdelt"))
    @eval begin
        function ($a)(f::FITSFile)

            # fits_close_file() is called during garbage collection, but file
            # may already be closed by user, so we need to check if it is open.
            if f.ptr != C_NULL
                status = Ref{Cint}(0)
                ccall(($b,libcfitsio), Cint,
                      (Ptr{Void},Ref{Cint}),
                      f.ptr, status)
                fits_assert_ok(status[])
                f.ptr = C_NULL
            end
        end
    end
end

close(f::FITSFile) = fits_close_file(f)

function fits_file_name(f::FITSFile)
    value = Vector{UInt8}(1025)
    status = Ref{Cint}(0)
    ccall((:ffflnm,libcfitsio), Cint,
          (Ptr{Void},Ptr{UInt8},Ref{Cint}),
          f.ptr, value, status)
    fits_assert_ok(status[])
    unsafe_string(pointer(value))
end

function fits_file_mode(f::FITSFile)
    result = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall(("ffflmd", libcfitsio), Cint, (Ptr{Void}, Ref{Cint}, Ref{Cint}),
          f.ptr, result, status)
    fits_assert_ok(status[])
    result[]
end

# -----------------------------------------------------------------------------
# header access functions

function fits_get_hdrspace(f::FITSFile)
    keysexist = Ref{Cint}(0)
    morekeys = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffghsp,libcfitsio), Cint,
        (Ptr{Void},Ref{Cint},Ref{Cint},Ref{Cint}),
        f.ptr, keysexist, morekeys, status)
    fits_assert_ok(status[])
    (keysexist[], morekeys[])
end

function fits_read_keyword(f::FITSFile, keyname::Compat.ASCIIString)
    value = Vector{UInt8}(71)
    comment = Vector{UInt8}(71)
    status = Ref{Cint}(0)
    ccall((:ffgkey,libcfitsio), Cint,
        (Ptr{Void},Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ref{Cint}),
        f.ptr, keyname, value, comment, status)
    fits_assert_ok(status[])
    unsafe_string(pointer(value)), unsafe_string(pointer(comment))
end

function fits_read_record(f::FITSFile, keynum::Integer)
    card = Vector{UInt8}(81)
    status = Ref{Cint}(0)
    ccall((:ffgrec,libcfitsio), Cint,
        (Ptr{Void},Cint,Ptr{UInt8},Ref{Cint}),
        f.ptr, keynum, card, status)
    fits_assert_ok(status[])
    unsafe_string(pointer(card))
end

function fits_read_keyn(f::FITSFile, keynum::Integer)
    keyname = Vector{UInt8}(9)
    value = Vector{UInt8}(71)
    comment = Vector{UInt8}(71)
    status = Ref{Cint}(0)
    ccall((:ffgkyn,libcfitsio), Cint,
        (Ptr{Void},Cint,Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ref{Cint}),
        f.ptr, keynum, keyname, value, comment, status)
    fits_assert_ok(status[])
    (unsafe_string(pointer(keyname)), unsafe_string(pointer(value)),
     unsafe_string(pointer(comment)))
end

function fits_write_key(f::FITSFile, keyname::Compat.ASCIIString,
                        value::@compat(Union{AbstractFloat,Compat.ASCIIString}),
                        comment::Compat.ASCIIString)
    cvalue = isa(value,Compat.ASCIIString) ?  value :
             isa(value,Bool) ? Cint[value] : [value]
    status = Ref{Cint}(0)
    ccall((:ffpky,libcfitsio), Cint,
        (Ptr{Void},Cint,Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ref{Cint}),
        f.ptr, cfitsio_typecode(typeof(value)), keyname,
        cvalue, comment, status)
    fits_assert_ok(status[])
end

function fits_write_comment(f::FITSFile, comment::Compat.ASCIIString)
    status = Ref{Cint}(0)
    ccall((:ffpcom, libcfitsio), Cint, (Ptr{Void}, Ptr{UInt8}, Ref{Cint}),
          f.ptr, comment, status)
    fits_assert_ok(status[])
end

function fits_write_history(f::FITSFile, history::Compat.ASCIIString)
    status = Ref{Cint}(0)
    ccall((:ffphis, libcfitsio), Cint, (Ptr{Void}, Ptr{UInt8}, Ref{Cint}),
          f.ptr, history, status)
    fits_assert_ok(status[])
end

# update key: if already present, update it, otherwise add it.
for (a,T,S) in (("ffukys", :(Compat.ASCIIString), :(Ptr{UInt8})),
                ("ffukyl", :Bool,        :Cint),
                ("ffukyj", :Integer,     :Int64))
    @eval begin
        function fits_update_key(f::FITSFile, key::Compat.ASCIIString, value::$T,
                                 comment::@compat(Union{Compat.ASCIIString, Ptr{Void}})=C_NULL)
            status = Ref{Cint}(0)
            ccall(($a, libcfitsio), Cint,
                  (Ptr{Void}, Ptr{UInt8}, $S, Ptr{UInt8}, Ref{Cint}),
                  f.ptr, key, value, comment, status)
            fits_assert_ok(status[])
        end
    end
end

function fits_update_key(f::FITSFile, key::Compat.ASCIIString, value::AbstractFloat,
                         comment::@compat(Union{Compat.ASCIIString, Ptr{Void}})=C_NULL)
    status = Ref{Cint}(0)
    ccall(("ffukyd", libcfitsio), Cint,
          (Ptr{Void}, Ptr{UInt8}, Cdouble, Cint, Ptr{UInt8}, Ref{Cint}),
          f.ptr, key, value, -15, comment, status)
    fits_assert_ok(status[])
end

function fits_update_key(f::FITSFile, key::Compat.ASCIIString, value::@compat(Void),
                         comment::@compat(Union{Compat.ASCIIString, Ptr{Void}})=C_NULL)
    status = Ref{Cint}(0)
    ccall(("ffukyu", libcfitsio), Cint,
          (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}, Ref{Cint}),
          f.ptr, key, comment, status)
    fits_assert_ok(status[])
end

function fits_write_record(f::FITSFile, card::Compat.ASCIIString)
    status = Ref{Cint}(0)
    ccall((:ffprec,libcfitsio), Cint,
        (Ptr{Void},Ptr{UInt8},Ref{Cint}),
        f.ptr, card, status)
    fits_assert_ok(status[])
end

function fits_delete_record(f::FITSFile, keynum::Integer)
    status = Ref{Cint}(0)
    ccall((:ffdrec,libcfitsio), Cint,
        (Ptr{Void},Cint,Ref{Cint}),
        f.ptr, keynum, status)
    fits_assert_ok(status[])
end

function fits_delete_key(f::FITSFile, keyname::Compat.ASCIIString)
    status = Ref{Cint}(0)
    ccall((:ffdkey,libcfitsio), Cint,
        (Ptr{Void},Ptr{UInt8},Ref{Cint}),
        f.ptr, keyname, status)
    fits_assert_ok(status[])
end

function fits_hdr2str(f::FITSFile, nocomments::Bool=false)
    status = Ref{Cint}(0)
    header = Ref{Ptr{UInt8}}()
    nkeys = Ref{Cint}(0)
    ccall((:ffhdr2str, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Ptr{UInt8}}, Cint,
           Ptr{Ptr{UInt8}}, Ref{Cint}, Ref{Cint}),
          f.ptr, nocomments, &C_NULL, 0, header, nkeys, status)
    result = unsafe_string(header[])

    # free header pointer allocated by cfitsio (result is a copy)
    ccall((:fffree, libcfitsio), Ref{Cint},
          (Ptr{UInt8}, Ref{Cint}),
          header[], status)
    fits_assert_ok(status[])
    result
end

# -----------------------------------------------------------------------------
# HDU info functions and moving the current HDU

function hdu_int_to_type(hdu_type_int)
    if hdu_type_int == 0
        return :image_hdu
    elseif hdu_type_int == 1
        return :ascii_table
    elseif hdu_type_int == 2
        return :binary_table
    end

    :unknown
end

for (a,b) in ((:fits_movabs_hdu,"ffmahd"),
              (:fits_movrel_hdu,"ffmrhd"))
    @eval begin
        function ($a)(f::FITSFile, hduNum::Integer)
            hdu_type = Ref{Cint}(0)
            status = Ref{Cint}(0)
            ccall(($b,libcfitsio), Cint,
                  (Ptr{Void}, Cint, Ref{Cint}, Ref{Cint}),
                  f.ptr, hduNum, hdu_type, status)
            fits_assert_ok(status[])
            hdu_int_to_type(hdu_type[])
        end
    end
end

function fits_movnam_hdu(f::FITSFile, extname::Compat.ASCIIString, extver::Integer=0,
                         hdu_type::Integer=-1)
    status = Ref{Cint}(0)
    ccall((:ffmnhd,libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{UInt8}, Cint, Ref{Cint}),
          f.ptr, hdu_type, extname, extver, status)
    fits_assert_ok(status[])
end

function fits_get_hdu_num(f::FITSFile)
    hdunum = Ref{Cint}(0)
    ccall((:ffghdn,libcfitsio), Cint,
          (Ptr{Void}, Ref{Cint}),
          f.ptr, hdunum)
    hdunum[]
end

function fits_get_hdu_type(f::FITSFile)
    hdutype = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffghdt, libcfitsio), Cint,
          (Ptr{Void}, Ref{Cint}, Ref{Cint}),
          f.ptr, hdutype, status)
    fits_assert_ok(status[])
    hdu_int_to_type(hdutype[])
end

# -----------------------------------------------------------------------------
# image HDU functions

for (a, b) in ((:fits_get_img_type,      "ffgidt"),
               (:fits_get_img_equivtype, "ffgiet"),
               (:fits_get_img_dim,       "ffgidm"))
    @eval function ($a)(f::FITSFile)
        result = Ref{Cint}(0)
        status = Ref{Cint}(0)
        ccall(($b, libcfitsio), Cint, (Ptr{Void}, Ref{Cint}, Ref{Cint}),
              f.ptr, result, status)
        fits_assert_ok(status[])
        result[]
    end
end

function fits_create_img{T, S<:Integer}(f::FITSFile, ::Type{T},
                                        naxes::Vector{S})
    status = Ref{Cint}(0)
    ccall((:ffcrimll, libcfitsio), Cint,
          (Ptr{Void}, Cint, Cint, Ptr{Int64}, Ref{Cint}),
          f.ptr, bitpix_from_type(T), length(naxes),
          convert(Vector{Int64}, naxes), status)
    fits_assert_ok(status[])
end

function fits_write_pix{S<:Integer,T}(f::FITSFile, fpixel::Vector{S},
                                      nelements::Integer, data::Array{T})
    status = Ref{Cint}(0)
    ccall((:ffppxll, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Int64}, Int64, Ptr{Void}, Ref{Cint}),
          f.ptr, cfitsio_typecode(T), convert(Vector{Int64}, fpixel),
          nelements, data, status)
    fits_assert_ok(status[])
end

function fits_write_pix(f::FITSFile, data::Array)
    fits_write_pix(f, ones(Int64, length(size(data))), length(data), data)
end

function fits_read_pix{S<:Integer,T}(f::FITSFile, fpixel::Vector{S},
                                     nelements::Int, nullval::T,
                                     data::Array{T})
    anynull = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffgpxvll, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Int64}, Int64, Ptr{Void}, Ptr{Void},
           Ref{Cint}, Ref{Cint}),
          f.ptr, cfitsio_typecode(T), convert(Vector{Int64}, fpixel),
          nelements, &nullval, data, anynull, status)
    fits_assert_ok(status[])
    anynull[]
end

function fits_read_pix{S<:Integer,T}(f::FITSFile, fpixel::Vector{S},
                                     nelements::Int, data::Array{T})
    anynull = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffgpxvll, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Int64}, Int64, Ptr{Void}, Ptr{Void},
           Ref{Cint}, Ref{Cint}),
          f.ptr, cfitsio_typecode(T), convert(Vector{Int64}, fpixel),
          nelements, C_NULL, data, anynull, status)
    fits_assert_ok(status[])
    anynull[]
end

function fits_read_pix(f::FITSFile, data::Array)
    fits_read_pix(f, ones(Int64,length(size(data))), length(data), data)
end

function fits_read_subset{S1<:Integer,S2<:Integer,S3<:Integer,T}(
             f::FITSFile, fpixel::Vector{S1}, lpixel::Vector{S2},
             inc::Vector{S3}, data::Array{T})
    anynull = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffgsv, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Void},
           Ptr{Void}, Ref{Cint}, Ref{Cint}),
          f.ptr, cfitsio_typecode(T),
          convert(Vector{Clong}, fpixel),
          convert(Vector{Clong}, lpixel),
          convert(Vector{Clong}, inc),
          C_NULL, data, anynull, status)
    fits_assert_ok(status[])
    anynull[]
end

function fits_copy_image_section(fin::FITSFile, fout::FITSFile,
                                 section::Compat.ASCIIString)
    status = Ref{Cint}(0)
    ccall((:fits_copy_image_section,libcfitsio), Cint,
          (Ptr{Void}, Ptr{Void}, Ptr{UInt8}, Ref{Cint}),
          fin.ptr, fout.ptr, section, status)
    fits_assert_ok(status[])
end

# -----------------------------------------------------------------------------
# ASCII/binary table HDU functions

# The three fields are: ttype, tform, tunit (CFITSIO's terminology)
typealias ColumnDef @compat Tuple{Compat.ASCIIString, Compat.ASCIIString, Compat.ASCIIString}

for (a,b) in ((:fits_create_binary_tbl, 2),
              (:fits_create_ascii_tbl,  1))
    @eval begin
        function ($a)(f::FITSFile, numrows::Integer,
                      coldefs::Array{ColumnDef}, extname::Compat.ASCIIString)

            # get length and convert coldefs to three arrays of Ptr{Uint8}
            ntype = length(coldefs)
            ttype = [pointer(x[1]) for x in coldefs]
            tform = [pointer(x[2]) for x in coldefs]
            tunit = [pointer(x[3]) for x in coldefs]
            status = Ref{Cint}(0)

            ccall(("ffcrtb", libcfitsio), Cint,
                  (Ptr{Void}, Cint, Int64, Cint,
                   Ptr{Ptr{UInt8}}, Ptr{Ptr{UInt8}},
                   Ptr{Ptr{UInt8}}, Ptr{UInt8}, Ref{Cint}),
                  f.ptr, $b, numrows, ntype,
                  ttype, tform, tunit, extname,
                  status)
            fits_assert_ok(status[])
        end
    end
end

for (a,b,T) in ((:fits_get_num_cols,  "ffgncl",  :Cint),
                (:fits_get_num_hdus,  "ffthdu",  :Cint),
                (:fits_get_rowsize,   "ffgrsz",  :Clong))
    @eval begin
        function ($a)(f::FITSFile)
            result = Ref{$T}(0)
            status = Ref{Cint}(0)
            ccall(($b,libcfitsio), Cint,
                  (Ptr{Void}, Ref{$T}, Ref{Cint}),
                  f.ptr, result, status)
            fits_assert_ok(status[])
            result[]
        end
    end
end

function fits_get_colnum(f::FITSFile, tmplt::Compat.ASCIIString)
    result = Ref{Cint}(0)
    status = Ref{Cint}(0)

    # Second argument is case-sensitivity of search: 0 = case-insensitive
    #                                                1 = case-sensitive
    ccall(("ffgcno", libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{UInt8}, Ref{Cint}, Ref{Cint}),
          f.ptr, 0, tmplt, result, status)
    fits_assert_ok(status[])
    return result[]
end

# The following block are all functions that have separate variants for Clong
# and 64-bit integers in cfitsio. Rather than providing both of these, we
# provide only one according to the native integer type on the platform.
if  promote_type(Int, Clong) == Clong
    T = Clong
    ffgtdm = "ffgtdm"
    ffgnrw = "ffgnrw"
    ffptdm = "ffptdm"
    ffgtcl = "ffgtcl"
    ffeqty = "ffeqty"
    ffgdes = "ffgdes"
    ffgisz = "ffgisz"
else
    T = Int64
    ffgtdm = "ffgtdmll"
    ffgnrw = "ffgnrwll"
    ffptdm = "ffptdmll"
    ffgtcl = "ffgtclll"
    ffeqty = "ffeqtyll"
    ffgdes = "ffgdesll"
    ffgisz = "ffgiszll"
end
@eval begin
    function fits_get_coltype(ff::FITSFile, colnum::Integer)
        typecode = Ref{Cint}(0)
        repcnt = Ref{$T}(0)
        width = Ref{$T}(0)
        status = Ref{Cint}(0)
        ccall(($ffgtcl,libcfitsio), Cint,
              (Ptr{Void}, Cint, Ref{Cint}, Ref{$T}, Ref{$T}, Ref{Cint}),
              ff.ptr, colnum, typecode, repcnt, width, status)
        fits_assert_ok(status[])
        return @compat Int(typecode[]), Int(repcnt[]), Int(width[])
    end

    function fits_get_eqcoltype(ff::FITSFile, colnum::Integer)
        typecode = Ref{Cint}(0)
        repcnt = Ref{$T}(0)
        width = Ref{$T}(0)
        status = Ref{Cint}(0)
        ccall(($ffeqty,libcfitsio), Cint,
              (Ptr{Void}, Cint, Ref{Cint}, Ref{$T}, Ref{$T}, Ref{Cint}),
              ff.ptr, colnum, typecode, repcnt, width, status)
        fits_assert_ok(status[])
        return @compat Int(typecode[]), Int(repcnt[]), Int(width[])
    end

    function fits_get_img_size(f::FITSFile)
        ndim = fits_get_img_dim(f)
        naxes = Vector{$T}(ndim)
        status = Ref{Cint}(0)
        ccall(($ffgisz, libcfitsio), Cint,
              (Ptr{Void}, Cint, Ptr{$T}, Ref{Cint}),
              f.ptr, ndim, naxes, status)
        fits_assert_ok(status[])
        naxes
    end

    function fits_get_num_rows(f::FITSFile)
        result = Ref{$T}(0)
        status = Ref{Cint}(0)
        ccall(($ffgnrw, libcfitsio), Cint,
              (Ptr{Void}, Ref{$T}, Ref{Cint}),
              f.ptr, result, status)
        fits_assert_ok(status[])
        return @compat Int(result[])
    end

    # `fits_read_tdim` returns the dimensions of a table column in a
    # binary table. Normally this information is given by the TDIMn
    # keyword, but if this keyword is not present then this routine
    # returns `[r]` with `r` equals to the repeat count in the TFORM
    # keyword.
    function fits_read_tdim(ff::FITSFile, colnum::Integer)
        naxes = Vector{$T}(99)  # 99 is the maximum allowed number of axes
        naxis = Ref{Cint}(0)
        status = Ref{Cint}(0)
        ccall(($ffgtdm,libcfitsio), Cint,
              (Ptr{Void}, Cint, Cint, Ref{Cint}, Ptr{$T}, Ref{Cint}),
              ff.ptr, colnum, length(naxes), naxis, naxes, status)
        fits_assert_ok(status[])
        return naxes[1:naxis[]]
    end

    function fits_write_tdim(ff::FITSFile, colnum::Integer,
                                 naxes::Array{$T})
        status = Ref{Cint}(0)
        ccall(($ffptdm, libcfitsio), Cint,
              (Ptr{Void}, Cint, Cint, Ptr{$T}, Ref{Cint}),
              ff.ptr, colnum, length(naxes), naxes, status)
        fits_assert_ok(status[])
    end

    function fits_read_descript(f::FITSFile, colnum::Integer, rownum::Integer)
        repeat = Ref{$T}(0)
        offset = Ref{$T}(0)
        status = Ref{Cint}(0)
        ccall(($ffgdes, libcfitsio), Cint,
              (Ptr{Void}, Cint, Int64, Ref{$T}, Ref{$T}, Ref{Cint}),
              f.ptr, colnum, rownum, repeat, offset, status)
        fits_assert_ok(status[])
        return @compat Int(repeat[]), Int(offset[])
    end
end

function fits_read_col(f::FITSFile,
                       colnum::Integer,
                       firstrow::Integer,
                       firstelem::Integer,
                       data::Array{Compat.ASCIIString})

    # get width: number of characters in each string
    typecode, repcount, width = fits_get_eqcoltype(f, colnum)

    # ensure that data are strings, otherwise cfitsio will try to write
    # formatted strings, which have widths given by fits_get_col_display_width
    # not by the repeat value from fits_get_coltype.
    abs(typecode) == 16 || error("not a string column")

    # create an array of character buffers of the correct width
    buffers = [Vector{UInt8}(width) for i in 1:length(data)]

    # Call the CFITSIO function
    anynull = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffgcvs, libcfitsio), Cint,
          (Ptr{Void}, Cint, Int64, Int64, Int64,
           Ptr{UInt8}, Ptr{Ptr{UInt8}}, Ref{Cint}, Ref{Cint}),
          f.ptr, colnum, firstrow, firstelem, length(data),
          " ", buffers, anynull, status)
    fits_assert_ok(status[])

    # Create strings out of the buffers, terminating at null characters.
    # Note that `Compat.ASCIIString(x)` does not copy the buffer x.
    for i in 1:length(data)
        zeropos = search(buffers[i], 0x00)
        data[i] = (zeropos >= 1) ? Compat.ASCIIString(buffers[i][1:(zeropos-1)]) :
                                   Compat.ASCIIString(buffers[i])
    end
end

function fits_read_col{T}(f::FITSFile,
                          colnum::Integer,
                          firstrow::Integer,
                          firstelem::Integer,
                          data::Array{T})
    anynull = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffgcv,libcfitsio), Cint,
          (Ptr{Void}, Cint, Cint, Int64, Int64, Int64,
           Ptr{Void}, Ptr{Void}, Ref{Cint}, Ref{Cint}),
          f.ptr, cfitsio_typecode(T), colnum,
          firstrow, firstelem, length(data),
          C_NULL, data, anynull, status)
    fits_assert_ok(status[])
end

function fits_write_col(f::FITSFile,
                        colnum::Integer,
                        firstrow::Integer,
                        firstelem::Integer,
                        data::Array{Compat.ASCIIString})
    status = Ref{Cint}(0)
    ccall((:ffpcls, libcfitsio), Cint,
          (Ptr{Void}, Cint, Int64, Int64, Int64,
           Ptr{Ptr{UInt8}}, Ref{Cint}),
          f.ptr, colnum, firstrow, firstelem, length(data),
          data, status)
    fits_assert_ok(status[])
end

function fits_write_col{T}(f::FITSFile,
                           colnum::Integer,
                           firstrow::Integer,
                           firstelem::Integer,
                           data::Array{T})
    status = Ref{Cint}(0)
    ccall((:ffpcl, libcfitsio), Cint,
          (Ptr{Void}, Cint, Cint, Int64, Int64, Int64,
           Ptr{Void}, Ref{Cint}),
          f.ptr, cfitsio_typecode(T), colnum,
          firstrow, firstelem, length(data),
          data, status)
    fits_assert_ok(status[])
end

for (a,b) in ((:fits_insert_rows, "ffirow"),
              (:fits_delete_rows, "ffdrow"))
    @eval begin
        function ($a)(f::FITSFile, firstrow::Integer, nrows::Integer)
            status = Ref{Cint}(0)
            ccall(($b,libcfitsio), Cint,
                  (Ptr{Void}, Int64, Int64, Ref{Cint}),
                  f.ptr, firstrow, nrows, status)
            fits_assert_ok(status[])
        end
    end
end

# -----------------------------------------------------------------------------
# deprecations

import Base: @deprecate

# Deprecated in v0.7

@deprecate fits_get_num_rowsll fits_get_num_rows

# Deprecated in v0.6

@deprecate fits_get_col_repeat(f::FITSFile, colnum::Integer) fits_get_coltype(f, colnum)[2, 3]

@deprecate fits_read_col{T}(f::FITSFile, ::Type{T}, colnum::Integer, firstrow::Integer, firstelem::Integer, data::Array{T}) fits_read_col(f, colnum, firstrow, firstelem, data)

@deprecate fits_write_col{T}(f::FITSFile, ::Type{T}, colnum::Integer, firstrow::Integer, firstelem::Integer, data::Array{T}) fits_write_col(f, colnum, firstrow, firstelem, data)

end # module


text:
using Pkg
Pkg.activate(".")
using EquivalenceClassesConstructor
using Printf

const nBose = 80
const nFermi = 80
const shift = 0

function symm_pred(l::Tuple{Int64,Int64,Int64}, r::Tuple{Int64,Int64,Int64})
    res = (l .== .-r) || (l[1] == r[1] && l[2] == r[3] && l[3] == r[2]) || # complex conjugation and time reversal
    ((l[1] == r[2] + r[3]) && l[2] == r[2] && (l[3] == r[1] + r[2])) ||    # single crossing symmetry (creation)
    ((r[1] == l[2] + l[3]) && r[2] == l[2] && (r[3] == l[1] + l[2])) ||    #     ---  "  ---
    ((l[1] == r[3] - r[2]) && (l[2] == r[1] + r[3]) && l[3] == r[3]) ||    # single crossing symmetry (annihilation)
    ((r[1] == l[3] - l[2]) && (r[2] == l[1] + l[3]) && r[3] == l[3]) ||    #     ---  "  ---
    (r[1] == -l[1] && (l[2] == r[1] + r[3]) && (l[3] == r[1] + r[2])) ||   # double crossing
    (l[1] == -r[1] && (r[2] == l[1] + l[3]) && (r[3] == l[1] + l[2]))     #     ---  "  ---
    return res
end

@inline function symm_map(t::Tuple{Int64,Int64,Int64})::Array{Tuple{Int64,Int64,Int64},1}
    res = [.- t .- (0,1,1), (t[1], t[3], t[2]),        # complex conjugation and time reversal
           (t[3]-t[2], t[2], t[1]+t[2]),    # single crossing (creation)
           (t[2]-t[3], t[1] + t[3], t[3]), #,  # single crossing (annihilation)
           (-t[1],t[1]+t[3], t[1] + t[2])]  # double crossing
end

println("Constructing Array")
const freqList = [(i,j,k) for i in (-nBose:nBose) for j in (-nFermi:nFermi-1) .- trunc(Int64,shift*i/2) for k in (-nFermi:nFermi-1) .- trunc(Int64,shift*i/2)]
const len_freq = (2*nBose+1)*(2*nFermi)^2
const mm = Mapping(symm_map)
println("Starting Computation 2")
maxF = nBose + 2*nFermi + 5
headerstr= @sprintf("  %26d  \n", maxF)
@time cl = find_classes(mm, freqList, vl_len=len_freq)
@time em = ExpandMap(cl)
#write_fixed_width("freqList.dat", em, sorted=true, header_add=headerstr)
write_JLD("freqList.jld", em_2, rm_2, cl_2)


text:
using PyPlot, MRIReco

##### simple example ####

N = 256
I = shepp_logan(N)

# simulation parameters
params = Dict{Symbol, Any}()
params[:simulation] = "fast"
params[:trajName] = "Radial"
params[:numProfiles] = floor(Int64, pi/2*N)
params[:numSamplingPerProfile] = 2*N

# do simulation
acqData = simulation(I, params)

# reco parameters
params = Dict{Symbol, Any}()
params[:reco] = "direct"
params[:reconSize] = (N,N)
Ireco = reconstruction(acqData, params)

# use PyPlot for interactive display
figure(1)
clf()
subplot(1,2,1)
title("Phantom")
imshow(abs.(I))
subplot(1,2,2)
title("Reconstruction")
imshow(abs.(Ireco[:,:,1,1,1]))

# export images
filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/simpleReco.png")
exportImage(filename, abs.(Ireco[:,:,1,1,1]))
filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/phantom.png")
exportImage(filename, abs.(I[:,:,1,1,1]))


text:
"""
    Word{T,A}

Word corresponding to token `t` and assignment `a`.

# Examples
```
# Generate single words
Word(1,Int16(2))
Word("a",Int16(2))

# Generate a corpus (Int64,full loop)
c = []
k = 3 # number of topics
ndoc = 10 # number of documents
nlex = 100 # number of terms in lexicon
nread = 1000 # number of tokens per document
lexicon = [1:nlex;]
topics = [1:k;]
for i in 1:ndoc
    x = rand(lexicon,1000)
    y = rand(topics,1000)
    arr = Word.(x,y)
    push!(c,arr)
end

# Generate a corpus (Int8 topic, map)
corp = map(x->Word.(rand(1:100,x),Int8.(rand(Categorical([0.5,0.3,0.2]),x))),[10,50,100])
```
"""
struct Word{T,A}
  t::T # the index of the token in the lexicon
  a::A # the latent topic associated with the token
end

"""
  getCorp(corp)

Get the parametric representation of a corpus.

# Examples
```
# Generate a corpus (Int8 topic, map)
corp = map(x->Word.(rand(1:100,x),Int8.(rand(Categorical([0.5,0.3,0.2]),x))),[10,50,100])
c = getCorp(c)
```
"""
getCorp(corp) = Array{Array{Word{typeof(corp[1][1].t),typeof(corp[1][1].a)},1},1}(corp)

"""
  mat2docind(mat)

Given a matrix of counts, get the indices of the resulting row-wise doc vector for each lexicon term.

# Examples
```
mat2docind(mat)
```
"""
function mat2docind(mat)
  ends = cumsum(mat,dims=2)
  starts = hcat(ones(Int,size(mat)[1]),(cumsum(mat,dims=2) .+ 1)[:,1:end-1])
  [[starts[i],ends[i]] for i in CartesianIndices(mat)]
end

"""
  vec2procind(vlen,nproc)

Partition the indices of a vector of length `vlen` into `nproc` approximately equal groups.

# Examples
```
vec2procind(vlen,nproc)
```
"""
function vec2procind(vlen,nproc)
  pvec = mod.(1:vlen,nproc).+1
  [findall(x->x==i,pvec) for i in unique(pvec)]
end

"""
  initterms!(data,Tdoc,Adoc,datind,docind,k)

For each term corresponding to `datind` in a row (`data`) from a doc x terms matrix (as in gene expression data), draw initial assignments from `[1:k;]` for tokens (corresponding to `docind`) as `Adoc` and update the token indices in `Tdoc`.

# Examples
```
initterms!(data,Tdoc,Adoc,datind,docind,k)
```
"""
function initterms!(data,Tdoc,Adoc,datind,docind,k)
  for i in 1:length(datind)
    # println("initializing term:",datind[i])
    Tdoc[docind[i][1]:docind[i][2]] .= datind[i]
    Adoc[docind[i][1]:docind[i][2]] .= rand([1:k;],data[datind[i]])
  end
  # println("init complete")
end

"""
  initcorp(counts,inds,procinds,doclengths,k,nwork)

Initialize a corpus over the matrix `counts` using `nwork` available workers.

# Examples
```
initcorp(counts,inds,procinds,doclengths,k,nwork)
```
"""
function initcorp(counts,inds,procinds,doclengths,k,nwork)
  corp = []
  for d in 1:size(counts)[1]
    println("initializing doc $d")
    Tdoc = SharedArray(zeros(Int32,doclengths[d]))
    Adoc = SharedArray(zeros(Int32,doclengths[d]))
    # @sync begin
    for i in 1:nwork
        # p = workers()[i]
        # initterms!(counts[d,:],Tdoc,Adoc,procinds[i],inds[d,procinds[i]],k)
        initterms!(counts[d,:],Tdoc,Adoc,procinds[i],inds[d,procinds[i]],k)
        # @async remotecall_wait(initterms!,p,counts[d,:],Tdoc,Adoc,procinds[i],inds[d,procinds[i]],k)
    end
    # end
    push!(corp,[Word(Tdoc[i],Adoc[i]) for i in 1:doclengths[d]])
  end
  corp
end

"""
  initcorpBeta!(corp,beta)

Randomly initialize the word assignments in a corpus, based on the topic word prior

# Examples
```
initcorpBeta!(corp,beta)
```
"""
function initcorpBeta(corp,beta)
  newcorp = deepcopy(corp)
  for i in 1:length(corp)
    for j in 1:length(corp[i])
      atype = typeof(corp[i][j].a)
      betanorm = beta[corp[i][j].t,:] ./ sum(beta[corp[i][j].t,:])
      global newcorp[i][j] = Word(corp[i][j].t,atype(rand(Categorical(betanorm))))
    end
  end
  newcorp
end

"""
  initcorpChain!(corp,chain_i)

Initialize the word assignments in a corpus, based on the chain id

# Examples
```
initcorpChain!(corp,chain_i)
```
"""
function initcorpChain(corp,chain_i,k)
  newcorp = deepcopy(corp)
  for i in 1:length(corp)
    for j in 1:length(corp[i])
      atype = typeof(corp[i][j].a)
      global newcorp[i][j] = Word(corp[i][j].t,atype(mod(chain_i,k)+1))
    end
  end
  newcorp
end

"""
  initcorpChain!(corp,chain_i)

Uniform random initialization

# Examples
```
initcorpUnif!(corp,chain_i)
```
"""
function initcorpUnif(corp,k)
  newcorp = deepcopy(corp)
  for i in 1:length(corp)
    for j in 1:length(corp[i])
      atype = typeof(corp[i][j].a)
      global newcorp[i][j] = Word(corp[i][j].t,atype(rand([1:k;])))
    end
  end
  newcorp
end

"""
  getSynthData(nbulk,ncells,nreads,zeta,K)

Create `nbulk` pseudobulk samples from `ncells` simulated cells.
Generate n x `ncells` reads where n ~ Poisson(λ = `nreads`).
Allocate these bulk reads based to a queue based on fidelity matrix `zeta`.

# Examples
```
getSynthData(nbulk,ncells,nreads,zeta,K)
```
"""
function getSynthData(nbulk,ncells,nreads,zeta,K)
  println("done")
end


text:
module TestLogging

import Historic
using Logging
using Test

using ..Samples: Samples, Default

function with_logging(f, recordmodule, logger = true)
    recordmodule.set_logger(logger)
    try
        Base.invokelatest(f)
    finally
        recordmodule.disable_logging()
    end
end

function test_default()
    Historic.clear(Default)
    main = Test.TestLogger()
    with_logging(Default, main) do
        Samples.default_simple()
        Samples.default_noyield()
        Samples.default_withdata(1)
        Samples.default_withdata(2)
        Samples.default_throw()
    end
    prefix = Historic.defaultprefix()
    @test main.logs[1].message == "$(prefix)simple"
    @test main.logs[2].message == "$(prefix)withdata"
    @test main.logs[3].message == "$(prefix)withdata"
    @test main.logs[3].message == "$(prefix)withdata"
    return main
end

end  # module


text:
# DCT operators: joDCT

## helper module
module joDCT_etc
    using JOLI: jo_convert
    using FFTW
    ### planned
    function apply_dct(pln::FFTW.DCTPlan,v::Vector{vdt},ms::Dims,rdt::DataType) where vdt<:Union{AbstractFloat,Complex}
        rv=reshape(v,ms)
        rv=pln*rv
        rv=vec(rv)
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    function apply_dct(pln::FFTW.DCTPlan,v::Matrix{vdt},ms::Dims,rdt::DataType) where vdt<:Union{AbstractFloat,Complex}
        lms=length(ms)
        nvc=size(v,2)
        msc=(ms...,nvc)
        rv=reshape(v,msc)
        pf=plan_dct(rv,[1:lms...])
        rv=pf*rv
        rv=reshape(rv,(prod(ms),nvc))
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    function apply_idct(pln::FFTW.DCTPlan,v::Vector{vdt},ms::Dims,rdt::DataType) where vdt<:Union{AbstractFloat,Complex}
        rv=reshape(v,ms)
        rv=pln*rv
        rv=vec(rv)
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    function apply_idct(pln::FFTW.DCTPlan,v::Matrix{vdt},ms::Dims,rdt::DataType) where vdt<:Union{AbstractFloat,Complex}
        lms=length(ms)
        nvc=size(v,2)
        msc=(ms...,nvc)
        rv=reshape(v,msc)
        pf=plan_idct(rv,[1:lms...])
        rv=pf*rv
        rv=reshape(rv,(prod(ms),nvc))
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    ### not planned
    function apply_dct(v::Vector{vdt},ms::Dims,rdt::DataType) where vdt<:Union{AbstractFloat,Complex}
        rv=reshape(v,ms)
        rv=dct(rv)
        rv=vec(rv)
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    function apply_dct(v::Matrix{vdt},ms::Dims,rdt::DataType) where vdt<:Union{AbstractFloat,Complex}
        lms=length(ms)
        nvc=size(v,2)
        msc=(ms...,nvc)
        rv=reshape(v,msc)
        pf=plan_dct(rv,[1:lms...])
        rv=pf*rv
        rv=reshape(rv,(prod(ms),nvc))
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    function apply_idct(v::Vector{vdt},ms::Dims,rdt::DataType) where vdt<:Union{AbstractFloat,Complex}
        rv=reshape(v,ms)
        rv=idct(rv)
        rv=vec(rv)
        rv=jo_convert(rdt,rv,false)
        return rv
    end
    function apply_idct(v::Matrix{vdt},ms::Dims,rdt::DataType) where vdt<:Union{AbstractFloat,Complex}
        lms=length(ms)
        nvc=size(v,2)
        msc=(ms...,nvc)
        rv=reshape(v,msc)
        pf=plan_idct(rv,[1:lms...])
        rv=pf*rv
        rv=reshape(rv,(prod(ms),nvc))
        rv=jo_convert(rdt,rv,false)
        return rv
    end
end
using .joDCT_etc

export joDCT
"""
    julia> op = joDCT(m[,n[,...]];[planned::Bool=true][,DDT=joFloat][,RDT=DDT])

Multi-dimensional DCT transform over fast dimension(s)

# Signature

    function joDCT(ms::Integer...;
        planned::Bool=true,
        DDT::DataType=joFloat,RDT::DataType=DDT,
        name::String="joDCT")

# Arguments

- `m`: dimension
- optional
    - `n,...`: more dimensions
- keywords
    - `planned`: use DCT plan
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Notes

- multidimensional image must be vectorized
- if you intend to use joDCT in remote* calls, you have to either set planned=false or create the operator on the worker
- joDCT is always planned if applied to multi-vector

# Examples

1D DCT

    joDCT(m)

1D DCT without the precomputed plan

    joDCT(m; planned=false)

2D DCT

    joDCT(m,n)

examples with DDT/RDT

    joDCT(m; DDT=Float32)
    joDCT(m; DDT=Float32,RDT=Float64)

"""
function joDCT(ms::Integer...;
    planned::Bool=true,
    DDT::DataType=joFloat,RDT::DataType=DDT,
    name::String="joDCT")

    if planned
        pf=plan_dct(zeros(ms))
        ipf=plan_idct(zeros(ms))
        return joLinearFunction_A(prod(ms),prod(ms),
            v1->joDCT_etc.apply_dct(pf,v1,ms,RDT),
            v2->joDCT_etc.apply_idct(ipf,v2,ms,DDT),
            v3->joDCT_etc.apply_idct(ipf,v3,ms,DDT),
            v4->joDCT_etc.apply_dct(pf,v4,ms,RDT),
            DDT,RDT;fMVok=true,iMVok=true,
            name=name*"_p"
            )
    else
        return joLinearFunction_A(prod(ms),prod(ms),
            v1->joDCT_etc.apply_dct(v1,ms,RDT),
            v2->joDCT_etc.apply_idct(v2,ms,DDT),
            v3->joDCT_etc.apply_idct(v3,ms,DDT),
            v4->joDCT_etc.apply_dct(v4,ms,RDT),
            DDT,RDT;fMVok=true,iMVok=true,
            name=name)
    end
end


text:
# Build the product of an MDP/POMDP and a deterministic finite state automata

struct ProductState{S, Q}
    s::S
    q::Q
end

function Base.:(==)(a::ProductState{S, Q}, b::ProductState{S, Q}) where {S, Q}
    return a.s == b.s && a.q == b.q
end

function Base.hash(s::ProductState{S, Q}, h::UInt64) where {S, Q}
    return hash(s.s, hash(s.q, h))
end

mutable struct ProductMDP{S, A, Q <: Int, R <: AbstractAutomata, F<:AbstractFloat} <: MDP{ProductState{S, Q}, A}
    problem::MDP{S, A}
    automata::R
    accepting_states::Set{ProductState{S, Q}}
    sink_state::ProductState{S, Q}
    discount::F
end

function ProductMDP(mdp::MDP{S}, automata::AbstractAutomata, sink_state::ProductState{S,Q}, discount::AbstractFloat) where {S,Q}
    ProductMDP(mdp, automata, Set{ProductState{S, Q}}(), sink_state, discount)
end

mutable struct ProductPOMDP{S, A, O, Q <: Int, R <: AbstractAutomata, F<:AbstractFloat} <: POMDP{ProductState{S, Q}, A, O}
    problem::POMDP{S, A, O}
    automata::R
    accepting_states::Set{ProductState{S, Q}}
    sink_state::ProductState{S, Q}
    discount::F
end

function ProductPOMDP(pomdp::POMDP{S}, automata::AbstractAutomata, sink_state::ProductState{S,Q}, discount::AbstractFloat) where {S,Q}
    ProductPOMDP(pomdp, automata, Set{ProductState{S, Q}}(), sink_state, discount)
end

# should be implemented by the problem writer
"""
Returns the labels associated with state s
For each state, it should return a list of atomic proposition that evaluate to true, all the other propositions are assumed false.
labels(mdp::M, s, a) where {M <: Union{MDP,POMDP}}
"""
function labels end

# returns the set of accepting states
function accepting_states!(mdp::M; verbose::Bool=false) where {M <: Union{ProductMDP, ProductPOMDP}}
    MECs = maximal_end_components(mdp, verbose=verbose)
    state_space = states(mdp)
    mdp.accepting_states = Set{statetype(mdp)}()
    fin_inf_sets = get_rabin_acceptance(mdp.automata)
    verbose ? println("Extracting accepting states from MECs ... \n") : nothing
    for ec in MECs
        ec_states = Set(state_space[i].q for i in ec)
        for (fin_q, inf_q) in fin_inf_sets
            if !isempty(intersect(ec_states, inf_q)) && isempty(intersect(ec_states, fin_q))
                for si in ec
                    push!(mdp.accepting_states, state_space[si])
                end
            end
        end
    end
    verbose ? println("Accepting states computed. \n") : nothing
    return mdp.accepting_states
end

# # inherit functions from the MDP api

function POMDPs.reward(mdp::Union{ProductMDP, ProductPOMDP}, s::ProductState{S, Q}, a::A, sp::ProductState{S, Q}) where {S, A, Q}
    if sp ∈ mdp.accepting_states
        return 1.
    end
    return 0.0
end

# POMDPs.reward(mdp::Union{ProductMDP, ProductPOMDP}, s::ProductState{S, Q}, a::A) where {S, A, Q} = reward(mdp, s, a, s)

function POMDPs.isterminal(mdp::Union{ProductMDP, ProductPOMDP}, s::ProductState{S, Q}) where {S, Q}
    if s ∈ mdp.accepting_states || s == mdp.sink_state
        return true
    end
    return false
end

POMDPs.discount(problem::Union{ProductMDP, ProductPOMDP}) = problem.discount

# in the product MDP, some transitions are "undefined" because the automata does not allow them.
# the transitions does not necessarily sums up to one!
function POMDPs.transition(problem::Union{ProductMDP, ProductPOMDP}, state::ProductState{S, Q}, action::A) where {S, A, Q}
    d = transition(problem.problem, state.s, action) # regular mdp transition
    if state ∈ problem.accepting_states || state == problem.sink_state
        return SparseCat{Vector{statetype(problem)}, Vector{Float64}}([problem.sink_state], [1.0])
    end
    l = labels(problem.problem, state.s, action)
    qp = nextstate(problem.automata, state.q, l)
    if qp != nothing
        new_vals = Vector{ProductState{S,Q}}(undef, length(support(d)))
        new_probs = Vector{Float64}(undef, length(support(d)))
        for (i, sp) in enumerate(support(d))
            new_vals[i] = ProductState(sp,qp)
            new_probs[i] = pdf(d, sp)
        end
        return SparseCat{Vector{ProductState{S,Q}}, Vector{Float64}}(new_vals, new_probs)
    else
        return SparseCat{Vector{statetype(problem)}, Vector{Float64}}([problem.sink_state], [1.0])
    end
end

function POMDPs.initialstate(problem::Union{ProductMDP, ProductPOMDP})
    b0 = initialstate(problem.problem)
    new_probs = Float64[]
    new_vals = Vector{statetype(problem)}()
    q0 = get_init_state_number(problem.automata)
    for (s0, p) in weighted_iterator(b0)
        push!(new_vals, ProductState(s0, q0))
        push!(new_probs, p)
    end
    normalize!(new_probs, 1)
    return SparseCat{Vector{statetype(problem)}, Vector{Float64}}(new_vals, new_probs)
end

function POMDPs.states(problem::Union{ProductMDP, ProductPOMDP})
    S = statetype(problem.problem)
    Q = eltype(problem.automata.states)
    state_space = ProductState{S, Q}[]
    for s in ordered_states(problem.problem)
        for q in 1:num_states(problem.automata)
            push!(state_space, ProductState(s, q))
        end
    end
    push!(state_space, problem.sink_state)
    return state_space
end

POMDPs.actions(problem::Union{ProductMDP, ProductPOMDP}) = actions(problem.problem)

function POMDPs.stateindex(problem::Union{ProductMDP, ProductPOMDP}, s::ProductState{S,Q}) where {S,Q}
    if s == problem.sink_state
        return length(states(problem))
    end
    si = stateindex(problem.problem, s.s)
    qi = s.q
    return LinearIndices((num_states(problem.automata), length(states(problem.problem))))[qi, si]
end

POMDPs.statetype(p::Union{ProductMDP, ProductPOMDP}) = ProductState{statetype(p.problem), typeof(p.sink_state.q)}

POMDPs.actiontype(p::Union{ProductMDP, ProductPOMDP}) = actiontype(p.problem)

POMDPs.actionindex(p::Union{ProductMDP, ProductPOMDP}, a) = actionindex(p.problem, a)

POMDPs.convert_a(T::Type{V}, a, p::Union{ProductMDP, ProductPOMDP}) where V<:AbstractArray = convert_a(T, a, p.problem)
POMDPs.convert_a(T::Type{A}, vec::V, p::Union{ProductMDP, ProductPOMDP}) where {A,V<:AbstractArray} = convert_a(T, vec, p.problem)

function POMDPs.convert_s(T::Type{Vector{Float64}}, s::ProductState{S,Int64}, p::Union{ProductMDP, ProductPOMDP}) where S
    v_mdp = convert_s(T, s.s, p.problem) # convert mdp state
    v_autom = zeros(num_states(p.automata))
    v_autom = 1.0
    return cat(1, v_mdp, v_autom)
end

function POMDPs.convert_s(::Type{ProductState{S,Int64}}, vec::Vector{Float64}, p::Union{ProductMDP, ProductPOMDP}) where {S}
    v_mdp = vec[1:end-num_states(p.automata)]
    v_autom = vec[end-num_states(p.automata)+1:end]
    s = convert_s(S, v_mdp, p.problem)
    q = findfirst(v_autom)
    return ProductState(s, q)
end

# ## POMDP Only

POMDPs.observation(p::ProductPOMDP, s::ProductState{S, Q}) where {S,Q} = observation(p.problem, s.s)
POMDPs.observation(p::ProductPOMDP, a::A, s::ProductState{S, Q}) where {S,Q,A} = observation(p.problem, a, s.s)
POMDPs.observation(p::ProductPOMDP, s::ProductState{S, Q}, a::A, sp::ProductState{S, Q}) where {S,Q,A}= observation(p.problem, s.s, a, sp.s)
POMDPs.observations(p::ProductPOMDP) = observations(p.problem)
POMDPs.obsindex(p::ProductPOMDP, o::O) where O = obsindex(p.problem, o)
POMDPs.convert_o(T::Type{V}, o, p::ProductPOMDP) where V<:AbstractArray = convert_o(T, o, p.problem)
POMDPs.convert_o(T::Type{O}, vec::V, p::ProductPOMDP) where {O,V<:AbstractArray} = convert_o(T, vec, p.problem)


text:
module TidalFluxStats

using TidalFluxQuantities

export tidalbalance,
    tidalaverage,
    segment,
    superimpose,
    adjustbias,
    adjustbias_average,
    adjustbias_discharge,
    prism,
    adjustbias_prism,
    Despike

include("averages.jl")
include("superimpose.jl")
include("bias.jl")
include("prism.jl")
include("despike.jl")
include("masks.jl")

end # module


text:
using OptimKit
using ..ITensorNetworks
using ..ITensorNetworks:
  PEPS, broadcast_add, broadcast_minus, broadcast_mul, broadcast_inner

function OptimKit.optimize(peps::PEPS, loss_w_grad; num_sweeps::Int, method="GD")
  @assert(method in ["GD", "LBFGS", "CG"])
  inner(x, peps1, peps2) = broadcast_inner(peps1, peps2)
  scale(peps, alpha) = broadcast_mul(alpha, peps)
  add(peps1, peps2, alpha) = broadcast_add(peps1, broadcast_mul(alpha, peps2))
  retract(peps1, peps2, alpha) = (add(peps1, peps2, alpha), peps2)
  linesearch = HagerZhangLineSearch(; c₁=0.1, c₂=0.9, verbosity=0)
  if method == "GD"
    alg = GradientDescent(num_sweeps, 1e-8, linesearch, 2)
  elseif method == "LBFGS"
    alg = LBFGS(16; maxiter=num_sweeps, gradtol=1e-8, linesearch=linesearch, verbosity=2)
  elseif method == "CG"
    alg = ConjugateGradient(;
      maxiter=num_sweeps, gradtol=1e-8, linesearch=linesearch, verbosity=2
    )
  end
  _, _, _, _, history = OptimKit.optimize(
    loss_w_grad, peps, alg; inner=inner, (scale!)=scale, (add!)=add, retract=retract
  )
  return history[:, 1]
end

function gd_error_tracker(
  peps,
  loss_w_grad,
  loss_w_grad_approx;
  stepsize::Float64,
  num_sweeps::Int,
  cutoff=1e-15,
  maxdim=100,
)
  for iter in 1:num_sweeps
    l, g = loss_w_grad(peps)
    l_approx, g_approx = loss_w_grad_approx(peps)
    g_diff = broadcast_minus(g, g_approx)
    g_diff_nrm = broadcast_inner(g_diff, g_diff)
    print("The gradient difference norm at iteraton $iter is $g_diff_nrm\n")
    print("The rayleigh quotient at iteraton $iter is $l\n")
    print("The approximate rayleigh quotient at iteraton $iter is $l_approx\n")
    peps = broadcast_minus(peps, broadcast_mul(stepsize, g))
  end
end

function backtracking_linesearch(beta::Float64, loss_w_grad, peps)
  stepsize = 1.0
  stepsize_lb = 0.0
  l, g = loss_w_grad(peps)
  update_peps = broadcast_minus(peps, broadcast_mul(stepsize, g))
  update_l, update_g = loss_w_grad(update_peps)
  gnrm_square = broadcast_inner(g, g)
  threshold = l - stepsize / 2.0 * broadcast_inner(g, g)
  print("threshold: $threshold, gradient norm: $gnrm_square \n")
  while update_l > threshold && stepsize >= stepsize_lb
    stepsize = stepsize * beta
    print("stepsize trial: $stepsize \n")
    update_peps = broadcast_minus(peps, broadcast_mul(stepsize, g))
    update_l, update_g = loss_w_grad(update_peps)
    threshold = l - stepsize / 2.0 * broadcast_inner(g, g)
    diffg = broadcast_minus(g, update_g)
    update_gnrm_square = broadcast_inner(diffg, diffg)
    print("threshold: $threshold, update_l: $update_l, update_g: $update_gnrm_square \n")
  end
  print("stepsize from backtracking_linesearch is $stepsize \n")
  return stepsize
end

function gradient_descent(peps::PEPS, loss_w_grad; stepsize::Float64, num_sweeps::Int)
  # gradient descent iterations
  losses = []
  for iter in 1:num_sweeps
    l, g = loss_w_grad(peps)
    print("The rayleigh quotient at iteraton $iter is $l\n")
    peps = broadcast_minus(peps, broadcast_mul(stepsize, g))
    push!(losses, l)
  end
  return losses
end

function gradient_descent(
  peps::PEPS, loss_w_grad, ::typeof(backtracking_linesearch); beta::Float64, num_sweeps::Int
)
  # gradient descent iterations
  losses = []
  for iter in 1:num_sweeps
    l, g = loss_w_grad(peps)
    print("The rayleigh quotient at iteraton $iter is $l\n")
    stepsize = backtracking_linesearch(beta, loss_w_grad, peps)
    peps = broadcast_minus(peps, broadcast_mul(stepsize, g))
    push!(losses, l)
  end
  return losses
end


text:
using StochasticDiffEq, Test, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_stiffquadito

Random.seed!(100)
prob = prob_sde_stiffquadito
prob = remake(prob;p=(1e5,2.))
alg = AutoSOSRA2(SKenCarp(), maxstiffstep=5, maxnonstiffstep=2, stiffalgfirst=true)
@test StochasticDiffEq.isadaptive(alg)
@time sol = solve(prob, alg);
@test typeof(alg.algs[sol.alg_choice[1]]) <: SKenCarp
@test length(unique(sol.alg_choice)) == 2

Random.seed!(100)
prob = prob_sde_stiffquadito
prob = remake(prob;p=(1e5,2.))
@time sol = solve(prob, AutoSOSRI2(ImplicitRKMil(), maxstiffstep=2, maxnonstiffstep=2, stiffalgfirst=true),
                        abstol=1e-4,reltol=1e-4);
@test length(unique(sol.alg_choice)) == 2


text:
# %% imports

using PyCall
gym = pyimport("gym")

include("../../Modules/Gridcells/StateSpikeGaussGrid.jl")
include("../../Modules/Networks/net_arch.jl")
include("../../Modules/Gridcells/GridGenerator.jl")

# %% sim functions

function count_spikes(arch::Array{Int64}, spiked::Array{Int64}, layer::Int64)
    count_layer = zeros(arch[layer])
    interval = [sum(arch[1:layer-1])+1, sum(arch[1:layer])]
    spikes_in = filter(x-> interval[1] <= x <= interval[2], spiked)
    spikes_in .-= sum(arch[1:layer-1])
    count_layer[spikes_in] .+= 1
    return count_layer
end

function left_right(count_spikes::Array{Float64}, n_left::Int64, n_right::Int64)
    sl = sum(count_spikes[1:n_left])
    sr = sum(count_spikes[n_left+1:end])
    (sl > sr) ? action = 0 : action = 1
    return action, sl, sr
end

# %% fitness function play 100 episodes and check les contraintes a chaque épisode

function play_episode(net::Network, n_steps::Int64 = 200)
    env = gym.make("CartPole-v1")
    obs = env.reset()
    win_size = 10 # ms
    grid_cells = grid_gen()
    spike_hist = Array{Int64}[]
    score_hist = [[],[]]
    decision  = 0
    rng_action = Random.MersenneTwister(0)
    @inbounds for j in 1:n_steps
        count_tot = zeros(net.arch[end])
        spikes = input_spikes(obs, grid_cells)
        influence = 2
        @inbounds for msec in 1:win_size
            (msec==10) ? train=true : train=false
            spiked = step!(net, spikes, train)
            if msec >= length(net.arch)-1
                count_tot .+= influence .* count_spikes(net.arch, spiked, length(net.arch))
                influence *= 0.95
                push!(spike_hist,spiked)
            end
        end
        n_left = div(net.arch[end],2)
        n_right = net.arch[end] - n_left
        action, sl, sr = left_right(count_tot, n_left, n_right)
        if sl == sr
            action = rand(rng_action, 0:1)
            decision += 0.001
        else
            decision += 0.01
        end
        push!(score_hist[1], sl)
        push!(score_hist[2], sr)
        obs_new, _, _, _ = env.step(action)
        obs = obs_new
    end
    env.close()
    env = nothing
    Base.GC.gc()
    return score_hist, spike_hist, decision
end

# CONSTRAINTS
function constraints(score::Array{Array{Any,1},1}, spikes::Array{Array{Int64}}, net::Network)
    mw_min = 0.0 #FILL
    mw_max = 1.0#FILL
    s_min = 1000#FILL
    s_max = 15000#FILL
    m_min = -5.0#FILL
    m_max = 5.0 #FILL
    # constraint #1 weight average
    weights = net.s.e
    for layer in 1:length(net.arch)-1
        mw = mean(weights[layer])
        if (mw < mw_min) | (mw > mw_max)
            println("-----> weights mean out of bounds")
            return true
        end
    end
    # constraint #2 layer activity
    for layer in 1:length(net.arch)
        count_layer = zeros(net.arch[layer])
        for i in 1:length(spikes)
            count_layer .+= count_spikes(net.arch, spikes[i], layer)
        end
        spikes_layer = sum(count_layer)
        if (spikes_layer < s_min) | (spikes_layer > s_max)
            println("-----> spike activity out of bounds")
            return true
        end
    end
    # constraint #3 mean left right
    mlr = mean(score[1] .- score[2])
    if (mlr > m_max) | (mlr < m_min)
        println("-----> left right mean out of bounds")
        return true
    end
    return false
end

# plays 100 episodes and returns the fitness ( and tells wether to drop batch )
function fitness(episodes::Int64, net::Network)
    #score_hist = [[0.,0.],[0.,0.]]
    drop = false
    fit = 0
    for episode in 1:episodes
        score, spikes, decision = play_episode(net)
        drop = constraints(score, spikes, net)
        #push!(score_hist[1],score[1]...)
        #push!(score_hist[2],score[2]...)
        fit += decision
        if drop
            break
        end

    end
    #fit = std(score_hist[1] .- score_hist[2]) + pond
    return fit
end

# %% objective function ( fitness + contraintes ) ( create new network and play 100 episodes with fitness function )

function rescale(genes::Array{Float64}, tunekeys::Array{String})
    cfg = YAML.load(open("/Users/titou/Documents/PFE/RL-STDP/Julia/CartPole/CMAES/cfg_cmaes.yml"))
    N = length(genes)
    scaled_genes = zeros(N)
    for idx in 1:N
        down = cfg[tunekeys[idx]*"_d"]
        up = cfg[tunekeys[idx]*"_u"]
        scaled_genes[idx] = ((genes[idx]+1)*0.5)*(up-down) + down
    end
    return scaled_genes
end

function descale(params::Array{Float64}, tunekeys::Array{String})
    cfg = YAML.load(open("/Users/titou/Documents/PFE/RL-STDP/Julia/CartPole/CMAES/cfg_cmaes.yml"))
    N = length(tunekeys)
    descaled_genes = zeros(N)
    for idx in 1:N
        down = cfg[tunekeys[idx]*"_d"]
        up = cfg[tunekeys[idx]*"_u"]
        descaled_genes[idx] = 2*(params[idx]-down)/(up-down) - 1
    end
    return descaled_genes
end

function generate(genes::Array{Float64})
    tunekeys = ["vthresh", "ae", "b", "c", "de","imin", "imax",
                "m", "std", "apost", "apre"]
    params = Dict(tunekeys .=> rescale(genes,tunekeys))
    keys = ["tda", "wmax" ,"theta", "ttheta", "ai", "di", "wei", "wie"]
    values = [0.995, 1, 0.0, 0.998, 0.1, 2.0, 0.0, 0.0]
    for idx in 1:length(keys)
        params[keys[idx]] = values[idx]
    end
    return params
end

function objective(genes::Array{Float64})
    params = generate(genes)
    n = Network([120,60,80],params)
    fit = fitness(50,n)
    @show fit
    -fit
end

# %% tests
#
# objective([-50.0,0.02,0.2,-65.0,8.0,0.5,0.1,0.0,0.998,1.0,1.5])

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
# ps3000_get_streaming_values.jl

export ps3000_get_streaming_values!, StreamingBufferAggregation

immutable StreamingBuffer
  a_max :: Array{Int16,1}
  a_mim :: Array{Int16,1}
  b_max :: Array{Int16,1}
  b_mim :: Array{Int16,1}
  c_max :: Array{Int16,1}
  c_mim :: Array{Int16,1}
  d_max :: Array{Int16,1}
  d_mim :: Array{Int16,1}
  function StreamingBuffer(length)
    new(Array(Int16,length),Array(Int16,length),Array(Int16,length),Array(Int16,length),
        Array(Int16,length),Array(Int16,length),Array(Int16,length),Array(Int16,length))
  end
end

function ps3000_get_streaming_values!(handle       :: Int16,
                                      st           :: Float64,
                                      buffer       :: StreamingBuffer,
                                      no_of_values :: UInt32,
                                      no_of_samples_per_aggreate :: UInt32)
  @assert length(buffer.buffer_a_max)>=no_of_values  # make sure buffer is big enough
  start_time = Ref{Float64}(st)   # start, relative to trigger, time in ns
  overflow = Ref{Int16}           # bitfield, see doc
  triggered_at = Ref{UInt32}      # index trigger occured at
  triggered = Ref{UInt32}         # was there a valid trigger?
  number_written_to_buffer = ccall((:ps3000_get_streaming_values, ps3000driver),UInt32,
    (Int16, Ref{Cdouble},
    Ptr{Int16}, Ptr{Int16},
    Ptr{Int16}, Ptr{Int16},
    Ptr{Int16}, Ptr{Int16},
    Ptr{Int16}, Ptr{Int16},
    Ref{Int16}, Ref{UInt32}, Ref{Int16}, UInt32, UInt32),
    handle, start_time,
    buffer.a_max, buffer.a_min,
    buffer.b_max, buffer.b_min,
    buffer.c_max, buffer.c_min,
    buffer.d_max, buffer.d_min,
    overflow, triggered_at, triggered, no_of_values, no_of_samples_aggreate)
  if number_written_to_buffer == 0
    error("ps3000 Error: Parameter Out Of Range")
  end
  return (number_written_to_buffer,overflow[],triggered_at[],triggered[]==0?false:true)
end


text:
# /* ----------------------------------------------------------------------
#    compute dEidRj
# ------------------------------------------------------------------------- */

function compute_deidrj(dedr::Vector{AbstractFloat}, jelem::Int)

    for k = 1:3
        dedr[k] = 0.0;
    end

    for j = 0:snap.twojmax #(int j = 0; j <= twojmax; j++) {
        jju = idxu_block[j+1];

        mb = 0
        ## Note, unlike many other while loops, this is strictly less than (see original code)
        while 2*mb < j
            for ma = 0:j

                dudr_r = dulist_r[jju+1];
                dudr_i = dulist_i[jju+1];
                jjjmambyarray_r = ylist_r[(jelem-1)*idxu_max+jju+1];
                jjjmambyarray_i = ylist_i[(jelem-1)*idxu_max+jju+1];

                for k = 1:3
                    dedr[k] +=
                    dudr_r[k] * jjjmambyarray_r +
                    dudr_i[k] * jjjmambyarray_i;
                end
                jju+=1;
            end
            mb += 1
        end # loop over ma mb

        # // For j even, handle middle column

        if (j%2 == 0)
            mb = floor(Int, j/2);
            ma = 0
            # Also strictly <
            while ma < mb
                dudr_r = dulist_r[jju+1, :];
                dudr_i = dulist_i[jju+1, :];
                jjjmambyarray_r = ylist_r[(jelem-1)*idxu_max+jju+1];
                jjjmambyarray_i = ylist_i[(jelem-1)*idxu_max+jju+1];

                for k=1:3
                    dedr[k] +=
                        dudr_r[k] * jjjmambyarray_r +
                        dudr_i[k] * jjjmambyarray_i;

                end
                jju+=1;
                ma+=1;
            end

            dudr_r = dulist_r[jju+1];
            dudr_i = dulist_i[jju+1];
            jjjmambyarray_r = ylist_r[(jelem-1)*idxu_max+jju+1];
            jjjmambyarray_i = ylist_i[(jelem-1)*idxu_max+jju+1];

            for k = 1:3
                dedr[k] +=
                    (dudr_r[k] * jjjmambyarray_r +
                    dudr_i[k] * jjjmambyarray_i)*0.5;
                # // jju++;
            end
        end # if j even

    end # loop over j

    for k = 1:3
        dedr[k] *= 2.0;
    end
end


text:
# This file illustrates calculation of energy spectral density and total enery of a rectangular pulse.

using FFTW
using Plots; plotlyjs()

# Define pulse
T = 1 / 10
p(t) = 0 ≤ t ≤ T ? 1 : 0

# Construct the signal
fs = 1000
ts = 1 / fs
l = 2^10
t = (0 : l - 1) * ts
x = p.(t)

# Compute fft
xf = fft(x)
psd = abs.(xf).^2
ff = (0 : 1 / (l - 1) : 1) * fs

# Compute energies
Δf = fs / (l - 1)
Δt = ts
Et = sum(abs.(x).^2) * Δt
Ef = sum(psd) * Δf / (fs^2)
@show Et, Ef

# Plots
plt = plot(layout=2)
plot!(t[1:500], x[1:500], subplot=1)
plot!(ff, psd, subplot=2)
display(plt)


text:
using Printf
using Logging
using Plots

using Oceananigans
using Oceananigans.Advection
using Oceananigans.TimeSteppers

using Oceananigans.Grids: xnodes

ENV["GKSwstype"] = "100"

#####
##### Initial conditions and analytic solutions
#####

# Convert x to periodic range in [-L/2, L/2] assuming u = 1.
@inline x′(x, t, L, U) = mod(x + L/2 - U * t, L) - L/2

# Analytic solutions
@inline ϕ_Gaussian(x, t; L, U, a=1, c=1/8) = a * exp(-x′(x, t, L, U)^2 / (2c^2))
@inline ϕ_Square(x, t; L, U, w=0.15) = -w <= x′(x, t, L, U) <= w ? 1.0 : 0.0

ic_name(::typeof(ϕ_Gaussian)) = "Gaussian"
ic_name(::typeof(ϕ_Square))   = "Square"

#####
##### Experiment functions
#####

function setup_model(N, L, U, ϕₐ, time_stepper, advection_scheme)
    topology = (Periodic, Flat, Flat)
    domain = (x=(-L/2, L/2), y=(0, 1), z=(0, 1))
    grid = RegularCartesianGrid(topology=topology, size=(N, 1, 1), halo=(9, 9, 9); domain...)

    model = IncompressibleModel(
               grid = grid,
        timestepper = time_stepper,
          advection = advection_scheme,
            tracers = :c,
           buoyancy = nothing,
            closure = IsotropicDiffusivity(ν=0, κ=0)
    )

    set!(model, u = U, v = (x, y, z) -> ϕₐ(x, 0; L=L, U=U), c = (x, y, z) -> ϕₐ(x, 0; L=L, U=U))

    return model
end

function short_name(ts)
    ts == :QuasiAdamsBashforth2 && return "QAB2"
    ts == :RungeKutta3 && return "RK3"
end

weno_order(::WENO{K}) where K = 2K-1

function create_animation(N, L, CFL, ϕₐ, time_stepper, advection_scheme; U=1.0, T=2.0)
    model = setup_model(N, L, U, ϕₐ, time_stepper, advection_scheme)

    v, c = model.velocities.v, model.tracers.c
    x = xnodes(c)
    Δt = CFL * model.grid.Δx / abs(U)
    Nt = ceil(Int, T/Δt)

    function every(n)
          0 < n <= 128 && return 1
        128 < n <= 256 && return 2
        256 < n <= 512 && return 4
        512 < n        && return 8
    end

    scheme_name = advection_scheme isa WENO ? "WENO$(weno_order(advection_scheme))" : typeof(advection_scheme)

    anim_filename = @sprintf("%s_%s_%s_N%d_CFL%.2f_U%+d.gif", ic_name(ϕₐ), time_stepper, scheme_name, N, CFL, U)

    anim = @animate for iter in 0:Nt
        iter % 10 == 0 && @info "$anim_filename, iter = $iter/$Nt"

        ϕ_analytic = ϕₐ.(x, model.clock.time; L=L, U=U)

        title = @sprintf("%s %s N=%d CFL=%.2f", short_name(time_stepper), scheme_name, N, CFL)
        plot(x, ϕ_analytic, lw=2, label="analytic", title=title, xlims=(-L/2, L/2), ylims=(-0.2, 1.2), dpi=200)
        plot!(x, interior(v)[:], ls=:solid, lw=2, label="Oceananigans v")
        plot!(x, interior(c)[:], ls=:dot,   lw=2, label="Oceananigans c")

        if time_stepper == :QuasiAdamsBashforth2
            time_step!(model, Δt, euler = iter == 0)
        else
            time_step!(model, Δt)
        end
    end every every(Nt)

    gif(anim, anim_filename, fps=15)

    return model
end

#####
##### Run 1D periodic advection experiments
#####

L = 1
ϕs = (ϕ_Gaussian, ϕ_Square)
time_steppers = (:RungeKutta3,)
advection_schemes = (CenteredSecondOrder(), CenteredFourthOrder(), UpwindBiasedThirdOrder(), UpwindBiasedFifthOrder(), WENO5())
Ns = [16, 64]
CFLs = (0.5, 1.7)
Us = [+1, -1]

for ϕ in ϕs, ts in time_steppers, scheme in advection_schemes, N in Ns, CFL in CFLs, U in Us
    scheme_name = scheme isa WENO ? "WENO$(weno_order(scheme))" : typeof(scheme)
    @info @sprintf("Creating two-revolution animation [%s, %s, %s, N=%d, CFL=%.2f, U=%+d]...", ic_name(ϕ), ts, scheme_name, N, CFL, U)
    create_animation(N, L, CFL, ϕ, ts, scheme, U=U)
end


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

struct PowerBasisCurve{P,S,N,M} <: Spline{P,S,N,M}
    controlpolygon::Array{S,2} #not exactly the right name, should be coefficients instead, but consistent with the other splines
    #data stored like this:
    # a0x a1x ... aM+1x
    # a0y a1y ... aM+1x
    # and so on for however many spatial dimensions

    function PowerBasisCurve{P,S,N,M}(coefficients::Array{S,2}) where {P,S,N,M}
        @assert (N, M + 1) == size(coefficients)

        return new{P,S,N,M}(copy(coefficients))
    end
end
export PowerBasisCurve

function coefficients(curve::PowerBasisCurve{P,S,N,M}, spatialindex) where {P,S,N,M}
    return curve.controlpolygon[spatialindex, :]
end

function PowerBasisCurve{P,S,N,M}(curve::BezierCurve{P,S,N,M}) where {P,S,N,M}
    controlpolygon = curve.controlpolygon
    coefficients = Array{S,2}(undef, N, M + 1) # array stores coefficients for all N spatial dimensions

    #this is probably not the most efficient way to convert from Bernstein to Power basis. Optimize later if necessary.
    for k in 0:M
        for i in k:M
            @. coefficients[:, i + 1] += (-1)^(i - k) * binomial(M, i) * binomial(i, k) * controlpolygon[k + 1]
        end
    end
    return PowerBasisCurve{P,S,N,M}(coefficients)
end

function beziertopowerbasis(k, n)
    coefficient = 0
    for i in k:n
        @. coefficients[:, i + 1] += (-1)^(i - k) * binomial(n, i) * binomial(i, k)
    end
end

function point(curve::PowerBasisCurve, u::T) where {T<:Number}
    _, n = size(curve.controlpolygon)
    sum = curve.controlpolygon[:, n]

    for i in (n - 1):-1:1
        @. sum = u * sum + curve.controlpolygon[:, i]
    end

    return sum
end


text:
# Unit tests for api/exceptions.jl

@testset "api/exceptions.jl" begin
    target = Atomistic.UnimplementedError(:func, 42)

    @test target isa Exception
    @test target isa Atomistic.UnimplementedError{Int64}

    message = sprint(showerror, target)
    @test contains(message, "func")
    @test contains(message, "Int64")
end