text:
update_g!(d, state, method) = nothing
function update_g!(d, state, method::M) where M<:Union{FirstOrderOptimizer, Newton}
    # Update the function value and gradient
    value_gradient!(d, state.x)
    if M <: FirstOrderOptimizer #only for methods that support manifold optimization
        project_tangent!(method.manifold, gradient(d), state.x)
    end
end
update_fg!(d, state, method) = nothing
update_fg!(d, state, method::ZerothOrderOptimizer) = value!(d, state.x)
function update_fg!(d, state, method::M) where M<:Union{FirstOrderOptimizer, Newton}
    value_gradient!(d, state.x)
    if M <: FirstOrderOptimizer #only for methods that support manifold optimization
        project_tangent!(method.manifold, gradient(d), state.x)
    end
end

# Update the Hessian
update_h!(d, state, method) = nothing
update_h!(d, state, method::SecondOrderOptimizer) = hessian!(d, state.x)

after_while!(d, state, method, options) = nothing

function initial_convergence(d, state, method::AbstractOptimizer, initial_x, options)
    gradient!(d, initial_x)
    norm(gradient(d), Inf) < options.g_abstol
end
initial_convergence(d, state, method::ZerothOrderOptimizer, initial_x, options) = false

function optimize(d::D, initial_x::Tx, method::M,
                  options::Options{T, TCallback} = Options(;default_options(method)...),
                  state = initial_state(method, options, d, initial_x)) where {D<:AbstractObjective, M<:AbstractOptimizer, Tx <: AbstractArray, T, TCallback}
    if length(initial_x) == 1 && typeof(method) <: NelderMead
        error("You cannot use NelderMead for univariate problems. Alternatively, use either interval bound univariate optimization, or another method such as BFGS or Newton.")
    end

    t0 = time() # Initial time stamp used to control early stopping by options.time_limit
    tr = OptimizationTrace{typeof(value(d)), typeof(method)}()
    tracing = options.store_trace || options.show_trace || options.extended_trace || options.callback != nothing
    stopped, stopped_by_callback, stopped_by_time_limit = false, false, false
    f_limit_reached, g_limit_reached, h_limit_reached = false, false, false
    x_converged, f_converged, f_increased, counter_f_tol = false, false, false, 0

    g_converged = initial_convergence(d, state, method, initial_x, options)
    converged = g_converged
    # prepare iteration counter (used to make "initial state" trace entry)
    iteration = 0

    options.show_trace && print_header(method)
    _time = time()
    trace!(tr, d, state, iteration, method, options, _time-t0)
    ls_success::Bool = true
    while !converged && !stopped && iteration < options.iterations
        iteration += 1

        ls_success = !update_state!(d, state, method)
        if !ls_success
            break # it returns true if it's forced by something in update! to stop (eg dx_dg == 0.0 in BFGS, or linesearch errors)
        end
        update_g!(d, state, method) # TODO: Should this be `update_fg!`?

        x_converged, f_converged,
        g_converged, f_increased = assess_convergence(state, d, options)
        # For some problems it may be useful to require `f_converged` to be hit multiple times
        # TODO: Do the same for x_tol?
        counter_f_tol = f_converged ? counter_f_tol+1 : 0
        converged = x_converged || g_converged || (counter_f_tol > options.successive_f_tol)

        if !(converged && method isa Newton)
            update_h!(d, state, method) # only relevant if not converged
        end
        if tracing
            # update trace; callbacks can stop routine early by returning true
            stopped_by_callback = trace!(tr, d, state, iteration, method, options, time()-t0)
        end

        # Check time_limit; if none is provided it is NaN and the comparison
        # will always return false.
        _time = time()
        stopped_by_time_limit = _time-t0 > options.time_limit
        f_limit_reached = options.f_calls_limit > 0 && f_calls(d) >= options.f_calls_limit ? true : false
        g_limit_reached = options.g_calls_limit > 0 && g_calls(d) >= options.g_calls_limit ? true : false
        h_limit_reached = options.h_calls_limit > 0 && h_calls(d) >= options.h_calls_limit ? true : false

        if (f_increased && !options.allow_f_increases) || stopped_by_callback ||
            stopped_by_time_limit || f_limit_reached || g_limit_reached || h_limit_reached
            stopped = true
        end
    end # while

    after_while!(d, state, method, options)

    # we can just check minimum, as we've earlier enforced same types/eltypes
    # in variables besides the option settings
    Tf = typeof(value(d))
    f_incr_pick = f_increased && !options.allow_f_increases
    stopped_by =(f_limit_reached=f_limit_reached,
                 g_limit_reached=g_limit_reached,
                 h_limit_reached=h_limit_reached,
                 time_limit=stopped_by_time_limit,
                 callback=stopped_by_callback,
                 f_increased=f_incr_pick)
    return MultivariateOptimizationResults{typeof(method),T,Tx,typeof(x_abschange(state)),Tf,typeof(tr), Bool, typeof(stopped_by)}(method,
                                        initial_x,
                                        pick_best_x(f_incr_pick, state),
                                        pick_best_f(f_incr_pick, state, d),
                                        iteration,
                                        iteration == options.iterations,
                                        x_converged,
                                        Tf(options.x_abstol),
                                        Tf(options.x_reltol),
                                        x_abschange(state),
                                        x_relchange(state),
                                        f_converged,
                                        Tf(options.f_abstol),
                                        Tf(options.f_reltol),
                                        f_abschange(d, state),
                                        f_relchange(d, state),
                                        g_converged,
                                        Tf(options.g_abstol),
                                        g_residual(d),
                                        f_increased,
                                        tr,
                                        f_calls(d),
                                        g_calls(d),
                                        h_calls(d),
                                        ls_success,
                                        options.time_limit,
                                        _time-t0,
                                        stopped_by,
                                        )
end


text:
using Base.Test
using Mimi
using DataFrames

@defcomp foo1 begin
    index1 = Index()

    par1 = Parameter()
    par2::Bool = Parameter(index=[time,index1], description="description par 1")
    par3 = Parameter(index=[time])

    var1 = Variable()
    var2 = Variable(index=[time])
    var3 = Variable(index=[time,index1])
    var4::Bool = Variable(index=[3])
    var5 = Variable(index=[index1,4])
end

x1 = foo1(Float64, Dict{Symbol, Int}(:time=>10, :index1=>3))

@test x1.Dimensions.index1.start == 1
@test x1.Dimensions.index1.stop == 3

# Check variable types
@test isa(x1.Variables.var1, Float64)
@test isa(x1.Variables.var2, Array{Float64,1})
@test isa(x1.Variables.var3, Array{Float64,2})
@test isa(x1.Variables.var4, Array{Bool,1})
@test isa(x1.Variables.var5, Array{Float64,2})

# Check variable sizes
@test size(x1.Variables.var2,1)==10
@test size(x1.Variables.var3)==(10,3)
@test size(x1.Variables.var4,1)==3
@test size(x1.Variables.var5)==(3,4)

resetvariables(x1)

# Check all variables are defaulted
@test isnan(x1.Variables.var1)

@test_throws MethodError x1.Parameters.par1 = Array(Float64, 1,2)

x1.Parameters.par1 = 5.0
@test x1.Parameters.par1 == 5.0

m = Model()
setindex(m, :time, 20)
setindex(m, :index1, 5)
addcomponent(m, foo1)

@test in(:var1, variables(m, :foo1))

@test isa(getdataframe(m, :foo1, :var3), DataFrame)

ext: jl


text:
using Sobol
using Base.Test

# test integrand from Joe and Kuo paper (integrates to 1)
function testfunc(x::Vector{Float64})
    f = 1.0
    for j = 1:length(x)
        cj = j^0.3333333333333333333
        f *= (abs(4*x[j] - 2) + cj) / (1 + cj)
    end
    return f
end

function testint(N, n)
    s = SobolSeq(N)
    x = Array(Float64, N)
    skip(s, n)
    sum = 0.0
    for j = 1:n
        sum += testfunc(next!(s, x))
    end
    sum / n
end

# test results from Joe and Kuo (2003)
N = [50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1111]
n = [1009, 1997, 4001, 8009, 16001, 32003, 64007, 128021]
JoeKuo = [0.9743 0.9544 0.9838 0.9982 0.9995 0.9976 0.9969 0.9975;
          1.0157 0.9288 0.9649 0.9866 0.9985 0.9817 0.9955 0.9995;
          1.1109 0.8877 0.9070 1.0312 1.0216 0.9772 0.9744 0.9925;
          0.9770 0.8715 0.9305 0.9546 1.0060 0.9712 0.9770 0.9843;
          1.1780 0.8908 0.9710 0.9527 1.0203 0.9754 0.9779 0.9725;
          1.0567 0.8548 0.9700 0.9704 0.9768 0.9634 0.9745 0.9670;
          0.8735 0.7674 1.0130 0.9874 0.9677 0.9416 0.9580 0.9564;
          0.8752 0.8241 1.0067 1.0078 0.9987 0.9236 0.9517 0.9475;
          0.9139 0.8822 1.0926 1.0560 1.0210 0.9454 0.9330 0.9359;
          0.7779 1.0225 1.0593 1.0991 1.0470 0.9115 0.9688 0.9615;
          0.8348 1.0443 1.0710 1.0862 1.0483 0.9188 0.9735 0.9683;
          0.7524 0.9139 0.9592 0.9703 1.0116 0.9124 1.0411 0.9904]

for i = 1:length(N)
    println("testing integrals of dimension $(N[i])")
    for j = 1:length(n)
        println("  ... with $(n[j]) samples")
        t = testint(N[i], n[j])
        @test_approx_eq_eps t JoeKuo[i,j] 1e-4
    end
end


text:
export decode

mutable struct Decoder
    worddict::Dict
    chardict::Dict
    tagdict::Dict
    nn
end

struct Sample
    w
    batchdims_w
    c
    batchdims_c
    t
end

macro logtime()
    return :( Dates.format(now(), "yyyy-mm-dd HH:MM:SS") )
end

function create_batch(samples::Vector{Sample}, batchsize::Int)
    batches = Sample[]
    for i = 1:batchsize:length(samples)
        range = i:min(i+batchsize-1,length(samples))
        s = samples[range]
        w = cat(2, map(x -> x.w, s)...)
        c = cat(2, map(x -> x.c, s)...)
        batchdims_w = cat(1, map(x -> x.batchdims_w, s)...)
        batchdims_c = cat(1, map(x -> x.batchdims_c, s)...)
        t = s[1].t == nothing ? nothing : cat(1, map(x -> x.t, s)...)

        w = BACKEND(w)
        batchdims_w = sort(batchdims_w, rev=true)
        c = BACKEND(c)
        batchdims_c = BACKEND(batchdims_c)
        t = BACKEND(t)

        push!(batches, Sample(w, batchdims_w, c, batchdims_c, t))
    end
    batches
end

