text:
module ColorMaps

export ColorMap, Gray, Brew, RGBArray, Distinguishable, SparseDistinguishable, write
import Images: grayim, save, ImageCmap
import Colors: RGB, distinguishable_colors, colormap
import ColorBrewer: palette

abstract ColorMap

type Gray <: ColorMap
    invert::Bool
    Gray(;invert = false) = new(invert)
end

type RGBArray <: ColorMap
    colors::Array{RGB{Float64},1}
    function RGBArray(colors; invert=false)
        if invert
            colors = flipdim(colors, 1)
        end
        new(colors)
    end
end

function Brew(name::AbstractString, number::Integer; invert = false)
    a = RGB{Float64}[convert(RGB{Float64}, c) for c in palette(name, number)]
    RGBArray(a, invert=invert)
end

function Distinguishable(n::Integer; invert=false)
    RGBArray(convert(Array{RGB{Float64},1}, distinguishable_colors(n)), invert=invert)
end

function SparseDistinguishable(n, nonzeroidx; zerocolor = RGB{Float64}(1.,1.,1.))
    m = RGB{Float64}[zerocolor for i = 1:n]
    sm = ColorMaps.Distinguishable(length(nonzeroidx))
    for i = 1:length(sm.colors)
        m[nonzeroidx[i]] = sm.colors[i]
    end
    ColorMaps.RGBArray(m)
end

function Named(name::AbstractString = "Jet", levels=255)
    if isequal(name, "Jet")
        return RGBArray([RGB(min(max(min(4*i/256-1.5,-4*i/256+4.5),0),1), min(max(min(4*i/256-0.5,-4*i/256+3.5),0),1),min(max(min(4*i/256+0.5,-4*i/256+2.5),0),1)) for i = 1:255])
    else
        return RGBArray(colormap(name, levels))
    end
end

function Base.write(colormap::ColorMaps.Gray, data, filename)
    if colormap.invert
        inverteddata = 1. - data
        save(filename, grayim(rotr90(flipdim(inverteddata, 1))))
    else
        save(filename, grayim(rotr90(flipdim(data, 1))))
    end
end

Base.write(colormap::ColorMap, data, filename) = error("Not supported")

function Base.write(colormap::ColorMaps.RGBArray, data, filename)
    img = ImageCmap(round(UInt8, 1.+(length(colormap.colors)-1).*(data)), colormap.colors)
    save(filename, img)
end

end


text:
type LRN <: Op; n; alpha; beta; k; end
#LRN(x,y;n=5,alpha=1e-4,beta=0.75,k=2.0,o...)=(LRN(n,alpha,beta,k),x,y)
LRN(x,y;o...)=error(:LRN_NOT_TESTED)
ninputs(::LRN)=1
overwrites(::LRN)=false
back_reads_x(::LRN)=true
back_reads_y(::LRN)=true
infersize(::LRN,xdims,ydims)=infersize(Sigm(),xdims,ydims)
forw(l::LRN, x, y; o...)=
    (cudnnLRNCrossChannelForward(x,y; n=l.n, alpha=l.alpha, beta=l.beta, k=l.k); gpusync(); y)
back(l::LRN, dy, dx; x=nothing, y=nothing, o...)=
    (dx!=nothing && cudnnLRNCrossChannelBackward(y,dy,x,dx;n=l.n, alpha=l.alpha, beta=l.beta, k=l.k); gpusync(); dx)

ext: jl


text:
module TransformVariables

using ArgCheck: @argcheck
using DocStringExtensions: SIGNATURES, TYPEDEF
import ForwardDiff
using LinearAlgebra: UpperTriangular, logabsdet
using MacroTools: @capture
using Parameters: @unpack
using Random: AbstractRNG, GLOBAL_RNG

include("utilities.jl")
include("generic.jl")
include("scalar.jl")
include("special_arrays.jl")
include("aggregation.jl")
include("custom.jl")

end # module


text:
using DagScheduler
using FoundationDB
using Base.Test

start_client()
atexit() do
    stop_client()
end
yield()

const waits = Float64[]
macro until(cond)
    quote
        t1 = time()
        while !($(esc(cond)))
            #yield()
            sleep(0.05)
        end
        push!(waits, time()-t1)
    end
end

@testset "fdb queue" begin
    open(FDBCluster()) do cluster
        open(FDBDatabase(cluster)) do db
            ncallbacks = 0
            sharemode = DagScheduler.ShareMode(1)
            ts = DagScheduler.FdbMeta.FdbTaskStore(db, "A", sharemode)
            DagScheduler.FdbMeta.set_callback(ts, (x)->(ncallbacks += 1))

            DagScheduler.FdbMeta.clear(ts)
            DagScheduler.FdbMeta.init(ts)
            DagScheduler.FdbMeta.start_processing_events(ts)

            @test isempty(ts.taskids)
            @test isempty(ts.taskprops)

            test_tids = DagScheduler.TaskIdType[1, 2, 3, 4, 5]
            test_annotated_tids = DagScheduler.TaskIdType[11, 12, 13, 14, 15]

            DagScheduler.FdbMeta.new_task(ts, test_tids[1], test_annotated_tids[1])
            @until (length(ts.taskids) == 1) && (sharemode.ncreated == 1)
            @test length(ts.taskids) == 1
            @test length(ts.taskprops) == 1
            @test ts.taskids[1] == test_tids[1]
            key, tid_annotated, reservation = ts.taskprops[test_tids[1]]
            @test tid_annotated == test_annotated_tids[1]
            @test reservation == 0
            @test sharemode.ncreated == 1
            @test sharemode.ndeleted == 0
            @test sharemode.nshared == 1

            oldkey = key
            DagScheduler.FdbMeta.new_task(ts, test_tids[1], test_annotated_tids[1])
            sleep(0.3)
            key, tid_annotated, reservation = ts.taskprops[test_tids[1]]
            @test key == oldkey

            DagScheduler.FdbMeta.reserve_task(ts, test_tids[1], 3)
            @until length(ts.taskids) == 0
            @test length(ts.taskids) == 0
            @test length(ts.taskprops) == 1
            key, tid_annotated, reservation = ts.taskprops[test_tids[1]]
            @test tid_annotated == test_annotated_tids[1]
            @test reservation == 3
            @test sharemode.ncreated == 1
            @test sharemode.ndeleted == 1
            @test sharemode.nshared == 0

            oldkey = key
            DagScheduler.FdbMeta.new_task(ts, test_tids[1], test_annotated_tids[1])
            sleep(0.3)
            key, tid_annotated, reservation = ts.taskprops[test_tids[1]]
            @test key == oldkey

            DagScheduler.FdbMeta.new_task(ts, test_tids[2], test_annotated_tids[2])
            @until length(ts.taskids) == 1
            @test length(ts.taskids) == 1
            @test length(ts.taskprops) == 2
            @test ts.taskids[1] == test_tids[2]
            key, tid_annotated, reservation = ts.taskprops[test_tids[2]]
            @test tid_annotated == test_annotated_tids[2]
            @test reservation == 0
            @test sharemode.ncreated == 2
            @test sharemode.ndeleted == 1
            @test sharemode.nshared == 1

            DagScheduler.FdbMeta.unreserve_task(ts, test_tids[1])
            @until length(ts.taskids) == 2
            @test length(ts.taskids) == 2
            @test length(ts.taskprops) == 2
            @test ts.taskids[1] == test_tids[2]
            @test ts.taskids[2] == test_tids[1]
            key, tid_annotated, reservation = ts.taskprops[test_tids[1]]
            @test tid_annotated == test_annotated_tids[1]
            @test reservation == 0
            @test sharemode.ncreated == 3
            @test sharemode.ndeleted == 1
            @test sharemode.nshared == 2

            DagScheduler.FdbMeta.reserve_task(ts, test_tids[1], 4)
            DagScheduler.FdbMeta.reserve_task(ts, test_tids[2], 4)
            @until length(ts.taskids) == 0
            @test length(ts.taskids) == 0
            @test length(ts.taskprops) == 2
            key, tid_annotated, reservation = ts.taskprops[test_tids[1]]
            @test tid_annotated == test_annotated_tids[1]
            @test reservation == 4
            key, tid_annotated, reservation = ts.taskprops[test_tids[2]]
            @test tid_annotated == test_annotated_tids[2]
            @test reservation == 4
            @test sharemode.ncreated == 3
            @test sharemode.ndeleted == 3
            @test sharemode.nshared == 0

            DagScheduler.FdbMeta.finish_task(ts, test_tids[1])
            DagScheduler.FdbMeta.finish_task(ts, test_tids[2])
            DagScheduler.FdbMeta.new_task(ts, test_tids[3:end], test_annotated_tids[3:end])
            @until length(ts.taskids) == 3
            @test length(ts.taskids) == 3
            @test length(ts.taskprops) == 5
            @test sharemode.ncreated == 6
            @test sharemode.ndeleted == 3
            @test sharemode.nshared == 3
            DagScheduler.FdbMeta.reserve_task(ts, test_tids[3], 4)
            DagScheduler.FdbMeta.reserve_task(ts, test_tids[4], 4)
            DagScheduler.FdbMeta.reserve_task(ts, test_tids[5], 4)
            @until length(ts.taskids) == 0
            @test length(ts.taskids) == 0
            @test length(ts.taskprops) == 5
            @test sharemode.ncreated == 6
            @test sharemode.ndeleted == 6
            @test sharemode.nshared == 0

            DagScheduler.FdbMeta.new_task(ts, test_tids[3:end], test_annotated_tids[3:end]; allow_duplicate=true)
            @until length(ts.taskids) == 3
            @test length(ts.taskprops) == 5
            @test sharemode.ncreated == 9

            DagScheduler.FdbMeta.reserve_task(ts, test_tids[3], 4)
            DagScheduler.FdbMeta.reserve_task(ts, test_tids[4], 4)
            DagScheduler.FdbMeta.reserve_task(ts, test_tids[5], 4)
            @until length(ts.taskids) == 0

            DagScheduler.FdbMeta.finish_task(ts, test_tids[3])
            DagScheduler.FdbMeta.finish_task(ts, test_tids[4])
            DagScheduler.FdbMeta.finish_task(ts, test_tids[5])
            DagScheduler.FdbMeta.prune_finished_tasks(ts)
            @test isempty(ts.taskids)
            @test isempty(ts.taskprops)
            @test sharemode.ncreated == 9
            @test sharemode.ndeleted == 9
            @test sharemode.nshared == 0

            DagScheduler.FdbMeta.stop_processing_events(ts)
            @test ts.eventprocessor == nothing

            @test ncallbacks > 0
        end
    end
    println("avg wait trigger time: $(mean(waits))")
end # testset fdb queue


text:
using OrdinaryDiffEq, Test, DiffEqBase
using ForwardDiff, Printf
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear,
                            prob_ode_bigfloatlinear, prob_ode_bigfloat2Dlinear
# use `PRINT_TESTS = true` to print the tests, including results
const PRINT_TESTS = false
print_results(x) = if PRINT_TESTS; @printf("%s \n", x) end

# points and storage arrays used in the interpolation tests
const interpolation_points = 0:1//2^(4):1
const interpolation_results_1d = fill(zero(typeof(prob_ode_linear.u0)), length(interpolation_points))
const interpolation_results_2d = Vector{typeof(prob_ode_2Dlinear.u0)}(undef, length(interpolation_points))
for idx in eachindex(interpolation_results_2d)
  interpolation_results_2d[idx] = zero(prob_ode_2Dlinear.u0)
end

f_linear_inplace = (du,u,p,t) -> begin @. du = 1.01 * u end
(::typeof(f_linear_inplace))(::Type{Val{:analytic}}, u0, p, t) = exp(1.01*t)*u0
prob_ode_linear_inplace = ODEProblem(f_linear_inplace, [0.5], (0.,1.))
const interpolation_results_1d_inplace = Vector{typeof(prob_ode_linear_inplace.u0)}(undef, length(interpolation_points))
for idx in eachindex(interpolation_results_1d_inplace)
  interpolation_results_1d_inplace[idx] = zero(prob_ode_linear_inplace.u0)
