text:
using ITensors
using Test

# Wrap an ITensor with pairs of primed and
# unprimed indices to pass to davidson
struct ITensorMap
  A::ITensor
end
Base.eltype(M::ITensorMap) = eltype(M.A)
Base.size(M::ITensorMap) = dim(IndexSet(inds(M.A;plev=0)...))
(M::ITensorMap)(v::ITensor) = noprime(M.A*v)

@testset "Complex davidson" begin
  d = 10
  i = Index(d,"i")
  A = randomITensor(ComplexF64,i,prime(i))
  A = mapprime(A*mapprime(dag(A),0,2),2,1)
  M = ITensorMap(A)

  v = randomITensor(i)
  λ,v = davidson(M,v;maxiter=10)
  @test M(v)≈λ*v

  v = randomITensor(ComplexF64, i)
  λ,v = davidson(M,v;maxiter=10)
  @test M(v)≈λ*v

end

nothing


text:
"""
    PolyharmonicSpline <: Kernel

TODO
"""
struct PolyharmonicSpline <: Kernel
    # TODO
end

ext: jl


text:
using Test, Avro, UUIDs, Dates, StructTypes, JSON3, Tables, SentinelArrays

@testset "Avro.jl" begin

# missing
buf = Avro.write(missing)
@test isempty(buf)

# booleans
buf = Avro.write(true)
@test buf[1] == 0x01
buf = Avro.write(false)
@test buf[1] == 0x00

# integers
buf = Avro.write(1)
@test buf[1] == 0x02
buf = Avro.write(63)
@test buf[1] == 0x7e
buf = Avro.write(64)
@test buf == [0x80, 0x01]

x = typemax(UInt8)
@test Avro.read(Avro.write(x), UInt8) === x

buf = Avro.write(-1)
@test buf[1] == 0x01
buf = Avro.write(-63)
@test buf[1] == 0x7d
buf = Avro.write(-64)
@test buf[1] == 0x7f
buf = Avro.write(-65)
@test buf == [0x81, 0x01]

for i = typemin(Int16):typemax(Int16)
    @test i == Avro.read(Avro.write(i), Int)
end

# floats
for x in (-0.0001, 0.0, -0.0, 1.0, floatmin(Float32), floatmax(Float32), floatmin(Float64), floatmax(Float64))
    @test x === Avro.read(Avro.write(x), typeof(x))
end

# bytes
x = [UInt8(y) for y in "hey there stranger"]
@test x == Avro.read(Avro.write(x), typeof(x))

x = UInt8[]
@test x == Avro.read(Avro.write(x), typeof(x))

# strings
x = "hey there stranger"
@test x == Avro.read(Avro.write(x), typeof(x))

x = ""
@test x == Avro.read(Avro.write(x), typeof(x))

# array
x = [1, 2, 3, 4, 5]
@test x == Avro.read(Avro.write(x), typeof(x))

# array of strings
x = ["hey", "there", "stranger"]
@test x == Avro.read(Avro.write(x), typeof(x))

x = String[]
@test x == Avro.read(Avro.write(x), typeof(x))

# fixed
x = tuple(b"hey"...)
@test x == Avro.read(Avro.write(x), typeof(x))

x = ()
@test x == Avro.read(Avro.write(x), typeof(x))

# maps
x = Dict("hey" => 1, "there" => 2, "stranger" => 3)
@test x == Avro.read(Avro.write(x), typeof(x))

x = Dict{String, Int}()
@test x == Avro.read(Avro.write(x), typeof(x))

# enums
x = Avro.Enum{(:hey, :there, :stranger)}(0)
@test x == Avro.read(Avro.write(x), typeof(x))

# unions
x = 1
@test x == Avro.read(Avro.write(x; schema=Union{Int, String}), Union{Int, String})

# records
x = (a=1, b=3.4, c="hey")
@test x == Avro.read(Avro.write(x), typeof(x))

r = Avro.read(Avro.write(x), Avro.Record{(:a, :b, :c), Tuple{Int, Float64, String}})
@test r.a == 1
@test r.b == 3.4
@test r.c == "hey"

struct Person
    id::Int
    name::String
end

StructTypes.StructType(::Type{Person}) = StructTypes.Struct()

x = Person(10, "Valentin")
@test x == Avro.read(Avro.write(x), typeof(x))

x = [Person(1, "meg"), Person(2, "jo"), Person(3, "beth"), Person(4, "amy")]
@test x == Avro.read(Avro.write(x), typeof(x))

# logical
x = Avro.Decimal{0, 4}(Int128(1))
@test x == Avro.read(Avro.write(x), typeof(x))

x = UUID(rand(UInt128))
@test x == Avro.read(Avro.write(x), typeof(x))

x = Dates.today()
@test x == Avro.read(Avro.write(x), typeof(x))

x = Time(Dates.now())
@test x == Avro.read(Avro.write(x), typeof(x))

x = Dates.now()
@test x == Avro.read(Avro.write(x), typeof(x))

x = Avro.Duration(1, 2, 3)
@test x == Avro.read(Avro.write(x), typeof(x))

# combinations
cases = [
    # arrays
    [missing, missing, missing],
    [true, false, true],
    [1.2, 3.4, 5.6],
    [Vector{UInt8}("hey"), Vector{UInt8}("there"), Vector{UInt8}("stranger")],
    [Avro.Enum{(:a, :b, :c)}(0), Avro.Enum{(:a, :b, :c)}(1), Avro.Enum{(:a, :b, :c)}(2)],
    [[1, 2], [3, 4, 5], [6, 7, 8, 9]],
    [Dict(:a => Float32(1)), Dict(:b => Float32(2)), Dict(:c => Float32(3))],
    [(a=Date(2021, 1, 1), b=true), (a=Date(2021, 1, 2), b=false), (a=Date(2021, 1, 3), b=true)],
    Union{Missing, UUID, Avro.Duration}[UUID(rand(UInt128)), missing, Avro.Duration(4, 5, 6)],
    Union{Missing, Int32, Vector{UInt8}, Date, UUID, Dict{String, NamedTuple{(:a,), Tuple{Union{Int64, Float32}}}}, Vector{Union{NamedTuple{(:a,), Tuple{Int64}}, Avro.Enum{(:a, :b)}}}}[missing, Int32(4), Vector{UInt8}("hey"), Date(2021, 2, 1), UUID(rand(UInt128)), Dict{String, NamedTuple{(:a,), Tuple{Union{Int64, Float32}}}}("a" => (a=Int64(1),), "b" => (a=Float32(3.14),)), Union{NamedTuple{(:a,), Tuple{Int64}}, Avro.Enum{(:a, :b)}}[(a=1001,), Avro.Enum{(:a, :b)}(0), Avro.Enum{(:a, :b)}(1)]],
    # maps
    Dict("a" => missing),
    Dict("a" => true, "b" => false),
    Dict("a" => 1.2, "b" => 3.4),
    Dict("a" => Vector{UInt8}("hey"), "b" => Vector{UInt8}("there")),
    Dict("a" => Avro.Enum{(:a, :b, :c)}(0), "b" => Avro.Enum{(:a, :b, :c)}(1)),
    Dict("a" => [1, 2], "b" => [3, 4, 5]),
    Dict("a" => Dict(:a => Float32(1)), "b" => Dict(:b => Float32(2))),
    Dict("a" => (a=Time(1, 2, 3), b=UUID(rand(UInt128))), "b" => (a=Time(4, 5, 6), b=UUID(rand(UInt128)))),
    Dict{String, Union{Missing, DateTime, Avro.Decimal{1, 4}}}("a" => missing, "b" => Dates.now(), "c" => Avro.Decimal{1, 4}(12345)),
    # records
    (a=missing,),
    (a=[missing, missing],),
    (a=Dict("a" => missing),),
    (a=Dict("a" => [missing, missing]),),
    (a=true, b=false),
    (a=[true, false], b=[false, false]),
    (a=Dict("a" => true),),
    (a=Dict("a" => [true, false], "b" => [false, false]),),
    (a=1.2,),
    (a=[1.2, 3.4],),
    (a=Dict("a" => 4.5),),
    (a=Dict("a" => [6.7, 8.9]),),
    (a=Vector{UInt8}("hey"), b=Avro.Enum{(:a, :b)}(0), c=[1, 2]),
    (a=Dict('a' => Float32(1.2)),),
    (a=(a=Date(2021, 1, 1), b=Time(1, 2, 3), c=UUID(rand(UInt128))),),
    NamedTuple{(:a,), Tuple{Union{Missing, UUID, Dict{String, Float32}}}}((Dict("a" => Float32(1.4)),)),
]

for case in cases
    @test isequal(case, Avro.read(Avro.write(case), typeof(case)))
    # write out schema, object to file from julia
    # sch = Avro.schematype(typeof(case))
    # JSON3.write("schema.avsc", sch)
    # Avro.write("x.avro", case)
    # js = """
    #     let avro = require('avro-js');
    #     let fs = require('fs');
    #     // read schema, object from file in node
    #     let sch = avro.parse('./schema.avsc');
    #     let x = sch.fromBuffer(fs.readFileSync('./x.avro'));
    #     // write out schema, object to file from node
    #     fs.writeFileSync('jsschema.avsc', sch.getSchema());
    #     fs.writeFileSync('jsx.avro', sch.toBuffer(x));
    # """
    # run(`node -e $js`)
    # # read in schema, object from file in julia
    # sch2 = Avro.parseschema("jsschema.avsc"))
    # case2 = Avro.read(Base.read("jsx.avro"), sch)
    # @test sch == sch2
    # @test isequal(case, case2)
end

nt = (a=1, b=2, c=3)
rt = [nt, nt, nt]
for comp in (:deflate, :bzip2, :xz, :zstd)
    io = Avro.tobuffer(rt; compress=comp)
    tbl = Avro.readtable(io)
    @test length(tbl) == 3
    @test tbl[1].a == nt.a
    @test tbl[1].b == nt.b
    @test tbl[1].c == nt.c
end

nt = (a=[1, 2, 3], b=[4.0, 5.0, 6.0], c=["7", "8", "9"])
io = Avro.tobuffer(nt)
tbl = Avro.readtable(io)
@test length(tbl) == 3
@test tbl.sch == Tables.Schema((:a, :b, :c), (Int, Float64, String))

rt = [(a=1, b=4.0, c="7"), (a=2.0, b=missing, c="8"), (a=3, b=6.0, c="9")]
io = Avro.tobuffer(rt)
tbl = Avro.readtable(io)
@test length(tbl) == 3

rt = [
    (a=1, b=2, c=3),
    (b=4.0, c=missing, d=5),
    (a=6, d=7),
    (a=8, b=9, c=10, d=missing),
    (d=11, c=10, b=9, a=8)
]

dct = Tables.dictcolumntable(rt)
io = Avro.tobuffer(dct)
tbl = Avro.readtable(io)
@test length(tbl) == 5

end

@testset "rowwise" begin
    rows = [
        (a=missing, b=2., c=3, d=4),
        (a=1, b=4.0, c=missing, d=5),
        (a=1, b=4.0, c=missing, d="5"),
    ]
    sch = Tables.Schema([:a, :b, :c, :d], [Union{Missing, Int64}, Float64, Union{Missing, Int64}, Union{Int64, String}])

    mktempdir() do p
        pth = joinpath(p,"tmp.avro")
        Avro.recordwriter(sch, pth) do rw
            Avro.writerecord(rw, rows[1])
            Avro.writerecord(rw, rows[2])
            Avro.writerecord(rw, rows[3])
        end

        Avro.recordreader(pth) do rr
            @test all(collect(skipmissing(first(rr))) .== collect(skipmissing(rows[1])))
            @test all(collect(skipmissing(first(rr))) .== collect(skipmissing(rows[1])))
            res= []
            for r in rr
                push!(res, NamedTuple(r))
            end
            @test all(collect.(skipmissing.(res)).==collect.(skipmissing.(rows)))
        end
        tbl = Avro.readtable(pth)
        @test length(tbl) == 3
        @test tbl.sch == sch

        pth = joinpath(p,"tmp.avro")
        Avro.recordwriter(sch, pth; compress=:bzip2) do rw
            Avro.writerecord(rw, rows[1])
        end
        tbl = Avro.readtable(pth)
        f = first(tbl)
        @test f.b == 2.
    end
end

# using CSV, Dates, Tables, Test
# const dir = joinpath(dirname(pathof(CSV)), "..", "test", "testfiles")
# include(joinpath(dirname(pathof(CSV)), "..", "test", "testfiles.jl"))
# for (i, test) in enumerate(testfiles)
#     file, kwargs, expected_sz, expected_sch, testfunc = test
#     println("testing $file, i = $i")
#     f = CSV.File(file isa IO ? file : joinpath(dir, file); kwargs...)
#     buf = Avro.tobuffer(f)
#     tbl = Avro.readtable(buf)
#     @test isequal(columntable(f), columntable(tbl))
# end


text:
import BenchPerfConfigSweeps
import DisplayAs
using DataFrames
using VegaLite

include("utils.jl")

sweepresult = BenchPerfConfigSweeps.load(joinpath(@__DIR__, "build"))

df_raw = DataFrame(sweepresult)

let
    global df = select(df_raw, Not(:trial))
    df[:, :time_ns] = map(t -> minimum(t.benchmark).time, df_raw.trial)
    # df[:, :evals] = map(t -> t.benchmark.params.evals, df_raw.trial)
    df[:, :L1_miss_percent] = map(df_raw.trial) do t
        t.perf["L1-dcache-load-misses"] / t.perf["L1-dcache-loads"] * 100
    end
    df[:, :LL_miss_percent] = map(df_raw.trial) do t
        get(t.perf, "LLC-load-misses", missing) / get(t.perf, "LLC-loads", missing) * 100
    end
    nrepeats = nrepeats_from_n.(df.n)
    bytes = df.n .* 2 * sizeof(Int)
    global throughput_unit = "GiB/thread/sec"
    df[:, :throughput] = bytes .* nrepeats ./ 2^30 ./ (df.time_ns / 1e9)
    df[:, :MiB] = bytes ./ 2^20
    df[:, :total_MiB] = df.nthreads .* bytes ./ 2^20
    df
end

resultdir = joinpath(@__DIR__, "result")
mkpath(resultdir)

saveresult(; plots...) = saveresult(:png, :svg; plots...)
function saveresult(exts::Symbol...; plots...)
    for (k, v) in plots
        for e in exts
            save(joinpath(resultdir, "$k.$e"), v)
        end
    end
end
nothing  # hide