function Decoder(config::Dict, iolog)
    words = h5read(config["wordvec_file"], "words")
    wordembeds = h5read(config["wordvec_file"], "vectors")
    worddict = Dict(words[i] => i for i=1:length(words))
    chardict, tagdict = initvocab(config["train_file"])
    charembeds = Normal(0,0.01)(Float32, 20, length(chardict))
    traindata = readdata(config["train_file"], worddict, chardict, tagdict)
    testdata = readdata(config["test_file"], worddict, chardict, tagdict)
    nepochs = haskey(config, "nepochs") ? config["nepochs"] : 1
    nlayers = haskey(config, "nlayers") ? config["nlayers"] : 1
    winsize_c = haskey(config, "winsize_c") ? config["winsize_c"] : 2
    winsize_w = haskey(config, "winsize_w") ? config["winsize_w"] : 5
    droprate = haskey(config, "droprate") ? config["droprate"] : 0.2

    procname = @sprintf("NER.Convolution[%s]", haskey(config, "jobid") ? config["jobid"] : "-")
    @printf(iolog, "%s %s training - train:%d test:%d words:%d, chars:%d tags:%d\n", @logtime, procname,
            length(traindata), length(testdata), length(worddict), length(chardict), length(tagdict))
    @printf(iolog, "%s %s arguments - nepochs:%d, nlayers:%d droprate:%f winsize_c:%d winsize_w:%d\n", @logtime, procname,
            nepochs, nlayers, droprate, winsize_c, winsize_w)
    flush(iolog)

    nn = NN(wordembeds, charembeds, length(tagdict);
            nlayers=nlayers, winsize_c=winsize_c, winsize_w=winsize_w, droprate=droprate)

    info("#Training examples:\t$(length(traindata))")
    info("#Testing examples:\t$(length(testdata))")
    info("#Words1:\t$(length(worddict))")
    info("#Chars:\t$(length(chardict))")
    info("#Tags:\t$(length(tagdict))")

    batchsize = config["batchsize"]
    testdata = create_batch(testdata, batchsize)

    opt = SGD()
    for epoch = 1:nepochs
        println("Epoch:\t$epoch")

        #opt.rate = LEARN_RATE / BATCHSIZE
        opt.rate = config["learning_rate"] * batchsize / sqrt(batchsize) / (1 + 0.05*(epoch-1))
        println("Learning rate: $(opt.rate)")

        @printf(iolog, "%s %s begin epoch %d - learnrate:%.5f\n", @logtime, procname,
                epoch, opt.rate)
        flush(iolog)

        shuffle!(traindata)
        batches = create_batch(traindata, batchsize)
        prog = Progress(length(batches))
        loss = 0.0
        for i in 1:length(batches)
            s = batches[i]
            z = nn(s, true)
            loss += sum(z.data)
            params = gradient!(z)
            foreach(opt, params)
            ProgressMeter.next!(prog)
        end
        loss /= length(batches)
        println("Loss:\t$loss")

        # test
        println("Testing...")
        preds = Int[]
        golds = Int[]
        for s in testdata
            z = nn(s, false)
            append!(preds, z)
            append!(golds, s.t)
        end
        length(preds) == length(golds) || throw("Length mismatch: $(length(preds)), $(length(golds))")

        preds = BIOES.decode(preds, tagdict)
        golds = BIOES.decode(golds, tagdict)
        prec, recall, fval = fscore(golds, preds)

        @printf(iolog, "%s %s end epoch %d - loss:%.4e fval:%.5f prec:%.5f recall:%.5f\n", @logtime, procname,
                epoch, loss, fval, prec, recall)
        flush(iolog)
        println()
    end

    @printf(iolog, "%s %s training complete\n", @logtime, procname)

    Decoder(worddict, chardict, tagdict, nn)
end

function initvocab(path::String)
    chardict = Dict{String,Int}()
    tagdict = Dict{String,Int}()
    lines = open(readlines, path)
    for line in lines
        isempty(line) && continue
        items = Vector{String}(split(line,"\t"))
        word = strip(items[1])
        chars = Vector{Char}(word)
        for c in chars
            c = string(c)
            if haskey(chardict, c)
                chardict[c] += 1
            else
                chardict[c] = 1
            end
        end

        tag = strip(items[2])
        haskey(tagdict,tag) || (tagdict[tag] = length(tagdict)+1)
    end

    chars = String[]
    for (k,v) in chardict
        v >= 3 && push!(chars,k)
    end
    chardict = Dict(chars[i] => i for i=1:length(chars))
    chardict["UNKNOWN"] = length(chardict) + 1

    chardict, tagdict
end

function decode(dec::Decoder, config::Dict)
    testdata = readdata(config["test_file"], dec.worddict, dec.chardict, dec.tagdict)
    testdata = create_batch(testdata, 10)
    id2tag = Array{String}(length(dec.tagdict))
    for (k, v) in dec.tagdict
        id2tag[v] = k
    end

    preds = Int[]
    for x in testdata
        y = dec.nn(x, false)
        append!(preds, y)
    end

    lines = open(readlines, config["test_file"])
    i = 1
    for line in lines
        if !isempty(strip(line))
            tag = id2tag[preds[i]]
            println("$line\t$tag")
            i += 1
        else
            println("")
        end
    end
end

function readdata(path::String, worddict::Dict, chardict::Dict, tagdict::Dict)
    samples = Sample[]
    words, tags = String[], String[]
    unkword = worddict["UNKNOWN"]
    unkchar = chardict["UNKNOWN"]

    lines = open(readlines, path)
    push!(lines, "")
    for i = 1:length(lines)
        line = lines[i]
        if isempty(line)
            isempty(words) && continue
            wordids = Int[]
            charids = Int[]
            batchdims_c = Int[]
            for w in words
                #w0 = replace(lowercase(w), r"[0-9]", '0')
                id = get(worddict, lowercase(w), unkword)
                push!(wordids, id)

                chars = Vector{Char}(w)
                ids = map(chars) do c
                    get(chardict, string(c), unkchar)
                end
                append!(charids, ids)
                push!(batchdims_c, length(ids))
            end
            batchdims_w = [length(words)]
            w = reshape(wordids, 1, length(wordids))
            c = reshape(charids, 1, length(charids))
            t = isempty(tags) ? nothing : map(t -> tagdict[t], tags)
            push!(samples, Sample(w,batchdims_w,c,batchdims_c,t))
            empty!(words)
            empty!(tags)
        else
            items = Vector{String}(split(line,"\t"))
            word = strip(items[1])
            @assert !isempty(word)
            push!(words, word)
            if length(items) >= 2
                tag = strip(items[2])
                push!(tags, tag)
            end
        end
    end
    samples
end

function fscore(golds::Vector{T}, preds::Vector{T}) where T
    set = intersect(Set(golds), Set(preds))
    count = length(set)
    prec = round(count/length(preds), 5)
    recall = round(count/length(golds), 5)
    fval = round(2*recall*prec/(recall+prec), 5)
    println("Prec:\t$prec")
    println("Recall:\t$recall")
    println("Fscore:\t$fval")
    prec, recall, fval
end


text:
# 载入共同的组件
include("common.jl")

# 载入空域掩膜算法(doi: 10.1109/JDT.2015.2479646)
include("sdm.jl")
# 载入冗余散斑降噪算法(doi: 10.1364/AO.390500)
include("rse.jl")
# 载入低维重建法(doi: 10.1364/AO.414773)
include("ldr.jl")


text:
@testset "plot_mesh.jl" begin

end


text:
using ITensors
N = 100
sites = siteinds("S=1",N)
ampo = AutoMPO()
for j=1:N-1
    ampo .+= ("Sz",j,"Sz",j+1)
    ampo .+= (0.5,"S+",j,"S-",j+1)
    ampo .+= (0.5,"S-",j,"S+",j+1)
end
H = MPO(ampo,sites)
psi0 = randomMPS(sites,10)
sweeps = Sweeps(5)
maxdim!(sweeps, 10,20,100,100,200)
cutoff!(sweeps, 1E-11)
energy, psi = dmrg(H,psi0, sweeps)


text:
module computer

export Context, io, run

using DelimitedFiles, Match

mutable struct Context
    program::Array{Int64,1}
    instruction::Int64
    relative_base::Int64
    debug::Bool
    buffer::String
end

function io()
    return (Channel{Int64}(Inf), Channel{Int64}(Inf))
end

function run(program::Array{Int64,1}, debug=false)
    stdin, stdout = io()
    run(program, stdin, stdout, debug)
end

function run(program::Array{Int64,1}, input::Int64, debug=false)
    stdin, stdout = io()
    put!(stdin, input)
    run(program, stdin, stdout, debug)
end

function run(program::Array{Int64,1}, stdin::Channel{Int64}, stdout::Channel{Int64}, debug=false)
    c = Context(program, 1, 0, debug, "")
    while execute(c, stdin, stdout); end
    return program
end

function execute(context::Context, stdin::Channel{Int64}, stdout::Channel{Int64})
    dbg(context, "@@")
    dbg(context, lpad(context.program[context.instruction], 8))
    op = rem(context.program[context.instruction], 100)
    @match op begin
        # ADD
        1 => begin
            dbg(context, rpad("  ADD", 8))
            x = rval(context, 1)
            y = rval(context, 2)
            wval(context, 3, x + y)
            context.instruction += 4
        end
        # MULTIPLY
        2 => begin
            dbg(context, rpad("  MULT", 8))
            x = rval(context, 1)
            y = rval(context, 2)
            wval(context, 3, x * y)
            context.instruction += 4
        end
        # INPUT
        3 => begin
            dbg(context, rpad("  STDIN", 8))
            val = take!(stdin)
            wval(context, 1, val)
            context.instruction += 2
        end
        # OUTPUT
        4 => begin
            dbg(context, rpad("  STDOUT", 8))
            val = rval(context, 1)
            put!(stdout, val)
            context.instruction += 2
        end
        # JUMP IF TRUE
        5 => begin
            dbg(context, rpad("  JMPT", 8))
            val = rval(context, 1)
            dest = rval(context, 2)
            val != 0 ?
                context.instruction = dest + 1 :
                context.instruction += 3
        end
        # JUMP IF FALSE
        6 => begin
            dbg(context, rpad("  JMPF", 8))
            val = rval(context, 1)
            dest = rval(context, 2)
            val == 0 ?
                context.instruction = dest + 1 :
                context.instruction += 3
        end
        # LESS THAN
        7 => begin
            dbg(context, rpad("  CMPLT", 8))
            x = rval(context, 1)
            y = rval(context, 2)
            x < y ? wval(context, 3, 1) : wval(context, 3, 0)
            context.instruction += 4
        end
        # EQUALS
        8 => begin
            dbg(context, rpad("  CMPEQ", 8))
            x = rval(context, 1)
            y = rval(context, 2)
            x == y ? wval(context, 3, 1) : wval(context, 3, 0)
            context.instruction += 4
        end
        # ADJUST RELATIVE BASE
        9 => begin
            dbg(context, rpad("  REL", 8))
            x = rval(context, 1)
            context.relative_base += x
            context.instruction += 2
        end
        # HALT
        99 => begin
            dbg(context, rpad("  HALT", 8))
            context.debug && println(context.buffer)
            context.buffer = ""
            return false
        end
        _ => throw(ErrorException("unknown op code: $op"))
    end
    context.debug && println(context.buffer)
    context.buffer = ""
    return true