end

const deriv_test_points = range(0, stop=1, length=5)

# perform the regression tests
# NOTE: If you want to add new tests (for new algorithms), you have to run the
#       commands below to get numerical values for `tol_ode_linear` and
#       `tol_ode_2Dlinear`.
function regression_test(alg, tol_ode_linear, tol_ode_2Dlinear; test_diff1 = false, nth_der = 1, dertol=1e-6)
  PRINT_TESTS && println("\n", alg)

  sol = solve(prob_ode_linear, alg, dt=1//2^(2), dense=true)
  sol(interpolation_results_1d, interpolation_points)
  sol(interpolation_points[1])
  sol2 = solve(prob_ode_linear, alg, dt=1//2^(4), dense=true, adaptive=false)
  for i in eachindex(sol2)
    print_results( @test maximum(abs.(sol2[i] - interpolation_results_1d[i])) < tol_ode_linear )
  end
  for N in 1:nth_der
    # prevent CI error
    nth_der != 1 && @show N
    if test_diff1
      sol(interpolation_results_1d, interpolation_points, Val{N})
      der = sol(interpolation_points[1], Val{N})
      @test interpolation_results_1d[1] ≈ der
      for t in deriv_test_points
        deriv = sol(t, Val{N})
        @test deriv ≈ ForwardDiff.derivative(t -> sol(t, Val{N-1}), t) rtol=dertol
      end
    end
  end

  sol = solve(prob_ode_linear_inplace, alg, dt=1//2^(2), dense=true)
  sol(interpolation_results_1d_inplace, interpolation_points, idxs=1:1)
  sol(interpolation_results_1d_inplace, interpolation_points)
  sol(interpolation_points[1], idxs=1:1)
  sol(interpolation_points[1])
  for N in 1:nth_der
    # prevent CI error
    nth_der != 1 && @show N
    if test_diff1
      sol(interpolation_results_1d_inplace, interpolation_points, Val{N}, idxs=1:1)
      sol(interpolation_results_1d_inplace, interpolation_points, Val{N})
      sol(interpolation_points[1], Val{N}, idxs=1:1)
      der = sol(interpolation_points[1], Val{N})
      @test interpolation_results_1d_inplace[1] ≈ der
      for t in deriv_test_points
        deriv = sol(t, Val{N}, idxs=1)
        @test deriv ≈ ForwardDiff.derivative(t -> sol(t, Val{N-1}; idxs=1), t) rtol=dertol
      end
    end
  end

  sol  = solve(prob_ode_2Dlinear, alg, dt=1//2^(2), dense=true)
  sol(interpolation_results_2d,  interpolation_points)
  sol(interpolation_points[1])
  sol2 = solve(prob_ode_2Dlinear, alg, dt=1//2^(4), dense=true, adaptive=false)
  for i in eachindex(sol2)
    print_results( @test maximum(maximum.(abs.(sol2[i] - interpolation_results_2d[i]))) < tol_ode_2Dlinear)
  end
end

# Some extra tests using Euler()
prob = prob_ode_linear
sol  = solve(prob, Euler(), dt=1//2^(2), dense=true)
interpd_1d = sol(0:1//2^(4):1)
sol2 = solve(prob, Euler(), dt=1//2^(4), dense=true)
sol3 = solve(prob, Euler(), dt=1//2^(5), dense=true)

prob = prob_ode_2Dlinear
sol  = solve(prob, Euler(), dt=1//2^(2), dense=true)
interpd = sol(0:1//2^(4):1)

interpd_idxs = sol(0:1//2^(4):1,idxs=1:2:5)

@test minimum([interpd_idxs[i] == interpd[i][1:2:5] for i in 1:length(interpd)])

interpd_single = sol(0:1//2^(4):1,idxs=1)

@test typeof(interpd_single.u) <: Vector{Float64}

@test typeof(sol(0.5,idxs=1)) <: Float64

A = rand(4,2)
sol(A,0.777)
A == sol(0.777)
sol2 = solve(prob, Euler(), dt=1//2^(4), dense=true)

@test maximum(map((x)->maximum(abs.(x)),sol2 - interpd)) < .2

sol(interpd, 0:1//2^(4):1)

@test maximum(map((x)->maximum(abs.(x)),sol2 - interpd)) < .2

sol  = solve(prob, Euler(), dt=1//2^(2), dense=false)

@test !(sol(0.6)[4,2] ≈ 0)

# Euler
regression_test(Euler(), 0.2, 0.2)

# Midpoint
regression_test(Midpoint(), 1.5e-2, 2.3e-2)

println("SSPRKs")

# SSPRK22
regression_test(SSPRK22(), 1.5e-2, 2.5e-2; test_diff1 = true, nth_der = 2, dertol = 1e-15)

# SSPRK33
regression_test(SSPRK33(), 7.5e-4, 1.5e-3; test_diff1 = true, nth_der = 2, dertol = 1e-15)

# SSPRK53
regression_test(SSPRK53(), 2.5e-4, 4.0e-4; test_diff1 = true)

# SSPRK63
regression_test(SSPRK63(), 1.5e-4, 3.0e-4; test_diff1 = true)

# SSPRK73
regression_test(SSPRK73(), 9.0e-5, 2.0e-4; test_diff1 = true)

# SSPRK83
regression_test(SSPRK83(), 6.5e-5, 1.5e-4; test_diff1 = true)

# SSPRK432
regression_test(SSPRK432(), 4.0e-4, 8.0e-4; test_diff1 = true, nth_der = 2, dertol = 1e-13)

# SSPRK932
regression_test(SSPRK932(), 6.0e-5, 1.0e-4; test_diff1 = true)

# SSPRK54
regression_test(SSPRK54(), 3.5e-5, 5.5e-5)

# SSPRK104
regression_test(SSPRK104(), 1.5e-5, 3e-5)

println("RKs")

# RK4
regression_test(RK4(), 4.5e-5, 1e-4)

# CarpenterKennedy2N54
regression_test(CarpenterKennedy2N54(), 3.0e-5, 5.0e-5)

# DP5
regression_test(DP5(), 5e-6, 1e-5; test_diff1 = true, nth_der = 4, dertol = 1e-14)

# BS3
regression_test(BS3(), 5e-4, 8e-4)

# OwrenZen3
regression_test(OwrenZen3(), 1.5e-4, 2.5e-4; test_diff1 = true, nth_der = 3, dertol = 1e-9)

# OwrenZen4
regression_test(OwrenZen4(), 6.5e-6, 1.5e-5; test_diff1 = true, nth_der = 4, dertol = 1e-10)

# OwrenZen5
regression_test(OwrenZen5(), 1.5e-6, 2.5e-6; test_diff1 = true, nth_der = 5, dertol = 1e-8)

# Tsit5
regression_test(Tsit5(), 2e-6, 4e-6; test_diff1 = true, nth_der = 4, dertol = 1e-6)

# TanYam7
regression_test(TanYam7(), 4e-4, 6e-4)

# TsitPap8
regression_test(TsitPap8(), 1e-3, 3e-3)

# Feagin10
regression_test(Feagin10(), 6e-4, 9e-4)

# BS5
regression_test(BS5(), 4e-8, 6e-8; test_diff1 = true, nth_der = 1, dertol = 1e-12)
regression_test(BS5(lazy=false), 4e-8, 6e-8; test_diff1 = true, nth_der = 1, dertol = 1e-12)

prob = prob_ode_linear
sol  = solve(prob, BS5(), dt=1//2^(1), dense=true, adaptive=false)
interpd_1d_long = sol(0:1//2^(7):1)
sol2 = solve(prob, BS5(), dt=1//2^(7), dense=true, adaptive=false)
print_results( @test maximum(map((x)->maximum(abs.(x)),sol2 - interpd_1d_long)) < 2e-7 )

# DP8
regression_test(DP8(), 2e-7, 3e-7; test_diff1 = true, nth_der = 1, dertol = 1e-15)

prob = prob_ode_linear
sol  = solve(prob, DP8(), dt=1//2^(2), dense=true)
sol(interpd_1d_long, 0:1//2^(7):1) # inplace update
sol2 = solve(prob, DP8(), dt=1//2^(7), dense=true, adaptive=false)
print_results( @test maximum(map((x)->maximum(abs.(x)),sol2 - interpd_1d_long)) < 2e-7 )

println("Verns")

# Vern6
regression_test(Vern6(), 7e-8, 7e-8; test_diff1 = true, nth_der = 1, dertol = 1e-9)
regression_test(Vern6(lazy=false), 7e-8, 7e-8; test_diff1 = true, nth_der = 1, dertol = 1e-9)

prob = remake(prob_ode_bigfloatlinear;u0=big(0.5))

sol  = solve(prob, Vern6(), dt=1//2^(2), dense=true)
interpd_1d_big = sol(0:1//2^(7):1)
sol2 = solve(prob, Vern6(), dt=1//2^(7), dense=true, adaptive=false)
print_results( @test maximum(map((x)->maximum(abs.(x)),sol2[:] - interpd_1d_big)) < 5e-8 )

prob_ode_bigfloatveclinear = ODEProblem((u,p,t)->p*u,[big(0.5)],(0.0,1.0),big(1.01))
prob = prob_ode_bigfloatveclinear
sol  = solve(prob, Vern6(), dt=1//2^(2), dense=true)
interpd_big = sol(0:1//2^(4):1)
sol2 = solve(prob, Vern6(), dt=1//2^(4), dense=true, adaptive=false)
print_results( @test maximum(map((x)->maximum(abs.(x)),sol2 - interpd_big)) < 5e-8 )

# Vern7
regression_test(Vern7(), 3e-9, 5e-9; test_diff1 = true, nth_der = 1, dertol = 1e-10)
regression_test(Vern7(lazy=false), 3e-9, 5e-9; test_diff1 = true, nth_der = 1, dertol = 1e-10)

# Vern8
regression_test(Vern8(), 3e-8, 5e-8; test_diff1 = true, nth_der = 1, dertol = 1e-7)
regression_test(Vern8(lazy=false), 3e-8, 5e-8; test_diff1 = true, nth_der = 1, dertol = 1e-7)

# Vern9
regression_test(Vern9(), 1e-9, 2e-9; test_diff1 = true, nth_der = 4, dertol = 2e-2)
regression_test(Vern9(lazy=false), 1e-9, 2e-9; test_diff1 = true, nth_der = 4, dertol = 2e-2)

println("Rosenbrocks")

# Rosenbrock23
regression_test(Rosenbrock23(), 3e-3, 6e-3; test_diff1 = true, nth_der = 1, dertol = 1e-14)

# Rosenbrock32
regression_test(Rosenbrock32(), 4e-4, 6e-4; test_diff1 = true, nth_der = 1, dertol = 1e-14)

# Rodas4
regression_test(Rodas4(), 8.5e-6, 2e-5)

# GenericImplicitEuler
regression_test(GenericImplicitEuler(), 0.2, 0.354)

# GenericTrapezoid
regression_test(GenericTrapezoid(), 7e-3, 1.4e-2)

# ExplicitRK
regression_test(ExplicitRK(), 7e-5, 2e-4)

prob = prob_ode_linear
sol  = solve(prob, ExplicitRK(), dt=1//2^(2), dense=true)
# inplace interp of solution
sol(interpd_1d_long,0:1//2^(7):1)
sol2 = solve(prob, ExplicitRK(), dt=1//2^(7), dense=true, adaptive=false)
print_results( @test maximum(map((x)->maximum(abs.(x)),sol2 - interpd_1d_long)) < 6e-5 )

prob = prob_ode_2Dlinear
sol  = solve(prob, ExplicitRK(), dt=1//2^(2), dense=true)
sol(interpd, 0:1//2^(4):1)
sol2 = solve(prob, ExplicitRK(), dt=1//2^(4), dense=true, adaptive=false)
print_results( @test maximum(map((x)->maximum(abs.(x)),sol2 - interpd)) < 2e-4 )


text:
# workspace()
cd("/home/astupidbear/.julia/v0.4/JuLearn/examples/rnn/")
using JuLearn,Plots,ArrayFire,JLD
@Eval lr=1e-2 batchsize=100

include("predict.jl")
include("dataProvider.jl")

inputBase,outputBase=dataProvider();
T=10
net=rnn([:d,:h,:h,:l],[1 1;1 10;10 1;1 1],T,:regression);

# backprop
err=backprop(net,Array(inputBase),Array(outputBase),nIter=30)
p=plot(err);savefig(p,"error")

# prediction
input=inputBase[:,100+(0:T÷2-1),1]
target=outputBase[:,100+(0:T-1),1]'
prediction=predict(net,Array(input))
prediction=map(i->Array(prediction[i])[1],1:T)
p=plot([target prediction],label=["target" "prediction"]);
savefig(p,"prediction.png");run(`xdg-open prediction.png`)


text:
# https://adventofcode.com/2020/day/25

"""
Single transformation step
"""
function transform(value::Integer, subject_number::Integer)
    return rem(value * subject_number, 20201227)
end

"""
Using the `subject_number` and `loop_size`, calculate the `key`
"""
function get_key(subject_number::Integer, loop_size::Integer)
    value = 1
    for _ in 1:loop_size
        value = transform(value, subject_number)
    end
    return value
end

# example: subject number -> public key
@assert get_key(7, 8) == 5764801
@assert get_key(7, 11) == 17807724
# example: public key -> encryption key
@assert get_key(17807724, 8) == 14897079
@assert get_key(5764801, 11) == 14897079

"""
Using the `subject_number` and `public_key`, guess the `loop_size`
"""
function get_loop_size(subject_number::Integer, public_key::Integer)
    loop_size = 0
    value = 1
    while value != public_key
        loop_size += 1
        value = transform(value, subject_number)
    end
    return loop_size
end

# example: (subject number, public key) -> loop size
@assert get_loop_size(7, 5764801) == 8
@assert get_loop_size(7, 17807724) == 11
# example: (public key, encryption key) -> loop size
@assert get_loop_size(17807724, 14897079) == 8
@assert get_loop_size(5764801, 14897079) == 11

function part1(card_key, door_key)
    card_loop_size = get_loop_size(7, card_key)
    door_loop_size = get_loop_size(7, door_key)

    card_encryption_key = get_key(card_key, door_loop_size)
    door_encryption_key = get_key(door_key, card_loop_size)

    @assert card_encryption_key == door_encryption_key

    return card_encryption_key
end

@assert part1(5764801, 17807724) == 14897079

test = (14205034, 18047856)
println("Part 1: $(result1 = part1(test...))")
# println("Part 2: $(result2 = part2(test...))")

@assert result1 == 297257
# @assert result2 == 418819514477

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
"""
Provides a collection of projection operators and error functionals.
"""

###### Projections ######

"Fourier modulus projection"
function pM(g, sqrtI)
    ang = angle.(fft(g))
    res = real(ifft(sqrtI .* exp.(im .* ang)))
end

"""
Incomplete (Elser) modulus projection

A variant of the modulus projection that leaves 0'th frequency
invariant or sets it to 0 if it is negative. This projection is
used, e.g., in
> "Benchmark problems for phase retrieval", V. Elser, T.-Y. Lan & T. Bendory
"""
function pM0(g, sqrtI)
    goalAbs = sqrtI
    goalAbs[1,1] = maximum([sum(g), 0.0])
    return real(ifft(goalAbs .* exp.(im .* angle.(fft(g)))))
end

"Positivity projection"
function pP(g)
    res = real(g .* (g .>= 0))
end

"Soft sparsity projection"
function pTa(g, alpha=0)
    res = real(g .* (g .>= alpha))
end

"""
Hard sparsity projection

Set all but suppSize largest elements of g to zero.
"""
function pTn(g, suppSize)
    whereLargest = partialsortperm(g[:], 1:suppSize, rev=true)
    res = zeros(size(g))
    res[whereLargest] = g[:][whereLargest]
    return res
end

###### Error functionals ######
"Modulus energy functional"
function eM(g, sqrtI)
    res = 0.5 * sum((abs.(fft(g)) - sqrtI).^2) / length(g)
end

"Incomplete modulus energy functional"
function eM0(g, sqrtI)
        res = 0.5 * sum((abs.(fft(g)[2:end]) - sqrtI[2:end]).^2) / length(g)
end

"Positive energy functional"
function eP(g)
    res = 0.5 * sum(real(g .* (g .< 0)).^2)
end

"Soft sparsity energy functional"
function eTa(g, alpha=0)
    res = 0.5 * sum((g - pTa(g, alpha)).^2)
end

"Hard sparsity energy functional"
function eTn(g, suppSize)
    res = 0.5 * sum((g - pTn(g, suppSize)).^2)
end

max_stars_repo_name: tsipenyuk/FPRS

max_issues_repo_name: tsipenyuk/FPRS

max_forks_repo_name: tsipenyuk/FPRS


text:
module Dwcl
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "1.0.1"

include("dashwindowcloselistener.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "dwcl",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "dwcl.min.js",
    external_url = "https://unpkg.com/dwcl@1.0.1/dwcl/dwcl.min.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dwcl.min.js.map",
    external_url = "https://unpkg.com/dwcl@1.0.1/dwcl/dwcl.min.js.map",
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end

max_stars_repo_name: raymon-ai/dwcl

max_issues_repo_name: raymon-ai/dwcl

max_forks_repo_name: raymon-ai/dwcl


text:
"""
A pure julia wrapper of R_tryEval.
"""
function tryEval{S<:Sxp}(expr::Ptr{S}, env::Ptr{EnvSxp})
    disable_sigint() do
        status = Array(Cint,1)
        protect(expr)
        protect(env)
        val = ccall((:R_tryEval,libR),UnknownSxpPtr,(Ptr{S},Ptr{EnvSxp},Ptr{Cint}),expr,env,status)
        unprotect(2)
        val, status[1]
    end
end

"""
Evaluate an R symbol or language object (i.e. a function call) in an R
try/catch block, returning a Sxp pointer.
"""
function reval_p{S<:Sxp}(expr::Ptr{S}, env::Ptr{EnvSxp})
    val, status = tryEval(expr, env)
    flush_print_buffer(STDOUT)
    if status !=0
        error("RCall.jl: ", takebuf_string(errorBuffer))
    elseif nb_available(errorBuffer) != 0
        warn("RCall.jl: ", takebuf_string(errorBuffer))
    end
    sexp(val)
end

"""
Evaluate an R expression array iteratively.
"""
function reval_p(expr::Ptr{ExprSxp}, env::Ptr{EnvSxp})
    local val           # the value of the last expression is returned
    protect(expr)
    try
        for e in expr
            val = reval_p(e,env)
        end
    finally
        unprotect(1)
    end
    val
end

reval_p{S<:Sxp}(expr::Ptr{S}, env::RObject{EnvSxp}) = reval_p(expr,sexp(env))
reval_p{S<:Sxp}(s::Ptr{S}) = reval_p(s,Const.GlobalEnv)

"""
Evaluate an R symbol or language object (i.e. a function call) in an R
try/catch block, returning an RObject.
"""
reval(r::RObject, env=Const.GlobalEnv) = reval(sexp(r), env)
function reval{S<:Sxp}(s::Ptr{S}, env=Const.GlobalEnv)
    ss = protect(s)
    local val
    try
        val = RObject(reval_p(ss,sexp(env)))
    finally
        unprotect(1)
    end
    val
end
reval(str::Union{AbstractString,Symbol}, env=Const.GlobalEnv) =
    reval(rparse_p(str))

"A pure julia wrapper of R_ParseVector"
function parseVector{S<:Sxp}(st::Ptr{StrSxp}, sf::Ptr{S}=sexp(Const.NilValue))
    protect(st)
    status = Array(Cint,1)
    val = ccall((:R_ParseVector,libR),UnknownSxpPtr,
                (Ptr{StrSxp},Cint,Ptr{Cint},UnknownSxpPtr),
                st,-1,status,sf)
    unprotect(1)
    val, status[1]
end

"Get the R parser error msg for the previous parsing result."
function getParseErrorMsg()
    Compat.unsafe_string(cglobal((:R_ParseErrorMsg, libR), UInt8))
end

"Parse a string as an R expression, returning a Sxp pointer."
function rparse_p(st::Ptr{StrSxp})
    val, status = parseVector(st)
    if status == 2 || status == 3
        error("RCall.jl: ", getParseErrorMsg())
    elseif status == 4
        throw(EOFError())
    end
    sexp(val)
end
rparse_p(st::AbstractString) = rparse_p(sexp(st))

"Parse a string as an R expression, returning an RObject."
rparse(st::AbstractString) = RObject(rparse_p(st))

"""
Evaluate and convert the result of a string as an R expression.
"""
rcopy(str::AbstractString) = rcopy(reval_p(rparse_p(str)))
rcopy(sym::Symbol) = rcopy(reval_p(sexp(sym)))
rcopy{T}(::Type{T}, str::AbstractString) = rcopy(T, reval_p(rparse_p(str)))
rcopy{T}(::Type{T}, sym::Symbol) = rcopy(T, reval_p(sexp(sym)))


text:
"""
    BoxLeastSquares.model(::BLSPeriodogram; kwargs...)

Create a transit model using the data and best-fitting parameters from the given BLS periodogram. Any keyword parameters can be overriden.
"""
function model(bls::BLSPeriodogram; kwargs...)
    # compute depth
    pars = params(bls)
    return model(bls.t, bls.y, bls.yerr; pars..., kwargs...)
end

"""
    BoxLeastSquares.model(t, y, [yerr]; period, duration, t0)

Evaluate the transit model on the given time grid. If `yerr` is not provided, it will default to 1. The following transit parameters must be set:
* `period` orbital period in the same units as `t`
* `duration` the transit duration in the same units as `t`
* `t0` the transit time (middle of transit) in the same units as `t`

If you are using `Unitful.jl`, the unit conversions will be made automatically.
"""
function model(t, y, yerr=fill!(similar(y), one(eltype(y))); period, duration, t0, kwargs...)
    # compute depth
    hp = 0.5 * period
    y_in = y_out = zero(typeof(first(y)))
    ivar_in = ivar_out = zero(typeof(inv(first(yerr)^2)))
    @inbounds for idx in eachindex(t, y, yerr)
        transiting = abs((t[idx] - t0 + hp) % period - hp) < 0.5 * duration
        iv = inv(yerr[idx]^2)
        χ = y[idx] * iv
        if transiting
            y_in += χ
            ivar_in += iv
        else
            y_out += χ
            ivar_out += iv
        end
    end
    y_in /= ivar_in
    y_out /= ivar_out
    y_model = similar(y)
    @inbounds for idx in eachindex(t)
        transiting = abs((t[idx] - t0 + hp) % period - hp) < 0.5 * duration
        y_model[idx] = ifelse(transiting, y_in, y_out)
    end
    return y_model
end


text:
@testset "Core" begin
        path = "test/data/case6ww.m"
        network = PowerFlowNetwork(path, "MATPOWERM")
        @test nbus(network) == 6
        @test nbranch(network) == 11
        @test ngen(network) == 3

        @test has_continuous_index(network) == true
        network.bus = network.bus[1:end .!= 2, :]
        @test has_continuous_index(network) == false

        @test is_disjoint(network) == false
        network.branch = DataFrame(zeros(0, 21), ["SRC", "DST", "BR_R", "BR_X", "BR_B", "RATE_A", "RATE_B", "RATE_C",
                                                  "TAP", "SHIFT", "BR_STATUS", "ANGMIN", "ANGMAX", "PF", "QF", "PT",
                                                  "QT", "MU_SF", "MU_ST", "MU_ANGMIN", "MU_ANGMAX"])
        @test is_disjoint(network) == true

        path = "test/data/C2FEN02312/scenario_3/"
        network = PowerFlowNetwork(path, "RAWGO")
        raw_path = joinpath(path, "case.raw")
        @test nbranch_unique(network) == nbranch_unique(raw_path)
        merge_duplicate_branch!(network)
        @test nbranch_unique(network) == nbranch_unique(raw_path)
        @test nbranch(network) == nbranch_unique(network)
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using VisualParticipationAnalyticsServers
using Test

@testset "VisualParticipationAnalyticsServers.jl" begin
end


text:
module CGS

export C, GRAV, BOLTZMANN, STEFANBOLTZMANN, MHYDROGEN, YEAR_SEC, AU, PC, MSUN, PLANCK
export RSUN, MEARTH, REARTH, MJUPITER, RJUPITER, AJUPITER, KILOGRAM, METER, MMOON, RMOON, TSUN

PLANCK = 6.62606983e-27
GRAV = 6.67384e-8
BOLTZMANN = 1.3806488e-16
STEFANBOLTZMANN = 5.670373e-5
MHYDROGEN = 1.00794*1.66053886e-24
YEAR_SEC = 3.15569e7
C = 29979245800e0

#astrophysical constants
AU = 1.49598e13
PC = 3.08568025e18

MSUN = 1.98892e33
RSUN = 6.955e10
TSUN = 5777.0

MEARTH = 5.9742e27
REARTH = 6.3781e8
MMOON = 7.342e25
RMOON = 1737.1e5

MJUPITER = 1.8986e30
RJUPITER = 6.9911e9
AJUPITER = 4.84143144246472090 # AU

KILOGRAM = 1e3
METER = 1e2

end


text:
using PMDPs
using Test
using POMDPSimulators

using StaticArrays
using POMDPs
using MCTS, DiscreteValueIteration
using StatsBase
using Plots
using POMDPSimulators
using D3Trees
using POMDPPolicies
using POMDPLinter
using Random
using DataFrames
using POMDPSimulators

mdp_vi = PMDPs.create_PMDP(PMDPe)
mdp_mc = PMDPs.create_PMDP(PMDPg)

policy = PMDPs.get_VI_policy(mdp_vi)
planner = PMDPs.get_MCTS_planner(mdp_mc)

rng = MersenneTwister(1234)

# s0 = rand(rng, initialstate(mdp_mc))

# function run_sim(mdp::PMDP, policy::Policy; rng_seed=1234)
#     rng = MersenneTwister(rng_seed)
#     hr = HistoryRecorder(max_steps=100, capture_exception=true, rng=rng)
#     h = simulate(hr, mdp, policy)
#     collect(eachstep(h, "s, a, r, user_budget"))
#     # sum(h[:r])
# end

# rng = MersenneTwister(rand_seed)
# hr = HistoryRecorder(max_steps=100, capture_exception=false, rng=rng)
# h_mc = simulate(hr, mdp_mc, planner)
# collect(eachstep(h, "s, a, r, info"))

rng_seed = 1
max_steps = mdp_mc.T+1

h_mc = run_sim(mdp_mc, planner; max_steps = max_steps, rng_seed = rng_seed)
h_vi = run_sim(mdp_mc, policy; max_steps = max_steps, rng_seed = rng_seed)

hindsight = PMDPs.LP.MILP_hindsight_pricing(mdp_mc, h_mc; optimization_goal="revenue", verbose=false)
flatrate = PMDPs.flatrate_pricing(mdp_mc, h_mc)

hindsight, flatrate, get_stats(h_mc), get_stats(h_vi)


text:
# these types and functions integrate the Boltzmann hierarchy through time

abstract type PerturbationIntegrator end
struct BasicNewtonian <: PerturbationIntegrator end

# a container for everything needed to integrate a hierarchy at wavenumber k
struct Hierarchy{T<:Real, PI<:PerturbationIntegrator, CP<:AbstractCosmoParams{T},
                 BG<:AbstractBackground, IH<:AbstractIonizationHistory, Tk<:Real}
    integrator::PI
    par::CP
    bg::BG
    ih::IH
    k::Tk
    ℓᵧ::Int  # Boltzmann hierarchy cutoff, i.e. Seljak & Zaldarriaga
end
Hierarchy(integrator::PerturbationIntegrator, par::AbstractCosmoParams, bg::AbstractBackground,
    ih::AbstractIonizationHistory, k::Real, ℓᵧ=8) = Hierarchy(integrator, par, bg, ih, k, ℓᵧ)

function boltsolve(hierarchy::Hierarchy{T}, ode_alg=Rodas5(); reltol=1e-10) where T
    xᵢ = first(hierarchy.bg.x_grid)
    u₀ = initial_conditions(xᵢ, hierarchy)
    prob = ODEProblem{true}(hierarchy!, u₀, (xᵢ , zero(T)), hierarchy)
    sol = solve(prob, ode_alg, reltol=reltol, saveat=hierarchy.bg.x_grid, dense=false)
    return sol
end

# basic Newtonian gauge: establish the order of perturbative variables in the ODE solve
function unpack(u, hierarchy::Hierarchy{T, BasicNewtonian}) where T
    ℓᵧ = hierarchy.ℓᵧ
    ℓ_ν = 10 #Callin06, for now
    Θ = OffsetVector(view(u, 1:(ℓᵧ+1)), 0:ℓᵧ)  # indexed 0 through ℓᵧ
    Θᵖ = OffsetVector(view(u, (ℓᵧ+2):(2ℓᵧ+2)), 0:ℓᵧ)  # indexed 0 through ℓᵧ
    𝒩 = OffsetVector(view(u, (2(ℓᵧ+1) + 1):(2(ℓᵧ+1)+(ℓ_ν+1))) , 0:ℓ_ν)  # indexed 0 through ℓ_ν
    Φ, δ, v, δ_b, v_b = view(u, (2(ℓᵧ+1)+(ℓ_ν+1)+1):(2(ℓᵧ+1)+(ℓ_ν+1)+5)) #getting a little messy...
    return Θ, Θᵖ, 𝒩, Φ, δ, v, δ_b, v_b
end

# BasicNewtonian comes from Callin+06 and the Dodelson textbook (dispatches on hierarchy.integrator)
function hierarchy!(du, u, hierarchy::Hierarchy{T, BasicNewtonian}, x) where T
    # compute cosmological quantities at time x, and do some unpacking
    k, ℓᵧ, par, bg, ih = hierarchy.k, hierarchy.ℓᵧ, hierarchy.par, hierarchy.bg, hierarchy.ih
    Ω_r, Ω_b, Ω_m, N_ν, H₀² = par.Ω_r, par.Ω_b, par.Ω_m, par.N_ν, bg.H₀^2 #add N_ν≡N_eff
    ℋₓ, ℋₓ′, ηₓ, τₓ′, τₓ′′ = bg.ℋ(x), bg.ℋ′(x), bg.η(x), ih.τ′(x), ih.τ′′(x)
    a = x2a(x)
    R = 4Ω_r / (3Ω_b * a)
    Ω_ν =  7N_ν/8 *(4/11)^(4/3) *Ω_r
    ℓ_ν = 10 #again, for now

    Θ, Θᵖ, 𝒩, Φ, δ, v, δ_b, v_b = unpack(u, hierarchy)  # the Θ, Θᵖ, 𝒩 are views (see unpack)
    Θ′, Θᵖ′, 𝒩′,_, _, _, _, _ = unpack(du, hierarchy)  # will be sweetened by .. syntax in 1.6

    # metric perturbations
    #Ψ = -Φ - 12H₀² / k^2 / a^2 * (Ω_r * Θ[2])
    Ψ = -Φ - 12H₀² / k^2 / a^2 * (Ω_r * Θ[2] + Ω_ν * 𝒩[2]) #add rel quadrupole
    Φ′ = Ψ - k^2 / (3ℋₓ^2) * Φ + H₀² / (2ℋₓ^2) * (
        Ω_m * a^(-1) * δ + Ω_b * a^(-1) * δ_b + 4Ω_r * a^(-2) * Θ[0]
        + 4Ω_ν * a^(-2) * 𝒩[0]) #add rel monopole on this line

    # matter
    δ′ = k / ℋₓ * v - 3Φ′
    v′ = -v - k / ℋₓ * Ψ
    δ_b′ = k / ℋₓ * v_b - 3Φ′
    v_b′ = -v_b - k / ℋₓ * Ψ + τₓ′ * R * (3Θ[1] + v_b)

    # relativistic neutrinos (massless)
    𝒩′[0] = -k / ℋₓ * 𝒩[1] - Φ′
    𝒩′[1] = k/(3ℋₓ)*𝒩[0] - 2*k/(3ℋₓ)*𝒩[2] + k/(3ℋₓ)*Ψ
    for ℓ in 2:(ℓ_ν-1) #ℓ_ν same as ℓᵧ for massless nu for now
        𝒩′[ℓ] =  k / ((2ℓ+1) * ℋₓ) *( ℓ*𝒩[ℓ-1] - (ℓ+1)*𝒩[ℓ+1])
    end
    #truncation
    𝒩′[ℓ_ν] =  k / ℋₓ  * 𝒩[ℓ_ν-1] - (ℓ_ν+1)/(ℋₓ *ηₓ) *𝒩[ℓ_ν]#Callin 06

    # WIP: massive nu
    # # neutrinos (massive, MB 57) - change convention
    # #units not yet right
    # Ψ_ν′[0] = -q*k / ϵ * Ψ_ν[1] - Φ′ *dnlnf0dlnq #FIXME dnln, def Psi,IC, Einstein int
    # Ψ_ν′[1] = q*k/(3 ϵ)*(𝒩[0] - 2*Ψ_ν[2]) - This k*ϵ / (3*q)*Ψ *dnlnf0dlnq
    # I think can't mutate the u variable...
    # Ψ_ν[ℓ_νm] = (2*ℓ_νm+1)*ϵ/(q*k*ηₓ)*Ψ_ν[ℓ_νm] - Ψ_ν[ℓ_νm-1] #truncation of MB (51)
    # for ℓ in 2:(ℓ_νm-1) #ℓ_νm should be smaller than massless case
    #     Ψ_ν′[ℓ] =  q*k / ((2ℓ+1) * ϵ) *( ℓ*Ψ_ν[ℓ-1] - (ℓ+1)*Ψ_ν[ℓ+1])
    # end

    # photons
    Π = Θ[2] + Θᵖ[2] + Θᵖ[0]
    Θ′[0] = -k / ℋₓ * Θ[1] - Φ′
    Θ′[1] = k / (3ℋₓ) * Θ[0] - 2k / (3ℋₓ) * Θ[2] + k / (3ℋₓ) * Ψ + τₓ′ * (Θ[1] + v_b/3)
    for ℓ in 2:(ℓᵧ-1)
        Θ′[ℓ] = ℓ * k / ((2ℓ+1) * ℋₓ) * Θ[ℓ-1] -
            (ℓ+1) * k / ((2ℓ+1) * ℋₓ) * Θ[ℓ+1] + τₓ′ * (Θ[ℓ] - Π * δ_kron(ℓ, 2) / 10)
    end

    # polarized photons
    Θᵖ′[0] = -k / ℋₓ * Θᵖ[1] + τₓ′ * (Θᵖ[0] - Π / 2)
    for ℓ in 1:(ℓᵧ-1)
        Θᵖ′[ℓ] = ℓ * k / ((2ℓ+1) * ℋₓ) * Θᵖ[ℓ-1] -
            (ℓ+1) * k / ((2ℓ+1) * ℋₓ) * Θᵖ[ℓ+1] + τₓ′ * (Θᵖ[ℓ] - Π * δ_kron(ℓ, 2) / 10)
    end

    # photon boundary conditions: diffusion damping
    Θ′[ℓᵧ] = k / ℋₓ * Θ[ℓᵧ-1] - (ℓᵧ + 1) / (ℋₓ * ηₓ) + τₓ′ * Θ[ℓᵧ]
    Θᵖ′[ℓᵧ] = k / ℋₓ * Θᵖ[ℓᵧ-1] - (ℓᵧ + 1) / (ℋₓ * ηₓ) + τₓ′ * Θᵖ[ℓᵧ]

    du[2(ℓᵧ+1)+(ℓ_ν+1)+1:2(ℓᵧ+1)+(ℓ_ν+1)+5] .= Φ′, δ′, v′, δ_b′, v_b′  # put non-photon perturbations back in
    return nothing
end

# BasicNewtonian Integrator (dispatches on hierarchy.integrator)
function initial_conditions(xᵢ, hierarchy::Hierarchy{T, BasicNewtonian}) where T
    k, ℓᵧ, par, bg, ih = hierarchy.k, hierarchy.ℓᵧ, hierarchy.par, hierarchy.bg, hierarchy.ih
    ℓ_ν = 10 #again, for now
    u = zeros(T, 2(ℓᵧ+1)+(ℓ_ν+1)+5)
    ℋₓ, ℋₓ′, ηₓ, τₓ′, τₓ′′ = bg.ℋ(xᵢ), bg.ℋ′(xᵢ), bg.η(xᵢ), ih.τ′(xᵢ), ih.τ′′(xᵢ)
    Θ, Θᵖ, 𝒩, Φ, δ, v, δ_b, v_b = unpack(u, hierarchy)  # the Θ, Θᵖ are mutable views (see unpack)
    H₀²,aᵢ² = bg.H₀^2,exp(xᵢ)^2

    # metric and matter perturbations
    Φ = 1.0
    δ = 3Φ / 2
    δ_b = δ
    v = k / (2ℋₓ) * Φ
    v_b = v

    # photon hierarchy
    Θ[0] = Φ / 2
    Θ[1] = -k * Φ / (6ℋₓ)
    Θ[2] = -8k / (15ℋₓ * τₓ′) * Θ[1]
    Θᵖ[0] = (5/4) * Θ[2]
    Θᵖ[1] = -k / (4ℋₓ * τₓ′) * Θ[2]
    Θᵖ[2] = (1/4) * Θ[2]
    for ℓ in 3:ℓᵧ
        Θ[ℓ] = -ℓ/(2ℓ+1) * k/(ℋₓ * τₓ′) * Θ[ℓ-1]
        Θᵖ[ℓ] = -ℓ/(2ℓ+1) * k/(ℋₓ * τₓ′) * Θᵖ[ℓ-1]
    end

    # neutrino hierarchy
    # for now we assume xᵢ is before neutrinos decouple
    f_ν = 1/(1 + 1/(7par.N_ν/8 *(4/11)^(4/3)))
    𝒩[0] = Θ[0]
    𝒩[1] = Θ[1]
    𝒩[2] = - (k^2 *aᵢ²*Φ) / (12H₀²) * 1 / (1 + 5f_ν/2) #Callin06 (71)
    for ℓ in 3:ℓ_ν
        𝒩[ℓ] = k/((2ℓ+1)ℋₓ) * 𝒩[ℓ-1] #approximation of Callin06 (72)
    end

    #WIP: massive nu
    #FIXME: nonrelativistic transition for massive species, needs to go in bg
    #^this will have to wait for m_ν to be added to pars
    #below notation is not right yet, starting from MB
    #x_nr = m_ν/5.3e-4 -1 #m_ν in eV (PDG26-pg3)
    # same as photons for 0,1
    # # massive #FIXME ingegrate the q moments, get dlnf0dlnq,define ϵ
    # σ_ν= (k*ηₓ)^2 *Ψ/ 15
    # Ψ_ν[0] = -δ_ν *dlnf0dlnq
    # Ψ_ν[1] = -ϵ/(3*q*k) θ_ν *dlnf0dlnq  #change θ to -k/ℋₓ v
    # Ψ_ν[2] = -σ_ν/2  *dlnf0dlnq
    #ignore ℓ>2, small

    u[(2(ℓᵧ+1)+(ℓ_ν+1)+1):(2(ℓᵧ+1)+(ℓ_ν+1)+5)] .= Φ, δ, v, δ_b, v_b  # write u with our variables
    return u
end

# TODO: this could be extended to any Newtonian gauge integrator if we specify the
# Bardeen potential Ψ and its derivative ψ′ for an integrator, or we saved them
function source_function(du, u, hierarchy::Hierarchy{T, BasicNewtonian}, x) where T
    # compute some quantities
    k, ℓᵧ, par, bg, ih = hierarchy.k, hierarchy.ℓᵧ, hierarchy.par, hierarchy.bg, hierarchy.ih
    H₀² = bg.H₀^2
    ℋₓ, ℋₓ′, ℋₓ′′ = bg.ℋ(x), bg.ℋ′(x), bg.ℋ′′(x)
    τₓ, τₓ′, τₓ′′ = ih.τ(x), ih.τ′(x), ih.τ′′(x)
    g̃ₓ, g̃ₓ′, g̃ₓ′′ = ih.g̃(x), ih.g̃′(x), ih.g̃′′(x)
    a = x2a(x)

    Θ, Θᵖ, Φ, δ, v, δ_b, v_b = unpack(u, hierarchy)  # the Θ, Θᵖ are mutable views (see unpack)
    Θ′, Θᵖ′, Φ′, δ′, v′, δ_b′, v_b′ = unpack(du, hierarchy)

    # recalulate these since we didn't save them
    Ψ = -Φ - 12H₀² / k^2 / a^2 * par.Ω_r * Θ[2]
    Ψ′ = -Φ′ - 12H₀² / k^2 / a^2 * par.Ω_r * (Θ′[2] - 2 * Θ[2])
    Π = Θ[2] + Θᵖ[2] + Θᵖ[0]
    Π′ = Θ′[2] + Θᵖ′[2] + Θᵖ′[0]

    term1 =  g̃ₓ * (Θ[0] + Ψ + Π/4) + exp(-τₓ) * (Ψ′ - Φ′)
    term2 = (-1/k) * (ℋₓ′ * g̃ₓ * v_b + ℋₓ * g̃ₓ′ * v_b + ℋₓ * g̃ₓ * v_b′)
    Π′′ = 2k / (5ℋₓ) * (-ℋₓ′ / ℋₓ * Θ[1] + Θ′[1]) + (3/10) * (τₓ′′ * Π + τₓ′ * Π′) -
        3k / (5ℋₓ) * (-ℋₓ′ / ℋₓ * (Θ[3] + Θᵖ[1] + Θᵖ[3]) + (Θ′[3] + Θᵖ′[1] + Θᵖ′[3]))
    term3 = (3/(4k^2)) * (
        (ℋₓ′^2 + ℋₓ * ℋₓ′′) * g̃ₓ * Π + 3 * ℋₓ * ℋₓ′ * (g̃ₓ′ * Π + g̃ₓ * Π′) +
        ℋₓ^2 * (g̃ₓ′′ * Π + 2g̃ₓ′ * Π′ + g̃ₓ * Π′′))
    return term1 + term2 + term3
end


text:
@testset "test/sub.jl" begin
    @test true
end

module SubSub
using ReTest
__revise_mode__ = :eval

@testset "SubSub" begin
    @test true
end

end # SubSub

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
using ObjFileBase, COFF

timeline, modules = replay("/home/kfischer/.local/share/rr/wine64-12");
function step_to_syscall(number)
  succ = true
  while succ
      session = current_session(timeline)
      frame = icxx"$session->current_trace_frame();"
      event = icxx"$frame.event();"
      succ = icxx"$(RR.step!(session)).status != rr::REPLAY_EXITED;"
      # If we're exiting the connect call
      if icxx"$event.is_syscall_event();" &&
          icxx"$event.Syscall().state == rr::EXITING_SYSCALL;" &&
          icxx"$event.Syscall().number == $number;";
            break
     end
   end
end

sock_fd = 0
# Run until connection to the wine server
while true
  step_to_syscall(42);
  regs = icxx"$(current_task(current_session(timeline)))->regs();";
  path = String(Gallium.load(current_task(current_session(timeline)),
    Gallium.RemotePtr{UInt8}(Gallium.get_dwarf(regs, :rsi)+2), 108))
  sock_fd = Gallium.get_dwarf(regs, :rdi)
  @show path[1:6]
  if String(path[1:6]) ==  "socket"
      break
  end
end

# Run until we receive the request fd
while true
  step_to_syscall(47);
  regs = icxx"$(current_task(current_session(timeline)))->regs();";
  if sock_fd == Gallium.get_dwarf(regs, :rdi)
      break
  end
end
step_to_syscall(72);
request_fd = Gallium.get_dwarf(icxx"$(current_task(current_session(timeline)))->regs();", :rdi)

using Wine
while true
  step_to_syscall(1);
  regs = icxx"$(current_task(current_session(timeline)))->regs();";
  if request_fd == Gallium.get_dwarf(regs, :rdi)
      try
        global req
        regs = icxx"$(current_task(current_session(timeline)))->regs();";
        req = Wine.load_request(current_task(current_session(timeline)), Gallium.get_dwarf(regs, :rsi))
      catch err
        @show err
        isa(err, ErrorException) && err.msg == "unknown request" && continue
      end
      break
  end
end

regs = icxx"$(current_task(current_session(timeline)))->regs();";
req = Wine.load_request(current_task(current_session(timeline)), Gallium.get_dwarf(regs, :rsi))
peb = req.entry

function Gallium.Win64DyldMoudles.get_peb_addr(vm)
    return peb
end

RR.continue!(timeline)

function reload_modules()
    global modules
    imageh = RR.read_exe(current_session(timeline))
    modules = Gallium.GlibcDyldModules.load_library_map(current_task(current_session(timeline)), imageh)
    win64modules = Gallium.Win64DyldMoudles.load_library_map(current_task(current_session(timeline)))
    merge!(modules, win64modules)
    nothing
end


text:
using ForwardDiff
using BenchmarkTools
using StaticArrays
using FiniteDiff
using LinearAlgebra
using SparseArrays

push!(LOAD_PATH, "./src")
using CommonUtils
using ExplicitJacobians
push!(LOAD_PATH, "./examples/EntropyStableEuler")
using EntropyStableEuler

F(x,y) = (x^2+x*y+y^2)/6
# F(x,y) = 1.0/(x+y)
# F(x,y) = logmean(x,y)
dF(x,y) = ForwardDiff.derivative(y->F(x,y),y)

K = 3
Q = randn(K,K)
Q = Q-transpose(Q)

function ff!(out,u,F,QTr)
    for i = 1:size(QTr,1)
        val = zero(eltype(u))
        ui = u[i]
        for j = 1:size(QTr,2)
            uj = u[j]
            val += QTr[j,i]*F(ui,uj)
        end
        out[i] = val
    end
end

# explicit Jacobian
function jac!(dfdu,u,dF,QTr)
    for i = 1:size(QTr,1)
        ui = u[i]
        for j = 1:size(QTr,2)
            uj = u[j]
            dfdu[i,j] = QTr[j,i]*dF(ui,uj)
        end
    end
    sum_df = sum(dfdu,dims=1)
    for i = 1:size(QTr,1)
        dfdu[i,i] += sum_df[i]
    end
end

u = 1 .+ rand(K)
cache = FiniteDiff.JacobianCache(u)
output = zeros(size(Q))
f!(out,u) = ff!(out,u,F,Q)
function ff(u,Q,F)
    uL,uR = meshgrid(u)
    return sum(Q.*F.(uL,uR),dims=2)
end
f(u) = ff(u,Q,F)

# @btime ForwardDiff.jacobian($f,$u) # AD jacobian
# @btime FiniteDiff.finite_difference_jacobian!($output,$f!,$u,$cache) # FiniteDiff.jl
# @btime jac!($output,$u,$dF,$Q) # explicit Jacobian
# out = similar(u); @btime f!($out,$u) # flux eval

function LF(uL,uR)
        return (@. .5*max(abs(uL),abs(uR))*(uL-uR))
        # return (@. (uL-uR))
end
B = randn(K,K)
B = B'*B
# B = [0 1;1 0]
d!(out,u) = ff!(out,u,LF,B)
FiniteDiff.finite_difference_jacobian!(output,d!,u,cache)

uR,uL = meshgrid(u)
dxLF(uL,uR) = ForwardDiff.derivative(uL->LF(uL,uR),uL)
dyLF(uL,uR) = ForwardDiff.derivative(uR->LF(uL,uR),uR)
jacx = -B.*transpose(dxLF.(uL,uR)) + diagm(vec(sum(B.* transpose(dxLF.(uL,uR)),dims=1)))
jacy = B.*dyLF.(uL,uR) - diagm(vec(sum(B.*dyLF.(uL,uR),dims=1)))

dxLFJ(uL,uR) = ForwardDiff.jacobian(uL->LF(uL,uR),uL)
dyLFJ(uL,uR) = ForwardDiff.jacobian(uR->LF(uL,uR),uR)
jacx2 = Matrix(hadamard_jacobian(sparse(B), dxLFJ, [u]))
jacy2 = Matrix(hadamard_jacobian(sparse(B), dyLFJ, [u]))

@show norm(output - jacx), norm(jacx - jacx2)
@show norm(output - jacy), norm(jacy - jacy2)

# F(x,y)  = SVector{2}(x[1]*y[2],y[1]+x[2])
# dF(x,y) = ForwardDiff.jacobian(y->F(x,y),y)
# uL = @SVector [.1 .+ rand(10000),.1 .+ rand(10000)]
# uR = @SVector [.1 .+ rand(10000),.1 .+ rand(10000)]
#
# function feval(uL,uR,F)
#     for i = 1:length(uL)
#         uLi = getindex.(uL,i)
#         uRi = getindex.(uR,i)
#         F(uLi,uRi)
#     end
# end
#
# @btime feval($uL,$uR,$F)
# @btime feval($uL,$uR,$dF)

max_stars_repo_name: jlchan/ESDG-CNS

max_issues_repo_name: jlchan/ESDG-CNS

max_forks_repo_name: jlchan/ESDG-CNS


text:
module Wicked

### Imports ###

using VT100
import VT100: Cell

### Exports ###

export Panel, PanelIO

### Types ###

"""
    Panel

A representation of a drawable terminal surface.
"""
struct Panel
  data::Matrix{Cell}
end

"""
    PanelIO <: IO

An IO object that can parse terminal input (including control codes) and store
them represented as a `Matrix` of `VT100.Cell`s. See `VT100.ScreenEmulator`
for details on the backend's capabilities.

Note: `PanelIO` epresents sizes and indices in (height, width) format, whereas
`VT100.ScreenEmulator` takes sizes as (width, height).
"""
struct PanelIO <: IO
  screen::ScreenEmulator
end
PanelIO(height, width) = PanelIO(ScreenEmulator(width, height))
function Panel(io::PanelIO)
  panel = Panel(undef, size(io)...)
  for h in 1:size(io)[1]
    for w in 1:size(io)[2]
      panel[h,w] = io[h,w]
    end
  end
  panel
end

### Methods ###

Base.size(panel::Panel) = size(panel.data)
Base.size(io::PanelIO) =
  (io.screen.ViewPortSize.height, io.screen.ViewPortSize.width)

function Base.getindex(io::PanelIO, row::Int, col::Int)
  se = io.screen
  height, width = size(io)
  @assert 1 <= row <= height "Invalid row index: $row"
  @assert 1 <= col <= width "Invalid col index: $col"
  if length(se.lines.data) >= row
    if length(se.lines[row].data) >= col
      return se.lines[row][col]
    end
  end
  return Cell(' ')
end

function Base.print(io::PanelIO, x::Union{SubString{String}, String})
  iob = IOBuffer(modified_parse(string(x)))
  parseall!(io.screen, iob)
end
# TODO: Just manually copy cells into io's buffer
function Base.print(io::PanelIO, panel::Panel)
  iob = IOBuffer()
  for cell in panel.data
    cell_dump!(iob, cell)
  end
  print(io, String(take!(iob)))
end
function Base.write(io::PanelIO, x::UInt8)
  # TODO: Manually handle the case where x is '\n'
  parseall!(io.screen, IOBuffer(modified_parse(string(Char(x)))))
  return 1
end

# TODO: This probably shouldn't reset the cursor to top-left
# TODO: since this might be used for display in the REPL
function Base.show(io::IO, panel::Panel)
  print(io, "\e[1;1H")
  iob = IOBuffer()
  panel_dump!(iob, panel)
  print(io, String(take!(iob)))
end

"Appends a terminal-parseable representation of `panel` into `iob`"
function panel_dump!(iob::IOBuffer, panel::Panel)
  height, width = size(panel.data)
  for h in 1:height
    print(iob, "\e[s")
    for w in 1:width
      cell_dump!(iob, panel.data[h,w])
      print(iob, "\e[0m")
    end
    print(iob, "\e[u\e[B")
  end
end

"Appends a terminal-parseable representation of `cell` into `iob`"
function cell_dump!(iob::IOBuffer, cell::Cell)
  # FIXME: More attrs
  if cell.attrs & VT100.Attributes.Blink > 0
    print(iob, "\e[5m")
  end
  if cell.flags & VT100.Flags.FG_IS_RGB > 0
    r = Int(floor(cell.fg_rgb.r*255))
    b = Int(floor(cell.fg_rgb.g*255))
    g = Int(floor(cell.fg_rgb.b*255))
    print(iob, "\e[38;2;$(r);$(g);$(b)m")
  elseif cell.flags & VT100.Flags.FG_IS_256 > 0
    if 0 <= cell.fg <= 7
      print(iob, "\e[3$(string(cell.fg))m")
    elseif 8 <= cell.fg <= 15
      print(iob, "\e[9$(string(cell.fg-8))m")
    else
      print(iob, "\e[38;5;$(string(cell.fg))m")
    end
  else
    print(iob, "\e[39m")
  end

  if cell.flags & VT100.Flags.BG_IS_RGB > 0
    r = Int(floor(cell.bg_rgb.r*255))
    b = Int(floor(cell.bg_rgb.g*255))
    g = Int(floor(cell.bg_rgb.b*255))
    print(iob, "\e[48;2;$(r);$(g);$(b)m")
  elseif cell.flags & VT100.Flags.BG_IS_256 > 0
    if 0 <= cell.bg <= 7
      print(iob, "\e[4$(string(cell.bg))m")
    elseif 8 <= cell.bg <= 15
      print(iob, "\e[10$(string(cell.bg-8))m")
    else
      print(iob, "\e[48;5;$(string(cell.bg))m")
    end
  else
    print(iob, "\e[49m")
  end

  print(iob, cell.content)
end

"""Returns a new `Panel` resized to the desired height and width. If a
dimension is extended, fills those `Cell`s with `fill_value`."""
function resize_panel(panel::Panel, new_size::Tuple{Int,Int},
    fill_value=Cell(' '))
  old_size = size(panel)
  old_size == new_size && return panel
  new_panel = fill(fill_value, new_size)
  shared_region = CartesianIndices(min.(old_size, new_size))
  copy!(new_panel, shared_region, panel, shared_region)
end

"Prepends carriage returns to newlines in a string."
# TODO: Add method for ::IO argument
function modified_parse(str::String)
  retstr = ""
  for char in str
    if char == '\n'
      retstr *= '\r'
    end
    retstr *= char
  end
  retstr
end

end # module


text:
"""
Adaptors for using tensor-network codes within the
[Qecsim](https://github.com/qecsim/Qecsim.jl) framework.

Code, [`QecsimTNCode`](@ref), and decoder, [`QecsimTNDecoder`](@ref), implementations are
provided to be used with Qecsim error models for QEC simulations using the Qecsim `App`
module.
"""
module QecsimAdaptors

#imports
using ..TensorNetworkCodes: TensorNetworkCode, num_qubits
using ..TNDecode: basic_contract, mps_contract, tn_decode
using Qecsim

#exports
export QecsimTNCode, QecsimTNDecoder

"""
    QecsimTNCode <: Qecsim.StabilizerCode

    QecsimTNCode(code::TensorNetworkCode; distance=missing, label=nothing)

Qecsim stabilizer code implementation based on a tensor-network code.

Pauli operators such as `stabilizers` and `logicals` are converted to Qecsim's format on
construction. These are available via Qecsim API methods of the same name. The Qecsim `nkd`
method returns calculated values for `n` (number of physical qubits) and `k` (number of
logical qubits) and the given `distance` for `d`, if provided, otherwise `missing`. The
Qecsim `label` method returns the given `label`, if provided, otherwise the label is
constructed from `[n,k,d]`.

Public fields:

    tn_code::TensorNetworkCode  # the given tensor-network code

# Examples
```jldoctest
julia> using TensorNetworkCodes.QecsimAdaptors

julia> using Qecsim: label, nkd, validate

julia> tn_code = TensorNetworkCode(five_qubit_code());

julia> qs_code = QecsimTNCode(tn_code; distance=3);

julia> label(qs_code)
"QecsimTNCode: [5,1,3]"

julia> nkd(qs_code)
(5, 1, 3)

julia> validate(qs_code)  # no error indicates success
```
"""
struct QecsimTNCode <: StabilizerCode
    _stabilizers::BitMatrix
    _logical_xs::BitMatrix
    _logical_zs::BitMatrix
    _nkd::Tuple{Int,Int,Union{Int,Missing}}
    _label::String
    tn_code::TensorNetworkCode
end
Model.stabilizers(code::QecsimTNCode) = code._stabilizers
Model.logical_xs(code::QecsimTNCode) = code._logical_xs
Model.logical_zs(code::QecsimTNCode) = code._logical_zs
Model.nkd(code::QecsimTNCode) = code._nkd
Model.label(code::QecsimTNCode) = code._label
function QecsimTNCode(code::TensorNetworkCode; distance=missing, label=nothing)
    # derived defaults
    dnkd = (num_qubits(code), length(code.logicals)÷2, distance)
    dlabel = isnothing(label) ? "QecsimTNCode: [$(dnkd[1]),$(dnkd[2]),$(dnkd[3])]" : label
    return QecsimTNCode(
        _tnpauli_to_bsf(code.stabilizers),
        _tnpauli_to_bsf(code.logicals[1:2:end]),
        _tnpauli_to_bsf(code.logicals[2:2:end]),
        dnkd,
        dlabel,
        code
    )
end

"""
    _tnpauli_to_bsf(pauli::AbstractVector{Int}) -> BitVector
    _tnpauli_to_bsf(paulis::AbstractVector{<:AbstractVector{Int}}) -> BitMatrix

Return the binary symplectic representation of Pauli given in integer vector representation.
"""
function _tnpauli_to_bsf(p::AbstractVector{Int})
    return vcat((x-> x ∈ (1, 2)).(p), (x-> x ∈ (2, 3)).(p))
end
function _tnpauli_to_bsf(ps::AbstractVector{<:AbstractVector{Int}})
    return transpose(reduce(hcat, _tnpauli_to_bsf.(ps)))
end
"""
    _bsf_to_tnpauli(b::AbstractVector{Bool}) -> Vector{Int}
    _bsf_to_tnpauli(bs::AbstractMatrix{Bool}) -> Vector{Vector{Int}}

Return the integer vector representation of Pauli given in binary symplectic representation.
"""
function _bsf_to_tnpauli(b::AbstractVector{Bool})
    return map(x -> (0,1,3,2)[x+1], b[1:end÷2] + 2b[end÷2+1:end])
end
function _bsf_to_tnpauli(bs::AbstractMatrix{Bool})
    return _bsf_to_tnpauli.(b for b in eachrow(bs))
end

@doc raw"""
    QecsimTNDecoder <: Decoder

    QecsimTNDecoder(chi::Union{Nothing,Integer}=nothing)

Qesim decoder implementation based on the [TNDecode module](@ref). A null value of `chi`
corresponds to exact contraction, otherwise chi defines the bond dimension used in MPS/MPO
contraction.

An `ArgumentError` is thrown if chi is not null or positive. The Qecsim `label` method
returns a label constructed from `chi`. The Qecsim `decode` method expects the keyword
argument `p` for error probability, which defaults to `0.1`, and returns a `DecodeResult`
object with `custom_values` containing the success probability. The success probability is
defined as the ratio of the probability of the mostly-likely logical coset to the sum of the
probabilities of all logical cosets, for example:

    DecodeResult.custom_values = [[0.92]] # e.g. where 0.92 is the success probability

The success probability is added to `custom_values` as a vector so that as Qecsim
`App.qec_run` aggregates success probabilities over many runs it extends the vector rather
than simply summing over values.

!!! note
    Currently the [TNDecode module](@ref) only supports codes that can be laid out on a
    square lattice.

# Examples
```jldoctest; filter = r"0\.92498\d*"
julia> using TensorNetworkCodes.QecsimAdaptors

julia> using Qecsim, Qecsim.GenericModels

julia> using Random:MersenneTwister  # use RNG for reproducible example

julia> code = QecsimTNCode(rotated_surface_code(3); distance=3);

julia> error_model = DepolarizingErrorModel();

julia> decoder = QecsimTNDecoder(4);

julia> label(decoder)
"QecsimTNDecoder (chi=4)"

julia> result = qec_run_once(code, error_model, decoder, 0.1, MersenneTwister(11))
RunResult{Vector{Vector{Float64}}}(true, 1, Bool[0, 0], [[0.9249813981321253]])

julia> success_flag = result.success
true

julia> success_probability = result.custom_values[1][1]
0.9249813981321253
```
"""
struct QecsimTNDecoder <: Decoder
    chi::Union{Nothing,Int}
    function QecsimTNDecoder(chi=nothing)
        if !(isnothing(chi) || chi > 0)
            throw(ArgumentError("QecsimTNDecoder supports `nothing` or positive `chi`."))
        end
        return new(chi)
    end
end
function Model.label(decoder::QecsimTNDecoder)
    chi = decoder.chi
    return isnothing(chi) ? "QecsimTNDecoder" : "QecsimTNDecoder (chi=$(chi))"
end
function Model.decode(
    decoder::QecsimTNDecoder, code::QecsimTNCode, syndrome::AbstractVector{Bool};
    p=0.1, kwargs...
)
    chi = decoder.chi
    tn_syndrome::Vector{Int} = syndrome  # promote to Vector{Int}
    tn_contract_fn = isnothing(chi) ? basic_contract() : mps_contract(chi)
    tn_recovery, predicted_success_rate = tn_decode(
        code.tn_code, tn_syndrome, p;
        contract_fn=tn_contract_fn
    )
    recovery = _tnpauli_to_bsf(tn_recovery)
    # add predicted_success_rate as vector so all values retained (not summed) over runs
    return DecodeResult(recovery=recovery, custom_values=[[predicted_success_rate]])
end

end


text:
export Sheep, Wolf, Grass

mutable struct Sheep <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    energy::Float64
end

mutable struct Wolf <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    energy::Float64
end

mutable struct Grass <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    fully_grown::Bool
    countdown::Int
end

"""
``` julia
predator_prey(;
    n_sheep = 100,
    n_wolves = 50,
    dims = (20, 20),
    regrowth_time = 30,
    Δenergy_sheep = 4,
    Δenergy_wolf = 20,
    sheep_reproduce = 0.04,
    wolf_reproduce = 0.05,
)
```
Same as in [Model of predator-prey dynamics](@ref).

To access the `Sheep`, `Wolf` and `Grass` types, simply call
``` julia
using Agents.Models: Sheep, Wolf, Grass
```
"""
function predator_prey(;
    n_sheep = 100,
    n_wolves = 50,
    dims = (20, 20),
    regrowth_time = 30,
    Δenergy_sheep = 4,
    Δenergy_wolf = 20,
    sheep_reproduce = 0.04,
    wolf_reproduce = 0.05,
)
    space = GridSpace(dims, periodic = false)
    properties = Dict(
        :Δenergy_wolf => Δenergy_wolf,
        :Δenergy_sheep => Δenergy_sheep,
        :regrowth_time => regrowth_time,
        :sheep_reproduce => sheep_reproduce,
        :wolf_reproduce => wolf_reproduce,
    )
    model = ABM(
        Union{Sheep,Wolf,Grass},
        space,
        scheduler = by_type(true, true),
        warn = false,
        properties = properties,
    )
    id = 0
    for _ in 1:n_sheep
        id += 1
        energy = rand(1:(Δenergy_sheep*2)) - 1
        sheep = Sheep(id, (0, 0), energy)
        add_agent!(sheep, model)
    end
    for _ in 1:n_wolves
        id += 1
        energy = rand(1:(Δenergy_wolf*2)) - 1
        wolf = Wolf(id, (0, 0), energy)
        add_agent!(wolf, model)
    end
    for p in positions(model)
        id += 1
        fully_grown = rand(Bool)
        countdown = fully_grown ? regrowth_time : rand(1:regrowth_time) - 1
        grass = Grass(id, (0, 0), fully_grown, countdown)
        add_agent!(grass, p, model)
    end
    return model, predator_prey_agent_step!, dummystep
end

function predator_prey_agent_step!(sheep::Sheep, model)
    move!(sheep, model)
    sheep.energy -= 1
    agents = collect(agents_in_position(sheep.pos, model))
    dinner = filter!(x -> isa(x, Grass), agents)
    eat!(sheep, dinner, model)
    if sheep.energy < 0
        kill_agent!(sheep, model)
        return
    end
    if rand() <= model.sheep_reproduce
        reproduce!(sheep, model)
    end
end

function predator_prey_agent_step!(wolf::Wolf, model)
    move!(wolf, model)
    wolf.energy -= 1
    agents = collect(agents_in_position(wolf.pos, model))
    dinner = filter!(x -> isa(x, Sheep), agents)
    eat!(wolf, dinner, model)
    if wolf.energy < 0
        kill_agent!(wolf, model)
        return
    end
    if rand() <= model.wolf_reproduce
        reproduce!(wolf, model)
    end
end

function predator_prey_agent_step!(grass::Grass, model)
    if !grass.fully_grown
        if grass.countdown <= 0
            grass.fully_grown = true
            grass.countdown = model.regrowth_time
        else
            grass.countdown -= 1
        end
    end
end

function move!(agent, model)
    neighbors = nearby_positions(agent, model)
    position = rand(collect(neighbors))
    move_agent!(agent, position, model)
end

function eat!(sheep::Sheep, grass_array, model)
    isempty(grass_array) && return
    grass = grass_array[1]
    if grass.fully_grown
        sheep.energy += model.Δenergy_sheep
        grass.fully_grown = false
    end
end

function eat!(wolf::Wolf, sheep, model)
    if !isempty(sheep)
        dinner = rand(sheep)
        kill_agent!(dinner, model)
        wolf.energy += model.Δenergy_wolf
    end
end

function reproduce!(agent, model)
    agent.energy /= 2
    id = nextid(model)
    A = typeof(agent)
    offspring = A(id, agent.pos, agent.energy)
    add_agent_pos!(offspring, model)
    return
end


text:
mutable struct GEK{X,Y,L,U,P,T,M,B,S,R} <: AbstractSurrogate
    x::X
    y::Y
    lb::L
    ub::U
    p::P
    theta::T
    mu::M
    b::B
    sigma::S
    inverse_of_R::R
 end

function _calc_gek_coeffs(x,y,p::Number,theta::Number)
    nd1 = length(y) #2n
    n = length(x)
    R = zeros(eltype(x[1]), nd1, nd1)

    #top left
    @inbounds for i = 1:n
        for j = 1:n
            R[i,j] = exp(-theta*abs(x[i]-x[j])^p)
        end
    end
    #top right
    @inbounds for i = 1:n
        for j = n+1:nd1
            R[i,j] = 2*theta*(x[i] - x[j-n])*exp(-theta*abs(x[i]-x[j-n])^p)
        end
    end
    #bottom left
    @inbounds for i = n+1:nd1
        for j = 1:n
            R[i,j] = -2*theta*(x[i-n] - x[j])*exp(-theta*abs(x[i-n]-x[j])^p)
        end
    end
    #bottom right
    @inbounds for i = n+1:nd1
        for j = n+1:nd1
            R[i,j] = -4*theta*(x[i-n] - x[j-n])^2*exp(-theta*abs(x[i-n]-x[j-n])^p)
        end
    end
    one = ones(eltype(x[1]),nd1,1)
    for i = n+1:nd1
        one[i] = zero(eltype(x[1]))
    end
    one_t = one'
    inverse_of_R = inv(R)
    mu = (one_t*inverse_of_R*y)/(one_t*inverse_of_R*one)
    b = inverse_of_R*(y-one*mu)
    sigma = ((y-one*mu)' * inverse_of_R * (y - one*mu))/n
    mu[1], b, sigma[1],inverse_of_R
end

function std_error_at_point(k::GEK,val::Number)
    phi(z) = exp(-(abs(z))^k.p)
    nd1 = length(k.y)
    n = length(k.x)
    r = zeros(eltype(k.x[1]),nd1,1)
    @inbounds for i = 1:n
        r[i] = phi(val - k.x[i])
    end
    one = ones(eltype(k.x[1]),nd1,1)
    @inbounds for i = n+1:nd1
        one[i] = zero(eltype(k.x[1]))
    end
    one_t = one'
    a = r'*k.inverse_of_R*r
    a = a[1]
    b = one_t*k.inverse_of_R*one
    b = b[1]
    mean_squared_error = k.sigma*(1 - a + (1-a)^2/(b))
    return sqrt(abs(mean_squared_error))
end

function (k::GEK)(val::Number)
    phi = z -> exp(-(abs(z))^k.p)
    n = length(k.x)
    prediction = zero(eltype(k.x[1]))
    @inbounds for i = 1:n
        prediction = prediction + k.b[i]*phi(val-k.x[i])
    end
    prediction = k.mu + prediction
    return prediction
end

function GEK(x,y,lb::Number,ub::Number; p=1.0,theta=1.0)
    if length(x) != length(unique(x))
        println("There exists a repetion in the samples, cannot build Kriging.")
        return
    end
    mu,b,sigma,inverse_of_R = _calc_gek_coeffs(x,y,p,theta)
    return GEK(x,y,lb,ub,p,theta,mu,b,sigma,inverse_of_R)
end

function _calc_gek_coeffs(x,y,p,theta)
    nd = length(y)
    n = length(x)
    d = length(x[1])
    R = zeros(eltype(x[1]), nd, nd)

    #top left
    @inbounds for i = 1:n
        for j = 1:n
            R[i,j] = prod([exp(-theta[l]*(x[i][l] - x[j][l])) for l = 1:d])
        end
    end

    #top right
    @inbounds for i = 1:n
        jr = 1
        for j = n+1:d:nd
            for l = 1:d
                R[i,j+l-1] = +2*theta[l]*(x[i][l]-x[jr][l])*R[i,jr]
            end
            jr = jr+1
        end
    end

    #bottom left
    @inbounds for j = 1:n
        ir = 1
        for i = n+1:d:nd
            for l = 1:d
                R[i+l-1,j] = -2*theta[l]*(x[ir][l]-x[j][l])*R[ir,j]
            end
            ir = ir + 1
        end
    end

    #bottom right
    ir = 1
    @inbounds for i = n+1:d:nd
        for j = n+1:d:nd
            jr = 1
            for l = 1:d
                for k  = 1:d
                    R[i+l-1,j+k-1] = -4*theta[l]*theta[k]*(x[ir][l]-x[jr][l])*(x[ir][k]-x[jr][k])*R[ir,jr]
                end
            end
            jr = jr + 1
        end
        ir = ir + + 1
    end

    one = ones(eltype(x[1]),nd,1)
    for i = n+1:nd
        one[i] = zero(eltype(x[1]))
    end
    one_t = one'
    inverse_of_R = inv(R)
    mu = (one_t*inverse_of_R*y)/(one_t*inverse_of_R*one)
    b = inverse_of_R*(y-one*mu)
    sigma = ((y-one*mu)' * inverse_of_R * (y - one*mu))/n
    mu[1], b, sigma[1],inverse_of_R

end

function std_error_at_point(k::GEK,val)
    nd1 = length(k.y)
    n = length(k.x)
    d = length(k.x[1])
    r = zeros(eltype(k.x[1]),nd1,1)
    @inbounds for i = 1:n
        sum = zero(eltype(k.x[1]))
        for l = 1:d
            sum = sum + k.theta[l]*norm(val[l]-k.x[i][l])^(k.p[l])
        end
        r[i] = exp(-sum)
    end
    one = ones(eltype(k.x[1]),nd1,1)
    @inbounds for i = n+1:nd1
        one[i] = zero(eltype(k.x[1]))
    end
    one_t = one'
    a = r'*k.inverse_of_R*r
    a = a[1]
    b = one_t*k.inverse_of_R*one
    b = b[1]
    mean_squared_error = k.sigma*(1 - a + (1-a)^2/(b))
    return sqrt(abs(mean_squared_error))
end

function (k::GEK)(val)
    d = length(val)
    n = length(k.x)
    return k.mu + sum(k.b[i] * exp( -sum(k.theta[j]*norm(val[j]-k.x[i][j])^k.p[j] for j = 1:d) ) for i = 1:n)
end

function GEK(x,y,lb,ub; p=collect(one.(x[1])),theta=collect(one.(x[1])))
    if length(x) != length(unique(x))
        println("There exists a repetition in the samples, cannot build Kriging.")
        return
    end
    mu,b,sigma,inverse_of_R = _calc_gek_coeffs(x,y,p,theta)
    GEK(x,y,lb,ub,p,theta,mu,b,sigma,inverse_of_R)
end

function add_point!(k::GEK,new_x,new_y)
    if new_x in k.x
        println("Adding a sample that already exists, cannot build Kriging.")
        return
    end
    if (length(new_x) == 1 && length(new_x[1]) == 1) || ( length(new_x) > 1 && length(new_x[1]) == 1 && length(k.theta)>1)
        n = length(k.x)
        k.x = insert!(k.x,n+1,new_x)
        k.y = insert!(k.y,n+1,new_y)
    else
        n = length(k.x)
        k.x = insert!(k.x,n+1,new_x)
        k.y = insert!(k.y,n+1,new_y)
    end
    k.mu,k.b,k.sigma,k.inverse_of_R = _calc_gek_coeffs(k.x,k.y,k.p,k.theta)
    nothing
end


text:
"""
Inner product (element-wise product) between two matrices
"""
function innerprod(A::AbstractMatrix{Float64},B::AbstractMatrix{Float64})::Float64
    if size(A) != size(B)
        throw(DimensionMismatch("size of A not equal to size of B"))
    end
    return sum(A.*B)/size(A)[1]
end

"""
Return a copy of the off-diagonal part of matrix
"""
function get_offdiag(A::AbstractMatrix{Float64})::AbstractMatrix{Float64}
    return A-Diagonal(A)
end

"""
Create center finite difference matrix
"""
function get_cfd_mat(ord::Int,n::Int)::SparseMatrixCSC{Float64,Int}
    if ord == 2
        arr = circshift(vcat(Vector{Float64}([1,-2,1]),zeros(Float64,n-3)),-1)
    elseif ord == 4
        arr = circshift(vcat(Vector{Float64}([-1/12,4/3,-5/2,4/3,-1/12]),zeros(Float64,n-5)),-2)
    elseif ord == 8
        arr = circshift(vcat(Vector{Float64}([1/90,-3/20,3/2,-49/18,3/2,-3/20,1/90]),zeros(Float64,n-7)),-3)
    else
        throw("unimplemented")
    end

    A = Matrix{Float64}(undef,n,n)
    for i in 1:n
        A[i,:] = arr
        arr = circshift(arr,1)
    end

    return A
end

function get_stagger_grids(ord::Int,n::Int,a::Vector{Float64})::Matrix{Float64}
    if length(a) != n
        throw("size of vector a not equal to n")
    end

    if ord == 2
        coefs = Float64[  1  0;
                         -1 -1;
                          0  1;]

        h = 1.0/n
        A = zeros(Float64,n,n)

        for i in 1:n
            arr = circshift(a,2-i)[1:2]
            A[i,:] = circshift(vcat(coefs*arr,zeros(Float64,n-3)),i-2)
        end

    elseif ord == 4
        coefs = Float64[ 1/576      0      0      0;
                         -3/64  -3/64      0      0;
                          3/64  81/64   3/64      0;
                        -1/576 -81/64 -81/64 -1/576;
                             0   3/64  81/64   3/64;
                             0      0  -3/64  -3/64;
                             0      0      0  1/576;]

        h = 1.0/n
        A = zeros(Float64,n,n)

        for i in 1:n
            arr = circshift(a,3-i)[1:4]
            A[i,:] = circshift(vcat(coefs*arr,zeros(Float64,n-7)),i-4)
        end

    elseif ord == 6
        coefs = Float64[ 9/409600           0          0          0           0        0;
                         -5/16384    -5/16384          0          0           0        0;
                          45/8192  625/147456    45/8192          0           0        0;
                         -45/8192   -625/8192  -625/8192   -45/8192           0        0;
                          5/16384    625/8192  5625/4096   625/8192     5/16384        0;
                        -9/409600 -625/147456 -5625/4096 -5625/4096 -625/147456 -9/409600;
                                0     5/16384   625/8192  5625/4096    625/8192   5/16384;
                                0           0   -45/8192  -625/8192   -625/8192  -45/8192;
                                0           0          0    45/8192  625/147456   45/8192;
                                0           0          0          0    -5/16384  -5/16384;
                                0           0          0          0           0  9/409600;]

        h = 1.0/n
        A = zeros(Float64,n,n)

        for i in 1:n
            arr = circshift(a,4-i)[1:6]
            A[i,:] = circshift(vcat(coefs*arr,zeros(Float64,n-11)),i-6)
        end
    else
        throw("unimplemented")
    end
    return A
end

"""
Check if a matrix is strictly diagonally dominant
"""
function is_diagally_dominant(A::AbstractMatrix{Float64})
    diff = abs.(diag(A)).-collect(Iterators.flatten(sum(abs.(get_offdiag(A)),dims=1)))
    return minimum(diff)
end

"""
Get 2D meshgrid
"""
function meshgrid(x,y)
    return (reshape(repeat(x,outer=length(y)),(length(x),length(y))),
        reshape(repeat(y,inner=length(x)),(length(x),length(y))))
end

"""
Get index for repartition of girds
"""
function get_repartiton_index(n::Int)::Vector{Int}
    idx = transpose(reshape(collect([1:n*n;]),(n,n)))
    return vcat(collect(Iterators.flatten(idx[1:n÷2,1:n÷2])),
        collect(Iterators.flatten(idx[n÷2+1:n,1:n÷2])),
        collect(Iterators.flatten(idx[1:n÷2,n÷2+1:n])),
        collect(Iterators.flatten(idx[n÷2+1:n,n÷2+1:n]))
    )
end

max_stars_repo_name: xinshuoyang/FactNSF

max_issues_repo_name: xinshuoyang/FactNSF

max_forks_repo_name: xinshuoyang/FactNSF


text:
get_theo_A(g::Grad; off_val=0) = begin
	A = g.A
	if length(A) == 1
		aux = [off_val; 0; A; A; 0]
	else
		aux = [off_val; 0; A; 0]
	end
	if sum(abs.(g.A)) == 0
		aux *= off_val
	end
	aux
end
get_theo_t(g::Grad) = begin
	NT, T, NA = length(g.T), g.T, length(g.A)
	if sum(T) != 0
		if NA == 1 && NT == 1
			trf = [0; T]
		elseif NA>1 && NT == 1
			trf = cumsum([0; T/(NA-1).*ones(NA-1)])
		elseif NA>1 && NT>1
			trf = cumsum([0; T.*ones(NT)])
		end
		t = g.delay .+ g.rise .+ trf
	else
		t = [g.delay+g.rise; g.delay+g.rise]
	end
	aux = [0; g.delay; t; g.delay+g.rise+sum(g.T)+g.fall]
	aux
end
get_theo_A(r::RF; off_val=0, max_rf_samples=Inf) = begin
	A = [r.A r.A]'; A = A[:]
	if sum(abs.(r.A)) == 0
		aux = off_val * [1; 1; ones(length(A)); 1]
	else
		#Calculating frequency waveform
		if length(r.A) != 1
			tf = range(0,sum(r.T),length(r.A))
		else
			tf = [0]
		end
		freq = exp.(1im*2π*r.Δf.*tf)
		freq = transpose([freq freq])[:]
		A = A .* freq
		A[abs.(A).<=1e-14] .= 0 #Remove small values
		#Output
		aux = [off_val; 0; A; 0]
	end
	#Subsample
	if length(aux) > max_rf_samples
		n = floor(Int, length(aux) / max_rf_samples)
		aux[1:n:end]
	else
		aux
	end
end
get_theo_t(r::RF; max_rf_samples=Inf) = begin
	if sum(r.T) != 0
		NA, T, NT = length(r.A), r.T, length(r.T)
		dT = T / NA * NT
		trf = cumsum([0; dT.*ones(NA)])
		t = r.delay .+ trf
	else
		t = [r.delay; r.delay]
	end
	t = [t t]'; t = t[:]
	t = [0; t]
	#Subsample
	if length(t) > max_rf_samples
		n = floor(Int, length(t) / max_rf_samples)
		t[1:n:end]
	else
		t
	end
end
get_theo_A(d::ADC; off_val=0) = begin
	N = [1 1]'; N = N[:]
	if d.N == 0
		aux = off_val * [1; 1; ones(length(N)); 1]
	else
		aux = [off_val; 0; N; 0]
	end
	aux
end
get_theo_t(d::ADC) = begin
	t = [d.delay; d.delay+d.T]
	t = [t t]'; t = t[1:end-1]
	[0; t; d.delay+d.T]
end
get_theo_Gi(seq, idx) = begin
	ΔT, N = durs(seq), length(seq)
	T0 = cumsum([0; ΔT],dims=1)
	t = vcat([get_theo_t(seq.GR[idx,i]) .+ T0[i] for i=1:N]...)
	G = vcat([get_theo_A(seq.GR[idx,i]) for i=1:N]...)
	#Removing duplicated points
	#TODO: do this properly. As it is now it generates a bug for slew rates that are too high
	Interpolations.deduplicate_knots!(t; move_knots=true)
	(t,G)
end