plt_throughput_cache_miss_vs_different_input_sizes =
    df |> [
        @vlplot(
            mark = {:line, point = true},
            x = {:MiB, scale = {type = :log}},
            y = {:throughput, title = throughput_unit},
            color = {:nthreads, type = :nominal},
        )
        @vlplot(
            mark = {:line, point = true},
            x = {:total_MiB, scale = {type = :log}},
            y = {:throughput, title = throughput_unit},
            color = {:nthreads, type = :nominal},
        )
    ]
saveresult(; plt_throughput_cache_miss_vs_different_input_sizes)
plt_throughput_cache_miss_vs_different_input_sizes
plt_throughput_cache_miss_vs_different_input_sizes |> DisplayAs.PNG

plt_throughput_cache_miss_vs_per_thread_size =
    df |> [
        @vlplot(
            mark = {:line, point = true},
            x = {:MiB, scale = {type = :log}},
            y = {:throughput, title = throughput_unit},
            color = {:nthreads, type = :nominal},
        )
        @vlplot(
            mark = {:line, point = true},
            x = {:MiB, scale = {type = :log}},
            y = {:L1_miss_percent, title = "L1 cache miss [%]"},
            color = {:nthreads, type = :nominal},
        )
        @vlplot(
            mark = {:line, point = true},
            x = {:MiB, scale = {type = :log}},
            y = {:LL_miss_percent, title = "LL cache miss [%]"},
            color = {:nthreads, type = :nominal},
        )
    ]
saveresult(; plt_throughput_cache_miss_vs_per_thread_size)
plt_throughput_cache_miss_vs_per_thread_size
plt_throughput_cache_miss_vs_per_thread_size |> DisplayAs.PNG
#-

plt_throughput_cache_miss_vs_total_size =
    df |> [
        @vlplot(
            mark = {:line, point = true},
            x = {:total_MiB, scale = {type = :log}},
            y = {:throughput, title = throughput_unit},
            color = {:nthreads, type = :nominal},
        )
        @vlplot(
            mark = {:line, point = true},
            x = {:total_MiB, scale = {type = :log}},
            y = {:L1_miss_percent, title = "L1 cache miss [%]"},
            color = {:nthreads, type = :nominal},
        )
        @vlplot(
            mark = {:line, point = true},
            x = {:total_MiB, scale = {type = :log}},
            y = {:LL_miss_percent, title = "LL cache miss [%]"},
            color = {:nthreads, type = :nominal},
        )
    ]
saveresult(; plt_throughput_cache_miss_vs_total_size)
plt_throughput_cache_miss_vs_total_size
plt_throughput_cache_miss_vs_total_size |> DisplayAs.PNG
#-

plt_throughput_cache_miss_vs_total_size  #src


text:
function download_pkg(pkg)
    Pkg.develop(PackageSpec(path="https://github.com/korbinian90/$pkg.jl"))
    Pkg.instantiate()
end

function findartifactpath(pth, name)
    for d in readdir(pth; join=true)
        if any(occursin.(lowercase(name), lowercase.(readdir(joinpath(d, "logs")))))
            return d
        end
    end
end

pathof(app) = normpath(homedir(), ".julia/dev", app)

function clean_app(path, app_name)
    # remove all dlls in mkl artifact but keep
    # mkl_core.1.dll and mkl_rt.1.dll

    artifact_path = joinpath(path, "share", "julia", "artifacts")
    if !ispath(artifact_path)
        artifact_path = joinpath(path, "artifacts")
    end

    mkl_path = findartifactpath(artifact_path, "mkl")
    if isdir(joinpath(mkl_path, "bin"))
        for f in readdir(joinpath(mkl_path, "bin"); join=true)
            if !(occursin("mkl_core.1.dll", f) || occursin("mkl_rt.1.dll", f))
                rm(f)
            end
        end
    end
    if isdir(joinpath(mkl_path, "lib"))
        for f in readdir(joinpath(mkl_path, "lib"); join=true)
            if !(occursin("libmkl_core.so", f) || occursin("libmkl_rt.so", f))
                rm(f)
            end
        end
    end

    try
        test_romeo(path, app_name)
    catch
        @warn("Artifact cleaning failed! Please recompile clearswi with the option `clean=false`. The artifacts folder will be very large but some of them might not needed and can be manually removed.")
    end
end

function copy_matlab(path)
    cp(joinpath(dirname(@__DIR__), "matlab"), joinpath(path, "matlab"))
end


text:
# written by Jeffrey Sarnoff, Feb 2020.
# the constants
# ---------------

#=
   One table is split into subsections of 16 entries each.
   This keeps the scan time minimal, even after accounting
   for the conditionals that select the proper subsection.

   These are Float8 values offset by 1/2 way to the next value,
   this lets `findfirst` return an appropriately rounded value.

   The `F8offsetN` tuples are used with normal values,
   values that are finite and are not subnormal. Values
   that round to zero are handled before these are used.
=#

const F8offset1 = (Float32[
   0.2578125, 0.2734375, 0.2890625, 0.3046875, 0.3203125, 0.3359375,
   0.3515625, 0.3671875, 0.3828125, 0.3984375, 0.4140625, 0.4296875,
   0.4453125, 0.4609375, 0.4765625, 0.4921875]...,)

const F8offset2 = (Float32[
    0.515625, 0.546875, 0.578125, 0.609375, 0.640625, 0.671875,
    0.703125, 0.734375, 0.765625, 0.796875, 0.828125, 0.859375,
    0.890625, 0.921875, 0.953125, 0.984375]...,);

const F8offset3 = (Float32[
    1.03125, 1.09375, 1.15625, 1.21875, 1.28125, 1.34375, 1.40625,
    1.46875, 1.53125, 1.59375, 1.65625, 1.71875, 1.78125, 1.84375,
    1.90625, 1.96875]...,);

const F8offset4 = (Float32[
    2.0625, 2.1875, 2.3125, 2.4375, 2.5625, 2.6875, 2.8125,
    2.9375, 3.0625, 3.1875, 3.3125, 3.4375, 3.5625, 3.6875,
    3.8125, 3.9375]...,);

const F8offset5 = (Float32[
    4.125, 4.375, 4.625, 4.875, 5.125, 5.375, 5.625, 5.875, 6.125,
    6.375,  6.625, 6.875, 7.125, 7.375, 7.625, 7.875]...,);

const F8offset6 = (Float32[
    8.25, 8.75, 9.25, 9.75, 10.25, 10.75, 11.25, 11.75, 12.25, 12.75,
    13.25, 13.75, 14.25, 14.75, 15.25, 15.75]...,);
#=
     There is one table used with subnormal values.
      It is derived from the actual values of each
      subnormal quantity, shifted up halfway to the
      next subnormal.  This lets scanning also round.
      An initial (anchor) value is prepended, that value
      is half of the smallest subnormal.

     A corresponding table of UInt8 values also is used.
=#

const F8offset_subnormal = (
    0.0078125f0, 0.0234375f0, 0.0390625f0, 0.0546875f0, 0.0703125f0,
    0.0859375f0, 0.1015625f0, 0.1171875f0, 0.1328125f0, 0.1484375f0,
    0.1640625f0, 0.1796825f0, 0.1953125f0, 0.2109375f0, 0.2265625f0,
    0.2421875f0)

const U8subnormal = (collect(UInt8.(0:15))...,)

# some named constants to clarify the source text

const roundsto_floatmax8 = 15.25f0
const roundsto_zero8 = 0.0078125f0
const roundsto_subnormal = 0.2421875f0
const floatmaxplus8 = 15.75f0 # floatmax(Float8) + floatmin(Float8)

const UNaN8 = 0x78
const UInf8 = 0x70
const UFloatmax8 = 0x6f
const UFloatmin8 = 0x10
const UZero8 = 0x00

#  the functions
# ----------------

function Float8(x::Float32)
    # s, absx = signbit(x), abs(x)
    ui = toUInt8(x)
    return reinterpret(Float8, ui)
end

function toUInt8(x::Float32)
    s, absx = signbit(x), abs(x)
    isnan(absx) && return s ? UNaN8|0x80 : UNaN8
    if absx >= roundsto_floatmax8
        if absx > floatmaxplus8
            return s ? UInf8|0x80 : UInf8
        else
            return s ? UFloatmax8|0x80 : UFloatmax8
        end
    end
    if absx < roundsto_zero8
        return s ? UZero8|0x80 : UZero8
    elseif absx < roundsto_subnormal
        return subnormal8(s, absx)
    end
    absx = min(15.5f0, max(0.25f0, absx))
    return normal8(s, absx)
end

@inline function subnormal8(s::Core.Bool, absx::Float32)
    idx = findfirst(a->absx <= a, F8offset_subnormal)
    return s ? U8subnormal[idx] | 0x80 : U8subnormal[idx]
end

function normal8(s::Core.Bool, absx::Float32)
   if absx <= 1.96875f0
     if absx <= 0.4921875f0
          idx = UInt8(15+firstof16lte(absx, F8offset1))
          return s ? idx|0x80 : idx
     elseif absx <= 0.984375f0
          idx = UInt8(15+16+firstof16lte(absx, F8offset2))
          return s ? idx|0x80 : idx
     else
         idx = UInt8(15+32+firstof16lte(absx, F8offset3))
         return s ? idx|0x80 : idx
     end
   else
     if absx <= 3.9375f0
         idx = UInt8(15+32+16+firstof16lte(absx, F8offset4))
         return s ? idx|0x80 : idx
     elseif absx <= 7.875f0
         idx = UInt8(15+64+firstof16lte(absx, F8offset5))
         return s ? idx|0x80 : idx
     else
         idx = UInt8(15+64+16+firstof16lte(absx, F8offset6))
         return s ? idx|0x80 : idx
     end
   end
end

function firstof16lte(needle, haystack)
    for idx = 1:16
        if needle <= haystack[idx]
            return idx
        end
    end
    error("should not be reached")
end

# """Old version, slower."""
# function normal8(s::Bool, absx::Float32)
#    if absx <= 0.4921875f0
#         idx = UInt8(15+findfirst(a->absx <= a, F8offset1))
#         return s ? idx|0x80 : idx
#    elseif absx <= 0.984375f0
#         idx = UInt8(15+16+findfirst(a->absx <= a, F8offset2))
#         return s ? idx|0x80 : idx
#    elseif absx <= 1.96875f0
#        idx = UInt8(15+32+findfirst(a->absx <= a, F8offset3))
#        return s ? idx|0x80 : idx
#    elseif absx <= 3.9375f0
#        idx = UInt8(15+32+16+findfirst(a->absx <= a, F8offset4))
#        return s ? idx|0x80 : idx
#    elseif absx <= 7.875f0
#        idx = UInt8(15+64+findfirst(a->absx <= a, F8offset5))
#        return s ? idx|0x80 : idx
#    elseif absx <= 15.75f0
#        idx = UInt8(15+64+16+findfirst(a->absx <= a, F8offset6))
#        return s ? idx|0x80 : idx
#    else
#        throw(DomainError(absx,"not normal for Float8"))
#    end
# end


text:
using Gnuplot
x = 0:0.01:2π
@gp xlab = "x" ylab = "sin"
for i in 1:50
    @gp :- i x sin.(x .+ i / 10.0) "w l lc 'black' notit"
end
save(term="gif animate size 600,400 delay 1", output="basicAnim.gif")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "rdma_core"
version = v"33.1.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/linux-rdma/rdma-core/releases/download/v33.1/rdma-core-33.1.tar.gz", "d179b102bec551ce62265ed463d1095fb2ae9baff604261ad63327fcd20650e5")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd rdma-core-*
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DNO_MAN_PAGES=ON
make -j
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("powerpc64le", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="musl"),
    Platform("x86_64", "linux"; libc="musl"),
    Platform("aarch64", "linux"; libc="musl")
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libibverbs", :libibverbs),
    LibraryProduct("librdmacm", :lbrdmacm)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="libnl_jll", uuid="7c700954-19d3-5208-81e2-8fa5fe7c0bd8"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"5.2.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using .BlackBoxOptim: bboptimize, best_candidate, best_fitness

export BBO
Base.@kwdef struct BBO <: ComradeSoss.AbstractOptimizer
    tracemode::Symbol = :compact
    maxevals::Int = 10_000
end

function optimize(opt::BBO, lj::Soss.ConditionalModel)
    tc = ascube(lj)
    lower = transform(tc, zeros(dimension(tc)))
    upper = transform(tc, 0.999*ones(dimension(tc)))

    upflat, _ = ParameterHandling.flatten(upper)
    lowflat, unflatten = ParameterHandling.flatten(lower)

    bounds = [(lowflat[i], upflat[i]) for i in eachindex(upflat)]
    println(bounds)
    f(x) = -logdensity(lj, unflatten(x))
    res = bboptimize(f, SearchRange=bounds, MaxFuncEvals=opt.maxevals, TraceMode=opt.tracemode)
    return unflatten(best_candidate(res)), -best_fitness(res)
end


text:
"""
emulate function of
'LinearAlgebraicRepresentation' package implemented by Alberto Paoluzzi.

We are using functions:
	characteristicMatrix()
	cuboidGrid()
and type definitions:
	Points
	Cells
	Cell
	LARmodel
	Chains
"""

# TODO better implementation of characteristicMatrix. Probably can use just my own function