end

function rval(context::Context, parameter::Int)
    program = context.program
    instruction = context.instruction
    rel = context.relative_base
    op = program[instruction]
    address = instruction + parameter
    mode = rem(div(op, ^(10, parameter + 1)), 10)
    dbg(context, lpad(address, 8))
    dbg(context, rpad(glyph(mode), 2))
    @match mode begin
        0 => begin
            idx = read(context, instruction + parameter)
            return read(context, idx + 1)
        end
        1 => return read(context, instruction + parameter)
        2 => begin
            offset = read(context, instruction + parameter)
            return read(context, rel + offset + 1)
        end
        _ => throw(ErrorException("unknown parameter mode: $mode"))
    end
end

function wval(context::Context, parameter::Int, value::Int64)
    program = context.program
    instruction = context.instruction
    rel = context.relative_base
    op = program[instruction]
    address = instruction + parameter
    mode = rem(div(op, ^(10, parameter + 1)), 10)
    dbg(context, lpad(address, 8))
    dbg(context, rpad(glyph(mode), 2))
    @match mode begin
        0 => begin
            idx = read(context, instruction + parameter)
            write(context, idx + 1, value)
        end
        1 => throw(ErrorException("invalid parameter mode: $mode"))
        2 => begin
            offset = read(context, instruction + parameter)
            write(context, rel + offset + 1, value)
        end
        _ => throw(ErrorException("unknown parameter mode: $mode"))
    end
end

function read(context::Context, index::Int64)
    return index > length(context.program) ? 0 : context.program[index]
end

function write(context::Context, index::Int64, value::Int64)
    if index > length(context.program)
        mem = zeros(Int64, index - length(context.program))
        context.program = cat(context.program, mem, dims=(1))
    end
    context.program[index] = value
end

function dbg(context, string)
    if context.debug
        context.buffer *= string
    end
end

function glyph(mode::Int)
    @match mode begin
        0 => return "x"
        1 => return "i"
        2 => return "r"
    end
end

function load_program(path::String)
    return vec(readdlm(path, ',', Int64))
end

end


text:
doc"""
    Skellam(μ1, μ2)

A *Skellam distribution* describes the difference between two independent [`Poisson`](:func:`Poisson`) variables, respectively with rate `μ1` and `μ2`.

$P(X = k) = e^{-(\mu_1 + \mu_2)} \left( \frac{\mu_1}{\mu_2} \right)^{k/2} I_k(2 \sqrt{\mu_1 \mu_2}) \quad \text{for integer } k$

where $I_k$ is the modified Bessel function of the first kind.

```julia
Skellam(mu1, mu2)   # Skellam distribution for the difference between two Poisson variables,
                    # respectively with expected values mu1 and mu2.

params(d)           # Get the parameters, i.e. (mu1, mu2)
```

External links:

* [Skellam distribution on Wikipedia](http://en.wikipedia.org/wiki/Skellam_distribution)
"""
immutable Skellam <: DiscreteUnivariateDistribution
    μ1::Float64
    μ2::Float64

    function Skellam(μ1::Real, μ2::Real)
        @check_args(Skellam, μ1 > zero(μ1) && μ2 > zero(μ2))
        new(μ1, μ2)
    end
    Skellam(μ::Real) = Skellam(μ, μ)
    Skellam() = new(1.0, 1.0)
end

@distr_support Skellam -Inf Inf

### Parameters

params(d::Skellam) = (d.μ1, d.μ2)

### Statistics

mean(d::Skellam) = d.μ1 - d.μ2

var(d::Skellam) = d.μ1 + d.μ2

skewness(d::Skellam) = mean(d) / (var(d)^1.5)

kurtosis(d::Skellam) = 1.0 / var(d)

### Evaluation

function logpdf(d::Skellam, x::Int)
    μ1, μ2 = params(d)
    - (μ1 + μ2) + (x / 2.0) * log(μ1 / μ2) + log(besseli(x, 2.0 * sqrt(μ1) * sqrt(μ2)))
end

pdf(d::Skellam, x::Int) = exp(logpdf(d, x))

function mgf(d::Skellam, t::Real)
    μ1, μ2 = params(d)
    exp(-(μ1 + μ2) + μ1 * exp(t) + μ2 * exp(-t))
end

function cf(d::Skellam, t::Real)
    μ1, μ2 = params(d)
    exp(-(μ1 + μ2) + μ1 * cis(t) + μ2 * cis(-t))
end

### Sampling

rand(d::Skellam) = rand(Poisson(d.μ1)) - rand(Poisson(d.μ2))


text:
# make a dataframe of the relative errors on a broad set of test cases
using DataFrames
df = DataFrame(err = Float64[], dualErr=Float64[], fnc=String[], bc=WaveletBoundary[], cwt = ContinuousWavelets.ContWaveClass[], β = Int[], aveLen = Int[], invMeth = ContinuousWavelets.InverseType[], N = Int[])
@testset "Inversion" begin
    bcs =(PerBoundary(),SymBoundary(), ZPBoundary())
    cwts = (cDb2,dog1,dog2,paul1,morl)
    βs =(1, 2)
    averagingLengths=(-1,0,1)
    typesOfTestFunctions = ["Doppler" "Bumps" "Blocks" "HeaviSine" "just Core"]
    inversionMethods = [NaiveDelta(), PenroseDelta(), DualFrames()]
    ns = (128, 1382,2039)
    @testset "length $n, with type $wave, bc $bc, β=$β, ave=$(ave) ex=$(testF)" for n in ns, wave in cwts, bc in bcs, β in βs, ave in averagingLengths, testF in typesOfTestFunctions, inverseType in inversionMethods
        if testF =="just Core"
            x̂ = zeros(ComplexF64, n>>1+1); x̂[n>>2 .+ (-10:10)] = 5*(11.5 .- abs.(range(-10/2,10/2,length=21))) .* randn(ComplexF64, 21)
            x = irfft(x̂, n) # test function that is definitely solidly in the "safe" range
        else
            x = testfunction(n,testF)
        end
        if inverseType isa NaiveDelta
            errTol = 3 # naiveDelta really isn't very good at all
        elseif inverseType isa DualFrames && wave isa Morlet
            errTol = 1 # generally good unless the frames are particularly bad at the end
        else
            errTol = .01 # everything else should have error less than 1%
        end

        wav = wavelet(wave, β=β,boundary=bc,averagingLength = ave)
        res = 3
        with_logger(ConsoleLogger(stderr, Logging.Error)) do
            res = cwt(x, wav)
            xRecon = real.(ContinuousWavelets.icwt(res, wav, inverseType))
            err = norm(xRecon-x)/norm(x)
            dualWav,dualErr = ContinuousWavelets.getDualCoverage(n, wav, DualFrames());
            push!(df, (err, dualErr, testF, bc, wave, β, ave, inverseType, n))
            #@test  err < errTol
        end
    end
end
#wave = morl; bc = PerBoundary(); β=2; ave = -1; n = 1382; testF = "just Core"; inverseType = DualFrames()
# put the data frame in order of increasing error
sort!(df, :err)
# see the non-naive entries where the error is greater than 50%
show(df[(df.err .≥ .5) .& (map(x->x!=NaiveDelta(), df.invMeth)), :], allrows=true)
sort!(df[(df.err .≥ .5) .& (map(x->x!=NaiveDelta(), df.invMeth)), :], :dualErr)
show(last(df[(df.err .≥ .5) .& (map(x->x!=NaiveDelta(), df.invMeth)), :],50), allrows=true)
show(last(df[(df.N .== 128) .& (map(x->x==ZPBoundary(), df.bc)), :], 50), allrows=true)
show(last(df[(map(x->x==ZPBoundary(), df.bc)), :], 50), allrows=true)

inverseType=PenroseDelta()