module Lar
	using SparseArrays
	using LinearAlgebra
		"""
		Points = Array{Number,2}

	Alias declation of LAR-specific data structure.
	Dense `Array{Number,2,1}` ``M x N`` to store the position of *vertices* (0-cells)
	of a *cellular complex*. The number of rows ``M`` is the dimension
	of the embedding space. The number of columns ``N`` is the number of vertices.
	"""
	const Points = Matrix

	const Cells = Array{Array{Int,1},1}
    const LARmodel = Tuple{Points,Array{Cells,1}}
	const Cell = SparseVector{Int8, Int}

	"""
		Chain = SparseArrays.SparseVector{Int8,Int}

	Alias declation of LAR-specific data structure.
	Binary `SparseVector` to store the coordinates of a `chain` of `N-cells`. It is
	`nnz=1` with `value=1` for the coordinates of an *elementary N-chain*, constituted by
	a single *N-chain*.
	"""
	const Chain = SparseArrays.SparseVector{Int8,Int}

	"""
		ChainOp = SparseArrays.SparseMatrixCSC{Int8,Int}

	Alias declation of LAR-specific data structure.
	`SparseMatrix` in *Compressed Sparse Column* format, contains the coordinate
	representation of an operator between linear spaces of `P-chains`.
	Operators ``P-Boundary : P-Chain -> (P-1)-Chain``
	and ``P-Coboundary : P-Chain -> (P+1)-Chain`` are typically stored as
	`ChainOp` with elements in ``{-1,0,1}`` or in ``{0,1}``, for
	*signed* and *unsigned* operators, respectively.
	"""
	const ChainOp = SparseArrays.SparseMatrixCSC{Int8,Int}

	"""
		ChainComplex = Array{ChainOp,1}

	Alias declation of LAR-specific data structure. It is a
	1-dimensional `Array` of `ChainOp` that provides storage for either the
	*chain of boundaries* (from `D` to `0`) or the transposed *chain of coboundaries*
	(from `0` to `D`), with `D` the dimension of the embedding space, which may be either
	``R^2`` or ``R^3``.
	"""
	const ChainComplex = Array{ChainOp,1}

	"""
		grid_0(n::Int)::Array{Int64,2}
	Generate a *uniform 0D cellular complex*.
	The `grid_0` function generates a 0-dimensional uniform complex embedding ``n+1`` equally-spaced  0-cells (at *unit interval* boundaries). It returns by columns the cells of this 0-complex as `Array{Int64,2}.
	#	Example
	```julia
	julia> grid_0(10)
	# output
	1×11 Array{Int64,2}:
	 0  1  2  3  4  5  6  7  8  9  10
	```
	"""
	function grid_0(n::Int)::Array{Int64,2}
	    return hcat([[i] for i in range(0, length=n+1)]...)
	end

	"""
		grid_1(n::Int)::Array{Int64,2}
	Generate a *uniform 1D cellular complex*.
	The `grid_1` function generates a 0-dimensional uniform complex embedding ``n+1`` equally-spaced  1-cells (*unit intervals*). It returns by columns the cells of this 1-complex as `Array{Int64,2}`.
	#	Example
	```julia
	julia> grid_1(10)
	# output
	2×10 Array{Int64,2}:
	 0  1  2  3  4  5  6  7  8   9
	 1  2  3  4  5  6  7  8  9  10
	```
	"""
	function grid_1(n)
	    return hcat([[i,i+1] for i in range(0, length=n)]...)
	end

	"""
		larGrid(n::Int)(d::Int)::Array{Int64,2}
	Generate either a *uniform 0D cellular complex* or a *uniform 1D cellular complex*.
	A `larGrid` function is given to generate the LAR representation of the cells of either a 0- or a 1-dimensional complex, depending on the value of the `d` parameter, to take values in the set ``{0,1}``, and providing the *order* of the output complex.
	#	Example
	```julia
	julia> larGrid(10)(0)
	# output
	1×11 Array{Int64,2}:
	 0  1  2  3  4  5  6  7  8  9  10
	julia> larGrid(10)(1)
	# output
	2×10 Array{Int64,2}:
	 0  1  2  3  4  5  6  7  8   9
	 1  2  3  4  5  6  7  8  9  10
	```
	"""
	function larGrid(n::Int)
	    function larGrid1(d::Int)::Array{Int64,2}
	        if d==0
	         return grid_0(n)
	        elseif d==1
	         return grid_1(n)
	        end
	    end
	    return larGrid1
	end

	"""
		cart(args::Array{Array{Any,1},1})::Array{Tuple,1}

	Cartesian product of collections given in the
	 unary `Array` argument. Return an `Array` of `Tuple`. The number ``n`` of output `Tuple` is equal to the *product of sizes* of input `args`
	#	Example
	```julia
	julia> cart([[1,2,3],["a","b"],[11,12]])
	# output
	12-element Array{Tuple{Any,Any,Any},1}:
	 (1, "a", 11)
	 (1, "a", 12)
	 (1, "b", 11)
	 (1, "b", 12)
	 (2, "a", 11)
	 (2, "a", 12)
	 (2, "b", 11)
	 (2, "b", 12)
	 (3, "a", 11)
	 (3, "a", 12)
	 (3, "b", 11)
	 (3, "b", 12)
	```
	"""
	function cart(args)::Array{Tuple,1}
	   return sort(vcat(collect(Iterators.product(args...))...))
	end

	"""
		larVertProd(vertLists::Array{Points,1})::Points
	Generate the integer *coordinates of vertices* (0-cells) of a *multidimensional grid*.
	*Grid n-vertices* are produced by the `larVertProd` function, via Cartesian product of vertices of ``n`` 0-dimensional arguments (vertex arrays in `vertLists`), orderly corresponding to ``x_1, x_2, ..., x_n`` coordinates in the output points ``(x_1, x_2,...,x_n)`` in ``R^n``.
	#	Example
	```julia
	julia> larVertProd([ larGrid(3)(0), larGrid(4)(0) ])
	# output
	2×20 Array{Int64,2}:
	 0  0  0  0  0  1  1  1  1  1  2  2  2  2  2  3  3  3  3  3
	 0  1  2  3  4  0  1  2  3  4  0  1  2  3  4  0  1  2  3  4
	```
	"""
	function larVertProd(vertLists::Array{Array{Int64,2},1})::Array{Int64,2}
	   coords = [[x[1] for x in v] for v in Lar.cart(vertLists)]
	   return sortslices(hcat(coords...), dims=2)
	end
	function larVertProd(vertLists::Array{Array{Float64,2},1})::Array{Float64,2}
	   coords = [[x[1] for x in v] for v in Lar.cart(vertLists)]
	   return sortslices(hcat(coords...), dims=2)
	end

	"""
		larCuboids( shape, filled=false )::Union( Cells, Array{Cells,1} )

	Multi-dimensional generator function.
	Generate either a solid *``d``-grid* of unit *``d``-cuboids* in ``d``-dimensional space, or the array of ``p``-skeletons (``0 <=p<= d``), depending on the Boolean variable `filled`. ``0``-cuboids are points, ``1``-cuboids are segments, , ``2``-cuboids are squares,  ``3``-cuboids are cubes, etc. The `shape=[a,b,c]` value determines the number ``a x b x c`` of ``d``-cells. Notice that `d = length(shape)`
	"""
	function larCuboids( shape, filled=false )
	   vertGrid = larImageVerts(shape)
	   gridMap = larGridSkeleton(shape)
	   if ! filled
	      cells = gridMap(length(shape))
	   else
	      skeletonIds = 0:length(shape)
	      cells = [ gridMap(id) for id in skeletonIds ]
	   end
	   return convert(Array{Float64,2},vertGrid), cells
	end

	cuboidGrid = larCuboids

	"""
		larCellProd(cellLists::Array{Cells,1})::Cells
	Generation of *grid cells* by *Cartesian product* of 0/1-complexes.
	The *output complex* is generated by the product of *any number* of either 0- or 1-dimensional cell complexes. The product of ``d`` 1-complexes generates *solid ``d``-cells*, while the product of ``n`` 0-complexes and ``d-n`` 1-complexes (``n < d``) generates *non-solid ``(d-n)``-cells*, properly embedded in ``d``-space, i.e. with vertices having ``d`` coordinates.
	# Examples
	To understand the *generation of cuboidal grids* from products of 0- or 1-dimensional complexes, below we show a simple example of 2D grids embedded in ``R^3``.
	In particular, `v1 = [0. 1. 2. 3.]` and `v0 = [0. 1. 2.]` are two 2-arrays of 1D vertices, `c1 = [[0,1],[1,2],[2,3]]` and `c0 = [[0],[1],[2]]` are the LAR representation of one *``1``-complex* and one *``0``-complex*, respectively. The solid 2-complex named `grid2D` is generated in 2D as follows:
	```julia
	julia> v1 = [0. 1. 2. 3.]
	1×4 Array{Float64,2}:
	 0.0  1.0  2.0  3.0
	julia> c1 = [[0,1],[1,2],[2,3]]
	3-element Array{Array{Int64,1},1}:
	 [0, 1]
	 [1, 2]
	 [2, 3]
	julia> grid2D = larVertProd([v1,v1]),larCellProd([c1,c1])
	([0.0 0.0 … 3.0 3.0; 0.0 1.0 … 2.0 3.0], Array{Int64,1}[[1, 2, 5, 6], [2, 3, 6, 7], [3, 4, 7, 8], [5, 6, 9, 10], [6, 7, 10, 11], [7, 8, 11, 12], [9, 10, 13, 14], [10, 11, 14, 15], [11, 12, 15, 16]])
	```
	whereas a *non-solid* ``2``-complex in ``3D`` is generated as:
	```julia
	julia> v1, c1 = [0. 1. 2. 3.],[[0,1],[1,2],[2,3]]
	([0.0 1.0 2.0 3.0], Array{Int64,1}[[0, 1], [1, 2], [2, 3]])
	julia> v0, c0 = [0. 1. 2.], [[0],[1],[2]]
	([0.0 1.0 2.0], Array{Int64,1}[[0], [1], [2]])
	julia> vertGrid = larVertProd([v1, v1, v0])
	3×48 Array{Float64,2}:
	 0.0  0.0  0.0  0.0  0.0  0.0  …  3.0  3.0  3.0  3.0  3.0  3.0  3.0  3.0  3.0
	 0.0  0.0  0.0  1.0  1.0  1.0  …  1.0  1.0  1.0  2.0  2.0  2.0  3.0  3.0  3.0
	 0.0  1.0  2.0  0.0  1.0  2.0  …  0.0  1.0  2.0  0.0  1.0  2.0  0.0  1.0  2.0
	julia> cellGrid = larCellProd([c1, c1, c0])
	27-element Array{Array{Int64,1},1}:
	 [1, 4, 13, 16]
	 [2, 5, 14, 17]
	 ...  ... ...
	 [32, 35, 44, 47]
	 [33, 36, 45, 48]
	julia> grid3D = vertGrid,cellGrid
	([0.0 0.0 … 3.0 3.0; 0.0 0.0 … 3.0 3.0; 0.0 1.0 … 1.0 2.0], Array{Int64,1}[[1, 4, 13, 16], [2, 5, 14, 17], … [32, 35, 44, 47], [33, 36, 45, 48]])
	julia> using Plasm
	julia> Plasm.view(grid3D)
	```
	"""
	function larCellProd(cellLists::Array{Cells,1})::Cells
	   shapes = [length(item) for item in cellLists]
	   subscripts = cart([collect(range(0, length=shape)) for shape in shapes])
	   indices = [collect(tuple) .+ 1 for tuple in subscripts]

	   jointCells = [cart([cells[k] for (k,cells) in zip(index,cellLists)])
	   				for index in indices]
	   convertIt = index2addr([ (length(cellLists[k][1]) > 1) ? shape .+ 1 : shape
	      for (k,shape) in enumerate(shapes) ])
	   [vcat(map(convertIt, map(collect,jointCells[j]))...) for j in 1:length(jointCells)]
	end

	"""
		filterByOrder( n::Int )Array{Array{Array{Int8,1},1},1}
	Filter the `array` of codes  (`Boolean` `String`) of *``n`` bits* depending on their integer value (*order*).
	# Example
	```julia
	julia> filterByOrder(3)
	# output
	4-element Array{Array{Array{Int8,1},1},1}:
	 Array{Int8,1}[Int8[0, 0, 0]]
	 Array{Int8,1}[Int8[0, 0, 1], Int8[0, 1, 0], Int8[1, 0, 0]]
	 Array{Int8,1}[Int8[0, 1, 1], Int8[1, 0, 1], Int8[1, 1, 0]]
	 Array{Int8,1}[Int8[1, 1, 1]]
	```"""
	function filterByOrder(n::Int)Array{Array{Array{Int8,1},1},1}
	   terms = [[parse(Int8,bit) for bit in collect(term)] for term in Lar.binaryRange(n)]
	   return [[term for term in terms if sum(term) == k] for k in 0:n]
	end

	"""
		larGridSkeleton( shape::Array{Int,1} )( d::Int )::Cells

	Produce the `d`-dimensional skeleton (set of `d`-cells) of a cuboidal grid of given `shape`.
	# Example
	A `shape=[1,1,1]` parameter refers to a *grid* with a single step on the three axes, i.e. to a single *3D unit cube*. Below all *skeletons* of such simplest grid are generated.
	```julia
	julia> Lar.larGridSkeleton([1,1,1])(0)
	# output
	8-element Array{Array{Int64,1},1}:
	[[1], [2], [3], [4], [5], [6], [7], [8]]
	julia> Lar.larGridSkeleton([1,1,1])(1)
	# output
	12-element Array{Array{Int64,1},1}:
	[[1,2],[3,4],[5,6],[7,8],[1,3],[2,4],[5,7],[6,8],[1,5],[2,6],[3,7],[4,8]]
	julia> Lar.larGridSkeleton([1,1,1])(2)
	# output
	6-element Array{Array{Int64,1},1}:
	[[1,2,3,4], [5,6,7,8], [1,2,5,6], [3,4,7,8], [1,3,5,7], [2,4,6,8]]
	julia> Lar.larGridSkeleton([1,1,1])(3)
	# output
	1-element Array{Array{Int64,1},1}:
	 [1, 2, 3, 4, 5, 6, 7, 8]
	```
	"""
	function larGridSkeleton(shape)
	    n = length(shape)
	    function larGridSkeleton0( d::Int )::Cells

	    	@assert d<=n

	        components = filterByOrder(n)[d .+ 1]
	        apply(fun,a) = fun(a)
			componentCellLists = [ [map(f,x)  for (f,x) in  zip( [larGrid(dim)
				for dim in shape], convert(Array{Int64,1},component) ) ]
					for component in components ]
	        colList(arr) = [arr[:,k]  for k in 1:size(arr,2)]
	        out = [ larCellProd(map(colList,cellLists)) for cellLists in componentCellLists ]
	        return vcat(out...)
	    end
	    return larGridSkeleton0
	end

	"""
		larImageVerts(shape::Array{Int,1})::Array{Int64,2}
	Linearize the *grid of integer vertices*, given the `shape` of a *cuboidal grid* (typically an *image*).
	# Examples
	```julia
	julia> larImageVerts([1024,1024])
	# output
	2×1050625 Array{Int64,2}:
	 0  0  0  0  0  0  0  0  0  0   0   0 … 1024  1024  1024  1024  1024  1024  1024  1024
	 0  1  2  3  4  5  6  7  8  9  10  11 … 1017  1018  1019  1020  1021  1022  1023  1024

	julia> larImageVerts([1,1,1])
	# output
	3×8 Array{Int64,2}:
	 0  0  0  0  1  1  1  1
	 0  0  1  1  0  0  1  1
	 0  1  0  1  0  1  0  1
	```
	"""
	function larImageVerts( shape::Array{Int,1} )::Array{Int64,2}
	   vertexDomain(n) = hcat([k for k in 0:n-1]...)
	   vertLists = [vertexDomain(k+1) for k in shape]
	   vertGrid = larVertProd(vertLists)
	   return vertGrid
	end

	"""
    binaryRange(n)
	Generate the first `n` binary numbers in string padded for max `2^n` length
	"""
	function binaryRange(n)
	    return string.(range(0, length=2^n), base=2, pad=n)
	end

	"""
		index2addr(shape::Array{Int64,1})(multiIndex)::Int
	*Multi-index to address* transformation. Multi-index is a *generalization* of the concept of an integer index to an *ordered tuple of indices*.
	The second-order utility function `index2addr`  transforms a `shape` list for a *multidimensional array* into a function that, when applied to a *multindex array*, i.e. to a list of integer `Tuple` within the `shape`'s bounds, returns the *integer addresses* of the corresponding array components within the *linear storage* of the multidimensional array.
	#	Example
	Notice that in the example below, there are ``3 x 6`` different *multi-index values* for the variable `index`, generated by `cart([ 0:2, 0:5 ])`.
	```julia
	julia> [index2addr([3,6])(collect(index)) for index in cart([ 0:2, 0:5 ])]'
	# output
	1×18 RowVector{Int64,Array{Int64,1}}:
	 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18

	julia> index2addr([3,6])([0,0])
	# output
	1
	julia> index2addr([3,6])([2,5])
	# output
	18
	```
	"""
	function index2addr( shape::Array{Int64,2} )
	    n = length(shape)
	    theShape = append!(shape[2:end],1)
	    weights = [prod(theShape[k:end]) for k in range(1, length=n)]

	    function index2addr0( multiIndex::Array{Int,1} )::Int
	        return dot(collect(multiIndex), weights) + 1
	    end

	    return index2addr0
	end

	"""
		index2addr(shape::Array{Int64,2})(multiIndex::Array{Int,1})::Int
	Multi-index to address transformation. Partial function allowing for using both horizontal and vertical vectors for `shape` parameter. Notice that multi-indices are used here as *coordinates of grid points*, hence they start from tuples of zeros. Accordingly, the translation formula for multi-index to address transformation is *0-based*.
	"""
	function index2addr(shape::Array{Int64,1})
	   index2addr(hcat(shape...))
	end

	"""
			characteristicMatrix( FV::Cells )::ChainOp
		Binary matrix representing by rows the `p`-cells of a cellular complex.
		The input parameter must be of `Cells` type. Return a sparse binary matrix,
		providing the basis of a ``Chain`` space of given dimension. Notice that the
		number of columns is equal to the number of vertices (0-cells).
		# Example
		```julia
		V,(VV,EV,FV,CV) = cuboid([1.,1.,1.], true);
		julia> Matrix(characteristicMatrix(FV))
		6×8 Array{Int8,2}:
		1  1  1  1  0  0  0  0
		0  0  0  0  1  1  1  1
		1  1  0  0  1  1  0  0
		0  0  1  1  0  0  1  1
		1  0  1  0  1  0  1  0
		0  1  0  1  0  1  0  1
		julia> Matrix(characteristicMatrix(CV))
		1×8 Array{Int8,2}:
		1  1  1  1  1  1  1  1
		julia> Matrix(characteristicMatrix(EV))
		12×8 Array{Int8,2}:
		1  1  0  0  0  0  0  0
		0  0  1  1  0  0  0  0
		0  0  0  0  1  1  0  0
		0  0  0  0  0  0  1  1
		1  0  1  0  0  0  0  0
		0  1  0  1  0  0  0  0
		0  0  0  0  1  0  1  0
		0  0  0  0  0  1  0  1
		1  0  0  0  1  0  0  0
		0  1  0  0  0  1  0  0
		0  0  1  0  0  0  1  0
		0  0  0  1  0  0  0  1
		```
	"""
	function characteristicMatrix( FV::Cells )::ChainOp
		I,J,V = Int64[],Int64[],Int8[]
		for f=1:length(FV)
			for k in FV[f]
			push!(I,f)
			push!(J,k)
			push!(V,1)
			end
		end
		M_2 = sparse(I,J,V)
		return M_2
	end

	# """
	# cuboidGrid( shape, filled=false )::Union( Cells, Array{Cells,1} )
	# Multi-dimensional generator function.
	# Generate either a solid *``d``-grid* of unit *``d``-cuboids* in ``d``-dimensional space, or the array of ``p``-skeletons (``0 <=p<= d``), depending on the Boolean variable `filled`. ``0``-cuboids are points, ``1``-cuboids are segments, , ``2``-cuboids are squares,  ``3``-cuboids are cubes, etc. The `shape=[a,b,c]` value determines the number ``a x b x c`` of ``d``-cells. Notice that `d = length(shape)`
	# """
	# function cuboidGrid( shape, filled=false )
	#    vertGrid = larImageVerts(shape)
	#    gridMap = larGridSkeleton(shape)
	#    if ! filled
	#       cells = gridMap(length(shape))
	#    else
	#       skeletonIds = 0:length(shape)
	#       cells = [ gridMap(id) for id in skeletonIds ]
	#    end
	#    return convert(Array{Float64,2},vertGrid), cells
	# end
	#

	"""
    lar2obj(V::Points, cc::ChainComplex)
Triangulated OBJ string representation of the model passed as input.
Use this function to export LAR models into OBJ
# Example
```julia
	julia> cube_1 = ([0 0 0 0 1 1 1 1; 0 0 1 1 0 0 1 1; 0 1 0 1 0 1 0 1],
	[[1,2,3,4],[5,6,7,8],[1,2,5,6],[3,4,7,8],[1,3,5,7],[2,4,6,8]],
	[[1,2],[3,4],[5,6],[7,8],[1,3],[2,4],[5,7],[6,8],[1,5],[2,6],[3,7],[4,8]] )
	julia> cube_2 = Lar.Struct([Lar.t(0,0,0.5), Lar.r(0,0,pi/3), cube_1])
	julia> V, FV, EV = Lar.struct2lar(Lar.Struct([ cube_1, cube_2 ]))
	julia> V, bases, coboundaries = Lar.chaincomplex(V,FV,EV)
	julia> (EV, FV, CV), (copEV, copFE, copCF) = bases, coboundaries
	julia> FV # bases[2]
	18-element Array{Array{Int64,1},1}:
	 [1, 3, 4, 6]
	 [2, 3, 5, 6]
	 [7, 8, 9, 10]
	 [1, 2, 3, 7, 8]
	 [4, 6, 9, 10, 11, 12]
	 [5, 6, 11, 12]
	 [1, 4, 7, 9]
	 [2, 5, 11, 13]
	 [2, 8, 10, 11, 13]
	 [2, 3, 14, 15, 16]
	 [11, 12, 13, 17]
	 [11, 12, 13, 18, 19, 20]
	 [2, 3, 13, 17]
	 [2, 13, 14, 18]
	 [15, 16, 19, 20]
	 [3, 6, 12, 15, 19]
	 [3, 6, 12, 17]
	 [14, 16, 18, 20]
	julia> CV # bases[3]
	3-element Array{Array{Int64,1},1}:
	 [2, 3, 5, 6, 11, 12, 13, 14, 15, 16, 18, 19, 20]
	 [2, 3, 5, 6, 11, 12, 13, 17]
	 [1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 17]
	julia> copEV # coboundaries[1]
	34×20 SparseMatrixCSC{Int8,Int64} with 68 stored entries: ...
	julia> copFE # coboundaries[2]
	18×34 SparseMatrixCSC{Int8,Int64} with 80 stored entries: ...
	julia> copCF # coboundaries[3]
	4×18 SparseMatrixCSC{Int8,Int64} with 36 stored entries: ...
	objs = Lar.lar2obj(V'::Lar.Points, [coboundaries...])
	open("./two_cubes.obj", "w") do f
    	write(f, objs)
	end
```
"""
function lar2obj(V::Lar.Points, cc::Lar.ChainComplex)
    # copEV, copFE, copCF = cc
	copEV = cc[1]
	copFE = cc[2]
	if size(V,2) > 3
		V = convert(Lar.Points, V') # out V by rows
	end
    obj = ""
    for v in 1:size(V, 1)
        obj = string(obj, "v ",
    	round(V[v, 1], digits=6), " ",
    	round(V[v, 2], digits=6), " ",
    	round(V[v, 3], digits=6), "\n")
    end

    # print("Triangulating")
    triangulated_faces = Lar.triangulate(V, cc[1:2])
    # println("DONE")

    for c in 1:copCF.m
        obj = string(obj, "\ng cell", c, "\n")
        for f in copCF[c, :].nzind
            triangles = triangulated_faces[f]
            for tri in triangles
                #t = copCF[c, f] > 0 ? tri : tri[end:-1:1]
				t = tri
                obj = string(obj, "f ", t[1], " ", t[2], " ", t[3], "\n")
            end
        end
    end

    return obj
end

function lar2obj(V::Lar.Points, FVtri::Array, filename::String)
	objlines = lar2obj(V, FVtri)
	open(filename, "w") do f
		for line in objlines
	    	write(f, line)
		end
	end
	return objlines
end

function lar2obj(V::Lar.Points, FVtri::Array)
    # copEV, copFE, copCF = cc
	# copEV = cc[1]
	# copFE = cc[2]
	lines = []
	if size(V,2) > 3
		V = convert(Lar.Points, V') # out V by rows
	end
    obj = ""
    for v in 1:size(V, 1)
		push!(lines,
	        string("v ",
	    	round(V[v, 1], digits=6), " ",
	    	round(V[v, 2], digits=6), " ",
	    	round(V[v, 3], digits=6), "\n")
		)
        # obj = string(obj, "v ",
    	# round(V[v, 1], digits=6), " ",
    	# round(V[v, 2], digits=6), " ",
    	# round(V[v, 3], digits=6), "\n")
    end

    # print("Triangulating")
    # triangulated_faces = Lar.triangulate(V, cc[1:2])
    # println("DONE")
	c = "object1"
    # for c in 1:copCF.m
		push!(lines,
        	string("\ng cell", c, "\n")
		)
        # for f in copCF[c, :].nzind
            # triangles = triangulated_faces[f]
			triangles = FVtri
            for tri in triangles
                #t = copCF[c, f] > 0 ? tri : tri[end:-1:1]
				t = tri
				push!(lines,
                	string("f ", t[1], " ", t[2], " ", t[3], "\n")
					)
            end
        # end
    # end

    return lines
end

end


text:
module Primitives

export AGen, AFacet

using ..Common
import ..Common.nconic

global numgenerators = 0

immutable AGen{T} <: Generator
    x::Vector{T}
    id::Int

    function AGen(x)
        global numgenerators
        g = new(T[x...], numgenerators += 1)
        g
    end
end

nconic(g::AGen) = length(g.x)

Base.getindex(g::AGen, k)  = g.x[k]
Base.show(io::IO, g::AGen) = print(io, "g$(g.id)")

end # module


text:
# main tests

function safe_mapslices(op, A, region)
    newregion = intersect(region, 1:ndims(A))
    return isempty(newregion) ? A : mapslices(op, A, newregion)
end
safe_sum{T}(A::Array{T}, region) = safe_mapslices(sum, A, region)
safe_prod{T}(A::Array{T}, region) = safe_mapslices(prod, A, region)
safe_maximum{T}(A::Array{T}, region) = safe_mapslices(maximum, A, region)
safe_minimum{T}(A::Array{T}, region) = safe_mapslices(minimum, A, region)
safe_sumabs{T}(A::Array{T}, region) = safe_mapslices(sum, abs(A), region)
safe_sumabs2{T}(A::Array{T}, region) = safe_mapslices(sum, abs2(A), region)

Areduc = rand(3, 4, 5, 6)
for region in {
        1, 2, 3, 4, 5, (1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4),
    (1, 2, 3), (1, 3, 4), (2, 3, 4), (1, 2, 3, 4)}

    r = fill(NaN, Base.reduced_dims(size(Areduc), region))
    @test_approx_eq sum!(r, Areduc) safe_sum(Areduc, region)
    @test_approx_eq prod!(r, Areduc) safe_prod(Areduc, region)
    @test_approx_eq maximum!(r, Areduc) safe_maximum(Areduc, region)
    @test_approx_eq minimum!(r, Areduc) safe_minimum(Areduc, region)
    @test_approx_eq Base.sumabs!(r, Areduc) safe_sumabs(Areduc, region)
    @test_approx_eq Base.sumabs2!(r, Areduc) safe_sumabs2(Areduc, region)

    @test_approx_eq sum(Areduc, region) safe_sum(Areduc, region)
    @test_approx_eq prod(Areduc, region) safe_prod(Areduc, region)
    @test_approx_eq maximum(Areduc, region) safe_maximum(Areduc, region)
    @test_approx_eq minimum(Areduc, region) safe_minimum(Areduc, region)
    @test_approx_eq Base.sumabs(Areduc, region) safe_sumabs(Areduc, region)
    @test_approx_eq Base.sumabs2(Areduc, region) safe_sumabs2(Areduc, region)
end

# Test reduction along first dimension; this is special-cased for
# size(A, 1) >= 16
Breduc = rand(64, 3)
r = fill(NaN, Base.reduced_dims(size(Breduc), 1))
@test_approx_eq sum!(r, Breduc) safe_sum(Breduc, 1)
@test_approx_eq Base.sumabs!(r, Breduc) safe_sumabs(Breduc, 1)
@test_approx_eq Base.sumabs2!(r, Breduc) safe_sumabs2(Breduc, 1)
@test_approx_eq sum(Breduc, 1) safe_sum(Breduc, 1)
@test_approx_eq Base.sumabs(Breduc, 1) safe_sumabs(Breduc, 1)
@test_approx_eq Base.sumabs2(Breduc, 1) safe_sumabs2(Breduc, 1)

@test reducedim((a,b) -> a|b, [true false; false false], 1, false) == [true false]
R = reducedim((a,b) -> a+b, [1 2; 3 4], 2, 0.0)
@test eltype(R) == Float64
@test_approx_eq R [3,7]
@test reducedim((a,b) -> a+b, [1 2; 3 4], 1, 0) == [4 6]

# inferred return types
rt = Base.return_types(reducedim, (Function, Array{Float64, 3}, Int, Float64))
@test length(rt) == 1 && rt[1] == Array{Float64, 3}

## findmin/findmax
A = [1.0 3.0 6.0;
     5.0 2.0 4.0]
@test findmin(A, (1,)) == ([1.0 2.0 4.0], [1 4 6])
@test findmin(A, (2,)) == (reshape([1.0,2.0], 2, 1), reshape([1,4], 2, 1))
@test findmin(A, (1,2)) == (fill(1.0,1,1),fill(1,1,1))
@test findmax(A, (1,)) == ([5.0 3.0 6.0], [2 3 5])
@test findmax(A, (2,)) == (reshape([6.0,5.0], 2, 1), reshape([5,2], 2, 1))
@test findmax(A, (1,2)) == (fill(6.0,1,1),fill(5,1,1))

# issue #6672
@test sum(Real[1 2 3; 4 5.3 7.1], 2) == reshape([6, 16.4], 2, 1)
@test std(FloatingPoint[1,2,3], 1) == [1.0]
@test sum({1 2;3 4},1) == [4 6]
@test sum(Vector{Int}[[1,2],[4,3]], 1)[1] == [5,5]

max_stars_repo_name: vchuravy/julia

max_issues_repo_name: vchuravy/julia

max_forks_repo_name: vchuravy/julia


text:
module UndirectedGraphBackend

    include("./UndirectedGraphFunctions.jl")
    include("../backend_Models.jl")
    export execute_backend
    using .UndirectedGraphFunctions
    using .BackendModels

    function execute_backend(n_qregs::Int, gates, model)
        backend_gates = call_backend_gates(gates)
        parsed_gates = parse_gates(backend_gates)
        _gates = convert(Array{UG_Gate,1}, parsed_gates)
        if model.measure_all
            input_state = zeros(Int, n_qregs)
            output_state = zeros(Int, n_qregs)
            state_vec = []
            buckets = create_buckets(n_qregs, _gates, input_state, output_state)
            result = contract_graph(buckets)
            push!(state_vec, result)
            while true
                for i in 1:n_qregs
                    if output_state[n_qregs - (i-1)] == 0
                        output_state[n_qregs - (i-1)] = 1
                        output_state[n_qregs - (i-1) + 1:n_qregs] .= 0
                        break
                    end
                end
                buckets = create_buckets(n_qregs, _gates, input_state, output_state)
                result = contract_graph(buckets)
                push!(state_vec, result)
                if sum(output_state) >= n_qregs
                    break
                end
            end
            return state_vec
        else
            if length(model.output_basis) != n_qregs
                error("Measured state length must be equal to qubit number.")
            end
            input_state = zeros(Int, n_qregs)
            output_state = model.output_basis
            buckets = create_buckets(n_qregs, _gates, input_state, output_state)
            result = contract_graph(buckets)
            return result
        end

    end

    function call_backend_gates(gates)
        backend_gates = []
        for gate in gates
            name = gate._name
            fields = fieldnames(typeof(gate))
            field_names = []
            for i in fields
                push!(field_names, string(i))
            end
            given_property = [getproperty(gate, property) for property in fields]
            code = gate._name * "("
            for prop in given_property[1:end-1]
                code *= string(prop)
                code *= ", "
            end
            code *= ")"
            new_gate = eval(Meta.parse(code))
            push!(backend_gates, new_gate)
        end
        return backend_gates
    end

    function parse_gates(backend_gates)
        parsed_gates =[]
        for gate in backend_gates
            _gates = gate_parser(gate)
            for _gate in _gates
                push!(parsed_gates, _gate)
            end
        end
        return parsed_gates
    end

end

max_stars_repo_name: terasakisatoshi/blueqat-julia


text:
using Plots

plotlyjs()

x=[-1, -1, 1, 1, -1, -1, 1, 1]
y=[-1, 1, 1, -1, -1, 1, 1, -1]
z=[-1, -1, -1, -1, 1, 1, 1, 1]

i= [7, 0, 0, 0, 4, 4, 6, 1, 4, 0, 3, 6]
j= [3, 4, 1, 2, 5, 6, 5, 2, 0, 1, 6, 3]
k= [0, 7, 2, 3, 6, 7, 1, 6, 5, 5, 7, 2]

fig1 = scatter(x, y, z; xlabel = "x", ylabel = "y", zlabel = "z", legend = :none, margin = 2 * Plots.mm, framestyle = :grid)
display(fig1)

fig2 = mesh3d(x, y, z; connections = (i, j, k), xlabel = "x", ylabel = "y", zlabel = "z", legend = :none, margin = 2 * Plots.mm, framestyle = :none)
display(fig2)


text:
module correlations

using ..operators
using ..timeevolution
using ..metrics
using ..steadystate

"""
Calculate two time correlation values :math:`\\langle A(t) B(0) \\rangle`

The calculation is done by multiplying the initial density operator
with :math:`B` performing a time evolution according to a master equation
and then calculating the expectation value :math:`\\mathrm{Tr} \\{ A \\rho \\}`

Arguments
---------

tspan
    Points of time at which the correlation should be calculated.
rho0
    Initial density operator.
H
    Operator specifying the Hamiltonian.
J
    Vector of jump operators.
op1
    Operator at time t.
op2
    Operator at time t=0.

Keyword Arguments
-----------------

Gamma
    Vector or matrix specifying the coefficients for the jump operators.
Jdagger (optional)
    Vector containing the hermitian conjugates of the jump operators. If they
    are not given they are calculated automatically.
kwargs
    Further arguments are passed on to the ode solver.
"""
function correlation(tspan::Vector{Float64}, rho0::DenseOperator, H::Operator, J::Vector,
                     op1::Operator, op2::Operator;
                     Gamma::Union{Real, Vector, Matrix}=ones(Float64, length(J)),
                     Jdagger::Vector=map(dagger, J),
                     tmp::DenseOperator=deepcopy(rho0),
                     kwargs...)
    exp_values = Complex128[]
    function fout(t, rho)
        push!(exp_values, expect(op1, rho))
    end
    timeevolution.master(tspan, op2*rho0, H, J; Gamma=Gamma, Jdagger=Jdagger,
                        tmp=tmp, fout=fout, kwargs...)
    return exp_values
end

"""
Calculate two time correlation values :math:`\\langle A(t) B(0) \\rangle`

The calculation is done by multiplying the initial density operator
with :math:`B` performing a time evolution according to a master equation
and then calculating the expectation value :math:`\\mathrm{Tr} \\{ A \\rho \\}`.
The points of time are chosen automatically from the ode solver and the final
time is determined by the steady state termination criterion specified in
:func:`steadystate.master`.

Arguments
---------

rho0
    Initial density operator.
H
    Operator specifying the Hamiltonian.
J
    Vector of jump operators.
op1
    Operator at time t.
op2
    Operator at time t=0.

Keyword Arguments
-----------------

eps
    Tracedistance used as termination criterion.
h0
    Initial time step used in the time evolution.
Gamma
    Vector or matrix specifying the coefficients for the jump operators.
Jdagger (optional)
    Vector containing the hermitian conjugates of the jump operators. If they
    are not given they are calculated automatically.
kwargs
    Further arguments are passed on to the ode solver.
"""
function correlation(rho0::DenseOperator, H::Operator, J::Vector,
                     op1::Operator, op2::Operator;
                     eps::Float64=1e-4, h0=10.,
                     Gamma::Union{Real, Vector, Matrix}=ones(Float64, length(J)),
                     Jdagger::Vector=map(dagger, J),
                     tmp::DenseOperator=deepcopy(rho0),
                     kwargs...)
    op2rho0 = op2*rho0
    tout = Float64[0.]
    exp_values = Complex128[expect(op1, op2rho0)]
    function fout(t, rho)
        push!(tout, t)
        push!(exp_values, expect(op1, rho))
    end
    steadystate.master(H, J; rho0=op2rho0, eps=eps, h0=h0, fout=fout,
                       Gamma=Gamma, Jdagger=Jdagger, tmp=tmp, kwargs...)
    return tout, exp_values
end

function correlationspectrum(omega_samplepoints::Vector{Float64},
                H::Operator, J::Vector, op::Operator;
                eps::Float64=1e-4,
                rho_ss::DenseOperator=steadystate.master(H, J; eps=eps),
                kwargs...)
    domega = minimum(diff(omega_samplepoints))
    dt = 2*pi/(omega_samplepoints[end] - omega_samplepoints[1])
    T = 2*pi/domega
    tspan = [0.:dt:T;]
    exp_values = correlation(tspan, rho_ss, H, J, dagger(op), op, kwargs...)
    # dtmin = minimum(diff(tspan))
    # T = tspan[end] - tspan[1]
    # domega = 2*pi/T
    # omega_min = -pi/dtmin
    # omega_max = pi/dtmin
    # omega_samplepoints = Float64[omega_min:domega:omega_max-domega/2;]
    S = Float64[]
    for omega=omega_samplepoints
        y = exp(1im*omega*tspan).*exp_values/pi
        I = 0.im
        for j=1:length(tspan)-1
            I += (tspan[j+1] - tspan[j])*(y[j+1] + y[j])
        end
        I = I/2
        push!(S, real(I))
    end
    return omega_samplepoints, S
end

function correlationspectrum(H::Operator, J::Vector, op::Operator;
                eps::Float64=1e-4, h0=10.,
                rho_ss::DenseOperator=steadystate.master(H, J; eps=eps),
                kwargs...)
    tspan, exp_values = correlation(rho_ss, H, J, dagger(op), op, eps=eps, h0=h0, kwargs...)
    dtmin = minimum(diff(tspan))
    T = tspan[end] - tspan[1]
    tspan = Float64[0.:dtmin:T;]
    return correlationspectrum(tspan, H, J, op; eps=eps, rho_ss=rho_ss, kwargs...)
    # domega = 1./T
    # omega_min = -pi/dtmin
    # omega_max = pi/dtmin
    # omega_samplepoints = Float64[omega_min:domega:omega_max;]
    # S = Float64[]
    # for omega=omega_samplepoints
    #     y = exp(1im*omega*tspan).*exp_values/pi
    #     I = 0.im
    #     for j=1:length(tspan)-1
    #         I += (tspan[j+1] - tspan[j])*(y[j+1] + y[j])
    #     end
    #     I = I/2
    #     push!(S, real(I))
    # end
    # return omega_samplepoints, S
end

end # module


text:
while true
	x = "43" - 1
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Proj
using Test
using StaticArrays: SA

# Base.isapprox doesn't work on tuples
is_approx(a, b) = all(isapprox(c[1], c[2]) for c in zip(a, b))

@testset "Proj" begin
    include("libproj.jl")
    include("applications.jl")

end # testset "Proj"


text:
@polly function kernel_floyd_warshall(path)
    for k = 1:size(path,1), i = 1:size(path,1), j = 1:size(path,1)
        path[i,j] = path[i,j] < path[i,k] + path[k,j] ? path[i,j] : path[i,k] + path[k,j]
    end
end

let
    n = 2800

    path = zeros(Float32, n, n)

    for i = 1:n, j = 1:n
        path[i,j] = i*j%7+1
        if ((i+j-2)%13 == 0 || (i+j-2)%7==0 || (i+j-2)%11 == 0)
            path[i,j] = 999
        end
    end

    SUITE["floyd-warshall"] = @benchmarkable kernel_floyd_warshall(path) setup = (path = copy($path))
end


text:
export     AbstractDrillingCost,
    AbstractDrillingCost_TimeFE,
    DrillingCost_TimeFE,
    DrillingCost_TimeFE_rigrate,
    DrillingCost_constant,
    DrillingCost_dgt1,
    DrillingCost_TimeFE_costdiffs,
    DrillingCost_TimeFE_rig_costdiffs,
    DrillingCost_DoubleTimeFE

@inline flowdψ(::AbstractDrillingCost, d, obs, theta, sim) = azero(theta)
@inline function flow!(grad, u::AbstractPayoffFunction, d, obs, θ, sim)
    dflow!(u, grad, d, obs, θ, sim)
    return flow(u, d, obs, θ, sim)
end

@inline function flow!(u::AbstractPayoffFunction, grad, d, obs, θ, sim)
    dflow!(u, grad, d, obs, θ, sim)
    return flow(u, d, obs, θ, sim)
end

"Abstract Type for Costs w Fixed Effects"
abstract type AbstractDrillingCost_TimeFE <: AbstractDrillingCost end
(x::Type{T})(a, b) where {T<:AbstractDrillingCost_TimeFE} = T(UnitRange(Float64(a), Float64(b)))
yearrange(x::AbstractDrillingCost_TimeFE) = x.yearrange
@inline start(x::AbstractDrillingCost_TimeFE) = first(yearrange(x))
@inline stop(x::AbstractDrillingCost_TimeFE) = last(yearrange(x))
@inline startstop(x::AbstractDrillingCost_TimeFE) = start(x), stop(x)
@inline time_idx( x::AbstractDrillingCost_TimeFE, t::Number) = Int(clamp(t, start(x), stop(x)) - start(x) + 1)
@inline time_idx(x, obs) = time_idx(x,year(obs))

# -------------------------------------------
# Drilling Cost
# -------------------------------------------

"Single drilling cost"
struct DrillingCost_constant <: AbstractDrillingCost end
@inline _nparm(x::DrillingCost_constant) = 1
@inline function flow!(grad, ::DrillingCost_constant, d, obs, θ, sim, dograd::Bool)
    if dograd
        grad[1] = d
    end
    return d*θ[1]
end
coefnames(::DrillingCost_constant) = ["\\alpha_c",]

"Drilling cost changes if `d>1`"
struct DrillingCost_dgt1 <: AbstractDrillingCost end
@inline _nparm(x::DrillingCost_dgt1) = 2
@inline function flow!(grad, ::DrillingCost_dgt1, d, obs, θ, sim, dograd::Bool)
    T = eltype(θ)
    if d == 0
        dograd && fill!(grad, 0)
        u = zero(T)
    else
        dgt1 = d>1
        u = d*θ[1+dgt1]
        if dograd
            grad[1+dgt1] = d
            grad[2-dgt1] = 0
        end
    end
    return u::T
end
coefnames(::DrillingCost_dgt1) = ["\\alpha_{d=1}", "\\alpha_{d>1}"]

"Time FE for 2008-2012"
struct DrillingCost_TimeFE <: AbstractDrillingCost_TimeFE
    yearrange::UnitRange{Float64}
    # start::Float64
    # stop::Float64
end
@inline _nparm(x::DrillingCost_TimeFE) = 1 + length(yearrange(x))
@inline function flow!(grad, u::DrillingCost_TimeFE, d, obs, θ, sim, dograd::Bool)
    T = eltype(θ)
    dograd && fill!(grad, 0)
    d == 0 && return zero(T)

    tidx = time_idx(u,obs)
    dgt1 = d > 1

    r = d*(θ[tidx] + dgt1*θ[end])
    if dograd
        grad[tidx] = d
        grad[end]  = d*dgt1
    end
    return r::T
end
function coefnames(x::DrillingCost_TimeFE)
    cfs = ["\\alpha_{c,$y}" for y in start(x):stop(x)]
    return vcat(cfs, "\\alpha_{d>1}")
end

"Time FE for 2008-2012 with shifters for (D==0,d>1), (D>1,d==1), (D>1,d>1)"
struct DrillingCost_TimeFE_costdiffs <: AbstractDrillingCost_TimeFE
    yearrange::UnitRange{Float64}
    # start::Float64
    # stop::Float64
end
@inline _nparm(x::DrillingCost_TimeFE_costdiffs) = 3 + length(yearrange(x))
@inline function flow!(grad, u::DrillingCost_TimeFE_costdiffs, d, obs, θ, sim, dograd::Bool)
    T = eltype(θ)
    dograd && fill!(grad, 0)
    d == 0 && return zero(T)

    tidx = time_idx(u,obs)
    Dgt0 = _Dgt0(obs)
    dgt1 = d > 1

    r = d * ( θ[tidx] + (!Dgt0*dgt1)*θ[end-2] + Dgt0*θ[end-1+dgt1] )
    if dograd
        grad[tidx]       = d
        grad[end-2]      = d * !Dgt0*dgt1
        grad[end-1+dgt1] = d * Dgt0
    end
    return r::T
end

function coefnames(x::DrillingCost_TimeFE_costdiffs)
    cfs = ["\\alpha_{c,$y}" for y in start(x):stop(x)]
    return vcat(cfs, "\\alpha_{D=0,d>1}", "\\alpha_{D>0,d=1}", "\\alpha_{D>0,d>1}")
end

"Time FE for 2008-2012 with bonkers flexibility"
struct DrillingCost_DoubleTimeFE <: AbstractDrillingCost_TimeFE
    yearrange::UnitRange{Float64}
end
@inline _nparm(x::DrillingCost_DoubleTimeFE) = 2*(1+length(yearrange(x)))
@inline function flow!(grad, u::DrillingCost_DoubleTimeFE, d, obs, θ, sim, dograd::Bool)
    T = eltype(θ)
    dograd && fill!(grad, 0)
    d == 0 && return zero(T)

    n = length(yearrange(u))
    tidx = time_idx(u,obs)
    Dgt0 = _Dgt0(obs)
    dgt1 = d > 1

    coef_tidx = tidx + Dgt0*n
    dgt1idx = 2*n + 1 + Dgt0

    r = d * ( θ[coef_tidx] + dgt1*θ[dgt1idx] )

    if dograd
        grad[coef_tidx]  = d
        grad[dgt1idx]    = d * dgt1
    end
    return r::T
end

function coefnames(x::DrillingCost_DoubleTimeFE)
    cfs0 = ["\\alpha_{c,D=0,$y}" for y in start(x):stop(x)]
    cfs1 = ["\\alpha_{c,D>0,$y}" for y in start(x):stop(x)]
    return vcat(cfs0, cfs1, "\\alpha_{D=0,d>1}", "\\alpha_{D>0,d>1}")
end

"Time FE w rig rates for 2008-2012"
struct DrillingCost_TimeFE_rigrate <: AbstractDrillingCost_TimeFE
    yearrange::UnitRange{Float64}
    # start::Float64
    # stop::Float64
end
@inline _nparm(x::DrillingCost_TimeFE_rigrate) = 2 + length(yearrange(x))
@inline function flow!(grad, u::DrillingCost_TimeFE_rigrate, d, obs, θ, sim, dograd::Bool)
    T = eltype(θ)
    dograd && fill!(grad, 0)
    d == 0 && return zero(T)

    tidx = time_idx(u,obs)
    rig = rigrate(obs)
    Dgt0 = _Dgt0(obs)
    dgt1 = d > 1

    r = d*(θ[tidx] + dgt1*θ[end-1] + θ[end]*rig)
    if dograd
        grad[tidx]  = d
        grad[end-1] = d*dgt1
        grad[end]   = d*rig
    end
    return r
end

function coefnames(x::DrillingCost_TimeFE_rigrate)
    cfs = ["\\alpha_{c,$y}" for y in start(x):stop(x)]
    return vcat(cfs, "\\alpha_{d>1}", "\\alpha_{rig}")
end

"Time FE for 2008-2012 with shifters for (D==0,d>1), (D>1,d==1), (D>1,d>1)"
struct DrillingCost_TimeFE_rig_costdiffs <: AbstractDrillingCost_TimeFE
    yearrange::UnitRange{Float64}
    # start::Float64
    # stop::Float64
end
@inline _nparm(x::DrillingCost_TimeFE_rig_costdiffs) = 4 + length(yearrange(x))
@inline function flow!(grad, u::DrillingCost_TimeFE_rig_costdiffs, d, obs, θ, sim, dograd::Bool)
    dograd && fill!(grad, 0)
    d == 0 && return zero(T)

    tidx = time_idx(u,obs)
    rig = rigrate(obs)
    Dgt0 = _Dgt0(obs)
    dgt1 = d > 1

    r = d * ( θ[tidx] + (!Dgt0*dgt1)*θ[end-3] + Dgt0*θ[end-2+dgt1] + rig)
    if dograd
        grad[tidx]       = d
        grad[end-3]      = d * (!Dgt0*dgt1)
        grad[end-2+dgt1] = d * Dgt0
        grad[end]        = d * rig
    end
    return r
end

function coefnames(x::DrillingCost_TimeFE_rig_costdiffs)
    cfs = ["\\alpha_{c,$y}" for y in start(x):stop(x)]
    return vcat(cfs, "\\alpha_{D=0,d>1}", "\\alpha_{D>0,d=1}", "\\alpha_{D>0,d>1}", "\\alpha_{rig}")
end


text:
"""
    Functions to copy structures (DataSource, GridStructure)
"""

function copy(ds::DataSource)
    bus = copy(ds.bus)
    link = copy(ds.link)
    return DataSource(bus, link)
end

function copy(gs::GridStructure)
    scenario = gs.scenario
    node_linksin = copy(gs.node_linksin)
    node_linksout = copy(gs.node_linksout)
    return GridStructure(scenario, node_linksin, node_linksout)
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
function UpdatePosition!(sys::System,integrator::NoIntegrator)
end

function UpdateVelocity!(sys::System,integrator::NoIntegrator)
end

function UpdatePosition!(sys::System,integrator::VVIntegrator)
    #sys.coords .+= sys.dt.*sys.vels.+integrator.hdtsqr .* sys.forces ./ sys.masses
    UpdatePosition_particle!.(sys.atoms.coords,sys.atoms.vels,sys.atoms.forces,sys.atoms.masses,sys.dt,integrator)
end

function UpdatePosition_particle!(  coords::Array{Float64,1},
                                    vels::Array{Float64,1},
                                    forces::Array{Float64,1},
                                    masses::Float64,
                                    dt::Float64,
                                    integrator::VVIntegrator)
    coords[1] += dt* vels[1] + integrator.hdtsqr * forces[1] / masses
    coords[2] += dt* vels[2] + integrator.hdtsqr * forces[2] / masses
    coords[3] += dt* vels[3] + integrator.hdtsqr * forces[3] / masses
end

function UpdateVelocity!(sys::System,integrator::VVIntegrator,comm::MPI.Comm)
    #sys.vels .+= 0.5 .* sys.dt .* sys.forces ./ sys.masses
    UpdateVelocity_particle!.(sys.atoms.vels,sys.atoms.forces,sys.atoms.masses,sys.dt,integrator)
end

function UpdateVelocity_particle!(  vels::Array{Float64,1},
                                    forces::Array{Float64,1},
                                    masses::Float64,
                                    dt::Float64,
                                    integrator::VVIntegrator)
    vels[1] += 0.5 * dt * forces[1] / masses
    vels[2] += 0.5 * dt * forces[2] / masses
    vels[3] += 0.5 * dt * forces[3] / masses
end


text:
using MAT
import sbp
datafile=ARGS[1]

vars = matread(datafile)

nx = vars["x"]
p = vars["p"]
v = vars["v"]
qp = vars["q_p"]
qv = vars["q_v"]
Fp = vars["Fp"]
Fv = vars["Fv"]
qFp = vars["q_Fp"]
qFv = vars["q_Fv"]

header = [ "Nx", "p", "v",
           "rate p", "rate v",
          "Fp", "rate Fp",
          "Fv", "rate Fv"]
variables = [
             nx,
             round.(log10.(p), sigdigits=4),
             round.(qp, sigdigits=4),
             round.(log10.(v), sigdigits=4),
             round.(qv, sigdigits=4),
             round.(log10.(Fp), sigdigits=4),
             round.(qFp, sigdigits=4),
             round.(log10.(Fv), sigdigits=4),
             round.(qFv, sigdigits=4),
             ]
print(join(sbp.Table.latex(header, variables)))

max_stars_repo_name: ooreilly/sbp

max_issues_repo_name: ooreilly/sbp

max_forks_repo_name: ooreilly/sbp


text:
#...length of percentage and ETA string with days is 29 characters
tty_width(desc) = max(0, displaysize()[2] - (length(desc) + 29))

function printvalues!(p::AbstractProgress, showvalues; color = false)
    length(showvalues) == 0 && return
    maxwidth = maximum(Int[length(string(name)) for (name, _) in showvalues])
    for (name, value) in showvalues
        msg = "\n  " * rpad(string(name) * ": ", maxwidth+2+1) * string(value)
        (color == false) ? print(p.output, msg) : print_with_color(color, p.output, msg)
    end
    p.numprintedvalues = length(showvalues)
end

function move_cursor_up_while_clearing_lines(io, numlinesup)
    for _ in 1:numlinesup
        print(io, "\u1b[1G\u1b[K\u1b[A")
    end
end

function printover(io::IO, s::AbstractString, color::Symbol = :color_normal)
    if isdefined(Main, :IJulia) || isdefined(Main, :ESS) || isdefined(Main, :Atom)
        print(io, "\r" * s)
    else
        print(io, "\u1b[1G")   # go to first column
        print_with_color(color, io, s)
        print(io, "\u1b[K")    # clear the rest of the line
    end
end

function durationstring(nsec)
    days = div(nsec, 60*60*24)
    r = nsec - 60*60*24*days
    hours = div(r,60*60)
    r = r - 60*60*hours
    minutes = div(r, 60)
    seconds = r - 60*minutes

    hhmmss = @sprintf "%u:%02u:%02u" hours minutes seconds
    if days>9
        return @sprintf "%.2f days" nsec/(60*60*24)
    elseif days>0
        return @sprintf "%u days, %s" days hhmmss
    end
    hhmmss
end


text:
"""
A multilayer perceptron.
"""
module MultilayerPerceptron

include("activation.jl")
include("loss.jl")
include("layer.jl")
include("util.jl")

export MLP, fit!, predict

struct MLP
    layers::Vector{Layer}
end

function MLP(layers::Vector{Tuple{Int, Int}})
    layers = [Layer(sz) for sz in layers]
    MLP(layers)
end

function fit!(
    mlp::MLP,
    train_x::AbstractMatrix{Float64}, train_y::AbstractMatrix,
    ;
    loss_fn::Function = mse,
    learning_rate::Float64 = 0.001,
    batch_size::Int = 1000,
    n_epochs::Int = 10,
    reporting_interval = 1000,
)
    (_, n) = size(train_x)
    batch_iter = Iterators.partition(1:n, batch_size)
    n_batches = 0
    ewma_loss = nothing  # Exponentially weighted moving average
    α = min(8 / reporting_interval, 1.0)  # EWMA α
    for epoch in 1:n_epochs
        for cols in batch_iter
            batch_x = @view train_x[:, cols]
            batch_y = @view train_y[:, cols]
            for l in mlp.layers
                batch_x = forward!(l, batch_x)
            end
            batch_ŷ = batch_x

            (loss, ∂e_over_∂y) = loss_fn(batch_ŷ, batch_y)
            ewma_loss = isnothing(ewma_loss) ? loss : loss * α + ewma_loss * (1.0 - α)
            if (n_batches += 1) % reporting_interval == 0
                println("#$n_batches: $ewma_loss")
            end

            for l in Iterators.reverse(mlp.layers)
                ∂e_over_∂y = backward!(l, ∂e_over_∂y, learning_rate)
            end
        end
    end
end

function predict(mlp::MLP, test_x::AbstractMatrix{Float64})::Matrix{Float64}
    x = test_x
    for l in mlp.layers
        x = forward!(l, x; training=false)
    end
    x
end

end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
"""
@knet Op add(x1,x2; alpha=1,beta=1) performs element-wise broadcasting
addition.  alpha and beta are used to scale x1 and x2 respectively.
The result computed is y = alpha x1 + beta x2.  The size of the output
y always matches the size of x2.  Broadcasting is performed as follows
(following Julia convention, size below refers to the tuple of
dimensions, fastest changing first, and ndims refers to the number of
dimensions):

- If x1 and x2 are the same size their elements are added.

- If x1 and x2 have the same ndims but different sizes, each dimension
  of the x1 array must match the corresponding dimension of the x2
  array or must be equal to 1. In the latter case, the same value from
  the x1 array for those dimensions will be used to blend into the x2
  array.  (At least this is what it says in the NVIDIA cudnnAddTensor
  doc but unfortunately not all combinations work as promised.  I sent
  them a bug report and wrote my own kernels to fill in the gaps for
  now -- Oct 16, 2015.)
  Example: (5,4,1,1)+(5,4,3,2)=>(5,4,3,2)

- If x1 and x2 have different ndims, x1 is assumed to be missing its
  rightmost dims and those are assumed to be 1.  Example:
  (5,4)+(5,4,3,2)=>(5,4,1,1)+(5,4,3,2)=>(5,4,3,2).

- The one exception to the last rule is when x1 has ndims=1 and x2 has
  ndims > 1.  In this case x1 is matched to the next to last dimension
  of x2.  Example: (3,)+(5,4,3,2)=>(1,1,3,1)+(5,4,3,2)=>(5,4,3,2).

"""
type Add <: Op; alpha; beta; Add(;alpha=1,beta=1,o...)=new(alpha,beta); end

ninputs(::Add)=2
canoverwrite(::Add)=true
back_reads_x(::Add)=false
back_reads_y(::Add)=false

function forw(a::Add, x1, x2, y; o...)
    @assert x2 == nothing || size(y) == size(x2)
    if x1!=nothing && x2!=nothing # we use nothing to represent the zero array
        baddforw!(a.alpha,x1,a.beta,x2,y)
    elseif x2 != nothing
        y===x2 || copysync!(y, x2)
        a.beta==1 || scale!(y, a.beta)
    elseif x1 != nothing && size(y)==size(x1)
        y===x1 || copysync!(y, x1)
        a.alpha==1 || scale!(y, a.alpha)
    else
        y=nothing
    end
    gpusync()
    return y
end

function back(a::Add, dy, dx1, dx2; o...)
    if dx2 != nothing
        size(dx2) == size(dy) || throw(DimensionMismatch("The size of the output must match x2"))
        dx2 === dy  || copysync!(dx2, dy)
        a.beta == 1 || scale!(a.beta, dx2)
    end
    if dx1 != nothing
        baddback!(dy, dx1)
        a.alpha == 1 || scale!(a.alpha, dx1)
    end
    gpusync()
end

addforw0!(alpha,a,beta,b,c)=(a===c||b===c||copysync!(c,b))
baddforw0!(alpha,a,beta,b,c)=(b===c||copysync!(c,b))
baddback0!(alpha,dy,db)=fillsync!(db,0)

function baddforw!(alpha,a,beta,b,c)
    size(c) == size(b) || throw(DimensionMismatch("The size of the output must match the second input"))
    size(a) == size(b) || (a = reshape_to_match(a,b; CUDNN_ADD_SAME_C=true))
    if size(a)==size(b)
        addforw3!(alpha,a,beta,b,c)
    elseif cudnnAddTensorCompatible(a,b)
        baddforw1!(alpha,a,beta,b,c)
    else
        baddforw2!(alpha,a,beta,b,c)
    end
end

@gpu function addforw1!{T}(alpha::Number,a::CudaArray{T},beta::Number,b::CudaArray{T},c::CudaArray{T})
    c===b || copysync!(c,b)
    cudnnAddTensor(a,c; alpha=T(alpha), beta=T(beta))
    gpusync(); return c
end

@gpu function addforw2!{T}(alpha::Number,a::CudaArray{T},beta::Number,b::CudaArray{T},c::CudaArray{T})
    CUBLAS.geam!('N','N',T(alpha),a,T(beta),b,c)
    gpusync(); return c
end

function addforw3!(alpha,a,beta,b,c)
    if c===b
        beta != 1 && scale!(beta, c)
        axpy!(alpha, a, c)
    elseif c===a
        alpha != 1 && scale!(alpha, c)
        axpy!(beta, b, c)
    else
        copysync!(c,b)
        beta != 1 && scale!(beta, c)
        axpy!(alpha, a, c)
    end
    gpusync(); return c
end

function baddforw!(alpha::Number,a::Array,beta::Number,b::Array,c::Array)
    broadcast!(+,c,(alpha==1 ? a : alpha*a),(beta==1 ? b : beta*b))
end

@gpu function baddforw1!{T}(alpha::Number,a::CudaArray{T},beta::Number,b::CudaArray{T},c::CudaArray{T}) # mnist2d:6.40
    c===b || copysync!(c,b)
    cudnnAddTensor(a,c; alpha=T(alpha), beta=T(beta))
    gpusync(); return c
end

function cudnnAddTensorCompatible(a,b)
    sa = size(a)
    sb = size(b)
    sa == sb && return true
    n = length(sa)
    length(sb) == n || return false
    sa == ntuple(i->(i==n-1 ? sb[i] : 1), n) && return true
    sa == ntuple(i->(i==n ? 1 : sb[i]), n) && return true
    sa == ntuple(i->(i>=n-1 ? 1 : sb[i]), n) && return true
    return false
end

@gpu function baddforw2!{T}(alpha::Number,a::CudaArray{T},beta::Number,b::CudaArray{T},c::CudaArray{T}) # mnist2d:7.55
    for i=1:ndims(a)
        size(a,i)==1 || size(a,i)==size(b,i) ||
        throw(DimensionMismatch("Each dimension of x1 must match x2 or be 1."))
    end
    ndims(b) <= 8 || error("add kernel supports dimensions up to 8")
    T <: Float32 ? ccall((:addforw32,libknet),Void,(Cint,Cfloat,Ptr{Cint},Ptr{Cfloat},Cfloat,Ptr{Cint},Ptr{Cfloat},Ptr{Cfloat}),ndims(a),T(alpha),cudadims(a),a,T(beta),cudadims(b),b,c) :
    T <: Float64 ? ccall((:addforw64,libknet),Void,(Cint,Cdouble,Ptr{Cint},Ptr{Cdouble},Cdouble,Ptr{Cint},Ptr{Cdouble},Ptr{Cdouble}),ndims(a),T(alpha),cudadims(a),a,T(beta),cudadims(b),b,c) :
    error("$T not supported")
    gpusync(); return c
end

# TODO: this does not cover all forms of db for cpu:
function baddback!(dy::Array, db::Vector)
    c=ndims(dy)-1
    fillsync!(db, zero(eltype(db)))
    @inbounds for i=1:length(dy)
        db[ind2sub(size(dy),i)[c]] += dy[i]
    end
    return db
end

function baddback!(dy, db)
    size(db) == size(dy) || (db = reshape_to_match(db,dy; CUDNN_ADD_SAME_C=true))
    if size(db) == size(dy)
        db === dy   || copysync!(db, dy)
    elseif cudnnConvolutionBackwardBiasCompatible(dy,db)
        baddback1!(dy, db)
    elseif ndims(dy) == 2
        baddback2!(dy, db)
    else
        baddback3!(dy, db)
    end
end

function cudnnConvolutionBackwardBiasCompatible(dy,db)
    ndims(dy)==ndims(db)==1 && return (size(dy)==size(db))
    size(db) == ntuple(i->(i==ndims(dy)-1 ? size(dy,i) : 1), ndims(dy))
end

@gpu function baddback1!{T}(dy::CudaArray{T}, db::CudaArray{T})
    cudnnConvolutionBackwardBias(dy, db)
    gpusync(); return db
end

@gpu function baddback2!{T}(dy::CudaArray{T}, db::CudaArray{T})
    ndims(db) == ndims(dy) || throw(DimensionMismatch())
    ndims(db)==2 || throw(DimensionMismatch())
    if size(db,1)==size(dy,1) && size(db,2)==1
        tmp = fillsync!(similar(dy, (size(dy,2),1)),1)
        A_mul_B!(db,dy,tmp)
        free(tmp)
    elseif size(db,2)==size(dy,2) && size(db,1)==1
        tmp = fillsync!(similar(dy, (1,size(dy,1))),1)
        A_mul_B!(db,tmp,dy)
        free(tmp)
    elseif size(db)==(1,1)
        baddback3!(dy,db)
    else
        throw(DimensionMismatch())
    end
    gpusync(); return db
end

@gpu function baddback3!{T}(dy::CudaArray{T}, db::CudaArray{T})
    ndims(db) == ndims(dy) || throw(DimensionMismatch())
    ndims(db) <= 8 || error("add kernel supports dimensions up to 8")
    for i=1:ndims(db)
        size(db,i)==1 || size(db,i)==size(dy,i) ||
        throw(DimensionMismatch("Each dimension of x1 must match x2 or be 1."))
    end
    fillsync!(db,0)
    T <: Float32 ? ccall((:addback32,libknet),Void,(Cint,Ptr{Cint},Ptr{Cfloat},Ptr{Cint},Ptr{Cfloat}),ndims(dy),cudadims(dy),dy,cudadims(db),db) :
    T <: Float64 ? ccall((:addback64,libknet),Void,(Cint,Ptr{Cint},Ptr{Cdouble},Ptr{Cint},Ptr{Cdouble}),ndims(dy),cudadims(dy),dy,cudadims(db),db) :
    error("$T not supported")
    gpusync(); return db
end

"reshape_to_match(x1,x2) reshapes x1 to match x2."
function reshape_to_match(x1,x2; CUDNN_ADD_SAME_C=false)
    if CUDNN_ADD_SAME_C && ndims(x1)==1
        newsize = ntuple(ndims(x2)) do i
            i != ndims(x2)-1 ? 1 :
            size(x2,i) == size(x1,1) ? size(x1,1) :
            throw(DimensionMismatch("$(size(x1)) $(size(x2))"))
        end
    else
        newsize = ntuple(ndims(x2)) do i
            i > ndims(x1) ? 1 :
            size(x1,i)==1 ? 1 :
            size(x1,i)==size(x2,i) ? size(x1,i) :
            throw(DimensionMismatch("$(size(x1)) $(size(x2))"))
        end
    end
    reshape(x1, newsize)
end

function infersize(a::Add, x1, x2, y)
    if x1==x2==y==nothing
        nothing
    elseif x2==y==nothing
        (x1, nothing, nothing)
    elseif x2==nothing
        infersize(a, x1, y, y)
    elseif y==nothing
        infersize(a, x1, x2, x2)
    else
        ydims = map(x2, y) do xdim, ydim
            xdim == ydim ? xdim :
            xdim == 0 ? ydim :
            ydim == 0 ? xdim :
            throw(DimensionMismatch("$y=add($x1,$x2)"))
        end
        x1==nothing && (return (ydims, ydims, ydims))
        xdims = [x1...]; ydims = [ydims...]
        if length(xdims) == 1
            xdims[1] == ydims[end-1] ? nothing :
            xdims[1] == 0 ? xdims[1]=ydims[end-1] :
            ydims[end-1]==0 ? ydims[end-1]=xdims[1] :
            throw(DimensionMismatch("$y=add($x1,$x2)"))
        elseif length(x1) <= length(ydims)
            for i=1:length(xdims)
                xdims[i] == ydims[i] ? continue :
                xdims[i] == 1 ? continue :
                xdims[i] == 0 ? continue :
                ydims[i] == 0 ? (ydims[i] = xdims[i]) :
                throw(DimensionMismatch("$y=add($x1,$x2)"))
            end
        else
            throw(DimensionMismatch("$y=add($x1,$x2)"))
        end
        xdims = tuple(xdims...); ydims = tuple(ydims...)
        return (xdims, ydims, ydims)
    end
end

### DEAD CODE:

# biasforw(b::Vector, x::Array, y::Array)=(c=ndims(x)-1; for i=1:length(y); y[i] = x[i] + b[ind2sub(size(x),i)[c]]; end; y)
# biasforw(b::Vector, x::Vector, y::Vector)=(for i=1:length(y); y[i] = x[i] + b[i]; end; y)
# @gpu biasforw(b::CudaArray, x::CudaArray, y::CudaArray)=(y===x||copysync!(y,x);cudnnAddTensor(b, y; mode=CUDNN_ADD_SAME_C); gpusync(); y)

# function sizeafterbias(x1,x2)
#     i1 = x1[1]
#     i2 = x2[end-1]
#     i1 == 0 && (i1=i2)
#     i2 == 0 && (i2=i1)
#     i1 == i2 || error()
#     tuple(x2[1:end-2]..., i2, x2[end])
# end

# function commonsize(x1,x2)
#     map(x1, x2) do i1,i2
#         i1 == 0 && (i1=i2)
#         i2 == 0 && (i2=i1)
#         i1 == i2 || error()
#         i1
#     end
# end

    # elseif biassize(dy, dx1)
    #     biasback(dy, dx1)
    #     a.alpha == 1 || scale!(a.alpha, dx1)
    # else
    #     # TODO: implement more broadcasting back.
    #     error("Don't know how to do back pass with $(size(dy))=$(size(dx1))+$(size(dx2))")
# biassize(dy,db)=(size(db,1)==size(dy, ndims(dy)==1 ? 1 : ndims(dy)-1) && all([size(db,i)==1 for i=2:ndims(db)]))

# baddback!(alpha::Number, dy::Array, db::Vector)=(c=ndims(dy)-1; fillsync!(db, zero(eltype(db))); for i=1:length(dy); db[ind2sub(size(dy),i)[c]] += dy[i]; end; alpha==1||scale!(alpha,db))
# baddback!(alpha::Number, dy::Vector, db::Vector)=(for i=1:length(dy); db[i]=dy[i]; end; alpha==1||scale!(alpha,db))

# TODO: averaging? keep both arr and avg inside par, set output to one without copy?
# should be handled between net and par not by the op.
# DONE: handle nothings?  -- net is handling them.
# DONE: handle scalar input for adding a constant. -- axpb will handle this
# DONE: back

# add(x1,x2,y; alpha=1, beta=1)=(Add(alpha,beta),x1,x2,y)
# _KENV[:add] = Add
# _KENV[:+] = Add


text:
# code for node updating via optimization

typealias explicit_opts Union{Val{:l_bfgs}, Val{:cg}}
function update!{D, S <: explicit_opts}(n::RandomNode{D}, m::VBModel, ::Type{S})
    # use const here so that the closure below can be optimized
    const fac_list = [f for (f, _) in m.graph[n]]

    const x0 = unroll_pars(n)
    # TODO: get lower and upper bounds from each variable in x from n

    # make an objective function that sets the parameters of n to x
    # and sums the values of all factors containing n
    function objfun(x)
        update_pars!(n, x)
        val = 0.0
        for f in fac_list
            val += value(f)
        end
        -val  # Optim minimizes, so objfun = -ELBO
    end

    # use ForwardDiff to get the gradient; autodiff in optimize uses
    # ReverseDiff, but ForwardDiff only requires args <: Real, but
    # ReverseDiff needs args <: Number, which Distributions doesn't allow
    gradf! = (out, x) -> ForwardDiff.gradient!(out, objfun, x)

    # define mutating gradient; store gradient in storage array
    objgrad!(x, storage) = gradf!(storage, x)

    if S.parameters[1] == :l_bfgs
        method = LBFGS()
    elseif S.parameters[1] == :cg
        method = ConjugateGradient()
    end

    # try optimization; if it fails, set parameters back to initial guess
    try
        optimize(objfun, objgrad!, x0, method=method)
    catch
        finalpars = unroll_pars(n)
        println("Optimization failed at pars:\n$finalpars")
        update_pars!(n, x0)
    end
end

function update_pars!(n::RandomNode, x)
    par_sizes = get_par_sizes(n[1])
    npars = mapreduce(prod, +, par_sizes)

    ctr = 0
    for i in eachindex(n.data)
        n[i] = reroll_pars(n[i], par_sizes, x[ctr + 1: ctr + npars])
        ctr += npars
    end
end


text:
# Based on https://github.com/JuliaLang/julia/blob/master/base/complex.jl (MIT license)

import Base: real, imag, conj, abs, +, -, *, /, exp, convert, promote_rule, isfinite

struct ComplexNumber{T<:Real} <: Number
    re::T
    im::T
end

const jm = ComplexNumber(0, 1)

ComplexNumber(re::Real, im::Real) = ComplexNumber(promote(re, im)...)
ComplexNumber(re::Real) = ComplexNumber(re, zero(re))

convert(::Type{ComplexNumber{T}}, x::Real) where {T<:Real} = ComplexNumber{T}(x, 0)
convert(::Type{ComplexNumber{T}}, z::ComplexNumber) where {T<:Real} = ComplexNumber{T}(real(z), imag(z))

promote_rule(::Type{ComplexNumber{T}}, ::Type{S}) where {T<:Real, S<:Real} = ComplexNumber{promote_type(T, S)}
promote_rule(::Type{ComplexNumber{T}}, ::Type{ComplexNumber{S}}) where {T<:Real, S<:Real} = ComplexNumber{promote_type(T, S)}

real(z::ComplexNumber) = z.re
imag(z::ComplexNumber) = z.im

reim(z::ComplexNumber) = z.re, z.im

conj(z::ComplexNumber) = ComplexNumber(real(z), -imag(z))
abs(z::ComplexNumber)  = hypot(real(z), imag(z))

isfinite(z::ComplexNumber) = isfinite(real(z)) && isfinite(imag(z))

+(u::ComplexNumber, v::ComplexNumber) = ComplexNumber(real(u) + real(v), imag(u) + imag(v))
-(u::ComplexNumber, v::ComplexNumber) = ComplexNumber(real(u) - real(v), imag(u) - imag(v))
*(u::ComplexNumber, v::ComplexNumber) = ComplexNumber(real(u) * real(v) - imag(u) * imag(v), real(u) * imag(v) + imag(u) * real(v))
/(u::ComplexNumber, v::ComplexNumber) = ComplexNumber((real(u) * real(v) + imag(u) * imag(v)) / (real(v)^2 + imag(v)^2), (imag(u) * real(v) - real(u) * imag(v)) / (real(v)^2 + imag(v)^2))
-(z::ComplexNumber) = ComplexNumber(-real(z), -imag(z))

function exp(z::ComplexNumber)
    z_re, z_im = reim(z)
    e_re = exp(z_re)
    ComplexNumber(e_re * cos(z_im), e_re * sin(z_im))
end


text:
using Revise, Parameters, Logging, Test
using CUDA
# include("/home/jakub/Documents/NuclearMedEval/src/utils/CUDAAtomicUtils.jl")

# includet("/home/jakub/Documents/NuclearMedEval/src/utils/BitWiseUtils.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/structs/BasicStructs.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/utils/CUDAGpuUtils.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/utils/IterationUtils.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/utils/ReductionUtils.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/utils/MemoryUtils.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/distanceMetrics/Housdorff/verB/MetaDataUtils.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/distanceMetrics/Housdorff/verB/ResultListUtils.jl")

# includet("/home/jakub/Documents/NuclearMedEval/src/overLap/MainOverlap.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/PairCounting/RandIndex.jl")

# includet("/home/jakub/Documents/NuclearMedEval/src/Probabilistic/ProbabilisticMetrics.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/InformationTheorhetic/InformationTheorhetic.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/volume/VolumeMetric.jl")

# includet("/home/jakub/Documents/NuclearMedEval/src/distanceMetrics/MeansMahalinobis.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/distanceMetrics/Housdorff/verB/HFUtils.jl")

# includet("/home/jakub/Documents/NuclearMedEval/src/distanceMetrics/Housdorff/verB/PrepareArrtoBool.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/distanceMetrics/Housdorff/verB/WorkQueueUtils.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/distanceMetrics/Housdorff/verB/ProcessMainDataVerB.jl")
# includet("/home/jakub/Documents/NuclearMedEval/test/GPUtestUtils.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/distanceMetrics/Housdorff/verB/ScanForDuplicates.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/distanceMetrics/Housdorff/verB/MetadataAnalyzePass.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/distanceMetrics/Housdorff/verB/MainLoopKernel.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/distanceMetrics/Housdorff/verB/Housdorff.jl")

# includet("/home/jakub/Documents/NuclearMedEval/src/kernels/TpfpfnKernel.jl")
# includet("/home/jakub/Documents/NuclearMedEval/src/kernels/InterClassCorrKernel.jl")

includet("C:/GitHub/GitHub/NuclearMedEval/src/utils/CUDAAtomicUtils.jl")

includet("C:/GitHub/GitHub/NuclearMedEval/src/utils/BitWiseUtils.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/structs/BasicStructs.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/utils/CUDAGpuUtils.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/utils/IterationUtils.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/utils/ReductionUtils.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/utils/MemoryUtils.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/Housdorff/verB/MetaDataUtils.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/Housdorff/verB/ResultListUtils.jl")

includet("C:/GitHub/GitHub/NuclearMedEval/src/overLap/MainOverlap.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/PairCounting/RandIndex.jl")

includet("C:/GitHub/GitHub/NuclearMedEval/src/Probabilistic/ProbabilisticMetrics.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/InformationTheorhetic/InformationTheorhetic.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/volume/VolumeMetric.jl")

includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/MeansMahalinobis.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/Housdorff/verB/HFUtils.jl")

includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/Housdorff/verB/PrepareArrtoBool.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/Housdorff/verB/WorkQueueUtils.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/Housdorff/verB/ProcessMainDataVerB.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/test/GPUtestUtils.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/Housdorff/verB/ScanForDuplicates.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/Housdorff/verB/MetadataAnalyzePass.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/Housdorff/verB/MainLoopKernel.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/distanceMetrics/Housdorff/verB/Housdorff.jl")

includet("C:/GitHub/GitHub/NuclearMedEval/src/kernels/TpfpfnKernel.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/kernels/InterClassCorrKernel.jl")
includet("C:/GitHub/GitHub/NuclearMedEval/src/MainAbstractions.jl")


text:
using CoordinateTransformations, Rotations, TestImages, ImageCore, StaticArrays, OffsetArrays, Interpolations, LinearAlgebra
using Test, ReferenceTests

refambs = detect_ambiguities(CoordinateTransformations, Base, Core)
using ImageTransformations
ambs = detect_ambiguities(ImageTransformations, CoordinateTransformations, Base, Core)
@test isempty(setdiff(ambs, refambs))

function typestring(::Type{T}) where T   # from https://github.com/JuliaImages/ImageCore.jl/pull/133
    buf = IOBuffer()
    show(buf, T)
    String(take!(buf))
end

# helper function to compare NaN
nearlysame(x, y) = x ≈ y || (isnan(x) & isnan(y))
nearlysame(A::AbstractArray, B::AbstractArray) = all(map(nearlysame, A, B))

tests = [
    "autorange.jl",
    "resizing.jl",
    "interpolations.jl",
    "warp.jl",
    "deprecated.jl" # test deprecations in the last
]

@testset "ImageTransformations" begin
for t in tests
    @testset "$t" begin
        include(t)
    end
end
end

nothing


text:
const MAX_IMPLEMENTATION_NAME = 64
const MAX_KERNEL_NAME = 256
const MAX_LOG_MESSAGE_LEN = 1024
const MAX_REFERENCE_NAME = 64
const VENDOR_MASK = 0xfff00000
const TYPE_MASK = 0x000fff00
const LIBRARY_MASK = 0x000ff000
const KERNEL_MASK = 0x00000fff
const ATTRIBUTE_ID_MASK = 0x000000ff
const ENUM_TYPE_MASK = 0x000ff000
const ENUM_MASK = 0x000000fff
const FMT_REF = "%p"
const FMT_SIZE = "%zu"
const SCALE_UNITY = UInt32(1024)
const SCALE_PYRAMID_HALF = Float32(0.5)
const SCALE_PYRAMID_ORB = Float32(0.8408964)
const OPENVX_KHR_NN = "vx_khr_nn"
const LIBRARY_KHR_NN_EXTENSION = 0x01
const AGO_TARGET_AFFINITY_CPU = 0x0010
const AGO_TARGET_AFFINITY_GPU = 0x0020
const AGO_MAX_PARAMS = 32
const AGO_MERGE_RULE_MAX_FIND = 4
const AGO_MERGE_RULE_MAX_REPLACE = 4
const AGO_MERGE_RULE_SOLITARY_FLAG = 0x20
const AGO_TARGET_AFFINITY_GPU_INFO_DEVICE_MASK = Float32(0x00)
const AGO_TARGET_AFFINITY_GPU_INFO_SVM_MASK = 0xf0
const AGO_TARGET_AFFINITY_GPU_INFO_SVM_ENABLE = 0x10
const AGO_TARGET_AFFINITY_GPU_INFO_SVM_AS_CLMEM = 0x20
const AGO_TARGET_AFFINITY_GPU_INFO_SVM_NO_FGS = 0x40
const MAX_STRING_BUFFER_SIZE_AMD = 256
const OPENVX_KHR_XML = "vx_khr_xml"
const OPENVX_KHR_IX = "vx_khr_ix"
const OPENVX_KHR_TILING = "vx_khr_tiling"
const MAX_TILING_PLANES = 4
const ENUM_IX_USE = 0x18

## Macro-dependent constants

const VERSION_1_0 = VERSION_MAJOR(1) | VERSION_MINOR(0)
const VERSION_1_1 = VERSION_MAJOR(1) | VERSION_MINOR(1)
const VERSION_1_2 = VERSION_MAJOR(1) | VERSION_MINOR(2)
const VERSION = VERSION_1_2

const KERNEL_ATTRIBUTE_LOCAL_DATA_PTR = ATTRIBUTE_BASE(ID_KHRONOS , TYPE_KERNEL) + 0x4
const GRAPH_ATTRIBUTE_STATUS = ATTRIBUTE_BASE(ID_KHRONOS , TYPE_GRAPH) + 0x1
const IMAGE_SIZE = ATTRIBUTE_BASE(ID_KHRONOS , TYPE_IMAGE) + 0x6
const META_FORMAT_ATTRIBUTE_DELTA_RECTANGLE = ATTRIBUTE_BASE(ID_KHRONOS , TYPE_META_FORMAT) + 0x0
const HINT_SERIALIZE = ENUM_BASE(ID_KHRONOS , ENUM_HINT) + 0x0
const KERNEL_INVALID = KERNEL_BASE(ID_KHRONOS , LIBRARY_KHR_BASE) + 0x0
const THRESHOLD_THRESHOLD_VALUE = ATTRIBUTE_BASE(ID_KHRONOS , TYPE_THRESHOLD) + 0x1
const THRESHOLD_THRESHOLD_LOWER = ATTRIBUTE_BASE(ID_KHRONOS , TYPE_THRESHOLD) + 0x2
const THRESHOLD_THRESHOLD_UPPER = ATTRIBUTE_BASE(ID_KHRONOS , TYPE_THRESHOLD) + 0x3
const THRESHOLD_TRUE_VALUE = ATTRIBUTE_BASE(ID_KHRONOS , TYPE_THRESHOLD) + 0x4
const THRESHOLD_FALSE_VALUE = ATTRIBUTE_BASE(ID_KHRONOS , TYPE_THRESHOLD) + 0x5
const THRESHOLD_DATA_TYPE = ATTRIBUTE_BASE(ID_KHRONOS , TYPE_THRESHOLD) + 0x6
const NN_ACTIVATION_LEAKY_RELU = ENUM_BASE(ID_AMD , ENUM_NN_ACTIVATION_FUNCTION_TYPE) + 0x9
const IX_USE_APPLICATION_CREATE = ENUM_BASE(ID_KHRONOS , ENUM_IX_USE) + 0x0
const IX_USE_EXPORT_VALUES = ENUM_BASE(ID_KHRONOS , ENUM_IX_USE) + 0x1
const IX_USE_NO_EXPORT_VALUES = ENUM_BASE(ID_KHRONOS , ENUM_IX_USE) + 0x2

## Extra macros -- implement as needed/skipped for now

# vxFormatArrayPointer ( ptr , index , stride ) ( & ( ( ( vx_uint8 * ) ( ptr ) ) [ ( index ) * ( stride ) ] ) )
# vxArrayItem ( type , ptr , index , stride ) ( * ( type * ) ( vxFormatArrayPointer ( ( ptr ) , ( index ) , ( stride ) ) ) )
# vxImageHeight ( ptile ) ( ( ptile ) ) -> image . height )
# vxImageWidth ( ptile ) ( ( ptile ) ) -> image . width )
# vxTileX ( ptile ) ( ( ptile ) -> tile_x )
# vxTileY ( ptile ) ( ( ptile ) -> tile_y )
# vxTileWidth ( ptile , index ) ( ( ptile ) -> addr [ index ] . dim_x )
# vxTileHeight ( ptile , index ) ( ( ptile ) -> addr [ index ] . dim_y )
# vxTileBlockHeight ( ptile ) ( ( ptile ) -> tile_block . height )
# vxTileBlockWidth ( ptile ) ( ( ptile ) -> tile_block . width )
# vxNeighborhoodLeft ( ptile ) ( ( ptile ) -> neighborhood . left )
# vxNeighborhoodRight ( ptile ) ( ( ptile ) -> neighborhood . right )
# vxNeighborhoodTop ( ptile ) ( ( ptile ) -> neighborhood . top )
# vxNeighborhoodBottom ( ptile ) ( ( ptile ) -> neighborhood . bottom )
# vxImageOffset ( ptile , i , x , y , ox , oy ) ( ( ptile ) -> addr [ i ] . stride_y * ( vx_int32 ) ( ( ( vx_int32 ) ( ( oy ) + ( y ) ) * ( vx_int32 ) ( ptile ) -> addr [ i ] . scale_y ) / ( vx_int32 ) SCALE_UNITY ) ) + ( ( ptile ) -> addr [ i ] . stride_x * ( vx_int32 ) ( ( ( vx_int32 ) ( ( ox ) + ( x ) ) * ( vx_int32 ) ( ptile ) -> addr [ i ] . scale_x ) / ( vx_int32 ) SCALE_UNITY ) )
# vxImagePixel ( type , ptile , i , x , y , ox , oy ) * ( ( type * ) ( & ( ( vx_uint8 * ) ( ptile ) -> base [ i ] ) [ vxImageOffset ( ptile , i , x , y , ox , oy ) ] ) )


text:
const AbstractScalarOrVec{T} = Union{T,AbstractVector{T}} where {T<:AbstractFloat}

# Support of passing a vector of RNGs

Base.rand(rng::AbstractVector{<:AbstractRNG}) = rand.(rng)

Base.randn(rng::AbstractVector{<:AbstractRNG}) = randn.(rng)

function Base.rand(rng::AbstractVector{<:AbstractRNG}, T, n_chains::Int)
    @argcheck length(rng) == n_chains
    return rand.(rng, T)
end

function Base.randn(rng::AbstractVector{<:AbstractRNG}, T, dim::Int, n_chains::Int)
    @argcheck length(rng) == n_chains
    return cat(randn.(rng, T, dim)...; dims=2)
end

# Sample from Categorical distributions

randcat_logp(rng::AbstractRNG, unnorm_ℓp::AbstractVector) =
    randcat(rng, exp.(unnorm_ℓp .- logsumexp(unnorm_ℓp)))

function randcat(rng::AbstractRNG, p::AbstractVector{T}) where {T}
    u = rand(rng, T)
    c = zero(eltype(p))
    i = 0
    while c < u
        c += p[i+=1]
    end
    return max(i, 1)
end

randcat_logp(
    rng::Union{AbstractRNG, AbstractVector{<:AbstractRNG}},
    unnorm_ℓP::AbstractMatrix
) = randcat(rng, exp.(unnorm_ℓP .- logsumexp(unnorm_ℓP; dims=2)))

function randcat(
    rng::Union{AbstractRNG, AbstractVector{<:AbstractRNG}},
    P::AbstractMatrix{T}
) where {T}
    u = rand(rng, T, size(P, 1))
    C = cumsum(P; dims=2)
    is = convert.(Int, vec(sum(C .< u; dims=2)))
    return max.(is, 1)
end