βp = computeNonNegativeDualWeights(Ŵ)
heatmap(abs.(res)')
wav = wavelet(wave, β=β,boundary=bc,averagingLength = ave)
Ŵ,ω = computeWavelets(n,wav)
plot(Ŵ)
dual,dualErrr = ContinuousWavelets.getDualCoverage(n,wav,inverseType)
plot(dual)
Ŵdual = conj.(Ŵ ./ [norm(ŵ)^2 for ŵ in eachslice(Ŵ,dims=1)])
[norm(ŵ)^2 for ŵ in eachslice(Ŵ,dims=1)]
plot(Ŵ,legend=false)
plot(Ŵdual[:,1:end-0],legend=false)
res = cwt(x, wav)
xRecon = real.(ContinuousWavelets.icwt(res,wav,inverseType))
plot([x xRecon .*norm(x)/norm(xRecon)],labels=["original" "reconstructed"])
plot([x xRecon],labels=["original" "reconstructed"])
norm(xRecon), norm(x)
norm(x-xRecon .*norm(x)/norm(xRecon))
norm(x-xRecon)/norm(x)
plot(imag(xRecon))
plot(abs.(xRecon))
plot(abs.(fft(xRecon)))
plot([x xRecon],labels=["original" "reconstructed"])
plot([x x-xRecon .*norm(x)/norm(xRecon)],labels=["original" "error"])
plot([x x-xRecon],labels=["original" "error"])
plot([abs.(rfft(x)) abs.(rfft(xRecon))],labels=["original" "reconstruction"])
plot([abs.(rfft([x; reverse(x)])) abs.(rfft([xRecon; reverse(xRecon)]))],labels=["original" "reconstruction"])
plot([abs.(rfft(x)) abs.(rfft(x-xRecon .*norm(x)/norm(xRecon)))],labels=["original" "error"])
plot(plot(x), plot(xRecon))
plot(plot(abs.(rfft(x))), plot(abs.(rfft(xRecon))))
plot(plot(abs.(rfft([x; reverse(x)]))), plot(abs.(rfft(xRecon))))
plot([x xRecon])
plot([x x-xRecon .*norm(x)/norm(xRecon)])


text:
function _lattice_primitives(::Bravais, φ, constants::NTuple{2})
    sinφ, cosφ = sincos(φ)
    return constants[1]*Point(cosφ,sinφ), constants[2]*Point(-sinφ,cosφ)
end
function _lattice_primitives(::Cartesian, φ, constants::NTuple{2})
    sinφ, cosφ = sincos(φ[1])
    return constants[1]*Point{Cartesian}(cosφ,sinφ), constants[2]*Point{Cartesian}(-sinφ,cosφ)
end
function _lattice_primitives(::Polar, φ, constants::NTuple{2})
    return Point{Polar}(constants[1],φ[1]), Point{Polar}(constants[2],φ[1]+π/2)
end
function _lattice_primitives(::Spherical, φ, constants::NTuple{2})
    return Point{Spherical}(constants[1],φ[1]), Point{Spherical}(constants[2],φ[1]+π/2)
end

Base.getindex(lat::Lattice{2},i::CartesianIndex) = getindex(lat,i[1],i[2])
Base.getindex(lat::Lattice{2},ij::Tuple) = getindex(lat,ij[1],ij[2])
Base.getindex(lat::Lattice{2},i,j) = map((i,j)->lat[i,j],i,j)

Base.getindex(lat::Lattice{2,Bravais}   ,i::Real, j::Real)  = lat.origin + i*lat.a*lat.e1 + j*lat.b*lat.e2
Base.getindex(lat::Lattice{2,Cartesian} ,i::Real, j::Real)  = lat.origin + i*lat.dx*lat.e1 + j*lat.dy*lat.e2
Base.getindex(lat::Lattice{2,Polar}     ,i::Real, j::Real)  = lat.origin + Point{Polar}(lat.r0 + i*lat.dr, j*lat.dϕ+lat.ϕ0)
Base.getindex(lat::Lattice{2,Spherical}     ,i::Real, j::Real)  = lat.origin + Point{Spherical}(lat.r0 + i*lat.dr, j*lat.dθ+lat.θ0)

function latticeindex(lat::Lattice{2,C}, p::Point{2}) where C<:Union{Bravais,Cartesian}
    p = p-lat.origin
    x = lat.e1[1]*lat.e2[1] + lat.e1[2]*lat.e2[2]
    u1 = (lat.e1 - x*lat.e2)/(1-x^2)
    u2 = (lat.e2 - x*lat.e1)/(1-x^2)
    if C<:Bravais
        return (u1[1]*p.x+u1[2]*p.y)/lat.a, (u2[1]*p.x+u2[2]*p.y)/lat.b
    else
        return (u1[1]*p.x+u1[2]*p.y)/lat.dx, (u2[1]*p.x+u2[2]*p.y)/lat.dy
    end
end

function latticeindex(lat::Lattice{2,Polar}, p::Point{2})
    p = p-lat.origin
    return p.r/lat.dr - 1/2, (p.ϕ-lat.ϕ0)/lat.dϕ
end

function latticeindex(lat::Lattice{2,Spherical}, p::Point{2})
    p = Cartesian(p-lat.origin)
    return p.r/lat.dr - 1/2, (π/2 - p.ϕ-lat.θ0)/lat.dθ
end

################################################################################
# Bravais
function Base.getproperty(lat::Lattice{2,Bravais}, sym::Symbol)
    if sym == :a
        return getfield(lat,:constants)[1]
    elseif sym == :b
        return getfield(lat,:constants)[2]
    elseif sym == :e1
        return getfield(lat,:e)[1]
    elseif sym == :e2
        return getfield(lat,:e)[2]
    else
        return getfield(lat,sym)
    end
end

function Base.propertynames(::Lattice{2,Bravais}, private=false)
    if private
        return fieldnames(Lattice)
    else
        return (:a, :b, :origin, :e1, :e2)
    end
end

################################################################################
# Cartesian
function Base.getproperty(lat::Lattice{2,Cartesian}, sym::Symbol)
    if sym == :dx
        return getfield(lat,:constants)[1]
    elseif sym == :dy
        return getfield(lat,:constants)[2]
    elseif sym == :x0
        return getfield(lat,:origin)[1]
    elseif sym == :y0
        return getfield(lat,:origin)[2]
    elseif sym == :e1
        return getfield(lat,:e)[1]
    elseif sym == :e2
        return getfield(lat,:e)[2]
    elseif Base.sym_in(sym,(:ϕ,:φ,:phi,:𝜑,:angle))
        return getfield(lat,:e)[1].ϕ
    else
        return getfield(lat,sym)
    end
end

function Base.propertynames(::Lattice{2,Cartesian},private=false)
    if private
        return fieldnames(Lattice)
    else
        return (:dx, :dy, :e1, :e2, :origin, :ϕ)
    end
end

################################################################################
# Polar
function Base.getproperty(lat::Lattice{2,Polar}, sym::Symbol)
    if sym == :dr
        return getfield(lat,:constants)[1]
    elseif Base.sym_in(sym,(:dϕ,:dφ,:dphi))
        return getfield(lat,:constants)[2]
    elseif sym == :r0
        return getfield(lat,:constants)[1]/2
    elseif Base.sym_in(sym,(:ϕ0,:φ0,:phi0,:ϕ,:φ,:𝜑,:𝜑0,:phi,:angle))
        return getfield(lat,:e)[1].ϕ
    else
        return getfield(lat,sym)
    end
end

function Base.propertynames(::Lattice{2,Polar}, private=false)
    if private
        return fieldnames(Lattice)
    else
        return (:dr, :dϕ, :r0, :ϕ0, :origin)
    end
end

################################################################################
# Spherical
function Base.getproperty(lat::Lattice{2,Spherical}, sym::Symbol)
    if sym == :dr
        return getfield(lat,:constants)[1]
    elseif Base.sym_in(sym,(:dθ,:dϑ,:d𝜃,:dtheta))
        return getfield(lat,:constants)[2]
    elseif sym == :r0
        return getfield(lat,:constants)[1]/2
    elseif Base.sym_in(sym,(:θ0,:ϑ0,:theta0,:θ,:𝜃,:𝜃0,:ϑ,:angle))
        return acsc(getfield(lat,:e)[1][2]/getfield(lat,:e)[1][1])
    else
        return getfield(lat,sym)
    end
end

function Base.propertynames(::Lattice{2,Spherical}, private=false)
    if private
        return fieldnames(Lattice)
    else
        return (:dr, :dθ, :r0, :θ0, :origin)
    end
end

################################################################################
# PRETTY PRINTING
import ..PRINTED_COLOR_NUMBER
import ..PRINTED_COLOR_WARN
import ..PRINTED_COLOR_DARK

function Base.show(io::IO, lat::Lattice{2,TYPE}) where TYPE
    get(io,:tabbed2,false) ? print(io,"\t\t") : nothing
    if TYPE <: Bravais
        printstyled(io,"2D Bravais Lattice ",color=PRINTED_COLOR_DARK)
        print(io,"(a=")
        printstyled(io, fmt("3.4f",lat.a),color=PRINTED_COLOR_NUMBER)
        print(io,", b=")
        printstyled(io, fmt("3.4f",lat.b),color=PRINTED_COLOR_NUMBER)
    elseif TYPE <: Cartesian
        printstyled(io,"2D Cartesian Lattice ",color=PRINTED_COLOR_DARK)
        print(io,"(dx=")
        printstyled(io, fmt("3.4f",lat.dx),color=PRINTED_COLOR_NUMBER)
        print(io,", dy=")
        printstyled(io, fmt("3.4f",lat.dy),color=PRINTED_COLOR_NUMBER)
    elseif TYPE <: Polar
        printstyled(io,"2D Polar Lattice ",color=PRINTED_COLOR_DARK)
        print(io,"(dr=")
        printstyled(io, fmt("3.4f",lat.dr),color=PRINTED_COLOR_NUMBER)
        print(io,", dϕ=")
        printstyled(io, fmt("3.4f",lat.dϕ),color=PRINTED_COLOR_NUMBER)
    elseif TYPE <: Spherical
        printstyled(io,"2D Spherical Lattice (w/ azimuthal symmetry) ",color=PRINTED_COLOR_DARK)
        print(io,"(dr=")
        printstyled(io, fmt("3.4f",lat.dr),color=PRINTED_COLOR_NUMBER)
        print(io,", dθ=")
        printstyled(io, fmt("3.4f",lat.dθ),color=PRINTED_COLOR_NUMBER)
    end
    print(io, ", origin=")
    printstyled(io, lat.origin, color=PRINTED_COLOR_NUMBER)
    if !iszero(lat.angle)
        print(io,", angle=")
        printstyled(io, fmt("3.2f",180*lat.angle[1]/π),"°",color=PRINTED_COLOR_NUMBER)
    end
    print(io,")")
end

################################################################################
# PLOTTING

import ..SHAPE_COLOR
import ..SHAPE_FILL_ALPHA

@recipe function f(lat::Lattice{2,Cartesian},n1=-2,n2=2)
    (lat, n1:n2, n1:n2)
end
@recipe function f(lat::Lattice{2,Cartesian},nx::AbstractVector,ny::AbstractVector)
    aspectratio --> 1
    lat[[i for i ∈ nx, j ∈ ny],[j for i ∈ nx, j ∈ ny]][:]
end
@recipe function f(lat::Lattice{2,Polar},nr=5,nϕ=floor(Int,2π/lat.dϕ))
    (lat,0:nr,0:nϕ)
end
@recipe function f(lat::Lattice{2,Polar},nr::AbstractVector,nϕ::AbstractVector)
    aspectratio --> 1
    lat[[i for i ∈ nr, j ∈ nϕ],[j for i ∈ nr, j ∈ nϕ]][:]
end
@recipe function f(lat::Lattice{2,Spherical},nr=5,nθ=floor(Int,2π/lat.dθ))
    (lat,0:nr,0:nθ)
end
@recipe function f(lat::Lattice{2,Spherical},nr::AbstractVector,nθ::AbstractVector)
    aspectratio --> 1
    lat[[i for i ∈ nr, j ∈ nθ],[j for i ∈ nr, j ∈ nθ]][:]
end
@recipe function f(lat::Lattice{2,Bravais},n1=-2,n2=2)
    (lat, n1:n2, n1:n2)
end
@recipe function f(lat::Lattice{2,Bravais},na::AbstractVector,nb::AbstractVector)
    @series begin
        aspectratio --> 1
        lat[[i for i ∈ na, j ∈ nb],[j for i ∈ na, j ∈ nb]][:]
    end
    @series begin
        alpha --> 0
        seriestype --> :path
        fillcolor --> SHAPE_COLOR
        fillrange --> mean(lat[na,nb]).y
        fillalpha --> SHAPE_FILL_ALPHA
        aspect_ratio --> 1
        legend --> false
        na0 = round(mean(na))
        nb0 = round(mean(nb))
        [lat[na0,nb0],lat[na0+1,nb0],lat[na0+1,nb0+1],lat[na0,nb0+1],lat[na0,nb0]]
    end
end

# R = @SMatrix [  cosφ    -sinφ;
                # sinφ    cosφ]
# Rᵀ = @SMatrix [  cosφ    sinφ;
#                 -sinφ    cosφ]
# return (Point(R[:,1].data),Point(R[:,2].data)), R, Rᵀ


text:
struct PoolOracle <: Oracle
    labels::Vector{Symbol}
end

function PoolOracle(data, labels, params::Dict{Symbol, Any})
    return PoolOracle(labels)
end

function ask_oracle(oracle::PoolOracle, query_ids::Vector{Int})::Vector{Symbol}
    return oracle.labels[query_ids]
end


text:
#!/usr/bin/env julia

using HDF5

coolName = "MB286.mcool"

# Read data from a cool file:
#h5open(coolName) do cin
#	data2 = read(cin,"/resolutions/1000/pixels/bin1_id")
#	println(length(data2))
#end

# You can create objects for different levels of the hdf5
# like so:
h5open(coolName) do cin
	onek=cin["resolutions/1000"]
	pixels = onek["pixels"]
	data3 = read(pixels,"bin1_id")
	println(length(data3))
end

h5open(coolName) do cin
	onek=cin["resolutions/1000"]
	keyexists = haskey(attributes(onek),"bin-size")
	binsize=read(attributes(onek)["bin-size"])
	println(binsize)
end

attrs = h5readattr(coolName,"resolutions/1000")
println(attrs)

#data1 = h5read(coolName,"/resolutions/1000/pixels/bin1_id")
#println(length(data1))

#println(cin["resolutions"])

# bin_size=h5.getIntArrayAttribute("/","bin-size")[0]
# binsizename=attributes(cin["resolutions"]["1000"])["bin-size"]

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
include("get_sepsis_score.jl")

using DelimitedFiles

function load_challenge_data(file)
    (data, header) = readdlm(file, '|', header=true)

    # Ignore SepsisLabel column if present.
    if header[end] == "SepsisLabel"
        header = header[1:end-1]
        data = data[:, 1:end-1]
    end

    return data
end

function save_challenge_predictions(file, scores, labels)
    f = open(file, "w")
    write(f, "PredictedProbability|PredictedLabel\n")
    writedlm(f, hcat(scores, labels), '|')
    close(f)
end

function driver(input_directory, output_directory)
    # Find files.
    files = []
    for f in readdir(input_directory)
        if isfile(joinpath(input_directory, f)) && !startswith(lowercase(f), ".") && endswith(lowercase(f), "psv")
            push!(files, f)
        end
    end

    if !isdir(output_directory)
        mkdir(output_directory)
    end

    # Load model.
    println("Loading sepsis model...")
    model = load_sepsis_model()

    # Iterate over files.
    println("Predicting sepsis labels...")
    num_files = length(files)
    for (i, f) in enumerate(files)
        println("    $i/$num_files...")

        # Load data.
        input_file = joinpath(input_directory, f)
        data = load_challenge_data(input_file)

        # Make predictions.
        num_rows = size(data, 1)
        scores = zeros(Float64, num_rows)
        labels = zeros(Int, num_rows)
        for t = 1:num_rows
            scores[t], labels[t] = get_sepsis_score(data[1:t, :], model)
        end

        # Save results.
        output_file = joinpath(output_directory, f)
        save_challenge_predictions(output_file, scores, labels)
    end

    println("Done.")
end

# Parse arguments.
if length(ARGS) != 2
    error("Include the input and output directories as arguments, e.g., julia driver.jl input output.")
end

driver(ARGS[1], ARGS[2])


text:
@testset "Attribute" begin

@testset "Constructor" begin
    attribute = Attribute( "Attribute" )
    @test all( [attribute.name == "Attribute",
        isempty( attribute.possibleValues ), isempty( attribute.initValues ),
        isempty( attribute.initValueWeights )] )
end  # @testset "Constructor"

attribute = Attribute( "Attribute" )

@testset "function addPossibleAttributeValue!" begin
    addPossibleAttributeValue!( attribute, "bim" )
    @test attribute.possibleValues == ["bim"]
    addPossibleAttributeValue!( attribute, "bam", "boom", "bim" )
    @test all( ["bim" ∈ attribute.possibleValues,
        "boom" ∈ attribute.possibleValues,
        length( attribute.possibleValues ) == 3] )
    @test !addPossibleAttributeValue!( attribute, "bim", "boom" )
end  # @testset "function addPossibleAttributeValue!"

@testset "function removePossibleAttributeValue!" begin
    removePossibleAttributeValue!( attribute, "bing", "bim", "bang" )
    @test ( length( attribute.possibleValues ) == 2 ) &&
        ( "bim" ∉ attribute.possibleValues )
    @test !removePossibleAttributeValue!( attribute, "boop", "beep" )
end  # @testset "function removePossibleAttributeValue!"

@testset "function clearPossibleAttributeValues!" begin
    @test clearPossibleAttributeValues!( attribute )
    @test isempty( attribute.possibleValues )
end  # @testset "function clearPossibleAttributeValues!"

@testset "function setPossibleAttributeValues!" begin
    setPossibleAttributeValues!( attribute, ["beep", "boop"] )
    @test attribute.possibleValues == ["beep", "boop"]
    setPossibleAttributeValues!( attribute, ["foo", "bar"] )
    @test attribute.possibleValues == ["foo", "bar"]
end  # @testset "function setPossibleAttributeValues!"

@testset "function addInitialAttributeValue!" begin
    addInitialAttributeValue!( attribute, ("foo", 2.0) )
    @test ( attribute.initValues == ["foo"] ) &&
        ( attribute.initValueWeights == [2.0] )
    addInitialAttributeValue!( attribute, ("beep", 5.0), ("bar", -3.0),
        ("foo", 4.0) )
    @test all( [attribute.initValues == ["foo", "beep"],
        attribute.initValueWeights == [4.0, 5.0],
        "beep" ∈ attribute.possibleValues] )
    addInitialAttributeValue!( attribute, "boop", 1 )
    @test ( "boop" ∈ attribute.initValues ) &&
        ( "boop" ∈ attribute.possibleValues )
    @test !addInitialAttributeValue!( attribute, ("foo", -2.0), ("foo", 1.5) )
end  # @testset "function addInitialAttributeValue!"

@testset "function removeInitialAttributeValue!" begin
    removeInitialAttributeValue!( attribute, "beep" )
    @test ( "beep" ∉ attribute.initValues ) &&
        ( "beep" ∈ attribute.possibleValues )
    removeInitialAttributeValue!( attribute, "beep", "boop", "flub" )
    @test attribute.initValues == ["foo"]
    @test !removeInitialAttributeValue!( attribute, "beep" )
end  # @testset "function removeInitialAttributeValue!"

@testset "function clearInitialAttributeValues!" begin
    @test clearInitialAttributeValues!( attribute )
    @test isempty( attribute.initValues ) &&
        isempty( attribute.initValueWeights )
end  # @testset "function clearInitialAttributeValues!"

@testset "function setInitialAttributeValues!" begin
    setInitialAttributeValues!( attribute, Dict( "foo" => 3.0, "boop" => 2.0 ) )
    @test all( [length( attribute.initValues ) == 2,
        length( attribute.initValueWeights ) == 2, "foo" ∈ attribute.initValues,
        "boop" ∈ attribute.initValues] )
    setInitialAttributeValues!( attribute,
        Dict( "boop" => 2.0, "boom" => 3.0 ) )
    @test "boom" ∈ attribute.possibleValues
    @test !setInitialAttributeValues!( attribute,
        Dict( "foo" => -3.0, "bar" => -4.0 ) )
    @test length( attribute.initValues ) == 2
    setInitialAttributeValues!( attribute, ("foo", 4.0), ("beep", -3.0),
        ("boop", 1.0) )
    @test ( length( attribute.initValues ) == 2 ) &&
        ( ["foo", "boop"] ⊆ attribute.initValues )
    @test !setInitialAttributeValues!( attribute, ("foo", 2.0), ("boop", 6.0),
        ("foo", 1.0) )
    setInitialAttributeValues!( attribute, ["foo", "beep", "boop"],
        [2.0, 3.0, 1.0] )
    @test ( length( attribute.initValues ) == 3 ) &&
        ( ["foo", "beep", "boop"] ⊆ attribute.initValues )
    @test !setInitialAttributeValues!( attribute, ["foo"], [2.0, 5.0] )
    @test !setInitialAttributeValues!( attribute, ["foo", "boop", "foo"],
        [3.0, 4.0, 2.0] )
end  # @testset "function setInitialAttributeValues!"

end  # @testset "Attribute"


text:
module PackedReadWrite

export @reflect, read_struct, write_struct

"""
    @reflect(type)
implements methods for `fieldtypes` and `fieldnames` as constants for this `type`.
This makes other methods that use e.g. `fieldtypes(type)` (like `read_struct` does) very fast as it can
now totally compile.

`type` must be a struct type.
"""
macro reflect(type)

    if !Base.eval(__module__, :(isstructtype($type)))
        error("Can only extend field... methods for struct types.")
    end

    types = Base.eval(__module__, :(fieldtypes($type)))
    names = Base.eval(__module__, :(fieldnames($type)))
    num_fields = length(types)

    constructor_ok = Base.eval(__module__, :(hasmethod($type, $types)))

    Base.eval(__module__, quote
        Base.fieldtypes(::Type{T}) where T <: $type = $types
        Base.fieldnames(::Type{T}) where T <: $type = $names
        Base.fieldcount(::Type{T}) where T <: $type = $num_fields

        PackedReadWrite.hasfullconstructor(::Type{T}) where T <: $type = $constructor_ok
    end)

    return
end

function hasfullconstructor(::Type{T})::Bool where T
    hasmethod(T, fieldtypes(T))
end

"""
    read_struct(io::IO, ::Type{T})
Reads a struct type `T` from `io`, field by field. In this version, it requires for `T`
to have a constructor accepting the fields in the same order as they are defined in.

Optionally call `@reflect T` in your module to speed up (>10x) this method for type instances of type `T`.
"""
function read_struct(io::IO, ::Type{T}) where T

    if !(isstructtype(T))
        error("`$T` is not a struct type.")
    end

    types = fieldtypes(T)

    if !(hasfullconstructor(T))
        error("No valid constructor found for $type.")
    end

    T(read.(Ref(io), types)...)
end

"""
    write_struct(io::IO, sample::T) where T
Writes a instance of struct type `T` to `io`, field-by-field.

Optionally call `@reflect T` in your module to speed up (>10x) this method for type instances of type `T`.
"""
function write_struct(io::IO, sample::T) where T

    if !(isstructtype(T))
        error("`$T` is not a struct type.")
    end

    fields = fieldnames(T)

    sum(write.(Ref(io), getfield.(Ref(sample), fields)))
end

end # module


text:
using ParallelAccelerator

function hello()
    iterations = 100
    A = zeros(iterations)
    #initiiialize to value
    A = map(x->2,A)
    A = map(x->x+1,A)
    #end
    return A
end

function main()
    tic()
    B = hello()
    println("SELFTIMED ", toq())

    println(B[end])

end

main()

max_stars_repo_name: malekbr/julia-to-tiramisu

max_issues_repo_name: malekbr/julia-to-tiramisu

max_forks_repo_name: malekbr/julia-to-tiramisu


text:
#=##############################################################################
# DESCRIPTION
    Simulation interface connecting vehicles and maneuvers.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2019
  * License   : MIT
=###############################################################################

################################################################################
# SIMULATION TYPE
################################################################################
"""
    `Simulation(vehicle, maneuver, Vref::Real, RPMref::R, ttot:R)`

Simulation interface. This type carries the simulation's options and connects
vehicle and maneuver together.
"""
mutable struct Simulation{V<:AbstractVehicle, M<:AbstractManeuver, R<:Real}
    # USER INPUTS: Simulation setup
    vehicle::V              # Vehicle
    maneuver::M             # Maneuver to be performed
    Vref::R                 # Reference velocity in this maneuver
    RPMref::R               # Reference RPM in this maneuver
    ttot::R                 # Total time in which to perform the maneuver

    # OPTION USER INPUTS
    Vinit::Any              # Initial vehicle velocity
    Winit::Any              # Initial vehicle angular velocity

    # INTERNAL PROPERTIES: Runtime parameters
    t::R                    # Dimensional time
    nt::Int                 # Current time step number

    Simulation{V, M, R}(
                            vehicle, maneuver, Vref, RPMref, ttot;
                            Vinit=nothing, Winit=nothing,
                            t=zero(R), nt=-1
                        ) where {V, M, R} = _check(vehicle, maneuver) ? new(
                            vehicle, maneuver, Vref, RPMref, ttot,
                            Vinit, Winit,
                            t, nt
                        ) : nothing
end

# Implicit V and M constructor
Simulation(v::AbstractVehicle, m::AbstractManeuver, n::Real, args...; optargs...
             ) = Simulation{typeof(v), typeof(m), typeof(n)}(v, m, n, args...;
                                                                     optargs...)

# Solvers of concrete AbstractVehicle implementations
for header_name in ["unsteady", "quasisteady"]
  include("FLOWUnsteady_simulation_types_"*header_name*".jl")
end

##### FUNCTIONS  ###############################################################

"""
    `nextstep_kinematic(self::Simulation, dt::Real)`

Takes a kinematic time step `dt` where the new velocity and angular velocity
is calculated and the vehicle is translated and rotated according to it. It also
updates the tilt angle and RPM of every system.

NOTE: No solver is run in this process, rather it uses the current aerodynamic
solution to calculate acceleration and moment.
"""
function nextstep_kinematic(self::Simulation, dt::Real)

    if self.nt==-1 # Setup step (first step of simulation)

        # Initial vehicle linear and angular velocity
        if self.Vinit!=nothing; set_V(self.vehicle, self.Vinit); end;
        if self.Winit!=nothing; set_W(self.vehicle, self.Winit); end;

    else
        # Linear velocity increment
        #  dV = calc_dV(self.maneuver, self.vehicle, self.t, dt, self.ttot, self.Vref)
          dV = calc_dV(self.maneuver, dt, self.Vref, self.t, self.nt)

        # Angular velocity increment
        #dW = calc_dW(self.maneuver, self.vehicle, self.t, dt, self.ttot)
        dW = calc_dW(self.maneuver, dt, self.t, self.nt)

        # Update vehicle linear and angular velocity
        add_dV(self.vehicle, dV)
        add_dW(self.vehicle, dW)

        # Translates and rotates the vehicle according to current velocity
        nextstep_kinematic(self.vehicle, dt)

        # Rotate tilting systems
        angles = get_angles(self.maneuver, self.t/self.ttot)
        tilt_systems(self.vehicle, angles)

        self.t += dt
    end

    self.nt += 1
end

function save_vtk(self::Simulation, filename; optargs...)
    return save_vtk(self.vehicle, filename; num=self.nt, optargs...)
end

"""
Precalculations before calling the solver.
"""
function precalculations(self::Simulation, Vinf::Function,
                                pfield::vpm.ParticleField, t::Real, dt::Real)
    # Rotate rotors
    rotate_rotors(self, dt)

    # Calculate kinematic velocities
    precalculations(self.vehicle, Vinf, pfield, t, dt, self.nt)

    # # Bring the rotors back to their initial positions. This is needed to avoid
    # # double accounting for the RPM when the solver is run
    # NOTE: I've moved this to after shedding the wake and before running the
    # solver, otherwise the particles shedding wouldn't see the kinematic
    # velocity due to rotor RPM
    # rotate_rotors(self, -dt)

    return nothing
end

function rotate_rotors(self::Simulation{V, M, R}, dt::Real
                        ) where {V<:AbstractVLMVehicle, M<:AbstractManeuver, R}

    for (si, rotors) in enumerate(self.vehicle.rotor_systems)

        # RPM of this rotor system
        RPM = self.RPMref*get_RPM(self.maneuver, si, self.t/self.ttot)

        for rotor in rotors

            # Save solution fields
            sol1 = deepcopy(rotor._wingsystem.sol)
            sol2 = deepcopy(rotor.sol)
            Vinf = rotor._wingsystem.Vinf

            # Rotate rotor
            vlm.rotate(rotor, 360*RPM/60*dt)

            # Force regeneration of horseshoes
            if "Gamma" in keys(sol1)
                vlm.setVinf(rotor, Vinf)
                vlm.setRPM(rotor, RPM)
                vlm.getHorseshoe(rotor, 1)
            end

            # Paste solution fields back
            for (key,val) in sol1; vlm._addsolution(rotor._wingsystem, key, val); end;
            for (key,val) in sol2; rotor.sol[key] = val; end;

        end
    end
    return nothing
end


text:
# Objective #
# ========================== #
struct FunctionizeDecisionObjectiveBridge{T} <: MOIB.Objective.AbstractBridge end

function MOIB.Objective.bridge_objective(::Type{FunctionizeDecisionObjectiveBridge{T}}, model::MOI.ModelLike,
                                         f::SingleDecision) where T
    F = AffineDecisionFunction{T}
    MOI.set(model, MOI.ObjectiveFunction{F}(), convert(F, f))
    return FunctionizeDecisionObjectiveBridge{T}()
end

function MOIB.Objective.supports_objective_function(
    ::Type{<:FunctionizeDecisionObjectiveBridge}, ::Type{SingleDecision})
    return true
end
MOIB.added_constrained_variable_types(::Type{<:FunctionizeDecisionObjectiveBridge}) = Tuple{DataType}[]
function MOIB.added_constraint_types(::Type{<:FunctionizeDecisionObjectiveBridge})
    return Tuple{DataType, DataType}[]
end
function MOIB.set_objective_function_type(::Type{FunctionizeDecisionObjectiveBridge{T}}) where T
    return AffineDecisionFunction{T}
end

# Attributes, Bridge acting as a model
function MOI.get(bridge::FunctionizeDecisionObjectiveBridge, ::MOI.NumberOfVariables)
    return 0
end
function MOI.get(bridge::FunctionizeDecisionObjectiveBridge, ::MOI.ListOfVariableIndices)
    return MOI.VariableIndex[]
end

# No variables or constraints are created in this bridge so there is nothing to
# delete.
function MOI.delete(model::MOI.ModelLike, bridge::FunctionizeDecisionObjectiveBridge) end

function MOI.set(::MOI.ModelLike, ::MOI.ObjectiveSense,
                 ::FunctionizeDecisionObjectiveBridge, ::MOI.OptimizationSense)
    # `FunctionizeDecisionObjectiveBridge` is sense agnostic, therefore, we don't need to change
    # anything.
    return nothing
end
function MOI.get(model::MOI.ModelLike,
                 attr::MOIB.ObjectiveFunctionValue{SingleDecision},
                 bridge::FunctionizeDecisionObjectiveBridge{T}) where T
    F = AffineDecisionFunction{T}
    return MOI.get(model, MOIB.ObjectiveFunctionValue{F}(attr.result_index))
end
function MOI.get(model::MOI.ModelLike, attr::MOI.ObjectiveFunction{SingleDecision},
                 bridge::FunctionizeDecisionObjectiveBridge{T}) where T
    F = AffineDecisionFunction{T}
    func = MOI.get(model, MOI.ObjectiveFunction{F}())
    return convert(SingleDecision, func)
end


text:
using PtFEM

ProjDir = dirname(@__FILE__)

data = Dict(
  # Beam(ndim, nst, nxe, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Beam(2, 1, 100, 1, :x, Line(2, 1), false),
  :properties => [1.0 100.0;],
  :x_coords => 0.0:0.01:1.0,
  :support => [
    (1, [0 1]),
    #(51, [0 1]),
    (101, [0 1])
    ],
  :limit => 100,
  :tol => 0.00001
)

data |> display
println()

@time m = p46(data)
println()

if VERSION.minor > 5
  println("\nThe buckling load = $(m[1])")
  m[4][1] = 0.0
  println("\nThe buckling mode ($(m[2]) iterations):\n")
  for i in 1:m[3]
    println("$i    $(m[4][m[5][:, i]+1])")
  end
  println()
else
  println("\nThe buckling load = $(m[1]) ($(m[2]) iterations):\n")
  buckling_df = DataFrame(
    translation = m[4][m[5][1,:]+1],
    rotation = m[4][m[5][2,:]+1]
  )
  display(buckling_df)
  println()

  using Plots
  gr(size=(400,600))

  p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
  titles = ["p46.2a translation", "p46.2a rotation"]
  p[1] = plot(convert(Array, buckling_df[:translation]),
    ylim=(-0.1, 0.1), xlabel="node",
    ylabel="y translation [m]", color=:blue,
    marker=(:circle,1,0.1,stroke(1,:black)),
    title=titles[1], leg=false)
  p[2] = plot(convert(Array, buckling_df[:rotation]),
    ylim=(-1.0, 1.0), xlabel="node",
    ylabel="rotation [radians]", color=:red,
    marker=(:circle,1,0.1,stroke(1,:black)),
    title=titles[2], leg=false)

  plot(p..., layout=(2, 1))
  savefig(ProjDir*"/Ex46.2a.png")

end


text:
using MusicXML
using Test

cd(@__DIR__)

@testset "parsing example" begin
    include("../examples/parsing.jl")
    @test scorepartwise isa MX.ScorePartwise
    @test scoreparts isa Vector{MX.ScorePart}
    @test parts isa Vector{MX.Part}
end

@testset "creating example" begin
    include("../examples/creating.jl")
    @test isfile("../examples/myscore.musicxml")
end

@testset "grace note" begin
    include("../examples/grace.jl")
end


text:
using Hyperscript, Markdown, Test
using JSServe, Observables
using JSServe: Session, evaljs, linkjs, div
using JSServe: onjs, JSString, Dependency, jsrender
using JSServe: @js_str, js_type, pointer_identity, uuid, serialize_js, SerializationContext, serialize_binary
using JSServe.DOM
using JSServe.HTTP
using Electron
using ElectronDisplay
using URIParser
using Random
using Hyperscript: children
using JSServe.MsgPack
using JSServe.CodecZlib
using Test
using JSServe: jsrender

include("ElectronTests.jl")

function wait_on_test_observable()
    global test_observable
    test_channel = Channel{Dict{String, Any}}(1)
    f = on(test_observable) do value
        put!(test_channel, value)
    end
    val = take!(test_channel)
    off(test_observable, f)
    return val
end

"""
    test_value(app, statement)
Executes statemen (js code, or julia function with 0 args),
And waits on `test_observable` to push a new value!
Returns new value from `test_observable`
"""
function test_value(app, statement)
    # First start waiting on the test communication channel
    # We do this async before scheduling the js, since otherwise there is a
    # chance, that the event gets triggered before we have a chance to wait for it
    # which would make use wait forever
    val_t = @async wait_on_test_observable()
    # eval our js expression that is supposed to write something to test_observable
    if statement isa JSServe.JSCode
        JSServe.evaljs(app.session, statement)
    else
        statement()
    end
    fetch(val_t) # fetch the value!
end

JSServe.JSSERVE_CONFIGURATION.listen_port[] = 8555

@testset "JSServe" begin
    @testset "serialization" begin; include("serialization.jl"); end
    @testset "checkbox" begin; include("checkbox.jl"); end
    @testset "various" begin; include("various.jl"); end
    @testset "markdown" begin; include("markdown.jl"); end
    @testset "basics" begin; include("basics.jl"); end
end


text:
mutable struct Sample
	name::String                # Avoid utf8-characters that could make us create weird filenames causing problems later.
	fastq::Vector{String}       # List of fastq-files belonging to this sample.
	fastq2::Vector{String}      # Mates if paired-end, otherwise empty.
	fastqLocal::Vector{String}  # Local path to fastq-files. Downloaded from Synapse if needed.
	fastqLocal2::Vector{String} # Mates if paired-end, otherwise empty.
	referenceName::String       # Filename of reference genome.
	referencePath::String       # Path to reference genome.
	referencePathLocal::String  # Local path to reference genome. Downloaded from Synapse if needed.

	bam::String			     # Output BAM file (sorted).
	consensus::String		     # Output consensus file.
end
Sample(name,fastq::Vector{String},fastq2::Vector{String},fastqLocal::Vector{String},fastqLocal2::Vector{String}) =
	Sample(name,fastq,fastq2,fastqLocal,fastqLocal2,"","","","","")
Sample(name,fastq::Vector{String},fastq2::Vector{String}) = Sample(name,fastq,fastq2,String[],String[])
Sample(name,fastq::String,fastq2::String) = Sample(name,[fastq],[fastq2])
Sample(name,fastq::Vector{String}) = Sample(name,fastq,String[])
Sample(name,fastq::String) = Sample(name,[fastq])

# utility function for filling sample array from table (Matrix or DataFrame)
function find_samples(table, namePrefix::String="";
	                  name=nothing,
	                  fastq=nothing, fastq2=nothing,
	                  fastqLocal=nothing, fastqLocal2=nothing,
	                  separator=';',
	                  log=devnull)
	if namePrefix != "" && namePrefix[1] != '_'
		namePrefix = namePrefix * "_"
	end

	@assert name!=nothing  "name column must be specfified"
	@assert fastq!=nothing "fastq column must be specfified"

	samples = Sample[]

	for i=1:size(table,1)
		n = string(namePrefix, table[i,name])

		f1  = convert(Vector{String},split(table[i,fastq],separator,keep=false))
		f2  = fastq2 != nothing ? convert(Vector{String},split(table[i,fastq2],separator,keep=false)) : String[]
		f1L = fastqLocal  != nothing ? convert(Vector{String},split(table[i,fastqLocal], separator,keep=false)) : String[]
		f2L = fastqLocal2 != nothing ? convert(Vector{String},split(table[i,fastqLocal2],separator,keep=false)) : String[]

		push!(samples, Sample(n, f1, f2, f1L, f2L))
	end
	samples
end

# utility function that tries to identify sample names in folder were each sample might have multiple fastq files
function find_samples(syn, fastqPath::String, runName::String, namePrefix::String="", pattern::Regex=r".+(?=_L\d+_R1_\d+.fastq.gz$)"; log=devnull)
	if namePrefix != "" && namePrefix[1] != '_'
		namePrefix = namePrefix * "_"
	end

	sampleFolder = childpath(syn, fastqPath, runName)
	if isempty(sampleFolder)
		error("Could not find \"$runName\" in \"$fastqPath\".")
	end
	filePaths, fileNames = listfiles(syn, sampleFolder) # filePaths are local paths or Synapse IDs

	# find files matching pattern
	# matches = map( x->match(pattern,x), fileNames )
	# mask = falses(matches)
	# map!( x->x!=nothing, mask, matches );
	matches = match.(pattern,fileNames)
	mask = matches.!=nothing

	# log warning for non-matching files
	for f in fileNames[.~mask]
		printifwarning(log, "File \"$f\" was ignored.")
	end

	# remove files that do not match pattern
	fileNames = fileNames[mask]
	filePaths = filePaths[mask]
	matches = convert(Vector{RegexMatch}, matches[mask])

	# extract the matching part of the regex as strings
	matchingNames = Vector{String}(undef,size(matches))
	map!( x->x.match, matchingNames, matches )

	# put all files with the same match together (and add the run name to the sample name)
	samples = [ Sample("$namePrefix$u", filePaths[matchingNames.==u]) for u in unique(matchingNames) ]

	for s in samples
		n = length(s.fastq)
		println(log, "Found sample \"$(s.name)\" with $n fastq files.")
	end
	samples
end

# find paths (Synapse or local) and sample IDs (i.e. base file names) of aligned samples (.bam)
function find_aligned(syn, path::AbstractString, runName::AbstractString)
	paths, names = listaligned(syn, path, runName)
	names = map(x->x[1:end-4], names) # keep only matches and remove ".bam" ending
	paths, names
end

# To avoid extending Base.ismatch
_ismatch(r::Regex, s::AbstractString) = match(r,s)!=nothing
_ismatch(f::Function, s::AbstractString) = f(s)

# each sample name should match exactly one of the patterns ref[i][1] ∀i=1...n and will be assigned the corresponding reference name ref[i][2]
function assign_reference!(samples::Array{Sample,1}, refs::Vector; log=devnull)
	patterns = [r[1] for r in refs]
	refNames = [r[2] for r in refs]

	M = falses(length(refs))
	for (i,s) in enumerate(samples)
		broadcast!(_ismatch,M,patterns,[s.name])
		nbrMatches = sum(M)

		@assert nbrMatches>=1 "Sample \"$(s.name)\" did not match any reference"
		@assert nbrMatches<=1 "Sample \"$(s.name)\" matches references $(refNames[M])"

		index = findfirst(M)
		s.referenceName      = refNames[index]
		s.referencePath      = refs[index][3]
		s.referencePathLocal = refs[index][4]
		println(log,"Assigned reference \"$(s.referenceName)\" to sample \"$(s.name)\".")
	end
end


text:
# Julia wrapper for header: /usr/include/libavutil/aes_ctr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_aes_ctr_alloc,
    av_aes_ctr_init,
    av_aes_ctr_free,
    av_aes_ctr_crypt,
    av_aes_ctr_get_iv,
    av_aes_ctr_set_random_iv,
    av_aes_ctr_set_iv,
    av_aes_ctr_set_full_iv,
    av_aes_ctr_increment_iv

function av_aes_ctr_alloc()
    ccall((:av_aes_ctr_alloc, libavutil), Ptr{AVAESCTR}, ())
end

function av_aes_ctr_init(a, key)
    ccall((:av_aes_ctr_init, libavutil), Cint, (Ptr{AVAESCTR}, Ptr{UInt8}), a, key)
end

function av_aes_ctr_free(a)
    ccall((:av_aes_ctr_free, libavutil), Cvoid, (Ptr{AVAESCTR},), a)
end

function av_aes_ctr_crypt(a, dst, src, size::Integer)
    ccall((:av_aes_ctr_crypt, libavutil), Cvoid, (Ptr{AVAESCTR}, Ptr{UInt8}, Ptr{UInt8}, Cint), a, dst, src, size)
end

function av_aes_ctr_get_iv(a)
    ccall((:av_aes_ctr_get_iv, libavutil), Ptr{UInt8}, (Ptr{AVAESCTR},), a)
end

function av_aes_ctr_set_random_iv(a)
    ccall((:av_aes_ctr_set_random_iv, libavutil), Cvoid, (Ptr{AVAESCTR},), a)
end

function av_aes_ctr_set_iv(a, iv)
    ccall((:av_aes_ctr_set_iv, libavutil), Cvoid, (Ptr{AVAESCTR}, Ptr{UInt8}), a, iv)
end

function av_aes_ctr_set_full_iv(a, iv)
    ccall((:av_aes_ctr_set_full_iv, libavutil), Cvoid, (Ptr{AVAESCTR}, Ptr{UInt8}), a, iv)
end

function av_aes_ctr_increment_iv(a)
    ccall((:av_aes_ctr_increment_iv, libavutil), Cvoid, (Ptr{AVAESCTR},), a)
end


text:
# WARNING: I still have doubts this works correctly, use with caution

mutable struct BND
    ptr::Ptr{Void}
end

function BND(xtd::TD, mode::Cuint)
    d = Ref{Cptr}(0)
    @cuda(cudnn, cudnnCreateTensorDescriptor, (Ptr{Cptr},), d)
    @cuda(cudnn, cudnnDeriveBNTensorDescriptor,
          (Ref{Void}, Cptr, Cuint),
          d[], xtd.ptr, mode)
    return BND(d[])
end

function batchnorm_param_size(x::CuArray, mode::UInt32)
    xsz = size(x)
    if mode == CUDNN_BATCHNORM_PER_ACTIVATION
        psz = (1, xsz[3], xsz[2], xsz[1])  # 1xCxHxW
    elseif mode == CUDNN_BATCHNORM_SPATIAL
        psz = (1, xsz[3], 1, 1)            # 1xCx1x1
    else
        error("Mode $mode is not supported")
    end
    return psz
end

mutable struct BatchNormState{T}
    bnScale::CuArray{T,4}
    bnBias::CuArray{T,4}
    resultRunningMean::CuArray{T,4}
    resultRunningVariance::CuArray{T,4}
    resultSaveMean::CuArray{T,4}
    resultSaveInvVariance::CuArray{T,4}
end

function BatchNormState(x::CuArray{T,4}, mode=CUDNN_BATCHNORM_SPATIAL) where T
    psz = batchnorm_param_size(x, mode)
    bnScale = CuArray(randn(T, psz))
    bnBias = CuArray(randn(T, psz))
    resultRunningMean = CuArray(zeros(T, psz))
    resultRunningVariance = CuArray(zeros(T, psz))
    resultSaveMean = CuArray(zeros(T, psz))
    resultSaveInvVariance = CuArray(zeros(T, psz))
    return BatchNormState(bnScale, bnBias, resultRunningMean, resultRunningVariance,
                          resultSaveMean, resultSaveInvVariance)
end

function batchnorm_train!(y::CuArray{T,4}, x::CuArray{T,4}, s::BatchNormState;
                          handle=cudnnhandle(), alpha=1, beta=0,
                          exponentialAverageFactor=T(1), mode=CUDNN_BATCHNORM_SPATIAL,
                          epsilon=CUDNN_BN_MIN_EPSILON) where T
    xtd = TD(x)
    ytd = TD(y)
    bnScaleBiasMeanVarDesc = BND(xtd, mode)
    @cuda(cudnn, cudnnBatchNormalizationForwardTraining,
          (Cptr, UInt32, Cptr, Cptr, Cptr, Cptr, Cptr, Cptr,
           Cptr, Cptr, Cptr, Cdouble,
           Cptr, Cptr, Cdouble,
           Cptr, Cptr),
          handle, mode, Ref(T(alpha)), Ref(T(beta)), xtd.ptr, x.ptr, ytd.ptr, y.ptr,
          bnScaleBiasMeanVarDesc.ptr, s.bnScale.ptr, s.bnBias.ptr, exponentialAverageFactor,
          s.resultRunningMean.ptr, s.resultRunningVariance.ptr, epsilon,
          s.resultSaveMean.ptr, s.resultSaveInvVariance.ptr)

end

function batchnorm_train(x::CuArray{T,4}, s::BatchNormState; opts...) where T
    y = similar(x)
    batchnorm_train!(y, x, s;  opts...)
    return y
end

function batchnorm_infer!(y::CuArray{T,4}, x::CuArray{T,4}, s::BatchNormState;
                          handle=cudnnhandle(), alpha=1, beta=0,
                          exponentialAverageFactor=T(1), mode=CUDNN_BATCHNORM_SPATIAL,
                          epsilon=CUDNN_BN_MIN_EPSILON) where T
    xtd = TD(x)
    ytd = TD(y)
    xsz = size(x)
    bnScaleBiasMeanVarDesc = BND(xtd, mode)
    estimatedMean = s.resultRunningMean
    estimatedVariance = s.resultRunningVariance
    @cuda(cudnn, cudnnBatchNormalizationForwardInference,
          (Cptr, UInt32, Cptr, Cptr, Cptr, Cptr, Cptr, Cptr,
           Cptr, Cptr, Cptr,
           Cptr, Cptr, Cdouble),
          handle, mode, Ref(T(alpha)), Ref(T(beta)), xtd.ptr, x.ptr, ytd.ptr, y.ptr,
          bnScaleBiasMeanVarDesc.ptr, s.bnScale.ptr, s.bnBias.ptr,
          estimatedMean.ptr, estimatedVariance.ptr, epsilon)
end

function batchnorm_infer(x::CuArray{T,4}, s::BatchNormState; opts...) where T
    y = similar(x)
    batchnorm_infer!(y, x, s; opts...)
    return y
end

function batchnorm_grad!(dx::CuArray{T,4}, x::CuArray{T,4}, dy::CuArray{T,4}, s::BatchNormState;
                         handle=cudnnhandle(), alpha_data=1, beta_data=0,
                         alpha_param=1, beta_param=0,
                         exponentialAverageFactor=T(1), mode=CUDNN_BATCHNORM_SPATIAL,
                         epsilon=CUDNN_BN_MIN_EPSILON) where T
    xtd = TD(x)
    dytd = TD(dy)
    dxtd = TD(dx)
    xsz = size(x)
    bnScaleBiasMeanVarDesc = BND(xtd, mode)
    # should we update bnScale & bnBias manually or cuDNN does it automatically?
    resultBnScaleDiff = similar(s.bnScale)
    resultBnBiasDiff = similar(s.bnScale)
    savedMean = s.resultSaveMean
    savedInvVariance = s.resultSaveInvVariance
    @cuda(cudnn, cudnnBatchNormalizationBackward,
          (Cptr, UInt32, Cptr, Cptr, Cptr, Cptr,
           Cptr, Cptr, Cptr, Cptr, Cptr, Cptr,
           Cptr, Cptr, Cptr, Cptr,
           Cdouble, Cptr, Cptr),
          handle, mode, Ref(T(alpha_data)), Ref(T(beta_data)), Ref(T(alpha_param)), Ref(T(beta_param)),
          xtd.ptr, x.ptr, dytd.ptr, dy.ptr, dxtd.ptr, dx.ptr,
          bnScaleBiasMeanVarDesc.ptr, s.bnScale.ptr, resultBnScaleDiff.ptr, resultBnBiasDiff.ptr,
          epsilon, savedMean.ptr, savedInvVariance.ptr)
end

function batchnorm_grad(x::CuArray{T,4}, dy::CuArray{T,4}, s::BatchNormState; opts...) where T
    dx = similar(x)
    batchnorm_grad!(dx, x, dy, s; opts...)
    return dx
end


text:
# Generates a random cone problem to test the MPB interface perf

rows = 10000
cols = 10000

srand(1988)
A = sprand(rows, cols, 0.01)
b = rand(rows)
c = rand(cols)

coneset = [:NonNeg, :NonPos, :SOC]
var_cones = {}
vars_coned = 0
while true
    (vars_coned >= cols) && break
    if vars_coned == cols - 1
        push!(var_cones, (:NonNeg, cols))
        break
    end
    cone = coneset[rand(1:3)]
    n_var = rand(2:5)
    end_var = min(vars_coned+n_var, cols)
    push!(var_cones, (cone, vars_coned+1:end_var) )
    vars_coned += n_var
end

coneset = [:Zero, :NonNeg, :NonPos, :SOC]
row_cones = {}
rows_coned = 0
while true
    (rows_coned >= rows) && break
    if rows_coned == rows - 1
        push!(row_cones, (:NonNeg, rows))
        break
    end
    cone = coneset[rand(1:4)]
    n_row = rand(2:5)
    end_row = min(rows_coned+n_row, rows)
    push!(row_cones, (cone, rows_coned+1:end_row) )
    rows_coned += n_row
end

using ECOS
m = ECOS.ECOSMathProgModel()

println("First run")
@time    ECOS.loadconicproblem!(m, c, A, b, row_cones, var_cones)
Base.gc()

println("Second run")
@profile ECOS.loadconicproblem!(m, c, A, b, row_cones, var_cones)
Profile.print(format=:flat)
Base.gc()

println("Third run")
@time    ECOS.loadconicproblem!(m, c, A, b, row_cones, var_cones)


text:
export Shamanskii_BFGS

function Shamanskii_BFGS(nlp_stop,
					     PData  :: PDataFact,
					     dₙ 	:: Vector,
					     gt 	:: Vector)
	# printstyled("on est ici ↓ \n", color = :bold)
	T = eltype(gt)
	nlp_at_x = nlp_stop.current_state
	∇fₖ = nlp_at_x.gx
	x = nlp_at_x.x
	# Bₖ = Symmetric(hess(nlp_stop.pb, x), :L)
	# @show nlp_at_x.Hx
	Bₖ = nlp_at_x.Hx
	# printstyled("on a Bₖ \n", color = :bold)

	sₖ = dₙ
	# printstyled("on a sₖ \n", color = :bold)
	xtemp = x .+ sₖ
	# printstyled("on a xtemp \n", color = :bold)
	∇fₖ₊₁ = grad(nlp_stop.pb, xtemp)
	# printstyled("on a fk+1 \n", color = :bold)
	yₖ = ∇fₖ₊₁ .- ∇fₖ

	Yₖ = yₖ * yₖ'
	yts = yₖ' * sₖ
	Yₖ = Yₖ ./ yts
	# printstyled("on a Yₖ \n", color = :bold)

	Bs = Bₖ * sₖ
	BssBₖ = Bs * Bs'
	sBs = Bs' * sₖ
	BssBₖ = (1.0 ./ sBs) .* BssBₖ

	Bₖ₊₁ = Bₖ .+ Yₖ .- BssBₖ
	# LDLTB = ldl(Bₖ₊₁)
	L, D, pp, ρ, ncomp = ldlt_symm(Bₖ₊₁, 'r')
    # printstyled("on a L, D,... \n", color = :bold)
	D = Hermitian(D)
    X = eigen(D)
	# printstyled("on a X \n", color = :bold)
    Δ = X.values
	Q =  X.vectors
    l_m, = findmin(Δ)
	# printstyled("on a l_m \n", color = :bold)
    ĝ = L \ (∇fₖ₊₁[pp])
	# printstyled("on a g chapeau \n", color = :bold)
    g̃ = Q' * ĝ
	# printstyled("on a g tilde \n", color = :bold)
    n_g = norm(∇fₖ₊₁)
    λ =  max(-l_m, 0.0)
	# printstyled("on a λ \n", color = :bold)
	ϵ2 = sqrt(eps(T))
	Γ = max.(abs.(Δ), ϵ2)

	d̃ = L \ ∇fₖ₊₁[pp]
   	d̂ = L' \ (Q * (Q' * d̃ ./ Γ))
   	dₛ = -d̂[invperm(pp)]
	# xt = xtemp + dₛ

	dₕₒ = dₙ .+ dₛ

	return dₕₒ
end


text:
struct ZulipClient
    ep::String
    headers::Vector{Pair{String, String}}
end

const ZulipGlobal = Ref(ZulipClient("", []))

function ZulipClient(; email = "", apikey = "", ep = "https://julialang.zulipchat.com", api_version = "v1")
    if isempty(email) || isempty(apikey)
        throw(ArgumentError("Arguments email and apikey should not be empty."))
    end

    key = base64encode(email * ":" * apikey)
    headers = ["Authorization" => "Basic " * key, "Content-Type" => "application/x-www-form-urlencoded"]
    endpoint = ep * "/api/" * api_version * "/"

    client = ZulipClient(endpoint, headers)

    return client
end

function global_zulip!(; email = "", apikey = "", ep = "https://julialang.zulipchat.com", api_version = "v1")
    client = ZulipClient(; email = email, apikey = apikey, ep = ep, api_version = api_version)
    ZulipGlobal[] = client
end

function query(client::ZulipClient, apimethod, params; method = "POST")
    params = HTTP.URIs.escapeuri(params)
    url = client.ep * apimethod
    JSON3.read(HTTP.request(method, url, client.headers, params).body)
end

function sendMessage(client::ZulipClient = ZulipGlobal[]; params...)
    query(client, "messages", Dict(params))
end

updateMessage(msg_id; params...) = updateMessage(ZulipGlobal[], msg_id; params...)
function updateMessage(client::ZulipClient, msg_id; params...)
    query(client, "messages/" * string(msg_id), Dict(params), method = "PATCH")
end

deleteMessage(msg_id; params...) = deleteMessage(ZulipGlobal[], msg_id; params...)
function deleteMessage(client::ZulipClient, msg_id; params...)
    query(client, "messages/" * string(msg_id), Dict(params), method = "DELETE")
end


text:
@diffrule NNlib.conv2d(x, w) w NNlib.conv2d_grad_w(x, w, ds;)
@diffrule NNlib.conv2d(x, w) x NNlib.conv2d_grad_x(x, w, ds)

@diffrule NNlib.conv2d(x, w; _opts...) w NNlib.conv2d_grad_w(x, w, ds; _opts...)
@diffrule NNlib.conv2d(x, w; _opts...) x NNlib.conv2d_grad_x(x, w, ds; _opts...)

@diffrule NNlib.pool(x) x NNlib.pool_grad(x, pool(x), ds)
@diffrule NNlib.pool(x; _opts...) x NNlib.pool_grad(x, pool(x), ds; _opts...)

@diffrule NNlib.σ(x::Number) x (NNlib.σ(x) .* (1 .- NNlib.σ(x)) .* ds)
@diffrule NNlib.sigmoid(x::Number) x (NNlib.σ(x) .* (1 .- NNlib.σ(x)) .* ds)

@diffrule NNlib.relu(x::Number) x (x .> 0) .* ds

# logistic(x) = 1 ./ (1 + exp.(-x))
# @diffrule logistic(x::Number) x (logistic(x) .* (1 .- logistic(x)) .* ds)

# softplus(x) = log(exp(x) + 1)
# @diffrule softplus(x::Number) x logistic(x) .* ds

@diffrule CUDAnative.log(x::Real) x ds / x
@diffrule CUDAnative.log(x::AbstractArray) x ds ./ x